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