bugfix: Use FindBin for correction matrix (otherwise bin mismatch for pT which has...
[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
15 #include "AliAnalysisTaskSE.h"
16 #include "AliAnalysisManager.h"
17
18 #include "AliESDVertex.h"
19 #include "AliESDEvent.h"
20 #include "AliESDInputHandler.h"
21 #include "AliAODEvent.h"
22 #include "AliAODTrack.h"
23 #include "AliAODInputHandler.h"
24 #include "AliAODMCParticle.h" 
25 #include "AliCollisionGeometry.h"
26 #include "AliGenEventHeader.h"
27 #include "AliMCEventHandler.h"
28 #include "AliMCEvent.h"
29 #include "AliStack.h"
30 #include "AliESDtrackCuts.h"
31 #include "AliEventplane.h"
32 #include "AliTHn.h"    
33 #include "AliLog.h"
34 #include "AliAnalysisUtils.h"
35
36 #include "AliEventPoolManager.h"           
37
38 #include "AliPID.h"                
39 #include "AliPIDResponse.h"        
40 #include "AliPIDCombined.h"        
41
42 #include "AliAnalysisTaskBFPsi.h"
43 #include "AliBalancePsi.h"
44 #include "AliAnalysisTaskTriggeredBF.h"
45
46
47 // Analysis task for the BF vs Psi code
48 // Authors: Panos.Christakoglou@nikhef.nl
49
50 using std::cout;
51 using std::endl;
52
53 ClassImp(AliAnalysisTaskBFPsi)
54
55 //________________________________________________________________________
56 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) 
57 : AliAnalysisTaskSE(name),
58   fDebugLevel(kFALSE),
59   fArrayMC(0), //+++++++++++++
60   fBalance(0),
61   fRunShuffling(kFALSE),
62   fShuffledBalance(0),
63   fRunMixing(kFALSE),
64   fRunMixingEventPlane(kFALSE),
65   fMixingTracks(50000),
66   fMixedBalance(0),
67   fPoolMgr(0),
68   fList(0),
69   fListBF(0),
70   fListBFS(0),
71   fListBFM(0),
72   fHistListPIDQA(0),
73   fHistEventStats(0),
74   fHistCentStats(0),
75   fHistCentStatsUsed(0),
76   fHistTriggerStats(0),
77   fHistTrackStats(0),
78   fHistVx(0),
79   fHistVy(0),
80   fHistVz(0),
81   fHistTPCvsVZEROMultiplicity(0),
82   fHistVZEROSignal(0),
83   fHistEventPlane(0),
84   fHistClus(0),
85   fHistDCA(0),
86   fHistChi2(0),
87   fHistPt(0),
88   fHistEta(0),
89   fHistRapidity(0),
90   fHistPhi(0),
91   fHistEtaPhiPos(0),             
92   fHistEtaPhiNeg(0), 
93   fHistPhiBefore(0),
94   fHistPhiAfter(0),
95   fHistPhiPos(0),
96   fHistPhiNeg(0),
97   fHistV0M(0),
98   fHistRefTracks(0),
99   fHistdEdxVsPTPCbeforePID(NULL),
100   fHistBetavsPTOFbeforePID(NULL), 
101   fHistProbTPCvsPtbeforePID(NULL), 
102   fHistProbTOFvsPtbeforePID(NULL), 
103   fHistProbTPCTOFvsPtbeforePID(NULL),
104   fHistNSigmaTPCvsPtbeforePID(NULL), 
105   fHistNSigmaTOFvsPtbeforePID(NULL), 
106   fHistBetaVsdEdXbeforePID(NULL), //+++++++ 
107   fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++
108   fHistdEdxVsPTPCafterPID(NULL),
109   fHistBetavsPTOFafterPID(NULL), 
110   fHistProbTPCvsPtafterPID(NULL), 
111   fHistProbTOFvsPtafterPID(NULL), 
112   fHistProbTPCTOFvsPtafterPID(NULL),
113   fHistNSigmaTPCvsPtafterPID(NULL), 
114   fHistNSigmaTOFvsPtafterPID(NULL),  
115   fHistBetaVsdEdXafterPID(NULL), //+++++++ 
116   fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++
117   fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++
118   fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++
119   fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++
120   fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++
121   fCentralityArrayBinsForCorrections(kCENTRALITY),
122   fCentralityWeights(0x0),
123   fPIDResponse(0x0),
124   fPIDCombined(0x0),
125   fParticleOfInterest(kPion),
126   fPidDetectorConfig(kTPCTOF),
127   fUsePID(kFALSE),
128   fUsePIDnSigma(kTRUE),
129   fUsePIDPropabilities(kFALSE), 
130   fPIDNSigma(3.),
131   fMinAcceptedPIDProbability(0.8),
132   fElectronRejection(kFALSE),
133   fElectronOnlyRejection(kFALSE),
134   fElectronRejectionNSigma(-1.),
135   fElectronRejectionMinPt(0.),
136   fElectronRejectionMaxPt(1000.),
137   fESDtrackCuts(0),
138   fCentralityEstimator("V0M"),
139   fUseCentrality(kFALSE),
140   fCentralityPercentileMin(0.), 
141   fCentralityPercentileMax(5.),
142   fImpactParameterMin(0.),
143   fImpactParameterMax(20.),
144   fMultiplicityEstimator("V0A"),
145   fUseMultiplicity(kFALSE),
146   fNumberOfAcceptedTracksMin(0),
147   fNumberOfAcceptedTracksMax(10000),
148   fHistNumberOfAcceptedTracks(0),
149   fHistMultiplicity(0),
150   fUseOfflineTrigger(kFALSE),
151   fCheckFirstEventInChunk(kFALSE),
152   fCheckPileUp(kFALSE),
153   fCheckPrimaryFlagAOD(kFALSE),
154   fUseMCforKinematics(kFALSE),
155   fVxMax(0.3),
156   fVyMax(0.3),
157   fVzMax(10.),
158   fnAODtrackCutBit(128),
159   fPtMin(0.3),
160   fPtMax(1.5),
161   fEtaMin(-0.8),
162   fEtaMax(-0.8),
163   fPhiMin(0.),
164   fPhiMax(360.),
165   fDCAxyCut(-1),
166   fDCAzCut(-1),
167   fTPCchi2Cut(-1),
168   fNClustersTPCCut(-1),
169   fAcceptanceParameterization(0),
170   fDifferentialV2(0),
171   fUseFlowAfterBurner(kFALSE),
172   fExcludeResonancesInMC(kFALSE),
173   fExcludeElectronsInMC(kFALSE),
174   fUseMCPdgCode(kFALSE),
175   fPDGCodeToBeAnalyzed(-1),
176   fEventClass("EventPlane"), 
177   fCustomBinning(""),
178   fHistVZEROAGainEqualizationMap(0),
179   fHistVZEROCGainEqualizationMap(0),
180   fHistVZEROChannelGainEqualizationMap(0) {
181   // Constructor
182   // Define input and output slots here
183   // Input slot #0 works with a TChain
184
185   //======================================================correction
186   for (Int_t i=0; i<kCENTRALITY; i++){
187     fHistCorrectionPlus[i] = NULL; 
188     fHistCorrectionMinus[i] = NULL; 
189     fCentralityArrayForCorrections[i] = -1.;
190   }
191   //=====================================================correction
192
193   DefineInput(0, TChain::Class());
194   // Output slot #0 writes into a TH1 container
195   DefineOutput(1, TList::Class());
196   DefineOutput(2, TList::Class());
197   DefineOutput(3, TList::Class());
198   DefineOutput(4, TList::Class());
199   DefineOutput(5, TList::Class());
200 }
201
202 //________________________________________________________________________
203 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
204
205   // delete fBalance; 
206   // delete fShuffledBalance; 
207   // delete fList;
208   // delete fListBF; 
209   // delete fListBFS;
210
211   // delete fHistEventStats; 
212   // delete fHistTrackStats; 
213   // delete fHistVx; 
214   // delete fHistVy; 
215   // delete fHistVz; 
216
217   // delete fHistClus;
218   // delete fHistDCA;
219   // delete fHistChi2;
220   // delete fHistPt;
221   // delete fHistEta;
222   // delete fHistPhi;
223   // delete fHistEtaPhiPos;                      
224   // delete fHistEtaPhiNeg;
225   // delete fHistV0M;
226 }
227
228 //________________________________________________________________________
229 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
230   // Create histograms
231   // Called once
232
233   // global switch disabling the reference 
234   // (to avoid "Replacing existing TH1" if several wagons are created in train)
235   Bool_t oldStatus = TH1::AddDirectoryStatus();
236   TH1::AddDirectory(kFALSE);
237
238   if(!fBalance) {
239     fBalance = new AliBalancePsi();
240     fBalance->SetAnalysisLevel("ESD");
241     fBalance->SetEventClass(fEventClass);
242     //fBalance->SetNumberOfBins(-1,16);
243     //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
244   }
245   if(fRunShuffling) {
246     if(!fShuffledBalance) {
247       fShuffledBalance = new AliBalancePsi();
248       fShuffledBalance->SetAnalysisLevel("ESD");
249       fShuffledBalance->SetEventClass(fEventClass);
250       //fShuffledBalance->SetNumberOfBins(-1,16);
251       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
252     }
253   }
254   if(fRunMixing) {
255     if(!fMixedBalance) {
256       fMixedBalance = new AliBalancePsi();
257       fMixedBalance->SetAnalysisLevel("ESD");
258       fMixedBalance->SetEventClass(fEventClass);
259     }
260   }
261
262   //QA list
263   fList = new TList();
264   fList->SetName("listQA");
265   fList->SetOwner();
266
267   //Balance Function list
268   fListBF = new TList();
269   fListBF->SetName("listBF");
270   fListBF->SetOwner();
271
272   if(fRunShuffling) {
273     fListBFS = new TList();
274     fListBFS->SetName("listBFShuffled");
275     fListBFS->SetOwner();
276   }
277
278   if(fRunMixing) {
279     fListBFM = new TList();
280     fListBFM->SetName("listTriggeredBFMixed");
281     fListBFM->SetOwner();
282   }
283
284   //PID QA list
285   if(fUsePID || fElectronRejection) {
286     fHistListPIDQA = new TList();
287     fHistListPIDQA->SetName("listQAPID");
288     fHistListPIDQA->SetOwner();
289   }
290
291   //Event stats.
292   TString gCutName[7] = {"Total","Offline trigger",
293                          "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
294   fHistEventStats = new TH2F("fHistEventStats",
295                              "Event statistics;;Centrality percentile;N_{events}",
296                              7,0.5,7.5,220,-5,105);
297   for(Int_t i = 1; i <= 7; i++)
298     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
299   fList->Add(fHistEventStats);
300
301   TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
302   fHistCentStats = new TH2F("fHistCentStats",
303                              "Centrality statistics;;Cent percentile",
304                             13,-0.5,12.5,220,-5,105);
305   for(Int_t i = 1; i <= 13; i++){
306     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
307     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
308   }
309   fList->Add(fHistCentStats);
310
311   fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
312   fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
313   fList->Add(fHistCentStatsUsed);
314
315   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
316   fList->Add(fHistTriggerStats);
317
318   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
319   fList->Add(fHistTrackStats);
320
321   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
322   fList->Add(fHistNumberOfAcceptedTracks);
323
324   fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
325   fList->Add(fHistMultiplicity);
326
327   // Vertex distributions
328   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
329   fList->Add(fHistVx);
330   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
331   fList->Add(fHistVy);
332   fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
333   fList->Add(fHistVz);
334
335   //TPC vs VZERO multiplicity
336   fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
337   if(fMultiplicityEstimator == "V0A") 
338     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
339   else if(fMultiplicityEstimator == "V0C") 
340     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
341   else 
342     fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
343   fList->Add(fHistTPCvsVZEROMultiplicity);
344
345   fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
346   fList->Add(fHistVZEROSignal);
347
348   //Event plane
349   fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
350   fList->Add(fHistEventPlane);
351
352   // QA histograms
353   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
354   fList->Add(fHistClus);
355   fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
356   fList->Add(fHistChi2);
357   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); 
358   fList->Add(fHistDCA);
359   fHistPt   = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
360   fList->Add(fHistPt);
361   fHistEta  = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
362   fList->Add(fHistEta);
363   fHistRapidity  = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
364   fList->Add(fHistRapidity);
365   fHistPhi  = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
366   fList->Add(fHistPhi);
367   fHistEtaPhiPos  = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);                      
368   fList->Add(fHistEtaPhiPos);                    
369   fHistEtaPhiNeg  = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105);              
370   fList->Add(fHistEtaPhiNeg);
371   fHistPhiBefore  = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
372   fList->Add(fHistPhiBefore);
373   fHistPhiAfter  = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
374   fList->Add(fHistPhiAfter);
375   fHistPhiPos  = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
376   fList->Add(fHistPhiPos);
377   fHistPhiNeg  = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
378   fList->Add(fHistPhiNeg);
379   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
380   fList->Add(fHistV0M);
381   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
382   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
383   for(Int_t i = 1; i <= 6; i++)
384     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
385   fList->Add(fHistRefTracks);
386
387   // Balance function histograms
388   // Initialize histograms if not done yet (including the custom binning)
389   if(!fBalance->GetHistNp()){
390     AliInfo("Histograms not yet initialized! --> Will be done now");
391     fBalance->SetCustomBinning(fCustomBinning);
392     fBalance->InitHistograms();
393   }
394
395   if(fRunShuffling) {
396     if(!fShuffledBalance->GetHistNp()) {
397       AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
398       fShuffledBalance->SetCustomBinning(fCustomBinning);
399       fShuffledBalance->InitHistograms();
400     }
401   }
402
403   if(fRunMixing) {
404     if(!fMixedBalance->GetHistNp()) {
405       AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
406       fMixedBalance->SetCustomBinning(fCustomBinning);
407       fMixedBalance->InitHistograms();
408     }
409   }
410
411   // QA histograms for different cuts
412   fList->Add(fBalance->GetQAHistHBTbefore());
413   fList->Add(fBalance->GetQAHistHBTafter());
414   fList->Add(fBalance->GetQAHistConversionbefore());
415   fList->Add(fBalance->GetQAHistConversionafter());
416   fList->Add(fBalance->GetQAHistPsiMinusPhi());
417   fList->Add(fBalance->GetQAHistResonancesBefore());
418   fList->Add(fBalance->GetQAHistResonancesRho());
419   fList->Add(fBalance->GetQAHistResonancesK0());
420   fList->Add(fBalance->GetQAHistResonancesLambda());
421   fList->Add(fBalance->GetQAHistQbefore());
422   fList->Add(fBalance->GetQAHistQafter());
423
424   //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
425   fListBF->Add(fBalance->GetHistNp());
426   fListBF->Add(fBalance->GetHistNn());
427   fListBF->Add(fBalance->GetHistNpn());
428   fListBF->Add(fBalance->GetHistNnn());
429   fListBF->Add(fBalance->GetHistNpp());
430   fListBF->Add(fBalance->GetHistNnp());
431
432   if(fRunShuffling) {
433     fListBFS->Add(fShuffledBalance->GetHistNp());
434     fListBFS->Add(fShuffledBalance->GetHistNn());
435     fListBFS->Add(fShuffledBalance->GetHistNpn());
436     fListBFS->Add(fShuffledBalance->GetHistNnn());
437     fListBFS->Add(fShuffledBalance->GetHistNpp());
438     fListBFS->Add(fShuffledBalance->GetHistNnp());
439   }  
440
441   if(fRunMixing) {
442     fListBFM->Add(fMixedBalance->GetHistNp());
443     fListBFM->Add(fMixedBalance->GetHistNn());
444     fListBFM->Add(fMixedBalance->GetHistNpn());
445     fListBFM->Add(fMixedBalance->GetHistNnn());
446     fListBFM->Add(fMixedBalance->GetHistNpp());
447     fListBFM->Add(fMixedBalance->GetHistNnp());
448   }
449   //}
450
451
452   // Event Mixing
453   if(fRunMixing){
454     Int_t trackDepth = fMixingTracks; 
455     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
456     
457     // centrality bins
458     Double_t* centbins = NULL;
459     Int_t nCentralityBins;
460     if(fBalance->IsUseVertexBinning()){
461       centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
462     }
463     else{
464       centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
465     }
466     
467     // multiplicity bins
468     Double_t* multbins = NULL;
469     Int_t nMultiplicityBins;
470     multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
471     
472     // Zvtx bins
473     Double_t* vtxbins = NULL; 
474     Int_t nVertexBins;
475     if(fBalance->IsUseVertexBinning()){
476       vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
477     }
478     else{
479       vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
480     }
481
482     // Event plane angle (Psi) bins
483     Double_t* psibins = NULL;
484     Int_t nPsiBins; 
485     psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
486
487     
488     // run the event mixing also in bins of event plane (statistics!)
489     if(fRunMixingEventPlane){
490       if(fEventClass=="Multiplicity"){
491         if(multbins && vtxbins && psibins){
492           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
493         }
494       }
495       else{
496         if(centbins && vtxbins && psibins){
497           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
498         }
499       }
500     }
501     else{
502       if(fEventClass=="Multiplicity"){
503         if(multbins && vtxbins){
504           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
505         }
506       }
507       else{
508         if(centbins && vtxbins){
509           fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
510         }
511       }
512     }
513
514     // check pool manager
515     if(!fPoolMgr){
516       AliError("Event Mixing required, but Pool Manager not initialized...");
517       return;
518     }
519   }
520   
521   if(fESDtrackCuts) fList->Add(fESDtrackCuts);
522
523   //====================PID========================//
524   if(fUsePID) {
525     fPIDCombined = new AliPIDCombined();
526     fPIDCombined->SetDefaultTPCPriors();
527
528     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); 
529     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
530     
531     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); 
532     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); 
533     
534     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); 
535     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); 
536     
537     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); 
538     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
539
540     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); 
541     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
542     
543     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); 
544     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
545     
546     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); 
547     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); 
548
549     fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); 
550     fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
551     
552     fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500); 
553     fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
554     
555     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); 
556     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
557     
558     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); 
559     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); 
560     
561     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); 
562     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
563   
564     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); 
565     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); 
566     
567     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); 
568     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
569
570     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); 
571     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
572     
573     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); 
574     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
575
576     fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); 
577     fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
578
579     fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500); 
580     fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
581   }
582
583   // for electron rejection only TPC nsigma histograms
584   if(fElectronRejection) {
585  
586     fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000); 
587     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
588     
589     fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500); 
590     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
591     
592     fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000); 
593     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
594
595     fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500); 
596     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron); 
597   }
598   //====================PID========================//
599
600   // Post output data.
601   PostData(1, fList);
602   PostData(2, fListBF);
603   if(fRunShuffling) PostData(3, fListBFS);
604   if(fRunMixing) PostData(4, fListBFM);
605   if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA);       //PID
606
607   AliInfo("Finished setting up the Output");
608
609   TH1::AddDirectory(oldStatus);
610 }
611
612
613 //________________________________________________________________________
614 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, 
615                                               Int_t nCentralityBins, 
616                                               Double_t *centralityArrayForCorrections) {
617   //Open files that will be used for correction
618   fCentralityArrayBinsForCorrections = nCentralityBins;
619   for (Int_t i=0; i<nCentralityBins; i++)
620     fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
621
622   // No file specified -> run without corrections
623   if(!filename.Contains(".root")) {
624     AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
625     return;
626   }
627
628   //Open the input file
629   TFile *f = TFile::Open(filename);
630   if(!f->IsOpen()) {
631     AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
632     return;
633   }
634     
635   //TString listEffName = "";
636   for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {    
637     //Printf("iCent %d:",iCent);    
638     TString histoName = "fHistCorrectionPlus";
639     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
640     fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
641     if(!fHistCorrectionPlus[iCent]) {
642       AliError(Form("fHist %s not found!!!",histoName.Data()));
643       return;
644     }
645     
646     histoName = "fHistCorrectionMinus";
647     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
648     fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); 
649     if(!fHistCorrectionMinus[iCent]) {
650       AliError(Form("fHist %s not found!!!",histoName.Data()));
651       return; 
652     }
653   }//loop over centralities: ONLY the PbPb case is covered
654 }
655
656 //________________________________________________________________________
657 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
658   // Main loop
659   // Called for each event
660
661   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
662   Int_t gNumberOfAcceptedTracks = 0;
663   Double_t lMultiplicityVar     = -999.; //-1
664   Double_t gReactionPlane       = -1.; 
665   Float_t bSign = 0.;
666   
667   // get the event (for generator level: MCEvent())
668   AliVEvent* eventMain = NULL;
669   if(gAnalysisLevel == "MC") {
670     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); 
671   }
672   else{
673     eventMain = dynamic_cast<AliVEvent*>(InputEvent());     
674     // for HBT like cuts need magnetic field sign
675     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
676   }
677   if(!eventMain) {
678     AliError("eventMain not available");
679     return;
680   }
681   
682   // PID Response task active?
683   if(fUsePID || fElectronRejection) {
684     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
685     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
686   }
687  
688   // check event cuts and fill event histograms
689   if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ 
690     return;
691   }
692   // get the reaction plane
693   if(fEventClass != "Multiplicity") {
694     gReactionPlane = GetEventPlane(eventMain);
695     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
696     if(gReactionPlane < 0){
697       return;
698     }
699   }
700   
701   // get the accepted tracks in main event
702   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
703   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
704
705   //multiplicity cut (used in pp)
706   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
707
708   // store charges of all accepted tracks,shuffle and reassign(two extra loops)
709   TObjArray* tracksShuffled = NULL;
710   if(fRunShuffling){
711     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
712   }
713   
714   // Event mixing 
715   if (fRunMixing)
716     {
717       // 1. First get an event pool corresponding in mult (cent) and
718       //    zvertex to the current event. Once initialized, the pool
719       //    should contain nMix (reduced) events. This routine does not
720       //    pre-scan the chain. The first several events of every chain
721       //    will be skipped until the needed pools are filled to the
722       //    specified depth. If the pool categories are not too rare, this
723       //    should not be a problem. If they are rare, you could lose`
724       //    statistics.
725       
726       // 2. Collect the whole pool's content of tracks into one TObjArray
727       //    (bgTracks), which is effectively a single background super-event.
728       
729       // 3. The reduced and bgTracks arrays must both be passed into
730       //    FillCorrelations(). Also nMix should be passed in, so a weight
731       //    of 1./nMix can be applied.
732       
733       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
734       
735       if (!pool){
736         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
737       }
738       else{
739         
740         //pool->SetDebug(1);
741
742         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ 
743           
744           
745           Int_t nMix = pool->GetCurrentNEvents();
746           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
747           
748           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
749           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
750           //if (pool->IsReady())
751           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
752           
753           // Fill mixed-event histos here  
754           for (Int_t jMix=0; jMix<nMix; jMix++) 
755             {
756               TObjArray* tracksMixed = pool->GetEvent(jMix);
757               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
758             }
759         }
760         
761         // Update the Event pool
762         pool->UpdatePool(tracksMain);
763         //pool->PrintInfo();
764         
765       }//pool NULL check  
766     }//run mixing
767   
768   // calculate balance function
769   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
770   
771   // calculate shuffled balance function
772   if(fRunShuffling && tracksShuffled != NULL) {
773     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
774   }
775 }      
776
777 //________________________________________________________________________
778 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
779   // Checks the Event cuts
780   // Fills Event statistics histograms
781   
782   Bool_t isSelectedMain = kTRUE;
783   Float_t gRefMultiplicity = -1.;
784   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
785
786   AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
787
788   fHistEventStats->Fill(1,gRefMultiplicity); //all events
789
790   // check first event in chunk (is not needed for new reconstructions)
791   if(fCheckFirstEventInChunk){
792     AliAnalysisUtils ut;
793     if(ut.IsFirstEventInChunk(event)) 
794       return -1.;
795     fHistEventStats->Fill(6,gRefMultiplicity); 
796   }
797   // check for pile-up event
798   if(fCheckPileUp){
799     AliAnalysisUtils ut;
800     ut.SetUseMVPlpSelection(kTRUE);
801     ut.SetUseOutOfBunchPileUp(kTRUE);
802     if(ut.IsPileUpEvent(event))
803       return -1.;
804     fHistEventStats->Fill(7,gRefMultiplicity); 
805   }
806
807   // Event trigger bits
808   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
809   if(fUseOfflineTrigger)
810     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
811   
812   if(isSelectedMain) {
813     fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
814  
815     // Event Vertex MC
816     if(gAnalysisLevel == "MC") {
817       if(!event) {
818         AliError("mcEvent not available");
819         return 0x0;
820       }
821       
822       if(mcevent){
823         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
824         if(header){  
825           TArrayF gVertexArray;
826           header->PrimaryVertex(gVertexArray);
827           //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
828           //gVertexArray.At(0),
829           //gVertexArray.At(1),
830           //gVertexArray.At(2));
831           fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
832           if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
833             if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
834               if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
835                 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
836
837                 // get the reference multiplicty or centrality
838                 gRefMultiplicity = GetRefMultiOrCentrality(event);
839
840                 fHistVx->Fill(gVertexArray.At(0));
841                 fHistVy->Fill(gVertexArray.At(1));
842                 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
843                 
844                 // take only events inside centrality class
845                 if(fUseCentrality) {
846                   if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
847                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
848                     return gRefMultiplicity;        
849                   }//centrality class
850                 }
851                 // take events only within the same multiplicity class
852                 else if(fUseMultiplicity) {
853                   if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
854                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
855                     return gRefMultiplicity;
856                   }
857                 }//multiplicity range
858               }//Vz cut
859             }//Vy cut
860           }//Vx cut
861         }//header    
862       }//MC event object
863     }//MC
864     
865     // Event Vertex AOD, ESD, ESDMC
866     else{
867       const AliVVertex *vertex = event->GetPrimaryVertex();
868       
869       if(vertex) {
870         Double32_t fCov[6];
871         vertex->GetCovarianceMatrix(fCov);
872         if(vertex->GetNContributors() > 0) {
873           if(fCov[5] != 0) {
874             fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
875             if(TMath::Abs(vertex->GetX()) < fVxMax) {
876               if(TMath::Abs(vertex->GetY()) < fVyMax) {
877                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
878                   fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
879
880                   // get the reference multiplicty or centrality
881                   gRefMultiplicity = GetRefMultiOrCentrality(event);
882                   
883                   fHistVx->Fill(vertex->GetX());
884                   fHistVy->Fill(vertex->GetY());
885                   fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
886                   
887                   // take only events inside centrality class
888                   if(fUseCentrality) {
889                     if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
890
891                       // centrality weighting (optional for 2011 if central and semicentral triggers are used)
892                       if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
893                         AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
894                         return -1;
895                       }
896                       
897                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
898                       return gRefMultiplicity;          
899                     }//centrality class
900                   }
901                   // take events only within the same multiplicity class
902                   else if(fUseMultiplicity) {
903                     //if(fDebugLevel) 
904                     //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
905                     //fNumberOfAcceptedTracksMax,gRefMultiplicity);
906
907                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
908                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
909                       return gRefMultiplicity;
910                     }
911                   }//multiplicity range
912                 }//Vz cut
913               }//Vy cut
914             }//Vx cut
915           }//proper vertex resolution
916         }//proper number of contributors
917       }//vertex object valid
918     }//triggered event 
919   }//AOD,ESD,ESDMC
920   
921   // in all other cases return -1 (event not accepted)
922   return -1;
923 }
924
925
926 //________________________________________________________________________
927 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
928     // Checks the Event cuts
929     // Fills Event statistics histograms
930   
931   Float_t gCentrality = -1.;
932   Double_t gMultiplicity = -1.;
933   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
934
935
936   // calculate centrality always (not only in centrality mode)
937   if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header  //++++++++++++++
938     AliAODHeader *header = (AliAODHeader*) event->GetHeader();
939     if(header){
940       gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
941
942       // QA for centrality estimators
943       fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
944       fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
945       fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
946       fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
947       fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
948       fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); 
949       fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
950       fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
951       fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
952       fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
953       fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
954       fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
955       fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
956       
957       // Centrality estimator USED   ++++++++++++++++++++++++++++++
958       fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
959       
960       // centrality QA (V0M)
961       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
962       
963       // centrality QA (reference tracks)
964       fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
965       fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
966       fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
967       fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
968       fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
969       fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
970       fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
971       fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
972       fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
973
974     }//AOD header
975   }//AOD
976   
977   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
978     AliCentrality *centrality = event->GetCentrality();
979     gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
980
981     // QA for centrality estimators
982     fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
983     fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
984     fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
985     fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
986     fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
987     fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
988     fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
989     fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
990     fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
991     fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
992     fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
993     fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
994     fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
995     
996     // Centrality estimator USED   ++++++++++++++++++++++++++++++
997     fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
998     
999     // centrality QA (V0M)
1000     fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1001   }//ESD
1002
1003   else if(gAnalysisLevel == "MC"){
1004     Double_t gImpactParameter = 0.;
1005     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1006     if(gMCEvent){
1007       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());      
1008       if(headerH){
1009         gImpactParameter = headerH->ImpactParameter();
1010         gCentrality      = gImpactParameter;
1011       }//MC header
1012     }//MC event cast
1013   }//MC
1014
1015   else{
1016     gCentrality = -1.;
1017   }
1018   
1019   // calculate reference multiplicity always (not only in multiplicity mode)
1020   if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1021     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1022     if(gESDEvent){
1023       gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1024       fHistMultiplicity->Fill(gMultiplicity);
1025     }//AliESDevent cast
1026   }//ESD mode
1027
1028   else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
1029     AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1030     if ((fMultiplicityEstimator == "V0M")||
1031         (fMultiplicityEstimator == "V0A")||
1032         (fMultiplicityEstimator == "V0C") ||
1033         (fMultiplicityEstimator == "TPC")) {
1034       gMultiplicity = GetReferenceMultiplicityFromAOD(event);
1035       if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
1036     }
1037     else {
1038       if(header)
1039         gMultiplicity = header->GetRefMultiplicity();
1040       if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1041     }
1042     fHistMultiplicity->Fill(gMultiplicity);
1043   }//AOD mode
1044   else if(gAnalysisLevel == "MC") {
1045     AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
1046     //Calculating the multiplicity as the number of charged primaries
1047     //within \pm 0.8 in eta and pT > 0.1 GeV/c
1048     for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
1049       AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
1050       if (!track) {
1051         AliError(Form("Could not receive particle %d", iParticle));
1052         continue;
1053       }
1054       
1055       //exclude non stable particles
1056       if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1057       
1058       //++++++++++++++++
1059       if (fMultiplicityEstimator == "V0M") {
1060         if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) 
1061           continue;}
1062       else if (fMultiplicityEstimator == "V0A") {
1063         if(track->Eta() > 5.1 || track->Eta() < 2.8)  continue;}
1064       else if (fMultiplicityEstimator == "V0C") {
1065         if(track->Eta() > -1.7 || track->Eta() < -3.7)  continue;}
1066       else if (fMultiplicityEstimator == "TPC") {
1067         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;
1068         if(track->Pt() < fPtMin || track->Pt() > fPtMax)  continue;
1069       }
1070       else{
1071         if(track->Pt() < fPtMin || track->Pt() > fPtMax)  continue;
1072         if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;
1073       }
1074       //++++++++++++++++
1075       
1076       if(track->Charge() == 0) continue;
1077       
1078       gMultiplicity += 1;
1079     }//loop over primaries
1080     fHistMultiplicity->Fill(gMultiplicity);
1081   }//MC mode
1082   else{
1083     gMultiplicity = -1;
1084   }
1085   
1086
1087   // decide what should be returned only here
1088   Double_t lReturnVal = -100;
1089   if(fEventClass=="Multiplicity"){
1090     lReturnVal = gMultiplicity;
1091   }else if(fEventClass=="Centrality"){
1092     lReturnVal = gCentrality;
1093   }
1094   return lReturnVal;
1095 }
1096
1097 //________________________________________________________________________
1098 Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
1099   //Function that returns the reference multiplicity from AODs (data or reco MC)
1100   //Different ref. mult. implemented: V0M, V0A, V0C, TPC
1101   Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
1102   Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
1103
1104   AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1105   if(!header) {
1106     Printf("ERROR: AOD header not available");
1107     return -999;
1108   }
1109   Int_t gRunNumber = header->GetRunNumber();
1110
1111   // Loop over tracks in event
1112   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1113     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1114     if (!aodTrack) {
1115       AliError(Form("Could not receive track %d", iTracks));
1116       continue;
1117     }
1118     
1119     // AOD track cuts
1120     if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1121             
1122     if(aodTrack->Charge() == 0) continue;
1123     // Kinematics cuts from ESD track cuts
1124     if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax)      continue;
1125     if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax)  continue;
1126     
1127     //=================PID (so far only for electron rejection)==========================//
1128     if(fElectronRejection) {
1129       // get the electron nsigma
1130       Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1131       
1132       // check only for given momentum range
1133       if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
1134         //look only at electron nsigma
1135         if(!fElectronOnlyRejection) {
1136           //Make the decision based on the n-sigma of electrons
1137           if(nSigma < fElectronRejectionNSigma) continue;
1138         }
1139         else {
1140           Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1141           Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1142           Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1143           
1144           //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1145           if(nSigma < fElectronRejectionNSigma
1146              && nSigmaPions   > fElectronRejectionNSigma
1147              && nSigmaKaons   > fElectronRejectionNSigma
1148              && nSigmaProtons > fElectronRejectionNSigma ) continue;
1149         }
1150       }
1151     }//electron rejection
1152     
1153     gRefMultiplicityTPC += 1.0;
1154   }// track loop
1155   
1156   //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
1157   for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
1158     fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
1159     
1160     if(iChannel < 32) 
1161       gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1162     else if(iChannel >= 32) 
1163       gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1164   }//loop over PMTs
1165   
1166   //Equalization of gain
1167   Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1168   if(gFactorA != 0)
1169     gRefMultiplicityVZEROA /= gFactorA;
1170   Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1171   if(gFactorC != 0)
1172     gRefMultiplicityVZEROC /= gFactorC;
1173   if((gFactorA != 0)&&(gFactorC != 0)) 
1174     gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1175   
1176   if(fDebugLevel) 
1177     Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1178
1179   fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1180
1181   if(fMultiplicityEstimator == "TPC") 
1182     gRefMultiplicity = gRefMultiplicityTPC;
1183   else if(fMultiplicityEstimator == "V0M")
1184     gRefMultiplicity = gRefMultiplicityVZERO;
1185   else if(fMultiplicityEstimator == "V0A")
1186     gRefMultiplicity = gRefMultiplicityVZEROA;
1187   else if(fMultiplicityEstimator == "V0C")
1188     gRefMultiplicity = gRefMultiplicityVZEROC;
1189   
1190   return gRefMultiplicity;
1191 }
1192
1193 //________________________________________________________________________
1194 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1195   // Get the event plane
1196
1197   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1198
1199   Float_t gVZEROEventPlane    = -10.;
1200   Float_t gReactionPlane      = -10.;
1201   Double_t qxTot = 0.0, qyTot = 0.0;
1202
1203   //MC: from reaction plane
1204   if(gAnalysisLevel == "MC"){
1205     if(!event) {
1206       AliError("mcEvent not available");
1207       return 0x0;
1208     }
1209
1210     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1211     if(gMCEvent){
1212       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    
1213       if (headerH) {
1214         gReactionPlane = headerH->ReactionPlaneAngle();
1215         //gReactionPlane *= TMath::RadToDeg();
1216       }//MC header
1217     }//MC event cast
1218   }//MC
1219   
1220   // AOD,ESD,ESDMC: from VZERO Event Plane
1221   else{
1222    
1223     AliEventplane *ep = event->GetEventplane();
1224     if(ep){ 
1225       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1226       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1227       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1228       gReactionPlane = gVZEROEventPlane;
1229     }
1230   }//AOD,ESD,ESDMC
1231
1232   return gReactionPlane;
1233 }
1234
1235 //________________________________________________________________________
1236 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, 
1237                                                                 Double_t vPhi, 
1238                                                                 Double_t vPt, 
1239                                                                 Short_t vCharge, 
1240                                                                 Double_t gCentrality) {
1241   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) 
1242
1243   Double_t correction = 1.;
1244   Int_t gCentralityInt = -1;
1245
1246   for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1247     if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1248       gCentralityInt = i;
1249       break;
1250     }
1251   }  
1252
1253   // centrality not in array --> no correction
1254   if(gCentralityInt < 0){
1255     correction = 1.;
1256   }
1257   else{
1258     
1259     //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1260     
1261     if(fHistCorrectionPlus[gCentralityInt]){
1262       if (vCharge > 0) {
1263         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1264         //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);  
1265       }
1266       if (vCharge < 0) {
1267         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1268         //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt); 
1269       }
1270     }
1271     else {
1272       correction = 1.;
1273     }
1274   }//centrality in array
1275   
1276   if (correction == 0.) { 
1277     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); 
1278     correction = 1.; 
1279   } 
1280   
1281   return correction;
1282 }
1283
1284 //________________________________________________________________________
1285 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
1286   // Returns TObjArray with tracks after all track cuts (only for AOD!)
1287   // Fills QA histograms
1288
1289   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1290
1291   //output TObjArray holding all good tracks
1292   TObjArray* tracksAccepted = new TObjArray;
1293   tracksAccepted->SetOwner(kTRUE);
1294
1295   Short_t vCharge;
1296   Double_t vEta;
1297   Double_t vY;
1298   Double_t vPhi;
1299   Double_t vPt;
1300
1301
1302   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1303     // Loop over tracks in event
1304     
1305     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1306       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1307       if (!aodTrack) {
1308         AliError(Form("Could not receive track %d", iTracks));
1309         continue;
1310       }
1311       
1312       // AOD track cuts
1313       
1314       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
1315       // take only TPC only tracks 
1316       //fHistTrackStats->Fill(aodTrack->GetFilterMap());
1317       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1318         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1319       }
1320
1321       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1322
1323
1324       // additional check on kPrimary flag
1325       if(fCheckPrimaryFlagAOD){
1326         if(aodTrack->GetType() != AliAODTrack::kPrimary)
1327           continue;
1328       }
1329
1330      
1331       vCharge = aodTrack->Charge();
1332       vEta    = aodTrack->Eta();
1333       vY      = aodTrack->Y();
1334       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1335       vPt     = aodTrack->Pt();
1336       
1337       //===========================PID (so far only for electron rejection)===============================//                
1338       if(fElectronRejection) {
1339
1340         // get the electron nsigma
1341         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1342         
1343         //Fill QA before the PID
1344         fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1345         fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1346         //end of QA-before pid
1347         
1348         // check only for given momentum range
1349         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1350                   
1351           //look only at electron nsigma
1352           if(!fElectronOnlyRejection){
1353             
1354             //Make the decision based on the n-sigma of electrons
1355             if(nSigma < fElectronRejectionNSigma) continue;
1356           }
1357           else{
1358             
1359             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1360             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1361             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1362             
1363             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1364             if(nSigma < fElectronRejectionNSigma
1365                && nSigmaPions   > fElectronRejectionNSigma
1366                && nSigmaKaons   > fElectronRejectionNSigma
1367                && nSigmaProtons > fElectronRejectionNSigma ) continue;
1368           }
1369         }
1370   
1371         //Fill QA after the PID
1372         fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1373         fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1374         
1375       }
1376       //===========================end of PID (so far only for electron rejection)===============================//
1377
1378       //+++++++++++++++++++++++++++++//
1379       //===========================PID===============================//             
1380       if(fUsePID) {
1381         Double_t prob[AliPID::kSPECIES]={0.};
1382         Double_t probTPC[AliPID::kSPECIES]={0.};
1383         Double_t probTOF[AliPID::kSPECIES]={0.};
1384         Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1385         
1386         Double_t nSigma = 0.;
1387         Double_t nSigmaTPC = 0.; //++++
1388         Double_t nSigmaTOF = 0.; //++++
1389         Double_t nSigmaTPCTOF = 0.; //++++
1390         UInt_t detUsedTPC = 0;
1391         UInt_t detUsedTOF = 0;
1392         UInt_t detUsedTPCTOF = 0;
1393         
1394         nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1395         nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1396         nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
1397
1398         //Decide what detector configuration we want to use
1399         switch(fPidDetectorConfig) {
1400         case kTPCpid:
1401           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1402           nSigma = nSigmaTPC;
1403           detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1404           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1405             prob[iSpecies] = probTPC[iSpecies];
1406           break;
1407         case kTOFpid:
1408           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1409           nSigma = nSigmaTOF;
1410           detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1411           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1412             prob[iSpecies] = probTOF[iSpecies];
1413           break;
1414         case kTPCTOF:
1415           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1416           nSigma = nSigmaTPCTOF;
1417           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
1418           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1419             prob[iSpecies] = probTPCTOF[iSpecies];
1420           break;
1421         default:
1422           break;
1423         }//end switch: define detector mask
1424         
1425         //Filling the PID QA
1426         Double_t tofTime = -999., length = 999., tof = -999.;
1427         Double_t c = TMath::C()*1.E-9;// m/ns
1428         Double_t beta = -999.;
1429         if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
1430              (aodTrack->IsOn(AliAODTrack::kTIME))  ) { 
1431           tofTime = aodTrack->GetTOFsignal();//in ps
1432           length = aodTrack->GetIntegratedLength();
1433           tof = tofTime*1E-3; // ns     
1434           
1435           if (tof <= 0) {
1436             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1437             continue;
1438           }
1439           if (length <= 0){
1440             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1441             continue;
1442           }
1443           
1444           length = length*0.01; // in meters
1445           tof = tof*c;
1446           beta = length/tof;
1447           
1448           fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1449           fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1450           fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1451         }//TOF signal 
1452         
1453         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1454         fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); 
1455         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
1456         
1457         fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1458
1459         //combined TPC&TOF
1460         fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++      
1461         fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
1462         Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
1463         
1464         //end of QA-before pid
1465         
1466         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1467           //Make the decision based on the n-sigma
1468           if(fUsePIDnSigma) {
1469             if(nSigma > fPIDNSigma) continue;  
1470             
1471             fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
1472             fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
1473             fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
1474           }
1475           //Make the decision based on the bayesian
1476           else if(fUsePIDPropabilities) {
1477             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1478             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      
1479           
1480             fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1481             fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); 
1482             fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1483           
1484           }
1485             
1486           //Fill QA after the PID
1487           fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1488           fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1489           fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++      
1490         }
1491       }
1492       //===========================PID===============================//
1493       //+++++++++++++++++++++++++++++//
1494
1495
1496       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)
1497       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
1498       
1499       
1500       // Kinematics cuts from ESD track cuts
1501       if( vPt < fPtMin || vPt > fPtMax)      continue;
1502       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1503       
1504       // Extra DCA cuts (for systematic studies [!= -1])
1505       if( fDCAxyCut != -1 && fDCAzCut != -1){
1506         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1507           continue;  // 2D cut
1508         }
1509       }
1510       
1511       // Extra TPC cuts (for systematic studies [!= -1])
1512       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1513         continue;
1514       }
1515       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1516         continue;
1517       }
1518       
1519       // fill QA histograms
1520       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1521       fHistDCA->Fill(dcaZ,dcaXY);
1522       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1523       fHistPt->Fill(vPt,gCentrality);
1524       fHistEta->Fill(vEta,gCentrality);
1525       fHistRapidity->Fill(vY,gCentrality);
1526       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1527       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1528       fHistPhi->Fill(vPhi,gCentrality);
1529       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  
1530       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1531       
1532       //=======================================correction
1533       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1534       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1535       
1536       // add the track to the TObjArray
1537       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1538     }//track loop
1539   }// AOD analysis
1540
1541   //==============================================================================================================
1542   else if(gAnalysisLevel == "MCAOD") {
1543     
1544     AliMCEvent* mcEvent = MCEvent();
1545     if (!mcEvent) {
1546       AliError("ERROR: Could not retrieve MC event");
1547     }
1548     else{
1549       
1550       for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1551         AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); 
1552         if (!aodTrack) {
1553           AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1554           continue;
1555         }
1556         
1557         if(!aodTrack->IsPhysicalPrimary()) continue;   
1558         
1559         vCharge = aodTrack->Charge();
1560         vEta    = aodTrack->Eta();
1561         vY      = aodTrack->Y();
1562         vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1563         vPt     = aodTrack->Pt();
1564         
1565         // Kinematics cuts from ESD track cuts
1566         if( vPt < fPtMin || vPt > fPtMax)      continue;
1567         if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1568         
1569         // Remove neutral tracks
1570         if( vCharge == 0 ) continue;
1571         
1572         //Exclude resonances
1573         if(fExcludeResonancesInMC) {
1574           
1575           Bool_t kExcludeParticle = kFALSE;
1576           Int_t gMotherIndex = aodTrack->GetMother();
1577           if(gMotherIndex != -1) {
1578             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1579             if(motherTrack) {
1580               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1581               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1582               //if(pdgCodeOfMother == 113) {
1583               if(pdgCodeOfMother == 113  // rho0
1584                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1585                  // || pdgCodeOfMother == 221  // eta
1586                  // || pdgCodeOfMother == 331  // eta'
1587                  // || pdgCodeOfMother == 223  // omega
1588                  // || pdgCodeOfMother == 333  // phi
1589                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0
1590                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*
1591                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*
1592                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1593                  || pdgCodeOfMother == 111  // pi0 Dalitz
1594                  || pdgCodeOfMother == 22  // photon
1595                  ) {
1596                 kExcludeParticle = kTRUE;
1597               }
1598             }
1599           }
1600           
1601           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1602           if(kExcludeParticle) continue;
1603         }
1604
1605         //Exclude electrons with PDG
1606         if(fExcludeElectronsInMC) {
1607           
1608           if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1609           
1610         }
1611         
1612         // fill QA histograms
1613         fHistPt->Fill(vPt,gCentrality);
1614         fHistEta->Fill(vEta,gCentrality);
1615         fHistRapidity->Fill(vY,gCentrality);
1616         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1617         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1618         fHistPhi->Fill(vPhi,gCentrality);
1619         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                
1620         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1621         
1622         //=======================================correction
1623         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1624         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   
1625         
1626         // add the track to the TObjArray
1627         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1628       }//aodTracks
1629     }//MC event
1630   }//MCAOD
1631   //==============================================================================================================
1632
1633   //==============================================================================================================
1634   else if(gAnalysisLevel == "MCAODrec") {
1635     
1636     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1637     if (!fAOD) {
1638       printf("ERROR: fAOD not available\n");
1639       return;
1640       }*/
1641     
1642     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); 
1643     if (!fArrayMC) { 
1644        AliError("No array of MC particles found !!!");
1645     }
1646
1647     AliMCEvent* mcEvent = MCEvent();
1648     if (!mcEvent) {
1649        AliError("ERROR: Could not retrieve MC event");
1650        return tracksAccepted;
1651     }
1652      
1653     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1654       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1655       if (!aodTrack) {
1656         AliError(Form("Could not receive track %d", iTracks));
1657         continue;
1658       }
1659       
1660       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1661         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1662       }
1663       if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1664       
1665       vCharge = aodTrack->Charge();
1666       vEta    = aodTrack->Eta();
1667       vY      = aodTrack->Y();
1668       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();
1669       vPt     = aodTrack->Pt();
1670       
1671       //===========================use MC information for Kinematics===============================//               
1672       if(fUseMCforKinematics){
1673
1674         Int_t label = TMath::Abs(aodTrack->GetLabel());
1675         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1676
1677         if(AODmcTrack){
1678           vCharge = AODmcTrack->Charge();
1679           vEta    = AODmcTrack->Eta();
1680           vY      = AODmcTrack->Y();
1681           vPhi    = AODmcTrack->Phi();// * TMath::RadToDeg();
1682           vPt     = AODmcTrack->Pt();
1683         }
1684         else{
1685           AliDebug(1, "no MC particle for this track"); 
1686           continue;
1687         }
1688       }
1689       //===========================end of use MC information for Kinematics========================//               
1690
1691
1692       //===========================PID (so far only for electron rejection)===============================//                
1693       if(fElectronRejection) {
1694
1695         // get the electron nsigma
1696         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1697
1698         //Fill QA before the PID
1699         fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1700         fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1701         //end of QA-before pid
1702         
1703         // check only for given momentum range
1704         if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1705                   
1706           //look only at electron nsigma
1707           if(!fElectronOnlyRejection){
1708             
1709             //Make the decision based on the n-sigma of electrons
1710             if(nSigma < fElectronRejectionNSigma) continue;
1711           }
1712           else{
1713             
1714             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1715             Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1716             Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1717             
1718             //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1719             if(nSigma < fElectronRejectionNSigma
1720                && nSigmaPions   > fElectronRejectionNSigma
1721                && nSigmaKaons   > fElectronRejectionNSigma
1722                && nSigmaProtons > fElectronRejectionNSigma ) continue;
1723           }
1724         }
1725   
1726         //Fill QA after the PID
1727         fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1728         fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); 
1729         
1730       }
1731       //===========================end of PID (so far only for electron rejection)===============================//
1732       
1733       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)
1734       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
1735             
1736       // Kinematics cuts from ESD track cuts
1737       if( vPt < fPtMin || vPt > fPtMax)      continue;
1738       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1739       
1740       // Extra DCA cuts (for systematic studies [!= -1])
1741       if( fDCAxyCut != -1 && fDCAzCut != -1){
1742         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1743           continue;  // 2D cut
1744         }
1745       }
1746       
1747       // Extra TPC cuts (for systematic studies [!= -1])
1748       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1749         continue;
1750       }
1751       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1752         continue;
1753       }
1754
1755       //Exclude resonances
1756       if(fExcludeResonancesInMC) {
1757         
1758         Bool_t kExcludeParticle = kFALSE;
1759
1760         Int_t label = TMath::Abs(aodTrack->GetLabel());
1761         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1762       
1763         if (AODmcTrack){ 
1764           //if (AODmcTrack->IsPhysicalPrimary()){
1765           
1766           Int_t gMotherIndex = AODmcTrack->GetMother();
1767           if(gMotherIndex != -1) {
1768             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1769             if(motherTrack) {
1770               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1771               if(pdgCodeOfMother == 113  // rho0
1772                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1773                  // || pdgCodeOfMother == 221  // eta
1774                  // || pdgCodeOfMother == 331  // eta'
1775                  // || pdgCodeOfMother == 223  // omega
1776                  // || pdgCodeOfMother == 333  // phi
1777                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0
1778                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*
1779                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*
1780                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1781                  || pdgCodeOfMother == 111  // pi0 Dalitz
1782                  || pdgCodeOfMother == 22   // photon
1783                  ) {
1784                 kExcludeParticle = kTRUE;
1785               }
1786             }
1787           }
1788         }       
1789         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1790         if(kExcludeParticle) continue;
1791       }
1792
1793       //Exclude electrons with PDG
1794       if(fExcludeElectronsInMC) {
1795         
1796         Int_t label = TMath::Abs(aodTrack->GetLabel());
1797         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1798         
1799         if (AODmcTrack){ 
1800           if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1801         }
1802       }
1803       
1804       // fill QA histograms
1805       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1806       fHistDCA->Fill(dcaZ,dcaXY);
1807       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1808       fHistPt->Fill(vPt,gCentrality);
1809       fHistEta->Fill(vEta,gCentrality);
1810       fHistRapidity->Fill(vY,gCentrality);
1811       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1812       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1813       fHistPhi->Fill(vPhi,gCentrality);
1814       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  
1815       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1816       
1817       //=======================================correction
1818       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
1819       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1820       
1821       // add the track to the TObjArray
1822       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  
1823     }//track loop
1824   }//MCAODrec
1825   //==============================================================================================================
1826
1827   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1828
1829     AliESDtrack *trackTPC   = NULL;
1830     AliMCParticle *mcTrack   = NULL;
1831
1832     AliMCEvent*  mcEvent     = NULL;
1833     
1834     // for MC ESDs use also MC information
1835     if(gAnalysisLevel == "MCESD"){
1836       mcEvent = MCEvent(); 
1837       if (!mcEvent) {
1838         AliError("mcEvent not available");
1839         return tracksAccepted;
1840       }
1841     }
1842     
1843     // Loop over tracks in event
1844     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1845       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
1846       if (!track) {
1847         AliError(Form("Could not receive track %d", iTracks));
1848         continue;
1849       } 
1850
1851       // for MC ESDs use also MC information --> MC track not used further???
1852       if(gAnalysisLevel == "MCESD"){
1853         Int_t label = TMath::Abs(track->GetLabel());
1854         if(label > mcEvent->GetNumberOfTracks()) continue;
1855         if(label > mcEvent->GetNumberOfPrimaries()) continue;
1856         
1857         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1858         if(!mcTrack) continue;
1859       }
1860
1861       // take only TPC only tracks
1862       trackTPC   = new AliESDtrack();
1863       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1864       
1865       //ESD track cuts
1866       if(fESDtrackCuts) 
1867         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1868       
1869       // fill QA histograms
1870       Float_t b[2];
1871       Float_t bCov[3];
1872       trackTPC->GetImpactParameters(b,bCov);
1873       if (bCov[0]<=0 || bCov[2]<=0) {
1874         AliDebug(1, "Estimated b resolution lower or equal zero!");
1875         bCov[0]=0; bCov[2]=0;
1876       }
1877       
1878       Int_t nClustersTPC = -1;
1879       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone
1880       //nClustersTPC = track->GetTPCclusters(0);   // global track
1881       Float_t chi2PerClusterTPC = -1;
1882       if (nClustersTPC!=0) {
1883         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone
1884         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track
1885       }
1886       
1887       //===========================PID===============================//             
1888       if(fUsePID) {
1889         Double_t prob[AliPID::kSPECIES]={0.};
1890         Double_t probTPC[AliPID::kSPECIES]={0.};
1891         Double_t probTOF[AliPID::kSPECIES]={0.};
1892         Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1893         
1894         Double_t nSigma = 0.;
1895         UInt_t detUsedTPC = 0;
1896         UInt_t detUsedTOF = 0;
1897         UInt_t detUsedTPCTOF = 0;
1898         
1899         //Decide what detector configuration we want to use
1900         switch(fPidDetectorConfig) {
1901         case kTPCpid:
1902           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1903           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
1904           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
1905           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1906             prob[iSpecies] = probTPC[iSpecies];
1907           break;
1908         case kTOFpid:
1909           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1910           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
1911           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
1912           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1913             prob[iSpecies] = probTOF[iSpecies];
1914           break;
1915         case kTPCTOF:
1916           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1917           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
1918           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1919             prob[iSpecies] = probTPCTOF[iSpecies];
1920           break;
1921         default:
1922           break;
1923         }//end switch: define detector mask
1924         
1925         //Filling the PID QA
1926         Double_t tofTime = -999., length = 999., tof = -999.;
1927         Double_t c = TMath::C()*1.E-9;// m/ns
1928         Double_t beta = -999.;
1929         Double_t  nSigmaTOFForParticleOfInterest = -999.;
1930         if ( (track->IsOn(AliESDtrack::kTOFin)) &&
1931              (track->IsOn(AliESDtrack::kTIME))  ) { 
1932           tofTime = track->GetTOFsignal();//in ps
1933           length = track->GetIntegratedLength();
1934           tof = tofTime*1E-3; // ns     
1935           
1936           if (tof <= 0) {
1937             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1938             continue;
1939           }
1940           if (length <= 0){
1941             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1942             continue;
1943           }
1944           
1945           length = length*0.01; // in meters
1946           tof = tof*c;
1947           beta = length/tof;
1948           
1949           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
1950           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
1951           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
1952           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
1953         }//TOF signal 
1954         
1955         
1956         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
1957         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
1958         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); 
1959         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); 
1960         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
1961         //end of QA-before pid
1962         
1963         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1964           //Make the decision based on the n-sigma
1965           if(fUsePIDnSigma) {
1966             if(nSigma > fPIDNSigma) continue;}
1967           
1968           //Make the decision based on the bayesian
1969           else if(fUsePIDPropabilities) {
1970             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1971             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      
1972           }
1973           
1974           //Fill QA after the PID
1975           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
1976           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
1977           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
1978           
1979           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
1980           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); 
1981           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
1982           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); 
1983         }
1984       }
1985       //===========================PID===============================//
1986       vCharge = trackTPC->Charge();
1987       vY      = trackTPC->Y();
1988       vEta    = trackTPC->Eta();
1989       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();
1990       vPt     = trackTPC->Pt();
1991       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
1992       fHistDCA->Fill(b[1],b[0]);
1993       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
1994       fHistPt->Fill(vPt,gCentrality);
1995       fHistEta->Fill(vEta,gCentrality);
1996       fHistPhi->Fill(vPhi,gCentrality);
1997       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1998       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1999       fHistRapidity->Fill(vY,gCentrality);
2000       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2001       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2002       
2003       //=======================================correction
2004       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
2005       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2006
2007       // add the track to the TObjArray
2008       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   
2009
2010       delete trackTPC;
2011     }//track loop
2012   }// ESD analysis
2013
2014   else if(gAnalysisLevel == "MC"){
2015     if(!event) {
2016       AliError("mcEvent not available");
2017       return 0x0;
2018     }
2019
2020     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2021     if(gMCEvent) {
2022       // Loop over tracks in event
2023       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
2024         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
2025         if (!track) {
2026           AliError(Form("Could not receive particle %d", iTracks));
2027           continue;
2028         }
2029         
2030         //exclude non stable particles
2031         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2032         
2033         vCharge = track->Charge();
2034         vEta    = track->Eta();
2035         vPt     = track->Pt();
2036         vY      = track->Y();
2037         
2038         if( vPt < fPtMin || vPt > fPtMax)      
2039           continue;
2040         if (!fUsePID) {
2041           if( vEta < fEtaMin || vEta > fEtaMax)  continue;
2042         }
2043         else if (fUsePID){
2044           if( vY < fEtaMin || vY > fEtaMax)  continue;
2045         }
2046
2047         // Remove neutral tracks
2048         if( vCharge == 0 ) continue;
2049         
2050         //analyze one set of particles
2051         if(fUseMCPdgCode) {
2052           TParticle *particle = track->Particle();
2053           if(!particle) continue;
2054           
2055           Int_t gPdgCode = particle->GetPdgCode();
2056           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) 
2057             continue;
2058         }
2059         
2060         //Use the acceptance parameterization
2061         if(fAcceptanceParameterization) {
2062           Double_t gRandomNumber = gRandom->Rndm();
2063           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) 
2064             continue;
2065         }
2066         
2067         //Exclude resonances
2068         if(fExcludeResonancesInMC) {
2069           TParticle *particle = track->Particle();
2070           if(!particle) continue;
2071           
2072           Bool_t kExcludeParticle = kFALSE;
2073           Int_t gMotherIndex = particle->GetFirstMother();
2074           if(gMotherIndex != -1) {
2075             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2076             if(motherTrack) {
2077               TParticle *motherParticle = motherTrack->Particle();
2078               if(motherParticle) {
2079                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2080                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2081                 if(pdgCodeOfMother == 113  // rho0
2082                    || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2083                    // || pdgCodeOfMother == 221  // eta
2084                    // || pdgCodeOfMother == 331  // eta'
2085                    // || pdgCodeOfMother == 223  // omega
2086                    // || pdgCodeOfMother == 333  // phi
2087                    || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0
2088                    // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*
2089                    // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*
2090                    || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2091                    || pdgCodeOfMother == 111  // pi0 Dalitz
2092                    ) {
2093                   kExcludeParticle = kTRUE;
2094                 }
2095               }
2096             }
2097           }
2098           
2099           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2100           if(kExcludeParticle) continue;
2101         }
2102
2103         //Exclude electrons with PDG
2104         if(fExcludeElectronsInMC) {
2105           
2106           TParticle *particle = track->Particle();
2107           
2108           if (particle){ 
2109             if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2110           }
2111         }
2112       
2113         vPhi    = track->Phi();
2114         //Printf("phi (before): %lf",vPhi);
2115         
2116         fHistPt->Fill(vPt,gCentrality);
2117         fHistEta->Fill(vEta,gCentrality);
2118         fHistPhi->Fill(vPhi,gCentrality);
2119         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2120         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2121         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2122         fHistRapidity->Fill(vY,gCentrality);
2123         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2124         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2125         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2126         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2127         
2128         //Flow after burner
2129         if(fUseFlowAfterBurner) {
2130           Double_t precisionPhi = 0.001;
2131           Int_t maxNumberOfIterations = 100;
2132           
2133           Double_t phi0 = vPhi;
2134           Double_t gV2 = fDifferentialV2->Eval(vPt);
2135           
2136           for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2137             Double_t phiprev = vPhi;
2138             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2139             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); 
2140             vPhi -= fl/fp;
2141             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2142           }
2143           //Printf("phi (after): %lf\n",vPhi);
2144           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2145           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2146           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2147           
2148           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2149           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2150           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2151           
2152         }
2153         
2154         //vPhi *= TMath::RadToDeg();
2155         
2156       //=======================================correction
2157       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  
2158       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2159
2160         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); 
2161       } //track loop
2162     }//MC event object
2163   }//MC
2164   
2165   return tracksAccepted;  
2166 }
2167
2168 //________________________________________________________________________
2169 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2170   // Clones TObjArray and returns it with tracks after shuffling the charges
2171
2172   TObjArray* tracksShuffled = new TObjArray;
2173   tracksShuffled->SetOwner(kTRUE);
2174
2175   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks 
2176
2177   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2178   {
2179     AliVParticle* track = (AliVParticle*) tracks->At(i);
2180     chargeVector->push_back(track->Charge());
2181   }  
2182  
2183   random_shuffle(chargeVector->begin(), chargeVector->end());
2184   
2185   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2186     AliVParticle* track = (AliVParticle*) tracks->At(i);
2187     //==============================correction
2188     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2189     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2190     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2191   }
2192
2193   delete chargeVector;
2194    
2195   return tracksShuffled;
2196 }
2197
2198 //________________________________________________________________________
2199 void  AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2200                                                     const char* lhcPeriod) {
2201   //Function to setup the VZERO gain equalization
2202     //============Get the equilization map============//
2203   TFile *calibrationFile = TFile::Open(filename);
2204   if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2205     Printf("No calibration file found!!!");
2206     return;
2207   }
2208
2209   TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2210   if(!list) {
2211     Printf("Calibration TList not found!!!");
2212     return;
2213   }
2214
2215   fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2216   if(!fHistVZEROAGainEqualizationMap) {
2217     Printf("VZERO-A calibration object not found!!!");
2218     return;
2219   }
2220   fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2221   if(!fHistVZEROCGainEqualizationMap) {
2222     Printf("VZERO-C calibration object not found!!!");
2223     return;
2224   }
2225
2226   fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2227   if(!fHistVZEROChannelGainEqualizationMap) {
2228     Printf("VZERO channel calibration object not found!!!");
2229     return;
2230   }
2231 }
2232
2233 //________________________________________________________________________
2234 Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, 
2235                                                             Int_t channel) {
2236   //
2237   if(!fHistVZEROAGainEqualizationMap) return 1.0;
2238
2239   for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2240     Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2241     if(gRunNumber == run)
2242       return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2243   }
2244
2245   return 1.0;
2246 }
2247
2248 //________________________________________________________________________
2249 Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, 
2250                                                      const char* side) {
2251   //
2252   if(!fHistVZEROAGainEqualizationMap) return 1.0;
2253
2254   TString gVZEROSide = side;
2255   for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2256     Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2257     //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2258     if(gRunNumber == run) {
2259       if(gVZEROSide == "A") 
2260         return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2261       else if(gVZEROSide == "C") 
2262         return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2263     }
2264   }
2265
2266   return 1.0;
2267 }
2268
2269 //____________________________________________________________________
2270 Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2271 {
2272   // copied from AliAnalysisTaskPhiCorrelations
2273   //
2274   // rejects "randomly" events such that the centrality gets flat
2275   // uses fCentralityWeights histogram
2276
2277   // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2278   
2279   Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2280   Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2281   
2282   Bool_t result = kFALSE;
2283   if (centralityDigits < weight) 
2284     result = kTRUE;
2285   
2286   AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2287   
2288   return result;
2289 }
2290
2291 //________________________________________________________________________
2292 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){
2293   //Printf("END BF");
2294
2295   if (!fBalance) {
2296     AliError("fBalance not available");
2297     return;
2298   }  
2299   if(fRunShuffling) {
2300     if (!fShuffledBalance) {
2301       AliError("fShuffledBalance not available");
2302       return;
2303     }
2304   }
2305
2306 }
2307
2308 //________________________________________________________________________
2309 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2310   // Draw result to the screen
2311   // Called once at the end of the query
2312
2313   // not implemented ...
2314
2315 }
2316