4 #include "TLorentzVector.h"
5 #include "TGraphErrors.h"
16 #include "AliAnalysisTaskSE.h"
17 #include "AliAnalysisManager.h"
19 #include "AliESDVertex.h"
20 #include "AliESDEvent.h"
21 #include "AliESDInputHandler.h"
22 #include "AliAODEvent.h"
23 #include "AliAODTrack.h"
24 #include "AliAODInputHandler.h"
25 #include "AliAODMCParticle.h"
26 #include "AliCollisionGeometry.h"
27 #include "AliGenEventHeader.h"
28 #include "AliMCEventHandler.h"
29 #include "AliMCEvent.h"
31 #include "AliESDtrackCuts.h"
32 #include "AliEventplane.h"
35 #include "AliAnalysisUtils.h"
37 #include "AliEventPoolManager.h"
40 #include "AliPIDResponse.h"
41 #include "AliPIDCombined.h"
43 #include "AliAnalysisTaskBFPsi.h"
44 #include "AliBalancePsi.h"
45 #include "AliAnalysisTaskTriggeredBF.h"
48 // Analysis task for the BF vs Psi code
49 // Authors: Panos.Christakoglou@nikhef.nl
54 ClassImp(AliAnalysisTaskBFPsi)
56 //________________________________________________________________________
57 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
58 : AliAnalysisTaskSE(name),
62 fRunShuffling(kFALSE),
65 fRunMixingEventPlane(kFALSE),
76 fHistCentStatsUsed(0),
82 fHistTPCvsVZEROMultiplicity(0),
100 fHistdEdxVsPTPCbeforePID(NULL),
101 fHistBetavsPTOFbeforePID(NULL),
102 fHistProbTPCvsPtbeforePID(NULL),
103 fHistProbTOFvsPtbeforePID(NULL),
104 fHistProbTPCTOFvsPtbeforePID(NULL),
105 fHistNSigmaTPCvsPtbeforePID(NULL),
106 fHistNSigmaTOFvsPtbeforePID(NULL),
107 fHistBetaVsdEdXbeforePID(NULL),
108 fHistNSigmaTPCTOFvsPtbeforePID(NULL),
109 fHistNSigmaTPCTOFPbefPID(NULL),
110 fHistdEdxVsPTPCafterPID(NULL),
111 fHistBetavsPTOFafterPID(NULL),
112 fHistProbTPCvsPtafterPID(NULL),
113 fHistProbTOFvsPtafterPID(NULL),
114 fHistProbTPCTOFvsPtafterPID(NULL),
115 fHistNSigmaTPCvsPtafterPID(NULL),
116 fHistNSigmaTOFvsPtafterPID(NULL),
117 fHistBetaVsdEdXafterPID(NULL),
118 fHistNSigmaTPCTOFvsPtafterPID(NULL),
119 fHistNSigmaTPCTOFPafterPID(NULL),
120 fHistdEdxVsPTPCbeforePIDelectron(NULL),
121 fHistNSigmaTPCvsPtbeforePIDelectron(NULL),
122 fHistdEdxVsPTPCafterPIDelectron(NULL),
123 fHistNSigmaTPCvsPtafterPIDelectron(NULL),
124 fCentralityArrayBinsForCorrections(kCENTRALITY),
125 fCentralityWeights(0x0),
128 fParticleOfInterest(kPion),
129 fPidDetectorConfig(kTPCTOF),
131 fUsePIDnSigma(kTRUE),
132 fUsePIDPropabilities(kFALSE),
134 fMinAcceptedPIDProbability(0.8),
135 fElectronRejection(kFALSE),
136 fElectronOnlyRejection(kFALSE),
137 fElectronRejectionNSigma(-1.),
138 fElectronRejectionMinPt(0.),
139 fElectronRejectionMaxPt(1000.),
141 fCentralityEstimator("V0M"),
142 fUseCentrality(kFALSE),
143 fCentralityPercentileMin(0.),
144 fCentralityPercentileMax(5.),
145 fImpactParameterMin(0.),
146 fImpactParameterMax(20.),
147 fMultiplicityEstimator("V0A"),
148 fUseMultiplicity(kFALSE),
149 fNumberOfAcceptedTracksMin(0),
150 fNumberOfAcceptedTracksMax(10000),
151 fHistNumberOfAcceptedTracks(0),
152 fHistMultiplicity(0),
153 fUseOfflineTrigger(kFALSE),
154 fCheckFirstEventInChunk(kFALSE),
155 fCheckPileUp(kFALSE),
156 fCheckPrimaryFlagAOD(kFALSE),
157 fUseMCforKinematics(kFALSE),
161 fnAODtrackCutBit(128),
171 fNClustersTPCCut(-1),
173 fAcceptanceParameterization(0),
175 fUseFlowAfterBurner(kFALSE),
176 fExcludeResonancesInMC(kFALSE),
177 fExcludeElectronsInMC(kFALSE),
178 fUseMCPdgCode(kFALSE),
179 fPDGCodeToBeAnalyzed(-1),
180 fEventClass("EventPlane"),
182 fHistVZEROAGainEqualizationMap(0),
183 fHistVZEROCGainEqualizationMap(0),
184 fHistVZEROChannelGainEqualizationMap(0) {
186 // Define input and output slots here
187 // Input slot #0 works with a TChain
189 //======================================================correction
190 for (Int_t i=0; i<kCENTRALITY; i++){
191 fHistCorrectionPlus[i] = NULL;
192 fHistCorrectionMinus[i] = NULL;
193 fCentralityArrayForCorrections[i] = -1.;
195 //=====================================================correction
197 DefineInput(0, TChain::Class());
198 // Output slot #0 writes into a TH1 container
199 DefineOutput(1, TList::Class());
200 DefineOutput(2, TList::Class());
201 DefineOutput(3, TList::Class());
202 DefineOutput(4, TList::Class());
203 DefineOutput(5, TList::Class());
206 //________________________________________________________________________
207 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
210 // delete fShuffledBalance;
215 // delete fHistEventStats;
216 // delete fHistTrackStats;
227 // delete fHistEtaPhiPos;
228 // delete fHistEtaPhiNeg;
232 //________________________________________________________________________
233 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
237 // global switch disabling the reference
238 // (to avoid "Replacing existing TH1" if several wagons are created in train)
239 Bool_t oldStatus = TH1::AddDirectoryStatus();
240 TH1::AddDirectory(kFALSE);
243 fBalance = new AliBalancePsi();
244 fBalance->SetAnalysisLevel("ESD");
245 fBalance->SetEventClass(fEventClass);
246 //fBalance->SetNumberOfBins(-1,16);
247 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
250 if(!fShuffledBalance) {
251 fShuffledBalance = new AliBalancePsi();
252 fShuffledBalance->SetAnalysisLevel("ESD");
253 fShuffledBalance->SetEventClass(fEventClass);
254 //fShuffledBalance->SetNumberOfBins(-1,16);
255 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
260 fMixedBalance = new AliBalancePsi();
261 fMixedBalance->SetAnalysisLevel("ESD");
262 fMixedBalance->SetEventClass(fEventClass);
268 fList->SetName("listQA");
271 //Balance Function list
272 fListBF = new TList();
273 fListBF->SetName("listBF");
277 fListBFS = new TList();
278 fListBFS->SetName("listBFShuffled");
279 fListBFS->SetOwner();
283 fListBFM = new TList();
284 fListBFM->SetName("listTriggeredBFMixed");
285 fListBFM->SetOwner();
289 if(fUsePID || fElectronRejection) {
290 fHistListPIDQA = new TList();
291 fHistListPIDQA->SetName("listQAPID");
292 fHistListPIDQA->SetOwner();
296 TString gCutName[7] = {"Total","Offline trigger",
297 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
298 fHistEventStats = new TH2F("fHistEventStats",
299 "Event statistics;;Centrality percentile;N_{events}",
300 7,0.5,7.5,220,-5,105);
301 for(Int_t i = 1; i <= 7; i++)
302 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
303 fList->Add(fHistEventStats);
305 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
306 fHistCentStats = new TH2F("fHistCentStats",
307 "Centrality statistics;;Cent percentile",
308 13,-0.5,12.5,220,-5,105);
309 for(Int_t i = 1; i <= 13; i++){
310 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
311 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
313 fList->Add(fHistCentStats);
315 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
316 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
317 fList->Add(fHistCentStatsUsed);
319 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
320 fList->Add(fHistTriggerStats);
322 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
323 fList->Add(fHistTrackStats);
325 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
326 fList->Add(fHistNumberOfAcceptedTracks);
328 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
329 fList->Add(fHistMultiplicity);
331 // Vertex distributions
332 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
334 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
336 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
339 //TPC vs VZERO multiplicity
340 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
341 if(fMultiplicityEstimator == "V0A")
342 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
343 else if(fMultiplicityEstimator == "V0C")
344 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
346 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
347 fList->Add(fHistTPCvsVZEROMultiplicity);
349 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
350 fList->Add(fHistVZEROSignal);
353 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
354 fList->Add(fHistEventPlane);
357 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
358 fList->Add(fHistClus);
359 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
360 fList->Add(fHistChi2);
361 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
362 fList->Add(fHistDCA);
363 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
365 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
366 fList->Add(fHistEta);
367 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
368 fList->Add(fHistRapidity);
369 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
370 fList->Add(fHistPhi);
371 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);
372 fList->Add(fHistEtaPhiPos);
373 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);
374 fList->Add(fHistEtaPhiNeg);
375 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
376 fList->Add(fHistPhiBefore);
377 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
378 fList->Add(fHistPhiAfter);
379 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
380 fList->Add(fHistPhiPos);
381 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
382 fList->Add(fHistPhiNeg);
383 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
384 fList->Add(fHistV0M);
385 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
386 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
387 for(Int_t i = 1; i <= 6; i++)
388 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
389 fList->Add(fHistRefTracks);
391 // Balance function histograms
392 // Initialize histograms if not done yet (including the custom binning)
393 if(!fBalance->GetHistNp()){
394 AliInfo("Histograms not yet initialized! --> Will be done now");
395 fBalance->SetCustomBinning(fCustomBinning);
396 fBalance->InitHistograms();
400 if(!fShuffledBalance->GetHistNp()) {
401 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
402 fShuffledBalance->SetCustomBinning(fCustomBinning);
403 fShuffledBalance->InitHistograms();
408 if(!fMixedBalance->GetHistNp()) {
409 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
410 fMixedBalance->SetCustomBinning(fCustomBinning);
411 fMixedBalance->InitHistograms();
415 // QA histograms for different cuts
416 fList->Add(fBalance->GetQAHistHBTbefore());
417 fList->Add(fBalance->GetQAHistHBTafter());
418 fList->Add(fBalance->GetQAHistConversionbefore());
419 fList->Add(fBalance->GetQAHistConversionafter());
420 fList->Add(fBalance->GetQAHistPsiMinusPhi());
421 fList->Add(fBalance->GetQAHistResonancesBefore());
422 fList->Add(fBalance->GetQAHistResonancesRho());
423 fList->Add(fBalance->GetQAHistResonancesK0());
424 fList->Add(fBalance->GetQAHistResonancesLambda());
425 fList->Add(fBalance->GetQAHistQbefore());
426 fList->Add(fBalance->GetQAHistQafter());
428 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
429 fListBF->Add(fBalance->GetHistNp());
430 fListBF->Add(fBalance->GetHistNn());
431 fListBF->Add(fBalance->GetHistNpn());
432 fListBF->Add(fBalance->GetHistNnn());
433 fListBF->Add(fBalance->GetHistNpp());
434 fListBF->Add(fBalance->GetHistNnp());
437 fListBFS->Add(fShuffledBalance->GetHistNp());
438 fListBFS->Add(fShuffledBalance->GetHistNn());
439 fListBFS->Add(fShuffledBalance->GetHistNpn());
440 fListBFS->Add(fShuffledBalance->GetHistNnn());
441 fListBFS->Add(fShuffledBalance->GetHistNpp());
442 fListBFS->Add(fShuffledBalance->GetHistNnp());
446 fListBFM->Add(fMixedBalance->GetHistNp());
447 fListBFM->Add(fMixedBalance->GetHistNn());
448 fListBFM->Add(fMixedBalance->GetHistNpn());
449 fListBFM->Add(fMixedBalance->GetHistNnn());
450 fListBFM->Add(fMixedBalance->GetHistNpp());
451 fListBFM->Add(fMixedBalance->GetHistNnp());
458 Int_t trackDepth = fMixingTracks;
459 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
462 Double_t* centbins = NULL;
463 Int_t nCentralityBins;
464 if(fBalance->IsUseVertexBinning()){
465 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
468 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
472 Double_t* multbins = NULL;
473 Int_t nMultiplicityBins;
474 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
477 Double_t* vtxbins = NULL;
479 if(fBalance->IsUseVertexBinning()){
480 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
483 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
486 // Event plane angle (Psi) bins
487 Double_t* psibins = NULL;
489 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
492 // run the event mixing also in bins of event plane (statistics!)
493 if(fRunMixingEventPlane){
494 if(fEventClass=="Multiplicity"){
495 if(multbins && vtxbins && psibins){
496 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
500 if(centbins && vtxbins && psibins){
501 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
506 if(fEventClass=="Multiplicity"){
507 if(multbins && vtxbins){
508 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
512 if(centbins && vtxbins){
513 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
518 if(centbins) delete [] centbins;
519 if(multbins) delete [] multbins;
520 if(vtxbins) delete [] vtxbins;
521 if(psibins) delete [] psibins;
523 // set minimum values for track depth, fraction, and number of events
524 fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
526 // check pool manager
528 AliError("Event Mixing required, but Pool Manager not initialized...");
533 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
535 //====================PID========================//
537 fPIDCombined = new AliPIDCombined();
538 fPIDCombined->SetDefaultTPCPriors();
540 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
541 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
543 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
544 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
546 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
547 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
549 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
550 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
552 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
553 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
555 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, -25, 25);
556 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
558 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, -25, 25);
559 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
561 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
562 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID);
564 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, -25, 25);
565 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID);
567 //+++++++++++++++++//
569 const Int_t pBins = 36;
570 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};
572 const Int_t nSigmaBins = 250;
573 Double_t nArrayS[nSigmaBins+1];
574 for (Int_t i = 0; i <= nSigmaBins; i++){
575 nArrayS[i]=i-125; //i+1
576 //Printf("nS: %lf - i: %d", nSigmaArray[i], i);
579 fHistNSigmaTPCTOFPbefPID = new TH3D ("fHistNSigmaTPCTOFPbefPID","fHistNSigmaTPCTOFPbefPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP);
580 fHistListPIDQA->Add(fHistNSigmaTPCTOFPbefPID);
583 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
584 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
586 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
587 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
589 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
590 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
592 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
593 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
595 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
596 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
598 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, -25, 25);
599 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
601 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, -25, 25);
602 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
604 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
605 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID);
607 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, -25, 25);
608 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID);
610 fHistNSigmaTPCTOFPafterPID = new TH3D ("fHistNSigmaTPCTOFPafterPID","fHistNSigmaTPCTOFPafterPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP);
611 fHistListPIDQA->Add(fHistNSigmaTPCTOFPafterPID); //++++++++++++++
614 // for electron rejection only TPC nsigma histograms
615 if(fElectronRejection) {
617 fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
618 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
620 fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
621 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
623 fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
624 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
626 fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
627 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
629 //====================PID========================//
633 PostData(2, fListBF);
634 if(fRunShuffling) PostData(3, fListBFS);
635 if(fRunMixing) PostData(4, fListBFM);
636 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
638 AliInfo("Finished setting up the Output");
640 TH1::AddDirectory(oldStatus);
645 //________________________________________________________________________
646 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
647 Int_t nCentralityBins,
648 Double_t *centralityArrayForCorrections) {
649 //Open files that will be used for correction
650 fCentralityArrayBinsForCorrections = nCentralityBins;
651 for (Int_t i=0; i<nCentralityBins; i++)
652 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
654 // No file specified -> run without corrections
655 if(!filename.Contains(".root")) {
656 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
660 //Open the input file
661 TFile *f = TFile::Open(filename);
663 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
667 //TString listEffName = "";
668 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
669 //Printf("iCent %d:",iCent);
670 TString histoName = "fHistCorrectionPlus";
671 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
672 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
673 if(!fHistCorrectionPlus[iCent]) {
674 AliError(Form("fHist %s not found!!!",histoName.Data()));
678 histoName = "fHistCorrectionMinus";
679 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
680 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
681 if(!fHistCorrectionMinus[iCent]) {
682 AliError(Form("fHist %s not found!!!",histoName.Data()));
685 }//loop over centralities: ONLY the PbPb case is covered
688 //________________________________________________________________________
689 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
691 // Called for each event
693 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
694 Int_t gNumberOfAcceptedTracks = 0;
695 Double_t lMultiplicityVar = -999.; //-1
696 Double_t gReactionPlane = -1.;
699 // get the event (for generator level: MCEvent())
700 AliVEvent* eventMain = NULL;
701 if(gAnalysisLevel == "MC") {
702 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
705 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
706 // for HBT like cuts need magnetic field sign
707 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
710 AliError("eventMain not available");
714 // PID Response task active?
715 if(fUsePID || fElectronRejection) {
716 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
717 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
720 // check event cuts and fill event histograms
721 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
724 // get the reaction plane
725 if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") {
726 gReactionPlane = GetEventPlane(eventMain);
727 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
728 if(gReactionPlane < 0){
733 // get the accepted tracks in main event
734 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
735 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
737 //multiplicity cut (used in pp)
738 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
740 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
741 TObjArray* tracksShuffled = NULL;
743 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
749 // 1. First get an event pool corresponding in mult (cent) and
750 // zvertex to the current event. Once initialized, the pool
751 // should contain nMix (reduced) events. This routine does not
752 // pre-scan the chain. The first several events of every chain
753 // will be skipped until the needed pools are filled to the
754 // specified depth. If the pool categories are not too rare, this
755 // should not be a problem. If they are rare, you could lose`
758 // 2. Collect the whole pool's content of tracks into one TObjArray
759 // (bgTracks), which is effectively a single background super-event.
761 // 3. The reduced and bgTracks arrays must both be passed into
762 // FillCorrelations(). Also nMix should be passed in, so a weight
763 // of 1./nMix can be applied.
765 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
768 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
774 if (pool->IsReady()){
777 Int_t nMix = pool->GetCurrentNEvents();
778 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
780 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
781 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
782 //if (pool->IsReady())
783 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
785 // Fill mixed-event histos here
786 for (Int_t jMix=0; jMix<nMix; jMix++)
788 TObjArray* tracksMixed = pool->GetEvent(jMix);
789 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
793 // Update the Event pool
794 pool->UpdatePool(tracksMain);
800 // calculate balance function
801 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
803 // calculate shuffled balance function
804 if(fRunShuffling && tracksShuffled != NULL) {
805 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
809 //________________________________________________________________________
810 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
811 // Checks the Event cuts
812 // Fills Event statistics histograms
814 Bool_t isSelectedMain = kTRUE;
815 Float_t gRefMultiplicity = -1.;
816 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
818 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
820 fHistEventStats->Fill(1,gRefMultiplicity); //all events
822 // check first event in chunk (is not needed for new reconstructions)
823 if(fCheckFirstEventInChunk){
825 if(ut.IsFirstEventInChunk(event))
827 fHistEventStats->Fill(6,gRefMultiplicity);
829 // check for pile-up event
832 ut.SetUseMVPlpSelection(kTRUE);
833 ut.SetUseOutOfBunchPileUp(kTRUE);
834 if(ut.IsPileUpEvent(event))
836 fHistEventStats->Fill(7,gRefMultiplicity);
839 // Event trigger bits
840 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
841 if(fUseOfflineTrigger)
842 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
845 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
848 if(gAnalysisLevel == "MC") {
850 AliError("mcEvent not available");
855 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
857 TArrayF gVertexArray;
858 header->PrimaryVertex(gVertexArray);
859 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
860 //gVertexArray.At(0),
861 //gVertexArray.At(1),
862 //gVertexArray.At(2));
863 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
864 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
865 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
866 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
867 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
869 // get the reference multiplicty or centrality
870 gRefMultiplicity = GetRefMultiOrCentrality(event);
872 fHistVx->Fill(gVertexArray.At(0));
873 fHistVy->Fill(gVertexArray.At(1));
874 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
876 // take only events inside centrality class
878 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
879 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
880 return gRefMultiplicity;
883 // take events only within the same multiplicity class
884 else if(fUseMultiplicity) {
885 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
886 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
887 return gRefMultiplicity;
889 }//multiplicity range
897 // Event Vertex AOD, ESD, ESDMC
899 const AliVVertex *vertex = event->GetPrimaryVertex();
903 vertex->GetCovarianceMatrix(fCov);
904 if(vertex->GetNContributors() > 0) {
906 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
907 if(TMath::Abs(vertex->GetX()) < fVxMax) {
908 if(TMath::Abs(vertex->GetY()) < fVyMax) {
909 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
910 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
912 // get the reference multiplicty or centrality
913 gRefMultiplicity = GetRefMultiOrCentrality(event);
915 fHistVx->Fill(vertex->GetX());
916 fHistVy->Fill(vertex->GetY());
917 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
919 // take only events inside centrality class
921 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
923 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
924 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
925 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
929 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
930 return gRefMultiplicity;
933 // take events only within the same multiplicity class
934 else if(fUseMultiplicity) {
936 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
937 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
939 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
940 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
941 return gRefMultiplicity;
943 }//multiplicity range
947 }//proper vertex resolution
948 }//proper number of contributors
949 }//vertex object valid
953 // in all other cases return -1 (event not accepted)
958 //________________________________________________________________________
959 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
960 // Checks the Event cuts
961 // Fills Event statistics histograms
963 Float_t gCentrality = -1.;
964 Double_t gMultiplicity = -1.;
965 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
968 // calculate centrality always (not only in centrality mode)
969 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
970 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
972 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
974 // QA for centrality estimators
975 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
976 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
977 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
978 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
979 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
980 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
981 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
982 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
983 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
984 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
985 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
986 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
987 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
989 // Centrality estimator USED ++++++++++++++++++++++++++++++
990 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
992 // centrality QA (V0M)
993 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
995 // centrality QA (reference tracks)
996 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
997 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
998 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
999 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
1000 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
1001 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
1002 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
1003 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
1004 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
1009 // calculate centrality always (not only in centrality mode)
1010 else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround
1012 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1014 gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0;
1017 fHistCentStatsUsed->Fill(0.,gCentrality);
1022 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
1023 AliCentrality *centrality = event->GetCentrality();
1024 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
1026 // QA for centrality estimators
1027 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
1028 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
1029 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
1030 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
1031 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
1032 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
1033 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
1034 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
1035 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
1036 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
1037 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
1038 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
1039 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
1041 // Centrality estimator USED ++++++++++++++++++++++++++++++
1042 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
1044 // centrality QA (V0M)
1045 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1048 else if(gAnalysisLevel == "MC"){
1049 Double_t gImpactParameter = 0.;
1050 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1052 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1054 gImpactParameter = headerH->ImpactParameter();
1055 gCentrality = gImpactParameter;
1064 // calculate reference multiplicity always (not only in multiplicity mode)
1065 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1066 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1068 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1069 fHistMultiplicity->Fill(gMultiplicity);
1073 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
1074 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1075 if ((fMultiplicityEstimator == "V0M")||
1076 (fMultiplicityEstimator == "V0A")||
1077 (fMultiplicityEstimator == "V0C") ||
1078 (fMultiplicityEstimator == "TPC")) {
1079 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
1080 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
1084 gMultiplicity = header->GetRefMultiplicity();
1085 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1087 fHistMultiplicity->Fill(gMultiplicity);
1089 else if(gAnalysisLevel == "MC") {
1090 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
1091 //Calculating the multiplicity as the number of charged primaries
1092 //within \pm 0.8 in eta and pT > 0.1 GeV/c
1093 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
1094 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
1096 AliError(Form("Could not receive particle %d", iParticle));
1100 //exclude non stable particles
1101 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1104 if (fMultiplicityEstimator == "V0M") {
1105 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
1107 else if (fMultiplicityEstimator == "V0A") {
1108 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
1109 else if (fMultiplicityEstimator == "V0C") {
1110 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
1111 else if (fMultiplicityEstimator == "TPC") {
1112 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1113 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1116 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1117 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1121 if(track->Charge() == 0) continue;
1124 }//loop over primaries
1125 fHistMultiplicity->Fill(gMultiplicity);
1132 // decide what should be returned only here
1133 Double_t lReturnVal = -100;
1134 if(fEventClass=="Multiplicity"){
1135 lReturnVal = gMultiplicity;
1136 }else if(fEventClass=="Centrality"){
1137 lReturnVal = gCentrality;
1142 //________________________________________________________________________
1143 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
1144 //Function that returns the reference multiplicity from AODs (data or reco MC)
1145 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
1146 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
1147 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
1149 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1151 Printf("ERROR: AOD header not available");
1154 Int_t gRunNumber = header->GetRunNumber();
1156 // Loop over tracks in event
1157 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1158 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1160 AliError(Form("Could not receive track %d", iTracks));
1165 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1167 if(aodTrack->Charge() == 0) continue;
1168 // Kinematics cuts from ESD track cuts
1169 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
1170 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
1172 //=================PID (so far only for electron rejection)==========================//
1173 if(fElectronRejection) {
1174 // get the electron nsigma
1175 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1177 // check only for given momentum range
1178 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
1179 //look only at electron nsigma
1180 if(!fElectronOnlyRejection) {
1181 //Make the decision based on the n-sigma of electrons
1182 if(nSigma < fElectronRejectionNSigma) continue;
1185 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1186 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1187 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1189 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1190 if(nSigma < fElectronRejectionNSigma
1191 && nSigmaPions > fElectronRejectionNSigma
1192 && nSigmaKaons > fElectronRejectionNSigma
1193 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1196 }//electron rejection
1198 gRefMultiplicityTPC += 1.0;
1201 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
1202 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
1203 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
1206 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1207 else if(iChannel >= 32)
1208 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1211 //Equalization of gain
1212 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1214 gRefMultiplicityVZEROA /= gFactorA;
1215 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1217 gRefMultiplicityVZEROC /= gFactorC;
1218 if((gFactorA != 0)&&(gFactorC != 0))
1219 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1222 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1224 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1226 if(fMultiplicityEstimator == "TPC")
1227 gRefMultiplicity = gRefMultiplicityTPC;
1228 else if(fMultiplicityEstimator == "V0M")
1229 gRefMultiplicity = gRefMultiplicityVZERO;
1230 else if(fMultiplicityEstimator == "V0A")
1231 gRefMultiplicity = gRefMultiplicityVZEROA;
1232 else if(fMultiplicityEstimator == "V0C")
1233 gRefMultiplicity = gRefMultiplicityVZEROC;
1235 return gRefMultiplicity;
1238 //________________________________________________________________________
1239 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1240 // Get the event plane
1242 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1244 Float_t gVZEROEventPlane = -10.;
1245 Float_t gReactionPlane = -10.;
1246 Double_t qxTot = 0.0, qyTot = 0.0;
1248 //MC: from reaction plane
1249 if(gAnalysisLevel == "MC"){
1251 AliError("mcEvent not available");
1255 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1257 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1259 gReactionPlane = headerH->ReactionPlaneAngle();
1260 //gReactionPlane *= TMath::RadToDeg();
1265 // AOD,ESD,ESDMC: from VZERO Event Plane
1268 AliEventplane *ep = event->GetEventplane();
1270 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1271 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1272 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1273 gReactionPlane = gVZEROEventPlane;
1277 return gReactionPlane;
1280 //________________________________________________________________________
1281 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
1285 Double_t gCentrality) {
1286 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
1288 Double_t correction = 1.;
1289 Int_t gCentralityInt = -1;
1291 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1292 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1298 // centrality not in array --> no correction
1299 if(gCentralityInt < 0){
1304 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1306 if(fHistCorrectionPlus[gCentralityInt]){
1308 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1309 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
1312 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1313 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
1319 }//centrality in array
1321 if (correction == 0.) {
1322 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
1329 //________________________________________________________________________
1330 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
1331 // Returns TObjArray with tracks after all track cuts (only for AOD!)
1332 // Fills QA histograms
1334 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1336 //output TObjArray holding all good tracks
1337 TObjArray* tracksAccepted = new TObjArray;
1338 tracksAccepted->SetOwner(kTRUE);
1347 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1348 // Loop over tracks in event
1350 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1351 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1353 AliError(Form("Could not receive track %d", iTracks));
1359 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
1360 // take only TPC only tracks
1361 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
1362 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1363 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1366 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1369 // additional check on kPrimary flag
1370 if(fCheckPrimaryFlagAOD){
1371 if(aodTrack->GetType() != AliAODTrack::kPrimary)
1376 vCharge = aodTrack->Charge();
1377 vEta = aodTrack->Eta();
1379 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1380 vPt = aodTrack->Pt();
1382 //===========================PID (so far only for electron rejection)===============================//
1383 if(fElectronRejection) {
1385 // get the electron nsigma
1386 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1388 //Fill QA before the PID
1389 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1390 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1391 //end of QA-before pid
1393 // check only for given momentum range
1394 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1396 //look only at electron nsigma
1397 if(!fElectronOnlyRejection){
1399 //Make the decision based on the n-sigma of electrons
1400 if(nSigma < fElectronRejectionNSigma) continue;
1404 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1405 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1406 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1408 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1409 if(nSigma < fElectronRejectionNSigma
1410 && nSigmaPions > fElectronRejectionNSigma
1411 && nSigmaKaons > fElectronRejectionNSigma
1412 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1416 //Fill QA after the PID
1417 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1418 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1421 //===========================end of PID (so far only for electron rejection)===============================//
1423 //+++++++++++++++++++++++++++++//
1424 //===========================PID===============================//
1426 Double_t prob[AliPID::kSPECIES]={0.};
1427 Double_t probTPC[AliPID::kSPECIES]={0.};
1428 Double_t probTOF[AliPID::kSPECIES]={0.};
1429 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1431 Double_t nSigma = 0.;
1432 Double_t nSigmaTPC = 0.;
1433 Double_t nSigmaTOF = 0.;
1434 Double_t nSigmaTPCTOF = 0.;
1435 UInt_t detUsedTPC = 0;
1436 UInt_t detUsedTOF = 0;
1437 UInt_t detUsedTPCTOF = 0;
1439 nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
1440 nSigmaTOF = fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest);
1441 nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
1442 if (nSigmaTOF == 999 || nSigmaTOF == -999){
1443 nSigmaTPCTOF = nSigmaTPC;
1446 //Decide what detector configuration we want to use
1447 switch(fPidDetectorConfig) {
1449 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1451 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1452 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1453 prob[iSpecies] = probTPC[iSpecies];
1456 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1458 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1459 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1460 prob[iSpecies] = probTOF[iSpecies];
1463 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1464 nSigma = nSigmaTPCTOF;
1465 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
1466 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1467 prob[iSpecies] = probTPCTOF[iSpecies];
1471 }//end switch: define detector mask
1473 //Filling the PID QA
1474 Double_t tofTime = -999., length = 999., tof = -999.;
1475 Double_t c = TMath::C()*1.E-9;// m/ns
1476 Double_t beta = -999.;
1478 Float_t probMis = fPIDResponse->GetTOFMismatchProbability(aodTrack);
1479 if (probMis < 0.01) { //if u want to reduce mismatch using also TPC
1481 if ((aodTrack->IsOn(AliAODTrack::kITSin)) && (aodTrack->IsOn(AliAODTrack::kTOFpid)) ) { //leonardo's analysis
1483 tofTime = aodTrack->GetTOFsignal();//in ps
1484 length = aodTrack->GetIntegratedLength();
1485 tof = tofTime*1E-3; // ns
1487 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1491 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1494 length = length*0.01; // in meters
1498 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1499 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1500 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1503 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->GetTPCmomentum()*aodTrack->Charge(),aodTrack->GetTPCsignal()); //aodTrack->P()*aodTrack->Charge()
1504 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1505 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
1506 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1509 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta);
1510 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
1511 fHistNSigmaTPCTOFPbefPID ->Fill(nSigmaTPC,nSigmaTOF,aodTrack->P()); //++++++++++++++
1512 //Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
1513 //end of QA-before pid
1515 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1516 //Make the decision based on the n-sigma
1518 if(nSigma > fPIDNSigma || nSigma < -fPIDNSigma) continue;
1520 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
1521 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
1522 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
1523 fHistNSigmaTPCTOFPafterPID ->Fill(nSigmaTPC,nSigmaTOF,aodTrack->P()); //++++++++++++++
1525 //Make the decision based on the bayesian
1526 else if(fUsePIDPropabilities) {
1527 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1528 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
1530 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1531 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1532 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1535 //Fill QA after the PID
1536 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1537 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1538 fHistBetaVsdEdXafterPID ->Fill(aodTrack->GetTPCsignal(),beta);
1542 //===========================PID===============================//
1543 //+++++++++++++++++++++++++++++//
1546 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1547 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1550 // Kinematics cuts from ESD track cuts
1551 if( vPt < fPtMin || vPt > fPtMax) continue;
1552 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1554 // Extra DCA cuts (for systematic studies [!= -1])
1555 if( fDCAxyCut != -1 && fDCAzCut != -1){
1556 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1561 // Extra TPC cuts (for systematic studies [!= -1])
1562 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1565 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1569 // Extra cut on shared clusters
1570 if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1574 // fill QA histograms
1575 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1576 fHistDCA->Fill(dcaZ,dcaXY);
1577 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1578 fHistPt->Fill(vPt,gCentrality);
1579 fHistEta->Fill(vEta,gCentrality);
1580 fHistRapidity->Fill(vY,gCentrality);
1581 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1582 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1583 fHistPhi->Fill(vPhi,gCentrality);
1584 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1585 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1587 //=======================================correction
1588 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1589 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1591 // add the track to the TObjArray
1592 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1598 else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing)
1599 // Loop over tracks in event
1601 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1602 AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks));
1604 AliError(Form("Could not receive track %d", iTracks));
1608 // AOD track cuts (not needed)
1609 //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1611 vCharge = aodTrack->Charge();
1612 vEta = aodTrack->Eta();
1614 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1615 vPt = aodTrack->Pt();
1618 // Kinematics cuts from ESD track cuts
1619 if( vPt < fPtMin || vPt > fPtMax) continue;
1620 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1623 // fill QA histograms
1624 fHistPt->Fill(vPt,gCentrality);
1625 fHistEta->Fill(vEta,gCentrality);
1626 fHistRapidity->Fill(vY,gCentrality);
1627 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1628 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1629 fHistPhi->Fill(vPhi,gCentrality);
1630 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1631 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1633 //=======================================correction
1634 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1635 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1637 // add the track to the TObjArray
1638 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1640 }// AOD nano analysis
1643 //==============================================================================================================
1644 else if(gAnalysisLevel == "MCAOD") {
1646 AliMCEvent* mcEvent = MCEvent();
1648 AliError("ERROR: Could not retrieve MC event");
1652 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1653 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
1655 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1659 if(!aodTrack->IsPhysicalPrimary()) continue;
1661 vCharge = aodTrack->Charge();
1662 vEta = aodTrack->Eta();
1664 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1665 vPt = aodTrack->Pt();
1667 // Kinematics cuts from ESD track cuts
1668 if( vPt < fPtMin || vPt > fPtMax) continue;
1669 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1671 // Remove neutral tracks
1672 if( vCharge == 0 ) continue;
1674 //Exclude resonances
1675 if(fExcludeResonancesInMC) {
1677 Bool_t kExcludeParticle = kFALSE;
1678 Int_t gMotherIndex = aodTrack->GetMother();
1679 if(gMotherIndex != -1) {
1680 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1682 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1683 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1684 //if(pdgCodeOfMother == 113) {
1685 if(pdgCodeOfMother == 113 // rho0
1686 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1687 // || pdgCodeOfMother == 221 // eta
1688 // || pdgCodeOfMother == 331 // eta'
1689 // || pdgCodeOfMother == 223 // omega
1690 // || pdgCodeOfMother == 333 // phi
1691 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1692 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1693 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1694 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1695 || pdgCodeOfMother == 111 // pi0 Dalitz
1696 || pdgCodeOfMother == 22 // photon
1698 kExcludeParticle = kTRUE;
1703 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1704 if(kExcludeParticle) continue;
1707 //Exclude electrons with PDG
1708 if(fExcludeElectronsInMC) {
1710 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1714 // fill QA histograms
1715 fHistPt->Fill(vPt,gCentrality);
1716 fHistEta->Fill(vEta,gCentrality);
1717 fHistRapidity->Fill(vY,gCentrality);
1718 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1719 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1720 fHistPhi->Fill(vPhi,gCentrality);
1721 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1722 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1724 //=======================================correction
1725 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1726 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1728 // add the track to the TObjArray
1729 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1733 //==============================================================================================================
1735 //==============================================================================================================
1736 else if(gAnalysisLevel == "MCAODrec") {
1738 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1740 printf("ERROR: fAOD not available\n");
1744 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
1746 AliError("No array of MC particles found !!!");
1749 AliMCEvent* mcEvent = MCEvent();
1751 AliError("ERROR: Could not retrieve MC event");
1752 return tracksAccepted;
1755 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1756 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1758 AliError(Form("Could not receive track %d", iTracks));
1762 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1763 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1765 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1767 vCharge = aodTrack->Charge();
1768 vEta = aodTrack->Eta();
1770 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1771 vPt = aodTrack->Pt();
1773 //===========================use MC information for Kinematics===============================//
1774 if(fUseMCforKinematics){
1776 Int_t label = TMath::Abs(aodTrack->GetLabel());
1777 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1780 vCharge = AODmcTrack->Charge();
1781 vEta = AODmcTrack->Eta();
1782 vY = AODmcTrack->Y();
1783 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
1784 vPt = AODmcTrack->Pt();
1787 AliDebug(1, "no MC particle for this track");
1791 //===========================end of use MC information for Kinematics========================//
1794 //===========================PID (so far only for electron rejection)===============================//
1795 if(fElectronRejection) {
1797 // get the electron nsigma
1798 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1800 //Fill QA before the PID
1801 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1802 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1803 //end of QA-before pid
1805 // check only for given momentum range
1806 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1808 //look only at electron nsigma
1809 if(!fElectronOnlyRejection){
1811 //Make the decision based on the n-sigma of electrons
1812 if(nSigma < fElectronRejectionNSigma) continue;
1816 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1817 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1818 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1820 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1821 if(nSigma < fElectronRejectionNSigma
1822 && nSigmaPions > fElectronRejectionNSigma
1823 && nSigmaKaons > fElectronRejectionNSigma
1824 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1828 //Fill QA after the PID
1829 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1830 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1833 //===========================end of PID (so far only for electron rejection)===============================//
1835 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1836 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1838 // Kinematics cuts from ESD track cuts
1839 if( vPt < fPtMin || vPt > fPtMax) continue;
1840 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1842 // Extra DCA cuts (for systematic studies [!= -1])
1843 if( fDCAxyCut != -1 && fDCAzCut != -1){
1844 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1849 // Extra TPC cuts (for systematic studies [!= -1])
1850 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1853 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1857 // Extra cut on shared clusters
1858 if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1862 //Exclude resonances
1863 if(fExcludeResonancesInMC) {
1865 Bool_t kExcludeParticle = kFALSE;
1867 Int_t label = TMath::Abs(aodTrack->GetLabel());
1868 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1871 //if (AODmcTrack->IsPhysicalPrimary()){
1873 Int_t gMotherIndex = AODmcTrack->GetMother();
1874 if(gMotherIndex != -1) {
1875 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1877 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1878 if(pdgCodeOfMother == 113 // rho0
1879 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1880 // || pdgCodeOfMother == 221 // eta
1881 // || pdgCodeOfMother == 331 // eta'
1882 // || pdgCodeOfMother == 223 // omega
1883 // || pdgCodeOfMother == 333 // phi
1884 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1885 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1886 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1887 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1888 || pdgCodeOfMother == 111 // pi0 Dalitz
1889 || pdgCodeOfMother == 22 // photon
1891 kExcludeParticle = kTRUE;
1896 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1897 if(kExcludeParticle) continue;
1900 //Exclude electrons with PDG
1901 if(fExcludeElectronsInMC) {
1903 Int_t label = TMath::Abs(aodTrack->GetLabel());
1904 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1907 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1911 // fill QA histograms
1912 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1913 fHistDCA->Fill(dcaZ,dcaXY);
1914 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1915 fHistPt->Fill(vPt,gCentrality);
1916 fHistEta->Fill(vEta,gCentrality);
1917 fHistRapidity->Fill(vY,gCentrality);
1918 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1919 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1920 fHistPhi->Fill(vPhi,gCentrality);
1921 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1922 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1924 //=======================================correction
1925 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1926 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1928 // add the track to the TObjArray
1929 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1932 //==============================================================================================================
1934 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1936 AliESDtrack *trackTPC = NULL;
1937 AliMCParticle *mcTrack = NULL;
1939 AliMCEvent* mcEvent = NULL;
1941 // for MC ESDs use also MC information
1942 if(gAnalysisLevel == "MCESD"){
1943 mcEvent = MCEvent();
1945 AliError("mcEvent not available");
1946 return tracksAccepted;
1950 // Loop over tracks in event
1951 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1952 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
1954 AliError(Form("Could not receive track %d", iTracks));
1958 // for MC ESDs use also MC information --> MC track not used further???
1959 if(gAnalysisLevel == "MCESD"){
1960 Int_t label = TMath::Abs(track->GetLabel());
1961 if(label > mcEvent->GetNumberOfTracks()) continue;
1962 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1964 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1965 if(!mcTrack) continue;
1968 // take only TPC only tracks
1969 trackTPC = new AliESDtrack();
1970 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1974 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1976 // fill QA histograms
1979 trackTPC->GetImpactParameters(b,bCov);
1980 if (bCov[0]<=0 || bCov[2]<=0) {
1981 AliDebug(1, "Estimated b resolution lower or equal zero!");
1982 bCov[0]=0; bCov[2]=0;
1985 Int_t nClustersTPC = -1;
1986 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1987 //nClustersTPC = track->GetTPCclusters(0); // global track
1988 Float_t chi2PerClusterTPC = -1;
1989 if (nClustersTPC!=0) {
1990 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1991 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1994 //===========================PID===============================//
1996 Double_t prob[AliPID::kSPECIES]={0.};
1997 Double_t probTPC[AliPID::kSPECIES]={0.};
1998 Double_t probTOF[AliPID::kSPECIES]={0.};
1999 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
2001 Double_t nSigma = 0.;
2002 UInt_t detUsedTPC = 0;
2003 UInt_t detUsedTOF = 0;
2004 UInt_t detUsedTPCTOF = 0;
2006 //Decide what detector configuration we want to use
2007 switch(fPidDetectorConfig) {
2009 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
2010 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
2011 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
2012 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
2013 prob[iSpecies] = probTPC[iSpecies];
2016 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
2017 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
2018 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
2019 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
2020 prob[iSpecies] = probTOF[iSpecies];
2023 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
2024 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
2025 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
2026 prob[iSpecies] = probTPCTOF[iSpecies];
2030 }//end switch: define detector mask
2032 //Filling the PID QA
2033 Double_t tofTime = -999., length = 999., tof = -999.;
2034 Double_t c = TMath::C()*1.E-9;// m/ns
2035 Double_t beta = -999.;
2036 Double_t nSigmaTOFForParticleOfInterest = -999.;
2037 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
2038 (track->IsOn(AliESDtrack::kTIME)) ) {
2039 tofTime = track->GetTOFsignal();//in ps
2040 length = track->GetIntegratedLength();
2041 tof = tofTime*1E-3; // ns
2044 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
2048 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
2052 length = length*0.01; // in meters
2056 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
2057 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
2058 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2059 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2063 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
2064 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2065 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2066 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2067 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2068 //end of QA-before pid
2070 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
2071 //Make the decision based on the n-sigma
2073 if(nSigma > fPIDNSigma) continue;}
2075 //Make the decision based on the bayesian
2076 else if(fUsePIDPropabilities) {
2077 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
2078 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
2081 //Fill QA after the PID
2082 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
2083 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2084 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2086 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2087 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2088 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2089 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2092 //===========================PID===============================//
2093 vCharge = trackTPC->Charge();
2095 vEta = trackTPC->Eta();
2096 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
2097 vPt = trackTPC->Pt();
2098 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
2099 fHistDCA->Fill(b[1],b[0]);
2100 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
2101 fHistPt->Fill(vPt,gCentrality);
2102 fHistEta->Fill(vEta,gCentrality);
2103 fHistPhi->Fill(vPhi,gCentrality);
2104 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2105 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2106 fHistRapidity->Fill(vY,gCentrality);
2107 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2108 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2110 //=======================================correction
2111 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2112 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2114 // add the track to the TObjArray
2115 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2121 else if(gAnalysisLevel == "MC"){
2123 AliError("mcEvent not available");
2127 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2129 // Loop over tracks in event
2130 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
2131 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
2133 AliError(Form("Could not receive particle %d", iTracks));
2137 //exclude non stable particles
2138 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2140 vCharge = track->Charge();
2141 vEta = track->Eta();
2145 if( vPt < fPtMin || vPt > fPtMax)
2148 if( vEta < fEtaMin || vEta > fEtaMax) continue;
2151 if( vY < fEtaMin || vY > fEtaMax) continue;
2154 // Remove neutral tracks
2155 if( vCharge == 0 ) continue;
2157 //analyze one set of particles
2159 TParticle *particle = track->Particle();
2160 if(!particle) continue;
2162 Int_t gPdgCode = particle->GetPdgCode();
2163 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
2167 //Use the acceptance parameterization
2168 if(fAcceptanceParameterization) {
2169 Double_t gRandomNumber = gRandom->Rndm();
2170 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
2174 //Exclude resonances
2175 if(fExcludeResonancesInMC) {
2176 TParticle *particle = track->Particle();
2177 if(!particle) continue;
2179 Bool_t kExcludeParticle = kFALSE;
2180 Int_t gMotherIndex = particle->GetFirstMother();
2181 if(gMotherIndex != -1) {
2182 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2184 TParticle *motherParticle = motherTrack->Particle();
2185 if(motherParticle) {
2186 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2187 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2188 if(pdgCodeOfMother == 113 // rho0
2189 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2190 // || pdgCodeOfMother == 221 // eta
2191 // || pdgCodeOfMother == 331 // eta'
2192 // || pdgCodeOfMother == 223 // omega
2193 // || pdgCodeOfMother == 333 // phi
2194 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
2195 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
2196 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
2197 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2198 || pdgCodeOfMother == 111 // pi0 Dalitz
2200 kExcludeParticle = kTRUE;
2206 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2207 if(kExcludeParticle) continue;
2210 //Exclude electrons with PDG
2211 if(fExcludeElectronsInMC) {
2213 TParticle *particle = track->Particle();
2216 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2220 vPhi = track->Phi();
2221 //Printf("phi (before): %lf",vPhi);
2223 fHistPt->Fill(vPt,gCentrality);
2224 fHistEta->Fill(vEta,gCentrality);
2225 fHistPhi->Fill(vPhi,gCentrality);
2226 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2227 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2228 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2229 fHistRapidity->Fill(vY,gCentrality);
2230 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2231 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2232 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2233 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2236 if(fUseFlowAfterBurner) {
2237 Double_t precisionPhi = 0.001;
2238 Int_t maxNumberOfIterations = 100;
2240 Double_t phi0 = vPhi;
2241 Double_t gV2 = fDifferentialV2->Eval(vPt);
2243 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2244 Double_t phiprev = vPhi;
2245 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2246 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2248 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2250 //Printf("phi (after): %lf\n",vPhi);
2251 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2252 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2253 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2255 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2256 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2257 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2261 //vPhi *= TMath::RadToDeg();
2263 //=======================================correction
2264 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2265 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2267 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2272 return tracksAccepted;
2275 //________________________________________________________________________
2276 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2277 // Clones TObjArray and returns it with tracks after shuffling the charges
2279 TObjArray* tracksShuffled = new TObjArray;
2280 tracksShuffled->SetOwner(kTRUE);
2282 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
2284 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2286 AliVParticle* track = (AliVParticle*) tracks->At(i);
2287 chargeVector->push_back(track->Charge());
2290 random_shuffle(chargeVector->begin(), chargeVector->end());
2292 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2293 AliVParticle* track = (AliVParticle*) tracks->At(i);
2294 //==============================correction
2295 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2296 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2297 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2300 delete chargeVector;
2302 return tracksShuffled;
2305 //________________________________________________________________________
2306 void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2307 const char* lhcPeriod) {
2308 //Function to setup the VZERO gain equalization
2309 //============Get the equilization map============//
2310 TFile *calibrationFile = TFile::Open(filename);
2311 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2312 Printf("No calibration file found!!!");
2316 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2318 Printf("Calibration TList not found!!!");
2322 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2323 if(!fHistVZEROAGainEqualizationMap) {
2324 Printf("VZERO-A calibration object not found!!!");
2327 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2328 if(!fHistVZEROCGainEqualizationMap) {
2329 Printf("VZERO-C calibration object not found!!!");
2333 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2334 if(!fHistVZEROChannelGainEqualizationMap) {
2335 Printf("VZERO channel calibration object not found!!!");
2340 //________________________________________________________________________
2341 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
2344 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2346 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2347 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2348 if(gRunNumber == run)
2349 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2355 //________________________________________________________________________
2356 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
2359 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2361 TString gVZEROSide = side;
2362 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2363 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2364 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2365 if(gRunNumber == run) {
2366 if(gVZEROSide == "A")
2367 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2368 else if(gVZEROSide == "C")
2369 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2376 //____________________________________________________________________
2377 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2379 // copied from AliAnalysisTaskPhiCorrelations
2381 // rejects "randomly" events such that the centrality gets flat
2382 // uses fCentralityWeights histogram
2384 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2386 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2387 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2389 Bool_t result = kFALSE;
2390 if (centralityDigits < weight)
2393 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2398 //________________________________________________________________________
2399 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
2403 AliError("fBalance not available");
2407 if (!fShuffledBalance) {
2408 AliError("fShuffledBalance not available");
2415 //________________________________________________________________________
2416 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2417 // Draw result to the screen
2418 // Called once at the end of the query
2420 // not implemented ...