diff --git a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
index b0f5f48513a951dc5742bdd96d1a3a8a36775689..a2491996b6bdf210fc53967f8dc5abf92312c369 100644
--- a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
+++ b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
@@ -6289,80 +6289,6 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
}
}
- //2.a - confirm primary track hypothesis by TRD
- if (nTrdHits > 0) {
- for (Int_t i = 0; i < nTofHits; i++) {
- Int_t j = iStsMin[i][1]; // index of STS hit in second plane
- if (j < 0) continue; // no STS assigned
- CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
- if (NULL == pTofHit) continue;
- if (pTofHit->GetZ() == 0) continue; // skip fake beam counter
- CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
- Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
- Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
- Double_t dDz = pTofHit->GetZ() - pStsHit->GetZ();
- LOG(DEBUG) << "Check for TRD hits between STS " << j << " and TOF " << i;
-
- for (Int_t l = 0; l < nTrdHits; l++) {
- CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(l);
- // calculate expected position in Trd layer
- Double_t dXexp =
- pStsHit->GetX() + dDx * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
- Double_t dYexp =
- pStsHit->GetY() + dDy * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
- Double_t dDtrans =
- TMath::Sqrt(TMath::Power(dXexp - pTrdHit->GetX(), 2)
- + TMath::Power(dYexp - pTrdHit->GetY(), 2));
- UInt_t iTrdLayer = CbmTrdAddress::GetLayerId(pTrdHit->GetAddress());
- LOG(DEBUG) << "Inspect TRD hit " << l << " in "
- << Form("Module 0x%08x, layer %d",
- pTrdHit->GetAddress(),
- CbmTrdAddress::GetLayerId(pTrdHit->GetAddress()))
- << " at z= " << pTrdHit->GetZ() << " dD = " << dDtrans
- << " < " << fdDistTRD;
- fhDTRDprim->Fill(dDtrans);
- if (dDtrans < fdDistTRD
- && dDtrans
- < dTrdDistMin
- [i]
- [iTrdLayer]) { // check if acceptable and take best match
- Int_t iMul = iTRD[i].size();
- if (dTrdDistMin[i][iTrdLayer] < 1.E3) { // modify previous entry
- //find old entry in vector
- Int_t ll = 0;
- for (; ll < iMul; ll++)
- if (CbmTrdAddress::GetLayerId(
- ((CbmTrdHit*) fTrdHits->At(iTRD[i][ll]))->GetAddress())
- == iTrdLayer)
- break;
- iTRD[i][ll] = l;
- } else { //add hit
- dTrdDistMin[i][iTrdLayer] = dDtrans;
- iTRD[i].resize(iMul + 1);
- iTRD[i][iMul] = l;
- }
- LOG(DEBUG) << "assign TrdHit " << l << " to TofHit " << i
- << " in layer " << iTrdLayer << " with d = " << dDtrans
- << ", TrdMul" << iMul
- << ", dEdx = " << pTrdHit->GetELoss();
- }
- }
- }
- //2.b - monitor TRD dEdx
- for (Int_t i = 0; i < nTofHits; i++) {
- Int_t iMul = iTRD[i].size();
- if (iMul > 0) {
- Double_t ddEdx = 0.;
- for (Int_t l = 0; l < iMul; l++) {
- CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(iTRD[i][l]);
- ddEdx += pTrdHit->GetELoss();
- }
- ddEdx /= (Double_t) iMul;
- fhdEdxMul->Fill((Double_t) iMul, ddEdx);
- }
- }
- }
-
//3. find secondary proton candidate
for (Int_t i = 0; i < nTofHits; i++) {
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
@@ -6566,97 +6492,6 @@ void CbmHadronAnalysis::ReconstructSecondaries() {
DX[i].SetXYZ(dDx, dDy, dDz);
pion_cand++;
}
- //4.a - confirm secondary track hypothesis by TRD
- if (nTrdHits > 0) {
- Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
- Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
- Double_t dDz = pTofHit->GetZ() - pStsHit->GetZ();
- for (Int_t l = 0; l < nTrdHits; l++) {
- CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(l);
- // calculate expected position in Trd layer
- Double_t dXexp =
- pStsHit->GetX() + dDx * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
- Double_t dYexp =
- pStsHit->GetY() + dDy * (pTrdHit->GetZ() - pStsHit->GetZ()) / dDz;
- Double_t dDtrans =
- TMath::Sqrt(TMath::Power(dXexp - pTrdHit->GetX(), 2)
- + TMath::Power(dYexp - pTrdHit->GetY(), 2));
- UInt_t iTrdLayer = CbmTrdAddress::GetLayerId(pTrdHit->GetAddress());
- fhDTRDsec->Fill(dDtrans);
- LOG(DEBUG) << "Inspect sec. TRD hit " << l << " in "
- << Form("Module 0x%08x, layer %d",
- pTrdHit->GetAddress(),
- CbmTrdAddress::GetLayerId(pTrdHit->GetAddress()))
- << " at z= " << pTrdHit->GetZ() << " dD = " << dDtrans
- << " < " << fdDistTRD;
- if (
- dDtrans < fdDistTRD
- && dDtrans
- < dTrdDistMin
- [i]
- [iTrdLayer]) { // check if acceptable and take best match
- Int_t iMul = iTRD[i].size();
- if (dTrdDistMin[i][iTrdLayer] < 1.E3) { // modify previous entry
- //find old entry in vector
- Int_t ll = 0;
- for (; ll < iMul; ll++)
- if (CbmTrdAddress::GetLayerId(
- ((CbmTrdHit*) fTrdHits->At(iTRD[i][ll]))->GetAddress())
- == iTrdLayer)
- break;
- iTRD[i][ll] = l;
- } else { //add hit
- dTrdDistMin[i][iTrdLayer] = dDtrans;
- iTRD[i].resize(iMul + 1);
- iTRD[i][iMul] = l;
- }
- LOG(DEBUG) << "assign TrdHit " << l << " to TofHit " << i
- << " in layer " << iTrdLayer << " with d = " << dDtrans
- << ", TrdMul" << iMul
- << ", dEdx = " << pTrdHit->GetELoss();
- }
- }
- //4.b - monitor TRD dEdx
- Int_t iMul = iTRD[i].size();
- if (iMul > 0) {
- Double_t ddEdx = 0.;
- for (Int_t l = 0; l < iMul; l++) {
- CbmTrdHit* pTrdHit = (CbmTrdHit*) fTrdHits->At(iTRD[i][l]);
- ddEdx += pTrdHit->GetELoss();
- }
- ddEdx /= (Double_t) iMul;
- fhdEdxMulsec->Fill((Double_t) iMul, ddEdx);
- }
- } // (nTrdHits > 0) end
- if ((Double_t) iTRD[i].size() >= fdTRDHmulMin) {
- Double_t dDx = pStsHit->GetX() - pSts2Hit->GetX();
- Double_t dDy = pStsHit->GetY() - pSts2Hit->GetY();
- Double_t dDz = pStsHit->GetZ() - pSts2Hit->GetZ();
- Double_t dDd = TMath::Sqrt(dDx * dDx + dDy * dDy + dDz * dDz);
- Double_t vel =
- pTofHit->GetR()
- / pTofHit->GetTime(); // approximation, ignoring decay kinematics
- Double_t bet = vel / clight;
- if (bet > 0.9999) continue; //bet=0.9999;
- Double_t m = secMass[0]; // assume pion
- Double_t pmag =
- m * bet / TMath::Sqrt(1. - bet * bet); // natural units
- Double_t pz = pmag * dDz / dDd;
- Double_t px = pmag * dDx / dDd;
- Double_t py = pmag * dDy / dDd;
- Double_t E = TMath::Sqrt(pmag * pmag + m * m);
- P[i].SetPxPyPzE(px, py, pz, E);
- X[i].SetXYZT(pTofHit->GetX(),
- pTofHit->GetY(),
- pTofHit->GetZ(),
- pTofHit->GetTime());
- LOG(DEBUG) << "Init pion LV at ind " << i << " with beta = " << bet
- << ", minv = " << P[i].M() << ", tof " << X[i].T()
- << ", TRDHmul " << iTRD[i].size();
- X0[i].SetXYZ(pSts2Hit->GetX(), pSts2Hit->GetY(), pSts2Hit->GetZ());
- DX[i].SetXYZ(dDx, dDy, dDz);
- pion_cand++;
- }
}
} //if( dStsDistMin[i] > dDistPrimLim) { // Sts hit not in the primary class
} //for (Int_t i=0; i<nTofHits; i++)