]>
Commit | Line | Data |
---|---|---|
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 | 37 | UInt_t requestedStationMask = 0; |
38 | Bool_t request2ChInSameSt45 = kFALSE; | |
39 | Double_t sigmaCut = -1.; | |
40 | ||
41 | Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, | |
42 | Bool_t useLabel, TH1F &hFractionOfConnectedClusters); | |
66085093 | 43 | |
44 | //----------------------------------------------------------------------- | |
45 | void 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 | ||
50 | //Reset ROOT and connect tree file | |
51 | gROOT->Reset(); | |
52 | ||
53 | // File for histograms and histogram booking | |
54 | TFile *histoFile = new TFile("Fakes.root", "RECREATE"); | |
55 | ||
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.); | |
58 | ||
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); | |
65 | ||
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.); | |
81 | ||
a99c3449 | 82 | // link to reconstructed and simulated tracks |
83 | AliMUONRecoCheck rc(esdFileName, SimDir); | |
66085093 | 84 | |
a99c3449 | 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; | |
66085093 | 90 | |
a99c3449 | 91 | // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used |
f202486b | 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(); | |
66085093 | 97 | |
98 | // initialize global counters | |
23035434 | 99 | Int_t nReconstructibleTracks = 0; |
100 | Int_t nReconstructedTracks = 0; | |
101 | Int_t nEventsWithTrackReconstructedYet = 0; | |
66085093 | 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; | |
23035434 | 107 | Int_t nTotConnectedTracks = 0; |
66085093 | 108 | Int_t nTotAdditionalTracks = 0; |
23035434 | 109 | Bool_t trackReconstructedYet; |
110 | TArrayI eventsWithTrackReconstructedYet(10); | |
111 | TArrayI eventsWithFake(10); | |
112 | TArrayI eventsWithAdditionalFake(10); | |
66085093 | 113 | |
114 | // Loop over ESD events | |
115 | FirstEvent = TMath::Max(0, FirstEvent); | |
a99c3449 | 116 | LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1; |
66085093 | 117 | for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) { |
118 | ||
a99c3449 | 119 | // get reconstructed and simulated tracks |
120 | AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE); | |
66085093 | 121 | AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent); |
a99c3449 | 122 | if (!muonTrackStore || !trackRefStore) continue; |
66085093 | 123 | |
23035434 | 124 | // count the number of reconstructible tracks |
125 | TIter next(trackRefStore->CreateIterator()); | |
126 | AliMUONTrack* trackRef; | |
127 | while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) { | |
f202486b | 128 | if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++; |
23035434 | 129 | } |
130 | ||
66085093 | 131 | // loop over ESD tracks |
132 | Int_t nTrackerTracks = 0; | |
23035434 | 133 | trackReconstructedYet = kFALSE; |
66085093 | 134 | AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore(); |
a99c3449 | 135 | const AliESDEvent* esd = rc.GetESDEvent(); |
66085093 | 136 | Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; |
137 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
138 | ||
a99c3449 | 139 | AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack); |
66085093 | 140 | |
141 | // skip ghosts | |
a99c3449 | 142 | if (!esdTrack->ContainTrackerData()) continue; |
66085093 | 143 | nTrackerTracks++; |
144 | ||
a99c3449 | 145 | // find the corresponding MUON track |
146 | AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID()); | |
147 | ||
66085093 | 148 | // get track info |
a99c3449 | 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(); | |
66085093 | 155 | |
156 | // fill global histograms | |
157 | hNumberOfClusters->Fill(nClusters); | |
158 | hChi2PerDof->Fill(normalizedChi2); | |
159 | hP->Fill(p); | |
160 | hPt->Fill(pT); | |
161 | hEta->Fill(eta); | |
162 | hPhi->Fill(phi); | |
163 | ||
164 | // try to match the reconstructed track with a simulated one | |
f202486b | 165 | Int_t nMatchClusters = 0; |
166 | AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut); | |
66085093 | 167 | |
168 | // take actions according to matching result | |
169 | if (matchedTrackRef) { | |
170 | ||
171 | // global counter | |
172 | nTotMatchedTracks++; | |
f202486b | 173 | if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) { |
23035434 | 174 | trackReconstructedYet = kTRUE; |
175 | nTotTracksReconstructedYet++; | |
176 | } | |
66085093 | 177 | |
178 | // fill histograms | |
f202486b | 179 | hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters)); |
66085093 | 180 | hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters()); |
181 | hNumberOfClustersM->Fill(nClusters); | |
182 | hChi2PerDofM->Fill(normalizedChi2); | |
183 | hPM->Fill(p); | |
184 | hPtM->Fill(pT); | |
185 | hEtaM->Fill(eta); | |
186 | hPhiM->Fill(phi); | |
187 | ||
188 | // remove already matched trackRefs | |
189 | trackRefStore->Remove(*matchedTrackRef); | |
190 | ||
191 | } else { | |
192 | ||
193 | // global counter | |
194 | nTotFakeTracks++; | |
195 | ||
196 | // fill histograms | |
197 | hNumberOfClustersF->Fill(nClusters); | |
198 | hChi2PerDofF->Fill(normalizedChi2); | |
199 | hPF->Fill(p); | |
200 | hPtF->Fill(pT); | |
201 | hEtaF->Fill(eta); | |
202 | hPhiF->Fill(phi); | |
203 | ||
204 | // store fake tracks | |
a99c3449 | 205 | fakeTrackStore->Add(*muonTrack); |
66085093 | 206 | |
207 | } | |
208 | ||
209 | } // end of loop over ESD tracks | |
210 | ||
211 | // fill histograms | |
212 | hNumberOfTracks->Fill(nTrackerTracks); | |
23035434 | 213 | nReconstructedTracks += nTrackerTracks; |
214 | if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent; | |
66085093 | 215 | |
216 | // count the number the additional fake tracks | |
217 | if (fakeTrackStore->GetSize() > 0) { | |
218 | ||
219 | // remove the most connected fake tracks | |
f202486b | 220 | Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, useLabel, *hFractionOfConnectedClusters); |
66085093 | 221 | |
222 | // remove the remaining free reconstructible tracks | |
223 | Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks; | |
224 | ||
225 | // fill histograms | |
23035434 | 226 | eventsWithFake[nEventsWithFake] = iEvent; |
66085093 | 227 | nEventsWithFake++; |
228 | if (nAdditionalTracks > 0) { | |
23035434 | 229 | eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent; |
66085093 | 230 | nEventsWithAdditionalFake++; |
231 | nTotAdditionalTracks += nAdditionalTracks; | |
232 | hNumberOfAdditionalTracks->Fill(nAdditionalTracks); | |
233 | } | |
234 | ||
235 | } | |
236 | ||
237 | delete fakeTrackStore; | |
238 | ||
239 | } // end of loop over events | |
23035434 | 240 | |
241 | // total number of connected tracks | |
242 | nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries(); | |
66085093 | 243 | |
244 | // plot results | |
245 | TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600); | |
246 | cFakesSummary.Divide(3,2); | |
247 | cFakesSummary.cd(1); | |
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); | |
257 | cFakesSummary.cd(2); | |
258 | cFakesSummary.GetPad(2)->SetLogy(); | |
259 | hChi2PerDof->Draw(); | |
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); | |
267 | cFakesSummary.cd(3); | |
268 | cFakesSummary.GetPad(3)->SetLogy(); | |
269 | hP->Draw(); | |
270 | hP->SetMinimum(0.5); | |
271 | hPM->Draw("same"); | |
272 | hPM->SetLineColor(4); | |
273 | hPF->Draw("same"); | |
274 | hPF->SetLineColor(2); | |
275 | hPF->SetFillColor(2); | |
276 | hPF->SetFillStyle(3017); | |
277 | cFakesSummary.cd(4); | |
278 | cFakesSummary.GetPad(4)->SetLogy(); | |
279 | hPt->Draw(); | |
280 | hPt->SetMinimum(0.5); | |
281 | hPtM->Draw("same"); | |
282 | hPtM->SetLineColor(4); | |
283 | hPtF->Draw("same"); | |
284 | hPtF->SetLineColor(2); | |
285 | hPtF->SetFillColor(2); | |
286 | hPtF->SetFillStyle(3017); | |
287 | cFakesSummary.cd(5); | |
288 | cFakesSummary.GetPad(5)->SetLogy(); | |
289 | hEta->Draw(); | |
290 | hEta->SetMinimum(0.5); | |
291 | hEtaM->Draw("same"); | |
292 | hEtaM->SetLineColor(4); | |
293 | hEtaF->Draw("same"); | |
294 | hEtaF->SetLineColor(2); | |
295 | hEtaF->SetFillColor(2); | |
296 | hEtaF->SetFillStyle(3017); | |
297 | cFakesSummary.cd(6); | |
298 | cFakesSummary.GetPad(6)->SetLogy(); | |
299 | hPhi->Draw(); | |
300 | hPhi->SetMinimum(0.5); | |
301 | hPhiM->Draw("same"); | |
302 | hPhiM->SetLineColor(4); | |
303 | hPhiF->Draw("same"); | |
304 | hPhiF->SetLineColor(2); | |
305 | hPhiF->SetFillColor(2); | |
306 | hPhiF->SetFillStyle(3017); | |
307 | ||
308 | // save results | |
309 | histoFile->cd(); | |
310 | histoFile->Write(); | |
311 | cFakesSummary.Write(); | |
312 | histoFile->Close(); | |
313 | ||
314 | // print results | |
315 | cout << endl; | |
23035434 | 316 | cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl; |
317 | cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl; | |
66085093 | 318 | cout << "- Number of matched tracks: " << nTotMatchedTracks << endl; |
23035434 | 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]; | |
324 | } | |
325 | cout << "))" << endl; | |
326 | } else cout << ")" << endl; | |
66085093 | 327 | cout << "- Number of fake tracks: " << nTotFakeTracks << endl; |
23035434 | 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]; | |
335 | } | |
336 | cout << ")" << endl; | |
337 | } else cout << endl; | |
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]; | |
343 | } | |
344 | cout << "))" << endl; | |
345 | } else cout << ")" << endl; | |
66085093 | 346 | cout << endl; |
347 | cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl; | |
348 | cout << endl; | |
349 | ||
66085093 | 350 | } |
351 | ||
352 | //----------------------------------------------------------------------- | |
f202486b | 353 | Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, |
354 | Bool_t useLabel, TH1F &hFractionOfConnectedClusters) | |
66085093 | 355 | { |
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 | |
360 | ||
361 | Int_t nFreeMissingTracks = 0; | |
362 | ||
363 | // loop over trackRefs | |
364 | TIter next(trackRefStore.CreateIterator()); | |
365 | AliMUONTrack* trackRef; | |
366 | while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) { | |
367 | ||
368 | // skip not reconstructible trackRefs | |
f202486b | 369 | if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue; |
66085093 | 370 | |
371 | Int_t label = trackRef->GetUniqueID(); | |
372 | ||
373 | // look for the most connected fake track | |
374 | AliMUONTrack *connectedFake = 0x0; | |
a99c3449 | 375 | Double_t fractionOfConnectedClusters = 0.; |
66085093 | 376 | TIter next2(fakeTrackStore.CreateIterator()); |
377 | AliMUONTrack* fakeTrack; | |
378 | while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) { | |
379 | ||
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]; | |
f202486b | 388 | nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack); |
66085093 | 389 | } |
390 | ||
391 | // skip non-connected fake tracks | |
392 | if (nConnectedClusters == 0) continue; | |
393 | ||
394 | // check if it is the most connected fake track | |
a99c3449 | 395 | Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters()); |
66085093 | 396 | if (f > fractionOfConnectedClusters) { |
397 | connectedFake = fakeTrack; | |
398 | fractionOfConnectedClusters = f; | |
399 | } | |
400 | ||
401 | } | |
402 | ||
403 | // remove the most connected fake track | |
404 | if (connectedFake) { | |
405 | hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters); | |
406 | fakeTrackStore.Remove(*connectedFake); | |
407 | } else nFreeMissingTracks++; | |
408 | ||
409 | } | |
410 | ||
411 | return nFreeMissingTracks; | |
412 | ||
413 | } | |
414 |