move momentum correction to track candidates and code clean up

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx
index 0b0a70fdc0d81d9bd615d92d46beb71935c0ef10..247aa7d8966f736cd2271b19f502059388320080 100644 (file)
--- a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx
+++ b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx
@@ -98,7 +98,7 @@ AliAnalysisTaskB2::AliAnalysisTaskB2()
 , fESDevent(0)
 , fOutputContainer(0)
 , fHistoMap(0)
-, fESDtrackCuts(0)
+, fTrackCuts(0)
 , fLnID(0)
 , fMaxNSigmaITS(3)
 , fMaxNSigmaTPC(3)
@@ -163,7 +163,7 @@ AliAnalysisTaskB2::AliAnalysisTaskB2(const char* name)
 , fESDevent(0)
 , fOutputContainer(0)
 , fHistoMap(0)
-, fESDtrackCuts(0)
+, fTrackCuts(0)
 , fLnID(0)
 , fMaxNSigmaITS(3)
 , fMaxNSigmaTPC(3)
@@ -238,14 +238,14 @@ void AliAnalysisTaskB2::CreateOutputObjects()
 //
 // Create output objects
 //
-       if(fHistoMap == 0) exit(1); // should be created somewhere else
+       if(fHistoMap == 0) AliFatal("no histogram map"); // should be created somewhere else
        
        fOutputContainer = new TList();
        fOutputContainer->SetOwner(kTRUE);
        
        TObjString* key;
        TIter iter(fHistoMap->GetMap());
-       while((key = (TObjString*)iter.Next())) fOutputContainer->Add((TH1*)fHistoMap->Get(key));
+       while((key = dynamic_cast<TObjString*>(iter.Next()))) fOutputContainer->Add((TH1*)fHistoMap->Get(key));
        
        PostData(1, fOutputContainer);
 }
@@ -257,7 +257,7 @@ AliAnalysisTaskB2::~AliAnalysisTaskB2()
 //
        delete fHistoMap;
        delete fOutputContainer;
-       delete fESDtrackCuts;
+       delete fTrackCuts;
        delete fLnID;
        
        delete fTrigAna;
@@ -287,17 +287,9 @@ void AliAnalysisTaskB2::Exec(Option_t* )
                return;
        }
        
-       if(fESDtrackCuts == 0)
-       {
-               this->Error("Exec", "ESD track cuts not set");
-               return;
-       }
+       if(fTrackCuts == 0) AliFatal("track cuts not set");
        
-       if(fLnID == 0)
-       {
-               this->Error("Exec", "PID not set");
-               return;
-       }
+       if(fLnID == 0) AliFatal("PID not set");
        
        if(fSimulation)
        {
@@ -310,7 +302,7 @@ void AliAnalysisTaskB2::Exec(Option_t* )
                if(fMCevent == 0) return;
        }
        
-       // --------- multiplicity and centrality ------------------
+       // multiplicity and centrality
        
        fNtrk = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent, AliESDtrackCuts::kTrackletsITSTPC, fMaxEta);
        if(fSimulation) fNch = this->GetChargedMultiplicity(fMaxEta);
@@ -337,7 +329,7 @@ void AliAnalysisTaskB2::Exec(Option_t* )
                ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Scaled_Mult"))->Fill(v0ScaMult);
        }
        
-       // ----------------- trigger --------------------
+       // trigger
        
        AliInputEventHandler* eventH = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
        if(eventH == 0)
@@ -360,7 +352,7 @@ void AliAnalysisTaskB2::Exec(Option_t* )
                if(fNtrkMultTrigger) fTriggerFired = ( fTriggerFired && fMultTriggerFired );
        }
        
-       // --------------------- vertex --------------
+       // vertex
        
        fGoodVertex  = kFALSE;
        
@@ -385,11 +377,11 @@ void AliAnalysisTaskB2::Exec(Option_t* )
        if( (vtx->GetY() <= fMinVy) || (vtx->GetY() >= fMaxVy) ) fGoodVertex=kFALSE;
        if( (vtx->GetZ() <= fMinVz) || (vtx->GetZ() >= fMaxVz) ) fGoodVertex=kFALSE;
        
-       // -------------------- pile up ------------------
+       // pile up
        
        fPileUpEvent = fESDevent->IsPileupFromSPDInMultBins();
        
-       // ---------------- event stats ----------------
+       // event stats
        
        ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(0); // number of events
        
@@ -424,13 +416,9 @@ void AliAnalysisTaskB2::Exec(Option_t* )
                }
        }
        
-       // ------------- Particles and tracks for this event ---------------
+       // Particles and tracks for this event
        
-       Int_t nParticles = 0;
-       if(fSimulation)
-       {
-               nParticles = this->GetParticles();
-       }
+       if(fSimulation) this->GetParticles();
        
        this->GetTracks();
        
@@ -462,7 +450,7 @@ Int_t AliAnalysisTaskB2::GetParticles()
                
                if(pid != fPartCode) continue;
                
-               // ------- physical primaries ------------
+               // physical primaries
                
                if(!stack->IsPhysicalPrimary(i)) continue;
                
@@ -475,7 +463,7 @@ Int_t AliAnalysisTaskB2::GetParticles()
                Double_t genPhi = iParticle->Phi();
                Double_t genEta = iParticle->Eta();
                
-               // ------- multiplicity and centrality -------------
+               // multiplicity and centrality
                
                if( fHeavyIons && !fCentTriggerFired) continue;
                if( fNtrkMultTrigger && !fMultTriggerFired) continue;
@@ -488,14 +476,14 @@ Int_t AliAnalysisTaskB2::GetParticles()
                ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_PtY"))->Fill(genY, genPt);
                ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_EtaY"))->Fill(genY, genEta);
                
-               // ------ is within phase space? ----------
+               // is within phase space?
                
                if(TMath::Abs(genY) >= fMaxY) continue;
                
                ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_P"))->Fill(genP);
                ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_Pt"))->Fill(genPt);
                
-               // ------- is from a triggering event? (same as rec.) --------
+               // is from a triggering event? (same as rec.)
                
                if(!fTriggerFired)
                {
@@ -505,7 +493,7 @@ Int_t AliAnalysisTaskB2::GetParticles()
                
                ((TH1D*)fHistoMap->Get(particle + "_Gen_Trig_Prim_Pt"))->Fill(genPt);
                
-               // ------- is from a triggering event with good vertex? -------------
+               // is from a triggering event with good vertex?
                
                if(!fESDevent->GetPrimaryVertex()->GetStatus())
                {
@@ -518,7 +506,7 @@ Int_t AliAnalysisTaskB2::GetParticles()
                ((TH1D*)fHistoMap->Get(particle + "_Gen_Nch"))->Fill(fNch);
                ((TH1D*)fHistoMap->Get(particle + "_Gen_Vtx_Prim_Pt"))->Fill(genPt);
                
-               // ------ is within the geometrical acceptance? ------------
+               // is within the geometrical acceptance?
                
                if(TMath::Abs(genEta) >= fMaxEta) continue;
                
@@ -541,13 +529,13 @@ Int_t AliAnalysisTaskB2::GetTracks()
        
        Int_t nTracks = 0;
        
-       // --------- trigger, vertex and pile-up -----------
+       // trigger, vertex and pile-up
        
        if(!fTriggerFired) return 0;
        if(!fGoodVertex)   return 0;
        if(fPileUpEvent)   return 0;
        
-       // ------------ this is a 'good' event -------------
+       // this is a 'good' event
        
        ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(2); // analyzed events
        
@@ -574,7 +562,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Scaled_Mult"))->Fill(v0ScaMult);
        }
        
-       // ------- track loop --------
+       // track loop
        
        for(Int_t i = 0; i < fESDevent->GetNumberOfTracks(); ++i)
        {
@@ -582,77 +570,41 @@ Int_t AliAnalysisTaskB2::GetTracks()
                
                if(!iTrack) continue;
                
-               Float_t sign = 1;
                TString particle = fSpecies;
-               if(iTrack->GetSign()<0 )
+               if(iTrack->GetSign() < 0) particle.Prepend("Anti");
+               
+               // impact parameters
+               Float_t dcaxy, dcaz;
+               iTrack->GetImpactParameters(dcaxy, dcaz);
+               
+               if(iTrack->GetSign() < 0) // in case of asymmetry
                {
-                       particle.Prepend("Anti");
-                       sign = -1;
+                       dcaxy = -dcaxy;
+                       dcaz  = -dcaz;
                }
                
-               // -------- track cuts ------------------------
+               Double_t nSigmaVtx = fTrackCuts->GetSigmaToVertex(iTrack);
+               
+               // momentum at DCA
+               Double_t p  = iTrack->GetP();
+               Double_t pt = iTrack->Pt();
+               Double_t pz = iTrack->Pz();
+               
+               // track cuts
                
                Double_t eta = iTrack->Eta();
                Double_t phi = this->GetPhi(iTrack);
                
                ((TH2D*)fHistoMap->Get(particle + "_Before_Phi_Eta"))->Fill(eta,phi);
                
-               if(!fESDtrackCuts->AcceptTrack(iTrack)) continue;  // with next track
+               if(!fTrackCuts->AcceptTrack(iTrack)) continue;  // with next track
                if(fTOFmatch && !this->AcceptTOFtrack(iTrack)) continue; // with next track
                
-               // --------------------- end track cuts -----------------------
-               
-               ((TH2D*)fHistoMap->Get(particle + "_After_Phi_Eta"))->Fill(eta, phi);
+               // end track cuts
                
                ++nTracks;
                
-               // charge
-               
-               Double_t z = 1;
-               if(fPartCode>AliPID::kTriton)  z = 2;
-               
-               // impact parameters
-               
-               Float_t dcaxy, dcaz;
-               iTrack->GetImpactParameters(dcaxy, dcaz);
-               
-               dcaxy *= sign;
-               dcaz  *= sign;
-               
-               Double_t nSigmaVtx = fESDtrackCuts->GetSigmaToVertex(iTrack);
-               
-               // momentum
-               
-               Double_t p       = iTrack->GetP()*z; // p at DCA
-               Double_t pt      = iTrack->Pt()*z; // pt at DCA
-               Double_t y       = this->GetRapidity(iTrack, fPartCode);
-               Double_t pITS    = this->GetITSmomentum(iTrack);
-               Double_t pTPC    = iTrack->GetTPCmomentum();
-               Double_t pTOF    = this->GetTOFmomentum(iTrack);
-               Double_t dEdxITS = iTrack->GetITSsignal();
-               Double_t dEdxTPC = iTrack->GetTPCsignal();
-               Int_t nPointsITS = this->GetITSnPointsPID(iTrack);
-               Int_t nPointsTPC = iTrack->GetTPCsignalN();
-               
-               if(fMomentumCorrection)
-               {
-                       pt += this->GetMomentumCorrection(pt);
-                       p   = TMath::Sqrt(pt*pt + iTrack->Pz()*iTrack->Pz());
-                       y   = this->GetRapidity(p, iTrack->Pz(), fPartCode);
-               }
-               
-               Double_t beta = 0;
-               Double_t mass = 0;
-               Double_t m2   = 0;
-               Double_t dm2  = -100;
-               Double_t t    = 0;
-               Double_t dt   = -1000;
-               
-               Double_t simPt  = 0;
-               Double_t simPhi = 0;
-               Double_t simY   = 0;
-               
-               // --------- track cuts ------------
+               ((TH2D*)fHistoMap->Get(particle + "_After_Phi_Eta"))->Fill(eta, phi);
                
                ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAxy"))->Fill(dcaxy);
                ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_DCAz"))->Fill(dcaz);
@@ -666,27 +618,37 @@ Int_t AliAnalysisTaskB2::GetTracks()
                ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2PerCls"))->Fill(iTrack->GetTPCchi2()/iTrack->GetTPCNcls());
                ((TH1D*)fHistoMap->Get(particle + "_TrackCuts_TPCchi2Global"))->Fill(iTrack->GetChi2TPCConstrainedVsGlobal(fESDevent->GetPrimaryVertex()));
                
-               // -------------
+               // detector signals
+               
+               Double_t pITS    = this->GetITSmomentum(iTrack);
+               Double_t pTPC    = iTrack->GetTPCmomentum();
+               Double_t pTOF    = this->GetTOFmomentum(iTrack);
+               Double_t dEdxITS = iTrack->GetITSsignal();
+               Double_t dEdxTPC = iTrack->GetTPCsignal();
+               Int_t nPointsITS = this->GetITSnPointsPID(iTrack);
+               Int_t nPointsTPC = iTrack->GetTPCsignalN();
+               Double_t beta    = 0;
+               Double_t mass    = 0;
+               Double_t m2      = 0;
                
                ((TH2D*)fHistoMap->Get(particle + "_ITS_dEdx_P"))->Fill(pITS, dEdxITS);
                ((TH2D*)fHistoMap->Get(particle + "_TPC_dEdx_P"))->Fill(pTPC, dEdxTPC);
                
-               if(this->TOFmatch(iTrack))
+               if(fTOFmatch)
                {
                        beta = this->GetBeta(iTrack);
-                       m2   = this->GetMassSquared(iTrack);
+                       m2   = this->GetMassSquared(p, beta);
                        mass = TMath::Sqrt(TMath::Abs(m2));
-                       dm2  = this->GetM2Difference(beta, pTOF, AliPID::ParticleMass(fPartCode));
-                       t    = this->GetTimeOfFlight(iTrack)*1.e-3; // ns
-                       dt   = t - 1.e-3*this->GetExpectedTime(iTrack, AliPID::ParticleMass(fPartCode));
                        
                        ((TH2D*)fHistoMap->Get(particle + "_TOF_Beta_P"))->Fill(pTOF, beta);
                        ((TH2D*)fHistoMap->Get(particle + "_TOF_Mass_P"))->Fill(pTOF, mass);
                }
                
-               // -----------------------------------------------
-               //  for pid and efficiency
-               // -----------------------------------------------
+               // for pid and efficiency
+               
+               Double_t simPt  = 0;
+               Double_t simPhi = 0;
+               Double_t simY   = 0;
                
                Int_t simpid = -1;
                TParticle* iParticle = 0;
@@ -710,7 +672,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                                
                                ((TH2D*)fHistoMap->Get(simparticle + "_Sim_PtY"))->Fill(simY, simPt);
                                
-                               if(TMath::Abs(y) < fMaxY)
+                               if(TMath::Abs(simY) < fMaxY)
                                {
                                        ((TH2D*)fHistoMap->Get(simparticle + "_Response_Matrix"))->Fill(pt, simPt);
                                        
@@ -737,7 +699,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                                                        
                                                        ((TH2D*)fHistoMap->Get(simparticle + "_Prim_DiffPt_RecPt"))->Fill(pt,simPt-pt);
                                                        
-                                                       if( this->TOFmatch(iTrack) )
+                                                       if(fTOFmatch)
                                                        {
                                                                ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_P"))->Fill(pTOF, m2);
                                                                ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_Pt"))->Fill(pt, m2);
@@ -749,7 +711,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                                                        
                                                        ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
                                                }
-                                               else // if(this->IsFromMaterial(iParticle)
+                                               else
                                                {
                                                        ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Mat_Pt"))->Fill(simPt);
                                                        
@@ -785,8 +747,8 @@ Int_t AliAnalysisTaskB2::GetTracks()
                
                if( fSimulation )
                {
-                       Int_t sim = simpid > AliPID::kProton ? simpid - offset : simpid;
-                       Int_t rec = pid > AliPID::kProton ? pid - offset : pid;
+                       Int_t sim = (simpid > AliPID::kProton) ? simpid - offset : simpid;
+                       Int_t rec = (pid > AliPID::kProton) ? pid - offset : pid;
                        
                        if((sim > -1) && (rec > -1)) // pid performance
                        {
@@ -802,12 +764,14 @@ Int_t AliAnalysisTaskB2::GetTracks()
                // for bayes iteration
                if(pid != -1)
                {
-                       Int_t iBin = pid > AliPID::kProton ? pid - offset : pid;
+                       Int_t iBin = (pid > AliPID::kProton) ? pid - offset : pid;
                        ((TH1D*)fHistoMap->Get(fSpecies + "_Stats_PID"))->Fill(iBin);
                }
                
                if(pid != fPartCode) continue;
                
+               // track candidates
+               
                Bool_t goodPid = 0;
                
                if(fSimulation)
@@ -822,7 +786,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                ((TH2D*)fHistoMap->Get(particle + "_PID_ITSdEdx_P"))->Fill(pITS, dEdxITS);
                ((TH2D*)fHistoMap->Get(particle + "_PID_TPCdEdx_P"))->Fill(pTPC, dEdxTPC);
                
-               if(this->TOFmatch(iTrack))
+               if(fTOFmatch)
                {
                        ((TH2D*)fHistoMap->Get(particle + "_PID_Beta_P"))->Fill(pTOF, beta);
                        ((TH2D*)fHistoMap->Get(particle + "_PID_Mass_P"))->Fill(pTOF, mass);
@@ -832,9 +796,30 @@ Int_t AliAnalysisTaskB2::GetTracks()
                        }
                }
                
-               // ---------------------------------
+               //
                //  results in |y| < fMaxY
-               // ---------------------------------
+               //
+               
+               if(fMomentumCorrection)
+               {
+                       pt += this->GetMomentumCorrection(pt);
+                       p   = TMath::Sqrt(pt*pt + pz*pz);
+                       
+                       if(fTOFmatch)
+                       {
+                               m2   = this->GetMassSquared(p, beta);
+                               mass = TMath::Sqrt(TMath::Abs(m2));
+                       }
+               }
+               
+               if(fPartCode > AliPID::kTriton)
+               {
+                       p  *= 2.;
+                       pt *= 2.;
+                       pz *= 2.;
+               }
+               
+               Double_t y = this->GetRapidity(p, pz, AliPID::ParticleMass(fPartCode));
                
                ((TH1D*)fHistoMap->Get(particle + "_PID_Y"))->Fill(y);
                ((TH2D*)fHistoMap->Get(particle + "_PID_Pt_Y"))->Fill(y, pt);
@@ -859,8 +844,12 @@ Int_t AliAnalysisTaskB2::GetTracks()
                        if(iTrack->IsOn(AliESDtrack::kTOFout)) ((TH1D*)fHistoMap->Get(particle + "_PID_TOFin_TOFout_Pt"))->Fill(pt);
                }
                
-               if( this->TOFmatch(iTrack) )
+               if(fTOFmatch)
                {
+                       Double_t dm2 = this->GetM2Difference(beta, pTOF, AliPID::ParticleMass(fPartCode));
+                       Double_t t   = this->GetTimeOfFlight(iTrack)*1.e-3; // ns
+                       Double_t dt  = t - this->GetExpectedTime(iTrack, AliPID::ParticleMass(fPartCode))*1.e-3;
+                       
                        ((TH2D*)fHistoMap->Get(particle + "_PID_M2_Pt"))->Fill(pt, m2);
                        ((TH2D*)fHistoMap->Get(particle + "_PID_DM2_Pt"))->Fill(pt, dm2);
                        ((TH2D*)fHistoMap->Get(particle + "_PID_Time_Pt"))->Fill(pt, t);
@@ -874,14 +863,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                
                if( fTOFmatch && (fLnID->GetPidProcedure() > AliLnID::kMaxLikelihood))
                {
-                       if( pt < 1.8 )
-                       {
-                               if ((m2 < fMinM2) || (m2 >= fMaxM2)) m2match = kFALSE;
-                       }
-                       else
-                       {
-                               if((dt < -0.6) || (dt >= 1)) m2match = kFALSE;
-                       }
+                       if ((m2 < fMinM2) || (m2 >= fMaxM2)) m2match = kFALSE;
                }
                
                if(m2match)
@@ -898,7 +880,7 @@ Int_t AliAnalysisTaskB2::GetTracks()
                        {
                                ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Pt"))->Fill(simPt);
                                
-                               if( this->TOFmatch(iTrack) )
+                               if(fTOFmatch)
                                {
                                        ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_M2_Pt"))->Fill(pt, m2);
                                }
@@ -989,20 +971,12 @@ Bool_t AliAnalysisTaskB2::IsMB(UInt_t triggerBits) const
        return ( (triggerBits&AliVEvent::kMB) == AliVEvent::kMB );
 }
 
-Bool_t AliAnalysisTaskB2::TOFmatch(const AliESDtrack* trk) const
-{
-//
-// Check for TOF match signal
-//
-       return ( trk->IsOn(AliESDtrack::kTOFout) && trk->IsOn(AliESDtrack::kTIME) );
-}
-
 Bool_t AliAnalysisTaskB2::AcceptTOFtrack(const AliESDtrack* trk) const
 {
 //
 // Additional checks for TOF match signal
 //
-       if( !this->TOFmatch(trk) ) return kFALSE;
+       if( !fTOFmatch ) return kFALSE;
        if( trk->GetIntegratedLength() < 350) return kFALSE;
        if( trk->GetTOFsignal() < 1e-6) return kFALSE;
        
@@ -1098,22 +1072,6 @@ Double_t AliAnalysisTaskB2::GetTheta(const AliESDtrack* trk) const
        return (pz == 0) ? TMath::PiOver2() : TMath::ACos(pz/p);
 }
 
-Double_t AliAnalysisTaskB2::GetRapidity(const AliESDtrack* trk, Int_t pid) const
-{
-//
-// Rapidity assuming the given particle hypothesis
-// and using the momentum at the DCA
-//
-       Double_t m  = AliPID::ParticleMass(pid);
-       Double_t p  = (pid>AliPID::kTriton) ? 2.*trk->GetP() : trk->GetP();
-       Double_t e  = TMath::Sqrt(p*p + m*m);
-       Double_t pz = (pid>AliPID::kTriton) ? 2.*trk->Pz() : trk->Pz();
-       
-       if(e <= pz) return 1.e+16;
-       
-       return 0.5*TMath::Log( (e+pz)/(e-pz) );
-}
-
 Double_t AliAnalysisTaskB2::GetITSmomentum(const AliESDtrack* trk) const
 {
 //
@@ -1163,23 +1121,19 @@ Double_t AliAnalysisTaskB2::GetExpectedTime(const AliESDtrack* trk, Double_t m)
        return l/beta/2.99792458e-2;
 }
 
-Double_t AliAnalysisTaskB2::GetMassSquared(const AliESDtrack* trk) const
+Double_t AliAnalysisTaskB2::GetMassSquared(Double_t p, Double_t beta) const
 {
 //
 // Mass squared
 //
-       Double_t p = (fPartCode>AliPID::kTriton) ? 2.*this->GetTOFmomentum(trk) : this->GetTOFmomentum(trk);
-       Double_t beta = this->GetBeta(trk);
-       
        return p*p*(1./(beta*beta) - 1.);
 }
 
-Double_t AliAnalysisTaskB2::GetM2Difference(Double_t beta, Double_t momentum, Double_t m) const
+Double_t AliAnalysisTaskB2::GetM2Difference(Double_t beta, Double_t p, Double_t m) const
 {
 //
 // Mass squared difference
 //
-       Double_t p = (fPartCode>AliPID::kTriton) ? 2.*momentum : momentum;
        Double_t expBeta2 = p*p/(p*p+m*m);
        
        return p*p*(1./(beta*beta)-1./expBeta2);
@@ -1315,7 +1269,7 @@ Double_t AliAnalysisTaskB2::GetMomentumCorrection(Double_t ptrec) const
 Double_t AliAnalysisTaskB2::GetRapidity(Double_t p, Double_t pz, Int_t pid) const
 {
 //
-// Rapidity (for momentum correction)
+// Rapidity
 //
        Double_t m  = AliPID::ParticleMass(pid);
        Double_t e  = TMath::Sqrt(p*p + m*m);

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h
index 7c5b2a167d255fbd792eb6e450827e30c32ff3e9..c8b492c795c72e35745a54db209813bd75c3349e 100644 (file)
--- a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h
+++ b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h
@@ -8,7 +8,6 @@
 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
 
 #include <AliAnalysisTask.h>
-#include <AliPID.h>
 
 class AliESDtrack;
 class AliMCEvent;
@@ -61,7 +60,7 @@ class AliAnalysisTaskB2: public AliAnalysisTask
        
        void SetHistogramMap(AliLnHistoMap* map) { fHistoMap = map; }
        
-       void SetESDtrackCuts(AliESDtrackCuts* esdTrackCuts) {fESDtrackCuts = esdTrackCuts; }
+       void SetESDtrackCuts(AliESDtrackCuts* esdTrackCuts) {fTrackCuts = esdTrackCuts; }
        
        void SetPID(AliLnID* lnID) { fLnID = lnID; }
        void SetMaxNSigmaITS(Double_t max) { fMaxNSigmaITS = max; }
@@ -89,8 +88,6 @@ class AliAnalysisTaskB2: public AliAnalysisTask
        Bool_t IsFastOnly(UInt_t triggerBits) const;
        Bool_t IsMB(UInt_t triggerBits) const;
        
-       Bool_t TOFmatch(const AliESDtrack* trk) const;
-       
        TParticle* GetParticle(const AliESDtrack* trk) const;
        
        Int_t GetChargedMultiplicity(Double_t maxEta) const;
@@ -104,12 +101,11 @@ class AliAnalysisTaskB2: public AliAnalysisTask
        
        Double_t GetPhi(const AliESDtrack* trk) const;
        Double_t GetTheta(const AliESDtrack* trk) const;
-       Double_t GetRapidity(const AliESDtrack* trk, Int_t pid) const;
        Double_t GetRapidity(Double_t p, Double_t pz, Int_t pid) const;
        Double_t GetITSmomentum(const AliESDtrack* trk) const;
        Double_t GetTOFmomentum(const AliESDtrack* trk) const;
        Double_t GetBeta(const AliESDtrack* trk) const;
-       Double_t GetMassSquared(const AliESDtrack* trk) const;
+       Double_t GetMassSquared(Double_t p, Double_t beta) const;
        Double_t GetTimeOfFlight(const AliESDtrack* trk) const;
        Double_t GetITSchi2PerCluster(const AliESDtrack* trk) const;
        Int_t    GetITSnClusters(const AliESDtrack* trk) const;
@@ -172,7 +168,7 @@ class AliAnalysisTaskB2: public AliAnalysisTask
        
        TList* fOutputContainer; // output container
        AliLnHistoMap*  fHistoMap; // histogram map (defined somewhere else)
-       AliESDtrackCuts* fESDtrackCuts; // track cuts (defined somewhere else)
+       AliESDtrackCuts* fTrackCuts; // track cuts (defined somewhere else)
        AliLnID* fLnID; // PID for light nuclei (defined somewhere else)
        
        Double_t fMaxNSigmaITS; // maximum number of sigmas to dEdx in the ITS

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx
index 3a139ba066c977a3236de2f9be4b71e5771aff49..9daceb54910a6aa005031ef7c7e8850210298a73 100644 (file)
--- a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx
+++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx
@@ -24,6 +24,7 @@
 #include <AliITSPIDResponse.h>
 #include <AliTPCPIDResponse.h>
 #include <TPDGCode.h>
+#include <AliAODMCParticle.h>
 #include "AliLnID.h"
 
 ClassImp(AliLnID)
@@ -178,10 +179,26 @@ Int_t AliLnID::GetPID(const TParticle* p) const
 {
 //
 // Montecarlo PID
+//
+       return this->GetPID(p->GetPdgCode());
+}
+
+Int_t AliLnID::GetPID(const AliAODMCParticle* p) const
+{
+//
+// Montecarlo PID
+//
+       return this->GetPID(p->GetPdgCode());
+}
+
+Int_t AliLnID::GetPID(Int_t pdgCode) const
+{
+//
+// Montecarlo PID
 //
        enum { kDeuteron=1000010020, kTriton=1000010030, kHelium3=1000020030, kAlpha=1000020040 };
        
-       switch(TMath::Abs(p->GetPdgCode()))
+       switch(TMath::Abs(pdgCode))
        {
                case kElectron:  return AliPID::kElectron;
                case kMuonMinus: return AliPID::kMuon;
@@ -210,6 +227,10 @@ Int_t AliLnID::GetPID(Int_t pidCode, Double_t pITS, Double_t dEdxITS, Int_t nPoi
        {
                return this->GetMaxLikelihoodPID(pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta);
        }
+       else if(fPidProcedure == kITS)
+       {
+               return this->GetITSpid(pidCode, pITS, dEdxITS, nPointsITS, nSigITS);
+       }
        else if(fPidProcedure == kTPC)
        {
                return this->GetTPCpid(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigTPC);
@@ -226,6 +247,17 @@ Int_t AliLnID::GetPID(Int_t pidCode, Double_t pITS, Double_t dEdxITS, Int_t nPoi
        return -1;
 }
 
+Int_t AliLnID::GetITSpid(Int_t pidCode, Double_t pITS, Double_t dEdxITS, Double_t nPointsITS, Double_t nSigmaITS) const
+{
+//
+// Check if particle with the given pid code is within
+// +/- nSigma around the expected dEdx in the ITS
+//
+       if( this->GetITSmatch(pidCode, pITS, dEdxITS, static_cast<Int_t>(nPointsITS), nSigmaITS)) return pidCode;
+       
+       return -1;
+}
+
 Int_t AliLnID::GetTPCpid(Int_t pidCode, Double_t pTPC, Double_t dEdxTPC, Double_t nPointsTPC, Double_t nSigmaTPC) const
 {
 //

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h
index eae19a4d30db48b6dac5b97e9fe92ad029642265..c5bfdc40ff131249f4f2dec856ddc73cecef5f7e 100644 (file)
--- a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h
+++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h
@@ -11,6 +11,7 @@
 #include <AliPID.h>
 
 class TParticle;
+class AliAODMCParticle;
 class AliITSPIDResponse;
 class AliTPCPIDResponse;
 
@@ -24,9 +25,13 @@ class AliLnID: public TObject
        virtual ~AliLnID();
        
        Int_t GetPID(const TParticle* p) const;
+       Int_t GetPID(const AliAODMCParticle* p) const;
+       Int_t GetPID(Int_t pdgCode) const;
        
        Int_t GetPID(Int_t partCode, Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta, Double_t nSigITS=3, Double_t nSigTPC=3, Double_t nSigTOF=3) const;
        
+       Int_t GetITSpid(Int_t partCode, Double_t pITS, Double_t dEdx, Double_t nPoints, Double_t nSigma=3) const;
+       
        Int_t GetTPCpid(Int_t partCode, Double_t pTPC, Double_t dEdx, Double_t nPoints, Double_t nSigma=3) const;
        
        Int_t GetITSTPCpid(Int_t partCode, Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t nSigmaITS, Double_t pTPC, Double_t dEdxTPC, Double_t nPointsTPC, Double_t nSigmaTPC) const;
@@ -66,7 +71,7 @@ class AliLnID: public TObject
        void SetTPCChargeCorrection(Double_t zexp) { fZexp = zexp; }
        
        enum { kSPECIES = 9 };
-       enum { kBayes=0, kMaxLikelihood, kTPC, kITSTPC, kTPCTOF };
+       enum { kBayes=0, kMaxLikelihood, kITS, kTPC, kITSTPC, kTPCTOF };
        
   private: