Formatting changes.
[u/mrichter/AliRoot.git] / MUON / MUONmassPlot_NewIO.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 "TParticle.h"
9 #include "TTree.h"
10
11 // STEER includes
12 #include "AliRun.h"
13 #include "AliRunLoader.h"
14 #include "AliHeader.h"
15 #include "AliLoader.h"
16 #include "AliStack.h"
17
18 // MUON includes
19 #include "AliMUON.h"
20 #include "AliMUONData.h"
21 #include "AliMUONHit.h"
22 #include "AliMUONConstants.h"
23 #include "AliMUONDigit.h"
24 #include "AliMUONRawCluster.h"
25 #include "AliMUONGlobalTrigger.h"
26 #include "AliMUONLocalTrigger.h"
27 #include "AliMUONTrack.h"
28 #include "AliMUONTrackParam.h"
29 #include "AliESDMuonTrack.h"
30 #endif
31 //
32 // Macro MUONmassPlot.C for new I/O
33 // Ch. Finck, Subatech, Jan. 2004
34 //
35
36 // macro to make invariant mass plots
37 // for combinations of 2 muons with opposite charges,
38 // from root file "MUON.tracks.root" containing the result of track reconstruction.
39 // Histograms are stored on the "MUONmassPlot.root" file.
40 // introducing TLorentzVector for parameter calculations (Pt, P,rap,etc...)
41 // using Invariant Mass for rapidity.
42
43 // Arguments:
44 //   FirstEvent (default 0)
45 //   LastEvent (default 0)
46 //   ResType (default 553)
47 //      553 for Upsilon, anything else for J/Psi
48 //   Chi2Cut (default 100)
49 //      to keep only tracks with chi2 per d.o.f. < Chi2Cut
50 //   PtCutMin (default 1)
51 //      to keep only tracks with transverse momentum > PtCutMin
52 //   PtCutMax (default 10000)
53 //      to keep only tracks with transverse momentum < PtCutMax
54 //   massMin (default 9.17 for Upsilon) 
55 //      &  massMax (default 9.77 for Upsilon) 
56 //         to calculate the reconstruction efficiency for resonances with invariant mass
57 //         massMin < mass < massMax.
58
59 // Add parameters and histograms for analysis 
60
61 void MUONmassPlot(char* filename="galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, 
62                   Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
63                   Float_t massMin = 9.17,Float_t massMax = 9.77)
64 {
65   cout << "MUONmassPlot " << endl;
66   cout << "FirstEvent " << FirstEvent << endl;
67   cout << "LastEvent " << LastEvent << endl;
68   cout << "ResType " << ResType << endl;
69   cout << "Chi2Cut " << Chi2Cut << endl;
70   cout << "PtCutMin " << PtCutMin << endl;
71   cout << "PtCutMax " << PtCutMax << endl;
72   cout << "massMin " << massMin << endl;
73   cout << "massMax " << massMax << endl;
74
75
76   //Reset ROOT and connect tree file
77   gROOT->Reset();
78
79
80   // File for histograms and histogram booking
81   TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
82   TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
83   TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
84   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
85   TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
86   TH1F *hInvMassRes;
87
88   if (ResType == 553) {
89     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
90   } else {
91     hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
92   }
93
94   TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
95   TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
96   TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
97   TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
98
99
100   // settings
101   Int_t EventInMass = 0;
102   Float_t muonMass = 0.105658389;
103 //   Float_t UpsilonMass = 9.46037;
104 //   Float_t JPsiMass = 3.097;
105
106   Double_t bendingSlope, nonBendingSlope, pYZ;
107   Double_t fPxRec1, fPyRec1, fPzRec1, fZRec1, fE1;
108   Double_t fPxRec2, fPyRec2, fPzRec2, fZRec2, fE2;
109   Int_t fCharge, fCharge2;
110
111   Int_t ntrackhits, nevents;
112   Double_t fitfmin;
113
114   TClonesArray * recTracksArray;
115   TLorentzVector fV1, fV2, fVtot;
116   
117   // Creating Run Loader and openning file containing Hits
118   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
119   if (RunLoader == 0x0) {
120     printf(">>> Error : Error Opening %s file \n",filename);
121     return;
122   }
123   
124   AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
125   MUONLoader->LoadTracks("READ");
126
127   // Creating MUON data container
128   AliMUONData muondata(MUONLoader,"MUON","MUON");
129   
130   nevents = RunLoader->GetNumberOfEvents();
131   
132   AliMUONTrack * rectrack;
133   AliMUONTrackParam *trackParam;
134         
135   // Loop over events
136   for (Int_t ievent = FirstEvent; ievent <= TMath::Min(LastEvent, nevents - 1); ievent++) {
137
138     // get current event
139     RunLoader->GetEvent(ievent);
140     
141     muondata.SetTreeAddress("RT");
142     muondata.GetRecTracks();
143     recTracksArray = muondata.RecTracks();
144
145     Int_t nrectracks = (Int_t) recTracksArray->GetEntriesFast(); //
146
147     printf("\n Nb of events analysed: %d\r",ievent);
148     //   cout << " number of tracks: " << nrectracks  <<endl;
149   
150     // loop over all reconstructed tracks (also first track of combination)
151      for (Int_t irectracks = 0; irectracks <  nrectracks;  irectracks++) {
152
153       rectrack = (AliMUONTrack*) recTracksArray->At(irectracks);
154
155       trackParam = rectrack->GetTrackParamAtVertex();
156       bendingSlope   = trackParam->GetBendingSlope();
157       nonBendingSlope = trackParam->GetNonBendingSlope();
158
159       pYZ     = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
160       fPzRec1  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
161       fPxRec1  = fPzRec1 * nonBendingSlope;
162       fPyRec1  = fPzRec1 * bendingSlope;
163       fZRec1   = trackParam->GetZ();
164       fCharge = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
165
166       fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1);
167       fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1);
168
169       ntrackhits = rectrack->GetNTrackHits();
170       fitfmin = rectrack->GetFitFMin();
171
172       // transverse momentum
173       Float_t pt1 = fV1.Pt();
174
175       // total momentum
176       Float_t p1 = fV1.P();
177
178       // Rapidity
179       Float_t rapMuon1 = fV1.Rapidity();
180
181       // chi2 per d.o.f.
182       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
183 //      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f charge %d\n", 
184 //           fPxRec1, fPyRec1, fPzRec1, ntrackhits, ch1, fCharge);
185
186       // condition for good track (Chi2Cut and PtCut)
187
188       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
189
190         // fill histos hPtMuon and hChi2PerDof
191         hPtMuon->Fill(pt1);
192         hPMuon->Fill(p1);
193         hChi2PerDof->Fill(ch1);
194         hRapMuon->Fill(rapMuon1);
195
196         // loop over second track of combination
197         for (Int_t irectracks2 = irectracks + 1; irectracks2 < nrectracks; irectracks2++) {
198
199           rectrack = (AliMUONTrack*) recTracksArray->At(irectracks2);
200          
201           trackParam = rectrack->GetTrackParamAtVertex();
202           bendingSlope    = trackParam->GetBendingSlope();
203           nonBendingSlope = trackParam->GetNonBendingSlope();
204           
205           pYZ      = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
206           fPzRec2  = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
207           fPxRec2  = fPzRec2 * nonBendingSlope;
208           fPyRec2  = fPzRec2 * bendingSlope;
209           fZRec2   = trackParam->GetZ();
210           fCharge2 = Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
211
212           fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
213           fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
214
215           ntrackhits = rectrack->GetNTrackHits();
216           fitfmin = rectrack->GetFitFMin();
217
218           // transverse momentum
219           Float_t pt2 = fV2.Pt();
220
221           // chi2 per d.o.f.
222           Float_t ch2 = fitfmin  / (2.0 * ntrackhits - 5);
223
224           // condition for good track (Chi2Cut and PtCut)
225           if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
226
227             // condition for opposite charges
228             if ((fCharge * fCharge2) == -1) {
229
230               // invariant mass
231               fVtot = fV1 + fV2;
232               Float_t invMass = fVtot.M();
233                     
234               // fill histos hInvMassAll and hInvMassRes
235               hInvMassAll->Fill(invMass);
236               hInvMassRes->Fill(invMass);
237
238               if (invMass > massMin && invMass < massMax) {
239                 EventInMass++;
240                 hRapResonance->Fill(fVtot.Rapidity());
241                 hPtResonance->Fill(fVtot.Pt());
242               }
243
244             } // if (fCharge * fCharge2) == -1)
245           } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
246         } //  for (Int_t irectracks2 = irectracks + 1; irectracks2 < irectracks; irectracks2++)
247       } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
248     } // for (Int_t irectracks = 0; irectracks < nrectracks; irectracks++)
249
250     hNumberOfTrack->Fill(nrectracks);
251   } // for (Int_t ievent = FirstEvent;
252
253   histoFile->Write();
254   histoFile->Close();
255
256   cout << "MUONmassPlot " << endl;
257   cout << "FirstEvent " << FirstEvent << endl;
258   cout << "LastEvent " << LastEvent << endl;
259   cout << "ResType " << ResType << endl;
260   cout << "Chi2Cut " << Chi2Cut << endl;
261   cout << "PtCutMin " << PtCutMin << endl;
262   cout << "PtCutMax " << PtCutMax << endl;
263   cout << "massMin " << massMin << endl;
264   cout << "massMax " << massMax << endl;
265   cout << "EventInMass " << EventInMass << endl;
266 }
267