Adding Eulogios analysis
authorakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Dec 2012 13:36:34 +0000 (13:36 +0000)
committerakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 17 Dec 2012 13:36:34 +0000 (13:36 +0000)
PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnID.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/AddTaskB2.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/BetheBlochParams.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/CreateHistograms.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/PriorProbabilities.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/TrackCuts.C [new file with mode: 0644]

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.cxx
new file mode 100644 (file)
index 0000000..a210682
--- /dev/null
@@ -0,0 +1,1232 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h b/PWGLF/SPECTRA/Nuclei/B2/AliAnalysisTaskB2.h
new file mode 100644 (file)
index 0000000..871427c
--- /dev/null
@@ -0,0 +1,175 @@
+#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
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.cxx
new file mode 100644 (file)
index 0000000..4e27e7f
--- /dev/null
@@ -0,0 +1,138 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnHistoMap.h
new file mode 100644 (file)
index 0000000..cf4c3c5
--- /dev/null
@@ -0,0 +1,48 @@
+#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
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.cxx
new file mode 100644 (file)
index 0000000..217a8a6
--- /dev/null
@@ -0,0 +1,673 @@
+/**************************************************************************
+ * 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(&param[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(&param[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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnID.h
new file mode 100644 (file)
index 0000000..9bd08a7
--- /dev/null
@@ -0,0 +1,88 @@
+#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
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/AddTaskB2.C b/PWGLF/SPECTRA/Nuclei/B2/macros/AddTaskB2.C
new file mode 100644 (file)
index 0000000..e65b26b
--- /dev/null
@@ -0,0 +1,172 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/BetheBlochParams.C b/PWGLF/SPECTRA/Nuclei/B2/macros/BetheBlochParams.C
new file mode 100644 (file)
index 0000000..fa6819e
--- /dev/null
@@ -0,0 +1,83 @@
+/**************************************************************************
+ * 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);
+       }
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/CreateHistograms.C b/PWGLF/SPECTRA/Nuclei/B2/macros/CreateHistograms.C
new file mode 100644 (file)
index 0000000..09315e8
--- /dev/null
@@ -0,0 +1,365 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/PriorProbabilities.C b/PWGLF/SPECTRA/Nuclei/B2/macros/PriorProbabilities.C
new file mode 100644 (file)
index 0000000..bf9d2aa
--- /dev/null
@@ -0,0 +1,374 @@
+/**************************************************************************
+ * 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;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/TrackCuts.C b/PWGLF/SPECTRA/Nuclei/B2/macros/TrackCuts.C
new file mode 100644 (file)
index 0000000..8707567
--- /dev/null
@@ -0,0 +1,127 @@
+/**************************************************************************
+ * 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;
+}