2D Invariant Mass plots
[u/mrichter/AliRoot.git] / MUON / MUONmassPlot_ESD.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 // ROOT includes
3 #include "TBranch.h"
4 #include "TClonesArray.h"
5 #include "TLorentzVector.h"
6 #include "TFile.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TParticle.h"
10 #include "TTree.h"
11 #include <Riostream.h>
12
13 // STEER includes
14 #include "AliRun.h"
15 #include "AliRunLoader.h"
16 #include "AliHeader.h"
17 #include "AliLoader.h"
18 #include "AliStack.h"
19 #include "AliESD.h"
20
21 // MUON includes
22 #include "AliESDMuonTrack.h"
23 #endif
24 //
25 // Macro MUONmassPlot.C for ESD
26 // Ch. Finck, Subatech, April. 2004
27 //
28
29 // macro to make invariant mass plots
30 // for combinations of 2 muons with opposite charges,
31 // from root file "MUON.tracks.root" containing the result of track reconstruction.
32 // Histograms are stored on the "MUONmassPlot.root" file.
33 // introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
34 // using Invariant Mass for rapidity.
35
36 // Arguments:
37 //   FirstEvent (default 0)
38 //   LastEvent (default 0)
39 //   ResType (default 553)
40 //      553 for Upsilon, anything else for J/Psi
41 //   Chi2Cut (default 100)
42 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
43 //   PtCutMin (default 1)
44 //      to keep only tracks with transverse momentum > PtCutMin
45 //   PtCutMax (default 10000)
46 //      to keep only tracks with transverse momentum < PtCutMax
47 //   massMin (default 9.17 for Upsilon) 
48 //      &  massMax (default 9.77 for Upsilon) 
49 //         to calculate the reconstruction efficiency for resonances with invariant mass
50 //         massMin < mass < massMax.
51
52 // Add parameters and histograms for analysis 
53
54 Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 10000,
55                   char* esdFileName = "AliESDs.root", Int_t ResType = 553, 
56                   Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
57                   Float_t massMin = 9.17,Float_t massMax = 9.77)
58 {
59   cout << "MUONmassPlot " << endl;
60   cout << "FirstEvent " << FirstEvent << endl;
61   cout << "LastEvent " << LastEvent << endl;
62   cout << "ResType " << ResType << endl;
63   cout << "Chi2Cut " << Chi2Cut << endl;
64   cout << "PtCutMin " << PtCutMin << endl;
65   cout << "PtCutMax " << PtCutMax << endl;
66   cout << "massMin " << massMin << endl;
67   cout << "massMax " << massMax << endl;
68
69  
70   //Reset ROOT and connect tree file
71   gROOT->Reset();
72
73
74   // File for histograms and histogram booking
75   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
76   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
77   TH1F *hPtMuonPlus = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
78   TH1F *hPtMuonMinus = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
79   TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
80   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
81   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
82   TH1F *hInvMassBg = new TH1F("hInvMassBg", "Mu+Mu- invariant mass BG(GeV/c2)", 480, 0., 12.);
83 TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
84 TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.);
85 TH1F *hInvMassRes;
86
87   if (ResType == 553) {
88     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
89   } else {
90     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
91   }
92
93   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
94   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
95   TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
96   TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
97   TH2F *hThetaPhiPlus = new TH2F("hThetaPhiPlus", "Theta vs Phi +", 760, -190., 190., 400, 160., 180.);
98   TH2F *hThetaPhiMinus = new TH2F("hThetaPhiMinus", "Theta vs Phi -", 760, -190., 190., 400, 160., 180.);
99
100
101   // settings
102   Int_t EventInMass = 0;
103   Float_t muonMass = 0.105658389;
104 //   Float_t UpsilonMass = 9.46037;
105 //   Float_t JPsiMass = 3.097;
106
107   Double_t thetaX, thetaY, pYZ;
108   Double_t fPxRec1, fPyRec1, fPzRec1, fE1;
109   Double_t fPxRec2, fPyRec2, fPzRec2, fE2;
110   Int_t fCharge, fCharge2;
111
112   Int_t ntrackhits, nevents;
113   Double_t fitfmin;
114
115  
116   TLorentzVector fV1, fV2, fVtot;
117   
118   // open run loader and load gAlice, kinematics and header
119   AliRunLoader* runLoader = AliRunLoader::Open(filename);
120   if (!runLoader) {
121     Error("MUONmass_ESD", "getting run loader from file %s failed", 
122             filename);
123     return kFALSE;
124   }
125
126   runLoader->LoadgAlice();
127   gAlice = runLoader->GetAliRun();
128   if (!gAlice) {
129     Error("MUONmass_ESD", "no galice object found");
130     return kFALSE;
131   }
132   
133
134   // open the ESD file
135   TFile* esdFile = TFile::Open(esdFileName);
136   if (!esdFile || !esdFile->IsOpen()) {
137     Error("MUONmass_ESD", "opening ESD file %s failed", esdFileName);
138     return kFALSE;
139   }
140   
141   runLoader->LoadHeader();
142   nevents = runLoader->GetNumberOfEvents();
143         
144   // Loop over events
145   for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) {
146
147     // get current event
148     runLoader->GetEvent(iEvent);
149     
150     // get the event summary data
151     char esdName[256]; 
152     sprintf(esdName, "ESD%d", iEvent);
153     AliESD* esd = (AliESD*) esdFile->Get(esdName);
154     if (!esd) {
155       Error("MUONmass_ESD", "no ESD object found for event %d", iEvent);
156       return kFALSE;
157     }
158
159     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; //
160
161     //    printf("\n Nb of events analysed: %d\r",iEvent);
162     //   cout << " number of tracks: " << nrectracks  <<endl;
163   
164     // loop over all reconstructed tracks (also first track of combination)
165     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
166
167       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
168
169       thetaX = muonTrack->GetThetaX();
170       thetaY = muonTrack->GetThetaY();
171
172       pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
173       fPzRec1  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
174       fPxRec1  = fPzRec1 * TMath::Tan(thetaX);
175       fPyRec1  = fPzRec1 * TMath::Tan(thetaY);
176       fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
177
178       fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
179       fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
180
181       ntrackhits = muonTrack->GetNHit();
182       fitfmin    = muonTrack->GetChi2();
183
184       // transverse momentum
185       Float_t pt1 = fV1.Pt();
186
187       // total momentum
188       Float_t p1 = fV1.P();
189
190       // Rapidity
191       Float_t rapMuon1 = fV1.Rapidity();
192
193       // chi2 per d.o.f.
194       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
195 //      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n", 
196 //           fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
197
198       // condition for good track (Chi2Cut and PtCut)
199
200       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
201
202         // fill histos hPtMuon and hChi2PerDof
203         hPtMuon->Fill(pt1);
204         hPMuon->Fill(p1);
205         hChi2PerDof->Fill(ch1);
206         hRapMuon->Fill(rapMuon1);
207         if (fCharge > 0) {
208           hPtMuonPlus->Fill(pt1);
209           hThetaPhiPlus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
210         } else {
211           hPtMuonMinus->Fill(pt1);
212           hThetaPhiMinus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415);
213         }
214         // loop over second track of combination
215         for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
216           
217           AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2);
218
219           thetaX = muonTrack->GetThetaX();
220           thetaY = muonTrack->GetThetaY();
221
222           pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
223           fPzRec2  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
224           fPxRec2  = fPzRec2 * TMath::Tan(thetaX);
225           fPyRec2  = fPzRec2 * TMath::Tan(thetaY);
226           fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
227
228           fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
229           fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
230
231           ntrackhits = muonTrack->GetNHit();
232           fitfmin    = muonTrack->GetChi2();
233
234           // transverse momentum
235           Float_t pt2 = fV2.Pt();
236
237           // chi2 per d.o.f.
238           Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
239
240           // condition for good track (Chi2Cut and PtCut)
241           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
242
243             // condition for opposite charges
244             if ((fCharge * fCharge2) == -1) {
245
246               // invariant mass
247               fVtot = fV1 + fV2;
248               Float_t invMass = fVtot.M();
249                     
250               // fill histos hInvMassAll and hInvMassRes
251               hInvMassAll->Fill(invMass);
252               hInvMassRes->Fill(invMass);
253               hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt());
254               if (invMass > massMin && invMass < massMax) {
255                 EventInMass++;
256                 hRapResonance->Fill(fVtot.Rapidity());
257                 hPtResonance->Fill(fVtot.Pt());
258               }
259
260             } // if (fCharge * fCharge2) == -1)
261           } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
262         } //  for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
263       } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
264     } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
265
266     hNumberOfTrack->Fill(nTracks);
267     esdFile->Delete();
268   } // for (Int_t iEvent = FirstEvent;
269
270 // Loop over events for bg event
271
272   Double_t thetaPlus,  phiPlus;
273   Double_t thetaMinus, phiMinus;
274   Float_t PtMinus, PtPlus;
275   
276   for (Int_t iEvent = 0; iEvent < hInvMassAll->Integral(); iEvent++) {
277
278     hThetaPhiPlus->GetRandom2(phiPlus, thetaPlus);
279     hThetaPhiMinus->GetRandom2(phiMinus,thetaMinus);
280     PtPlus = hPtMuonPlus->GetRandom();
281     PtMinus = hPtMuonMinus->GetRandom();
282
283     fPxRec1  = PtPlus * TMath::Cos(TMath::Pi()/180.*phiPlus);
284     fPyRec1  = PtPlus * TMath::Sin(TMath::Pi()/180.*phiPlus);
285     fPzRec1  = PtPlus / TMath::Tan(TMath::Pi()/180.*thetaPlus);
286
287     fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
288     fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
289
290     fPxRec2  = PtMinus * TMath::Cos(TMath::Pi()/180.*phiMinus);
291     fPyRec2  = PtMinus * TMath::Sin(TMath::Pi()/180.*phiMinus);
292     fPzRec2  = PtMinus / TMath::Tan(TMath::Pi()/180.*thetaMinus);
293
294     fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
295     fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
296
297     // invariant mass
298     fVtot = fV1 + fV2;
299       
300     // fill histos hInvMassAll and hInvMassRes
301     hInvMassBg->Fill(fVtot.M());
302     hInvMassBgk_vs_Pt->Fill( fVtot.M(), fVtot.Pt() );
303   }
304
305   histoFile->Write();
306   histoFile->Close();
307
308   cout << "MUONmassPlot " << endl;
309   cout << "FirstEvent " << FirstEvent << endl;
310   cout << "LastEvent " << LastEvent << endl;
311   cout << "ResType " << ResType << endl;
312   cout << "Chi2Cut " << Chi2Cut << endl;
313   cout << "PtCutMin " << PtCutMin << endl;
314   cout << "PtCutMax " << PtCutMax << endl;
315   cout << "massMin " << massMin << endl;
316   cout << "massMax " << massMax << endl;
317   cout << "EventInMass " << EventInMass << endl;
318
319   return kTRUE;
320 }
321