Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskMuonCollisionMultiplicity.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /// Compute the number of Muons tracks as a function of the SPD tracklets multiplicity
17 /// Compare with Monte Carlo tracks
18 /// Author Matthieu LENHARDT - SUBATECH, Nantes
19
20
21 //PWG3/muon includes
22 #include "AliAnalysisTaskMuonCollisionMultiplicity.h"
23
24 //STEER includes
25 #include "AliAODEvent.h"
26 #include "AliAODVertex.h"
27 #include "AliESDEvent.h"
28 #include "AliESDVertex.h"
29 #include "AliAODTrack.h"
30 #include "AliESDMuonTrack.h"
31 #include "AliAODTracklets.h"
32 #include "AliAODInputHandler.h"
33 #include "AliAnalysisManager.h"
34 #include "AliAODDimuon.h"
35 #include "AliMCEvent.h"
36 #include "AliMCParticle.h"
37 #include "AliMultiplicity.h"
38
39 //ROOT includes
40 #include <TH2D.h>
41 #include <THnSparse.h>
42 #include <TChain.h>
43 #include <TList.h>
44 #include <TArrayD.h>
45 #include <Riostream.h>
46 #include <TParticle.h>
47 #include <TLorentzVector.h>
48
49 //___________________________________________________
50 AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity()
51   :
52   AliAnalysisTaskSE(),
53   fIsInit(kFALSE),
54   fAOD(0x0),
55   fESD(0x0),
56   fEtaCut(0),
57   fTrackletMultiplicity(0),
58   fTriggerList(0),
59   fSingleMuonList(0),
60   fDimuonList(0),
61   fMonteCarloList(0)
62 {
63   ///Default Constructor
64 }
65
66
67 //___________________________________________________
68 AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity(const AliAnalysisTaskMuonCollisionMultiplicity& src)
69   :
70   AliAnalysisTaskSE(),
71   fIsInit(kFALSE),
72   fAOD(0x0),
73   fESD(0x0),
74   fEtaCut(0),
75   fTrackletMultiplicity(0),
76   fTriggerList(0),
77   fSingleMuonList(0),
78   fDimuonList(0),
79   fMonteCarloList(0)
80 {
81   /// copy ctor
82   src.Copy(*this);
83
84 }
85
86 //___________________________________________________
87 AliAnalysisTaskMuonCollisionMultiplicity& AliAnalysisTaskMuonCollisionMultiplicity::operator=(const AliAnalysisTaskMuonCollisionMultiplicity& src)
88 {
89   /// assignement operator
90   if ( this != &src ) 
91   {
92     src.Copy(*this);
93   }
94   return *this;
95 }
96
97
98
99
100 //___________________________________________________
101 AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity(const Char_t *name)
102   :
103   AliAnalysisTaskSE(name),
104   fIsInit(kFALSE),
105   fAOD(0x0),
106   fESD(0x0),
107   fEtaCut(0),
108   fTrackletMultiplicity(0),
109   fTriggerList(0),
110   fSingleMuonList(0),
111   fDimuonList(0),
112   fMonteCarloList(0)
113 {
114   // Define Inputs and outputs
115   //DefineInput(0, TChain::Class());
116   //DefineInput(1, TChain::Class());
117   DefineOutput(0, TList::Class());
118   DefineOutput(1, TList::Class());
119   DefineOutput(2, TList::Class());
120   DefineOutput(3, TList::Class());
121   DefineOutput(4, TList::Class());
122 }
123
124
125 //______________________________________________________________________________
126 AliAnalysisTaskMuonCollisionMultiplicity::~AliAnalysisTaskMuonCollisionMultiplicity()
127 {
128 // Destructor.
129   delete fAOD;
130   delete fESD;
131   delete fTriggerList;
132   delete fSingleMuonList;
133   delete fDimuonList;
134   delete fMonteCarloList;
135 }
136
137
138
139 //________________________________________________________________________
140 void AliAnalysisTaskMuonCollisionMultiplicity::UserCreateOutputObjects()
141 {
142   // Initialise the object, and open the file.
143   if (!fIsInit)
144     Init();
145   OpenFile(0);
146 }
147
148
149
150
151 //________________________________________________________________________
152 void AliAnalysisTaskMuonCollisionMultiplicity::UserExec(Option_t */*option*/)
153 {
154   // Execute the analysis task
155   fAOD = 0x0;
156   fESD = 0x0;
157
158   if (!fIsInit)
159     Init();
160
161   fAOD = dynamic_cast<AliAODEvent *> (InputEvent());
162   if (!fAOD)
163     fESD = dynamic_cast<AliESDEvent *> (InputEvent());
164   
165  
166   if (fAOD)
167     CheckEventAOD();
168
169   if (fESD)
170     CheckEventESD();
171
172   PostData(1, fTriggerList);
173   PostData(2, fSingleMuonList);
174   PostData(3, fDimuonList);
175   PostData(4, fMonteCarloList);
176 }
177
178
179
180 //________________________________________________________________________
181 void AliAnalysisTaskMuonCollisionMultiplicity::NotifyRun()
182 {
183   // Notify run
184 }
185
186
187 //________________________________________________________________________
188 void AliAnalysisTaskMuonCollisionMultiplicity::FinishTaskOutput()
189 {
190   // Finish the task
191 }
192
193
194 //__________________________________________________________________________
195 Bool_t AliAnalysisTaskMuonCollisionMultiplicity::CheckEventAOD()
196 {
197   // Check if the AOD event pass the cuts
198   AliAODVertex *vertex = fAOD->GetPrimaryVertex();
199   
200   if (!vertex)
201      return kFALSE;
202
203   if (vertex->GetNContributors() < 1)
204     return kFALSE;
205   
206   ComputeMultiplicity();
207   
208   // Variables use to determine the type of trigger :
209   // 0 for minimum bias : CINT1B, CINT1-B or MB1
210   // 1 for muon events : CMUS1B, CMUS1-B or MULow
211   // -1 for everything else
212   TString trigger = fAOD->GetFiredTriggerClasses();
213   Int_t triggerClass = -1; 
214
215   if (trigger.Contains("CINT1B") || trigger.Contains("CINT1-B") || trigger.Contains("MB1"))
216     triggerClass = 0;
217
218   if (trigger.Contains("CMUS1B") || trigger.Contains("CMUS1-B") || trigger.Contains("MULow"))
219     triggerClass = 1;
220
221   if (triggerClass >= 0)
222     FillHistosAOD(triggerClass);
223   
224   return kTRUE;
225 }
226
227
228
229 //__________________________________________________________________________
230 Bool_t AliAnalysisTaskMuonCollisionMultiplicity::CheckEventESD()
231 {
232   // Check if the ESD event pass the cuts
233   const AliESDVertex *vertex = fESD->GetPrimaryVertex();
234   
235   if (!vertex)
236     return kFALSE;
237
238   if (vertex->GetNContributors() < 1)
239     return kFALSE;
240
241   ComputeMultiplicity();
242
243   // Variables use to determine the type of trigger :
244   // 0 for minimum bias : CINT1B, CINT1-B or MB1
245   // 1 for muon events : CMUS1B, CMUS1-B or MULow
246   // -1 for everything else
247   TString trigger = fESD->GetFiredTriggerClasses();
248   Int_t triggerClass = -1; 
249
250   if (trigger.Contains("CINT1B") || trigger.Contains("CINT1-B") || trigger.Contains("MB1") || trigger.Contains("CMBAC-B") || trigger.Contains("CMBACS2-B"))
251     triggerClass = 0;
252
253   if (trigger.Contains("CMUS1B") || trigger.Contains("CMUS1-B") || trigger.Contains("MULow"))
254     triggerClass = 1;
255
256   if (triggerClass >= 0)
257     FillHistosESD(triggerClass);
258
259   if (fMCEvent)
260     FillHistosMC();
261
262   return kTRUE;
263 }
264
265
266
267
268 //________________________________________________________________________
269 void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosAOD(Int_t triggerClass)
270 {
271   // Fill histos for AOD events
272   Int_t nTracks = fAOD->GetNTracks();
273   Int_t nDimuons = fAOD->GetNDimuons();
274   
275   // Fill histos
276   Double_t vertexPosition = fAOD->GetPrimaryVertex()->GetZ();
277   Double_t pileUp = !(fAOD->IsPileupFromSPD());
278   
279   Double_t valuesTrigger[3] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp};
280   ((THnSparseD *)fTriggerList->At(triggerClass))->Fill(valuesTrigger);
281   
282
283   // Loop on the muons tracks
284   for (Int_t ii = 0; ii < nTracks; ii++)
285     if (IsUsableMuon(fAOD->GetTrack(ii)))
286       {
287         Double_t matchTrigger = fAOD->GetTrack(ii)->GetMatchTrigger();
288         if (matchTrigger > 1.0)
289           matchTrigger = 1.0;                                            // We don't care what type of trigger it is (low or high pT)
290         
291         Double_t thetaAbs = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetTrack(ii)->GetRAtAbsorberEnd()/505.0);
292         Double_t eta = fAOD->GetTrack(ii)->Eta();
293         Double_t p = fAOD->GetTrack(ii)->P();
294         Double_t pT = fAOD->GetTrack(ii)->Pt();
295         // For the p used in the pDCA, we want the mean between the p before the muon go through the absorber (p corrected) and after (p uncorrected)
296         // However p uncorrected is not saved in the AODs
297         // Instead we define p uncorrected as p corrected minus the mean p lost when a muon go through the absorber
298         // p lost for muons (it depends of theta_abs) :
299         // 2.0 < theta_abs < 3.0 ---> 2.98 GeV
300         // 3.0 < theta_abs < 10.0 --->  2.4 GeV
301         // No correction applied otherwise
302         Double_t pDCA = p*fAOD->GetTrack(ii)->DCA();
303         if (2.0 < thetaAbs && thetaAbs < 3.0)
304           pDCA = (p-2.98/2.0) * fAOD->GetTrack(ii)->DCA();
305         if (3.0 < thetaAbs && thetaAbs < 10.0)
306           pDCA = (p-2.4/2.0) * fAOD->GetTrack(ii)->DCA();
307         
308         Double_t valuesMuon[9] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp, matchTrigger, thetaAbs, eta, pDCA, static_cast<Double_t>(p), static_cast<Double_t>(pT)};
309         ((THnSparseD *)fSingleMuonList->At(triggerClass))->Fill(valuesMuon);
310       }
311   
312   // Loop on Dimuons
313   for (Int_t ii = 0; ii < nDimuons; ii++)
314     if (fAOD->GetDimuon(ii)->Charge() == 0.0)
315       if (IsUsableMuon(fAOD->GetDimuon(ii)->GetMu(0)))
316         if (IsUsableMuon(fAOD->GetDimuon(ii)->GetMu(1)))
317           {
318             Double_t matchTrigger1 = fAOD->GetDimuon(ii)->GetMu(0)->GetMatchTrigger();
319             if (matchTrigger1 > 0.0)
320               matchTrigger1 = 1.0;
321             Double_t matchTrigger2 = fAOD->GetDimuon(ii)->GetMu(1)->GetMatchTrigger();
322             if (matchTrigger2 > 0.0)
323               matchTrigger2 = 1.0;
324             
325             Double_t nMatchTrigger = matchTrigger1 + matchTrigger2;
326             
327             Double_t thetaAbs1 = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetDimuon(ii)->GetMu(0)->GetRAtAbsorberEnd()/505.0);
328             Double_t thetaAbs2 = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetDimuon(ii)->GetMu(1)->GetRAtAbsorberEnd()/505.0);
329             Double_t eta1 = fAOD->GetDimuon(ii)->GetMu(0)->Eta();
330             Double_t eta2 = fAOD->GetDimuon(ii)->GetMu(1)->Eta();
331             
332             Double_t p1 = fAOD->GetDimuon(ii)->GetMu(0)->P();
333             Double_t p2 = fAOD->GetDimuon(ii)->GetMu(1)->P();
334             // See the explanation on how the pDCA is computed in the single muon loop
335             Double_t pDCA1 = p1*fAOD->GetDimuon(ii)->GetMu(0)->DCA();
336             if (2.0 < thetaAbs1 && thetaAbs1 < 3.0)
337               pDCA1 = (p1-2.98/2.0) * fAOD->GetDimuon(ii)->GetMu(0)->DCA();
338             if (3.0 < thetaAbs1 && thetaAbs1 < 10.0)
339               pDCA1 = (p1-2.4/2.0) * fAOD->GetDimuon(ii)->GetMu(0)->DCA();
340             Double_t pDCA2 = p2*fAOD->GetDimuon(ii)->GetMu(1)->DCA();
341             if (2.0 < thetaAbs2 && thetaAbs2 < 3.0)
342               pDCA2 = (p2-2.98/2.0) * fAOD->GetDimuon(ii)->GetMu(1)->DCA();
343             if (3.0 < thetaAbs2 && thetaAbs2 < 10.0)
344               pDCA2 = (p2-2.4/2.0) * fAOD->GetDimuon(ii)->GetMu(1)->DCA();
345             
346             Double_t y = fAOD->GetDimuon(ii)->Y();
347             Double_t p = fAOD->GetDimuon(ii)->P();
348             Double_t pT = fAOD->GetDimuon(ii)->Pt();
349             Double_t M = fAOD->GetDimuon(ii)->M();
350             
351             Double_t valuesDimuon[18] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp, matchTrigger1, matchTrigger2, static_cast<Double_t>(nMatchTrigger), thetaAbs1, thetaAbs2,
352                                          eta1, eta2, pDCA1, pDCA2, static_cast<Double_t>(p1), p2, y, p, pT, M};
353             ((THnSparseD *)fDimuonList->At(triggerClass))->Fill(valuesDimuon);
354           }
355 }
356
357
358
359 void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosESD(Int_t triggerClass)
360 {
361   // Fill the histos for ESD events 
362   Int_t nTracks = fESD->GetNumberOfMuonTracks();
363
364   Double_t vertexPosition = fESD->GetPrimaryVertex()->GetZ();
365   Double_t pileUp = !(fESD->IsPileupFromSPD());
366   
367   Double_t valuesTrigger[3] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp};
368   ((THnSparseD *)fTriggerList->At(triggerClass))->Fill(valuesTrigger);
369   
370
371   // Loop on the muons tracks
372   for (Int_t ii = 0; ii < nTracks; ii++)
373     if (IsUsableMuon(fESD->GetMuonTrack(ii)))
374       {
375         Double_t matchTrigger1 = fESD->GetMuonTrack(ii)->GetMatchTrigger();
376         if (matchTrigger1 > 1.0)
377           matchTrigger1 = 1.0;                                            // We don't care what type of trigger it is (low or high pT)
378         
379         Double_t thetaAbs1 = (180.0 / TMath::Pi()) * TMath::ATan(fESD->GetMuonTrack(ii)->GetRAtAbsorberEnd()/505.0);
380         Double_t eta1 = fESD->GetMuonTrack(ii)->Eta();
381         Double_t p1 = fESD->GetMuonTrack(ii)->P();
382         Double_t pT1 = fESD->GetMuonTrack(ii)->Pt();
383
384         // For the p used in the pDCA, we want the mean between the p before the muon go through the absorber (p corrected) and after (p uncorrected)
385         // These are respectively AliESDMuonTrack::P() and AliESDMuonTrack::PUncorrected()
386         Double_t pUncor1 = fESD->GetMuonTrack(ii)->PUncorrected();
387         Double_t pDCA1 = (p1+pUncor1) * fESD->GetMuonTrack(ii)->GetDCA() / 2.0;
388         
389         Double_t valuesMuon[9] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp, matchTrigger1, thetaAbs1, eta1, pDCA1, p1, pT1};
390         ((THnSparseD *)fSingleMuonList->At(triggerClass))->Fill(valuesMuon);
391         
392         // Second loop on muons, to fill the dimuons histos
393         for (Int_t jj = ii+1; jj < nTracks; jj++)
394           if (IsUsableMuon(fESD->GetMuonTrack(jj)))
395             if (fESD->GetMuonTrack(ii)->Charge() + fESD->GetMuonTrack(jj)->Charge() == 0.0)
396               {
397                 Double_t matchTrigger2 = fESD->GetMuonTrack(jj)->GetMatchTrigger();
398                 if (matchTrigger2 > 0.0)
399                   matchTrigger2 = 1.0;
400                 
401                 Double_t nMatchTrigger = matchTrigger1 + matchTrigger2;
402
403                 Double_t thetaAbs2 = (180.0 / TMath::Pi()) * TMath::ATan(fESD->GetMuonTrack(jj)->GetRAtAbsorberEnd()/505.0);
404                 Double_t eta2 = fESD->GetMuonTrack(jj)->Eta();
405                 Double_t p2 = fESD->GetMuonTrack(jj)->P();
406
407                 // For the p used in the pDCA, we want the mean between the p before the muon go through the absorber (p corrected) and after (p uncorrected)
408                 // These are respectively AliESDMuonTrack::P() and AliESDMuonTrack::PUncorrected()
409                 Double_t pUncor2 = fESD->GetMuonTrack(jj)->PUncorrected();
410                 Double_t pDCA2 = (p2+pUncor2) * fESD->GetMuonTrack(jj)->GetDCA() / 2.0;
411                 
412                 // To compute the p, pT and M of the dimuon, we need a TLorentz vector of the dimuon
413                 Double_t E = fESD->GetMuonTrack(ii)->E() + fESD->GetMuonTrack(jj)->E();
414                 Double_t pX = fESD->GetMuonTrack(ii)->Px() + fESD->GetMuonTrack(jj)->Px();
415                 Double_t pY = fESD->GetMuonTrack(ii)->Py() + fESD->GetMuonTrack(jj)->Py();
416                 Double_t pZ = fESD->GetMuonTrack(ii)->Pz() + fESD->GetMuonTrack(jj)->Pz();
417                 TLorentzVector *dimuonVector = new TLorentzVector(pX, pY, pZ, E);
418                 dimuonVector->SetPxPyPzE(pX, pY, pZ, E);
419                 
420                 Double_t y = dimuonVector->Rapidity();
421                 Double_t p = dimuonVector->P();
422                 Double_t pT = TMath::Sqrt(pX*pX + pY*pY);
423                 Double_t M = dimuonVector->M();
424                 
425                 Double_t valuesDimuon[18] = {static_cast<Double_t>(fTrackletMultiplicity), vertexPosition, pileUp, matchTrigger1, matchTrigger2, static_cast<Double_t>(nMatchTrigger), thetaAbs1, thetaAbs2,
426                                              eta1, eta2, pDCA1, pDCA2, p1, p2, y, p, pT, M};
427                 ((THnSparseD *)fDimuonList->At(triggerClass))->Fill(valuesDimuon);
428                 delete dimuonVector;
429               }
430       }
431   
432
433   // Since this is the ESD, we are also going to fill a the V0amp vs multiplicity histos
434   if (triggerClass == 0)
435     {
436       AliESDVZERO *v0 = fESD->GetVZEROData();
437       Float_t multV0 = 0;
438       Float_t multV0A = 0;
439       Float_t multV0C = 0;
440
441       for (Int_t ii = 0; ii < 64; ii++)
442         {
443           multV0 += v0->GetMultiplicity(ii);
444           if (ii < 32)
445             {
446               multV0C += v0->GetMultiplicityV0C(ii);
447               multV0A += v0->GetMultiplicityV0A(ii);
448             }
449         }
450
451       ((TH2D *)fTriggerList->At(2))->Fill(fTrackletMultiplicity, multV0);
452       ((TH2D *)fTriggerList->At(3))->Fill(fTrackletMultiplicity, multV0A);
453       ((TH2D *)fTriggerList->At(4))->Fill(fTrackletMultiplicity, multV0C);
454     }
455
456
457 }
458
459
460 //________________________________________________________________________
461 void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosMC()
462 {
463   // Fill the histo of the correlation between MC tracks and ESD tracks
464
465   Int_t multiplicityGenerated = 0;
466   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
467
468   for (Int_t nn = 0; nn < fMCEvent->GetNumberOfTracks(); nn++)
469     {
470       AliMCParticle *particle = (AliMCParticle *) fMCEvent->GetTrack(nn);
471       Bool_t isGoodMult = kTRUE;
472       
473       if (particle->Particle()->GetStatusCode() != 1)
474         isGoodMult = kFALSE;
475       
476       if (particle->Charge() == 0)
477         isGoodMult = kFALSE;
478
479       if (TMath::Abs(particle->Eta()) > 1.6)
480         isGoodMult = kFALSE;
481
482       // Check if the particle is a pion, kaon, proton, electron or muon
483       if (TMath::Abs(particle->PdgCode()) != 211 && TMath::Abs(particle->PdgCode()) != 321 && TMath::Abs(particle->PdgCode()) != 2212 &&
484           TMath::Abs(particle->PdgCode()) != 11 && TMath::Abs(particle->PdgCode()) != 13)
485         isGoodMult = kFALSE;
486
487       // Check if the distance to vertex is inferior to 1 cm
488       Double_t distanceToVertex = TMath::Sqrt((particle->Xv() - vertex->GetXv())*(particle->Xv() - vertex->GetXv()) + 
489                                               (particle->Yv() - vertex->GetYv())*(particle->Yv() - vertex->GetYv()) + 
490                                               (particle->Zv() - vertex->GetZv())*(particle->Zv() - vertex->GetZv()));
491       if (distanceToVertex > 1.0)
492         isGoodMult = kFALSE;
493         
494       if (isGoodMult)
495         multiplicityGenerated += 1;
496     }
497
498   ((TH2D *)fMonteCarloList->At(0))->Fill(multiplicityGenerated, fTrackletMultiplicity);
499 }
500
501
502
503 //________________________________________________________________________
504 void AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
505 {
506   // Compute the collision multiplicity based on AOD or ESD tracklets
507
508   Int_t multiplicity = 0;
509   
510   if (fAOD)
511     {
512       AliAODTracklets *tracklets = fAOD->GetTracklets();
513       Int_t nTracklets = tracklets->GetNumberOfTracklets();
514       for (Int_t nn = 0; nn < nTracklets; nn++)
515         {
516           Double_t theta = tracklets->GetTheta(nn);
517           Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
518           
519           if (TMath::Abs(eta) < fEtaCut)
520             multiplicity += 1;
521         }
522     }
523
524
525   if (fESD)
526     {
527       // In ESDs, we use the EstimateMultiplicity function
528       Int_t multTracklets = 0;
529       Int_t multTPC = 0;
530       Int_t multITSSA = 0;
531       fESD->EstimateMultiplicity(multTracklets, multTPC, multITSSA, fEtaCut);
532       multiplicity = multTracklets;
533     }
534
535   fTrackletMultiplicity =  multiplicity;
536 }
537
538
539
540 //________________________________________________________________________
541 Bool_t AliAnalysisTaskMuonCollisionMultiplicity::IsUsableMuon(AliAODTrack *track)
542 {
543   // Check if the track is a usable muon track
544   // Cuts applied :
545   // - is it a muon track?
546   // - does it have a pT > 0.0?
547
548   Bool_t isGood = kFALSE;
549
550   if (!track->IsMuonTrack())
551     return isGood;
552
553   if (!(track->Pt() > 0.0))
554     return isGood;
555
556   isGood = kTRUE;
557   return isGood;
558 }
559
560
561
562 //________________________________________________________________________
563 Bool_t AliAnalysisTaskMuonCollisionMultiplicity::IsUsableMuon(AliESDMuonTrack *track)
564 {
565   // Check if the track is a usable muon track
566   // Cuts applied :
567   // - is it a muon track?
568   // - does it have a pT > 0.0?
569
570   Bool_t isGood = kFALSE;
571
572   if (!track->ContainTrackerData())
573     return isGood;
574
575   if (!(track->Pt() > 0.0))
576     return isGood;
577
578   isGood = kTRUE;
579   return isGood;
580 }
581
582
583
584 //________________________________________________________________________
585 void AliAnalysisTaskMuonCollisionMultiplicity::Init()
586 {
587   // Initialize the object
588
589   fTriggerList = new TList();
590   fSingleMuonList = new TList();
591   fDimuonList = new TList();
592   fMonteCarloList = new TList();
593
594   fTriggerList->SetOwner();
595   fSingleMuonList->SetOwner();
596   fDimuonList->SetOwner();
597   fMonteCarloList->SetOwner();
598
599
600
601   // Trigger histos
602   // dimension 0 : multiplicity of the event
603   // dimension 1 : z vertex of the event
604   // dimension 2 : is it an event without pile up (0 for no, 1 for yes)?
605   Int_t nBinsTrigger[3] =       {  150,    60,   2};
606   Double_t minRangeTrigger[3] = {  0.0, -30.0, 0.0};
607   Double_t maxRangeTrigger[3] = {150.0,  30.0, 2.0};
608   THnSparseD *CINT1B = new THnSparseD ("CINT1B", "CINT1B", 3, nBinsTrigger, minRangeTrigger, maxRangeTrigger);
609   THnSparseD *CMUS1B = new THnSparseD ("CMUS1B", "CMUS1B", 3, nBinsTrigger, minRangeTrigger, maxRangeTrigger);
610   CINT1B->Sumw2();
611   CMUS1B->Sumw2();
612
613   TH2D *CompSPDV0 = new TH2D ("CompSPDV0", "CompSPDV0", 150, 0.0, 150.0, 2000, 0.0, 2000.0);
614   CompSPDV0->Sumw2();
615   TH2D *CompSPDV0A = new TH2D ("CompSPDV0A", "CompSPDV0A", 150, 0.0, 150.0, 1000, 0.0, 1000.0);
616   CompSPDV0A->Sumw2();
617   TH2D *CompSPDV0C = new TH2D ("CompSPDV0C", "CompSPDV0C", 150, 0.0, 150.0, 1000, 0.0, 1000.0);
618   CompSPDV0C->Sumw2();
619
620   fTriggerList->AddAt(CINT1B, 0);
621   fTriggerList->AddAt(CMUS1B, 1);
622   fTriggerList->AddAt(CompSPDV0, 2);
623   fTriggerList->AddAt(CompSPDV0A, 3);
624   fTriggerList->AddAt(CompSPDV0C, 4);
625
626
627
628
629   // Muons histos
630   // dimension 0 : multiplicity of the event
631   // dimension 1 : z vertex of the event
632   // dimension 2 : is it an event without pile up (0 for no, 1 for yes)?
633   // dimension 3 : does the muon match the trigger (0 for no, 1 for yes)?
634   // dimension 4 : theta_abs of the muon
635   // dimension 5 : eta of the muon
636   // dimension 6 : p DCA of the muon
637   // dimension 7 : p of the muon
638   // dimension 8 : pT of the muon
639
640   Int_t nBinsMuon[9] =       {  150,    60,   2,   2,  110,   35,   300,   300,  300};
641   Double_t minRangeMuon[9] = {  0.0, -30.0, 0.0, 0.0,  0.0, -5.0,   0.0,   0.0,  0.0};
642   Double_t maxRangeMuon[9] = {150.0,  30.0, 2.0, 2.0, 11.0, -1.5, 300.0, 150.0, 30.0};
643
644   THnSparseD *muonCINT1B = new THnSparseD("muonCINT1B", "muonCINT1B", 9, nBinsMuon, minRangeMuon, maxRangeMuon);
645   THnSparseD *muonCMUS1B = new THnSparseD("muonCMUS1B", "muonCMUS1B", 9, nBinsMuon, minRangeMuon, maxRangeMuon);
646   muonCINT1B->Sumw2();
647   muonCMUS1B->Sumw2();
648
649   fSingleMuonList->AddAt(muonCINT1B, 0);
650   fSingleMuonList->AddAt(muonCMUS1B, 1);
651
652
653   // Dimuons histos
654   // dimension 0  : multiplicity of the event
655   // dimension 1  : z range (0 for no, 1 for yes)?
656   // dimension 2  : is it an event without pile up (0 for no, 1 for yes)?
657   // dimension 3  : does the first muon match the trigger (0 for no, 1 for yes)?
658   // dimension 4  : does the second muon match the trigger (0 for no, 1 for yes)?
659   // dimension 5  : number of muons matching the trigger in the dimuon
660   // dimension 6  : theta_abs of the first muon
661   // dimension 7  : theta_abs of the second muon
662   // dimension 8  : eta of the first muon
663   // dimension 9  : eta of the second muon
664   // dimension 10 : p DCA of the first muon
665   // dimension 11 : p DCA of the second muon
666   // dimension 12 : p of the first muon
667   // dimension 13 : p of the second muon
668   // dimension 14 : y of the dimuon
669   // dimension 15 : p of the dimuon
670   // dimension 16 : pT of the dimuon
671   // dimension 17 : invariant mass of the dimuon
672
673   Int_t nBinsDimuon[18] =       {  150,    60,   2,   2,   2,   3,  110,  110,   35,   35,   200,   200,   300,   300,   35,   300,  300,  375};
674   Double_t minRangeDimuon[18] = {  0.0, -30.0, 0.0, 0.0, 0.0, 0.0,  0.0,  0.0, -5.0, -5.0,   0.0,   0.0,   0.0,   0.0, -5.0,   0.0,  0.0,  0.0};
675   Double_t maxRangeDimuon[18] = {150.0,  30.0, 2.0, 2.0, 2.0, 3.0, 11.0, 11.0, -1.5, -1.5, 600.0, 600.0, 150.0, 150.0, -1.5, 300.0, 30.0, 15.0};
676
677   THnSparseD *dimuonCINT1B = new THnSparseD("dimuonCINT1B", "dimuonCINT1B", 18, nBinsDimuon, minRangeDimuon, maxRangeDimuon);
678   THnSparseD *dimuonCMUS1B = new THnSparseD("dimuonCMUS1B", "dimuonCMUS1B", 18, nBinsDimuon, minRangeDimuon, maxRangeDimuon);
679
680   dimuonCINT1B->Sumw2();
681   dimuonCMUS1B->Sumw2();
682
683   fDimuonList->AddAt(dimuonCINT1B, 0);
684   fDimuonList->AddAt(dimuonCMUS1B, 1);
685
686   // MonteCarlo Histo
687   TH2D *correlGenerReco = new TH2D("correlGenerReco", "correlGenerReco", 250, 0.0, 250.0, 250, 0.0, 250.0);
688   correlGenerReco->GetXaxis()->SetTitle("N ch gener");
689   correlGenerReco->GetYaxis()->SetTitle("N reco tracklets");
690
691   correlGenerReco->Sumw2();
692
693   fMonteCarloList->AddAt(correlGenerReco, 0);
694
695   fIsInit = kTRUE;
696 }
697
698
699
700
701
702 //________________________________________________________________________
703 void AliAnalysisTaskMuonCollisionMultiplicity::Terminate(Option_t */*option*/)
704 {
705 //Terminate analysis
706   
707   fTriggerList = (TList *) GetOutputData(0);
708   fSingleMuonList = (TList *) GetOutputData(1);
709   fDimuonList = (TList *) GetOutputData(2);
710 }