fix finding of pad neighbours; remove methods to write them in OCDB
[u/mrichter/AliRoot.git] / MUON / MUONFakes.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 // ROOT includes
3 #include <TFile.h>
4 #include <TH1.h>
5 #include <TCanvas.h>
6 #include <Riostream.h>
7 #include <TROOT.h>
8 #include <TObjArray.h>
9 #include <TArrayI.h>
10
11 // STEER includes
12 #include "AliLog.h"
13 #include "AliESDEvent.h"
14 #include "AliESDMuonTrack.h"
15 #include "AliCDBManager.h"
16   
17 // MUON includes
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"
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
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);
43
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")
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   
82   // link to reconstructed and simulated tracks
83   AliMUONRecoCheck rc(esdFileName, SimDir);
84   
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;
90   
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();
97   
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);
113   
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++) {
118     
119     // get reconstructed and simulated tracks
120     AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
121     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
122     if (!muonTrackStore || !trackRefStore) continue;
123     
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++;
129     }
130     
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++) {
138       
139       AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
140       
141       // skip ghosts
142       if (!esdTrack->ContainTrackerData()) continue;
143       nTrackerTracks++;
144       
145       // find the corresponding MUON track
146       AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
147       
148       // get track info
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();
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
165       Int_t nMatchClusters = 0;
166       AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
167       
168       // take actions according to matching result
169       if (matchedTrackRef) {
170         
171         // global counter
172         nTotMatchedTracks++;
173         if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
174           trackReconstructedYet = kTRUE;
175           nTotTracksReconstructedYet++;
176         }
177         
178         // fill histograms
179         hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
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
205         fakeTrackStore->Add(*muonTrack);
206         
207       }
208       
209     } // end of loop over ESD tracks
210     
211     // fill histograms
212     hNumberOfTracks->Fill(nTrackerTracks);
213     nReconstructedTracks += nTrackerTracks;
214     if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
215     
216     // count the number the additional fake tracks
217     if (fakeTrackStore->GetSize() > 0) {
218       
219       // remove the most connected fake tracks
220       Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore,  useLabel, *hFractionOfConnectedClusters);
221       
222       // remove the remaining free reconstructible tracks
223       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
224       
225       // fill histograms
226       eventsWithFake[nEventsWithFake] = iEvent;
227       nEventsWithFake++;
228       if (nAdditionalTracks > 0) {
229         eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
230         nEventsWithAdditionalFake++;
231         nTotAdditionalTracks += nAdditionalTracks;
232         hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
233       }
234       
235     }
236     
237     delete fakeTrackStore;
238     
239   } // end of loop over events
240
241   // total number of connected tracks
242   nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
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;
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];
324     }
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];
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;
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   
350 }
351
352 //-----------------------------------------------------------------------
353 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
354                            Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
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
369     if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
370     
371     Int_t label = trackRef->GetUniqueID();
372     
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()) ) ) {
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];
388         nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
389       }
390       
391       // skip non-connected fake tracks
392       if (nConnectedClusters == 0) continue;
393       
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;
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