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