1 #if !defined(__CINT__) || defined(__MAKECINT__)
9 #include <TClonesArray.h>
10 #include <TGeoGlobalMagField.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"
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"
36 /// \author Ph. Pillot, Subatech, March. 2009
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
42 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
43 TTree* GetESDTree(TFile *esdFile);
44 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
45 Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
46 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
47 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
48 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
49 Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
51 //-----------------------------------------------------------------------
52 void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
53 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
56 //Reset ROOT and connect tree file
59 // File for histograms and histogram booking
60 TFile *histoFile = new TFile("Fakes.root", "RECREATE");
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.);
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);
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.);
88 // prepare for analysis
89 AliMUONRecoParam *recoParam = 0x0;
90 Double_t sigmaCut = -1;
91 Prepare(recoParam, sigmaCut);
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);
99 // link to simulated tracks
100 AliMUONRecoCheck rc(esdFileName, SimDir);
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;
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++) {
115 // get the ESD of current event
116 esdTree->GetEvent(iEvent);
118 Error("CheckESD", "no ESD object found for event %d", iEvent);
122 // convert TrackRef to MUON tracks
123 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
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++) {
131 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
134 if (!muonTrack->ContainTrackerData()) continue;
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();
145 // fill global histograms
146 hNumberOfClusters->Fill(nClusters);
147 hChi2PerDof->Fill(normalizedChi2);
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);
157 // take actions according to matching result
158 if (matchedTrackRef) {
162 if (!IsRecontructible(*matchedTrackRef,*recoParam)) nTotTracksReconstructedYet++;
165 hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
166 hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
167 hNumberOfClustersM->Fill(nClusters);
168 hChi2PerDofM->Fill(normalizedChi2);
174 // remove already matched trackRefs
175 trackRefStore->Remove(*matchedTrackRef);
183 hNumberOfClustersF->Fill(nClusters);
184 hChi2PerDofF->Fill(normalizedChi2);
191 AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
195 } // end of loop over ESD tracks
198 hNumberOfTracks->Fill(nTrackerTracks);
200 // count the number the additional fake tracks
201 if (fakeTrackStore->GetSize() > 0) {
203 // remove the most connected fake tracks
204 Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
205 useLabel, sigmaCut, *hFractionOfConnectedClusters);
207 // remove the remaining free reconstructible tracks
208 Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
212 if (nAdditionalTracks > 0) {
213 nEventsWithAdditionalFake++;
214 nTotAdditionalTracks += nAdditionalTracks;
215 hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
220 delete fakeTrackStore;
222 } // end of loop over events
225 TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
226 cFakesSummary.Divide(3,2);
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);
238 cFakesSummary.GetPad(2)->SetLogy();
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);
248 cFakesSummary.GetPad(3)->SetLogy();
252 hPM->SetLineColor(4);
254 hPF->SetLineColor(2);
255 hPF->SetFillColor(2);
256 hPF->SetFillStyle(3017);
258 cFakesSummary.GetPad(4)->SetLogy();
260 hPt->SetMinimum(0.5);
262 hPtM->SetLineColor(4);
264 hPtF->SetLineColor(2);
265 hPtF->SetFillColor(2);
266 hPtF->SetFillStyle(3017);
268 cFakesSummary.GetPad(5)->SetLogy();
270 hEta->SetMinimum(0.5);
272 hEtaM->SetLineColor(4);
274 hEtaF->SetLineColor(2);
275 hEtaF->SetFillColor(2);
276 hEtaF->SetFillStyle(3017);
278 cFakesSummary.GetPad(6)->SetLogy();
280 hPhi->SetMinimum(0.5);
282 hPhiM->SetLineColor(4);
284 hPhiF->SetLineColor(2);
285 hPhiF->SetFillColor(2);
286 hPhiF->SetFillStyle(3017);
291 cFakesSummary.Write();
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;
303 cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
313 //-----------------------------------------------------------------------
314 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
316 /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
318 // prepare OCDB access
319 AliCDBManager* man = AliCDBManager::Instance();
320 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
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);
330 // set the magnetic field for track extrapolations
331 AliMUONTrackExtrap::SetField();
333 // Load initial reconstruction parameters from OCDB
334 AliCDBPath path("MUON","Calib","RecoParam");
335 AliCDBEntry *entry=man->Get(path.GetPath());
337 recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
339 AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
342 printf("Couldn't find RecoParam object in OCDB: create default one");
343 recoParam = AliMUONRecoParam::GetLowFluxParam();
346 Info("MUONFakes", "\n recontruction parameters:");
347 recoParam->Print("FULL");
348 AliMUONESDInterface::ResetTracker(recoParam);
350 // sigma cut to associate clusters with TrackRefs in case the label are not used
351 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
355 //-----------------------------------------------------------------------
356 TTree* GetESDTree(TFile *esdFile)
358 /// Check that the file is properly open
359 /// Return pointer to the ESD Tree
361 if (!esdFile || !esdFile->IsOpen()) {
362 Error("GetESDTree", "opening ESD file failed");
366 TTree* tree = (TTree*) esdFile->Get("esdTree");
368 Error("GetESDTree", "no ESD tree found");
376 //-----------------------------------------------------------------------
377 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
379 /// Try to match 2 tracks
381 Bool_t compTrack[10];
382 Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
383 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
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
392 //-----------------------------------------------------------------------
393 Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
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);
401 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
404 currentCh = trackParam->GetClusterPtr()->GetChamberId();
406 clusterInSt[currentCh/2] = kTRUE;
408 if (currentCh > 5 && currentCh != previousCh) {
410 previousCh = currentCh;
413 trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
416 return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
417 clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
421 //-----------------------------------------------------------------------
422 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
423 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
425 /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
427 AliMUONTrack *matchedTrackRef = 0x0;
428 fractionOfMatchCluster = 0.;
430 if (useLabel) { // by using the MC label
432 // get the corresponding simulated track if any
433 Int_t label = muonTrack.GetLabel();
434 matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
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();
442 if (cluster->GetLabel() == label) nMatchClusters++;
443 cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
446 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
449 } else { // by comparing cluster/TrackRef positions
451 // convert ESD track to MUON track
453 AliMUONESDInterface::ESDToMUON(muonTrack,track);
455 // look for the corresponding simulated track if any
456 TIter next(trackRefStore.CreateIterator());
457 AliMUONTrack* trackRef;
458 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
460 // check compatibility
462 if (TrackMatched(track, *trackRef, f, sigmaCut)) {
463 matchedTrackRef = trackRef;
464 fractionOfMatchCluster = f;
472 return matchedTrackRef;
476 //-----------------------------------------------------------------------
477 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
478 Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
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
485 Int_t nFreeMissingTracks = 0;
487 // loop over trackRefs
488 TIter next(trackRefStore.CreateIterator());
489 AliMUONTrack* trackRef;
490 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
492 // skip not reconstructible trackRefs
493 if (!IsRecontructible(*trackRef,recoParam)) continue;
495 Int_t label = trackRef->GetUniqueID();
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()) ) ) {
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);
515 // skip non-connected fake tracks
516 if (nConnectedClusters == 0) continue;
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;
527 // remove the most connected fake track
529 hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
530 fakeTrackStore.Remove(*connectedFake);
531 } else nFreeMissingTracks++;
535 return nFreeMissingTracks;