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