renaming function to avoid library conflict (discovered in test/generators/TUHKMgen)
[u/mrichter/AliRoot.git] / MUON / MUONFakes.C
CommitLineData
66085093 1#if !defined(__CINT__) || defined(__MAKECINT__)
2// ROOT includes
66085093 3#include <TFile.h>
4#include <TH1.h>
5#include <TCanvas.h>
6#include <Riostream.h>
7#include <TROOT.h>
ce8cd162 8#include <TObjArray.h>
23035434 9#include <TArrayI.h>
66085093 10
11// STEER includes
12#include "AliLog.h"
66085093 13#include "AliESDEvent.h"
14#include "AliESDMuonTrack.h"
66085093 15#include "AliCDBManager.h"
16
17// MUON includes
a99c3449 18#include "AliMUONCDB.h"
66085093 19#include "AliMUONTrack.h"
20#include "AliMUONVTrackStore.h"
21#include "AliMUONTrackParam.h"
66085093 22#include "AliMUONESDInterface.h"
23#include "AliMUONRecoCheck.h"
24#include "AliMUONVCluster.h"
25#include "AliMUONRecoParam.h"
26#endif
27
28/// \ingroup macros
29/// \file MUONFakes.C
30///
31/// \author Ph. Pillot, Subatech, March. 2009
32///
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
36
f202486b 37UInt_t requestedStationMask = 0;
38Bool_t request2ChInSameSt45 = kFALSE;
39Double_t sigmaCut = -1.;
40
41Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
42 Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
66085093 43
44//-----------------------------------------------------------------------
45void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
a99c3449 46 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
47 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
66085093 48{
49
19a4c30d 50 // Disable printout of AliMCEvent
51 AliLog::SetClassDebugLevel("AliMCEvent",-1);
52
66085093 53 //Reset ROOT and connect tree file
54 gROOT->Reset();
55
56 // File for histograms and histogram booking
57 TFile *histoFile = new TFile("Fakes.root", "RECREATE");
58
59 TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
60 TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
61
62 TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
63 TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
64 TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
65 TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
66 TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
67 TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
68
69 TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
70 TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
71 TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
72 TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
73 TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
74 TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
75 TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
76 TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
77 TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
78 TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
79 TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
80 TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
81 TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
82 TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
83 TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
19a4c30d 84 TH1F *hNumberOfChamberHit = new TH1F("hNumberOfChamberHit","nb of chambers hit /track",20,0.,20.);
85 TH1F *hNumberOfChamberHitM = new TH1F("hNumberOfChamberHitM","nb of chambers hit /matched track",20,0.,20.);
86 TH1F *hNumberOfChamberHitF = new TH1F("hNumberOfChamberHitF","nb of chambers hit /fake track",20,0.,20.);
66085093 87
a99c3449 88 // link to reconstructed and simulated tracks
89 AliMUONRecoCheck rc(esdFileName, SimDir);
66085093 90
a99c3449 91 // load necessary data from OCDB
92 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
93 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
94 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
95 if (!recoParam) return;
66085093 96
19a4c30d 97 // get sigma cut from recoParam to associate clusters with TrackRefs in case the labels are not used
f202486b 98 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
99 // compute the mask of requested stations from recoParam
100 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
101 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
102 request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
66085093 103
104 // initialize global counters
23035434 105 Int_t nReconstructibleTracks = 0;
106 Int_t nReconstructedTracks = 0;
107 Int_t nEventsWithTrackReconstructedYet = 0;
66085093 108 Int_t nEventsWithFake = 0;
109 Int_t nEventsWithAdditionalFake = 0;
110 Int_t nTotMatchedTracks = 0;
111 Int_t nTotTracksReconstructedYet = 0;
112 Int_t nTotFakeTracks = 0;
23035434 113 Int_t nTotConnectedTracks = 0;
66085093 114 Int_t nTotAdditionalTracks = 0;
23035434 115 Bool_t trackReconstructedYet;
116 TArrayI eventsWithTrackReconstructedYet(10);
117 TArrayI eventsWithFake(10);
118 TArrayI eventsWithAdditionalFake(10);
66085093 119
120 // Loop over ESD events
121 FirstEvent = TMath::Max(0, FirstEvent);
a99c3449 122 LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
66085093 123 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
124
a99c3449 125 // get reconstructed and simulated tracks
126 AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
66085093 127 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
a99c3449 128 if (!muonTrackStore || !trackRefStore) continue;
66085093 129
23035434 130 // count the number of reconstructible tracks
131 TIter next(trackRefStore->CreateIterator());
132 AliMUONTrack* trackRef;
133 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
f202486b 134 if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
23035434 135 }
136
66085093 137 // loop over ESD tracks
138 Int_t nTrackerTracks = 0;
23035434 139 trackReconstructedYet = kFALSE;
66085093 140 AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
a99c3449 141 const AliESDEvent* esd = rc.GetESDEvent();
66085093 142 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
143 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
144
a99c3449 145 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
66085093 146
147 // skip ghosts
a99c3449 148 if (!esdTrack->ContainTrackerData()) continue;
66085093 149 nTrackerTracks++;
150
a99c3449 151 // find the corresponding MUON track
152 AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
153
66085093 154 // get track info
a99c3449 155 Int_t nClusters = esdTrack->GetNClusters();
156 Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
157 Double_t p = esdTrack->P();
158 Double_t pT = esdTrack->Pt();
159 Double_t eta = esdTrack->Eta();
160 Double_t phi = esdTrack->Phi();
19a4c30d 161 Int_t nChamberHit = 0;
162 for (Int_t ich=0; ich<10; ich++) if (esdTrack->IsInMuonClusterMap(ich)) nChamberHit++;
66085093 163
164 // fill global histograms
165 hNumberOfClusters->Fill(nClusters);
166 hChi2PerDof->Fill(normalizedChi2);
167 hP->Fill(p);
168 hPt->Fill(pT);
169 hEta->Fill(eta);
170 hPhi->Fill(phi);
19a4c30d 171 hNumberOfChamberHit->Fill(nChamberHit);
66085093 172
173 // try to match the reconstructed track with a simulated one
f202486b 174 Int_t nMatchClusters = 0;
175 AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
66085093 176
177 // take actions according to matching result
178 if (matchedTrackRef) {
179
180 // global counter
181 nTotMatchedTracks++;
f202486b 182 if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
23035434 183 trackReconstructedYet = kTRUE;
184 nTotTracksReconstructedYet++;
185 }
66085093 186
187 // fill histograms
f202486b 188 hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
66085093 189 hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
190 hNumberOfClustersM->Fill(nClusters);
191 hChi2PerDofM->Fill(normalizedChi2);
192 hPM->Fill(p);
193 hPtM->Fill(pT);
194 hEtaM->Fill(eta);
195 hPhiM->Fill(phi);
19a4c30d 196 hNumberOfChamberHitM->Fill(nChamberHit);
66085093 197
198 // remove already matched trackRefs
199 trackRefStore->Remove(*matchedTrackRef);
200
201 } else {
202
203 // global counter
204 nTotFakeTracks++;
205
206 // fill histograms
207 hNumberOfClustersF->Fill(nClusters);
208 hChi2PerDofF->Fill(normalizedChi2);
209 hPF->Fill(p);
210 hPtF->Fill(pT);
211 hEtaF->Fill(eta);
212 hPhiF->Fill(phi);
19a4c30d 213 hNumberOfChamberHitF->Fill(nChamberHit);
66085093 214
215 // store fake tracks
a99c3449 216 fakeTrackStore->Add(*muonTrack);
66085093 217
218 }
219
220 } // end of loop over ESD tracks
221
222 // fill histograms
223 hNumberOfTracks->Fill(nTrackerTracks);
23035434 224 nReconstructedTracks += nTrackerTracks;
19a4c30d 225 if (trackReconstructedYet) {
226 if (nEventsWithTrackReconstructedYet >= eventsWithTrackReconstructedYet.GetSize())
227 eventsWithTrackReconstructedYet.Set(10*eventsWithTrackReconstructedYet.GetSize());
228 eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
229 }
66085093 230
231 // count the number the additional fake tracks
232 if (fakeTrackStore->GetSize() > 0) {
233
234 // remove the most connected fake tracks
f202486b 235 Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, useLabel, *hFractionOfConnectedClusters);
66085093 236
237 // remove the remaining free reconstructible tracks
238 Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
239
240 // fill histograms
19a4c30d 241 if (nEventsWithFake >= eventsWithFake.GetSize()) eventsWithFake.Set(10*eventsWithFake.GetSize());
242 eventsWithFake[nEventsWithFake++] = iEvent;
66085093 243 if (nAdditionalTracks > 0) {
19a4c30d 244 if (nEventsWithAdditionalFake >= eventsWithAdditionalFake.GetSize()) eventsWithAdditionalFake.Set(10*eventsWithAdditionalFake.GetSize());
245 eventsWithAdditionalFake[nEventsWithAdditionalFake++] = iEvent;
66085093 246 nTotAdditionalTracks += nAdditionalTracks;
247 hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
248 }
249
250 }
251
252 delete fakeTrackStore;
253
254 } // end of loop over events
23035434 255
256 // total number of connected tracks
257 nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
66085093 258
259 // plot results
260 TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
261 cFakesSummary.Divide(3,2);
262 cFakesSummary.cd(1);
19a4c30d 263 TVirtualPad* pad1 = cFakesSummary.GetPad(1);
264 pad1->Divide(0,2);
265 pad1->cd(1);
266 pad1->GetPad(1)->SetLogy();
66085093 267 hNumberOfClusters->Draw();
268 hNumberOfClusters->SetMinimum(0.5);
269 hNumberOfClustersM->Draw("same");
270 hNumberOfClustersM->SetLineColor(4);
271 hNumberOfClustersF->Draw("same");
272 hNumberOfClustersF->SetLineColor(2);
273 hNumberOfClustersF->SetFillColor(2);
274 hNumberOfClustersF->SetFillStyle(3017);
19a4c30d 275 pad1->cd(2);
276 pad1->GetPad(2)->SetLogy();
277 hNumberOfChamberHit->Draw();
278 hNumberOfChamberHit->SetMinimum(0.5);
279 hNumberOfChamberHitM->Draw("same");
280 hNumberOfChamberHitM->SetLineColor(4);
281 hNumberOfChamberHitF->Draw("same");
282 hNumberOfChamberHitF->SetLineColor(2);
283 hNumberOfChamberHitF->SetFillColor(2);
284 hNumberOfChamberHitF->SetFillStyle(3017);
66085093 285 cFakesSummary.cd(2);
286 cFakesSummary.GetPad(2)->SetLogy();
287 hChi2PerDof->Draw();
288 hChi2PerDof->SetMinimum(0.5);
289 hChi2PerDofM->Draw("same");
290 hChi2PerDofM->SetLineColor(4);
291 hChi2PerDofF->Draw("same");
292 hChi2PerDofF->SetLineColor(2);
293 hChi2PerDofF->SetFillColor(2);
294 hChi2PerDofF->SetFillStyle(3017);
295 cFakesSummary.cd(3);
296 cFakesSummary.GetPad(3)->SetLogy();
297 hP->Draw();
298 hP->SetMinimum(0.5);
299 hPM->Draw("same");
300 hPM->SetLineColor(4);
301 hPF->Draw("same");
302 hPF->SetLineColor(2);
303 hPF->SetFillColor(2);
304 hPF->SetFillStyle(3017);
305 cFakesSummary.cd(4);
306 cFakesSummary.GetPad(4)->SetLogy();
307 hPt->Draw();
308 hPt->SetMinimum(0.5);
309 hPtM->Draw("same");
310 hPtM->SetLineColor(4);
311 hPtF->Draw("same");
312 hPtF->SetLineColor(2);
313 hPtF->SetFillColor(2);
314 hPtF->SetFillStyle(3017);
315 cFakesSummary.cd(5);
316 cFakesSummary.GetPad(5)->SetLogy();
317 hEta->Draw();
318 hEta->SetMinimum(0.5);
319 hEtaM->Draw("same");
320 hEtaM->SetLineColor(4);
321 hEtaF->Draw("same");
322 hEtaF->SetLineColor(2);
323 hEtaF->SetFillColor(2);
324 hEtaF->SetFillStyle(3017);
325 cFakesSummary.cd(6);
326 cFakesSummary.GetPad(6)->SetLogy();
327 hPhi->Draw();
328 hPhi->SetMinimum(0.5);
329 hPhiM->Draw("same");
330 hPhiM->SetLineColor(4);
331 hPhiF->Draw("same");
332 hPhiF->SetLineColor(2);
333 hPhiF->SetFillColor(2);
334 hPhiF->SetFillStyle(3017);
335
336 // save results
337 histoFile->cd();
338 histoFile->Write();
339 cFakesSummary.Write();
340 histoFile->Close();
341
342 // print results
343 cout << endl;
23035434 344 cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
345 cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
66085093 346 cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
23035434 347 cout << " (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
348 if (nTotTracksReconstructedYet > 0) {
349 for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
350 if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
351 else cout << ", " << eventsWithTrackReconstructedYet[i];
352 }
353 cout << "))" << endl;
354 } else cout << ")" << endl;
66085093 355 cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
23035434 356 cout << " (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
357 cout << " (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
358 cout << "- Number of events with fake track(s): " << nEventsWithFake;
359 if (nEventsWithFake > 0) {
360 for(Int_t i=0; i<nEventsWithFake; i++){
361 if (i==0) cout << " (eventID = " << eventsWithFake[i];
362 else cout << ", " << eventsWithFake[i];
363 }
364 cout << ")" << endl;
365 } else cout << endl;
366 cout << " (including " << nEventsWithAdditionalFake << " events with additional track(s)";
367 if (nEventsWithAdditionalFake > 0) {
368 for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
369 if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
370 else cout << ", " << eventsWithAdditionalFake[i];
371 }
372 cout << "))" << endl;
373 } else cout << ")" << endl;
66085093 374 cout << endl;
375 cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
376 cout << endl;
377
66085093 378}
379
380//-----------------------------------------------------------------------
f202486b 381Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
382 Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
66085093 383{
384 /// loop over reconstructible TrackRef not associated with reconstructed track:
385 /// for each of them, find and remove the most connected the fake track, if any,
386 /// and fill the histograms with the fraction of connected clusters.
387 /// Return the number of reconstructible track not connected to any fake
388
389 Int_t nFreeMissingTracks = 0;
390
391 // loop over trackRefs
392 TIter next(trackRefStore.CreateIterator());
393 AliMUONTrack* trackRef;
394 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
395
396 // skip not reconstructible trackRefs
f202486b 397 if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
66085093 398
399 Int_t label = trackRef->GetUniqueID();
400
401 // look for the most connected fake track
402 AliMUONTrack *connectedFake = 0x0;
a99c3449 403 Double_t fractionOfConnectedClusters = 0.;
66085093 404 TIter next2(fakeTrackStore.CreateIterator());
405 AliMUONTrack* fakeTrack;
406 while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
407
408 // get the number of connected clusters
409 Int_t nConnectedClusters = 0;
410 if (useLabel) { // by using the MC label
411 for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
412 if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
413 nConnectedClusters++;
414 } else { // by comparing cluster/TrackRef positions
415 Bool_t compTrack[10];
f202486b 416 nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
66085093 417 }
418
419 // skip non-connected fake tracks
420 if (nConnectedClusters == 0) continue;
421
422 // check if it is the most connected fake track
a99c3449 423 Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
66085093 424 if (f > fractionOfConnectedClusters) {
425 connectedFake = fakeTrack;
426 fractionOfConnectedClusters = f;
427 }
428
429 }
430
431 // remove the most connected fake track
432 if (connectedFake) {
433 hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
434 fakeTrackStore.Remove(*connectedFake);
435 } else nFreeMissingTracks++;
436
437 }
438
439 return nFreeMissingTracks;
440
441}
442