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