// ------------------------------------------------------------------------------- // ----- PndDchPreFitterTR source file ----- // ----- Created 26.02.2008 by A. Wronska, P.Hawranek ----- // ------------------------------------------------------------------------------- #include #include #include "TClonesArray.h" #include "TVector3.h" #include "TVector2.h" #include "TString.h" #include "TGeoMatrix.h" #include "TH1F.h" #include "TMath.h" #include "TFile.h" #include "FairRootManager.h" #include "FairRuntimeDb.h" #include "FairBaseParSet.h" #include "FairRunAna.h" #include "FairRun.h" #include "FairTrackParam.h" #include "PndDchPreFitterTR.h" #include "PndDchHit.h" #include "PndDchPoint.h" #include "PndDchCylinderHit.h" #include "PndTrackCand.h" #include "PndTrackCandHit.h" #include "PndDchStructure.h" #include "PndDchMapper.h" using std::cout; using std::endl; PndDchPreFitterTR::PndDchPreFitterTR() : FairTask("Dch Prefitter for Forward Spectrometer") { fInHitsAfterXZ = 0; fInHitsInXZ = 0; fInHitsBeforeXZ = 0; fInHitsYZ = 0; fDetIdList = 0; } PndDchPreFitterTR::~PndDchPreFitterTR() { if (0 != fInHitsBeforeXZ) { fInHitsBeforeXZ->Delete(); delete fInHitsBeforeXZ; } if (0 != fInHitsAfterXZ) { fInHitsAfterXZ->Delete(); delete fInHitsAfterXZ; } if (0 != fInHitsInXZ) { fInHitsInXZ->Delete(); delete fInHitsInXZ; } if (0 != fInHitsYZ) { fInHitsYZ->Delete(); delete fInHitsYZ; } if( 0 != fDetIdList){ delete fDetIdList; } } // Public method Init InitStatus PndDchPreFitterTR::Init() { cout << "PndDchPreFitterTR::Init()... " << endl; // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); TFile *infile = ioman->GetInFile(); if ( ! ioman ) { cout << "-E- PndDchPreFitterTR::Init():\n\t " << "RootManager not instantiated!" << endl; return kFATAL; } // TGeoManager *geoMan = (TGeoManager*) infile->Get("FAIRGeom"); TGeoManager *geoMan =gGeoManager; fStructure = PndDchStructure::Instance(geoMan); // Get input point array - only to get real value of momentum! fPointArray = (TClonesArray*) ioman->GetObject("PndDchPoint"); if ( ! fPointArray ) { cout << "-W- PndDchPreFitterTR::Init(): " << "No input PndDchPoint array!" << endl; return kERROR; } // Get input hit array fInHitArray = (TClonesArray*) ioman->GetObject("PndDchCylinderHit"); if ( ! fInHitArray ) { cout << "-W- PndDchPreFitterTR::Init(): " << "No input PndDchCylinderHit array!" << endl; return kERROR; } fDetIdList = fStructure->GetDetectorIDList(); fDetIdListSize = fDetIdList->GetSize(); // Get input track array fTrackArray = (TClonesArray*) ioman->GetObject("DCHTrackCand"); if ( !fTrackArray ) { cout << "-W- PndDchPreFitterTR::Init(): No DCHTrackCand array!" << endl; return kERROR; } FairRunAna *fRun= FairRunAna::Instance(); Double_t fBeamMom = 0.; if(fRun){ FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); FairBaseParSet* par=(FairBaseParSet*) (rtdb->findContainer("FairBaseParSet")); fBeamMom= par->GetBeamMom(); std::cout<<" Beam momentum read out from DB = "<GetEntries(); PndTrackCand* track = 0; for(Int_t iTrack = 0; iTrack < nTracks; iTrack++) { track = (PndTrackCand*) fTrackArray->At(iTrack); fInHitsBeforeXZ->Delete(); fInHitsAfterXZ->Delete(); fInHitsInXZ->Delete(); fInHitsYZ->Delete(); TVector2* tmp; std::map mapa; mapa = GetHitPointsInChambers(track); std::map::iterator iter; if(fVerbose>0){ cout<<"Cross Points in Chambers"<first<second.Print(); } cout<<"Cross Points in Chambers END"<second.X(),mapa.find(3)->second.Z()); tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(3)->second.Y(),mapa.find(3)->second.Z()); } if (mapa.find(4) != mapa.end()) { tmp = new((*fInHitsBeforeXZ)[nBeforeXZ++]) TVector2(mapa.find(4)->second.X(),mapa.find(4)->second.Z()); tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(4)->second.Y(),mapa.find(4)->second.Z()); } if (mapa.find(5) != mapa.end()) { tmp = new((*fInHitsInXZ)[nInXZ++]) TVector2(mapa.find(5)->second.X(),mapa.find(5)->second.Z()); tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(5)->second.Y(),mapa.find(5)->second.Z()); } if (mapa.find(6) != mapa.end()) { tmp = new((*fInHitsInXZ)[nInXZ++]) TVector2(mapa.find(6)->second.X(),mapa.find(6)->second.Z()); tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(6)->second.Y(),mapa.find(6)->second.Z()); } if (mapa.find(7) != mapa.end()) { tmp = new((*fInHitsAfterXZ)[nAfterXZ++]) TVector2(mapa.find(7)->second.X(),mapa.find(7)->second.Z()); // tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(7)->second.Y(),mapa.find(7)->second.Z()); } if (mapa.find(8) != mapa.end()) { tmp = new((*fInHitsAfterXZ)[nAfterXZ++]) TVector2(mapa.find(8)->second.X(),mapa.find(8)->second.Z()); // tmp = new((*fInHitsYZ)[nYZ++]) TVector2(mapa.find(8)->second.Y(),mapa.find(8)->second.Z()); } Bool_t isRecoTrackBeforeXZ = GetLineParameters(fInHitsBeforeXZ,fTrackBeforeXZ); if(isRecoTrackBeforeXZ){ Double_t xIn = (fZFieldBegin-fTrackBeforeXZ.Y())/fTrackBeforeXZ.X(); //artificial dch inhit for entering magnetic field TVector2 *pos = new((*fInHitsInXZ)[nInXZ++]) TVector2(xIn,fZFieldBegin); } Bool_t isRecoTrackAfterXZ = GetLineParameters(fInHitsAfterXZ,fTrackAfterXZ); if(isRecoTrackAfterXZ){ Double_t xOut = (fZFieldEnd-fTrackAfterXZ.Y())/fTrackAfterXZ.X(); //artificial dch inhit for exiting the field TVector2 *pos = new((*fInHitsInXZ)[nInXZ++]) TVector2(xOut,fZFieldEnd); } Bool_t isRecoTrackYZ = GetLineParameters(fInHitsYZ,fTrackYZ); Bool_t isRecoTrackInXZ = GetCircleParameters(fInHitsInXZ,fTrackInXZ); TVector3 startMomentum(0., 0., 0.); TVector3 startPosition(0., 0., 0.); if(!(isRecoTrackYZ && isRecoTrackBeforeXZ && isRecoTrackInXZ) ){ cout<<"Tracking failed (YZ BeforeXZ InXZ AfterXZ) " <GetSortedHit(0); Int_t idx = candHit.GetHitId(); PndDchCylinderHit* chit = (PndDchCylinderHit*)fInHitArray->At(idx); GetInHitAtZ(chit->GetWireZcoordGlobal(),startPosition); Int_t chargeSign = GetChargeSign(); const TMatrixFSym* covMatrix = new TMatrixFSym(15); FairTrackParam parset(startPosition.X(),startPosition.Y(),startPosition.Z(), startMomentum.X()/ startMomentum.Z(), startMomentum.Y()/ startMomentum.Z(), chargeSign/ startMomentum.Mag(), *covMatrix); // track->SetParamFirst(parset); track->setTrackSeed(startPosition, TVector3(startMomentum.X()/ startMomentum.Z(), startMomentum.Y()/ startMomentum.Z(), 1.), chargeSign/startMomentum.Mag()); } //end of loop over tracks cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"< fZFieldBegin && zpos < fZFieldEnd) { Double_t sgn=1; if(fTrackInXZ.Z()>zpos) sgn=-1; Double_t tmp=pow(fTrackInXZ.Z(),2)-pow(zpos-fTrackInXZ.Y(),2); if (tmp<0) return kFALSE; xpos=sgn*sqrt(tmp)+fTrackInXZ.X(); //after dipole } else if( fTrackAfterXZ.X()!=0){ xpos =(zpos-fTrackAfterXZ.Y())/fTrackAfterXZ.X(); } else return kFALSE; inhit.SetX(xpos); return kTRUE; } Int_t PndDchPreFitterTR::GetChargeSign(){ Double_t zReal = ((TVector2*)fInHitsInXZ->At( fInHitsInXZ->GetEntries()-1))->Y(); Double_t xReal = ((TVector2*)fInHitsInXZ->At( fInHitsInXZ->GetEntries()-1))->X(); Double_t xExtrapolated = (zReal-fTrackBeforeXZ.Y())/fTrackBeforeXZ.X(); if(xReal>xExtrapolated) return -1; return 1; } Bool_t PndDchPreFitterTR::GetMomentum(TVector3& momentum ) { Double_t magXZ = 0.3*fField*fTrackInXZ(2)/100.; // std::cout<<"magXZ= "<GetEntries(); if(nInHits<2) { cout<<" nHits="<At(0); p2=(TVector2*)fInHits->At(1); if(TMath::Abs(p2->X()-p1->X())<1e-10) { Error("PndDchPreFitterTR::GetLineParameters","2 points paralel with x1=%f and x2=%f",p2->X(),p1->X()); return kFALSE; } Double_t slope = (p2->Y()-p1->Y())/(p2->X()-p1->X()); Double_t offset = p2->Y()-slope*p2->X(); result.Set(slope, offset); return kTRUE; } for (Int_t iInHit=0; iInHitAt(iInHit); if ((inhit->X()!=0) && (inhit->Y()!=0)) { x += inhit->X(); y += inhit->Y(); xx += inhit->X()*inhit->X(); xy += inhit->X()*inhit->Y(); n++; } } WG = x*x - n*xx; WA = x*y - n*xy; WB = x*xy - xx*y; if (WG==0){ return kFALSE; } //result = (slope, offset) result.Set(WA/WG,WB/WG); return kTRUE; } Bool_t PndDchPreFitterTR::GetCrossPoint(PndDchCylinderHit *chit1, PndDchCylinderHit *chit2, TVector2 &point){ TVector2 wire1End1 = chit1->GetWireEnd1(); TVector2 wire1End2 = chit1->GetWireEnd2(); TVector2 wire2End1 = chit2->GetWireEnd1(); TVector2 wire2End2 = chit2->GetWireEnd2(); TVector2 dir1 = (wire1End2-wire1End1).Unit(); TVector2 dir2 = (wire2End2-wire2End1).Unit(); if((dir1-dir2).Mod()<1e-6 || (dir1+dir2).Mod()<1e-6) { return kFALSE; //lines are parallel } if(TMath::Abs((dir1+dir2).X())<1e-6){ // cout<<"The same X"< PndDchPreFitterTR::GetHitPointsInChambers(PndTrackCand* tr) { Int_t nHits = tr->GetNHits(); std::map chamberHitMap; // contains map chamberNo<->average of hit wires std::map chamberCountsMap; // contains map chamber No<->number of hits std::map::iterator iter; std::map::iterator iterCounts; Int_t ch1, ch2; for(Int_t i = 0; iGetSortedHit(i); Int_t idx1 = candHit1.GetHitId(); PndDchCylinderHit* chit1 = (PndDchCylinderHit*)fInHitArray->At(idx1); ch1 = fStructure->InWhichChamber(chit1->GetWireZcoordGlobal()); for(Int_t j = i; jGetSortedHit(j); Int_t idx2 = candHit2.GetHitId(); PndDchCylinderHit* chit2 = (PndDchCylinderHit*)fInHitArray->At(idx2); ch2 = fStructure->InWhichChamber(chit2->GetWireZcoordGlobal()); if(ch1 != ch2) continue; TVector2 wi1, wi2, point; if (GetCrossPoint(chit1,chit2,point)) { TVector3 point3d(point.X(), point.Y(), (chit1->GetWireZcoordGlobal()+chit2->GetWireZcoordGlobal())/2.); iter = chamberHitMap.find(ch1); if( iter == chamberHitMap.end() ) { chamberHitMap[ch1] = point3d; } else { point3d = point3d + iter->second; chamberHitMap[ch1] = point3d; } ++chamberCountsMap[ch1]; } } } for( iter = chamberHitMap.begin(); iter != chamberHitMap.end(); ++iter ) { iterCounts = chamberCountsMap.find(iter->first); iter->second = iter->second *(1./(Double_t)iterCounts->second); } return chamberHitMap; } ClassImp(PndDchPreFitterTR)