Version number ++
[u/mrichter/AliRoot.git] / 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 }