--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// Analysis task for B2
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TString.h>
+
+#include <AliLog.h>
+
+#include <AliAnalysisTask.h>
+#include <AliAnalysisManager.h>
+
+#include <AliESD.h>
+#include <AliESDEvent.h>
+#include <AliESDInputHandler.h>
+#include <AliESDtrackCuts.h>
+#include <AliExternalTrackParam.h>
+
+#include <AliMCEvent.h>
+#include <AliMCEventHandler.h>
+#include <AliMCVertex.h>
+#include <AliStack.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TMCProcess.h>
+
+#include <AliTriggerAnalysis.h> // for offline signals
+#include <AliCentrality.h>
+#include <AliESDUtils.h>
+
+#include <AliESDpid.h>
+
+#include <TH1D.h>
+#include <TH2D.h>
+
+#include "AliLnID.h"
+#include "AliLnHistoMap.h"
+#include "AliAnalysisTaskB2.h"
+
+ClassImp(AliAnalysisTaskB2)
+
+AliAnalysisTaskB2::AliAnalysisTaskB2()
+: AliAnalysisTask()
+, fSpecies("Proton")
+, fPartCode(AliPID::kProton)
+
+, fHeavyIons(0)
+, fSimulation(0)
+
+, fMultTrigger(0)
+, fCentTrigger(0)
+, fTriggerFired(0)
+, fGoodVertex(0)
+, fPileUpEvent(0)
+
+, fV0AND(0)
+, fNoFastOnly(0)
+, fMinKNOmult(-10)
+, fMaxKNOmult(100000)
+, fMinCentrality(0)
+, fMaxCentrality(100)
+
+, fNtrk(0)
+, fMeanNtrk(1)
+, fKNOmult(1)
+
+, fMinVx(-1)
+, fMaxVx(1)
+, fMinVy(-1)
+, fMaxVy(1)
+, fMinVz(-10)
+, fMaxVz(10)
+
+, fMinDCAxy(-1)
+, fMaxDCAxy(1)
+, fMinDCAz(-2)
+, fMaxDCAz(2)
+, fMaxNSigma(3)
+
+, fMinEta(-0.8)
+, fMaxEta(0.8)
+, fMinY(-0.5)
+, fMaxY(0.5)
+
+, fMCevent(0)
+, fESDevent(0)
+
+, fHistoMap(0)
+, fESDtrackCuts(0)
+, fLnID(0)
+
+, fMaxNSigmaITS(3)
+, fMaxNSigmaTPC(3)
+, fMaxNSigmaTOF(3)
+
+, fTrigAna(0)
+, fESDpid(0)
+, fIsPidOwner(0)
+, fTimeZeroType(AliESDpid::kTOF_T0)
+, fTOFmatch(0)
+
+{
+//
+// Default constructor
+//
+ AliLog::SetGlobalLogLevel(AliLog::kFatal);
+
+ fTrigAna = new AliTriggerAnalysis();
+}
+
+AliAnalysisTaskB2::AliAnalysisTaskB2(const char* name)
+: AliAnalysisTask(name,"")
+, fSpecies("Proton")
+, fPartCode(AliPID::kProton)
+
+, fHeavyIons(0)
+, fSimulation(0)
+
+, fMultTrigger(0)
+, fCentTrigger(0)
+, fTriggerFired(0)
+, fGoodVertex(0)
+, fPileUpEvent(0)
+
+, fV0AND(0)
+, fNoFastOnly(0)
+, fMinKNOmult(-10)
+, fMaxKNOmult(100000)
+, fMinCentrality(-1)
+, fMaxCentrality(100)
+
+, fNtrk(0)
+, fMeanNtrk(1)
+, fKNOmult(1)
+
+, fMinVx(-1)
+, fMaxVx(1)
+, fMinVy(-1)
+, fMaxVy(1)
+, fMinVz(-10)
+, fMaxVz(10)
+
+, fMinDCAxy(-1)
+, fMaxDCAxy(1)
+, fMinDCAz(-2)
+, fMaxDCAz(2)
+, fMaxNSigma(3)
+
+, fMinEta(-0.8)
+, fMaxEta(0.8)
+, fMinY(-0.5)
+, fMaxY(0.5)
+
+, fMCevent(0)
+, fESDevent(0)
+
+, fHistoMap(0)
+, fESDtrackCuts(0)
+, fLnID(0)
+
+, fMaxNSigmaITS(3)
+, fMaxNSigmaTPC(3)
+, fMaxNSigmaTOF(3)
+
+, fTrigAna(0)
+, fESDpid(0)
+, fIsPidOwner(0)
+, fTimeZeroType(AliESDpid::kTOF_T0)
+, fTOFmatch(0)
+
+{
+//
+// Constructor
+//
+ DefineInput(0, TChain::Class());
+ DefineOutput(0, AliLnHistoMap::Class());
+
+ //kFatal, kError, kWarning, kInfo, kDebug, kMaxType
+ AliLog::SetGlobalLogLevel(AliLog::kFatal);
+
+ fTrigAna = new AliTriggerAnalysis();
+}
+
+void AliAnalysisTaskB2::ConnectInputData(Option_t *)
+{
+//
+// Connect input data
+//
+ using namespace std;
+
+ TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+ if(!tree)
+ {
+ cerr << "ERROR: Could not read chain from input slot 0" << endl;
+ return;
+ }
+
+ // Get the pointer to the existing analysis manager via the static access method.
+ AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ cerr << "ERROR: could not get analysis manager" << endl;
+ return;
+ }
+
+ // cast type AliVEventHandler
+ AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(mgr->GetInputEventHandler());
+ if (!esdH)
+ {
+ cerr << "ERROR: could not get ESDInputHandler" << endl;
+ return;
+ }
+
+ // Get pointer to esd event from input handler
+ fESDevent = esdH->GetEvent();
+
+ // PID object for TOF
+ fESDpid = esdH->GetESDpid();
+
+ if(!fESDpid)
+ { //in case of no Tender attached
+ fESDpid = new AliESDpid();
+ fIsPidOwner = kTRUE;
+ }
+
+ if(!fSimulation) return;
+
+ AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (mgr->GetMCtruthEventHandler());
+ if (!mcH)
+ {
+ cerr << "ERROR: could not get AliMCEventHandler" << endl;
+ return;
+ }
+
+ fMCevent = mcH->MCEvent();
+ if (!fMCevent)
+ {
+ cerr << "ERROR: could not get MC fLnEvent" << endl;
+ return;
+ }
+}
+
+void AliAnalysisTaskB2::CreateOutputObjects()
+{
+//
+// Create output objects
+//
+ if(fHistoMap == 0) exit(1); // should be created somewhere else
+
+ PostData(0, fHistoMap);
+}
+
+AliAnalysisTaskB2::~AliAnalysisTaskB2()
+{
+//
+// Default destructor
+//
+ delete fHistoMap;
+ delete fESDtrackCuts;
+ delete fLnID;
+
+ delete fTrigAna;
+
+ if(fIsPidOwner) delete fESDpid;
+}
+
+void AliAnalysisTaskB2::SetParticleSpecies(const TString& species)
+{
+//
+// set the particle species and the AliPID code
+//
+ fSpecies = species;
+ fPartCode = this->GetPidCode(species);
+}
+
+void AliAnalysisTaskB2::Exec(Option_t* )
+{
+//
+// Execute analysis for the current event
+//
+ using namespace std;
+
+ if (fESDevent == 0)
+ {
+ cerr << "ERROR: AliESDEvent not available" << endl;
+ return;
+ }
+
+ if(fESDtrackCuts == 0)
+ {
+ cerr << "ERROR: ESD track cuts not set" << endl;
+ return;
+ }
+
+ if(fLnID == 0)
+ {
+ cerr << "ERROR: PID not set" << endl;
+ return;
+ }
+
+ fMultTrigger = kFALSE;
+ fCentTrigger = kFALSE;
+ fTriggerFired = kFALSE;
+ fGoodVertex = kFALSE;
+ fPileUpEvent = kFALSE;
+
+ // --------- multiplicity and centrality ------------------
+
+ fNtrk = AliESDtrackCuts::GetReferenceMultiplicity(fESDevent);
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Event_Ntrk"))->Fill(fNtrk);
+
+ fKNOmult = fNtrk/fMeanNtrk;
+
+ if( (fKNOmult >= fMinKNOmult) && (fKNOmult < fMaxKNOmult) ) fMultTrigger = kTRUE;
+
+ if(fHeavyIons)
+ {
+ AliCentrality* esdCent = fESDevent->GetCentrality();
+
+ Float_t centrality = esdCent->GetCentralityPercentile("V0M");
+
+ if((centrality >= fMinCentrality) && (centrality < fMaxCentrality)) fCentTrigger = kTRUE;
+
+ Float_t v0ScaMult;
+ Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Mult"))->Fill(v0Mult);
+ ((TH1D*)fHistoMap->Get(fSpecies + "_V0_Scaled_Mult"))->Fill(v0ScaMult);
+ }
+
+ // ----------------- trigger --------------------
+
+ UInt_t triggerBits = (dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+
+ if(fHeavyIons)
+ {
+ fTriggerFired = ( this->IsMB(triggerBits) && fCentTrigger );
+ }
+ else
+ {
+ fTriggerFired = this->IsMB(triggerBits);
+ if(fNoFastOnly) fTriggerFired = ( fTriggerFired && !this->IsFastOnly(triggerBits) );
+ if(fV0AND) fTriggerFired = ( fTriggerFired && this->IsV0AND() );
+
+ fTriggerFired = ( fTriggerFired && fMultTrigger );
+ }
+
+ // --------------------- vertex --------------
+
+ const AliESDVertex* vtx = fESDevent->GetPrimaryVertex(); // best primary vertex
+
+ if(vtx->GetStatus())
+ {
+ fGoodVertex=kTRUE;
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
+ }
+
+ if(fESDevent->GetPrimaryVertex()->IsFromVertexerZ())
+ {
+ if(fESDevent->GetPrimaryVertex()->GetDispersion() > 0.02) fGoodVertex=kFALSE;
+ if(fESDevent->GetPrimaryVertex()->GetZRes() > 0.25 ) fGoodVertex=kFALSE;
+ }
+
+ if( (vtx->GetX() <= fMinVx) || (vtx->GetX() >= fMaxVx) ) fGoodVertex=kFALSE;
+ if( (vtx->GetY() <= fMinVy) || (vtx->GetY() >= fMaxVy) ) fGoodVertex=kFALSE;
+ if( (vtx->GetZ() <= fMinVz) || (vtx->GetZ() >= fMaxVz) ) fGoodVertex=kFALSE;
+
+ // -------------------- pile up ------------------
+
+ if(fESDevent->IsPileupFromSPDInMultBins()) fPileUpEvent = kTRUE;
+
+ // ---------------- event stats ----------------
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(0); // number of events
+
+ if(fTriggerFired)
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(1); // triggered events
+
+ if(vtx->GetStatus())
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(3); // valid vertex within a triggered event
+
+ if((vtx->GetZ() > fMinVz) && (vtx->GetZ() < fMaxVz))
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(4); // Vz
+
+ if( (vtx->GetX() > fMinVx) && (vtx->GetX() < fMaxVx)
+ && (vtx->GetY() > fMinVy) && (vtx->GetY() < fMaxVy))
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(5); // Vx, Vy
+
+ if(fGoodVertex)
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(6); // VertexerZ
+
+ if(!fPileUpEvent)
+ {
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(7); // no pile-up
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // ------------- Particles and tracks for this event ---------------
+
+ Int_t nParticles = 0;
+ if(fSimulation)
+ {
+ nParticles = this->GetParticles();
+ }
+
+ this->GetTracks();
+
+ // Post the data (input/output slots #0 already used by the base class)
+ PostData(0, fHistoMap);
+}
+
+Int_t AliAnalysisTaskB2::GetParticles()
+{
+//
+// Get particles from current event
+//
+ int nParticles = 0;
+
+ AliStack* stack = fMCevent->Stack();
+ if (!stack)
+ {
+ AliDebug(AliLog::kWarning, "stack not available");
+ return 0;
+ }
+
+ for (Int_t i = 0; i < fMCevent->GetNumberOfTracks(); ++i)
+ {
+ TParticle* iParticle = stack->Particle(i);
+
+ if(!iParticle) continue;
+
+ Int_t pid = fLnID->GetPID(iParticle);
+
+ if(pid != fPartCode) continue;
+
+ // ------- physical primaries ------------
+
+ if(!stack->IsPhysicalPrimary(i)) continue;
+
+ TString particle = fSpecies;
+ if(iParticle->GetPDG()->Charge() < 0) particle.Prepend("Anti");
+
+ Double_t genP = iParticle->P();
+ Double_t genPt = iParticle->Pt();
+ Double_t genY = iParticle->Y();
+ Double_t genPhi = iParticle->Phi();
+ Double_t genEta = iParticle->Eta();
+
+ // --------- all events -----------
+
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_P"))->Fill(genP);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Pt"))->Fill(genPt);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Y"))->Fill(genY);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Phi"))->Fill(genPhi);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Prim_Eta"))->Fill(genEta);
+ ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_PtY"))->Fill(genY, genPt);
+ ((TH2D*)fHistoMap->Get(particle + "_Gen_Prim_EtaY"))->Fill(genY, genEta);
+
+ // ------- multiplicity and centrality -------------
+
+ if( fHeavyIons && !fCentTrigger) continue;
+ if(!fHeavyIons && !fMultTrigger) continue;
+
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Mult_Prim_PtY"))->Fill(genY,genPt);
+
+ // ------ is within phase space? ----------
+
+ if(TMath::Abs(genY) >= fMaxY) continue;
+
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_PhS_Prim_Pt"))->Fill(genPt);
+
+ // ------- is from a triggered event? (same as rec.) --------
+
+ if(!fTriggerFired)
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_NoTrig_Prim_Pt"))->Fill(genPt);
+ continue;
+ }
+
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Trig_Prim_Pt"))->Fill(genPt);
+
+ // ------- is from an triggered event with good vertex? -------------
+
+ if(!fESDevent->GetPrimaryVertex()->GetStatus())
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_NoVtx_Prim_Pt"))->Fill(genPt);
+ }
+
+ if(!fGoodVertex) continue;
+ if(fPileUpEvent) continue;
+
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Vtx_Prim_Pt"))->Fill(genPt);
+
+ // ------ is within the geometrical acceptance? ------------
+
+ if(TMath::Abs(iParticle->Eta()) >= fMaxEta) continue;
+
+ ((TH2D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_PtY"))->Fill(genY,genPt);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_P"))->Fill(genP);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Pt"))->Fill(genPt);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Phi"))->Fill(genPhi);
+ ((TH1D*)fHistoMap->Get(particle + "_Gen_Acc_Prim_Y"))->Fill(genY);
+ }
+
+ return nParticles;
+}
+
+Int_t AliAnalysisTaskB2::GetTracks()
+{
+//
+// Get tracks from current event
+//
+ using namespace std;
+
+ Int_t nTracks = 0;
+
+ // --------- trigger, vertex and pile-up -----------
+
+ if(!fTriggerFired) return 0;
+ if(!fGoodVertex) return 0;
+ if(fPileUpEvent) return 0;
+
+ // ------------ this is a 'good' event -------------
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats"))->Fill(2); // analyzed events
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_Event_Ntrk"))->Fill(fNtrk);
+
+ if(fSimulation)
+ {
+ Double_t nch = this->GetChargedMultiplicity(0.5);
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Event_Nch_Ntrk"))->Fill(fNtrk, nch);
+ }
+
+ const AliESDVertex* vtx = fESDevent->GetPrimaryVertex();
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YX"))->Fill(vtx->GetX(), vtx->GetY());
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_YZ"))->Fill(vtx->GetZ(), vtx->GetY());
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Ana_Vertex_XZ"))->Fill(vtx->GetZ(), vtx->GetX());
+
+ if(fHeavyIons)
+ {
+ Float_t v0ScaMult;
+ Float_t v0Mult = AliESDUtils::GetCorrV0(fESDevent,v0ScaMult);
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Mult"))->Fill(v0Mult);
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Ana_V0_Scaled_Mult"))->Fill(v0ScaMult);
+ }
+
+ // ------- track loop --------
+
+ for(Int_t i = 0; i < fESDevent->GetNumberOfTracks(); ++i)
+ {
+ AliESDtrack* iTrack = fESDevent->GetTrack(i);
+
+ if(!iTrack) continue;
+
+ // -------- track cuts ------------------------
+
+ Double_t theta = this->GetTheta(iTrack);
+ Double_t phi = this->GetPhi(iTrack);
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Before_Phi_Theta"))->Fill(theta,phi);
+
+ if(!fESDtrackCuts->AcceptTrack(iTrack)) continue; // with next track
+ if(fTOFmatch && !this->AcceptTOFtrack(iTrack)) continue; // with next track
+
+ // --------------------- end track cuts -----------------------
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_After_Phi_Theta"))->Fill(theta, phi);
+
+ ++nTracks;
+
+ // charge
+
+ Double_t z = 1;
+ if(fPartCode>AliPID::kTriton) z = 2;
+
+ Float_t sign = 1;
+ TString particle = fSpecies;
+ if(iTrack->GetSign()<0 )
+ {
+ particle.Prepend("Anti");
+ sign = -1;
+ }
+
+ // 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();
+ Double_t pt = iTrack->Pt(); // 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();
+
+ Double_t beta = 0;
+ Double_t mass = 0;
+ Double_t m2 = 0;
+
+ Double_t simPt = 0;
+ Double_t simPhi = 0;
+ Double_t simY = 0;
+
+ // --------- track cuts ------------
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_DCAxy"))->Fill(dcaxy);
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_DCAz"))->Fill(dcaz);
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_NSigma"))->Fill(nSigmaVtx);
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_ITSchi2PerCls"))->Fill(this->GetITSchi2PerCluster(iTrack));
+
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_TPCncls"))->Fill(iTrack->GetTPCNcls());
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_TPCclsOverF"))->Fill((Double_t)iTrack->GetTPCNcls()/(Double_t)iTrack->GetTPCNclsF());
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_TPCxRows"))->Fill(iTrack->GetTPCCrossedRows());
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_TPCchi2PerCls"))->Fill(iTrack->GetTPCchi2()/iTrack->GetTPCNcls());
+ ((TH1D*)fHistoMap->Get(fSpecies + "_TrackCuts_TPCchi2Global"))->Fill(iTrack->GetChi2TPCConstrainedVsGlobal(fESDevent->GetPrimaryVertex()));
+
+ // -------------
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_ITS_dEdx_P"))->Fill(pITS, dEdxITS);
+ ((TH2D*)fHistoMap->Get(fSpecies + "_TPC_dEdx_P"))->Fill(pTPC, dEdxTPC);
+
+ if(this->TOFmatch(iTrack))
+ {
+ beta = this->GetBeta(iTrack);
+ m2 = this->GetMassSquare(iTrack);
+ mass = TMath::Sqrt(TMath::Abs(m2));
+
+ ((TH2D*)fHistoMap->Get(fSpecies + "_TOF_Beta_P"))->Fill(pTOF, beta);
+ ((TH2D*)fHistoMap->Get(fSpecies + "_TOF_Mass_P"))->Fill(pTOF, mass);
+ }
+
+ // -----------------------------------------------
+ // for pid and efficiency
+ // -----------------------------------------------
+
+ Int_t simpid = -1;
+ TParticle* iParticle = 0;
+
+ if(fSimulation)
+ {
+ iParticle = this->GetParticle(iTrack);
+
+ if(iParticle == 0) continue;
+
+ simPt = iParticle->Pt();
+ simPhi = iParticle->Phi();
+ simY = iParticle->Y();
+
+ simpid = fLnID->GetPID(iParticle);
+
+ if(simpid == fPartCode)
+ {
+ TString simparticle = fSpecies;
+ if(this->GetSign(iParticle)<0) simparticle.Prepend("Anti");
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_PtY"))->Fill(simY, simPt);
+
+ if(TMath::Abs(y) < fMaxY)
+ {
+ ((TH2D*)fHistoMap->Get(simparticle + "_Response_Matrix"))->Fill(pt, simPt);
+
+ // for pid
+ if(!this->IsFakeTrack(iTrack))
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Pt"))->Fill(simPt);
+
+ if(this->IsPhysicalPrimary(iParticle)) // the efficiency is calculated on the primaries
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Pt"))->Fill(simPt);
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Y"))->Fill(simY);
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Prim_Phi"))->Fill(simPhi);
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Prim_Response_Matrix"))->Fill(pt, simPt);
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Prim_DiffPt_RecPt"))->Fill(pt,simPt-pt);
+
+ if( this->TOFmatch(iTrack) )
+ {
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_P"))->Fill(pTOF, m2);
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Prim_M2_Pt"))->Fill(pt, m2);
+ }
+
+ }
+ else
+ {
+ if(this->IsFromWeakDecay(iParticle))
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_Pt"))->Fill(simPt);
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ }
+ else // if(this->IsFromMaterial(iParticle)
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Mat_Pt"))->Fill(simPt);
+
+ ((TH2D*)fHistoMap->Get(simparticle + "_Sim_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ }
+ }
+ }
+ else // fake tracks
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Pt"))->Fill(simPt);
+ if(this->IsPhysicalPrimary(iParticle))
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Prim_Pt"))->Fill(simPt);
+ }
+ else if(this->IsFromWeakDecay(iParticle))
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Fdwn_Pt"))->Fill(simPt);
+ }
+ else if(this->IsFromMaterial(iParticle))
+ {
+ ((TH1D*)fHistoMap->Get(simparticle + "_Sim_Fake_Mat_Pt"))->Fill(simPt);
+ }
+ }
+ }
+ }
+ }
+
+ // get the pid
+ Int_t pid = fLnID->GetPID( fPartCode, pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, fMaxNSigmaITS, fMaxNSigmaTPC, fMaxNSigmaTOF);
+
+ // pid table for prior probabilities (only Bayes)
+
+ Int_t offset = AliPID::kDeuteron - 5;
+
+ if( fSimulation )
+ {
+ 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
+ {
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, rec);
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(9, rec);
+ }
+ if(sim > -1)
+ {
+ ((TH2D*)fHistoMap->Get(fSpecies + "_Stats_PID_Table"))->Fill(sim, 9);
+ }
+ }
+
+ // for bayes iteration
+ if(pid != -1)
+ {
+ Int_t iBin = pid > AliPID::kProton ? pid - offset : pid;
+ ((TH1D*)fHistoMap->Get(fSpecies + "_Stats_PID"))->Fill(iBin);
+ }
+
+ if(pid != fPartCode) continue;
+
+ Bool_t goodPid = 0;
+
+ if(fSimulation)
+ {
+ goodPid = ( simpid == pid );
+ }
+
+ // pid performance
+ ((TH2D*)fHistoMap->Get(particle + "_PID_ITSdEdx_P"))->Fill(p*z, dEdxITS);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_TPCdEdx_P"))->Fill(pTPC*z, dEdxTPC);
+
+ if(this->TOFmatch(iTrack))
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_PID_Beta_P"))->Fill(pTOF*z, beta);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_Mass_P"))->Fill(pTOF*z, mass);
+ if(fSimulation && goodPid)
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Mass"))->Fill(mass);
+ }
+ }
+
+ // ---------------------------------
+ // results in |y| < fMaxY
+ // ---------------------------------
+
+ ((TH1D*)fHistoMap->Get(particle + "_PID_Y"))->Fill(y);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_Pt_Y"))->Fill(y, pt);
+
+ if(TMath::Abs(y) >= fMaxY) continue;
+
+ ((TH1D*)fHistoMap->Get(particle + "_PID_Pt"))->Fill(pt);
+ ((TH1D*)fHistoMap->Get(particle + "_PID_Phi"))->Fill(phi);
+
+ if( this->TOFmatch(iTrack) )
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_PID_M2_Pt"))->Fill(pt, m2);
+ }
+
+ // secondaries
+ if(fTOFmatch && fLnID->GetPidProcedure() > AliLnID::kMaxLikelihood)
+ {
+ if(fLnID->GetM2match(fPartCode, pTOF, m2, 3))
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_PID_DCAxy_Pt"))->Fill(pt, dcaxy);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_DCAz_Pt"))->Fill(pt, dcaz);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_NSigma_Pt"))->Fill(pt, nSigmaVtx);
+ }
+ }
+ else
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_PID_DCAxy_Pt"))->Fill(pt, dcaxy);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_DCAz_Pt"))->Fill(pt, dcaz);
+ ((TH2D*)fHistoMap->Get(particle + "_PID_NSigma_Pt"))->Fill(pt, nSigmaVtx);
+ }
+
+ if(fSimulation && goodPid)
+ {
+ // for unfolding and pid contamination
+ if(!this->IsFakeTrack(iTrack))
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Pt"))->Fill(simPt);
+
+ if( this->TOFmatch(iTrack) )
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_M2_Pt"))->Fill(pt,m2);
+ }
+
+ if(this->IsPhysicalPrimary(iParticle))
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Prim_Pt"))->Fill(simPt);
+
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_DCAz_Pt"))->Fill(simPt,dcaz);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Prim_NSigma_Pt"))->Fill(simPt, nSigmaVtx);
+ }
+ else
+ {
+ if(this->IsFromWeakDecay(iParticle))
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_DCAz_Pt"))->Fill(simPt,dcaz);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fdwn_NSigma_Pt"))->Fill(simPt, nSigmaVtx);
+ }
+ else
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_DCAz_Pt"))->Fill(simPt,dcaz);
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Mat_NSigma_Pt"))->Fill(simPt, nSigmaVtx);
+ }
+ }
+ }
+ else // fake tracks
+ {
+ ((TH1D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Pt"))->Fill(simPt);
+
+ if(this->IsPhysicalPrimary(iParticle))
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Prim_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ }
+ else if(this->IsFromWeakDecay(iParticle))
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ }
+ else
+ {
+ ((TH2D*)fHistoMap->Get(particle + "_Sim_PID_Fake_Mat_DCAxy_Pt"))->Fill(simPt,dcaxy);
+ }
+ }
+ }
+ }
+
+ return nTracks;
+}
+
+void AliAnalysisTaskB2::Terminate(Option_t* )
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+}
+
+Bool_t AliAnalysisTaskB2::IsV0AND() const
+{
+//
+// signals in both V0A and V0C
+//
+ return ( fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0A) &&
+ fTrigAna->IsOfflineTriggerFired(fESDevent, AliTriggerAnalysis::kV0C) );
+}
+
+Bool_t AliAnalysisTaskB2::IsFastOnly(UInt_t triggerBits) const
+{
+//
+// kFastOnly trigger
+//
+ return ( (triggerBits&AliVEvent::kFastOnly) == AliVEvent::kFastOnly );
+}
+
+Bool_t AliAnalysisTaskB2::IsMB(UInt_t triggerBits) const
+{
+//
+// MB event
+//
+ 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( trk->GetIntegratedLength() < 350) return kFALSE;
+ if( trk->GetTOFsignal() < 1e-6) return kFALSE;
+
+ return kTRUE;
+}
+
+TParticle* AliAnalysisTaskB2::GetParticle(const AliESDtrack* trk) const
+{
+//
+// Particle that left the track
+//
+ AliStack* stack = fMCevent->Stack();
+
+ Int_t label = TMath::Abs(trk->GetLabel()); // if negative then it shares points from other tracks
+ if( label >= fMCevent->GetNumberOfTracks() ) return 0;
+
+ return stack->Particle(label);
+}
+
+Bool_t AliAnalysisTaskB2::IsFakeTrack(const AliESDtrack* trk) const
+{
+//
+// Check if the track shares some clusters with different particles
+//
+ return ( trk->GetLabel() < 0 );
+}
+
+Bool_t AliAnalysisTaskB2::IsPhysicalPrimary(const TParticle* prt) const
+{
+//
+// Check if the particle is physical primary
+//
+ AliStack* stack = fMCevent->Stack();
+ Int_t index = stack->Particles()->IndexOf(prt);
+
+ return stack->IsPhysicalPrimary(index);
+}
+
+Bool_t AliAnalysisTaskB2::IsFromMaterial(const TParticle* prt) const
+{
+//
+// Check if the particle is originated at the materials
+//
+ AliStack* stack = fMCevent->Stack();
+ Int_t index = stack->Particles()->IndexOf(prt);
+
+ return stack->IsSecondaryFromMaterial(index);
+}
+
+Bool_t AliAnalysisTaskB2::IsFromWeakDecay(const TParticle* prt) const
+{
+//
+// Check if the particle comes from a weak decay
+//
+ AliStack* stack = fMCevent->Stack();
+ Int_t index = stack->Particles()->IndexOf(prt);
+
+ return stack->IsSecondaryFromWeakDecay(index);
+}
+
+Double_t AliAnalysisTaskB2::GetSign(TParticle* prt) const
+{
+//
+// Sign of the particle
+//
+ TParticlePDG* pdg = prt->GetPDG();
+
+ if(pdg != 0) return pdg->Charge();
+
+ return 0;
+}
+
+Double_t AliAnalysisTaskB2::GetPhi(const AliESDtrack* trk) const
+{
+//
+// Azimuthal angle [0,2pi) using the pt at the DCA point
+//
+ Double_t px = trk->Px();
+ Double_t py = trk->Py();
+
+ return TMath::Pi()+TMath::ATan2(-py, -px);
+}
+
+Double_t AliAnalysisTaskB2::GetTheta(const AliESDtrack* trk) const
+{
+//
+// Polar angle using the pt at the DCA point
+//
+ Double_t p = trk->GetP();
+ Double_t pz = trk->Pz();
+
+ 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 = trk->GetP();
+ Double_t e = TMath::Sqrt(p*p + m*m);
+ Double_t 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
+{
+//
+// Momentum for ITS pid
+//
+ Double_t pDCA = trk->GetP();
+ Double_t pTPC = trk->GetTPCmomentum();
+
+ return (pDCA+pTPC)/2.;
+}
+
+Double_t AliAnalysisTaskB2::GetTOFmomentum(const AliESDtrack* trk) const
+{
+//
+// Momentum for TOF pid
+//
+ Double_t pDCA = trk->GetP();
+
+ const AliExternalTrackParam* param = trk->GetOuterParam();
+ Double_t pOut = param ? param->GetP() : trk->GetP();
+
+ return (pDCA+pOut)/2.;
+}
+
+Double_t AliAnalysisTaskB2::GetBeta(const AliESDtrack* trk) const
+{
+//
+// Velocity
+//
+ Double_t t = this->GetTimeOfFlight(trk); // ps
+ Double_t l = trk->GetIntegratedLength(); // cm
+
+ if(t <= 0) return 1.e6; // 1M times the speed of light ;)
+
+ return (l/t)/2.99792458e-2;
+}
+
+Double_t AliAnalysisTaskB2::GetMassSquare(const AliESDtrack* trk) const
+{
+//
+// Square mass
+//
+ Double_t p = this->GetTOFmomentum(trk);
+ Double_t beta = this->GetBeta(trk);
+
+ return p*p*(1./(beta*beta) - 1.);
+}
+
+Int_t AliAnalysisTaskB2::GetChargedMultiplicity(double etaMax) const
+{
+//
+// Charged particle multiplicity using ALICE physical primary definition
+//
+ AliStack* stack = fMCevent->Stack();
+
+ int nch = 0;
+ //for (Int_t i=0; i < stack->GetNprimary(); ++i)
+ for (Int_t i=0; i < stack->GetNtrack(); ++i)
+ {
+ if(!stack->IsPhysicalPrimary(i)) continue;
+
+ TParticle* iParticle = stack->Particle(i);
+
+ if(TMath::Abs(iParticle->Eta()) >= etaMax) continue;
+ //if (iParticle->Pt() < -1) continue;
+
+ TParticlePDG* iPDG = iParticle->GetPDG(); // There are some particles with no PDG
+ if (iPDG && iPDG->Charge() == 0) continue;
+
+ ++nch;
+ }
+
+ return nch;
+}
+
+Double_t AliAnalysisTaskB2::GetTimeOfFlight(const AliESDtrack* trk) const
+{
+//
+// Time of flight associated to the track.
+// Adapted from ANALYSIS/AliAnalysisTaskESDfilter.cxx
+//
+ if(!fESDevent->GetTOFHeader())
+ { //protection in case the pass2 LHC10b,c,d have been processed without tender.
+ Float_t t0spread[10];
+ Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!!
+ for (Int_t j=0; j<10; j++) t0spread[j] = (TMath::Sqrt(fESDevent->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
+ fESDpid->GetTOFResponse().SetT0resolution(t0spread);
+ fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
+
+ fESDpid->SetTOFResponse(fESDevent, (AliESDpid::EStartTimeType_t)fTimeZeroType);
+ }
+
+ if(fESDevent->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(fESDevent, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender.
+
+ Double_t timeZero = fESDpid->GetTOFResponse().GetStartTime(trk->P());
+
+ return trk->GetTOFsignal()-timeZero;
+}
+
+Int_t AliAnalysisTaskB2::GetITSnClusters(const AliESDtrack* trk) const
+{
+//
+// ITS number of clusters
+//
+ UChar_t map = trk->GetITSClusterMap();
+
+ Int_t npoints=0;
+ for(Int_t j=0; j<6; j++)
+ {
+ if(map&(1<<j)) ++npoints;
+ }
+
+ return npoints;
+}
+
+Double_t AliAnalysisTaskB2::GetITSchi2PerCluster(const AliESDtrack* trk) const
+{
+//
+// ITS chi2 per number of clusters
+//
+ Double_t chi2 = trk->GetITSchi2();
+ Int_t ncls = this->GetITSnClusters(trk);
+
+ Double_t chi2ncls = (ncls==0) ? 1.e10 : chi2/ncls;
+
+ return chi2ncls;
+}
+
+Int_t AliAnalysisTaskB2::GetITSnPointsPID(const AliESDtrack* trk) const
+{
+//
+// ITS number of points for PID
+//
+ UChar_t map = trk->GetITSClusterMap();
+
+ Int_t npoints = 0;
+ for(Int_t j=2; j<6; j++)
+ {
+ if(map&(1<<j)) ++npoints;
+ }
+
+ return npoints;
+}
+
+Int_t AliAnalysisTaskB2::GetPidCode(const TString& species) const
+{
+//
+// Return AliPID code of the given species
+//
+ TString name = species;
+ name.ToLower();
+
+ if(name == "electron") return AliPID::kElectron;
+ if(name == "muon") return AliPID::kMuon;
+ if(name == "pion") return AliPID::kPion;
+ if(name == "kaon") return AliPID::kKaon;
+ if(name == "proton") return AliPID::kProton;
+ if(name == "deuteron") return AliPID::kDeuteron;
+ if(name == "triton") return AliPID::kTriton;
+ if(name == "he3") return AliPID::kHe3;
+ if(name == "alpha") return AliPID::kAlpha;
+
+ return -1;
+}
--- /dev/null
+#ifndef ALIANALYSISTASKB2_H
+#define ALIANALYSISTASKB2_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Analysis task for B2
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <AliAnalysisTask.h>
+#include <AliPID.h>
+
+class AliESDtrack;
+class AliMCEvent;
+class AliESDEvent;
+class TString;
+class AliESDtrackCuts;
+class AliLnHistoMap;
+class AliLnID;
+class TParticle;
+
+class AliAnalysisTaskB2: public AliAnalysisTask
+{
+ public:
+ AliAnalysisTaskB2();
+ AliAnalysisTaskB2(const char* name);
+
+ virtual ~AliAnalysisTaskB2();
+
+ virtual void ConnectInputData(Option_t *);
+ virtual void CreateOutputObjects();
+ virtual void Exec(Option_t* option);
+ virtual void Terminate(Option_t *);
+
+ void SetV0ANDtrigger(Bool_t flag=1) { fV0AND = flag; }
+ void SetNoFastOnlyTrigger(Bool_t flag=1) { fNoFastOnly = flag; }
+ void SetKNOmultInterval(Double_t min, Double_t max) { fMinKNOmult = min; fMaxKNOmult = max; }
+ void SetCentralityInterval(Double_t min, Double_t max) { fMinCentrality = min; fMaxCentrality = max; }
+
+ void SetMeanNtrk(Double_t mean) { fMeanNtrk = mean; }
+
+ void SetVertexXInterval(Double_t min, Double_t max) { fMinVx = min; fMaxVx = max; }
+ void SetVertexYInterval(Double_t min, Double_t max) { fMinVy = min; fMaxVy = max; }
+ void SetVertexZInterval(Double_t min, Double_t max) { fMinVz = min; fMaxVz = max; }
+
+ void SetDCAxyInterval(Double_t min, Double_t max) { fMinDCAxy = min; fMaxDCAxy = max; }
+ void SetDCAzInterval(Double_t min, Double_t max) { fMinDCAz = min; fMaxDCAz = max; }
+ void SetMaxNSigmaToVertex(Double_t max ) { fMaxNSigma = max; }
+
+ void SetEtaInterval(Double_t min, Double_t max) { fMinEta = min; fMaxEta = max; }
+ void SetRapidityInterval(Double_t min, Double_t max) { fMinY = min; fMaxY = max; }
+
+ void SetSimulation(Bool_t flag = kTRUE) { fSimulation = flag; }
+ void SetHeavyIons(Bool_t flag = kTRUE) { fHeavyIons = flag; }
+
+ void SetHistogramMap(AliLnHistoMap* map) { fHistoMap = map; }
+
+ void SetESDtrackCuts(AliESDtrackCuts* esdTrackCuts) {fESDtrackCuts = esdTrackCuts; }
+
+ void SetPID(AliLnID* lnID) { fLnID = lnID; }
+ void SetMaxNSigmaITS(Double_t max) { fMaxNSigmaITS = max; }
+ void SetMaxNSigmaTPC(Double_t max) { fMaxNSigmaTPC = max; }
+ void SetMaxNSigmaTOF(Double_t max) { fMaxNSigmaTOF = max; }
+
+ void SetTOFmatch(Bool_t flag = kTRUE) { fTOFmatch = flag; }
+ Bool_t AcceptTOFtrack(const AliESDtrack* trk) const;
+
+ void SetParticleSpecies(const TString& species);
+
+ private:
+
+ AliAnalysisTaskB2(const AliAnalysisTaskB2& other);
+ AliAnalysisTaskB2& operator=(const AliAnalysisTaskB2& other);
+
+ AliLnHistoMap* CreateHistograms();
+
+ Int_t GetParticles();
+ Int_t GetTracks();
+
+ Bool_t IsV0AND() const;
+ 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 maxEta) const;
+
+ Bool_t IsFakeTrack(const AliESDtrack* trk) const;
+ Bool_t IsPhysicalPrimary(const TParticle* prt) const;
+ Bool_t IsFromMaterial(const TParticle* prt) const;
+ Bool_t IsFromWeakDecay(const TParticle* prt) const;
+
+ Double_t GetSign(TParticle* prt) const;
+
+ 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 GetITSmomentum(const AliESDtrack* trk) const;
+ Double_t GetTOFmomentum(const AliESDtrack* trk) const;
+ Double_t GetBeta(const AliESDtrack* trk) const;
+ Double_t GetMassSquare(const AliESDtrack* trk) const;
+ Double_t GetTimeOfFlight(const AliESDtrack* trk) const;
+ Double_t GetITSchi2PerCluster(const AliESDtrack* trk) const;
+ Int_t GetITSnClusters(const AliESDtrack* trk) const;
+ Int_t GetITSnPointsPID(const AliESDtrack* trk) const;
+
+ Int_t GetPidCode(const TString& species) const;
+
+ private:
+
+ TString fSpecies; // particle species for the analysis
+ Int_t fPartCode; // particle species code
+
+ Bool_t fHeavyIons; // analysis of heavy ions data
+ Bool_t fSimulation; // analysis of MC simulation
+
+ Bool_t fMultTrigger; //! track multiplicity trigger flag
+ Bool_t fCentTrigger; //! centrality trigger flag
+ Bool_t fTriggerFired; //! trigger flag
+ Bool_t fGoodVertex; //! good vertex flag
+ Bool_t fPileUpEvent; //! pile-up flag
+
+ Bool_t fV0AND; // V0AND trigger flag
+ Bool_t fNoFastOnly; // No kFastOnly trigger flag
+ Double_t fMinKNOmult; // minimum KNO track multiplicity scaling
+ Double_t fMaxKNOmult; // maximum KNO track multiplicity scaling
+ Double_t fMinCentrality; // minimum centrality for HI
+ Double_t fMaxCentrality; // maximum centrality for HI
+
+ Double_t fNtrk; //! current track multipicity
+ Double_t fMeanNtrk; // average track multiplicity
+ Double_t fKNOmult; //! KNO track multiplicity scaling
+
+ Double_t fMinVx; // vertex low X value
+ Double_t fMaxVx; // vertex high X value
+ Double_t fMinVy; // vertex low Y value
+ Double_t fMaxVy; // vertex high Y value
+ Double_t fMinVz; // vertex low Z value
+ Double_t fMaxVz; // vertex high Z value
+
+ Double_t fMinDCAxy; // minimum DCAxy
+ Double_t fMaxDCAxy; // maximum DCAxy
+ Double_t fMinDCAz; // minimum DCAz
+ Double_t fMaxDCAz; // maximum DCAz
+ Double_t fMaxNSigma; // maximum number of sigmas to primary vertex
+
+ Double_t fMinEta; // minimum pseudorapidity
+ Double_t fMaxEta; // maximum pseudorapidity
+ Double_t fMinY; // minimum rapidity
+ Double_t fMaxY; // maximum rapidity
+
+ AliMCEvent* fMCevent; //! monte carlo event
+ AliESDEvent* fESDevent; //! ESD event
+
+ AliLnHistoMap* fHistoMap; // histogram map (defined somewhere else)
+ AliESDtrackCuts* fESDtrackCuts; // 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
+ Double_t fMaxNSigmaTPC; // maximum number of sigmas to dEdx in the TPC
+ Double_t fMaxNSigmaTOF; // maximum number of sigmas to dEdx in the TOF
+
+ class AliTriggerAnalysis* fTrigAna; //! to access trigger information
+ class AliESDpid* fESDpid; //! ESD pid
+ Bool_t fIsPidOwner; // whether we own fESDpid
+ Int_t fTimeZeroType; // time zero type
+
+ Bool_t fTOFmatch; // TOF match flag
+
+ ClassDef(AliAnalysisTaskB2, 1)
+};
+
+#endif // ALIANALYSISTASKB2_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// class for handling histograms
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TString.h>
+#include <TMap.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TObjString.h>
+#include "AliLnHistoMap.h"
+
+ClassImp(AliLnHistoMap)
+
+AliLnHistoMap::AliLnHistoMap()
+: TObject()
+, fHistoMap(0)
+{
+//
+// Default constructor
+//
+ fHistoMap = new TMap();
+ fHistoMap->SetOwner(kTRUE);
+}
+
+AliLnHistoMap::~AliLnHistoMap()
+{
+//
+// Default destructor
+//
+ delete fHistoMap;
+
+}
+
+Int_t AliLnHistoMap::Write(const char *name, Int_t option, Int_t bsize) const
+{
+//
+// write only the histograms
+//
+ Int_t nbytes = 0;
+
+ TObjString* key;
+ TIter hIter(fHistoMap);
+ while( (key = (TObjString*)hIter.Next()) )
+ {
+ if(fHistoMap->GetValue(key)->InheritsFrom("TH1"))
+ {
+ nbytes += ((TH1*)fHistoMap->GetValue(key))->Write(name,option,bsize);
+ }
+ }
+
+ return nbytes;
+}
+
+Int_t AliLnHistoMap::Write(const char *name, Int_t option, Int_t bsize)
+{
+// write only the histograms
+//
+ return ((const AliLnHistoMap*)this)->Write(name,option,bsize);
+}
+
+TObject* AliLnHistoMap::Add(const TString& keyname, TObject* value)
+{
+//
+// Add a TObject to the histogram map
+//
+ if(fHistoMap->Contains(keyname.Data()))
+ {
+ std::cerr << "WARNING: object " << keyname << " already exists" << std::endl;
+ return 0;
+ }
+
+ TObjString* key = new TObjString(keyname.Data());
+
+ fHistoMap->Add((TObject*)key, value);
+ return value;
+}
+
+TH1D* AliLnHistoMap::Add(const TString& name, Int_t nbins, Double_t xmin, Double_t xmax, const TString& title, const TString& xlabel, const TString& ylabel)
+{
+//
+// Add a TH1D histogram to the histogram map
+//
+ TH1D* value = 0;
+ if(fHistoMap->Contains(name.Data()))
+ {
+ std::cerr << "WARNING: histogram " << name << " already exists" << std::endl;
+ }
+ else
+ {
+ TObjString* key = new TObjString(name.Data());
+ value = new TH1D(name.Data(),title.Data(),nbins,xmin,xmax);
+ value->SetXTitle(xlabel.Data());
+ value->SetYTitle(ylabel.Data());
+
+ fHistoMap->Add((TObject*)key, (TObject*)value);
+ }
+
+ return value;
+}
+
+TH2D* AliLnHistoMap::Add(const TString& name, Int_t xbins, Double_t xmin, Double_t xmax, Int_t ybins, Double_t ymin, Double_t ymax, const TString& title, const TString& xlabel, const TString& ylabel)
+{
+//
+// Add a TH2D histogram to the output map
+//
+ TH2D* value = 0;
+ if(fHistoMap->Contains(name.Data()))
+ {
+ std::cerr << "WARNING: histogram " << name << " already exists" << std::endl;
+ }
+ else
+ {
+ TObjString* key = new TObjString(name.Data());
+
+ value = new TH2D(name.Data(),title.Data(),xbins,xmin,xmax,ybins,ymin,ymax);
+ value->SetXTitle(xlabel.Data());
+ value->SetYTitle(ylabel.Data());
+
+ fHistoMap->Add((TObject*)key, (TObject*)value);
+ }
+
+ return value;
+}
--- /dev/null
+#ifndef ALILNHISTOMAP_H
+#define ALILNHISTOMAP_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// class for handling histograms
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TMap.h>
+
+class TString;
+class TH1D;
+class TH2D;
+
+class AliLnHistoMap: public TObject
+{
+ public:
+ AliLnHistoMap();
+
+ virtual ~AliLnHistoMap();
+
+ virtual Int_t Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0);
+ virtual Int_t Write(const char* name = 0, Int_t option = 0, Int_t bufsize = 0) const;
+
+ TObject* Get(const TString& keyname) const { return fHistoMap->GetValue(keyname.Data()); }
+ TObject* Get(const char* keyname) const { return fHistoMap->GetValue(keyname); }
+
+ TObject* Add(const TString& keyname, TObject* value);
+
+ TH1D* Add(const TString& name, Int_t nbins, Double_t xmin, Double_t xmax, const TString& title="", const TString& xlabel="", const TString& ylabel="");
+
+ TH2D* Add(const TString& name, Int_t xbins, Double_t xmin, Double_t xmax, Int_t ybins, Double_t ymin, Double_t ymax, const TString& title="", const TString& xlabel="", const TString& ylabel="");
+
+ private:
+
+ AliLnHistoMap(const AliLnHistoMap& other);
+ AliLnHistoMap& operator=(const AliLnHistoMap& other);
+
+ private:
+
+ TMap* fHistoMap; //-> histogram map
+
+ ClassDef(AliLnHistoMap, 1)
+};
+
+#endif // ALILNHISTOMAP_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// class for light nuclei identification
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <AliExternalTrackParam.h>
+#include <TParticle.h>
+#include <AliESDtrack.h>
+#include <AliPID.h>
+#include <AliITSPIDResponse.h>
+#include <AliTPCPIDResponse.h>
+#include <TPDGCode.h>
+#include "AliLnID.h"
+
+ClassImp(AliLnID)
+
+AliLnID::AliLnID()
+: TObject()
+, fPidProcedure(kBayes)
+, fkNumSpecies(AliLnID::kSPECIES)
+, fRange(5.)
+, fITSpid(0)
+, fTPCpid(0)
+{
+//
+// Default constructor
+//
+ fSpecies[0] = AliPID::kElectron;
+ fSpecies[1] = AliPID::kMuon;
+ fSpecies[2] = AliPID::kPion;
+ fSpecies[3] = AliPID::kKaon;
+ fSpecies[4] = AliPID::kProton;
+ fSpecies[5] = AliPID::kDeuteron;
+ fSpecies[6] = AliPID::kTriton;
+ fSpecies[7] = AliPID::kHe3;
+ fSpecies[8] = AliPID::kAlpha;
+
+ // equal prior probabilities for all particle species
+ for(Int_t i=0; i<fkNumSpecies; ++i) fPrior[i]=1./fkNumSpecies;
+
+ // ALEPH Bethe Bloch parameters for ITS
+ Double_t param[] = { 0.13, 15.77/0.95, 4.95, 0.312, 2.14, 0.82};
+ fITSpid = new AliITSPIDResponse(¶m[0]);
+
+ fTPCpid = new AliTPCPIDResponse();
+ fTPCpid->SetSigma(0.06,0);
+}
+
+AliLnID::AliLnID(const AliLnID& other)
+: TObject(other)
+, fPidProcedure(other.fPidProcedure)
+, fkNumSpecies(other.fkNumSpecies)
+, fRange(other.fRange)
+, fITSpid(0)
+, fTPCpid(0)
+{
+//
+// Copy constructor
+//
+ for(Int_t i=0; i < fkNumSpecies; ++i)
+ {
+ fSpecies[i] = other.fSpecies[i];
+ fPrior[i] = other.fPrior[i];
+ }
+
+ fITSpid = new AliITSPIDResponse(*(other.fITSpid));
+ fTPCpid = new AliTPCPIDResponse(*(other.fTPCpid));
+}
+
+AliLnID& AliLnID::operator=(const AliLnID& other)
+{
+//
+// Assignment operator
+//
+ if(&other == this) return *this; // check for self-assignment
+
+ fPidProcedure = other.fPidProcedure;
+
+ // deallocate memory
+ delete fITSpid;
+ delete fTPCpid;
+
+ // copy
+ TObject::operator=(other);
+
+ for(Int_t i=0; i < fkNumSpecies; ++i)
+ {
+ fSpecies[i] = other.fSpecies[i];
+ fPrior[i] = other.fPrior[i];
+ }
+
+ fITSpid = new AliITSPIDResponse(*(other.fITSpid));
+ fTPCpid = new AliTPCPIDResponse(*(other.fTPCpid));
+
+ return *this;
+}
+
+void AliLnID::SetPriorProbabilities(Double_t e, Double_t mu, Double_t pi, Double_t k, Double_t p, Double_t d, Double_t t, Double_t he3, Double_t alpha)
+{
+//
+// Set prior probabilities
+//
+ fPrior[0] = e; // electron
+ fPrior[1] = mu; // muon
+ fPrior[2] = pi; // pion
+ fPrior[3] = k; // kaon
+ fPrior[4] = p; // proton
+ fPrior[5] = d; // deuteron
+ fPrior[6] = t; // triton
+ fPrior[7] = he3; // he3
+ fPrior[8] = alpha; // alpha
+}
+
+void AliLnID::SetPriorProbabilities(const Double_t* prob)
+{
+//
+// Set prior probabilities
+//
+ for(Int_t i=0; i < fkNumSpecies; ++i) fPrior[i] = prob[i];
+}
+
+void AliLnID::SetITSBetheBlochParams(const Double_t* par, Double_t res)
+{
+//
+// Set ALEPH Bethe Bloch parameters for ITS
+//
+ delete fITSpid;
+ Double_t param[] = { res, par[0], par[1], par[2], par[3], par[4]};
+ fITSpid = new AliITSPIDResponse(¶m[0]);
+}
+
+void AliLnID::SetTPCBetheBlochParams(const Double_t* par, Double_t mip, Double_t res)
+{
+//
+// Set ALEPH Bethe Bloch parameters for TPC
+//
+ fTPCpid->SetMip(mip);
+ fTPCpid->SetSigma(res,0);
+ fTPCpid->SetBetheBlochParameters(par[0], par[1], par[2], par[3], par[4]);
+}
+
+void AliLnID::SetTPCBetheBlochParams(Double_t c0, Double_t c1, Double_t c2, Double_t c3, Double_t c4, Double_t mip, Double_t res)
+{
+//
+// Set ALEPH Bethe Bloch parameters for TPC
+//
+ fTPCpid->SetMip(mip);
+ fTPCpid->SetSigma(res,0);
+ fTPCpid->SetBetheBlochParameters(c0, c1, c2, c3, c4);
+}
+
+AliLnID::~AliLnID()
+{
+//
+// Default destructor
+//
+ delete fITSpid;
+ delete fTPCpid;
+}
+
+Int_t AliLnID::GetPID(const TParticle* p) const
+{
+//
+// Montecarlo PID
+//
+ enum { kDeuteron=1000010020, kTriton=1000010030, kHelium3=1000020030, kAlpha=1000020040 };
+
+ switch(TMath::Abs(p->GetPdgCode()))
+ {
+ case kElectron: return AliPID::kElectron;
+ case kMuonMinus: return AliPID::kMuon;
+ case kPiPlus: return AliPID::kPion;
+ case kKPlus: return AliPID::kKaon;
+ case kProton: return AliPID::kProton;
+ case kDeuteron: return AliPID::kDeuteron;
+ case kTriton: return AliPID::kTriton;
+ case kHelium3: return AliPID::kHe3;
+ case kAlpha: return AliPID::kAlpha;
+ }
+
+ return -1;
+}
+
+Int_t AliLnID::GetPID(Int_t pidCode, 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, Double_t nSigTPC, Double_t nSigTOF) const
+{
+//
+// PID according to the seletected procedure
+//
+ if(fPidProcedure == kBayes)
+ {
+ return this->GetBayesPID(pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta);
+ }
+ else if(fPidProcedure == kMaxLikelihood)
+ {
+ return this->GetMaxLikelihoodPID(pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta);
+ }
+ else if(fPidProcedure == kTPC)
+ {
+ return this->GetTPCpid(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigTPC);
+ }
+ else if(fPidProcedure == kITSTPC)
+ {
+ return this->GetITSTPCpid(pidCode, pITS, dEdxITS, nPointsITS, nSigITS, pTPC, dEdxTPC, nPointsTPC, nSigTPC);
+ }
+ else if(fPidProcedure == kTPCTOF)
+ {
+ return this->GetTPCTOFpid(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigTPC, pTOF, beta, nSigTOF);
+ }
+
+ return -1;
+}
+
+Int_t AliLnID::GetTPCpid(Int_t pidCode, Double_t pTPC, Double_t dEdxTPC, Double_t nPointsTPC, Double_t nSigmaTPC) const
+{
+//
+// Check if particle with the given pid code is within
+// +/- nSigma around the expected dEdx in the TPC
+//
+ if( this->GetTPCmatch(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigmaTPC)) return pidCode;
+
+ return -1;
+}
+
+Int_t AliLnID::GetTPCTOFpid(Int_t pidCode, Double_t pTPC, Double_t dEdxTPC, Double_t nPointsTPC, Double_t nSigmaTPC, Double_t pTOF, Double_t beta, Double_t nSigmaTOF) const
+{
+//
+// Check if particle with the given pid code is within
+// +/- nSigmaTPC around the expected dEdx in the TPC
+// AND +/- nSigmaTOF around the expected beta in the TOF (when available)
+//
+ if( beta > 0 )
+ {
+ if( this->GetTPCmatch(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigmaTPC)
+ && this->GetTOFmatch(pidCode, pTOF, beta, nSigmaTOF))
+ {
+ return pidCode;
+ }
+ }
+ else
+ {
+ if( this->GetTPCmatch(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigmaTPC))
+ {
+ return pidCode;
+ }
+ }
+
+ return -1;
+}
+
+Int_t AliLnID::GetITSTPCpid(Int_t pidCode, 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
+{
+//
+// Check if particle with the given pid code is within
+// +/- nSigmaITS around the expected dEdx in the ITS
+// AND +/- nSigmaTPC around the expected dEdx in the TPC
+//
+ if( this->GetITSmatch(pidCode, pITS, dEdxITS, nPointsITS, nSigmaITS)
+ && this->GetTPCmatch(pidCode, pTPC, dEdxTPC, nPointsTPC, nSigmaTPC))
+ {
+ return pidCode;
+ }
+
+ return -1;
+}
+
+Int_t AliLnID::GetIndexOfMaxValue(const Double_t* w) const
+{
+//
+// Index with maximum value in the array of size fkNumSpecies
+//
+ Double_t wmax= 0;
+ Int_t imax = -1;
+ for(Int_t i=0; i < fkNumSpecies; ++i)
+ {
+ if(w[i] > wmax)
+ {
+ wmax = w[i];
+ imax = i;
+ }
+ }
+
+ return imax;
+}
+
+Bool_t AliLnID::GetLikelihood(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* r) const
+{
+//
+// Fill r array with the combined likelihood for ITS, TPC and TOF when possible
+// return 1 if success
+//
+ Double_t its[fkNumSpecies];
+ Double_t tpc[fkNumSpecies];
+ Double_t tof[fkNumSpecies];
+
+ Bool_t itsPid = this->GetITSlikelihood(pITS, dEdxITS, nPointsITS, its);
+ Bool_t tpcPid = this->GetTPClikelihood(pTPC, dEdxTPC, nPointsTPC, tpc);
+ Bool_t tofPid = this->GetTOFlikelihood(pTOF, beta, tof);
+
+ if(!itsPid && !tpcPid && !tofPid) return kFALSE;
+
+ // Combine detector responses
+
+ for(Int_t i=0; i<fkNumSpecies; ++i) r[i]=1.;
+
+ if(itsPid)
+ {
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] *= its[i];
+ }
+ if(tpcPid)
+ {
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] *= tpc[i];
+ }
+ if(tofPid)
+ {
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] *= tof[i];
+ }
+
+ return kTRUE;
+}
+
+Int_t AliLnID::GetMaxLikelihoodPID(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta) const
+{
+//
+// Maximum likelihood principle
+//
+ Double_t r[fkNumSpecies];
+
+ if(!this->GetLikelihood(pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, r)) return -1;
+
+ Int_t imax = this->GetIndexOfMaxValue(r);
+
+ if(imax < 0) return -1;
+
+ return fSpecies[imax];
+}
+
+Int_t AliLnID::GetBayesPID(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta) const
+{
+//
+// Bayesian inference
+//
+ Double_t r[fkNumSpecies];
+
+ if(!this->GetLikelihood(pITS, dEdxITS, nPointsITS, pTPC, dEdxTPC, nPointsTPC, pTOF, beta, r)) return -1;
+
+ // Bayes' rule
+
+ Double_t w[fkNumSpecies]; // NOTE: no need to normalize
+ for(Int_t i = 0; i < fkNumSpecies; ++i) w[i] = r[i]*fPrior[i];
+
+ Int_t imax = this->GetIndexOfMaxValue(w);
+
+ if(imax < 0) return -1;
+
+ return fSpecies[imax];
+}
+
+Bool_t AliLnID::GetITSlikelihood(Double_t pITS, Double_t dEdx, Int_t nPointsITS, Double_t* r) const
+{
+//
+// Probability of dEdx in the ITS for each particle species.
+// Adapted from STEER/ESD/AliESDpid.cxx
+// (truncated mean method)
+//
+ if( pITS <= 0) return kFALSE;
+ if( dEdx < 1.) return kFALSE;
+ if( nPointsITS < 3) return kFALSE;
+
+ Bool_t mismatch = kTRUE;
+
+ Double_t sum = 0;
+ for (Int_t i=0; i<fkNumSpecies; ++i)
+ {
+ Double_t mass = AliPID::ParticleMass(fSpecies[i]);
+ Double_t p = fSpecies[i] > AliPID::kTriton ? 2.*pITS : pITS; // correct by Z
+ Double_t expDedx = fITSpid->BetheAleph(p, mass);
+ Double_t sigma = fITSpid->GetResolution(expDedx, nPointsITS);
+
+ //correct by Z^2
+ if(fSpecies[i] > AliPID::kTriton)
+ {
+ expDedx *= 4.;
+ sigma *= 4.;
+ }
+
+ if (TMath::Abs(dEdx - expDedx) > fRange*sigma)
+ {
+ r[i]=TMath::Exp(-0.5*fRange*fRange)/sigma;
+ }
+ else
+ {
+ r[i]=TMath::Exp(-0.5*(dEdx-expDedx)*(dEdx-expDedx)/(sigma*sigma))/sigma;
+ mismatch = kFALSE;
+ }
+
+ sum += r[i];
+ }
+
+ if(sum <= 0. || mismatch)
+ {
+ return kFALSE;
+ }
+
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] /= sum;
+
+ return kTRUE;
+}
+
+Bool_t AliLnID::GetTPClikelihood(Double_t pTPC, Double_t dEdx, Int_t nPointsTPC, Double_t* r) const
+{
+//
+// Probability of dEdx for each particle species in the TPC.
+// Adapted from STEER/ESD/AliESDpid.cxx
+//
+ if( pTPC <= 0) return kFALSE;
+ if( dEdx < 1.) return kFALSE;
+
+ Bool_t mismatch = kTRUE;
+
+ Double_t sum = 0;
+ for(Int_t i=0; i < fkNumSpecies; ++i)
+ {
+ Double_t p = fSpecies[i] > AliPID::kTriton ? 2.*pTPC : pTPC; // correct by Z
+ Double_t expDedx = fTPCpid->GetExpectedSignal(p, fSpecies[i]);
+ Double_t sigma = fTPCpid->GetExpectedSigma(p, nPointsTPC, fSpecies[i]);
+
+ if(fSpecies[i] > AliPID::kTriton) // correct by Z^2
+ {
+ expDedx *= 4.;
+ sigma *= 4.;
+ }
+
+ if(TMath::Abs(dEdx - expDedx) > fRange*sigma)
+ {
+ r[i] = TMath::Exp(-0.5*fRange*fRange)/sigma;
+ }
+ else
+ {
+ r[i] = TMath::Exp(-0.5*(dEdx-expDedx)*(dEdx-expDedx)/(sigma*sigma))/sigma;
+ mismatch=kFALSE;
+ }
+
+ sum += r[i];
+ }
+
+ if(sum <= 0. || mismatch)
+ {
+ return kFALSE;
+ }
+
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] /= sum;
+
+ return kTRUE;
+}
+
+Bool_t AliLnID::GetTOFlikelihood(Double_t pTOF, Double_t beta, Double_t* r) const
+{
+//
+// Probability of beta for each particle species
+//
+ if( pTOF <= 0) return kFALSE;
+ if( beta <= 0) return kFALSE;
+
+ Double_t mismatch = kTRUE;
+
+ Double_t sum = 0;
+ for(Int_t i=0; i < fkNumSpecies; ++i)
+ {
+ Double_t mass = AliPID::ParticleMass(fSpecies[i]);
+ Double_t p = fSpecies[i] > AliPID::kTriton ? 2.*pTOF : pTOF; // correct by Z
+ Double_t expBeta = this->Beta(p,mass);
+ Double_t sigma = this->GetBetaExpectedSigma(p,mass);
+
+ if(TMath::Abs(beta-expBeta) > fRange*sigma)
+ {
+ r[i] = TMath::Exp(-0.5*fRange*fRange)/sigma;
+ }
+ else
+ {
+ r[i] = TMath::Exp(-0.5*(beta-expBeta)*(beta-expBeta)/(sigma*sigma))/sigma;
+ mismatch = kFALSE;
+ }
+
+ sum += r[i];
+ }
+
+ if(sum <= 0. || mismatch)
+ {
+ return kFALSE;
+ }
+
+ for(Int_t i=0; i < fkNumSpecies; ++i) r[i] /= sum;
+
+ return kTRUE;
+}
+
+Double_t AliLnID::Beta(Double_t p, Double_t m) const
+{
+//
+// Expected beta for mass hypothesis m
+//
+ return p/TMath::Sqrt(p*p+m*m);
+}
+
+Double_t AliLnID::GetBetaExpectedSigma(Double_t p, Double_t m) const
+{
+//
+// Expected sigma for the given mass hypothesis
+//
+ const Double_t kC0 = 0.0131203;
+ const Double_t kC1 = 0.00670148;
+
+ Double_t kappa = kC0*m*m/(p*(p*p+m*m)) + kC1;
+ return kappa*Beta(p,m);
+}
+
+Double_t AliLnID::GetM2ExpectedSigma(Double_t p, Double_t m2) const
+{
+//
+// Expected deuteron M2 width for the given momentum
+//
+ Double_t c0 = 6.76363e-03;
+ Double_t c1 = 7.04984e-02;
+
+ return c0/(p*p) + c1*(p*p/m2 + 1.);
+}
+
+Bool_t AliLnID::GetITSmatch(Int_t pid, Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t nSigma) const
+{
+//
+// Check if the signal is less than nSigma from the expected ITS dEdx
+//
+ Double_t mass = AliPID::ParticleMass(pid);
+ Double_t p = pid > AliPID::kTriton ? 2.*pITS : pITS; // correct by Z
+ Double_t expDedx = fITSpid->BetheAleph(p, mass);
+ Double_t sigma = fITSpid->GetResolution(expDedx, nPointsITS);
+
+ if(pid > AliPID::kTriton) //correct by Z^2
+ {
+ expDedx *= 4.;
+ sigma *= 4.;
+ }
+
+ if(TMath::Abs(dEdxITS-expDedx) < nSigma*sigma)
+ {
+ return kTRUE;
+ }
+
+ return kFALSE;
+}
+
+Bool_t AliLnID::GetTPCmatch(Int_t pid, Double_t pTPC, Double_t dEdxTPC, Double_t nPointsTPC, Double_t nSigma) const
+{
+//
+// Check if the signal is less than nSigma from the expected TPC dEdx
+//
+ Double_t p = pid > AliPID::kTriton ? 2.*pTPC : pTPC; // correct by Z
+ Double_t expDedx = fTPCpid->GetExpectedSignal(p, (AliPID::EParticleType)pid);
+ Double_t sigma = fTPCpid->GetExpectedSigma(p, nPointsTPC, (AliPID::EParticleType)pid);
+
+ if(pid > AliPID::kTriton) // correct by Z^2
+ {
+ expDedx *= 4.;
+ sigma *= 4.;
+ }
+
+ if(TMath::Abs(dEdxTPC-expDedx) < nSigma*sigma)
+ {
+ return kTRUE;
+ }
+
+ return kFALSE;
+}
+
+Bool_t AliLnID::GetTOFmatch(Int_t pid, Double_t pTOF, Double_t beta, Double_t nSigma) const
+{
+//
+// Check if the signal is less than nSigma from the expected velocity
+//
+ Double_t mass = AliPID::ParticleMass(pid);
+ Double_t p = pid > AliPID::kTriton ? 2.*pTOF : pTOF;
+ Double_t expBeta = this->Beta(p, mass);
+ Double_t sigma = this->GetBetaExpectedSigma(p, mass);
+
+ if(TMath::Abs(beta-expBeta) < nSigma*sigma)
+ {
+ return kTRUE;
+ }
+
+ return kFALSE;
+}
+
+Bool_t AliLnID::GetM2match(Int_t pid, Double_t pTOF, Double_t m2, Double_t nSigma) const
+{
+//
+// Check if the signal is less than nSigma from the expected m2
+//
+ Double_t p = pid > AliPID::kTriton ? 2.*pTOF : pTOF; // correct by Z
+ Double_t expM2 = AliPID::ParticleMass(pid); expM2 *= expM2;
+ Double_t sigma = this->GetM2ExpectedSigma(p, expM2);
+
+ if(TMath::Abs(m2-expM2) < nSigma*sigma)
+ {
+ return kTRUE;
+ }
+
+ return kFALSE;
+}
+
+Bool_t AliLnID::IsITSTPCmismatch(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t nSigma) const
+{
+//
+// Check track TPC mismatch with ITS
+//
+ for(Int_t i=0; i<fkNumSpecies; ++i)
+ {
+ if(this->GetTPCmatch(fSpecies[i], pTPC, dEdxTPC, nPointsTPC, nSigma) &&
+ this->GetITSmatch(fSpecies[i], pITS, dEdxITS, nPointsITS, 5.))
+ {
+ return kFALSE;
+ }
+ }
+
+ return kTRUE;
+}
+
+Bool_t AliLnID::IsITSTOFmismatch(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTOF, Double_t beta, Double_t nSigma) const
+{
+//
+// Check track TOF mismatch with ITS
+//
+ for(Int_t i=0; i<fkNumSpecies; ++i)
+ {
+ if(this->GetTOFmatch(fSpecies[i], pTOF, beta, nSigma) &&
+ this->GetITSmatch(fSpecies[i], pITS, dEdxITS, nPointsITS, 3.))
+ {
+ return kFALSE;
+ }
+ }
+
+ return kTRUE;
+}
+
+Bool_t AliLnID::IsTPCTOFmismatch(Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta, Double_t nSigma) const
+{
+//
+// Check track TOF mismatch with TPC
+//
+ for(Int_t i=0; i<fkNumSpecies; ++i)
+ {
+ if(this->GetTOFmatch(fSpecies[i], pTOF, beta, nSigma) &&
+ this->GetTPCmatch(fSpecies[i], pTPC, dEdxTPC, nPointsTPC, 3.))
+ {
+ return kFALSE;
+ }
+ }
+
+ return kTRUE;
+}
--- /dev/null
+#ifndef ALILNID_H
+#define ALILNID_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// class for light nuclei identification
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <AliPID.h>
+
+class AliLnID: public TObject
+{
+ public:
+ AliLnID();
+ AliLnID(const AliLnID& other);
+ AliLnID& operator=(const AliLnID& other);
+
+ virtual ~AliLnID();
+
+ Int_t GetPID(const class TParticle* p) 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 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;
+
+ Int_t GetTPCTOFpid(Int_t partCode, Double_t pTPC, Double_t dEdx, Double_t nPoints, Double_t nSigmaTPC, Double_t pTOF, Double_t beta, Double_t nSigmaTOF) const;
+
+ Int_t GetMaxLikelihoodPID(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta) const;
+
+ Int_t GetBayesPID(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta) const;
+
+ Bool_t GetITSlikelihood(Double_t p, Double_t dEdx, Int_t nPoints, Double_t* r) const;
+ Bool_t GetTPClikelihood(Double_t p, Double_t dEdx, Int_t nPoints, Double_t* r) const;
+ Bool_t GetTOFlikelihood(Double_t p, Double_t beta, Double_t* r) const;
+
+ Double_t GetBetaExpectedSigma(Double_t p, Double_t mass) const;
+ Double_t GetM2ExpectedSigma(Double_t p, Double_t mass) const;
+
+ Bool_t GetITSmatch(Int_t pid, Double_t p, Double_t dEdx, Int_t nPoints, Double_t nSigma=3.) const;
+ Bool_t GetTPCmatch(Int_t pid, Double_t pTPC, Double_t dEdx, Double_t nPoints, Double_t nSigma=3.) const;
+ Bool_t GetTOFmatch(Int_t pid, Double_t pTOF, Double_t beta, Double_t nSigma=3.) const;
+ Bool_t GetM2match(Int_t pid, Double_t pTOF, Double_t m2, Double_t nSigma=3.) const;
+
+ Int_t GetPidProcedure() const { return fPidProcedure; }
+
+ Bool_t IsITSTPCmismatch(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t nSigma=5.) const;
+ Bool_t IsITSTOFmismatch(Double_t pITS, Double_t dEdxITS, Int_t nPointsITS, Double_t pTOF, Double_t beta, Double_t nSigma=5.) const;
+ Bool_t IsTPCTOFmismatch(Double_t pTPC, Double_t dEdxTPC, Int_t nPointsTPC, Double_t pTOF, Double_t beta, Double_t nSigma=5.) const;
+
+ void SetPriorProbabilities(const Double_t* prob);
+ void SetPriorProbabilities(Double_t e, Double_t mu, Double_t pi, Double_t k, Double_t p, Double_t d, Double_t t, Double_t he3, Double_t alpha);
+
+ void SetPidProcedure(Int_t proc) { fPidProcedure = proc; }
+
+ void SetITSBetheBlochParams(const Double_t* param, Double_t res=0.13);
+
+ void SetTPCBetheBlochParams(const Double_t* param, Double_t mip=1., Double_t res=0.06);
+ void SetTPCBetheBlochParams(Double_t c0, Double_t c1, Double_t c2, Double_t c3, Double_t c4, Double_t mip=1., Double_t res=0.06);
+
+ enum { kSPECIES = 9 };
+ enum { kBayes, kMaxLikelihood, kTPC, kITSTPC, kTPCTOF };
+
+ private:
+
+ Double_t Beta(Double_t p, Double_t m) const;
+ Bool_t GetLikelihood(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* r) const;
+ Int_t GetIndexOfMaxValue(const Double_t* w) const;
+
+ private:
+
+ Int_t fPidProcedure; // PID procedure code
+ const Int_t fkNumSpecies; // number of particle species
+ AliPID::EParticleType fSpecies[kSPECIES]; // particle species known by the pid
+
+ Double_t fPrior[kSPECIES]; // prior probabilities
+ Double_t fRange; // number of sigmas to the expected values
+
+ class AliITSPIDResponse* fITSpid; // ITS likelihood
+ class AliTPCPIDResponse* fTPCpid; // TPC likelihood
+
+ ClassDef(AliLnID, 1)
+};
+
+#endif // ALILNID_H
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// Add task B2
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+Double_t GetMeanNtrk(const TString& period);
+
+AliAnalysisTaskB2* AddTaskB2(const TString& species,
+ const TString& outputfile,
+ const TString& trksel,
+ Int_t pidProc,
+ const TString& periodname,
+ Bool_t simulation=kFALSE,
+ Bool_t heavyIons=kFALSE,
+ Double_t maxDCAxy=1,
+ Double_t maxDCAz=2,
+ Double_t minKNOmult=-10,
+ Double_t maxKNOmult=10000,
+ Bool_t V0AND=kFALSE,
+ Double_t minCentrality=0,
+ Double_t maxCentrality=20)
+{
+//
+// Create, configure and add the analysis task to the analysis manager
+//
+ using namespace std;
+
+ // sample config
+
+ const Double_t kMaxVx = 1.;
+ const Double_t kMaxVy = 1.;
+ const Double_t kMaxVz = 10.;
+
+ const Double_t kMaxY = 0.5;
+ const Double_t kMaxEta = 0.8;
+ const Double_t kMaxNSigma = 3.;
+
+ const Int_t kMaxNSigmaITS = 3;
+ const Int_t kMaxNSigmaTPC = 3;
+ const Int_t kMaxNSigmaTOF = 3;
+ const Int_t kMinTPCnCls = 70;
+
+ TString period = periodname;
+ period.ToLower();
+
+ // Get pointer to the existing analysis manager
+
+ AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ cerr << "AddTaskB2: no analysis manager to connect to" << endl;
+ return 0;
+ }
+
+ if (!mgr->GetInputEventHandler())
+ {
+ cerr << "AddTaskB2: this task requires an input event handler" << endl;
+ return 0;
+ }
+
+ // Create and configure the task
+
+ AliAnalysisTaskB2* task = new AliAnalysisTaskB2(Form("B2_%s",species.Data()));
+
+ task->SetParticleSpecies(species);
+ task->SetSimulation(simulation);
+ task->SetHeavyIons(heavyIons);
+ task->SetV0ANDtrigger(V0AND);
+
+ task->SetMaxNSigmaITS(kMaxNSigmaITS);
+ task->SetMaxNSigmaTPC(kMaxNSigmaTPC);
+ task->SetMaxNSigmaTOF(kMaxNSigmaTOF);
+
+ task->SetMeanNtrk(GetMeanNtrk(period));
+ task->SetKNOmultInterval(minKNOmult, maxKNOmult);
+ task->SetVertexXInterval(-kMaxVx, kMaxVx);
+ task->SetVertexYInterval(-kMaxVy, kMaxVy);
+ task->SetVertexZInterval(-kMaxVz, kMaxVz);
+
+ task->SetCentralityInterval(minCentrality, maxCentrality);
+
+ if(period=="lhc11a_wsdd" || period=="lhc11a_wosdd")
+ {
+ task->SetNoFastOnlyTrigger();
+ }
+
+ // histograms
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/CreateHistograms.C");
+
+ AliLnHistoMap* hMap = CreateHistograms(species, simulation, maxDCAxy, kMaxEta, kMaxY, heavyIons);
+
+ task->SetHistogramMap(hMap);
+
+ // track selection criteria
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/TrackCuts.C");
+
+ AliESDtrackCuts* trkCuts = TrackCuts(task, trksel, maxDCAxy, maxDCAz, kMaxNSigma, kMinTPCnCls, kMaxEta);
+ task->SetESDtrackCuts(trkCuts);
+
+ // PID
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/BetheBlochParams.C");
+
+ Double_t bethe[5];
+ BetheBlochParams(bethe, period);
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/PriorProbabilities.C");
+
+ Double_t prob[9];
+ PriorProbabilities(prob, period, trksel);
+
+ AliLnID* lnID = new AliLnID();
+
+ lnID->SetTPCBetheBlochParams(bethe);
+ lnID->SetPriorProbabilities(prob);
+ lnID->SetPidProcedure(pidProc);
+
+ task->SetPID(lnID);
+
+ // Add task to the manager
+
+ mgr->AddTask(task);
+
+ // input and output containers
+
+ AliAnalysisDataContainer* output = mgr->CreateContainer(outputfile.Data(), AliLnHistoMap::Class(), AliAnalysisManager::kOutputContainer, outputfile.Data());
+
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 0, output);
+
+ return task;
+}
+
+Double_t GetMeanNtrk(const TString& period)
+{
+//
+// average track multiplicity <Ntrk> for the given period
+//
+ if(period =="lhc10c900") return 3.70158; // pass3
+ if(period =="lhc10b_pass2") return 9.86258;
+ if(period =="lhc10c_pass2") return 9.61402;
+ if(period =="lhc10b") return 5.96104; // pass3
+ if(period =="lhc10c") return 5.94719; // pass3
+ if(period =="lhc10d") return 5.82333; // pass2
+ if(period =="lhc10e") return 5.89367; // pass2
+ if(period =="lhc11a_wosdd") return 4.28597; // pass3
+ if(period =="lhc11a_wsdd") return 4.57960; // pass4
+
+ // MC
+ if(period =="lhc10e13") return 3.17763;
+ if(period =="lhc10f6a") return 4.41362;
+ if(period =="lhc10e21") return 4.74991;
+ if(period =="lhc11e3a_plus_wosdd") return 3.37669;
+ if(period =="lhc11e3a_plus_wosdd") return 3.55973;
+
+ return 1;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// macro for Bethe-Bloch parameters
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+void SetParameters(Double_t* param, Double_t c0, Double_t c1, Double_t c2, Double_t c3, Double_t c4)
+{
+//
+// assign values to an array
+//
+ param[0] = c0;
+ param[1] = c1;
+ param[2] = c2;
+ param[3] = c3;
+ param[4] = c4;
+}
+
+void BetheBlochParams(Double_t* param, const TString& periodname)
+{
+//
+// TPC Bethe Bloch ALEPH paramaters for different beam periods
+//
+ TString period = periodname;
+ period.ToLower();
+
+ if(period == "lhc10b_pass2")
+ {
+ SetParameters(param, 3.22422, 10.9345, 1.26309e-05, 2.26343, 2.43587);
+ }
+ else if(period == "lhc10b" || period == "lhc10b_pass3")
+ {
+ SetParameters(param, 5.24531, 5.82813, 0.000431522, 2.47198, 1.38539);
+
+ }
+ else if(period == "lhc10c900") // pass3
+ {
+ SetParameters(param, 1.4857, 22.9345, 1.77678e-11, 2.26456, 4.44306);
+ }
+ else if(period == "lhc10c_pass2")
+ {
+ SetParameters(param, 1.49726, 24.5879, 2.76442e-11, 2.15661, 4.91248);
+ }
+ else if(period == "lhc10c" || period == "lhc10c_pass3")
+ {
+ SetParameters(param, 7.41249, 5.13753, 0.000738319, 2.55495, 1.33519);
+ }
+ else if(period == "lhc10d" || period == "lhc10e" || period == "lhc10d_pass2" || period == "lhc10e_pass2")
+ {
+ SetParameters(param, 1.59526, 24.6438, 3.5082e-11, 2.18753, 3.7487);
+ }
+ else if(period == "lhc10h") // heavy ions
+ {
+ SetParameters(param, 2.77047, 14.6777, 5.62959e-08, 2.30422, 2.35623);
+ }
+ else if(period == "lhc11a" || period == "lhc11a_wsdd" || period == "lhc11a_wosdd") // pp 2.76 TeV
+ {
+ SetParameters(param, 4.94865, 8.29784, 9.95186e-05, 2.22417, 1.51139);
+ }
+
+ // simulation
+
+ else if(period == "lhc10e12" || period == "lhc10e13" || period == "lhc10f6a" || period == "lhc10f6" || period == "lhc10e21")
+ {
+ SetParameters(param, 1.98509, 16.9132, 2.27954e-10, 2.1544, 3.94486);
+ }
+ else
+ {
+ SetParameters(param, 4.4194, 7.50931, 1.34e-05, 2.22085, 1.80461);
+ }
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// macro for creating histograms
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+AliLnHistoMap* CreateHistograms(const TString& species, Bool_t simulation, Double_t maxDCAxy, Double_t maxEta, Double_t maxY, Bool_t heavyIons)
+{
+//
+// Define an create the histograms for the analysis
+//
+ AliLnHistoMap* hMap = new AliLnHistoMap();
+
+ Int_t A = 1;
+ if(species == "Deuteron") A = 2;
+ else if(species == "Triton") A = 3;
+ else if(species == "He3") A = 3;
+ else if(species == "Alpha") A = 4;
+
+ // pt
+ Int_t ptBins = 100;
+ Double_t ptMin = 0.*A;
+ Double_t ptMax = 10.*A;
+
+ // eta and rapidity
+ Int_t etaBins = 300;
+ Double_t etaMin = -1.5;
+ Double_t etaMax = 1.5;
+
+ // phi
+ Int_t phiBins = 180;
+ Double_t phiMin = 0.;
+ Double_t phiMax = 2.*TMath::Pi();
+
+ // DCA distribution bin width = 0.001 cm
+ Int_t dcaxyBins = 2.*maxDCAxy/0.001;
+ Double_t dcaxyMin = -maxDCAxy;
+ Double_t dcaxyMax = maxDCAxy;
+
+ Int_t dcazBins = 1300;
+ Double_t dcazMin = -3.25;
+ Double_t dcazMax = 3.25;
+
+ // m2 bins
+ Int_t m2Bins = 200;
+ Double_t m2Min = 0;
+ Double_t m2Max = 10;
+
+ // track multiplicity
+ Int_t ntrkBins = 200;
+ Double_t ntrkMin = 0;
+ Double_t ntrkMax = 200;
+
+ if(heavyIons)
+ {
+ ntrkBins = 2000;
+ ntrkMin = 0;
+ ntrkMax = 2000;
+ }
+
+ // stats
+
+ TH1D* hStats = new TH1D(species + "_Stats", "Stats", 8, 0, 8);
+
+ hStats->GetXaxis()->SetBinLabel(1,"Events");
+ hStats->GetXaxis()->SetBinLabel(2,"Trig");
+ hStats->GetXaxis()->SetBinLabel(3,"Ana");
+ hStats->GetXaxis()->SetBinLabel(4,"Vtx");
+ hStats->GetXaxis()->SetBinLabel(5,"Vz");
+ hStats->GetXaxis()->SetBinLabel(6,"Vxy");
+ hStats->GetXaxis()->SetBinLabel(7,"VertexerZ");
+ hStats->GetXaxis()->SetBinLabel(8,"NoPileUp");
+
+ hStats->SetStats(0);
+ hMap->Add( species + "_Stats", (TObject*)hStats);
+
+ if(simulation)
+ {
+ // PID table
+
+ TH2D* hPidTable = new TH2D(species + "_Stats_PID_Table", "Particle table (Row: PID, Col: MC)", 10,0,10,10,0,10);
+
+ hPidTable->GetXaxis()->SetBinLabel(1,"e");
+ hPidTable->GetXaxis()->SetBinLabel(2,"#mu");
+ hPidTable->GetXaxis()->SetBinLabel(3,"#pi");
+ hPidTable->GetXaxis()->SetBinLabel(4,"K");
+ hPidTable->GetXaxis()->SetBinLabel(5,"p");
+ hPidTable->GetXaxis()->SetBinLabel(6,"d");
+ hPidTable->GetXaxis()->SetBinLabel(7,"t");
+ hPidTable->GetXaxis()->SetBinLabel(8,"h");
+ hPidTable->GetXaxis()->SetBinLabel(9,"#alpha");
+ hPidTable->GetXaxis()->SetBinLabel(10,"sum");
+
+ hPidTable->GetYaxis()->SetBinLabel(1,"e");
+ hPidTable->GetYaxis()->SetBinLabel(2,"#mu");
+ hPidTable->GetYaxis()->SetBinLabel(3,"#pi");
+ hPidTable->GetYaxis()->SetBinLabel(4,"K");
+ hPidTable->GetYaxis()->SetBinLabel(5,"p");
+ hPidTable->GetYaxis()->SetBinLabel(6,"d");
+ hPidTable->GetYaxis()->SetBinLabel(7,"t");
+ hPidTable->GetYaxis()->SetBinLabel(8,"h");
+ hPidTable->GetYaxis()->SetBinLabel(9,"#alpha");
+ hPidTable->GetYaxis()->SetBinLabel(10,"sum");
+
+ hPidTable->SetStats(0);
+ hMap->Add( species + "_Stats_PID_Table", (TObject*)hPidTable);
+ }
+
+ TH1D* hStatsPid = new TH1D(species + "_Stats_PID", "Stats PID", 9, 0, 9);
+ hStatsPid->GetXaxis()->SetBinLabel(1,"e");
+ hStatsPid->GetXaxis()->SetBinLabel(2,"\\mu");
+ hStatsPid->GetXaxis()->SetBinLabel(3,"\\pi");
+ hStatsPid->GetXaxis()->SetBinLabel(4,"K");
+ hStatsPid->GetXaxis()->SetBinLabel(5,"p");
+ hStatsPid->GetXaxis()->SetBinLabel(6,"d");
+ hStatsPid->GetXaxis()->SetBinLabel(7,"t");
+ hStatsPid->GetXaxis()->SetBinLabel(8,"h");
+ hStatsPid->GetXaxis()->SetBinLabel(9,"\\alpha");
+ hStatsPid->SetStats(0);
+ hMap->Add( species + "_Stats_PID", (TObject*)hStatsPid);
+
+ // detector signals
+
+ hMap->Add( species + "_ITS_dEdx_P", 1000, 0.001, 10., 1500, 0.001, 1500., "", "p/Z (GeV/c)", "dE/dx");
+ hMap->Add( species + "_TPC_dEdx_P", 1000, 0.001, 10., 1500, 0.001, 1500., "", "p/Z (GeV/c)", "dE/dx");
+ hMap->Add( species + "_TOF_Beta_P", 1000, 0.001, 10., 1200, 0.001, 1.2, "", "p/Z (GeV/c)", "#beta");
+ hMap->Add( species + "_TOF_Mass_P", 1000, 0.001, 10., 500, 0.01, 5., "", "p/Z (GeV/c)", "Mass (GeV/c^{2})");
+
+ // TrackCuts
+
+ hMap->Add( species + "_TrackCuts_DCAxy", 200, -4, 4, "Track cuts", "DCA_{xy} (cm)");
+ hMap->Add( species + "_TrackCuts_DCAz", 200, -6, 6, "Track cuts", "DCA_{z} (cm)");
+ hMap->Add( species + "_TrackCuts_NSigma", 200, 0, 10, "Track cuts", "N_{#sigma}");
+
+ hMap->Add( species + "_TrackCuts_ITSchi2PerCls", 200, 0, 40, "Track cuts", "ITS #chi^{2} / Cluster");
+
+ hMap->Add( species + "_TrackCuts_TPCncls", 200, 0, 200, "Track cuts", "TPC clusters");
+ hMap->Add( species + "_TrackCuts_TPCclsOverF", 200, 0, 2, "Track cuts", "TPC clusters/Findable");
+ hMap->Add( species + "_TrackCuts_TPCxRows", 200, 0, 200, "Track cuts", "TPC crossed rows");
+ hMap->Add( species + "_TrackCuts_TPCchi2PerCls", 200, 0, 20, "Track cuts", "TPC #chi^{2} / Cluster");
+ hMap->Add( species + "_TrackCuts_TPCchi2Global", 200, -2, 38, "Track cuts", "#chi^{2} of constrained TPC vs global track");
+
+ hMap->Add( species + "_Before_Phi_Theta", 180, 0, TMath::Pi(), 360, 0, 2.*TMath::Pi(), "Global tracks (before track cuts)", "#theta (rad)", "#phi (rad)");
+
+ hMap->Add( species + "_After_Phi_Theta", 180, 0, TMath::Pi(), 360, 0, 2.*TMath::Pi(), "Global tracks (after track cuts)", "#theta (rad)", "#phi (rad)");
+
+ // multiplicity
+
+ hMap->Add( species + "_Event_Ntrk", ntrkBins, ntrkMin, ntrkMax, "Track multiplicity (all events)", "N_{trk}", "Events");
+ hMap->Add( species + "_Ana_Event_Ntrk", ntrkBins, ntrkMin, ntrkMax, "Track multiplicity", "N_{trk}", "Events");
+
+ if(simulation)
+ {
+ hMap->Add( species + "_Ana_Event_Nch_Ntrk", 200, 0., 200., 200, 0., 200., "Track multiplicity", "N_{trk}", "N_{ch}");
+ }
+
+ hMap->Add( species + "_Vertex_XZ", 500, -50., 50., 800, -1., 1., "Primary Vertex", "Z (cm)", "X (cm)");
+ hMap->Add( species + "_Vertex_YZ", 500, -50., 50., 800, -1., 1., "Primary Vertex", "Z (cm)", "Y (cm)");
+ hMap->Add( species + "_Vertex_YX", 800, -1., 1., 800, -1., 1., "Primary Vertex", "X (cm)", "Y (cm)");
+
+ hMap->Add( species + "_Ana_Vertex_XZ", 500, -50., 50., 800, -1., 1., "Primary Vertex", "Z (cm)", "X (cm)");
+ hMap->Add( species + "_Ana_Vertex_YZ", 500, -50., 50., 800, -1., 1., "Primary Vertex", "Z (cm)", "Y (cm)");
+ hMap->Add( species + "_Ana_Vertex_YX", 800, -1., 1., 800, -1., 1., "Primary Vertex", "X (cm)", "Y (cm)");
+
+ if(heavyIons)
+ {
+ hMap->Add( species + "_V0_Mult", 1000, 0, 100000., "Corrected V0 Multiplicity", "V0M", "Events");
+ hMap->Add( species + "_V0_Scaled_Mult", 1000, 0, 10000., "Corrected Scaled V0 Multiplicity", "V0SM", "Events");
+
+ hMap->Add( species + "_Ana_V0_Mult", 1000, 0, 100000., "Corrected V0 Multiplicity (after centrality sel.)", "V0M", "Events");
+ hMap->Add( species + "_Ana_V0_Scaled_Mult", 1000, 0, 10000., "Corrected Scaled V0 Multiplicity (after centrality sel.)", "V0SM", "Events");
+ }
+
+ // positive and negative
+
+ TString particle[] = { Form("Anti%s",species.Data()), species };
+
+ for(Int_t i=0; i<2; ++i)
+ {
+ // pid with TOF
+
+ hMap->Add( particle[i] + "_PID_M2_Pt", ptBins, ptMin, ptMax, m2Bins, m2Min, m2Max, Form("%s candidate",particle[i].Data()), "p_{T} (GeV/c)", "m^{2} (GeV^{2}/c^{4})");
+
+ // as a function of momentum
+ hMap->Add( particle[i] + "_PID_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("%s candidate",particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_PID_DCAz_Pt", ptBins, ptMin, ptMax, dcazBins, dcazMin, dcazMax, Form("%s candidate",particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{z} (cm)");
+
+ hMap->Add( particle[i] + "_PID_NSigma_Pt", ptBins, ptMin, ptMax, 200, 0., 10., Form("%s candidate",particle[i].Data()), "p_{T} (GeV/c)", "N\\sigma");
+
+ if(simulation)
+ {
+ hMap->Add( particle[i] + "_Sim_PID_M2_Pt", ptBins, ptMin, ptMax, m2Bins, m2Min, m2Max, Form("%s after PID", particle[i].Data()), "p_{T} (GeV/c)", "m^{2} (GeV^{2}/c^{4})");
+
+ hMap->Add( particle[i] + "_Sim_Prim_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Primary %s", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_Fdwn_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Feed-down %s", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_Mat_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("%s from materials", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ // primaries
+ hMap->Add( particle[i] + "_Sim_PID_Prim_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Primary %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Prim_DCAz_Pt", ptBins, ptMin, ptMax, dcazBins, dcazMin, dcazMax, Form("Primary %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{z} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Prim_NSigma_Pt", ptBins, ptMin, ptMax, 200, 0., 10., Form("Primary %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "N\\sigma");
+
+ // feed-down
+ hMap->Add( particle[i] + "_Sim_PID_Fdwn_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Feed-down %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Fdwn_DCAz_Pt", ptBins, ptMin, ptMax, dcazBins, dcazMin, dcazMax, Form("Feed-down %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{z} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Fdwn_NSigma_Pt", ptBins, ptMin, ptMax, 200, 0., 10., Form("Feed-down %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "N\\sigma");
+
+ // secondaries from materials
+ hMap->Add( particle[i] + "_Sim_PID_Mat_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("%s from materials after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Mat_DCAz_Pt", ptBins, ptMin, ptMax, dcazBins, dcazMin, dcazMax, Form("%s from materials after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{z} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Mat_NSigma_Pt", ptBins, ptMin, ptMax, 200, 0., 10., Form("%s from materials after PID", particle[i].Data()), "p_{T} (GeV/c)", "N\\sigma");
+
+ // fake tracks
+ hMap->Add( particle[i] + "_Sim_PID_Fake_Prim_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Fake Primary %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Fake_Fdwn_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Fake feed-down %s after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Fake_Mat_DCAxy_Pt", ptBins, ptMin, ptMax, dcaxyBins, dcaxyMin, dcaxyMax, Form("Fake %s from materials after PID", particle[i].Data()), "p_{T} (GeV/c)", "sign #times DCA_{xy} (cm)");
+ }
+
+ // identification
+
+ hMap->Add( particle[i] + "_PID_ITSdEdx_P", 1000, 1.e-3, 10., 1500, 1.e-3, 1500., Form("%s candidate",particle[i].Data()), "p (GeV/c)", "dE/dx");
+
+ hMap->Add( particle[i] + "_PID_TPCdEdx_P", 1000, 1.e-3, 10., 1500, 1.e-3, 1500., Form("%s candidate",particle[i].Data()), "p_{TPC} (GeV/c)", "dE/dx");
+
+ hMap->Add( particle[i] + "_PID_Beta_P", 1000, 1.e-3, 10., 1200, 1.e-3, 1.2, Form("%s candidate",particle[i].Data()), "p_{TOF} (GeV/c)", "\\beta");
+
+ hMap->Add( particle[i] + "_PID_Mass_P", 1000, 1.e-3, 10., 500, 0.01, 5., Form("%s candidate",particle[i].Data()), "p_{TOF} (GeV/c)", "m (GeV/c^{2})");
+
+ if(simulation)
+ {
+ hMap->Add( particle[i] + "_Sim_PID_Mass", 500, 0.01, 5., Form("%s after PID", particle[i].Data()), "m (GeV/c^{2})");
+ }
+
+ // results
+
+ hMap->Add( particle[i] + "_PID_Pt_Y", etaBins, etaMin, etaMax, ptBins, ptMin, ptMax, Form("%s candidate",particle[i].Data()), "y", "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_PID_Pt", ptBins, ptMin, ptMax, Form("%s candidate (|y| < %0.1f, 0 < \\phi < 2\\pi)", particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( Form("%s_PID_Y",particle[i].Data()), etaBins, etaMin, etaMax, Form("%s candidate (p_{T} > 0, 0 < \\phi < 2\\pi)",particle[i].Data()), "y");
+
+ hMap->Add( particle[i] + "_PID_Phi", phiBins, phiMin, phiMax, Form("%s candidate (p_{T} > 0, |y| < %0.1f)",particle[i].Data(),maxY), "\\phi (rad)");
+
+ if(simulation)
+ {
+ hMap->Add( particle[i] + "_Sim_Pt", ptBins, ptMin, ptMax, Form("%s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Prim_Phi", phiBins, phiMin, phiMax, Form("Primary %s (p_{T} > 0, |y| < %0.1f)",particle[i].Data(),maxY), "\\phi (rad)");
+
+ hMap->Add( particle[i] + "_Sim_Prim_Y", etaBins, etaMin, etaMax, Form("Primary %s (p_{T} > 0, 0 < \\phi < 2\\pi)",particle[i].Data()), "y");
+
+ hMap->Add( particle[i] + "_Sim_Fdwn_Pt", ptBins, ptMin, ptMax, Form("Feed-down %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Mat_Pt", ptBins, ptMin, ptMax, Form("%s from materials (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s after PID (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Pt", ptBins, ptMin, ptMax, Form("%s after PID (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Prim_R", 200, 0., 200., Form("%s after PID (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "R_{xyz} (cm)");
+
+ hMap->Add( particle[i] + "_Sim_PtY", etaBins, etaMin, etaMax, ptBins, ptMin, ptMax, particle[i].Data(), "y", "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Prim_M2_P", ptBins, ptMin, ptMax, m2Bins, m2Min, m2Max, Form("Primary %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p (GeV/c)", "m^{2} (GeV^{2}/c^{4})");
+
+ hMap->Add( particle[i] + "_Sim_Prim_M2_Pt", ptBins, ptMin, ptMax, 100, 0., 10., Form("Primary %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)", "m^{2} (GeV^{2}/c^{4})");
+
+ // unfolding
+
+ hMap->Add( particle[i] + "_Response_Matrix", ptBins, ptMin, ptMax, ptBins, ptMin, ptMax, particle[i].Data(), "Measured p_{T} (GeV/c)", "True p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Prim_Response_Matrix", ptBins, ptMin, ptMax, ptBins, ptMin, ptMax, Form("Primary %s",particle[i].Data()), "Measured p_{T} (GeV/c)", "True p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Prim_DiffPt_RecPt", ptBins, ptMin, ptMax, 1000, -0.5, 0.5, Form("Primary %s",particle[i].Data()), "p_{T}^{rec} (GeV/c)", "p_{T}^{gen}-p_{T}^{rec} (GeV/c)");
+
+ // fake tracks
+ hMap->Add( particle[i] + "_Sim_Fake_Pt", ptBins, ptMin, ptMax, Form("Fake %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Fake_Prim_Pt", ptBins, ptMin, ptMax, Form("Fake Primary %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Fake_Fdwn_Pt", ptBins, ptMin, ptMax, Form("Fake Feed-down %s (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_Fake_Mat_Pt", ptBins, ptMin, ptMax, Form("Fake %s from materials (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Sim_PID_Fake_Pt", ptBins, ptMin, ptMax, Form("Fake %s after PID (|y| < %0.1f, 0 < \\phi < 2\\pi)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+ }
+ }
+
+ // generation
+
+ if(simulation)
+ {
+ for(int i=0; i<2; ++i)
+ {
+ hMap->Add( particle[i] + "_Gen_Prim_P", ptBins, ptMin, ptMax, Form("Primary %s", particle[i].Data()), "p (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s", particle[i].Data()), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Prim_Y", 600, -12, 12, Form("Primary %s", particle[i].Data()), "y");
+
+ hMap->Add( particle[i] + "_Gen_Prim_Eta", 600, -12, 12, Form("Primary %s", particle[i].Data()), "#eta");
+
+ hMap->Add( particle[i] + "_Gen_Prim_Phi", phiBins, phiMin, phiMax, Form("Primary %s", particle[i].Data()), "\\phi (rad)");
+
+ hMap->Add( particle[i] + "_Gen_Prim_PtY", etaBins, etaMin, etaMax, ptBins, ptMin, ptMax, Form("Primary %s", particle[i].Data()), "y", "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Prim_EtaY", etaBins, etaMin, etaMax, etaBins, etaMin, etaMax, Form("Primary %s", particle[i].Data()), "y", "#eta");
+
+ // after multiplicity cut
+
+ hMap->Add( particle[i] + "_Gen_Mult_Prim_PtY", etaBins, etaMin, etaMax, ptBins, ptMin, ptMax, Form("Primary %s", particle[i].Data()), "y", "p_{T} (GeV/c)");
+
+ // in the phase space within the acceptance
+
+ hMap->Add( particle[i] + "_Gen_PhS_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (all events, |y| < %0.1f)",particle[i].Data(),maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Trig_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (triggered events, |y| < %0.1f)", particle[i].Data(), maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Vtx_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (good vertex, |y| < %0.1f)", particle[i].Data(), maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_NoTrig_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (no triggered events, |y| < %0.1f)", particle[i].Data(), maxY), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_NoVtx_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (no reconstructed vertex, |y| < %0.1f)", particle[i].Data(), maxY), "p_{T} (GeV/c)");
+
+ // within the acceptance
+
+ hMap->Add( particle[i] + "_Gen_Acc_Prim_P", ptBins, ptMin, ptMax, Form("Primary %s (|y| < %0.1f and |\\eta| < %0.1f)",particle[i].Data(),maxY,maxEta), "p (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Acc_Prim_Pt", ptBins, ptMin, ptMax, Form("Primary %s (|y| < %0.1f and |\\eta| < %0.1f)",particle[i].Data(),maxY,maxEta), "p_{T} (GeV/c)");
+
+ hMap->Add( particle[i] + "_Gen_Acc_Prim_Y", etaBins, etaMin, etaMax, Form("Primary %s (|\\eta| < %0.1f)",particle[i].Data(),maxEta), "y");
+
+ hMap->Add( particle[i] + "_Gen_Acc_Prim_Phi", phiBins, phiMin, phiMax, Form("Primary %s (|y|< %0.1f and |\\eta| < %0.1f)",particle[i].Data(),maxY,maxEta), "\\phi (rad)");
+
+ hMap->Add( particle[i] + "_Gen_Acc_Prim_PtY", etaBins, etaMin, etaMax, ptBins, ptMin, ptMax, Form("Primary %s",particle[i].Data()), "y", "p_{T} (GeV/c)");
+ }
+ }
+
+ return hMap;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// macro for prior probabilities
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+void SetParameters(double* prob, double c0, double c1, double c2, double c3, double c4, double c5, double c6, double c7, double c8)
+{
+//
+// assign values to an array
+//
+ prob[0] = c0;
+ prob[1] = c1;
+ prob[2] = c2;
+ prob[3] = c3;
+ prob[4] = c4;
+ prob[5] = c5;
+ prob[6] = c6;
+ prob[7] = c7;
+ prob[8] = c8;
+}
+
+Bool_t PriorProbabilities(Double_t* prob, const TString& periodname, const TString& tracksel)
+{
+//
+// Prior probabilities for the different samples and track selection criteria.
+// Obtained from several iterations of the analysis program
+//
+ TString period = periodname;
+ period.ToLower();
+
+ TString trksel = tracksel;
+ trksel.ToLower();
+
+ if(period == "lhc10c900") // pp 900 GeV
+ {
+ if (trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0229385, 2.13852e-06, 0.868889, 0.0596281, 0.0484369, 0.00010232, 2.54977e-06, 0, 0);
+ }
+ else if (trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00670787, 1.03214e-06, 0.887492, 0.0627365, 0.0429736, 8.72157e-05, 2.06428e-06, 0, 0);
+ }
+ else if(trksel=="its_tpc_nsigma_spd")
+ {
+ SetParameters(prob, 0.00598811, 0.000283397, 0.889894, 0.0655886, 0.0382077, 3.76161e-05, 8.51043e-08, 0, 0);
+ }
+ else if(trksel == "its_tpc_tof_nsigma")
+ {
+ SetParameters(prob, 0.00968656, 0.00230198, 0.844375, 0.084224, 0.059334, 7.8127e-05, 7.89161e-07, 0, 0);
+ }
+ else if(trksel == "its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0123784, 0.00368129, 0.830011, 0.0795686, 0.0742296, 0.000127461, 2.99204e-06, 1.19682e-07, 0);
+ }
+ else // its_tpc_nsigma
+ {
+ SetParameters(prob, 0.0181447, 1.09966e-06, 0.880063, 0.0635145, 0.0382359, 4.06874e-05, 1.83277e-07, 0, 0);
+ }
+ }
+
+ else if( period.Contains("lhc10b")
+ || period.Contains("lhc10c")
+ || period.Contains("lhc10d")
+ || period.Contains("lhc10e") ) // pp 7 TeV
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0250199, 0.00330477, 0.856905, 0.0621491, 0.0524821, 0.000136071, 3.44981e-06, 6.86137e-08, 0);
+ }
+ else if(trksel=="its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00883871, 0.00306025, 0.877528, 0.0646921, 0.0457766, 0.000102032, 2.23315e-06, 3.71013e-08, 0);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0155752, 0.0130319, 0.813626, 0.0808387, 0.0767456, 0.000177436, 4.32771e-06, 3.93428e-07, 0);
+ }
+ else // its_tpc_nsigma
+ {
+ SetParameters(prob, 0.0207349, 0.000587565, 0.871052, 0.0673557, 0.0402021, 6.74649e-05, 0, 0, 0);
+ }
+ }
+
+ else if(period == "lhc10h") // heavy ions 2.76 TeV, 20% centrality
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.032994, 0.012492, 0.824044, 0.0809969, 0.0492257, 0.000242368, 5.01681e-06, 3.34454e-07, 0);
+ }
+ else if(trksel=="its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.0262206, 0.00398172, 0.845586, 0.081878, 0.0421142, 0.000215343, 3.92153e-06, 0, 0);
+ }
+ else if(trksel=="its_tpc_nsigma")
+ {
+ SetParameters(prob, 0.0279673, 0.00329778, 0.845802, 0.085995, 0.0367557, 0.000179264, 2.62155e-06, 2.49671e-07, 1.24836e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca")
+ {
+ SetParameters(prob, 0.017615, 0.00752739, 0.784779, 0.109672, 0.080031, 0.000370991, 3.95272e-06, 8.47012e-07, 0);
+ }
+ else if(trksel=="its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0130896, 0.0079917, 0.797236, 0.112833, 0.0685078, 0.000338538, 3.07762e-06, 0, 0);
+ }
+ else if(trksel=="its_tpc_tof_nsigma")
+ {
+ SetParameters(prob, 0.0157654, 0.00721323, 0.79879, 0.117158, 0.0607624, 0.000308557, 1.26199e-06, 6.30996e-07, 0);
+ }
+ else
+ {
+ SetParameters(prob, 0.032994, 0.012492, 0.824044, 0.0809969, 0.0492257, 0.000242368, 5.01681e-06, 3.34454e-07, 3.34454e-07);
+ }
+ }
+
+ else if(period == "lhc11a_wosdd") // pp 2.76 TeV without SDD pass3
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0172658, 0.0179012, 0.857385, 0.0590938, 0.0482499, 0.000104021, 8.05072e-07, 0, 0);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00808143, 0.016243, 0.873576, 0.0593675, 0.0426588, 7.27021e-05, 6.45822e-07, 0, 0);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0113244, 0.0275349, 0.790324, 0.0899881, 0.0806648, 0.000161212, 2.45189e-06, 0, 0);
+ }
+ else
+ {
+ SetParameters(prob, 0.0132164, 0.0209664, 0.864317, 0.0624972, 0.0389603, 4.30391e-05, 0, 0, 0);
+ }
+ }
+
+ else if(period == "lhc11a_wsdd") // pp 2.76 TeV without SDD pass3
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0198007, 0.0130987, 0.857079, 0.0604184, 0.0494965, 0.000104874, 1.5572e-06, 0, 0);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00808143, 0.016243, 0.873576, 0.0593675, 0.0426588, 7.27021e-05, 6.45822e-07, 0, 0);
+ }
+ else if(trksel == "its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0206518, 0.00922775, 0.859646, 0.0595307, 0.0508038, 0.000137122, 3.02079e-06, 0, 0);
+ }
+ else
+ {
+ SetParameters(prob, 0.0150916, 0.0122592, 0.868485, 0.0649622, 0.0391601, 4.18627e-05, 0, 0, 0);
+ }
+ }
+
+ // --------- simulation ----------
+
+ else if(period == "lhc12a5a")
+ {
+ if(trksel == "its_tpc_dca")
+ {
+ SetParameters(prob, 0.00269435, 0.00142106, 0.0831012, 0.00522487, 0.371055, 0.332289, 0.0602906, 0.0830945, 0.0608297);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.000903954, 0.000740461, 0.0787924, 0.0052352, 0.377698, 0.331812, 0.0598412, 0.0834821, 0.0614951);
+ }
+ else if(trksel == "its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.000699209, 0.000628928, 0.0488789, 0.00274458, 0.378142, 0.35414, 0.0643921, 0.0830996, 0.0672741);
+ }
+ else
+ {
+ SetParameters(prob, 0.00177397, 0.000490646, 0.0745835, 0.00552979, 0.383725, 0.331756, 0.0592882, 0.0827403, 0.0601128);
+ }
+ }
+
+ else if(period == "lhc12a5bb")
+ {
+ if(trksel == "its_tpc_dca")
+ {
+ SetParameters(prob, 0.00370442, 0.00158498, 0.108206, 0.00792265, 0.375296, 0.33667, 0.0610851, 0.0715099, 0.034021);
+ }
+ else if(trksel == "its_tpc_tof_dca")
+ {
+ SetParameters(prob, 0.00112885, 0.000832592, 0.0672591, 0.00469105, 0.385956, 0.361319, 0.0655949, 0.0767322, 0.0364865);
+ }
+ }
+
+ else if(period == "lhc12a5bc")
+ {
+ if(trksel == "its_tpc_dca")
+ {
+ SetParameters(prob, 0.00370442, 0.00158498, 0.108206, 0.00792265, 0.375296, 0.33667, 0.0610851, 0.0715099, 0.034021);
+ }
+ else if(trksel == "its_tpc_tof_dca")
+ {
+ SetParameters(prob, 0.00112885, 0.000832592, 0.0672591, 0.00469105, 0.385956, 0.361319, 0.0655949, 0.0767322, 0.0364865);
+ }
+ }
+
+ else if(period == "lhc12a5bd")
+ {
+ if(trksel == "its_tpc_dca")
+ {
+ SetParameters(prob, 0.00389124, 0.00181114, 0.10725, 0.00753775, 0.361361, 0.322078, 0.0584554, 0.0785471, 0.0590672);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00129393, 0.000963743, 0.102737, 0.00760901, 0.36841, 0.322256, 0.0581189, 0.0791033, 0.0595082);
+ }
+ else if(trksel == "its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.00122986, 0.000908924, 0.0696766, 0.00462655, 0.385381, 0.360127, 0.0653381, 0.0762508, 0.0364611);
+ }
+ else
+ {
+ SetParameters(prob, 0.00265834, 0.000648637, 0.0985374, 0.00801093, 0.373786, 0.322076, 0.0576089, 0.0782477, 0.0584268);
+ }
+ }
+
+ else if(period == "lhc12a5be")
+ {
+ if(trksel == "its_tpc_tof_dca")
+ {
+ SetParameters(prob, 0.00114523, 0.000795738, 0.0666132, 0.00465889, 0.386032, 0.36171, 0.0656093, 0.0767938, 0.0366419);
+ }
+ else
+ {
+ SetParameters(prob, 0.00114523, 0.000795738, 0.0666132, 0.00465889, 0.386032, 0.36171, 0.0656093, 0.0767938, 0.0366419);
+ }
+ }
+
+ else if(period == "lhc12a5c_wosdd" || period == "lhc12a5c_wsdd")
+ {
+ if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00140969, 0.00100911, 0.0901532, 0.00627527, 0.375878, 0.333126, 0.0601881, 0.0752193, 0.0567413);
+ }
+ else if(trksel == "its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.00254534, 0.00138188, 0.0920481, 0.00642957, 0.381436, 0.344603, 0.0627622, 0.0736729, 0.0351211);
+ }
+ else
+ {
+ SetParameters(prob, 0.00254534, 0.00138188, 0.0920481, 0.00642957, 0.381436, 0.344603, 0.0627622, 0.0736729, 0.0351211);
+ }
+ }
+
+ else if(period == "lhc10e12")
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0274572, 0.0132276, 0.868083, 0.0579849, 0.0332055, 2.94442e-05, 1.1225e-05, 0, 8.63465e-07);
+ }
+ else if(trksel=="its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00902306, 0.00879859, 0.891372, 0.0600273, 0.0307466, 2.43343e-05, 7.44928e-06, 0, 4.96619e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca")
+ {
+ SetParameters(prob, 0.0113533, 0.0104639, 0.877447, 0.0577014, 0.0429803, 3.77285e-05, 1.40853e-05, 0, 2.51523e-06);
+ }
+ else
+ {
+ SetParameters(prob, 0.0228473, 0.00644153, 0.881383, 0.061429, 0.0278978, 1.21828e-06, 2.81142e-07, 0, 0);
+ }
+ }
+
+ else if(period == "lhc10e13")
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0275037, 0.0129471, 0.864088, 0.0566767, 0.03875, 2.46576e-05, 9.52441e-06, 0, 5.29134e-07);
+ }
+ else if(trksel=="its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.00881856, 0.00853512, 0.887761, 0.0587557, 0.03611, 1.46288e-05, 5.01096e-06, 0, 4.0411e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0118565, 0.0102582, 0.869917, 0.0579405, 0.0499923, 2.71412e-05, 8.44393e-06, 0, 0);
+ }
+ else
+ {
+ SetParameters(prob, 0.0227488, 0.00644463, 0.877399, 0.0598542, 0.033552, 1.08815e-06, 3.43625e-07, 0, 0);
+ }
+ }
+
+ else if(period == "lhc10f6a")
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0301784, 0.0130575, 0.849464, 0.0631854, 0.0440839, 1.65739e-05, 8.28697e-06, 0, 5.52465e-06);
+ }
+ else if(trksel=="its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.0101971, 0.00853634, 0.875408, 0.0657018, 0.0401396, 1.18943e-05, 4.23147e-06, 0, 8.53405e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0149055, 0.0105359, 0.852053, 0.0663013, 0.0561727, 2.17209e-05, 9.18959e-06, 0, 8.35418e-07);
+ }
+ else
+ {
+ SetParameters(prob, 0.0248513, 0.00652446, 0.863878, 0.0670812, 0.0376654, 3.27516e-05, 2.45637e-05, 0, 5.4586e-06);
+ }
+ }
+
+ else if(period == "lhc11e3a_wsdd") // from the grid without SDD
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0288269, 0.0143436, 0.85568, 0.0600814, 0.0410309, 2.04805e-05, 1.56616e-05, 0, 1.20474e-06);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.0104567, 0.00998477, 0.878832, 0.0622996, 0.0384064, 1.43839e-05, 5.60782e-06, 0, 7.60383e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0129161, 0.0117642, 0.852733, 0.0671011, 0.0554566, 2.17717e-05, 5.44293e-06, 0, 1.36073e-06);
+ }
+ else
+ {
+ SetParameters(prob, 0.0249858, 0.00773001, 0.87308, 0.0601418, 0.0340623, 0, 0, 0, 0);
+ }
+ }
+
+ else if(period == "lhc11e3a_wosdd") // from the grid without SDD
+ {
+ if(trksel=="its_tpc_dca")
+ {
+ SetParameters(prob, 0.0222021, 0.0131403, 0.862641, 0.0613467, 0.0406294, 2.3342e-05, 1.55613e-05, 0, 1.29678e-06);
+ }
+ else if(trksel == "its_tpc_dca_spd")
+ {
+ SetParameters(prob, 0.0104567, 0.00998477, 0.878832, 0.0622996, 0.0384064, 1.43839e-05, 5.60782e-06, 0, 7.60383e-07);
+ }
+ else if(trksel=="its_tpc_tof_dca" || trksel == "its_tpc_tof_dca_spd")
+ {
+ SetParameters(prob, 0.0118237, 0.0110965, 0.855035, 0.0672175, 0.0547956, 2.49547e-05, 6.23867e-06, 0, 6.93185e-07);
+ }
+ else
+ {
+ SetParameters(prob, 0.0186686, 0.00736997, 0.875823, 0.062905, 0.0352281, 4.98982e-06, 0, 0, 0);
+ }
+ }
+
+ else
+ {
+ std::cout << "Setting equal probabilities to all particle species: " << periodname << ", " << trksel << std::endl;
+
+ SetParameters(prob, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.11);
+
+ return kFALSE;
+ }
+
+ return kTRUE;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// macro for a predefined set of track cuts
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+void AddNSigmaCuts(AliESDtrackCuts* trkCuts, Float_t maxNSigma=3)
+{
+//
+// nsigma cuts
+//
+ trkCuts->SetRequireSigmaToVertex(kTRUE);
+ trkCuts->SetMaxNsigmaToVertex(maxNSigma);
+}
+
+void AddDCACuts(AliESDtrackCuts* trkCuts, Float_t maxDCAxy=1.5, Float_t maxDCAz=2)
+{
+//
+// DCA cuts
+//
+ trkCuts->SetMaxDCAToVertexXY(maxDCAxy);
+ trkCuts->SetMaxDCAToVertexZ(maxDCAz);
+}
+
+AliESDtrackCuts* TrackCuts(AliAnalysisTaskB2* task, const TString& trackselection, Double_t maxDCAxy, Double_t maxDCAz, Double_t maxNSigma, Int_t minTPCnCls, Double_t maxEta)
+{
+//
+// Create an AliESDtrackCuts from a predefined set
+//
+ AliESDtrackCuts* trkCuts = new AliESDtrackCuts("AliESDtrackCuts");
+
+ trkCuts->SetEtaRange(-maxEta, maxEta);
+ trkCuts->SetPtRange(0.15, 100.);
+ trkCuts->SetAcceptKinkDaughters(kFALSE);
+
+ // ITS
+ trkCuts->SetRequireITSRefit(kTRUE);
+ trkCuts->SetMinNClustersITS(2);
+ trkCuts->SetMaxChi2PerClusterITS(36);
+
+ // TPC
+ trkCuts->SetRequireTPCRefit(kTRUE);
+ trkCuts->SetMinNClustersTPC(minTPCnCls);
+ trkCuts->SetMaxChi2PerClusterTPC(4.);
+
+ TString tracksel = trackselection;
+ tracksel.ToLower();
+
+ if(tracksel == "its_tpc_nsigma")
+ {
+ AddNSigmaCuts(trkCuts, maxNSigma);
+ }
+
+ else if(tracksel == "its_tpc_nsigma_spd")
+ {
+ AddNSigmaCuts(trkCuts, maxNSigma);
+ trkCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
+ }
+
+ else if(tracksel == "its_tpc_dca")
+ {
+ AddDCACuts(trkCuts,maxDCAxy,maxDCAz);
+ }
+
+ else if(tracksel == "its_tpc_dca_spd")
+ {
+ AddDCACuts(trkCuts,maxDCAxy,maxDCAz);
+ trkCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
+ }
+
+ else if(tracksel == "its_tpc_tof_nsigma")
+ {
+ AddNSigmaCuts(trkCuts, maxNSigma);
+ task->SetTOFmatch(1);
+ }
+
+ else if(tracksel == "its_tpc_tof_nsigma_spd")
+ {
+ AddNSigmaCuts(trkCuts, maxNSigma);
+ task->SetTOFmatch(1);
+ trkCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
+ }
+
+ else if(tracksel == "its_tpc_tof_dca")
+ {
+ AddDCACuts(trkCuts,maxDCAxy,maxDCAz);
+ task->SetTOFmatch(1);
+ }
+
+ else if(tracksel == "its_tpc_tof_dca_spd")
+ {
+ AddDCACuts(trkCuts,maxDCAxy,maxDCAz);
+ task->SetTOFmatch(1);
+ trkCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
+ }
+
+ else if(tracksel == "std_its_tpc_2009") // pp data 2009
+ {
+ delete trkCuts;
+ return AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(1);
+ }
+
+ else if(tracksel == "std_its_tpc_2010") // pp data 2010
+ {
+ delete trkCuts;
+ return AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
+ }
+
+ else
+ {
+ std::cerr << "Warning: no track selection criteria selected" << std::endl;
+ }
+
+ return trkCuts;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff