]>
Commit | Line | Data |
---|---|---|
30b3fe0f | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | // ROOT includes | |
3 | #include <TTree.h> | |
4 | #include <TFile.h> | |
5 | #include <TH1.h> | |
6 | #include <TCanvas.h> | |
7 | #include <Riostream.h> | |
8 | #include <TROOT.h> | |
9 | #include <TClonesArray.h> | |
10 | #include <TGeoGlobalMagField.h> | |
11 | #include <TLorentzVector.h> | |
12 | ||
13 | // STEER includes | |
14 | #include "AliLog.h" | |
15 | #include "AliMagF.h" | |
16 | #include "AliESDEvent.h" | |
17 | #include "AliESDMuonTrack.h" | |
18 | #include "AliESDMuonCluster.h" | |
19 | #include "AliCDBPath.h" | |
20 | #include "AliCDBEntry.h" | |
21 | #include "AliCDBManager.h" | |
22 | ||
23 | // MUON includes | |
24 | #include "AliMUONTrack.h" | |
25 | #include "AliMUONVTrackStore.h" | |
26 | #include "AliMUONTrackParam.h" | |
27 | #include "AliMUONTrackExtrap.h" | |
28 | #include "AliMUONESDInterface.h" | |
29 | #include "AliMUONRecoCheck.h" | |
30 | #include "AliMUONVCluster.h" | |
31 | #include "AliMUONRecoParam.h" | |
32 | #endif | |
33 | ||
34 | /// \ingroup macros | |
35 | /// \file DIMUONFakes.C | |
36 | /// | |
37 | /// \author Ph. Pillot, Subatech, March. 2009 | |
38 | /// | |
39 | /// Macro to study the effects of fake tracks on the dimuon spectra | |
40 | /// Results are saved in the root file DiFakes.root | |
41 | /// Results are relevent provided that you use the same recoParams as for the reconstruction | |
42 | ||
43 | void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut); | |
44 | TTree* GetESDTree(TFile *esdFile); | |
45 | Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut); | |
46 | AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore, | |
47 | Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut); | |
48 | ||
49 | //----------------------------------------------------------------------- | |
50 | void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1, | |
51 | const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/") | |
52 | { | |
53 | ||
54 | //Reset ROOT and connect tree file | |
55 | gROOT->Reset(); | |
56 | ||
57 | // File for histograms and histogram booking | |
58 | TFile *histoFile = new TFile("DiFakes.root", "RECREATE"); | |
59 | ||
23035434 | 60 | TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.); |
30b3fe0f | 61 | TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.); |
62 | TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.); | |
23035434 | 63 | TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.); |
30b3fe0f | 64 | TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.); |
65 | TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.); | |
23035434 | 66 | TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.); |
30b3fe0f | 67 | TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.); |
68 | TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.); | |
23035434 | 69 | TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0); |
30b3fe0f | 70 | TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0); |
71 | TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0); | |
23035434 | 72 | TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0); |
30b3fe0f | 73 | TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0); |
74 | TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0); | |
23035434 | 75 | TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.); |
30b3fe0f | 76 | TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.); |
77 | TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.); | |
78 | ||
79 | // prepare for analysis | |
80 | AliMUONRecoParam *recoParam = 0x0; | |
81 | Double_t sigmaCut = -1; | |
82 | Prepare(recoParam, sigmaCut); | |
83 | ||
84 | // link to reconstructed tracks | |
85 | TFile* esdFile = TFile::Open(esdFileName); | |
86 | TTree* esdTree = GetESDTree(esdFile); | |
87 | AliESDEvent* esd = new AliESDEvent(); | |
88 | esd->ReadFromTree(esdTree); | |
89 | ||
90 | // link to simulated tracks | |
91 | AliMUONRecoCheck rc(esdFileName, SimDir); | |
92 | ||
93 | TLorentzVector vMu1, vMu2, vDiMu; | |
94 | ||
95 | // Loop over ESD events | |
96 | FirstEvent = TMath::Max(0, FirstEvent); | |
97 | LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1; | |
98 | for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) { | |
99 | ||
100 | // get the ESD of current event | |
101 | esdTree->GetEvent(iEvent); | |
102 | if (!esd) { | |
103 | Error("CheckESD", "no ESD object found for event %d", iEvent); | |
104 | return; | |
105 | } | |
106 | ||
107 | // convert TrackRef to MUON tracks | |
108 | AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent); | |
109 | ||
110 | // loop over ESD tracks and flag them | |
111 | Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; | |
112 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
113 | ||
114 | AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); | |
115 | ||
116 | // skip ghosts | |
117 | if (!muonTrack->ContainTrackerData()) continue; | |
118 | ||
119 | // try to match the reconstructed track with a simulated one | |
120 | Float_t fractionOfMatchCluster = 0.; | |
121 | AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut); | |
122 | ||
123 | // take actions according to matching result | |
124 | if (matchedTrackRef) { | |
125 | ||
126 | // flag matched tracks | |
127 | muonTrack->SetLabel(matchedTrackRef->GetUniqueID()); | |
128 | ||
129 | // remove already matched trackRefs | |
130 | trackRefStore->Remove(*matchedTrackRef); | |
131 | ||
132 | } else { | |
133 | ||
134 | // flag fake tracks | |
135 | muonTrack->SetLabel(-1); | |
136 | ||
137 | } | |
138 | ||
139 | } | |
140 | ||
141 | // double loop over ESD tracks, build pairs and fill histograms according to their label | |
142 | for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) { | |
143 | AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1); | |
144 | ||
145 | // skip ghosts | |
146 | if (!muonTrack1->ContainTrackerData()) continue; | |
147 | ||
148 | // get track info | |
149 | Short_t charge1 = muonTrack1->Charge(); | |
150 | Int_t label1 = muonTrack1->GetLabel(); | |
151 | muonTrack1->LorentzP(vMu1); | |
152 | ||
153 | for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) { | |
154 | AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2); | |
155 | ||
156 | // skip ghosts | |
157 | if (!muonTrack2->ContainTrackerData()) continue; | |
158 | ||
159 | // keep only opposite sign pairs | |
160 | Short_t charge2 = muonTrack2->Charge(); | |
161 | if (charge1*charge2 > 0) continue; | |
162 | ||
163 | // get track info | |
164 | Int_t label2 = muonTrack2->GetLabel(); | |
165 | muonTrack2->LorentzP(vMu2); | |
166 | ||
167 | // compute kinematics of the pair | |
168 | vDiMu = vMu1 + vMu2; | |
169 | Float_t mass = vDiMu.M(); | |
170 | Float_t p = vDiMu.P(); | |
171 | Float_t pt = vDiMu.Pt(); | |
172 | Float_t y = vDiMu.Rapidity(); | |
173 | Float_t eta = vDiMu.Eta(); | |
23035434 | 174 | Float_t phi = vDiMu.Phi(); |
30b3fe0f | 175 | if (phi < 0) phi += 2.*TMath::Pi(); |
176 | ||
177 | // fill global histograms | |
178 | hMass->Fill(mass); | |
179 | hP->Fill(p); | |
180 | hPt->Fill(pt); | |
181 | hY->Fill(y); | |
182 | hEta->Fill(eta); | |
183 | hPhi->Fill(phi); | |
184 | ||
185 | // fill histograms according to labels | |
186 | if (label1 >= 0 && label2 >= 0) { | |
187 | ||
188 | hMassM->Fill(mass); | |
189 | hPM->Fill(p); | |
190 | hPtM->Fill(pt); | |
191 | hYM->Fill(y); | |
192 | hEtaM->Fill(eta); | |
193 | hPhiM->Fill(phi); | |
194 | ||
195 | } else { | |
196 | ||
197 | hMassF->Fill(mass); | |
198 | hPF->Fill(p); | |
199 | hPtF->Fill(pt); | |
200 | hYF->Fill(y); | |
201 | hEtaF->Fill(eta); | |
202 | hPhiF->Fill(phi); | |
203 | ||
204 | } | |
205 | ||
206 | } | |
207 | ||
208 | } | |
209 | ||
210 | } // end of loop over events | |
211 | ||
212 | // plot results | |
213 | TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600); | |
214 | cDiFakesSummary.Divide(3,2); | |
215 | cDiFakesSummary.cd(1); | |
216 | cDiFakesSummary.GetPad(1)->SetLogy(); | |
217 | hMass->Draw(); | |
218 | hMass->SetMinimum(0.5); | |
219 | hMassM->Draw("same"); | |
220 | hMassM->SetLineColor(4); | |
221 | hMassF->Draw("same"); | |
222 | hMassF->SetLineColor(2); | |
223 | hMassF->SetFillColor(2); | |
224 | hMassF->SetFillStyle(3017); | |
225 | cDiFakesSummary.cd(2); | |
226 | cDiFakesSummary.GetPad(3)->SetLogy(); | |
227 | hP->Draw(); | |
228 | hP->SetMinimum(0.5); | |
229 | hPM->Draw("same"); | |
230 | hPM->SetLineColor(4); | |
231 | hPF->Draw("same"); | |
232 | hPF->SetLineColor(2); | |
233 | hPF->SetFillColor(2); | |
234 | hPF->SetFillStyle(3017); | |
235 | cDiFakesSummary.cd(3); | |
236 | cDiFakesSummary.GetPad(4)->SetLogy(); | |
237 | hPt->Draw(); | |
238 | hPt->SetMinimum(0.5); | |
239 | hPtM->Draw("same"); | |
240 | hPtM->SetLineColor(4); | |
241 | hPtF->Draw("same"); | |
242 | hPtF->SetLineColor(2); | |
243 | hPtF->SetFillColor(2); | |
244 | hPtF->SetFillStyle(3017); | |
245 | cDiFakesSummary.cd(4); | |
246 | cDiFakesSummary.GetPad(2)->SetLogy(); | |
247 | hY->Draw(); | |
248 | hY->SetMinimum(0.5); | |
249 | hYM->Draw("same"); | |
250 | hYM->SetLineColor(4); | |
251 | hYF->Draw("same"); | |
252 | hYF->SetLineColor(2); | |
253 | hYF->SetFillColor(2); | |
254 | hYF->SetFillStyle(3017); | |
255 | cDiFakesSummary.cd(5); | |
256 | cDiFakesSummary.GetPad(5)->SetLogy(); | |
257 | hEta->Draw(); | |
258 | hEta->SetMinimum(0.5); | |
259 | hEtaM->Draw("same"); | |
260 | hEtaM->SetLineColor(4); | |
261 | hEtaF->Draw("same"); | |
262 | hEtaF->SetLineColor(2); | |
263 | hEtaF->SetFillColor(2); | |
264 | hEtaF->SetFillStyle(3017); | |
265 | cDiFakesSummary.cd(6); | |
266 | cDiFakesSummary.GetPad(6)->SetLogy(); | |
267 | hPhi->Draw(); | |
268 | hPhi->SetMinimum(0.5); | |
269 | hPhiM->Draw("same"); | |
270 | hPhiM->SetLineColor(4); | |
271 | hPhiF->Draw("same"); | |
272 | hPhiF->SetLineColor(2); | |
273 | hPhiF->SetFillColor(2); | |
274 | hPhiF->SetFillStyle(3017); | |
275 | ||
276 | // save results | |
277 | histoFile->cd(); | |
278 | histoFile->Write(); | |
279 | cDiFakesSummary.Write(); | |
280 | histoFile->Close(); | |
281 | ||
282 | // clear memory | |
283 | delete esd; | |
284 | esdFile->Close(); | |
285 | delete recoParam; | |
286 | ||
287 | } | |
288 | ||
289 | //----------------------------------------------------------------------- | |
290 | void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut) | |
291 | { | |
292 | /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef | |
293 | ||
294 | // prepare OCDB access | |
295 | AliCDBManager* man = AliCDBManager::Instance(); | |
296 | man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); | |
297 | man->SetRun(0); | |
298 | ||
299 | // set mag field | |
300 | // waiting for mag field in CDB | |
301 | if (!TGeoGlobalMagField::Instance()->GetField()) { | |
302 | printf("Loading field map...\n"); | |
4642ac4b | 303 | AliMagF* field = new AliMagF("Maps","Maps",1.,1.,AliMagF::k5kG); |
30b3fe0f | 304 | TGeoGlobalMagField::Instance()->SetField(field); |
305 | } | |
306 | // set the magnetic field for track extrapolations | |
307 | AliMUONTrackExtrap::SetField(); | |
308 | ||
309 | // Load initial reconstruction parameters from OCDB | |
310 | AliCDBPath path("MUON","Calib","RecoParam"); | |
311 | AliCDBEntry *entry=man->Get(path.GetPath()); | |
312 | if(entry) { | |
313 | recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject()); | |
314 | entry->SetOwner(0); | |
315 | AliCDBManager::Instance()->UnloadFromCache(path.GetPath()); | |
316 | } | |
317 | if (!recoParam) { | |
318 | printf("Couldn't find RecoParam object in OCDB: create default one"); | |
319 | recoParam = AliMUONRecoParam::GetLowFluxParam(); | |
320 | } | |
321 | ||
322 | Info("MUONFakes", "\n recontruction parameters:"); | |
323 | recoParam->Print("FULL"); | |
324 | AliMUONESDInterface::ResetTracker(recoParam); | |
325 | ||
326 | // sigma cut to associate clusters with TrackRefs in case the label are not used | |
327 | sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); | |
328 | ||
329 | } | |
330 | ||
331 | //----------------------------------------------------------------------- | |
332 | TTree* GetESDTree(TFile *esdFile) | |
333 | { | |
334 | /// Check that the file is properly open | |
335 | /// Return pointer to the ESD Tree | |
336 | ||
337 | if (!esdFile || !esdFile->IsOpen()) { | |
338 | Error("GetESDTree", "opening ESD file failed"); | |
339 | exit(-1); | |
340 | } | |
341 | ||
342 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
343 | if (!tree) { | |
344 | Error("GetESDTree", "no ESD tree found"); | |
345 | exit(-1); | |
346 | } | |
347 | ||
348 | return tree; | |
349 | ||
350 | } | |
351 | ||
352 | //----------------------------------------------------------------------- | |
353 | Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut) | |
354 | { | |
355 | /// Try to match 2 tracks | |
356 | ||
357 | Bool_t compTrack[10]; | |
358 | Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack); | |
359 | fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters()); | |
360 | ||
361 | if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 | |
362 | (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 | |
363 | fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched | |
364 | else return kFALSE; | |
365 | ||
366 | } | |
367 | ||
368 | //----------------------------------------------------------------------- | |
369 | AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore, | |
370 | Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut) | |
371 | { | |
372 | /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters | |
373 | ||
374 | AliMUONTrack *matchedTrackRef = 0x0; | |
375 | fractionOfMatchCluster = 0.; | |
376 | ||
377 | if (useLabel) { // by using the MC label | |
378 | ||
379 | // get the corresponding simulated track if any | |
380 | Int_t label = muonTrack.GetLabel(); | |
381 | matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label); | |
382 | ||
383 | // get the fraction of matched clusters | |
384 | if (matchedTrackRef) { | |
385 | Int_t nMatchClusters = 0; | |
386 | if (muonTrack.ClustersStored()) { | |
387 | AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First(); | |
388 | while (cluster) { | |
389 | if (cluster->GetLabel() == label) nMatchClusters++; | |
390 | cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster); | |
391 | } | |
392 | } | |
393 | fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters()); | |
394 | } | |
395 | ||
396 | } else { // by comparing cluster/TrackRef positions | |
397 | ||
398 | // convert ESD track to MUON track | |
399 | AliMUONTrack track; | |
400 | AliMUONESDInterface::ESDToMUON(muonTrack,track); | |
401 | ||
402 | // look for the corresponding simulated track if any | |
403 | TIter next(trackRefStore.CreateIterator()); | |
404 | AliMUONTrack* trackRef; | |
405 | while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) { | |
406 | ||
407 | // check compatibility | |
408 | Float_t f = 0.; | |
409 | if (TrackMatched(track, *trackRef, f, sigmaCut)) { | |
410 | matchedTrackRef = trackRef; | |
411 | fractionOfMatchCluster = f; | |
412 | break; | |
413 | } | |
414 | ||
415 | } | |
416 | ||
417 | } | |
418 | ||
419 | return matchedTrackRef; | |
420 | ||
421 | } | |
422 |