1 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <TClonesArray.h>
13 #include "AliESDEvent.h"
14 #include "AliESDMuonTrack.h"
15 #include "AliCDBManager.h"
18 #include "AliMUONCDB.h"
19 #include "AliMUONTrack.h"
20 #include "AliMUONVTrackStore.h"
21 #include "AliMUONTrackParam.h"
22 #include "AliMUONESDInterface.h"
23 #include "AliMUONRecoCheck.h"
24 #include "AliMUONVCluster.h"
25 #include "AliMUONRecoParam.h"
31 /// \author Ph. Pillot, Subatech, March. 2009
33 /// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
34 /// Results are saved in the root file Fakes.root
35 /// Results are relevent provided that you use the same recoParams as for the reconstruction
37 UInt_t requestedStationMask = 0;
38 Bool_t request2ChInSameSt45 = kFALSE;
39 Double_t sigmaCut = -1.;
41 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
42 Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
44 //-----------------------------------------------------------------------
45 void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
46 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
47 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
50 //Reset ROOT and connect tree file
53 // File for histograms and histogram booking
54 TFile *histoFile = new TFile("Fakes.root", "RECREATE");
56 TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
57 TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
59 TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
60 TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
61 TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
62 TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
63 TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
64 TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
66 TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
67 TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
68 TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
69 TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
70 TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
71 TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
72 TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
73 TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
74 TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
75 TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
76 TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
77 TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
78 TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
79 TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
80 TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
82 // link to reconstructed and simulated tracks
83 AliMUONRecoCheck rc(esdFileName, SimDir);
85 // load necessary data from OCDB
86 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
87 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
88 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
89 if (!recoParam) return;
91 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
92 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
93 // compute the mask of requested stations from recoParam
94 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
95 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
96 request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
98 // initialize global counters
99 Int_t nReconstructibleTracks = 0;
100 Int_t nReconstructedTracks = 0;
101 Int_t nEventsWithTrackReconstructedYet = 0;
102 Int_t nEventsWithFake = 0;
103 Int_t nEventsWithAdditionalFake = 0;
104 Int_t nTotMatchedTracks = 0;
105 Int_t nTotTracksReconstructedYet = 0;
106 Int_t nTotFakeTracks = 0;
107 Int_t nTotConnectedTracks = 0;
108 Int_t nTotAdditionalTracks = 0;
109 Bool_t trackReconstructedYet;
110 TArrayI eventsWithTrackReconstructedYet(10);
111 TArrayI eventsWithFake(10);
112 TArrayI eventsWithAdditionalFake(10);
114 // Loop over ESD events
115 FirstEvent = TMath::Max(0, FirstEvent);
116 LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
117 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
119 // get reconstructed and simulated tracks
120 AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
121 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
122 if (!muonTrackStore || !trackRefStore) continue;
124 // count the number of reconstructible tracks
125 TIter next(trackRefStore->CreateIterator());
126 AliMUONTrack* trackRef;
127 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
128 if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
131 // loop over ESD tracks
132 Int_t nTrackerTracks = 0;
133 trackReconstructedYet = kFALSE;
134 AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
135 const AliESDEvent* esd = rc.GetESDEvent();
136 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
137 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
139 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
142 if (!esdTrack->ContainTrackerData()) continue;
145 // find the corresponding MUON track
146 AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
149 Int_t nClusters = esdTrack->GetNClusters();
150 Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
151 Double_t p = esdTrack->P();
152 Double_t pT = esdTrack->Pt();
153 Double_t eta = esdTrack->Eta();
154 Double_t phi = esdTrack->Phi();
156 // fill global histograms
157 hNumberOfClusters->Fill(nClusters);
158 hChi2PerDof->Fill(normalizedChi2);
164 // try to match the reconstructed track with a simulated one
165 Int_t nMatchClusters = 0;
166 AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
168 // take actions according to matching result
169 if (matchedTrackRef) {
173 if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
174 trackReconstructedYet = kTRUE;
175 nTotTracksReconstructedYet++;
179 hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
180 hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
181 hNumberOfClustersM->Fill(nClusters);
182 hChi2PerDofM->Fill(normalizedChi2);
188 // remove already matched trackRefs
189 trackRefStore->Remove(*matchedTrackRef);
197 hNumberOfClustersF->Fill(nClusters);
198 hChi2PerDofF->Fill(normalizedChi2);
205 fakeTrackStore->Add(*muonTrack);
209 } // end of loop over ESD tracks
212 hNumberOfTracks->Fill(nTrackerTracks);
213 nReconstructedTracks += nTrackerTracks;
214 if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
216 // count the number the additional fake tracks
217 if (fakeTrackStore->GetSize() > 0) {
219 // remove the most connected fake tracks
220 Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, useLabel, *hFractionOfConnectedClusters);
222 // remove the remaining free reconstructible tracks
223 Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
226 eventsWithFake[nEventsWithFake] = iEvent;
228 if (nAdditionalTracks > 0) {
229 eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
230 nEventsWithAdditionalFake++;
231 nTotAdditionalTracks += nAdditionalTracks;
232 hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
237 delete fakeTrackStore;
239 } // end of loop over events
241 // total number of connected tracks
242 nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
245 TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
246 cFakesSummary.Divide(3,2);
248 cFakesSummary.GetPad(1)->SetLogy();
249 hNumberOfClusters->Draw();
250 hNumberOfClusters->SetMinimum(0.5);
251 hNumberOfClustersM->Draw("same");
252 hNumberOfClustersM->SetLineColor(4);
253 hNumberOfClustersF->Draw("same");
254 hNumberOfClustersF->SetLineColor(2);
255 hNumberOfClustersF->SetFillColor(2);
256 hNumberOfClustersF->SetFillStyle(3017);
258 cFakesSummary.GetPad(2)->SetLogy();
260 hChi2PerDof->SetMinimum(0.5);
261 hChi2PerDofM->Draw("same");
262 hChi2PerDofM->SetLineColor(4);
263 hChi2PerDofF->Draw("same");
264 hChi2PerDofF->SetLineColor(2);
265 hChi2PerDofF->SetFillColor(2);
266 hChi2PerDofF->SetFillStyle(3017);
268 cFakesSummary.GetPad(3)->SetLogy();
272 hPM->SetLineColor(4);
274 hPF->SetLineColor(2);
275 hPF->SetFillColor(2);
276 hPF->SetFillStyle(3017);
278 cFakesSummary.GetPad(4)->SetLogy();
280 hPt->SetMinimum(0.5);
282 hPtM->SetLineColor(4);
284 hPtF->SetLineColor(2);
285 hPtF->SetFillColor(2);
286 hPtF->SetFillStyle(3017);
288 cFakesSummary.GetPad(5)->SetLogy();
290 hEta->SetMinimum(0.5);
292 hEtaM->SetLineColor(4);
294 hEtaF->SetLineColor(2);
295 hEtaF->SetFillColor(2);
296 hEtaF->SetFillStyle(3017);
298 cFakesSummary.GetPad(6)->SetLogy();
300 hPhi->SetMinimum(0.5);
302 hPhiM->SetLineColor(4);
304 hPhiF->SetLineColor(2);
305 hPhiF->SetFillColor(2);
306 hPhiF->SetFillStyle(3017);
311 cFakesSummary.Write();
316 cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
317 cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
318 cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
319 cout << " (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
320 if (nTotTracksReconstructedYet > 0) {
321 for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
322 if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
323 else cout << ", " << eventsWithTrackReconstructedYet[i];
325 cout << "))" << endl;
326 } else cout << ")" << endl;
327 cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
328 cout << " (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
329 cout << " (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
330 cout << "- Number of events with fake track(s): " << nEventsWithFake;
331 if (nEventsWithFake > 0) {
332 for(Int_t i=0; i<nEventsWithFake; i++){
333 if (i==0) cout << " (eventID = " << eventsWithFake[i];
334 else cout << ", " << eventsWithFake[i];
338 cout << " (including " << nEventsWithAdditionalFake << " events with additional track(s)";
339 if (nEventsWithAdditionalFake > 0) {
340 for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
341 if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
342 else cout << ", " << eventsWithAdditionalFake[i];
344 cout << "))" << endl;
345 } else cout << ")" << endl;
347 cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
352 //-----------------------------------------------------------------------
353 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
354 Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
356 /// loop over reconstructible TrackRef not associated with reconstructed track:
357 /// for each of them, find and remove the most connected the fake track, if any,
358 /// and fill the histograms with the fraction of connected clusters.
359 /// Return the number of reconstructible track not connected to any fake
361 Int_t nFreeMissingTracks = 0;
363 // loop over trackRefs
364 TIter next(trackRefStore.CreateIterator());
365 AliMUONTrack* trackRef;
366 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
368 // skip not reconstructible trackRefs
369 if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
371 Int_t label = trackRef->GetUniqueID();
373 // look for the most connected fake track
374 AliMUONTrack *connectedFake = 0x0;
375 Double_t fractionOfConnectedClusters = 0.;
376 TIter next2(fakeTrackStore.CreateIterator());
377 AliMUONTrack* fakeTrack;
378 while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
380 // get the number of connected clusters
381 Int_t nConnectedClusters = 0;
382 if (useLabel) { // by using the MC label
383 for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
384 if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
385 nConnectedClusters++;
386 } else { // by comparing cluster/TrackRef positions
387 Bool_t compTrack[10];
388 nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
391 // skip non-connected fake tracks
392 if (nConnectedClusters == 0) continue;
394 // check if it is the most connected fake track
395 Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
396 if (f > fractionOfConnectedClusters) {
397 connectedFake = fakeTrack;
398 fractionOfConnectedClusters = f;
403 // remove the most connected fake track
405 hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
406 fakeTrackStore.Remove(*connectedFake);
407 } else nFreeMissingTracks++;
411 return nFreeMissingTracks;