6b18204ff5405426d9f99561aef8984ba938fde8
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / MuonHadron / AliAnalysisTaskMuonHadronCorrelations.cxx
1 #include "AliLog.h"
2 #include "TH1D.h"
3 #include "TH2D.h"
4 #include "TTree.h"
5 #include "TAxis.h"
6 #include "AliAODEvent.h"
7 #include "AliAODTrack.h"
8 #include "TMath.h"
9 #include "TString.h"
10
11 #include "AliAnalysisTaskMuonHadronCorrelations.h"
12 #include "AliAnalysisManager.h"
13 #include "AliInputEventHandler.h"
14 #include "AliEventPoolManager.h"
15
16 ClassImp(AliAnalysisTaskMuonHadronCorrelations)
17
18 //====================================================================================================================================================
19
20 AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations() : 
21   AliAnalysisTaskSE(), 
22   fAOD(0x0),
23   fPoolMgr(0x0),
24   fTrackCB(0x0),
25   fTrackMA(0x0),
26   fFilterBitCentralBarrel(0),
27   fMaxEtaCentralBarrel(1.0),
28   fMaxChi2Muon(9999999999.), 
29   fMinRAbsMuon(0), 
30   fMaxRAbsMuon(9999999999.),
31   fTriggerMatchLevelMuon(0),
32   fNbinsCent(1), 
33   fNbinsPt(1),
34   fCentAxis(0x0), 
35   fPtAxis(0x0),
36   fHistV0Multiplicity(0x0), 
37   fHistITSMultiplicity(0x0),
38   fHistCentrality(0x0),
39   fCentMethod(0),
40   fOutputList(0x0)
41 {
42
43   // Default constructor
44
45   for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
46     for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
47       for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
48         fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = NULL;
49         fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
50       }
51     }
52     fHistNTracksCB_vs_NTracksMA[iCent]  = NULL;
53     fHistNTracksCB_vs_NTracksMAmixed[iCent]  = NULL;
54     fHistSingleMuonsEtaVsPt[iCent]   = NULL;
55     fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
56   }  
57   
58 }
59
60
61 //====================================================================================================================================================
62
63 AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations(const char *name) : 
64   AliAnalysisTaskSE(name), 
65   fAOD(0x0),
66   fPoolMgr(0x0),
67   fTrackCB(0x0),
68   fTrackMA(0x0),
69   fFilterBitCentralBarrel(0),
70   fMaxEtaCentralBarrel(1.0),
71   fMaxChi2Muon(9999999999.), 
72   fMinRAbsMuon(0), 
73   fMaxRAbsMuon(9999999999.),
74   fTriggerMatchLevelMuon(0),
75   fNbinsCent(1), 
76   fNbinsPt(1),
77   fCentAxis(0x0), 
78   fPtAxis(0x0),
79   fHistV0Multiplicity(0x0), 
80   fHistITSMultiplicity(0x0),
81   fHistCentrality(0x0),
82   fCentMethod(0),
83   fOutputList(0x0)
84 {
85
86   // Constructor
87
88   for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
89     for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
90       for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
91         fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = NULL;
92         fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
93       }
94     }
95     fHistNTracksCB_vs_NTracksMA[iCent]  = NULL;
96     fHistNTracksCB_vs_NTracksMAmixed[iCent]  = NULL;
97     fHistSingleMuonsEtaVsPt[iCent]   = NULL;
98     fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
99   }  
100   
101   // Define input and output slots here
102   DefineOutput(1, TList::Class());
103   
104 }
105
106 //====================================================================================================================================================
107
108 AliAnalysisTaskMuonHadronCorrelations::~AliAnalysisTaskMuonHadronCorrelations() {
109   
110   delete fCentAxis;
111   delete fPtAxis;
112
113   if (fOutputList  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
114     delete fOutputList;
115 }
116
117 //====================================================================================================================================================
118
119 void AliAnalysisTaskMuonHadronCorrelations::UserCreateOutputObjects() {
120   
121   fOutputList = new TList();
122   fOutputList->SetOwner(kTRUE);
123
124   for (Int_t iCent=0; iCent<fNbinsCent; iCent++) {
125     for (Int_t iPtBinCB=0; iPtBinCB<fNbinsPt; iPtBinCB++) {
126       for (Int_t iPtBinMA=0; iPtBinMA<fNbinsPt; iPtBinMA++) {
127         
128         fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = new TH1D(Form("fHistDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA), 
129                                                                Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f",
130                                                                     Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
131                                                                     Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
132                                                                     fPtAxis->GetBinLowEdge(iPtBinCB+1),
133                                                                     fPtAxis->GetBinUpEdge(iPtBinCB+1),
134                                                                     fPtAxis->GetBinLowEdge(iPtBinMA+1),
135                                                                     fPtAxis->GetBinUpEdge(iPtBinMA+1)),
136                                                                100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
137         
138         fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA]    = new TH1D(Form("fHistDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA), 
139                                                                   Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f MIXED",
140                                                                        Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
141                                                                        Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
142                                                                        fPtAxis->GetBinLowEdge(iPtBinCB+1),
143                                                                        fPtAxis->GetBinUpEdge(iPtBinCB+1),
144                                                                        fPtAxis->GetBinLowEdge(iPtBinMA+1),
145                                                                        fPtAxis->GetBinUpEdge(iPtBinMA+1)),
146                                                                   100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
147         
148         fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    -> SetXTitle("#Delta#varphi  [degrees]");
149         fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> SetXTitle("#Delta#varphi  [degrees]");
150
151         fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    -> Sumw2();
152         fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> Sumw2();
153         
154         fOutputList -> Add(fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]);
155         fOutputList -> Add(fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA]);
156
157       }
158     }
159
160     fHistNTracksCB_vs_NTracksMA[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMA_Cent%02d",iCent),
161                                                   Form("%d-%d %%",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
162                                                   100, 0, 500, 20, 0, 20);
163     fHistNTracksCB_vs_NTracksMA[iCent] -> SetXTitle("N_{tracks} Central Barrel");
164     fHistNTracksCB_vs_NTracksMA[iCent] -> SetYTitle("N_{tracks} Muon Arm");
165     fHistNTracksCB_vs_NTracksMA[iCent] -> Sumw2();
166
167     fHistNTracksCB_vs_NTracksMAmixed[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMAmixed_Cent%02d",iCent),
168                                                   Form("%d-%d %% MIXED",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
169                                                   100, 0, 500, 20, 0, 20);
170     fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetXTitle("N_{tracks} Central Barrel");
171     fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetYTitle("N_{tracks} Muon Arm");
172     fHistNTracksCB_vs_NTracksMAmixed[iCent] -> Sumw2();
173
174     fOutputList -> Add(fHistNTracksCB_vs_NTracksMA[iCent]);
175     fOutputList -> Add(fHistNTracksCB_vs_NTracksMAmixed[iCent]);
176
177     fHistSingleMuonsEtaVsPt[iCent]   = new TH2D(Form("fHistSingleMuonsEtaVsPt_Cent%02d",iCent),
178                                                "#eta vs p_{T} for single muons",
179                                                100, -4.5, -2., 100, 0., 10.);
180     fOutputList->Add(fHistSingleMuonsEtaVsPt[iCent]);
181     fHistSingleMuonsEtaVsRAbs[iCent] = new TH2D(Form("fHistSingleMuonsEtaVsRAbs_Cent%02d",iCent),
182                                                "#eta vs R_{Abs} for single muons",
183                                                100, -4.5, -2., 100, 0, 100);
184     fOutputList->Add(fHistSingleMuonsEtaVsRAbs[iCent]);
185     
186   }
187   
188   fHistV0Multiplicity = new TH1D("fHistV0Multiplicity", "V0 Multiplicity", 500, 0, 1000);
189   fHistV0Multiplicity -> SetXTitle("Multiplicity");
190   fHistV0Multiplicity -> Sumw2();
191
192   fHistITSMultiplicity = new TH1D("fHistITSMultiplicity", "ITS Multiplicity", 500, 0, 500);
193   fHistITSMultiplicity -> SetXTitle("N_{Clusters}");
194   fHistITSMultiplicity -> Sumw2();
195
196   fHistCentrality = new TH1D("fHistCentrality", Form("%s Centrality",fCentMethod.Data()), 300, -100, 200);
197   fHistCentrality -> SetXTitle("Centrality  [%]");
198   fHistCentrality -> Sumw2();
199
200   fOutputList -> Add(fHistV0Multiplicity);
201   fOutputList -> Add(fHistITSMultiplicity);
202   fOutputList -> Add(fHistCentrality);
203
204   const Int_t kNZvtxBins  = 10;
205   // bins for further buffers are shifted by 100 cm
206   Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10 };
207   Int_t nZvtxBins  = kNZvtxBins;
208   Double_t* zvtxbin = vertexBins;
209
210   fPoolMgr = new AliEventPoolManager(1000, 20000, fNbinsCent, (Double_t*)fCentAxis->GetXbins()->GetArray(), nZvtxBins, zvtxbin);
211
212   PostData(1, fOutputList); 
213
214 }
215
216 //====================================================================================================================================================
217
218 void AliAnalysisTaskMuonHadronCorrelations::UserExec(Option_t *) {
219
220   fAOD = dynamic_cast<AliAODEvent *>(InputEvent());  
221   if (!fAOD) return;  
222   
223   // Trigger selection
224   if (!(IsTriggerFired())) return;
225   
226   fHistV0Multiplicity  -> Fill(GetV0Multiplicity());
227   fHistITSMultiplicity -> Fill(GetITSMultiplicity());
228
229   Int_t centBin = GetCentBin();
230   if (centBin<0) return;
231   Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
232   fHistCentrality->Fill(percentile);
233
234   const AliVVertex* vertex = fAOD->GetPrimaryVertex();
235   Double_t zVtx = vertex->GetZ();
236   if (TMath::Abs(zVtx) > 10.) return;
237   
238   TObjArray *tracksMuonArm = GetAcceptedTracksMuonArm(fAOD,centBin);
239   if (tracksMuonArm->GetEntriesFast() == 0) {
240     delete tracksMuonArm;
241     return;
242   }
243   TObjArray *tracksCentralBarrel = GetAcceptedTracksCentralBarrel(fAOD);
244   
245   fHistNTracksCB_vs_NTracksMA[centBin] -> Fill(tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries());
246
247   AliDebug(1, Form("Single Event analysis : nTracksCB = %4d, nTracksMA = %4d", tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries()));
248
249   // Same event  
250   for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
251     fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
252     for (Int_t iTrCB=0; iTrCB<tracksCentralBarrel->GetEntriesFast(); iTrCB++) {
253       fTrackCB = (AliAODTrack*) tracksCentralBarrel -> At(iTrCB);
254       FillHistograms(centBin, kSingleEvent);
255     }
256   }
257
258   // Mixed event
259   {
260     AliEventPool* pool = fPoolMgr->GetEventPool(percentile, zVtx);
261     //pool->PrintInfo();
262     if (pool->IsReady() || pool->NTracksInPool() > 2000 || pool->GetCurrentNEvents() >= 5) 
263       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {
264         TObjArray *mixedTracks = pool->GetEvent(jMix);
265         fHistNTracksCB_vs_NTracksMAmixed[centBin]->Fill(mixedTracks->GetEntriesFast(), tracksMuonArm->GetEntriesFast());
266         for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
267           fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
268           for (Int_t iTrCB=0; iTrCB<mixedTracks->GetEntriesFast(); iTrCB++) {
269             fTrackCB = (AliAODTrack*) mixedTracks -> At(iTrCB);
270             FillHistograms(centBin, kMixedEvent);
271           }
272         }
273       }
274     pool->UpdatePool(tracksCentralBarrel);
275   }
276
277   delete tracksMuonArm;
278
279   PostData(1, fOutputList); 
280
281 }
282
283 //====================================================================================================================================================
284
285 void AliAnalysisTaskMuonHadronCorrelations::FillHistograms(Int_t centrality, Int_t option) {
286
287   Int_t ptBinTrackCB = fPtAxis -> FindBin(fTrackCB->Pt());
288   Int_t ptBinTrackMA = fPtAxis -> FindBin(fTrackMA->Pt());
289
290   if (ptBinTrackCB<1 || ptBinTrackCB>fNbinsPt || ptBinTrackMA<1 || ptBinTrackMA>fNbinsPt) return;
291
292   Double_t deltaPhi = fTrackCB->Phi() - fTrackMA->Phi();
293   if (deltaPhi >  1.5*TMath::Pi()) deltaPhi -= TMath::TwoPi();
294   if (deltaPhi < -0.5*TMath::Pi()) deltaPhi += TMath::TwoPi();
295
296   if (option==kSingleEvent)     fHistDeltaPhi[centrality][ptBinTrackCB-1][ptBinTrackMA-1]    -> Fill(TMath::RadToDeg()*deltaPhi);
297   else if (option==kMixedEvent) fHistDeltaPhiMix[centrality][ptBinTrackCB-1][ptBinTrackMA-1] -> Fill(TMath::RadToDeg()*deltaPhi);
298
299 }
300
301 //====================================================================================================================================================
302
303 Bool_t AliAnalysisTaskMuonHadronCorrelations::IsTriggerFired() {
304   
305   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); 
306
307   return isSelected;
308 }
309
310 //====================================================================================================================================================
311
312 Float_t AliAnalysisTaskMuonHadronCorrelations::GetV0Multiplicity() {
313   
314   Float_t multiplicity=0;
315   for (Int_t iChannel=0; iChannel<64; iChannel++) multiplicity += fAOD->GetVZEROData()->GetMultiplicity(iChannel);
316   return multiplicity;
317
318 }
319
320 //====================================================================================================================================================
321
322 TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksCentralBarrel(AliAODEvent *aodEvent) {
323
324   // fills the array of central barrel tracks that pass the cuts
325
326   TObjArray *tracks = new TObjArray;
327   tracks->SetOwner(kTRUE);
328
329   Int_t nTracks = aodEvent->GetNTracks();
330
331   AliAODTrack *track = 0;
332
333   for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
334     track = aodEvent->GetTrack(iTrack);
335     if (track->TestFilterBit(fFilterBitCentralBarrel) && TMath::Abs(track->Eta())<fMaxEtaCentralBarrel) {
336       tracks->Add(new AliAODTrack(*track));
337     }
338   }
339
340   return tracks;
341
342 }
343
344 //====================================================================================================================================================
345
346 TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksMuonArm(AliAODEvent *aodEvent, Int_t centBin) {
347
348   // fills the array of muon tracks that pass the cuts
349
350   TObjArray *tracks = new TObjArray;
351   tracks->SetOwner(kFALSE);
352
353   Int_t nTracks = aodEvent->GetNTracks();
354
355   AliAODTrack *track = 0;
356   
357   for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
358     track = aodEvent->GetTrack(iTrack);
359     if (track->IsMuonTrack() && track->GetMatchTrigger()>=fTriggerMatchLevelMuon) {
360       if (track->Chi2perNDF() < 5.) {
361         // histos
362         fHistSingleMuonsEtaVsPt[centBin]->Fill(track->Eta(),track->Pt());
363         fHistSingleMuonsEtaVsRAbs[centBin]->Fill(track->Eta(),track->GetRAtAbsorberEnd());
364         if (track->Eta() > -4. && track->Eta() < -2.5) {
365           Double_t rabs = track->GetRAtAbsorberEnd();
366           if (rabs > 17.6 && rabs < 89.5) {
367             tracks->Add(new AliAODTrack(*track));
368           }
369         }
370       }
371     }
372   }
373
374   return tracks;
375
376 }
377
378 //====================================================================================================================================================
379
380 void AliAnalysisTaskMuonHadronCorrelations::SetCentBinning(Int_t nBins, Double_t *limits) {
381
382   if (nBins>fNMaxBinsCentrality) {
383     AliInfo(Form("WARNING : only %d centrality bins (out of the %d proposed) will be considered",fNMaxBinsCentrality,nBins));
384     nBins = fNMaxBinsCentrality;
385   }
386   if (nBins<=0) {
387     AliInfo("WARNING : at least one centrality bin must be considered");
388     nBins = 1;
389   }
390   
391   fNbinsCent = nBins;
392   fCentAxis  = new TAxis(fNbinsCent, limits);
393
394 }
395
396 //====================================================================================================================================================
397
398 void AliAnalysisTaskMuonHadronCorrelations::SetPtBinning(Int_t nBins, Double_t *limits) {
399
400   if (nBins>fNMaxBinsPt) {
401     AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsPt,nBins));
402     nBins = fNMaxBinsPt;
403   }
404   if (nBins<=0) {
405     AliInfo("WARNING : at least one pt bin must be considered");
406     nBins = 1;
407   }
408   
409   fNbinsPt = nBins;
410   fPtAxis  = new TAxis(fNbinsPt, limits);
411
412 }
413
414 //====================================================================================================================================================
415
416 Int_t AliAnalysisTaskMuonHadronCorrelations::GetCentBin() {
417
418   Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
419
420   Int_t bin = fCentAxis->FindBin(percentile) - 1;
421   if (bin >= fNbinsCent) bin = -1;
422   return bin;
423   
424 }
425
426 //====================================================================================================================================================
427
428 Double_t AliAnalysisTaskMuonHadronCorrelations::GetITSMultiplicity() {
429
430   Double_t multiplicity = fAOD->GetHeader()->GetNumberOfITSClusters(1);
431
432   return multiplicity;
433
434 }
435
436 //====================================================================================================================================================
437
438 void AliAnalysisTaskMuonHadronCorrelations::Terminate(Option_t *) {
439
440   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
441   if (!fOutputList) {
442     printf("ERROR: Output list not available\n");
443     return;
444   }  
445
446 }
447
448 //====================================================================================================================================================
449