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