--- /dev/null
+/* *************************************************************************
+ * AliEbyE Analysis for Particle Ratio Fluctuation *
+ * Deepika Rathee | Satyajit Jena *
+ * drathee@cern.ch | sjena@cern.ch *
+ * Date: Wed Jul 9 18:38:30 CEST 2014 *
+ * New approch to find particle ratio to reduce memory *
+ * (Test Only) *
+ ***************************************************************************/
+
+AliAnalysisTask *AddEbyEPidRatioTask(const Char_t *name = "TPC_NuDyn", Bool_t isModeDist = 1,
+ Bool_t isModeEff = 0, Bool_t isModeDCA = 0, Bool_t isModeQA = 0,
+ Bool_t isCreateCSC = 0, Bool_t isModeAOD = 0, Bool_t isSetExt = 0,
+ Int_t aodFilterBit = 1024,
+ Int_t modeCSC = 0, Int_t modeCuts = 0, Int_t modePID =-1,
+ Float_t gMinPt = 0.3, Float_t gMaxPt = 2.5,
+ Float_t gMinPtForTof = 0.21, Float_t gMaxPtForTPClow = 0.69,
+ Float_t gMinPtEff = 0.3, Float_t gMaxPtEff = 2.5,
+ Float_t gSigmaITS = 4.0, Float_t gSigmaTOF = 4.0,
+ Float_t gSigmaTPC = 4.0, Float_t gSigmaTPClow = 3.0) {
+
+ TString sName(name);
+
+ if (isCreateCSC && !isModeEff) {
+ Error("AddTaskNetParticle", "Creating CrossSectionCorrection needs 'isModeEff' to be set.");
+ return NULL;
+ }
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ Error("AddTaskNetParticle", "No analysis manager found.");
+ return NULL;
+ }
+
+ Bool_t isMC = (mgr->GetMCtruthEventHandler() != NULL);
+ if (isMC)
+ Info("AddTaskNetParticle", "This task has MC.");
+
+ AliEbyEPidRatioTask *task = new AliEbyEPidRatioTask("EbyEPidRatio");
+ if (!task) {
+ Error("EbyEPidRatio", "Task could not be created.");
+ return NULL;
+ }
+
+ if (isMC)
+ task->SetIsMC();
+ 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)
+ if (isModeAOD) {
+ task->SetIsAOD(1); // => 1 = AOD | 0 = ESD (default)
+ task->SetTrackFilterBit(aodFilterBit);
+ }
+ if (isModeQA)
+ task->SetModeQACreation(1); // => 1 = on | 0 = off (default)
+
+ Float_t minPt, maxPt, minPtEff, maxPtEff, minPtForTOF;
+ Float_t nSigmaITS, nSigmaTPC, nSigmaTPClow, nSigmaTOF, maxPtForTPClow;
+ Float_t etaMax, etaMaxEff;
+ Int_t pidStrategy;
+
+ minPtForTOF = 0.69;
+ maxPtForTPClow = 0.69;
+ minPt = 0.5;
+ maxPt = 2.0;
+ minPtEff = 0.3;
+ maxPtEff = 2.5;
+
+ etaMax = 0.8;
+ etaMaxEff = 0.8;
+ nSigmaITS = 4.0;
+ nSigmaTPC = 4.0;
+ nSigmaTPClow = 3.0;
+ nSigmaTOF = 4.0;
+
+ if (isCreateCSC && modeCSC == 1)
+ minPtForTOF = maxPtEff;
+
+ if (isSetExt) {
+ minPt = gMinPt;
+ maxPt = gMaxPt;
+ minPtForTOF = gMinPtForTof;
+ maxPtForTPClow = gMaxPtForTPClow;
+ minPtEff = gMinPtEff;
+ maxPtEff = gMaxPtEff;
+
+ nSigmaITS = gSigmaITS;
+ nSigmaTPC = gSigmaTPC;
+ nSigmaTPClow = gSigmaTPClow;
+ nSigmaTOF = gSigmaTOF;
+ }
+
+ if (modePID == -1) { // default
+ pidStrategy = 7; // 7: ITS + TPC + TOF (using minPtForTOF)
+ if (modeCuts == 1)
+ pidStrategy = 5; // 5: TPC + TOF (using minPtForTOF)
+ }
+ else
+ pidStrategy = modePID;
+
+ AliEbyEPidRatioHelper *helper = new AliEbyEPidRatioHelper;
+ if (!helper) {
+ Error("AddTaskNetParticle", "Helper could not be created.");
+ delete task;
+ return NULL;
+ }
+
+
+ task->SetESDTrackCutMode(modeCuts); // => 0 = normal | 1 = LF
+ task->SetEtaMax(etaMax); // eta cut
+ task->SetEtaMaxEff(etaMaxEff); // eta cut for efficiency
+ task->SetPtRange(minPt, maxPt); // pt cut range for the analysis
+ task->SetPtRangeEff(minPtEff, maxPtEff); // pt cut range for the correction / efficiency / contamination creation
+ helper->SetVertexZMax(10.);
+ helper->SetCentralityBinMax(7);
+ helper->SetRapidityMax(0.5);
+ helper->SetMinTrackLengthMC(70.);
+ helper->SetNSigmaMaxCdd(0.); // 3. || ->> Turn off sigmaDCA cuts for now
+ helper->SetNSigmaMaxCzz(0.); // 3. || ->> Turn off sigmaDCA cuts for now
+ helper->SetPhiRange(0., 3.88); // Only used if requested in task - default is TwoPi
+ helper->SetPIDStrategy(pidStrategy);
+ helper->SetNSigmaMaxITS(nSigmaITS);
+ helper->SetNSigmaMaxTPC(nSigmaTPC);
+ helper->SetNSigmaMaxTPClow(nSigmaTPClow);
+ helper->SetNSigmaMaxTOF(nSigmaTOF);
+ helper->SetMinPtForTOFRequired(minPtForTOF);
+ helper->SetMaxPtForTPClow(maxPtForTPClow);
+ task->SetNetParticleHelper(helper);
+ mgr->AddTask(task);
+
+ AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+ TString outputFileName = Form("%s:%s", AliAnalysisManager::GetCommonFileName(), name);
+ TString outputQAFileName = Form("%s:%s", AliAnalysisManager::GetCommonFileName(), name);
+
+ AliAnalysisDataContainer *coutput = mgr->CreateContainer(name, TList::Class(), AliAnalysisManager::kOutputContainer, outputFileName);
+ AliAnalysisDataContainer *coutputEff = mgr->CreateContainer(Form("%s_eff", name), TList::Class(), AliAnalysisManager::kOutputContainer, outputFileName);
+ AliAnalysisDataContainer *coutputCont = mgr->CreateContainer(Form("%s_cont", name), TList::Class(), AliAnalysisManager::kOutputContainer, outputFileName);
+ AliAnalysisDataContainer *coutputDca = mgr->CreateContainer(Form("%s_dca", name), TList::Class(), AliAnalysisManager::kOutputContainer, outputFileName);
+
+ AliAnalysisDataContainer *coutputQA = mgr->CreateContainer(Form("%sQA", name), TList::Class(), AliAnalysisManager::kOutputContainer, outputQAFileName);
+
+ 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;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "TMath.h"
+#include "TAxis.h"
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliESDtrackCuts.h"
+
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+#include "AliEbyEPidRatioBase.h"
+
+using namespace std;
+ClassImp(AliEbyEPidRatioBase)
+//________________________________________________________________________
+AliEbyEPidRatioBase::AliEbyEPidRatioBase() :
+ TNamed(),
+ fHelper(NULL),
+
+ fESD(NULL),
+ fESDTrackCuts(NULL),
+ fAOD(NULL),
+ fArrayMC(NULL),
+ fAODtrackCutBit(1024),
+ fIsMC(kFALSE),
+ fMCEvent(NULL),
+ fStack(NULL),
+
+ fCentralityBin(-1.),
+ fNTracks(0) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioBase",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioBase::AliEbyEPidRatioBase(const Char_t* name, const Char_t* title) :
+ TNamed(name, title),
+ fHelper(NULL),
+
+ fESD(NULL),
+ fESDTrackCuts(NULL),
+ fAOD(NULL),
+ fArrayMC(NULL),
+ fAODtrackCutBit(1024),
+ fIsMC(kFALSE),
+ fMCEvent(NULL),
+ fStack(NULL),
+
+ fCentralityBin(-1.),
+ fNTracks(0) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioBase",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioBase::~AliEbyEPidRatioBase() {
+ // Destructor
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioBase::Initialize(AliEbyEPidRatioHelper* helper, AliESDtrackCuts* cuts) {
+ fHelper = helper;
+ fESDTrackCuts = (cuts) ? cuts : helper->GetESDTrackCuts();
+ fIsMC = helper->GetIsMC();
+ fAODtrackCutBit = helper->GetAODtrackCutBit();
+ Init();
+ CreateHistograms();
+
+ return;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioBase::SetupEvent() {
+ // -- Setup event
+
+ ResetEvent();
+
+ // -- Get ESD objects
+ if (dynamic_cast<AliESDInputHandler*>(fHelper->GetInputEventHandler())) {
+ fESD = dynamic_cast<AliESDEvent*>(fHelper->GetInputEventHandler()->GetEvent());
+ fNTracks = fESD->GetNumberOfTracks();
+ }
+
+ // -- Get AOD objects
+ else if (dynamic_cast<AliAODInputHandler*>(fHelper->GetInputEventHandler())) {
+ fAOD = dynamic_cast<AliAODEvent*>(fHelper->GetInputEventHandler()->GetEvent());
+ fNTracks = fAOD->GetNumberOfTracks();
+
+ if (fIsMC) {
+ fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!fArrayMC)
+ AliFatal("No array of MC particles found !!!"); // MW no AliFatal use return values
+ }
+ }
+
+ if (fIsMC) {
+ fMCEvent = fHelper->GetMCEvent();
+ if (fMCEvent)
+ fStack = fMCEvent->Stack();
+ }
+
+ fCentralityBin = fHelper->GetCentralityBin();
+
+ return Setup();
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioBase::ResetEvent() {
+ // -- Reset ESD Event
+ fESD = NULL;
+ // -- Reset AOD Event
+ fAOD = NULL;
+ // -- Reset MC Event
+ if (fIsMC)
+ fMCEvent = NULL;
+ // -- Reset in class
+ Reset();
+ return;
+}
+
+
--- /dev/null
+#ifndef ALIEBYEPIDRATIOBASE_H
+#define ALIEBYEPIDRATIOBASE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "THnSparse.h"
+#include "AliEbyEPidRatioHelper.h"
+
+class AliESDEvent;
+class AliESDInputHandler;
+class AliMCEvent;
+class AliStack;
+class AliAODInputHandler;
+
+class AliEbyEPidRatioBase : public TNamed {
+
+ public:
+
+ AliEbyEPidRatioBase();
+ AliEbyEPidRatioBase(const Char_t* name, const Char_t* title);
+ virtual ~AliEbyEPidRatioBase();
+
+ void Initialize(AliEbyEPidRatioHelper* helper, AliESDtrackCuts* cuts = NULL);
+ Int_t SetupEvent();
+ void ResetEvent();
+ virtual void Process() = 0;
+
+ private:
+
+ AliEbyEPidRatioBase(const AliEbyEPidRatioBase&); // not implemented
+ AliEbyEPidRatioBase& operator=(const AliEbyEPidRatioBase&); // not implemented
+
+ protected:
+ virtual void Init() {};
+ virtual void CreateHistograms() {};
+ virtual void Reset() {};
+ virtual Int_t Setup() { return 0;};
+ AliEbyEPidRatioHelper *fHelper; //! Ptr to helper class
+ AliESDEvent *fESD; //! ESD object
+ AliESDtrackCuts *fESDTrackCuts; //! ESD cuts
+ AliAODEvent *fAOD; //! AOD object
+ TClonesArray *fArrayMC; //! array of MC particles
+ Int_t fAODtrackCutBit; // Track filter bit for AOD tracks
+ Bool_t fIsMC; // Is MC event
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+ AliStack *fStack; //! Ptr to stack
+ Float_t fCentralityBin; // Centrality of current event
+ Int_t fNTracks; // N Tracks in the current event
+
+ ClassDef(AliEbyEPidRatioBase, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "TMath.h"
+#include "TAxis.h"
+
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliESDtrackCuts.h"
+
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+
+#include "AliEbyEPidRatioDCA.h"
+
+using namespace std;
+
+ClassImp(AliEbyEPidRatioDCA)
+
+//________________________________________________________________________
+AliEbyEPidRatioDCA::AliEbyEPidRatioDCA() :
+ AliEbyEPidRatioBase("DCA", "DCA"),
+ fESDTrackCutsBkg(NULL),
+ fHnDCA(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioDCA",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioDCA::~AliEbyEPidRatioDCA() {
+ // Destructor
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioDCA::Process() {
+ Float_t etaRange[2];
+ fESDTrackCutsBkg->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCutsBkg->GetPtRange(ptRange[0],ptRange[1]);
+
+ // -- Track Loop
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
+
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted - ESD
+ if (fESD && !fESDTrackCutsBkg->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
+ continue;
+
+ // -- Check if accepted - AOD
+ if (fAOD){
+ AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
+
+ if (!trackAOD) {
+ AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
+ continue;
+ }
+ if (!trackAOD->TestFilterBit(fAODtrackCutBit))
+ continue;
+
+ // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
+ if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
+ continue;
+ }
+
+ Int_t gPdgCode = 0;
+ Int_t iPid = 0;
+ Double_t pid[3];
+ if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion))) { iPid = 1; gPdgCode = 211;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon))) { iPid = 2; gPdgCode = 321;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){ iPid = 3; gPdgCode = 2212;}
+ else iPid = 0;
+
+
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
+ continue;
+
+ Bool_t isDCArAccepted = fHelper->IsTrackAcceptedDCA(track);
+
+ // -- Check if accepted with thighter DCA cuts
+ // ?!?!? How to mimic this in AODs?
+ if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
+ isDCArAccepted = kFALSE;
+
+ // -- Check for contamination
+ Int_t contIdx = (fIsMC) ? GetContIdxTrack(TMath::Abs(track->GetLabel()), track->Charge(), gPdgCode) : 1;
+
+ // -- Get DCAs (dca_r, dca_z, sigma_xy, sigma_xy_z, sigma_z)
+ Float_t dca[2], cov[3]; //
+ if (fESD)
+ (dynamic_cast<AliESDtrack*>(track))->GetImpactParameters(dca, cov);
+ else
+ dca[0] = 1.;
+
+ // -- Fill THnSparse
+
+ if(iPid != 0) {
+ Double_t hnDCA[10] = {fCentralityBin,
+ track->Eta(),
+ yP,
+ track->Phi(),
+ track->Pt(),
+ track->Charge(),
+ contIdx,
+ isDCArAccepted,
+ dca[0],
+ 0};
+ fHnDCA->Fill(hnDCA);
+ }
+
+ Double_t hnDCA[10] = {fCentralityBin,
+ track->Eta(),
+ yP,
+ track->Phi(),
+ track->Pt(),
+ track->Charge(),
+ contIdx,
+ isDCArAccepted,
+ dca[0],
+ iPid};
+ fHnDCA->Fill(hnDCA);
+
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioDCA::CreateHistograms() {
+ Int_t binHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, AliEbyEPidRatioHelper::fgkfHistNBinsEta, // cent | etaRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsRap, AliEbyEPidRatioHelper::fgkfHistNBinsPhi, // yRec | phiRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsPt, AliEbyEPidRatioHelper::fgkfHistNBinsSign, // ptRec | sign
+ 4, 2, 77,4}; // contPart| DCArAccepted | DCAr
+
+ Double_t minHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[0], AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[0], AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
+ 0.5, -0.5, -3.,-0.5};
+
+ Double_t maxHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[1], AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[1], AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
+ 4.5, 1.5, 3., 3.5};
+
+
+ fHnDCA = new THnSparseD("hnDCA", "cent:etaRec:yRec:phiRec:ptRec:sign:contPart:contSign:DCArAccepted:DCAr", 10, 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_{Rec}"); // eta [-0.9, 0.9]
+ fHnDCA->GetAxis(2)->SetTitle("#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
+ fHnDCA->GetAxis(3)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
+ fHnDCA->GetAxis(4)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pT [ 0.2, 2.6]
+ 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 [-3, 3]
+ fHnDCA->GetAxis(9)->SetTitle("PID"); // 0 | 1 | 2 | 3
+
+ fHelper->BinLogAxis(fHnDCA, 4, fESDTrackCuts);
+ fHelper->BinLogAxis(fHnDCA, 4, fESDTrackCutsBkg);
+
+ // -- Set binning for DCAr
+ Double_t binsDCAr[77] = {-3.,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0.,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3.};
+ fHnDCA->GetAxis(8)->Set(76, binsDCAr);
+
+ // ------------------------------------------------------------------
+
+ return;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioDCA::GetContIdxTrack(Int_t label, Int_t sign, Int_t gPdgCode) {
+ Int_t contIdx = -1;
+
+ AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
+ if (!particle)
+ return contIdx;
+
+ Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
+ Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
+ Bool_t isSecondaryFromMaterial = (fESD) ? fStack->IsSecondaryFromMaterial(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
+
+ if (isPhysicalPrimary) {
+ if (gPdgCode == 0) {
+ // -- Check if correctly identified
+ if (particle->PdgCode() == (sign*gPdgCode))
+ contIdx = 1;
+ // -- MissIdentification
+ else
+ contIdx = 2;
+ }
+ else
+ contIdx = 1;
+ }
+
+ // -- Check if secondaries from material or weak decay
+ else if(isSecondaryFromWeakDecay)
+ contIdx = 3;
+ else if (isSecondaryFromMaterial)
+ contIdx = 4;
+ else
+ contIdx = -1;
+
+ return contIdx;
+}
+
+
+
--- /dev/null
+#ifndef ALIEBYEPIDRSTIODCA_H
+#define ALIEBYEPIDRSTIODCA_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "THnSparse.h"
+#include "AliEbyEPidRatioBase.h"
+
+class AliEbyEPidRatioDCA : public AliEbyEPidRatioBase {
+
+ public:
+
+ AliEbyEPidRatioDCA();
+ virtual ~AliEbyEPidRatioDCA();
+
+ virtual void Process();
+ void SetESDTrackCutsBkg(AliESDtrackCuts *p) {fESDTrackCutsBkg = p;}
+ THnSparseD* GetHnDCA() {return fHnDCA;}
+
+ private:
+
+ AliEbyEPidRatioDCA(const AliEbyEPidRatioDCA&); // not implemented
+ AliEbyEPidRatioDCA& operator=(const AliEbyEPidRatioDCA&); // not implemented
+
+ virtual void CreateHistograms();
+ Int_t GetContIdxTrack(Int_t label, Int_t sign, Int_t gPdgCode);
+ AliESDtrackCuts *fESDTrackCutsBkg; //! ESD cuts
+ THnSparseD *fHnDCA; // THnSparseD contamination DCA
+
+ ClassDef(AliEbyEPidRatioDCA, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "TMath.h"
+#include "TAxis.h"
+
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliESDtrackCuts.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+
+#include "AliEbyEPidRatioEffCont.h"
+
+using namespace std;
+
+
+ClassImp(AliEbyEPidRatioEffCont)
+
+//________________________________________________________________________
+AliEbyEPidRatioEffCont::AliEbyEPidRatioEffCont() :
+ AliEbyEPidRatioBase("EffCont", "EffCont"),
+ fLabelsRec(NULL),
+ fHnEff(NULL),
+ fHnCont(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioEffCont",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioEffCont::~AliEbyEPidRatioEffCont() {
+ // Destructor
+
+ for (Int_t ii = 0; ii < 2; ++ii) {
+ for (Int_t kk = 0; kk < 4; ++kk)
+ if (fLabelsRec[ii][kk]) delete[] fLabelsRec[ii][kk];
+ if (fLabelsRec[ii]) delete[] fLabelsRec[ii];
+ }
+ if (fLabelsRec) delete[] fLabelsRec;
+
+
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::Process() {
+ // -- Process event
+
+ // -- Setup (clean, create and fill) MC labels
+ FillMCLabels();
+
+ // -- Fill MC histograms for efficiency studies
+ FillMCEffHist();
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::Init() {
+ // -- Init eventwise
+
+ fLabelsRec = new Int_t**[2];
+ for (Int_t ii = 0 ; ii < 2; ++ii) {
+ fLabelsRec[ii] = new Int_t*[4];
+ for (Int_t kk = 0 ; kk < 4; ++kk)
+ fLabelsRec[ii][kk] = NULL;
+ }
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::CreateHistograms() {
+ // Copied from NetParticle class
+ Int_t binHnEff[20] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, AliEbyEPidRatioHelper::fgkfHistNBinsEta, // cent | etaMC
+ AliEbyEPidRatioHelper::fgkfHistNBinsRap, AliEbyEPidRatioHelper::fgkfHistNBinsPhi, // yMC | phiMC
+ AliEbyEPidRatioHelper::fgkfHistNBinsPt, AliEbyEPidRatioHelper::fgkfHistNBinsSign, // ptMC | signMC
+ 2, 2 , 2 , // findable | recStatus | pidStatus
+ AliEbyEPidRatioHelper::fgkfHistNBinsEta, AliEbyEPidRatioHelper::fgkfHistNBinsRap, // etaRec | yRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsPhi, AliEbyEPidRatioHelper::fgkfHistNBinsPt, // phiRec | ptRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsSign, // signRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsEta, AliEbyEPidRatioHelper::fgkfHistNBinsRap, // deltaEta | deltaY
+ 2*AliEbyEPidRatioHelper::fgkfHistNBinsPhi+1, 2*AliEbyEPidRatioHelper::fgkfHistNBinsPt+1, // deltaPhi | deltaPt
+ AliEbyEPidRatioHelper::fgkfHistNBinsSign+2, 4}; // deltaSign
+
+
+ Double_t minHnEff[20] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[0], AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[0], AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
+ -0.5, -0.5, -0.5,
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[0], AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[0], AliEbyEPidRatioHelper::fgkfHistRangePt[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[0], AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
+ -1.*AliEbyEPidRatioHelper::fgkfHistRangePhi[1], -1.*AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[0]-1., -0.5};
+
+ Double_t maxHnEff[20] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[1], AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[1], AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
+ 1.5, 1.5, 1.5,
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[1], AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[1], AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[1], AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[1], AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[1]+1., -3.5};
+
+ fHnEff = new THnSparseF("hnEff", "cent:etaMC:yMC:phiMC:ptMC:signMC:findable:recStatus:pidStatus:etaRec:yRec:phiRec:ptRec:signRec:deltaEta:deltaY:deltaPhi:deltaPt:deltaSign:PID",
+ 20, 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.2, 2.3]
+
+ 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("#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
+ fHnEff->GetAxis(11)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
+ fHnEff->GetAxis(12)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.2, 2.3]
+ fHnEff->GetAxis(13)->SetTitle("sign"); // -1 | 0 | +1
+
+ fHnEff->GetAxis(14)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9, 0.9]
+ fHnEff->GetAxis(15)->SetTitle("#it{y}_{MC}-#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
+ fHnEff->GetAxis(16)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ -2Pi , 2Pi]
+ fHnEff->GetAxis(17)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ -2.3, 2.3]
+ fHnEff->GetAxis(18)->SetTitle("sign_{MC}-sign_{Rec}"); // -2 | 0 | +2
+ fHnEff->GetAxis(19)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
+
+ fHelper->BinLogAxis(fHnEff, 4);
+ fHelper->BinLogAxis(fHnEff, 12);
+
+ /*
+ >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Copied from NetParticle
+
+ 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 THnSparse - Cont
+ // ------------------------------------------------------------------
+
+ Int_t binHnCont[18] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, AliEbyEPidRatioHelper::fgkfHistNBinsEta, // cent | etaMC
+ AliEbyEPidRatioHelper::fgkfHistNBinsRap, AliEbyEPidRatioHelper::fgkfHistNBinsPhi, // yMC | phiMC
+ AliEbyEPidRatioHelper::fgkfHistNBinsPt, AliEbyEPidRatioHelper::fgkfHistNBinsSign, // ptMC | signMC=contSign
+ 8, // contPart |
+ AliEbyEPidRatioHelper::fgkfHistNBinsEta, AliEbyEPidRatioHelper::fgkfHistNBinsRap, // etaRec | yRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsPhi, AliEbyEPidRatioHelper::fgkfHistNBinsPt, // phiRec | ptRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsSign, // signRec
+ AliEbyEPidRatioHelper::fgkfHistNBinsEta, AliEbyEPidRatioHelper::fgkfHistNBinsRap, // deltaEta | deltaY
+ 2*AliEbyEPidRatioHelper::fgkfHistNBinsPhi+1, 2*AliEbyEPidRatioHelper::fgkfHistNBinsPt+1, // deltaPhi | deltaPt
+ AliEbyEPidRatioHelper::fgkfHistNBinsSign+2, 4}; // deltaSign| pid 0-3
+
+ Double_t minHnCont[18] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[0], AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[0], AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
+ 0.5,
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[0], AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[0], AliEbyEPidRatioHelper::fgkfHistRangePt[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[0], AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
+ -1.*AliEbyEPidRatioHelper::fgkfHistRangePhi[1], -1.*AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[0]-1., -0.5};
+
+
+ Double_t maxHnCont[18] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[1], AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[1], AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
+ 8.5,
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[1], AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[1], AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[1], AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePhi[1], AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[1]+1., -3.5};
+
+ fHnCont = new THnSparseF("hnCont", "cent:etaMC:yMC:phiMC:ptMC:signMC:contPart:etaRec:yRec:phiRec:ptRec:signRec:deltaEta:deltaY:deltaPhi:deltaPt:deltaSign:PID",
+ 18, 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.2,2.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("#eta_{Rec}"); // eta [-0.9, 0.9]
+ fHnCont->GetAxis(8)->SetTitle("#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
+ fHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)"); // phi [ 0. , 2Pi]
+ fHnCont->GetAxis(10)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ 0.2, 2.3]
+ fHnCont->GetAxis(11)->SetTitle("sign"); // -1 | 0 | +1
+
+ fHnCont->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}"); // eta [-0.9, 0.9]
+ fHnCont->GetAxis(13)->SetTitle("#it{y}_{MC}-#it{y}_{Rec}"); // rapidity [-0.5, 0.5]
+ fHnCont->GetAxis(14)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); // phi [ -2Pi , 2Pi]
+ fHnCont->GetAxis(15)->SetTitle("#it{p}_{T,MC}-#it{p}_{T,Rec} (GeV/#it{c})"); // pt [ -2.3, 2.3]
+ fHnCont->GetAxis(16)->SetTitle("sign_{MC}-sign_{Rec}"); // -2 | 0 | +2
+ fHnCont->GetAxis(17)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}"); // 0 | 1 | 2 | 3
+
+ fHelper->BinLogAxis(fHnCont, 4);
+ fHelper->BinLogAxis(fHnCont, 10);
+
+ return;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioEffCont::Setup() {
+ // -- Setup eventwise
+
+ // -- Create label arrays
+ for(Int_t i = 0; i < 4; i++) {
+ fLabelsRec[0][i] = new Int_t[fNTracks];
+ if(!fLabelsRec[0][i]) {
+ AliError("Cannot create fLabelsRec[0]");
+ return -1;
+ }
+
+ fLabelsRec[1][i] = new Int_t[fNTracks];
+ if(!fLabelsRec[1][i]) {
+ AliError("Cannot create fLabelsRec[1] for PID");
+ return -1;
+ }
+
+ for(Int_t ii = 0; ii < fNTracks; ++ii) {
+ fLabelsRec[0][i][ii] = 0;
+ fLabelsRec[1][i][ii] = 0;
+ }
+ }
+ return 0;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::Reset() {
+ // -- Reset eventwise
+ for(Int_t i = 0; i < 4; i++) {
+ for (Int_t ii = 0; ii < 2 ; ++ii) {
+ if (fLabelsRec[ii][i])
+ delete[] fLabelsRec[ii][i];
+ fLabelsRec[ii][i] = NULL;
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::FillMCLabels() {
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+
+ // -- Track Loop
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
+
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted - ESD
+ if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
+ continue;
+
+ // -- Check if accepted - AOD
+ if (fAOD){
+ AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
+
+ if (!trackAOD) {
+ AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
+ continue;
+ }
+ if (!trackAOD->TestFilterBit(fAODtrackCutBit))
+ continue;
+
+ // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
+ if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
+ continue;
+ }
+
+ Int_t gPdgCode = 0;
+
+ Int_t iPid = 0;
+ Double_t pid[3];
+ if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion))) { iPid = 1; gPdgCode = 211;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon))) { iPid = 2; gPdgCode = 321;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){ iPid = 3; gPdgCode = 2212;}
+ else iPid = 0;
+
+
+ Double_t yP;
+ if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
+ continue;
+
+ if (!fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ Int_t label = TMath::Abs(track->GetLabel());
+
+ // -- Fill Label of all reconstructed
+ if(iPid != 0) fLabelsRec[0][0][idxTrack] = label;
+ fLabelsRec[0][iPid][idxTrack] = label;
+
+ // -- Fill Label of all reconstructed && recPid_TPC+TOF
+ if(iPid != 0) fLabelsRec[1][0][idxTrack] = label;
+ fLabelsRec[1][iPid][idxTrack] = label;
+
+ // -- Check for contamination and fill contamination THnSparse
+ CheckContTrack(track, iPid, gPdgCode);
+ if(iPid != 0) CheckContTrack(track, 0, 0);
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::CheckContTrack(AliVTrack *track, Int_t iPid, Int_t gPdgCode) {
+ Int_t label = TMath::Abs(track->GetLabel());
+ Float_t signRec = track->Charge();
+
+
+ AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
+ if (!particle)
+ return;
+
+ Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
+ if (iPid == 0) {
+ if (particle->PdgCode() == (signRec*gPdgCode))
+ if (isPhysicalPrimary)
+ return;
+ }
+ else {
+ if (isPhysicalPrimary)
+ return;
+ }
+
+ // -- Check if secondaries from material or weak decay
+ Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
+ Bool_t isSecondaryFromMaterial = (fESD) ? fStack->IsSecondaryFromMaterial(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
+
+ // -- Get PDG Charge of contaminating particle
+ Float_t signMC = 0.;
+ if (particle->Charge() == 0.) signMC = 0.;
+ else if (particle->Charge() < 0.) signMC = -1.;
+ else if (particle->Charge() > 0.) signMC = 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->PdgCode()) == 211) contPart = 1; // pion
+ else if (TMath::Abs(particle->PdgCode()) == 321) contPart = 2; // kaon
+ else if (TMath::Abs(particle->PdgCode()) == 2212) contPart = 3; // proton
+ else if (TMath::Abs(particle->PdgCode()) == 11) contPart = 4; // electron
+ else if (TMath::Abs(particle->PdgCode()) == 13) contPart = 5; // muon
+ else contPart = 6; // other
+ }
+
+ // -- Get Reconstructed y
+ // yRec = y for identified particles | yRec = eta for charged particles
+ Double_t yRec = 0.;
+ fHelper->IsTrackAcceptedRapidity(track, yRec, iPid);
+
+ Double_t deltaPhi = particle->Phi()-track->Phi();
+ if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
+ if (deltaPhi < 0)
+ deltaPhi += TMath::TwoPi();
+ else
+ deltaPhi -= TMath::TwoPi();
+ }
+
+ Double_t hnCont[18] = {fCentralityBin,
+ particle->Eta(),
+ particle->Y(),
+ particle->Phi(),
+ particle->Pt(),
+ signMC,
+ contPart,
+ track->Eta(),
+ yRec,
+ track->Phi(),
+ track->Pt(),
+ signRec,
+ particle->Eta()-track->Eta(),
+ particle->Y()-yRec,
+ deltaPhi,
+ particle->Pt()-track->Pt(),
+ signMC-signRec,
+ iPid};
+ fHnCont->Fill(hnCont);
+
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioEffCont::FillMCEffHist() {
+ // Fill efficiency THnSparse for ESDs
+
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Int_t nPart = (fESD) ? fStack->GetNprimary() : fArrayMC->GetEntriesFast();
+
+ for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+ AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : static_cast<AliVParticle*>(fArrayMC->At(idxMC));
+
+ // -- Check basic MC properties -> charged physical primary
+ if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
+ continue;
+
+ Int_t iPid = 0;
+ Int_t gPdgCode = 0;
+ if ( particle->PdgCode() == 211 ) { iPid = 1; gPdgCode = 211;}
+ if ( particle->PdgCode() == 321 ) { iPid = 1; gPdgCode = 321;}
+ if ( particle->PdgCode() == 2212 ) { iPid = 1; gPdgCode = 2212;}
+ else {iPid = 0; gPdgCode = 0;}
+
+ // -- Check if accepted in rapidity window -- for identified particles
+ Double_t yMC;
+ if (iPid != 0 && !fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
+ continue;
+
+ // -- Check if accepted in eta window -- for charged particles
+ if (iPid == 0 && TMath::Abs(particle->Eta()) > etaRange[1])
+ continue;
+
+ // -- Check if probeParticle / anti-probeParticle
+ // > skip check if PID is not required
+ if (iPid == 0 && TMath::Abs(particle->PdgCode()) != gPdgCode)
+ continue;
+
+ // -- Get sign of particle
+ Float_t signMC = (particle->PdgCode() < 0) ? -1. : 1.;
+
+ // -- Get if particle is findable --- not availible for AODs yet
+ Float_t findable = (fESD) ? Float_t(fHelper->IsParticleFindable(idxMC)) : 1.;
+
+ // -- 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.;
+ Double_t yRec = 0.;
+ Float_t signRec = 0.;
+
+ // -- Loop over all labels
+ for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
+ if (idxMC == fLabelsRec[0][iPid][idxRec]) {
+ recStatus = 1.;
+
+ if (idxMC == fLabelsRec[1][iPid][idxRec])
+ recPid = 1.;
+
+ AliVTrack *track = NULL;
+ if(fESD)
+ track = fESD->GetTrack(idxRec);
+ else if(fAOD)
+ track = fAOD->GetTrack(idxRec);
+
+ if (track) {
+ // if no track present (which should not happen)
+ // -> pt = 0. , which is not in the looked at range
+
+ // -- Get Reconstructed values
+ etaRec = track->Eta();
+ phiRec = track->Phi();
+ ptRec = track->Pt();
+ signRec = track->Charge();
+ fHelper->IsTrackAcceptedRapidity(track, yRec, iPid); // yRec = y for identified particles | yRec = eta for charged particles
+ }
+ break;
+ }
+ } // for (Int_t idxRec=0; idxRec < fNTracks; ++idxRec) {
+
+ Double_t deltaPhi = particle->Phi()-phiRec;
+ if (TMath::Abs(deltaPhi) > TMath::TwoPi()) {
+ if (deltaPhi < 0)
+ deltaPhi += TMath::TwoPi();
+ else
+ deltaPhi -= TMath::TwoPi();
+ }
+
+ if(iPid != 0) {
+ Double_t hnEff[20] = {fCentralityBin,
+ particle->Eta(),
+ particle->Y(),
+ particle->Phi(),
+ particle->Pt(),
+ signMC,
+ findable,
+ recStatus,
+ recPid,
+ etaRec,
+ yRec,
+ phiRec,
+ ptRec,
+ signRec,
+ particle->Eta()-etaRec,
+ particle->Y()-yRec,
+ deltaPhi,
+ particle->Pt()-ptRec,
+ signMC-signRec,
+ 0};
+ fHnEff->Fill(hnEff);
+ }
+
+ Double_t hnEff[20] = {fCentralityBin,
+ particle->Eta(),
+ particle->Y(),
+ particle->Phi(),
+ particle->Pt(),
+ signMC,
+ findable,
+ recStatus,
+ recPid,
+ etaRec,
+ yRec,
+ phiRec,
+ ptRec,
+ signRec,
+ particle->Eta()-etaRec,
+ particle->Y()-yRec,
+ deltaPhi,
+ particle->Pt()-ptRec,
+ signMC-signRec,
+ iPid};
+ fHnEff->Fill(hnEff);
+
+
+ } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+
+ return;
+}
--- /dev/null
+#ifndef ALIEBYEPIDRSTIOEFFCONT_H
+#define ALIEBYEPIDRSTIOEFFCONT_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+class AliVTrack;
+
+#include "THnSparse.h"
+
+#include "AliEbyEPidRatioBase.h"
+
+class AliEbyEPidRatioEffCont: public AliEbyEPidRatioBase {
+
+ public:
+
+ AliEbyEPidRatioEffCont();
+ virtual ~AliEbyEPidRatioEffCont();
+ virtual void Process();
+ THnSparseF* GetHnEff() {return fHnEff;}
+ THnSparseF* GetHnCont() {return fHnCont;}
+
+ private:
+
+ AliEbyEPidRatioEffCont(const AliEbyEPidRatioEffCont&); // not implemented
+ AliEbyEPidRatioEffCont& operator=(const AliEbyEPidRatioEffCont&); // not implemented
+
+
+ virtual void Init();
+ virtual void CreateHistograms();
+ virtual void Reset();
+ virtual Int_t Setup();
+ void FillMCLabels();
+ void FillMCEffHist();
+ void CheckContTrack(AliVTrack* track, Int_t iPid, Int_t gPdgCode);
+
+ Int_t ***fLabelsRec; //! 2x nTracks large array with labels for MC particles
+ THnSparseF *fHnEff; // THnSparseF efficiency
+ THnSparseF *fHnCont; // THnSparseF contamination
+
+ ClassDef(AliEbyEPidRatioEffCont, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#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 "AliAODInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+#include "AliCentrality.h"
+#include "AliTracker.h"
+
+#include "AliEbyEPidRatioHelper.h"
+
+using namespace std;
+
+
+ClassImp(AliEbyEPidRatioHelper)
+//________________________________________________________________________
+AliEbyEPidRatioHelper::AliEbyEPidRatioHelper() :
+ fModeDistCreation(0),
+
+ fInputEventHandler(NULL),
+ fPIDResponse(NULL),
+ fESD(NULL),
+ fESDTrackCuts(NULL),
+ fAOD(NULL),
+ fAODtrackCutBit(1024),
+ fIsMC(kFALSE),
+ fMCEvent(NULL),
+ fStack(NULL),
+
+ fCentralityBin(-1),
+ fCentralityPercentile(-1.),
+
+ fCentralityBinMax(7),
+ fVertexZMax(10.),
+ fRapidityMax(0.5),
+ fPhiMin(0.),
+ fPhiMax(TMath::TwoPi()),
+ fMinTrackLengthMC(70.),
+ fNSigmaMaxCdd(3.),
+ fNSigmaMaxCzz(3.),
+
+ fPIDStrategy(0),
+ fNSigmaMaxITS(4.),
+ fNSigmaMaxTPC(4.),
+ fNSigmaMaxTPClow(3.),
+ fNSigmaMaxTOF(4.),
+ fMinPtForTOFRequired(0.69),
+ fMaxPtForTPClow(0.69),
+
+ fHEventStat0(NULL),
+ fHEventStat1(NULL),
+ fHEventStatMax(6),
+
+ fHTriggerStat(NULL),
+ fNTriggers(5),
+
+ fHCentralityStat(NULL),
+ fNCentralityBins(10),
+
+ fRandom(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioHelper",10);
+}
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthRap = 0.075;
+const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthPt = 0.3; // 0.08 // 300 MeV // was 80 MeV
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangeCent[] = {-0.5, 10.5};
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsCent = 11 ;
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangeEta[] = {-0.9, 0.9};
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsEta = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeEta[1] -
+ AliEbyEPidRatioHelper::fgkfHistRangeEta[0]) /
+ AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangeRap[] = {-0.5, 0.5};
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsRap = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeRap[1] - AliEbyEPidRatioHelper::fgkfHistRangeRap[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangePhi[] = {0.0, TMath::TwoPi()};
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsPhi = 42 ;
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangePt[] = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsPt = Int_t((AliEbyEPidRatioHelper::fgkfHistRangePt[1] - AliEbyEPidRatioHelper::fgkfHistRangePt[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthPt);
+
+const Float_t AliEbyEPidRatioHelper::fgkfHistRangeSign[] = {-1.5, 1.5};
+const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsSign = 3;
+
+const Char_t* AliEbyEPidRatioHelper::fgkEventNames[] = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
+const Char_t* AliEbyEPidRatioHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
+const Char_t* AliEbyEPidRatioHelper::fgkTriggerNames[] = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
+const Char_t* AliEbyEPidRatioHelper::fgkCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%","50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
+
+const Char_t* AliEbyEPidRatioHelper::fgkPidName[4] = {"Nch","Npi","Nka","Npr"};
+const Char_t* AliEbyEPidRatioHelper::fgkPidLatex[4][2] = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
+const Char_t* AliEbyEPidRatioHelper::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
+
+
+
+//________________________________________________________________________
+AliEbyEPidRatioHelper::~AliEbyEPidRatioHelper() {
+ // Destructor
+
+ if (fRandom)
+ delete fRandom;
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioHelper::SetPhiRange(Float_t f1, Float_t f2) {
+ // -- Set phi range and adopt to phi-histogram
+
+ fPhiMin = f1;
+ fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
+
+ Float_t phiMin = fPhiMin;
+ Float_t phiMax = fPhiMax;
+
+ // -- Update Ranges
+ Float_t binWidth = (AliEbyEPidRatioHelper::fgkfHistRangePhi[1] - AliEbyEPidRatioHelper::fgkfHistRangePhi[0]) /
+ Float_t(AliEbyEPidRatioHelper::fgkfHistNBinsPhi);
+
+ Float_t lowEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0] - binWidth;
+ Float_t highEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0];
+
+ for (Int_t ii = 1; ii <= AliEbyEPidRatioHelper::fgkfHistNBinsPhi; ++ii) {
+ lowEdge += binWidth;
+ highEdge += binWidth;
+
+ if (phiMin >= lowEdge && phiMin < highEdge )
+ phiMin = lowEdge;
+ if (phiMax > lowEdge && phiMax <= highEdge )
+ phiMax = highEdge;
+ }
+
+ printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
+ fPhiMin = phiMin;
+ fPhiMax = phiMax;
+}
+
+
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Int_t trackCutBit, Int_t modeDistCreation) {
+ // -- Initialize helper
+
+ Int_t iResult = 0;
+
+ // -- ESD track cuts
+ fESDTrackCuts = cuts;
+
+ // -- Is MC
+ fIsMC = isMC;
+
+ // -- AOD track filter bit
+ fAODtrackCutBit = trackCutBit;
+
+ // -- mode Distribution creation
+ fModeDistCreation = modeDistCreation;
+
+ // -- Setup event cut statistics
+ InitializeEventStats();
+
+ // -- Setup trigger statistics
+ InitializeTriggerStats();
+
+ // -- Setup centrality statistics
+ InitializeCentralityStats();
+
+ // -- PRINT PID Strategy
+ // 0 : TPC(TPClow+TPCHigh)
+ // 1 : ITS
+ // 2 : TOF
+ // 3 : ITS+TPC(TPClow+TPCHigh)
+ // 4 : TPC(TPClow+TPCHigh)+TOF
+ // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
+ // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
+ // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
+ // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
+ printf(">>>> PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);
+
+ // -- Initialize random number generator
+ fRandom = new TRandom3();
+ fRandom->SetSeed();
+
+ return iResult;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
+ // -- Setup Event
+
+ // -- Get ESD objects
+ if(esdHandler){
+ fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
+ fESD = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
+ if (!fESD) {
+ AliError("ESD event handler not available");
+ return -1;
+ }
+ }
+
+ // -- Get AOD objects
+ else if(aodHandler){
+ fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
+ fAOD = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
+ if (!fAOD) {
+ AliError("AOD event handler not available");
+ return -1;
+ }
+ }
+
+ // -- Get Common objects
+ fPIDResponse = fInputEventHandler->GetPIDResponse();
+
+ // -- 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 = NULL;
+
+ if(esdHandler)
+ centrality = fESD->GetCentrality();
+ else if(aodHandler)
+ centrality = fAOD->GetHeader()->GetCentralityP();
+
+ if (!centrality) {
+ AliError("Centrality not available");
+ return -1;
+ }
+
+ 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");
+
+
+ return 0;
+}
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::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 ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
+ if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
+ if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
+ if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
+ if ((fInputEventHandler->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 AliEbyEPidRatioHelper::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 = NULL;
+ const AliAODVertex* vtxAOD = NULL;
+ if (fESD){
+ vtxESD = fESD->GetPrimaryVertexTracks();
+ if (!vtxESD)
+ aEventCuts[iCut] = 1;
+ }
+ else if (fAOD){
+ vtxAOD = fAOD->GetPrimaryVertex();
+ if (!vtxAOD)
+ aEventCuts[iCut] = 1;
+ }
+
+ // -- 3 - Vertex z outside cut window
+ ++iCut;
+ if (vtxESD){
+ if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
+ aEventCuts[iCut] = 1;
+ }
+ else if(vtxAOD){
+ if(TMath::Abs(vtxAOD->GetZ()) > 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;
+
+ //cout << isRejected << endl;
+
+ return isRejected;
+}
+
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
+ // -- Check if MC particle is accepted for basic parameters
+
+ if (!particle)
+ return kFALSE;
+
+ // -- check if charged
+ if (particle->Charge() == 0.0)
+ return kFALSE;
+
+ // -- check if physical primary - ESD
+ if (fESD) {
+ if(!fStack->IsPhysicalPrimary(idxMC))
+ return kFALSE;
+ }
+ // -- check if physical primary - AOD
+ else {
+ if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
+ return kFALSE;
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
+ // -- Check if MC particle is accepted for basic parameters
+
+ if (!particle)
+ return kFALSE;
+
+ // -- check if charged
+ if (particle->Charge() != 0.0)
+ return kFALSE;
+
+ // -- check if physical primary - ESD
+ if (fESD) {
+ if(!fStack->IsPhysicalPrimary(idxMC))
+ return kFALSE;
+ }
+ // -- check if physical primary - AOD
+ else {
+ if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
+ return kFALSE;
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP, Int_t gCurPid) {
+ // -- Check if particle is accepted
+ // > in rapidity
+ // > if no pid : return kTRUE, yP = eta
+ // > return 0 if not accepted
+
+ if (gCurPid == 0) {
+ yP = particle->Eta();
+ return kTRUE;
+ }
+
+ Double_t mP;
+ if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
+ if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
+ if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
+
+ // -- 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 AliEbyEPidRatioHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
+ // -- Check if particle is accepted
+ // > in phi
+ // > return 0 if not accepted
+
+ if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
+ return kTRUE;
+ else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::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);
+}
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
+ // -- Check if track is accepted
+ // > for basic parameters
+
+ if (!track)
+ return kFALSE;
+
+ if (track->Charge() == 0)
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP, Int_t gCurPid) {
+ if (gCurPid == 0) {
+ yP = track->Eta();
+ return kTRUE;
+ }
+
+ Double_t mP;
+ if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
+ if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
+ if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
+
+ // -- Calculate rapidities and kinematics
+ Double_t pvec[3];
+ track->GetPxPyPz(pvec);
+
+ Double_t p = track->P();
+ 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 AliEbyEPidRatioHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
+ Bool_t isAccepted = kTRUE;
+
+ if (!fESD)
+ return isAccepted;
+
+ AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
+ if (!track)
+ return kFALSE;
+
+ // -- 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 AliEbyEPidRatioHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid, AliPID::EParticleType gCurPid) {
+
+ Bool_t isAcceptedITS = kFALSE;
+ Bool_t isAcceptedTPC = kFALSE;
+ Bool_t isAcceptedTPClow = kFALSE;
+ Bool_t isAcceptedTOF = kFALSE;
+ Bool_t isAccepted = kFALSE;
+
+ // -- In case not PID is used
+ if (gCurPid == 0) {
+ pid[0] = 10.;
+ pid[1] = 10.;
+ pid[2] = 10.;
+ return kTRUE;
+ }
+
+ if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, gCurPid, pid[0]) == AliPIDResponse::kDetPidOk) {
+ if (TMath::Abs(pid[0]) < fNSigmaMaxITS)
+ isAcceptedITS = kTRUE;
+ }
+
+ if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, gCurPid, pid[1]) == AliPIDResponse::kDetPidOk) {
+ if (TMath::Abs(pid[1]) < fNSigmaMaxTPC)
+ isAcceptedTPC = kTRUE;
+ if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow)
+ isAcceptedTPClow = kTRUE;
+ if (track->Pt() < fMaxPtForTPClow)
+ isAcceptedTPC = isAcceptedTPClow;
+ }
+
+ Bool_t hasPIDTOF = kFALSE;
+ if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, gCurPid, pid[2]) == AliPIDResponse::kDetPidOk) {
+ hasPIDTOF = kTRUE;
+ if (TMath::Abs(pid[2]) < fNSigmaMaxTOF)
+ isAcceptedTOF = kTRUE;
+ }
+
+ // -- Check TOF missmatch for MC
+
+ //if (ESD)
+ if (fIsMC && isAcceptedTOF) {
+ Int_t tofLabel[3];
+ // AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
+ // TODO add code for AOD
+
+ (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
+
+ Bool_t hasMatchTOF = kTRUE;
+ if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0)
+ hasMatchTOF = kFALSE;
+
+ TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
+ if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
+ hasMatchTOF = kTRUE;
+
+ isAcceptedTOF = hasMatchTOF;
+ }
+
+ // 0 : TPC(TPClow+TPCHigh)
+ // 1 : ITS
+ // 2 : TOF
+ // 3 : ITS+TPC(TPClow+TPCHigh)
+ // 4 : TPC(TPClow+TPCHigh)+TOF
+ // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
+ // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
+ // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
+ // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
+ if (fPIDStrategy == 0) { // TPC PID
+ isAccepted = isAcceptedTPC;
+ }
+ else if (fPIDStrategy == 1) { // ITS PID
+ isAccepted = isAcceptedITS;
+ }
+ else if (fPIDStrategy == 2) { // TOF PID
+ isAccepted = isAcceptedTOF;
+ }
+ else if (fPIDStrategy == 3) { // TPC+ITS PID
+ isAccepted = isAcceptedTPC && isAcceptedITS;
+ }
+ else if (fPIDStrategy == 4) { // TPC+TOF PID
+ isAccepted = isAcceptedTPC && isAcceptedTOF;
+ }
+ else if (fPIDStrategy == 5) { // TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
+ if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
+ isAcceptedTOF = kTRUE;
+ isAccepted = isAcceptedTPC && isAcceptedTOF;
+ }
+ else if (fPIDStrategy == 6) { // ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
+ isAccepted = isAcceptedTPC && isAcceptedITS;
+ if (hasPIDTOF)
+ isAccepted = isAccepted && isAcceptedTOF;
+ }
+ else if (fPIDStrategy == 7) { // ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
+ if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
+ isAcceptedTOF = kTRUE;
+ isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
+ }
+ else if (fPIDStrategy == 8) { // ITS+TPC+TOF PID
+ isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
+ }
+
+ return isAccepted;
+}
+
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPhi(AliVTrack *track) {
+ // -- Check if track is accepted
+ // > in phi
+ // > return 0 if not accepted
+
+ if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
+ return kTRUE;
+ else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
+ AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
+
+ // -- Make logarithmic binning
+ TAxis *axis = hn->GetAxis(axisNumber);
+ Int_t nBins = axis->GetNbins();
+
+ Double_t from = axis->GetXmin();
+ Double_t to = axis->GetXmax();
+ Double_t *newBins = new Double_t[nBins + 1];
+
+ newBins[0] = from;
+ Double_t factor = TMath::Power(to/from, 1./nBins);
+
+ for (int ii = 1; ii <= nBins; ii++)
+ newBins[ii] = factor * newBins[ii-1];
+
+ axis->Set(nBins, newBins);
+
+ delete [] newBins;
+
+ // -- Update Ranges
+ // ------------------
+ Float_t ptRange[2];
+ Float_t oldPtRange[2];
+ esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+ esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
+
+ Float_t minPtForTOFRequired = fMinPtForTOFRequired;
+ Float_t maxPtForTPClow = fMaxPtForTPClow;
+
+ // -- Update minPt Cut
+ Int_t bin = axis->FindBin(ptRange[0]+10e-6);
+ ptRange[0] = axis->GetBinLowEdge(bin);
+
+ // -- Update maxPt Cut
+ bin = axis->FindBin(ptRange[1]-10e-6);
+ ptRange[1] = axis->GetBinUpEdge(bin);
+
+ if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
+ printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
+ esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
+ }
+
+ // -- Update MinPtForTOFRequired
+ bin = axis->FindBin(minPtForTOFRequired-10e-6);
+ minPtForTOFRequired = axis->GetBinUpEdge(bin);
+
+ if (minPtForTOFRequired != fMinPtForTOFRequired) {
+ printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
+ fMinPtForTOFRequired = minPtForTOFRequired;
+ }
+
+ // -- Update MaxPtForTPClow
+ bin = axis->FindBin(maxPtForTPClow-10e-6);
+ maxPtForTPClow = axis->GetBinUpEdge(bin);
+
+ if (maxPtForTPClow != fMaxPtForTPClow) {
+ printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
+ fMaxPtForTPClow = maxPtForTPClow;
+ }
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioHelper::InitializeEventStats() {
+ // -- Initialize event statistics histograms
+
+ fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
+ fHEventStat1 = new TH1F("hEventStat1","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, AliEbyEPidRatioHelper::fgkEventNames[ii]);
+ fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
+ }
+
+ fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
+ fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioHelper::InitializeTriggerStats() {
+ // -- Initialize trigger statistics histograms
+
+ fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
+
+ for ( Int_t ii=0; ii < fNTriggers; ii++ )
+ fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkTriggerNames[ii]);
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioHelper::InitializeCentralityStats() {
+ // -- Initialize trigger statistics histograms
+
+ fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events",
+ fNCentralityBins,-0.5,fNCentralityBins-0.5);
+
+ for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
+ fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkCentralityNames[ii]);
+}
+//________________________________________________________________________
+Bool_t AliEbyEPidRatioHelper::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
+#ifndef ALIEBYEPIDRATIOHELPER_H
+#define ALIEBYEPIDRATIOHELPER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "THnBase.h"
+#include "THn.h"
+#include "TH1F.h"
+#include "TF1.h"
+#include "TProfile2D.h"
+#include "TRandom3.h"
+
+class AliESDtrack;
+class AliMCEvent;
+class AliStack;
+class AliPIDResponse;
+class AliESDtrackCuts;
+class AliInputEventHandler;
+class AliESDInputHandler;
+class AliAODInputHandler;
+class AliAODEvent;
+class AliAODTrack;
+class AliAODMCParticle;
+class AliMCParticle;
+
+class AliEbyEPidRatioHelper : public TNamed {
+
+ public:
+
+ AliEbyEPidRatioHelper();
+ virtual ~AliEbyEPidRatioHelper();
+
+ 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 SetPIDStrategy(Int_t i) {fPIDStrategy = i;}
+ void SetNSigmaMaxITS(Float_t f) {fNSigmaMaxITS = f;}
+ void SetNSigmaMaxTPC(Float_t f) {fNSigmaMaxTPC = f;}
+ void SetNSigmaMaxTPClow(Float_t f) {fNSigmaMaxTPClow = f;}
+ void SetNSigmaMaxTOF(Float_t f) {fNSigmaMaxTOF = f;}
+ void SetMinPtForTOFRequired(Float_t f) {fMinPtForTOFRequired = f;}
+ void SetMaxPtForTPClow(Float_t f) {fMaxPtForTPClow = f;}
+
+ void SetPhiRange(Float_t f1, Float_t f2);
+
+ 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;}
+
+ Float_t GetMinPtForTOFRequired() {return fMinPtForTOFRequired;}
+ Float_t GetMaxPtForTPClow() {return fMaxPtForTPClow;}
+ Float_t GetRapidityMax() {return fRapidityMax;}
+ Float_t GetPhiMin() {return fPhiMin;}
+ Float_t GetPhiMax() {return fPhiMax;}
+ AliESDtrackCuts* GetESDTrackCuts() {return fESDTrackCuts;}
+ Bool_t GetIsMC() {return fIsMC;}
+ Int_t GetAODtrackCutBit() {return fAODtrackCutBit;}
+ AliInputEventHandler* GetInputEventHandler() {return fInputEventHandler;}
+ AliMCEvent* GetMCEvent() {return fMCEvent;}
+
+ /** Initialize Helper */
+ Int_t Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Int_t trackCutBit, Int_t modeDistCreation);
+
+ /** Setup Event */
+ Int_t SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent);
+
+ /** Check if event is triggred */
+ Bool_t IsEventTriggered();
+
+ /** Fill event cut statistics */
+ Bool_t IsEventRejected();
+
+ /** Check if charged MC particle is accepted for basic parameters */
+ Bool_t IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC);
+
+ /** Check if neutral MC particle is accepted for basic parameters */
+ Bool_t IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC);
+
+ /** Check if MC particle is accepted for Rapidity */
+ Bool_t IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP, Int_t gCurPid);
+
+ /** Check if MC particle is accepted for Phi */
+ Bool_t IsParticleAcceptedPhi(AliVParticle *particle);
+
+ /** Check if MC particle is findable tracks */
+ Bool_t IsParticleFindable(Int_t label);
+
+ /** Check if track is accepted for basic parameters */
+ Bool_t IsTrackAcceptedBasicCharged(AliVTrack *track);
+
+ /** Check if track is accepted for Rapidity */
+ Bool_t IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP, Int_t gCurPid);
+
+ /** Check if track is accepted for DCA */
+ Bool_t IsTrackAcceptedDCA(AliVTrack *track);
+
+ /** Check if track is accepted for PID */
+ Bool_t IsTrackAcceptedPID(AliVTrack *track, Double_t *pid, AliPID::EParticleType gCurPid);
+
+ /** Check if trackis accepted for Phi */
+ Bool_t IsTrackAcceptedPhi(AliVTrack *track);
+
+ /** Method for the correct logarithmic binning of histograms
+ * and Update MinPtForTOFRequired, using the pT log-scale
+ */
+ void BinLogAxis(const THnBase *h, Int_t axisNumber, AliESDtrackCuts* cuts = NULL);
+
+ static const Float_t fgkfHistBinWitdthRap; // Histogram std bin width for rapidity/eta
+ static const Float_t fgkfHistBinWitdthPt; // Histogram std bin width for pt
+
+ static const Float_t fgkfHistRangeCent[]; // Histogram range for centrality
+ static const Int_t fgkfHistNBinsCent; // Histogram N bins for centrality
+ static const Float_t fgkfHistRangeEta[]; // Histogram range for eta
+ static const Int_t fgkfHistNBinsEta; // Histogram N bins for eta
+ static const Float_t fgkfHistRangeRap[]; // Histogram range for rapidity
+ static const Int_t fgkfHistNBinsRap; // Histogram N bins for rapidity
+ static const Float_t fgkfHistRangePhi[]; // Histogram range for phi
+ static const Int_t fgkfHistNBinsPhi; // Histogram N bins for phi
+ static const Float_t fgkfHistRangePt[]; // Histogram range for pt
+ static const Int_t fgkfHistNBinsPt; // Histogram N bins for pt
+ static const Float_t fgkfHistRangeSign[]; // Histogram range for sign
+ static const Int_t fgkfHistNBinsSign; // Histogram N bins for sign
+
+ static const Char_t* fgkEventNames[]; // Event names
+ static const Char_t* fgkCentralityMaxNames[]; // Centrality names
+ static const Char_t* fgkTriggerNames[]; // Trigger names
+ static const Char_t* fgkCentralityNames[]; // Centrality names
+ static const Char_t* fgkPidName[4];
+ static const Char_t* fgkPidLatex[4][2];
+ static const Char_t* fgkPidTitles[4][2];
+
+ private:
+
+ AliEbyEPidRatioHelper(const AliEbyEPidRatioHelper&); // not implemented
+ AliEbyEPidRatioHelper& operator=(const AliEbyEPidRatioHelper&); // not implemented
+
+ void InitializeEventStats();
+ void InitializeTriggerStats();
+ void InitializeCentralityStats();
+ Bool_t FillEventStats(Int_t *aEventCuts);
+
+ Int_t fModeDistCreation; // Dist creation mode : 1 = on | 0 = off
+ AliInputEventHandler *fInputEventHandler; //! Ptr to input event handler (ESD or AOD)
+ AliPIDResponse *fPIDResponse; //! Ptr to PID response Object
+ AliESDEvent *fESD; //! Ptr to ESD event
+ AliESDtrackCuts *fESDTrackCuts; //! Ptr to ESD cuts
+ AliAODEvent *fAOD; //! Ptr to AOD event
+ Int_t fAODtrackCutBit; // Track filter bit for AOD tracks
+ Bool_t fIsMC; // Is MC 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 fPhiMin; // Phi min cut
+ Float_t fPhiMax; // Phi max 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.
+
+ Int_t fPIDStrategy; // PID Strategy to be used
+ Float_t fNSigmaMaxITS; // N Sigma for ITS PID
+ Float_t fNSigmaMaxTPC; // N Sigma for TPC PID
+ Float_t fNSigmaMaxTPClow; // N Sigma for TPC PID lower part
+ Float_t fNSigmaMaxTOF; // N Sigma for TOF PID
+ Float_t fMinPtForTOFRequired; // Min pt from where TOF is required
+ Float_t fMaxPtForTPClow; // Max pt until TPClow is used
+ 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
+ TRandom3 *fRandom; // Random generator
+
+ ClassDef(AliEbyEPidRatioHelper, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "TMath.h"
+#include "TAxis.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TH2D.h"
+#include "TH3D.h"
+
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
+#include "AliESDEvent.h"
+#include "AliESDtrackCuts.h"
+
+#include "AliCentrality.h"
+#include "AliTracker.h"
+
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+
+#include "AliEbyEPidRatioPhy.h"
+
+using namespace std;
+
+ClassImp(AliEbyEPidRatioPhy)
+//________________________________________________________________________
+AliEbyEPidRatioPhy::AliEbyEPidRatioPhy() :
+ AliEbyEPidRatioBase("Dist", "Dist"),
+ fOutList(NULL),
+ fOrder(8),
+ fNNp(6),
+ fNp(NULL),
+ fNMCNp(7),
+ fMCNp(NULL),
+ fRedFactp(NULL) {
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioPhy",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioPhy::~AliEbyEPidRatioPhy() {
+ // Destructor
+
+ for (Int_t ii = 0; ii < fNNp; ++ii) {
+ for (Int_t kk = 0; kk < 2; ++kk)
+ if (fNp[ii][kk]) delete[] fNp[ii][kk];
+ if (fNp[ii]) delete[] fNp[ii];
+ }
+ if (fNp) delete[] fNp;
+
+ for (Int_t ii = 0; ii < fNMCNp; ++ii) {
+ for (Int_t kk = 0; kk < 2; ++kk)
+ if (fMCNp[ii][kk]) delete[] fMCNp[ii][kk];
+ if (fMCNp[ii]) delete[] fMCNp[ii];
+ }
+ if (fMCNp) delete[] fMCNp;
+
+ for (Int_t ii = 0; ii <= fOrder; ++ii)
+ if (fRedFactp[ii]) delete[] fRedFactp[ii];
+ if (fRedFactp) delete[] fRedFactp;
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::Process() {
+ ProcessTracks();
+ if (fIsMC)
+ ProcessParticles();
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::Init() {
+ fNp = new Int_t**[fNNp];
+ for (Int_t ii = 0 ; ii < fNNp; ++ii) {
+ fNp[ii] = new Int_t*[4];
+ for (Int_t kk = 0 ; kk < 4; ++kk)
+ fNp[ii][kk] = new Int_t[2];
+ }
+
+ fMCNp = new Int_t**[fNMCNp];
+ for (Int_t ii = 0 ; ii < fNMCNp; ++ii) {
+ fMCNp[ii] = new Int_t*[4];
+ for (Int_t kk = 0 ; kk < 4; ++kk)
+ fMCNp[ii][kk] = new Int_t[2];
+ }
+
+ fRedFactp = new Double_t*[fOrder+1];
+ for (Int_t ii = 0 ; ii <= fOrder; ++ii)
+ fRedFactp[ii] = new Double_t[2];
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::Reset() {
+ for (Int_t ii = 0; ii < fNNp; ++ii)
+ for (Int_t kk = 0; kk < 4; ++kk)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fNp[ii][kk][jj] = 0;
+ for (Int_t ii = 0; ii < fNMCNp; ++ii)
+ for (Int_t kk = 0; kk < 4; ++kk)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fMCNp[ii][kk][jj] = 0;
+ for (Int_t ii = 0; ii <= fOrder; ++ii)
+ for (Int_t jj = 0; jj < 2; ++jj)
+ fRedFactp[ii][jj] = 1.;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::CreateHistograms() {
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+ TString sTitle("");
+ AddHistSetCent("Phy", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
+ AddHistSetCent("PhyTPC", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
+ AddHistSetCent("PhyTOF", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
+
+#if USE_PHI
+ AddHistSetCent("Phyphi", Form("%s,#it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+ AddHistSetCent("PhyTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+ AddHistSetCent("PhyTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
+
+ if (fIsMC) {
+ TString sMCTitle("");
+
+ AddHistSetCent("MC", Form("%s", sTitle.Data()));
+ AddHistSetCent("MCpt", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1]));
+ AddHistSetCent("MCTPC", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
+ AddHistSetCent("MCTOF", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
+
+#if USE_PHI
+ AddHistSetCent("MCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+ AddHistSetCent("MCTPCphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+ AddHistSetCent("MCTOFphi",Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]",sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
+
+ }
+
+ return;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioPhy::ProcessTracks() {
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+ // -- Track Loop
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+ AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+ // -- Check if accepted - ESD
+ if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
+ continue;
+ // -- Check if accepted - AOD
+ if (fAOD){
+ AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
+
+ if (!trackAOD) {
+ AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
+ continue;
+ }
+ if (!trackAOD->TestFilterBit(fAODtrackCutBit))
+ continue;
+ if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
+ continue;
+ }
+
+ if (!fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ Int_t iPid = 0;
+ Double_t pid[3];
+ if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion))) iPid = 1;
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon))) iPid = 2;
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))) iPid = 3;
+ else iPid = 0;
+
+ Double_t yP;
+ if (iPid != 0 && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
+ continue;
+
+ Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
+ // -- in pt Range
+ fNp[0][0][idxPart] += 1;
+ if(iPid != 0) fNp[0][iPid][idxPart] += 1;
+ // -- in TPC pt Range
+ if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
+ fNp[1][0][idxPart] += 1;
+ if(iPid != 0) fNp[1][iPid][idxPart] += 1;
+ }
+ // -- in TPC+TOF pt Range
+ if (track->Pt() > fHelper->GetMinPtForTOFRequired()){
+ fNp[2][0][idxPart] += 1;
+ if(iPid != 0) fNp[2][iPid][idxPart] += 1;
+ }
+
+#if USE_PHI
+ if(!fHelper->IsTrackAcceptedPhi(track))
+ continue;
+
+ // -- in pt Range
+ fNp[3][iPid][idxPart] += 1;
+ if(iPid != 0) fNp[3][0][idxPart] += 1;
+ // -- in TPC pt Range
+ if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) {
+ fNp[4][0][idxPart] += 1;
+ if(iPid != 0)fNp[4][iPid][idxPart] += 1;
+ }
+ // -- in TPC+TOF pt Range
+ if (track->Pt() > fHelper->GetMinPtForTOFRequired()) {
+ fNp[5][0][idxPart] += 1;
+ if(iPid != 0) fNp[5][iPid][idxPart] += 1;
+ }
+#endif
+ } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
+
+ FillHistSetCent("Phy", 0, kFALSE);
+ FillHistSetCent("PhyTPC", 1, kFALSE);
+ FillHistSetCent("PhyTOF", 2, kFALSE);
+
+#if USE_PHI
+ FillHistSetCent("Phyphi", 3, kFALSE);
+ FillHistSetCent("PhyTPCphi", 4, kFALSE);
+ FillHistSetCent("PhyTOFphi", 5, kFALSE);
+
+ #endif
+ /* Printf("<<<<<<<<<< Inside Loop >>>>>>>>>>");
+ for (Int_t id = 0; id < 6; ++id) {
+ printf(" == %2d == ", id);
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ printf("%7d %7d ", fNp[id][iPid][0] , fNp[id][iPid][1]);
+ }
+ Printf("");
+ } */
+
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioPhy::ProcessParticles() {
+ // -- Process primary particles from the stack and fill histograms
+
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+
+ for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
+ AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : NULL;
+
+ if (!particle)
+ continue;
+ if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
+ continue;
+
+ Int_t iPid = 0;
+ if (TMath::Abs(particle->PdgCode()) == 211) iPid = 1; // pion
+ else if (TMath::Abs(particle->PdgCode()) == 321) iPid = 2; // kaon
+ else if (TMath::Abs(particle->PdgCode()) == 2212) iPid = 3; // proton
+ else iPid = 0;
+
+ Double_t yMC;
+ if ((iPid != 0) && !fHelper->IsParticleAcceptedRapidity(particle, yMC, iPid))
+ continue;
+
+ // -- Check eta window -- for charged particles
+ if ((iPid == 0) && TMath::Abs(particle->Eta()) > etaRange[1])
+ continue;
+
+ Int_t idxPart = (particle->PdgCode() < 0) ? 0 : 1;
+
+ fMCNp[0][0][idxPart] += 1.;
+ if(iPid != 0) fMCNp[0][iPid][idxPart] += 1.;
+
+ // -- Check main pt window
+ if (!(particle->Pt() > ptRange[0] && particle->Pt() <= ptRange[1]))
+ continue;
+
+ // -- in pt Range
+ fMCNp[1][0][idxPart] += 1.;
+ if(iPid != 0)fMCNp[1][iPid][idxPart] += 1.;
+
+ // -- in TPC pt Range
+ if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
+ fMCNp[2][0][idxPart] += 1;
+ if(iPid != 0)fMCNp[2][iPid][idxPart] += 1;
+ }
+ // -- in TPC+TOF pt Range
+ if (particle->Pt() > fHelper->GetMinPtForTOFRequired()) {
+ fMCNp[3][0][idxPart] += 1;
+ if(iPid != 0) fMCNp[3][iPid][idxPart] += 1;
+ }
+
+#if USE_PHI
+ if(!fHelper->IsParticleAcceptedPhi(particle))
+ continue;
+ // idxPhi = 1;
+
+ // -- in pt Range
+ fMCNp[4][0][idxPart] += 1;
+ if(iPid != 0)fMCNp[4][iPid][idxPart] += 1;
+
+ // -- in TPC pt Range
+ if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) {
+ fMCNp[5][0][idxPart] += 1;
+ if(iPid != 0)fMCNp[5][iPid][idxPart] += 1;
+ }
+
+ // -- in TPC+TOF pt Range
+ if (particle->Pt() > fHelper->GetMinPtForTOFRequired()){
+ fMCNp[6][0][idxPart] += 1;
+ if(iPid != 0)fMCNp[6][iPid][idxPart] += 1;
+ }
+#endif
+ } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+
+ FillHistSetCent("MC", 0, kTRUE);
+ FillHistSetCent("MCpt", 1, kTRUE);
+ FillHistSetCent("MCTPC", 2, kTRUE);
+ FillHistSetCent("MCTOF", 3, kTRUE);
+
+#if USE_PHI
+ FillHistSetCent("MCphi", 4, kTRUE);
+ FillHistSetCent("MCTPCphi", 5, kTRUE);
+ FillHistSetCent("MCTOFphi", 6, kTRUE);
+
+
+#endif
+
+ return 0;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::AddHistSetCent(const Char_t *name, const Char_t *title) {
+ TString sName(name);
+ TString sTitle(title);
+ TList *list[4];
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ fOutList->Add(new TList);
+ list[iPid] = static_cast<TList*>(fOutList->Last());
+ list[iPid]->SetName(Form("f%s_%s", name,AliEbyEPidRatioHelper::fgkPidName[iPid]));
+ list[iPid]->SetOwner(kTRUE);
+ TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
+
+ sTitle = (iPid != 0 ) ? Form("|y| < %.1f", fHelper->GetRapidityMax()) : Form(" |#eta|<%.1f", etaRange[1]);
+
+ for (Int_t idx = 1; idx <= fOrder; ++idx) {
+ list[iPid]->Add(new TProfile(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
+ Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+ 100,-0.5,99.5));
+ }
+
+ for (Int_t ii = 0; ii <= fOrder; ++ii) {
+ for (Int_t kk = 0; kk <= fOrder; ++kk) {
+ list[iPid]->Add(new TProfile(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
+ Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
+ 100,-0.5,99.5));
+ }
+ }
+
+ }
+
+ fOutList->Add(new TList);
+ TList *list_nu = static_cast<TList*>(fOutList->Last());
+ list_nu->SetName(Form("f%s_nu", name));
+ list_nu->SetOwner(kTRUE);
+
+ for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
+ list_nu->Add(new TProfile(Form("fProf%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),100,-0.5,99.5));
+ }
+
+ Int_t nBinsCent = AliEbyEPidRatioHelper::fgkfHistNBinsCent;
+ Double_t centBinRange[] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeCent[1]};
+
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ TString sNetTitle(Form("%s - %s", AliEbyEPidRatioHelper::fgkPidLatex[iPid][1], AliEbyEPidRatioHelper::fgkPidLatex[iPid][0]));
+ sTitle = (iPid != 0 ) ? Form(" |y|<%.1f", fHelper->GetRapidityMax()) : Form(" |#eta| < %.1f", etaRange[1]);
+
+ for (Int_t idx = 1; idx <= fOrder; ++idx) {
+ list[iPid]->Add(new TProfile(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid],name, idx),
+ Form("(%s)^{%d} : %s;Centrality(11);(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+ nBinsCent, centBinRange[0], centBinRange[1]));
+ }
+
+ for (Int_t ii = 0; ii <= fOrder; ++ii) {
+ for (Int_t kk = 0; kk <= fOrder; ++kk) {
+ list[iPid]->Add(new TProfile(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk),
+ Form("f_{%02d%02d} : %s;Centrality(11);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
+ nBinsCent, centBinRange[0], centBinRange[1]));
+ }
+ }
+
+ }
+
+
+ for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
+ list_nu->Add(new TProfile(Form("fProfBin%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),nBinsCent, centBinRange[0], centBinRange[1]));
+ }
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioPhy::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC) {
+ /*printf(" !! %2d !! ", idx);
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ printf("%7d %7d ", fNp[idx][iPid][0] , fNp[idx][iPid][1]);
+ }
+ Printf("");*/
+
+ Int_t ***np = (isMC) ? fMCNp : fNp;
+
+ /* printf(" ** %2d ** ", idx);
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ printf("%7d %7d ", np[idx][iPid][0] , np[idx][iPid][1]);
+ }
+ Printf("");*/
+
+ TList *list[4];
+ Float_t centralityBin = fHelper->GetCentralityBin();
+ Float_t centralityPer = fHelper->GetCentralityPercentile();//fHelper->GetCentralityBin();
+
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+
+ list[iPid] = static_cast<TList*>(fOutList->FindObject(Form("f%s_%s",name, AliEbyEPidRatioHelper::fgkPidName[iPid])));
+
+ Int_t deltaNp = np[idx][iPid][1]-np[idx][iPid][0];
+
+ Double_t delta = 1.;
+ for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
+ delta *= deltaNp;
+ (static_cast<TProfile*>(list[iPid]->FindObject(Form("fProfBin%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityBin, delta);
+ (static_cast<TProfile*>(list[iPid]->FindObject(Form("fProf%s%sNet%dM", AliEbyEPidRatioHelper::fgkPidName[iPid], name, idxOrder))))->Fill(centralityPer, delta);
+ }
+
+ for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
+ fRedFactp[idxOrder][0] = 1.;
+ fRedFactp[idxOrder][1] = 1.;
+ }
+
+ for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
+ fRedFactp[idxOrder][0] = fRedFactp[idxOrder-1][0] * Double_t(np[idx][iPid][0]-(idxOrder-1));
+ fRedFactp[idxOrder][1] = fRedFactp[idxOrder-1][1] * Double_t(np[idx][iPid][1]-(idxOrder-1));
+ }
+
+ for (Int_t ii = 0; ii <= fOrder; ++ii) { // ii -> p -> n1
+ for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
+ Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0]; // n1 *n2 -> p * pbar
+ (static_cast<TProfile*>(list[iPid]->FindObject(Form("fProfBin%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityBin, fik);
+ (static_cast<TProfile*>(list[iPid]->FindObject(Form("fProf%s%sNetF%02d%02d", AliEbyEPidRatioHelper::fgkPidName[iPid], name, ii, kk))))->Fill(centralityPer, fik);
+ }
+ }
+ }
+ Double_t a[6][4]; Double_t b[22];
+ for (Int_t iPid = 0; iPid < 4; ++iPid) {
+ a[0][iPid] = np[idx][iPid][1]+np[idx][iPid][0]; // 0 n+ + n-
+ a[1][iPid] = np[idx][iPid][1]; // 1 n+
+ a[2][iPid] = np[idx][iPid][0]; // 2 n-
+ a[3][iPid] = np[idx][iPid][1]*np[idx][iPid][0]; // 3 n+ . n-
+ a[4][iPid] = np[idx][iPid][1]*(np[idx][iPid][1]-1); // 4 n+ (n+ - 1)
+ a[5][iPid] = np[idx][iPid][0]*(np[idx][iPid][0]-1); // 5 n- (n- - 1)
+ }
+
+ b[0] = a[0][0]*a[0][2]; // 24 N K
+ b[1] = a[0][1]*a[0][2]; // 25 Pi K
+ b[2] = a[1][1]*a[1][2]; // 26 pi+ k+
+ b[3] = a[1][1]*a[2][2]; // 27 pi+ k-
+ b[4] = a[2][1]*a[1][2]; // 28 pi- k+
+ b[5] = a[2][1]*a[2][2]; // 29 pi- k-
+
+ b[6] = a[0][0]*a[0][3]; // 30 N P
+ b[7] = a[0][2]*a[0][3]; // 31 K P
+ b[8] = a[1][2]*a[1][3]; // 32 k+ p+
+ b[9] = a[1][2]*a[2][3]; // 33 k+ p-
+ b[10] = a[2][2]*a[1][3]; // 34 k- p+
+ b[11] = a[2][2]*a[2][3]; // 35 k- p-
+
+ b[12] = a[0][0]*a[0][1]; // 36 N Pi
+ b[13] = a[0][3]*a[0][1]; // 37 P Pi
+ b[14] = a[1][3]*a[1][1]; // 38 p+ pi+
+ b[15] = a[1][3]*a[2][1]; // 39 p+ pi-
+ b[16] = a[2][3]*a[1][1]; // 40 p- pi+
+ b[17] = a[2][3]*a[2][1]; // 41 p- pi-
+
+ b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
+ b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
+ b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
+ b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
+ TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
+ Int_t k = 0;
+ for (Int_t i = 0; i < 6; i++) {
+ for (Int_t j = 0; j < 4; j++) {
+ (static_cast<TProfile*>(list_nu->FindObject(Form("fProfBin%sNu%02d", name,k))))->Fill(centralityBin,a[i][j]);
+ (static_cast<TProfile*>(list_nu->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(centralityPer,a[i][j]);
+ k++;
+ }
+ }
+
+ for (Int_t j = 0; j < 22; j++) {
+ (static_cast<TProfile*>(list_nu->FindObject(Form("fProfBin%sNu%02d", name,j+23))))->Fill(centralityBin,b[j]);
+ (static_cast<TProfile*>(list_nu->FindObject(Form("fProf%sNu%02d", name,j+23))))->Fill(centralityPer,b[j]);
+ }
+
+ return;
+}
--- /dev/null
+#ifndef ALIEBYEPIDRSTIOPHY_H
+#define ALIEBYEPIDRSTIOPHY_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#include "THnSparse.h"
+#include "TList.h"
+
+#include "AliEbyEPidRatioBase.h"
+
+class AliEbyEPidRatioPhy : public AliEbyEPidRatioBase {
+
+ public:
+
+ AliEbyEPidRatioPhy();
+ virtual ~AliEbyEPidRatioPhy();
+ virtual void Process();
+ void SetOutList(TList* l) {fOutList = l;}
+
+ private:
+
+ AliEbyEPidRatioPhy(const AliEbyEPidRatioPhy&);
+ AliEbyEPidRatioPhy& operator=(const AliEbyEPidRatioPhy&);
+ virtual void Init();
+ virtual void Reset();
+ virtual void CreateHistograms();
+ Int_t ProcessTracks();
+ Int_t ProcessParticles();
+ void ResetHistSet();
+ void AddHistSetCent(const Char_t *name, const Char_t *title);
+ void FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC);
+
+ TList *fOutList; //! Output data container
+ Int_t fOrder; // Max order of higher order distributions
+ Int_t fNNp; // N sets of arrays of particle/anti-particle counts
+ Int_t ***fNp; // Array of particle/anti-particle counts
+ Int_t fNMCNp; // N sets of arrays of MC particle/anti-particle counts
+ Int_t ***fMCNp; // Array of MC particle/anti-particle counts
+ Double_t **fRedFactp; // Array of particle/anti-particle reduced factorial
+
+ ClassDef(AliEbyEPidRatioPhy, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#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 "AliAODInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliEbyEPidRatioQA.h"
+
+using namespace std;
+
+ClassImp(AliEbyEPidRatioQA)
+//________________________________________________________________________
+AliEbyEPidRatioQA::AliEbyEPidRatioQA() :
+ AliEbyEPidRatioBase("QA", "QA"),
+
+ fHnQA(NULL) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioQA",10);
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioQA::~AliEbyEPidRatioQA() {
+ // Destructor
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioQA::CreateHistograms() {
+
+ Int_t binHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, AliEbyEPidRatioHelper::fgkfHistNBinsEta, // cent | eta
+ AliEbyEPidRatioHelper::fgkfHistNBinsRap, AliEbyEPidRatioHelper::fgkfHistNBinsPhi, // y | phi
+ AliEbyEPidRatioHelper::fgkfHistNBinsPt, AliEbyEPidRatioHelper::fgkfHistNBinsPt, // pt | pInner
+ AliEbyEPidRatioHelper::fgkfHistNBinsSign, 500, // sign | TPCsignal |
+ 50, 50, 50, // nSigmaITS | nSigmaTPC | nSigmaTOF
+ 50, 50, 50, 50, // DCAr | DCAz | nSigmaDCAr | nSigmaDCAz
+ 3,4}; // MCisProbe
+
+ Double_t minHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[0], AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[0], AliEbyEPidRatioHelper::fgkfHistRangePt[0],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 30,
+ -10., -10., -10.,
+ -10., -10., -10., -10.,
+ -0.5,-0.5};
+
+ Double_t maxHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeRap[1], AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
+ AliEbyEPidRatioHelper::fgkfHistRangePt[1], AliEbyEPidRatioHelper::fgkfHistRangePt[1],
+ AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 500,
+ 10., 10., 10.,
+ 10., 10., 10., 10.,
+ 2.5,3.5};
+
+
+ fHnQA = new THnSparseF("hnQA", "cent:eta:y:phi:pt:pInner:sign:TPCsignal:nSigmaITS:nSigmaTPC:nSigmaTOF:DCAr:DCAz:nSigmaDCAr:nSigmaDCAz:MCisProbe",
+ 17, binHnQA, minHnQA, maxHnQA);
+
+ fHnQA->Sumw2();
+
+ fHnQA->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
+ fHnQA->GetAxis(1)->SetTitle("#eta"); // eta [-0.9, 0.9]
+ fHnQA->GetAxis(2)->SetTitle("#it{y}"); // rapidity [-0.5, 0.5]
+ fHnQA->GetAxis(3)->SetTitle("#varphi"); // phi [ 0. , 2Pi]
+ fHnQA->GetAxis(4)->SetTitle("#it{p}_{T} (GeV/#it{c})"); // pt [ 0.46, 2.22]
+ fHnQA->GetAxis(5)->SetTitle("#it{p}_{Inner} (GeV/#it{c})"); // pInner [ 0.1, 3.0]
+ fHnQA->GetAxis(6)->SetTitle("sign"); // -1 | 0 | +1
+
+ fHnQA->GetAxis(7)->SetTitle("TPC signal"); // TPCsignal [30, 500]
+ fHnQA->GetAxis(8)->SetTitle("n #sigma ITS"); // nSigma ITS [-10, 10]
+ fHnQA->GetAxis(9)->SetTitle("n #sigma TPC"); // nSigma TPC [-10, 10]
+ fHnQA->GetAxis(10)->SetTitle("n #sigma TOF"); // nSigma TOF [-10, 10]
+
+ fHnQA->GetAxis(11)->SetTitle("DCAr"); // DCAr [-10, 10]
+ fHnQA->GetAxis(12)->SetTitle("DCAz"); // DCAz [-10, 10]
+ fHnQA->GetAxis(13)->SetTitle("n #sigma #sqrt(Cdd)/DCAr"); // nSigma DCAr [-10, 10]
+ fHnQA->GetAxis(14)->SetTitle("n #sigma #sqrt(Czz)/DCAz"); // nSigma DCAz [-10, 10]
+ fHnQA->GetAxis(15)->SetTitle("MCisProbe"); // 0 | 1 (isProbeParticle) | 2 (isProbeParticle wrong sign)
+ fHnQA->GetAxis(16)->SetTitle("PID"); // 0 | 1 | 2 | 3
+
+ fHelper->BinLogAxis(fHnQA, 4);
+ fHelper->BinLogAxis(fHnQA, 5);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioQA::Process() {
+ Float_t etaRange[2];
+ fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+ Float_t ptRange[2];
+ fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+
+ // -- Track Loop
+ for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
+
+ // -- Check if track is accepted for basic parameters
+ if (!fHelper->IsTrackAcceptedBasicCharged(track))
+ continue;
+
+ // -- Check if accepted - ESD
+ if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
+ continue;
+
+ // -- Check if accepted - AOD
+ if (fAOD){
+ AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
+
+ if (!trackAOD) {
+ AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
+ continue;
+ }
+ if (!trackAOD->TestFilterBit(fAODtrackCutBit))
+ continue;
+
+ // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
+ if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
+ continue;
+ }
+
+ Int_t gPdgCode = 0;
+
+ Int_t iPid = 0;
+ Double_t pid[3];
+ if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion))) { iPid = 1; gPdgCode = 211;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon))) { iPid = 2; gPdgCode = 321;}
+ else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){ iPid = 3; gPdgCode = 2212;}
+ else iPid = 0;
+
+ Double_t yP;
+ if ((iPid != 0) && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
+ continue;
+
+ if (!fHelper->IsTrackAcceptedDCA(track))
+ continue;
+
+ // -- Check if probe particle
+ Int_t isProbeParticle = 0;
+ if (fIsMC) {
+ Int_t label = TMath::Abs(track->GetLabel());
+
+ AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
+ if (particle) {
+ if (TMath::Abs(particle->PdgCode()) == gPdgCode) {
+ ++isProbeParticle;
+ if (particle->PdgCode() != (track->Charge()*gPdgCode))
+ ++isProbeParticle;
+ }
+ }
+ }
+
+ // -- Get dca r/z
+ Float_t dca[] = {0.,0.}; // dca_xy, dca_z,
+ Float_t cov[] = {0.,0.,0.}; // sigma_xy, sigma_xy_z, sigma_z
+ if (fESD)
+ (static_cast<AliESDtrack*>(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;
+
+ if (iPid == 0)
+ yP = track->Eta();
+
+
+ // cout << pid[0] << " " << pid[1] << " " << pid[2] << yP << " " << iPid << endl;
+
+ if (iPid != 0) {
+ Double_t aTrack[17] = {
+ Double_t(fCentralityBin), // 0 centrality
+ track->Eta(), // 1 eta
+ yP, // 2 rapidity
+ track->Phi(), // 3 phi
+ track->Pt(), // 4 pt
+ track->GetTPCmomentum(), // 5 pInner
+ track->Charge(), // 6 sign
+ track->GetTPCsignal(), // 7 TPC dE/dx
+ pid[0], // 8 n Sigma ITS
+ pid[1], // 9 n Sigma TPC
+ pid[2], // 10 n Sigma TOF
+ dca[0], // 11 dca r
+ dca[1], // 12 dca z
+ dcaRoverCdd, // 13 sqrt(cov[dd])
+ dcaZoverCzz, // 14 sqrt(cov[zz])
+ isProbeParticle, // 15 isProbe
+ 0 // 16 PID
+ };
+
+ fHnQA->Fill(aTrack);
+ }
+
+ Double_t aTrack[17] = {
+ Double_t(fCentralityBin), // 0 centrality
+ track->Eta(), // 1 eta
+ yP, // 2 rapidity
+ track->Phi(), // 3 phi
+ track->Pt(), // 4 pt
+ track->GetTPCmomentum(), // 5 pInner
+ track->Charge(), // 6 sign
+ track->GetTPCsignal(), // 7 TPC dE/dx
+ pid[0], // 8 n Sigma ITS
+ pid[1], // 9 n Sigma TPC
+ pid[2], // 10 n Sigma TOF
+ dca[0], // 11 dca r
+ dca[1], // 12 dca z
+ dcaRoverCdd, // 13 sqrt(cov[dd])
+ dcaZoverCzz, // 14 sqrt(cov[zz])
+ isProbeParticle, // 15 isProbe
+ iPid // 16 PID
+ };
+
+ fHnQA->Fill(aTrack);
+
+ } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
+
+ return;
+}
--- /dev/null
+#ifndef ALIEBYEPIDRATIOQA_H
+#define ALIEBYEPIDRATIOQA_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+
+#include "THnSparse.h"
+#include "AliEbyEPidRatioBase.h"
+
+class AliEbyEPidRatioQA : public AliEbyEPidRatioBase {
+
+ public:
+
+ AliEbyEPidRatioQA();
+ virtual ~AliEbyEPidRatioQA();
+ virtual void Process();
+ THnSparseF* GetHnQA() {return fHnQA;}
+
+ private:
+
+ AliEbyEPidRatioQA(const AliEbyEPidRatioQA&); // not implemented
+ AliEbyEPidRatioQA& operator=(const AliEbyEPidRatioQA&); // not implemented
+ virtual void CreateHistograms();
+
+ THnSparseF *fHnQA; //! THnSparseF : tracks for QA
+
+ ClassDef(AliEbyEPidRatioQA, 1);
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: ALICE Offline. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+#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 "AliEbyEPidRatioTask.h"
+#include "AliGenEventHeader.h"
+#include "AliCentrality.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+
+using namespace std;
+ClassImp(AliEbyEPidRatioTask)
+//________________________________________________________________________
+AliEbyEPidRatioTask::AliEbyEPidRatioTask(const char *name) :
+ AliAnalysisTaskSE(name),
+ fHelper(NULL),
+ fEffCont(NULL),
+ fDCA(NULL),
+ fDist(NULL),
+ fQA(NULL),
+
+ fOutList(NULL),
+ fOutListEff(NULL),
+ fOutListCont(NULL),
+ fOutListDCA(NULL),
+ fOutListQA(NULL),
+
+ fESD(NULL),
+ fESDHandler(NULL),
+
+ fESDTrackCutsBase(NULL),
+ fESDTrackCuts(NULL),
+ fESDTrackCutsBkg(NULL),
+ fESDTrackCutsEff(NULL),
+
+ fAOD(NULL),
+ fAODHandler(NULL),
+
+ fIsMC(kFALSE),
+ fIsAOD(kFALSE),
+ fESDTrackCutMode(0),
+ fModeEffCreation(0),
+ fModeDCACreation(0),
+ fModeDistCreation(0),
+ fModeQACreation(0),
+
+ fMCEvent(NULL),
+ fMCStack(NULL),
+
+ fEtaMax(0.8),
+ fEtaMaxEff(0.9),
+ fPtRange(),
+ fPtRangeEff(),
+
+ fAODtrackCutBit(1024) {
+ // Constructor
+
+ AliLog::SetClassDebugLevel("AliEbyEPidRatioTask",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());
+}
+
+//________________________________________________________________________
+AliEbyEPidRatioTask::~AliEbyEPidRatioTask() {
+ // Destructor
+
+ if (fESDTrackCutsBase) delete fESDTrackCutsBase;
+ if (fESDTrackCuts) delete fESDTrackCuts;
+ if (fESDTrackCutsBkg) delete fESDTrackCutsBkg;
+ if (fESDTrackCutsEff) delete fESDTrackCutsEff;
+
+ if (fEffCont) delete fEffCont;
+ if (fDCA) delete fDCA;
+ if (fDist) delete fDist;
+ if (fQA) delete fQA;
+ if (fHelper) delete fHelper;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Public Methods
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+void AliEbyEPidRatioTask::UserCreateOutputObjects() {
+ // Create histograms
+
+ // ------------------------------------------------------------------
+ // -- Create Output Lists
+ // ------------------------------------------------------------------
+ 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) ;
+
+ Initialize();
+
+ fOutList->Add(new TList);
+ TList *list = static_cast<TList*>(fOutList->Last());
+ list->SetName(Form("fStat"));
+ list->SetOwner(kTRUE);
+
+ list->Add(fHelper->GetHEventStat0());
+ list->Add(fHelper->GetHEventStat1());
+ list->Add(fHelper->GetHTriggerStat());
+ list->Add(fHelper->GetHCentralityStat());
+
+ if ((fIsAOD||fIsMC) && fModeEffCreation == 1) {
+ fOutListEff->Add(fEffCont->GetHnEff());
+ fOutListCont->Add(fEffCont->GetHnCont());
+ }
+
+ if (fModeDCACreation == 1)
+ fOutListDCA->Add(fDCA->GetHnDCA());
+
+ if (fModeQACreation == 1)
+ fOutListQA->Add(fQA->GetHnQA());
+
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1,fOutList);
+ PostData(2,fOutListEff);
+ PostData(3,fOutListCont);
+ PostData(4,fOutListDCA);
+ PostData(5,fOutListQA);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioTask::UserExec(Option_t *) {
+
+ if (SetupEvent() < 0) {
+ PostData(1,fOutList);
+ PostData(2,fOutListEff);
+ PostData(3,fOutListCont);
+ PostData(4,fOutListDCA);
+ PostData(5,fOutListQA);
+ return;
+ }
+
+
+ if ((fIsMC||fIsAOD) && fModeEffCreation == 1)
+ fEffCont->Process();
+
+ if (fModeDCACreation == 1)
+ fDCA->Process();
+
+ if (fModeDistCreation == 1)
+ fDist->Process();
+
+ if (fModeQACreation == 1)
+ fQA->Process();
+
+ PostData(1,fOutList);
+ PostData(2,fOutListEff);
+ PostData(3,fOutListCont);
+ PostData(4,fOutListDCA);
+ PostData(5,fOutListQA);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliEbyEPidRatioTask::Terminate(Option_t *){
+ // Terminate
+}
+
+Int_t AliEbyEPidRatioTask::Initialize() {
+ // Initialize event
+
+ // ------------------------------------------------------------------
+ // -- ESD TrackCuts
+ // ------------------------------------------------------------------
+ TString sModeName("");
+
+ // -- Create ESD track cuts
+ // --------------------------
+ fESDTrackCutsBase = new AliESDtrackCuts;
+
+ if (fESDTrackCutMode == 0) {
+ fESDTrackCutsBase->SetMinNCrossedRowsTPC(70); // TPC
+ fESDTrackCutsBase->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); // TPC
+ }
+ else if (fESDTrackCutMode == 1) {
+ fESDTrackCutsBase->SetMinNClustersTPC(70); // TPC 2010
+ }
+
+ fESDTrackCutsBase->SetMaxChi2PerClusterTPC(4); // TPC 2010
+ fESDTrackCutsBase->SetAcceptKinkDaughters(kFALSE); // TPC 2010
+ fESDTrackCutsBase->SetRequireTPCRefit(kTRUE); // TPC 2010
+
+ if (fESDTrackCutMode == 0) {
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff); // ITS
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kOff); // ITS
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSSD,AliESDtrackCuts::kOff); // ITS
+ }
+ else if (fESDTrackCutMode == 1) {
+ fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); // ITS 2010
+ // fESDTrackCutsBase->SetMinNClustersITS(4);
+ }
+
+ fESDTrackCutsBase->SetRequireITSRefit(kTRUE); // ITS 2010
+ fESDTrackCutsBase->SetMaxChi2PerClusterITS(36); // ITS 2010
+
+ fESDTrackCutsBase->SetDCAToVertex2D(kFALSE); // VertexConstrained 2010
+ fESDTrackCutsBase->SetRequireSigmaToVertex(kFALSE); // VertexConstrained 2010
+ fESDTrackCutsBase->SetMaxDCAToVertexZ(2); // VertexConstrained 2010
+
+ fESDTrackCutsBase->SetEtaRange(-1.*fEtaMax, fEtaMax); // Acceptance
+ fESDTrackCutsBase->SetPtRange(fPtRange[0],fPtRange[1]); // Acceptance
+
+ // -- Mode : standard cuts
+ if (fESDTrackCutMode == 0)
+ sModeName = "Std";
+ // -- Mode : for comparison to LF
+ else if (fESDTrackCutMode == 1)
+ sModeName = "LF";
+ // -- Mode : Default
+ else
+ sModeName = "Base";
+
+ fESDTrackCutsBase->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
+
+ // -- Create ESD track cuts -> Base + DCA
+ // ------------------------------
+ fESDTrackCuts = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
+ fESDTrackCuts->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
+ if (fESDTrackCutMode == 0)
+ fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); // 2010 VertexConstrained -> 7*(0.0026+0.0050/pt^1.01)
+ // fESDTrackCuts->SetMaxDCAToVertexXY(0.3);
+ else if (fESDTrackCutMode == 1)
+ fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); // 2010 VertexConstrained -> 7*(0.0026+0.0050/pt^1.01)
+
+ // fESDTrackCuts->SetMaxChi2TPCConstrainedGlobal(36); // golden cut off
+
+ // -- Create ESD BKG track cuts -> Base + Acceptance(Eff)
+ // ------------------------------
+ fESDTrackCutsBkg = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
+ fESDTrackCutsBkg->SetName(Form("NetParticleCuts2010_%s_Bkg",sModeName.Data()));
+ fESDTrackCutsBkg->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]); // Acceptance
+ fESDTrackCutsBkg->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff); // Acceptance
+
+ // -- Create ESD Eff track cuts -> Base + DCA + Acceptance(Eff)
+ // ------------------------------
+ fESDTrackCutsEff = static_cast<AliESDtrackCuts*>(fESDTrackCuts->Clone());
+ fESDTrackCutsEff->SetName(Form("NetParticleCuts2010_%s_Eff",sModeName.Data()));
+ fESDTrackCutsEff->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]); // Acceptance
+ fESDTrackCutsEff->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff); // Acceptance
+
+ // ------------------------------------------------------------------
+ // -- Initialize Helper
+ // ------------------------------------------------------------------
+ if (fHelper->Initialize(fESDTrackCutsEff, fIsMC, fAODtrackCutBit, fModeDistCreation))
+ return -1;
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize Efficiency/Contamination
+ // ------------------------------------------------------------------
+ if ((fIsMC||fIsAOD) && fModeEffCreation == 1) {
+ fEffCont = new AliEbyEPidRatioEffCont;
+ fEffCont->Initialize(fHelper);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize DCA Determination
+ // ------------------------------------------------------------------
+ if (fModeDCACreation == 1) {
+ fDCA = new AliEbyEPidRatioDCA;
+ fDCA->SetESDTrackCutsBkg(fESDTrackCutsBkg);
+ fDCA->Initialize(fHelper);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize Phy Determination
+ // ------------------------------------------------------------------
+ if (fModeDistCreation == 1) {
+ fDist = new AliEbyEPidRatioPhy;
+ fDist->SetOutList(fOutList);
+ fDist->Initialize(fHelper, fESDTrackCuts);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Create / Initialize QA Determination
+ // ------------------------------------------------------------------
+ if (fModeQACreation == 1) {
+ fQA = new AliEbyEPidRatioQA();
+ fQA->Initialize(fHelper);
+ }
+
+ // ------------------------------------------------------------------
+ // -- Reset Event
+ // ------------------------------------------------------------------
+ ResetEvent();
+
+ return 0;
+}
+
+/*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioTask::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 (!fIsAOD && SetupESDEvent() < 0) {
+ AliError("Setup ESD Event failed");
+ return -1;
+ }
+
+ // -- AOD Event
+ // ------------------------------------------------------------------
+ if (fIsAOD && SetupAODEvent() < 0) {
+ AliError("Setup AOD 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 / QA classes
+ // ------------------------------------------------------------------
+ fHelper->SetupEvent(fESDHandler, fAODHandler, fMCEvent);
+
+
+ if (fModeEffCreation && (fIsMC || fIsAOD) )
+ fEffCont->SetupEvent();
+
+ if (fModeDCACreation == 1)
+ fDCA->SetupEvent();
+
+ if (fModeDistCreation == 1)
+ fDist->SetupEvent();
+
+ if (fModeQACreation == 1)
+ fQA->SetupEvent();
+
+
+ // -- Evaluate Event cuts
+ // ------------------------------------------------------------------
+ return fHelper->IsEventRejected() ? -2 : 0;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioTask::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 AliEbyEPidRatioTask::SetupAODEvent() {
+ // -- Setup AOD Event
+ // > return 0 for success
+ // > return -1 for failed setup
+
+ fAODHandler= dynamic_cast<AliAODInputHandler*>
+ (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ if (!fAODHandler) {
+ AliError("Could not get AOD input handler");
+ return -1;
+ }
+
+ fAOD = fAODHandler->GetEvent();
+ if (!fAOD) {
+ AliError("Could not get AOD event");
+ return -1;
+ }
+
+ // -- Check PID response
+ // ------------------------------------------------------------------
+ if (!fAODHandler->GetPIDResponse()) {
+ AliError("Could not get PID response");
+ return -1;
+ }
+
+ // -- Check Vertex
+ // ------------------------------------------------------------------
+ if (!fAOD->GetPrimaryVertex()) {
+ AliError("Could not get primary vertex");
+ return -1;
+ }
+
+ // -- Check Centrality
+ // ------------------------------------------------------------------
+ if (!fAOD->GetHeader()->GetCentralityP()) {
+ AliError("Could not get centrality");
+ return -1;
+ }
+
+ return 0;
+}
+
+//________________________________________________________________________
+Int_t AliEbyEPidRatioTask::SetupMCEvent() {
+ AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>
+ (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+
+ if (!mcH) {
+ AliError("MC event handler not available");
+ return -1;
+ }
+
+ 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 AliEbyEPidRatioTask::ResetEvent() {
+ // -- Reset event
+
+ // -- Reset ESD Event
+ fESD = NULL;
+
+ // -- Reset AOD Event
+ fAOD = NULL;
+
+ // -- Reset MC Event
+ if (fIsMC)
+ fMCEvent = NULL;
+
+ // -- Reset Dist Creation
+ if (fModeDistCreation == 1)
+ fDist->ResetEvent();
+
+ return;
+}
+
+
--- /dev/null
+#ifndef ALIEBYEPIDRATIOTASK_H
+#define ALIEBYEPIDRATIOTASK_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+//=========================================================================//
+// AliEbyE Analysis for Particle Ratio Fluctuation //
+// Deepika Rathee | Satyajit Jena //
+// drathee@cern.ch | sjena@cern.ch //
+// Date: Wed Jul 9 18:38:30 CEST 2014 //
+// New approch to find particle ratio to reduce memory //
+// (Test Only) //
+//=========================================================================//
+
+
+#include "AliAnalysisTaskSE.h"
+#include "TList.h"
+#include "THnSparse.h"
+
+#include "AliEbyEPidRatioHelper.h"
+#include "AliEbyEPidRatioEffCont.h"
+#include "AliEbyEPidRatioDCA.h"
+#include "AliEbyEPidRatioPhy.h"
+#include "AliEbyEPidRatioQA.h"
+
+class AliESDEvent;
+class AliESDInputHandler;
+class AliAODEvent;
+class AliAODInputHandler;
+class AliMCEvent;
+class AliStack;
+class AliPIDResponse;
+class AliESDResponse;
+
+class AliEbyEPidRatioTask : public AliAnalysisTaskSE {
+
+ public:
+
+ AliEbyEPidRatioTask(const char *name = "AliEbyEPidRatioTask");
+ virtual ~AliEbyEPidRatioTask();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * 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 SetIsAOD(Bool_t b) {fIsAOD = 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 SetModeQACreation(Int_t i) {fModeQACreation = i;}
+
+ void SetEtaMax(Float_t f) {fEtaMax = f;}
+ void SetEtaMaxEff(Float_t f) {fEtaMaxEff = 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;}
+
+ void SetTrackFilterBit(Int_t i) {fAODtrackCutBit = i;}
+
+ void SetNetParticleHelper(AliEbyEPidRatioHelper *helper) {
+ if (fHelper)
+ delete fHelper;
+ fHelper = helper;
+ }
+
+ ///////////////////////////////////////////////////////////////////////////////////
+
+ private:
+
+ AliEbyEPidRatioTask(const AliEbyEPidRatioTask&); // not implemented
+ AliEbyEPidRatioTask& operator=(const AliEbyEPidRatioTask&); // not implemented
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Setup/Reset Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ /** Setup Event */
+ Int_t SetupEvent();
+
+ /** Setup ESD Event */
+ Int_t SetupESDEvent();
+
+ /** Setup AOD Event */
+ Int_t SetupAODEvent();
+
+ /** Setup MC Event */
+ Int_t SetupMCEvent();
+
+ /** Reset Event */
+ void ResetEvent();
+
+ /*
+ * ---------------------------------------------------------------------------------
+ * Members - private
+ * ---------------------------------------------------------------------------------
+ */
+
+ AliEbyEPidRatioHelper *fHelper; // Helper class
+ AliEbyEPidRatioEffCont *fEffCont; //! Efficiency and Contamination class
+ AliEbyEPidRatioDCA *fDCA; //! DCA class
+ AliEbyEPidRatioPhy *fDist; //! Distributions class
+ AliEbyEPidRatioQA *fQA; //! QA 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 output data container - QA
+
+ // --- 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
+
+ // --- AOD only ----------------------------------------------------------
+ AliAODEvent *fAOD; //! Ptr to AOD event
+ AliAODInputHandler *fAODHandler; //! Ptr to AOD Handler
+
+ // --- Flags -------------------------------------------------------------
+ Bool_t fIsMC; // Is MC event
+ Bool_t fIsAOD; // analysis mode : 0 = ESDs | 1 = AODs
+ 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
+ Int_t fModeQACreation; // QA creation mode : 1 = on | 0 = off
+
+ // --- MC only -----------------------------------------------------------
+ AliMCEvent *fMCEvent; //! Ptr to MC event
+ AliStack *fMCStack; //! Ptr to MC stack
+ // -----------------------------------------------------------------------
+ Float_t fEtaMax; // Max, absolut eta
+ Float_t fEtaMaxEff; // Max, absolut eta for efficiency
+ Float_t fPtRange[2]; // Array of pt [min,max]
+ Float_t fPtRangeEff[2]; // Array of pt [min,max] for efficiency
+ // -----------------------------------------------------------------------
+ Int_t fAODtrackCutBit; // Track filter bit for AOD tracks
+ // =======================================================================
+
+
+ ClassDef(AliEbyEPidRatioTask, 1);
+};
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff