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