// ------------------------------------------------------------------------- // ----- PndDchPrepareKalmanTracks2 source file ----- // ----- Created 15.05.2008 by A. Wronska ----- // ----- based on the recotasks/demo code by S.Neubert ----- // ------------------------------------------------------------------------- // Panda Headers ---------------------- #include "CbmRootManager.h" #include "Track.h" #include "TrackCand.h" #include "CbmMCTrack.h" #include "LSLTrackRep.h" #include "GeaneTrackRep.h" #include "Kalman.h" #include "FitterExceptions.h" #include "CbmGeanePro.h" #include "CbmTrackParP.h" #include "PndDchPrepareKalmanTracks2.h" #include "PndDchPoint.h" #include "PndDchTrack.h" #include "PndDchTrackMatch.h" #include "PndDchCylinderHit.h" // C/C++ Headers ---------------------- #include #include // ROOT Class Headers -------- #include "TClonesArray.h" PndDchPrepareKalmanTracks2::PndDchPrepareKalmanTracks2() : CbmTask("Translation of PndDchTracks to Tracks"), fPersistence(kFALSE), fUseGeane(kFALSE) { } PndDchPrepareKalmanTracks2::~PndDchPrepareKalmanTracks2() { } InitStatus PndDchPrepareKalmanTracks2::Init() { //Get ROOT Manager CbmRootManager* ioman= CbmRootManager::Instance(); if(ioman==0) { Error("PndDchPrepareKalmanTracks2::Init","RootManager not instantiated!"); return kERROR; } AddHitBranch(1,"PndDchCylinderHit"); // open hit arrays std::map::iterator iter=fHitBranchNameMap.begin(); while(iter!=fHitBranchNameMap.end()){ TClonesArray* ar=(TClonesArray*) ioman->GetObject(iter->second); if(ar==0){ Error("DemoPRTask::Init","hit-array %s not found!",iter->second.Data()); } else{ fHitBranchMap[iter->first] = ar; } ++iter; }//end loops over hit types fCHitArray = fHitBranchMap[1]; // open MCTruth array fMcArray=(TClonesArray*) ioman->GetObject("MCTrack"); if(fMcArray==0){ Error("PndDchPrepareKalmanTracks2::Init","mctrack-array not found!"); return kERROR; } // open point array fDchPointArray=(TClonesArray*) ioman->GetObject("PndDchPoint"); if(fDchPointArray==0){ Error("PndDchPrepareKalmanTracks2::Init","dchpoint-array not found!"); return kERROR; } // open input array of PndDchTracks fDchTrackArray = (TClonesArray*) ioman->GetObject("PndDchTrack"); if(fDchTrackArray==0){ Error("PndDchPrepareKalmanTracks2::Init","PndDchTrack array not found!"); return kERROR; } std::cout<<"DchTrack array found "<< fDchTrackArray<GetObject("PndDchTrackMatch"); if(fDchTrackMatchArray==0){ Error("PndDchPrepareKalmanTracks2::Init","PndDchTrackMatch array not found!"); return kERROR; } std::cout<<"DchTrackMatch array found "<< fDchTrackMatchArray<Register("Track","GenFit",fTrackArray,fPersistence); std::cout<<"Track array created "<< fTrackArray<Delete(); std::map candmap; Int_t nuOfTracks = fDchTrackArray->GetEntriesFast(); for (Int_t id=0; idAt(id); Int_t tridx = dchtrmatch->GetRecTrackID(); PndDchTrack* dchtrack = (PndDchTrack*) fDchTrackArray->At(tridx); Int_t nuOfChits = dchtrack->GetNofDchCylinderHits(); std::cout<<"PndDchPrepareKalmanTracks2::Exec(): I found here "<GetDchCylinderHitIndex(nuhit); candmap[id]->addHit(1,globalCHitNu); } } // ------------------------------------------------------------------------ // try to find some starting values // loop over tracks std::map::iterator candIter=candmap.begin(); while(candIter!=candmap.end()){ TrackCand* cand=candIter->second; if(cand->getNHits()<10){ ++candIter; continue; } // Get MCTrack PndDchTrackMatch* dchTrackMatch = (PndDchTrackMatch*) fDchTrackMatchArray->At(candIter->first); Int_t mcTrId = dchTrackMatch->GetMCTrackID(); CbmMCTrack* mc=(CbmMCTrack*)fMcArray->At(mcTrId); if(mc==0){ Error("PndDchPrepareKalmanTracks2::Exec","MCTrack Id=&i not found!",mcTrId); ++candIter; continue; } int pdg=mc->GetPdgCode(); double q=TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/3.; Int_t pointidx = 0; TVector3 pos; TVector3 mom; while(pointidxGetEntries()){ PndDchPoint* pnt=(PndDchPoint*)fDchPointArray->At(pointidx); if(pnt->GetTrackID()==mcTrId){ pnt->Position(pos); pnt->Momentum(mom); break; } } // pos=mc->GetStartVertex(); Double_t startPosAccuracy = 0.5; TVector3 poserr(startPosAccuracy,startPosAccuracy,startPosAccuracy); TVector3 startMomAccuracy(0.1,0.1,0.1); TVector3 mcMom = mc->GetMomentum(); // mom.SetXYZ(gRandom->Gaus(mom.X(),mom.X()*startMomAccuracy.X()), // gRandom->Gaus(mom.Y(),mom.Y()*startMomAccuracy.Y()), // gRandom->Gaus(mom.Z(),mom.Z()*startMomAccuracy.Z())); TVector3 momerr(mom.X()*startMomAccuracy.X(), mom.Y()*startMomAccuracy.Y(), mom.Z()*startMomAccuracy.Z()); TVector3 u(1.,0.,0.); TVector3 v(0.,1.,0.); // create track-representation object and initialize with start values AbsTrackRep* rep=0; if(fUseGeane){ DetPlane pl(pos,u,v); GeaneTrackRep *grep=new GeaneTrackRep(fGeanePro,pl,mom,poserr,momerr,q,pdg); grep->setPropDir(1); // propagate in flight direction! std::cout<<" ^^^^^^^^^^^^^I prepare the following GeaneTrackRep:"<Print(); std::cout<<" ^^^^^^^^^^^^^End of GeaneTrackRep printout"<Print(); std::cout<<" ^^^^^^^^^^^^^End of LSLTrackRep printout"<Print(); // create track object Track* trk=new((*fTrackArray)[fTrackArray->GetEntriesFast()]) Track(rep); trk->setCandidate(*cand); // here the candidate is copied! ++candIter; }// end loop over tracks std::cout<GetEntriesFast()<<" tracks created"<