Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / MUONmassPlot.C
CommitLineData
c2a319d4 1// macro to make invariant mass plots
2// for combinations of 2 muons with opposite charges,
3// from root file "MUONtrackReco.root" containing the result of track reconstruction,
4// generated by the macro "MUONrecoNtuple.C".
5// Histograms are stored on the "MUONmassPlot.root" file.
6// A model for macros using the Ntuple in the file "MUONtrackReco.root"
7// may be found in "MUONtrackRecoModel.C":
8// it has been obtained by reading with Root a file "MUONtrackReco.root",
9// and executing the command:
10// MUONtrackReco->MakeCode("MUONtrackRecoModel.C")
11
12// Arguments:
13// FirstEvent (default 0)
14// LastEvent (default 0)
15// ResType (default 553)
16// 553 for Upsilon, anything else for J/Psi
c2a319d4 17// Chi2Cut (default 100)
18// to keep only tracks with chi2 per d.o.f. < Chi2Cut
8904b7e3 19// PtCutMin (default 1)
20// to keep only tracks with transverse momentum > PtCutMin
21// PtCutMax (default 10000)
22// to keep only tracks with transverse momentum < PtCutMax
23// massMin (default 9.17 for Upsilon)
24// & massMax (default 9.77 for Upsilon)
25// to calculate the reconstruction efficiency for resonances with invariant mass
26// massMin < mass < massMax.
27
28// compile MUONrecoNtuple.C and load shared library in this macro
29// Add parameters and histograms for analysis
30
31void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77)
c2a319d4 32{
8904b7e3 33 cout << "MUONmassPlot " << endl;
34 cout << "FirstEvent " << FirstEvent << endl;
35 cout << "LastEvent " << LastEvent << endl;
36 cout << "ResType " << ResType << endl;
37// cout << "Nsig " << Nsig << endl;
38 cout << "Chi2Cut " << Chi2Cut << endl;
39 cout << "PtCutMin " << PtCutMin << endl;
40 cout << "PtCutMax " << PtCutMax << endl;
41 cout << "massMin " << massMin << endl;
42 cout << "massMax " << massMax << endl;
43
c2a319d4 44
45 //////////////////////////////////////////////////////////
46 // This file has been automatically generated
47 // (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
48 // from TTree MUONtrackReco/MUONtrackReco
49 // found on file: MUONtrackReco.root
50 //////////////////////////////////////////////////////////
51
52
53 //Reset ROOT and connect tree file
54 gROOT->Reset();
8904b7e3 55// Dynamically link some shared libs
56
57 if (gClassTable->GetID("AliRun") < 0) {
58 gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER");
5b64e914 59// gROOT->LoadMacro("loadlibs.C");
60// loadlibs();
8904b7e3 61// compile MUONrecoNtuple and load shared library
62 gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf");
63 }
64
c2a319d4 65 TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
66 if (!f) {
67 f = new TFile("MUONtrackReco.root");
68 }
69 TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
eff57104 70 MUONtrackReco->SetMakeClass(1);
c2a319d4 71
72//Declaration of leaves types
73 Int_t fEvent;
74 UInt_t fUniqueID;
75 UInt_t fBits;
76 Int_t Tracks_;
8904b7e3 77 Int_t Tracks_fCharge[500];
78 Float_t Tracks_fPxRec[500];
79 Float_t Tracks_fPyRec[500];
80 Float_t Tracks_fPzRec[500];
81 Float_t Tracks_fZRec[500];
82 Float_t Tracks_fZRec1[500];
83 Int_t Tracks_fNHits[500];
84 Float_t Tracks_fChi2[500];
85 Float_t Tracks_fPxGen[500];
86 Float_t Tracks_fPyGen[500];
87 Float_t Tracks_fPzGen[500];
88 UInt_t Tracks_fUniqueID[500];
89 UInt_t Tracks_fBits[500];
c2a319d4 90
91 //Set branch addresses
92 //MUONtrackReco->SetBranchAddress("Header",&Header);
93 MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
94 MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
95 MUONtrackReco->SetBranchAddress("fBits",&fBits);
eff57104 96// MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
97 MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
c2a319d4 98 MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
99 MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
100 MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
101 MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
102 MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
103 MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
104 MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
105 MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
106 MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
107 MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
108 MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
109 MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
110 MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
111
112// This is the loop skeleton
113// To read only selected branches, Insert statements like:
114// MUONtrackReco->SetBranchStatus("*",0); // disable all branches
115// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
116
117 Int_t nentries = MUONtrackReco->GetEntries();
118
119 Int_t nbytes = 0;
120 // for (Int_t i=0; i<nentries;i++) {
121 // nbytes += MUONtrackReco->GetEntry(i);
122 // }
123
124 /////////////////////////////////////////////////////////////////
125 // Here comes the specialized part for MUONmassPlot
126 /////////////////////////////////////////////////////////////////
127
128 // File for histograms and histogram booking
129 TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
130 TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
8904b7e3 131 TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
c2a319d4 132 TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
8904b7e3 133 TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
134 if (ResType == 553) {
135 TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
136 } else {
137 TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
138 }
c2a319d4 139
8904b7e3 140 TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
5b64e914 141 TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5.,-2);
142 TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
8904b7e3 143 TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
144 Int_t EventInMass = 0;
145 Float_t muonMass = 0.105658389;
146 Float_t UpsilonMass = 9.46037;
147 Float_t JPsiMass = 3.097;
c2a319d4 148 // Loop over events
149 for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
150 // get current event
8904b7e3 151 Int_t NumberOfTrack = 0;
152 cout <<" event: " << event <<endl;
153 cout << " number of tracks: " << Tracks_ <<endl;
c2a319d4 154 nbytes += MUONtrackReco->GetEntry(event);
155 // loop over all reconstructed tracks (also first track of combination)
156 for (Int_t t1 = 0; t1 < Tracks_; t1++) {
157 // transverse momentum
158 Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
8904b7e3 159
160 // total momentum
161 Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] );
162 // Rapidity
163 Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass);
164
c2a319d4 165 // chi2 per d.o.f.
166 Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
8904b7e3 167 printf(" px %f py %f pz %f NHits %d Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1);
c2a319d4 168 // condition for good track (Chi2Cut and PtCut)
8904b7e3 169 if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
c2a319d4 170 // fill histos hPtMuon and hChi2PerDof
8904b7e3 171 NumberOfTrack++;
c2a319d4 172 hPtMuon->Fill(pt1);
8904b7e3 173 hPMuon->Fill(p1);
c2a319d4 174 hChi2PerDof->Fill(ch1);
8904b7e3 175 hRapMuon->Fill(rapMuon1);
c2a319d4 176 // loop over second track of combination
177 for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
178 // transverse momentum
179 Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
180 // chi2 per d.o.f.
181 Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
182 // condition for good track (Chi2Cut and PtCut)
8904b7e3 183 if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) {
c2a319d4 184 // condition for opposite charges
185 if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
186 // invariant mass
8904b7e3 187 Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass);
c2a319d4 188 // fill histos hInvMassAll and hInvMassRes
189 hInvMassAll->Fill(invMass);
190 hInvMassRes->Fill(invMass);
8904b7e3 191 Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2];
192 Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2];
193 Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2];
194 if (invMass > massMin && invMass < massMax) {
195 EventInMass++;
196 Float_t rapRes = 0;
197 if (ResType == 553) {
198 rapRes = rapParticle(pxRes,pyRes,pzRes,UpsilonMass);
199 } else {
200 rapRes = rapParticle(pxRes,pyRes,pzRes,JPsiMass);
201 }
202 hRapResonance->Fill(rapRes);
203 hPtResonance->Fill(TMath::Sqrt(pxRes*pxRes+pyRes*pyRes));
204 }
205
c2a319d4 206 } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
8904b7e3 207 } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
c2a319d4 208 } // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
8904b7e3 209 } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
c2a319d4 210 } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
8904b7e3 211 hNumberOfTrack->Fill(NumberOfTrack);
c2a319d4 212 } // for (Int_t event = FirstEvent;
213
214 histoFile->Write();
215 histoFile->Close();
8904b7e3 216
217 cout << "MUONmassPlot " << endl;
218 cout << "FirstEvent " << FirstEvent << endl;
219 cout << "LastEvent " << LastEvent << endl;
220 cout << "ResType " << ResType << endl;
221// cout << "Nsig " << Nsig << endl;
222 cout << "Chi2Cut " << Chi2Cut << endl;
223 cout << "PtCutMin " << PtCutMin << endl;
224 cout << "PtCutMax " << PtCutMax << endl;
225 cout << "massMin " << massMin << endl;
226 cout << "massMax " << massMax << endl;
227 cout << "EventInMass " << EventInMass << endl;
c2a319d4 228}
229
8904b7e3 230
231Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
c2a319d4 232{
8904b7e3 233 // return invariant mass for particle 1 & 2 whose masses are equal to muonMass
234
c2a319d4 235 Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
236 Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
237 return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
238}
8904b7e3 239
240Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
241{
242 // return rapidity for particle
243
244 // Particle energy
245 Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
246 // Rapidity
247 Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
248 return Rapidity;
249}