Adding the possibility to switch off the centrality usage (pp case) + addition of...
[u/mrichter/AliRoot.git] / PWG2 / EBYE / AliAnalysisTaskBF.cxx
CommitLineData
3feee083 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\r
9#include "AliAnalysisTaskSE.h"\r
10#include "AliAnalysisManager.h"\r
11\r
12#include "AliESDVertex.h"\r
13#include "AliESDEvent.h"\r
14#include "AliESDInputHandler.h"\r
15#include "AliAODEvent.h"\r
16#include "AliAODTrack.h"\r
17#include "AliAODInputHandler.h"\r
18#include "AliMCEventHandler.h"\r
19#include "AliMCEvent.h"\r
20#include "AliStack.h"\r
21#include "AliESDtrackCuts.h"\r
22\r
23#include "AliAnalysisTaskBF.h"\r
24#include "AliBalance.h"\r
25\r
26\r
27// Analysis task for the BF code\r
28// Authors: Panos Cristakoglou@cern.ch\r
29\r
30ClassImp(AliAnalysisTaskBF)\r
31\r
32//________________________________________________________________________\r
33AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name) \r
34: AliAnalysisTaskSE(name), \r
35 fBalance(0),\r
10a747d2 36 fRunShuffling(kFALSE),\r
1900bf4e 37 fShuffledBalance(0),\r
3feee083 38 fList(0),\r
39 fListBF(0),\r
6432ac6a 40 fListBFS(0),\r
3feee083 41 fHistEventStats(0),\r
a0ce3c86 42 fHistCentStats(0),\r
a839b0a3 43 fHistTriggerStats(0),\r
3feee083 44 fHistTrackStats(0),\r
45 fHistVx(0),\r
46 fHistVy(0),\r
47 fHistVz(0),\r
48 fHistClus(0),\r
49 fHistDCA(0),\r
50 fHistChi2(0),\r
51 fHistPt(0),\r
52 fHistEta(0),\r
53 fHistPhi(0),\r
54 fHistV0M(0),\r
a0ce3c86 55 fHistRefTracks(0),\r
3feee083 56 fESDtrackCuts(0),\r
32a4f715 57 fCentralityEstimator("V0M"),\r
c735aca6 58 fUseCentrality(kFALSE),\r
3feee083 59 fCentralityPercentileMin(0.), \r
60 fCentralityPercentileMax(5.),\r
61 fUseOfflineTrigger(kFALSE),\r
62 fVxMax(0.3),\r
63 fVyMax(0.3),\r
64 fVzMax(10.),\r
6432ac6a 65 nAODtrackCutBit(128),\r
3feee083 66 fPtMin(0.3),\r
67 fPtMax(1.5),\r
68 fEtaMin(-0.8),\r
1900bf4e 69 fEtaMax(-0.8),\r
a839b0a3 70 fDCAxyCut(-1),\r
71 fDCAzCut(-1),\r
72 fTPCchi2Cut(-1),\r
73 fNClustersTPCCut(-1){\r
3feee083 74 // Constructor\r
3feee083 75\r
76 // Define input and output slots here\r
77 // Input slot #0 works with a TChain\r
78 DefineInput(0, TChain::Class());\r
79 // Output slot #0 writes into a TH1 container\r
6432ac6a 80 DefineOutput(1, TList::Class());\r
81 DefineOutput(2, TList::Class());\r
3feee083 82 DefineOutput(3, TList::Class());\r
83}\r
84\r
6432ac6a 85//________________________________________________________________________\r
86AliAnalysisTaskBF::~AliAnalysisTaskBF() {\r
87\r
88 // delete fBalance; \r
89 // delete fShuffledBalance; \r
90 // delete fList;\r
91 // delete fListBF; \r
92 // delete fListBFS;\r
93\r
94 // delete fHistEventStats; \r
95 // delete fHistTrackStats; \r
96 // delete fHistVx; \r
97 // delete fHistVy; \r
98 // delete fHistVz; \r
99\r
100 // delete fHistClus;\r
101 // delete fHistDCA;\r
102 // delete fHistChi2;\r
103 // delete fHistPt;\r
104 // delete fHistEta;\r
105 // delete fHistPhi;\r
106 // delete fHistV0M;\r
107}\r
3feee083 108\r
109//________________________________________________________________________\r
110void AliAnalysisTaskBF::UserCreateOutputObjects() {\r
111 // Create histograms\r
112 // Called once\r
113 if(!fBalance) {\r
114 fBalance = new AliBalance();\r
115 fBalance->SetAnalysisLevel("ESD");\r
6432ac6a 116 //fBalance->SetNumberOfBins(-1,16);\r
117 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
3feee083 118 }\r
10a747d2 119 if(fRunShuffling) {\r
120 if(!fShuffledBalance) {\r
121 fShuffledBalance = new AliBalance();\r
122 fShuffledBalance->SetAnalysisLevel("ESD");\r
123 //fShuffledBalance->SetNumberOfBins(-1,16);\r
124 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
125 }\r
1900bf4e 126 }\r
3feee083 127\r
128 //QA list\r
129 fList = new TList();\r
130 fList->SetName("listQA");\r
131 fList->SetOwner();\r
132\r
133 //Balance Function list\r
134 fListBF = new TList();\r
135 fListBF->SetName("listBF");\r
136 fListBF->SetOwner();\r
137\r
10a747d2 138 if(fRunShuffling) {\r
139 fListBFS = new TList();\r
140 fListBFS->SetName("listBFShuffled");\r
141 fListBFS->SetOwner();\r
142 }\r
6432ac6a 143\r
3feee083 144 //Event stats.\r
145 TString gCutName[4] = {"Total","Offline trigger",\r
146 "Vertex","Analyzed"};\r
147 fHistEventStats = new TH1F("fHistEventStats",\r
148 "Event statistics;;N_{events}",\r
149 4,0.5,4.5);\r
150 for(Int_t i = 1; i <= 4; i++)\r
151 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
152 fList->Add(fHistEventStats);\r
153\r
a0ce3c86 154 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
155 fHistCentStats = new TH2F("fHistCentStats",\r
156 "Centrality statistics;;Cent percentile",\r
157 9,-0.5,8.5,220,-5,105);\r
158 for(Int_t i = 1; i <= 9; i++)\r
159 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
160 fList->Add(fHistCentStats);\r
161\r
a839b0a3 162 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
163 fList->Add(fHistTriggerStats);\r
164\r
a0ce3c86 165 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
3feee083 166 fList->Add(fHistTrackStats);\r
167\r
168 // Vertex distributions\r
169 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
170 fList->Add(fHistVx);\r
171 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
172 fList->Add(fHistVy);\r
173 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
174 fList->Add(fHistVz);\r
175\r
176 // QA histograms\r
177 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
178 fList->Add(fHistClus);\r
179 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
180 fList->Add(fHistChi2);\r
181 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
182 fList->Add(fHistDCA);\r
183 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
184 fList->Add(fHistPt);\r
185 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
186 fList->Add(fHistEta);\r
187 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
188 fList->Add(fHistPhi);\r
189 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
190 fList->Add(fHistV0M);\r
a0ce3c86 191 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
192 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
193 for(Int_t i = 1; i <= 6; i++)\r
194 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
195 fList->Add(fHistRefTracks);\r
3feee083 196\r
197\r
6432ac6a 198 // Balance function histograms\r
3feee083 199\r
6432ac6a 200 // Initialize histograms if not done yet\r
10a747d2 201 if(!fBalance->GetHistNp(0)){\r
6432ac6a 202 AliWarning("Histograms not yet initialized! --> Will be done now");\r
203 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
204 fBalance->InitHistograms();\r
10a747d2 205 }\r
206\r
207 if(fRunShuffling) {\r
208 if(!fShuffledBalance->GetHistNp(0)) {\r
209 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
210 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
211 fShuffledBalance->InitHistograms();\r
212 }\r
3feee083 213 }\r
6432ac6a 214\r
215 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
216 fListBF->Add(fBalance->GetHistNp(a));\r
217 fListBF->Add(fBalance->GetHistNn(a));\r
218 fListBF->Add(fBalance->GetHistNpn(a));\r
219 fListBF->Add(fBalance->GetHistNnn(a));\r
220 fListBF->Add(fBalance->GetHistNpp(a));\r
221 fListBF->Add(fBalance->GetHistNnp(a));\r
222\r
10a747d2 223 if(fRunShuffling) {\r
224 fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
225 fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
226 fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
227 fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
228 fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
229 fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
230 } \r
231 }\r
6432ac6a 232\r
233 if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
3feee083 234\r
235 // Post output data.\r
6432ac6a 236 PostData(1, fList);\r
237 PostData(2, fListBF);\r
10a747d2 238 if(fRunShuffling) PostData(3, fListBFS);\r
a839b0a3 239}\r
3feee083 240\r
241//________________________________________________________________________\r
242void AliAnalysisTaskBF::UserExec(Option_t *) {\r
243 // Main loop\r
244 // Called for each event\r
245 TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
246\r
1900bf4e 247 TObjArray *array = new TObjArray();\r
a839b0a3 248 AliESDtrack *track_TPC = NULL;\r
249\r
1900bf4e 250\r
251 // vector holding the charges of all tracks\r
252 vector<Int_t> chargeVectorShuffle; // this will be shuffled\r
a839b0a3 253 vector<Int_t> chargeVector; // original charge \r
3feee083 254 \r
255 //ESD analysis\r
256 if(gAnalysisLevel == "ESD") {\r
257 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
258 if (!gESD) {\r
259 Printf("ERROR: gESD not available");\r
260 return;\r
261 }\r
262\r
a839b0a3 263 // store offline trigger bits\r
264 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
265\r
6432ac6a 266 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
3feee083 267 fHistEventStats->Fill(1); //all events\r
268 Bool_t isSelected = kTRUE;\r
269 if(fUseOfflineTrigger)\r
270 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
271 if(isSelected) {\r
272 fHistEventStats->Fill(2); //triggered events\r
273\r
c735aca6 274 if(fUseCentrality) {\r
275 //Centrality stuff\r
276 AliCentrality *centrality = gESD->GetCentrality();\r
277 //Int_t nCentrality = 0;\r
278 //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
279 //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
280 \r
281 // take only events inside centrality class\r
282 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
283 fCentralityPercentileMax,\r
284 fCentralityEstimator.Data()))\r
285 return;\r
286 \r
3feee083 287 // centrality QA (V0M)\r
288 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
c735aca6 289 }\r
3feee083 290 \r
c735aca6 291 const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
292 if(vertex) {\r
293 if(vertex->GetNContributors() > 0) {\r
294 if(vertex->GetZRes() != 0) {\r
295 fHistEventStats->Fill(3); //events with a proper vertex\r
296 if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
297 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
298 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
299 fHistEventStats->Fill(4); //analayzed events\r
300 fHistVx->Fill(vertex->GetXv());\r
301 fHistVy->Fill(vertex->GetYv());\r
302 fHistVz->Fill(vertex->GetZv());\r
303 \r
304 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
305 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
306 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
307 if (!track) {\r
308 Printf("ERROR: Could not receive track %d", iTracks);\r
309 continue;\r
310 } \r
3feee083 311 \r
c735aca6 312 // take only TPC only tracks\r
313 track_TPC = new AliESDtrack();\r
314 if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
315 \r
316 //ESD track cuts\r
317 if(fESDtrackCuts) \r
318 if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
319 \r
320 // fill QA histograms\r
321 Float_t b[2];\r
322 Float_t bCov[3];\r
323 track_TPC->GetImpactParameters(b,bCov);\r
324 if (bCov[0]<=0 || bCov[2]<=0) {\r
325 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
326 bCov[0]=0; bCov[2]=0;\r
327 }\r
328 \r
329 Int_t nClustersTPC = -1;\r
330 nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
331 //nClustersTPC = track->GetTPCclusters(0); // global track\r
332 Float_t chi2PerClusterTPC = -1;\r
333 if (nClustersTPC!=0) {\r
334 chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
335 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
336 }\r
337 \r
338 fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
339 fHistDCA->Fill(b[1],b[0]);\r
340 fHistChi2->Fill(chi2PerClusterTPC);\r
341 fHistPt->Fill(track_TPC->Pt());\r
342 fHistEta->Fill(track_TPC->Eta());\r
343 fHistPhi->Fill(track_TPC->Phi()*TMath::RadToDeg());\r
344 \r
345 // fill BF array\r
346 array->Add(track_TPC);\r
347 \r
348 // fill charge vector\r
349 chargeVector.push_back(track_TPC->Charge());\r
350 if(fRunShuffling){\r
351 chargeVectorShuffle.push_back(track_TPC->Charge());\r
352 }\r
353 \r
354 delete track_TPC;\r
355 \r
356 } //track loop\r
357 }//Vz cut\r
358 }//Vy cut\r
359 }//Vx cut\r
360 }//proper vertex resolution\r
361 }//proper number of contributors\r
362 }//vertex object valid\r
3feee083 363 }//triggered event \r
364 }//ESD analysis\r
365 \r
3feee083 366 //AOD analysis (vertex and track cuts also here!!!!)\r
367 else if(gAnalysisLevel == "AOD") {\r
368 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
369 if(!gAOD) {\r
370 Printf("ERROR: gAOD not available");\r
371 return;\r
372 }\r
373\r
a839b0a3 374 AliAODHeader *aodHeader = gAOD->GetHeader();\r
375\r
376 // store offline trigger bits\r
377 fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
378\r
6432ac6a 379 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
3feee083 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
3feee083 386 \r
387 //Centrality stuff (centrality in AOD header)\r
c735aca6 388 if(fUseCentrality) {\r
389 Float_t fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
390 // cout<<fCentralityEstimator.Data()<<" = "<<fCentrality<<" , others are V0M = "\r
391 // << aodHeader->GetCentralityP()->GetCentralityPercentile("V0M")\r
392 // <<" FMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("FMD")\r
393 // <<" TRK = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TRK")\r
394 // <<" TKL = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKL")\r
395 // <<" CL0 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL0")\r
396 // <<" CL1 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL1")\r
397 // <<" V0MvsFMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")\r
398 // <<" TKLvsV0M = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")\r
399 // <<" ZEMvsZDC = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")\r
400 // <<endl;\r
401 \r
402 // QA for centrality estimators\r
403 fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
404 fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
405 fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
406 fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
407 fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
408 fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
409 fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
410 fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
411 fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
412 \r
413 // take only events inside centrality class\r
414 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
415 return;\r
416 \r
3feee083 417 // centrality QA (V0M)\r
418 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
c735aca6 419 \r
a0ce3c86 420 // centrality QA (reference tracks)\r
421 fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
422 fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
423 fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
424 fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
425 fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
426 fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
427 fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
428 fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
429 fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
c735aca6 430 }\r
a0ce3c86 431\r
c735aca6 432 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
433 \r
434 if(vertex) {\r
435 Double32_t fCov[6];\r
436 vertex->GetCovarianceMatrix(fCov);\r
437 \r
438 if(vertex->GetNContributors() > 0) {\r
439 if(fCov[5] != 0) {\r
440 fHistEventStats->Fill(3); //events with a proper vertex\r
441 if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
442 if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
443 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
444 fHistEventStats->Fill(4); //analyzed events\r
445 fHistVx->Fill(vertex->GetX());\r
446 fHistVy->Fill(vertex->GetY());\r
447 fHistVz->Fill(vertex->GetZ());\r
448 \r
449 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
450 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
451 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
452 if (!aodTrack) {\r
453 Printf("ERROR: Could not receive track %d", iTracks);\r
454 continue;\r
455 }\r
456 \r
457 // AOD track cuts\r
458 \r
459 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
460 // take only TPC only tracks \r
461 fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
462 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
463 \r
464 Float_t pt = aodTrack->Pt();\r
465 Float_t eta = aodTrack->Eta();\r
466 \r
467 Float_t DCAxy = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
468 Float_t DCAz = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
469 \r
470 \r
471 // Kinematics cuts from ESD track cuts\r
472 if( pt < fPtMin || pt > fPtMax) continue;\r
473 if( eta < fEtaMin || eta > fEtaMax) continue;\r
474 \r
475 // Extra DCA cuts (for systematic studies [!= -1])\r
476 if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
477 if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
478 continue; // 2D cut\r
10a747d2 479 }\r
c735aca6 480 }\r
481 \r
482 // Extra TPC cuts (for systematic studies [!= -1])\r
483 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
484 continue;\r
485 }\r
486 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
487 continue;\r
488 }\r
489 \r
490 \r
491 // fill QA histograms\r
492 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
493 fHistDCA->Fill(DCAz,DCAxy);\r
494 fHistChi2->Fill(aodTrack->Chi2perNDF());\r
495 fHistPt->Fill(pt);\r
496 fHistEta->Fill(eta);\r
497 fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());\r
498 \r
499 // fill BF array\r
500 array->Add(aodTrack);\r
501 \r
502 // fill charge vector\r
503 chargeVector.push_back(aodTrack->Charge());\r
504 if(fRunShuffling) {\r
505 chargeVectorShuffle.push_back(aodTrack->Charge());\r
506 }\r
507 } //track loop\r
508 }//Vz cut\r
509 }//Vy cut\r
510 }//Vx cut\r
511 }//proper vertex resolution\r
512 }//proper number of contributors\r
513 }//vertex object valid\r
3feee083 514 }//triggered event \r
515 }//AOD analysis\r
516\r
517 //MC analysis\r
518 else if(gAnalysisLevel == "MC") {\r
519 \r
520 AliMCEvent* mcEvent = MCEvent(); \r
521 if (!mcEvent) {\r
522 Printf("ERROR: mcEvent not available");\r
523 return;\r
524 }\r
525 \r
526 Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
527 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
528 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
529 if (!track) {\r
530 Printf("ERROR: Could not receive particle %d", iTracks);\r
531 continue;\r
532 }\r
1900bf4e 533\r
10a747d2 534 if( track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r
535 if( track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r
1900bf4e 536\r
10a747d2 537 array->Add(track);\r
1900bf4e 538\r
a839b0a3 539 // fill charge vector\r
540 chargeVector.push_back(track->Charge());\r
10a747d2 541 if(fRunShuffling) {\r
10a747d2 542 chargeVectorShuffle.push_back(track->Charge());\r
543 }\r
a839b0a3 544\r
3feee083 545 } //track loop\r
546 }//MC analysis\r
547 \r
10a747d2 548 // calculate balance function\r
32a4f715 549 fBalance->CalculateBalance(array,chargeVector);\r
1900bf4e 550\r
10a747d2 551 if(fRunShuffling) {\r
552 // shuffle charges\r
553 random_shuffle( chargeVectorShuffle.begin(), chargeVectorShuffle.end() );\r
554 fShuffledBalance->CalculateBalance(array,chargeVectorShuffle);\r
555 }\r
3feee083 556 \r
557 delete array;\r
558 \r
6432ac6a 559\r
3feee083 560} \r
561\r
562//________________________________________________________________________\r
563void AliAnalysisTaskBF::FinishTaskOutput(){\r
564 //Printf("END BF");\r
565\r
3feee083 566 if (!fBalance) {\r
567 Printf("ERROR: fBalance not available");\r
568 return;\r
1900bf4e 569 } \r
10a747d2 570 if(fRunShuffling) {\r
571 if (!fShuffledBalance) {\r
572 Printf("ERROR: fShuffledBalance not available");\r
573 return;\r
574 }\r
3feee083 575 }\r
3feee083 576\r
3feee083 577}\r
578\r
579//________________________________________________________________________\r
580void AliAnalysisTaskBF::Terminate(Option_t *) {\r
581 // Draw result to the screen\r
582 // Called once at the end of the query\r
583\r
584 // not implemented ...\r
6432ac6a 585\r
3feee083 586}\r