]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/MUON/dep/AliAnalysisTaskMuonFakes.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGPP / MUON / dep / AliAnalysisTaskMuonFakes.cxx
CommitLineData
fa97374c 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
27de2dfb 16/* $Id$ */
17
fa97374c 18// ROOT includes
19#include <TObjArray.h>
20#include <TClonesArray.h>
21#include <TH1F.h>
22#include <TH2F.h>
23#include <TCanvas.h>
66c5256f 24#include <TArrayI.h>
fa97374c 25
26// STEER includes
27#include "AliLog.h"
28#include "AliESDEvent.h"
29#include "AliESDMuonTrack.h"
30#include "AliESDInputHandler.h"
31#include "AliMCEventHandler.h"
0ffd1c5a 32#include "AliCDBManager.h"
33#include "AliCentrality.h"
66c5256f 34#include "AliMCParticle.h"
35#include "AliMCEvent.h"
36#include "AliCounterCollection.h"
fa97374c 37
38// ANALYSIS includes
c121b0eb 39#include "AliAnalysisDataSlot.h"
40#include "AliAnalysisDataContainer.h"
fa97374c 41#include "AliAnalysisManager.h"
fa97374c 42
43// MUON includes
44#include "AliMUONCDB.h"
45#include "AliMUONRecoParam.h"
46#include "AliMUONRecoCheck.h"
47#include "AliMUONVCluster.h"
48#include "AliMUONVTrackStore.h"
c121b0eb 49#include "AliMUONVTriggerTrackStore.h"
fa97374c 50#include "AliMUONTrack.h"
51#include "AliMUONTrackParam.h"
52#include "AliMUONESDInterface.h"
66c5256f 53#include "AliMUONTriggerTrack.h"
fa97374c 54
55#include "AliAnalysisTaskMuonFakes.h"
fa97374c 56
e8c28d9e 57using std::cout;
58using std::endl;
59using std::flush;
60
fa97374c 61ClassImp(AliAnalysisTaskMuonFakes)
62
63//________________________________________________________________________
64AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes() :
65AliAnalysisTaskSE(),
66fList(0x0),
0ffd1c5a 67fList2(0x0),
fa97374c 68fCanvases(0x0),
69fTrackCounters(0x0),
70fFakeTrackCounters(0x0),
71fMatchedTrackCounters(0x0),
72fEventCounters(0x0),
0ffd1c5a 73fPairCounters(0x0),
fa97374c 74fCurrentFileName(""),
75fRequestedStationMask(0),
76fRequest2ChInSameSt45(kFALSE),
77fSigmaCut(-1.),
66c5256f 78fNEvents(0),
79fShowProgressBar(kFALSE),
fa97374c 80fUseLabel(kFALSE),
66c5256f 81fCombineMCId(kFALSE),
82fExternalSigmaCut(-1.),
0ffd1c5a 83fMatchTrig(kFALSE),
84fApplyAccCut(kFALSE),
66c5256f 85fChi2Cut(-1.),
86fPtCut(-1.),
87fRecoParamLocation(""),
88fDecayAsFake(kFALSE),
89fPrintDecayChain(kFALSE),
90fDisableDetailedCounters(kFALSE),
91fMuonTrackCuts(0x0)
fa97374c 92{
93 /// Default constructor.
94}
95
96//________________________________________________________________________
97AliAnalysisTaskMuonFakes::AliAnalysisTaskMuonFakes(const char *name) :
98AliAnalysisTaskSE(name),
99fList(0x0),
0ffd1c5a 100fList2(0x0),
fa97374c 101fCanvases(0x0),
102fTrackCounters(0x0),
103fFakeTrackCounters(0x0),
104fMatchedTrackCounters(0x0),
105fEventCounters(0x0),
0ffd1c5a 106fPairCounters(0x0),
fa97374c 107fCurrentFileName(""),
108fRequestedStationMask(0),
109fRequest2ChInSameSt45(kFALSE),
110fSigmaCut(-1.),
66c5256f 111fNEvents(0),
112fShowProgressBar(kFALSE),
fa97374c 113fUseLabel(kFALSE),
66c5256f 114fCombineMCId(kFALSE),
115fExternalSigmaCut(-1.),
0ffd1c5a 116fMatchTrig(kFALSE),
117fApplyAccCut(kFALSE),
66c5256f 118fChi2Cut(-1.),
119fPtCut(-1.),
120fRecoParamLocation("raw://"),
121fDecayAsFake(kFALSE),
122fPrintDecayChain(kFALSE),
123fDisableDetailedCounters(kFALSE),
124fMuonTrackCuts(0x0)
fa97374c 125{
126 /// Constructor.
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());
0ffd1c5a 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());
fa97374c 141}
142
143//________________________________________________________________________
144AliAnalysisTaskMuonFakes::~AliAnalysisTaskMuonFakes()
145{
146 /// Destructor.
147 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
148 delete fList;
0ffd1c5a 149 delete fList2;
fa97374c 150 delete fTrackCounters;
151 delete fFakeTrackCounters;
152 delete fMatchedTrackCounters;
153 delete fEventCounters;
0ffd1c5a 154 delete fPairCounters;
fa97374c 155 }
156 delete fCanvases;
66c5256f 157 delete fMuonTrackCuts;
fa97374c 158}
159
160//___________________________________________________________________________
161void AliAnalysisTaskMuonFakes::UserCreateOutputObjects()
162{
163 /// Create histograms and counters.
164
0ffd1c5a 165 // single track histograms
fa97374c 166 fList = new TObjArray(100);
167 fList->SetOwner();
168
66c5256f 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
fa97374c 234 // number of tracks
235 TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks", "nb of tracks /evt", 21, -0.5, 20.5);
66c5256f 236 fList->AddAtAndExpand(hNumberOfTracks, 2*kNhistTrack+kNumberOfTracks);
fa97374c 237 TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks", "nb of fake - nb of missing track", 21, -0.5, 20.5);
66c5256f 238 fList->AddAtAndExpand(hNumberOfAdditionalTracks, 2*kNhistTrack+kNumberOfAdditionalTracks);
239
240 // number of clusters MC / fraction of clusters
fa97374c 241 TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC", "nb of clusters /MC track", 21, -0.5, 20.5);
66c5256f 242 fList->AddAtAndExpand(hNumberOfClustersMC, 2*kNhistTrack+kNumberOfClustersMC);
fa97374c 243 TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters", "nb of matched clusters / nb of clusters", 110, 0., 1.1);
66c5256f 244 fList->AddAtAndExpand(hFractionOfMatchedClusters, 2*kNhistTrack+kFractionOfMatchedClusters);
fa97374c 245 TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters", "nb of connected clusters / nb of clusters in fake tracks", 110, 0., 1.1);
66c5256f 246 fList->AddAtAndExpand(hFractionOfConnectedClusters, 2*kNhistTrack+kFractionOfConnectedClusters);
0ffd1c5a 247
248 // track pair histograms
249 fList2 = new TObjArray(100);
250 fList2->SetOwner();
251
252 // physics quantities of opposite-sign track pairs
0ffd1c5a 253 TString nameSuffix[4] = {"", "M", "F1", "F2"};
66c5256f 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 }
0ffd1c5a 276 }
fa97374c 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
66c5256f 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
fa97374c 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
c121b0eb 288 fTrackCounters = new AliCounterCollection(GetOutputSlot(2)->GetContainer()->GetName());
66c5256f 289 fTrackCounters->AddRubric("track", "reconstructible/reconstructed/matched/matchedyet/decay/decayyet/fake/connected/additional");
fa97374c 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");
0ffd1c5a 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());
fa97374c 296 fTrackCounters->Init();
297
66c5256f 298 // detailled counters of decays and fake tracks:
c121b0eb 299 fFakeTrackCounters = new AliCounterCollection(GetOutputSlot(3)->GetContainer()->GetName());
66c5256f 300 fFakeTrackCounters->AddRubric("position", "matched/decay/decayyet/matchedyet/fake/connected/additional");
301 fFakeTrackCounters->AddRubric("label", "matched/decay/decayyet/matchedyet/fake/connected/additional");
fa97374c 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");
0ffd1c5a 308 fFakeTrackCounters->AddRubric("cent", centralityClasses.Data());
fa97374c 309 fFakeTrackCounters->Init();
310
311 // counters of tracks matched by position or by using MC labels
c121b0eb 312 fMatchedTrackCounters = new AliCounterCollection(GetOutputSlot(4)->GetContainer()->GetName());
66c5256f 313 fMatchedTrackCounters->AddRubric("position", "matched/decay/decayyet/matchedyet/fake");
314 fMatchedTrackCounters->AddRubric("label", "matched/decay/decayyet/matchedyet/fake/matchedother");
fa97374c 315 fMatchedTrackCounters->AddRubric("run", 1000000);
316 fMatchedTrackCounters->AddRubric("trig", "yes/no");
317 fMatchedTrackCounters->AddRubric("selected", "yes/no");
318 fMatchedTrackCounters->AddRubric("acc", "in/out");
0ffd1c5a 319 fMatchedTrackCounters->AddRubric("cent", centralityClasses.Data());
fa97374c 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
c121b0eb 329 fEventCounters = new AliCounterCollection(GetOutputSlot(5)->GetContainer()->GetName());
fa97374c 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");
0ffd1c5a 334 fEventCounters->AddRubric("cent", centralityClasses.Data());
fa97374c 335 fEventCounters->Init();
336
0ffd1c5a 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
fa97374c 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);
0ffd1c5a 361 PostData(6, fList2);
362 PostData(7, fPairCounters);
fa97374c 363}
364
365//________________________________________________________________________
366void AliAnalysisTaskMuonFakes::UserExec(Option_t *)
367{
368 /// Process event: looks for fakes...
369
370 // check that reconstructions parameters for that run have been properly set
371 if (fSigmaCut < 0) return;
372
66c5256f 373 if (fShowProgressBar && (++fNEvents)%100 == 0) cout<<"\rEvent processing... "<<fNEvents<<"\r"<<flush;
374
fa97374c 375 // check physics selection
376 TString selected = (fInputHandler && fInputHandler->IsEventSelected() != 0) ? "selected:yes" : "selected:no";
377
378 // current file name
66c5256f 379 if (fDisableDetailedCounters) fCurrentFileName = "any";
380 else {
381 fCurrentFileName = CurrentFileName();
382 fCurrentFileName.ReplaceAll("alien://","");
383 fCurrentFileName.ReplaceAll("/","\\");
384 fCurrentFileName.ReplaceAll(":",";");
385 }
fa97374c 386
387 // Load ESD event
388 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
2426f648 389 if (!esd) {
390 AliError("Cannot get input event");
391 return;
392 }
fa97374c 393
66c5256f 394 // event number in current file
395 TString eventNumberInFile = (fDisableDetailedCounters) ? "event:any" : Form("event:%d",esd->GetEventNumberInFile());
396
0ffd1c5a 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
fa97374c 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);
c121b0eb 420 AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
fa97374c 421 if (!muonTrackStore || !trackRefStore) return;
422
0ffd1c5a 423 // loop over trackRefs
424 Int_t nMuPlus[2] = {0, 0};
425 Int_t nMuMinus[2] = {0, 0};
fa97374c 426 TIter next(trackRefStore->CreateIterator());
427 AliMUONTrack* trackRef;
428 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
0ffd1c5a 429
430 // skip trackRefs that are not reconstructible
431 if (!trackRef->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) continue;
432
433 // trigger condition
66c5256f 434 AliMUONTriggerTrack *trigRef = static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(trackRef->GetUniqueID()));
435 Bool_t trigger = (trigRef && trigRef->GetPtCutLevel() > 0);
0ffd1c5a 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
fa97374c 446 }
447
0ffd1c5a 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
fa97374c 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};
66c5256f 458 AliMUONVTrackStore *usedTrackRefStore = AliMUONESDInterface::NewTrackStore();
fa97374c 459 AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
66c5256f 460 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
461 TArrayI mcLabels(nTracks);
0ffd1c5a 462 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
fa97374c 463
464 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
465
466 // skip ghosts
66c5256f 467 if (!IsSelected(*esdTrack)) continue;
fa97374c 468 containTrack[0] = kTRUE;
469
470 // trigger condition
471 Bool_t trigger = esdTrack->ContainTriggerData();
0ffd1c5a 472 TString trig = trigger ? "trig:yes" : "trig:no";
473 if (trigger) containTrack[1] = kTRUE;
fa97374c 474
475 // acceptance condition
0ffd1c5a 476 Double_t rAbs = esdTrack->GetRAtAbsorberEnd();
477 Double_t thetaTrackAbsEnd = TMath::ATan(rAbs/505.) * TMath::RadToDeg();
fa97374c 478 Double_t eta = esdTrack->Eta();
0ffd1c5a 479 Bool_t inAcc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5);
480 TString acc = inAcc ? "acc:in" : "acc:out";
fa97374c 481
482 // fill global counters
0ffd1c5a 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()));
fa97374c 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();
66c5256f 498 Double_t pU = esdTrack->PUncorrected();
499 Double_t pdca = 0.5*(p+pU)*dca;
fa97374c 500
501 // fill global histograms
66c5256f 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);
fa97374c 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);
66c5256f 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);
fa97374c 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);
66c5256f 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);
fa97374c 531
532 // fill global counters
66c5256f 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";
fa97374c 539 TString labelCase = "label:";
66c5256f 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
fa97374c 553 // take actions according to the matching result we are interested in
66c5256f 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
fa97374c 622 if (matchedTrackRef) {
623
624 // fill global counters
0ffd1c5a 625 fTrackCounters->Count(Form("track:matched/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
fa97374c 626
627 // track matched with a trackRef that is not reconstructible
66c5256f 628 if (isMatchedYet) {
fa97374c 629
630 containMatchedYetTrack[0] = kTRUE;
631 if (trigger) containMatchedYetTrack[1] = kTRUE;
632
633 // fill global counters
0ffd1c5a 634 fTrackCounters->Count(Form("track:matchedyet/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
66c5256f 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
fa97374c 641 }
642
643 // fill histograms
66c5256f 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);
0ffd1c5a 649
650 // flag matched tracks
66c5256f 651 mcLabels[iTrack] = matchedTrackRef->GetUniqueID();
fa97374c 652
66c5256f 653 // move already matched trackRefs
654 usedTrackRefStore->Add(*matchedTrackRef);
fa97374c 655 trackRefStore->Remove(*matchedTrackRef);
656
657 } else {
658
66c5256f 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
0ffd1c5a 685 }
686
66c5256f 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 }
fa97374c 707
708 }
709
710 } // end of loop over ESD tracks
711
712 // fill histogram and global counters
66c5256f 713 ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kNumberOfTracks))->Fill(nTrackerTracks);
0ffd1c5a 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()));
fa97374c 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
0ffd1c5a 725 Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, selected, centrality);
fa97374c 726
727 if (fakeTrackStore->GetSize() > 0) {
728
729 // fill global counters
0ffd1c5a 730 fEventCounters->Count(Form("event:notconnected/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
fa97374c 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
0ffd1c5a 743 if (containMatchedFake) fEventCounters->Count(Form("event:notconnected/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
fa97374c 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
66c5256f 751 ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kNumberOfAdditionalTracks))->Fill(nAdditionalTracks);
0ffd1c5a 752 fEventCounters->Count(Form("event:additional/run:%d/trig:any/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
fa97374c 753 if (!containUnmatchedFake) { // all matched
0ffd1c5a 754 fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
66c5256f 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);
fa97374c 758 } else if (!containMatchedFake) { // none matched
0ffd1c5a 759 fTrackCounters->Count(Form("track:additional/run:%d/trig:no/%s/acc:unknown/%s", fCurrentRunNumber, selected.Data(), centrality.Data()), nAdditionalTracks);
66c5256f 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);
fa97374c 762 } else { // mixed
0ffd1c5a 763 fEventCounters->Count(Form("event:additional/run:%d/trig:yes/%s/%s", fCurrentRunNumber, selected.Data(), centrality.Data()));
66c5256f 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);
fa97374c 767 }
768
769 }
770
771 }
772
773 }
774
0ffd1c5a 775 // clean memory
66c5256f 776 delete usedTrackRefStore;
fa97374c 777 delete fakeTrackStore;
778
0ffd1c5a 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
66c5256f 785 if (!IsSelected(*muonTrack1)) continue;
0ffd1c5a 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();
66c5256f 793 Int_t label1 = mcLabels[iTrack1];
0ffd1c5a 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
66c5256f 800 if (!IsSelected(*muonTrack2)) continue;
0ffd1c5a 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);
66c5256f 811 Int_t label2 = mcLabels[iTrack2];
0ffd1c5a 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
66c5256f 831 Bool_t inAcc = (acc1 && acc2 && y >= -4. && y <= -2.5);
0ffd1c5a 832 TString acc = inAcc ? "acc:in" : "acc:out";
833
834 // fill global histograms
66c5256f 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);
0ffd1c5a 838
839 TString pair = "pair:";
840
841 // fill histograms according to labels
842 if (label1 >= 0 && label2 >= 0) {
843
844 pair += "matched";
845
66c5256f 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);
0ffd1c5a 849
850 } else if (label1 >= 0 || label2 >= 0) {
851
852 pair += "1fake";
853
66c5256f 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);
0ffd1c5a 857
858 } else {
859
860 pair += "2fakes";
861
66c5256f 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);
0ffd1c5a 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
fa97374c 876 // Post final data
877 PostData(1, fList);
878 PostData(2, fTrackCounters);
879 PostData(3, fFakeTrackCounters);
880 PostData(4, fMatchedTrackCounters);
881 PostData(5, fEventCounters);
0ffd1c5a 882 PostData(6, fList2);
883 PostData(7, fPairCounters);
fa97374c 884}
885
886//________________________________________________________________________
887void AliAnalysisTaskMuonFakes::NotifyRun()
888{
889 /// Prepare processing of new run: load corresponding OCDB objects...
890
66c5256f 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
fa97374c 901 // load necessary data from OCDB
fa97374c 902 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
903 if (!recoParam) {
904 fRequestedStationMask = 0;
905 fRequest2ChInSameSt45 = kFALSE;
906 fSigmaCut = -1.;
907 AliError("--> skip this run");
908 return;
909 }
910
911 // compute the mask of requested stations from recoParam
912 fRequestedStationMask = 0;
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
66c5256f 918 // get sigma cut to associate clusters with TrackRefs from recoParam if not already set manually
919 if (fExternalSigmaCut > 0) fSigmaCut = fExternalSigmaCut;
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);
fa97374c 925}
926
927//________________________________________________________________________
928void AliAnalysisTaskMuonFakes::Terminate(Option_t *)
929{
930 /// Draw results to the screen and print statistics.
931
932 // recover output objects
933 fList = static_cast<TObjArray*> (GetOutputData(1));
0ffd1c5a 934 fList2 = static_cast<TObjArray*> (GetOutputData(6));
935 if (!fList || !fList2) return;
fa97374c 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));
0ffd1c5a 940 fPairCounters = static_cast<AliCounterCollection*> (GetOutputData(7));
fa97374c 941
66c5256f 942 TString extention = GetName();
943 extention.ReplaceAll("MUONFakes","");
944
fa97374c 945 // add canvas to compare histograms
66c5256f 946 fCanvases = new TObjArray(13);
fa97374c 947 fCanvases->SetOwner();
fa97374c 948
66c5256f 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
976 Int_t iHist1[9] = {kNumberOfClusters, kNumberOfChamberHit, kChi2PerDof, kDCA, kRAbs, kEta, kP, kPt, kPhi};
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
0ffd1c5a 1120 }
1121
fa97374c 1122 // print
1123 if (fTrackCounters && fFakeTrackCounters && fMatchedTrackCounters && fEventCounters) {
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");
66c5256f 1127 fFakeTrackCounters->Print("label/position/trig");
fa97374c 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
0ffd1c5a 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");
fa97374c 1140}
1141
66c5256f 1142
1143//________________________________________________________________________
1144Bool_t AliAnalysisTaskMuonFakes::IsSelected(AliESDMuonTrack &esdTrack)
1145{
1146 /// return kTRUE if the track pass the section criteria
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//________________________________________________________________________
1165void 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{
1169 /// fill global histograms at track level
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//________________________________________________________________________
1186void 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{
1189 /// fill global histograms at pair level
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
fa97374c 1198//________________________________________________________________________
0ffd1c5a 1199Int_t AliAnalysisTaskMuonFakes::RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
1200 TString &selected, TString &centrality)
fa97374c 1201{
1202 /// loop over reconstructible TrackRef not associated with reconstructed track:
1203 /// for each of them, find and remove the most connected the fake track, if any,
1204 /// and fill the histograms with the fraction of connected clusters.
1205 /// Return the number of reconstructible track not connected to any fake
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;
66c5256f 1228 if (fUseLabel || fCombineMCId) { // by using the MC label
fa97374c 1229 for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
1230 if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
1231 nConnectedClusters++;
66c5256f 1232 }
1233 if (!fUseLabel || fCombineMCId) { // by comparing cluster/TrackRef positions
fa97374c 1234 Bool_t compTrack[10];
66c5256f 1235 nConnectedClusters = TMath::Max(nConnectedClusters, fakeTrack->FindCompatibleClusters(*trackRef, fSigmaCut, compTrack));
fa97374c 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());
2426f648 1254 if (!esd) {
1255 AliError("Cannot get input event");
1256 return nFreeMissingTracks;
1257 }
fa97374c 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();
0ffd1c5a 1272 TString acc = (thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5) ? "acc:in" : "acc:out";
fa97374c 1273
1274 // fill histogram and counters
66c5256f 1275 TString eventNumberInFile = (fDisableDetailedCounters) ? "event:any" : Form("event:%d",esd->GetEventNumberInFile());
1276 ((TH1F*)fList->UncheckedAt(2*kNhistTrack+kFractionOfConnectedClusters))->Fill(fractionOfConnectedClusters);
0ffd1c5a 1277 fTrackCounters->Count(Form("track:connected/run:%d/%s/%s/%s/%s", fCurrentRunNumber, trig.Data(), selected.Data(), acc.Data(), centrality.Data()));
66c5256f 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()));
fa97374c 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
66c5256f 1292//________________________________________________________________________
1293Int_t AliAnalysisTaskMuonFakes::IsDecay(Int_t nClusters, Int_t *chId, Int_t *labels,
1294 Bool_t &isReconstructible, Int_t &lastCh) const
1295{
1296 /// Check whether this combination of clusters correspond to a decaying particle or not:
1297 /// More than 50% of clusters, including 1 before and 1 after the dipole, must be connected.
1298 /// - Return the MC label of the most downstream decay product or -1 if not a decay.
1299 /// - "isReconstructible" tells if the combination of matched clusters fulfil the reconstruction criteria.
1300 /// - As soon as we realized the decay chain cannot be tagged as reconstructible, we reject any chain ending
1301 /// on a chamber equal to or upstream "lastCh" (used to select the best chain in case of multiple choices).
1302 /// - "lastCh" is reset the most downstream chamber of the found decay chain if any.
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//________________________________________________________________________
1411void AliAnalysisTaskMuonFakes::AddCompatibleClusters(const AliMUONTrack &track, const AliMUONTrack &trackRef,
1412 TArrayI *labels, Int_t *nLabels) const
1413{
1414 /// Try to match clusters between track and trackRef and add the corresponding MC labels to the arrays
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//________________________________________________________________________
1452Int_t AliAnalysisTaskMuonFakes::IsDecayByLabel(const AliMUONTrack &track, Bool_t &isReconstructible,
1453 Int_t &lastCh) const
1454{
1455 /// Check whether this track correspond to a decaying particle by using cluster MC labels.
1456 /// "lastCh" contains the chamber Id of the most downstream chamber hit by the decay chain
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//________________________________________________________________________
1495Int_t AliAnalysisTaskMuonFakes::IsDecayByPosition(const AliMUONTrack &track, const AliMUONVTrackStore &trackRefStore,
1496 const AliMUONVTrackStore &usedTrackRefStore, Bool_t &isReconstructible,
1497 Int_t &lastCh) const
1498{
1499 /// Check whether this track correspond to a decaying particle by comparing clusters position
1500 /// All possible combinations of compatible clusters from every trackRefs are considered
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