#include "PndTorinoDetectorTrackFinderIdeal.h" #include "PndTorinoDetectorPoint.h" #include "PndTorinoDetectorHit.h" #include "PndTrackCand.h" #include "PndTrack.h" #include "PndMCTrack.h" #include "FairHit.h" #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "TClonesArray.h" #include "TDatabasePDG.h" #include "TRandom.h" #include using namespace std; PndTorinoDetectorTrackFinderIdeal::PndTorinoDetectorTrackFinderIdeal() : FairTask(), fTrackCandArray(new TClonesArray("PndTrackCand")), fTrackArray(new TClonesArray("PndTrack")), fMCTrackArray(NULL), fPointArray(NULL), fHitArray(NULL), fSmearMomentum(kTRUE), fErrpx(0.02), fErrpy(0.02), fErrpz(0.03) {} PndTorinoDetectorTrackFinderIdeal::~PndTorinoDetectorTrackFinderIdeal() { fMCTrackArray->Delete(); fPointArray->Delete(); fHitArray->Delete(); } // initialization InitStatus PndTorinoDetectorTrackFinderIdeal::Init() { // instance of FairRootManager FairRootManager *ioman = FairRootManager::Instance(); if (!ioman) { cout << "-E- RootManager not instantised!" << endl; return kFATAL; } // get input data ------------------------------------------ // get MC point array fPointArray = (TClonesArray*) ioman->GetObject("PndTorinoDetectorPoint"); if (!fPointArray) { cout << "-E- PndTorinoDetectorTrackFinderIdeal::Init: No MC point array!" << endl; return kERROR; } // get MC track array fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack"); if (!fMCTrackArray) { cout << "-E- PndTorinoDetectorTrackFinderIdeal::Init: No MC track array!" << endl; return kERROR; } // get hit array fHitArray = (TClonesArray*) ioman->GetObject("PndTorinoDetectorHit"); if (!fHitArray) { cout << "-E- PndTorinoDetectorTrackFinderIdeal::Init: No hit array!" << endl; return kERROR; } // register output ----------------------------------------- ioman->Register("PndTorinoDetectorTrackCand", "PndTorinoDetector", fTrackCandArray, kTRUE); ioman->Register("PndTorinoDetectorTrack", "PndTorinoDetector", fTrackArray, kTRUE); return kSUCCESS; } // Exec void PndTorinoDetectorTrackFinderIdeal::Exec(Option_t* opt) { // call this to clean up the output arrays fTrackCandArray->Delete(); fTrackArray->Delete(); // create the trackcand and track object PndTrackCand* trackcand = NULL; PndTrack* track = NULL; PndTorinoDetectorHit *hit = NULL; PndTorinoDetectorPoint* point = NULL; Int_t nMCtracks = fMCTrackArray->GetEntriesFast(); Int_t nhits = fHitArray->GetEntriesFast(); // create map & list Int_t mapTrack2Hit[nMCtracks][nhits]; Int_t listofhits[nMCtracks]; // initialization for(int itrk = 0; itrk < nMCtracks; itrk++) listofhits[itrk] = 0; for(int itrk = 0; itrk < nMCtracks; itrk++) for(int ihit = 0; ihit < nhits; ihit++) mapTrack2Hit[itrk][ihit] = 0; // loop over the hits ...................................................... // 1. get the associated point, via refindex // 2. get its MCtrack via trackID for(int ihit = 0; ihit < fHitArray->GetEntries(); ihit++) { hit = (PndTorinoDetectorHit*) fHitArray->At(ihit); if(!hit) continue; // 1. get the associated point, via refindex int refindex = hit->GetRefIndex(); point = (PndTorinoDetectorPoint*) fPointArray->At(refindex); if(!point) continue; // 2. get its MCtrack via trackID Int_t trackID = point->GetTrackID(); // match this hit to the MC track #trackID mapTrack2Hit[trackID][listofhits[trackID]] = ihit; // how many hits belong to the MC track #trackID? listofhits[trackID]++; } // ========================================================================= // ========================== set up PndTrack(Cand) ======================== // loop over the MC tracks for(int itrk = 0; itrk < fMCTrackArray->GetEntries(); itrk++) { // if the itrk MC track left no hits, continue; if(listofhits[itrk] == 0) continue; // else, create a new TRACK CAND ........................................ Int_t size = fTrackCandArray->GetEntriesFast(); trackcand = new ((*fTrackCandArray)[size]) PndTrackCand(); // add all the hits for(int ihit = 0; ihit < listofhits[itrk]; ihit++) { Int_t hitID = mapTrack2Hit[itrk][ihit]; // .............. <------------------------ detID -------------------------------> . hitID . rho . trackcand->AddHit(FairRootManager::Instance()->GetBranchId("PndTorinoDetectorHit"), hitID, hitID); // CHECK } // now the TRACK parameters ---------------------------------------------- // we need to have TVector3 firstPos, firstMom, firstPosErr, firstMomErr; TVector3 lastPos, lastMom, lastPosErr, lastMomErr; TVector3 dj, dk; // --------------------------------------------------- // take first and last hit values PndMCTrack *mctrack = (PndMCTrack*) fMCTrackArray->At(itrk); int pdgcode = mctrack->GetPdgCode(); // charge // get and instance of the PDG database TDatabasePDG *dbPDG = TDatabasePDG::Instance(); TParticlePDG *particle = dbPDG->GetParticle(pdgcode); if(!particle) continue; // charge in units of |e|/3 double q = particle->Charge()/3.; Int_t charge = (Int_t) q; Int_t hitcounter = trackcand->GetNHits(); // 1st hit ............................................................... // get cand hit PndTrackCandHit candhit = trackcand->GetSortedHit(0); Int_t hitID = candhit.GetHitId(); // get hit PndTorinoDetectorHit *firsthit = (PndTorinoDetectorHit*) fHitArray->At(hitID); if(!firsthit) continue; // get MC point PndTorinoDetectorPoint *firstpnt = (PndTorinoDetectorPoint*) fPointArray->At(firsthit->GetRefIndex()); if(!firstpnt) continue; // now, we can get from the first point, the position and momentum, from MC truth firstPos.SetXYZ(firsthit->GetX(), firsthit->GetY(), firsthit->GetZ()); firstMom.SetXYZ(firstpnt->GetPx(), firstpnt->GetPy(), firstpnt->GetPz()); // last hit ............................................................... // get cand hit candhit = trackcand->GetSortedHit(hitcounter - 1); hitID = candhit.GetHitId(); // get hit PndTorinoDetectorHit *lasthit = (PndTorinoDetectorHit*) fHitArray->At(hitID); // get MC point PndTorinoDetectorPoint *lastpnt = (PndTorinoDetectorPoint*) fPointArray->At(lasthit->GetRefIndex()); if(!lastpnt) continue; // now, we can get from the last point, the position and momentum, from MC truth lastPos.SetXYZ(lasthit->GetX(), lasthit->GetY(), lasthit->GetZ()); lastMom.SetXYZ(lastpnt->GetPx(), lastpnt->GetPy(), lastpnt->GetPz()); // out plane is the xy plane! ............................................... dj.SetXYZ(1., 0., 0.), dk.SetXYZ(0., 1., 0.); // smearing of the track momentum at start/last points with gaussians ....... if(fSmearMomentum) { double smear1_momx = fErrpx * firstMom.X(); double smear1_momy = fErrpy * firstMom.Y(); double smear1_momz = fErrpz * firstMom.Z(); double smeared1_momx = gRandom->Gaus(firstMom.X(), smear1_momx); double smeared1_momy = gRandom->Gaus(firstMom.Y(), smear1_momy); double smeared1_momz = gRandom->Gaus(firstMom.Z(), smear1_momz); firstMom.SetXYZ(smeared1_momx, smeared1_momy, smeared1_momz); double smear2_momx = fErrpx * lastMom.X(); double smear2_momy = fErrpy * lastMom.Y(); double smear2_momz = fErrpz * lastMom.Z(); double smeared2_momx = gRandom->Gaus(lastMom.X(), smear2_momx); double smeared2_momy = gRandom->Gaus(lastMom.Y(), smear2_momy); double smeared2_momz = gRandom->Gaus(lastMom.Z(), smear2_momz); lastMom.SetXYZ(smeared2_momx, smeared2_momy, smeared2_momz); // retrieve the errors ..................................................... firstPosErr.SetXYZ(firsthit->GetDx(), firsthit->GetDy(), firsthit->GetDz()); firstMomErr.SetXYZ(smear1_momx, smear1_momy, smear1_momz); lastPosErr.SetXYZ(lasthit->GetDx(), lasthit->GetDy(), lasthit->GetDz()); lastMomErr.SetXYZ(smear2_momx, smear2_momy, smear2_momz); } // create parameters @ 1st ................................................. FairTrackParP first(firstPos, firstMom, firstPosErr, firstMomErr, charge, firstPos, dj, dk); // create parameters @ last ................................................. FairTrackParP last(lastPos, lastMom, lastPosErr, lastMomErr, charge, lastPos, dj, dk); size = fTrackArray->GetEntriesFast(); // fill the track ........................................................... track = new((*fTrackArray)[size]) PndTrack(first, last, *trackcand); // just some print outs ..................................................... if(fVerbose > 0) { cout << "creating track for track cand " << size << " with hits: "; hitcounter = trackcand->GetNHits(); for(int ihit = 0; ihit < hitcounter; ihit++) { candhit = trackcand->GetSortedHit(ihit); cout << " " << candhit.GetHitId(); } cout << endl; } } } // ClassImp(ClassName) ClassImp(PndTorinoDetectorTrackFinderIdeal)