--- /dev/null
+//-*- Mode: C++ -*-
+
+#include "TMath.h"
+#include "TAxis.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+
+#include "AliESDtrackCuts.h"
+
+#include "AliAnalysisNetParticleDCA.h"
+
+using namespace std;
+
+// Task for NetParticle checks
+// Author: Jochen Thaeder <jochen@thaeder.de>
+
+ClassImp(AliAnalysisNetParticleDCA)
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+AliAnalysisNetParticleDCA::AliAnalysisNetParticleDCA() :
+ fHelper(NULL),
+
+ fESD(NULL),
+ fESDTrackCuts(NULL),
+ fESDTrackCutsBkg(NULL),
+
+ fStack(NULL),
+ fMCEvent(NULL),
+
+ fHnDCA(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliAnalysisNetParticleDCA",10);
+}
+
+//________________________________________________________________________
+AliAnalysisNetParticleDCA::~AliAnalysisNetParticleDCA() {
+ // Destructor
+
+}
+
+
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::Initialize(AliESDtrackCuts *cuts, AliESDtrackCuts *cutsBkg, AliAnalysisNetParticleHelper* helper) {
+
+ // -- Get Helper
+ // ---------------
+ fHelper = helper;
+
+ // -- ESD track cuts
+ // -------------------
+ fESDTrackCuts = cuts;
+ fESDTrackCutsBkg = cutsBkg;
+
+#if 0
+ // -- Get particle species / pdgCode
+ // -------------------------
+ fPdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies());
+
+ // -- MC Labels for efficiency
+ // -----------------------------
+ fLabelsRec = new Int_t*[2];
+ for (Int_t ii = 0; ii < 2 ; ++ii)
+ fLabelsRec[ii] = NULL;
+
+ // -- Create THnSparse Histograms
+ // --------------------------------
+ CreateHistograms();
+#endif
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDCA::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
+ // -- Setup event
+
+ // -- Reset of event
+ // -------------------
+ ResetEvent();
+
+ // -- Setup of event
+ // -------------------
+ fESD = esdHandler->GetEvent();
+
+ fMCEvent = mcEvent;
+ fStack = fMCEvent->Stack();
+#if 0
+ fCentralityBin = centBin;
+
+ // -- Create label arrays
+ // ------------------------
+ fLabelsRec[0] = new Int_t[fNTracks];
+ if(!fLabelsRec[0]) {
+ AliError("Cannot create fLabelsRec[0]");
+ return -1;
+ }
+
+ fLabelsRec[1] = new Int_t[fNTracks];
+ if(!fLabelsRec[1]) {
+ AliError("Cannot create fLabelsRec[1] for TPC");
+ return -1;
+ }
+
+ for(Int_t ii = 0; ii < fNTracks; ++ii) {
+ fLabelsRec[0][ii] = 0;
+ fLabelsRec[1][ii] = 0;
+ }
+#endif
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::ResetEvent() {
+ // -- Reset event
+ /*
+ for (Int_t ii = 0; ii < 2 ; ++ii) {
+ if (fLabelsRec[ii])
+ delete[] fLabelsRec[ii];
+ fLabelsRec[ii] = NULL;
+ }
+ */
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::Process() {
+ // -- Process event
+
+ // -- Setup (clean, create and fill) MC labels
+ // ---------------------------------------------
+ // FillMCLabels();
+
+ // -- Fill MC histograms for efficiency studies
+ // -----------------------------------------------
+ // FillMCEffHist();
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Private Methods
+ * ---------------------------------------------------------------------------------
+ */
+#if 0
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::CreateHistograms() {
+ // -- Create histograms
+
+ Double_t dCent[2] = {-0.5, 8.5};
+ Int_t iCent = 9;
+
+ Double_t dEta[2] = {-0.9, 0.9};
+ Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
+
+ Double_t dRap[2] = {-0.5, 0.5};
+ Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
+
+ Double_t dPhi[2] = {0.0, TMath::TwoPi()};
+ Int_t iPhi = 42;
+
+ Double_t dPt[2] = {0.1, 3.0};
+ Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
+
+ Double_t dSign[2] = {-1.5, 1.5};
+ Int_t iSign = 3;
+
+ // ------------------------------------------------------------------
+ // -- Create THnSparseF - DCA
+ // ------------------------------------------------------------------
+
+ // cent: etaMC: yMC: phiMC: ptMC: sign:contPart:DCArAccepted:DCAr
+ Int_t binHnDCA[9] = { iCent, iEta, iRap, iPhi, iPt, iSign, 4, 2, 50 };
+ Double_t minHnDCA[9] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], 0.5, -0.5, -10.};
+ Double_t maxHnDCA[9] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 4.5, 1.5, 10.};
+ fHnDCA = new THnSparseF("fHnDCA", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:DCArAccepted:DCAr", 9, binHnDCA, minHnDCA, maxHnDCA);
+
+ fHnDCA->Sumw2();
+ fHnDCA->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
+ fHnDCA->GetAxis(1)->SetTitle("#eta_{MC}"); // eta [-0.9,0.9]
+ fHnDCA->GetAxis(2)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
+ fHnDCA->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. ,2Pi]
+ fHnDCA->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pT [ 0.1,1.3]
+ fHnDCA->GetAxis(5)->SetTitle("sign"); // -1 | 0 | +1
+ fHnDCA->GetAxis(6)->SetTitle("contPart"); // 1 primary | 2 missId | 3 from WeakDecay | 4 p from Material
+ fHnDCA->GetAxis(7)->SetTitle("DCArAccepted"); // 0 not accepted | 1 accepted
+ fHnDCA->GetAxis(8)->SetTitle("DCAr"); // DCAr [-10,10]
+
+ // ------------------------------------------------------------------
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::FillMCDCA() {
+ // Fill MC labels
+ // Loop over ESD tracks and fill arrays with MC lables
+ // fLabelsRec[0] : all Tracks
+ // fLabelsRec[1] : all Tracks accepted by PID of TPC
+ // Check every accepted track if correctly identified
+ // otherwise check for contamination
+
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted
+ if (!fESDTrackCuts->AcceptTrack(track))
+ continue;
+
+ // -- Check if accepted in rapidity window
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP))
+ continue;
+
+ // -- Check if accepted with thighter DCA cuts
+ if (!fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ Int_t label = TMath::Abs(track->GetLabel());
+
+ // -- Fill Label of all reconstructed
+ fLabelsRec[0][idxTrack] = label;
+
+ // -- Check if accepted by PID from TPC or TPC+TOF
+ Double_t pid[2];
+ if (!fHelper->IsTrackAcceptedPID(track, pid))
+ continue;
+
+ // -- Fill Label of all reconstructed && recPid_TPC+TOF
+ fLabelsRec[1][idxTrack] = label;
+
+ // -- Check for contamination and fill contamination THnSparse
+ CheckContTrack(label, track->GetSign(), idxTrack);
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDCA::CheckDCATrack(Int_t label, Float_t sign, Float_t isDCArAccepted, Float_t dcar) {
+ // Check if DCA of contamination or correctly identified
+ // Fill contamination DCA THnSparse
+
+ TParticle* particle = fStack->Particle(label);
+ if (!particle)
+ return;
+
+ Bool_t isSecondaryFromWeakDecay = kFALSE;
+ Bool_t isSecondaryFromMaterial = kFALSE;
+
+ Int_t contPart = 0;
+
+ // -- Check if correctly identified
+ if (particle->GetPdgCode() == (sign*fPdgCode)) {
+
+ // Check if is physical primary -> all ok
+ if (fStack->IsPhysicalPrimary(label))
+ contPart = 1;
+ // -- Check if secondaries from weak decay
+ else if(isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label))
+ contPart = 3;
+ // -- Check if secondaries from material decay
+ else if (isSecondaryFromMaterial = fStack->IsSecondaryFromMaterial(label))
+ contPart = 4;
+ }
+ // -- MissIdentification
+ else
+ contPart = 2;
+
+ Double_t hnDCA[9] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, contPart, isDCArAccepted, dcar};
+ fHnDCA->Fill(hnContDCA);
+
+ return;
+}
+
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#ifndef ALIANALYSISNETPARTICLEDCA_H
+#define ALIANALYSISNETPARTICLEDCA_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Efficiency and Contaminations for NetParticle Distributions
+// Authors: Jochen Thaeder <jochen@thaeder.de>
+
+#include "THnSparse.h"
+
+#include "AliAnalysisNetParticleHelper.h"
+
+class AliESDEvent;
+class AliESDInputHandler;
+class AliMCEvent;
+class AliStack;
+
+class AliAnalysisNetParticleDCA : public TNamed {
+
+ public:
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleDCA();
+ virtual ~AliAnalysisNetParticleDCA();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize */
+ void Initialize(AliESDtrackCuts *cuts, AliESDtrackCuts *cutsBkg, AliAnalysisNetParticleHelper* helper);
+
+ /** Setup Event */
+ Int_t SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent);
+
+ /** Reset Event */
+ void ResetEvent();
+
+ /** Process Event */
+ void Process();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Getter
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Get Ptr to DCA THnSparse */
+ THnSparseF* GetHnDCA() {return fHnDCA;}
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliAnalysisNetParticleDCA(const AliAnalysisNetParticleDCA&); // not implemented
+ AliAnalysisNetParticleDCA& operator=(const AliAnalysisNetParticleDCA&); // not implemented
+
+
+#if 0
+ /*
+ * ---------------------------------------------------------------------------------
+ * Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Create the efficiency / contamination THnSparseF */
+ void CreateHistograms();
+
+ /** Fill DCA ThnSparse */
+ void FillDCA();
+
+ /** Check if particle is contamination */
+ void CheckDCATrack(Int_t label, Float_t sign, Int_t idxTrack);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+#endif
+ AliAnalysisNetParticleHelper *fHelper; //! Ptr to helper class
+
+ // --- ESD only ----------------------------------------------------------
+
+ AliESDEvent *fESD; //! ESD object
+ AliESDtrackCuts *fESDTrackCuts; //! ESD cuts
+ AliESDtrackCuts *fESDTrackCutsBkg; //! ESD cuts
+
+ // --- MC only -----------------------------------------------------------
+
+ AliStack *fStack; //! Ptr to stack
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+
+ // -----------------------------------------------------------------------
+
+ THnSparseF *fHnDCA; // THnSparseF contamination DCA
+
+ // -----------------------------------------------------------------------
+
+ ClassDef(AliAnalysisNetParticleDCA, 1);
+};
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#include "TMath.h"
+#include "TAxis.h"
+#include "TSystem.h"
+#include "TProfile.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TFile.h"
+#include "TPRegexp.h"
+
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDInputHandler.h"
+#include "AliESDpid.h"
+#include "AliCentrality.h"
+#include "AliTracker.h"
+
+#include "AliAnalysisNetParticleDistribution.h"
+
+using namespace std;
+
+// Task for NetParticle checks
+// Author: Jochen Thaeder <jochen@thaeder.de>
+
+ClassImp(AliAnalysisNetParticleDistribution)
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() :
+ fHelper(NULL),
+
+ fOutList(NULL),
+
+ fESDHandler(NULL),
+ fPIDResponse(NULL),
+ fESD(NULL),
+ fIsMC(kFALSE),
+ fMCEvent(NULL),
+ fStack(NULL),
+ fESDTrackCuts(NULL),
+ fEtaMax(0.9),
+ fPtRange(NULL),
+ fNp(NULL),
+ fNCorrNp(2),
+ fCorrNp(NULL),
+ fNMCNp(5),
+ fMCNp(NULL),
+ fNControlMCNp(5),
+ fControlMCNp(NULL),
+
+ fHnTrackUnCorr(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliAnalysisNetParticleDistribution",10);
+}
+
+//________________________________________________________________________
+AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() {
+ // Destructor
+
+ if (fNp) delete[] fNp;
+
+ for (Int_t ii = 0; ii < fNCorrNp; ++ii)
+ if (fCorrNp[ii]) delete[] fCorrNp[ii];
+ if (fCorrNp) delete[] fCorrNp;
+
+ for (Int_t ii = 0; ii < fNMCNp; ++ii)
+ if (fMCNp[ii]) delete[] fMCNp[ii];
+ if (fMCNp) delete[] fMCNp;
+
+ for (Int_t ii = 0; ii < fNControlMCNp; ++ii)
+ if (fControlMCNp[ii]) delete[] fControlMCNp[ii];
+ if (fControlMCNp) delete[] fControlMCNp;
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax) {
+ // -- Initialize
+
+ fHelper = helper;
+ fESDTrackCuts = cuts;
+ fIsMC = isMC;
+ fPtRange = ptRange;
+ fEtaMax = etaMax;
+
+ // ------------------------------------------------------------------
+ // -- N particles / N anti-particles
+ // ------------------------------------------------------------------
+ // Np : arr[particle]
+ // CorrNp : arr[corrSet][particle]
+ // MCNp : arr[corrSet][particle] - MC
+ // ControlMCNp : arr[corrSet][particle] - Control MC
+
+ fNp = new Float_t[2];
+
+ fCorrNp = new Float_t*[fNCorrNp];
+ for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
+ fCorrNp[ii] = new Float_t[2];
+
+ fMCNp = new Float_t*[fNMCNp];
+ for (Int_t ii = 0 ; ii < fNMCNp; ++ii)
+ fMCNp[ii] = new Float_t[2];
+
+ fControlMCNp = new Float_t*[fNControlMCNp];
+ for (Int_t ii = 0 ; ii < fNControlMCNp; ++ii)
+ fControlMCNp[ii] = new Float_t[2];
+
+ ResetEvent();
+
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::CreateHistograms(TList* outList) {
+ // -- Add histograms to outlist
+
+ fOutList = outList;
+
+ // ------------------------------------------------------------------
+ // -- Create net particle histograms
+ // ------------------------------------------------------------------
+
+ AddHistSet("fHDist", "Uncorrected");
+ AddHistSet("fHDistCorr0", "Corrected [without cross section correction]");
+ AddHistSet("fHDistCorr1", "Corrected [with cross section correction]");
+
+ if (fIsMC) {
+ AddHistSet("fHMCrapidity", "MC primary in |y| < 0.5");
+ AddHistSet("fHMCptMin", "MC primary in |y| + #it{p}_{T} > 0.1");
+ AddHistSet("fHMCpt", Form("MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
+ AddHistSet("fHMCeta", Form("MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
+ AddHistSet("fHMCetapt", Form("MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
+
+ AddHistSet("fHControlMCLambdarapidity", "Control MC Lambda primary in |y| < 0.5");
+ AddHistSet("fHControlMCLambdaptMin", "Control MC Lambda primary in |y| + #it{p}_{T} > 0.1");
+ AddHistSet("fHControlMCLambdapt", Form("Control MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
+ AddHistSet("fHControlMCLambdaeta", Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
+ AddHistSet("fHControlMCLambdaetapt", Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
+ }
+
+ // ------------------------------------------------------------------
+ // -- Get Probe Particle Container
+ // ------------------------------------------------------------------
+
+ Double_t dCent[2] = {-0.5, 8.5};
+ Int_t iCent = 9;
+
+ Double_t dEta[2] = {-0.9, 0.9}; // -> 37 bins
+ Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
+
+ Double_t dRap[2] = {-0.5, 0.5};
+ Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
+
+ Double_t dPhi[2] = {0.0, TMath::TwoPi()};
+ Int_t iPhi = 42;
+
+ Double_t dPt[2] = {0.1, 3.0};
+ Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
+
+ Double_t dSign[2] = {-1.5, 1.5};
+ Int_t iSign = 3;
+
+ // cent: pt: sign: eta: phi: y
+ Int_t binShort[6] = { iCent, iPt, iSign, iEta, iPhi, iRap};
+ Double_t minShort[6] = {dCent[0], dPt[0], dSign[0], dEta[0], dPhi[0], dRap[0]};
+ Double_t maxShort[6] = {dCent[1], dPt[1], dSign[1], dEta[1], dPhi[1], dRap[1]};
+
+ // -- UnCorrected
+ fOutList->Add(new THnSparseF("fHnTrackUnCorr", "cent:pt:sign:eta:phi:y", 6, binShort, minShort, maxShort));
+ fHnTrackUnCorr = static_cast<THnSparseF*>(fOutList->Last());
+ fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality");
+ fHnTrackUnCorr->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ fHnTrackUnCorr->GetAxis(2)->SetTitle("sign");
+ fHnTrackUnCorr->GetAxis(3)->SetTitle("#eta");
+ fHnTrackUnCorr->GetAxis(4)->SetTitle("#varphi");
+ fHnTrackUnCorr->GetAxis(5)->SetTitle("#it{y}");
+
+ fHelper->BinLogAxis(fHnTrackUnCorr, 1);
+
+ // ------------------------------------------------------------------
+
+ return;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
+ // -- Setup Event
+
+ ResetEvent();
+
+ // -- Get ESD objects
+ fESDHandler = esdHandler;
+ fPIDResponse = esdHandler->GetPIDResponse();
+ fESD = fESDHandler->GetEvent();
+
+ // -- Get MC objects
+ fMCEvent = mcEvent;
+ if (fMCEvent)
+ fStack = fMCEvent->Stack();
+
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::ResetEvent() {
+ // -- Reset event
+
+ // -- Reset ESD Event
+ fESD = NULL;
+
+ // -- Reset MC Event
+ if (fIsMC)
+ fMCEvent = NULL;
+
+ // -- Reset N particles/anti-particles
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fNp[jj] = 0.;
+
+ // -- Reset N corrected particles/anti-particles
+ for (Int_t ii = 0; ii < fNCorrNp; ++ii)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fCorrNp[ii][jj] = 0.;
+
+ // -- Reset N MC particles/anti-particles
+ for (Int_t ii = 0; ii < fNMCNp; ++ii)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fMCNp[ii][jj] = 0.;
+
+ // -- Reset N control MC particles/anti-particles
+ for (Int_t ii = 0; ii < fNControlMCNp; ++ii)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fControlMCNp[ii][jj] = 0.;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::Process() {
+ // -- Process NetParticle Distributions
+
+ // -- Fill ESD tracks
+ ProcessESDTracks();
+
+ // -- Fill MC truth particles
+ if (fIsMC) {
+ ProcessStackParticles();
+ ProcessStackControlParticles();
+ }
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() {
+ // -- Process ESD tracks and fill histograms
+
+ for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Check track cuts
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ // -- Check if charged track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted
+ if (!fESDTrackCuts->AcceptTrack(track))
+ continue;
+
+ // -- Check if accepted in rapidity window
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP))
+ continue;
+
+ // -- Check if accepted bt PID from TPC or TPC+TOF
+ Double_t pid[2];
+ if (!fHelper->IsTrackAcceptedPID(track, pid))
+ continue;
+
+ // -- Check if accepted with thighter DCA cuts
+ if (fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill Probe Particle
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ Double_t aTrack[6] = {
+ Double_t(fHelper->GetCentralityBin()), // 0 centrality
+ track->Pt(), // 1 pt
+ track->GetSign(), // 2 sign
+ track->Eta(), // 3 eta
+ track->Phi(), // 4 phi
+ yP // 5 rapidity
+ };
+
+ fHnTrackUnCorr->Fill(aTrack);
+
+ // -- Count particle / anti-particle
+ // ------------------------------------------------------------------
+ // idxPart = 0 -> anti particle
+ // idxPart = 1 -> particle
+
+ Int_t idxPart = (track->GetSign() < 0) ? 0 : 1;
+ fNp[idxPart] += 1.;
+
+ for (Int_t ii = 0; ii < fNCorrNp; ++ii)
+ fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);
+
+ } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill Particle Fluctuation Histograms
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ // -- Uncorrected
+ FillHistSet("fHDist", fNp);
+
+ // -- Corrected
+ for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
+ FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]);
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() {
+ // -- Process primary particles from the stack and fill histograms
+
+ Int_t pdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies());
+
+ for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
+ TParticle* particle = fStack->Particle(idxMC);
+ if (!particle)
+ continue;
+
+ // -- Check basic MC properties -> charged physical primary
+ if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
+ continue;
+
+ // -- Check if particle / anti-particle
+ if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
+ continue;
+
+ // -- Get particle : 0 anti-particle / 1 particle
+ Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
+
+ // >> NOW only anti-particle / particle
+ // >> With idxPart
+
+ // -- Check rapidity window
+ Double_t yMC;
+ if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
+ continue;
+
+ fMCNp[0][idxPart] += 1.; // -> MCrapidity
+
+ // -- Check acceptance
+ if (particle->Pt() > 0.1 )
+ fMCNp[1][idxPart] += 1.; // -> MCptMin
+
+ if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
+ fMCNp[2][idxPart] += 1.; // -> MCpt
+
+ if (TMath::Abs(particle->Eta()) <= fEtaMax)
+ fMCNp[3][idxPart] += 1.; // -> MCeta
+
+ if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
+ fMCNp[4][idxPart] += 1.; // -> MCetapt
+
+ } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill Particle Fluctuation Histograms
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ FillHistSet("fHMCrapidity", fMCNp[0]);
+ FillHistSet("fHMCptMin", fMCNp[1]);
+ FillHistSet("fHMCpt", fMCNp[2]);
+ FillHistSet("fHMCeta", fMCNp[3]);
+ FillHistSet("fHMCetapt", fMCNp[4]);
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleDistribution::ProcessStackControlParticles() {
+ // -- Process primary particles from the stack and fill histograms
+
+ Int_t pdgCode = fHelper->GetControlParticleCode();
+ Bool_t isNeutral = fHelper->IsControlParticleNeutral();
+ const Char_t* name = fHelper->GetControlParticleName().Data();
+
+ for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
+ TParticle* particle = fStack->Particle(idxMC);
+ if (!particle)
+ continue;
+
+ // -- Check basic MC properties -> neutral or charged physical primary
+ if (isNeutral) {
+ if (!fHelper->IsParticleAcceptedBasicNeutral(particle, idxMC))
+ continue;
+ }
+ else {
+ if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
+ continue;
+ }
+
+ // -- Check if particle / anti-particle
+ if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
+ continue;
+
+ // -- Get particle : 0 anti-particle / 1 particle
+ Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
+
+ // >> NOW only anti-particle / particle
+ // >> With idxPart
+
+ // -- Check rapidity window
+ Double_t yMC;
+ if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
+ continue;
+
+ fControlMCNp[0][idxPart] += 1.; // -> ControlMCrapidity
+
+ // -- Check acceptance
+ if (particle->Pt() > 0.1 )
+ fControlMCNp[1][idxPart] += 1.; // -> ControlMCptMin
+
+ if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
+ fControlMCNp[2][idxPart] += 1.; // -> ControlMCpt
+
+ if (TMath::Abs(particle->Eta()) <= fEtaMax)
+ fControlMCNp[3][idxPart] += 1.; // -> ControlMCeta
+
+ if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
+ fControlMCNp[4][idxPart] += 1.; // -> ControlMCetapt
+ } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill Particle Fluctuation Histograms
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ FillHistSet(Form("fHControlMC%srapidity", name), fControlMCNp[0], 0);
+ FillHistSet(Form("fHControlMC%sptMin", name), fControlMCNp[1], 1);
+ FillHistSet(Form("fHControlMC%spt", name), fControlMCNp[2], 2);
+ FillHistSet(Form("fHControlMC%seta", name), fControlMCNp[3], 3);
+ FillHistSet(Form("fHControlMC%setapt", name), fControlMCNp[4], 4);
+
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::AddHistSet(const Char_t *name, const Char_t *title) {
+ // -- Add histogram sets for particle and anti-particle
+
+ const Char_t* aPartNames[] = {"pbar", "p"};
+ const Char_t* aPartTitles[] = {"Anti-Proton", "Proton"};
+ const Char_t* aPartAxisTitles[] = {"#bar{p}", "p"};
+
+ const Char_t* aControlPartNames[] = {"lambdabar", "lambda"};
+ const Char_t* aControlPartTitles[] = {"Anti-Lambda", "Lambda"};
+
+ fOutList->Add(new TList);
+ TList *list = static_cast<TList*>(fOutList->Last());
+ list->SetName(name) ;
+ list->SetOwner(kTRUE);
+
+ TString sName(name);
+
+ TString sPtTitle("");
+ if (!sName.Contains("fHMC") || !sName.Contains("fHControlMC"))
+ sPtTitle += Form("#it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]);
+
+ for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+ list->Add(new TH2F(Form("%s%s", name, aPartNames[idxPart]),
+ Form("%s %s Dist %s;Centrality;N_{%s}", title, aPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]),
+ 24, -0.5, 11.5, 501, -0.5, 500.49));
+ } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+
+ list->Add(new TH2F(Form("%sNet%s", name, aPartNames[1]),
+ Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5, 201, -100.5, 100.49));
+
+ list->Add(new TProfile(Form("%sNet%sM", name, aPartNames[1]),
+ Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+ list->Add(new TProfile(Form("%sNet%s2M", name, aPartNames[1]),
+ Form("%s (Net %s)^{2} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{2}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+ list->Add(new TProfile(Form("%sNet%s3M", name, aPartNames[1]),
+ Form("%s (Net %s)^{3} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{3}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+ list->Add(new TProfile(Form("%sNet%s4M", name, aPartNames[1]),
+ Form("%s (Net %s)^{4} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{4}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+ list->Add(new TProfile(Form("%sNet%s5M", name, aPartNames[1]),
+ Form("%s (Net %s)^{5} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{5}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+ list->Add(new TProfile(Form("%sNet%s6M", name, aPartNames[1]),
+ Form("%s (Net %s)^{6} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{6}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
+
+ list->Add(new TH2F(Form("%sNet%sOverSum", name, aPartNames[1]),
+ Form("%s (Net %s)/ Sum Dist %s;Centrality;(N_{p} - N_{#bar{p}})/(N_{p} + N_{#bar{p}})", title, aPartTitles[1], sPtTitle.Data()),
+ 24, -0.5, 11.5, 801, -50.5, 50.49));
+
+ if (sName.Contains("fHControlMC")) {
+ for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+ list->Add(new TH2F(Form("%s%sOver%s", name, aPartNames[idxPart], aControlPartNames[idxPart]),
+ Form("%s %s / %s Dist %s;Centrality;N_{%s}/N_{Tracks}",
+ title, aPartTitles[idxPart], aControlPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]),
+ 24, -0.5, 11.5, 101, 0., 1.));
+ } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+ }
+
+ // if (sName.Contains("fHDist")) {
+ // list->Add(new TH3F(Form("%sNet%s_RecVsMC", name, aPartNames[1]),
+ // Form("%s Net %s Dist - Rec Vs MC - %s;Centrality;(N_{p} - N_{#bar{p}})_{rec};(N_{p} - N_{#bar{p}})_{MC}",
+ // title, aPartTitles[1], sPtTitle.Data()),
+ // 24, -0.5, 11.5, 201, -100.5, 100.49, 201, -100.5, 100.49));
+ // }
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx) {
+ // -- Add histogram sets for particle and anti-particle
+
+ TList *list = static_cast<TList*>(fOutList->FindObject(name));
+
+ Float_t centralityBin = fHelper->GetCentralityBin();
+
+ (static_cast<TH2F*>(list->FindObject(Form("%spbar", name))))->Fill(centralityBin, np[0]);
+ (static_cast<TH2F*>(list->FindObject(Form("%sp", name))))->Fill(centralityBin, np[1]);
+
+ Float_t sumNp = np[1]+np[0];
+ Float_t deltaNp = np[1]-np[0];
+ Float_t deltaNp2 = deltaNp * deltaNp;
+ Float_t deltaNp3 = deltaNp2 * deltaNp;
+
+ (static_cast<TH2F*>(list->FindObject(Form("%sNetp", name))))->Fill(centralityBin, deltaNp);
+
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetpM", name))))->Fill(centralityBin, deltaNp);
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetp2M", name))))->Fill(centralityBin, deltaNp2);
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetp3M", name))))->Fill(centralityBin, deltaNp2*deltaNp);
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetp4M", name))))->Fill(centralityBin, deltaNp2*deltaNp2);
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetp5M", name))))->Fill(centralityBin, deltaNp3*deltaNp2);
+ (static_cast<TProfile*>(list->FindObject(Form("%sNetp6M", name))))->Fill(centralityBin, deltaNp3*deltaNp3);
+
+ (static_cast<TH2F*>(list->FindObject(Form("%sNetpOverSum", name))))->Fill(centralityBin, deltaNp/sumNp);
+
+ TString sName(name);
+ if (sName.Contains("fHControlMC") && controlIdx >= 0) {
+ (static_cast<TH2F*>(list->FindObject(Form("%spbarOverlambdabar", name))))->Fill(centralityBin, np[0]/fControlMCNp[controlIdx][0]);
+ (static_cast<TH2F*>(list->FindObject(Form("%spOverlambda", name))))->Fill(centralityBin, np[1]/fControlMCNp[controlIdx][1]);
+ }
+
+ // if (sName.Contains("fHDist")) {
+ // Float_t deltaMCNp = fMCNp[controlIdx][1]-fMCNp[controldx][0];
+ // (static_cast<TH3F*>(list->FindObject(Form("%sNetp_RecVsMC", name))))->Fill(centralityBin, deltaNp, deltaMCNp);
+ // }
+
+ return;
+}
+
--- /dev/null
+//-*- Mode: C++ -*-
+
+#ifndef ALIANALYSISNETPARTICLEDISTRIBUTION_H
+#define ALIANALYSISNETPARTICLEDISTRIBUTION_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Helper Class for for NetParticle Distributions
+// Authors: Jochen Thaeder <jochen@thaeder.de>
+
+#include "THnSparse.h"
+#include "TH1F.h"
+#include "TF1.h"
+
+#include "AliAnalysisNetParticleHelper.h"
+
+class AliESDtrack;
+class AliMCEvent;
+class AliStack;
+class AliPIDResponse;
+class AliESDInputHandler;
+class AliESDtrackCuts;
+
+class AliAnalysisNetParticleDistribution : public TNamed {
+
+ public:
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleDistribution();
+ virtual ~AliAnalysisNetParticleDistribution();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize */
+ Int_t Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax);
+
+ /** Add histograms to outlist */
+ void CreateHistograms(TList *outList);
+
+ /** Setup Event */
+ Int_t SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent);
+
+ /** Resre Event */
+ void ResetEvent();
+
+ /** Process NetParticle Distributions */
+ Int_t Process();
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliAnalysisNetParticleDistribution(const AliAnalysisNetParticleDistribution&); // not implemented
+ AliAnalysisNetParticleDistribution& operator=(const AliAnalysisNetParticleDistribution&); // not implemented
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Process - Private
+ * ---------------------------------------------------------------------------------
+ */
+ /** Process ESD tracks and fill histograms */
+ Int_t ProcessESDTracks();
+
+ /** Process primary particles from the stack and fill histograms */
+ Int_t ProcessStackParticles();
+
+ /** Process control particles from the stack and fill histograms */
+ Int_t ProcessStackControlParticles();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Helper Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Add set of histograms */
+ void AddHistSet(const Char_t *name, const Char_t *title);
+
+ /** Fill set of histograms */
+ void FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx = -1);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleHelper *fHelper; //! Ptr to helper class
+
+ TList *fOutList; //! Output data container
+ // -----------------------------------------------------------------------
+ AliESDInputHandler *fESDHandler; //! Ptr to ESD handler
+ AliPIDResponse *fPIDResponse; //! Ptr to PID response Object
+ AliESDEvent *fESD; //! Ptr to ESD event
+
+ Bool_t fIsMC; // Is MC event
+
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+ AliStack *fStack; //! Ptr to stack
+
+ AliESDtrackCuts *fESDTrackCuts; //! ESD cuts
+ // -----------------------------------------------------------------------
+ Float_t fEtaMax; // Max, absolut eta
+ Float_t *fPtRange; // Array of pt [min,max]
+ // -----------------------------------------------------------------------
+ Float_t *fNp; // Array of particle/anti-particle counts
+
+ Int_t fNCorrNp; // N sets of arrays of corrected particle/anti-particle counts
+ Float_t **fCorrNp; // Array of corrected particle/anti-particle counts
+
+ Int_t fNMCNp; // N sets of arrays of MC particle/anti-particle counts
+ Float_t **fMCNp; // Array of MC particle/anti-particle counts
+
+ Int_t fNControlMCNp; // N sets of arrays of control MC particle/anti-particle counts
+ Float_t **fControlMCNp; // Array of control MC particle/anti-particle counts
+ // -----------------------------------------------------------------------
+ THnSparseF *fHnTrackUnCorr; // THnSparseF : uncorrected probe particles
+ // -----------------------------------------------------------------------
+
+ ClassDef(AliAnalysisNetParticleDistribution, 1);
+};
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#include "TMath.h"
+#include "TAxis.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+
+#include "AliESDtrackCuts.h"
+
+#include "AliAnalysisNetParticleEffCont.h"
+
+using namespace std;
+
+// Task for NetParticle checks
+// Author: Jochen Thaeder <jochen@thaeder.de>
+
+ClassImp(AliAnalysisNetParticleEffCont)
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+AliAnalysisNetParticleEffCont::AliAnalysisNetParticleEffCont() :
+ fHelper(NULL),
+ fPdgCode(2212),
+
+ fESD(NULL),
+ fESDTrackCuts(NULL),
+
+ fCentralityBin(-1.),
+ fNTracks(0),
+
+ fStack(NULL),
+ fMCEvent(NULL),
+
+ fLabelsRec(NULL),
+
+ fHnEff(NULL),
+ fHnCont(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliAnalysisNetParticleEffCont",10);
+}
+
+//________________________________________________________________________
+AliAnalysisNetParticleEffCont::~AliAnalysisNetParticleEffCont() {
+ // Destructor
+
+ if (fLabelsRec[0])
+ delete[] (fLabelsRec[0]);
+ if (fLabelsRec[1])
+ delete[] (fLabelsRec[1]);
+ if (fLabelsRec)
+ delete[] fLabelsRec;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::Initialize(AliESDtrackCuts *cuts , AliAnalysisNetParticleHelper* helper) {
+
+ // -- ESD track cuts
+ // -------------------
+ fESDTrackCuts = cuts;
+
+ // -- Get Helper
+ // ---------------
+ fHelper = helper;
+
+ // -- Get particle species / pdgCode
+ // -------------------------
+ fPdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies());
+
+ // -- MC Labels for efficiency
+ // -----------------------------
+ fLabelsRec = new Int_t*[2];
+ for (Int_t ii = 0; ii < 2 ; ++ii)
+ fLabelsRec[ii] = NULL;
+
+ // -- Create THnSparse Histograms
+ // --------------------------------
+ CreateHistograms();
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleEffCont::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
+ // -- Setup event
+
+ // -- Reset of event
+ // -------------------
+ ResetEvent();
+
+ // -- Setup of event
+ // -------------------
+ fESD = esdHandler->GetEvent();
+ fNTracks = fESD->GetNumberOfTracks();
+
+ fMCEvent = mcEvent;
+ fStack = fMCEvent->Stack();
+
+ fCentralityBin = fHelper->GetCentralityBin();
+
+ // -- Create label arrays
+ // ------------------------
+ fLabelsRec[0] = new Int_t[fNTracks];
+ if(!fLabelsRec[0]) {
+ AliError("Cannot create fLabelsRec[0]");
+ return -1;
+ }
+
+ fLabelsRec[1] = new Int_t[fNTracks];
+ if(!fLabelsRec[1]) {
+ AliError("Cannot create fLabelsRec[1] for TPC");
+ return -1;
+ }
+
+ for(Int_t ii = 0; ii < fNTracks; ++ii) {
+ fLabelsRec[0][ii] = 0;
+ fLabelsRec[1][ii] = 0;
+ }
+
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::ResetEvent() {
+ // -- Reset event
+
+ for (Int_t ii = 0; ii < 2 ; ++ii) {
+ if (fLabelsRec[ii])
+ delete[] fLabelsRec[ii];
+ fLabelsRec[ii] = NULL;
+ }
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::Process() {
+ // -- Process event
+
+ // -- Setup (clean, create and fill) MC labels
+ // ---------------------------------------------
+ FillMCLabels();
+
+ // -- Fill MC histograms for efficiency studies
+ // -----------------------------------------------
+ FillMCEffHist();
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Private Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::CreateHistograms() {
+ // -- Create histograms
+
+ Double_t dCent[2] = {-0.5, 8.5};
+ Int_t iCent = 9;
+
+ Double_t dEta[2] = {-0.9, 0.9};
+ Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
+
+ Double_t dRap[2] = {-0.5, 0.5};
+ Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
+
+ Double_t dPhi[2] = {0.0, TMath::TwoPi()};
+ Int_t iPhi = 42;
+
+ Double_t dPt[2] = {0.1, 3.0};
+ Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
+
+ Double_t dSign[2] = {-1.5, 1.5};
+ Int_t iSign = 3;
+
+ // ------------------------------------------------------------------
+ // -- Create THnSparseF - Eff
+ // ------------------------------------------------------------------
+
+ // cent: etaMC: yMC: phiMC: ptMC: sign:findable:recStatus:pidStatus: etaRec: phiRec: ptRec:deltaEta: deltaPhi: deltaPt
+ Int_t binHnEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2 , 2 , iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1};
+ Double_t minHnEff[15] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], -0.5, -0.5, -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
+ Double_t maxHnEff[15] = {dCent[1], dEta[1], dRap[0], dPhi[1], dPt[1], dSign[1], 1.5, 1.5, 1.5, dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1]};
+
+ fHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
+ 15, binHnEff, minHnEff, maxHnEff);
+
+ fHnEff->Sumw2();
+
+ fHnEff->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
+ fHnEff->GetAxis(1)->SetTitle("#eta_{MC}"); // eta [-0.9, 0.9]
+ fHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
+ fHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. , 2Pi]
+ fHnEff->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pt [ 0.1, 3.0]
+
+ fHnEff->GetAxis(5)->SetTitle("sign"); // -1 | 0 | +1
+ fHnEff->GetAxis(6)->SetTitle("findable"); // 0 not findable | 1 findable
+ fHnEff->GetAxis(7)->SetTitle("recStatus"); // 0 not reconstructed | 1 reconstructed
+ fHnEff->GetAxis(8)->SetTitle("recPid"); // 0 not accepted | 1 accepted
+
+ fHnEff->GetAxis(9)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
+ fHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
+ fHnEff->GetAxis(11)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.1, 3.0]
+
+ fHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9,0.9]
+ fHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ 0. ,2Pi]
+ fHnEff->GetAxis(14)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 3.0,3.0]
+
+ fHelper->BinLogAxis(fHnEff, 4);
+ fHelper->BinLogAxis(fHnEff, 11);
+
+ /*
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+ creation -> findable -> reconstructed -> pid_TPC+TOF
+ (1) (2) (3) (4)
+ || findable | recStatus | recPid
+ 1) all primary probeParticles_MC || - - -
+ 2) all findable primary probeParticles_MC || x - -
+ 3) all reconstructed primary probeParticles_MC || x x -
+ 4) all reconstructed primary probeParticles_MC & recPid_TPC+TOF || x x x
+
+ <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ */
+
+ // ------------------------------------------------------------------
+ // -- Create THnSparseF - Cont
+ // ------------------------------------------------------------------
+
+ // cent: etaMC: yMC: phiMC: ptMC: sign:contPart: contSign: etaRec: phiRec: ptRec:deltaEta: deltaPhi: deltaPt
+ Int_t binHnCont[14] = { iCent, iEta, iRap, iPhi, iPt, iSign, 8, iSign, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1};
+ Double_t minHnCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], 0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
+ Double_t maxHnCont[14] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1]};
+ fHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",
+ 14, binHnCont, minHnCont, maxHnCont);
+
+ fHnCont->Sumw2();
+
+ fHnCont->GetAxis(0)->SetTitle("centrality"); // 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
+ fHnCont->GetAxis(1)->SetTitle("#eta_{MC}"); // eta [-0.9,0.9]
+ fHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}"); // rapidity [-0.5, 0.5]
+ fHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); // phi [ 0. ,2Pi]
+ fHnCont->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})"); // pT [ 0.1,1.3]
+ fHnCont->GetAxis(5)->SetTitle("sign"); // -1 | 0 | +1
+
+ fHnCont->GetAxis(6)->SetTitle("contPart"); // 1 pi | 2 K | 3 p | 4 e | 5 mu | 6 other | 7 p from WeakDecay | 8 p from Material
+ fHnCont->GetAxis(7)->SetTitle("contSign"); // -1 | 0 | +1
+
+ fHnCont->GetAxis(8)->SetTitle("#eta_{Rec}"); // eta [-0.9, 0.9]
+ fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
+ fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.1, 3.0]
+
+ fHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9,0.9]
+ fHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ 0. ,2Pi]
+ fHnCont->GetAxis(13)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 3.0,3.0]
+
+ fHelper->BinLogAxis(fHnCont, 4);
+ fHelper->BinLogAxis(fHnCont, 10);
+
+ // ------------------------------------------------------------------
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::FillMCLabels() {
+ // Fill MC labels
+ // Loop over ESD tracks and fill arrays with MC lables
+ // fLabelsRec[0] : all Tracks
+ // fLabelsRec[1] : all Tracks accepted by PID of TPC
+ // Check every accepted track if correctly identified
+ // otherwise check for contamination
+
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted
+ if (!fESDTrackCuts->AcceptTrack(track))
+ continue;
+
+ // -- Check if accepted in rapidity window
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP))
+ continue;
+
+ // -- Check if accepted with thighter DCA cuts
+ if (!fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ Int_t label = TMath::Abs(track->GetLabel());
+
+ // -- Fill Label of all reconstructed
+ fLabelsRec[0][idxTrack] = label;
+
+ // -- Check if accepted by PID from TPC or TPC+TOF
+ Double_t pid[2];
+ if (!fHelper->IsTrackAcceptedPID(track, pid))
+ continue;
+
+ // -- Fill Label of all reconstructed && recPid_TPC+TOF
+ fLabelsRec[1][idxTrack] = label;
+
+ // -- Check for contamination and fill contamination THnSparse
+ CheckContTrack(label, track->GetSign(), idxTrack);
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack) {
+ // Check if particle is contamination or correctly identified
+ // Fill contamination THnSparse
+
+ TParticle* particle = fStack->Particle(label);
+ if (!particle)
+ return;
+
+ Bool_t isSecondaryFromWeakDecay = kFALSE;
+ Bool_t isSecondaryFromMaterial = kFALSE;
+
+ // -- Check if correctly identified
+ if (particle->GetPdgCode() == (sign*fPdgCode)) {
+
+ // Check if is physical primary -> all ok
+ if (fStack->IsPhysicalPrimary(label))
+ return;
+
+ // -- Check if secondaries from material or weak decay
+ isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label);
+ isSecondaryFromMaterial = fStack->IsSecondaryFromMaterial(label);
+ }
+
+ // -- Get MC pdg
+ Float_t contSign = 0.;
+ if (particle->GetPDG()->Charge() == 0.) contSign = 0.;
+ else if (particle->GetPDG()->Charge() < 0.) contSign = -1.;
+ else if (particle->GetPDG()->Charge() > 0.) contSign = 1.;
+
+ // -- Get contaminating particle
+ Float_t contPart = 0;
+ if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
+ else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
+ else {
+ if (TMath::Abs(particle->GetPdgCode()) == 211) contPart = 1; // pion
+ else if (TMath::Abs(particle->GetPdgCode()) == 321) contPart = 2; // kaon
+ else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
+ else if (TMath::Abs(particle->GetPdgCode()) == 11) contPart = 4; // electron
+ else if (TMath::Abs(particle->GetPdgCode()) == 13) contPart = 5; // muon
+ else contPart = 6; // other
+ }
+
+ // -- Get Reconstructed values
+ Float_t etaRec = 0.;
+ Float_t phiRec = 0.;
+ Float_t ptRec = 0.;
+
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+ if (track) {
+ // if no track present (which should not happen)
+ // -> pt = 0. , which is not in the looked at range
+
+ // -- Get Reconstructed eta/phi/pt
+ etaRec = track->Eta();
+ phiRec = track->Phi();
+ ptRec = track->Pt();
+ }
+
+ // -- Fill THnSparse
+ Double_t hnCont[14] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
+ contPart, contSign, etaRec, phiRec, ptRec,
+ particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
+ fHnCont->Fill(hnCont);
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleEffCont::FillMCEffHist() {
+ // Fill efficiency THnSparse
+
+ Int_t nPart = fStack->GetNprimary();
+
+ for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+ TParticle* particle = fStack->Particle(idxMC);
+
+ // -- Check basic MC properties -> charged physical primary
+ if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
+ continue;
+
+ // -- Check rapidity window
+ Double_t yMC;
+ if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
+ continue;
+
+ // -- Check if probeParticle / anti-probeParticle
+ if (TMath::Abs(particle->GetPdgCode()) != fPdgCode)
+ continue;
+
+ // -- Get sign of particle
+ Float_t sign = (particle->GetPdgCode() < 0) ? -1. : 1.;
+
+ // -- Get if particle is findable
+ Float_t findable = Float_t(fHelper->IsParticleFindable(idxMC));
+
+ // -- Get recStatus and pidStatus
+ Float_t recStatus = 0.;
+ Float_t recPid = 0.;
+
+ // -- Get Reconstructed values
+ Float_t etaRec = 0.;
+ Float_t phiRec = 0.;
+ Float_t ptRec = 0.;
+
+ // -- Loop over all labels
+ for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
+ if (idxMC == fLabelsRec[0][idxRec]) {
+ recStatus = 1.;
+
+ if (idxMC == fLabelsRec[1][idxRec])
+ recPid = 1.;
+
+ AliESDtrack *track = fESD->GetTrack(idxRec);
+ if (track) {
+ // if no track present (which should not happen)
+ // -> pt = 0. , which is not in the looked at range
+
+ // -- Get Reconstructed eta/phi/pt
+ etaRec = track->Eta();
+ phiRec = track->Phi();
+ ptRec = track->Pt();
+ }
+ break;
+ }
+ } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
+
+ // -- Fill THnSparse
+ Double_t hnEff[15] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign,
+ findable, recStatus, recPid, etaRec, phiRec, ptRec,
+ particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec};
+ fHnEff->Fill(hnEff);
+
+ } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+
+ return;
+}
--- /dev/null
+//-*- Mode: C++ -*-
+
+#ifndef ALIANALYSISNETPARTICLEEFFCONT_H
+#define ALIANALYSISNETPARTICLEEFFCONT_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Efficiency and Contaminations for NetParticle Distributions
+// Authors: Jochen Thaeder <jochen@thaeder.de>
+
+#include "THnSparse.h"
+
+#include "AliAnalysisNetParticleHelper.h"
+
+class AliESDEvent;
+class AliESDInputHandler;
+class AliMCEvent;
+
+class AliAnalysisNetParticleEffCont : public TNamed {
+
+ public:
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleEffCont();
+ virtual ~AliAnalysisNetParticleEffCont();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize */
+ void Initialize(AliESDtrackCuts *cuts, AliAnalysisNetParticleHelper* helper);
+
+ /** Setup Event */
+ Int_t SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent);
+
+ /** Reset Event */
+ void ResetEvent();
+
+ /** Process Event */
+ void Process();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Getter
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Get Ptr to efficiency THnSparse */
+ THnSparseF* GetHnEff() {return fHnEff;}
+
+ /** Get Ptr to contaminiation THnSparse */
+ THnSparseF* GetHnCont() {return fHnCont;}
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliAnalysisNetParticleEffCont(const AliAnalysisNetParticleEffCont&); // not implemented
+ AliAnalysisNetParticleEffCont& operator=(const AliAnalysisNetParticleEffCont&); // not implemented
+
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Create the efficiency / contamination THnSparseF */
+ void CreateHistograms();
+
+ /** Fill MC labels */
+ void FillMCLabels();
+
+ /** Fill efficiency THnSparse */
+ void FillMCEffHist();
+
+ /** Check if particle is contamination */
+ void CheckContTrack(Int_t label, Float_t sign, Int_t idxTrack);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleHelper *fHelper; //! Ptr to helper class
+
+ // -----------------------------------------------------------------------
+
+ Int_t fPdgCode; // PDG code of particle to be found
+
+ // --- ESD only ----------------------------------------------------------
+
+ AliESDEvent *fESD; //! ESD object
+ AliESDtrackCuts *fESDTrackCuts; //! ESD cuts
+
+ // -----------------------------------------------------------------------
+
+ Float_t fCentralityBin; // Centrality of current event
+ Int_t fNTracks; // N Tracks in the current event
+
+ // --- MC only -----------------------------------------------------------
+
+ AliStack *fStack; //! Ptr to stack
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+
+ Int_t **fLabelsRec; //! 2x nTracks large array with labels for MC particles
+
+ // -----------------------------------------------------------------------
+
+ THnSparseF *fHnEff; // THnSparseF efficiency
+ THnSparseF *fHnCont; // THnSparseF contamination
+
+ // -----------------------------------------------------------------------
+
+ ClassDef(AliAnalysisNetParticleEffCont, 1);
+};
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#include "TMath.h"
+#include "TAxis.h"
+#include "TSystem.h"
+#include "TFile.h"
+#include "TPRegexp.h"
+
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDInputHandler.h"
+#include "AliESDpid.h"
+#include "AliCentrality.h"
+#include "AliTracker.h"
+
+#include "AliAnalysisNetParticleHelper.h"
+
+using namespace std;
+
+// Task for NetParticle checks
+// Author: Jochen Thaeder <jochen@thaeder.de>
+
+ClassImp(AliAnalysisNetParticleHelper)
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+AliAnalysisNetParticleHelper::AliAnalysisNetParticleHelper() :
+ fESDHandler(NULL),
+ fPIDResponse(NULL),
+ fESD(NULL),
+ fMCEvent(NULL),
+ fStack(NULL),
+
+ fCentralityBin(-1),
+ fCentralityPercentile(-1.),
+
+ fCentralityBinMax(7),
+ fVertexZMax(10.),
+ fRapidityMax(0.5),
+ fMinTrackLengthMC(70.),
+ fNSigmaMaxCdd(3.),
+ fNSigmaMaxCzz(3.),
+
+ fParticleSpecies(AliPID::kProton),
+ fControlParticleCode(3122),
+ fControlParticleIsNeutral(kTRUE),
+ fControlParticleName("Lambda"),
+
+ fNSigmaMaxTPC(2.5),
+ fNSigmaMaxTOF(2.5),
+ fMinPtForTOFRequired(0.8),
+
+ fHEventStat0(NULL),
+ fHEventStat1(NULL),
+ fHEventStatMax(6),
+
+ fHTriggerStat(NULL),
+ fNTriggers(5),
+
+ fHCentralityStat(NULL),
+ fNCentralityBins(10),
+
+ fEtaCorrFunc(NULL),
+ fCorr0(NULL),
+ fCorr1(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
+}
+
+//________________________________________________________________________
+AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
+ // Destructor
+
+ for (Int_t jj = 0; jj < 2; ++jj) {
+ if (fCorr0[jj]) delete[] fCorr0[jj];
+ if (fCorr1[jj]) delete[] fCorr1[jj];
+ }
+ if (fCorr0) delete[] fCorr0;
+ if (fCorr1) delete[] fCorr1;
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleHelper::Initialize(Bool_t isMC) {
+ // -- Initialize helper
+
+ Int_t iResult = 0;
+
+ // -- Setup event cut statistics
+ InitializeEventStats();
+
+ // -- Setup trigger statistics
+ InitializeTriggerStats();
+
+ // -- Setup centrality statistics
+ InitializeCentralityStats();
+
+ // -- Load Eta correction function
+ iResult = InitializeEtaCorrection(isMC);
+
+ // -- Load Eta correction function
+ iResult = InitializeTrackbyTrackCorrection();
+
+ return iResult;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
+ // -- Setup Event
+
+ // -- Get ESD objects
+ fESDHandler = esdHandler;
+ fPIDResponse = esdHandler->GetPIDResponse();
+ fESD = fESDHandler->GetEvent();
+
+ // -- Get MC objects
+ fMCEvent = mcEvent;
+ if (fMCEvent)
+ fStack = fMCEvent->Stack();
+
+ // -- Get event centrality
+ // > 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
+ // > 0 1 2 3 4 5 6 7 8 9
+
+ AliCentrality *centrality = fESD->GetCentrality();
+
+ Int_t centBin = centrality->GetCentralityClass10("V0M");
+ if (centBin == 0)
+ fCentralityBin = centrality->GetCentralityClass5("V0M");
+ else if (centBin == 10 || centBin == -1.)
+ fCentralityBin = -1;
+ else if (centBin > 0 && centBin < fNCentralityBins)
+ fCentralityBin = centBin + 1;
+ else
+ fCentralityBin = -2;
+
+ // -- Stay within the max centrality bin
+ if (fCentralityBin >= fCentralityBinMax)
+ fCentralityBin = -2;
+
+ fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
+
+ // -- Update TPC pid with eta correction
+ UpdateEtaCorrectedTPCPid();
+
+ return 0;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Event / Trigger Statistics
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
+ // -- Check if Event is triggered and fill Trigger Histogram
+
+ Bool_t *aTriggerFired = new Bool_t[fNTriggers];
+ for (Int_t ii = 0; ii < fNTriggers; ++ii)
+ aTriggerFired[ii] = kFALSE;
+
+ if ((fESDHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
+ if ((fESDHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
+ if ((fESDHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
+ if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
+ if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEGA)) aTriggerFired[4] = kTRUE;
+
+ Bool_t isTriggered = kFALSE;
+
+ for (Int_t ii=0; ii<fNTriggers; ++ii) {
+ if(aTriggerFired[ii]) {
+ isTriggered = kTRUE;
+ fHTriggerStat->Fill(ii);
+ }
+ }
+
+ delete[] aTriggerFired;
+
+ return isTriggered;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
+ // -- Evaluate event statistics histograms
+
+ Int_t *aEventCuts = new Int_t[fHEventStatMax];
+ // set aEventCuts[ii] to 1 in case of reject
+
+ for (Int_t ii=0;ii<fHEventStatMax; ++ii)
+ aEventCuts[ii] = 0;
+
+ Int_t iCut = 0;
+
+ // -- 0 - Before Physics Selection
+ aEventCuts[iCut] = 0;
+
+ // -- 1 - No Trigger fired
+ ++iCut;
+ if (!IsEventTriggered())
+ aEventCuts[iCut] = 1;
+
+ // -- 2 - No Vertex
+ ++iCut;
+ const AliESDVertex* vtxESD = fESD->GetPrimaryVertexTracks();
+ if (!vtxESD)
+ aEventCuts[iCut] = 1;
+
+ // -- 3 - Vertex z outside cut window
+ ++iCut;
+ if (vtxESD){
+ if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
+ aEventCuts[iCut] = 1;
+ }
+ else
+ aEventCuts[iCut] = 1;
+
+ // -- 4 - Centrality = -1 (no centrality or not hadronic)
+ ++iCut;
+ if(fCentralityBin == -1.)
+ aEventCuts[iCut] = 1;
+
+ // -- 5 - Centrality < fCentralityMax
+ ++iCut;
+ if(fCentralityBin == -2.)
+ aEventCuts[iCut] = 1;
+
+ // -- Fill statistics / reject event
+ Bool_t isRejected = FillEventStats(aEventCuts);
+
+ // -- Cleanup
+ delete[] aEventCuts;
+
+ return isRejected;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Accept Particle Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(TParticle *particle, Int_t idxMC) {
+ // -- Check if MC particle is accepted for basic parameters
+
+ if (!particle)
+ return kFALSE;
+
+ // -- check if PDF code exists
+ if (!particle->GetPDG())
+ return kFALSE;
+
+ // -- check if charged
+ if (particle->GetPDG()->Charge() == 0.0)
+ return kFALSE;
+
+ // -- check if physical primary
+ if(!fStack->IsPhysicalPrimary(idxMC))
+ return kFALSE;
+
+ return kTRUE;
+}
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(TParticle *particle, Int_t idxMC) {
+ // -- Check if MC particle is accepted for basic parameters
+
+ if (!particle)
+ return kFALSE;
+
+ // -- check if PDF code exists
+ if (!particle->GetPDG())
+ return kFALSE;
+
+ // -- check if neutral
+ if (particle->GetPDG()->Charge() != 0.0)
+ return kFALSE;
+
+ // -- check if physical primary
+ if(!fStack->IsPhysicalPrimary(idxMC))
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(TParticle *particle, Double_t &yP) {
+ // -- Check if particle is accepted
+ // > in rapidity
+ // > return 0 if not accepted
+
+ Double_t mP = AliPID::ParticleMass(fParticleSpecies);
+
+ // -- Calculate rapidities and kinematics
+ Double_t p = particle->P();
+ Double_t pz = particle->Pz();
+
+ Double_t eP = TMath::Sqrt(p*p + mP*mP);
+ yP = 0.5 * TMath::Log((eP + pz) / (eP - pz));
+
+ // -- Check Rapidity window
+ if (TMath::Abs(yP) > fRapidityMax)
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
+ // -- Check if MC particle is findable tracks
+
+ AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
+ if(!mcParticle)
+ return kFALSE;
+
+ Int_t counter;
+ Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
+
+ return (tpcTrackLength > fMinTrackLengthMC);
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Accept Track Methods - public
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliESDtrack *track) {
+ // -- Check if track is accepted
+ // > for basic parameters
+
+ if (!track)
+ return kFALSE;
+
+ if (track->Charge() == 0)
+ return kFALSE;
+
+ // -- Get momentum for dEdx
+ if (!track->GetInnerParam())
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliESDtrack *track, Double_t &yP) {
+ // -- Check if track is accepted
+ // > in rapidity
+ // > return 0 if not accepted
+
+ Double_t mP = AliPID::ParticleMass(fParticleSpecies);
+
+ // -- Calculate rapidities and kinematics
+ Double_t pvec[3];
+ track->GetPxPyPz(pvec);
+
+ Double_t p = track->GetP();
+ Double_t eP = TMath::Sqrt(p*p + mP*mP);
+ yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
+
+ // -- Check Rapidity window
+ if (TMath::Abs(yP) > fRapidityMax)
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliESDtrack *track) {
+ // -- Check if track is accepted
+ // > for DCA, if both SPD layers have hits
+
+ Bool_t isAccepted = kTRUE;
+
+ // -- Get nHits SPD
+ if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
+
+ // -- Get DCA nSigmas
+ Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
+ track->GetImpactParameters(dca,cov);
+
+ Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
+ Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
+
+ if (fNSigmaMaxCdd != 0.) {
+ if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
+ isAccepted = kFALSE;
+ }
+
+ if (fNSigmaMaxCzz != 0.) {
+ if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
+ isAccepted = kFALSE;
+ }
+ }
+
+ return isAccepted;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliESDtrack *track, Double_t* pid) {
+ // -- Check if track is accepted
+ // > provides TPC and TOF nSigmas to the argument
+
+ Bool_t isAcceptedTPC = kFALSE;
+ Bool_t isAcceptedTOF = kTRUE;
+
+ // -- Get PID
+ pid[0] = fPIDResponse->NumberOfSigmasTPC(track, fParticleSpecies);
+ pid[1] = fPIDResponse->NumberOfSigmasTOF(track, fParticleSpecies);
+
+ // -- Check PID with TPC
+ if (TMath::Abs(pid[0]) < fNSigmaMaxTPC)
+ isAcceptedTPC = kTRUE;
+
+ // -- Check PID with TOF
+ if (isAcceptedTPC) {
+
+ // Has TOF PID availible
+ if (track->GetStatus() & AliVTrack::kTOFpid) {
+ if (TMath::Abs(pid[1]) < fNSigmaMaxTOF)
+ isAcceptedTOF = kTRUE;
+ else
+ isAcceptedTOF = kFALSE;
+ }
+
+ // Has no TOF PID availible
+ else {
+ if (track->Pt() > fMinPtForTOFRequired)
+ isAcceptedTOF = kFALSE;
+ else
+ isAcceptedTOF = kTRUE;
+ }
+ } // if (isAcceptedTPC) {
+
+ return (isAcceptedTPC && isAcceptedTOF);
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Helper Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisNetParticleHelper::UpdateEtaCorrectedTPCPid() {
+ // -- Update eta corrected TPC pid
+
+ if (!fEtaCorrFunc)
+ return;
+
+ for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+
+ // -- Check if charged track is accepted for basic parameters
+ if (!IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ Double_t newTPCSignal = track->GetTPCsignal() / fEtaCorrFunc->Eval(track->Eta());
+ track->SetTPCsignal(newTPCSignal, track->GetTPCsignalSigma(), track->GetTPCsignalN());
+ }
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisNetParticleHelper::GetTrackbyTrackCorrectionFactor(Double_t *aTrack, Int_t flag) {
+ // -- Get efficiency correctionf of particle dependent on (eta, phi, pt, centrality)
+
+ Int_t idxPart = (aTrack[2] < 0) ? 0 : 1;
+ THnSparseF* corrMatrix = (flag == 0) ? fCorr0[idxPart][fCentralityBin] : fCorr1[idxPart][fCentralityBin];
+
+ Double_t dimBin[3] = {aTrack[3], aTrack[4], aTrack[1]}; // eta, phi, pt
+
+ Double_t corr = corrMatrix->GetBinContent(corrMatrix->GetBin(dimBin));
+ if (corr == 0.) {
+ AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",
+ aTrack[3], aTrack[4], aTrack[1], fCentralityBin));
+ corr = 1.;
+ }
+
+ return corr;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleHelper::BinLogAxis(const THnSparseF *h, Int_t axisNumber) {
+ // -- Method for the correct logarithmic binning of histograms
+
+ TAxis *axis = h->GetAxis(axisNumber);
+ Int_t bins = axis->GetNbins();
+
+ Double_t from = axis->GetXmin();
+ Double_t to = axis->GetXmax();
+ Double_t *newBins = new Double_t[bins + 1];
+
+ newBins[0] = from;
+ Double_t factor = pow(to/from, 1./bins);
+
+ for (int i = 1; i <= bins; i++) {
+ newBins[i] = factor * newBins[i-1];
+ }
+ axis->Set(bins, newBins);
+ delete [] newBins;
+}
+
+///////////////////////////////////////////////////////////////////////////////////
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Initialize - Private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisNetParticleHelper::InitializeEventStats() {
+ // -- Initialize event statistics histograms
+
+ const Char_t* aEventNames[] = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
+ const Char_t* aCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
+
+ fHEventStat0 = new TH1F("fHEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
+ fHEventStat1 = new TH1F("fHEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
+
+ for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
+ fHEventStat0->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
+ fHEventStat1->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
+ }
+ fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
+ fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
+ // -- Initialize trigger statistics histograms
+
+ const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
+
+ fHTriggerStat = new TH1F("fHTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
+
+ for ( Int_t ii=0; ii < fNTriggers; ii++ )
+ fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
+ // -- Initialize trigger statistics histograms
+
+ const Char_t* aCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
+ "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
+
+ fHCentralityStat = new TH1F("fHCentralityStat","Centrality statistics;Centrality Bins;Events",
+ fNCentralityBins,-0.5,fNCentralityBins-0.5);
+
+ for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
+ fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, aCentralityNames[ii]);
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleHelper::InitializeEtaCorrection(Bool_t isMC) {
+ // -- Initialize eta correction maps for TPC pid
+
+ TFile fileEtaCorrMaps("${ALICE_ROOT}/PWGDQ/dielectron/files/EtaCorrMaps.root");
+ if (!fileEtaCorrMaps.IsOpen())
+ return -1;
+
+ TList *keyList = fileEtaCorrMaps.GetListOfKeys();
+
+ TString sList("");
+ sList = (isMC) ? "LHC11a10" : "LHC10h.pass2";
+
+ for (Int_t idx = 0; idx < keyList->GetEntries(); ++idx) {
+ TString keyName = keyList->At(idx)->GetName();
+ TPRegexp reg(keyName);
+ if (reg.MatchB(sList)){
+ AliInfo(Form("Using Eta Correction Function: %s",keyName.Data()));
+ fEtaCorrFunc = static_cast<TF1*>(fileEtaCorrMaps.Get(keyName.Data()));
+ return 0;
+ }
+ }
+
+ return -2;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisNetParticleHelper::InitializeTrackbyTrackCorrection() {
+ // -- Initialize track by track correction matrices
+
+ AliInfo("TODO ... get the correct name from particle"); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ TFile* corrFile = TFile::Open("${ALICE_ROOT}/PWGCF/EBYE/NetParticle/eff/effectiveCorrectionProton.root");
+ if (!corrFile) {
+ AliError("Track-by-track correction file can not be opened!");
+ return -1;
+ }
+
+ // -- correction - not cross section corrected
+ fCorr0 = new THnSparseF**[2];
+ fCorr0[0] = new THnSparseF*[fCentralityBinMax];
+ fCorr0[1] = new THnSparseF*[fCentralityBinMax];
+
+ for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
+ THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr0_Cent_%d", idxCent)));
+ THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr0_Cent_%d", idxCent)));
+
+ if (!sp0 || !sp1) {
+ AliError(Form("Effective correction objects 0 - idxCent %d can not be retrieved!",idxCent));
+ return -1;
+ }
+
+ fCorr0[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
+ fCorr0[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
+ }
+
+ // -- correction - ross section corrected
+ fCorr1 = new THnSparseF**[2];
+ fCorr1[0] = new THnSparseF*[fCentralityBinMax];
+ fCorr1[1] = new THnSparseF*[fCentralityBinMax];
+
+ for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
+ THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr1_Cent_%d", idxCent)));
+ THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr1_Cent_%d", idxCent)));
+
+ if (!sp0 || !sp1) {
+ AliError(Form("Effective correction objects 1 - idxCent %d can not be retrieved!",idxCent));
+ return -1;
+ }
+
+ fCorr1[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
+ fCorr1[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
+ }
+
+ return 0;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Event / Trigger Statistics - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
+ // -- Fill event / centrality statistics
+
+ Bool_t isRejected = kFALSE;
+
+ // -- Fill event statistics
+ for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
+ if (aEventCuts[idx])
+ isRejected = kTRUE;
+ else
+ fHEventStat0->Fill(idx);
+ }
+
+ for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
+ if (aEventCuts[idx])
+ break;
+ fHEventStat1->Fill(idx);
+ }
+
+ // -- Fill centrality statistics of accepted events
+ if (!isRejected)
+ fHCentralityStat->Fill(fCentralityBin);
+
+ return isRejected;
+}
--- /dev/null
+//-*- Mode: C++ -*-
+
+#ifndef ALIANALYSISNETPARTICLEHELPER_H
+#define ALIANALYSISNETPARTICLEHELPER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Helper Class for for NetParticle Distributions
+// Authors: Jochen Thaeder <jochen@thaeder.de>
+
+#include "THnSparse.h"
+#include "TParticle.h"
+#include "TH1F.h"
+#include "TF1.h"
+
+class AliESDtrack;
+class AliMCEvent;
+class AliStack;
+class AliPIDResponse;
+class AliESDInputHandler;
+
+class AliAnalysisNetParticleHelper : public TNamed {
+
+ public:
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleHelper();
+ virtual ~AliAnalysisNetParticleHelper();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Setter
+ * ---------------------------------------------------------------------------------
+ */
+
+ void SetCentralityBinMax(Int_t d) {fCentralityBinMax = d;}
+ void SetVertexZMax(Float_t f) {fVertexZMax = f;}
+ void SetRapidityMax(Float_t f) {fRapidityMax = f;}
+ void SetMinTrackLengthMC(Float_t f) {fMinTrackLengthMC = f;}
+ void SetNSigmaMaxCdd(Float_t f) {fNSigmaMaxCdd = f;}
+ void SetNSigmaMaxCzz(Float_t f) {fNSigmaMaxCzz = f;}
+
+ void SetParticleSpecies(AliPID::EParticleType pid) {fParticleSpecies = pid;}
+ void SetControlParticleSpecies(Int_t pdgCode, Bool_t isNeutral, TString name) {
+ fControlParticleCode = pdgCode;
+ fControlParticleIsNeutral = isNeutral;
+ fControlParticleName = name;
+ }
+
+ void SetNSigmaMaxTPC(Float_t f) {fNSigmaMaxTPC = f;}
+ void SetNSigmaMaxTOF(Float_t f) {fNSigmaMaxTOF = f;}
+ void SetMinPtForTOFRequired(Float_t f) {fMinPtForTOFRequired = f;}
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Getter
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliPID::EParticleType GetParticleSpecies(){return fParticleSpecies;}
+
+ TH1F* GetHEventStat0() {return fHEventStat0;}
+ TH1F* GetHEventStat1() {return fHEventStat1;}
+ TH1F* GetHTriggerStat() {return fHTriggerStat;}
+ TH1F* GetHCentralityStat() {return fHCentralityStat;}
+
+ Int_t GetCentralityBin() {return fCentralityBin;}
+ Float_t GetCentralityPercentile() {return fCentralityPercentile;}
+
+ Int_t GetControlParticleCode() {return fControlParticleCode;}
+ Bool_t IsControlParticleNeutral() {return fControlParticleIsNeutral;}
+ TString& GetControlParticleName() {return fControlParticleName;}
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize Helper */
+ Int_t Initialize(Bool_t isMC);
+
+ /** Setup Event */
+ Int_t SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Event / Trigger Statistics
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Check if event is triggred */
+ Bool_t IsEventTriggered();
+
+ /** Fill event cut statistics */
+ Bool_t IsEventRejected();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Accept Particle Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Check if charged MC particle is accepted for basic parameters */
+ Bool_t IsParticleAcceptedBasicCharged(TParticle *particle, Int_t idxMC);
+
+ /** Check if neutral MC particle is accepted for basic parameters */
+ Bool_t IsParticleAcceptedBasicNeutral(TParticle *particle, Int_t idxMC);
+
+ /** Check if MC particle is accepted for Rapidity */
+ Bool_t IsParticleAcceptedRapidity(TParticle *particle, Double_t &yP);
+
+ /** Check if MC particle is findable tracks */
+ Bool_t IsParticleFindable(Int_t label);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Accept Track Methods - public
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Check if track is accepted for basic parameters */
+ Bool_t IsTrackAcceptedBasicCharged(AliESDtrack *track);
+
+ /** Check if track is accepted for Rapidity */
+ Bool_t IsTrackAcceptedRapidity(AliESDtrack *track, Double_t &yP);
+
+ /** Check if track is accepted for DCA */
+ Bool_t IsTrackAcceptedDCA(AliESDtrack *track);
+
+ /** Check if track is accepted for PID */
+ Bool_t IsTrackAcceptedPID(AliESDtrack *track, Double_t *pid);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Helper Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Update eta corrected TPC pid */
+ void UpdateEtaCorrectedTPCPid();
+
+ /** Get efficiency correctionf of particle dependent on (eta, phi, pt, centrality) */
+ Double_t GetTrackbyTrackCorrectionFactor(Double_t *aTrack, Int_t flag);
+
+ /** Method for the correct logarithmic binning of histograms */
+ void BinLogAxis(const THnSparseF *h, Int_t axisNumber);
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliAnalysisNetParticleHelper(const AliAnalysisNetParticleHelper&); // not implemented
+ AliAnalysisNetParticleHelper& operator=(const AliAnalysisNetParticleHelper&); // not implemented
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Initialize - Private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize event cut statistics */
+ void InitializeEventStats();
+
+ /** Initialize trigger statistics */
+ void InitializeTriggerStats();
+
+ /** Initialize centrality statistics */
+ void InitializeCentralityStats();
+
+ /** Initialize eta correction maps for TPC pid */
+ Int_t InitializeEtaCorrection(Bool_t isMC);
+
+ /** Initialize track by track correction matrices */
+ Int_t InitializeTrackbyTrackCorrection();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Event / Trigger Statistics - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Fill event cut statistics */
+ Bool_t FillEventStats(Int_t *aEventCuts);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliESDInputHandler *fESDHandler; //! Ptr to ESD handler
+ AliPIDResponse *fPIDResponse; //! Ptr to PID response Object
+ AliESDEvent *fESD; //! Ptr to ESD event
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+ AliStack *fStack; //! Ptr to stack
+
+ // =======================================================================
+
+ Int_t fCentralityBin; // Centrality bin of current event within max centrality bin
+ Float_t fCentralityPercentile; // Centrality percentile of current event
+ // ----------------------------------------------------------------------
+ Int_t fCentralityBinMax; // Max centrality bin to be used
+ Float_t fVertexZMax; // VertexZ cut
+ Float_t fRapidityMax; // Rapidity cut
+ Float_t fMinTrackLengthMC; // Min track length for MC tracks
+ Float_t fNSigmaMaxCdd; // N Sigma for dcar / sqrt(cdd) - turn off with 0.
+ Float_t fNSigmaMaxCzz; // N Sigma for dcaz / sqrt(czz) - turn off with 0.
+ // -----------------------------------------------------------------------
+ AliPID::EParticleType fParticleSpecies; // Particle species on basis of AliPID
+ Int_t fControlParticleCode; // PDG code control particle
+ Bool_t fControlParticleIsNeutral; // Is control particle neutral
+ TString fControlParticleName; // Name of control particle
+ // -----------------------------------------------------------------------
+ Float_t fNSigmaMaxTPC; // N Sigma for TPC PID
+ Float_t fNSigmaMaxTOF; // N Sigma for TOF PID
+ Float_t fMinPtForTOFRequired; // Min pt from where TOF is required
+
+ // =======================================================================
+
+ TH1F *fHEventStat0; // Event cut statistics
+ TH1F *fHEventStat1; // Event cut statistics - incremental
+ Int_t fHEventStatMax; // Max N cuts to be included in HEventStat
+ // -----------------------------------------------------------------------
+ TH1F *fHTriggerStat; // Trigger statistics
+ Int_t fNTriggers; // N triggers used
+ // -----------------------------------------------------------------------
+ TH1F *fHCentralityStat; // Centrality statistics
+ Int_t fNCentralityBins; // N centrality bins used
+
+ // =======================================================================
+
+ TF1 *fEtaCorrFunc; //! Eta correction function for TPC dE/dx
+ THnSparseF ***fCorr0; // Correction matrices for particle / anti-particle
+ THnSparseF ***fCorr1; // Correction matrices [cross section corrected] matrices for particle/ anti-particle
+ // -----------------------------------------------------------------------
+
+ ClassDef(AliAnalysisNetParticleHelper, 1);
+};
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#include "TFile.h"
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TCanvas.h"
+#include "TStopwatch.h"
+#include "TMath.h"
+#include "THashList.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliTracker.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDpid.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliKineTrackCuts.h"
+#include "AliMCParticle.h"
+#include "AliESDVZERO.h"
+#include "AliAnalysisTaskNetParticle.h"
+#include "AliGenEventHeader.h"
+#include "AliCentrality.h"
+
+using namespace std;
+
+// Task for NetParticle checks
+// Author: Jochen Thaeder <jochen@thaeder.de>
+
+ClassImp(AliAnalysisTaskNetParticle)
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+AliAnalysisTaskNetParticle::AliAnalysisTaskNetParticle(const char *name) :
+ AliAnalysisTaskSE(name),
+ fHelper(new AliAnalysisNetParticleHelper),
+ fEffCont(NULL),
+ fDCA(NULL),
+ fDist(NULL),
+
+ fOutList(NULL),
+ fOutListEff(NULL),
+ fOutListCont(NULL),
+ fOutListDCA(NULL),
+ fOutListQA(NULL),
+
+ fESD(NULL),
+ fESDHandler(NULL),
+
+ fESDTrackCutsBase(new AliESDtrackCuts),
+ fESDTrackCuts(NULL),
+ fESDTrackCutsBkg(NULL),
+ fESDTrackCutsEff(NULL),
+
+ fIsMC(kFALSE),
+ fESDTrackCutMode(0),
+ fModeEffCreation(0),
+ fModeDCACreation(0),
+ fModeDistCreation(0),
+
+ fMCEvent(NULL),
+ fMCStack(NULL),
+
+ fHnQA(NULL),
+ fUseQATHnSparse(kFALSE),
+
+ fEtaMax(0.9),
+ fPtRange(new Float_t[2]),
+ fPtRangeEff(new Float_t[2]) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliAnalysisTaskNetParticle",10);
+
+ fPtRange[0] = 0.4;
+ fPtRange[1] = 0.8;
+ fPtRangeEff[0] = 0.2;
+ fPtRangeEff[1] = 1.6;
+
+ // -- Output slots
+ // -------------------------------------------------
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
+ DefineOutput(3, TList::Class());
+ DefineOutput(4, TList::Class());
+ DefineOutput(5, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskNetParticle::~AliAnalysisTaskNetParticle() {
+ // Destructor
+
+ if (fPtRange) delete[] fPtRange;
+ if (fPtRangeEff) delete[] fPtRangeEff;
+
+ if (fESDTrackCutsBase) delete fESDTrackCutsBase;
+ if (fESDTrackCuts) delete fESDTrackCuts;
+ if (fESDTrackCutsBkg) delete fESDTrackCutsBkg;
+ if (fESDTrackCutsEff) delete fESDTrackCutsEff;
+
+ if (fHelper) delete fHelper;
+ if (fEffCont) delete fEffCont;
+ if (fDCA) delete fDCA;
+ if (fDist) delete fDist;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliAnalysisTaskNetParticle::UserCreateOutputObjects() {
+ // Create histograms
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ fOutList = new TList;
+ fOutList->SetName(GetName()) ;
+ fOutList->SetOwner(kTRUE);
+
+ fOutListEff = new TList;
+ fOutListEff->SetName(Form("%s_eff",GetName()));
+ fOutListEff->SetOwner(kTRUE) ;
+
+ fOutListCont = new TList;
+ fOutListCont->SetName(Form("%s_cont",GetName()));
+ fOutListCont->SetOwner(kTRUE) ;
+
+ fOutListDCA = new TList;
+ fOutListDCA->SetName(Form("%s_dca",GetName()));
+ fOutListDCA->SetOwner(kTRUE) ;
+
+ fOutListQA = new TList;
+ fOutListQA->SetName(Form("%s_qa",GetName()));
+ fOutListQA->SetOwner(kTRUE) ;
+
+ // ------------------------------------------------------------------
+ // -- Get event / trigger statistics histograms
+ // ------------------------------------------------------------------
+
+ fOutList->Add(fHelper->GetHEventStat0());
+ fOutList->Add(fHelper->GetHEventStat1());
+ fOutList->Add(fHelper->GetHTriggerStat());
+ fOutList->Add(fHelper->GetHCentralityStat());
+
+ // ------------------------------------------------------------------
+ // -- Add histograms from distribution class
+ // ------------------------------------------------------------------
+ if (fModeDistCreation == 1)
+ fDist->CreateHistograms(fOutList);
+
+ // ------------------------------------------------------------------
+ // -- Add histograms from efficiency/contamination class
+ // ------------------------------------------------------------------
+ if (fIsMC && fModeEffCreation == 1) {
+ fOutListEff->Add(fEffCont->GetHnEff());
+ fOutListCont->Add(fEffCont->GetHnCont());
+ }
+
+ // ------------------------------------------------------------------
+ // -- Add histograms from DCA class
+ // ------------------------------------------------------------------
+ if (fModeDCACreation == 1)
+ fOutListDCA->Add(fDCA->GetHnDCA());
+
+ // ------------------------------------------------------------------
+ // -- Add Multiplicity Correlations
+ // ------------------------------------------------------------------
+ /*
+ fOutList->Add(new TH2F("fHMultCorrTPCref0", "Centrality vs Multiplicity (TPCref);Centrality;Multiplicity (TPCref)",
+ 24, -0.5, 11.49, 10001, 0.5, 10000.49));
+ fOutList->Add(new TH2F("fHMultCorrTracks0", "Centrality vs Multiplicity (Tracks);Centrality;Multiplicity (Tracks)",
+ 24, -0.5, 11.49, 10001, 0.5, 10000.49));
+ fOutList->Add(new TH2F("fHMultCorrESDTracks0", "Centrality vs Multiplicity (ESDTracks);Centrality;Multiplicity (ESDTracks)",
+ 24, -0.5, 11.49, 10001, 0.5, 10000.49));
+
+ fOutList->Add(new TH2F("fHMultCorrTPCref1", "Centrality vs Multiplicity (TPCref);Centrality;Multiplicity (TPCref)",
+ 101, -0.5, 100.49, 10001, 0.5, 10000.49));
+ fOutList->Add(new TH2F("fHMultCorrTracks1", "Centrality vs Multiplicity (Tracks);Centrality;Multiplicity (Tracks)",
+ 101, -0.5, 100.49, 10001, 0.5, 10000.49));
+ fOutList->Add(new TH2F("fHMultCorrESDTracks1", "Centrality vs Multiplicity (ESDTracks);Centrality;Multiplicity (ESDTracks)",
+ 101, -0.5, 100.49, 10001, 0.5, 10000.49));
+ */
+ // ------------------------------------------------------------------
+ // -- Create THnSparseF - QA
+ // ------------------------------------------------------------------
+ /*
+ Double_t dCent[2] = {-0.5, 8.5};
+ Int_t iCent = 9;
+
+ Double_t dEta[2] = {-0.9, 0.9}; // -> 37 bins
+ Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
+
+ Double_t dRap[2] = {-0.5, 0.5};
+ Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
+
+ Double_t dPhi[2] = {0.0, TMath::TwoPi()};
+ Int_t iPhi = 42;
+
+ Double_t dPt[2] = {0.1, 3.0};
+ Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
+
+ Double_t dSign[2] = {-1.5, 1.5};
+ Int_t iSign = 3;
+
+ // cent:isAccepted: pInner: pt: sign:tpcSignal:nSigmaTPC:nSigmaTOF: eta: phi: y: dcar: dcaz: nSigmaCdd: nSigmaCzz
+ Int_t bin[15] = { iCent, 2 , iPt, iPt, iSign, 500, 50 , 50 , iEta, iPhi, iRap, 50 , 50 , 50 , 50 };
+ Double_t min[15] = {dCent[0], -0.5, dPt[0], dPt[0], dSign[0], 30, -5.0, -5.0, dEta[0], dPhi[0], dRap[0], -10., -10., -10., -10.};
+ Double_t max[15] = {dCent[1], 1.5, dPt[1], dPt[1], dSign[1], 500, 5.0, 5.0, dEta[1], dPhi[1], dRap[1], 10., 10., 10., 10.};
+
+ if (fUseQATHnSparse) {
+ fOutListQA->Add(new THnSparseF("fHnQA", "cent:isAccepted:pInner:pt:sign:tpcSignal:nSigmaTPC:nSigmaTOF:eta:phi:u:dcar:dcaz:nSigmaCdd:nSigmaCzz",
+ 15, bin, min, max));
+
+ fHnQA = static_cast<THnSparseF*>(fOutListQA->Last());
+ fHnQA->Sumw2();
+
+ fHnQA->GetAxis(0)->SetTitle("centrality");
+ fHnQA->GetAxis(1)->SetTitle("isAccepted");
+ fHnQA->GetAxis(2)->SetTitle("#it{p}_{Inner} (GeV/#it{c})");
+ fHnQA->GetAxis(3)->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ fHnQA->GetAxis(4)->SetTitle("sign");
+
+ fHnQA->GetAxis(5)->SetTitle("TPC signal");
+ fHnQA->GetAxis(6)->SetTitle("n #sigma TPC");
+ fHnQA->GetAxis(7)->SetTitle("n #sigma TOF");
+
+ fHnQA->GetAxis(8)->SetTitle("#eta");
+ fHnQA->GetAxis(9)->SetTitle("#varphi");
+ fHnQA->GetAxis(10)->SetTitle("#it{y}");
+
+ fHnQA->GetAxis(11)->SetTitle("DCAr");
+ fHnQA->GetAxis(12)->SetTitle("DCAz");
+
+ fHnQA->GetAxis(13)->SetTitle("n #sigma #sqrt(Cdd)/DCAr");
+ fHnQA->GetAxis(14)->SetTitle("n #sigma #sqrt(Czz)/DCAz");
+
+ fHelper->BinLogAxis(fHnQA, 2);
+ fHelper->BinLogAxis(fHnQA, 3);
+ }
+ */
+ // ------------------------------------------------------------------
+
+ TH1::AddDirectory(oldStatus);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskNetParticle::UserExec(Option_t *) {
+ // Called for each event
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Setup Event
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ if (SetupEvent() < 0) {
+ PostData(1,fOutList);
+ PostData(2,fOutListEff);
+ PostData(3,fOutListCont);
+ PostData(4,fOutListDCA);
+ PostData(5,fOutListQA);
+ return;
+ }
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Process Efficiency / Contamination Determination
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ if (fIsMC && fModeEffCreation == 1)
+ fEffCont->Process();
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Process DCA Determination
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ if (fModeDCACreation == 1)
+ fDCA->Process();
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Process Distributions
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ if (fModeDistCreation == 1)
+ fDist->Process();
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Post output data
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ PostData(1,fOutList);
+ PostData(2,fOutListEff);
+ PostData(3,fOutListCont);
+ PostData(4,fOutListDCA);
+ PostData(5,fOutListQA);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskNetParticle::Terminate(Option_t *){
+ // Terminate
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskNetParticle::Initialize() {
+ // Initialize event
+
+ // ------------------------------------------------------------------
+ // -- ESD TrackCuts
+ // ------------------------------------------------------------------
+ TString sModeName("");
+
+ // -- Create ESD track cuts
+ // --------------------------
+ fESDTrackCutsBase->SetMinNCrossedRowsTPC(70); // TPC
+ fESDTrackCutsBase->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); // TPC
+
+ fESDTrackCutsBase->SetMaxChi2PerClusterTPC(4); // TPC
+ fESDTrackCutsBase->SetAcceptKinkDaughters(kFALSE); // TPC
+ fESDTrackCutsBase->SetRequireTPCRefit(kTRUE); // TPC
+
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff); // ITS
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kOff); // ITS
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSSD,AliESDtrackCuts::kOff); // ITS
+
+ fESDTrackCutsBase->SetDCAToVertex2D(kFALSE); // VertexConstrained
+ fESDTrackCutsBase->SetRequireSigmaToVertex(kFALSE); // VertexConstrained
+
+ fESDTrackCutsBase->SetEtaRange(-1.*fEtaMax, fEtaMax); // Acceptance
+ fESDTrackCutsBase->SetPtRange(fPtRange[0],fPtRange[1]); // Acceptance
+
+ // -- Mode : clean cuts -> small contamination
+ if (fESDTrackCutMode == 0) {
+ sModeName = "Clean";
+ fESDTrackCutsBase->SetRequireITSRefit(kTRUE); // ITS
+ fESDTrackCutsBase->SetMaxChi2PerClusterITS(36); // ITS
+ }
+ // -- Mode : dirty cuts -> high efficiency
+ else if (fESDTrackCutMode == 1) {
+ sModeName = "Dirty";
+ fESDTrackCutsBase->SetRequireITSRefit(kFALSE); // ITS
+ }
+ // -- Mode : Default
+ else {
+ sModeName = "Base";
+ }
+
+ fESDTrackCutsBase->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
+
+ // -- Create ESD BKG track cuts
+ // ------------------------------
+ fESDTrackCutsBkg = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
+ fESDTrackCutsBkg->SetName(Form("NetParticleCuts2010_%s_Bkg",sModeName.Data()));
+ fESDTrackCutsBkg->SetMaxDCAToVertexZ(10.); // VertexConstrained
+
+ // -- Create ESD track cuts
+ // ------------------------------
+ fESDTrackCuts = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
+ fESDTrackCuts->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
+ fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); // VertexConstrained -> 7*(0.0026+0.0050/pt^1.01)
+ fESDTrackCuts->SetMaxChi2TPCConstrainedGlobal(36); // VertexConstrained
+ fESDTrackCuts->SetMaxDCAToVertexZ(2); // VertexConstrained
+
+ // -- Create ESD Eff track cuts
+ // ------------------------------
+ fESDTrackCutsEff = static_cast<AliESDtrackCuts*>(fESDTrackCuts->Clone());
+ fESDTrackCutsEff->SetName(Form("NetParticleCuts2010_%s_Eff",sModeName.Data()));
+ fESDTrackCutsEff->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]); // Acceptance
+
+ // ------------------------------------------------------------------
+ // -- Initialize Helper
+ // ------------------------------------------------------------------
+ if (fHelper->Initialize(fIsMC))
+ return -1;
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize Efficiency/Contamination
+ // ------------------------------------------------------------------
+ if (fIsMC && fModeEffCreation == 1) {
+ fEffCont = new AliAnalysisNetParticleEffCont;
+ fEffCont->Initialize(fESDTrackCutsEff, fHelper);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize DCA Determination
+ // ------------------------------------------------------------------
+ if (fModeDCACreation == 1) {
+ fDCA = new AliAnalysisNetParticleDCA;
+ fDCA->Initialize(fESDTrackCuts, fESDTrackCutsBkg, fHelper);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize DCA Determination
+ // ------------------------------------------------------------------
+ if (fModeDistCreation == 1) {
+ fDist = new AliAnalysisNetParticleDistribution;
+ fDist->Initialize(fHelper, fESDTrackCuts, fIsMC, fPtRange, fEtaMax);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Reset Event
+ // ------------------------------------------------------------------
+ ResetEvent();
+
+ return 0;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskNetParticle::SetupEvent() {
+ // Setup Reading of event
+ // > return 0 for success / accepted event
+ // > return -1 for failed setup
+ // > return -2 for rejected event
+
+ ResetEvent();
+
+ // -- ESD Event
+ // ------------------------------------------------------------------
+ if (SetupESDEvent() < 0) {
+ AliError("Setup ESD Event failed");
+ return -1;
+ }
+
+ // -- Setup MC Event
+ // ------------------------------------------------------------------
+ if (fIsMC && SetupMCEvent() < 0) {
+ AliError("Setup MC Event failed");
+ return -1;
+ }
+
+ // -- Setup Event for Helper / EffCont / DCA / Dist classes
+ // ------------------------------------------------------------------
+ fHelper->SetupEvent(fESDHandler, fMCEvent);
+
+ if (fIsMC && fModeEffCreation)
+ fEffCont->SetupEvent(fESDHandler, fMCEvent);
+
+ if (fModeDCACreation == 1)
+ fDCA->SetupEvent(fESDHandler, fMCEvent);
+
+ if (fModeDistCreation == 1)
+ fDist->SetupEvent(fESDHandler, fMCEvent);
+
+ // -- Evaluate Event cuts
+ // ------------------------------------------------------------------
+ return fHelper->IsEventRejected() ? -2 : 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskNetParticle::SetupESDEvent() {
+ // -- Setup ESD Event
+ // > return 0 for success
+ // > return -1 for failed setup
+
+ fESDHandler= dynamic_cast<AliESDInputHandler*>
+ (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ if (!fESDHandler) {
+ AliError("Could not get ESD input handler");
+ return -1;
+ }
+
+ fESD = fESDHandler->GetEvent();
+ if (!fESD) {
+ AliError("Could not get ESD event");
+ return -1;
+ }
+
+ // -- Check PID response
+ // ------------------------------------------------------------------
+ if (!fESDHandler->GetPIDResponse()) {
+ AliError("Could not get PID response");
+ return -1;
+ }
+
+ // -- Check Vertex
+ // ------------------------------------------------------------------
+ if (!fESD->GetPrimaryVertexTracks()) {
+ AliError("Could not get vertex from tracks");
+ return -1;
+ }
+
+ // -- Check Centrality
+ // ------------------------------------------------------------------
+ if (!fESD->GetCentrality()) {
+ AliError("Could not get centrality");
+ return -1;
+ }
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskNetParticle::SetupMCEvent() {
+ // -- Setup MC Event
+ // > return 0 for success
+ // > return -1 for failed setup
+
+ AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>
+ (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+
+ fMCEvent = mcH->MCEvent();
+ if (!fMCEvent) {
+ AliError("MC event not available");
+ return -1;
+ }
+
+ // -- Get MC header
+ // ------------------------------------------------------------------
+ AliHeader* header = fMCEvent->Header();
+ if (!header) {
+ AliError("MC header not available");
+ return -1;
+ }
+
+ // -- Check Stack
+ // ------------------------------------------------------------------
+ fMCStack = fMCEvent->Stack();
+ if (!fMCStack) {
+ AliError("MC stack not available");
+ return -1;
+ }
+
+ // -- Check GenHeader
+ // ------------------------------------------------------------------
+ if (!header->GenEventHeader()) {
+ AliError("Could not retrieve genHeader from header");
+ return -1;
+ }
+
+ // -- Check primary vertex
+ // ------------------------------------------------------------------
+ if (!fMCEvent->GetPrimaryVertex()){
+ AliError("Could not get MC vertex");
+ return -1;
+ }
+
+ return 0;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskNetParticle::ResetEvent() {
+ // -- Reset event
+
+ // -- Reset ESD Event
+ fESD = NULL;
+
+ // -- Reset MC Event
+ if (fIsMC)
+ fMCEvent = NULL;
+
+ // -- Reset Dist Creation
+ if (fModeDistCreation == 1)
+ fDist->ResetEvent();
+
+ return;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Process Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+#if 0
+
+//________________________________________________________________________
+void AliAnalysisTaskNetParticle::ProcessESDTrackLoop() {
+ // -- Process ESD tracks and fill histograms
+
+ Int_t nTracks = 0;
+
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+ AliESDtrack *track = fESD->GetTrack(idxTrack);
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Check track cuts
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ // -- Check if charged track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted - BKG
+ if (!fESDTrackCutsBkg->AcceptTrack(track))
+ continue;
+
+ // -- Check if accepted in rapidity window
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP))
+ continue;
+
+ // -- Check if accepted bt PID from TPC or TPC+TOF
+ Double_t pid[2];
+ Bool_t isAcceptedPID = fHelper->IsTrackAcceptedPID(track, pid);
+
+ // -- Check if accepted with thighter DCA cuts
+ Bool_t isAcceptedDCA = fHelper->IsTrackAcceptedDCA(track);
+
+ // -- Check track cuts
+ Bool_t isAcceptedVertex = fESDTrackCuts->AcceptTrack(track);
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill tracks in QA THnSparseF
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ if (fUseQATHnSparse && fHnQA->GetEntries() < 5e5) {
+
+ // -- Get dca r/z
+ Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
+ track->GetImpactParameters(dca,cov);
+
+ Float_t dcaRoverCdd = ( TMath::Sqrt(cov[0]) != 0. ) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
+ Float_t dcaZoverCzz = ( TMath::Sqrt(cov[2]) != 0. ) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
+
+ Double_t aQA[15] = {
+ Double_t(fHelper->GetCentralityBin()), // 0 centrality
+ Double_t(isAcceptedVertex || isAcceptedDCA), // 1 isAccepted -> Vertex || isAcceptedDCA
+ track->GetInnerParam()->GetP(), // 2 p at InnerParam
+ track->Pt(), // 3 pt
+ track->GetSign(), // 4 sign
+ track->GetTPCsignal(), // 5 TPC dE/dx
+ pid[0], // 6 n Sigma TPC
+ pid[1], // 7 n Sigma TOF
+ track->Eta(), // 8 eta
+ track->Phi(), // 9 phi
+ yP, // 10 rapidity
+ dca[0], // 11 dca r
+ dca[1], // 12 dca z
+ dcaRoverCdd, // 13 sqrt(cov[dd])
+ dcaZoverCzz, // 14 sqrt(cov[zz])
+ };
+
+ fHnQA->Fill(aQA);
+ }
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Apply track cuts
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ if (!isAcceptedVertex || !isAcceptedDCA)
+ continue;
+
+ ++nTracks;
+
+ if (!isAcceptedPID)
+ continue;
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+ // -- Fill Particle Fluctuation Histograms
+ // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+
+ /*
+ Int_t tpcRef = fESDTrackCuts->GetReferenceMultiplicity(fESD,kTRUE);
+
+ (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTPCref0")))->Fill(fHelper->GetCentralityBin(), tpcRef);
+ (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTracks0")))->Fill(fHelper->GetCentralityBin(), nTracks);
+
+ Float_t centPercentile = fHelper->GetCentralityPercentile();
+ (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTPCref1")))->Fill(centPercentile, tpcRef);
+ (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTracks1")))->Fill(centPercentile, nTracks);
+*/
+ return;
+}
+
+#endif
--- /dev/null
+//-*- Mode: C++ -*-
+
+#ifndef ALIANALYSISTASKNETPARTICLE_H
+#define ALIANALYSISTASKNETPARTICLE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+// Task for NetParticle Distributions
+// Authors: Jochen Thaeder <jochen@thaeder.de>
+
+#include "AliAnalysisTaskSE.h"
+#include "TList.h"
+#include "THnSparse.h"
+
+#include "AliAnalysisNetParticleHelper.h"
+#include "AliAnalysisNetParticleEffCont.h"
+#include "AliAnalysisNetParticleDCA.h"
+#include "AliAnalysisNetParticleDistribution.h"
+
+class AliESDEvent;
+class AliESDInputHandler;
+class AliMCEvent;
+class AliStack;
+class AliPIDResponse;
+class AliESDResponse;
+
+class AliAnalysisTaskNetParticle : public AliAnalysisTaskSE {
+
+ public:
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Constructor / Destructor
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisTaskNetParticle(const char *name = "AliAnalysisTaskNetParticle");
+ virtual ~AliAnalysisTaskNetParticle();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Task Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ virtual void UserCreateOutputObjects();
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *);
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Initialize Task */
+ Int_t Initialize();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Setter
+ * ---------------------------------------------------------------------------------
+ */
+
+ void SetIsMC() {fIsMC = kTRUE;}
+ void SetUseQATHnSparse(Bool_t b) {fUseQATHnSparse = b;}
+
+ void SetESDTrackCutMode(Int_t i) {fESDTrackCutMode = i;}
+ void SetModeEffCreation(Int_t i) {fModeEffCreation = i;}
+ void SetModeDCACreation(Int_t i) {fModeDCACreation = i;}
+ void SetModeDistCreation(Int_t i) {fModeDistCreation = i;}
+
+ void SetEtaMax(Float_t f) {fEtaMax = f;}
+ void SetPtRange(Float_t f1, Float_t f2) {fPtRange[0] = f1; fPtRange[1] = f2;}
+ void SetPtRangeEff(Float_t f1, Float_t f2) {fPtRangeEff[0] = f1; fPtRangeEff[1] = f2;}
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Setter - Pass-Through
+ * ---------------------------------------------------------------------------------
+ */
+
+ void SetCentralityBinMax(Int_t d) {fHelper->SetCentralityBinMax(d);}
+ void SetVertexZMax(Float_t f) {fHelper->SetVertexZMax(f);}
+
+ void SetRapidityMax(Float_t f) {fHelper->SetRapidityMax(f);}
+ void SetMinTrackLengthMC(Float_t f) {fHelper->SetMinTrackLengthMC(f);}
+ void SetNSigmaMaxCdd(Float_t f) {fHelper->SetNSigmaMaxCdd(f);}
+ void SetNSigmaMaxCzz(Float_t f) {fHelper->SetNSigmaMaxCzz(f);}
+
+ void SetParticleSpecies(AliPID::EParticleType pid) {fHelper->SetParticleSpecies(pid);}
+ void SetControlParticleSpecies(Int_t pdgCode, Bool_t isNeutral, const Char_t *name) {
+ fHelper->SetControlParticleSpecies(pdgCode, isNeutral, name);
+ }
+
+ void SetNSigmaMaxTPC(Float_t f) {fHelper->SetNSigmaMaxTPC(f);}
+ void SetNSigmaMaxTOF(Float_t f) {fHelper->SetNSigmaMaxTOF(f);}
+ void SetMinPtForTOFRequired(Float_t f) {fHelper->SetMinPtForTOFRequired(f);}
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliAnalysisTaskNetParticle(const AliAnalysisTaskNetParticle&); // not implemented
+ AliAnalysisTaskNetParticle& operator=(const AliAnalysisTaskNetParticle&); // not implemented
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Setup Event */
+ Int_t SetupEvent();
+
+ /** Setup ESD Event */
+ Int_t SetupESDEvent();
+
+ /** Setup MC Event */
+ Int_t SetupMCEvent();
+
+ /** Reset Event */
+ void ResetEvent();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliAnalysisNetParticleHelper *fHelper; //! Helper class
+ AliAnalysisNetParticleEffCont *fEffCont; //! Efficiency and Contamination class
+ AliAnalysisNetParticleDCA *fDCA; //! DCA class
+ AliAnalysisNetParticleDistribution *fDist; //! Distributions class
+
+ // --- OutLists ----------------------------------------------------------
+
+ TList *fOutList; //! Ptr to output data container
+ TList *fOutListEff; //! Ptr to output data container - efficiency
+ TList *fOutListCont; //! Ptr to output data container - contamination
+ TList *fOutListDCA; //! Ptr to output data container - DCA
+ TList *fOutListQA; //! Ptr to QA output data container
+
+ // --- ESD only ----------------------------------------------------------
+
+ AliESDEvent *fESD; //! Ptr to ESD event
+ AliESDInputHandler *fESDHandler; //! Ptr to ESD Handler
+
+ AliESDtrackCuts *fESDTrackCutsBase; //! ESD cuts - base settings
+ AliESDtrackCuts *fESDTrackCuts; //! ESD cuts
+ AliESDtrackCuts *fESDTrackCutsBkg; //! ESD cuts for Bkg
+ AliESDtrackCuts *fESDTrackCutsEff; //! ESD cuts for efficiency determination -> larger pt Range
+
+ // --- Flags -------------------------------------------------------------
+
+ Bool_t fIsMC; // Is MC event
+ Int_t fESDTrackCutMode; // ESD track cut mode : 0 = clean | 1 dirty
+ Int_t fModeEffCreation ; // Correction creation mode : 1 = on | 0 = off
+ Int_t fModeDCACreation; // DCA creation mode : 1 = on | 0 = off
+ Int_t fModeDistCreation; // Dist creation mode : 1 = on | 0 = off
+
+ // --- MC only -----------------------------------------------------------
+
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+ AliStack *fMCStack; //! Ptr to MC stack
+
+ // -----------------------------------------------------------------------
+
+ THnSparseF *fHnQA; // THnSparseF : tracks for QA
+ Float_t fUseQATHnSparse; // Usage of THnSparse for QA
+
+ // -----------------------------------------------------------------------
+
+ Float_t fEtaMax; // Max, absolut eta
+ Float_t *fPtRange; // Array of pt [min,max]
+ Float_t *fPtRangeEff; // Array of pt [min,max] for efficiency
+
+ // -----------------------------------------------------------------------
+
+ ClassDef(AliAnalysisTaskNetParticle, 1);
+};
+
+#endif
--- /dev/null
+/* *********************************************************************************
+ * File : AddTaskNetParticle.C
+ * Author : Jochen Thaeder <jochen@thaeder.de>
+ * *********************************************************************************
+ * Configuring NetParticle Task:
+ * - ARGUMENTS :
+ * name -> Name of the task, containing partcile type :
+ * Currently : Proton, Pion, Kaon
+ * isModeDist -> Fill Distributions
+ * isModeEff -> Fill Efficiency/Contamination ThnSparse
+ * isModeDCA -> Fill DCA ThnSparse
+ * useQAThnSparse -> Fill QA ThnSparse
+ * isCreateCSC -> Prepare for CrossSectionCorrection
+ * - requires isModeEff to be set
+ * - Proton only
+ *
+ * - OUTPUT CONTAINER : #N = 5
+ * (1) - Standard Output, Distributions
+ * (2) - Efficiency ThnSparse
+ * (3) - Contamination ThnSparse
+ * (4) - DCA ThnSparse
+ * (5) - QA ThnSparse
+ *
+ ********************************************************************************* */
+
+AliAnalysisTask *AddTaskNetParticle(const Char_t * name = "jthaeder_NetProton",
+ Bool_t isModeDist, Bool_t isModeEff, Bool_t isModeDCA, Bool_t useQAThnSparse = kFALSE,
+ Bool_t isCreateCSC = kFALSE) {
+
+ TString sName(name);
+
+ if (isCreateCSC && !isModeEff) {
+ Error("AddTaskNetParticle", "Creating CrossSectionCorrection needs 'isModeEff' to be set.");
+ return NULL;
+ }
+
+ // ----------------------------------------------
+ // -- Get the current analysis manager
+ // ----------------------------------------------
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ Error("AddTaskNetParticle", "No analysis manager found.");
+ return NULL;
+ }
+
+ // ----------------------------------------------
+ // -- Check for MC
+ // ----------------------------------------------
+ Bool_t isMC = (mgr->GetMCtruthEventHandler() != NULL);
+ if (isMC)
+ Info("AddTaskNetParticle", "This task has MC.");
+
+ // ----------------------------------------------
+ // -- Create task configure it
+ // ----------------------------------------------
+ AliAnalysisTaskNetParticle *task = new AliAnalysisTaskNetParticle("AliAnalysisTaskNetParticle");
+ if (!task) {
+ Error("AddTaskNetParticle", "Task could not be created.");
+ return NULL;
+ }
+
+ if (isMC)
+ task->SetIsMC();
+
+ // -- Set particle type
+ // ----------------------------------------------
+ Float_t minPt, maxPt, minPtEff, maxPtEff, minPtForTOF;
+
+ if (sName.Contains("Proton")) {
+ task->SetParticleSpecies(AliPID::kProton);
+ task->SetControlParticleSpecies(3122, kTRUE, "Lambda");
+ minPt = 0.4; maxPt = 2.2;
+ minPtEff = 0.2; maxPtEff = 2.6;
+ minPtForTOF = 0.8;
+ if (isCreateCSC) {
+ minPtForTOF = maxPtEff;
+ }
+ }
+ else if (sName.Contains("Pion")) {
+ task->SetParticleSpecies(AliPID::kPion);
+ task->SetControlParticleSpecies(3122, kTRUE, "Lambda"); /// maybe something else ...
+ minPt = 0.3; maxPt = 0.9;
+ minPtEff = 0.2; maxPtEff = 1.2;
+ minPtForTOF = 0.8;
+ }
+ else if (sName.Contains("Kaon")) {
+ task->SetParticleSpecies(AliPID::kKaon);
+ task->SetControlParticleSpecies(3122, kTRUE, "Lambda"); /// maybe something else ...
+ minPt = 0.2; maxPt = 0.4;
+ minPtEff = 0.1; maxPtEff = 0.8;
+ minPtForTOF = 0.8;
+ }
+ else {
+ Error("AddTaskNetParticle", "Unknown Particle type.");
+ return NULL;
+ }
+
+ // -- Configure flags
+ // ----------------------------------------------
+ if (isModeEff)
+ task->SetModeEffCreation(1); // => 1 = on | 0 = off (default)
+ if (isModeDCA)
+ task->SetModeDCACreation(1); // => 1 = on | 0 = off (default)
+ if (isModeDist)
+ task->SetModeDistCreation(1); // => 1 = on | 0 = off (default)
+
+ // -- Enable QA plots
+ if (useQAThnSparse)
+ task->SetUseQATHnSparse(kTRUE);
+
+ // -- Configure cuts
+ // ----------------------------------------------
+
+ // -- Set cut flags ...
+ task->SetESDTrackCutMode(0); // => 0 = clean | 1 = dirty
+
+ // -- Set standard event cuts
+ task->SetVertexZMax(10.);
+ task->SetCentralityBinMax(7);
+
+ // -- Set track event cuts
+ task->SetRapidityMax(0.5);
+ task->SetMinTrackLengthMC(70.);
+ task->SetNSigmaMaxCdd(3.);
+ task->SetNSigmaMaxCzz(3.);
+
+ task->SetNSigmaMaxTPC(2.5);
+ task->SetNSigmaMaxTOF(2.5);
+ task->SetMinPtForTOFRequired(minPtForTOF);
+
+ // -- Set analysis ranges
+ task->SetEtaMax(0.9);
+ task->SetPtRange(minPt, maxPt); // pt cut range for the analysis
+ task->SetPtRangeEff(minPtEff, maxPtEff); // pt cut range for the correction / efficiency / contamination creation
+
+ // ----------------------------------------------
+ // -- Initialize and add task to the ANALYSIS manager
+ // ----------------------------------------------
+ if (!task->Initialize())
+ mgr->AddTask(task);
+ else {
+ Error("Initialize failed, not adding AliAnalysistaskNetParticle !!!");
+ delete task;
+ return NULL;
+ }
+
+ // ----------------------------------------------
+ // -- data containers - input
+ // ----------------------------------------------
+ AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+
+ // ----------------------------------------------
+ // -- data containers - output
+ // ----------------------------------------------
+ AliAnalysisDataContainer *coutput = mgr->CreateContainer(name, TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s.root", name));
+ AliAnalysisDataContainer *coutputEff = mgr->CreateContainer(Form("%s_eff", name), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s.root", name));
+ AliAnalysisDataContainer *coutputCont = mgr->CreateContainer(Form("%s_cont", name), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s.root", name));
+ AliAnalysisDataContainer *coutputDca = mgr->CreateContainer(Form("%s_dca", name), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s.root", name));
+
+ AliAnalysisDataContainer *coutputQA = mgr->CreateContainer(Form("%sQA", name), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%sQA.root", name));
+
+ mgr->ConnectInput (task, 0, cinput );
+ mgr->ConnectOutput (task, 1, coutput);
+ mgr->ConnectOutput (task, 2, coutputEff);
+ mgr->ConnectOutput (task, 3, coutputCont);
+ mgr->ConnectOutput (task, 4, coutputDca);
+ mgr->ConnectOutput (task, 5, coutputQA);
+
+ return task;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff