Identify decay chains + new opt. (track selec., comb. label/pos. match, ...): added...
[u/mrichter/AliRoot.git] / PWGPP / MUON / dep / AliAnalysisTaskMuonFakes.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // ROOT includes
19 #include <TObjArray.h>
20 #include <TClonesArray.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TCanvas.h>
24 #include <TArrayI.h>
25
26 // STEER includes
27 #include "AliLog.h"
28 #include "AliESDEvent.h"
29 #include "AliESDMuonTrack.h"
30 #include "AliESDInputHandler.h"
31 #include "AliMCEventHandler.h"
32 #include "AliCDBManager.h"
33 #include "AliCentrality.h"
34 #include "AliMCParticle.h"
35 #include "AliMCEvent.h"
36 #include "AliCounterCollection.h"
37
38 // ANALYSIS includes
39 #include "AliAnalysisDataSlot.h"
40 #include "AliAnalysisDataContainer.h"
41 #include "AliAnalysisManager.h"
42
43 // MUON includes
44 #include "AliMUONCDB.h"
45 #include "AliMUONRecoParam.h"
46 #include "AliMUONRecoCheck.h"
47 #include "AliMUONVCluster.h"
48 #include "AliMUONVTrackStore.h"
49 #include "AliMUONVTriggerTrackStore.h"
50 #include "AliMUONTrack.h"
51 #include "AliMUONTrackParam.h"
52 #include "AliMUONESDInterface.h"
53 #include "AliMUONTriggerTrack.h"
54
55 #include "AliAnalysisTaskMuonFakes.h"
56
57 ClassImp(AliAnalysisTaskMuonFakes)
58
59 //________________________________________________________________________
60 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes() :
61 AliAnalysisTaskSE(), 
62 fList(0x0),
63 fList2(0x0),
64 fCanvases(0x0),
65 fTrackCounters(0x0),
66 fFakeTrackCounters(0x0),
67 fMatchedTrackCounters(0x0),
68 fEventCounters(0x0),
69 fPairCounters(0x0),
70 fCurrentFileName(""),
71 fRequestedStationMask(0),
72 fRequest2ChInSameSt45(kFALSE),
73 fSigmaCut(-1.),
74 fNEvents(0),
75 fShowProgressBar(kFALSE),
76 fUseLabel(kFALSE),
77 fCombineMCId(kFALSE),
78 fExternalSigmaCut(-1.),
79 fMatchTrig(kFALSE),
80 fApplyAccCut(kFALSE),
81 fChi2Cut(-1.),
82 fPtCut(-1.),
83 fRecoParamLocation(""),
84 fDecayAsFake(kFALSE),
85 fPrintDecayChain(kFALSE),
86 fDisableDetailedCounters(kFALSE),
87 fMuonTrackCuts(0x0)
88 {
89   /// Default constructor.
90 }
91
92 //________________________________________________________________________
93 AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes(const char *name) :
94 AliAnalysisTaskSE(name), 
95 fList(0x0),
96 fList2(0x0),
97 fCanvases(0x0),
98 fTrackCounters(0x0),
99 fFakeTrackCounters(0x0),
100 fMatchedTrackCounters(0x0),
101 fEventCounters(0x0),
102 fPairCounters(0x0),
103 fCurrentFileName(""),
104 fRequestedStationMask(0),
105 fRequest2ChInSameSt45(kFALSE),
106 fSigmaCut(-1.),
107 fNEvents(0),
108 fShowProgressBar(kFALSE),
109 fUseLabel(kFALSE),
110 fCombineMCId(kFALSE),
111 fExternalSigmaCut(-1.),
112 fMatchTrig(kFALSE),
113 fApplyAccCut(kFALSE),
114 fChi2Cut(-1.),
115 fPtCut(-1.),
116 fRecoParamLocation("raw://"),
117 fDecayAsFake(kFALSE),
118 fPrintDecayChain(kFALSE),
119 fDisableDetailedCounters(kFALSE),
120 fMuonTrackCuts(0x0)
121 {
122   /// Constructor.
123   // Output slot #1 writes into a TObjArray container
124   DefineOutput(1,TObjArray::Class());
125   // Output slot #2 writes into an AliCounterCollection container
126   DefineOutput(2,AliCounterCollection::Class());
127   // Output slot #3 writes into an AliCounterCollection container
128   DefineOutput(3,AliCounterCollection::Class());
129   // Output slot #4 writes into an AliCounterCollection container
130   DefineOutput(4,AliCounterCollection::Class());
131   // Output slot #5 writes into an AliCounterCollection container
132   DefineOutput(5,AliCounterCollection::Class());
133   // Output slot #6 writes into a TObjArray container
134   DefineOutput(6,TObjArray::Class());
135   // Output slot #7 writes into an AliCounterCollection container
136   DefineOutput(7,AliCounterCollection::Class());
137 }
138
139 //________________________________________________________________________
140 AliAnalysisTaskMuonFakes::~AliAnalysisTaskMuonFakes()
141 {
142   /// Destructor.
143   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
144     delete fList;
145     delete fList2;
146     delete fTrackCounters;
147     delete fFakeTrackCounters;
148     delete fMatchedTrackCounters;
149     delete fEventCounters;
150     delete fPairCounters;
151   }
152   delete fCanvases;
153   delete fMuonTrackCuts;
154 }
155
156 //___________________________________________________________________________
157 void AliAnalysisTaskMuonFakes::UserCreateOutputObjects()
158 {
159   /// Create histograms and counters.
160   
161   // single track histograms
162   fList = new TObjArray(100);
163   fList->SetOwner();
164   
165   TH1F *h = 0x0;
166   TH2F *h2 = 0x0;
167   TString nameSuffix0[2] = {"", "S"};
168   TString nameSuffixT[6] = {"", "M", "MY", "D", "DY", "F"};
169   TString titlePrefix0[2] = {"", "selected "};
170   TString titlePrefixT[6] = {"", "matched ", "not reconstructible matched ", "decay ", "not reconstructible decay ", "fake "};
171   for (Int_t i = 0; i < 2; i++) {
172     
173     for (Int_t j = 0; j < 6; j++) {
174       
175       // number of clusters
176       h = new TH1F(Form("hNumberOfClusters%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
177                    Form("nb of clusters /%s%strack",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 21, -0.5, 20.5);
178       fList->AddAtAndExpand(h, kNumberOfClusters+i*kNhistTrack+j);
179       
180       // number of fired chambers
181       h = new TH1F(Form("hNumberOfChamberHit%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
182                    Form("nb of chambers hit /%s%strack",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 16, -0.5, 15.5);
183       fList->AddAtAndExpand(h, kNumberOfChamberHit+i*kNhistTrack+j);
184       
185       // chi2
186       h = new TH1F(Form("hChi2PerDof%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
187                    Form("%s%strack chi2/d.o.f.",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 200, 0., 20.);
188       fList->AddAtAndExpand(h, kChi2PerDof+i*kNhistTrack+j);
189       
190       // chi2 versus number of clusters
191       h2 = new TH2F(Form("hChi2PerDofVsNClusters%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
192                     Form("%s%strack chi2/d.o.f. versus nb of clusters",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 21, -0.5, 20.5, 100, 0., 20.);
193       fList->AddAtAndExpand(h2, kChi2PerDofVsNClusters+i*kNhistTrack+j);
194       
195       // chi2 versus number of fired chambers
196       h2 = new TH2F(Form("hChi2PerDofVsNChamberHit%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
197                     Form("%s%strack chi2/d.o.f. versus nb of fired chambers",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 16, -0.5, 15.5, 100, 0., 20.);
198       fList->AddAtAndExpand(h2, kChi2PerDofVsNChamberHit+i*kNhistTrack+j);
199       
200       // physics quantities
201       h = new TH1F(Form("hP%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
202                    Form("%s%strack P distribution (GeV/c)",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 100, 0., 200.);
203       fList->AddAtAndExpand(h, kP+i*kNhistTrack+j);
204       h = new TH1F(Form("hPt%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
205                    Form("%s%strack Pt distribution (GeV/c)",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 100, 0., 20.);
206       fList->AddAtAndExpand(h, kPt+i*kNhistTrack+j);
207       h = new TH1F(Form("hEta%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
208                    Form("%s%strack pseudo-rapidity distribution",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 200, -10., 0.);
209       fList->AddAtAndExpand(h , kEta+i*kNhistTrack+j);
210       h = new TH1F(Form("hPhi%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
211                    Form("%s%strack phi distribution",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 100, -1., 9.);
212       fList->AddAtAndExpand(h, kPhi+i*kNhistTrack+j);
213       h = new TH1F(Form("hDCA%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
214                    Form("%s%strack DCA distribution",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 250, 0., 500.);
215       fList->AddAtAndExpand(h, kDCA+i*kNhistTrack+j);
216       h = new TH1F(Form("hPDCA23%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
217                    Form("%s%strack P*DCA distribution in 2-3 deg",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 250, 0., 5000.);
218       fList->AddAtAndExpand(h, kPDCA23+i*kNhistTrack+j);
219       h = new TH1F(Form("hPDCA310%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
220                    Form("%s%strack P*DCA distribution in 3-10 deg",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 250, 0., 5000.);
221       fList->AddAtAndExpand(h, kPDCA310+i*kNhistTrack+j);
222       h = new TH1F(Form("hRAbs%s%s",nameSuffix0[i].Data(),nameSuffixT[j].Data()),
223                    Form("%s%strack R_{Abs} distribution",titlePrefix0[i].Data(),titlePrefixT[j].Data()), 200, 0., 100.);
224       fList->AddAtAndExpand(h, kRAbs+i*kNhistTrack+j);
225       
226     }
227     
228   }
229   
230   // number of tracks
231   TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks", "nb of tracks /evt", 21, -0.5, 20.5);
232   fList->AddAtAndExpand(hNumberOfTracks, 2*kNhistTrack+kNumberOfTracks);
233   TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks", "nb of fake - nb of missing track", 21, -0.5, 20.5);
234   fList->AddAtAndExpand(hNumberOfAdditionalTracks, 2*kNhistTrack+kNumberOfAdditionalTracks);
235   
236   // number of clusters MC / fraction of clusters
237   TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC", "nb of clusters /MC track", 21, -0.5, 20.5);
238   fList->AddAtAndExpand(hNumberOfClustersMC, 2*kNhistTrack+kNumberOfClustersMC);
239   TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters", "nb of matched clusters / nb of clusters", 110, 0., 1.1);
240   fList->AddAtAndExpand(hFractionOfMatchedClusters, 2*kNhistTrack+kFractionOfMatchedClusters);
241   TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters", "nb of connected clusters / nb of clusters in fake tracks", 110, 0., 1.1);
242   fList->AddAtAndExpand(hFractionOfConnectedClusters, 2*kNhistTrack+kFractionOfConnectedClusters);
243   
244   // track pair histograms
245   fList2 = new TObjArray(100);
246   fList2->SetOwner();
247   
248   // physics quantities of opposite-sign track pairs
249   TString nameSuffix[4] = {"", "M", "F1", "F2"};
250   TString titlePrefix[4] = {"dimuon ", "matched-matched pair ", "matched-fake pair ", "fake-fake pair "};
251   for (Int_t i = 0; i < 2; i++) {
252     for (Int_t j = 0; j < 4; j++) {
253       h = new TH1F(Form("h2Mass%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
254                    Form("%s%smass distribution (GeV/c^{2})",titlePrefix0[i].Data(),titlePrefix[j].Data()), 300, 0., 15.);
255       fList2->AddAtAndExpand(h, k2Mass+i*kNhistPair+j);
256       h = new TH1F(Form("h2P%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
257                    Form("%s%sP distribution (GeV/c)",titlePrefix0[i].Data(),titlePrefix[j].Data()), 100, 0., 200.);
258       fList2->AddAtAndExpand(h, k2P+i*kNhistPair+j);
259       h = new TH1F(Form("h2Pt%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
260                    Form("%s%sPt distribution (GeV/c)",titlePrefix0[i].Data(),titlePrefix[j].Data()), 100, 0., 20.);
261       fList2->AddAtAndExpand(h, k2Pt+i*kNhistPair+j);
262       h = new TH1F(Form("h2Y%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
263                    Form("%s%srapidity distribution",titlePrefix0[i].Data(),titlePrefix[j].Data()), 200, -10., 0.);
264       fList2->AddAtAndExpand(h, k2Y+i*kNhistPair+j);
265       h = new TH1F(Form("h2Eta%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
266                    Form("%s%spseudo-rapidity distribution",titlePrefix0[i].Data(),titlePrefix[j].Data()), 200, -10., 0.);
267       fList2->AddAtAndExpand(h, k2Eta+i*kNhistPair+j);
268       h = new TH1F(Form("h2Phi%s%s",nameSuffix0[i].Data(),nameSuffix[j].Data()),
269                    Form("%s%sphi distribution",titlePrefix0[i].Data(),titlePrefix[j].Data()), 100, -1., 9.);
270       fList2->AddAtAndExpand(h, k2Phi+i*kNhistPair+j);
271     }
272   }
273   
274   // global counters of tracks:
275   // - reconstructible = number of reconstructible tracks
276   // - reconstructed   = number of reconstructed tracks
277   // - matched         = number of reconstructed tracks matched with a simulated one (reconstructible or not)
278   // - matchedyet      = number of reconstructed tracks matched with a simulated one that is not reconstructible
279   // - decay           = number of reconstructed tracks matched with a decay chain (reconstructible or not)
280   // - decayyet        = number of reconstructed tracks matched with a decay chain that is not reconstructible
281   // - fake            = number of fake tracks
282   // - connected       = number of fake tracks connected to a reconstructible simulated track
283   // - additional      = number of additional (fake) tracks compared to the number of reconstructible ones
284   fTrackCounters = new AliCounterCollection(GetOutputSlot(2)->GetContainer()->GetName());
285   fTrackCounters->AddRubric("track", "reconstructible/reconstructed/matched/matchedyet/decay/decayyet/fake/connected/additional");
286   fTrackCounters->AddRubric("run", 1000000);
287   fTrackCounters->AddRubric("trig", "yes/no/unknown");
288   fTrackCounters->AddRubric("selected", "yes/no");
289   fTrackCounters->AddRubric("acc", "in/out/unknown");
290   TString centralityClasses = "5/10/15/20/25/30/35/40/45/50/55/60/65/70/75/80/85/90/95/100";
291   fTrackCounters->AddRubric("cent", centralityClasses.Data());
292   fTrackCounters->Init();
293   
294   // detailled counters of decays and fake tracks:
295   fFakeTrackCounters = new AliCounterCollection(GetOutputSlot(3)->GetContainer()->GetName());
296   fFakeTrackCounters->AddRubric("position", "matched/decay/decayyet/matchedyet/fake/connected/additional");
297   fFakeTrackCounters->AddRubric("label", "matched/decay/decayyet/matchedyet/fake/connected/additional");
298   fFakeTrackCounters->AddRubric("run", 1000000);
299   fFakeTrackCounters->AddRubric("file", 1000000);
300   fFakeTrackCounters->AddRubric("event", 1000000);
301   fFakeTrackCounters->AddRubric("trig", "yes/no/unknown");
302   fFakeTrackCounters->AddRubric("selected", "yes/no");
303   fFakeTrackCounters->AddRubric("acc", "in/out/unknown");
304   fFakeTrackCounters->AddRubric("cent", centralityClasses.Data());
305   fFakeTrackCounters->Init();
306   
307   // counters of tracks matched by position or by using MC labels
308   fMatchedTrackCounters = new AliCounterCollection(GetOutputSlot(4)->GetContainer()->GetName());
309   fMatchedTrackCounters->AddRubric("position", "matched/decay/decayyet/matchedyet/fake");
310   fMatchedTrackCounters->AddRubric("label", "matched/decay/decayyet/matchedyet/fake/matchedother");
311   fMatchedTrackCounters->AddRubric("run", 1000000);
312   fMatchedTrackCounters->AddRubric("trig", "yes/no");
313   fMatchedTrackCounters->AddRubric("selected", "yes/no");
314   fMatchedTrackCounters->AddRubric("acc", "in/out");
315   fMatchedTrackCounters->AddRubric("cent", centralityClasses.Data());
316   fMatchedTrackCounters->Init();
317   
318   // global counters of events
319   // - any             = total number of events with reconstructed tracks
320   // - fake            = number of events with fake track(s)
321   // - notconnected    = number of events with fake tracks that are not connected to a reconstructible simulated track
322   // - additional      = number of events with additional (fake) tracks compared to the number of reconstructible ones
323   // - matchedyet      = number of events with reconstructed tracks matched with a simulated one that is not reconstructible
324   // if trig = yes: only the tracks matched with the trigger are considered in the above logic
325   fEventCounters = new AliCounterCollection(GetOutputSlot(5)->GetContainer()->GetName());
326   fEventCounters->AddRubric("event", "any/fake/notconnected/additional/matchedyet");
327   fEventCounters->AddRubric("run", 1000000);
328   fEventCounters->AddRubric("trig", "any/yes");
329   fEventCounters->AddRubric("selected", "yes/no");
330   fEventCounters->AddRubric("cent", centralityClasses.Data());
331   fEventCounters->Init();
332   
333   // global counters of track pairs:
334   // - reconstructible = number of reconstructible track pairs
335   // - reconstructed   = number of reconstructed track pairs
336   // - matched         = number of reconstructed track pairs fully matched with a simulated one (reconstructible or not)
337   // - 1fake           = number of reconstructed track pairs made of one matched track and one fake
338   // - 2fakes          = number of reconstructed track pairs fully fake
339   fPairCounters = new AliCounterCollection(GetOutputSlot(7)->GetContainer()->GetName());
340   fPairCounters->AddRubric("pair", "reconstructible/reconstructed/matched/1fake/2fakes");
341   fPairCounters->AddRubric("run", 1000000);
342   fPairCounters->AddRubric("trig", "0/1/2");
343   fPairCounters->AddRubric("selected", "yes/no");
344   fPairCounters->AddRubric("acc", "in/out/unknown");
345   fPairCounters->AddRubric("cent", centralityClasses.Data());
346   fPairCounters->Init();
347   
348   // Disable printout of AliMCEvent
349   AliLog::SetClassDebugLevel("AliMCEvent",-1);
350   
351   // Post data at least once per task to ensure data synchronisation (required for merging)
352   PostData(1, fList);
353   PostData(2, fTrackCounters);
354   PostData(3, fFakeTrackCounters);
355   PostData(4, fMatchedTrackCounters);
356   PostData(5, fEventCounters);
357   PostData(6, fList2);
358   PostData(7, fPairCounters);
359 }
360
361 //________________________________________________________________________
362 void AliAnalysisTaskMuonFakes::UserExec(Option_t *)
363 {
364   /// Process event: looks for fakes...
365   
366   // check that reconstructions parameters for that run have been properly set
367   if (fSigmaCut < 0) return;
368   
369   if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
370   
371   // check physics selection
372   TString selected = (fInputHandler && fInputHandler->IsEventSelected() != 0) ? "selected:yes" : "selected:no";
373   
374   // current file name
375   if (fDisableDetailedCounters) fCurrentFileName = "any";
376   else {
377     fCurrentFileName = CurrentFileName();
378     fCurrentFileName.ReplaceAll("alien://","");
379     fCurrentFileName.ReplaceAll("/","\\");
380     fCurrentFileName.ReplaceAll(":",";");
381   }
382   
383   // Load ESD event
384   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
385   if (!esd) {
386     AliError("Cannot get input event");
387     return;
388   }      
389   
390   // event number in current file
391   TString eventNumberInFile = (fDisableDetailedCounters) ? "event:any" : Form("event:%d",esd->GetEventNumberInFile());
392   
393   // current centrality class
394   TString centrality = "cent:";
395   Double_t centralityValue = esd->GetCentrality()->GetCentralityPercentile("V0M");
396   TObjArray* centralylimits = fTrackCounters->GetKeyWords("cent").Tokenize(",");
397   TObjString* limit = 0x0;
398   TIter nextLimit(centralylimits);
399   while ((limit = static_cast<TObjString*>(nextLimit()))) {
400     if (centralityValue < limit->String().Atoi()) {
401       centrality += limit->GetName();
402       break;
403     }
404   }
405   if (!limit) centrality += static_cast<TObjString*>(centralylimits->Last())->GetName();
406   delete centralylimits;
407   
408   // Load MC event 
409   AliMCEventHandler *mcH = 0;
410   if(MCEvent()) mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
411   
412   // get reconstructed and simulated tracks
413   AliMUONRecoCheck rc(esd,mcH);
414   AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(-1, kFALSE);
415   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
416   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
417   if (!muonTrackStore || !trackRefStore) return;
418   
419   // loop over trackRefs
420   Int_t nMuPlus[2] = {0, 0};
421   Int_t nMuMinus[2] = {0, 0};
422   TIter next(trackRefStore->CreateIterator());
423   AliMUONTrack* trackRef;
424   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
425     
426     // skip trackRefs that are not reconstructible
427     if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
428     
429     // trigger condition
430     AliMUONTriggerTrack *trigRef = static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(trackRef->GetUniqueID()));
431     Bool_t trigger = (trigRef && trigRef->GetPtCutLevel() > 0);
432     Int_t iTrig = trigger ? 1 : 0;
433     TString trig = trigger ? "trig:yes" : "trig:no";
434     
435     // count muons
436     if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus[iTrig]++;
437     else nMuMinus[iTrig]++;
438     
439     // fill global counters
440     fTrackCounters->Count(Form("track:reconstructible/run:%d/%s/%s/acc:unknown/%s", fCurrentRunNumber, trig.Data(), selected.Data(), centrality.Data()));
441     
442   }
443   
444   // fill global counters
445   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:0/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[0]*nMuMinus[0]);
446   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:1/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[1]*nMuMinus[0]+nMuPlus[0]*nMuMinus[1]);
447   fPairCounters->Count(Form("pair:reconstructible/run:%d/trig:2/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nMuPlus[1]*nMuMinus[1]);
448   
449   // loop over ESD tracks
450   Int_t nTrackerTracks = 0;
451   Bool_t containTrack[2] = {kFALSE, kFALSE};
452   Bool_t containFakeTrack[2] = {kFALSE, kFALSE};
453   Bool_t containMatchedYetTrack[2] = {kFALSE, kFALSE};
454   AliMUONVTrackStore *usedTrackRefStore = AliMUONESDInterface::NewTrackStore();
455   AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
456   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
457   TArrayI mcLabels(nTracks);
458   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
459     
460     AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
461     
462     // skip ghosts
463     if (!IsSelected(*esdTrack)) continue;
464     containTrack[0] = kTRUE;
465     
466     // trigger condition
467     Bool_t trigger = esdTrack->ContainTriggerData();
468     TString trig = trigger ? "trig:yes" : "trig:no";
469     if (trigger) containTrack[1] = kTRUE;
470     
471     // acceptance condition
472     Double_t rAbs = esdTrack->GetRAtAbsorberEnd();
473     Double_t thetaTrackAbsEnd = TMath::ATan(rAbs/505.) * TMath::RadToDeg();
474     Double_t eta = esdTrack->Eta();
475     Bool_t inAcc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5);
476     TString acc = inAcc ? "acc:in" : "acc:out";
477     
478     // fill global counters
479     if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) nTrackerTracks++;
480     fTrackCounters->Count(Form("track:reconstructed/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
481     
482     // find the corresponding MUON track
483     AliMUONTrack* muonTrack = static_cast<AliMUONTrack*>(muonTrackStore->FindObject(esdTrack->GetUniqueID()));
484     
485     // get track info
486     Int_t nClusters = esdTrack->GetNClusters();
487     Int_t nChamberHit = 0;
488     for (Int_t ich=0; ich<10; ich++) if (esdTrack->IsInMuonClusterMap(ich)) nChamberHit++;
489     Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
490     Double_t p = esdTrack->P();
491     Double_t pT = esdTrack->Pt();
492     Double_t phi = esdTrack->Phi();
493     Double_t dca = esdTrack->GetDCA();
494     Double_t pU = esdTrack->PUncorrected();
495     Double_t pdca = 0.5*(p+pU)*dca;
496     
497     // fill global histograms
498     FillHistoTrack(0, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
499     if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc))
500       FillHistoTrack(kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
501     
502     // try to match, by position, the reconstructed track with a simulated one
503     Int_t nMatchClustersByPosition = 0;
504     AliMUONTrack* matchedTrackRefByPosition = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByPosition, kFALSE, fSigmaCut);
505     Bool_t isMatchedYetByPosition = kFALSE;
506     Bool_t isRecoDecayByPosition = kFALSE;
507     Int_t decayLabelByPosition = -1, lastChDecayByPosition = 0;
508     if (!matchedTrackRefByPosition || !matchedTrackRefByPosition->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) {
509       decayLabelByPosition = IsDecayByPosition(*muonTrack, *trackRefStore, *usedTrackRefStore, isRecoDecayByPosition, lastChDecayByPosition);
510       if (decayLabelByPosition >= 0) matchedTrackRefByPosition = 0x0;
511       else if (matchedTrackRefByPosition) isMatchedYetByPosition = kTRUE;
512     }
513     Bool_t isFakeByPosition = (!matchedTrackRefByPosition && decayLabelByPosition < 0);
514     
515     // try to match, by using MC labels, the reconstructed track with a simulated one
516     Int_t nMatchClustersByLabel = 0;
517     AliMUONTrack* matchedTrackRefByLabel = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClustersByLabel, kTRUE, fSigmaCut);
518     Bool_t isMatchedYetByLabel = kFALSE;
519     Bool_t isRecoDecayByLabel = kFALSE;
520     Int_t decayLabelByLabel = -1, lastChDecayByLabel = 0;
521     if (!matchedTrackRefByLabel || !matchedTrackRefByLabel->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) {
522       decayLabelByLabel = IsDecayByLabel(*muonTrack, isRecoDecayByLabel, lastChDecayByLabel);
523       if (decayLabelByLabel >= 0) matchedTrackRefByLabel = 0x0;
524       else if (matchedTrackRefByLabel) isMatchedYetByLabel = kTRUE;
525     }
526     Bool_t isFakeByLabel = (!matchedTrackRefByLabel && decayLabelByLabel < 0);
527     
528     // fill global counters
529     TString positionCase = "position:";
530     if (isMatchedYetByPosition) positionCase += "matchedyet";
531     else if (isRecoDecayByPosition) positionCase += "decay";
532     else if (decayLabelByPosition >= 0) positionCase += "decayyet";
533     else if (isFakeByPosition) positionCase += "fake";
534     else positionCase += "matched";
535     TString labelCase = "label:";
536     if (isMatchedYetByLabel) labelCase += "matchedyet";
537     else if (isRecoDecayByLabel) labelCase += "decay";
538     else if (decayLabelByLabel >= 0) labelCase += "decayyet";
539     else if (isFakeByLabel) labelCase += "fake";
540     else labelCase += "matched";
541     if (!matchedTrackRefByPosition || isMatchedYetByPosition || !matchedTrackRefByLabel || isMatchedYetByLabel)
542       fFakeTrackCounters->Count(Form("%s/%s/run:%d/file:%s/%s/%s/%s/%s/%s", positionCase.Data(), labelCase.Data(), fCurrentRunNumber,
543                                      fCurrentFileName.Data(), eventNumberInFile.Data(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
544     if (matchedTrackRefByLabel && matchedTrackRefByPosition &&
545         matchedTrackRefByLabel->GetUniqueID() != matchedTrackRefByPosition->GetUniqueID()) labelCase = "label:matchedother";
546     fMatchedTrackCounters->Count(Form("%s/%s/run:%d/%s/%s/%s/%s", positionCase.Data(), labelCase.Data(),
547                                       fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
548     
549     // take actions according to the matching result we are interested in
550     Int_t nMatchClusters = 0;
551     AliMUONTrack* matchedTrackRef = 0x0;
552     Bool_t isFake = kFALSE, isMatchedYet = kFALSE, isRecoDecay = kFALSE;
553     Int_t decayLabel = -1;
554     if (fCombineMCId) {
555       
556       // choose the best, or the only available, matched track
557       if (matchedTrackRefByPosition && matchedTrackRefByLabel && ((!isMatchedYetByPosition && !isMatchedYetByLabel) ||
558                                                                   (isMatchedYetByPosition && isMatchedYetByLabel))) {
559         
560         nMatchClusters = TMath::Max(nMatchClustersByPosition, nMatchClustersByLabel);
561         matchedTrackRef = (nMatchClusters == nMatchClustersByPosition) ? matchedTrackRefByPosition : matchedTrackRefByLabel;
562         isMatchedYet = isMatchedYetByPosition;
563         
564       } else if (matchedTrackRefByPosition && (!isMatchedYetByPosition || isFakeByLabel)) {
565         
566         nMatchClusters = nMatchClustersByPosition;
567         matchedTrackRef = matchedTrackRefByPosition;
568         isMatchedYet = isMatchedYetByPosition;
569         
570       } else if (matchedTrackRefByLabel && (!isMatchedYetByLabel || isFakeByPosition)) {
571         
572         nMatchClusters = nMatchClustersByLabel;
573         matchedTrackRef = matchedTrackRefByLabel;
574         isMatchedYet = isMatchedYetByLabel;
575         
576         // choose the best (even if it does not matter here), or the only available, decay chain
577       } else if (decayLabelByPosition >= 0 && decayLabelByLabel >= 0 && ((isRecoDecayByPosition && isRecoDecayByLabel) ||
578                                                                          (!isRecoDecayByPosition && !isRecoDecayByLabel))) {
579         
580         decayLabel = (lastChDecayByLabel > lastChDecayByPosition) ? decayLabelByLabel : decayLabelByPosition;
581         isRecoDecay = isRecoDecayByPosition;
582         
583       } else if (decayLabelByPosition >= 0 && (isRecoDecayByPosition || decayLabelByLabel < 0)) {
584         
585         decayLabel = decayLabelByPosition;
586         isRecoDecay = isRecoDecayByPosition;
587         
588       } else if (decayLabelByLabel >= 0) {
589         
590         decayLabel = decayLabelByLabel;
591         isRecoDecay = isRecoDecayByLabel;
592         
593         // no matched track and no decay chain... It must be fakes!
594       } else isFake = kTRUE;
595       
596     } else if (fUseLabel) {
597       
598       // choose the track matched by MC labels
599       nMatchClusters = nMatchClustersByLabel;
600       matchedTrackRef = matchedTrackRefByLabel;
601       isMatchedYet = isMatchedYetByLabel;
602       decayLabel = decayLabelByLabel;
603       isRecoDecay = isRecoDecayByLabel;
604       isFake = isFakeByLabel;
605       
606     } else {
607       
608       // choose the track matched by position
609       nMatchClusters = nMatchClustersByPosition;
610       matchedTrackRef = matchedTrackRefByPosition;
611       isMatchedYet = isMatchedYetByPosition;
612       decayLabel = decayLabelByPosition;
613       isRecoDecay = isRecoDecayByPosition;
614       isFake = isFakeByPosition;
615       
616     }
617     
618     if (matchedTrackRef) {
619       
620       // fill global counters
621       fTrackCounters->Count(Form("track:matched/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
622       
623       // track matched with a trackRef that is not reconstructible
624       if (isMatchedYet) {
625         
626         containMatchedYetTrack[0] = kTRUE;
627         if (trigger) containMatchedYetTrack[1] = kTRUE;
628         
629         // fill global counters
630         fTrackCounters->Count(Form("track:matchedyet/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
631         
632         // fill histograms
633         FillHistoTrack(2, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
634         if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) 
635           FillHistoTrack(2+kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
636         
637       }
638       
639       // fill histograms
640       if (nClusters > 0) ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kFractionOfMatchedClusters))->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
641       ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kNumberOfClustersMC))->Fill(matchedTrackRef->GetNClusters());
642       FillHistoTrack(1, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
643       if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) 
644         FillHistoTrack(1+kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
645       
646       // flag matched tracks
647       mcLabels[iTrack] = matchedTrackRef->GetUniqueID();
648       
649       // move already matched trackRefs
650       usedTrackRefStore->Add(*matchedTrackRef);
651       trackRefStore->Remove(*matchedTrackRef);
652       
653     } else {
654       
655       if (decayLabel >= 0) {
656         
657         // fill global counters
658         fTrackCounters->Count(Form("track:decay/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
659         
660         // track matched with a decay that has not be tagged reconstructible
661         if (!isRecoDecay) {
662           
663           // fill global counters
664           fTrackCounters->Count(Form("track:decayyet/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
665           
666           // fill histograms
667           FillHistoTrack(4, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
668           if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) 
669             FillHistoTrack(4+kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
670           
671         }
672         
673         // fill histograms
674         FillHistoTrack(3, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
675         if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) 
676           FillHistoTrack(3+kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
677         
678         // flag decay tracks
679         mcLabels[iTrack] = decayLabel;
680         
681       }
682       
683       if (isFake || fDecayAsFake) {
684         
685         containFakeTrack[0] = kTRUE;
686         if (trigger) containFakeTrack[1] = kTRUE;
687         
688         // fill global counters
689         fTrackCounters->Count(Form("track:fake/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
690         
691         // fill histograms
692         FillHistoTrack(5, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
693         if ((!fMatchTrig || trigger) && (!fApplyAccCut || inAcc)) 
694           FillHistoTrack(5+kNhistTrack, nClusters, nChamberHit, normalizedChi2, p, pT, eta, phi, dca, thetaTrackAbsEnd, pdca, rAbs);
695         
696         // flag fake tracks
697         mcLabels[iTrack] = -1;
698         
699         // store fake tracks
700         fakeTrackStore->Add(*muonTrack);
701         
702       }
703       
704     }
705     
706   } // end of loop over ESD tracks
707   
708   // fill histogram and global counters
709   ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kNumberOfTracks))->Fill(nTrackerTracks);
710   if (containTrack[0]) fEventCounters->Count(Form("event:any/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
711   if (containTrack[1]) fEventCounters->Count(Form("event:any/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
712   if (containFakeTrack[0]) fEventCounters->Count(Form("event:fake/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
713   if (containFakeTrack[1]) fEventCounters->Count(Form("event:fake/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
714   if (containMatchedYetTrack[0]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
715   if (containMatchedYetTrack[1]) fEventCounters->Count(Form("event:matchedyet/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
716   
717   // count the number of not connected and additional fake tracks
718   if (fakeTrackStore->GetSize() > 0) {
719     
720     // remove the most connected fake tracks
721     Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, selected, centrality);
722     
723     if (fakeTrackStore->GetSize() > 0) {
724       
725       // fill global counters
726       fEventCounters->Count(Form("event:notconnected/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
727       
728       // check status of remaining fakes with respect to the matching with trigger
729       Bool_t containMatchedFake = kFALSE;
730       Bool_t containUnmatchedFake = kFALSE;
731       AliMUONTrack* fakeTrack = 0x0;
732       TIter next3(fakeTrackStore->CreateIterator());
733       while ( ( fakeTrack = static_cast<AliMUONTrack*>(next3()) ) ) {
734         if (fakeTrack->GetMatchTrigger() > 0) containMatchedFake = kTRUE;
735         else containUnmatchedFake = kTRUE;
736       }
737       
738       // fill global counters
739       if (containMatchedFake) fEventCounters->Count(Form("event:notconnected/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
740       
741       // remove the remaining free reconstructible tracks
742       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
743       
744       if (nAdditionalTracks > 0) {
745         
746         // fill histogram and global counters
747         ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kNumberOfAdditionalTracks))->Fill(nAdditionalTracks);
748         fEventCounters->Count(Form("event:additional/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
749         if (!containUnmatchedFake) { // all matched
750           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
751           fTrackCounters->Count(Form("track:additional/run:%d/trig:yes/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
752           fFakeTrackCounters->Count(Form("position:additional/label:additional/run:%d/file:%s/%s/trig:yes/%s/acc:unknown/%s", fCurrentRunNumber,
753                                          fCurrentFileName.Data(), eventNumberInFile.Data(), selected.Data(), centrality.Data()), nAdditionalTracks);
754         } else if (!containMatchedFake) { // none matched
755           fTrackCounters->Count(Form("track:additional/run:%d/trig:no/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
756           fFakeTrackCounters->Count(Form("position:additional/label:additional/run:%d/file:%s/%s/trig:no/%s/acc:unknown/%s", fCurrentRunNumber,
757                                          fCurrentFileName.Data(), eventNumberInFile.Data(), selected.Data(), centrality.Data()), nAdditionalTracks);
758         } else { // mixed
759           fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
760           fTrackCounters->Count(Form("track:additional/run:%d/trig:unknown/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
761           fFakeTrackCounters->Count(Form("position:additional/label:additional/run:%d/file:%s/%s/trig:unknown/%s/acc:unknown/%s", fCurrentRunNumber,
762                                          fCurrentFileName.Data(), eventNumberInFile.Data(), selected.Data(), centrality.Data()), nAdditionalTracks);
763         }
764         
765       }
766       
767     }
768     
769   }
770   
771   // clean memory
772   delete usedTrackRefStore;
773   delete fakeTrackStore;
774   
775   // double loop over ESD tracks, build pairs and fill histograms and counters according to their label
776   TLorentzVector vMu1, vMu2, vDiMu;
777   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
778     AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
779     
780     // skip ghosts
781     if (!IsSelected(*muonTrack1)) continue;
782     
783     // get track info
784     Bool_t trigger1 = muonTrack1->ContainTriggerData();
785     Double_t thetaAbs1 = TMath::ATan(muonTrack1->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
786     Double_t eta1 = muonTrack1->Eta();
787     Bool_t acc1 = (thetaAbs1 >= 2. && thetaAbs1 <= 10. && eta1 >= -4. && eta1 <= -2.5);
788     Short_t charge1 = muonTrack1->Charge();
789     Int_t label1 = mcLabels[iTrack1];
790     muonTrack1->LorentzP(vMu1);
791     
792     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
793       AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
794       
795       // skip ghosts
796       if (!IsSelected(*muonTrack2)) continue;
797       
798       // keep only opposite sign pairs
799       Short_t charge2 = muonTrack2->Charge();
800       if (charge1*charge2 > 0) continue;
801       
802       // get track info
803       Bool_t trigger2 = muonTrack2->ContainTriggerData();
804       Double_t thetaAbs2 = TMath::ATan(muonTrack2->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
805       Double_t eta2 = muonTrack2->Eta();
806       Bool_t acc2 = (thetaAbs2 >= 2. && thetaAbs2 <= 10. && eta2 >= -4. && eta2 <= -2.5);
807       Int_t label2 = mcLabels[iTrack2];
808       muonTrack2->LorentzP(vMu2);
809       
810       // compute kinematics of the pair
811       vDiMu = vMu1 + vMu2;
812       Float_t mass = vDiMu.M();
813       Float_t p = vDiMu.P();
814       Float_t pt = vDiMu.Pt();
815       Float_t y = vDiMu.Rapidity();
816       Float_t eta = vDiMu.Eta();
817       Float_t phi = vDiMu.Phi();
818       if (phi < 0) phi += 2.*TMath::Pi();
819       
820       // trigger condition
821       TString trig = "trig:";
822       if (trigger1 && trigger2) trig += "2";
823       else if (trigger1 || trigger2) trig += "1";
824       else trig += "0";
825       
826       // acceptance condition
827       Bool_t inAcc = (acc1 && acc2 && y >= -4. && y <= -2.5);
828       TString acc = inAcc ? "acc:in" : "acc:out";
829       
830       // fill global histograms
831       FillHistoPair(0, mass, p, pt, y, eta, phi);
832       if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc))
833         FillHistoPair(kNhistPair, mass, p, pt, y, eta, phi);
834       
835       TString pair = "pair:";
836       
837       // fill histograms according to labels
838       if (label1 >= 0 && label2 >= 0) {
839         
840         pair += "matched";
841         
842         FillHistoPair(1, mass, p, pt, y, eta, phi);
843         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc))
844           FillHistoPair(1+kNhistPair, mass, p, pt, y, eta, phi);
845         
846       } else if (label1 >= 0 || label2 >= 0) {
847         
848         pair += "1fake";
849         
850         FillHistoPair(2, mass, p, pt, y, eta, phi);
851         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc))
852           FillHistoPair(2+kNhistPair, mass, p, pt, y, eta, phi);
853         
854       } else {
855         
856         pair += "2fakes";
857         
858         FillHistoPair(3, mass, p, pt, y, eta, phi);
859         if ((!fMatchTrig || (trigger1 && trigger2)) && (!fApplyAccCut || inAcc))
860           FillHistoPair(3+kNhistPair, mass, p, pt, y, eta, phi);
861         
862       }
863       
864       // fill global counters
865       fPairCounters->Count(Form("pair:reconstructed/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
866       fPairCounters->Count(Form("%s/run:%d/%s/%s/%s/%s", pair.Data(), fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
867       
868     }
869     
870   }
871   
872   // Post final data
873   PostData(1, fList);
874   PostData(2, fTrackCounters);
875   PostData(3, fFakeTrackCounters);
876   PostData(4, fMatchedTrackCounters);
877   PostData(5, fEventCounters);
878   PostData(6, fList2);
879   PostData(7, fPairCounters);
880 }
881
882 //________________________________________________________________________
883 void AliAnalysisTaskMuonFakes::NotifyRun()
884 {
885   /// Prepare processing of new run: load corresponding OCDB objects...
886   
887   // load OCDB objects only once
888   if (fSigmaCut > 0) return;
889   
890   // set OCDB location
891   AliCDBManager* cdbm = AliCDBManager::Instance();
892   if (cdbm->IsDefaultStorageSet()) printf("FakeTask: CDB default storage already set!\n");
893   else cdbm->SetDefaultStorage(fRecoParamLocation.Data());
894   if (cdbm->GetRun() > -1) printf("FakeTask: run number already set!\n");
895   else cdbm->SetRun(fCurrentRunNumber);
896   
897   // load necessary data from OCDB
898   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
899   if (!recoParam) {
900     fRequestedStationMask = 0;
901     fRequest2ChInSameSt45 = kFALSE;
902     fSigmaCut = -1.;
903     AliError("--> skip this run");
904     return;
905   }
906   
907   // compute the mask of requested stations from recoParam
908   fRequestedStationMask = 0;
909   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i );
910   
911   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
912   fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
913   
914   // get sigma cut to associate clusters with TrackRefs from recoParam if not already set manually
915   if (fExternalSigmaCut > 0) fSigmaCut = fExternalSigmaCut;
916   else if (recoParam->ImproveTracks()) fSigmaCut = recoParam->GetSigmaCutForImprovement();
917   else fSigmaCut = recoParam->GetSigmaCutForTracking();
918   
919   // get the trackCuts for this run
920   if (fMuonTrackCuts) fMuonTrackCuts->SetRun(fInputHandler);
921 }
922
923 //________________________________________________________________________
924 void AliAnalysisTaskMuonFakes::Terminate(Option_t *)
925 {
926   /// Draw results to the screen and print statistics.
927   
928   // recover output objects
929   fList = static_cast<TObjArray*> (GetOutputData(1));
930   fList2 = static_cast<TObjArray*> (GetOutputData(6));
931   if (!fList || !fList2) return;
932   fTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(2));
933   fFakeTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(3));
934   fMatchedTrackCounters = static_cast<AliCounterCollection*> (GetOutputData(4));
935   fEventCounters = static_cast<AliCounterCollection*> (GetOutputData(5));
936   fPairCounters = static_cast<AliCounterCollection*> (GetOutputData(7));
937   
938   TString extention = GetName();
939   extention.ReplaceAll("MUONFakes","");
940   
941   // add canvas to compare histograms
942   fCanvases = new TObjArray(13);
943   fCanvases->SetOwner();
944   
945   TString nameSuffix[2] = {"", "S"};
946   TString titleSuffix[2] = {"", "selected "};
947   for (Int_t j = 0; j < 2; j++) {
948     
949     TCanvas *cFakesSummary11 = new TCanvas(Form("cTracks11%s_%s",nameSuffix[j].Data(),extention.Data()),
950                                            Form("distributions of %stracks (%s)",titleSuffix[j].Data(),extention.Data()),900,900);
951     fCanvases->AddAtAndExpand(cFakesSummary11, 0+7*j);
952     TCanvas *cFakesSummary12 = new TCanvas(Form("cTracks12%s_%s",nameSuffix[j].Data(),extention.Data()),
953                                            Form("detailled distributions of %stracks (%s)",titleSuffix[j].Data(),extention.Data()),900,900);
954     fCanvases->AddAtAndExpand(cFakesSummary12, 1+7*j);
955     TCanvas *cFakesSummary13 = new TCanvas(Form("cTracks13%s_%s",nameSuffix[j].Data(),extention.Data()),
956                                            Form("p*DCA distributions of %stracks (%s)",titleSuffix[j].Data(),extention.Data()),600,300);
957     fCanvases->AddAtAndExpand(cFakesSummary13, 2+7*j);
958     TCanvas *cFakesSummary14 = new TCanvas(Form("cTracks14%s_%s",nameSuffix[j].Data(),extention.Data()),
959                                            Form("detailled p*DCA distributions of %stracks (%s)",titleSuffix[j].Data(),extention.Data()),600,300);
960     fCanvases->AddAtAndExpand(cFakesSummary14, 3+7*j);
961     TCanvas *cFakesSummary21 = new TCanvas(Form("cTracks21%s_%s",nameSuffix[j].Data(),extention.Data()),
962                                           Form("correlations at the %strack level (%s)",titleSuffix[j].Data(),extention.Data()),1200,600);
963     fCanvases->AddAtAndExpand(cFakesSummary21, 4+7*j);
964     TCanvas *cFakesSummary22 = new TCanvas(Form("cTracks22%s_%s",nameSuffix[j].Data(),extention.Data()),
965                                           Form("detailled correlations at the %strack level (%s)",titleSuffix[j].Data(),extention.Data()),1200,600);
966     fCanvases->AddAtAndExpand(cFakesSummary22, 5+7*j);
967     TCanvas *cFakesSummary3 = new TCanvas(Form("cPairs%s_%s",nameSuffix[j].Data(),extention.Data()),
968                                           Form("distributions of %spairs (%s)",titleSuffix[j].Data(),extention.Data()),900,600);
969     fCanvases->AddAtAndExpand(cFakesSummary3, 6+7*j);
970     
971     // display
972     Int_t iHist1[9] = {kNumberOfClusters, kNumberOfChamberHit, kChi2PerDof, kDCA, kRAbs, kEta, kP, kPt, kPhi};
973     cFakesSummary11->Divide(3,3);
974     for (Int_t i=0; i<9; i++) {
975       cFakesSummary11->cd(i+1);
976       cFakesSummary11->GetPad(i+1)->SetLogy();
977       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack))->SetMinimum(0.5);
978       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack))->DrawCopy();
979       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+1))->SetLineColor(4);
980       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+1))->DrawCopy("sames");
981       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+3))->SetLineColor(kViolet-3);
982       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+3))->SetFillColor(kViolet-3);
983       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+3))->SetFillStyle(3018);
984       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+3))->DrawCopy("sames");
985       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->SetLineColor(2);
986       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->SetFillColor(2);
987       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->SetFillStyle(3017);
988       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->DrawCopy("sames");
989     }
990     
991     cFakesSummary12->Divide(3,3);
992     for (Int_t i=0; i<9; i++) {
993       cFakesSummary12->cd(i+1);
994       cFakesSummary12->GetPad(i+1)->SetLogy();
995       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack))->SetMinimum(0.5);
996       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack))->DrawCopy();
997       TH1F *hClone = (TH1F*) fList->UncheckedAt(iHist1[i]+j*kNhistTrack+1)->Clone();
998       hClone->Add(((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+2)), -1);
999       hClone->SetLineColor(4);
1000       hClone->DrawCopy("sames");
1001       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+2))->SetLineColor(7);
1002       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+2))->DrawCopy("sames");
1003       hClone = (TH1F*) fList->UncheckedAt(iHist1[i]+j*kNhistTrack+3)->Clone();
1004       hClone->Add(((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+4)), -1);
1005       hClone->SetLineColor(3);
1006       hClone->SetFillStyle(0);
1007       hClone->DrawCopy("sames");
1008       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+4))->SetLineColor(32);
1009       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+4))->DrawCopy("sames");
1010       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->SetLineColor(2);
1011       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->SetFillStyle(0);
1012       ((TH1F*)fList->UncheckedAt(iHist1[i]+j*kNhistTrack+5))->DrawCopy("sames");
1013     }
1014     
1015     Int_t iHist2[2] = {kPDCA23, kPDCA310};
1016     cFakesSummary13->Divide(2,1);
1017     for (Int_t i=0; i<2; i++) {
1018       cFakesSummary13->cd(i+1);
1019       cFakesSummary13->GetPad(i+1)->SetLogy();
1020       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack))->SetMinimum(0.5);
1021       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack))->DrawCopy();
1022       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+1))->SetLineColor(4);
1023       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+1))->DrawCopy("sames");
1024       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+3))->SetLineColor(kViolet-3);
1025       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+3))->SetFillColor(kViolet-3);
1026       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+3))->SetFillStyle(3018);
1027       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+3))->DrawCopy("sames");
1028       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->SetLineColor(2);
1029       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->SetFillColor(2);
1030       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->SetFillStyle(3017);
1031       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->DrawCopy("sames");
1032     }
1033     
1034     cFakesSummary14->Divide(2,1);
1035     for (Int_t i=0; i<2; i++) {
1036       cFakesSummary14->cd(i+1);
1037       cFakesSummary14->GetPad(i+1)->SetLogy();
1038       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack))->SetMinimum(0.5);
1039       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack))->DrawCopy();
1040       TH1F *hClone = (TH1F*) fList->UncheckedAt(iHist2[i]+j*kNhistTrack+1)->Clone();
1041       hClone->Add(((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+2)), -1);
1042       hClone->SetLineColor(4);
1043       hClone->DrawCopy("sames");
1044       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+2))->SetLineColor(7);
1045       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+2))->DrawCopy("sames");
1046       hClone = (TH1F*) fList->UncheckedAt(iHist2[i]+j*kNhistTrack+3)->Clone();
1047       hClone->Add(((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+4)), -1);
1048       hClone->SetLineColor(3);
1049       hClone->SetFillStyle(0);
1050       hClone->DrawCopy("sames");
1051       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+4))->SetLineColor(32);
1052       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+4))->DrawCopy("sames");
1053       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->SetLineColor(2);
1054       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->SetFillStyle(0);
1055       ((TH1F*)fList->UncheckedAt(iHist2[i]+j*kNhistTrack+5))->DrawCopy("sames");
1056     }
1057     
1058     Int_t iHist3[2] = {kChi2PerDofVsNClusters, kChi2PerDofVsNChamberHit};
1059     cFakesSummary21->Divide(2);
1060     for (Int_t i=0; i<2; i++) {
1061       cFakesSummary21->cd(i+1);
1062       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+1))->SetMarkerColor(4);
1063       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+1))->DrawCopy();
1064       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+3))->SetMarkerColor(kViolet-3);
1065       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+3))->SetMarkerStyle(6);
1066       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+3))->DrawCopy("sames");
1067       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->SetMarkerColor(2);
1068       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->SetMarkerStyle(7);
1069       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->DrawCopy("sames");
1070     }
1071     
1072     cFakesSummary22->Divide(2);
1073     for (Int_t i=0; i<2; i++) {
1074       cFakesSummary22->cd(i+1);
1075       TH2F *hClone = (TH2F*) fList->UncheckedAt(iHist3[i]+j*kNhistTrack+1)->Clone();
1076       hClone->Add(((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+2)), -1);
1077       hClone->SetMarkerColor(4);
1078       hClone->DrawCopy();
1079       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+2))->SetMarkerColor(7);
1080       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+2))->SetMarkerStyle(6);
1081       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+2))->DrawCopy("sames");
1082       hClone = (TH2F*) fList->UncheckedAt(iHist3[i]+j*kNhistTrack+3)->Clone();
1083       hClone->Add(((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+4)), -1);
1084       hClone->SetMarkerColor(kViolet-3);
1085       hClone->SetMarkerStyle(6);
1086       hClone->DrawCopy("sames");
1087       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+4))->SetLineColor(32);
1088       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+4))->SetMarkerStyle(6);
1089       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+4))->DrawCopy("sames");
1090       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->SetMarkerColor(2);
1091       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->SetMarkerStyle(7);
1092       ((TH2F*)fList->UncheckedAt(iHist3[i]+j*kNhistTrack+5))->DrawCopy("sames");
1093     }
1094     
1095     Int_t iHist4[6] = {k2Mass, k2P, k2Pt, k2Y, k2Eta, k2Phi};
1096     cFakesSummary3->Divide(3,2);
1097     for (Int_t i=0; i<6; i++) {
1098       cFakesSummary3->cd(i+1);
1099       cFakesSummary3->GetPad(i+1)->SetLogy();
1100       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair))->SetMinimum(0.5);
1101       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair))->DrawCopy();
1102       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+1))->SetLineColor(4);
1103       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+1))->DrawCopy("sames");
1104       TH1F* hClone = (TH1F*) fList2->UncheckedAt(iHist4[i]+j*kNhistPair+2)->Clone();
1105       hClone->Add((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+3));
1106       hClone->SetLineColor(2);
1107       hClone->SetFillColor(2);
1108       hClone->SetFillStyle(3017);
1109       hClone->DrawCopy("sames");
1110       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+3))->SetLineColor(6);
1111       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+3))->SetFillColor(6);
1112       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+3))->SetFillStyle(3018);
1113       ((TH1F*)fList2->UncheckedAt(iHist4[i]+j*kNhistPair+3))->DrawCopy("sames");
1114     }
1115     
1116   }
1117   
1118   // print
1119   if (fTrackCounters && fFakeTrackCounters && fMatchedTrackCounters && fEventCounters) {
1120     printf("\nGlobal statistics of reconstructed tracks matched or not with the trigger:\n");
1121     fTrackCounters->Print("track/trig");
1122     printf("\nGlobal statistics of pathological tracks matched or not with the trigger:\n");
1123     fFakeTrackCounters->Print("label/position/trig");
1124     printf("\nDetailled statistics of tracks matched per label vs position:\n");
1125     fMatchedTrackCounters->Print("label/position");
1126     printf("\nGlobal statistics of events containing pathological tracks:\n");
1127     fEventCounters->Print("event/trig");
1128   }
1129   
1130   if (fPairCounters) {
1131     printf("\nGlobal statistics of reconstructed track pairs matched or not with the trigger:\n");
1132     fPairCounters->Print("pair/trig");
1133   }
1134   
1135   printf("\nREMINDER: results are relevent provided that you use the same recoParams as for the reconstruction\n\n");
1136 }
1137
1138
1139 //________________________________________________________________________
1140 Bool_t AliAnalysisTaskMuonFakes::IsSelected(AliESDMuonTrack &esdTrack)
1141 {
1142   /// return kTRUE if the track pass the section criteria
1143   
1144   // make sure to skip ghosts
1145   if (!esdTrack.ContainTrackerData()) return kFALSE;
1146   
1147   // apply standard track cuts if any
1148   if (fMuonTrackCuts && !fMuonTrackCuts->IsSelected(&esdTrack)) return kFALSE;
1149   
1150   // apply specific chi2 cut if required
1151   if (fChi2Cut > 0. && esdTrack.GetNormalizedChi2() > fChi2Cut) return kFALSE;
1152   
1153   // apply specific pt cut if required
1154   if (fPtCut > 0. && esdTrack.Pt() < fPtCut) return kFALSE;
1155   
1156   return kTRUE;
1157 }
1158
1159
1160 //________________________________________________________________________
1161 void AliAnalysisTaskMuonFakes::FillHistoTrack(Int_t histShift, Int_t nClusters, Int_t nChamberHit, Double_t normalizedChi2,
1162                                               Double_t p, Double_t pT, Double_t eta, Double_t phi, Double_t dca,
1163                                               Double_t thetaTrackAbsEnd, Double_t pdca, Double_t rAbs)
1164 {
1165   /// fill global histograms at track level
1166   ((TH1F*)fList->UncheckedAt(kNumberOfClusters+histShift))->Fill(nClusters);
1167   ((TH1F*)fList->UncheckedAt(kNumberOfChamberHit+histShift))->Fill(nChamberHit);
1168   ((TH1F*)fList->UncheckedAt(kChi2PerDof+histShift))->Fill(normalizedChi2);
1169   ((TH1F*)fList->UncheckedAt(kP+histShift))->Fill(p);
1170   ((TH1F*)fList->UncheckedAt(kPt+histShift))->Fill(pT);
1171   ((TH1F*)fList->UncheckedAt(kEta+histShift))->Fill(eta);
1172   ((TH1F*)fList->UncheckedAt(kPhi+histShift))->Fill(phi);
1173   ((TH1F*)fList->UncheckedAt(kDCA+histShift))->Fill(dca);
1174   if (thetaTrackAbsEnd > 2 && thetaTrackAbsEnd <= 3) ((TH1F*)fList->UncheckedAt(kPDCA23+histShift))->Fill(pdca);
1175   else if (thetaTrackAbsEnd > 3 && thetaTrackAbsEnd < 10) ((TH1F*)fList->UncheckedAt(kPDCA310+histShift))->Fill(pdca);
1176   ((TH1F*)fList->UncheckedAt(kRAbs+histShift))->Fill(rAbs);
1177   ((TH2F*)fList->UncheckedAt(kChi2PerDofVsNClusters+histShift))->Fill(nClusters,normalizedChi2);
1178   ((TH2F*)fList->UncheckedAt(kChi2PerDofVsNChamberHit+histShift))->Fill(nChamberHit,normalizedChi2);
1179 }
1180
1181 //________________________________________________________________________
1182 void AliAnalysisTaskMuonFakes::FillHistoPair(Int_t histShift, Double_t mass, Double_t p, Double_t pt,
1183                                              Double_t y, Double_t eta, Double_t phi)
1184 {
1185   /// fill global histograms at pair level
1186   ((TH1F*)fList2->UncheckedAt(k2Mass+histShift))->Fill(mass);
1187   ((TH1F*)fList2->UncheckedAt(k2P+histShift))->Fill(p);
1188   ((TH1F*)fList2->UncheckedAt(k2Pt+histShift))->Fill(pt);
1189   ((TH1F*)fList2->UncheckedAt(k2Y+histShift))->Fill(y);
1190   ((TH1F*)fList2->UncheckedAt(k2Eta+histShift))->Fill(eta);
1191   ((TH1F*)fList2->UncheckedAt(k2Phi+histShift))->Fill(phi);
1192 }
1193
1194 //________________________________________________________________________
1195 Int_t AliAnalysisTaskMuonFakes::RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
1196                                                      TString &selected, TString &centrality)
1197 {
1198   /// loop over reconstructible TrackRef not associated with reconstructed track:
1199   /// for each of them, find and remove the most connected the fake track, if any,
1200   /// and fill the histograms with the fraction of connected clusters.
1201   /// Return the number of reconstructible track not connected to any fake
1202   
1203   Int_t nFreeMissingTracks = 0;
1204   
1205   // loop over trackRefs
1206   TIter next(trackRefStore.CreateIterator());
1207   AliMUONTrack* trackRef;
1208   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
1209     
1210     // skip not reconstructible trackRefs
1211     if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
1212     
1213     Int_t label = trackRef->GetUniqueID();
1214     
1215     // look for the most connected fake track
1216     AliMUONTrack *connectedFake = 0x0;
1217     Double_t fractionOfConnectedClusters = 0.;
1218     TIter next2(fakeTrackStore.CreateIterator());
1219     AliMUONTrack* fakeTrack;
1220     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
1221       
1222       // get the number of connected clusters
1223       Int_t nConnectedClusters = 0;
1224       if (fUseLabel || fCombineMCId) { // by using the MC label
1225         for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
1226           if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
1227             nConnectedClusters++;
1228       }
1229       if (!fUseLabel || fCombineMCId) { // by comparing cluster/TrackRef positions
1230         Bool_t compTrack[10];
1231         nConnectedClusters = TMath::Max(nConnectedClusters, fakeTrack->FindCompatibleClusters(*trackRef, fSigmaCut, compTrack));
1232       }
1233       
1234       // skip non-connected fake tracks
1235       if (nConnectedClusters == 0) continue;
1236       
1237       // check if it is the most connected fake track
1238       Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
1239       if (f > fractionOfConnectedClusters) {
1240         connectedFake = fakeTrack;
1241         fractionOfConnectedClusters = f;
1242       }
1243       
1244     }
1245     
1246     if (connectedFake) {
1247       
1248       // find the corresponding ESD MUON track
1249       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
1250       if (!esd) {
1251         AliError("Cannot get input event");
1252         return nFreeMissingTracks;
1253       }      
1254       TIter next3(static_cast<TClonesArray*>(esd->FindListObject("MuonTracks")));
1255       AliESDMuonTrack* esdTrack = 0x0;
1256       while ((esdTrack = static_cast<AliESDMuonTrack*>(next3())) && esdTrack->GetUniqueID() != connectedFake->GetUniqueID()) {}
1257       if (!esdTrack) {
1258         AliError("unable to find the corresponding ESD track???");
1259         continue;
1260       }
1261       
1262       // trigger condition
1263       TString trig = (esdTrack->ContainTriggerData()) ? "trig:yes" : "trig:no";
1264       
1265       // acceptance condition
1266       Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
1267       Double_t eta = esdTrack->Eta();
1268       TString acc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5) ? "acc:in" : "acc:out";
1269       
1270       // fill histogram and counters
1271       TString eventNumberInFile = (fDisableDetailedCounters) ? "event:any" : Form("event:%d",esd->GetEventNumberInFile());
1272       ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kFractionOfConnectedClusters))->Fill(fractionOfConnectedClusters);
1273       fTrackCounters->Count(Form("track:connected/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
1274       fFakeTrackCounters->Count(Form("position:connected/label:connected/run:%d/file:%s/%s/%s/%s/%s/%s", fCurrentRunNumber, fCurrentFileName.Data(),
1275                                      eventNumberInFile.Data(), trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
1276       
1277       // remove the most connected fake track
1278       fakeTrackStore.Remove(*connectedFake);
1279       
1280     } else nFreeMissingTracks++;
1281     
1282   }
1283   
1284   return nFreeMissingTracks;
1285   
1286 }
1287
1288 //________________________________________________________________________
1289 Int_t AliAnalysisTaskMuonFakes::IsDecay(Int_t nClusters, Int_t *chId, Int_t *labels,
1290                                         Bool_t &isReconstructible, Int_t &lastCh) const
1291 {
1292   /// Check whether this combination of clusters correspond to a decaying particle or not:
1293   /// More than 50% of clusters, including 1 before and 1 after the dipole, must be connected.
1294   /// - Return the MC label of the most downstream decay product or -1 if not a decay.
1295   /// - "isReconstructible" tells if the combination of matched clusters fulfil the reconstruction criteria.
1296   /// - As soon as we realized the decay chain cannot be tagged as reconstructible, we reject any chain ending
1297   ///   on a chamber equal to or upstream "lastCh" (used to select the best chain in case of multiple choices).
1298   /// - "lastCh" is reset the most downstream chamber of the found decay chain if any.
1299   
1300   Int_t halfCluster = nClusters/2;
1301   
1302   // loop over last clusters (if nClusters left < halfCluster the conditions cannot be fulfilled)
1303   Int_t firstLabel = -1, decayLabel = -1;
1304   isReconstructible = kFALSE;
1305   for (Int_t iCluster1 = nClusters-1; iCluster1 >= halfCluster; iCluster1--) {
1306     
1307     // if the last cluster is not on station 4 or 5 the conditions cannot be fulfilled
1308     if (chId[iCluster1] < 6) break;
1309     
1310     // skip clusters with no label or same label as at the begining of the previous step (already tested)
1311     if (labels[iCluster1] < 0 || labels[iCluster1] == firstLabel) continue;
1312     
1313     // is there any chance the hypothetical decay chain can be tagged reconstructible?
1314     Int_t stationId = chId[iCluster1]/2;
1315     Int_t stationMask = 1 << stationId;
1316     Int_t requestedStations = fRequestedStationMask >> stationId;
1317     Bool_t isValid = ((1 & requestedStations) == requestedStations);
1318     
1319     // if not: check whether we can find a better chain than already found
1320     if (!isValid && chId[iCluster1] <= lastCh) break;
1321     
1322     // count the number of fired chambers on stations 4 & 5
1323     Int_t nChHitInSt45[2] = {0, 0};
1324     nChHitInSt45[stationId-3] = 1;
1325     Int_t currentCh = chId[iCluster1];
1326     
1327     // get the ancestors
1328     TArrayI chainLabels(100);
1329     Int_t nParticles = 0;
1330     Int_t currentLabel = labels[iCluster1];
1331     do {
1332       chainLabels[nParticles++] = currentLabel;
1333       if (nParticles >= chainLabels.GetSize()) chainLabels.Set(2*chainLabels.GetSize());
1334       AliMCParticle* currentParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(currentLabel));
1335       currentLabel = (currentParticle) ? currentParticle->GetMother() : -1;
1336     } while (currentLabel >= 0);
1337     
1338     // Loop over prior clusters
1339     firstLabel = labels[iCluster1];
1340     Int_t nCompatibleLabel = 1;
1341     Int_t currentParticle = 0;
1342     for (Int_t iCluster2 = iCluster1-1; iCluster2 >= 0; iCluster2--) {
1343       
1344       // if the number of clusters left is not enough the conditions cannot be fulfilled
1345       if (iCluster2 < halfCluster-nCompatibleLabel) break;
1346       
1347       if (labels[iCluster2] < 0) continue;
1348       
1349       // check if the cluster belong to the same particle or one of its ancestors
1350       Bool_t matchFound = kFALSE;
1351       for (Int_t iParticle = currentParticle; iParticle < nParticles; iParticle++) {
1352         if (labels[iCluster2] == chainLabels[iParticle]) {
1353           currentParticle = iParticle;
1354           matchFound = kTRUE;
1355           break;
1356         }
1357       }
1358       if (matchFound) nCompatibleLabel++;
1359       else continue;
1360       
1361       // add this station to the mask
1362       stationId = chId[iCluster2]/2;
1363       stationMask |= 1 << stationId;
1364       
1365       // count the number of fired chamber on stations 4 & 5
1366       if (stationId > 2 && chId[iCluster2] < currentCh) {
1367         nChHitInSt45[stationId-3]++;
1368         currentCh = chId[iCluster2];
1369       }
1370       
1371       // check if we matched enough clusters to tag the track as a decay
1372       if (nCompatibleLabel <= halfCluster || chId[iCluster2] > 3 || chainLabels[currentParticle] == firstLabel) continue;
1373       
1374       // check if this chain is better than already found
1375       if (chId[iCluster1] > lastCh) {
1376         decayLabel = firstLabel;
1377         lastCh = chId[iCluster1];
1378       }
1379       
1380       // is there enough matched clusters on station 4 & 5 to make the track reconstructible?
1381       Bool_t isEnoughClOnSt45 = fRequest2ChInSameSt45 ? (nChHitInSt45[0] == 2 || nChHitInSt45[1] == 2)
1382                                                       : (nChHitInSt45[0]+nChHitInSt45[1] >= 2);
1383       
1384       // is there any chance the current decay chain can still be tagged reconstructible?
1385       requestedStations = fRequestedStationMask >> stationId;
1386       isValid = (((stationMask >> stationId) & requestedStations) == requestedStations &&
1387                  (chId[iCluster2] > 5 || isEnoughClOnSt45));
1388       
1389       // if not then we cannot do better with this trial
1390       if (!isValid) break;
1391       
1392       // take in priority the decay chain that can be tagged reconstructible
1393       if (((stationMask & fRequestedStationMask) == fRequestedStationMask) && isEnoughClOnSt45) {
1394         lastCh = chId[iCluster1];
1395         isReconstructible = kTRUE;
1396         return firstLabel;
1397       }
1398       
1399     }
1400     
1401   }
1402   
1403   return decayLabel;
1404 }
1405
1406 //________________________________________________________________________
1407 void AliAnalysisTaskMuonFakes::AddCompatibleClusters(const AliMUONTrack &track, const AliMUONTrack &trackRef,
1408                                                      TArrayI *labels, Int_t *nLabels) const
1409 {
1410   /// Try to match clusters between track and trackRef and add the corresponding MC labels to the arrays
1411   
1412   Double_t chi2Max = 2. * fSigmaCut * fSigmaCut; // 2 because 2 quantities in chi2
1413   
1414   // Loop over clusters of first track
1415   Int_t nCl1 = track.GetNClusters();
1416   for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) {
1417     AliMUONVCluster *cluster1 = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCl1))->GetClusterPtr();
1418     
1419     // Loop over clusters of second track
1420     Int_t nCl2 = trackRef.GetNClusters();
1421     for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) {
1422       AliMUONVCluster *cluster2 = static_cast<AliMUONTrackParam*>(trackRef.GetTrackParamAtCluster()->UncheckedAt(iCl2))->GetClusterPtr();
1423       
1424       // check DE Id
1425       if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue;
1426       
1427       // check local chi2
1428       Double_t dX = cluster1->GetX() - cluster2->GetX();
1429       Double_t dY = cluster1->GetY() - cluster2->GetY();
1430       Double_t chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
1431       if (chi2 > chi2Max) continue;
1432       
1433       // expand array if needed
1434       if (nLabels[iCl1] >= labels[iCl1].GetSize()) labels[iCl1].Set(2*labels[iCl1].GetSize());
1435       
1436       // save label
1437       labels[iCl1][nLabels[iCl1]] = static_cast<Int_t>(trackRef.GetUniqueID());
1438       nLabels[iCl1]++;
1439       break;
1440       
1441     }
1442     
1443   }
1444   
1445 }
1446
1447 //________________________________________________________________________
1448 Int_t AliAnalysisTaskMuonFakes::IsDecayByLabel(const AliMUONTrack &track, Bool_t &isReconstructible,
1449                                                Int_t &lastCh) const
1450 {
1451   /// Check whether this track correspond to a decaying particle by using cluster MC labels.
1452   /// "lastCh" contains the chamber Id of the most downstream chamber hit by the decay chain
1453   if (fPrintDecayChain) printf("\nBY LABEL\n");
1454   
1455   Int_t nClusters = track.GetNClusters();
1456   if (nClusters <= 0) return -1;
1457   Int_t *chId = new Int_t[nClusters];
1458   Int_t *labels = new Int_t[nClusters];
1459   
1460   // copy labels and chamber Ids
1461   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
1462     AliMUONVCluster* cluster = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1463     chId[iCluster] = cluster->GetChamberId();
1464     labels[iCluster] = cluster->GetMCLabel();
1465     if (fPrintDecayChain) {
1466       printf("ch%d: %d",chId[iCluster],labels[iCluster]);
1467       Int_t currentLabel = labels[iCluster];
1468       while (currentLabel >= 0) {
1469         AliMCParticle* currentParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(currentLabel));
1470         printf("(%s)",(currentParticle) ? currentParticle->Particle()->GetName() : "");
1471         if (currentLabel == labels[iCluster]) printf(" (");
1472         currentLabel = (currentParticle) ? currentParticle->GetMother() : -1;
1473         if (currentLabel >= 0) printf(" %d",currentLabel);
1474       }
1475       printf(" )\n");
1476     }
1477   }
1478   
1479   // look for decay
1480   lastCh = 0;
1481   Int_t decayLabel = IsDecay(nClusters, chId, labels, isReconstructible, lastCh);
1482   if (fPrintDecayChain) printf("---> decayLabel = %d (reco = %d / lastCh = %d)\n",decayLabel,isReconstructible,lastCh);
1483   
1484   delete[] chId;
1485   delete[] labels;
1486   
1487   return decayLabel;
1488 }
1489
1490 //________________________________________________________________________
1491 Int_t AliAnalysisTaskMuonFakes::IsDecayByPosition(const AliMUONTrack &track, const AliMUONVTrackStore &trackRefStore,
1492                                                   const AliMUONVTrackStore &usedTrackRefStore, Bool_t &isReconstructible,
1493                                                   Int_t &lastCh) const
1494 {
1495   /// Check whether this track correspond to a decaying particle by comparing clusters position
1496   /// All possible combinations of compatible clusters from every trackRefs are considered
1497   if (fPrintDecayChain) printf("\nBY POSITION\n");
1498   
1499   Int_t nClusters = track.GetNClusters();
1500   if (nClusters <= 0) return -1;
1501   Int_t *chId = new Int_t[nClusters];
1502   Int_t *nLabels = new Int_t[nClusters];
1503   TArrayI *labels = new TArrayI[nClusters];
1504   
1505   // copy chamber Ids
1506   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
1507     AliMUONVCluster* cluster = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr();
1508     chId[iCluster] = cluster->GetChamberId();
1509     nLabels[iCluster] = 0;
1510     labels[iCluster].Set(100);
1511   }
1512   
1513   // loop over trackRef store and add label of compatible clusters
1514   TIter next1(trackRefStore.CreateIterator());
1515   AliMUONTrack* trackRef;
1516   while ( ( trackRef = static_cast<AliMUONTrack*>(next1()) ) )
1517     AddCompatibleClusters(track, *trackRef, labels, nLabels);
1518   
1519   // loop over usedTrackRef store and add label of compatible clusters
1520   TIter next2(usedTrackRefStore.CreateIterator());
1521   while ( ( trackRef = static_cast<AliMUONTrack*>(next2()) ) )
1522     AddCompatibleClusters(track, *trackRef, labels, nLabels);
1523   
1524   // complete the arrays of labels with "-1" if no label was found for a given cluster
1525   for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
1526     if (nLabels[iCluster] == 0) {
1527       labels[iCluster][0] = -1;
1528       nLabels[iCluster]++;
1529     }
1530   }
1531   
1532   // loop over all possible combinations
1533   Int_t *iLabel = new Int_t[nClusters];
1534   memset(iLabel,0,nClusters*sizeof(Int_t));
1535   iLabel[nClusters-1] = -1;
1536   Int_t *currentLabels = new Int_t[nClusters];
1537   Int_t decayLabel = -1;
1538   lastCh = 0;
1539   isReconstructible = kFALSE;
1540   while (kTRUE) {
1541     
1542     // go to the next combination
1543     Int_t iCl = nClusters-1;
1544     while (++iLabel[iCl] >= nLabels[iCl] && iCl > 0) iLabel[iCl--] = 0;
1545     if (iLabel[iCl] >= nLabels[iCl]) break; // no more combination
1546     
1547     // copy labels
1548     if (fPrintDecayChain) printf("\n");
1549     for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
1550       currentLabels[iCluster] = labels[iCluster][iLabel[iCluster]];
1551       if (fPrintDecayChain) {
1552         printf("ch%d: %d",chId[iCluster],currentLabels[iCluster]);
1553         Int_t currentLabel = currentLabels[iCluster];
1554         while (currentLabel >= 0) {
1555           AliMCParticle* currentParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(currentLabel));
1556           printf("(%s)",(currentParticle) ? currentParticle->Particle()->GetName() : "");
1557           if (currentLabel == currentLabels[iCluster]) printf(" (");
1558           currentLabel = (currentParticle) ? currentParticle->GetMother() : -1;
1559           if (currentLabel >= 0) printf(" %d",currentLabel);
1560         }
1561         printf(" )\n");
1562       }
1563     }
1564     
1565     // look for decay
1566     Int_t currentDecayLabel = IsDecay(nClusters, chId, currentLabels, isReconstructible, lastCh);
1567     if (fPrintDecayChain) printf("---> decayLabel = %d (reco = %d / lastCh = %d)\n",currentDecayLabel,isReconstructible,lastCh);
1568     if (currentDecayLabel >= 0) {
1569       decayLabel = currentDecayLabel;
1570       if (isReconstructible) break;
1571     }
1572     
1573   }
1574   
1575   if (fPrintDecayChain) printf("------> decayLabel = %d (reco = %d / lastCh = %d)\n",decayLabel,isReconstructible,lastCh);
1576   
1577   delete[] chId;
1578   delete[] nLabels;
1579   delete[] labels;
1580   delete[] iLabel;
1581   delete[] currentLabels;
1582   
1583   return decayLabel;  
1584 }
1585