]>
Commit | Line | Data |
---|---|---|
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 | ||
31 | void 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"); | |
59 | gROOT->LoadMacro("loadlibs.C"); | |
60 | loadlibs(); | |
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); |
141 | TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,2.,4.5); | |
142 | TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,2.,4.5); | |
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 | |
231 | Float_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 | |
240 | Float_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 | } |