From 1d26da6d8e27b8543e7c65c22f19dff5775cbd93 Mon Sep 17 00:00:00 2001 From: shahoian Date: Fri, 17 Jan 2014 17:52:18 +0100 Subject: [PATCH 1/1] Modifications needed to use PID framework based mass during tracking and to include heavy nuclei into tracking. The ESDtrack PID probability arrays are made transient (kept for backward compatibility only) --- HMPID/AliHMPIDtrack.cxx | 2 +- ITS/AliITStrackV2.cxx | 4 +- ITS/AliITStrackerMI.cxx | 4 -- ITS/UPGRADE/AliITSUTrackHyp.cxx | 4 +- ITS/UPGRADE/AliITSUTrackerGlo.cxx | 5 +- STEER/ESD/AliESDpid.cxx | 26 +++++++++ STEER/ESD/AliESDpid.h | 4 ++ STEER/ESD/AliESDtrack.cxx | 64 +++++------------------ STEER/ESD/AliESDtrack.h | 17 +++--- STEER/ESD/AliKalmanTrack.cxx | 11 ++-- STEER/ESD/AliTrackerBase.cxx | 4 +- STEER/STEER/AliReconstruction.cxx | 29 ++++++---- STEER/STEER/AliReconstruction.h | 5 +- STEER/STEERBase/AliExternalTrackParam.cxx | 28 +++++++--- STEER/STEERBase/AliPID.cxx | 14 +++-- STEER/STEERBase/AliPID.h | 12 ++--- TOF/AliTOFtrack.cxx | 2 +- TPC/Rec/AliTPCtrack.cxx | 2 +- TRD/AliTRDtrackV1.cxx | 4 +- TRD/AliTRDtrackerV1.cxx | 2 +- 20 files changed, 133 insertions(+), 110 deletions(-) diff --git a/HMPID/AliHMPIDtrack.cxx b/HMPID/AliHMPIDtrack.cxx index 451a83c99ac..1f742105e39 100644 --- a/HMPID/AliHMPIDtrack.cxx +++ b/HMPID/AliHMPIDtrack.cxx @@ -41,7 +41,7 @@ AliHMPIDtrack::AliHMPIDtrack(const AliESDtrack& t):AliKalmanTrack() // SetLabel(t.GetLabel()); SetChi2(0.); - SetMass(t.GetMass()); + SetMass(t.GetMassForTracking()); Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance()); diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index ae71822d02f..86a2fe5166d 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -66,7 +66,7 @@ AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c): Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance()); SetLabel(t.GetITSLabel()); - SetMass(t.GetMass()); + SetMass(t.GetMassForTracking()); SetNumberOfClusters(t.GetITSclusters(fIndex)); if (t.GetStatus()&AliESDtrack::kTIME) { @@ -456,6 +456,7 @@ Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) { Double_t dx = x - xv, dy = par[0] - yv, dz = par[1] - zv; Double_t r2=dx*dx + dy*dy; Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt()); + if (GetMass()<0) p2 *= 4; // q=2 Double_t beta2=p2/(p2 + GetMass()*GetMass()); x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp())); Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0; @@ -690,6 +691,7 @@ Bool_t AliITStrackV2::ImproveKalman(Double_t xyz[3],Double_t ers[3], const Doubl double ms44t = p34*p34; // double p2=(1.+ par[3]*par[3])/(par[4]*par[4]); + if (GetMass()<0) p2 *= 4; // q=2 double beta2 = p2/(p2+GetMass()*GetMass()); double theta2t = 14.1*14.1/(beta2*p2*1e6) * (1. + par[3]*par[3]); // diff --git a/ITS/AliITStrackerMI.cxx b/ITS/AliITStrackerMI.cxx index e6877423890..15cc909281a 100644 --- a/ITS/AliITStrackerMI.cxx +++ b/ITS/AliITStrackerMI.cxx @@ -563,7 +563,6 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) { // temporary Int_t noesd = 0; {/* Read ESD tracks */ - Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass(); Int_t nentr=event->GetNumberOfTracks(); noesd=nentr; // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr); @@ -580,9 +579,6 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) { t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B. Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1)); - - // look at the ESD mass hypothesys ! - if (t->GetMass()<0.9*pimass) t->SetMass(pimass); // write expected q t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.)); diff --git a/ITS/UPGRADE/AliITSUTrackHyp.cxx b/ITS/UPGRADE/AliITSUTrackHyp.cxx index fcbf4ce9926..28b4bd18fc0 100644 --- a/ITS/UPGRADE/AliITSUTrackHyp.cxx +++ b/ITS/UPGRADE/AliITSUTrackHyp.cxx @@ -81,7 +81,7 @@ AliITSUTrackHyp::AliITSUTrackHyp(const AliESDtrack &src) { // copy c-tor. from ESD track: take only kinematics, mass and time integral AliExternalTrackParam::operator=((const AliExternalTrackParam&)src); - SetMass(src.GetMass()); + SetMass(src.GetMassForTracking()); if (src.IsOn(AliESDtrack::kTIME)) { StartTimeIntegral(); SetIntegratedLength(src.GetIntegratedLength()); @@ -130,7 +130,7 @@ AliITSUTrackHyp &AliITSUTrackHyp::operator=(const AliESDtrack &src) { // copy oparator from ESD track: take only kinematics, mass and time integral AliExternalTrackParam::operator=((const AliExternalTrackParam&)src); - SetMass(src.GetMass()); + SetMass(src.GetMassForTracking()); if (src.IsOn(AliESDtrack::kTIME)) { StartTimeIntegral(); SetIntegratedLength(src.GetIntegratedLength()); diff --git a/ITS/UPGRADE/AliITSUTrackerGlo.cxx b/ITS/UPGRADE/AliITSUTrackerGlo.cxx index c42c13cd1f1..fc843fcfe8b 100644 --- a/ITS/UPGRADE/AliITSUTrackerGlo.cxx +++ b/ITS/UPGRADE/AliITSUTrackerGlo.cxx @@ -263,7 +263,7 @@ Int_t AliITSUTrackerGlo::Clusters2Tracks(AliESDEvent *esdEv) AliLog::SetClassDebugLevel("AliITSUTrackerGlo",dbg ? 10:0); */ #ifdef _ITSU_DEBUG_ - AliDebug(1,Form("Processing track %d(esd%d) | M=%.3f Pt=%.3f | MCLabel: %d",itr,trID,fCurrESDtrack->GetMass(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS + AliDebug(1,Form("Processing track %d(esd%d) | M=%.3f Pt=%.3f | MCLabel: %d",itr,trID,fCurrESDtrack->GetMassForTracking(kTRUE),fCurrESDtrack->Pt(),fCurrESDtrMClb));//RS #endif FindTrack(fCurrESDtrack, trID); } @@ -680,8 +680,7 @@ AliITSUTrackHyp* AliITSUTrackerGlo::InitHypothesis(AliESDtrack *esdTr, Int_t esd // fCountProlongationTrials++; // - fCurrMass = esdTr->GetMass(); - if (fCurrMassGetMassForTracking(); // hyp = new AliITSUTrackHyp(fNLrActive); hyp->SetESDTrack(esdTr); diff --git a/STEER/ESD/AliESDpid.cxx b/STEER/ESD/AliESDpid.cxx index 89f6d92f8d4..66020a62214 100644 --- a/STEER/ESD/AliESDpid.cxx +++ b/STEER/ESD/AliESDpid.cxx @@ -467,3 +467,29 @@ Float_t AliESDpid::GetNumberOfSigmasTOFold(const AliVParticle *track, AliPID::EP Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type); return (tofTime - fTOFResponse.GetStartTime(vtrack->P()) - expTime)/fTOFResponse.GetExpectedSigma(vtrack->P(),expTime,AliPID::ParticleMassZ(type)); } + +//_________________________________________________________________________ +void AliESDpid::SetMassForTracking(AliESDtrack *esdtr) const +{ + // assign mass for tracking + // + int pid = AliPID::kPion; // this should be substituted by real most probable TPC pid (e,mu -> pion) or poin if no PID possible + // + if (pidAliPID::kSPECIESC-1) pid = AliPID::kPion; + // + esdtr->SetMassForTracking( AliPID::ParticleCharge(pid)==1 ? AliPID::ParticleMass(pid) : -AliPID::ParticleMass(pid)); + // +} + + +//_________________________________________________________________________ +void AliESDpid::MakePIDForTracking(AliESDEvent *event) const +{ + // assign masses using for tracking + Int_t nTrk=event->GetNumberOfTracks(); + for (Int_t iTrk=nTrk; iTrk--;) { + AliESDtrack *track = event->GetTrack(iTrk); + SetMassForTracking(track); + } + +} diff --git a/STEER/ESD/AliESDpid.h b/STEER/ESD/AliESDpid.h index 2c16383ea7e..6e1d5279e76 100644 --- a/STEER/ESD/AliESDpid.h +++ b/STEER/ESD/AliESDpid.h @@ -32,6 +32,8 @@ AliESDpid& operator=(const AliESDpid& a){if (this==&a) return *this; AliPIDRespo virtual ~AliESDpid() {} Int_t MakePID(AliESDEvent *event, Bool_t TPCOnly = kFALSE, Float_t timeZeroTOF=9999) const; + void MakePIDForTracking(AliESDEvent *event) const; + void MakeTPCPID(AliESDtrack *track) const; void MakeITSPID(AliESDtrack *track) const; void MakeTOFPID(AliESDtrack *track, Float_t /*timeZeroTOF*/) const; @@ -40,6 +42,8 @@ AliESDpid& operator=(const AliESDpid& a){if (this==&a) return *this; AliPIDRespo void MakeTRDPID(AliESDtrack *track) const; void CombinePID(AliESDtrack *track) const; + void SetMassForTracking(AliESDtrack *track) const; + // Float_t NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type) const {return AliPIDResponse::NumberOfSigmasTOF(track,type);} // Float_t GetNumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type, const Float_t timeZeroTOF) const; diff --git a/STEER/ESD/AliESDtrack.cxx b/STEER/ESD/AliESDtrack.cxx index 0e228f7bb13..da01af84eb4 100644 --- a/STEER/ESD/AliESDtrack.cxx +++ b/STEER/ESD/AliESDtrack.cxx @@ -185,6 +185,7 @@ AliESDtrack::AliESDtrack() : fHMPIDqn(0), fHMPIDcluIdx(-1), fCaloIndex(kEMCALNoMatch), + fMassForTracking(0.13957), fHMPIDtrkTheta(0), fHMPIDtrkPhi(0), fHMPIDsignal(0), @@ -259,14 +260,7 @@ AliESDtrack::AliESDtrack() : Int_t i; for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0; - for (i=0; i 100) return; for (Int_t i=0; iGetProcInfo(&procInfo); AliInfo(Form("Current memory usage %ld %ld", procInfo.fMemResident, procInfo.fMemVirtual)); + // PID + fESDpid = new AliESDpid(); + //QA //Initialize the QA and start of cycle if (fRunQA || fRunGlobalQA) @@ -1915,8 +1923,6 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) AliCodeTimerAuto("",0); - AliESDpid pid; - AliSysInfo::AddStamp(Form("StartEv_%d",iEvent), 0,0,iEvent); if (iEvent >= fRunLoader->GetNumberOfEvents()) { @@ -1961,7 +1967,7 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) if (reconstructor && fRecoParam.GetDetRecoParamArray(iDet)) { const AliDetectorRecoParam *par = fRecoParam.GetDetRecoParam(iDet); reconstructor->SetRecoParam(par); - reconstructor->GetPidSettings(&pid); + reconstructor->GetPidSettings(fESDpid); reconstructor->SetEventInfo(&fEventInfo); if (fRunQA) { AliQAManager::QAManager()->SetEventInfo(&fEventInfo) ; @@ -2122,7 +2128,7 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) // barrel tracking if (!fRunTracking.IsNull()) { - if (!RunTracking(fesd,pid)) { + if (!RunTracking(fesd,*fESDpid)) { if (fStopOnError) {CleanUp(); return kFALSE;} } } @@ -2166,7 +2172,7 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) ok = kFALSE; if (tpcTrack) ok = AliTracker:: - PropagateTrackToBxByBz(tpcTrack,kRadius,track->GetMass(),kMaxStep,kFALSE); + PropagateTrackToBxByBz(tpcTrack,kRadius,track->GetMassForTracking(),kMaxStep,kFALSE); if (ok) { Int_t n=trkArray.GetEntriesFast(); @@ -2178,7 +2184,7 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) if (track->IsOn(AliESDtrack::kITSrefit)) continue; AliTracker:: - PropagateTrackToBxByBz(track,kRadius,track->GetMass(),kMaxStep,kFALSE); + PropagateTrackToBxByBz(track,kRadius,track->GetMassForTracking(),kMaxStep,kFALSE); Double_t x[3]; track->GetXYZ(x); Double_t b[3]; AliTracker::GetBxByBz(x,b); track->RelateToVertexBxByBz(fesd->GetPrimaryVertexSPD(), b, kVeryBig); @@ -2303,10 +2309,10 @@ Bool_t AliReconstruction::ProcessEvent(Int_t iEvent) // AdC+FN if (fReconstructor[3]) - GetReconstructor(3)->FillEventTimeWithTOF(fesd,&pid); + GetReconstructor(3)->FillEventTimeWithTOF(fesd,fESDpid); // combined PID - pid.MakePID(fesd); + // fESDpid->MakePID(fesd); if (fFillTriggerESD) { if (!FillTriggerESD(fesd)) { @@ -2972,7 +2978,7 @@ Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd,AliESDpid &PID) // preliminary PID in TPC needed by the ITS tracker if (iDet == 1) { GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd); - PID.MakePID(esd,kTRUE); + PID.MakePIDForTracking(esd); AliSysInfo::AddStamp(Form("MakePID0%s_%d",fgkDetectorName[iDet],eventNr), iDet,4,eventNr); } } @@ -3027,7 +3033,7 @@ Bool_t AliReconstruction::RunTracking(AliESDEvent*& esd,AliESDpid &PID) if (iDet == 1) { //GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd); //AliESDpid::MakePID(esd); - PID.MakePID(esd,kTRUE); + PID.MakePIDForTracking(esd); AliSysInfo::AddStamp(Form("MakePID1%s_%d",fgkDetectorName[iDet],eventNr), iDet,4,eventNr); } @@ -3559,6 +3565,9 @@ void AliReconstruction::CleanUp() delete fParentRawReader; fParentRawReader=NULL; + delete fESDpid; + fESDpid = NULL; + if (ffile) { ffile->Close(); delete ffile; diff --git a/STEER/STEER/AliReconstruction.h b/STEER/STEER/AliReconstruction.h index 6946a67fe3d..f10ef53b9bd 100644 --- a/STEER/STEER/AliReconstruction.h +++ b/STEER/STEER/AliReconstruction.h @@ -363,6 +363,9 @@ private: // Plane Efficiency Evaluation Bool_t fRunPlaneEff ; // Evaluate Plane Efficiency + // PID + AliESDpid* fESDpid; // PID object + // New members needed in order to split Run method // into InitRun,RunEvent,FinishRun methods AliESDEvent* fesd; //! Pointer to the ESD event object @@ -407,7 +410,7 @@ private: Int_t fMaxVMEM; // max VMEM memory, MB static const char* fgkStopEvFName; // filename for stop.event stamp // - ClassDef(AliReconstruction, 46) // class for running the reconstruction + ClassDef(AliReconstruction, 47) // class for running the reconstruction }; #endif diff --git a/STEER/STEERBase/AliExternalTrackParam.cxx b/STEER/STEERBase/AliExternalTrackParam.cxx index 0f5bca4f1d5..950f1d51908 100644 --- a/STEER/STEERBase/AliExternalTrackParam.cxx +++ b/STEER/STEERBase/AliExternalTrackParam.cxx @@ -418,7 +418,7 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx // "xTimesRho" - is the product length*density (g/cm^2). // It should be passed as negative when propagating tracks // from the intreaction point to the outside of the central barrel. - // "mass" - the mass of this particle (GeV/c^2). + // "mass" - the mass of this particle (GeV/c^2). Negative mass means charge=2 particle // "dEdx" - mean enery loss (GeV/(g/cm^2) // "anglecorr" - switch for the angular correction //------------------------------------------------------------------ @@ -439,6 +439,7 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx } Double_t p=GetP(); + if (mass<0) p += p; // q=2 particle Double_t p2=p*p; Double_t beta2=p2/(p2 + mass*mass); @@ -454,6 +455,7 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx double lt = 1+0.038*TMath::Log(TMath::Abs(xOverX0)); if (lt>0) theta2 *= lt*lt; } + if (mass<0) theta2 *= 4; // q=2 particle if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE; cC22 = theta2*((1.-fP2)*(1.+fP2))*(1. + fP3*fP3); cC33 = theta2*(1. + fP3*fP3)*(1. + fP3*fP3); @@ -467,7 +469,6 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx Double_t dE=dEdx*xTimesRho; Double_t e=TMath::Sqrt(p2 + mass*mass); if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much! - //cP4 = (1.- e/p2*dE); if ( (1.+ dE/p2*(dE + 2*e)) < 0. ) return kFALSE; cP4 = 1./TMath::Sqrt(1.+ dE/p2*(dE + 2*e)); //A precise formula by Ruben ! if (TMath::Abs(fP4*cP4)>100.) return kFALSE; //Do not track below 10 MeV/c @@ -502,16 +503,20 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterial // "xTimesRho" - is the product length*density (g/cm^2). // It should be passed as negative when propagating tracks // from the intreaction point to the outside of the central barrel. - // "mass" - the mass of this particle (GeV/c^2). + // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2 // "anglecorr" - switch for the angular correction // "Bethe" - function calculating the energy loss (GeV/(g/cm^2)) //------------------------------------------------------------------ Double_t bg=GetP()/mass; - if (bg=kSPECIESC) return 0; - return TMath::Nint(fgkParticleMass[iType]/fgkParticleMassZ[iType]); - } static const char* ParticleName(Int_t iType) {return fgkParticleName[iType];}; static const char* ParticleShortName(Int_t iType) @@ -100,12 +99,13 @@ class AliPID : public TObject { static /*const*/ Float_t fgkParticleMass[kSPECIESCN+1]; // particle masses static /*const*/ Float_t fgkParticleMassZ[kSPECIESCN+1]; // particle masses/charge + static /*const*/ Char_t fgkParticleCharge[kSPECIESCN+1]; // particle charge (in e units!) static const char* fgkParticleName[kSPECIESCN+1]; // particle names static const char* fgkParticleShortName[kSPECIESCN+1]; // particle names static const char* fgkParticleLatexName[kSPECIESCN+1]; // particle names static const Int_t fgkParticleCode[kSPECIESCN+1]; // particle codes - ClassDef(AliPID, 3) // particle id probability densities + ClassDef(AliPID, 4) // particle id probability densities }; diff --git a/TOF/AliTOFtrack.cxx b/TOF/AliTOFtrack.cxx index 64bd058b1fc..eee16aaf28b 100644 --- a/TOF/AliTOFtrack.cxx +++ b/TOF/AliTOFtrack.cxx @@ -66,7 +66,7 @@ AliTOFtrack::AliTOFtrack(const AliESDtrack& t) : // SetLabel(t.GetLabel()); SetChi2(0.); - SetMass(t.GetMass()); + SetMass(t.GetMassForTracking()); Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance()); diff --git a/TPC/Rec/AliTPCtrack.cxx b/TPC/Rec/AliTPCtrack.cxx index 5cd0ec68f50..5398c2b5bad 100644 --- a/TPC/Rec/AliTPCtrack.cxx +++ b/TPC/Rec/AliTPCtrack.cxx @@ -142,7 +142,7 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t, TTreeSRedirector *pcstream) : const Double_t kmaxC[4]={10,10,0.1,0.1}; // cuts on the rms /fP0,fP1,fP2,fP3 SetNumberOfClusters(t.GetTPCclusters(fIndex)); SetLabel(t.GetLabel()); - SetMass(t.GetMass()); + SetMass(t.GetMassForTracking()); for (Int_t i=0; i<4;i++) fPoints[i]=0.; for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.; for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0; diff --git a/TRD/AliTRDtrackV1.cxx b/TRD/AliTRDtrackV1.cxx index a13ddaa294a..aceb18107cc 100644 --- a/TRD/AliTRDtrackV1.cxx +++ b/TRD/AliTRDtrackV1.cxx @@ -132,7 +132,7 @@ AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack() SetLabel(t.GetLabel()); SetChi2(0.0); - SetMass(t.GetMass(kTRUE)); + SetMass(t.GetMassForTracking()); AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls()); Int_t ti[]={-1, -1, -1, -1, -1, -1}; t.GetTRDtracklets(&ti[0]); for(int ip=0; ip=3){ const char *pn[] = {"rho", "x/X0", "", "", "L", "", "Nb"}; printf("D-AliTRDtrackerV1::PropagateTo(): x[%6.2f] bg[%6.2f]\n", xpos, bg); -- 2.39.3