Macro to study the effects of fake tracks on the dimuon spectra.
[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.C
CommitLineData
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
43void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
44TTree* GetESDTree(TFile *esdFile);
45Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
46AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
47 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
48
49//-----------------------------------------------------------------------
50void 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
60 TH1F *hMass = new TH1F("hMass", "Muon mass distribution (GeV/c^{2})", 100, 0., 12.);
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.);
63 TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
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.);
66 TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
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.);
69 TH1F *hY = new TH1F("hY"," Muon rapidity distribution",100,-10,0);
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);
72 TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
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);
75 TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
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();
174 Float_t phi = vMu1.Phi();
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//-----------------------------------------------------------------------
290void 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");
303 AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
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//-----------------------------------------------------------------------
332TTree* 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//-----------------------------------------------------------------------
353Bool_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//-----------------------------------------------------------------------
369AliMUONTrack* 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