Update event mixing
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / MuonHadron / AliAnalysisTaskDiMuonCorrelations.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 "AliAnalysisTaskDiMuonCorrelations.h"
12 #include "AliAnalysisManager.h"
13 #include "AliInputEventHandler.h"
14 #include "AliEventPoolManager.h"
15
16 ClassImp(AliAnalysisTaskDiMuonCorrelations)
17
18 //====================================================================================================================================================
19
20 AliAnalysisTaskDiMuonCorrelations::AliAnalysisTaskDiMuonCorrelations() : 
21   AliAnalysisTaskSE(), 
22   fAOD(0x0),
23   fPoolMgr(0x0),
24   fMaxChi2Muon(5.), 
25   fMinEtaMuon(-4.0), 
26   fMaxEtaMuon(-2.5),
27   fTriggerMatchLevelMuon(0),
28   fNbinsCent(1), 
29   fNbinsPt(1),
30   fCentAxis(0x0), 
31   fPtAxis(0x0),
32   fEtaAxis(0x0),
33   fHistV0Multiplicity(0x0), 
34   fHistITSMultiplicity(0x0),
35   fHistCentrality(0x0),
36   fHistEvStat(0x0),
37   fCentMethod(0),
38   fOutputList(0x0)
39 {
40
41   // Default constructor
42
43   fMuonTrack[0] = NULL;
44   fMuonTrack[1] = NULL;
45
46   for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
47     for (Int_t iPtBinMuon1=0; iPtBinMuon1<fNMaxBinsPt; iPtBinMuon1++) {
48       for (Int_t iPtBinMuon2=0; iPtBinMuon2<fNMaxBinsPt; iPtBinMuon2++) {
49         fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]       = NULL;
50         fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2]    = NULL;
51         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    = NULL;
52         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] = NULL;
53       }
54     }
55     fHistNMuons_vs_NMuons[iCent]         = NULL;
56     fHistNMuons_vs_NMuons_Mixed[iCent]   = NULL;
57     fHistTracksEtavsEta[iCent]           = NULL;
58     fHistTracksEtavsEta_Mixed[iCent]     = NULL;
59     fHistSingleMuonsPt[iCent]            = NULL;
60     fHistSingleMuonsPt_Mixed[iCent]      = NULL;
61     fHistSingleMuonsEtaPt[iCent]         = NULL;
62     fHistSingleMuonsEtaPt_Mixed[iCent]   = NULL;
63   }  
64   
65 }
66
67
68 //====================================================================================================================================================
69
70 AliAnalysisTaskDiMuonCorrelations::AliAnalysisTaskDiMuonCorrelations(const char *name) : 
71   AliAnalysisTaskSE(name), 
72   fAOD(0x0),
73   fPoolMgr(0x0),
74   fMaxChi2Muon(5.), 
75   fMinEtaMuon(-4.0), 
76   fMaxEtaMuon(-2.5),
77   fTriggerMatchLevelMuon(0),
78   fNbinsCent(1), 
79   fNbinsPt(1),
80   fCentAxis(0x0), 
81   fPtAxis(0x0),
82   fEtaAxis(0x0),
83   fHistV0Multiplicity(0x0), 
84   fHistITSMultiplicity(0x0),
85   fHistCentrality(0x0),
86   fHistEvStat(0x0),
87   fCentMethod(0),
88   fOutputList(0x0)
89 {
90
91   // Constructor
92
93   fMuonTrack[0] = NULL;
94   fMuonTrack[1] = NULL;
95
96   for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
97     for (Int_t iPtBinMuon1=0; iPtBinMuon1<fNMaxBinsPt; iPtBinMuon1++) {
98       for (Int_t iPtBinMuon2=0; iPtBinMuon2<fNMaxBinsPt; iPtBinMuon2++) {
99         fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    = NULL;
100         fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] = NULL;
101         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    = NULL;
102         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] = NULL;
103       }
104     }
105     fHistNMuons_vs_NMuons[iCent]         = NULL;
106     fHistNMuons_vs_NMuons_Mixed[iCent]   = NULL;
107     fHistTracksEtavsEta[iCent]           = NULL;
108     fHistTracksEtavsEta_Mixed[iCent]     = NULL;
109     fHistSingleMuonsPt[iCent]            = NULL;
110     fHistSingleMuonsPt_Mixed[iCent]      = NULL;
111     fHistSingleMuonsEtaPt[iCent]         = NULL;
112     fHistSingleMuonsEtaPt_Mixed[iCent]   = NULL;
113   }  
114   
115   // Define input and output slots here
116   DefineOutput(1, TList::Class());
117   
118 }
119
120 //====================================================================================================================================================
121
122 AliAnalysisTaskDiMuonCorrelations::~AliAnalysisTaskDiMuonCorrelations() {
123   
124   delete fCentAxis;
125   delete fPtAxis;
126   delete fEtaAxis;
127
128   if (fOutputList  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
129
130 }
131
132 //====================================================================================================================================================
133
134 void AliAnalysisTaskDiMuonCorrelations::UserCreateOutputObjects() {
135   
136   fOutputList = new TList();
137   fOutputList->SetOwner(kTRUE);
138
139   for (Int_t iCent=0; iCent<fNbinsCent; iCent++) {
140     for (Int_t iPtBinMuon1=0; iPtBinMuon1<fNbinsPt; iPtBinMuon1++) {
141       for (Int_t iPtBinMuon2=0; iPtBinMuon2<fNbinsPt; iPtBinMuon2++) {
142         
143         fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    = new TH1D(Form("fHistDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinMuon1,iPtBinMuon2), 
144                                                                      Form("%d-%d %%, %3.1f<p_{T}^{#mu1}<%3.1f, %3.1f<p_{T}^{#mu2}<%3.1f",
145                                                                           Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
146                                                                           Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
147                                                                           fPtAxis->GetBinLowEdge(iPtBinMuon1+1),
148                                                                           fPtAxis->GetBinUpEdge(iPtBinMuon1+1),
149                                                                           fPtAxis->GetBinLowEdge(iPtBinMuon2+1),
150                                                                           fPtAxis->GetBinUpEdge(iPtBinMuon2+1)),
151                                                                      100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
152         
153         fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2]    = new TH1D(Form("fHistDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinMuon1,iPtBinMuon2), 
154                                                                         Form("%d-%d %%, %3.1f<p_{T}^{#mu1}<%3.1f, %3.1f<p_{T}^{#mu2}<%3.1f MIXED",
155                                                                              Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
156                                                                              Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
157                                                                              fPtAxis->GetBinLowEdge(iPtBinMuon1+1),
158                                                                              fPtAxis->GetBinUpEdge(iPtBinMuon1+1),
159                                                                              fPtAxis->GetBinLowEdge(iPtBinMuon2+1),
160                                                                              fPtAxis->GetBinUpEdge(iPtBinMuon2+1)),
161                                                                         100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
162         
163         fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    -> SetXTitle("#Delta#varphi  [degrees]");
164         fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] -> SetXTitle("#Delta#varphi  [degrees]");
165         
166         fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    -> Sumw2();
167         fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] -> Sumw2();
168         
169         fOutputList -> Add(fHistDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]);
170         fOutputList -> Add(fHistDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2]);
171         
172         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    = new TH2D(Form("fHistEtaDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinMuon1,iPtBinMuon2), 
173                                                                         Form("%d-%d %%, %3.1f<p_{T}^{#mu1}<%3.1f, %3.1f<p_{T}^{#mu2}<%3.1f",
174                                                                              Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
175                                                                              Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
176                                                                              fPtAxis->GetBinLowEdge(iPtBinMuon1+1),
177                                                                              fPtAxis->GetBinUpEdge(iPtBinMuon1+1),
178                                                                              fPtAxis->GetBinLowEdge(iPtBinMuon2+1),
179                                                                              fPtAxis->GetBinUpEdge(iPtBinMuon2+1)),
180                                                                         100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi(),
181                                                                         fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
182         
183         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2]    = new TH2D(Form("fHistEtaDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinMuon1,iPtBinMuon2), 
184                                                                            Form("%d-%d %%, %3.1f<p_{T}^{#mu1}<%3.1f, %3.1f<p_{T}^{#mu2}<%3.1f MIXED",
185                                                                                 Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
186                                                                                 Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
187                                                                                 fPtAxis->GetBinLowEdge(iPtBinMuon1+1),
188                                                                                 fPtAxis->GetBinUpEdge(iPtBinMuon1+1),
189                                                                                 fPtAxis->GetBinLowEdge(iPtBinMuon2+1),
190                                                                                 fPtAxis->GetBinUpEdge(iPtBinMuon2+1)),
191                                                                            100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi(),
192                                                                            fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
193         
194         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    -> SetXTitle("#Delta#varphi  [degrees]");
195         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] -> SetXTitle("#Delta#varphi  [degrees]");
196         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    -> SetYTitle("#eta");
197         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] -> SetYTitle("#eta");
198         
199         fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]    -> Sumw2();
200         fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2] -> Sumw2();
201         
202         fOutputList -> Add(fHistEtaDeltaPhi[iCent][iPtBinMuon1][iPtBinMuon2]);
203         fOutputList -> Add(fHistEtaDeltaPhiMix[iCent][iPtBinMuon1][iPtBinMuon2]);
204
205       }
206     }
207
208     fHistNMuons_vs_NMuons[iCent] = new TH2D(Form("fHistNMuons_vs_NMuons_Cent%02d",iCent),
209                                                   Form("%d-%d %%",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
210                                                   100, 0, 500, 20, 0, 20);
211     fHistNMuons_vs_NMuons[iCent] -> SetXTitle("N_{tracks} Muon Arm");
212     fHistNMuons_vs_NMuons[iCent] -> SetYTitle("N_{tracks} Muon Arm");
213     fHistNMuons_vs_NMuons[iCent] -> Sumw2();
214
215     fHistNMuons_vs_NMuons_Mixed[iCent] = new TH2D(Form("fHistNMuons_vs_NMuons_Mixed_Cent%02d",iCent),
216                                                   Form("%d-%d %% MIXED",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
217                                                   100, 0, 500, 20, 0, 20);
218     fHistNMuons_vs_NMuons_Mixed[iCent] -> SetXTitle("N_{tracks} Muon Arm");
219     fHistNMuons_vs_NMuons_Mixed[iCent] -> SetYTitle("N_{tracks} Muon Arm");
220     fHistNMuons_vs_NMuons_Mixed[iCent] -> Sumw2();
221
222     fHistTracksEtavsEta[iCent]       = new TH2D(Form("fHistTracksEtavsEta_%02d",iCent),       "#eta 1st muon vs #eta 2nd muon", 100,-4.5,-2.,100,-1.5,1.5);
223     fHistTracksEtavsEta_Mixed[iCent] = new TH2D(Form("fHistTracksEtavsEta_Mixed_%02d",iCent), "#eta 1st muon vs #eta 2nd muon", 100,-4.5,-2.,100,-1.5,1.5);
224
225     fOutputList -> Add(fHistNMuons_vs_NMuons[iCent]);
226     fOutputList -> Add(fHistNMuons_vs_NMuons_Mixed[iCent]);
227     fOutputList -> Add(fHistTracksEtavsEta[iCent]);
228     fOutputList -> Add(fHistTracksEtavsEta_Mixed[iCent]);
229
230     fHistSingleMuonsPt[iCent]       = new TH1D(Form("fHistSingleMuonPt_Cent%02d",iCent),      "p_{T} for single muons", fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray());
231     fHistSingleMuonsPt_Mixed[iCent] = new TH1D(Form("fHistSingleMuonPtmixed_Cent%02d",iCent), "p_{T} for single muons", fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray());
232     fOutputList -> Add(fHistSingleMuonsPt[iCent]);
233     fOutputList -> Add(fHistSingleMuonsPt_Mixed[iCent]);
234
235     fHistSingleMuonsEtaPt[iCent] = new TH2D(Form("fHistSingleMuonEtaPt_Cent%02d",iCent), "#eta vs p_{T} for single muons",
236                                             fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray(),
237                                             fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
238     fHistSingleMuonsEtaPt_Mixed[iCent] = new TH2D(Form("fHistSingleMuonEtaPtmixed_Cent%02d",iCent), "#eta vs p_{T} for single muons",
239                                                   fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray(),
240                                                   fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
241     fOutputList -> Add(fHistSingleMuonsEtaPt[iCent]);
242     fOutputList -> Add(fHistSingleMuonsEtaPt_Mixed[iCent]);
243
244   }
245   
246   fHistV0Multiplicity = new TH1D("fHistV0Multiplicity", "V0 Multiplicity", 500, 0, 1000);
247   fHistV0Multiplicity -> SetXTitle("Multiplicity");
248   fHistV0Multiplicity -> Sumw2();
249
250   fHistITSMultiplicity = new TH1D("fHistITSMultiplicity", "ITS Multiplicity", 500, 0, 500);
251   fHistITSMultiplicity -> SetXTitle("N_{Clusters}");
252   fHistITSMultiplicity -> Sumw2();
253
254   fHistCentrality = new TH1D("fHistCentrality", Form("%s Centrality",fCentMethod.Data()), 300, -100, 200);
255   fHistCentrality -> SetXTitle("Centrality  [%]");
256   fHistCentrality -> Sumw2();
257
258   fOutputList -> Add(fHistV0Multiplicity);
259   fOutputList -> Add(fHistITSMultiplicity);
260   fOutputList -> Add(fHistCentrality);
261
262   fHistEvStat = new TH1D("fHistEvStat","Event cuts statistics",20,-0.5,19.5);
263   fHistEvStat->SetXTitle("Cut index");
264   fOutputList->Add(fHistEvStat);
265
266   const Int_t kNZvtxBins  = 10;
267   // bins for further buffers are shifted by 100 cm
268   Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10 };
269   Int_t nZvtxBins  = kNZvtxBins;
270   Double_t* zvtxbin = vertexBins;
271
272   fPoolMgr = new AliEventPoolManager(1000, 20000, fNbinsCent, (Double_t*)fCentAxis->GetXbins()->GetArray(), nZvtxBins, zvtxbin);
273
274   PostData(1, fOutputList); 
275
276 }
277
278 //====================================================================================================================================================
279
280 void AliAnalysisTaskDiMuonCorrelations::UserExec(Option_t *) {
281
282   fAOD = dynamic_cast<AliAODEvent *>(InputEvent());  
283   if (!fAOD) return;  
284
285   Int_t cutIndex = 0;
286   fHistEvStat->Fill(cutIndex++);
287   // Trigger selection
288   if (!(IsTriggerFired())) return;
289   fHistEvStat->Fill(cutIndex++);
290   
291   fHistV0Multiplicity  -> Fill(GetV0Multiplicity());
292   fHistITSMultiplicity -> Fill(GetITSMultiplicity());
293
294   Int_t centBin = GetCentBin();
295   if (centBin<0) return;
296   fHistEvStat->Fill(cutIndex++);
297   Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
298   fHistCentrality->Fill(percentile);
299
300   // Vertex selection
301   const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
302   if (!trkVtx || trkVtx->GetNContributors()<=0) return;
303   TString vtxTtl = trkVtx->GetTitle();
304   if (!vtxTtl.Contains("VertexerTracks")) return;
305   fHistEvStat->Fill(cutIndex++);
306   Double_t zvtx = trkVtx->GetZ();
307   const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
308   if (spdVtx->GetNContributors()<=0) return;
309   TString vtxTyp = spdVtx->GetTitle();
310   Double_t cov[6]={0};
311   spdVtx->GetCovarianceMatrix(cov);
312   Double_t zRes = TMath::Sqrt(cov[5]);
313   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
314   if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
315   fHistEvStat->Fill(cutIndex++);
316
317   if (TMath::Abs(zvtx) > 10.) return;
318   fHistEvStat->Fill(cutIndex++);
319
320   TObjArray *tracksMuonArm = GetAcceptedTracksMuonArm(fAOD);
321   if (tracksMuonArm->GetEntriesFast() < 2) {
322     delete tracksMuonArm;
323     return;
324   }
325   fHistEvStat->Fill(cutIndex++);
326   
327   fHistNMuons_vs_NMuons[centBin] -> Fill(tracksMuonArm->GetEntries(), tracksMuonArm->GetEntries());
328
329   AliDebug(1, Form("Single Event analysis : nTracksMuonArm = %4d", tracksMuonArm->GetEntries()));
330
331   // Same event  
332   for (Int_t iTrMuon1=0; iTrMuon1<tracksMuonArm->GetEntriesFast(); iTrMuon1++) {
333     fMuonTrack[0] = (AliAODTrack*) tracksMuonArm->At(iTrMuon1);
334     fHistSingleMuonsPt[centBin]->Fill(fMuonTrack[0]->Pt());
335     fHistSingleMuonsEtaPt[centBin]->Fill(fMuonTrack[0]->Pt(),fMuonTrack[0]->Eta());
336     for (Int_t iTrMuon2=iTrMuon1+1; iTrMuon2<tracksMuonArm->GetEntriesFast(); iTrMuon2++) {
337       fMuonTrack[1] = (AliAODTrack*) tracksMuonArm -> At(iTrMuon2);
338       FillHistograms(centBin, kSingleEvent);
339     }
340   }
341
342   // Mixed event
343   AliEventPool* pool = fPoolMgr->GetEventPool(percentile, zvtx);
344   //pool->PrintInfo();
345   if (pool->IsReady() || pool->NTracksInPool() > 2000 || pool->GetCurrentNEvents() >= 5) {
346     for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {
347       TObjArray *mixedTracks = pool->GetEvent(jMix);
348       fHistNMuons_vs_NMuons_Mixed[centBin]->Fill(mixedTracks->GetEntriesFast(), tracksMuonArm->GetEntriesFast());
349       for (Int_t iTrMuon1=0; iTrMuon1<tracksMuonArm->GetEntriesFast(); iTrMuon1++) {
350         fMuonTrack[0] = (AliAODTrack*) tracksMuonArm->At(iTrMuon1);
351         fHistSingleMuonsPt_Mixed[centBin]->Fill(fMuonTrack[0]->Pt());
352         fHistSingleMuonsEtaPt_Mixed[centBin]->Fill(fMuonTrack[0]->Pt(),fMuonTrack[0]->Eta());
353         for (Int_t iTrMuon2=0; iTrMuon2<mixedTracks->GetEntriesFast(); iTrMuon2++) {
354           fMuonTrack[1] = (AliAODTrack*) mixedTracks -> At(iTrMuon2);
355           FillHistograms(centBin, kMixedEvent);
356         }
357       }
358     }
359   }
360   pool->UpdatePool(tracksMuonArm);
361
362   delete tracksMuonArm;
363
364   PostData(1, fOutputList); 
365
366 }
367
368 //====================================================================================================================================================
369
370 void AliAnalysisTaskDiMuonCorrelations::FillHistograms(Int_t centrality, Int_t option) {
371
372   Int_t ptBinTrackMuon1 = fPtAxis -> FindBin(fMuonTrack[0]->Pt());
373   Int_t ptBinTrackMuon2 = fPtAxis -> FindBin(fMuonTrack[1]->Pt());
374
375   if (ptBinTrackMuon1<1 || ptBinTrackMuon1>fNbinsPt || ptBinTrackMuon2<1 || ptBinTrackMuon2>fNbinsPt) return;
376
377   Double_t deltaPhi = fMuonTrack[0]->Phi() - fMuonTrack[1]->Phi();
378   if (deltaPhi >  1.5*TMath::Pi()) deltaPhi -= TMath::TwoPi();
379   if (deltaPhi < -0.5*TMath::Pi()) deltaPhi += TMath::TwoPi();
380
381   if (option==kSingleEvent)     fHistDeltaPhi[centrality][ptBinTrackMuon1-1][ptBinTrackMuon2-1]    -> Fill(TMath::RadToDeg()*deltaPhi);
382   else if (option==kMixedEvent) fHistDeltaPhiMix[centrality][ptBinTrackMuon1-1][ptBinTrackMuon2-1] -> Fill(TMath::RadToDeg()*deltaPhi);
383
384   if (option==kSingleEvent)     fHistEtaDeltaPhi[centrality][ptBinTrackMuon1-1][ptBinTrackMuon2-1]    -> Fill(TMath::RadToDeg()*deltaPhi,fMuonTrack[0]->Eta());
385   else if (option==kMixedEvent) fHistEtaDeltaPhiMix[centrality][ptBinTrackMuon1-1][ptBinTrackMuon2-1] -> Fill(TMath::RadToDeg()*deltaPhi,fMuonTrack[0]->Eta());
386
387   if (option==kSingleEvent)     fHistTracksEtavsEta[centrality]->Fill(fMuonTrack[0]->Eta(),fMuonTrack[1]->Eta());
388   else if (option==kMixedEvent) fHistTracksEtavsEta_Mixed[centrality]->Fill(fMuonTrack[0]->Eta(),fMuonTrack[1]->Eta());
389
390 }
391
392 //====================================================================================================================================================
393
394 Bool_t AliAnalysisTaskDiMuonCorrelations::IsTriggerFired() {
395   
396   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); 
397
398   return isSelected;
399 }
400
401 //====================================================================================================================================================
402
403 Float_t AliAnalysisTaskDiMuonCorrelations::GetV0Multiplicity() {
404   
405   Float_t multiplicity=0;
406   for (Int_t iChannel=0; iChannel<64; iChannel++) multiplicity += fAOD->GetVZEROData()->GetMultiplicity(iChannel);
407   return multiplicity;
408
409 }
410
411 //====================================================================================================================================================
412
413 TObjArray* AliAnalysisTaskDiMuonCorrelations::GetAcceptedTracksMuonArm(AliAODEvent *aodEvent) {
414
415   // fills the array of muon tracks that pass the cuts
416
417   TObjArray *tracks = new TObjArray;
418   tracks->SetOwner(kFALSE);
419
420   Int_t nTracks = aodEvent->GetNTracks();
421
422   AliAODTrack *track = 0;
423   
424   for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
425     track = aodEvent->GetTrack(iTrack);
426     if (track->IsMuonTrack() && track->GetMatchTrigger()>=fTriggerMatchLevelMuon) {
427       if (track->Chi2perNDF() < fMaxChi2Muon) {
428         if (track->Eta() > fMinEtaMuon && track->Eta() < fMaxEtaMuon) {
429           tracks->Add(new AliAODTrack(*track));
430         }
431       }
432     }
433   }
434
435   return tracks;
436
437 }
438
439 //====================================================================================================================================================
440
441 void AliAnalysisTaskDiMuonCorrelations::SetCentBinning(Int_t nBins, Double_t *limits) {
442
443   if (nBins>fNMaxBinsCentrality) {
444     AliInfo(Form("WARNING : only %d centrality bins (out of the %d proposed) will be considered",fNMaxBinsCentrality,nBins));
445     nBins = fNMaxBinsCentrality;
446   }
447   if (nBins<=0) {
448     AliInfo("WARNING : at least one centrality bin must be considered");
449     nBins = 1;
450   }
451   
452   fNbinsCent = nBins;
453   fCentAxis  = new TAxis(fNbinsCent, limits);
454
455 }
456
457 //====================================================================================================================================================
458
459 void AliAnalysisTaskDiMuonCorrelations::SetPtBinning(Int_t nBins, Double_t *limits) {
460
461   if (nBins>fNMaxBinsPt) {
462     AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsPt,nBins));
463     nBins = fNMaxBinsPt;
464   }
465   if (nBins<=0) {
466     AliInfo("WARNING : at least one pt bin must be considered");
467     nBins = 1;
468   }
469   
470   fNbinsPt = nBins;
471   fPtAxis  = new TAxis(fNbinsPt, limits);
472
473 }
474
475 //====================================================================================================================================================
476
477 void AliAnalysisTaskDiMuonCorrelations::SetEtaBinning(Int_t nBins, Double_t *limits) {
478
479   if (nBins>fNMaxBinsEta) {
480     AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsEta,nBins));
481     nBins = fNMaxBinsEta;
482   }
483   if (nBins<=0) {
484     AliInfo("WARNING : at least one pt bin must be considered");
485     nBins = 1;
486   }
487   
488   fEtaAxis  = new TAxis(nBins, limits);
489
490 }
491
492 //====================================================================================================================================================
493
494 Int_t AliAnalysisTaskDiMuonCorrelations::GetCentBin() {
495
496   Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
497
498   Int_t bin = fCentAxis->FindBin(percentile) - 1;
499   if (bin >= fNbinsCent) bin = -1;
500   return bin;
501   
502 }
503
504 //====================================================================================================================================================
505
506 Double_t AliAnalysisTaskDiMuonCorrelations::GetITSMultiplicity() {
507
508   Double_t multiplicity = fAOD->GetHeader()->GetNumberOfITSClusters(1);
509
510   return multiplicity;
511
512 }
513
514 //====================================================================================================================================================
515
516 void AliAnalysisTaskDiMuonCorrelations::Terminate(Option_t *) {
517
518   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
519   if (!fOutputList) {
520     printf("ERROR: Output list not available\n");
521     return;
522   }  
523
524 }
525
526 //====================================================================================================================================================
527