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