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