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 }