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