fixing psi in MC from the header
[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 "TH2D.h"                  \r
9 #include "TH3D.h"\r
10 #include "TArrayF.h"\r
11 #include "TF1.h"\r
12 #include "TRandom.h"\r
13 \r
14 #include "AliAnalysisTaskSE.h"\r
15 #include "AliAnalysisManager.h"\r
16 \r
17 #include "AliESDVertex.h"\r
18 #include "AliESDEvent.h"\r
19 #include "AliESDInputHandler.h"\r
20 #include "AliAODEvent.h"\r
21 #include "AliAODTrack.h"\r
22 #include "AliAODInputHandler.h"\r
23 #include "AliGenEventHeader.h"\r
24 #include "AliGenHijingEventHeader.h"\r
25 #include "AliMCEventHandler.h"\r
26 #include "AliMCEvent.h"\r
27 #include "AliStack.h"\r
28 #include "AliESDtrackCuts.h"\r
29 #include "AliEventplane.h"\r
30 \r
31 #include "AliPID.h"                \r
32 #include "AliPIDResponse.h"        \r
33 #include "AliPIDCombined.h"        \r
34 \r
35 #include "AliAnalysisTaskBFPsi.h"\r
36 #include "AliBalancePsi.h"\r
37 \r
38 \r
39 // Analysis task for the BF vs Psi code\r
40 // Authors: Panos.Christakoglou@nikhef.nl\r
41 \r
42 ClassImp(AliAnalysisTaskBFPsi)\r
43 \r
44 //________________________________________________________________________\r
45 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
46 : AliAnalysisTaskSE(name), \r
47   fBalance(0),\r
48   fRunShuffling(kFALSE),\r
49   fShuffledBalance(0),\r
50   fList(0),\r
51   fListBF(0),\r
52   fListBFS(0),\r
53   fHistListPIDQA(0),\r
54   fHistEventStats(0),\r
55   fHistCentStats(0),\r
56   fHistTriggerStats(0),\r
57   fHistTrackStats(0),\r
58   fHistVx(0),\r
59   fHistVy(0),\r
60   fHistVz(0),\r
61   fHistEventPlane(0),\r
62   fHistClus(0),\r
63   fHistDCA(0),\r
64   fHistChi2(0),\r
65   fHistPt(0),\r
66   fHistEta(0),\r
67   fHistPhi(0),\r
68   fHistPhiBefore(0),\r
69   fHistPhiAfter(0),\r
70   fHistV0M(0),\r
71   fHistRefTracks(0),\r
72   fHistdEdxVsPTPCbeforePID(NULL),\r
73   fHistBetavsPTOFbeforePID(NULL), \r
74   fHistProbTPCvsPtbeforePID(NULL), \r
75   fHistProbTOFvsPtbeforePID(NULL), \r
76   fHistProbTPCTOFvsPtbeforePID(NULL),\r
77   fHistNSigmaTPCvsPtbeforePID(NULL), \r
78   fHistNSigmaTOFvsPtbeforePID(NULL), \r
79   fHistdEdxVsPTPCafterPID(NULL),\r
80   fHistBetavsPTOFafterPID(NULL), \r
81   fHistProbTPCvsPtafterPID(NULL), \r
82   fHistProbTOFvsPtafterPID(NULL), \r
83   fHistProbTPCTOFvsPtafterPID(NULL),\r
84   fHistNSigmaTPCvsPtafterPID(NULL), \r
85   fHistNSigmaTOFvsPtafterPID(NULL),  \r
86   fPIDResponse(0x0),\r
87   fPIDCombined(0x0),\r
88   fParticleOfInterest(kPion),\r
89   fPidDetectorConfig(kTPCTOF),\r
90   fUsePID(kFALSE),\r
91   fUsePIDnSigma(kTRUE),\r
92   fUsePIDPropabilities(kFALSE), \r
93   fPIDNSigma(3.),\r
94   fMinAcceptedPIDProbability(0.8),\r
95   fESDtrackCuts(0),\r
96   fCentralityEstimator("V0M"),\r
97   fUseCentrality(kFALSE),\r
98   fCentralityPercentileMin(0.), \r
99   fCentralityPercentileMax(5.),\r
100   fImpactParameterMin(0.),\r
101   fImpactParameterMax(20.),\r
102   fUseMultiplicity(kFALSE),\r
103   fNumberOfAcceptedTracksMin(0),\r
104   fNumberOfAcceptedTracksMax(10000),\r
105   fHistNumberOfAcceptedTracks(0),\r
106   fUseOfflineTrigger(kFALSE),\r
107   fVxMax(0.3),\r
108   fVyMax(0.3),\r
109   fVzMax(10.),\r
110   nAODtrackCutBit(128),\r
111   fPtMin(0.3),\r
112   fPtMax(1.5),\r
113   fEtaMin(-0.8),\r
114   fEtaMax(-0.8),\r
115   fDCAxyCut(-1),\r
116   fDCAzCut(-1),\r
117   fTPCchi2Cut(-1),\r
118   fNClustersTPCCut(-1),\r
119   fAcceptanceParameterization(0),\r
120   fDifferentialV2(0),\r
121   fUseFlowAfterBurner(kFALSE),\r
122   fExcludeResonancesInMC(kFALSE),\r
123   fUseMCPdgCode(kFALSE),\r
124   fPDGCodeToBeAnalyzed(-1) {\r
125   // Constructor\r
126   // Define input and output slots here\r
127   // Input slot #0 works with a TChain\r
128   DefineInput(0, TChain::Class());\r
129   // Output slot #0 writes into a TH1 container\r
130   DefineOutput(1, TList::Class());\r
131   DefineOutput(2, TList::Class());\r
132   DefineOutput(3, TList::Class());\r
133   DefineOutput(4, TList::Class());\r
134 }\r
135 \r
136 //________________________________________________________________________\r
137 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
138 \r
139   // delete fBalance; \r
140   // delete fShuffledBalance; \r
141   // delete fList;\r
142   // delete fListBF; \r
143   // delete fListBFS;\r
144 \r
145   // delete fHistEventStats; \r
146   // delete fHistTrackStats; \r
147   // delete fHistVx; \r
148   // delete fHistVy; \r
149   // delete fHistVz; \r
150 \r
151   // delete fHistClus;\r
152   // delete fHistDCA;\r
153   // delete fHistChi2;\r
154   // delete fHistPt;\r
155   // delete fHistEta;\r
156   // delete fHistPhi;\r
157   // delete fHistV0M;\r
158 }\r
159 \r
160 //________________________________________________________________________\r
161 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
162   // Create histograms\r
163   // Called once\r
164   if(!fBalance) {\r
165     fBalance = new AliBalancePsi();\r
166     fBalance->SetAnalysisLevel("ESD");\r
167     //fBalance->SetNumberOfBins(-1,16);\r
168     fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
169   }\r
170   if(fRunShuffling) {\r
171     if(!fShuffledBalance) {\r
172       fShuffledBalance = new AliBalancePsi();\r
173       fShuffledBalance->SetAnalysisLevel("ESD");\r
174       //fShuffledBalance->SetNumberOfBins(-1,16);\r
175       fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
176     }\r
177   }\r
178 \r
179   //QA list\r
180   fList = new TList();\r
181   fList->SetName("listQA");\r
182   fList->SetOwner();\r
183 \r
184   //Balance Function list\r
185   fListBF = new TList();\r
186   fListBF->SetName("listBF");\r
187   fListBF->SetOwner();\r
188 \r
189   if(fRunShuffling) {\r
190     fListBFS = new TList();\r
191     fListBFS->SetName("listBFShuffled");\r
192     fListBFS->SetOwner();\r
193   }\r
194 \r
195   //PID QA list\r
196   if(fUsePID) {\r
197     fHistListPIDQA = new TList();\r
198     fHistListPIDQA->SetName("listQAPID");\r
199     fHistListPIDQA->SetOwner();\r
200   }\r
201 \r
202   //Event stats.\r
203   TString gCutName[4] = {"Total","Offline trigger",\r
204                          "Vertex","Analyzed"};\r
205   fHistEventStats = new TH1F("fHistEventStats",\r
206                              "Event statistics;;N_{events}",\r
207                              4,0.5,4.5);\r
208   for(Int_t i = 1; i <= 4; i++)\r
209     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
210   fList->Add(fHistEventStats);\r
211 \r
212   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
213   fHistCentStats = new TH2F("fHistCentStats",\r
214                              "Centrality statistics;;Cent percentile",\r
215                             9,-0.5,8.5,220,-5,105);\r
216   for(Int_t i = 1; i <= 9; i++)\r
217     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
218   fList->Add(fHistCentStats);\r
219 \r
220   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
221   fList->Add(fHistTriggerStats);\r
222 \r
223   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
224   fList->Add(fHistTrackStats);\r
225 \r
226   fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
227   fList->Add(fHistNumberOfAcceptedTracks);\r
228 \r
229   // Vertex distributions\r
230   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
231   fList->Add(fHistVx);\r
232   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
233   fList->Add(fHistVy);\r
234   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
235   fList->Add(fHistVz);\r
236 \r
237   //Event plane\r
238   fHistEventPlane = new TH1F("fHistEventPlane",";#Psi_{2} [deg.];Counts",100,0,360.);\r
239   fList->Add(fHistEventPlane);\r
240 \r
241   // QA histograms\r
242   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
243   fList->Add(fHistClus);\r
244   fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
245   fList->Add(fHistChi2);\r
246   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
247   fList->Add(fHistDCA);\r
248   fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
249   fList->Add(fHistPt);\r
250   fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
251   fList->Add(fHistEta);\r
252   fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
253   fList->Add(fHistPhi);\r
254   fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
255   fList->Add(fHistPhiBefore);\r
256   fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
257   fList->Add(fHistPhiAfter);\r
258   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
259   fList->Add(fHistV0M);\r
260   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
261   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
262   for(Int_t i = 1; i <= 6; i++)\r
263     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
264   fList->Add(fHistRefTracks);\r
265 \r
266   // Balance function histograms\r
267   // Initialize histograms if not done yet\r
268   if(!fBalance->GetHistNp(0)){\r
269     AliWarning("Histograms not yet initialized! --> Will be done now");\r
270     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
271     fBalance->InitHistograms();\r
272   }\r
273 \r
274   if(fRunShuffling) {\r
275     if(!fShuffledBalance->GetHistNp(0)) {\r
276       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
277       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
278       fShuffledBalance->InitHistograms();\r
279     }\r
280   }\r
281 \r
282   for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
283     fListBF->Add(fBalance->GetHistNp(a));\r
284     fListBF->Add(fBalance->GetHistNn(a));\r
285     fListBF->Add(fBalance->GetHistNpn(a));\r
286     fListBF->Add(fBalance->GetHistNnn(a));\r
287     fListBF->Add(fBalance->GetHistNpp(a));\r
288     fListBF->Add(fBalance->GetHistNnp(a));\r
289 \r
290     if(fRunShuffling) {\r
291       fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
292       fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
293       fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
294       fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
295       fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
296       fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
297     }  \r
298   }\r
299 \r
300   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
301 \r
302   //====================PID========================//\r
303   if(fUsePID) {\r
304     fPIDCombined = new AliPIDCombined();\r
305     fPIDCombined->SetDefaultTPCPriors();\r
306 \r
307     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
308     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
309     \r
310     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
311     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
312     \r
313     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
314     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
315     \r
316     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
317     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
318 \r
319     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
320     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
321     \r
322     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
323     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
324     \r
325     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
326     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
327     \r
328     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
329     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
330     \r
331     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
332     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
333     \r
334     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
335     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
336   \r
337     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
338     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
339     \r
340     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
341     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
342 \r
343     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
344     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
345     \r
346     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
347     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
348   }\r
349   //====================PID========================//\r
350 \r
351   // Post output data.\r
352   PostData(1, fList);\r
353   PostData(2, fListBF);\r
354   if(fRunShuffling) PostData(3, fListBFS);\r
355   if(fUsePID) PostData(4, fHistListPIDQA);       //PID\r
356 }\r
357 \r
358 //________________________________________________________________________\r
359 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
360   // Main loop\r
361   // Called for each event\r
362   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
363 \r
364   AliESDtrack *track_TPC   = NULL;\r
365 \r
366   Int_t gNumberOfAcceptedTracks = 0;\r
367   Float_t fCentrality           = 0.;\r
368   Double_t gReactionPlane       = 0.;\r
369 \r
370   // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
371   vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
372   vector<Double_t> *chargeVector[9];          // original charge\r
373   for(Int_t i = 0; i < 9; i++){\r
374     chargeVectorShuffle[i] = new vector<Double_t>;\r
375     chargeVector[i]        = new vector<Double_t>;\r
376   }\r
377 \r
378   Double_t v_charge;\r
379   Double_t v_y;\r
380   Double_t v_eta;\r
381   Double_t v_phi;\r
382   Double_t v_p[3];\r
383   Double_t v_pt;\r
384   Double_t v_E;\r
385 \r
386   if(fUsePID) {\r
387     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
388     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
389   }\r
390  \r
391   //ESD analysis\r
392   if(gAnalysisLevel == "ESD") {\r
393     AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
394     if (!gESD) {\r
395       Printf("ERROR: gESD not available");\r
396       return;\r
397     }\r
398 \r
399     // store offline trigger bits\r
400     fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
401 \r
402     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
403     fHistEventStats->Fill(1); //all events\r
404     Bool_t isSelected = kTRUE;\r
405     if(fUseOfflineTrigger)\r
406       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
407     if(isSelected) {\r
408       fHistEventStats->Fill(2); //triggered events\r
409 \r
410       if(fUseCentrality) {\r
411         //Centrality stuff\r
412         AliCentrality *centrality = gESD->GetCentrality();\r
413         \r
414         fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
415         \r
416         // take only events inside centrality class\r
417         if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
418                                                  fCentralityPercentileMax,\r
419                                                  fCentralityEstimator.Data()))\r
420           return;\r
421         \r
422         // centrality QA (V0M)\r
423         fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
424       }\r
425         \r
426       const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
427       if(vertex) {\r
428         if(vertex->GetNContributors() > 0) {\r
429           if(vertex->GetZRes() != 0) {\r
430             fHistEventStats->Fill(3); //events with a proper vertex\r
431             if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
432               if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
433                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
434                   fHistEventStats->Fill(4); //analayzed events\r
435                   fHistVx->Fill(vertex->GetXv());\r
436                   fHistVy->Fill(vertex->GetYv());\r
437                   fHistVz->Fill(vertex->GetZv());\r
438                   \r
439                   //========Get the VZERO event plane========//\r
440                   Double_t gVZEROEventPlane = -10.0;\r
441                   Double_t qxTot = 0.0, qyTot = 0.0;\r
442                   AliEventplane *ep = gESD->GetEventplane();\r
443                   if(ep) \r
444                     gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
445                   if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
446                   fHistEventPlane->Fill(gVZEROEventPlane);\r
447                   gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
448                   //========Get the VZERO event plane========//\r
449 \r
450                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
451                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
452                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
453                     if (!track) {\r
454                       Printf("ERROR: Could not receive track %d", iTracks);\r
455                       continue;\r
456                     }   \r
457                     \r
458                     // take only TPC only tracks\r
459                     track_TPC   = new AliESDtrack();\r
460                     if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
461                     \r
462                     //ESD track cuts\r
463                     if(fESDtrackCuts) \r
464                       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
465                     \r
466                     // fill QA histograms\r
467                     Float_t b[2];\r
468                     Float_t bCov[3];\r
469                     track_TPC->GetImpactParameters(b,bCov);\r
470                     if (bCov[0]<=0 || bCov[2]<=0) {\r
471                       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
472                       bCov[0]=0; bCov[2]=0;\r
473                     }\r
474                     \r
475                     Int_t nClustersTPC = -1;\r
476                     nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
477                     //nClustersTPC = track->GetTPCclusters(0);   // global track\r
478                     Float_t chi2PerClusterTPC = -1;\r
479                     if (nClustersTPC!=0) {\r
480                       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
481                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
482                     }\r
483 \r
484                     //===========================PID===============================//               \r
485                     if(fUsePID) {\r
486                       Double_t prob[AliPID::kSPECIES]={0.};\r
487                       Double_t probTPC[AliPID::kSPECIES]={0.};\r
488                       Double_t probTOF[AliPID::kSPECIES]={0.};\r
489                       Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
490 \r
491                       Double_t nSigma = 0.;\r
492                       UInt_t detUsedTPC = 0;\r
493                       UInt_t detUsedTOF = 0;\r
494                       UInt_t detUsedTPCTOF = 0;\r
495 \r
496                       //Decide what detector configuration we want to use\r
497                       switch(fPidDetectorConfig) {\r
498                       case kTPCpid:\r
499                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
500                         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
501                         detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
502                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
503                           prob[iSpecies] = probTPC[iSpecies];\r
504                         break;\r
505                       case kTOFpid:\r
506                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
507                         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
508                         detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
509                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
510                           prob[iSpecies] = probTOF[iSpecies];\r
511                         break;\r
512                       case kTPCTOF:\r
513                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
514                         detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
515                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
516                           prob[iSpecies] = probTPCTOF[iSpecies];\r
517                         break;\r
518                       default:\r
519                         break;\r
520                       }//end switch: define detector mask\r
521                       \r
522                       //Filling the PID QA\r
523                       Double_t tofTime = -999., length = 999., tof = -999.;\r
524                       Double_t c = TMath::C()*1.E-9;// m/ns\r
525                       Double_t beta = -999.;\r
526                       Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
527                       if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
528                            (track->IsOn(AliESDtrack::kTIME))  ) { \r
529                         tofTime = track->GetTOFsignal();//in ps\r
530                         length = track->GetIntegratedLength();\r
531                         tof = tofTime*1E-3; // ns       \r
532                         \r
533                         if (tof <= 0) {\r
534                           //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
535                           continue;\r
536                         }\r
537                         if (length <= 0){\r
538                           //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
539                           continue;\r
540                         }\r
541                         \r
542                         length = length*0.01; // in meters\r
543                         tof = tof*c;\r
544                         beta = length/tof;\r
545                         \r
546                         nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
547                         fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
548                         fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
549                         fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
550                       }//TOF signal \r
551                       \r
552                       \r
553                       Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
554                       fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
555                       fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
556                       fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
557                       fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
558                       //end of QA-before pid\r
559                       \r
560                       if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
561                         //Make the decision based on the n-sigma\r
562                         if(fUsePIDnSigma) {\r
563                           if(nSigma > fPIDNSigma) continue;}\r
564                         \r
565                         //Make the decision based on the bayesian\r
566                         else if(fUsePIDPropabilities) {\r
567                           if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
568                           if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
569                         }\r
570                         \r
571                         //Fill QA after the PID\r
572                         fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
573                         fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
574                         fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
575                         \r
576                         fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
577                         fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
578                         fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
579                         fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
580                       }\r
581                       \r
582                       PostData(4, fHistListPIDQA);\r
583                     }\r
584                     //===========================PID===============================//\r
585                     v_charge = track_TPC->Charge();\r
586                     v_y      = track_TPC->Y();\r
587                     v_eta    = track_TPC->Eta();\r
588                     v_phi    = track_TPC->Phi() * TMath::RadToDeg();\r
589                     v_E      = track_TPC->E();\r
590                     v_pt     = track_TPC->Pt();\r
591                     track_TPC->PxPyPz(v_p);\r
592                     fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
593                     fHistDCA->Fill(b[1],b[0]);\r
594                     fHistChi2->Fill(chi2PerClusterTPC);\r
595                     fHistPt->Fill(v_pt);\r
596                     fHistEta->Fill(v_eta);\r
597                     fHistPhi->Fill(v_phi);\r
598 \r
599                     // fill charge vector\r
600                     chargeVector[0]->push_back(v_charge);\r
601                     chargeVector[1]->push_back(v_y);\r
602                     chargeVector[2]->push_back(v_eta);\r
603                     chargeVector[3]->push_back(v_phi);\r
604                     chargeVector[4]->push_back(v_p[0]);\r
605                     chargeVector[5]->push_back(v_p[1]);\r
606                     chargeVector[6]->push_back(v_p[2]);\r
607                     chargeVector[7]->push_back(v_pt);\r
608                     chargeVector[8]->push_back(v_E);\r
609 \r
610                     if(fRunShuffling) {\r
611                       chargeVectorShuffle[0]->push_back(v_charge);\r
612                       chargeVectorShuffle[1]->push_back(v_y);\r
613                       chargeVectorShuffle[2]->push_back(v_eta);\r
614                       chargeVectorShuffle[3]->push_back(v_phi);\r
615                       chargeVectorShuffle[4]->push_back(v_p[0]);\r
616                       chargeVectorShuffle[5]->push_back(v_p[1]);\r
617                       chargeVectorShuffle[6]->push_back(v_p[2]);\r
618                       chargeVectorShuffle[7]->push_back(v_pt);\r
619                       chargeVectorShuffle[8]->push_back(v_E);\r
620                     }\r
621                     \r
622                     delete track_TPC;\r
623                     \r
624                   } //track loop\r
625                 }//Vz cut\r
626               }//Vy cut\r
627             }//Vx cut\r
628           }//proper vertex resolution\r
629         }//proper number of contributors\r
630       }//vertex object valid\r
631     }//triggered event \r
632   }//ESD analysis\r
633   \r
634   //AOD analysis (vertex and track cuts also here!!!!)\r
635   else if(gAnalysisLevel == "AOD") {\r
636     AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
637     if(!gAOD) {\r
638       Printf("ERROR: gAOD not available");\r
639       return;\r
640     }\r
641 \r
642     AliAODHeader *aodHeader = gAOD->GetHeader();\r
643 \r
644     // store offline trigger bits\r
645     fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
646     \r
647     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
648     fHistEventStats->Fill(1); //all events\r
649     Bool_t isSelected = kTRUE;\r
650     if(fUseOfflineTrigger)\r
651       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
652     if(isSelected) {\r
653       fHistEventStats->Fill(2); //triggered events\r
654                   \r
655       //Centrality stuff (centrality in AOD header)\r
656       if(fUseCentrality) {\r
657         fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
658         \r
659         // QA for centrality estimators\r
660         fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
661         fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
662         fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
663         fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
664         fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
665         fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
666         fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
667         fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
668         fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
669         \r
670         // take only events inside centrality class\r
671         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
672           return;\r
673         \r
674         // centrality QA (V0M)\r
675         fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
676         \r
677         // centrality QA (reference tracks)\r
678         fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
679         fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
680         fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
681         fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
682         fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
683         fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
684         fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
685         fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
686         fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
687       }\r
688 \r
689       const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
690       \r
691       if(vertex) {\r
692         Double32_t fCov[6];\r
693         vertex->GetCovarianceMatrix(fCov);\r
694         \r
695         if(vertex->GetNContributors() > 0) {\r
696           if(fCov[5] != 0) {\r
697             fHistEventStats->Fill(3); //events with a proper vertex\r
698             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
699               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
700                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
701                   fHistEventStats->Fill(4); //analyzed events\r
702                   fHistVx->Fill(vertex->GetX());\r
703                   fHistVy->Fill(vertex->GetY());\r
704                   fHistVz->Fill(vertex->GetZ());\r
705         \r
706                   //========Get the VZERO event plane========//\r
707                   Double_t gVZEROEventPlane = -10.0;\r
708                   Double_t qxTot = 0.0, qyTot = 0.0;\r
709                   AliEventplane *ep = gAOD->GetEventplane();\r
710                   if(ep) \r
711                     gVZEROEventPlane = ep->CalculateVZEROEventPlane(gAOD,10,2,qxTot,qyTot);\r
712                   if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
713                   fHistEventPlane->Fill(gVZEROEventPlane);\r
714                   gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
715                   //========Get the VZERO event plane========//\r
716 \r
717                   //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
718                   for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
719                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
720                     if (!aodTrack) {\r
721                       Printf("ERROR: Could not receive track %d", iTracks);\r
722                       continue;\r
723                     }\r
724                     \r
725                     // AOD track cuts\r
726                     \r
727                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
728                     // take only TPC only tracks \r
729                     fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
730                     if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
731                     \r
732                     v_charge = aodTrack->Charge();\r
733                     v_y      = aodTrack->Y();\r
734                     v_eta    = aodTrack->Eta();\r
735                     v_phi    = aodTrack->Phi() * TMath::RadToDeg();\r
736                     v_E      = aodTrack->E();\r
737                     v_pt     = aodTrack->Pt();\r
738                     aodTrack->PxPyPz(v_p);\r
739                     \r
740                     Float_t DCAxy = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
741                     Float_t DCAz  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
742                     \r
743                     \r
744                     // Kinematics cuts from ESD track cuts\r
745                     if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
746                     if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
747                     \r
748                     // Extra DCA cuts (for systematic studies [!= -1])\r
749                     if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
750                       if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
751                         continue;  // 2D cut\r
752                       }\r
753                     }\r
754                     \r
755                     // Extra TPC cuts (for systematic studies [!= -1])\r
756                     if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
757                       continue;\r
758                     }\r
759                     if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
760                       continue;\r
761                     }\r
762                     \r
763                     // fill QA histograms\r
764                     fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
765                     fHistDCA->Fill(DCAz,DCAxy);\r
766                     fHistChi2->Fill(aodTrack->Chi2perNDF());\r
767                     fHistPt->Fill(v_pt);\r
768                     fHistEta->Fill(v_eta);\r
769                     fHistPhi->Fill(v_phi);\r
770 \r
771                     // fill charge vector\r
772                     chargeVector[0]->push_back(v_charge);\r
773                     chargeVector[1]->push_back(v_y);\r
774                     chargeVector[2]->push_back(v_eta);\r
775                     chargeVector[3]->push_back(v_phi);\r
776                     chargeVector[4]->push_back(v_p[0]);\r
777                     chargeVector[5]->push_back(v_p[1]);\r
778                     chargeVector[6]->push_back(v_p[2]);\r
779                     chargeVector[7]->push_back(v_pt);\r
780                     chargeVector[8]->push_back(v_E);\r
781 \r
782                     if(fRunShuffling) {\r
783                       chargeVectorShuffle[0]->push_back(v_charge);\r
784                       chargeVectorShuffle[1]->push_back(v_y);\r
785                       chargeVectorShuffle[2]->push_back(v_eta);\r
786                       chargeVectorShuffle[3]->push_back(v_phi);\r
787                       chargeVectorShuffle[4]->push_back(v_p[0]);\r
788                       chargeVectorShuffle[5]->push_back(v_p[1]);\r
789                       chargeVectorShuffle[6]->push_back(v_p[2]);\r
790                       chargeVectorShuffle[7]->push_back(v_pt);\r
791                       chargeVectorShuffle[8]->push_back(v_E);\r
792                     }\r
793                                     \r
794                     gNumberOfAcceptedTracks += 1;\r
795                     \r
796                   } //track loop\r
797                 }//Vz cut\r
798               }//Vy cut\r
799             }//Vx cut\r
800           }//proper vertex resolution\r
801         }//proper number of contributors\r
802       }//vertex object valid\r
803     }//triggered event \r
804   }//AOD analysis\r
805 \r
806   //MC-ESD analysis\r
807   if(gAnalysisLevel == "MCESD") {\r
808     AliMCEvent*  mcEvent = MCEvent(); \r
809     if (!mcEvent) {\r
810       Printf("ERROR: mcEvent not available");\r
811       return;\r
812     }\r
813 \r
814     AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
815     if (!gESD) {\r
816       Printf("ERROR: gESD not available");\r
817       return;\r
818     }\r
819 \r
820     // store offline trigger bits\r
821     fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
822 \r
823     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
824     fHistEventStats->Fill(1); //all events\r
825     Bool_t isSelected = kTRUE;\r
826     if(fUseOfflineTrigger)\r
827       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
828     if(isSelected) {\r
829       fHistEventStats->Fill(2); //triggered events\r
830 \r
831       if(fUseCentrality) {\r
832         //Centrality stuff\r
833         AliCentrality *centrality = gESD->GetCentrality();\r
834 \r
835         fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
836         \r
837         // take only events inside centrality class\r
838         if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
839                                                  fCentralityPercentileMax,\r
840                                                  fCentralityEstimator.Data()))\r
841           return;\r
842         \r
843         // centrality QA (V0M)\r
844         fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
845       }\r
846         \r
847       const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
848       if(vertex) {\r
849         if(vertex->GetNContributors() > 0) {\r
850           if(vertex->GetZRes() != 0) {\r
851             fHistEventStats->Fill(3); //events with a proper vertex\r
852             if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
853               if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
854                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
855                   fHistEventStats->Fill(4); //analayzed events\r
856                   fHistVx->Fill(vertex->GetXv());\r
857                   fHistVy->Fill(vertex->GetYv());\r
858                   fHistVz->Fill(vertex->GetZv());\r
859                   \r
860                   //========Get the VZERO event plane========//\r
861                   Double_t gVZEROEventPlane = -10.0;\r
862                   Double_t qxTot = 0.0, qyTot = 0.0;\r
863                   AliEventplane *ep = gESD->GetEventplane();\r
864                   if(ep) \r
865                     gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
866                   if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
867                   fHistEventPlane->Fill(gVZEROEventPlane);\r
868                   gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
869                   //========Get the VZERO event plane========//\r
870 \r
871                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
872                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
873                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
874                     if (!track) {\r
875                       Printf("ERROR: Could not receive track %d", iTracks);\r
876                       continue;\r
877                     }   \r
878                     \r
879                     Int_t label = TMath::Abs(track->GetLabel());\r
880                     if(label > mcEvent->GetNumberOfTracks()) continue;\r
881                     if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
882                     \r
883                     AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
884                     if(!mcTrack) continue;\r
885 \r
886                     // take only TPC only tracks\r
887                     track_TPC   = new AliESDtrack();\r
888                     if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
889                     \r
890                     //ESD track cuts\r
891                     if(fESDtrackCuts) \r
892                       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
893                     \r
894                     // fill QA histograms\r
895                     Float_t b[2];\r
896                     Float_t bCov[3];\r
897                     track_TPC->GetImpactParameters(b,bCov);\r
898                     if (bCov[0]<=0 || bCov[2]<=0) {\r
899                       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
900                       bCov[0]=0; bCov[2]=0;\r
901                     }\r
902                     \r
903                     Int_t nClustersTPC = -1;\r
904                     nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
905                     //nClustersTPC = track->GetTPCclusters(0);   // global track\r
906                     Float_t chi2PerClusterTPC = -1;\r
907                     if (nClustersTPC!=0) {\r
908                       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
909                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
910                     }\r
911                     \r
912                     v_charge = track_TPC->Charge();\r
913                     v_y      = track_TPC->Y();\r
914                     v_eta    = track_TPC->Eta();\r
915                     v_phi    = track_TPC->Phi() * TMath::RadToDeg();\r
916                     v_E      = track_TPC->E();\r
917                     v_pt     = track_TPC->Pt();\r
918                     track_TPC->PxPyPz(v_p);\r
919 \r
920                     fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
921                     fHistDCA->Fill(b[1],b[0]);\r
922                     fHistChi2->Fill(chi2PerClusterTPC);\r
923                     fHistPt->Fill(v_pt);\r
924                     fHistEta->Fill(v_eta);\r
925                     fHistPhi->Fill(v_phi);\r
926 \r
927                     // fill charge vector\r
928                     chargeVector[0]->push_back(v_charge);\r
929                     chargeVector[1]->push_back(v_y);\r
930                     chargeVector[2]->push_back(v_eta);\r
931                     chargeVector[3]->push_back(v_phi);\r
932                     chargeVector[4]->push_back(v_p[0]);\r
933                     chargeVector[5]->push_back(v_p[1]);\r
934                     chargeVector[6]->push_back(v_p[2]);\r
935                     chargeVector[7]->push_back(v_pt);\r
936                     chargeVector[8]->push_back(v_E);\r
937 \r
938                     if(fRunShuffling) {\r
939                       chargeVectorShuffle[0]->push_back(v_charge);\r
940                       chargeVectorShuffle[1]->push_back(v_y);\r
941                       chargeVectorShuffle[2]->push_back(v_eta);\r
942                       chargeVectorShuffle[3]->push_back(v_phi);\r
943                       chargeVectorShuffle[4]->push_back(v_p[0]);\r
944                       chargeVectorShuffle[5]->push_back(v_p[1]);\r
945                       chargeVectorShuffle[6]->push_back(v_p[2]);\r
946                       chargeVectorShuffle[7]->push_back(v_pt);\r
947                       chargeVectorShuffle[8]->push_back(v_E);\r
948                     }\r
949                     \r
950                     delete track_TPC;\r
951                     gNumberOfAcceptedTracks += 1;\r
952                     \r
953                   } //track loop\r
954                 }//Vz cut\r
955               }//Vy cut\r
956             }//Vx cut\r
957           }//proper vertex resolution\r
958         }//proper number of contributors\r
959       }//vertex object valid\r
960     }//triggered event \r
961   }//MC-ESD analysis\r
962 \r
963   //MC analysis\r
964   else if(gAnalysisLevel == "MC") {\r
965     AliMCEvent*  mcEvent = MCEvent(); \r
966     if (!mcEvent) {\r
967       Printf("ERROR: mcEvent not available");\r
968       return;\r
969     }\r
970     fHistEventStats->Fill(1); //total events\r
971     fHistEventStats->Fill(2); //offline trigger\r
972 \r
973     Double_t gImpactParameter = 0.;\r
974     if(fUseCentrality) {\r
975       //Get the MC header\r
976       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
977       if (headerH) {\r
978         //Printf("=====================================================");\r
979         //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
980         //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
981         //Printf("=====================================================");\r
982         gReactionPlane = headerH->ReactionPlaneAngle();\r
983         gImpactParameter = headerH->ImpactParameter();\r
984         fCentrality = gImpactParameter;\r
985       }\r
986       fCentrality = gImpactParameter;\r
987       fHistEventPlane->Fill(gReactionPlane*TMath::RadToDeg());\r
988 \r
989       // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)\r
990       if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
991         return;\r
992     }\r
993     \r
994     AliGenEventHeader *header = mcEvent->GenEventHeader();\r
995     if(!header) return;\r
996     \r
997     TArrayF gVertexArray;\r
998     header->PrimaryVertex(gVertexArray);\r
999     //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
1000     //gVertexArray.At(0),\r
1001     //gVertexArray.At(1),\r
1002     //gVertexArray.At(2));\r
1003     fHistEventStats->Fill(3); //events with a proper vertex\r
1004     if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
1005       if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
1006         if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
1007           fHistEventStats->Fill(4); //analayzed events\r
1008           fHistVx->Fill(gVertexArray.At(0));\r
1009           fHistVy->Fill(gVertexArray.At(1));\r
1010           fHistVz->Fill(gVertexArray.At(2));\r
1011           \r
1012           Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
1013           for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
1014             AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
1015             if (!track) {\r
1016               Printf("ERROR: Could not receive particle %d", iTracks);\r
1017               continue;\r
1018             }\r
1019             \r
1020             //exclude non stable particles\r
1021             if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;\r
1022 \r
1023             v_eta    = track->Eta();\r
1024             v_pt     = track->Pt();\r
1025             \r
1026             if( v_pt < fPtMin || v_pt > fPtMax)      \r
1027               continue;\r
1028             if( v_eta < fEtaMin || v_eta > fEtaMax)  \r
1029               continue;\r
1030             \r
1031             //analyze one set of particles\r
1032             if(fUseMCPdgCode) {\r
1033               TParticle *particle = track->Particle();\r
1034               if(!particle) continue;\r
1035               \r
1036               Int_t gPdgCode = particle->GetPdgCode();\r
1037               if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1038                 continue;\r
1039             }\r
1040             \r
1041             //Use the acceptance parameterization\r
1042             if(fAcceptanceParameterization) {\r
1043               Double_t gRandomNumber = gRandom->Rndm();\r
1044               if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1045                 continue;\r
1046             }\r
1047             \r
1048             //Exclude resonances\r
1049             if(fExcludeResonancesInMC) {\r
1050               TParticle *particle = track->Particle();\r
1051               if(!particle) continue;\r
1052               \r
1053               Bool_t kExcludeParticle = kFALSE;\r
1054               Int_t gMotherIndex = particle->GetFirstMother();\r
1055               if(gMotherIndex != -1) {\r
1056                 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1057                 if(motherTrack) {\r
1058                   TParticle *motherParticle = motherTrack->Particle();\r
1059                   if(motherParticle) {\r
1060                     Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1061                     //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1062                     if(pdgCodeOfMother == 113) {\r
1063                       kExcludeParticle = kTRUE;\r
1064                     }\r
1065                   }\r
1066                 }\r
1067               }\r
1068               \r
1069               //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1070               if(kExcludeParticle) continue;\r
1071             }\r
1072 \r
1073             v_charge = track->Charge();\r
1074             v_y      = track->Y();\r
1075             v_phi    = track->Phi();\r
1076             v_E      = track->E();\r
1077             track->PxPyPz(v_p);\r
1078             //Printf("phi (before): %lf",v_phi);\r
1079 \r
1080             fHistPt->Fill(v_pt);\r
1081             fHistEta->Fill(v_eta);\r
1082             fHistPhi->Fill(v_phi);\r
1083 \r
1084             //Flow after burner\r
1085             if(fUseFlowAfterBurner) {\r
1086               Double_t precisionPhi = 0.001;\r
1087               Int_t maxNumberOfIterations = 100;\r
1088 \r
1089               Double_t phi0 = v_phi;\r
1090               Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
1091 \r
1092               for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1093                 Double_t phiprev = v_phi;\r
1094                 Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
1095                 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
1096                 v_phi -= fl/fp;\r
1097                 if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
1098               }\r
1099               //Printf("phi (after): %lf\n",v_phi);\r
1100                       Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
1101               if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
1102               fHistPhiBefore->Fill(v_DeltaphiBefore);\r
1103 \r
1104               Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
1105               if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
1106               fHistPhiAfter->Fill(v_DeltaphiAfter);\r
1107             }\r
1108             \r
1109             v_phi *= TMath::RadToDeg();\r
1110 \r
1111             // fill charge vector\r
1112             chargeVector[0]->push_back(v_charge);\r
1113             chargeVector[1]->push_back(v_y);\r
1114             chargeVector[2]->push_back(v_eta);\r
1115             chargeVector[3]->push_back(v_phi);\r
1116             chargeVector[4]->push_back(v_p[0]);\r
1117             chargeVector[5]->push_back(v_p[1]);\r
1118             chargeVector[6]->push_back(v_p[2]);\r
1119             chargeVector[7]->push_back(v_pt);\r
1120             chargeVector[8]->push_back(v_E);\r
1121             \r
1122             if(fRunShuffling) {\r
1123               chargeVectorShuffle[0]->push_back(v_charge);\r
1124               chargeVectorShuffle[1]->push_back(v_y);\r
1125               chargeVectorShuffle[2]->push_back(v_eta);\r
1126               chargeVectorShuffle[3]->push_back(v_phi);\r
1127               chargeVectorShuffle[4]->push_back(v_p[0]);\r
1128               chargeVectorShuffle[5]->push_back(v_p[1]);\r
1129               chargeVectorShuffle[6]->push_back(v_p[2]);\r
1130               chargeVectorShuffle[7]->push_back(v_pt);\r
1131               chargeVectorShuffle[8]->push_back(v_E);\r
1132             }\r
1133             gNumberOfAcceptedTracks += 1;\r
1134                     \r
1135           } //track loop\r
1136         }//Vz cut\r
1137       }//Vy cut\r
1138     }//Vx cut\r
1139   }//MC analysis\r
1140   \r
1141   //multiplicity cut (used in pp)\r
1142   if(fUseMultiplicity) {\r
1143     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
1144       return;\r
1145   }\r
1146   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks);\r
1147   \r
1148   // calculate balance function\r
1149   if(fUseMultiplicity) \r
1150     fBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane*TMath::RadToDeg(),chargeVector);\r
1151   else                 \r
1152     fBalance->CalculateBalance(fCentrality,gReactionPlane*TMath::RadToDeg(),chargeVector);\r
1153 \r
1154   if(fRunShuffling) {\r
1155     // shuffle charges\r
1156     random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
1157     if(fUseMultiplicity) \r
1158       fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVectorShuffle);\r
1159     else                 \r
1160       fShuffledBalance->CalculateBalance(fCentrality,gReactionPlane,chargeVectorShuffle);\r
1161   }\r
1162 }      \r
1163 \r
1164 //________________________________________________________________________\r
1165 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1166   //Printf("END BF");\r
1167 \r
1168   if (!fBalance) {\r
1169     Printf("ERROR: fBalance not available");\r
1170     return;\r
1171   }  \r
1172   if(fRunShuffling) {\r
1173     if (!fShuffledBalance) {\r
1174       Printf("ERROR: fShuffledBalance not available");\r
1175       return;\r
1176     }\r
1177   }\r
1178 \r
1179 }\r
1180 \r
1181 //________________________________________________________________________\r
1182 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1183   // Draw result to the screen\r
1184   // Called once at the end of the query\r
1185 \r
1186   // not implemented ...\r
1187 \r
1188 }\r