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