Add Shuffled Event Analysis to Triggered BF analysis
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskTriggeredBF.cxx
CommitLineData
f8b2d882 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 "AliLog.h"\r
13\r
14#include "AliAnalysisTaskSE.h"\r
15#include "AliAnalysisManager.h"\r
16\r
17#include "AliESDVertex.h"\r
18#include "AliESDEvent.h"\r
19#include "AliESDInputHandler.h"\r
20#include "AliAODEvent.h"\r
21#include "AliAODTrack.h"\r
22#include "AliAODInputHandler.h"\r
23#include "AliGenEventHeader.h"\r
24#include "AliGenHijingEventHeader.h"\r
25#include "AliMCEventHandler.h"\r
26#include "AliMCEvent.h"\r
27#include "AliMixInputEventHandler.h"\r
28#include "AliStack.h"\r
29\r
30#include "TH2D.h" \r
31#include "AliTHn.h" \r
32\r
33#include "AliAnalysisTaskTriggeredBF.h"\r
34#include "AliBalanceTriggered.h"\r
35\r
36\r
37// Analysis task for the TriggeredBF code\r
38// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
39\r
40ClassImp(AliAnalysisTaskTriggeredBF)\r
41\r
42//________________________________________________________________________\r
43AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name) \r
44: AliAnalysisTaskSE(name), \r
45 fBalance(0),\r
46 fRunShuffling(kFALSE),\r
47 fShuffledBalance(0),\r
48 fList(0),\r
49 fListTriggeredBF(0),\r
50 fListTriggeredBFS(0),\r
51 fHistListPIDQA(0),\r
52 fHistEventStats(0),\r
53 fHistCentStats(0),\r
54 fHistTriggerStats(0),\r
55 fHistTrackStats(0),\r
56 fHistVx(0),\r
57 fHistVy(0),\r
58 fHistVz(0),\r
59 fHistClus(0),\r
60 fHistDCA(0),\r
61 fHistChi2(0),\r
62 fHistPt(0),\r
63 fHistEta(0),\r
64 fHistPhi(0),\r
65 fHistPhiBefore(0),\r
66 fHistPhiAfter(0),\r
67 fHistV0M(0),\r
68 fHistRefTracks(0),\r
69 fCentralityEstimator("V0M"),\r
70 fUseCentrality(kFALSE),\r
71 fCentralityPercentileMin(0.), \r
72 fCentralityPercentileMax(5.),\r
73 fImpactParameterMin(0.),\r
74 fImpactParameterMax(20.),\r
75 fUseMultiplicity(kFALSE),\r
76 fNumberOfAcceptedTracksMin(0),\r
77 fNumberOfAcceptedTracksMax(10000),\r
78 fHistNumberOfAcceptedTracks(0),\r
79 fUseOfflineTrigger(kFALSE),\r
80 fVxMax(0.3),\r
81 fVyMax(0.3),\r
82 fVzMax(10.),\r
83 nAODtrackCutBit(128),\r
84 fPtMin(0.3),\r
85 fPtMax(1.5),\r
86 fEtaMin(-0.8),\r
87 fEtaMax(-0.8),\r
88 fDCAxyCut(-1),\r
89 fDCAzCut(-1),\r
90 fTPCchi2Cut(-1),\r
91 fNClustersTPCCut(-1)\r
92{\r
93 // Constructor\r
94 // Define input and output slots here\r
95 // Input slot #0 works with a TChain\r
96 DefineInput(0, TChain::Class());\r
97 // Output slot #0 writes into a TH1 container\r
98 DefineOutput(1, TList::Class());\r
99 DefineOutput(2, TList::Class());\r
100 DefineOutput(3, TList::Class());\r
101}\r
102\r
103//________________________________________________________________________\r
104AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {\r
105\r
106 // Destructor\r
107\r
108}\r
109\r
110//________________________________________________________________________\r
111void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {\r
112 // Create histograms\r
113 // Called once\r
114 if(!fBalance) {\r
115 fBalance = new AliBalanceTriggered();\r
116 fBalance->SetAnalysisLevel("AOD");\r
117 }\r
118 if(fRunShuffling) {\r
119 if(!fShuffledBalance) {\r
120 fShuffledBalance = new AliBalanceTriggered();\r
121 fShuffledBalance->SetAnalysisLevel("AOD");\r
122 }\r
123 }\r
124\r
125 //QA list\r
126 fList = new TList();\r
127 fList->SetName("listQA");\r
128 fList->SetOwner();\r
129\r
130 //Balance Function list\r
131 fListTriggeredBF = new TList();\r
132 fListTriggeredBF->SetName("listTriggeredBF");\r
133 fListTriggeredBF->SetOwner();\r
134\r
135 if(fRunShuffling) {\r
136 fListTriggeredBFS = new TList();\r
137 fListTriggeredBFS->SetName("listTriggeredBFShuffled");\r
138 fListTriggeredBFS->SetOwner();\r
139 }\r
140 \r
141 \r
142 //Event stats.\r
143 TString gCutName[4] = {"Total","Offline trigger",\r
144 "Vertex","Analyzed"};\r
145 fHistEventStats = new TH1F("fHistEventStats",\r
146 "Event statistics;;N_{events}",\r
147 4,0.5,4.5);\r
148 for(Int_t i = 1; i <= 4; i++)\r
149 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
150 fList->Add(fHistEventStats);\r
151 \r
152 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
153 fHistCentStats = new TH2F("fHistCentStats",\r
154 "Centrality statistics;;Cent percentile",\r
155 9,-0.5,8.5,220,-5,105);\r
156 for(Int_t i = 1; i <= 9; i++)\r
157 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
158 fList->Add(fHistCentStats);\r
159 \r
160 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
161 fList->Add(fHistTriggerStats);\r
162 \r
163 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
164 fList->Add(fHistTrackStats);\r
165\r
166 fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
167 fList->Add(fHistNumberOfAcceptedTracks);\r
168\r
169 // Vertex distributions\r
170 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
171 fList->Add(fHistVx);\r
172 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
173 fList->Add(fHistVy);\r
174 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
175 fList->Add(fHistVz);\r
176\r
177 // QA histograms\r
178 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
179 fList->Add(fHistClus);\r
180 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
181 fList->Add(fHistChi2);\r
182 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
183 fList->Add(fHistDCA);\r
184 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
185 fList->Add(fHistPt);\r
186 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
187 fList->Add(fHistEta);\r
188 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
189 fList->Add(fHistPhi);\r
190 fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
191 fList->Add(fHistPhiBefore);\r
192 fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
193 fList->Add(fHistPhiAfter);\r
194 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
195 fList->Add(fHistV0M);\r
196 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
197 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
198 for(Int_t i = 1; i <= 6; i++)\r
199 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
200 fList->Add(fHistRefTracks);\r
201\r
202\r
203\r
204 // Balance function histograms\r
205 // Initialize histograms if not done yet\r
206 if(!fBalance->GetHistNp()){\r
207 AliWarning("Histograms not yet initialized! --> Will be done now");\r
208 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
209 fBalance->InitHistograms();\r
210 }\r
211\r
212 if(fRunShuffling) {\r
213 if(!fShuffledBalance->GetHistNp()) {\r
214 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
215 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
216 fShuffledBalance->InitHistograms();\r
217 }\r
218 }\r
219\r
220 fListTriggeredBF->Add(fBalance->GetHistNp());\r
221 fListTriggeredBF->Add(fBalance->GetHistNn());\r
222 fListTriggeredBF->Add(fBalance->GetHistNpn());\r
223 fListTriggeredBF->Add(fBalance->GetHistNnn());\r
224 fListTriggeredBF->Add(fBalance->GetHistNpp());\r
225 fListTriggeredBF->Add(fBalance->GetHistNnp());\r
226 \r
227 if(fRunShuffling) {\r
228 fListTriggeredBFS->Add(fShuffledBalance->GetHistNp());\r
229 fListTriggeredBFS->Add(fShuffledBalance->GetHistNn());\r
230 fListTriggeredBFS->Add(fShuffledBalance->GetHistNpn());\r
231 fListTriggeredBFS->Add(fShuffledBalance->GetHistNnn());\r
232 fListTriggeredBFS->Add(fShuffledBalance->GetHistNpp());\r
233 fListTriggeredBFS->Add(fShuffledBalance->GetHistNnp());\r
234 } \r
235\r
236\r
237 \r
238\r
239 // Post output data.\r
240 PostData(1, fList);\r
241 PostData(2, fListTriggeredBF);\r
242 if(fRunShuffling) PostData(3, fListTriggeredBFS);\r
243}\r
244\r
245//________________________________________________________________________\r
246void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {\r
247 // Main loop\r
248 // Called for each event\r
249\r
250 TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
251\r
252 Float_t fCentrality = 0.;\r
253 \r
254 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
255 vector<Double_t> *chargeVector[9]; // original charge\r
87961e2b 256 vector<Double_t> *chargeVectorShuffled[9]; // shuffled charge\r
257\r
f8b2d882 258 for(Int_t i = 0; i < 9; i++){\r
87961e2b 259 chargeVector[i] = new vector<Double_t>;\r
260 chargeVectorShuffled[i] = new vector<Double_t>;\r
f8b2d882 261 }\r
262 \r
263 Double_t v_charge;\r
264 Double_t v_y;\r
265 Double_t v_eta;\r
266 Double_t v_phi;\r
267 Double_t v_p[3];\r
268 Double_t v_pt;\r
269 Double_t v_E;\r
270\r
271 // ------------------------------------------------------------- \r
272 // AOD analysis (vertex and track cuts also here!!!!)\r
273 if(gAnalysisLevel == "AOD") {\r
274 AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(InputEvent()); \r
275 if(!aodEventMain) {\r
276 AliError("aodEventMain not available");\r
277 return;\r
278 }\r
279 \r
280 \r
281 AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\r
282 \r
283 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
284 fHistEventStats->Fill(1); //all events\r
285 \r
286 Bool_t isSelectedMain = kTRUE;\r
287 \r
288 if(fUseOfflineTrigger)\r
289 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
290 \r
291 if(isSelectedMain) {\r
292 fHistEventStats->Fill(2); //triggered events\r
293 \r
294 //Centrality stuff (centrality in AOD header)\r
295 if(fUseCentrality) {\r
296 fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
297 \r
298 // QA for centrality estimators\r
299 fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
300 fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
301 fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
302 fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
303 fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
304 fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
305 fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
306 fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
307 fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
308 \r
309 // take only events inside centrality class\r
310 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
311 return;\r
312 \r
313 // centrality QA (V0M)\r
314 fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
315 \r
316 // centrality QA (reference tracks)\r
317 fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
318 fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
319 fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
320 fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
321 fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
322 fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
323 fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
324 fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
325 fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
326 }\r
327 \r
328 const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
329 \r
330 if(vertexMain) {\r
331 Double32_t fCovMain[6];\r
332 vertexMain->GetCovarianceMatrix(fCovMain);\r
333 \r
334 if(vertexMain->GetNContributors() > 0) {\r
335 if(fCovMain[5] != 0) {\r
336 fHistEventStats->Fill(3); //events with a proper vertex\r
337 if(TMath::Abs(vertexMain->GetX()) < fVxMax) {\r
338 if(TMath::Abs(vertexMain->GetY()) < fVyMax) {\r
339 if(TMath::Abs(vertexMain->GetZ()) < fVzMax) {\r
340 fHistEventStats->Fill(4); //analyzed events\r
341 fHistVx->Fill(vertexMain->GetX());\r
342 fHistVy->Fill(vertexMain->GetY());\r
343 fHistVz->Fill(vertexMain->GetZ());\r
344 \r
345 // Loop over tracks in main event\r
346 for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
347 AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
348 if (!aodTrackMain) {\r
349 AliError(Form("Could not receive track %d", iTracksMain));\r
350 continue;\r
351 }\r
352 \r
353 // AOD track cuts\r
354 \r
355 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
356 // take only TPC only tracks \r
357 fHistTrackStats->Fill(aodTrackMain->GetFilterMap());\r
358 if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
359 \r
360 v_charge = aodTrackMain->Charge();\r
361 v_y = aodTrackMain->Y();\r
362 v_eta = aodTrackMain->Eta();\r
363 v_phi = aodTrackMain->Phi() * TMath::RadToDeg();\r
364 v_E = aodTrackMain->E();\r
365 v_pt = aodTrackMain->Pt();\r
366 aodTrackMain->PxPyPz(v_p);\r
367 \r
368 Float_t DCAxy = aodTrackMain->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
369 Float_t DCAz = aodTrackMain->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
370 \r
371 \r
372 // Kinematics cuts from ESD track cuts\r
373 if( v_pt < fPtMin || v_pt > fPtMax) continue;\r
374 if( v_eta < fEtaMin || v_eta > fEtaMax) continue;\r
375 \r
376 // Extra DCA cuts (for systematic studies [!= -1])\r
377 if( fDCAxyCut != -1 && fDCAzCut != -1){\r
378 if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
379 continue; // 2D cut\r
380 }\r
381 }\r
382 \r
383 // Extra TPC cuts (for systematic studies [!= -1])\r
384 if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
385 continue;\r
386 }\r
387 if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
388 continue;\r
389 }\r
390 \r
391 // fill QA histograms\r
392 fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
393 fHistDCA->Fill(DCAz,DCAxy);\r
394 fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
395 fHistPt->Fill(v_pt);\r
396 fHistEta->Fill(v_eta);\r
397 fHistPhi->Fill(v_phi);\r
398 \r
399 // fill charge vector\r
400 chargeVector[0]->push_back(v_charge);\r
401 chargeVector[1]->push_back(v_y);\r
402 chargeVector[2]->push_back(v_eta);\r
403 chargeVector[3]->push_back(v_phi);\r
404 chargeVector[4]->push_back(v_p[0]);\r
405 chargeVector[5]->push_back(v_p[1]);\r
406 chargeVector[6]->push_back(v_p[2]);\r
407 chargeVector[7]->push_back(v_pt);\r
408 chargeVector[8]->push_back(v_E);\r
87961e2b 409\r
410 if(fRunShuffling) {\r
411 chargeVectorShuffled[0]->push_back(v_charge);\r
412 chargeVectorShuffled[1]->push_back(v_y);\r
413 chargeVectorShuffled[2]->push_back(v_eta);\r
414 chargeVectorShuffled[3]->push_back(v_phi);\r
415 chargeVectorShuffled[4]->push_back(v_p[0]);\r
416 chargeVectorShuffled[5]->push_back(v_p[1]);\r
417 chargeVectorShuffled[6]->push_back(v_p[2]);\r
418 chargeVectorShuffled[7]->push_back(v_pt);\r
419 chargeVectorShuffled[8]->push_back(v_E);\r
420 }\r
f8b2d882 421 \r
422 } //track loop\r
423 \r
87961e2b 424 // calculate balance function\r
f8b2d882 425 fBalance->FillBalance(fCentrality,chargeVector);\r
87961e2b 426 \r
427 // calculate shuffled balance function\r
428 if(fRunShuffling) {\r
429 random_shuffle(chargeVectorShuffled[0]->begin(), chargeVectorShuffled[0]->end());\r
430 fShuffledBalance->FillBalance(fCentrality,chargeVectorShuffled);\r
431 }\r
432\r
f8b2d882 433 // clean charge vector afterwards\r
434 for(Int_t i = 0; i < 9; i++){ \r
435 chargeVector[i]->clear();\r
87961e2b 436 chargeVectorShuffled[i]->clear();\r
f8b2d882 437 }\r
438\r
439 }//Vz cut\r
440 }//Vy cut\r
441 }//Vx cut\r
442 }//proper vertex resolution\r
443 }//proper number of contributors\r
444 }//vertex object valid\r
445 }//triggered event \r
446 }//AOD analysis\r
447 else{\r
448 AliError("Triggered Balance Function analysis only for AODs!");\r
449 }\r
450} \r
451\r
452//________________________________________________________________________\r
453void AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r
454\r
455 if (!fBalance) {\r
456 AliError("fBalance not available");\r
457 return;\r
458 } \r
459 if(fRunShuffling) {\r
460 if (!fShuffledBalance) {\r
461 AliError("fShuffledBalance not available");\r
462 return;\r
463 }\r
464 }\r
465\r
466}\r
467\r
468//________________________________________________________________________\r
469void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {\r
470 // Called once at the end of the query\r
471\r
472 // not implemented ...\r
473\r
474}\r
475\r
476void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)\r
477{\r
478\r
479 // not yet done for event mixing!\r
480 return;\r
481\r
482}\r
483\r