#include "PndDetectorList.h" #include "PndPidCorrelator.h" #include "PndPidCandidate.h" #include "PndMCTrack.h" #include "PndTrack.h" #include "PndTrackID.h" #include "PndFtofHit.h" #include "PndFtofPoint.h" #include "FairTrackParH.h" #include "FairMCApplication.h" #include "FairRunAna.h" #include "FairRootManager.h" #include "FairRuntimeDb.h" Bool_t PndPidCorrelator::GetFtofInfo(FairTrackParH* helix, PndPidCandidate* pidCand) { if(!helix) { std::cerr << " PndPidCorrelator::GetFtofInfo: FairTrackParH NULL pointer parameter."< PndPidCorrelator::GetFtofInfo: pidCand NULL pointer parameter."<GetZ() < fCorrPar->GetZLastPlane()) { if (fVerbose>0) std::cout << "-W- PndPidCorrelator::GetFtofInfo: Skipping tracks not reaching the last FTS layer" << std::endl; return kFALSE; } if(helix->GetPz() <= 0.) { std::cout << "-W- PndPidCorrelator::GetFtofInfo: Skipping tracks going backward" << std::endl; return kFALSE; } PndFtofHit *tofHit = NULL; Int_t tofEntries = fFtofHit->GetEntriesFast(); Int_t tofIndex = -1; Float_t tofTof = 0., tofLength = -1000, tofGLength = -1000, tofTrackLength = -1000;; Float_t tofQuality = 1000000; //Float_t chi2 = 0; //[R.K. 01/2017] unused variable TVector3 vertex(0., 0., -10000.); TVector3 vertexrec(0., 0., -10000.); TVector3 momrec(0., 0., -10000.); TVector3 tofPos(0., 0., 0.); TVector3 momentum(0., 0., 0.); if (fGeanePro) // Overwrites vertex if Geane is used { // calculates track length from (0,0,0) to last point fGeanePropagator->SetPoint(TVector3(0,0,0)); fGeanePropagator->PropagateToPCA(1, -1); FairTrackParH *fRes= new FairTrackParH(); Bool_t rc = fGeanePropagator->Propagate(helix, fRes, fPidHyp*pidCand->GetCharge()); if (rc) { tofTrackLength = fGeanePropagator->GetLengthAtPCA(); vertexrec.SetXYZ(fRes->GetX(), fRes->GetY(), fRes->GetZ()); momrec.SetXYZ(fRes->GetPx(), fRes->GetPy(), fRes->GetPz()); } } for (Int_t tt = 0; ttAt(tt); if ( fIdeal && ( ((PndFtofPoint*)fFtofPoint->At(tofHit->GetRefIndex()))->GetTrackID() !=pidCand->GetMcIndex()) ) continue; tofHit->Position(tofPos); Float_t propX = helix->GetX() + (fCorrPar->GetFtofZ() - helix->GetZ()) * helix->GetPx() / helix->GetPz(); Float_t propY = helix->GetY() + (fCorrPar->GetFtofZ() - helix->GetZ()) * helix->GetPy() / helix->GetPz(); Float_t propZ = fCorrPar->GetFtofZ(); vertex.SetXYZ(propX, propY, propZ); tofGLength = (vertex-helix->GetPosition()).Mag(); if (fGeanePro) // Overwrites vertex if Geane is used { fGeanePropagator->SetPoint(tofPos); fGeanePropagator->PropagateToPCA(1, 1); FairTrackParH *fRes= new FairTrackParH(); Bool_t rc = fGeanePropagator->Propagate(helix, fRes, fPidHyp*pidCand->GetCharge()); if (!rc) continue; tofGLength = fGeanePropagator->GetLengthAtPCA(); vertex.SetXYZ(fRes->GetX(), fRes->GetY(), fRes->GetZ()); } Float_t dist = (tofPos-vertex).Mag2(); if ( tofQuality > dist) { tofIndex = tt; tofQuality = dist; tofTof = tofHit->GetTime(); tofLength = tofTrackLength+tofGLength; } if (fDebugMode) { Float_t ntuple[] = {static_cast(vertex.X()), static_cast(vertex.Y()), static_cast(vertex.Z()), static_cast(vertexrec.X()), static_cast(vertexrec.Y()), static_cast(vertexrec.Z()), static_cast(momrec.X()), static_cast(momrec.Y()), static_cast(momrec.Z()), static_cast(helix->GetMomentum().Mag()), static_cast(helix->GetQ()), static_cast(helix->GetMomentum().Theta()), static_cast(helix->GetZ()), static_cast(tofPos.X()), static_cast(tofPos.Y()), static_cast(tofPos.Z()), dist, tofLength, tofGLength, tofTrackLength}; // Float_t ntuple[] = {vertex.X(), vertex.Y(), vertex.Z(), // vertexrec.X(), vertexrec.Y(), vertexrec.Z(), // momrec.X(), momrec.Y(), momrec.Z(), // helix->GetMomentum().Mag(), helix->GetQ(), helix->GetMomentum().Theta(), helix->GetZ(), // tofPos.X(), tofPos.Y(), tofPos.Z(), // dist, tofLength, tofGLength,tofTrackLength}; ftofCorr->Fill(ntuple); } } if ( (tofQualityGetFTofCut()) || (fIdeal && tofIndex!=-1) ) { pidCand->SetTofQuality(tofQuality); pidCand->SetTofStopTime(tofTof); pidCand->SetTofTrackLength(tofLength); pidCand->SetTofIndex(tofIndex); if (tofLength>0.) { // mass^2 = p^2 * ( 1/beta^2 - 1 ) Float_t mass2 = helix->GetMomentum().Mag()*helix->GetMomentum().Mag()*(30.*30.*tofTof*tofTof/tofLength/tofLength-1.); pidCand->SetTofM2(mass2); } } return kTRUE; }