4 #include "TLorentzVector.h"
5 #include "TGraphErrors.h"
15 #include "AliAnalysisTaskSE.h"
16 #include "AliAnalysisManager.h"
18 #include "AliESDVertex.h"
19 #include "AliESDEvent.h"
20 #include "AliESDInputHandler.h"
21 #include "AliAODEvent.h"
22 #include "AliAODTrack.h"
23 #include "AliAODInputHandler.h"
24 #include "AliAODMCParticle.h"
25 #include "AliCollisionGeometry.h"
26 #include "AliGenEventHeader.h"
27 #include "AliMCEventHandler.h"
28 #include "AliMCEvent.h"
30 #include "AliESDtrackCuts.h"
31 #include "AliEventplane.h"
34 #include "AliAnalysisUtils.h"
36 #include "AliEventPoolManager.h"
39 #include "AliPIDResponse.h"
40 #include "AliPIDCombined.h"
42 #include "AliAnalysisTaskBFPsi.h"
43 #include "AliBalancePsi.h"
44 #include "AliAnalysisTaskTriggeredBF.h"
47 // Analysis task for the BF vs Psi code
48 // Authors: Panos.Christakoglou@nikhef.nl
53 ClassImp(AliAnalysisTaskBFPsi)
55 //________________________________________________________________________
56 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
57 : AliAnalysisTaskSE(name),
59 fArrayMC(0), //+++++++++++++
61 fRunShuffling(kFALSE),
64 fRunMixingEventPlane(kFALSE),
75 fHistCentStatsUsed(0),
81 fHistTPCvsVZEROMultiplicity(0),
99 fHistdEdxVsPTPCbeforePID(NULL),
100 fHistBetavsPTOFbeforePID(NULL),
101 fHistProbTPCvsPtbeforePID(NULL),
102 fHistProbTOFvsPtbeforePID(NULL),
103 fHistProbTPCTOFvsPtbeforePID(NULL),
104 fHistNSigmaTPCvsPtbeforePID(NULL),
105 fHistNSigmaTOFvsPtbeforePID(NULL),
106 fHistBetaVsdEdXbeforePID(NULL), //+++++++
107 fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++
108 fHistdEdxVsPTPCafterPID(NULL),
109 fHistBetavsPTOFafterPID(NULL),
110 fHistProbTPCvsPtafterPID(NULL),
111 fHistProbTOFvsPtafterPID(NULL),
112 fHistProbTPCTOFvsPtafterPID(NULL),
113 fHistNSigmaTPCvsPtafterPID(NULL),
114 fHistNSigmaTOFvsPtafterPID(NULL),
115 fHistBetaVsdEdXafterPID(NULL), //+++++++
116 fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++
117 fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++
118 fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++
119 fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++
120 fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++
121 fCentralityArrayBinsForCorrections(kCENTRALITY),
122 fCentralityWeights(0x0),
125 fParticleOfInterest(kPion),
126 fPidDetectorConfig(kTPCTOF),
128 fUsePIDnSigma(kTRUE),
129 fUsePIDPropabilities(kFALSE),
131 fMinAcceptedPIDProbability(0.8),
132 fElectronRejection(kFALSE),
133 fElectronOnlyRejection(kFALSE),
134 fElectronRejectionNSigma(-1.),
135 fElectronRejectionMinPt(0.),
136 fElectronRejectionMaxPt(1000.),
138 fCentralityEstimator("V0M"),
139 fUseCentrality(kFALSE),
140 fCentralityPercentileMin(0.),
141 fCentralityPercentileMax(5.),
142 fImpactParameterMin(0.),
143 fImpactParameterMax(20.),
144 fMultiplicityEstimator("V0A"),
145 fUseMultiplicity(kFALSE),
146 fNumberOfAcceptedTracksMin(0),
147 fNumberOfAcceptedTracksMax(10000),
148 fHistNumberOfAcceptedTracks(0),
149 fHistMultiplicity(0),
150 fUseOfflineTrigger(kFALSE),
151 fCheckFirstEventInChunk(kFALSE),
152 fCheckPileUp(kFALSE),
153 fCheckPrimaryFlagAOD(kFALSE),
154 fUseMCforKinematics(kFALSE),
158 fnAODtrackCutBit(128),
168 fNClustersTPCCut(-1),
169 fAcceptanceParameterization(0),
171 fUseFlowAfterBurner(kFALSE),
172 fExcludeResonancesInMC(kFALSE),
173 fExcludeElectronsInMC(kFALSE),
174 fUseMCPdgCode(kFALSE),
175 fPDGCodeToBeAnalyzed(-1),
176 fEventClass("EventPlane"),
178 fHistVZEROAGainEqualizationMap(0),
179 fHistVZEROCGainEqualizationMap(0),
180 fHistVZEROChannelGainEqualizationMap(0) {
182 // Define input and output slots here
183 // Input slot #0 works with a TChain
185 //======================================================correction
186 for (Int_t i=0; i<kCENTRALITY; i++){
187 fHistCorrectionPlus[i] = NULL;
188 fHistCorrectionMinus[i] = NULL;
189 fCentralityArrayForCorrections[i] = -1.;
191 //=====================================================correction
193 DefineInput(0, TChain::Class());
194 // Output slot #0 writes into a TH1 container
195 DefineOutput(1, TList::Class());
196 DefineOutput(2, TList::Class());
197 DefineOutput(3, TList::Class());
198 DefineOutput(4, TList::Class());
199 DefineOutput(5, TList::Class());
202 //________________________________________________________________________
203 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
206 // delete fShuffledBalance;
211 // delete fHistEventStats;
212 // delete fHistTrackStats;
223 // delete fHistEtaPhiPos;
224 // delete fHistEtaPhiNeg;
228 //________________________________________________________________________
229 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
233 // global switch disabling the reference
234 // (to avoid "Replacing existing TH1" if several wagons are created in train)
235 Bool_t oldStatus = TH1::AddDirectoryStatus();
236 TH1::AddDirectory(kFALSE);
239 fBalance = new AliBalancePsi();
240 fBalance->SetAnalysisLevel("ESD");
241 fBalance->SetEventClass(fEventClass);
242 //fBalance->SetNumberOfBins(-1,16);
243 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
246 if(!fShuffledBalance) {
247 fShuffledBalance = new AliBalancePsi();
248 fShuffledBalance->SetAnalysisLevel("ESD");
249 fShuffledBalance->SetEventClass(fEventClass);
250 //fShuffledBalance->SetNumberOfBins(-1,16);
251 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
256 fMixedBalance = new AliBalancePsi();
257 fMixedBalance->SetAnalysisLevel("ESD");
258 fMixedBalance->SetEventClass(fEventClass);
264 fList->SetName("listQA");
267 //Balance Function list
268 fListBF = new TList();
269 fListBF->SetName("listBF");
273 fListBFS = new TList();
274 fListBFS->SetName("listBFShuffled");
275 fListBFS->SetOwner();
279 fListBFM = new TList();
280 fListBFM->SetName("listTriggeredBFMixed");
281 fListBFM->SetOwner();
285 if(fUsePID || fElectronRejection) {
286 fHistListPIDQA = new TList();
287 fHistListPIDQA->SetName("listQAPID");
288 fHistListPIDQA->SetOwner();
292 TString gCutName[7] = {"Total","Offline trigger",
293 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
294 fHistEventStats = new TH2F("fHistEventStats",
295 "Event statistics;;Centrality percentile;N_{events}",
296 7,0.5,7.5,220,-5,105);
297 for(Int_t i = 1; i <= 7; i++)
298 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
299 fList->Add(fHistEventStats);
301 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
302 fHistCentStats = new TH2F("fHistCentStats",
303 "Centrality statistics;;Cent percentile",
304 13,-0.5,12.5,220,-5,105);
305 for(Int_t i = 1; i <= 13; i++){
306 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
307 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
309 fList->Add(fHistCentStats);
311 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
312 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
313 fList->Add(fHistCentStatsUsed);
315 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
316 fList->Add(fHistTriggerStats);
318 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
319 fList->Add(fHistTrackStats);
321 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
322 fList->Add(fHistNumberOfAcceptedTracks);
324 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
325 fList->Add(fHistMultiplicity);
327 // Vertex distributions
328 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
330 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
332 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
335 //TPC vs VZERO multiplicity
336 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
337 if(fMultiplicityEstimator == "V0A")
338 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
339 else if(fMultiplicityEstimator == "V0C")
340 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
342 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
343 fList->Add(fHistTPCvsVZEROMultiplicity);
345 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
346 fList->Add(fHistVZEROSignal);
349 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
350 fList->Add(fHistEventPlane);
353 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
354 fList->Add(fHistClus);
355 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
356 fList->Add(fHistChi2);
357 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
358 fList->Add(fHistDCA);
359 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
361 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
362 fList->Add(fHistEta);
363 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
364 fList->Add(fHistRapidity);
365 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
366 fList->Add(fHistPhi);
367 fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
368 fList->Add(fHistEtaPhiPos);
369 fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);
370 fList->Add(fHistEtaPhiNeg);
371 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
372 fList->Add(fHistPhiBefore);
373 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
374 fList->Add(fHistPhiAfter);
375 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
376 fList->Add(fHistPhiPos);
377 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
378 fList->Add(fHistPhiNeg);
379 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
380 fList->Add(fHistV0M);
381 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
382 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
383 for(Int_t i = 1; i <= 6; i++)
384 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
385 fList->Add(fHistRefTracks);
387 // Balance function histograms
388 // Initialize histograms if not done yet (including the custom binning)
389 if(!fBalance->GetHistNp()){
390 AliInfo("Histograms not yet initialized! --> Will be done now");
391 fBalance->SetCustomBinning(fCustomBinning);
392 fBalance->InitHistograms();
396 if(!fShuffledBalance->GetHistNp()) {
397 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
398 fShuffledBalance->SetCustomBinning(fCustomBinning);
399 fShuffledBalance->InitHistograms();
404 if(!fMixedBalance->GetHistNp()) {
405 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
406 fMixedBalance->SetCustomBinning(fCustomBinning);
407 fMixedBalance->InitHistograms();
411 // QA histograms for different cuts
412 fList->Add(fBalance->GetQAHistHBTbefore());
413 fList->Add(fBalance->GetQAHistHBTafter());
414 fList->Add(fBalance->GetQAHistConversionbefore());
415 fList->Add(fBalance->GetQAHistConversionafter());
416 fList->Add(fBalance->GetQAHistPsiMinusPhi());
417 fList->Add(fBalance->GetQAHistResonancesBefore());
418 fList->Add(fBalance->GetQAHistResonancesRho());
419 fList->Add(fBalance->GetQAHistResonancesK0());
420 fList->Add(fBalance->GetQAHistResonancesLambda());
421 fList->Add(fBalance->GetQAHistQbefore());
422 fList->Add(fBalance->GetQAHistQafter());
424 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
425 fListBF->Add(fBalance->GetHistNp());
426 fListBF->Add(fBalance->GetHistNn());
427 fListBF->Add(fBalance->GetHistNpn());
428 fListBF->Add(fBalance->GetHistNnn());
429 fListBF->Add(fBalance->GetHistNpp());
430 fListBF->Add(fBalance->GetHistNnp());
433 fListBFS->Add(fShuffledBalance->GetHistNp());
434 fListBFS->Add(fShuffledBalance->GetHistNn());
435 fListBFS->Add(fShuffledBalance->GetHistNpn());
436 fListBFS->Add(fShuffledBalance->GetHistNnn());
437 fListBFS->Add(fShuffledBalance->GetHistNpp());
438 fListBFS->Add(fShuffledBalance->GetHistNnp());
442 fListBFM->Add(fMixedBalance->GetHistNp());
443 fListBFM->Add(fMixedBalance->GetHistNn());
444 fListBFM->Add(fMixedBalance->GetHistNpn());
445 fListBFM->Add(fMixedBalance->GetHistNnn());
446 fListBFM->Add(fMixedBalance->GetHistNpp());
447 fListBFM->Add(fMixedBalance->GetHistNnp());
454 Int_t trackDepth = fMixingTracks;
455 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
458 Double_t* centbins = NULL;
459 Int_t nCentralityBins;
460 if(fBalance->IsUseVertexBinning()){
461 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
464 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
468 Double_t* multbins = NULL;
469 Int_t nMultiplicityBins;
470 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
473 Double_t* vtxbins = NULL;
475 if(fBalance->IsUseVertexBinning()){
476 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
479 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
482 // Event plane angle (Psi) bins
483 Double_t* psibins = NULL;
485 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
488 // run the event mixing also in bins of event plane (statistics!)
489 if(fRunMixingEventPlane){
490 if(fEventClass=="Multiplicity"){
491 if(multbins && vtxbins && psibins){
492 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
496 if(centbins && vtxbins && psibins){
497 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
502 if(fEventClass=="Multiplicity"){
503 if(multbins && vtxbins){
504 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
508 if(centbins && vtxbins){
509 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
514 if(centbins) delete [] centbins;
515 if(multbins) delete [] multbins;
516 if(vtxbins) delete [] vtxbins;
518 // check pool manager
520 AliError("Event Mixing required, but Pool Manager not initialized...");
525 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
527 //====================PID========================//
529 fPIDCombined = new AliPIDCombined();
530 fPIDCombined->SetDefaultTPCPriors();
532 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
533 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
535 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
536 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
538 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
539 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
541 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
542 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
544 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
545 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
547 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
548 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
550 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
551 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
553 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
554 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
556 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500);
557 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
559 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
560 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
562 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
563 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
565 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
566 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
568 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
569 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
571 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
572 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
574 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
575 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
577 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
578 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
580 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
581 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
583 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500);
584 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
587 // for electron rejection only TPC nsigma histograms
588 if(fElectronRejection) {
590 fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
591 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
593 fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
594 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
596 fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
597 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
599 fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
600 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
602 //====================PID========================//
606 PostData(2, fListBF);
607 if(fRunShuffling) PostData(3, fListBFS);
608 if(fRunMixing) PostData(4, fListBFM);
609 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
611 AliInfo("Finished setting up the Output");
613 TH1::AddDirectory(oldStatus);
618 //________________________________________________________________________
619 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
620 Int_t nCentralityBins,
621 Double_t *centralityArrayForCorrections) {
622 //Open files that will be used for correction
623 fCentralityArrayBinsForCorrections = nCentralityBins;
624 for (Int_t i=0; i<nCentralityBins; i++)
625 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
627 // No file specified -> run without corrections
628 if(!filename.Contains(".root")) {
629 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
633 //Open the input file
634 TFile *f = TFile::Open(filename);
636 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
640 //TString listEffName = "";
641 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
642 //Printf("iCent %d:",iCent);
643 TString histoName = "fHistCorrectionPlus";
644 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
645 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
646 if(!fHistCorrectionPlus[iCent]) {
647 AliError(Form("fHist %s not found!!!",histoName.Data()));
651 histoName = "fHistCorrectionMinus";
652 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
653 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
654 if(!fHistCorrectionMinus[iCent]) {
655 AliError(Form("fHist %s not found!!!",histoName.Data()));
658 }//loop over centralities: ONLY the PbPb case is covered
661 //________________________________________________________________________
662 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
664 // Called for each event
666 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
667 Int_t gNumberOfAcceptedTracks = 0;
668 Double_t lMultiplicityVar = -999.; //-1
669 Double_t gReactionPlane = -1.;
672 // get the event (for generator level: MCEvent())
673 AliVEvent* eventMain = NULL;
674 if(gAnalysisLevel == "MC") {
675 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
678 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
679 // for HBT like cuts need magnetic field sign
680 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
683 AliError("eventMain not available");
687 // PID Response task active?
688 if(fUsePID || fElectronRejection) {
689 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
690 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
693 // check event cuts and fill event histograms
694 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
697 // get the reaction plane
698 if(fEventClass != "Multiplicity") {
699 gReactionPlane = GetEventPlane(eventMain);
700 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
701 if(gReactionPlane < 0){
706 // get the accepted tracks in main event
707 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
708 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
710 //multiplicity cut (used in pp)
711 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
713 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
714 TObjArray* tracksShuffled = NULL;
716 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
722 // 1. First get an event pool corresponding in mult (cent) and
723 // zvertex to the current event. Once initialized, the pool
724 // should contain nMix (reduced) events. This routine does not
725 // pre-scan the chain. The first several events of every chain
726 // will be skipped until the needed pools are filled to the
727 // specified depth. If the pool categories are not too rare, this
728 // should not be a problem. If they are rare, you could lose`
731 // 2. Collect the whole pool's content of tracks into one TObjArray
732 // (bgTracks), which is effectively a single background super-event.
734 // 3. The reduced and bgTracks arrays must both be passed into
735 // FillCorrelations(). Also nMix should be passed in, so a weight
736 // of 1./nMix can be applied.
738 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
741 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
747 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
750 Int_t nMix = pool->GetCurrentNEvents();
751 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
753 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
754 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
755 //if (pool->IsReady())
756 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
758 // Fill mixed-event histos here
759 for (Int_t jMix=0; jMix<nMix; jMix++)
761 TObjArray* tracksMixed = pool->GetEvent(jMix);
762 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
766 // Update the Event pool
767 pool->UpdatePool(tracksMain);
773 // calculate balance function
774 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
776 // calculate shuffled balance function
777 if(fRunShuffling && tracksShuffled != NULL) {
778 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
782 //________________________________________________________________________
783 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
784 // Checks the Event cuts
785 // Fills Event statistics histograms
787 Bool_t isSelectedMain = kTRUE;
788 Float_t gRefMultiplicity = -1.;
789 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
791 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
793 fHistEventStats->Fill(1,gRefMultiplicity); //all events
795 // check first event in chunk (is not needed for new reconstructions)
796 if(fCheckFirstEventInChunk){
798 if(ut.IsFirstEventInChunk(event))
800 fHistEventStats->Fill(6,gRefMultiplicity);
802 // check for pile-up event
805 ut.SetUseMVPlpSelection(kTRUE);
806 ut.SetUseOutOfBunchPileUp(kTRUE);
807 if(ut.IsPileUpEvent(event))
809 fHistEventStats->Fill(7,gRefMultiplicity);
812 // Event trigger bits
813 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
814 if(fUseOfflineTrigger)
815 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
818 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
821 if(gAnalysisLevel == "MC") {
823 AliError("mcEvent not available");
828 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
830 TArrayF gVertexArray;
831 header->PrimaryVertex(gVertexArray);
832 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
833 //gVertexArray.At(0),
834 //gVertexArray.At(1),
835 //gVertexArray.At(2));
836 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
837 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
838 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
839 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
840 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
842 // get the reference multiplicty or centrality
843 gRefMultiplicity = GetRefMultiOrCentrality(event);
845 fHistVx->Fill(gVertexArray.At(0));
846 fHistVy->Fill(gVertexArray.At(1));
847 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
849 // take only events inside centrality class
851 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
852 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
853 return gRefMultiplicity;
856 // take events only within the same multiplicity class
857 else if(fUseMultiplicity) {
858 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
859 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
860 return gRefMultiplicity;
862 }//multiplicity range
870 // Event Vertex AOD, ESD, ESDMC
872 const AliVVertex *vertex = event->GetPrimaryVertex();
876 vertex->GetCovarianceMatrix(fCov);
877 if(vertex->GetNContributors() > 0) {
879 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
880 if(TMath::Abs(vertex->GetX()) < fVxMax) {
881 if(TMath::Abs(vertex->GetY()) < fVyMax) {
882 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
883 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
885 // get the reference multiplicty or centrality
886 gRefMultiplicity = GetRefMultiOrCentrality(event);
888 fHistVx->Fill(vertex->GetX());
889 fHistVy->Fill(vertex->GetY());
890 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
892 // take only events inside centrality class
894 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
896 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
897 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
898 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
902 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
903 return gRefMultiplicity;
906 // take events only within the same multiplicity class
907 else if(fUseMultiplicity) {
909 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
910 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
912 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
913 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
914 return gRefMultiplicity;
916 }//multiplicity range
920 }//proper vertex resolution
921 }//proper number of contributors
922 }//vertex object valid
926 // in all other cases return -1 (event not accepted)
931 //________________________________________________________________________
932 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
933 // Checks the Event cuts
934 // Fills Event statistics histograms
936 Float_t gCentrality = -1.;
937 Double_t gMultiplicity = -1.;
938 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
941 // calculate centrality always (not only in centrality mode)
942 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
943 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
945 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
947 // QA for centrality estimators
948 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
949 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
950 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
951 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
952 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
953 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
954 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
955 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
956 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
957 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
958 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
959 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
960 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
962 // Centrality estimator USED ++++++++++++++++++++++++++++++
963 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
965 // centrality QA (V0M)
966 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
968 // centrality QA (reference tracks)
969 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
970 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
971 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
972 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
973 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
974 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
975 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
976 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
977 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
982 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
983 AliCentrality *centrality = event->GetCentrality();
984 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
986 // QA for centrality estimators
987 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
988 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
989 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
990 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
991 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
992 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
993 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
994 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
995 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
996 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
997 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
998 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
999 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
1001 // Centrality estimator USED ++++++++++++++++++++++++++++++
1002 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
1004 // centrality QA (V0M)
1005 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1008 else if(gAnalysisLevel == "MC"){
1009 Double_t gImpactParameter = 0.;
1010 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1012 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1014 gImpactParameter = headerH->ImpactParameter();
1015 gCentrality = gImpactParameter;
1024 // calculate reference multiplicity always (not only in multiplicity mode)
1025 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1026 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1028 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1029 fHistMultiplicity->Fill(gMultiplicity);
1033 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
1034 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1035 if ((fMultiplicityEstimator == "V0M")||
1036 (fMultiplicityEstimator == "V0A")||
1037 (fMultiplicityEstimator == "V0C") ||
1038 (fMultiplicityEstimator == "TPC")) {
1039 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
1040 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
1044 gMultiplicity = header->GetRefMultiplicity();
1045 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1047 fHistMultiplicity->Fill(gMultiplicity);
1049 else if(gAnalysisLevel == "MC") {
1050 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
1051 //Calculating the multiplicity as the number of charged primaries
1052 //within \pm 0.8 in eta and pT > 0.1 GeV/c
1053 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
1054 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
1056 AliError(Form("Could not receive particle %d", iParticle));
1060 //exclude non stable particles
1061 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1064 if (fMultiplicityEstimator == "V0M") {
1065 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
1067 else if (fMultiplicityEstimator == "V0A") {
1068 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
1069 else if (fMultiplicityEstimator == "V0C") {
1070 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
1071 else if (fMultiplicityEstimator == "TPC") {
1072 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1073 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1076 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1077 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1081 if(track->Charge() == 0) continue;
1084 }//loop over primaries
1085 fHistMultiplicity->Fill(gMultiplicity);
1092 // decide what should be returned only here
1093 Double_t lReturnVal = -100;
1094 if(fEventClass=="Multiplicity"){
1095 lReturnVal = gMultiplicity;
1096 }else if(fEventClass=="Centrality"){
1097 lReturnVal = gCentrality;
1102 //________________________________________________________________________
1103 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
1104 //Function that returns the reference multiplicity from AODs (data or reco MC)
1105 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
1106 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
1107 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
1109 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1111 Printf("ERROR: AOD header not available");
1114 Int_t gRunNumber = header->GetRunNumber();
1116 // Loop over tracks in event
1117 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1118 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1120 AliError(Form("Could not receive track %d", iTracks));
1125 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1127 if(aodTrack->Charge() == 0) continue;
1128 // Kinematics cuts from ESD track cuts
1129 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
1130 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
1132 //=================PID (so far only for electron rejection)==========================//
1133 if(fElectronRejection) {
1134 // get the electron nsigma
1135 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1137 // check only for given momentum range
1138 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
1139 //look only at electron nsigma
1140 if(!fElectronOnlyRejection) {
1141 //Make the decision based on the n-sigma of electrons
1142 if(nSigma < fElectronRejectionNSigma) continue;
1145 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1146 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1147 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1149 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1150 if(nSigma < fElectronRejectionNSigma
1151 && nSigmaPions > fElectronRejectionNSigma
1152 && nSigmaKaons > fElectronRejectionNSigma
1153 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1156 }//electron rejection
1158 gRefMultiplicityTPC += 1.0;
1161 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
1162 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
1163 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
1166 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1167 else if(iChannel >= 32)
1168 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1171 //Equalization of gain
1172 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1174 gRefMultiplicityVZEROA /= gFactorA;
1175 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1177 gRefMultiplicityVZEROC /= gFactorC;
1178 if((gFactorA != 0)&&(gFactorC != 0))
1179 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1182 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1184 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1186 if(fMultiplicityEstimator == "TPC")
1187 gRefMultiplicity = gRefMultiplicityTPC;
1188 else if(fMultiplicityEstimator == "V0M")
1189 gRefMultiplicity = gRefMultiplicityVZERO;
1190 else if(fMultiplicityEstimator == "V0A")
1191 gRefMultiplicity = gRefMultiplicityVZEROA;
1192 else if(fMultiplicityEstimator == "V0C")
1193 gRefMultiplicity = gRefMultiplicityVZEROC;
1195 return gRefMultiplicity;
1198 //________________________________________________________________________
1199 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1200 // Get the event plane
1202 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1204 Float_t gVZEROEventPlane = -10.;
1205 Float_t gReactionPlane = -10.;
1206 Double_t qxTot = 0.0, qyTot = 0.0;
1208 //MC: from reaction plane
1209 if(gAnalysisLevel == "MC"){
1211 AliError("mcEvent not available");
1215 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1217 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1219 gReactionPlane = headerH->ReactionPlaneAngle();
1220 //gReactionPlane *= TMath::RadToDeg();
1225 // AOD,ESD,ESDMC: from VZERO Event Plane
1228 AliEventplane *ep = event->GetEventplane();
1230 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1231 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1232 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1233 gReactionPlane = gVZEROEventPlane;
1237 return gReactionPlane;
1240 //________________________________________________________________________
1241 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
1245 Double_t gCentrality) {
1246 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
1248 Double_t correction = 1.;
1249 Int_t gCentralityInt = -1;
1251 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1252 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1258 // centrality not in array --> no correction
1259 if(gCentralityInt < 0){
1264 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1266 if(fHistCorrectionPlus[gCentralityInt]){
1268 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1269 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
1272 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1273 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
1279 }//centrality in array
1281 if (correction == 0.) {
1282 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
1289 //________________________________________________________________________
1290 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
1291 // Returns TObjArray with tracks after all track cuts (only for AOD!)
1292 // Fills QA histograms
1294 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1296 //output TObjArray holding all good tracks
1297 TObjArray* tracksAccepted = new TObjArray;
1298 tracksAccepted->SetOwner(kTRUE);
1307 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1308 // Loop over tracks in event
1310 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1311 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1313 AliError(Form("Could not receive track %d", iTracks));
1319 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
1320 // take only TPC only tracks
1321 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
1322 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1323 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1326 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1329 // additional check on kPrimary flag
1330 if(fCheckPrimaryFlagAOD){
1331 if(aodTrack->GetType() != AliAODTrack::kPrimary)
1336 vCharge = aodTrack->Charge();
1337 vEta = aodTrack->Eta();
1339 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1340 vPt = aodTrack->Pt();
1342 //===========================PID (so far only for electron rejection)===============================//
1343 if(fElectronRejection) {
1345 // get the electron nsigma
1346 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1348 //Fill QA before the PID
1349 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1350 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1351 //end of QA-before pid
1353 // check only for given momentum range
1354 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1356 //look only at electron nsigma
1357 if(!fElectronOnlyRejection){
1359 //Make the decision based on the n-sigma of electrons
1360 if(nSigma < fElectronRejectionNSigma) continue;
1364 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1365 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1366 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1368 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1369 if(nSigma < fElectronRejectionNSigma
1370 && nSigmaPions > fElectronRejectionNSigma
1371 && nSigmaKaons > fElectronRejectionNSigma
1372 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1376 //Fill QA after the PID
1377 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1378 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1381 //===========================end of PID (so far only for electron rejection)===============================//
1383 //+++++++++++++++++++++++++++++//
1384 //===========================PID===============================//
1386 Double_t prob[AliPID::kSPECIES]={0.};
1387 Double_t probTPC[AliPID::kSPECIES]={0.};
1388 Double_t probTOF[AliPID::kSPECIES]={0.};
1389 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1391 Double_t nSigma = 0.;
1392 Double_t nSigmaTPC = 0.; //++++
1393 Double_t nSigmaTOF = 0.; //++++
1394 Double_t nSigmaTPCTOF = 0.; //++++
1395 UInt_t detUsedTPC = 0;
1396 UInt_t detUsedTOF = 0;
1397 UInt_t detUsedTPCTOF = 0;
1399 nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1400 nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1401 nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
1403 //Decide what detector configuration we want to use
1404 switch(fPidDetectorConfig) {
1406 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1408 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1409 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1410 prob[iSpecies] = probTPC[iSpecies];
1413 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1415 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1416 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1417 prob[iSpecies] = probTOF[iSpecies];
1420 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1421 nSigma = nSigmaTPCTOF;
1422 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
1423 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1424 prob[iSpecies] = probTPCTOF[iSpecies];
1428 }//end switch: define detector mask
1430 //Filling the PID QA
1431 Double_t tofTime = -999., length = 999., tof = -999.;
1432 Double_t c = TMath::C()*1.E-9;// m/ns
1433 Double_t beta = -999.;
1434 if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
1435 (aodTrack->IsOn(AliAODTrack::kTIME)) ) {
1436 tofTime = aodTrack->GetTOFsignal();//in ps
1437 length = aodTrack->GetIntegratedLength();
1438 tof = tofTime*1E-3; // ns
1441 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1445 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1449 length = length*0.01; // in meters
1453 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1454 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1455 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1458 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1459 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1460 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
1462 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1465 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1466 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
1467 Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
1469 //end of QA-before pid
1471 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1472 //Make the decision based on the n-sigma
1474 if(nSigma > fPIDNSigma) continue;
1476 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
1477 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
1478 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
1480 //Make the decision based on the bayesian
1481 else if(fUsePIDPropabilities) {
1482 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1483 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
1485 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1486 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1487 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1491 //Fill QA after the PID
1492 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1493 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1494 fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1497 //===========================PID===============================//
1498 //+++++++++++++++++++++++++++++//
1501 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1502 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1505 // Kinematics cuts from ESD track cuts
1506 if( vPt < fPtMin || vPt > fPtMax) continue;
1507 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1509 // Extra DCA cuts (for systematic studies [!= -1])
1510 if( fDCAxyCut != -1 && fDCAzCut != -1){
1511 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1516 // Extra TPC cuts (for systematic studies [!= -1])
1517 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1520 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1524 // fill QA histograms
1525 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1526 fHistDCA->Fill(dcaZ,dcaXY);
1527 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1528 fHistPt->Fill(vPt,gCentrality);
1529 fHistEta->Fill(vEta,gCentrality);
1530 fHistRapidity->Fill(vY,gCentrality);
1531 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1532 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1533 fHistPhi->Fill(vPhi,gCentrality);
1534 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1535 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1537 //=======================================correction
1538 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1539 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1541 // add the track to the TObjArray
1542 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1546 //==============================================================================================================
1547 else if(gAnalysisLevel == "MCAOD") {
1549 AliMCEvent* mcEvent = MCEvent();
1551 AliError("ERROR: Could not retrieve MC event");
1555 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1556 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
1558 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1562 if(!aodTrack->IsPhysicalPrimary()) continue;
1564 vCharge = aodTrack->Charge();
1565 vEta = aodTrack->Eta();
1567 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1568 vPt = aodTrack->Pt();
1570 // Kinematics cuts from ESD track cuts
1571 if( vPt < fPtMin || vPt > fPtMax) continue;
1572 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1574 // Remove neutral tracks
1575 if( vCharge == 0 ) continue;
1577 //Exclude resonances
1578 if(fExcludeResonancesInMC) {
1580 Bool_t kExcludeParticle = kFALSE;
1581 Int_t gMotherIndex = aodTrack->GetMother();
1582 if(gMotherIndex != -1) {
1583 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1585 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1586 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1587 //if(pdgCodeOfMother == 113) {
1588 if(pdgCodeOfMother == 113 // rho0
1589 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1590 // || pdgCodeOfMother == 221 // eta
1591 // || pdgCodeOfMother == 331 // eta'
1592 // || pdgCodeOfMother == 223 // omega
1593 // || pdgCodeOfMother == 333 // phi
1594 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1595 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1596 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1597 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1598 || pdgCodeOfMother == 111 // pi0 Dalitz
1599 || pdgCodeOfMother == 22 // photon
1601 kExcludeParticle = kTRUE;
1606 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1607 if(kExcludeParticle) continue;
1610 //Exclude electrons with PDG
1611 if(fExcludeElectronsInMC) {
1613 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1617 // fill QA histograms
1618 fHistPt->Fill(vPt,gCentrality);
1619 fHistEta->Fill(vEta,gCentrality);
1620 fHistRapidity->Fill(vY,gCentrality);
1621 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1622 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1623 fHistPhi->Fill(vPhi,gCentrality);
1624 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1625 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1627 //=======================================correction
1628 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1629 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1631 // add the track to the TObjArray
1632 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1636 //==============================================================================================================
1638 //==============================================================================================================
1639 else if(gAnalysisLevel == "MCAODrec") {
1641 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1643 printf("ERROR: fAOD not available\n");
1647 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
1649 AliError("No array of MC particles found !!!");
1652 AliMCEvent* mcEvent = MCEvent();
1654 AliError("ERROR: Could not retrieve MC event");
1655 return tracksAccepted;
1658 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1659 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1661 AliError(Form("Could not receive track %d", iTracks));
1665 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1666 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1668 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1670 vCharge = aodTrack->Charge();
1671 vEta = aodTrack->Eta();
1673 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1674 vPt = aodTrack->Pt();
1676 //===========================use MC information for Kinematics===============================//
1677 if(fUseMCforKinematics){
1679 Int_t label = TMath::Abs(aodTrack->GetLabel());
1680 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1683 vCharge = AODmcTrack->Charge();
1684 vEta = AODmcTrack->Eta();
1685 vY = AODmcTrack->Y();
1686 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
1687 vPt = AODmcTrack->Pt();
1690 AliDebug(1, "no MC particle for this track");
1694 //===========================end of use MC information for Kinematics========================//
1697 //===========================PID (so far only for electron rejection)===============================//
1698 if(fElectronRejection) {
1700 // get the electron nsigma
1701 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1703 //Fill QA before the PID
1704 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1705 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1706 //end of QA-before pid
1708 // check only for given momentum range
1709 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1711 //look only at electron nsigma
1712 if(!fElectronOnlyRejection){
1714 //Make the decision based on the n-sigma of electrons
1715 if(nSigma < fElectronRejectionNSigma) continue;
1719 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1720 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1721 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1723 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1724 if(nSigma < fElectronRejectionNSigma
1725 && nSigmaPions > fElectronRejectionNSigma
1726 && nSigmaKaons > fElectronRejectionNSigma
1727 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1731 //Fill QA after the PID
1732 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1733 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1736 //===========================end of PID (so far only for electron rejection)===============================//
1738 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1739 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1741 // Kinematics cuts from ESD track cuts
1742 if( vPt < fPtMin || vPt > fPtMax) continue;
1743 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1745 // Extra DCA cuts (for systematic studies [!= -1])
1746 if( fDCAxyCut != -1 && fDCAzCut != -1){
1747 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1752 // Extra TPC cuts (for systematic studies [!= -1])
1753 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1756 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1760 //Exclude resonances
1761 if(fExcludeResonancesInMC) {
1763 Bool_t kExcludeParticle = kFALSE;
1765 Int_t label = TMath::Abs(aodTrack->GetLabel());
1766 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1769 //if (AODmcTrack->IsPhysicalPrimary()){
1771 Int_t gMotherIndex = AODmcTrack->GetMother();
1772 if(gMotherIndex != -1) {
1773 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1775 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1776 if(pdgCodeOfMother == 113 // rho0
1777 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1778 // || pdgCodeOfMother == 221 // eta
1779 // || pdgCodeOfMother == 331 // eta'
1780 // || pdgCodeOfMother == 223 // omega
1781 // || pdgCodeOfMother == 333 // phi
1782 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1783 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1784 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1785 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1786 || pdgCodeOfMother == 111 // pi0 Dalitz
1787 || pdgCodeOfMother == 22 // photon
1789 kExcludeParticle = kTRUE;
1794 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1795 if(kExcludeParticle) continue;
1798 //Exclude electrons with PDG
1799 if(fExcludeElectronsInMC) {
1801 Int_t label = TMath::Abs(aodTrack->GetLabel());
1802 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1805 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1809 // fill QA histograms
1810 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1811 fHistDCA->Fill(dcaZ,dcaXY);
1812 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1813 fHistPt->Fill(vPt,gCentrality);
1814 fHistEta->Fill(vEta,gCentrality);
1815 fHistRapidity->Fill(vY,gCentrality);
1816 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1817 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1818 fHistPhi->Fill(vPhi,gCentrality);
1819 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1820 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1822 //=======================================correction
1823 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1824 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1826 // add the track to the TObjArray
1827 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1830 //==============================================================================================================
1832 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1834 AliESDtrack *trackTPC = NULL;
1835 AliMCParticle *mcTrack = NULL;
1837 AliMCEvent* mcEvent = NULL;
1839 // for MC ESDs use also MC information
1840 if(gAnalysisLevel == "MCESD"){
1841 mcEvent = MCEvent();
1843 AliError("mcEvent not available");
1844 return tracksAccepted;
1848 // Loop over tracks in event
1849 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1850 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
1852 AliError(Form("Could not receive track %d", iTracks));
1856 // for MC ESDs use also MC information --> MC track not used further???
1857 if(gAnalysisLevel == "MCESD"){
1858 Int_t label = TMath::Abs(track->GetLabel());
1859 if(label > mcEvent->GetNumberOfTracks()) continue;
1860 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1862 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1863 if(!mcTrack) continue;
1866 // take only TPC only tracks
1867 trackTPC = new AliESDtrack();
1868 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1872 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1874 // fill QA histograms
1877 trackTPC->GetImpactParameters(b,bCov);
1878 if (bCov[0]<=0 || bCov[2]<=0) {
1879 AliDebug(1, "Estimated b resolution lower or equal zero!");
1880 bCov[0]=0; bCov[2]=0;
1883 Int_t nClustersTPC = -1;
1884 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1885 //nClustersTPC = track->GetTPCclusters(0); // global track
1886 Float_t chi2PerClusterTPC = -1;
1887 if (nClustersTPC!=0) {
1888 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1889 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1892 //===========================PID===============================//
1894 Double_t prob[AliPID::kSPECIES]={0.};
1895 Double_t probTPC[AliPID::kSPECIES]={0.};
1896 Double_t probTOF[AliPID::kSPECIES]={0.};
1897 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1899 Double_t nSigma = 0.;
1900 UInt_t detUsedTPC = 0;
1901 UInt_t detUsedTOF = 0;
1902 UInt_t detUsedTPCTOF = 0;
1904 //Decide what detector configuration we want to use
1905 switch(fPidDetectorConfig) {
1907 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1908 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
1909 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
1910 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1911 prob[iSpecies] = probTPC[iSpecies];
1914 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1915 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
1916 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
1917 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1918 prob[iSpecies] = probTOF[iSpecies];
1921 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1922 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
1923 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1924 prob[iSpecies] = probTPCTOF[iSpecies];
1928 }//end switch: define detector mask
1930 //Filling the PID QA
1931 Double_t tofTime = -999., length = 999., tof = -999.;
1932 Double_t c = TMath::C()*1.E-9;// m/ns
1933 Double_t beta = -999.;
1934 Double_t nSigmaTOFForParticleOfInterest = -999.;
1935 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
1936 (track->IsOn(AliESDtrack::kTIME)) ) {
1937 tofTime = track->GetTOFsignal();//in ps
1938 length = track->GetIntegratedLength();
1939 tof = tofTime*1E-3; // ns
1942 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1946 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1950 length = length*0.01; // in meters
1954 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
1955 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
1956 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
1957 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
1961 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
1962 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
1963 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
1964 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
1965 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
1966 //end of QA-before pid
1968 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1969 //Make the decision based on the n-sigma
1971 if(nSigma > fPIDNSigma) continue;}
1973 //Make the decision based on the bayesian
1974 else if(fUsePIDPropabilities) {
1975 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1976 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
1979 //Fill QA after the PID
1980 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
1981 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
1982 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
1984 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
1985 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
1986 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
1987 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
1990 //===========================PID===============================//
1991 vCharge = trackTPC->Charge();
1993 vEta = trackTPC->Eta();
1994 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
1995 vPt = trackTPC->Pt();
1996 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
1997 fHistDCA->Fill(b[1],b[0]);
1998 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
1999 fHistPt->Fill(vPt,gCentrality);
2000 fHistEta->Fill(vEta,gCentrality);
2001 fHistPhi->Fill(vPhi,gCentrality);
2002 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2003 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2004 fHistRapidity->Fill(vY,gCentrality);
2005 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2006 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2008 //=======================================correction
2009 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2010 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2012 // add the track to the TObjArray
2013 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2019 else if(gAnalysisLevel == "MC"){
2021 AliError("mcEvent not available");
2025 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2027 // Loop over tracks in event
2028 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
2029 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
2031 AliError(Form("Could not receive particle %d", iTracks));
2035 //exclude non stable particles
2036 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2038 vCharge = track->Charge();
2039 vEta = track->Eta();
2043 if( vPt < fPtMin || vPt > fPtMax)
2046 if( vEta < fEtaMin || vEta > fEtaMax) continue;
2049 if( vY < fEtaMin || vY > fEtaMax) continue;
2052 // Remove neutral tracks
2053 if( vCharge == 0 ) continue;
2055 //analyze one set of particles
2057 TParticle *particle = track->Particle();
2058 if(!particle) continue;
2060 Int_t gPdgCode = particle->GetPdgCode();
2061 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
2065 //Use the acceptance parameterization
2066 if(fAcceptanceParameterization) {
2067 Double_t gRandomNumber = gRandom->Rndm();
2068 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
2072 //Exclude resonances
2073 if(fExcludeResonancesInMC) {
2074 TParticle *particle = track->Particle();
2075 if(!particle) continue;
2077 Bool_t kExcludeParticle = kFALSE;
2078 Int_t gMotherIndex = particle->GetFirstMother();
2079 if(gMotherIndex != -1) {
2080 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2082 TParticle *motherParticle = motherTrack->Particle();
2083 if(motherParticle) {
2084 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2085 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2086 if(pdgCodeOfMother == 113 // rho0
2087 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2088 // || pdgCodeOfMother == 221 // eta
2089 // || pdgCodeOfMother == 331 // eta'
2090 // || pdgCodeOfMother == 223 // omega
2091 // || pdgCodeOfMother == 333 // phi
2092 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
2093 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
2094 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
2095 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2096 || pdgCodeOfMother == 111 // pi0 Dalitz
2098 kExcludeParticle = kTRUE;
2104 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2105 if(kExcludeParticle) continue;
2108 //Exclude electrons with PDG
2109 if(fExcludeElectronsInMC) {
2111 TParticle *particle = track->Particle();
2114 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2118 vPhi = track->Phi();
2119 //Printf("phi (before): %lf",vPhi);
2121 fHistPt->Fill(vPt,gCentrality);
2122 fHistEta->Fill(vEta,gCentrality);
2123 fHistPhi->Fill(vPhi,gCentrality);
2124 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2125 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2126 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2127 fHistRapidity->Fill(vY,gCentrality);
2128 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2129 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2130 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2131 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2134 if(fUseFlowAfterBurner) {
2135 Double_t precisionPhi = 0.001;
2136 Int_t maxNumberOfIterations = 100;
2138 Double_t phi0 = vPhi;
2139 Double_t gV2 = fDifferentialV2->Eval(vPt);
2141 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2142 Double_t phiprev = vPhi;
2143 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2144 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2146 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2148 //Printf("phi (after): %lf\n",vPhi);
2149 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2150 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2151 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2153 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2154 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2155 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2159 //vPhi *= TMath::RadToDeg();
2161 //=======================================correction
2162 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2163 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2165 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2170 return tracksAccepted;
2173 //________________________________________________________________________
2174 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2175 // Clones TObjArray and returns it with tracks after shuffling the charges
2177 TObjArray* tracksShuffled = new TObjArray;
2178 tracksShuffled->SetOwner(kTRUE);
2180 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
2182 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2184 AliVParticle* track = (AliVParticle*) tracks->At(i);
2185 chargeVector->push_back(track->Charge());
2188 random_shuffle(chargeVector->begin(), chargeVector->end());
2190 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2191 AliVParticle* track = (AliVParticle*) tracks->At(i);
2192 //==============================correction
2193 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2194 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2195 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2198 delete chargeVector;
2200 return tracksShuffled;
2203 //________________________________________________________________________
2204 void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2205 const char* lhcPeriod) {
2206 //Function to setup the VZERO gain equalization
2207 //============Get the equilization map============//
2208 TFile *calibrationFile = TFile::Open(filename);
2209 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2210 Printf("No calibration file found!!!");
2214 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2216 Printf("Calibration TList not found!!!");
2220 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2221 if(!fHistVZEROAGainEqualizationMap) {
2222 Printf("VZERO-A calibration object not found!!!");
2225 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2226 if(!fHistVZEROCGainEqualizationMap) {
2227 Printf("VZERO-C calibration object not found!!!");
2231 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2232 if(!fHistVZEROChannelGainEqualizationMap) {
2233 Printf("VZERO channel calibration object not found!!!");
2238 //________________________________________________________________________
2239 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
2242 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2244 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2245 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2246 if(gRunNumber == run)
2247 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2253 //________________________________________________________________________
2254 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
2257 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2259 TString gVZEROSide = side;
2260 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2261 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2262 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2263 if(gRunNumber == run) {
2264 if(gVZEROSide == "A")
2265 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2266 else if(gVZEROSide == "C")
2267 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2274 //____________________________________________________________________
2275 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2277 // copied from AliAnalysisTaskPhiCorrelations
2279 // rejects "randomly" events such that the centrality gets flat
2280 // uses fCentralityWeights histogram
2282 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2284 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2285 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2287 Bool_t result = kFALSE;
2288 if (centralityDigits < weight)
2291 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2296 //________________________________________________________________________
2297 void AliAnalysisTaskBFPsi::FinishTaskOutput(){
2301 AliError("fBalance not available");
2305 if (!fShuffledBalance) {
2306 AliError("fShuffledBalance not available");
2313 //________________________________________________________________________
2314 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2315 // Draw result to the screen
2316 // Called once at the end of the query
2318 // not implemented ...