-#include "TChain.h"\r
-#include "TList.h"\r
-#include "TCanvas.h"\r
-#include "TLorentzVector.h"\r
-#include "TGraphErrors.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "TH3F.h" \r
-#include "TH2D.h" \r
-#include "TH3D.h"\r
-#include "TArrayF.h"\r
-#include "TF1.h"\r
-#include "TRandom.h"\r
-\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliAnalysisManager.h"\r
-\r
-#include "AliESDVertex.h"\r
-#include "AliESDEvent.h"\r
-#include "AliESDInputHandler.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODInputHandler.h"\r
-#include "AliAODMCParticle.h" \r
-#include "AliCollisionGeometry.h"\r
-#include "AliGenEventHeader.h"\r
-#include "AliMCEventHandler.h"\r
-#include "AliMCEvent.h"\r
-#include "AliStack.h"\r
-#include "AliESDtrackCuts.h"\r
-#include "AliEventplane.h"\r
-#include "AliTHn.h" \r
-\r
-#include "AliEventPoolManager.h" \r
-\r
-#include "AliPID.h" \r
-#include "AliPIDResponse.h" \r
-#include "AliPIDCombined.h" \r
-\r
-#include "AliAnalysisTaskBFPsi.h"\r
-#include "AliBalancePsi.h"\r
-#include "AliAnalysisTaskTriggeredBF.h"\r
-\r
-\r
-// Analysis task for the BF vs Psi code\r
-// Authors: Panos.Christakoglou@nikhef.nl\r
-\r
-using std::cout;\r
-using std::endl;\r
-\r
-ClassImp(AliAnalysisTaskBFPsi)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
-: AliAnalysisTaskSE(name), \r
- fBalance(0),\r
- fRunShuffling(kFALSE),\r
- fShuffledBalance(0),\r
- fRunMixing(kFALSE),\r
- fRunMixingEventPlane(kFALSE),\r
- fMixingTracks(50000),\r
- fMixedBalance(0),\r
- fPoolMgr(0),\r
- fList(0),\r
- fListBF(0),\r
- fListBFS(0),\r
- fListBFM(0),\r
- fHistListPIDQA(0),\r
- fHistEventStats(0),\r
- fHistCentStats(0),\r
- fHistTriggerStats(0),\r
- fHistTrackStats(0),\r
- fHistVx(0),\r
- fHistVy(0),\r
- fHistVz(0),\r
- fHistEventPlane(0),\r
- fHistClus(0),\r
- fHistDCA(0),\r
- fHistChi2(0),\r
- fHistPt(0),\r
- fHistEta(0),\r
- fHistRapidity(0),\r
- fHistPhi(0),\r
- fHistEtaPhiPos(0), \r
- fHistEtaPhiNeg(0), \r
- fHistPhiBefore(0),\r
- fHistPhiAfter(0),\r
- fHistPhiPos(0),\r
- fHistPhiNeg(0),\r
- fHistV0M(0),\r
- fHistRefTracks(0),\r
- fHistdEdxVsPTPCbeforePID(NULL),\r
- fHistBetavsPTOFbeforePID(NULL), \r
- fHistProbTPCvsPtbeforePID(NULL), \r
- fHistProbTOFvsPtbeforePID(NULL), \r
- fHistProbTPCTOFvsPtbeforePID(NULL),\r
- fHistNSigmaTPCvsPtbeforePID(NULL), \r
- fHistNSigmaTOFvsPtbeforePID(NULL), \r
- fHistdEdxVsPTPCafterPID(NULL),\r
- fHistBetavsPTOFafterPID(NULL), \r
- fHistProbTPCvsPtafterPID(NULL), \r
- fHistProbTOFvsPtafterPID(NULL), \r
- fHistProbTPCTOFvsPtafterPID(NULL),\r
- fHistNSigmaTPCvsPtafterPID(NULL), \r
- fHistNSigmaTOFvsPtafterPID(NULL), \r
- fCentralityArrayBinsForCorrections(kCENTRALITY),\r
- fPIDResponse(0x0),\r
- fPIDCombined(0x0),\r
- fParticleOfInterest(kPion),\r
- fPidDetectorConfig(kTPCTOF),\r
- fUsePID(kFALSE),\r
- fUsePIDnSigma(kTRUE),\r
- fUsePIDPropabilities(kFALSE), \r
- fPIDNSigma(3.),\r
- fMinAcceptedPIDProbability(0.8),\r
- fElectronRejection(kFALSE),\r
- fElectronRejectionNSigma(-1.),\r
- fESDtrackCuts(0),\r
- fCentralityEstimator("V0M"),\r
- fUseCentrality(kFALSE),\r
- fCentralityPercentileMin(0.), \r
- fCentralityPercentileMax(5.),\r
- fImpactParameterMin(0.),\r
- fImpactParameterMax(20.),\r
- fUseMultiplicity(kFALSE),\r
- fNumberOfAcceptedTracksMin(0),\r
- fNumberOfAcceptedTracksMax(10000),\r
- fHistNumberOfAcceptedTracks(0),\r
- fUseOfflineTrigger(kFALSE),\r
- fVxMax(0.3),\r
- fVyMax(0.3),\r
- fVzMax(10.),\r
- nAODtrackCutBit(128),\r
- fPtMin(0.3),\r
- fPtMax(1.5),\r
- fPtMinForCorrections(0.3),//=================================correction\r
- fPtMaxForCorrections(1.5),//=================================correction\r
- fPtBinForCorrections(36), //=================================correction\r
- fEtaMin(-0.8),\r
- fEtaMax(-0.8),\r
- fEtaMinForCorrections(-0.8),//=================================correction\r
- fEtaMaxForCorrections(0.8),//=================================correction\r
- fEtaBinForCorrections(16), //=================================correction\r
- fPhiMin(0.),\r
- fPhiMax(360.),\r
- fPhiMinForCorrections(0.),//=================================correction\r
- fPhiMaxForCorrections(360.),//=================================correction\r
- fPhiBinForCorrections(100), //=================================correction\r
- fDCAxyCut(-1),\r
- fDCAzCut(-1),\r
- fTPCchi2Cut(-1),\r
- fNClustersTPCCut(-1),\r
- fAcceptanceParameterization(0),\r
- fDifferentialV2(0),\r
- fUseFlowAfterBurner(kFALSE),\r
- fExcludeResonancesInMC(kFALSE),\r
- fUseMCPdgCode(kFALSE),\r
- fPDGCodeToBeAnalyzed(-1),\r
- fEventClass("EventPlane") \r
-{\r
- // Constructor\r
- // Define input and output slots here\r
- // Input slot #0 works with a TChain\r
-\r
- //======================================================correction\r
- for (Int_t i=0; i<kCENTRALITY; i++){\r
- fHistCorrectionPlus[i] = NULL; \r
- fHistCorrectionMinus[i] = NULL; \r
- fCentralityArrayForCorrections[i] = -1.;\r
- }\r
- //=====================================================correction\r
-\r
- DefineInput(0, TChain::Class());\r
- // Output slot #0 writes into a TH1 container\r
- DefineOutput(1, TList::Class());\r
- DefineOutput(2, TList::Class());\r
- DefineOutput(3, TList::Class());\r
- DefineOutput(4, TList::Class());\r
- DefineOutput(5, TList::Class());\r
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
-\r
- // delete fBalance; \r
- // delete fShuffledBalance; \r
- // delete fList;\r
- // delete fListBF; \r
- // delete fListBFS;\r
-\r
- // delete fHistEventStats; \r
- // delete fHistTrackStats; \r
- // delete fHistVx; \r
- // delete fHistVy; \r
- // delete fHistVz; \r
-\r
- // delete fHistClus;\r
- // delete fHistDCA;\r
- // delete fHistChi2;\r
- // delete fHistPt;\r
- // delete fHistEta;\r
- // delete fHistPhi;\r
- // delete fHistEtaPhiPos; \r
- // delete fHistEtaPhiNeg;\r
- // delete fHistV0M;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
- // Create histograms\r
- // Called once\r
-\r
- // global switch disabling the reference \r
- // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
- Bool_t oldStatus = TH1::AddDirectoryStatus();\r
- TH1::AddDirectory(kFALSE);\r
-\r
- if(!fBalance) {\r
- fBalance = new AliBalancePsi();\r
- fBalance->SetAnalysisLevel("ESD");\r
- fBalance->SetEventClass(fEventClass);\r
- //fBalance->SetNumberOfBins(-1,16);\r
- //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
- }\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance) {\r
- fShuffledBalance = new AliBalancePsi();\r
- fShuffledBalance->SetAnalysisLevel("ESD");\r
- fShuffledBalance->SetEventClass(fEventClass);\r
- //fShuffledBalance->SetNumberOfBins(-1,16);\r
- //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
- }\r
- }\r
- if(fRunMixing) {\r
- if(!fMixedBalance) {\r
- fMixedBalance = new AliBalancePsi();\r
- fMixedBalance->SetAnalysisLevel("ESD");\r
- fMixedBalance->SetEventClass(fEventClass);\r
- }\r
- }\r
-\r
- //QA list\r
- fList = new TList();\r
- fList->SetName("listQA");\r
- fList->SetOwner();\r
-\r
- //Balance Function list\r
- fListBF = new TList();\r
- fListBF->SetName("listBF");\r
- fListBF->SetOwner();\r
-\r
- if(fRunShuffling) {\r
- fListBFS = new TList();\r
- fListBFS->SetName("listBFShuffled");\r
- fListBFS->SetOwner();\r
- }\r
-\r
- if(fRunMixing) {\r
- fListBFM = new TList();\r
- fListBFM->SetName("listTriggeredBFMixed");\r
- fListBFM->SetOwner();\r
- }\r
-\r
- //PID QA list\r
- if(fUsePID || fElectronRejection) {\r
- fHistListPIDQA = new TList();\r
- fHistListPIDQA->SetName("listQAPID");\r
- fHistListPIDQA->SetOwner();\r
- }\r
-\r
- //Event stats.\r
- TString gCutName[5] = {"Total","Offline trigger",\r
- "Vertex","Analyzed","sel. Centrality"};\r
- fHistEventStats = new TH2F("fHistEventStats",\r
- "Event statistics;;Centrality percentile;N_{events}",\r
- 5,0.5,5.5,220,-5,105);\r
- for(Int_t i = 1; i <= 5; i++)\r
- fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
- fList->Add(fHistEventStats);\r
-\r
- TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
- fHistCentStats = new TH2F("fHistCentStats",\r
- "Centrality statistics;;Cent percentile",\r
- 9,-0.5,8.5,220,-5,105);\r
- for(Int_t i = 1; i <= 9; i++)\r
- fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
- fList->Add(fHistCentStats);\r
-\r
- fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
- fList->Add(fHistTriggerStats);\r
-\r
- fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
- fList->Add(fHistTrackStats);\r
-\r
- fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
- fList->Add(fHistNumberOfAcceptedTracks);\r
-\r
- // Vertex distributions\r
- fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
- fList->Add(fHistVx);\r
- fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
- fList->Add(fHistVy);\r
- fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r
- fList->Add(fHistVz);\r
-\r
- //Event plane\r
- fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r
- fList->Add(fHistEventPlane);\r
-\r
- // QA histograms\r
- fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
- fList->Add(fHistClus);\r
- fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);\r
- fList->Add(fHistChi2);\r
- fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
- fList->Add(fHistDCA);\r
- fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);\r
- fList->Add(fHistPt);\r
- fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);\r
- fList->Add(fHistEta);\r
- fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r
- fList->Add(fHistRapidity);\r
- fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
- fList->Add(fHistPhi);\r
- fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
- fList->Add(fHistEtaPhiPos); \r
- fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); \r
- fList->Add(fHistEtaPhiNeg);\r
- fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
- fList->Add(fHistPhiBefore);\r
- fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
- fList->Add(fHistPhiAfter);\r
- fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
- fList->Add(fHistPhiPos);\r
- fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r
- fList->Add(fHistPhiNeg);\r
- fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
- fList->Add(fHistV0M);\r
- TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
- fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
- for(Int_t i = 1; i <= 6; i++)\r
- fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
- fList->Add(fHistRefTracks);\r
-\r
- // QA histograms for different cuts\r
- fList->Add(fBalance->GetQAHistHBTbefore());\r
- fList->Add(fBalance->GetQAHistHBTafter());\r
- fList->Add(fBalance->GetQAHistConversionbefore());\r
- fList->Add(fBalance->GetQAHistConversionafter());\r
- fList->Add(fBalance->GetQAHistPsiMinusPhi());\r
- fList->Add(fBalance->GetQAHistResonancesBefore());\r
- fList->Add(fBalance->GetQAHistResonancesRho());\r
- fList->Add(fBalance->GetQAHistResonancesK0());\r
- fList->Add(fBalance->GetQAHistResonancesLambda());\r
- fList->Add(fBalance->GetQAHistQbefore());\r
- fList->Add(fBalance->GetQAHistQafter());\r
-\r
- // Balance function histograms\r
- // Initialize histograms if not done yet\r
- if(!fBalance->GetHistNp()){\r
- AliWarning("Histograms not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fBalance->InitHistograms();\r
- }\r
-\r
- if(fRunShuffling) {\r
- if(!fShuffledBalance->GetHistNp()) {\r
- AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
- AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
- fShuffledBalance->InitHistograms();\r
- }\r
- }\r
-\r
- //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
- fListBF->Add(fBalance->GetHistNp());\r
- fListBF->Add(fBalance->GetHistNn());\r
- fListBF->Add(fBalance->GetHistNpn());\r
- fListBF->Add(fBalance->GetHistNnn());\r
- fListBF->Add(fBalance->GetHistNpp());\r
- fListBF->Add(fBalance->GetHistNnp());\r
-\r
- if(fRunShuffling) {\r
- fListBFS->Add(fShuffledBalance->GetHistNp());\r
- fListBFS->Add(fShuffledBalance->GetHistNn());\r
- fListBFS->Add(fShuffledBalance->GetHistNpn());\r
- fListBFS->Add(fShuffledBalance->GetHistNnn());\r
- fListBFS->Add(fShuffledBalance->GetHistNpp());\r
- fListBFS->Add(fShuffledBalance->GetHistNnp());\r
- } \r
-\r
- if(fRunMixing) {\r
- fListBFM->Add(fMixedBalance->GetHistNp());\r
- fListBFM->Add(fMixedBalance->GetHistNn());\r
- fListBFM->Add(fMixedBalance->GetHistNpn());\r
- fListBFM->Add(fMixedBalance->GetHistNnn());\r
- fListBFM->Add(fMixedBalance->GetHistNpp());\r
- fListBFM->Add(fMixedBalance->GetHistNnp());\r
- }\r
- //}\r
-\r
-\r
- // Event Mixing\r
- if(fRunMixing){\r
- Int_t trackDepth = fMixingTracks; \r
- Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
- \r
- // centrality bins\r
- Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* centbins = centralityBins;\r
- Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
-\r
- // multiplicity bins\r
- Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* multbins = multiplicityBins;\r
- Int_t nMultiplicityBins = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
- \r
- // Zvtx bins\r
- Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* vtxbins = vertexBins;\r
- Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
- \r
- // Event plane angle (Psi) bins\r
- Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
- Double_t* psibins = psiBins;\r
- Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;\r
- \r
- // run the event mixing also in bins of event plane (statistics!)\r
- if(fRunMixingEventPlane){\r
- if(fEventClass=="Multiplicity"){\r
- fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
- }\r
- else{\r
- fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
- }\r
- }\r
- else{\r
- if(fEventClass=="Multiplicity"){\r
- fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
- }\r
- else{\r
- fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
- }\r
- }\r
- }\r
-\r
- if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
-\r
- //====================PID========================//\r
- if(fUsePID) {\r
- fPIDCombined = new AliPIDCombined();\r
- fPIDCombined->SetDefaultTPCPriors();\r
-\r
- fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
- \r
- fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
- \r
- fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
- \r
- fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
-\r
- fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
- \r
- fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
- \r
- fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
- \r
- fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
- \r
- fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
- fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
- \r
- fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
- \r
- fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
- fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
- \r
- fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
- fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
-\r
- fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
- \r
- fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
- }\r
-\r
- // for electron rejection only TPC nsigma histograms\r
- if(!fUsePID && fElectronRejection) {\r
- \r
- fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
- \r
- fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
- \r
- fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
- fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
-\r
- fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
- fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
- }\r
- //====================PID========================//\r
-\r
- // Post output data.\r
- PostData(1, fList);\r
- PostData(2, fListBF);\r
- if(fRunShuffling) PostData(3, fListBFS);\r
- if(fRunMixing) PostData(4, fListBFM);\r
- if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID\r
-\r
- TH1::AddDirectory(oldStatus);\r
-}\r
-\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, \r
- Int_t nCentralityBins, \r
- Double_t *centralityArrayForCorrections) {\r
- //Open files that will be used for correction\r
- fCentralityArrayBinsForCorrections = nCentralityBins;\r
- for (Int_t i=0; i<nCentralityBins; i++)\r
- fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];\r
-\r
- // No file specified -> run without corrections\r
- if(!filename.Contains(".root")) {\r
- AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));\r
- return;\r
- }\r
-\r
- //Open the input file\r
- TFile *f = TFile::Open(filename);\r
- if(!f->IsOpen()) {\r
- AliInfo(Form("File %s not found --> run without corrections",filename.Data()));\r
- return;\r
- }\r
- \r
- //TString listEffName = "";\r
- for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) { \r
- //Printf("iCent %d:",iCent); \r
- TString histoName = "fHistCorrectionPlus";\r
- histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
- fHistCorrectionPlus[iCent]= dynamic_cast<TH3D *>(f->Get(histoName.Data()));\r
- if(!fHistCorrectionPlus[iCent]) {\r
- AliError(Form("fHist %s not found!!!",histoName.Data()));\r
- return;\r
- }\r
- \r
- histoName = "fHistCorrectionMinus";\r
- histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
- fHistCorrectionMinus[iCent] = dynamic_cast<TH3D *>(f->Get(histoName.Data())); \r
- if(!fHistCorrectionMinus[iCent]) {\r
- AliError(Form("fHist %s not found!!!",histoName.Data()));\r
- return; \r
- }\r
- }//loop over centralities: ONLY the PbPb case is covered\r
-\r
- if(fHistCorrectionPlus[0]){\r
- fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();\r
- fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();\r
- fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();\r
- \r
- fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();\r
- fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();\r
- fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();\r
- \r
- fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();\r
- fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();\r
- fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();\r
- }\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
- // Main loop\r
- // Called for each event\r
-\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
- Int_t gNumberOfAcceptedTracks = 0;\r
- Double_t gCentrality = -1.;\r
- Double_t gReactionPlane = -1.; \r
- Float_t bSign = 0.;\r
- \r
- // get the event (for generator level: MCEvent())\r
- AliVEvent* eventMain = NULL;\r
- if(gAnalysisLevel == "MC") {\r
- eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
- }\r
- else{\r
- eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
- \r
- // for HBT like cuts need magnetic field sign\r
- bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
- }\r
- if(!eventMain) {\r
- AliError("eventMain not available");\r
- return;\r
- }\r
- \r
- // PID Response task active?\r
- if(fUsePID || fElectronRejection) {\r
- fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
- if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
- }\r
- \r
- // check event cuts and fill event histograms\r
- if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
- return;\r
- }\r
- \r
- //Compute Multiplicity or Centrality variable\r
- Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
-\r
- // get the reaction plane\r
- if(fEventClass != "Multiplicity") {\r
- gReactionPlane = GetEventPlane(eventMain);\r
- fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
- if(gReactionPlane < 0){\r
- return;\r
- }\r
- }\r
- \r
- // get the accepted tracks in main event\r
- TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
- gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
-\r
- //multiplicity cut (used in pp)\r
- fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
-\r
- // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
- TObjArray* tracksShuffled = NULL;\r
- if(fRunShuffling){\r
- tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
- }\r
- \r
- // Event mixing \r
- if (fRunMixing)\r
- {\r
- // 1. First get an event pool corresponding in mult (cent) and\r
- // zvertex to the current event. Once initialized, the pool\r
- // should contain nMix (reduced) events. This routine does not\r
- // pre-scan the chain. The first several events of every chain\r
- // will be skipped until the needed pools are filled to the\r
- // specified depth. If the pool categories are not too rare, this\r
- // should not be a problem. If they are rare, you could lose`\r
- // statistics.\r
- \r
- // 2. Collect the whole pool's content of tracks into one TObjArray\r
- // (bgTracks), which is effectively a single background super-event.\r
- \r
- // 3. The reduced and bgTracks arrays must both be passed into\r
- // FillCorrelations(). Also nMix should be passed in, so a weight\r
- // of 1./nMix can be applied.\r
- \r
- AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
- \r
- if (!pool){\r
- AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
- }\r
- else{\r
- \r
- //pool->SetDebug(1);\r
-\r
- if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
- \r
- \r
- Int_t nMix = pool->GetCurrentNEvents();\r
- //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
- \r
- //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
- //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
- //if (pool->IsReady())\r
- //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
- \r
- // Fill mixed-event histos here \r
- for (Int_t jMix=0; jMix<nMix; jMix++) \r
- {\r
- TObjArray* tracksMixed = pool->GetEvent(jMix);\r
- fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
- }\r
- }\r
- \r
- // Update the Event pool\r
- pool->UpdatePool(tracksMain);\r
- //pool->PrintInfo();\r
- \r
- }//pool NULL check \r
- }//run mixing\r
- \r
- // calculate balance function\r
- fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
- \r
- // calculate shuffled balance function\r
- if(fRunShuffling && tracksShuffled != NULL) {\r
- fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
- }\r
-} \r
-\r
-//________________________________________________________________________\r
-Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
- // Checks the Event cuts\r
- // Fills Event statistics histograms\r
- \r
- Bool_t isSelectedMain = kTRUE;\r
- Float_t gCentrality = -1.;\r
- Float_t gRefMultiplicity = -1.;\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- fHistEventStats->Fill(1,gCentrality); //all events\r
-\r
- // Event trigger bits\r
- fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
- if(fUseOfflineTrigger)\r
- isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
- \r
- if(isSelectedMain) {\r
- fHistEventStats->Fill(2,gCentrality); //triggered events\r
-\r
- //Centrality stuff \r
- if(fUseCentrality) {\r
- if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header\r
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
- if(header){\r
- gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));\r
- fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
- \r
- // centrality QA (V0M)\r
- fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
- \r
- // centrality QA (reference tracks)\r
- fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
- fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
- fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
- fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
- fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
- fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
- fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
- fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
- fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
- }//AOD header\r
- }//AOD\r
-\r
- else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
- AliCentrality *centrality = event->GetCentrality();\r
- gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // QA for centrality estimators\r
- fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
- fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("FMD"));\r
- fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("TRK"));\r
- fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("TKL"));\r
- fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("CL0"));\r
- fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("CL1"));\r
- fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
- fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
- fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
-\r
- // centrality QA (V0M)\r
- fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
- }//ESD\r
- else if(gAnalysisLevel == "MC"){\r
- Double_t gImpactParameter = 0.;\r
- if(dynamic_cast<AliMCEvent*>(event)){\r
- AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
- if(headerH){\r
- gImpactParameter = headerH->ImpactParameter();\r
- gCentrality = gImpactParameter;\r
- }//MC header\r
- }//MC event cast\r
- }//MC\r
- else{\r
- gCentrality = -1.;\r
- }\r
- }//centrality\r
-\r
- //Multiplicity stuff \r
- if(fUseMultiplicity)\r
- gRefMultiplicity = GetRefMultiOrCentrality(event);\r
-\r
- // Event Vertex MC\r
- if(gAnalysisLevel == "MC") {\r
- if(!event) {\r
- AliError("mcEvent not available");\r
- return 0x0;\r
- }\r
- \r
- if(dynamic_cast<AliMCEvent*>(event)){\r
- AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
- if(header){ \r
- TArrayF gVertexArray;\r
- header->PrimaryVertex(gVertexArray);\r
- //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
- //gVertexArray.At(0),\r
- //gVertexArray.At(1),\r
- //gVertexArray.At(2));\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
- else \r
- fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
- if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
- if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
- if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
- else \r
- fHistEventStats->Fill(4,gCentrality); //analyzed events\r
- fHistVx->Fill(gVertexArray.At(0));\r
- fHistVy->Fill(gVertexArray.At(1));\r
- fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
- \r
- // take only events inside centrality class\r
- if(fUseCentrality) {\r
- if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
- fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
- return gCentrality; \r
- }//centrality class\r
- }\r
- // take events only within the same multiplicity class\r
- else if(fUseMultiplicity) {\r
- if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
- fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
- return gRefMultiplicity;\r
- }\r
- }//multiplicity range\r
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//header \r
- }//MC event object\r
- }//MC\r
- \r
- // Event Vertex AOD, ESD, ESDMC\r
- else{\r
- const AliVVertex *vertex = event->GetPrimaryVertex();\r
- \r
- if(vertex) {\r
- Double32_t fCov[6];\r
- vertex->GetCovarianceMatrix(fCov);\r
- if(vertex->GetNContributors() > 0) {\r
- if(fCov[5] != 0) {\r
- if(fUseMultiplicity) \r
- fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
- else if(fUseCentrality)\r
- fHistEventStats->Fill(3,gCentrality); //proper vertex\r
- if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
- if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
- if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
- if(fUseMultiplicity) {\r
- //cout<<"Filling 4 for multiplicity..."<<endl;\r
- fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
- }\r
- else if(fUseCentrality) {\r
- //cout<<"Filling 4 for centrality..."<<endl;\r
- fHistEventStats->Fill(4,gCentrality); //analyzed events\r
- }\r
- fHistVx->Fill(vertex->GetX());\r
- fHistVy->Fill(vertex->GetY());\r
- fHistVz->Fill(vertex->GetZ(),gCentrality);\r
- \r
- // take only events inside centrality class\r
- if(fUseCentrality) {\r
- //cout<<"Centrality..."<<endl;\r
- if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
- fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
- return gCentrality; \r
- }//centrality class\r
- }\r
- // take events only within the same multiplicity class\r
- else if(fUseMultiplicity) {\r
- //cout<<"N(min): "<<fNumberOfAcceptedTracksMin<<" - N(max): "<<fNumberOfAcceptedTracksMax<<" - Nref: "<<gRefMultiplicity<<endl;\r
-\r
- if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
- fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
- return gRefMultiplicity;\r
- }\r
- }//multiplicity range\r
- }//Vz cut\r
- }//Vy cut\r
- }//Vx cut\r
- }//proper vertex resolution\r
- }//proper number of contributors\r
- }//vertex object valid\r
- }//triggered event \r
- }//AOD,ESD,ESDMC\r
- \r
- // in all other cases return -1 (event not accepted)\r
- return -1;\r
-}\r
-\r
-\r
-//________________________________________________________________________\r
-Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
- // Checks the Event cuts\r
- // Fills Event statistics histograms\r
- \r
- Float_t gCentrality = -1.;\r
- Double_t gMultiplicity = 0.;\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- if(fEventClass == "Centrality"){\r
- if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header\r
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
- if(header){\r
- gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
- }//AOD header\r
- }//AOD\r
- \r
- else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
- AliCentrality *centrality = event->GetCentrality();\r
- gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
- }//ESD\r
- else if(gAnalysisLevel == "MC"){\r
- Double_t gImpactParameter = 0.;\r
- if(dynamic_cast<AliMCEvent*>(event)){\r
- AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader()); \r
- if(headerH){\r
- gImpactParameter = headerH->ImpactParameter();\r
- gCentrality = gImpactParameter;\r
- }//MC header\r
- }//MC event cast\r
- }//MC\r
- else{\r
- gCentrality = -1.;\r
- }\r
- }//End if "Centrality"\r
- if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
- AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
- if(gESDEvent){\r
- gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
- }//AliESDevent cast\r
- }\r
- else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "AOD"){\r
- AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
- if(header){\r
- gMultiplicity = header->GetRefMultiplicity();\r
- }//AOD header\r
- }\r
- else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "MC") {\r
- AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
- //Calculating the multiplicity as the number of charged primaries\r
- //within \pm 0.8 in eta and pT > 0.1 GeV/c\r
- for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r
- AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r
- if (!track) {\r
- AliError(Form("Could not receive particle %d", iParticle));\r
- continue;\r
- }\r
- \r
- //exclude non stable particles\r
- if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
- if(track->Pt() < 0.1) continue;\r
- if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
- if(track->Charge() == 0) continue;\r
-\r
- gMultiplicity += 1;\r
- }//loop over primaries\r
- }//MC mode & multiplicity class\r
-\r
- Double_t lReturnVal = -100;\r
- if(fEventClass=="Multiplicity"){\r
- lReturnVal = gMultiplicity;\r
- }else if(fEventClass=="Centrality"){\r
- lReturnVal = gCentrality;\r
- }\r
- return lReturnVal;\r
-}\r
-\r
-//________________________________________________________________________\r
-Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
- // Get the event plane\r
-\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- Float_t gVZEROEventPlane = -10.;\r
- Float_t gReactionPlane = -10.;\r
- Double_t qxTot = 0.0, qyTot = 0.0;\r
-\r
- //MC: from reaction plane\r
- if(gAnalysisLevel == "MC"){\r
- if(!event) {\r
- AliError("mcEvent not available");\r
- return 0x0;\r
- }\r
-\r
- AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
- if(gMCEvent){\r
- AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r
- if (headerH) {\r
- gReactionPlane = headerH->ReactionPlaneAngle();\r
- //gReactionPlane *= TMath::RadToDeg();\r
- }//MC header\r
- }//MC event cast\r
- }//MC\r
- \r
- // AOD,ESD,ESDMC: from VZERO Event Plane\r
- else{\r
- \r
- AliEventplane *ep = event->GetEventplane();\r
- if(ep){ \r
- gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
- if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
- //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
- gReactionPlane = gVZEROEventPlane;\r
- }\r
- }//AOD,ESD,ESDMC\r
-\r
- return gReactionPlane;\r
-}\r
-\r
-//________________________________________________________________________\r
-Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
- Double_t vPhi, \r
- Double_t vPt, \r
- Short_t vCharge, \r
- Double_t gCentrality) {\r
- // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
-\r
- Double_t correction = 1.;\r
- Int_t binEta = 0, binPt = 0, binPhi = 0;\r
-\r
- //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
- if(fEtaBinForCorrections != 0) {\r
- Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
- if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
- binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r
- }\r
-\r
- if(fPtBinForCorrections != 0) {\r
- Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
- if(fPtMaxForCorrections != fPtMinForCorrections) \r
- binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r
- }\r
- \r
- if(fPhiBinForCorrections != 0) {\r
- Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
- if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
- binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r
- }\r
-\r
- Int_t gCentralityInt = 1;\r
- for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r
- if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r
- gCentralityInt = i;\r
- break;\r
- }\r
- }\r
- \r
- //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r
-\r
- if(fHistCorrectionPlus[gCentralityInt]){\r
- if (vCharge > 0) {\r
- correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
- //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt); \r
- }\r
- if (vCharge < 0) {\r
- correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
- //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r
- }\r
- }\r
- else {\r
- correction = 1.;\r
- }\r
- if (correction == 0.) { \r
- AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
- correction = 1.; \r
- } \r
- \r
- return correction;\r
-}\r
-\r
-//________________________________________________________________________\r
-TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
- // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
- // Fills QA histograms\r
-\r
- TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
- //output TObjArray holding all good tracks\r
- TObjArray* tracksAccepted = new TObjArray;\r
- tracksAccepted->SetOwner(kTRUE);\r
-\r
- Short_t vCharge;\r
- Double_t vEta;\r
- Double_t vY;\r
- Double_t vPhi;\r
- Double_t vPt;\r
-\r
-\r
- if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
- // Loop over tracks in event\r
- for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
- AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
- if (!aodTrack) {\r
- AliError(Form("Could not receive track %d", iTracks));\r
- continue;\r
- }\r
- \r
- // AOD track cuts\r
- \r
- // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
- // take only TPC only tracks \r
- //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
- for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
- fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
- }\r
- if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
-\r
- //===========================PID (so far only for electron rejection)===============================// \r
- if(fElectronRejection) {\r
- \r
- // get the electron nsigma\r
- Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
- \r
- //Fill QA before the PID\r
- fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
- //end of QA-before pid\r
- \r
- //Make the decision based on the n-sigma\r
- if(nSigma < fElectronRejectionNSigma) continue;\r
- \r
- //Fill QA after the PID\r
- fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
- fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
- }\r
- //===========================end of PID (so far only for electron rejection)===============================//\r
- \r
- vCharge = aodTrack->Charge();\r
- vEta = aodTrack->Eta();\r
- vY = aodTrack->Y();\r
- vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
- vPt = aodTrack->Pt();\r
- \r
- Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
- Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
- \r
- \r
- // Kinematics cuts from ESD track cuts\r
- if( vPt < fPtMin || vPt > fPtMax) continue;\r
- if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
- \r
- // Extra DCA cuts (for systematic studies [!= -1])\r
- if( fDCAxyCut != -1 && fDCAzCut != -1){\r
- if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
- continue; // 2D cut\r
- }\r
- }\r
- \r
- // Extra TPC cuts (for systematic studies [!= -1])\r
- if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
- continue;\r
- }\r
- if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
- continue;\r
- }\r
- \r
- // fill QA histograms\r
- fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
- fHistDCA->Fill(dcaZ,dcaXY);\r
- fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
- fHistPt->Fill(vPt,gCentrality);\r
- fHistEta->Fill(vEta,gCentrality);\r
- fHistRapidity->Fill(vY,gCentrality);\r
- if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
- else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
- fHistPhi->Fill(vPhi,gCentrality);\r
- if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
- else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
- \r
- //=======================================correction\r
- Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
- //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
- \r
- // add the track to the TObjArray\r
- tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
- }//track loop\r
- }// AOD analysis\r
-\r
- //==============================================================================================================\r
- else if(gAnalysisLevel == "MCAOD") {\r
- \r
- AliMCEvent* mcEvent = MCEvent();\r
- if (!mcEvent) {\r
- Printf("ERROR: Could not retrieve MC event");\r
- }\r
- \r
- for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
- AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r
- if (!aodTrack) {\r
- Printf("ERROR: Could not receive track %d (mc loop)", iTracks);\r
- continue;\r
- }\r
- \r
- if(!aodTrack->IsPhysicalPrimary()) continue; \r
-\r
- vCharge = aodTrack->Charge();\r
- vEta = aodTrack->Eta();\r
- vY = aodTrack->Y();\r
- vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
- vPt = aodTrack->Pt();\r
- \r
- // Kinematics cuts from ESD track cuts\r
- if( vPt < fPtMin || vPt > fPtMax) continue;\r
- if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
-\r
- // Remove neutral tracks\r
- if( vCharge == 0 ) continue;\r
-\r
- //Exclude resonances\r
- if(fExcludeResonancesInMC) {\r
- \r
- Bool_t kExcludeParticle = kFALSE;\r
- Int_t gMotherIndex = aodTrack->GetMother();\r
- if(gMotherIndex != -1) {\r
- AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
- if(motherTrack) {\r
- Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
- //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
- //if(pdgCodeOfMother == 113) {\r
- if(pdgCodeOfMother == 113 // rho0\r
- || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
- || pdgCodeOfMother == 221 // eta\r
- || pdgCodeOfMother == 331 // eta'\r
- || pdgCodeOfMother == 223 // omega\r
- || pdgCodeOfMother == 333 // phi\r
- || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r
- || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r
- || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r
- \r
- ) {\r
- kExcludeParticle = kTRUE;\r
- }\r
- }\r
- }\r
- \r
- //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
- if(kExcludeParticle) continue;\r
- }\r
-\r
- // fill QA histograms\r
- fHistPt->Fill(vPt,gCentrality);\r
- fHistEta->Fill(vEta,gCentrality);\r
- fHistRapidity->Fill(vY,gCentrality);\r
- if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
- else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
- fHistPhi->Fill(vPhi,gCentrality);\r
- if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r
- else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
- \r
- //=======================================correction\r
- Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
- //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); \r
- \r
- // add the track to the TObjArray\r
- tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
- }//aodTracks\r
- }//MCAOD\r
- //==============================================================================================================\r
-\r
- else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
-\r
- AliESDtrack *trackTPC = NULL;\r
- AliMCParticle *mcTrack = NULL;\r
-\r
- AliMCEvent* mcEvent = NULL;\r
- \r
- // for MC ESDs use also MC information\r
- if(gAnalysisLevel == "MCESD"){\r
- mcEvent = MCEvent(); \r
- if (!mcEvent) {\r
- AliError("mcEvent not available");\r
- return tracksAccepted;\r
- }\r
- }\r
- \r
- // Loop over tracks in event\r
- for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
- AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
- if (!track) {\r
- AliError(Form("Could not receive track %d", iTracks));\r
- continue;\r
- } \r
-\r
- // for MC ESDs use also MC information --> MC track not used further???\r
- if(gAnalysisLevel == "MCESD"){\r
- Int_t label = TMath::Abs(track->GetLabel());\r
- if(label > mcEvent->GetNumberOfTracks()) continue;\r
- if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
- \r
- mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
- if(!mcTrack) continue;\r
- }\r
-\r
- // take only TPC only tracks\r
- trackTPC = new AliESDtrack();\r
- if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
- \r
- //ESD track cuts\r
- if(fESDtrackCuts) \r
- if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
- \r
- // fill QA histograms\r
- Float_t b[2];\r
- Float_t bCov[3];\r
- trackTPC->GetImpactParameters(b,bCov);\r
- if (bCov[0]<=0 || bCov[2]<=0) {\r
- AliDebug(1, "Estimated b resolution lower or equal zero!");\r
- bCov[0]=0; bCov[2]=0;\r
- }\r
- \r
- Int_t nClustersTPC = -1;\r
- nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone\r
- //nClustersTPC = track->GetTPCclusters(0); // global track\r
- Float_t chi2PerClusterTPC = -1;\r
- if (nClustersTPC!=0) {\r
- chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
- //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
- }\r
- \r
- //===========================PID===============================// \r
- if(fUsePID) {\r
- Double_t prob[AliPID::kSPECIES]={0.};\r
- Double_t probTPC[AliPID::kSPECIES]={0.};\r
- Double_t probTOF[AliPID::kSPECIES]={0.};\r
- Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
- \r
- Double_t nSigma = 0.;\r
- UInt_t detUsedTPC = 0;\r
- UInt_t detUsedTOF = 0;\r
- UInt_t detUsedTPCTOF = 0;\r
- \r
- //Decide what detector configuration we want to use\r
- switch(fPidDetectorConfig) {\r
- case kTPCpid:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
- nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
- detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTPC[iSpecies];\r
- break;\r
- case kTOFpid:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
- nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
- detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTOF[iSpecies];\r
- break;\r
- case kTPCTOF:\r
- fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
- detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
- for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
- prob[iSpecies] = probTPCTOF[iSpecies];\r
- break;\r
- default:\r
- break;\r
- }//end switch: define detector mask\r
- \r
- //Filling the PID QA\r
- Double_t tofTime = -999., length = 999., tof = -999.;\r
- Double_t c = TMath::C()*1.E-9;// m/ns\r
- Double_t beta = -999.;\r
- Double_t nSigmaTOFForParticleOfInterest = -999.;\r
- if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
- (track->IsOn(AliESDtrack::kTIME)) ) { \r
- tofTime = track->GetTOFsignal();//in ps\r
- length = track->GetIntegratedLength();\r
- tof = tofTime*1E-3; // ns \r
- \r
- if (tof <= 0) {\r
- //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
- continue;\r
- }\r
- if (length <= 0){\r
- //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
- continue;\r
- }\r
- \r
- length = length*0.01; // in meters\r
- tof = tof*c;\r
- beta = length/tof;\r
- \r
- nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
- fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
- fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
- fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
- }//TOF signal \r
- \r
- \r
- Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
- fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
- fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
- fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
- fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
- //end of QA-before pid\r
- \r
- if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
- //Make the decision based on the n-sigma\r
- if(fUsePIDnSigma) {\r
- if(nSigma > fPIDNSigma) continue;}\r
- \r
- //Make the decision based on the bayesian\r
- else if(fUsePIDPropabilities) {\r
- if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
- if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
- }\r
- \r
- //Fill QA after the PID\r
- fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
- fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
- fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
- \r
- fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
- fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
- fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
- fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
- }\r
- }\r
- //===========================PID===============================//\r
- vCharge = trackTPC->Charge();\r
- vY = trackTPC->Y();\r
- vEta = trackTPC->Eta();\r
- vPhi = trackTPC->Phi();// * TMath::RadToDeg();\r
- vPt = trackTPC->Pt();\r
- fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
- fHistDCA->Fill(b[1],b[0]);\r
- fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
- fHistPt->Fill(vPt,gCentrality);\r
- fHistEta->Fill(vEta,gCentrality);\r
- fHistPhi->Fill(vPhi,gCentrality);\r
- if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
- else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
- fHistRapidity->Fill(vY,gCentrality);\r
- if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
- else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
- \r
- //=======================================correction\r
- Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
- //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
-\r
- // add the track to the TObjArray\r
- tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
-\r
- delete trackTPC;\r
- }//track loop\r
- }// ESD analysis\r
-\r
- else if(gAnalysisLevel == "MC"){\r
- if(!event) {\r
- AliError("mcEvent not available");\r
- return 0x0;\r
- }\r
-\r
- AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
- if(gMCEvent) {\r
- // Loop over tracks in event\r
- for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
- AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
- if (!track) {\r
- AliError(Form("Could not receive particle %d", iTracks));\r
- continue;\r
- }\r
- \r
- //exclude non stable particles\r
- if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
- \r
- vCharge = track->Charge();\r
- vEta = track->Eta();\r
- vPt = track->Pt();\r
- vY = track->Y();\r
- \r
- if( vPt < fPtMin || vPt > fPtMax) \r
- continue;\r
- if (!fUsePID) {\r
- if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
- }\r
- else if (fUsePID){\r
- if( vY < fEtaMin || vY > fEtaMax) continue;\r
- }\r
-\r
- // Remove neutral tracks\r
- if( vCharge == 0 ) continue;\r
- \r
- //analyze one set of particles\r
- if(fUseMCPdgCode) {\r
- TParticle *particle = track->Particle();\r
- if(!particle) continue;\r
- \r
- Int_t gPdgCode = particle->GetPdgCode();\r
- if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
- continue;\r
- }\r
- \r
- //Use the acceptance parameterization\r
- if(fAcceptanceParameterization) {\r
- Double_t gRandomNumber = gRandom->Rndm();\r
- if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
- continue;\r
- }\r
- \r
- //Exclude resonances\r
- if(fExcludeResonancesInMC) {\r
- TParticle *particle = track->Particle();\r
- if(!particle) continue;\r
- \r
- Bool_t kExcludeParticle = kFALSE;\r
- Int_t gMotherIndex = particle->GetFirstMother();\r
- if(gMotherIndex != -1) {\r
- AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
- if(motherTrack) {\r
- TParticle *motherParticle = motherTrack->Particle();\r
- if(motherParticle) {\r
- Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
- //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
- if(pdgCodeOfMother == 113) {\r
- kExcludeParticle = kTRUE;\r
- }\r
- }\r
- }\r
- }\r
- \r
- //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
- if(kExcludeParticle) continue;\r
- }\r
- \r
- vPhi = track->Phi();\r
- //Printf("phi (before): %lf",vPhi);\r
- \r
- fHistPt->Fill(vPt,gCentrality);\r
- fHistEta->Fill(vEta,gCentrality);\r
- fHistPhi->Fill(vPhi,gCentrality);\r
- if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
- else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
- //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
- fHistRapidity->Fill(vY,gCentrality);\r
- //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
- //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
- if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
- else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
- \r
- //Flow after burner\r
- if(fUseFlowAfterBurner) {\r
- Double_t precisionPhi = 0.001;\r
- Int_t maxNumberOfIterations = 100;\r
- \r
- Double_t phi0 = vPhi;\r
- Double_t gV2 = fDifferentialV2->Eval(vPt);\r
- \r
- for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
- Double_t phiprev = vPhi;\r
- Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
- Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
- vPhi -= fl/fp;\r
- if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
- }\r
- //Printf("phi (after): %lf\n",vPhi);\r
- Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
- if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
- fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
- \r
- Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
- if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
- fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
- \r
- }\r
- \r
- //vPhi *= TMath::RadToDeg();\r
- \r
- //=======================================correction\r
- Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
- //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
-\r
- tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
- } //track loop\r
- }//MC event object\r
- }//MC\r
- \r
- return tracksAccepted; \r
-}\r
-\r
-//________________________________________________________________________\r
-TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
- // Clones TObjArray and returns it with tracks after shuffling the charges\r
-\r
- TObjArray* tracksShuffled = new TObjArray;\r
- tracksShuffled->SetOwner(kTRUE);\r
-\r
- vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks \r
-\r
- for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
- {\r
- AliVParticle* track = (AliVParticle*) tracks->At(i);\r
- chargeVector->push_back(track->Charge());\r
- } \r
- \r
- random_shuffle(chargeVector->begin(), chargeVector->end());\r
- \r
- for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
- AliVParticle* track = (AliVParticle*) tracks->At(i);\r
- //==============================correction\r
- Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
- //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
- tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
- }\r
-\r
- delete chargeVector;\r
- \r
- return tracksShuffled;\r
-}\r
-\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
- //Printf("END BF");\r
-\r
- if (!fBalance) {\r
- AliError("fBalance not available");\r
- return;\r
- } \r
- if(fRunShuffling) {\r
- if (!fShuffledBalance) {\r
- AliError("fShuffledBalance not available");\r
- return;\r
- }\r
- }\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
- // Draw result to the screen\r
- // Called once at the end of the query\r
-\r
- // not implemented ...\r
-\r
-}\r
-\r
+#include "TChain.h"
+#include "TList.h"
+#include "TCanvas.h"
+#include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TH2D.h"
+#include "TH3D.h"
+#include "TArrayF.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TROOT.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+#include "AliCollisionGeometry.h"
+#include "AliGenEventHeader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+#include "AliEventplane.h"
+#include "AliTHn.h"
+#include "AliLog.h"
+#include "AliAnalysisUtils.h"
+
+#include "AliEventPoolManager.h"
+
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+
+#include "AliAnalysisTaskBFPsi.h"
+#include "AliBalancePsi.h"
+#include "AliAnalysisTaskTriggeredBF.h"
+
+
+// Analysis task for the BF vs Psi code
+// Authors: Panos.Christakoglou@nikhef.nl
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskBFPsi)
+
+//________________________________________________________________________
+AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
+: AliAnalysisTaskSE(name),
+ fDebugLevel(kFALSE),
+ fArrayMC(0),
+ fBalance(0),
+ fRunShuffling(kFALSE),
+ fShuffledBalance(0),
+ fRunMixing(kFALSE),
+ fRunMixingEventPlane(kFALSE),
+ fMixingTracks(50000),
+ fMixedBalance(0),
+ fPoolMgr(0),
+ fList(0),
+ fListBF(0),
+ fListBFS(0),
+ fListBFM(0),
+ fHistListPIDQA(0),
+ fHistEventStats(0),
+ fHistCentStats(0),
+ fHistCentStatsUsed(0),
+ fHistTriggerStats(0),
+ fHistTrackStats(0),
+ fHistVx(0),
+ fHistVy(0),
+ fHistVz(0),
+ fHistTPCvsVZEROMultiplicity(0),
+ fHistVZEROSignal(0),
+ fHistEventPlane(0),
+ fHistClus(0),
+ fHistDCA(0),
+ fHistChi2(0),
+ fHistPt(0),
+ fHistEta(0),
+ fHistRapidity(0),
+ fHistPhi(0),
+ fHistEtaPhiPos(0),
+ fHistEtaPhiNeg(0),
+ fHistPhiBefore(0),
+ fHistPhiAfter(0),
+ fHistPhiPos(0),
+ fHistPhiNeg(0),
+ fHistV0M(0),
+ fHistRefTracks(0),
+ fHistdEdxVsPTPCbeforePID(NULL),
+ fHistBetavsPTOFbeforePID(NULL),
+ fHistProbTPCvsPtbeforePID(NULL),
+ fHistProbTOFvsPtbeforePID(NULL),
+ fHistProbTPCTOFvsPtbeforePID(NULL),
+ fHistNSigmaTPCvsPtbeforePID(NULL),
+ fHistNSigmaTOFvsPtbeforePID(NULL),
+ fHistBetaVsdEdXbeforePID(NULL),
+ fHistNSigmaTPCTOFvsPtbeforePID(NULL),
+ fHistNSigmaTPCTOFPbefPID(NULL),
+ fHistdEdxVsPTPCafterPID(NULL),
+ fHistBetavsPTOFafterPID(NULL),
+ fHistProbTPCvsPtafterPID(NULL),
+ fHistProbTOFvsPtafterPID(NULL),
+ fHistProbTPCTOFvsPtafterPID(NULL),
+ fHistNSigmaTPCvsPtafterPID(NULL),
+ fHistNSigmaTOFvsPtafterPID(NULL),
+ fHistBetaVsdEdXafterPID(NULL),
+ fHistNSigmaTPCTOFvsPtafterPID(NULL),
+ fHistNSigmaTPCTOFPafterPID(NULL),
+ fHistdEdxVsPTPCbeforePIDelectron(NULL),
+ fHistNSigmaTPCvsPtbeforePIDelectron(NULL),
+ fHistdEdxVsPTPCafterPIDelectron(NULL),
+ fHistNSigmaTPCvsPtafterPIDelectron(NULL),
+ fCentralityArrayBinsForCorrections(kCENTRALITY),
+ fCentralityWeights(0x0),
+ fPIDResponse(0x0),
+ fPIDCombined(0x0),
+ fParticleOfInterest(kPion),
+ fPidDetectorConfig(kTPCTOF),
+ fUsePID(kFALSE),
+ fUsePIDnSigma(kTRUE),
+ fUsePIDPropabilities(kFALSE),
+ fPIDNSigma(3.),
+ fMinAcceptedPIDProbability(0.8),
+ fElectronRejection(kFALSE),
+ fElectronOnlyRejection(kFALSE),
+ fElectronRejectionNSigma(-1.),
+ fElectronRejectionMinPt(0.),
+ fElectronRejectionMaxPt(1000.),
+ fESDtrackCuts(0),
+ fCentralityEstimator("V0M"),
+ fUseCentrality(kFALSE),
+ fCentralityPercentileMin(0.),
+ fCentralityPercentileMax(5.),
+ fImpactParameterMin(0.),
+ fImpactParameterMax(20.),
+ fMultiplicityEstimator("V0A"),
+ fUseMultiplicity(kFALSE),
+ fNumberOfAcceptedTracksMin(0),
+ fNumberOfAcceptedTracksMax(10000),
+ fHistNumberOfAcceptedTracks(0),
+ fHistMultiplicity(0),
+ fUseOfflineTrigger(kFALSE),
+ fCheckFirstEventInChunk(kFALSE),
+ fCheckPileUp(kFALSE),
+ fCheckPrimaryFlagAOD(kFALSE),
+ fUseMCforKinematics(kFALSE),
+ fVxMax(0.3),
+ fVyMax(0.3),
+ fVzMax(10.),
+ fnAODtrackCutBit(128),
+ fPtMin(0.3),
+ fPtMax(1.5),
+ fEtaMin(-0.8),
+ fEtaMax(0.8),
+ fPhiMin(0.),
+ fPhiMax(360.),
+ fDCAxyCut(-1),
+ fDCAzCut(-1),
+ fTPCchi2Cut(-1),
+ fNClustersTPCCut(-1),
+ fTPCsharedCut(-1),
+ fAcceptanceParameterization(0),
+ fDifferentialV2(0),
+ fUseFlowAfterBurner(kFALSE),
+ fExcludeResonancesInMC(kFALSE),
+ fExcludeElectronsInMC(kFALSE),
+ fUseMCPdgCode(kFALSE),
+ fPDGCodeToBeAnalyzed(-1),
+ fEventClass("EventPlane"),
+ fCustomBinning(""),
+ fHistVZEROAGainEqualizationMap(0),
+ fHistVZEROCGainEqualizationMap(0),
+ fHistVZEROChannelGainEqualizationMap(0) {
+ // Constructor
+ // Define input and output slots here
+ // Input slot #0 works with a TChain
+
+ //======================================================correction
+ for (Int_t i=0; i<kCENTRALITY; i++){
+ fHistCorrectionPlus[i] = NULL;
+ fHistCorrectionMinus[i] = NULL;
+ fCentralityArrayForCorrections[i] = -1.;
+ }
+ //=====================================================correction
+
+ DefineInput(0, TChain::Class());
+ // Output slot #0 writes into a TH1 container
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
+ DefineOutput(3, TList::Class());
+ DefineOutput(4, TList::Class());
+ DefineOutput(5, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
+
+ // delete fBalance;
+ // delete fShuffledBalance;
+ // delete fList;
+ // delete fListBF;
+ // delete fListBFS;
+
+ // delete fHistEventStats;
+ // delete fHistTrackStats;
+ // delete fHistVx;
+ // delete fHistVy;
+ // delete fHistVz;
+
+ // delete fHistClus;
+ // delete fHistDCA;
+ // delete fHistChi2;
+ // delete fHistPt;
+ // delete fHistEta;
+ // delete fHistPhi;
+ // delete fHistEtaPhiPos;
+ // delete fHistEtaPhiNeg;
+ // delete fHistV0M;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
+ // Create histograms
+ // Called once
+
+ // global switch disabling the reference
+ // (to avoid "Replacing existing TH1" if several wagons are created in train)
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ if(!fBalance) {
+ fBalance = new AliBalancePsi();
+ fBalance->SetAnalysisLevel("ESD");
+ fBalance->SetEventClass(fEventClass);
+ //fBalance->SetNumberOfBins(-1,16);
+ //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
+ }
+ if(fRunShuffling) {
+ if(!fShuffledBalance) {
+ fShuffledBalance = new AliBalancePsi();
+ fShuffledBalance->SetAnalysisLevel("ESD");
+ fShuffledBalance->SetEventClass(fEventClass);
+ //fShuffledBalance->SetNumberOfBins(-1,16);
+ //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
+ }
+ }
+ if(fRunMixing) {
+ if(!fMixedBalance) {
+ fMixedBalance = new AliBalancePsi();
+ fMixedBalance->SetAnalysisLevel("ESD");
+ fMixedBalance->SetEventClass(fEventClass);
+ }
+ }
+
+ //QA list
+ fList = new TList();
+ fList->SetName("listQA");
+ fList->SetOwner();
+
+ //Balance Function list
+ fListBF = new TList();
+ fListBF->SetName("listBF");
+ fListBF->SetOwner();
+
+ if(fRunShuffling) {
+ fListBFS = new TList();
+ fListBFS->SetName("listBFShuffled");
+ fListBFS->SetOwner();
+ }
+
+ if(fRunMixing) {
+ fListBFM = new TList();
+ fListBFM->SetName("listTriggeredBFMixed");
+ fListBFM->SetOwner();
+ }
+
+ //PID QA list
+ if(fUsePID || fElectronRejection) {
+ fHistListPIDQA = new TList();
+ fHistListPIDQA->SetName("listQAPID");
+ fHistListPIDQA->SetOwner();
+ }
+
+ //Event stats.
+ TString gCutName[7] = {"Total","Offline trigger",
+ "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
+ fHistEventStats = new TH2F("fHistEventStats",
+ "Event statistics;;Centrality percentile;N_{events}",
+ 7,0.5,7.5,220,-5,105);
+ for(Int_t i = 1; i <= 7; i++)
+ fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+ fList->Add(fHistEventStats);
+
+ TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
+ fHistCentStats = new TH2F("fHistCentStats",
+ "Centrality statistics;;Cent percentile",
+ 13,-0.5,12.5,220,-5,105);
+ for(Int_t i = 1; i <= 13; i++){
+ fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
+ //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
+ }
+ fList->Add(fHistCentStats);
+
+ fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
+ fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
+ fList->Add(fHistCentStatsUsed);
+
+ fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
+ fList->Add(fHistTriggerStats);
+
+ fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
+ fList->Add(fHistTrackStats);
+
+ fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
+ fList->Add(fHistNumberOfAcceptedTracks);
+
+ fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
+ fList->Add(fHistMultiplicity);
+
+ // Vertex distributions
+ fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
+ fList->Add(fHistVx);
+ fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
+ fList->Add(fHistVy);
+ fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
+ fList->Add(fHistVz);
+
+ //TPC vs VZERO multiplicity
+ fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
+ if(fMultiplicityEstimator == "V0A")
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
+ else if(fMultiplicityEstimator == "V0C")
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
+ else
+ fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
+ fList->Add(fHistTPCvsVZEROMultiplicity);
+
+ fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
+ fList->Add(fHistVZEROSignal);
+
+ //Event plane
+ fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
+ fList->Add(fHistEventPlane);
+
+ // QA histograms
+ fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
+ fList->Add(fHistClus);
+ fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
+ fList->Add(fHistChi2);
+ fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
+ fList->Add(fHistDCA);
+ fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
+ fList->Add(fHistPt);
+ fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
+ fList->Add(fHistEta);
+ fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
+ fList->Add(fHistRapidity);
+ fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
+ fList->Add(fHistPhi);
+ fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",40,-1.6,1.6,72,0.,2.*TMath::Pi(),220,-5,105);
+ fList->Add(fHistEtaPhiPos);
+ fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",40,-1.6,1.6,72,0.,2.*TMath::Pi(),220,-5,105);
+ fList->Add(fHistEtaPhiNeg);
+ fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
+ fList->Add(fHistPhiBefore);
+ fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
+ fList->Add(fHistPhiAfter);
+ fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
+ fList->Add(fHistPhiPos);
+ fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
+ fList->Add(fHistPhiNeg);
+ fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
+ fList->Add(fHistV0M);
+ TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
+ fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
+ for(Int_t i = 1; i <= 6; i++)
+ fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
+ fList->Add(fHistRefTracks);
+
+ // Balance function histograms
+ // Initialize histograms if not done yet (including the custom binning)
+ if(!fBalance->GetHistNp()){
+ AliInfo("Histograms not yet initialized! --> Will be done now");
+ fBalance->SetCustomBinning(fCustomBinning);
+ fBalance->InitHistograms();
+ }
+
+ if(fRunShuffling) {
+ if(!fShuffledBalance->GetHistNp()) {
+ AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
+ fShuffledBalance->SetCustomBinning(fCustomBinning);
+ fShuffledBalance->InitHistograms();
+ }
+ }
+
+ if(fRunMixing) {
+ if(!fMixedBalance->GetHistNp()) {
+ AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
+ fMixedBalance->SetCustomBinning(fCustomBinning);
+ fMixedBalance->InitHistograms();
+ }
+ }
+
+ // QA histograms for different cuts
+ fList->Add(fBalance->GetQAHistHBTbefore());
+ fList->Add(fBalance->GetQAHistHBTafter());
+ fList->Add(fBalance->GetQAHistConversionbefore());
+ fList->Add(fBalance->GetQAHistConversionafter());
+ fList->Add(fBalance->GetQAHistPsiMinusPhi());
+ fList->Add(fBalance->GetQAHistResonancesBefore());
+ fList->Add(fBalance->GetQAHistResonancesRho());
+ fList->Add(fBalance->GetQAHistResonancesK0());
+ fList->Add(fBalance->GetQAHistResonancesLambda());
+ fList->Add(fBalance->GetQAHistQbefore());
+ fList->Add(fBalance->GetQAHistQafter());
+
+ //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
+ fListBF->Add(fBalance->GetHistNp());
+ fListBF->Add(fBalance->GetHistNn());
+ fListBF->Add(fBalance->GetHistNpn());
+ fListBF->Add(fBalance->GetHistNnn());
+ fListBF->Add(fBalance->GetHistNpp());
+ fListBF->Add(fBalance->GetHistNnp());
+
+ if(fRunShuffling) {
+ fListBFS->Add(fShuffledBalance->GetHistNp());
+ fListBFS->Add(fShuffledBalance->GetHistNn());
+ fListBFS->Add(fShuffledBalance->GetHistNpn());
+ fListBFS->Add(fShuffledBalance->GetHistNnn());
+ fListBFS->Add(fShuffledBalance->GetHistNpp());
+ fListBFS->Add(fShuffledBalance->GetHistNnp());
+ }
+
+ if(fRunMixing) {
+ fListBFM->Add(fMixedBalance->GetHistNp());
+ fListBFM->Add(fMixedBalance->GetHistNn());
+ fListBFM->Add(fMixedBalance->GetHistNpn());
+ fListBFM->Add(fMixedBalance->GetHistNnn());
+ fListBFM->Add(fMixedBalance->GetHistNpp());
+ fListBFM->Add(fMixedBalance->GetHistNnp());
+ }
+ //}
+
+
+ // Event Mixing
+ if(fRunMixing){
+ Int_t trackDepth = fMixingTracks;
+ Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
+
+ // centrality bins
+ Double_t* centbins = NULL;
+ Int_t nCentralityBins;
+ if(fBalance->IsUseVertexBinning()){
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
+ }
+ else{
+ centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
+ }
+
+ // multiplicity bins
+ Double_t* multbins = NULL;
+ Int_t nMultiplicityBins;
+ multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
+
+ // Zvtx bins
+ Double_t* vtxbins = NULL;
+ Int_t nVertexBins;
+ if(fBalance->IsUseVertexBinning()){
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
+ }
+ else{
+ vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
+ }
+
+ // Event plane angle (Psi) bins
+ Double_t* psibins = NULL;
+ Int_t nPsiBins;
+ psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
+
+
+ // run the event mixing also in bins of event plane (statistics!)
+ if(fRunMixingEventPlane){
+ if(fEventClass=="Multiplicity"){
+ if(multbins && vtxbins && psibins){
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
+ }
+ }
+ else{
+ if(centbins && vtxbins && psibins){
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
+ }
+ }
+ }
+ else{
+ if(fEventClass=="Multiplicity"){
+ if(multbins && vtxbins){
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
+ }
+ }
+ else{
+ if(centbins && vtxbins){
+ fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
+ }
+ }
+ }
+
+ if(centbins) delete [] centbins;
+ if(multbins) delete [] multbins;
+ if(vtxbins) delete [] vtxbins;
+ if(psibins) delete [] psibins;
+
+ // set minimum values for track depth, fraction, and number of events
+ fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
+
+ // check pool manager
+ if(!fPoolMgr){
+ AliError("Event Mixing required, but Pool Manager not initialized...");
+ return;
+ }
+ }
+
+ if(fESDtrackCuts) fList->Add(fESDtrackCuts);
+
+ //====================PID========================//
+ if(fUsePID) {
+ fPIDCombined = new AliPIDCombined();
+ fPIDCombined->SetDefaultTPCPriors();
+
+ fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
+
+ fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
+
+ fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
+
+ fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
+
+ fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
+
+ fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
+
+ fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
+
+ fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID);
+
+ fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID);
+
+ //+++++++++++++++++//
+ //p array
+ const Int_t pBins = 36;
+ Double_t nArrayP[pBins+1]={0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 15.0, 20.0};
+ //nSigma Array
+ const Int_t nSigmaBins = 250;
+ Double_t nArrayS[nSigmaBins+1];
+ for (Int_t i = 0; i <= nSigmaBins; i++){
+ nArrayS[i]=i-125; //i+1
+ //Printf("nS: %lf - i: %d", nSigmaArray[i], i);
+ }
+
+ fHistNSigmaTPCTOFPbefPID = new TH3D ("fHistNSigmaTPCTOFPbefPID","fHistNSigmaTPCTOFPbefPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP);
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFPbefPID);
+ //++++++++++++++//
+
+ fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
+
+ fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
+
+ fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
+ fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
+
+ fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
+ fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
+
+ fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
+ fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
+
+ fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
+
+ fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
+
+ fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
+ fHistListPIDQA->Add(fHistBetaVsdEdXafterPID);
+
+ fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, -25, 25);
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID);
+
+ fHistNSigmaTPCTOFPafterPID = new TH3D ("fHistNSigmaTPCTOFPafterPID","fHistNSigmaTPCTOFPafterPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP);
+ fHistListPIDQA->Add(fHistNSigmaTPCTOFPafterPID); //++++++++++++++
+ }
+
+ // for electron rejection only TPC nsigma histograms
+ if(fElectronRejection) {
+
+ fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
+
+ fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
+
+ fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
+ fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
+
+ fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
+ fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
+ }
+ //====================PID========================//
+
+ // Post output data.
+ PostData(1, fList);
+ PostData(2, fListBF);
+ if(fRunShuffling) PostData(3, fListBFS);
+ if(fRunMixing) PostData(4, fListBFM);
+ if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
+
+ AliInfo("Finished setting up the Output");
+
+ TH1::AddDirectory(oldStatus);
+
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
+ Int_t nCentralityBins,
+ Double_t *centralityArrayForCorrections) {
+ //Open files that will be used for correction
+ fCentralityArrayBinsForCorrections = nCentralityBins;
+ for (Int_t i=0; i<nCentralityBins; i++)
+ fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
+
+ // No file specified -> run without corrections
+ if(!filename.Contains(".root")) {
+ AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
+ return;
+ }
+
+ //Open the input file
+ TFile *f = TFile::Open(filename);
+ if(!f->IsOpen()) {
+ AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
+ return;
+ }
+
+ //TString listEffName = "";
+ for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
+ //Printf("iCent %d:",iCent);
+ TString histoName = "fHistCorrectionPlus";
+ histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
+ fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
+ if(!fHistCorrectionPlus[iCent]) {
+ AliError(Form("fHist %s not found!!!",histoName.Data()));
+ return;
+ }
+
+ histoName = "fHistCorrectionMinus";
+ histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
+ fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
+ if(!fHistCorrectionMinus[iCent]) {
+ AliError(Form("fHist %s not found!!!",histoName.Data()));
+ return;
+ }
+ }//loop over centralities: ONLY the PbPb case is covered
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
+ // Main loop
+ // Called for each event
+
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+ Int_t gNumberOfAcceptedTracks = 0;
+ Double_t lMultiplicityVar = -999.; //-1
+ Double_t gReactionPlane = -1.;
+ Float_t bSign = 0.;
+
+ // get the event (for generator level: MCEvent())
+ AliVEvent* eventMain = NULL;
+ if(gAnalysisLevel == "MC") {
+ eventMain = dynamic_cast<AliVEvent*>(MCEvent());
+ }
+ else{
+ eventMain = dynamic_cast<AliVEvent*>(InputEvent());
+ // for HBT like cuts need magnetic field sign
+ bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
+ }
+ if(!eventMain) {
+ AliError("eventMain not available");
+ return;
+ }
+
+ // PID Response task active?
+ if(fUsePID || fElectronRejection) {
+ fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
+ if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
+ }
+
+ // check event cuts and fill event histograms
+ if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
+ return;
+ }
+ // get the reaction plane
+ if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") {
+ gReactionPlane = GetEventPlane(eventMain);
+ fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
+ if(gReactionPlane < 0){
+ return;
+ }
+ }
+
+ // get the accepted tracks in main event
+ TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
+ gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
+
+ //multiplicity cut (used in pp)
+ fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
+
+ // store charges of all accepted tracks,shuffle and reassign(two extra loops)
+ TObjArray* tracksShuffled = NULL;
+ if(fRunShuffling){
+ tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
+ }
+
+ // Event mixing
+ if (fRunMixing)
+ {
+ // 1. First get an event pool corresponding in mult (cent) and
+ // zvertex to the current event. Once initialized, the pool
+ // should contain nMix (reduced) events. This routine does not
+ // pre-scan the chain. The first several events of every chain
+ // will be skipped until the needed pools are filled to the
+ // specified depth. If the pool categories are not too rare, this
+ // should not be a problem. If they are rare, you could lose`
+ // statistics.
+
+ // 2. Collect the whole pool's content of tracks into one TObjArray
+ // (bgTracks), which is effectively a single background super-event.
+
+ // 3. The reduced and bgTracks arrays must both be passed into
+ // FillCorrelations(). Also nMix should be passed in, so a weight
+ // of 1./nMix can be applied.
+
+ AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
+
+ if (!pool){
+ AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
+ }
+ else{
+
+ //pool->SetDebug(1);
+
+ if (pool->IsReady()){
+
+
+ Int_t nMix = pool->GetCurrentNEvents();
+ //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
+
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
+ //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
+ //if (pool->IsReady())
+ //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
+
+ // Fill mixed-event histos here
+ for (Int_t jMix=0; jMix<nMix; jMix++)
+ {
+ TObjArray* tracksMixed = pool->GetEvent(jMix);
+ fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
+ }
+ }
+
+ // Update the Event pool
+ pool->UpdatePool(tracksMain);
+ //pool->PrintInfo();
+
+ }//pool NULL check
+ }//run mixing
+
+ // calculate balance function
+ fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
+
+ // calculate shuffled balance function
+ if(fRunShuffling && tracksShuffled != NULL) {
+ fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
+ }
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
+ // Checks the Event cuts
+ // Fills Event statistics histograms
+
+ Bool_t isSelectedMain = kTRUE;
+ Float_t gRefMultiplicity = -1.;
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+ AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
+
+ fHistEventStats->Fill(1,gRefMultiplicity); //all events
+
+ // check first event in chunk (is not needed for new reconstructions)
+ if(fCheckFirstEventInChunk){
+ AliAnalysisUtils ut;
+ if(ut.IsFirstEventInChunk(event))
+ return -1.;
+ fHistEventStats->Fill(6,gRefMultiplicity);
+ }
+ // check for pile-up event
+ if(fCheckPileUp){
+ AliAnalysisUtils ut;
+ ut.SetUseMVPlpSelection(kTRUE);
+ ut.SetUseOutOfBunchPileUp(kTRUE);
+ if(ut.IsPileUpEvent(event))
+ return -1.;
+ fHistEventStats->Fill(7,gRefMultiplicity);
+ }
+
+ // Event trigger bits
+ fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
+ if(fUseOfflineTrigger)
+ isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+
+ if(isSelectedMain) {
+ fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
+
+ // Event Vertex MC
+ if(gAnalysisLevel == "MC") {
+ if(!event) {
+ AliError("mcEvent not available");
+ return 0x0;
+ }
+
+ if(mcevent){
+ AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
+ if(header){
+ TArrayF gVertexArray;
+ header->PrimaryVertex(gVertexArray);
+ //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
+ //gVertexArray.At(0),
+ //gVertexArray.At(1),
+ //gVertexArray.At(2));
+ fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
+ if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
+ if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
+ if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
+
+ // get the reference multiplicty or centrality
+ gRefMultiplicity = GetRefMultiOrCentrality(event);
+
+ fHistVx->Fill(gVertexArray.At(0));
+ fHistVy->Fill(gVertexArray.At(1));
+ fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
+
+ // take only events inside centrality class
+ if(fUseCentrality) {
+ if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
+ return gRefMultiplicity;
+ }//centrality class
+ }
+ // take events only within the same multiplicity class
+ else if(fUseMultiplicity) {
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
+ return gRefMultiplicity;
+ }
+ }//multiplicity range
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//header
+ }//MC event object
+ }//MC
+
+ // Event Vertex AOD, ESD, ESDMC
+ else{
+ const AliVVertex *vertex = event->GetPrimaryVertex();
+
+ if(vertex) {
+ Double32_t fCov[6];
+ vertex->GetCovarianceMatrix(fCov);
+ if(vertex->GetNContributors() > 0) {
+ if(fCov[5] != 0) {
+ fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
+ if(TMath::Abs(vertex->GetX()) < fVxMax) {
+ if(TMath::Abs(vertex->GetY()) < fVyMax) {
+ if(TMath::Abs(vertex->GetZ()) < fVzMax) {
+ fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
+
+ // get the reference multiplicty or centrality
+ gRefMultiplicity = GetRefMultiOrCentrality(event);
+
+ fHistVx->Fill(vertex->GetX());
+ fHistVy->Fill(vertex->GetY());
+ fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
+
+ // take only events inside centrality class
+ if(fUseCentrality) {
+ if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
+
+ // centrality weighting (optional for 2011 if central and semicentral triggers are used)
+ if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
+ AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
+ return -1;
+ }
+
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
+ return gRefMultiplicity;
+ }//centrality class
+ }
+ // take events only within the same multiplicity class
+ else if(fUseMultiplicity) {
+ //if(fDebugLevel)
+ //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
+ //fNumberOfAcceptedTracksMax,gRefMultiplicity);
+
+ if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
+ fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
+ return gRefMultiplicity;
+ }
+ }//multiplicity range
+ }//Vz cut
+ }//Vy cut
+ }//Vx cut
+ }//proper vertex resolution
+ }//proper number of contributors
+ }//vertex object valid
+ }//triggered event
+ }//AOD,ESD,ESDMC
+
+ // in all other cases return -1 (event not accepted)
+ return -1;
+}
+
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
+ // Checks the Event cuts
+ // Fills Event statistics histograms
+
+ Float_t gCentrality = -1.;
+ Double_t gMultiplicity = -1.;
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+
+ // calculate centrality always (not only in centrality mode)
+ if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();
+ if(header){
+ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // QA for centrality estimators
+ fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
+ fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
+ fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
+ fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
+ fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
+ fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
+ fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
+ fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
+ fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
+ fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
+ fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
+ fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
+ fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
+
+ // Centrality estimator USED ++++++++++++++++++++++++++++++
+ fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
+
+ // centrality QA (reference tracks)
+ fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
+ fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
+ fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
+ fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
+ fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
+ fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
+ fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
+ fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
+ fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
+
+ }//AOD header
+ }//AOD
+
+ // calculate centrality always (not only in centrality mode)
+ else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround
+
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();
+ if(header){
+ gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0;
+
+ // QA histogram
+ fHistCentStatsUsed->Fill(0.,gCentrality);
+
+ }//AOD header
+ }//AODnano
+
+ else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
+ AliCentrality *centrality = event->GetCentrality();
+ gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // QA for centrality estimators
+ fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
+ fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
+ fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
+ fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
+ fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
+ fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
+ fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
+ fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
+ fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
+ fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
+ fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
+ fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
+ fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
+
+ // Centrality estimator USED ++++++++++++++++++++++++++++++
+ fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
+
+ // centrality QA (V0M)
+ fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
+ }//ESD
+
+ else if(gAnalysisLevel == "MC"){
+ Double_t gImpactParameter = 0.;
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
+ if(gMCEvent){
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
+ if(headerH){
+ gImpactParameter = headerH->ImpactParameter();
+ gCentrality = gImpactParameter;
+ }//MC header
+ }//MC event cast
+ }//MC
+
+ else{
+ gCentrality = -1.;
+ }
+
+ // calculate reference multiplicity always (not only in multiplicity mode)
+ if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
+ AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
+ if(gESDEvent){
+ gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
+ fHistMultiplicity->Fill(gMultiplicity);
+ }//AliESDevent cast
+ }//ESD mode
+
+ else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
+ AliAODHeader *header = (AliAODHeader*) event->GetHeader();
+ if ((fMultiplicityEstimator == "V0M")||
+ (fMultiplicityEstimator == "V0A")||
+ (fMultiplicityEstimator == "V0C") ||
+ (fMultiplicityEstimator == "TPC")) {
+ gMultiplicity = GetReferenceMultiplicityFromAOD(event);
+ if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
+ }
+ else {
+ if(header)
+ gMultiplicity = header->GetRefMultiplicity();
+ if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
+ }
+ fHistMultiplicity->Fill(gMultiplicity);
+ }//AOD mode
+ else if(gAnalysisLevel == "MC") {
+ AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
+ //Calculating the multiplicity as the number of charged primaries
+ //within \pm 0.8 in eta and pT > 0.1 GeV/c
+ for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
+ if (!track) {
+ AliError(Form("Could not receive particle %d", iParticle));
+ continue;
+ }
+
+ //exclude non stable particles
+ if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
+
+ //++++++++++++++++
+ if (fMultiplicityEstimator == "V0M") {
+ if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
+ continue;}
+ else if (fMultiplicityEstimator == "V0A") {
+ if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
+ else if (fMultiplicityEstimator == "V0C") {
+ if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
+ else if (fMultiplicityEstimator == "TPC") {
+ if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
+ }
+ else{
+ if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
+ if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
+ }
+ //++++++++++++++++
+
+ if(track->Charge() == 0) continue;
+
+ gMultiplicity += 1;
+ }//loop over primaries
+ fHistMultiplicity->Fill(gMultiplicity);
+ }//MC mode
+ else{
+ gMultiplicity = -1;
+ }
+
+
+ // decide what should be returned only here
+ Double_t lReturnVal = -100;
+ if(fEventClass=="Multiplicity"){
+ lReturnVal = gMultiplicity;
+ }else if(fEventClass=="Centrality"){
+ lReturnVal = gCentrality;
+ }
+ return lReturnVal;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
+ //Function that returns the reference multiplicity from AODs (data or reco MC)
+ //Different ref. mult. implemented: V0M, V0A, V0C, TPC
+ Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
+ Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
+
+ AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
+ if(!header) {
+ Printf("ERROR: AOD header not available");
+ return -999;
+ }
+ Int_t gRunNumber = header->GetRunNumber();
+
+ // Loop over tracks in event
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
+ if (!aodTrack) {
+ AliError(Form("Could not receive track %d", iTracks));
+ continue;
+ }
+
+ // AOD track cuts
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
+
+ if(aodTrack->Charge() == 0) continue;
+ // Kinematics cuts from ESD track cuts
+ if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
+ if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
+
+ //=================PID (so far only for electron rejection)==========================//
+ if(fElectronRejection) {
+ // get the electron nsigma
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
+
+ // check only for given momentum range
+ if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
+ //look only at electron nsigma
+ if(!fElectronOnlyRejection) {
+ //Make the decision based on the n-sigma of electrons
+ if(nSigma < fElectronRejectionNSigma) continue;
+ }
+ else {
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
+
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
+ if(nSigma < fElectronRejectionNSigma
+ && nSigmaPions > fElectronRejectionNSigma
+ && nSigmaKaons > fElectronRejectionNSigma
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;
+ }
+ }
+ }//electron rejection
+
+ gRefMultiplicityTPC += 1.0;
+ }// track loop
+
+ //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
+ for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
+ fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
+
+ if(iChannel < 32)
+ gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
+ else if(iChannel >= 32)
+ gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
+ }//loop over PMTs
+
+ //Equalization of gain
+ Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
+ if(gFactorA != 0)
+ gRefMultiplicityVZEROA /= gFactorA;
+ Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
+ if(gFactorC != 0)
+ gRefMultiplicityVZEROC /= gFactorC;
+ if((gFactorA != 0)&&(gFactorC != 0))
+ gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
+
+ if(fDebugLevel)
+ Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
+
+ fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
+
+ if(fMultiplicityEstimator == "TPC")
+ gRefMultiplicity = gRefMultiplicityTPC;
+ else if(fMultiplicityEstimator == "V0M")
+ gRefMultiplicity = gRefMultiplicityVZERO;
+ else if(fMultiplicityEstimator == "V0A")
+ gRefMultiplicity = gRefMultiplicityVZEROA;
+ else if(fMultiplicityEstimator == "V0C")
+ gRefMultiplicity = gRefMultiplicityVZEROC;
+
+ return gRefMultiplicity;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
+ // Get the event plane
+
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+ Float_t gVZEROEventPlane = -10.;
+ Float_t gReactionPlane = -10.;
+ Double_t qxTot = 0.0, qyTot = 0.0;
+
+ //MC: from reaction plane
+ if(gAnalysisLevel == "MC"){
+ if(!event) {
+ AliError("mcEvent not available");
+ return 0x0;
+ }
+
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
+ if(gMCEvent){
+ AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
+ if (headerH) {
+ gReactionPlane = headerH->ReactionPlaneAngle();
+ //gReactionPlane *= TMath::RadToDeg();
+ }//MC header
+ }//MC event cast
+ }//MC
+
+ // AOD,ESD,ESDMC: from VZERO Event Plane
+ else{
+
+ AliEventplane *ep = event->GetEventplane();
+ if(ep){
+ gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
+ if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
+ //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
+ gReactionPlane = gVZEROEventPlane;
+ }
+ }//AOD,ESD,ESDMC
+
+ return gReactionPlane;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
+ Double_t vPhi,
+ Double_t vPt,
+ Short_t vCharge,
+ Double_t gCentrality) {
+ // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
+
+ Double_t correction = 1.;
+ Int_t gCentralityInt = -1;
+
+ for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
+ if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
+ gCentralityInt = i;
+ break;
+ }
+ }
+
+ // centrality not in array --> no correction
+ if(gCentralityInt < 0){
+ correction = 1.;
+ }
+ else{
+
+ //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
+
+ if(fHistCorrectionPlus[gCentralityInt]){
+ if (vCharge > 0) {
+ correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
+ //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
+ }
+ if (vCharge < 0) {
+ correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
+ //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
+ }
+ }
+ else {
+ correction = 1.;
+ }
+ }//centrality in array
+
+ if (correction == 0.) {
+ AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
+ correction = 1.;
+ }
+
+ return correction;
+}
+
+//________________________________________________________________________
+TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
+ // Returns TObjArray with tracks after all track cuts (only for AOD!)
+ // Fills QA histograms
+
+ TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+
+ //output TObjArray holding all good tracks
+ TObjArray* tracksAccepted = new TObjArray;
+ tracksAccepted->SetOwner(kTRUE);
+
+ Short_t vCharge;
+ Double_t vEta;
+ Double_t vY;
+ Double_t vPhi;
+ Double_t vPt;
+
+
+ if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
+ // Loop over tracks in event
+
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
+ if (!aodTrack) {
+ AliError(Form("Could not receive track %d", iTracks));
+ continue;
+ }
+
+ // AOD track cuts
+
+ // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
+ // take only TPC only tracks
+ //fHistTrackStats->Fill(aodTrack->GetFilterMap());
+ for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
+ fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
+ }
+
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
+
+
+ // additional check on kPrimary flag
+ if(fCheckPrimaryFlagAOD){
+ if(aodTrack->GetType() != AliAODTrack::kPrimary)
+ continue;
+ }
+
+
+ vCharge = aodTrack->Charge();
+ vEta = aodTrack->Eta();
+ vY = aodTrack->Y();
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();
+ vPt = aodTrack->Pt();
+
+ //===========================PID (so far only for electron rejection)===============================//
+ if(fElectronRejection) {
+
+ // get the electron nsigma
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
+
+ //Fill QA before the PID
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
+ //end of QA-before pid
+
+ // check only for given momentum range
+ if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
+
+ //look only at electron nsigma
+ if(!fElectronOnlyRejection){
+
+ //Make the decision based on the n-sigma of electrons
+ if(nSigma < fElectronRejectionNSigma) continue;
+ }
+ else{
+
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
+
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
+ if(nSigma < fElectronRejectionNSigma
+ && nSigmaPions > fElectronRejectionNSigma
+ && nSigmaKaons > fElectronRejectionNSigma
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;
+ }
+ }
+
+ //Fill QA after the PID
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
+
+ }
+ //===========================end of PID (so far only for electron rejection)===============================//
+
+ //+++++++++++++++++++++++++++++//
+ //===========================PID===============================//
+ if(fUsePID) {
+ Double_t prob[AliPID::kSPECIES]={0.};
+ Double_t probTPC[AliPID::kSPECIES]={0.};
+ Double_t probTOF[AliPID::kSPECIES]={0.};
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+
+ Double_t nSigma = 0.;
+ Double_t nSigmaTPC = 0.;
+ Double_t nSigmaTOF = 0.;
+ Double_t nSigmaTPCTOF = 0.;
+ UInt_t detUsedTPC = 0;
+ UInt_t detUsedTOF = 0;
+ UInt_t detUsedTPCTOF = 0;
+
+ nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
+ nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
+ nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
+ if (nSigmaTOF == 999 || nSigmaTOF == -999){
+ nSigmaTPCTOF = nSigmaTPC;
+ }
+
+ //Decide what detector configuration we want to use
+ switch(fPidDetectorConfig) {
+ case kTPCpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+ nSigma = nSigmaTPC;
+ detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPC[iSpecies];
+ break;
+ case kTOFpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+ nSigma = nSigmaTOF;
+ detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTOF[iSpecies];
+ break;
+ case kTPCTOF:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
+ nSigma = nSigmaTPCTOF;
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPCTOF[iSpecies];
+ break;
+ default:
+ break;
+ }//end switch: define detector mask
+
+ //Filling the PID QA
+ Double_t tofTime = -999., length = 999., tof = -999.;
+ Double_t c = TMath::C()*1.E-9;// m/ns
+ Double_t beta = -999.;
+
+ Float_t probMis = fPIDResponse->GetTOFMismatchProbability(aodTrack);
+ if (probMis < 0.01) { //if u want to reduce mismatch using also TPC
+
+ if ((aodTrack->IsOn(AliAODTrack::kITSin)) && (aodTrack->IsOn(AliAODTrack::kTOFpid)) ) { //leonardo's analysis
+
+ tofTime = aodTrack->GetTOFsignal();//in ps
+ length = aodTrack->GetIntegratedLength();
+ tof = tofTime*1E-3; // ns
+ if (tof <= 0) {
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
+ continue;
+ }
+ if (length <= 0){
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
+ continue;
+ }
+ length = length*0.01; // in meters
+ tof = tof*c;
+ beta = length/tof;
+
+ fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
+ fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
+ }//TOF signal
+
+ fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->GetTPCmomentum()*aodTrack->Charge(),aodTrack->GetTPCsignal()); //aodTrack->P()*aodTrack->Charge()
+ fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
+ fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
+
+ //combined TPC&TOF
+ fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta);
+ fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
+ fHistNSigmaTPCTOFPbefPID ->Fill(nSigmaTPC,nSigmaTOF,aodTrack->P()); //++++++++++++++
+ //Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
+ //end of QA-before pid
+
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
+ //Make the decision based on the n-sigma
+ if(fUsePIDnSigma) {
+ if(nSigma > fPIDNSigma || nSigma < -fPIDNSigma) continue;
+
+ fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
+ fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
+ fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
+ fHistNSigmaTPCTOFPafterPID ->Fill(nSigmaTPC,nSigmaTOF,aodTrack->P()); //++++++++++++++
+ }
+ //Make the decision based on the bayesian
+ else if(fUsePIDPropabilities) {
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
+
+ fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
+ fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
+ fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
+ }
+
+ //Fill QA after the PID
+ fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
+ fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
+ fHistBetaVsdEdXafterPID ->Fill(aodTrack->GetTPCsignal(),beta);
+ }
+ }
+ }
+ //===========================PID===============================//
+ //+++++++++++++++++++++++++++++//
+
+
+ Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
+ Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
+
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+ // Extra DCA cuts (for systematic studies [!= -1])
+ if( fDCAxyCut != -1 && fDCAzCut != -1){
+ if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
+ continue; // 2D cut
+ }
+ }
+
+ // Extra TPC cuts (for systematic studies [!= -1])
+ if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
+ continue;
+ }
+ if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
+ continue;
+ }
+
+ // Extra cut on shared clusters
+ if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
+ continue;
+ }
+
+ // fill QA histograms
+ fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
+ fHistDCA->Fill(dcaZ,dcaXY);
+ fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ // add the track to the TObjArray
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+ }//track loop
+ }// AOD analysis
+
+
+ // nano AODs
+ else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing)
+ // Loop over tracks in event
+
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+ AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks));
+ if (!aodTrack) {
+ AliError(Form("Could not receive track %d", iTracks));
+ continue;
+ }
+
+ // AOD track cuts (not needed)
+ //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
+
+ vCharge = aodTrack->Charge();
+ vEta = aodTrack->Eta();
+ vY = -999.;
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();
+ vPt = aodTrack->Pt();
+
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+
+ // fill QA histograms
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ // add the track to the TObjArray
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+ }//track loop
+ }// AOD nano analysis
+
+
+ //==============================================================================================================
+ else if(gAnalysisLevel == "MCAOD") {
+
+ AliMCEvent* mcEvent = MCEvent();
+ if (!mcEvent) {
+ AliError("ERROR: Could not retrieve MC event");
+ }
+ else{
+
+ for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
+ AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
+ if (!aodTrack) {
+ AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
+ continue;
+ }
+
+ if(!aodTrack->IsPhysicalPrimary()) continue;
+
+ vCharge = aodTrack->Charge();
+ vEta = aodTrack->Eta();
+ vY = aodTrack->Y();
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();
+ vPt = aodTrack->Pt();
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+ // Remove neutral tracks
+ if( vCharge == 0 ) continue;
+
+ //Exclude resonances
+ if(fExcludeResonancesInMC) {
+
+ Bool_t kExcludeParticle = kFALSE;
+ Int_t gMotherIndex = aodTrack->GetMother();
+ if(gMotherIndex != -1) {
+ AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+ if(motherTrack) {
+ Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+ //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
+ //if(pdgCodeOfMother == 113) {
+ if(pdgCodeOfMother == 113 // rho0
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
+ // || pdgCodeOfMother == 221 // eta
+ // || pdgCodeOfMother == 331 // eta'
+ // || pdgCodeOfMother == 223 // omega
+ // || pdgCodeOfMother == 333 // phi
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
+ || pdgCodeOfMother == 111 // pi0 Dalitz
+ || pdgCodeOfMother == 22 // photon
+ ) {
+ kExcludeParticle = kTRUE;
+ }
+ }
+ }
+
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
+ if(kExcludeParticle) continue;
+ }
+
+ //Exclude electrons with PDG
+ if(fExcludeElectronsInMC) {
+
+ if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
+
+ }
+
+ // fill QA histograms
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ // add the track to the TObjArray
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+ }//aodTracks
+ }//MC event
+ }//MCAOD
+ //==============================================================================================================
+
+ //==============================================================================================================
+ else if(gAnalysisLevel == "MCAODrec") {
+
+ /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+ if (!fAOD) {
+ printf("ERROR: fAOD not available\n");
+ return;
+ }*/
+
+ fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!fArrayMC) {
+ AliError("No array of MC particles found !!!");
+ }
+
+ AliMCEvent* mcEvent = MCEvent();
+ if (!mcEvent) {
+ AliError("ERROR: Could not retrieve MC event");
+ return tracksAccepted;
+ }
+
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
+ if (!aodTrack) {
+ AliError(Form("Could not receive track %d", iTracks));
+ continue;
+ }
+
+ for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
+ fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
+ }
+ if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
+
+ vCharge = aodTrack->Charge();
+ vEta = aodTrack->Eta();
+ vY = aodTrack->Y();
+ vPhi = aodTrack->Phi();// * TMath::RadToDeg();
+ vPt = aodTrack->Pt();
+
+ //===========================use MC information for Kinematics===============================//
+ if(fUseMCforKinematics){
+
+ Int_t label = TMath::Abs(aodTrack->GetLabel());
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+
+ if(AODmcTrack){
+ vCharge = AODmcTrack->Charge();
+ vEta = AODmcTrack->Eta();
+ vY = AODmcTrack->Y();
+ vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
+ vPt = AODmcTrack->Pt();
+ }
+ else{
+ AliDebug(1, "no MC particle for this track");
+ continue;
+ }
+ }
+ //===========================end of use MC information for Kinematics========================//
+
+
+ //===========================PID (so far only for electron rejection)===============================//
+ if(fElectronRejection) {
+
+ // get the electron nsigma
+ Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
+
+ //Fill QA before the PID
+ fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
+ fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
+ //end of QA-before pid
+
+ // check only for given momentum range
+ if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
+
+ //look only at electron nsigma
+ if(!fElectronOnlyRejection){
+
+ //Make the decision based on the n-sigma of electrons
+ if(nSigma < fElectronRejectionNSigma) continue;
+ }
+ else{
+
+ Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
+ Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
+ Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
+
+ //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
+ if(nSigma < fElectronRejectionNSigma
+ && nSigmaPions > fElectronRejectionNSigma
+ && nSigmaKaons > fElectronRejectionNSigma
+ && nSigmaProtons > fElectronRejectionNSigma ) continue;
+ }
+ }
+
+ //Fill QA after the PID
+ fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
+ fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
+
+ }
+ //===========================end of PID (so far only for electron rejection)===============================//
+
+ Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
+ Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
+
+ // Kinematics cuts from ESD track cuts
+ if( vPt < fPtMin || vPt > fPtMax) continue;
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+
+ // Extra DCA cuts (for systematic studies [!= -1])
+ if( fDCAxyCut != -1 && fDCAzCut != -1){
+ if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
+ continue; // 2D cut
+ }
+ }
+
+ // Extra TPC cuts (for systematic studies [!= -1])
+ if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
+ continue;
+ }
+ if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
+ continue;
+ }
+
+ // Extra cut on shared clusters
+ if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
+ continue;
+ }
+
+ //Exclude resonances
+ if(fExcludeResonancesInMC) {
+
+ Bool_t kExcludeParticle = kFALSE;
+
+ Int_t label = TMath::Abs(aodTrack->GetLabel());
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+
+ if (AODmcTrack){
+ //if (AODmcTrack->IsPhysicalPrimary()){
+
+ Int_t gMotherIndex = AODmcTrack->GetMother();
+ if(gMotherIndex != -1) {
+ AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+ if(motherTrack) {
+ Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+ if(pdgCodeOfMother == 113 // rho0
+ || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
+ // || pdgCodeOfMother == 221 // eta
+ // || pdgCodeOfMother == 331 // eta'
+ // || pdgCodeOfMother == 223 // omega
+ // || pdgCodeOfMother == 333 // phi
+ || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
+ // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
+ // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
+ || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
+ || pdgCodeOfMother == 111 // pi0 Dalitz
+ || pdgCodeOfMother == 22 // photon
+ ) {
+ kExcludeParticle = kTRUE;
+ }
+ }
+ }
+ }
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
+ if(kExcludeParticle) continue;
+ }
+
+ //Exclude electrons with PDG
+ if(fExcludeElectronsInMC) {
+
+ Int_t label = TMath::Abs(aodTrack->GetLabel());
+ AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+
+ if (AODmcTrack){
+ if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
+ }
+ }
+
+ // fill QA histograms
+ fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
+ fHistDCA->Fill(dcaZ,dcaXY);
+ fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ // add the track to the TObjArray
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+ }//track loop
+ }//MCAODrec
+ //==============================================================================================================
+
+ else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
+
+ AliESDtrack *trackTPC = NULL;
+ AliMCParticle *mcTrack = NULL;
+
+ AliMCEvent* mcEvent = NULL;
+
+ // for MC ESDs use also MC information
+ if(gAnalysisLevel == "MCESD"){
+ mcEvent = MCEvent();
+ if (!mcEvent) {
+ AliError("mcEvent not available");
+ return tracksAccepted;
+ }
+ }
+
+ // Loop over tracks in event
+ for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+ AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
+ if (!track) {
+ AliError(Form("Could not receive track %d", iTracks));
+ continue;
+ }
+
+ // for MC ESDs use also MC information --> MC track not used further???
+ if(gAnalysisLevel == "MCESD"){
+ Int_t label = TMath::Abs(track->GetLabel());
+ if(label > mcEvent->GetNumberOfTracks()) continue;
+ if(label > mcEvent->GetNumberOfPrimaries()) continue;
+
+ mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
+ if(!mcTrack) continue;
+ }
+
+ // take only TPC only tracks
+ trackTPC = new AliESDtrack();
+ if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
+
+ //ESD track cuts
+ if(fESDtrackCuts)
+ if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
+
+ // fill QA histograms
+ Float_t b[2];
+ Float_t bCov[3];
+ trackTPC->GetImpactParameters(b,bCov);
+ if (bCov[0]<=0 || bCov[2]<=0) {
+ AliDebug(1, "Estimated b resolution lower or equal zero!");
+ bCov[0]=0; bCov[2]=0;
+ }
+
+ Int_t nClustersTPC = -1;
+ nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
+ //nClustersTPC = track->GetTPCclusters(0); // global track
+ Float_t chi2PerClusterTPC = -1;
+ if (nClustersTPC!=0) {
+ chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
+ //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
+ }
+
+ //===========================PID===============================//
+ if(fUsePID) {
+ Double_t prob[AliPID::kSPECIES]={0.};
+ Double_t probTPC[AliPID::kSPECIES]={0.};
+ Double_t probTOF[AliPID::kSPECIES]={0.};
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+
+ Double_t nSigma = 0.;
+ UInt_t detUsedTPC = 0;
+ UInt_t detUsedTOF = 0;
+ UInt_t detUsedTPCTOF = 0;
+
+ //Decide what detector configuration we want to use
+ switch(fPidDetectorConfig) {
+ case kTPCpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
+ detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPC[iSpecies];
+ break;
+ case kTOFpid:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+ nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
+ detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTOF[iSpecies];
+ break;
+ case kTPCTOF:
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
+ detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
+ for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
+ prob[iSpecies] = probTPCTOF[iSpecies];
+ break;
+ default:
+ break;
+ }//end switch: define detector mask
+
+ //Filling the PID QA
+ Double_t tofTime = -999., length = 999., tof = -999.;
+ Double_t c = TMath::C()*1.E-9;// m/ns
+ Double_t beta = -999.;
+ Double_t nSigmaTOFForParticleOfInterest = -999.;
+ if ( (track->IsOn(AliESDtrack::kTOFin)) &&
+ (track->IsOn(AliESDtrack::kTIME)) ) {
+ tofTime = track->GetTOFsignal();//in ps
+ length = track->GetIntegratedLength();
+ tof = tofTime*1E-3; // ns
+
+ if (tof <= 0) {
+ //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
+ continue;
+ }
+ if (length <= 0){
+ //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
+ continue;
+ }
+
+ length = length*0.01; // in meters
+ tof = tof*c;
+ beta = length/tof;
+
+ nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
+ fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
+ fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
+ }//TOF signal
+
+
+ Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
+ fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
+ fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
+ fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
+ //end of QA-before pid
+
+ if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
+ //Make the decision based on the n-sigma
+ if(fUsePIDnSigma) {
+ if(nSigma > fPIDNSigma) continue;}
+
+ //Make the decision based on the bayesian
+ else if(fUsePIDPropabilities) {
+ if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
+ if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
+ }
+
+ //Fill QA after the PID
+ fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
+ fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
+ fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
+
+ fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
+ fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
+ fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
+ fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
+ }
+ }
+ //===========================PID===============================//
+ vCharge = trackTPC->Charge();
+ vY = trackTPC->Y();
+ vEta = trackTPC->Eta();
+ vPhi = trackTPC->Phi();// * TMath::RadToDeg();
+ vPt = trackTPC->Pt();
+ fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
+ fHistDCA->Fill(b[1],b[0]);
+ fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ // add the track to the TObjArray
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+
+ delete trackTPC;
+ }//track loop
+ }// ESD analysis
+
+ else if(gAnalysisLevel == "MC"){
+ if(!event) {
+ AliError("mcEvent not available");
+ return 0x0;
+ }
+
+ AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
+ if(gMCEvent) {
+ // Loop over tracks in event
+ for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
+ AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
+ if (!track) {
+ AliError(Form("Could not receive particle %d", iTracks));
+ continue;
+ }
+
+ //exclude non stable particles
+ if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
+
+ vCharge = track->Charge();
+ vEta = track->Eta();
+ vPt = track->Pt();
+ vY = track->Y();
+
+ if( vPt < fPtMin || vPt > fPtMax)
+ continue;
+ if (!fUsePID) {
+ if( vEta < fEtaMin || vEta > fEtaMax) continue;
+ }
+ else if (fUsePID){
+ if( vY < fEtaMin || vY > fEtaMax) continue;
+ }
+
+ // Remove neutral tracks
+ if( vCharge == 0 ) continue;
+
+ //analyze one set of particles
+ if(fUseMCPdgCode) {
+ TParticle *particle = track->Particle();
+ if(!particle) continue;
+
+ Int_t gPdgCode = particle->GetPdgCode();
+ if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
+ continue;
+ }
+
+ //Use the acceptance parameterization
+ if(fAcceptanceParameterization) {
+ Double_t gRandomNumber = gRandom->Rndm();
+ if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
+ continue;
+ }
+
+ //Exclude resonances
+ if(fExcludeResonancesInMC) {
+ TParticle *particle = track->Particle();
+ if(!particle) continue;
+
+ Bool_t kExcludeParticle = kFALSE;
+ Int_t gMotherIndex = particle->GetFirstMother();
+
+ if(gMotherIndex != -1) {
+ // AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
+ // if(motherTrack) {
+ // TParticle *motherParticle = motherTrack->Particle();
+ // if(motherParticle) {
+ // Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
+ // //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
+ // if(pdgCodeOfMother == 113 // rho0
+ // || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
+ // // || pdgCodeOfMother == 221 // eta
+ // // || pdgCodeOfMother == 331 // eta'
+ // // || pdgCodeOfMother == 223 // omega
+ // // || pdgCodeOfMother == 333 // phi
+ // || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
+ // // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
+ // // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
+ // || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
+ // || pdgCodeOfMother == 111 // pi0 Dalitz
+ // ) {
+ kExcludeParticle = kTRUE;
+ // }
+ // }
+ //}
+ }
+
+ //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
+ if(kExcludeParticle) continue;
+
+ }
+
+ //Exclude electrons with PDG
+ if(fExcludeElectronsInMC) {
+
+ TParticle *particle = track->Particle();
+
+ if (particle){
+ if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
+ }
+ }
+
+ vPhi = track->Phi();
+ //Printf("phi (before): %lf",vPhi);
+
+ fHistPt->Fill(vPt,gCentrality);
+ fHistEta->Fill(vEta,gCentrality);
+ fHistPhi->Fill(vPhi,gCentrality);
+ if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
+ else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
+ //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
+ fHistRapidity->Fill(vY,gCentrality);
+ //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
+ //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
+ if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
+ else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
+
+ //Flow after burner
+ if(fUseFlowAfterBurner) {
+ Double_t precisionPhi = 0.001;
+ Int_t maxNumberOfIterations = 100;
+
+ Double_t phi0 = vPhi;
+ Double_t gV2 = fDifferentialV2->Eval(vPt);
+
+ for (Int_t j = 0; j < maxNumberOfIterations; j++) {
+ Double_t phiprev = vPhi;
+ Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
+ Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
+ vPhi -= fl/fp;
+ if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
+ }
+ //Printf("phi (after): %lf\n",vPhi);
+ Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
+ if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
+ fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
+
+ Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
+ if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
+ fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
+
+ }
+
+ //vPhi *= TMath::RadToDeg();
+
+ //=======================================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+
+ tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
+ } //track loop
+ }//MC event object
+ }//MC
+
+ return tracksAccepted;
+}
+
+//________________________________________________________________________
+TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
+ // Clones TObjArray and returns it with tracks after shuffling the charges
+
+ TObjArray* tracksShuffled = new TObjArray;
+ tracksShuffled->SetOwner(kTRUE);
+
+ vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
+
+ for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
+ {
+ AliVParticle* track = (AliVParticle*) tracks->At(i);
+ chargeVector->push_back(track->Charge());
+ }
+
+ random_shuffle(chargeVector->begin(), chargeVector->end());
+
+ for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
+ AliVParticle* track = (AliVParticle*) tracks->At(i);
+ //==============================correction
+ Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
+ //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
+ tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
+ }
+
+ delete chargeVector;
+
+ return tracksShuffled;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
+ const char* lhcPeriod) {
+ //Function to setup the VZERO gain equalization
+ //============Get the equilization map============//
+ TFile *calibrationFile = TFile::Open(filename);
+ if((!calibrationFile)||(!calibrationFile->IsOpen())) {
+ Printf("No calibration file found!!!");
+ return;
+ }
+
+ TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
+ if(!list) {
+ Printf("Calibration TList not found!!!");
+ return;
+ }
+
+ fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
+ if(!fHistVZEROAGainEqualizationMap) {
+ Printf("VZERO-A calibration object not found!!!");
+ return;
+ }
+ fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
+ if(!fHistVZEROCGainEqualizationMap) {
+ Printf("VZERO-C calibration object not found!!!");
+ return;
+ }
+
+ fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
+ if(!fHistVZEROChannelGainEqualizationMap) {
+ Printf("VZERO channel calibration object not found!!!");
+ return;
+ }
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
+ Int_t channel) {
+ //
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;
+
+ for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
+ Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
+ if(gRunNumber == run)
+ return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
+ }
+
+ return 1.0;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
+ const char* side) {
+ //
+ if(!fHistVZEROAGainEqualizationMap) return 1.0;
+
+ TString gVZEROSide = side;
+ for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
+ Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
+ //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
+ if(gRunNumber == run) {
+ if(gVZEROSide == "A")
+ return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
+ else if(gVZEROSide == "C")
+ return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
+ }
+ }
+
+ return 1.0;
+}
+
+//____________________________________________________________________
+Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
+{
+ // copied from AliAnalysisTaskPhiCorrelations
+ //
+ // rejects "randomly" events such that the centrality gets flat
+ // uses fCentralityWeights histogram
+
+ // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
+
+ Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
+ Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
+
+ Bool_t result = kFALSE;
+ if (centralityDigits < weight)
+ result = kTRUE;
+
+ AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
+
+ return result;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::FinishTaskOutput(){
+ //Printf("END BF");
+
+ if (!fBalance) {
+ AliError("fBalance not available");
+ return;
+ }
+ if(fRunShuffling) {
+ if (!fShuffledBalance) {
+ AliError("fShuffledBalance not available");
+ return;
+ }
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
+ // Draw result to the screen
+ // Called once at the end of the query
+
+ // not implemented ...
+
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff