ANALYSIS/AliMuonAnalysis.cxx

   1 #include "AliMuonAnalysis.h"
   2 //________________________________
   3 ///////////////////////////////////////////////////////////
   4 //
   5 // class AliMuonAnalysis
   6 //
   7 // MUON Analysis
   8 //
   9 //
  10 //
  11 // finck@subatech.in2p3.fr
  12 //
  13 ///////////////////////////////////////////////////////////
  14 /*********************************************************/
  15
  16 #include <TString.h>
  17 #include <TParticle.h>
  18
  19 #include <AliStack.h>
  20 #include <AliAOD.h>
  21 #include <AliAODParticle.h>
  22 #include <AliAODParticleCut.h>
  23
  24 #include <AliESDtrack.h>
  25 #include <AliESD.h>
  26
  27 #include "TFile.h"
  28 #include "TH1.h"
  29 #include "TH2.h"
  30
  31 ClassImp(AliMuonAnalysis)
  32
  33 AliMuonAnalysis::AliMuonAnalysis():
  34  fPartCut(0x0)
  35 {
  36  //ctor
  37 }
  38
  39 /*********************************************************/
  40 AliMuonAnalysis::~AliMuonAnalysis()
  41 {
  42  //dtor
  43   delete fPartCut;
  44   delete fHistoFile;
  45   delete fHPtMuon;
  46   delete fHPtMuonPlus;
  47   delete fHPtMuonMinus;
  48   delete fHPMuon;
  49   delete fHInvMassAll;
  50   delete fHRapMuon;
  51   delete fHRapResonance;
  52   delete fHPtResonance;
  53   delete fHInvMassAllvsPt;
  54 }
  55 /*********************************************************/
  56
  57 Int_t AliMuonAnalysis::Init()
  58 {
  59   //Initilizes analysis
  60   Info("Init","Histo initialized for MUON Analysis");
  61
  62   fHistoFile         = new TFile("MUONmassPlot.root", "RECREATE");
  63   fHPtMuon           = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
  64   fHPMuon            = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
  65   fHPtMuonPlus       = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
  66   fHPtMuonMinus      = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
  67   fHInvMassAll       = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
  68   fHRapMuon          = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
  69   fHRapResonance     = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
  70   fHPtResonance      = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
  71   fHInvMassAllvsPt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
  72
  73   return 0;
  74 }
  75 /*********************************************************/
  76
  77 Int_t AliMuonAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
  78 {
  79   //
  80   // process the event
  81   //
  82   if (aodrec) {
  83     GetInvMass(aodrec);
  84     //    Info("ProcessEvent","Inv Mass Rec");
  85   }
  86
  87   if (aodsim) {
  88     //     Info("ProcessEvent","aodsim not implemented");
  89   }
  90
  91   return 0;
  92
  93 }
  94
  95 /*********************************************************/
  96
  97 Int_t AliMuonAnalysis::Finish()
  98 {
  99   //Finish analysis and writes results
 100   Info("Finish","Histo writing for MUON Analysis");
 101
 102   fHistoFile->Write();
 103   fHistoFile->Close();
 104
 105   return 0;
 106 }
 107 /*********************************************************/
 108
 109 void AliMuonAnalysis::GetInvMass(AliAOD* aod)
 110 {
 111   // get the invariant mass distribution
 112   // from the oad events
 113
 114   TLorentzVector lorV1, lorV2, lorVtot;
 115   Float_t massMin = 9.17;
 116   Float_t massMax = 9.77;
 117   Int_t charge1, charge2;
 118
 119 //returns flow parameters: v2 and event plane
 120   if (aod == 0x0) {
 121      Error("AliMuonAnalysis::GetInvMass","Pointer to AOD is NULL");
 122      return;
 123   }
 124
 125   Int_t nPart = aod->GetNumberOfParticles();
 126
 127   for (Int_t iPart1 = 0; iPart1 < nPart; iPart1++)  {
 128     AliAODParticle* aodPart1 =  (AliAODParticle*)aod->GetParticle(iPart1);
 129
 130     if (aodPart1 == 0x0) {
 131       Error("AliMuonAnalysis::GetInvMass","Cannot get particle %d", iPart1);
 132       continue;
 133     }
 134
 135     lorV1.SetPxPyPzE(aodPart1->Px(),
 136                          aodPart1->Py(),
 137                          aodPart1->Pz(),
 138                          aodPart1->E());
 139
 140     fHPtMuon->Fill(lorV1.Pt());
 141     fHPMuon->Fill(lorV1.P());
 142
 143     charge1 = TMath::Sign(1,aodPart1->GetPdgCode());
 144
 145     if (charge1 > 0) {
 146       fHPtMuonPlus->Fill(lorV1.Pt());
 147     } else {
 148       fHPtMuonMinus->Fill(lorV1.Pt());
 149     }
 150     fHRapMuon->Fill(lorV1.Rapidity());
 151     for (Int_t iPart2 = iPart1 + 1; iPart2 < nPart; iPart2++)  {
 152
 153       AliAODParticle* aodPart2 =  (AliAODParticle*)aod->GetParticle(iPart2);
 154
 155       lorV2.SetPxPyPzE(aodPart2->Px(),
 156                          aodPart2->Py(),
 157                          aodPart2->Pz(),
 158                          aodPart2->E());
 159
 160       charge2 = TMath::Sign(1,aodPart2->GetPdgCode());
 161
 162       if ((charge1 * charge2) == -1) {
 163
 164         lorVtot = lorV1;
 165         lorVtot += lorV2;
 166         Float_t invMass = lorVtot.M();
 167
 168         fHInvMassAll->Fill(invMass);
 169         fHInvMassAllvsPt->Fill(invMass,lorVtot.Pt());
 170
 171         if (invMass > massMin && invMass < massMax) {
 172           fHRapResonance->Fill(lorVtot.Rapidity());
 173           fHPtResonance->Fill(lorVtot.Pt());
 174         }
 175       }
 176
 177     }
 178   }
 179 }