exclude all but kaons in MC mode exclude resonances (for testing only)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskBFPsi.cxx
1 #include "TChain.h"
2 #include "TList.h"
3 #include "TCanvas.h"
4 #include "TLorentzVector.h"
5 #include "TGraphErrors.h"
6 #include "TH1F.h"
7 #include "TH2F.h"
8 #include "TH3F.h" 
9 #include "TH2D.h"                  
10 #include "TH3D.h"
11 #include "TArrayF.h"
12 #include "TF1.h"
13 #include "TRandom.h"
14 #include "TROOT.h"
15
16 #include "AliAnalysisTaskSE.h"
17 #include "AliAnalysisManager.h"
18
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"
30 #include "AliStack.h"
31 #include "AliESDtrackCuts.h"
32 #include "AliEventplane.h"
33 #include "AliTHn.h"    
34 #include "AliLog.h"
35 #include "AliAnalysisUtils.h"
36
37 #include "AliEventPoolManager.h"           
38
39 #include "AliPID.h"                
40 #include "AliPIDResponse.h"        
41 #include "AliPIDCombined.h"        
42
43 #include "AliAnalysisTaskBFPsi.h"
44 #include "AliBalancePsi.h"
45 #include "AliAnalysisTaskTriggeredBF.h"
46
47
48 // Analysis task for the BF vs Psi code
49 // Authors: Panos.Christakoglou@nikhef.nl
50
51 using std::cout;
52 using std::endl;
53
54 ClassImp(AliAnalysisTaskBFPsi)
55
56 //________________________________________________________________________
57 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) 
58 : AliAnalysisTaskSE(name),
59   fDebugLevel(kFALSE),
60   fArrayMC(0),
61   fBalance(0),
62   fRunShuffling(kFALSE),
63   fShuffledBalance(0),
64   fRunMixing(kFALSE),
65   fRunMixingEventPlane(kFALSE),
66   fMixingTracks(50000),
67   fMixedBalance(0),
68   fPoolMgr(0),
69   fList(0),
70   fListBF(0),
71   fListBFS(0),
72   fListBFM(0),
73   fHistListPIDQA(0),
74   fHistEventStats(0),
75   fHistCentStats(0),
76   fHistCentStatsUsed(0),
77   fHistTriggerStats(0),
78   fHistTrackStats(0),
79   fHistVx(0),
80   fHistVy(0),
81   fHistVz(0),
82   fHistTPCvsVZEROMultiplicity(0),
83   fHistVZEROSignal(0),
84   fHistEventPlane(0),
85   fHistClus(0),
86   fHistDCA(0),
87   fHistChi2(0),
88   fHistPt(0),
89   fHistEta(0),
90   fHistRapidity(0),
91   fHistPhi(0),
92   fHistEtaPhiPos(0),             
93   fHistEtaPhiNeg(0), 
94   fHistPhiBefore(0),
95   fHistPhiAfter(0),
96   fHistPhiPos(0),
97   fHistPhiNeg(0),
98   fHistV0M(0),
99   fHistRefTracks(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),
126   fPIDResponse(0x0),
127   fPIDCombined(0x0),
128   fParticleOfInterest(kPion),
129   fPidDetectorConfig(kTPCTOF),
130   fUsePID(kFALSE),
131   fUsePIDnSigma(kTRUE),
132   fUsePIDPropabilities(kFALSE), 
133   fPIDNSigma(3.),
134   fMinAcceptedPIDProbability(0.8),
135   fElectronRejection(kFALSE),
136   fElectronOnlyRejection(kFALSE),
137   fElectronRejectionNSigma(-1.),
138   fElectronRejectionMinPt(0.),
139   fElectronRejectionMaxPt(1000.),
140   fESDtrackCuts(0),
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),
158   fVxMax(0.3),
159   fVyMax(0.3),
160   fVzMax(10.),
161   fnAODtrackCutBit(128),
162   fPtMin(0.3),
163   fPtMax(1.5),
164   fEtaMin(-0.8),
165   fEtaMax(0.8),
166   fPhiMin(0.),
167   fPhiMax(360.),
168   fDCAxyCut(-1),
169   fDCAzCut(-1),
170   fTPCchi2Cut(-1),
171   fNClustersTPCCut(-1),
172   fTPCsharedCut(-1),
173   fAcceptanceParameterization(0),
174   fDifferentialV2(0),
175   fUseFlowAfterBurner(kFALSE),
176   fExcludeResonancesInMC(kFALSE),
177   fExcludeElectronsInMC(kFALSE),
178   fUseMCPdgCode(kFALSE),
179   fPDGCodeToBeAnalyzed(-1),
180   fEventClass("EventPlane"), 
181   fCustomBinning(""),
182   fHistVZEROAGainEqualizationMap(0),
183   fHistVZEROCGainEqualizationMap(0),
184   fHistVZEROChannelGainEqualizationMap(0) {
185   // Constructor
186   // Define input and output slots here
187   // Input slot #0 works with a TChain
188
189   //======================================================correction
190   for (Int_t i=0; i<kCENTRALITY; i++){
191     fHistCorrectionPlus[i] = NULL; 
192     fHistCorrectionMinus[i] = NULL; 
193     fCentralityArrayForCorrections[i] = -1.;
194   }
195   //=====================================================correction
196
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());
204 }
205
206 //________________________________________________________________________
207 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
208
209   // delete fBalance; 
210   // delete fShuffledBalance; 
211   // delete fList;
212   // delete fListBF; 
213   // delete fListBFS;
214
215   // delete fHistEventStats; 
216   // delete fHistTrackStats; 
217   // delete fHistVx; 
218   // delete fHistVy; 
219   // delete fHistVz; 
220
221   // delete fHistClus;
222   // delete fHistDCA;
223   // delete fHistChi2;
224   // delete fHistPt;
225   // delete fHistEta;
226   // delete fHistPhi;
227   // delete fHistEtaPhiPos;                      
228   // delete fHistEtaPhiNeg;
229   // delete fHistV0M;
230 }
231
232 //________________________________________________________________________
233 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
234   // Create histograms
235   // Called once
236
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);
241
242   if(!fBalance) {
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.);
248   }
249   if(fRunShuffling) {
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.);
256     }
257   }
258   if(fRunMixing) {
259     if(!fMixedBalance) {
260       fMixedBalance = new AliBalancePsi();
261       fMixedBalance->SetAnalysisLevel("ESD");
262       fMixedBalance->SetEventClass(fEventClass);
263     }
264   }
265
266   //QA list
267   fList = new TList();
268   fList->SetName("listQA");
269   fList->SetOwner();
270
271   //Balance Function list
272   fListBF = new TList();
273   fListBF->SetName("listBF");
274   fListBF->SetOwner();
275
276   if(fRunShuffling) {
277     fListBFS = new TList();
278     fListBFS->SetName("listBFShuffled");
279     fListBFS->SetOwner();
280   }
281
282   if(fRunMixing) {
283     fListBFM = new TList();
284     fListBFM->SetName("listTriggeredBFMixed");
285     fListBFM->SetOwner();
286   }
287
288   //PID QA list
289   if(fUsePID || fElectronRejection) {
290     fHistListPIDQA = new TList();
291     fHistListPIDQA->SetName("listQAPID");
292     fHistListPIDQA->SetOwner();
293   }
294
295   //Event stats.
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);
304
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());
312   }
313   fList->Add(fHistCentStats);
314
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);
318
319   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
320   fList->Add(fHistTriggerStats);
321
322   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
323   fList->Add(fHistTrackStats);
324
325   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
326   fList->Add(fHistNumberOfAcceptedTracks);
327
328   fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
329   fList->Add(fHistMultiplicity);
330
331   // Vertex distributions
332   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
333   fList->Add(fHistVx);
334   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
335   fList->Add(fHistVy);
336   fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
337   fList->Add(fHistVz);
338
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.)");
345   else 
346     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
347   fList->Add(fHistTPCvsVZEROMultiplicity);
348
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);
351
352   //Event plane
353   fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
354   fList->Add(fHistEventPlane);
355
356   // QA histograms
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);
364   fList->Add(fHistPt);
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);
390
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();
397   }
398
399   if(fRunShuffling) {
400     if(!fShuffledBalance->GetHistNp()) {
401       AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
402       fShuffledBalance->SetCustomBinning(fCustomBinning);
403       fShuffledBalance->InitHistograms();
404     }
405   }
406
407   if(fRunMixing) {
408     if(!fMixedBalance->GetHistNp()) {
409       AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
410       fMixedBalance->SetCustomBinning(fCustomBinning);
411       fMixedBalance->InitHistograms();
412     }
413   }
414
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());
427
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());
435
436   if(fRunShuffling) {
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());
443   }  
444
445   if(fRunMixing) {
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());
452   }
453   //}
454
455
456   // Event Mixing
457   if(fRunMixing){
458     Int_t trackDepth = fMixingTracks; 
459     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
460     
461     // centrality bins
462     Double_t* centbins = NULL;
463     Int_t nCentralityBins;
464     if(fBalance->IsUseVertexBinning()){
465       centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
466     }
467     else{
468       centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
469     }
470     
471     // multiplicity bins
472     Double_t* multbins = NULL;
473     Int_t nMultiplicityBins;
474     multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
475     
476     // Zvtx bins
477     Double_t* vtxbins = NULL; 
478     Int_t nVertexBins;
479     if(fBalance->IsUseVertexBinning()){
480       vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
481     }
482     else{
483       vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
484     }
485
486     // Event plane angle (Psi) bins
487     Double_t* psibins = NULL;
488     Int_t nPsiBins; 
489     psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
490
491   
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);
497         }
498       }
499       else{
500         if(centbins && vtxbins && psibins){
501           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
502         }
503       }
504     }
505     else{
506       if(fEventClass=="Multiplicity"){
507         if(multbins && vtxbins){
508           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
509         }
510       }
511       else{
512         if(centbins && vtxbins){
513           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
514         }
515       }
516     }
517     
518     if(centbins) delete [] centbins; 
519     if(multbins) delete [] multbins; 
520     if(vtxbins)  delete [] vtxbins; 
521     if(psibins)  delete [] psibins; 
522
523     // set minimum values for track depth, fraction, and number of events
524     fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
525     
526     // check pool manager
527     if(!fPoolMgr){
528       AliError("Event Mixing required, but Pool Manager not initialized...");
529       return;
530     }
531   }
532   
533   if(fESDtrackCuts) fList->Add(fESDtrackCuts);
534
535   //====================PID========================//
536   if(fUsePID) {
537     fPIDCombined = new AliPIDCombined();
538     fPIDCombined->SetDefaultTPCPriors();
539
540     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); 
541     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
542     
543     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); 
544     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); 
545     
546     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); 
547     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); 
548     
549     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); 
550     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
551
552     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); 
553     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
554     
555     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, -25, 25); 
556     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
557     
558     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, -25, 25); 
559     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); 
560
561     fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); 
562     fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID);
563     
564     fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, -25, 25); 
565     fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID);
566     
567     //+++++++++++++++++//
568     //p array
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};
571     //nSigma Array
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);
577     }
578  
579     fHistNSigmaTPCTOFPbefPID = new TH3D ("fHistNSigmaTPCTOFPbefPID","fHistNSigmaTPCTOFPbefPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP); 
580     fHistListPIDQA->Add(fHistNSigmaTPCTOFPbefPID); 
581     //++++++++++++++//
582
583     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); 
584     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
585     
586     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); 
587     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); 
588     
589     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); 
590     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
591   
592     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); 
593     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); 
594     
595     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); 
596     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
597
598     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, -25, 25); 
599     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
600     
601     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, -25, 25); 
602     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
603
604     fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); 
605     fHistListPIDQA->Add(fHistBetaVsdEdXafterPID);
606
607     fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, -25, 25); 
608     fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID);
609
610     fHistNSigmaTPCTOFPafterPID = new TH3D ("fHistNSigmaTPCTOFPafterPID","fHistNSigmaTPCTOFPafterPID;#sigma_{TPC};#sigma_{TOF};p_{T} (GeV/c)", nSigmaBins, nArrayS, nSigmaBins, nArrayS, pBins,nArrayP); 
611     fHistListPIDQA->Add(fHistNSigmaTPCTOFPafterPID); //++++++++++++++
612   }
613
614   // for electron rejection only TPC nsigma histograms
615   if(fElectronRejection) {
616  
617     fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000); 
618     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
619     
620     fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500); 
621     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
622     
623     fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000); 
624     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
625
626     fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500); 
627     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron); 
628   }
629   //====================PID========================//
630
631   // Post output data.
632   PostData(1, fList);
633   PostData(2, fListBF);
634   if(fRunShuffling) PostData(3, fListBFS);
635   if(fRunMixing) PostData(4, fListBFM);
636   if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA);       //PID
637
638   AliInfo("Finished setting up the Output");
639
640   TH1::AddDirectory(oldStatus);
641
642 }
643
644
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];
653
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()));
657     return;
658   }
659
660   //Open the input file
661   TFile *f = TFile::Open(filename);
662   if(!f->IsOpen()) {
663     AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
664     return;
665   }
666     
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()));
675       return;
676     }
677     
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()));
683       return; 
684     }
685   }//loop over centralities: ONLY the PbPb case is covered
686 }
687
688 //________________________________________________________________________
689 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
690   // Main loop
691   // Called for each event
692
693   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
694   Int_t gNumberOfAcceptedTracks = 0;
695   Double_t lMultiplicityVar     = -999.; //-1
696   Double_t gReactionPlane       = -1.; 
697   Float_t bSign = 0.;
698   
699   // get the event (for generator level: MCEvent())
700   AliVEvent* eventMain = NULL;
701   if(gAnalysisLevel == "MC") {
702     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); 
703   }
704   else{
705     eventMain = dynamic_cast<AliVEvent*>(InputEvent());     
706     // for HBT like cuts need magnetic field sign
707     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
708   }
709   if(!eventMain) {
710     AliError("eventMain not available");
711     return;
712   }
713   
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");
718   }
719  
720   // check event cuts and fill event histograms
721   if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ 
722     return;
723   }
724   // get the reaction plane
725   if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") {
726     gReactionPlane = GetEventPlane(eventMain);
727     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
728     if(gReactionPlane < 0){
729       return;
730     }
731   }
732   
733   // get the accepted tracks in main event
734   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
735   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
736
737   //multiplicity cut (used in pp)
738   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
739
740   // store charges of all accepted tracks,shuffle and reassign(two extra loops)
741   TObjArray* tracksShuffled = NULL;
742   if(fRunShuffling){
743     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
744   }
745   
746   // Event mixing 
747   if (fRunMixing)
748     {
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`
756       //    statistics.
757       
758       // 2. Collect the whole pool's content of tracks into one TObjArray
759       //    (bgTracks), which is effectively a single background super-event.
760       
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.
764       
765       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
766       
767       if (!pool){
768         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
769       }
770       else{
771         
772         //pool->SetDebug(1);
773
774         if (pool->IsReady()){ 
775           
776           
777           Int_t nMix = pool->GetCurrentNEvents();
778           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
779           
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);
784           
785           // Fill mixed-event histos here  
786           for (Int_t jMix=0; jMix<nMix; jMix++) 
787             {
788               TObjArray* tracksMixed = pool->GetEvent(jMix);
789               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
790             }
791         }
792         
793         // Update the Event pool
794         pool->UpdatePool(tracksMain);
795         //pool->PrintInfo();
796         
797       }//pool NULL check  
798     }//run mixing
799   
800   // calculate balance function
801   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
802   
803   // calculate shuffled balance function
804   if(fRunShuffling && tracksShuffled != NULL) {
805     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
806   }
807 }      
808
809 //________________________________________________________________________
810 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
811   // Checks the Event cuts
812   // Fills Event statistics histograms
813   
814   Bool_t isSelectedMain = kTRUE;
815   Float_t gRefMultiplicity = -1.;
816   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
817
818   AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
819
820   fHistEventStats->Fill(1,gRefMultiplicity); //all events
821
822   // check first event in chunk (is not needed for new reconstructions)
823   if(fCheckFirstEventInChunk){
824     AliAnalysisUtils ut;
825     if(ut.IsFirstEventInChunk(event)) 
826       return -1.;
827     fHistEventStats->Fill(6,gRefMultiplicity); 
828   }
829   // check for pile-up event
830   if(fCheckPileUp){
831     AliAnalysisUtils ut;
832     ut.SetUseMVPlpSelection(kTRUE);
833     ut.SetUseOutOfBunchPileUp(kTRUE);
834     if(ut.IsPileUpEvent(event))
835       return -1.;
836     fHistEventStats->Fill(7,gRefMultiplicity); 
837   }
838
839   // Event trigger bits
840   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
841   if(fUseOfflineTrigger)
842     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
843   
844   if(isSelectedMain) {
845     fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
846  
847     // Event Vertex MC
848     if(gAnalysisLevel == "MC") {
849       if(!event) {
850         AliError("mcEvent not available");
851         return 0x0;
852       }
853       
854       if(mcevent){
855         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
856         if(header){  
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
868
869                 // get the reference multiplicty or centrality
870                 gRefMultiplicity = GetRefMultiOrCentrality(event);
871
872                 fHistVx->Fill(gVertexArray.At(0));
873                 fHistVy->Fill(gVertexArray.At(1));
874                 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
875                 
876                 // take only events inside centrality class
877                 if(fUseCentrality) {
878                   if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
879                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
880                     return gRefMultiplicity;        
881                   }//centrality class
882                 }
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;
888                   }
889                 }//multiplicity range
890               }//Vz cut
891             }//Vy cut
892           }//Vx cut
893         }//header    
894       }//MC event object
895     }//MC
896     
897     // Event Vertex AOD, ESD, ESDMC
898     else{
899       const AliVVertex *vertex = event->GetPrimaryVertex();
900       
901       if(vertex) {
902         Double32_t fCov[6];
903         vertex->GetCovarianceMatrix(fCov);
904         if(vertex->GetNContributors() > 0) {
905           if(fCov[5] != 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
911
912                   // get the reference multiplicty or centrality
913                   gRefMultiplicity = GetRefMultiOrCentrality(event);
914                   
915                   fHistVx->Fill(vertex->GetX());
916                   fHistVy->Fill(vertex->GetY());
917                   fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
918                   
919                   // take only events inside centrality class
920                   if(fUseCentrality) {
921                     if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
922
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));
926                         return -1;
927                       }
928                       
929                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
930                       return gRefMultiplicity;          
931                     }//centrality class
932                   }
933                   // take events only within the same multiplicity class
934                   else if(fUseMultiplicity) {
935                     //if(fDebugLevel) 
936                     //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
937                     //fNumberOfAcceptedTracksMax,gRefMultiplicity);
938
939                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
940                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
941                       return gRefMultiplicity;
942                     }
943                   }//multiplicity range
944                 }//Vz cut
945               }//Vy cut
946             }//Vx cut
947           }//proper vertex resolution
948         }//proper number of contributors
949       }//vertex object valid
950     }//triggered event 
951   }//AOD,ESD,ESDMC
952   
953   // in all other cases return -1 (event not accepted)
954   return -1;
955 }
956
957
958 //________________________________________________________________________
959 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
960     // Checks the Event cuts
961     // Fills Event statistics histograms
962
963   Float_t gCentrality = -1.;
964   Double_t gMultiplicity = -1.;
965   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
966
967
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();
971     if(header){
972       gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
973
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"));
988       
989       // Centrality estimator USED   ++++++++++++++++++++++++++++++
990       fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
991       
992       // centrality QA (V0M)
993       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
994       
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));
1005
1006     }//AOD header
1007   }//AOD
1008
1009   // calculate centrality always (not only in centrality mode)
1010   else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround
1011     
1012     AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1013     if(header){
1014       gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0;
1015       
1016       // QA histogram
1017       fHistCentStatsUsed->Fill(0.,gCentrality);
1018
1019     }//AOD header
1020   }//AODnano
1021   
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());
1025
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"));
1040     
1041     // Centrality estimator USED   ++++++++++++++++++++++++++++++
1042     fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
1043     
1044     // centrality QA (V0M)
1045     fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1046   }//ESD
1047
1048   else if(gAnalysisLevel == "MC"){
1049     Double_t gImpactParameter = 0.;
1050     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1051     if(gMCEvent){
1052       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());      
1053       if(headerH){
1054         gImpactParameter = headerH->ImpactParameter();
1055         gCentrality      = gImpactParameter;
1056       }//MC header
1057     }//MC event cast
1058   }//MC
1059
1060   else{
1061     gCentrality = -1.;
1062   }
1063   
1064   // calculate reference multiplicity always (not only in multiplicity mode)
1065   if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1066     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1067     if(gESDEvent){
1068       gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1069       fHistMultiplicity->Fill(gMultiplicity);
1070     }//AliESDevent cast
1071   }//ESD mode
1072
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);
1081     }
1082     else {
1083       if(header)
1084         gMultiplicity = header->GetRefMultiplicity();
1085       if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1086     }
1087     fHistMultiplicity->Fill(gMultiplicity);
1088   }//AOD mode
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));
1095       if (!track) {
1096         AliError(Form("Could not receive particle %d", iParticle));
1097         continue;
1098       }
1099       
1100       //exclude non stable particles
1101       if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1102       
1103       //++++++++++++++++
1104       if (fMultiplicityEstimator == "V0M") {
1105         if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) 
1106           continue;}
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;
1114       }
1115       else{
1116         if(track->Pt() < fPtMin || track->Pt() > fPtMax)  continue;
1117         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;
1118       }
1119       //++++++++++++++++
1120       
1121       if(track->Charge() == 0) continue;
1122       
1123       gMultiplicity += 1;
1124     }//loop over primaries
1125     fHistMultiplicity->Fill(gMultiplicity);
1126   }//MC mode
1127   else{
1128     gMultiplicity = -1;
1129   }
1130   
1131
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;
1138   }
1139   return lReturnVal;
1140 }
1141
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.;
1148
1149   AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1150   if(!header) {
1151     Printf("ERROR: AOD header not available");
1152     return -999;
1153   }
1154   Int_t gRunNumber = header->GetRunNumber();
1155
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));
1159     if (!aodTrack) {
1160       AliError(Form("Could not receive track %d", iTracks));
1161       continue;
1162     }
1163     
1164     // AOD track cuts
1165     if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1166             
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;
1171     
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));
1176       
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;
1183         }
1184         else {
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));
1188           
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;
1194         }
1195       }
1196     }//electron rejection
1197     
1198     gRefMultiplicityTPC += 1.0;
1199   }// track loop
1200   
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));
1204     
1205     if(iChannel < 32) 
1206       gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1207     else if(iChannel >= 32) 
1208       gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1209   }//loop over PMTs
1210   
1211   //Equalization of gain
1212   Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1213   if(gFactorA != 0)
1214     gRefMultiplicityVZEROA /= gFactorA;
1215   Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1216   if(gFactorC != 0)
1217     gRefMultiplicityVZEROC /= gFactorC;
1218   if((gFactorA != 0)&&(gFactorC != 0)) 
1219     gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1220   
1221   if(fDebugLevel) 
1222     Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1223
1224   fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1225
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;
1234   
1235   return gRefMultiplicity;
1236 }
1237
1238 //________________________________________________________________________
1239 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1240   // Get the event plane
1241
1242   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1243
1244   Float_t gVZEROEventPlane    = -10.;
1245   Float_t gReactionPlane      = -10.;
1246   Double_t qxTot = 0.0, qyTot = 0.0;
1247
1248   //MC: from reaction plane
1249   if(gAnalysisLevel == "MC"){
1250     if(!event) {
1251       AliError("mcEvent not available");
1252       return 0x0;
1253     }
1254
1255     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1256     if(gMCEvent){
1257       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    
1258       if (headerH) {
1259         gReactionPlane = headerH->ReactionPlaneAngle();
1260         //gReactionPlane *= TMath::RadToDeg();
1261       }//MC header
1262     }//MC event cast
1263   }//MC
1264   
1265   // AOD,ESD,ESDMC: from VZERO Event Plane
1266   else{
1267    
1268     AliEventplane *ep = event->GetEventplane();
1269     if(ep){ 
1270       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1271       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1272       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1273       gReactionPlane = gVZEROEventPlane;
1274     }
1275   }//AOD,ESD,ESDMC
1276
1277   return gReactionPlane;
1278 }
1279
1280 //________________________________________________________________________
1281 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, 
1282                                                                 Double_t vPhi, 
1283                                                                 Double_t vPt, 
1284                                                                 Short_t vCharge, 
1285                                                                 Double_t gCentrality) {
1286   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) 
1287
1288   Double_t correction = 1.;
1289   Int_t gCentralityInt = -1;
1290
1291   for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1292     if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1293       gCentralityInt = i;
1294       break;
1295     }
1296   }  
1297
1298   // centrality not in array --> no correction
1299   if(gCentralityInt < 0){
1300     correction = 1.;
1301   }
1302   else{
1303     
1304     //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1305     
1306     if(fHistCorrectionPlus[gCentralityInt]){
1307       if (vCharge > 0) {
1308         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1309         //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);  
1310       }
1311       if (vCharge < 0) {
1312         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1313         //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt); 
1314       }
1315     }
1316     else {
1317       correction = 1.;
1318     }
1319   }//centrality in array
1320   
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)); 
1323     correction = 1.; 
1324   } 
1325   
1326   return correction;
1327 }
1328
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
1333
1334   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1335
1336   //output TObjArray holding all good tracks
1337   TObjArray* tracksAccepted = new TObjArray;
1338   tracksAccepted->SetOwner(kTRUE);
1339
1340   Short_t vCharge;
1341   Double_t vEta;
1342   Double_t vY;
1343   Double_t vPhi;
1344   Double_t vPt;
1345
1346
1347   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1348     // Loop over tracks in event
1349     
1350     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1351       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1352       if (!aodTrack) {
1353         AliError(Form("Could not receive track %d", iTracks));
1354         continue;
1355       }
1356       
1357       // AOD track cuts
1358       
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));
1364       }
1365
1366       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1367
1368
1369       // additional check on kPrimary flag
1370       if(fCheckPrimaryFlagAOD){
1371         if(aodTrack->GetType() != AliAODTrack::kPrimary)
1372           continue;
1373       }
1374
1375      
1376       vCharge = aodTrack->Charge();
1377       vEta    = aodTrack->Eta();
1378       vY      = aodTrack->Y();
1379       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1380       vPt     = aodTrack->Pt();
1381       
1382       //===========================PID (so far only for electron rejection)===============================//                
1383       if(fElectronRejection) {
1384
1385         // get the electron nsigma
1386         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1387         
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
1392         
1393         // check only for given momentum range
1394         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1395                   
1396           //look only at electron nsigma
1397           if(!fElectronOnlyRejection){
1398             
1399             //Make the decision based on the n-sigma of electrons
1400             if(nSigma < fElectronRejectionNSigma) continue;
1401           }
1402           else{
1403             
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));
1407             
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;
1413           }
1414         }
1415   
1416         //Fill QA after the PID
1417         fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1418         fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1419         
1420       }
1421       //===========================end of PID (so far only for electron rejection)===============================//
1422
1423       //+++++++++++++++++++++++++++++//
1424       //===========================PID===============================//             
1425       if(fUsePID) {
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.};
1430         
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;
1438         
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;
1444         }
1445
1446         //Decide what detector configuration we want to use
1447         switch(fPidDetectorConfig) {
1448         case kTPCpid:
1449           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1450           nSigma = nSigmaTPC;
1451           detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1452           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1453             prob[iSpecies] = probTPC[iSpecies];
1454           break;
1455         case kTOFpid:
1456           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1457           nSigma = nSigmaTOF;
1458           detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1459           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1460             prob[iSpecies] = probTOF[iSpecies];
1461           break;
1462         case kTPCTOF:
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];
1468           break;
1469         default:
1470           break;
1471         }//end switch: define detector mask
1472         
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.;
1477
1478         Float_t probMis = fPIDResponse->GetTOFMismatchProbability(aodTrack);
1479         if (probMis < 0.01) { //if u want to reduce mismatch using also TPC
1480
1481           if ((aodTrack->IsOn(AliAODTrack::kITSin)) &&  (aodTrack->IsOn(AliAODTrack::kTOFpid)) ) { //leonardo's analysis
1482             
1483             tofTime = aodTrack->GetTOFsignal();//in ps
1484             length = aodTrack->GetIntegratedLength();
1485             tof = tofTime*1E-3; // ns                 
1486             if (tof <= 0) {
1487               //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1488               continue;
1489             }
1490             if (length <= 0){
1491               //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1492               continue;
1493             }         
1494             length = length*0.01; // in meters
1495             tof = tof*c;
1496             beta = length/tof;
1497             
1498             fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1499             fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1500             fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1501           }//TOF signal 
1502           
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]);
1507           
1508           //combined TPC&TOF
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
1514           
1515           if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1516             //Make the decision based on the n-sigma
1517             if(fUsePIDnSigma) {
1518               if(nSigma > fPIDNSigma || nSigma < -fPIDNSigma) continue;  
1519               
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());  //++++++++++++++
1524             }
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;      
1529               
1530               fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1531               fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); 
1532               fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1533             }      
1534             
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);
1539           }
1540         }
1541       }
1542       //===========================PID===============================//
1543       //+++++++++++++++++++++++++++++//
1544
1545
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)
1548       
1549       
1550       // Kinematics cuts from ESD track cuts
1551       if( vPt < fPtMin || vPt > fPtMax)      continue;
1552       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1553       
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 ){
1557           continue;  // 2D cut
1558         }
1559       }
1560       
1561       // Extra TPC cuts (for systematic studies [!= -1])
1562       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1563         continue;
1564       }
1565       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1566         continue;
1567       }
1568
1569       // Extra cut on shared clusters
1570       if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1571         continue;
1572       }
1573       
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);
1586       
1587       //=======================================correction
1588       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1589       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1590       
1591       // add the track to the TObjArray
1592       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1593     }//track loop
1594   }// AOD analysis
1595
1596
1597   // nano AODs
1598   else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing)
1599     // Loop over tracks in event
1600     
1601     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1602       AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks));
1603       if (!aodTrack) {
1604         AliError(Form("Could not receive track %d", iTracks));
1605         continue;
1606       }
1607       
1608       // AOD track cuts (not needed)
1609       //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1610      
1611       vCharge = aodTrack->Charge();
1612       vEta    = aodTrack->Eta();
1613       vY      = -999.;
1614       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1615       vPt     = aodTrack->Pt();
1616            
1617       
1618       // Kinematics cuts from ESD track cuts
1619       if( vPt < fPtMin || vPt > fPtMax)      continue;
1620       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1621       
1622        
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);
1632       
1633       //=======================================correction
1634       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1635       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1636       
1637       // add the track to the TObjArray
1638       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1639     }//track loop
1640   }// AOD nano analysis
1641
1642
1643   //==============================================================================================================
1644   else if(gAnalysisLevel == "MCAOD") {
1645     
1646     AliMCEvent* mcEvent = MCEvent();
1647     if (!mcEvent) {
1648       AliError("ERROR: Could not retrieve MC event");
1649     }
1650     else{
1651       
1652       for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1653         AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); 
1654         if (!aodTrack) {
1655           AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1656           continue;
1657         }
1658         
1659         if(!aodTrack->IsPhysicalPrimary()) continue;   
1660         
1661         vCharge = aodTrack->Charge();
1662         vEta    = aodTrack->Eta();
1663         vY      = aodTrack->Y();
1664         vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1665         vPt     = aodTrack->Pt();
1666         
1667         // Kinematics cuts from ESD track cuts
1668         if( vPt < fPtMin || vPt > fPtMax)      continue;
1669         if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1670         
1671         // Remove neutral tracks
1672         if( vCharge == 0 ) continue;
1673         
1674         //Exclude resonances
1675         if(fExcludeResonancesInMC) {
1676           
1677           Bool_t kExcludeParticle = kFALSE;
1678           Int_t gMotherIndex = aodTrack->GetMother();
1679           if(gMotherIndex != -1) {
1680             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1681             if(motherTrack) {
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
1697                  ) {
1698                 kExcludeParticle = kTRUE;
1699               }
1700             }
1701           }
1702           
1703           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1704           if(kExcludeParticle) continue;
1705         }
1706
1707         //Exclude electrons with PDG
1708         if(fExcludeElectronsInMC) {
1709           
1710           if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1711           
1712         }
1713         
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);
1723         
1724         //=======================================correction
1725         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1726         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   
1727         
1728         // add the track to the TObjArray
1729         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1730       }//aodTracks
1731     }//MC event
1732   }//MCAOD
1733   //==============================================================================================================
1734
1735   //==============================================================================================================
1736   else if(gAnalysisLevel == "MCAODrec") {
1737     
1738     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1739     if (!fAOD) {
1740       printf("ERROR: fAOD not available\n");
1741       return;
1742       }*/
1743     
1744     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); 
1745     if (!fArrayMC) { 
1746        AliError("No array of MC particles found !!!");
1747     }
1748
1749     AliMCEvent* mcEvent = MCEvent();
1750     if (!mcEvent) {
1751        AliError("ERROR: Could not retrieve MC event");
1752        return tracksAccepted;
1753     }
1754      
1755     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1756       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1757       if (!aodTrack) {
1758         AliError(Form("Could not receive track %d", iTracks));
1759         continue;
1760       }
1761       
1762       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1763         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1764       }
1765       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1766       
1767       vCharge = aodTrack->Charge();
1768       vEta    = aodTrack->Eta();
1769       vY      = aodTrack->Y();
1770       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1771       vPt     = aodTrack->Pt();
1772       
1773       //===========================use MC information for Kinematics===============================//               
1774       if(fUseMCforKinematics){
1775
1776         Int_t label = TMath::Abs(aodTrack->GetLabel());
1777         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1778
1779         if(AODmcTrack){
1780           vCharge = AODmcTrack->Charge();
1781           vEta    = AODmcTrack->Eta();
1782           vY      = AODmcTrack->Y();
1783           vPhi    = AODmcTrack->Phi();// * TMath::RadToDeg();
1784           vPt     = AODmcTrack->Pt();
1785         }
1786         else{
1787           AliDebug(1, "no MC particle for this track"); 
1788           continue;
1789         }
1790       }
1791       //===========================end of use MC information for Kinematics========================//               
1792
1793
1794       //===========================PID (so far only for electron rejection)===============================//                
1795       if(fElectronRejection) {
1796
1797         // get the electron nsigma
1798         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1799
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
1804         
1805         // check only for given momentum range
1806         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1807                   
1808           //look only at electron nsigma
1809           if(!fElectronOnlyRejection){
1810             
1811             //Make the decision based on the n-sigma of electrons
1812             if(nSigma < fElectronRejectionNSigma) continue;
1813           }
1814           else{
1815             
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));
1819             
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;
1825           }
1826         }
1827   
1828         //Fill QA after the PID
1829         fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1830         fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1831         
1832       }
1833       //===========================end of PID (so far only for electron rejection)===============================//
1834       
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)
1837             
1838       // Kinematics cuts from ESD track cuts
1839       if( vPt < fPtMin || vPt > fPtMax)      continue;
1840       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1841       
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 ){
1845           continue;  // 2D cut
1846         }
1847       }
1848       
1849       // Extra TPC cuts (for systematic studies [!= -1])
1850       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1851         continue;
1852       }
1853       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1854         continue;
1855       }
1856
1857      // Extra cut on shared clusters
1858       if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1859         continue;
1860       }
1861
1862       //Exclude resonances
1863       if(fExcludeResonancesInMC) {
1864         
1865         Bool_t kExcludeParticle = kFALSE;
1866
1867         Int_t label = TMath::Abs(aodTrack->GetLabel());
1868         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1869       
1870         if (AODmcTrack){ 
1871           //if (AODmcTrack->IsPhysicalPrimary()){
1872           
1873           Int_t gMotherIndex = AODmcTrack->GetMother();
1874           if(gMotherIndex != -1) {
1875             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1876             if(motherTrack) {
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
1890                  ) {
1891                 kExcludeParticle = kTRUE;
1892               }
1893             }
1894           }
1895         }       
1896         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1897         if(kExcludeParticle) continue;
1898       }
1899
1900       //Exclude electrons with PDG
1901       if(fExcludeElectronsInMC) {
1902         
1903         Int_t label = TMath::Abs(aodTrack->GetLabel());
1904         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1905         
1906         if (AODmcTrack){ 
1907           if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1908         }
1909       }
1910       
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);
1923       
1924       //=======================================correction
1925       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1926       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1927       
1928       // add the track to the TObjArray
1929       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1930     }//track loop
1931   }//MCAODrec
1932   //==============================================================================================================
1933
1934   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1935
1936     AliESDtrack *trackTPC   = NULL;
1937     AliMCParticle *mcTrack   = NULL;
1938
1939     AliMCEvent*  mcEvent     = NULL;
1940     
1941     // for MC ESDs use also MC information
1942     if(gAnalysisLevel == "MCESD"){
1943       mcEvent = MCEvent(); 
1944       if (!mcEvent) {
1945         AliError("mcEvent not available");
1946         return tracksAccepted;
1947       }
1948     }
1949     
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));
1953       if (!track) {
1954         AliError(Form("Could not receive track %d", iTracks));
1955         continue;
1956       } 
1957
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;
1963         
1964         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1965         if(!mcTrack) continue;
1966       }
1967
1968       // take only TPC only tracks
1969       trackTPC   = new AliESDtrack();
1970       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1971       
1972       //ESD track cuts
1973       if(fESDtrackCuts) 
1974         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1975       
1976       // fill QA histograms
1977       Float_t b[2];
1978       Float_t bCov[3];
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;
1983       }
1984       
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
1992       }
1993       
1994       //===========================PID===============================//             
1995       if(fUsePID) {
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.};
2000         
2001         Double_t nSigma = 0.;
2002         UInt_t detUsedTPC = 0;
2003         UInt_t detUsedTOF = 0;
2004         UInt_t detUsedTPCTOF = 0;
2005         
2006         //Decide what detector configuration we want to use
2007         switch(fPidDetectorConfig) {
2008         case kTPCpid:
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];
2014           break;
2015         case kTOFpid:
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];
2021           break;
2022         case kTPCTOF:
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];
2027           break;
2028         default:
2029           break;
2030         }//end switch: define detector mask
2031         
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     
2042           
2043           if (tof <= 0) {
2044             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
2045             continue;
2046           }
2047           if (length <= 0){
2048             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
2049             continue;
2050           }
2051           
2052           length = length*0.01; // in meters
2053           tof = tof*c;
2054           beta = length/tof;
2055           
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);
2060         }//TOF signal 
2061         
2062         
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
2069         
2070         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
2071           //Make the decision based on the n-sigma
2072           if(fUsePIDnSigma) {
2073             if(nSigma > fPIDNSigma) continue;}
2074           
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;      
2079           }
2080           
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);
2085           
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); 
2090         }
2091       }
2092       //===========================PID===============================//
2093       vCharge = trackTPC->Charge();
2094       vY      = trackTPC->Y();
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);
2109       
2110       //=======================================correction
2111       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
2112       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2113
2114       // add the track to the TObjArray
2115       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   
2116
2117       delete trackTPC;
2118     }//track loop
2119   }// ESD analysis
2120
2121   else if(gAnalysisLevel == "MC"){
2122     if(!event) {
2123       AliError("mcEvent not available");
2124       return 0x0;
2125     }
2126
2127     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2128     if(gMCEvent) {
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));
2132         if (!track) {
2133           AliError(Form("Could not receive particle %d", iTracks));
2134           continue;
2135         }
2136         
2137         //exclude non stable particles
2138         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2139         
2140         vCharge = track->Charge();
2141         vEta    = track->Eta();
2142         vPt     = track->Pt();
2143         vY      = track->Y();
2144         
2145         if( vPt < fPtMin || vPt > fPtMax)      
2146           continue;
2147         if (!fUsePID) {
2148           if( vEta < fEtaMin || vEta > fEtaMax)  continue;
2149         }
2150         else if (fUsePID){
2151           if( vY < fEtaMin || vY > fEtaMax)  continue;
2152         }
2153
2154         // Remove neutral tracks
2155         if( vCharge == 0 ) continue;
2156         
2157         //analyze one set of particles
2158         if(fUseMCPdgCode) {
2159           TParticle *particle = track->Particle();
2160           if(!particle) continue;
2161           
2162           Int_t gPdgCode = particle->GetPdgCode();
2163           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) 
2164             continue;
2165         }
2166         
2167         //Use the acceptance parameterization
2168         if(fAcceptanceParameterization) {
2169           Double_t gRandomNumber = gRandom->Rndm();
2170           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) 
2171             continue;
2172         }
2173         
2174         //Exclude resonances
2175         if(fExcludeResonancesInMC) {
2176           TParticle *particle = track->Particle();
2177           if(!particle) continue;
2178           
2179           Bool_t kExcludeParticle = kFALSE;
2180           Int_t gMotherIndex = particle->GetFirstMother();
2181
2182           if(TMath::Abs(particle->GetPdgCode()) != 321) {
2183             // AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2184             // if(motherTrack) {
2185             //   TParticle *motherParticle = motherTrack->Particle();
2186             //   if(motherParticle) {
2187             //  Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2188             //  //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2189             //  if(pdgCodeOfMother == 113  // rho0
2190             //     || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2191             //     // || pdgCodeOfMother == 221  // eta
2192             //     // || pdgCodeOfMother == 331  // eta'
2193             //     // || pdgCodeOfMother == 223  // omega
2194             //     // || pdgCodeOfMother == 333  // phi
2195             //     || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0
2196             //     // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*
2197             //     // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*
2198             //     || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2199             //     || pdgCodeOfMother == 111  // pi0 Dalitz
2200             //     ) {
2201             kExcludeParticle = kTRUE;
2202             //  }
2203             //        }
2204             //}
2205           }
2206           
2207           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2208           if(kExcludeParticle) continue;
2209
2210         }
2211
2212         //Exclude electrons with PDG
2213         if(fExcludeElectronsInMC) {
2214           
2215           TParticle *particle = track->Particle();
2216           
2217           if (particle){ 
2218             if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2219           }
2220         }
2221       
2222         vPhi    = track->Phi();
2223         //Printf("phi (before): %lf",vPhi);
2224         
2225         fHistPt->Fill(vPt,gCentrality);
2226         fHistEta->Fill(vEta,gCentrality);
2227         fHistPhi->Fill(vPhi,gCentrality);
2228         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2229         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2230         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2231         fHistRapidity->Fill(vY,gCentrality);
2232         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2233         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2234         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2235         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2236         
2237         //Flow after burner
2238         if(fUseFlowAfterBurner) {
2239           Double_t precisionPhi = 0.001;
2240           Int_t maxNumberOfIterations = 100;
2241           
2242           Double_t phi0 = vPhi;
2243           Double_t gV2 = fDifferentialV2->Eval(vPt);
2244           
2245           for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2246             Double_t phiprev = vPhi;
2247             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2248             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); 
2249             vPhi -= fl/fp;
2250             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2251           }
2252           //Printf("phi (after): %lf\n",vPhi);
2253           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2254           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2255           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2256           
2257           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2258           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2259           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2260           
2261         }
2262         
2263         //vPhi *= TMath::RadToDeg();
2264         
2265       //=======================================correction
2266       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
2267       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2268
2269         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); 
2270       } //track loop
2271     }//MC event object
2272   }//MC
2273   
2274   return tracksAccepted;  
2275 }
2276
2277 //________________________________________________________________________
2278 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2279   // Clones TObjArray and returns it with tracks after shuffling the charges
2280
2281   TObjArray* tracksShuffled = new TObjArray;
2282   tracksShuffled->SetOwner(kTRUE);
2283
2284   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks 
2285
2286   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2287   {
2288     AliVParticle* track = (AliVParticle*) tracks->At(i);
2289     chargeVector->push_back(track->Charge());
2290   }  
2291  
2292   random_shuffle(chargeVector->begin(), chargeVector->end());
2293   
2294   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2295     AliVParticle* track = (AliVParticle*) tracks->At(i);
2296     //==============================correction
2297     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2298     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2299     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2300   }
2301
2302   delete chargeVector;
2303    
2304   return tracksShuffled;
2305 }
2306
2307 //________________________________________________________________________
2308 void  AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2309                                                     const char* lhcPeriod) {
2310   //Function to setup the VZERO gain equalization
2311     //============Get the equilization map============//
2312   TFile *calibrationFile = TFile::Open(filename);
2313   if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2314     Printf("No calibration file found!!!");
2315     return;
2316   }
2317
2318   TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2319   if(!list) {
2320     Printf("Calibration TList not found!!!");
2321     return;
2322   }
2323
2324   fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2325   if(!fHistVZEROAGainEqualizationMap) {
2326     Printf("VZERO-A calibration object not found!!!");
2327     return;
2328   }
2329   fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2330   if(!fHistVZEROCGainEqualizationMap) {
2331     Printf("VZERO-C calibration object not found!!!");
2332     return;
2333   }
2334
2335   fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2336   if(!fHistVZEROChannelGainEqualizationMap) {
2337     Printf("VZERO channel calibration object not found!!!");
2338     return;
2339   }
2340 }
2341
2342 //________________________________________________________________________
2343 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, 
2344                                                             Int_t channel) {
2345   //
2346   if(!fHistVZEROAGainEqualizationMap) return 1.0;
2347
2348   for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2349     Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2350     if(gRunNumber == run)
2351       return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2352   }
2353
2354   return 1.0;
2355 }
2356
2357 //________________________________________________________________________
2358 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, 
2359                                                      const char* side) {
2360   //
2361   if(!fHistVZEROAGainEqualizationMap) return 1.0;
2362
2363   TString gVZEROSide = side;
2364   for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2365     Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2366     //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2367     if(gRunNumber == run) {
2368       if(gVZEROSide == "A") 
2369         return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2370       else if(gVZEROSide == "C") 
2371         return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2372     }
2373   }
2374
2375   return 1.0;
2376 }
2377
2378 //____________________________________________________________________
2379 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2380 {
2381   // copied from AliAnalysisTaskPhiCorrelations
2382   //
2383   // rejects "randomly" events such that the centrality gets flat
2384   // uses fCentralityWeights histogram
2385
2386   // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2387   
2388   Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2389   Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2390   
2391   Bool_t result = kFALSE;
2392   if (centralityDigits < weight) 
2393     result = kTRUE;
2394   
2395   AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2396   
2397   return result;
2398 }
2399
2400 //________________________________________________________________________
2401 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){
2402   //Printf("END BF");
2403
2404   if (!fBalance) {
2405     AliError("fBalance not available");
2406     return;
2407   }  
2408   if(fRunShuffling) {
2409     if (!fShuffledBalance) {
2410       AliError("fShuffledBalance not available");
2411       return;
2412     }
2413   }
2414
2415 }
2416
2417 //________________________________________________________________________
2418 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2419   // Draw result to the screen
2420   // Called once at the end of the query
2421
2422   // not implemented ...
2423
2424 }
2425