]>
Commit | Line | Data |
---|---|---|
d949b15e | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | // Macro (upgraded version of MUONmassPlot_ESD.C) to make : | |
19 | // 1) Ntuple (Ktuple) containing Upsilon kinematics variables (from kinematics.root files) | |
20 | // 2) Ntuple (ESDtuple) containing Upsilon kinematics variables from reconstruction and | |
21 | // combinations of 2 muons with opposite charges, | |
22 | // 3) Some QA histograms | |
23 | // Ntuple are stored in the file MUONefficiency.root and ESD tree and QA histograms in AliESDs.root | |
24 | ||
25 | // Arguments: | |
26 | // FirstEvent (default 0) | |
27 | // LastEvent (default 0) | |
28 | // ResType (default 553) | |
29 | // 553 for Upsilon, anything else for J/Psi | |
30 | // Chi2Cut (default 100) | |
31 | // to keep only tracks with chi2 per d.o.f. < Chi2Cut | |
32 | // PtCutMin (default 1) | |
33 | // to keep only tracks with transverse momentum > PtCutMin | |
34 | // PtCutMax (default 10000) | |
35 | // to keep only tracks with transverse momentum < PtCutMax | |
36 | // massMin (default 9.17 for Upsilon) | |
37 | // & massMax (default 9.77 for Upsilon) | |
38 | // to calculate the reconstruction efficiency for resonances with invariant mass | |
39 | // massMin < mass < massMax. | |
40 | ||
41 | // Add parameters and histograms for analysis | |
42 | ||
43 | // Christophe Suire, IPN Orsay | |
44 | ||
45 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
46 | // ROOT includes | |
47 | #include "TTree.h" | |
48 | #include "TBranch.h" | |
49 | #include "TClonesArray.h" | |
50 | #include "TLorentzVector.h" | |
51 | #include "TFile.h" | |
52 | #include "TH1.h" | |
53 | #include "TH2.h" | |
54 | #include "TParticle.h" | |
55 | #include "TTree.h" | |
56 | #include "TString.h" | |
57 | #include <Riostream.h> | |
58 | ||
59 | // STEER includes | |
60 | #include "AliRun.h" | |
61 | #include "AliRunLoader.h" | |
62 | #include "AliHeader.h" | |
63 | #include "AliLoader.h" | |
64 | #include "AliStack.h" | |
65 | #include "AliMagF.h" | |
66 | #include "AliESD.h" | |
67 | ||
68 | // MUON includes | |
69 | #include "AliESDMuonTrack.h" | |
70 | #endif | |
71 | ||
72 | ||
73 | Bool_t MUONefficiency(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 11000000, | |
74 | char* esdFileName = "AliESDs.root", Int_t ResType = 553, | |
75 | Float_t Chi2Cut = 100., Float_t PtCutMin = 0., Float_t PtCutMax = 10000., | |
76 | Float_t massMin = 9.17,Float_t massMax = 9.77) | |
77 | { // MUONefficiency starts | |
78 | cout << "MUONmassPlot " << endl; | |
79 | cout << "FirstEvent " << FirstEvent << endl; | |
80 | cout << "LastEvent " << LastEvent << endl; | |
81 | cout << "ResType " << ResType << endl; | |
82 | cout << "Chi2Cut " << Chi2Cut << endl; | |
83 | cout << "PtCutMin " << PtCutMin << endl; | |
84 | cout << "PtCutMax " << PtCutMax << endl; | |
85 | cout << "massMin " << massMin << endl; | |
86 | cout << "massMax " << massMax << endl; | |
87 | ||
88 | ||
89 | //Reset ROOT and connect tree file | |
90 | gROOT->Reset(); | |
91 | ||
92 | // Printing Level | |
93 | Int_t PRINTLEVEL = 0 ; | |
94 | Int_t SELECT = 0 ; // not used | |
95 | ||
96 | //for kinematic, i.e. reference tracks | |
97 | TNtuple *Ktuple = new TNtuple("Ktuple","Kinematics NTuple","ev:npart:id:idmo:idgdmo:p:pt:y:theta:pseudorap:vx:vy:vz"); | |
98 | ||
99 | //for reconstruction | |
100 | TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.); | |
101 | TH1F *hPtMuonPlus = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.); | |
102 | TH1F *hPtMuonMinus = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.); | |
103 | TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.); | |
104 | TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.); | |
105 | TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.); | |
106 | TH1F *hInvMassBg = new TH1F("hInvMassBg", "Mu+Mu- invariant mass BG(GeV/c2)", 480, 0., 12.); | |
107 | TH2F *hInvMassAll_vs_Pt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.); | |
108 | TH2F *hInvMassBgk_vs_Pt = new TH2F("hInvMassBgk_vs_Pt","hInvMassBgk_vs_Pt",480,0.,12.,80,0.,20.); | |
109 | TH1F *hInvMassRes; | |
110 | if (ResType == 553) { | |
111 | hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.); | |
112 | } else { | |
113 | hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.); | |
114 | } | |
115 | ||
116 | TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5); | |
117 | TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2); | |
118 | TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2); | |
119 | TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.); | |
120 | TH2F *hThetaPhiPlus = new TH2F("hThetaPhiPlus", "Theta vs Phi +", 760, -190., 190., 400, 160., 180.); | |
121 | TH2F *hThetaPhiMinus = new TH2F("hThetaPhiMinus", "Theta vs Phi -", 760, -190., 190., 400, 160., 180.); | |
122 | ||
123 | TNtuple *ESDtuple = new TNtuple("ESDtuple","Reconstructed Mu+Mu- pairs and Upsilon","ev:tw:pt:y:theta:minv:pt1:y1:theta1:q1:trig1:pt2:y2:theta2:q2:trig2"); | |
124 | TNtuple *ESDtupleBck = new TNtuple("ESDtupleBck","Reconstructed Mu+Mu- pairs for Background","ev:pt:y:theta:minv:pt1:y1:theta1:pt2:y2:theta2"); | |
125 | ||
126 | ||
127 | // settings | |
128 | Int_t EventInMass = 0; | |
129 | Float_t muonMass = 0.105658389; | |
130 | Float_t UpsilonMass = 9.46037; | |
131 | Float_t JPsiMass = 3.097; | |
132 | ||
133 | Double_t thetaX, thetaY, pYZ; | |
134 | Double_t fPxRec1, fPyRec1, fPzRec1, fE1; | |
135 | Double_t fPxRec2, fPyRec2, fPzRec2, fE2; | |
136 | Int_t fCharge1, fCharge2; | |
137 | ||
138 | Int_t ntrackhits, nevents; | |
139 | Double_t fitfmin; | |
140 | ||
141 | TLorentzVector fV1, fV2, fVtot; | |
142 | ||
143 | // set off mag field | |
144 | AliMagF::SetReadField(kFALSE); | |
145 | ||
146 | // open run loader and load gAlice, kinematics and header | |
147 | AliRunLoader* runLoader = AliRunLoader::Open(filename); | |
148 | if (!runLoader) { | |
149 | Error("MUONefficiency", "getting run loader from file %s failed", filename); | |
150 | return kFALSE; | |
151 | } | |
152 | ||
153 | runLoader->LoadgAlice(); | |
154 | gAlice = runLoader->GetAliRun(); | |
155 | if (!gAlice) { | |
156 | Error("MUONefficiency", "no galice object found"); | |
157 | return kFALSE; | |
158 | } | |
159 | ||
160 | // open the ESD file | |
161 | TFile* esdFile = TFile::Open(esdFileName); | |
162 | if (!esdFile || !esdFile->IsOpen()) { | |
163 | Error("MUONefficiency", "opening ESD file %s failed", esdFileName); | |
164 | return kFALSE; | |
165 | } | |
166 | ||
167 | AliESD* esd = new AliESD(); | |
168 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
169 | if (!tree) { | |
170 | Error("CheckESD", "no ESD tree found"); | |
171 | return kFALSE; | |
172 | } | |
173 | tree->SetBranchAddress("ESD", &esd); | |
174 | ||
175 | runLoader->LoadHeader(); | |
176 | nevents = runLoader->GetNumberOfEvents(); | |
177 | ||
178 | // to access the particle Stack | |
179 | runLoader->LoadKinematics("READ"); | |
180 | ||
181 | TParticle *particle; | |
182 | Int_t track1Id = 0 ; | |
183 | Int_t track1PDGId = 0 ; | |
184 | Int_t track1MotherId = 0 ; | |
185 | Int_t track1MotherPDGId = 0 ; | |
186 | Int_t track1Trigger = 0 ; | |
187 | Float_t track1TriggerChi2 = 0 ; | |
188 | Int_t track2Id = 0 ; | |
189 | Int_t track2PDGId = 0 ; | |
190 | Int_t track2MotherId = 0 ; | |
191 | Int_t track2MotherPDGId = 0 ; | |
192 | Int_t track2Trigger = 0 ; | |
193 | Float_t track2TriggerChi2 = 0 ; | |
194 | ||
195 | ||
196 | // Loop over events | |
197 | for (Int_t iEvent = FirstEvent; iEvent <= TMath::Min(LastEvent, nevents - 1); iEvent++) { // Start event loop | |
198 | ||
199 | if (iEvent%1000 == 0 ) | |
200 | printf("\n Nb of events analysed: %d \n",iEvent); | |
201 | ||
202 | // get current event | |
203 | runLoader->GetEvent(iEvent); | |
204 | ||
205 | // get the stack and fill the kine tree | |
206 | AliStack *theStack = runLoader->Stack(); | |
207 | if (PRINTLEVEL > 0) theStack->DumpPStack (); | |
208 | ||
209 | Int_t nparticles = (Int_t)runLoader->TreeK()->GetEntries(); | |
210 | Int_t nprimarypart = theStack->GetNprimary(); | |
211 | Int_t ntracks = theStack->GetNtrack(); | |
212 | ||
213 | if (PRINTLEVEL || (iEvent%100==0)) printf("\n >>> Event %d \n",iEvent); | |
214 | if (PRINTLEVEL) cout << nprimarypart << " Particles generated (total is " << ntracks << ")"<< endl ; | |
215 | ||
216 | ||
217 | ||
218 | for(Int_t iparticle=0; iparticle<nparticles; iparticle++) { // Start loop over particles | |
219 | particle = theStack->Particle(iparticle); | |
220 | ||
221 | Int_t muId = particle->GetPdgCode(); | |
222 | Int_t muM = particle->GetFirstMother(); | |
223 | Int_t muGM = 0; | |
224 | Float_t muP = particle->P(); | |
225 | Float_t muPt = TMath::Sqrt(particle->Px()*particle->Px()+particle->Py()*particle->Py()); | |
226 | Float_t muY = 0.5*TMath::Log((particle->Energy()+particle->Pz()+1.e-13)/(particle->Energy()-particle->Pz()+1.e-13)); | |
227 | if (muM >= 0) { | |
228 | //cout << "in stack " << partM << endl ; | |
229 | TParticle *theMum = theStack->Particle(muM); | |
230 | muM = theMum->GetPdgCode(); | |
231 | //cout << "the Mum " << partM << endl ; | |
232 | ||
233 | muGM = theMum->GetFirstMother() ; | |
234 | if (muGM >= 0){ | |
235 | TParticle *grandMa = theStack->Particle(muGM); | |
236 | muGM = grandMa->GetPdgCode(); | |
237 | } | |
238 | else muGM=0; | |
239 | } | |
240 | else muM=0; | |
241 | ||
242 | Float_t muT = particle->Theta()*180/TMath::Pi(); | |
243 | Float_t muE = particle->Eta(); | |
244 | ||
245 | Float_t muVx = particle->Vx(); | |
246 | Float_t muVy = particle->Vy(); | |
247 | Float_t muVz = particle->Vz(); | |
248 | ||
249 | // If a write error occurs, the number of bytes returned is -1. | |
250 | // If no data are written, because e.g. the branch is disabled, | |
251 | // the number of bytes returned is 0. | |
252 | Int_t errCode = Ktuple->Fill(iEvent,nparticles,muId,muM,muGM,muP,muPt,muY,muT,muE,muVx,muVy,muVz); | |
253 | if (PRINTLEVEL || errCode < 1) printf("iEvent %d,nparticles %d,muId %d,muM %d,muGM %d,muP %.2f,muPt %.2f,muY %.2f,muT %.2f,muE %.2f,muVx %.2f,muVy %.2f,muVz %.2f \n", iEvent,nparticles,muId,muM,muGM,muP,muPt,muY,muT,muE,muVx,muVy,muVz); | |
254 | ||
255 | } // End loop over particles | |
256 | ||
257 | ||
258 | ||
259 | // get the event summary data | |
260 | tree->GetEvent(iEvent); | |
261 | if (!esd) { | |
262 | Error("CheckESD", "no ESD object found for event %d", iEvent); | |
263 | return kFALSE; | |
264 | } | |
265 | ||
266 | Int_t triggerWord = esd->GetTrigger(); | |
267 | Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; | |
268 | ||
269 | if (PRINTLEVEL > 0){ | |
270 | printf("\n Nb of events analysed: %d \n",iEvent); | |
271 | cout << " number of tracks: " << nTracks <<endl; | |
272 | } | |
273 | ||
274 | // loop over all reconstructed tracks (also first track of combination) | |
275 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
276 | ||
277 | AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); | |
278 | ||
279 | //if (PRINTLEVEL > 5) cout << "1st muonTrack->GetTrackID() : " << track1Id << endl; | |
280 | ||
281 | if(SELECT && track1Id) { | |
282 | particle = theStack->Particle(track1Id); | |
283 | track1PDGId = particle->GetPdgCode() ; | |
284 | track1MotherId = particle->GetFirstMother(); | |
285 | if (track1MotherId >=0 ) | |
286 | track1MotherPDGId = ((TParticle*) theStack->Particle(track1MotherId))->GetPdgCode(); | |
287 | if (PRINTLEVEL > 0) cout << "track1MotherPDGId = " << track1MotherPDGId << endl ; | |
288 | } | |
289 | ||
290 | // Trigger | |
291 | if (PRINTLEVEL > 5) cout << "MatchTrigger " << muonTrack->GetMatchTrigger() << " and Chi2 of matching tracks " << track1TriggerChi2 << endl ; | |
292 | track1Trigger = muonTrack->GetMatchTrigger(); | |
293 | if (track1Trigger) | |
294 | track1TriggerChi2 = muonTrack->GetChi2MatchTrigger(); | |
295 | else | |
296 | track1TriggerChi2 = 0. ; | |
297 | ||
298 | thetaX = muonTrack->GetThetaX(); | |
299 | thetaY = muonTrack->GetThetaY(); | |
300 | ||
301 | pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum()); | |
302 | fPzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY)); | |
303 | fPxRec1 = fPzRec1 * TMath::Tan(thetaX); | |
304 | fPyRec1 = fPzRec1 * TMath::Tan(thetaY); | |
305 | fCharge1 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); | |
306 | ||
307 | fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1); | |
308 | fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1); | |
309 | ||
310 | ntrackhits = muonTrack->GetNHit(); | |
311 | fitfmin = muonTrack->GetChi2(); | |
312 | ||
313 | // transverse momentum | |
314 | Float_t pt1 = fV1.Pt(); | |
315 | ||
316 | // total momentum | |
317 | Float_t p1 = fV1.P(); | |
318 | ||
319 | // Rapidity | |
320 | Float_t rapMuon1 = fV1.Rapidity(); | |
321 | ||
322 | // chi2 per d.o.f. | |
323 | ||
324 | Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5); | |
325 | if (PRINTLEVEL > 5 ) printf(" px %f py %f pz %f pt %f NHits %d Norm.chi2 %f charge %d\n",fPxRec1, fPyRec1, fPzRec1, pt1, ntrackhits, ch1, fCharge1); | |
326 | ||
327 | ||
328 | if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) { // condition for good track (Chi2Cut and PtCut) | |
329 | if (PRINTLEVEL > 8) cout << "inside pt and chi2 cuts " << endl ; | |
330 | ||
331 | // fill histos hPtMuon and hChi2PerDof | |
332 | hPtMuon->Fill(pt1); | |
333 | hPMuon->Fill(p1); | |
334 | hChi2PerDof->Fill(ch1); | |
335 | hRapMuon->Fill(rapMuon1); | |
336 | ||
337 | if (fCharge1 > 0) { | |
338 | hPtMuonPlus->Fill(pt1); | |
339 | hThetaPhiPlus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415); | |
340 | } else { | |
341 | hPtMuonMinus->Fill(pt1); | |
342 | hThetaPhiMinus->Fill(TMath::ATan2(fPyRec1,fPxRec1)*180./TMath::Pi(),TMath::ATan2(pt1,fPzRec1)*180./3.1415); | |
343 | } | |
344 | ||
345 | // loop over second track of combination | |
346 | for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) { | |
347 | ||
348 | AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2); | |
349 | //Int_t track2Id = muonTrack->GetTrackID(); | |
350 | //if (PRINTLEVEL > 5) cout << "2nd muonTrack->GetTrackID() : " << track2Id << endl; | |
351 | ||
352 | if(SELECT && track2Id) { | |
353 | particle = theStack->Particle(track2Id); | |
354 | track2PDGId = particle->GetPdgCode(); | |
355 | track2MotherId = particle->GetFirstMother(); | |
356 | if (track2MotherId >=0 ) | |
357 | track2MotherPDGId = ((TParticle*) theStack->Particle(track2MotherId))->GetPdgCode(); | |
358 | } | |
359 | ||
360 | track2Trigger = muonTrack->GetMatchTrigger(); | |
361 | if (track2Trigger) | |
362 | track2TriggerChi2 = muonTrack->GetChi2MatchTrigger(); | |
363 | else | |
364 | track2TriggerChi2 = 0. ; | |
365 | ||
366 | thetaX = muonTrack->GetThetaX(); | |
367 | thetaY = muonTrack->GetThetaY(); | |
368 | ||
369 | pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum()); | |
370 | fPzRec2 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY)); | |
371 | fPxRec2 = fPzRec2 * TMath::Tan(thetaX); | |
372 | fPyRec2 = fPzRec2 * TMath::Tan(thetaY); | |
373 | fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum())); | |
374 | ||
375 | fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2); | |
376 | fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2); | |
377 | ||
378 | ntrackhits = muonTrack->GetNHit(); | |
379 | fitfmin = muonTrack->GetChi2(); | |
380 | ||
381 | // transverse momentum | |
382 | Float_t pt2 = fV2.Pt(); | |
383 | ||
384 | // chi2 per d.o.f. | |
385 | Float_t ch2 = fitfmin / (2.0 * ntrackhits - 5); | |
386 | ||
387 | if (PRINTLEVEL > 5) cout << "track1MotherId : "<< track1MotherId << " track2MotherId : " << track2MotherId << endl ; | |
388 | if (PRINTLEVEL > 5) cout << "track1MotherPDGId : " << track1MotherPDGId << " track2MotherPDGId : " << track2MotherPDGId << endl ; | |
389 | ||
390 | // Select Condition | |
391 | if (!SELECT || (track2MotherId == track1MotherId && track2MotherPDGId == ResType && TMath::Abs(track1PDGId)==13 && TMath::Abs(track2PDGId)==13 )) { | |
392 | ||
393 | // condition for good track (Chi2Cut and PtCut) | |
394 | if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) { | |
395 | ||
396 | // condition for opposite charges | |
397 | if ((fCharge1 * fCharge2) == -1) { | |
398 | ||
399 | if (PRINTLEVEL > 8) cout << "---------> Now filling the Ntuple " << endl ; | |
400 | ||
401 | // invariant mass | |
402 | fVtot = fV1 + fV2; | |
403 | Float_t invMass = fVtot.M(); | |
404 | ||
405 | if (fCharge1 < 0){ //mu_minus is index 1 in the ntuple | |
406 | Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2}; | |
407 | ESDtuple->Fill(ESDFill); | |
408 | } | |
409 | else{ | |
410 | Float_t ESDFill[16] = {iEvent,triggerWord,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,invMass,fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fCharge2,track2TriggerChi2,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180,fCharge1,track1TriggerChi2}; | |
411 | ESDtuple->Fill(ESDFill); | |
412 | } | |
413 | ||
414 | // fill histos hInvMassAll and hInvMassRes | |
415 | hInvMassAll->Fill(invMass); | |
416 | hInvMassRes->Fill(invMass); | |
417 | hInvMassAll_vs_Pt->Fill(invMass,fVtot.Pt()); | |
418 | if (invMass > massMin && invMass < massMax) { | |
419 | EventInMass++; | |
420 | hRapResonance->Fill(fVtot.Rapidity()); | |
421 | hPtResonance->Fill(fVtot.Pt()); | |
422 | } | |
423 | ||
424 | } // if (fCharge1 * fCharge2) == -1) | |
425 | } // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) | |
426 | } // if (track2MotherId == track1MotherId && track2MotherPDGId == ResType) | |
427 | } // for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++) | |
428 | } // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) ) | |
429 | } // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++) | |
430 | ||
431 | hNumberOfTrack->Fill(nTracks); | |
432 | // esdFile->Delete(); | |
433 | ||
434 | } // End of event loop | |
435 | ||
436 | ||
437 | // Loop over events for bg event | |
438 | ||
439 | Double_t thetaPlus, phiPlus; | |
440 | Double_t thetaMinus, phiMinus; | |
441 | Float_t PtMinus, PtPlus; | |
442 | ||
443 | for (Int_t iEvent = 0; iEvent < hInvMassAll->Integral(); iEvent++) { // Loop over events for bg event | |
444 | // according to Christian a 3d histo phi-theta-pt would take better care | |
445 | // of all correlations | |
446 | ||
447 | hThetaPhiPlus->GetRandom2(phiPlus, thetaPlus); | |
448 | hThetaPhiMinus->GetRandom2(phiMinus,thetaMinus); | |
449 | PtPlus = hPtMuonPlus->GetRandom(); | |
450 | PtMinus = hPtMuonMinus->GetRandom(); | |
451 | ||
452 | fPxRec1 = PtPlus * TMath::Cos(TMath::Pi()/180.*phiPlus); | |
453 | fPyRec1 = PtPlus * TMath::Sin(TMath::Pi()/180.*phiPlus); | |
454 | fPzRec1 = PtPlus / TMath::Tan(TMath::Pi()/180.*thetaPlus); | |
455 | ||
456 | fE1 = TMath::Sqrt(muonMass * muonMass + fPxRec1 * fPxRec1 + fPyRec1 * fPyRec1 + fPzRec1 * fPzRec1); | |
457 | fV1.SetPxPyPzE(fPxRec1, fPyRec1, fPzRec1, fE1); | |
458 | ||
459 | fPxRec2 = PtMinus * TMath::Cos(TMath::Pi()/180.*phiMinus); | |
460 | fPyRec2 = PtMinus * TMath::Sin(TMath::Pi()/180.*phiMinus); | |
461 | fPzRec2 = PtMinus / TMath::Tan(TMath::Pi()/180.*thetaMinus); | |
462 | ||
463 | fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2); | |
464 | fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2); | |
465 | ||
466 | // invariant mass | |
467 | fVtot = fV1 + fV2; | |
468 | ||
469 | // fill histos hInvMassAll and hInvMassRes | |
470 | hInvMassBg->Fill(fVtot.M()); | |
471 | hInvMassBgk_vs_Pt->Fill( fVtot.M(), fVtot.Pt() ); | |
472 | ||
473 | // Ntuple for background... more convenient | |
474 | ESDtupleBck->Fill(iEvent,fVtot.Pt(),fVtot.Rapidity(),fVtot.Theta()/TMath::Pi()*180,fVtot.M(),fV2.Pt(),fV2.Rapidity(),fV2.Theta()/TMath::Pi()*180,fV1.Pt(),fV1.Rapidity(),fV1.Theta()/TMath::Pi()*180); | |
475 | ||
476 | } // End loop over events for background | |
477 | ||
478 | ||
479 | // File for histograms and histogram booking | |
480 | TString outfilename = "MUONefficiency.root"; | |
481 | TFile *histoFile = new TFile(outfilename.Data(), "RECREATE"); | |
482 | ||
483 | Ktuple->Write(); | |
484 | ESDtuple->Write(); | |
485 | //histoFile->Write(); | |
486 | ||
487 | histoFile->Close(); | |
488 | ||
489 | cout << "MUONefficiency " << endl; | |
490 | cout << "FirstEvent " << FirstEvent << endl; | |
491 | cout << "LastEvent " << LastEvent << endl; | |
492 | cout << "ResType " << ResType << endl; | |
493 | cout << "Chi2Cut " << Chi2Cut << endl; | |
494 | cout << "PtCutMin " << PtCutMin << endl; | |
495 | cout << "PtCutMax " << PtCutMax << endl; | |
496 | cout << "massMin " << massMin << endl; | |
497 | cout << "massMax " << massMax << endl; | |
498 | cout << "EventInMass " << EventInMass << endl; | |
499 | ||
500 | return kTRUE; | |
501 | } | |
502 |