Macro to study fake tracks by comparing reconstructed tracks with TrackRefs.
[u/mrichter/AliRoot.git] / MUON / MUONFakes.C
CommitLineData
66085093 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
12// STEER includes
13#include "AliLog.h"
14#include "AliMagF.h"
15#include "AliESDEvent.h"
16#include "AliESDMuonTrack.h"
17#include "AliESDMuonCluster.h"
18#include "AliCDBPath.h"
19#include "AliCDBEntry.h"
20#include "AliCDBManager.h"
21
22// MUON includes
23#include "AliMUONTrack.h"
24#include "AliMUONVTrackStore.h"
25#include "AliMUONTrackParam.h"
26#include "AliMUONTrackExtrap.h"
27#include "AliMUONESDInterface.h"
28#include "AliMUONRecoCheck.h"
29#include "AliMUONVCluster.h"
30#include "AliMUONRecoParam.h"
31#endif
32
33/// \ingroup macros
34/// \file MUONFakes.C
35///
36/// \author Ph. Pillot, Subatech, March. 2009
37///
38/// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
39/// Results are saved in the root file Fakes.root
40/// Results are relevent provided that you use the same recoParams as for the reconstruction
41
42void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
43TTree* GetESDTree(TFile *esdFile);
44Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
45Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
46AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
47 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
48Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
49 Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
50
51//-----------------------------------------------------------------------
52void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
53 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
54{
55
56 //Reset ROOT and connect tree file
57 gROOT->Reset();
58
59 // File for histograms and histogram booking
60 TFile *histoFile = new TFile("Fakes.root", "RECREATE");
61
62 TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
63 TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
64
65 TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
66 TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
67 TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
68 TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
69 TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
70 TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
71
72 TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
73 TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
74 TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
75 TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
76 TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
77 TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
78 TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
79 TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
80 TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
81 TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
82 TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
83 TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
84 TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
85 TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
86 TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
87
88 // prepare for analysis
89 AliMUONRecoParam *recoParam = 0x0;
90 Double_t sigmaCut = -1;
91 Prepare(recoParam, sigmaCut);
92
93 // link to reconstructed tracks
94 TFile* esdFile = TFile::Open(esdFileName);
95 TTree* esdTree = GetESDTree(esdFile);
96 AliESDEvent* esd = new AliESDEvent();
97 esd->ReadFromTree(esdTree);
98
99 // link to simulated tracks
100 AliMUONRecoCheck rc(esdFileName, SimDir);
101
102 // initialize global counters
103 Int_t nEventsWithFake = 0;
104 Int_t nEventsWithAdditionalFake = 0;
105 Int_t nTotMatchedTracks = 0;
106 Int_t nTotTracksReconstructedYet = 0;
107 Int_t nTotFakeTracks = 0;
108 Int_t nTotAdditionalTracks = 0;
109
110 // Loop over ESD events
111 FirstEvent = TMath::Max(0, FirstEvent);
112 LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
113 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
114
115 // get the ESD of current event
116 esdTree->GetEvent(iEvent);
117 if (!esd) {
118 Error("CheckESD", "no ESD object found for event %d", iEvent);
119 return;
120 }
121
122 // convert TrackRef to MUON tracks
123 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
124
125 // loop over ESD tracks
126 Int_t nTrackerTracks = 0;
127 AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
128 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
129 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
130
131 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
132
133 // skip ghosts
134 if (!muonTrack->ContainTrackerData()) continue;
135 nTrackerTracks++;
136
137 // get track info
138 Int_t nClusters = muonTrack->GetNClusters();
139 Double_t normalizedChi2 = muonTrack->GetChi2() / (2. * muonTrack->GetNHit() - 5);
140 Double_t p = muonTrack->P();
141 Double_t pT = muonTrack->Pt();
142 Double_t eta = muonTrack->Eta();
143 Double_t phi = muonTrack->Phi();
144
145 // fill global histograms
146 hNumberOfClusters->Fill(nClusters);
147 hChi2PerDof->Fill(normalizedChi2);
148 hP->Fill(p);
149 hPt->Fill(pT);
150 hEta->Fill(eta);
151 hPhi->Fill(phi);
152
153 // try to match the reconstructed track with a simulated one
154 Float_t fractionOfMatchCluster = 0.;
155 AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
156
157 // take actions according to matching result
158 if (matchedTrackRef) {
159
160 // global counter
161 nTotMatchedTracks++;
162 if (!IsRecontructible(*matchedTrackRef,*recoParam)) nTotTracksReconstructedYet++;
163
164 // fill histograms
165 hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
166 hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
167 hNumberOfClustersM->Fill(nClusters);
168 hChi2PerDofM->Fill(normalizedChi2);
169 hPM->Fill(p);
170 hPtM->Fill(pT);
171 hEtaM->Fill(eta);
172 hPhiM->Fill(phi);
173
174 // remove already matched trackRefs
175 trackRefStore->Remove(*matchedTrackRef);
176
177 } else {
178
179 // global counter
180 nTotFakeTracks++;
181
182 // fill histograms
183 hNumberOfClustersF->Fill(nClusters);
184 hChi2PerDofF->Fill(normalizedChi2);
185 hPF->Fill(p);
186 hPtF->Fill(pT);
187 hEtaF->Fill(eta);
188 hPhiF->Fill(phi);
189
190 // store fake tracks
191 AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
192
193 }
194
195 } // end of loop over ESD tracks
196
197 // fill histograms
198 hNumberOfTracks->Fill(nTrackerTracks);
199
200 // count the number the additional fake tracks
201 if (fakeTrackStore->GetSize() > 0) {
202
203 // remove the most connected fake tracks
204 Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
205 useLabel, sigmaCut, *hFractionOfConnectedClusters);
206
207 // remove the remaining free reconstructible tracks
208 Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
209
210 // fill histograms
211 nEventsWithFake++;
212 if (nAdditionalTracks > 0) {
213 nEventsWithAdditionalFake++;
214 nTotAdditionalTracks += nAdditionalTracks;
215 hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
216 }
217
218 }
219
220 delete fakeTrackStore;
221
222 } // end of loop over events
223
224 // plot results
225 TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
226 cFakesSummary.Divide(3,2);
227 cFakesSummary.cd(1);
228 cFakesSummary.GetPad(1)->SetLogy();
229 hNumberOfClusters->Draw();
230 hNumberOfClusters->SetMinimum(0.5);
231 hNumberOfClustersM->Draw("same");
232 hNumberOfClustersM->SetLineColor(4);
233 hNumberOfClustersF->Draw("same");
234 hNumberOfClustersF->SetLineColor(2);
235 hNumberOfClustersF->SetFillColor(2);
236 hNumberOfClustersF->SetFillStyle(3017);
237 cFakesSummary.cd(2);
238 cFakesSummary.GetPad(2)->SetLogy();
239 hChi2PerDof->Draw();
240 hChi2PerDof->SetMinimum(0.5);
241 hChi2PerDofM->Draw("same");
242 hChi2PerDofM->SetLineColor(4);
243 hChi2PerDofF->Draw("same");
244 hChi2PerDofF->SetLineColor(2);
245 hChi2PerDofF->SetFillColor(2);
246 hChi2PerDofF->SetFillStyle(3017);
247 cFakesSummary.cd(3);
248 cFakesSummary.GetPad(3)->SetLogy();
249 hP->Draw();
250 hP->SetMinimum(0.5);
251 hPM->Draw("same");
252 hPM->SetLineColor(4);
253 hPF->Draw("same");
254 hPF->SetLineColor(2);
255 hPF->SetFillColor(2);
256 hPF->SetFillStyle(3017);
257 cFakesSummary.cd(4);
258 cFakesSummary.GetPad(4)->SetLogy();
259 hPt->Draw();
260 hPt->SetMinimum(0.5);
261 hPtM->Draw("same");
262 hPtM->SetLineColor(4);
263 hPtF->Draw("same");
264 hPtF->SetLineColor(2);
265 hPtF->SetFillColor(2);
266 hPtF->SetFillStyle(3017);
267 cFakesSummary.cd(5);
268 cFakesSummary.GetPad(5)->SetLogy();
269 hEta->Draw();
270 hEta->SetMinimum(0.5);
271 hEtaM->Draw("same");
272 hEtaM->SetLineColor(4);
273 hEtaF->Draw("same");
274 hEtaF->SetLineColor(2);
275 hEtaF->SetFillColor(2);
276 hEtaF->SetFillStyle(3017);
277 cFakesSummary.cd(6);
278 cFakesSummary.GetPad(6)->SetLogy();
279 hPhi->Draw();
280 hPhi->SetMinimum(0.5);
281 hPhiM->Draw("same");
282 hPhiM->SetLineColor(4);
283 hPhiF->Draw("same");
284 hPhiF->SetLineColor(2);
285 hPhiF->SetFillColor(2);
286 hPhiF->SetFillStyle(3017);
287
288 // save results
289 histoFile->cd();
290 histoFile->Write();
291 cFakesSummary.Write();
292 histoFile->Close();
293
294 // print results
295 cout << endl;
296 cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
297 cout << " (including " << nTotTracksReconstructedYet << " tracks matched with a TrackRef that is not reconstructible)" << endl;
298 cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
299 cout << " (including " << nTotAdditionalTracks << " additional tracks (compared to the number of expected ones))" << endl;
300 cout << "- Number of events with fake track(s): " << nEventsWithFake << endl;
301 cout << " (including " << nEventsWithAdditionalFake << " events with additional tracks)" << endl;
302 cout << endl;
303 cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
304 cout << endl;
305
306 // clear memory
307 delete esd;
308 esdFile->Close();
309 delete recoParam;
310
311}
312
313//-----------------------------------------------------------------------
314void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
315{
316 /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
317
318 // prepare OCDB access
319 AliCDBManager* man = AliCDBManager::Instance();
320 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
321 man->SetRun(0);
322
323 // set mag field
324 // waiting for mag field in CDB
325 if (!TGeoGlobalMagField::Instance()->GetField()) {
326 printf("Loading field map...\n");
327 AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
328 TGeoGlobalMagField::Instance()->SetField(field);
329 }
330 // set the magnetic field for track extrapolations
331 AliMUONTrackExtrap::SetField();
332
333 // Load initial reconstruction parameters from OCDB
334 AliCDBPath path("MUON","Calib","RecoParam");
335 AliCDBEntry *entry=man->Get(path.GetPath());
336 if(entry) {
337 recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
338 entry->SetOwner(0);
339 AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
340 }
341 if (!recoParam) {
342 printf("Couldn't find RecoParam object in OCDB: create default one");
343 recoParam = AliMUONRecoParam::GetLowFluxParam();
344 }
345
346 Info("MUONFakes", "\n recontruction parameters:");
347 recoParam->Print("FULL");
348 AliMUONESDInterface::ResetTracker(recoParam);
349
350 // sigma cut to associate clusters with TrackRefs in case the label are not used
351 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
352
353}
354
355//-----------------------------------------------------------------------
356TTree* GetESDTree(TFile *esdFile)
357{
358 /// Check that the file is properly open
359 /// Return pointer to the ESD Tree
360
361 if (!esdFile || !esdFile->IsOpen()) {
362 Error("GetESDTree", "opening ESD file failed");
363 exit(-1);
364 }
365
366 TTree* tree = (TTree*) esdFile->Get("esdTree");
367 if (!tree) {
368 Error("GetESDTree", "no ESD tree found");
369 exit(-1);
370 }
371
372 return tree;
373
374}
375
376//-----------------------------------------------------------------------
377Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
378{
379 /// Try to match 2 tracks
380
381 Bool_t compTrack[10];
382 Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
383 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
384
385 if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
386 (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
387 fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched
388 else return kFALSE;
389
390}
391
392//-----------------------------------------------------------------------
393Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
394{
395 /// Check il the track is reconstructible
396 Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3;
397 Int_t currentCh, previousCh = -1, nChHitInSt45 = 0;
398 Bool_t clusterInSt[5];
399 for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt);
400
401 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
402 while (trackParam) {
403
404 currentCh = trackParam->GetClusterPtr()->GetChamberId();
405
406 clusterInSt[currentCh/2] = kTRUE;
407
408 if (currentCh > 5 && currentCh != previousCh) {
409 nChHitInSt45++;
410 previousCh = currentCh;
411 }
412
413 trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
414 }
415
416 return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
417 clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
418
419}
420
421//-----------------------------------------------------------------------
422AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
423 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
424{
425 /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
426
427 AliMUONTrack *matchedTrackRef = 0x0;
428 fractionOfMatchCluster = 0.;
429
430 if (useLabel) { // by using the MC label
431
432 // get the corresponding simulated track if any
433 Int_t label = muonTrack.GetLabel();
434 matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
435
436 // get the fraction of matched clusters
437 if (matchedTrackRef) {
438 Int_t nMatchClusters = 0;
439 if (muonTrack.ClustersStored()) {
440 AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
441 while (cluster) {
442 if (cluster->GetLabel() == label) nMatchClusters++;
443 cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
444 }
445 }
446 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
447 }
448
449 } else { // by comparing cluster/TrackRef positions
450
451 // convert ESD track to MUON track
452 AliMUONTrack track;
453 AliMUONESDInterface::ESDToMUON(muonTrack,track);
454
455 // look for the corresponding simulated track if any
456 TIter next(trackRefStore.CreateIterator());
457 AliMUONTrack* trackRef;
458 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
459
460 // check compatibility
461 Float_t f = 0.;
462 if (TrackMatched(track, *trackRef, f, sigmaCut)) {
463 matchedTrackRef = trackRef;
464 fractionOfMatchCluster = f;
465 break;
466 }
467
468 }
469
470 }
471
472 return matchedTrackRef;
473
474}
475
476//-----------------------------------------------------------------------
477Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
478 Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
479{
480 /// loop over reconstructible TrackRef not associated with reconstructed track:
481 /// for each of them, find and remove the most connected the fake track, if any,
482 /// and fill the histograms with the fraction of connected clusters.
483 /// Return the number of reconstructible track not connected to any fake
484
485 Int_t nFreeMissingTracks = 0;
486
487 // loop over trackRefs
488 TIter next(trackRefStore.CreateIterator());
489 AliMUONTrack* trackRef;
490 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
491
492 // skip not reconstructible trackRefs
493 if (!IsRecontructible(*trackRef,recoParam)) continue;
494
495 Int_t label = trackRef->GetUniqueID();
496
497 // look for the most connected fake track
498 AliMUONTrack *connectedFake = 0x0;
499 Float_t fractionOfConnectedClusters = 0.;
500 TIter next2(fakeTrackStore.CreateIterator());
501 AliMUONTrack* fakeTrack;
502 while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
503
504 // get the number of connected clusters
505 Int_t nConnectedClusters = 0;
506 if (useLabel) { // by using the MC label
507 for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
508 if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
509 nConnectedClusters++;
510 } else { // by comparing cluster/TrackRef positions
511 Bool_t compTrack[10];
512 nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack);
513 }
514
515 // skip non-connected fake tracks
516 if (nConnectedClusters == 0) continue;
517
518 // check if it is the most connected fake track
519 Float_t f = ((Float_t)nConnectedClusters) / ((Float_t)fakeTrack->GetNClusters());
520 if (f > fractionOfConnectedClusters) {
521 connectedFake = fakeTrack;
522 fractionOfConnectedClusters = f;
523 }
524
525 }
526
527 // remove the most connected fake track
528 if (connectedFake) {
529 hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
530 fakeTrackStore.Remove(*connectedFake);
531 } else nFreeMissingTracks++;
532
533 }
534
535 return nFreeMissingTracks;
536
537}
538