70ed78f1a8039fc1eba5db33f832473877f17c30
[u/mrichter/AliRoot.git] / MUON / MUONFakes.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 // ROOT includes
3 #include <TTree.h>
4 #include <TFile.h>
5 #include <TH1.h>
6 #include <TCanvas.h>
7 #include <Riostream.h>
8 #include <TROOT.h>
9 #include <TClonesArray.h>
10 #include <TGeoGlobalMagField.h>
11 #include <TArrayI.h>
12
13 // STEER includes
14 #include "AliLog.h"
15 #include "AliMagF.h"
16 #include "AliESDEvent.h"
17 #include "AliESDMuonTrack.h"
18 #include "AliESDMuonCluster.h"
19 #include "AliCDBPath.h"
20 #include "AliCDBEntry.h"
21 #include "AliCDBManager.h"
22   
23 // MUON includes
24 #include "AliMUONTrack.h"
25 #include "AliMUONVTrackStore.h"
26 #include "AliMUONTrackParam.h"
27 #include "AliMUONTrackExtrap.h"
28 #include "AliMUONESDInterface.h"
29 #include "AliMUONRecoCheck.h"
30 #include "AliMUONVCluster.h"
31 #include "AliMUONRecoParam.h"
32 #endif
33
34 /// \ingroup macros
35 /// \file MUONFakes.C
36 ///
37 /// \author Ph. Pillot, Subatech, March. 2009
38 ///
39 /// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
40 /// Results are saved in the root file Fakes.root
41 /// Results are relevent provided that you use the same recoParams as for the reconstruction
42
43 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
44 TTree* GetESDTree(TFile *esdFile);
45 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
46 Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
47 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
48                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
49 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
50                            Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
51
52 //-----------------------------------------------------------------------
53 void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
54                const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
55 {
56   
57   //Reset ROOT and connect tree file
58   gROOT->Reset();
59   
60   // File for histograms and histogram booking
61   TFile *histoFile = new TFile("Fakes.root", "RECREATE");
62   
63   TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
64   TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
65   
66   TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
67   TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
68   TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
69   TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
70   TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
71   TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
72   
73   TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
74   TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
75   TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
76   TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
77   TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
78   TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
79   TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
80   TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
81   TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
82   TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
83   TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
84   TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
85   TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
86   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
87   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
88   
89   // prepare for analysis
90   AliMUONRecoParam *recoParam = 0x0;
91   Double_t sigmaCut = -1;
92   Prepare(recoParam, sigmaCut);
93   
94   // link to reconstructed tracks
95   TFile* esdFile = TFile::Open(esdFileName);
96   TTree* esdTree = GetESDTree(esdFile);
97   AliESDEvent* esd = new AliESDEvent();
98   esd->ReadFromTree(esdTree);
99   
100   // link to simulated tracks
101   AliMUONRecoCheck rc(esdFileName, SimDir);
102   
103   // initialize global counters
104   Int_t nReconstructibleTracks = 0;
105   Int_t nReconstructedTracks = 0;
106   Int_t nEventsWithTrackReconstructedYet = 0;
107   Int_t nEventsWithFake = 0;
108   Int_t nEventsWithAdditionalFake = 0;
109   Int_t nTotMatchedTracks = 0;
110   Int_t nTotTracksReconstructedYet = 0;
111   Int_t nTotFakeTracks = 0;
112   Int_t nTotConnectedTracks = 0;
113   Int_t nTotAdditionalTracks = 0;
114   Bool_t trackReconstructedYet;
115   TArrayI eventsWithTrackReconstructedYet(10);
116   TArrayI eventsWithFake(10);
117   TArrayI eventsWithAdditionalFake(10);
118   
119   // Loop over ESD events
120   FirstEvent = TMath::Max(0, FirstEvent);
121   LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
122   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
123     
124     // get the ESD of current event
125     esdTree->GetEvent(iEvent);
126     if (!esd) {
127       Error("CheckESD", "no ESD object found for event %d", iEvent);
128       return;
129     }
130     
131     // convert TrackRef to MUON tracks
132     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
133     
134     // count the number of reconstructible tracks
135     TIter next(trackRefStore->CreateIterator());
136     AliMUONTrack* trackRef;
137     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
138       if (IsRecontructible(*trackRef,*recoParam)) nReconstructibleTracks++;
139     }
140     
141     // loop over ESD tracks
142     Int_t nTrackerTracks = 0;
143     trackReconstructedYet = kFALSE;
144     AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
145     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
146     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
147       
148       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
149       
150       // skip ghosts
151       if (!muonTrack->ContainTrackerData()) continue;
152       nTrackerTracks++;
153       
154       // get track info
155       Int_t nClusters = muonTrack->GetNClusters();
156       Double_t normalizedChi2 = muonTrack->GetChi2() / (2. * muonTrack->GetNHit() - 5);
157       Double_t p = muonTrack->P();
158       Double_t pT = muonTrack->Pt();
159       Double_t eta = muonTrack->Eta();
160       Double_t phi = muonTrack->Phi();
161       
162       // fill global histograms
163       hNumberOfClusters->Fill(nClusters);
164       hChi2PerDof->Fill(normalizedChi2);
165       hP->Fill(p);
166       hPt->Fill(pT);
167       hEta->Fill(eta);
168       hPhi->Fill(phi);
169       
170       // try to match the reconstructed track with a simulated one
171       Float_t fractionOfMatchCluster = 0.;
172       AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
173       
174       // take actions according to matching result
175       if (matchedTrackRef) {
176         
177         // global counter
178         nTotMatchedTracks++;
179         if (!IsRecontructible(*matchedTrackRef,*recoParam)) {
180           trackReconstructedYet = kTRUE;
181           nTotTracksReconstructedYet++;
182         }
183         
184         // fill histograms
185         hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
186         hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
187         hNumberOfClustersM->Fill(nClusters);
188         hChi2PerDofM->Fill(normalizedChi2);
189         hPM->Fill(p);
190         hPtM->Fill(pT);
191         hEtaM->Fill(eta);
192         hPhiM->Fill(phi);
193         
194         // remove already matched trackRefs
195         trackRefStore->Remove(*matchedTrackRef);
196         
197       } else {
198         
199         // global counter
200         nTotFakeTracks++;
201         
202         // fill histograms
203         hNumberOfClustersF->Fill(nClusters);
204         hChi2PerDofF->Fill(normalizedChi2);
205         hPF->Fill(p);
206         hPtF->Fill(pT);
207         hEtaF->Fill(eta);
208         hPhiF->Fill(phi);
209         
210         // store fake tracks
211         AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
212         
213       }
214       
215     } // end of loop over ESD tracks
216     
217     // fill histograms
218     hNumberOfTracks->Fill(nTrackerTracks);
219     nReconstructedTracks += nTrackerTracks;
220     if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
221     
222     // count the number the additional fake tracks
223     if (fakeTrackStore->GetSize() > 0) {
224       
225       // remove the most connected fake tracks
226       Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
227                                                       useLabel, sigmaCut, *hFractionOfConnectedClusters);
228       
229       // remove the remaining free reconstructible tracks
230       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
231       
232       // fill histograms
233       eventsWithFake[nEventsWithFake] = iEvent;
234       nEventsWithFake++;
235       if (nAdditionalTracks > 0) {
236         eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
237         nEventsWithAdditionalFake++;
238         nTotAdditionalTracks += nAdditionalTracks;
239         hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
240       }
241       
242     }
243     
244     delete fakeTrackStore;
245     
246   } // end of loop over events
247
248   // total number of connected tracks
249   nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
250   
251   // plot results
252   TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
253   cFakesSummary.Divide(3,2);
254   cFakesSummary.cd(1);
255   cFakesSummary.GetPad(1)->SetLogy();
256   hNumberOfClusters->Draw();
257   hNumberOfClusters->SetMinimum(0.5);
258   hNumberOfClustersM->Draw("same");
259   hNumberOfClustersM->SetLineColor(4);
260   hNumberOfClustersF->Draw("same");
261   hNumberOfClustersF->SetLineColor(2);
262   hNumberOfClustersF->SetFillColor(2);
263   hNumberOfClustersF->SetFillStyle(3017);
264   cFakesSummary.cd(2);
265   cFakesSummary.GetPad(2)->SetLogy();
266   hChi2PerDof->Draw();
267   hChi2PerDof->SetMinimum(0.5);
268   hChi2PerDofM->Draw("same");
269   hChi2PerDofM->SetLineColor(4);
270   hChi2PerDofF->Draw("same");
271   hChi2PerDofF->SetLineColor(2);
272   hChi2PerDofF->SetFillColor(2);
273   hChi2PerDofF->SetFillStyle(3017);
274   cFakesSummary.cd(3);
275   cFakesSummary.GetPad(3)->SetLogy();
276   hP->Draw();
277   hP->SetMinimum(0.5);
278   hPM->Draw("same");
279   hPM->SetLineColor(4);
280   hPF->Draw("same");
281   hPF->SetLineColor(2);
282   hPF->SetFillColor(2);
283   hPF->SetFillStyle(3017);
284   cFakesSummary.cd(4);
285   cFakesSummary.GetPad(4)->SetLogy();
286   hPt->Draw();
287   hPt->SetMinimum(0.5);
288   hPtM->Draw("same");
289   hPtM->SetLineColor(4);
290   hPtF->Draw("same");
291   hPtF->SetLineColor(2);
292   hPtF->SetFillColor(2);
293   hPtF->SetFillStyle(3017);
294   cFakesSummary.cd(5);
295   cFakesSummary.GetPad(5)->SetLogy();
296   hEta->Draw();
297   hEta->SetMinimum(0.5);
298   hEtaM->Draw("same");
299   hEtaM->SetLineColor(4);
300   hEtaF->Draw("same");
301   hEtaF->SetLineColor(2);
302   hEtaF->SetFillColor(2);
303   hEtaF->SetFillStyle(3017);
304   cFakesSummary.cd(6);
305   cFakesSummary.GetPad(6)->SetLogy();
306   hPhi->Draw();
307   hPhi->SetMinimum(0.5);
308   hPhiM->Draw("same");
309   hPhiM->SetLineColor(4);
310   hPhiF->Draw("same");
311   hPhiF->SetLineColor(2);
312   hPhiF->SetFillColor(2);
313   hPhiF->SetFillStyle(3017);
314   
315   // save results
316   histoFile->cd();
317   histoFile->Write();
318   cFakesSummary.Write();
319   histoFile->Close();
320   
321   // print results
322   cout << endl;
323   cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
324   cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
325   cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
326   cout << "  (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
327   if (nTotTracksReconstructedYet > 0) {
328     for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
329       if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
330       else cout << ", " << eventsWithTrackReconstructedYet[i];
331     }
332     cout << "))" << endl;
333   } else cout << ")" << endl;
334   cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
335   cout << "  (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
336   cout << "  (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
337   cout << "- Number of events with fake track(s): " << nEventsWithFake;
338   if (nEventsWithFake > 0) {
339     for(Int_t i=0; i<nEventsWithFake; i++){
340       if (i==0) cout << " (eventID = " << eventsWithFake[i];
341       else cout << ", " << eventsWithFake[i];
342     }
343     cout << ")" << endl;
344   } else cout << endl;
345   cout << "  (including " << nEventsWithAdditionalFake << " events with additional track(s)";
346   if (nEventsWithAdditionalFake > 0) {
347     for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
348       if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
349       else cout << ", " << eventsWithAdditionalFake[i];
350     }
351     cout << "))" << endl;
352   } else cout << ")" << endl;
353   cout << endl;
354   cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
355   cout << endl;
356   
357   // clear memory
358   delete esd;
359   esdFile->Close();
360   delete recoParam;
361   
362 }
363
364 //-----------------------------------------------------------------------
365 void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
366 {
367   /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
368   
369   // prepare OCDB access
370   AliCDBManager* man = AliCDBManager::Instance();
371   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
372   man->SetRun(0);
373   
374   // set  mag field 
375   // waiting for mag field in CDB 
376   if (!TGeoGlobalMagField::Instance()->GetField()) {
377     printf("Loading field map...\n");
378     AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
379     TGeoGlobalMagField::Instance()->SetField(field);
380   }
381   // set the magnetic field for track extrapolations
382   AliMUONTrackExtrap::SetField();
383   
384   // Load initial reconstruction parameters from OCDB
385   AliCDBPath path("MUON","Calib","RecoParam");
386   AliCDBEntry *entry=man->Get(path.GetPath());
387   if(entry) {
388     recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
389     entry->SetOwner(0);
390     AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
391   }
392   if (!recoParam) {
393     printf("Couldn't find RecoParam object in OCDB: create default one");
394     recoParam = AliMUONRecoParam::GetLowFluxParam();
395   }
396   
397   Info("MUONFakes", "\n recontruction parameters:");
398   recoParam->Print("FULL");
399   AliMUONESDInterface::ResetTracker(recoParam);
400   
401   // sigma cut to associate clusters with TrackRefs in case the label are not used
402   sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); 
403   
404 }
405
406 //-----------------------------------------------------------------------
407 TTree* GetESDTree(TFile *esdFile)
408 {
409   /// Check that the file is properly open
410   /// Return pointer to the ESD Tree
411   
412   if (!esdFile || !esdFile->IsOpen()) {
413     Error("GetESDTree", "opening ESD file failed");
414     exit(-1);
415   }
416   
417   TTree* tree = (TTree*) esdFile->Get("esdTree");
418   if (!tree) {
419     Error("GetESDTree", "no ESD tree found");
420     exit(-1);
421   }
422   
423   return tree;
424   
425 }
426
427 //-----------------------------------------------------------------------
428 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
429 {
430   /// Try to match 2 tracks
431   
432   Bool_t compTrack[10];
433   Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
434   fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
435   
436   if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
437       (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
438       fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
439   else return kFALSE;
440   
441 }
442
443 //-----------------------------------------------------------------------
444 Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
445 {
446   /// Check il the track is reconstructible
447   Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3;
448   Int_t currentCh, previousCh = -1, nChHitInSt45 = 0;
449   Bool_t clusterInSt[5];
450   for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt);
451   
452   AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
453   while (trackParam) {
454     
455     currentCh = trackParam->GetClusterPtr()->GetChamberId();
456     
457     clusterInSt[currentCh/2] = kTRUE;
458     
459     if (currentCh > 5 && currentCh != previousCh) {
460       nChHitInSt45++;
461       previousCh = currentCh;
462     }
463     
464     trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
465   }
466   
467   return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
468           clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
469   
470 }
471
472 //-----------------------------------------------------------------------
473 AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
474                                 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
475 {
476   /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
477   
478   AliMUONTrack *matchedTrackRef = 0x0;
479   fractionOfMatchCluster = 0.;
480   
481   if (useLabel) { // by using the MC label
482     
483     // get the corresponding simulated track if any
484     Int_t label = muonTrack.GetLabel();
485     matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
486     
487     // get the fraction of matched clusters
488     if (matchedTrackRef) {
489       Int_t nMatchClusters = 0;
490       if (muonTrack.ClustersStored()) {
491         AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
492         while (cluster) {
493           if (cluster->GetLabel() == label) nMatchClusters++;
494           cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
495         }
496       }
497       fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
498     }
499     
500   } else { // by comparing cluster/TrackRef positions
501     
502     // convert ESD track to MUON track
503     AliMUONTrack track;
504     AliMUONESDInterface::ESDToMUON(muonTrack,track);
505     
506     // look for the corresponding simulated track if any
507     TIter next(trackRefStore.CreateIterator());
508     AliMUONTrack* trackRef;
509     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
510       
511       // check compatibility
512       Float_t f = 0.;
513       if (TrackMatched(track, *trackRef, f, sigmaCut)) {
514         matchedTrackRef = trackRef;
515         fractionOfMatchCluster = f;
516         break;
517       }
518       
519     }
520     
521   }
522   
523   return matchedTrackRef;
524   
525 }
526
527 //-----------------------------------------------------------------------
528 Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
529                           Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
530 {
531   /// loop over reconstructible TrackRef not associated with reconstructed track:
532   /// for each of them, find and remove the most connected the fake track, if any,
533   /// and fill the histograms with the fraction of connected clusters.
534   /// Return the number of reconstructible track not connected to any fake
535   
536   Int_t nFreeMissingTracks = 0;
537   
538   // loop over trackRefs
539   TIter next(trackRefStore.CreateIterator());
540   AliMUONTrack* trackRef;
541   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
542     
543     // skip not reconstructible trackRefs
544     if (!IsRecontructible(*trackRef,recoParam)) continue;
545     
546     Int_t label = trackRef->GetUniqueID();
547     
548     // look for the most connected fake track
549     AliMUONTrack *connectedFake = 0x0;
550     Float_t fractionOfConnectedClusters = 0.;
551     TIter next2(fakeTrackStore.CreateIterator());
552     AliMUONTrack* fakeTrack;
553     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
554       
555       // get the number of connected clusters
556       Int_t nConnectedClusters = 0;
557       if (useLabel) { // by using the MC label
558         for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
559           if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
560             nConnectedClusters++;
561       } else { // by comparing cluster/TrackRef positions
562         Bool_t compTrack[10];
563         nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack);
564       }
565       
566       // skip non-connected fake tracks
567       if (nConnectedClusters == 0) continue;
568       
569       // check if it is the most connected fake track
570       Float_t f = ((Float_t)nConnectedClusters) / ((Float_t)fakeTrack->GetNClusters());
571       if (f > fractionOfConnectedClusters) {
572         connectedFake = fakeTrack;
573         fractionOfConnectedClusters = f;
574       }
575       
576     }
577     
578     // remove the most connected fake track
579     if (connectedFake) {
580       hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
581       fakeTrackStore.Remove(*connectedFake);
582     } else nFreeMissingTracks++;
583     
584   }
585   
586   return nFreeMissingTracks;
587   
588 }
589