Triggered Balance Function analysis (multidimensions)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskTriggeredBF.cxx
1 #include "TChain.h"\r
2 #include "TList.h"\r
3 #include "TCanvas.h"\r
4 #include "TLorentzVector.h"\r
5 #include "TGraphErrors.h"\r
6 #include "TH1F.h"\r
7 #include "TH2F.h"\r
8 #include "TArrayF.h"\r
9 #include "TF1.h"\r
10 #include "TRandom.h"\r
11 \r
12 #include "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
40 ClassImp(AliAnalysisTaskTriggeredBF)\r
41 \r
42 //________________________________________________________________________\r
43 AliAnalysisTaskTriggeredBF::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
104 AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {\r
105 \r
106   // Destructor\r
107 \r
108 }\r
109 \r
110 //________________________________________________________________________\r
111 void 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
246 void 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
256   for(Int_t i = 0; i < 9; i++){\r
257     chargeVector[i]        = new vector<Double_t>;\r
258   }\r
259   \r
260   Double_t v_charge;\r
261   Double_t v_y;\r
262   Double_t v_eta;\r
263   Double_t v_phi;\r
264   Double_t v_p[3];\r
265   Double_t v_pt;\r
266   Double_t v_E;\r
267 \r
268   // -------------------------------------------------------------                   \r
269   // AOD analysis (vertex and track cuts also here!!!!)\r
270   if(gAnalysisLevel == "AOD") {\r
271     AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(InputEvent()); \r
272     if(!aodEventMain) {\r
273       AliError("aodEventMain not available");\r
274       return;\r
275     }\r
276     \r
277     \r
278     AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\r
279     \r
280     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
281     fHistEventStats->Fill(1); //all events\r
282     \r
283     Bool_t isSelectedMain = kTRUE;\r
284     \r
285     if(fUseOfflineTrigger)\r
286       isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
287     \r
288     if(isSelectedMain) {\r
289       fHistEventStats->Fill(2); //triggered events\r
290       \r
291       //Centrality stuff (centrality in AOD header)\r
292       if(fUseCentrality) {\r
293         fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
294         \r
295         // QA for centrality estimators\r
296         fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
297         fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
298         fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
299         fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
300         fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
301         fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
302         fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
303         fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
304         fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
305         \r
306         // take only events inside centrality class\r
307         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
308           return;\r
309         \r
310         // centrality QA (V0M)\r
311         fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
312         \r
313         // centrality QA (reference tracks)\r
314         fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
315         fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
316         fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
317         fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
318         fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
319         fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
320         fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
321         fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
322         fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
323       }\r
324       \r
325       const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
326       \r
327       if(vertexMain) {\r
328         Double32_t fCovMain[6];\r
329         vertexMain->GetCovarianceMatrix(fCovMain);\r
330         \r
331         if(vertexMain->GetNContributors() > 0) {\r
332           if(fCovMain[5] != 0) {\r
333             fHistEventStats->Fill(3); //events with a proper vertex\r
334             if(TMath::Abs(vertexMain->GetX()) < fVxMax) {\r
335               if(TMath::Abs(vertexMain->GetY()) < fVyMax) {\r
336                 if(TMath::Abs(vertexMain->GetZ()) < fVzMax) {\r
337                   fHistEventStats->Fill(4); //analyzed events\r
338                   fHistVx->Fill(vertexMain->GetX());\r
339                   fHistVy->Fill(vertexMain->GetY());\r
340                   fHistVz->Fill(vertexMain->GetZ());\r
341                   \r
342                   // Loop over tracks in main event\r
343                   for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
344                     AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
345                     if (!aodTrackMain) {\r
346                       AliError(Form("Could not receive track %d", iTracksMain));\r
347                       continue;\r
348                     }\r
349                     \r
350                     // AOD track cuts\r
351                     \r
352                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
353                     // take only TPC only tracks \r
354                     fHistTrackStats->Fill(aodTrackMain->GetFilterMap());\r
355                     if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
356                     \r
357                     v_charge = aodTrackMain->Charge();\r
358                     v_y      = aodTrackMain->Y();\r
359                     v_eta    = aodTrackMain->Eta();\r
360                     v_phi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
361                     v_E      = aodTrackMain->E();\r
362                     v_pt     = aodTrackMain->Pt();\r
363                     aodTrackMain->PxPyPz(v_p);\r
364                     \r
365                     Float_t DCAxy = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
366                     Float_t DCAz  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
367                     \r
368                     \r
369                     // Kinematics cuts from ESD track cuts\r
370                     if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
371                     if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
372                     \r
373                     // Extra DCA cuts (for systematic studies [!= -1])\r
374                     if( fDCAxyCut != -1 && fDCAzCut != -1){\r
375                       if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
376                         continue;  // 2D cut\r
377                       }\r
378                     }\r
379                     \r
380                     // Extra TPC cuts (for systematic studies [!= -1])\r
381                     if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
382                       continue;\r
383                     }\r
384                     if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
385                       continue;\r
386                     }\r
387                     \r
388                     // fill QA histograms\r
389                     fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
390                     fHistDCA->Fill(DCAz,DCAxy);\r
391                     fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
392                     fHistPt->Fill(v_pt);\r
393                     fHistEta->Fill(v_eta);\r
394                     fHistPhi->Fill(v_phi);\r
395                     \r
396                     // fill charge vector\r
397                     chargeVector[0]->push_back(v_charge);\r
398                     chargeVector[1]->push_back(v_y);\r
399                     chargeVector[2]->push_back(v_eta);\r
400                     chargeVector[3]->push_back(v_phi);\r
401                     chargeVector[4]->push_back(v_p[0]);\r
402                     chargeVector[5]->push_back(v_p[1]);\r
403                     chargeVector[6]->push_back(v_p[2]);\r
404                     chargeVector[7]->push_back(v_pt);\r
405                     chargeVector[8]->push_back(v_E);\r
406                     \r
407                   } //track loop\r
408                   \r
409                     // calculate balance function\r
410                   fBalance->FillBalance(fCentrality,chargeVector);\r
411                   // clean charge vector afterwards\r
412                   for(Int_t i = 0; i < 9; i++){                \r
413                     chargeVector[i]->clear();\r
414                   }\r
415 \r
416                 }//Vz cut\r
417               }//Vy cut\r
418             }//Vx cut\r
419           }//proper vertex resolution\r
420         }//proper number of contributors\r
421       }//vertex object valid\r
422     }//triggered event \r
423   }//AOD analysis\r
424   else{\r
425     AliError("Triggered Balance Function analysis only for AODs!");\r
426   }\r
427 }     \r
428 \r
429 //________________________________________________________________________\r
430 void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r
431 \r
432   if (!fBalance) {\r
433     AliError("fBalance not available");\r
434     return;\r
435   }  \r
436   if(fRunShuffling) {\r
437     if (!fShuffledBalance) {\r
438       AliError("fShuffledBalance not available");\r
439       return;\r
440     }\r
441   }\r
442 \r
443 }\r
444 \r
445 //________________________________________________________________________\r
446 void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {\r
447   // Called once at the end of the query\r
448 \r
449   // not implemented ...\r
450 \r
451 }\r
452 \r
453 void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)\r
454 {\r
455 \r
456   // not yet done for event mixing!\r
457   return;\r
458 \r
459 }\r
460 \r