]>
Commit | Line | Data |
---|---|---|
63c5d0a6 | 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 | ||
e5b66c63 | 18 | // ROOT includes |
19 | #include <TArrayI.h> | |
20 | #include <TClonesArray.h> | |
21 | ||
b54df4f9 | 22 | // STEER includes |
63c5d0a6 | 23 | #include "AliESDEvent.h" |
63c5d0a6 | 24 | #include "AliESDMuonTrack.h" |
63c5d0a6 | 25 | #include "AliLog.h" |
63c5d0a6 | 26 | #include "AliMCEventHandler.h" |
b54df4f9 | 27 | #include "AliCDBManager.h" |
e5b66c63 | 28 | #include "AliMCParticle.h" |
29 | #include "AliMCEvent.h" | |
b54df4f9 | 30 | |
31 | // ANALYSIS includes | |
32 | #include "AliAnalysisManager.h" | |
63c5d0a6 | 33 | |
b54df4f9 | 34 | // MUON includes |
35 | #include "AliMUONCDB.h" | |
36 | #include "AliMUONRecoParam.h" | |
63c5d0a6 | 37 | #include "AliMUONRecoCheck.h" |
38 | #include "AliMUONESDInterface.h" | |
e5b66c63 | 39 | #include "AliMUONVTrackReconstructor.h" |
63c5d0a6 | 40 | #include "AliMUONTrack.h" |
e5b66c63 | 41 | #include "AliMUONTrackParam.h" |
63c5d0a6 | 42 | #include "AliMUONVTrackStore.h" |
b54df4f9 | 43 | #include "AliMUONTriggerTrack.h" |
44 | #include "AliMUONVTriggerTrackStore.h" | |
45 | #include "AliMUONLocalTrigger.h" | |
e5b66c63 | 46 | #include "AliMUONVCluster.h" |
63c5d0a6 | 47 | |
b54df4f9 | 48 | #include "AliAnalysisTaskESDMCLabelAddition.h" |
63c5d0a6 | 49 | |
b54df4f9 | 50 | ClassImp(AliAnalysisTaskESDMCLabelAddition) |
63c5d0a6 | 51 | |
52 | //---------------------------------------------------------------------- | |
53 | AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition(): | |
b54df4f9 | 54 | AliAnalysisTaskSE(), |
55 | fDefaultStorage(""), | |
e5b66c63 | 56 | fRequestedStationMask(0), |
57 | fRequest2ChInSameSt45(kFALSE), | |
58 | fExternalTrkSigmaCut(-1.), | |
b54df4f9 | 59 | fSigmaCut(-1.), |
e5b66c63 | 60 | fExternalTrgSigmaCut(-1.), |
b54df4f9 | 61 | fSigmaCutTrig(-1.) |
63c5d0a6 | 62 | { |
b54df4f9 | 63 | /// Default constructor |
63c5d0a6 | 64 | } |
65 | ||
66 | ||
67 | //---------------------------------------------------------------------- | |
68 | AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition(const char* name): | |
b54df4f9 | 69 | AliAnalysisTaskSE(name), |
70 | fDefaultStorage("raw://"), | |
e5b66c63 | 71 | fRequestedStationMask(0), |
72 | fRequest2ChInSameSt45(kFALSE), | |
73 | fExternalTrkSigmaCut(-1.), | |
b54df4f9 | 74 | fSigmaCut(-1.), |
e5b66c63 | 75 | fExternalTrgSigmaCut(-1.), |
b54df4f9 | 76 | fSigmaCutTrig(-1.) |
63c5d0a6 | 77 | { |
b54df4f9 | 78 | /// Constructor |
63c5d0a6 | 79 | } |
80 | ||
81 | ||
82 | //---------------------------------------------------------------------- | |
83 | void AliAnalysisTaskESDMCLabelAddition::UserCreateOutputObjects() | |
84 | { | |
b54df4f9 | 85 | /// Create output objects |
63c5d0a6 | 86 | } |
87 | ||
88 | ||
89 | //---------------------------------------------------------------------- | |
b54df4f9 | 90 | void AliAnalysisTaskESDMCLabelAddition::NotifyRun() |
63c5d0a6 | 91 | { |
b54df4f9 | 92 | /// Load OCDB inputs |
93 | ||
e5b66c63 | 94 | // load OCDB objects only once |
95 | if (fSigmaCut > 0) return; | |
96 | ||
b54df4f9 | 97 | // set OCDB location |
98 | AliCDBManager* cdbm = AliCDBManager::Instance(); | |
e5b66c63 | 99 | if (cdbm->IsDefaultStorageSet()) printf("MCLabelAddition: CDB default storage already set!\n"); |
100 | else cdbm->SetDefaultStorage(fDefaultStorage.Data()); | |
101 | if (cdbm->GetRun() > -1) printf("MCLabelAddition: run number already set!\n"); | |
102 | else cdbm->SetRun(fCurrentRunNumber); | |
b54df4f9 | 103 | |
104 | // load recoParam | |
e5b66c63 | 105 | const AliMUONRecoParam* recoParam = (AliMUONESDInterface::GetTracker()) |
106 | ? AliMUONESDInterface::GetTracker()->GetRecoParam() | |
107 | : AliMUONCDB::LoadRecoParam(); | |
108 | ||
b54df4f9 | 109 | if (!recoParam) { |
e5b66c63 | 110 | fRequestedStationMask = 0; |
111 | fRequest2ChInSameSt45 = kFALSE; | |
b54df4f9 | 112 | fSigmaCut = -1.; |
113 | fSigmaCutTrig = -1.; | |
114 | return; | |
115 | } | |
116 | ||
e5b66c63 | 117 | // compute the mask of requested stations from recoParam |
118 | fRequestedStationMask = 0; | |
119 | for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i ); | |
120 | ||
121 | // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible | |
122 | fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates(); | |
b54df4f9 | 123 | |
e5b66c63 | 124 | // get sigma cut to associate clusters with TrackRefs from recoParam if not already set manually |
125 | if (fExternalTrkSigmaCut > 0) fSigmaCut = fExternalTrkSigmaCut; | |
126 | else if (recoParam->ImproveTracks()) fSigmaCut = recoParam->GetSigmaCutForImprovement(); | |
127 | else fSigmaCut = recoParam->GetSigmaCutForTracking(); | |
128 | ||
129 | // get sigma cut to associate trigger to triggerable track from recoParam if not already set manually | |
130 | if (fExternalTrgSigmaCut > 0) fSigmaCutTrig = fExternalTrgSigmaCut; | |
131 | else fSigmaCutTrig = recoParam->GetSigmaCutForTrigger(); | |
b54df4f9 | 132 | |
63c5d0a6 | 133 | } |
134 | ||
135 | ||
136 | //---------------------------------------------------------------------- | |
137 | void AliAnalysisTaskESDMCLabelAddition::UserExec(Option_t */*option*/) | |
138 | { | |
b54df4f9 | 139 | /// Execute analysis for current event |
140 | ||
141 | AliDebug(1, Form("MCLabel Addition: Analysing event # %5d\n",(Int_t) Entry())); | |
142 | ||
e5b66c63 | 143 | // make sure necessary information from OCDB have been loaded |
b54df4f9 | 144 | if (fSigmaCut < 0) return; |
145 | ||
146 | /// Load ESD event | |
63c5d0a6 | 147 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()); |
2426f648 | 148 | if (!esd) { |
149 | AliError("Cannot get input event"); | |
150 | return; | |
151 | } | |
63c5d0a6 | 152 | |
153 | // Load MC event | |
b54df4f9 | 154 | AliMCEventHandler* mcH = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); |
155 | if ( ! mcH ) { | |
156 | AliError ("MCH event handler not found. Nothing done!"); | |
157 | return; | |
158 | } | |
159 | ||
63c5d0a6 | 160 | |
161 | // Get reference tracks | |
162 | AliMUONRecoCheck rc(esd,mcH); | |
163 | AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1); | |
b54df4f9 | 164 | AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1); |
63c5d0a6 | 165 | |
166 | // Loop over reconstructed tracks | |
167 | AliESDMuonTrack *esdTrack = 0x0; | |
168 | Int_t nMuTracks = esd->GetNumberOfMuonTracks(); | |
169 | for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) { | |
170 | ||
171 | esdTrack = esd->GetMuonTrack(nMuTrack); | |
172 | ||
b54df4f9 | 173 | // tracker tracks |
174 | if (esdTrack->ContainTrackerData()) { | |
175 | ||
176 | // convert ESD track to MUON track (without recomputing track parameters at each clusters) | |
177 | AliMUONTrack muonTrack; | |
178 | AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE); | |
179 | ||
e5b66c63 | 180 | // try to match, by position, the reconstructed track with a simulated one |
181 | Int_t nMatchClustersByPosition = 0; | |
182 | AliMUONTrack* matchedTrackRefByPosition = rc.FindCompatibleTrack(muonTrack, *trackRefStore, nMatchClustersByPosition, kFALSE, fSigmaCut); | |
183 | Bool_t isMatchedYetByPosition = kFALSE; | |
184 | Bool_t isRecoDecayByPosition = kFALSE; | |
185 | Int_t decayLabelByPosition = -1, lastChDecayByPosition = 0; | |
186 | if (!matchedTrackRefByPosition || !matchedTrackRefByPosition->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) { | |
187 | decayLabelByPosition = IsDecayByPosition(muonTrack, *trackRefStore, isRecoDecayByPosition, lastChDecayByPosition); | |
188 | if (decayLabelByPosition >= 0) matchedTrackRefByPosition = 0x0; | |
189 | else if (matchedTrackRefByPosition) isMatchedYetByPosition = kTRUE; | |
190 | } | |
191 | Bool_t isFakeByPosition = (!matchedTrackRefByPosition && decayLabelByPosition < 0); | |
b54df4f9 | 192 | |
e5b66c63 | 193 | // try to match, by using MC labels, the reconstructed track with a simulated one |
194 | Int_t nMatchClustersByLabel = 0; | |
195 | AliMUONTrack* matchedTrackRefByLabel = rc.FindCompatibleTrack(muonTrack, *trackRefStore, nMatchClustersByLabel, kTRUE, fSigmaCut); | |
196 | Bool_t isMatchedYetByLabel = kFALSE; | |
197 | Bool_t isRecoDecayByLabel = kFALSE; | |
198 | Int_t decayLabelByLabel = -1, lastChDecayByLabel = 0; | |
199 | if (!matchedTrackRefByLabel || !matchedTrackRefByLabel->IsValid(fRequestedStationMask, fRequest2ChInSameSt45)) { | |
200 | decayLabelByLabel = IsDecayByLabel(muonTrack, isRecoDecayByLabel, lastChDecayByLabel); | |
201 | if (decayLabelByLabel >= 0) matchedTrackRefByLabel = 0x0; | |
202 | else if (matchedTrackRefByLabel) isMatchedYetByLabel = kTRUE; | |
203 | } | |
204 | Bool_t isFakeByLabel = (!matchedTrackRefByLabel && decayLabelByLabel < 0); | |
205 | ||
206 | // choose the best, or the only available, matched track | |
207 | AliMUONTrack* matchedTrackRef = 0x0; | |
208 | Bool_t isMatchedYet = kFALSE, isRecoDecay = kFALSE; | |
209 | Int_t decayLabel = -1; | |
210 | if (matchedTrackRefByPosition && matchedTrackRefByLabel && ((!isMatchedYetByPosition && !isMatchedYetByLabel) || | |
211 | (isMatchedYetByPosition && isMatchedYetByLabel))) { | |
212 | ||
213 | Int_t nMatchClusters = TMath::Max(nMatchClustersByPosition, nMatchClustersByLabel); | |
214 | matchedTrackRef = (nMatchClusters == nMatchClustersByPosition) ? matchedTrackRefByPosition : matchedTrackRefByLabel; | |
215 | isMatchedYet = isMatchedYetByPosition; | |
216 | ||
217 | } else if (matchedTrackRefByPosition && (!isMatchedYetByPosition || isFakeByLabel)) { | |
218 | ||
219 | matchedTrackRef = matchedTrackRefByPosition; | |
220 | isMatchedYet = isMatchedYetByPosition; | |
221 | ||
222 | } else if (matchedTrackRefByLabel && (!isMatchedYetByLabel || isFakeByPosition)) { | |
223 | ||
224 | matchedTrackRef = matchedTrackRefByLabel; | |
225 | isMatchedYet = isMatchedYetByLabel; | |
226 | ||
227 | // choose the best, or the only available, decay chain | |
228 | } else if (decayLabelByPosition >= 0 && decayLabelByLabel >= 0 && ((isRecoDecayByPosition && isRecoDecayByLabel) || | |
229 | (!isRecoDecayByPosition && !isRecoDecayByLabel))) { | |
230 | ||
231 | decayLabel = (lastChDecayByLabel > lastChDecayByPosition) ? decayLabelByLabel : decayLabelByPosition; | |
232 | isRecoDecay = isRecoDecayByPosition; | |
233 | ||
234 | } else if (decayLabelByPosition >= 0 && (isRecoDecayByPosition || decayLabelByLabel < 0)) { | |
235 | ||
236 | decayLabel = decayLabelByPosition; | |
237 | isRecoDecay = isRecoDecayByPosition; | |
238 | ||
239 | } else if (decayLabelByLabel >= 0) { | |
240 | ||
241 | decayLabel = decayLabelByLabel; | |
242 | isRecoDecay = isRecoDecayByLabel; | |
243 | ||
244 | } | |
245 | ||
246 | // set the MC label, the decay flag (bit 22) and the not-reconstructible flag (bit 23) | |
247 | if (matchedTrackRef) { | |
248 | ||
249 | esdTrack->SetLabel(matchedTrackRef->GetUniqueID()); | |
250 | esdTrack->SetBit(BIT(22), kFALSE); | |
251 | esdTrack->SetBit(BIT(23), isMatchedYet); | |
252 | ||
253 | } else if (decayLabel >= 0) { | |
254 | ||
255 | esdTrack->SetLabel(decayLabel); | |
256 | esdTrack->SetBit(BIT(22), kTRUE); | |
257 | esdTrack->SetBit(BIT(23), !isRecoDecay); | |
258 | ||
259 | } else { | |
260 | ||
261 | esdTrack->SetLabel(-1); | |
262 | esdTrack->SetBit(BIT(22), kFALSE); | |
263 | esdTrack->SetBit(BIT(23), kFALSE); | |
264 | ||
265 | } | |
b54df4f9 | 266 | |
267 | } else { // ghosts | |
268 | ||
269 | // Convert ESD track to trigger track | |
270 | AliMUONLocalTrigger locTrg; | |
271 | AliMUONESDInterface::ESDToMUON(*esdTrack, locTrg); | |
272 | AliMUONTriggerTrack trigTrack; | |
273 | rc.TriggerToTrack(locTrg, trigTrack); | |
274 | ||
275 | // try to match the reconstructed track with a simulated one | |
276 | AliMUONTriggerTrack* matchedTrigTrackRef = rc.FindCompatibleTrack(trigTrack, *triggerTrackRefStore, fSigmaCutTrig); | |
277 | ||
278 | // set the MC label | |
279 | if (matchedTrigTrackRef) esdTrack->SetLabel(matchedTrigTrackRef->GetUniqueID()); | |
280 | else esdTrack->SetLabel(-1); | |
281 | ||
282 | } | |
63c5d0a6 | 283 | |
284 | } | |
285 | ||
286 | } | |
287 | ||
288 | ||
289 | //---------------------------------------------------------------------- | |
290 | void AliAnalysisTaskESDMCLabelAddition::Terminate(Option_t */*option*/) | |
291 | { | |
b54df4f9 | 292 | /// Terminate analysis |
d5e5d942 | 293 | AliDebug(2, "Terminate()"); |
63c5d0a6 | 294 | } |
295 | ||
e5b66c63 | 296 | //---------------------------------------------------------------------- |
297 | Int_t AliAnalysisTaskESDMCLabelAddition::IsDecay(Int_t nClusters, Int_t *chId, Int_t *labels, | |
298 | Bool_t &isReconstructible, Int_t &lastCh) const | |
299 | { | |
300 | /// Check whether this combination of clusters correspond to a decaying particle or not: | |
301 | /// More than 50% of clusters, including 1 before and 1 after the dipole, must be connected. | |
302 | /// - Return the MC label of the most downstream decay product or -1 if not a decay. | |
303 | /// - "isReconstructible" tells if the combination of matched clusters fulfil the reconstruction criteria. | |
304 | /// - As soon as we realized the decay chain cannot be tagged as reconstructible, we reject any chain ending | |
305 | /// on a chamber equal to or upstream "lastCh" (used to select the best chain in case of multiple choices). | |
306 | /// - "lastCh" is reset the most downstream chamber of the found decay chain if any. | |
307 | ||
308 | Int_t halfCluster = nClusters/2; | |
309 | ||
310 | // loop over last clusters (if nClusters left < halfCluster the conditions cannot be fulfilled) | |
311 | Int_t firstLabel = -1, decayLabel = -1; | |
312 | isReconstructible = kFALSE; | |
313 | for (Int_t iCluster1 = nClusters-1; iCluster1 >= halfCluster; iCluster1--) { | |
314 | ||
315 | // if the last cluster is not on station 4 or 5 the conditions cannot be fulfilled | |
316 | if (chId[iCluster1] < 6) break; | |
317 | ||
318 | // skip clusters with no label or same label as at the begining of the previous step (already tested) | |
319 | if (labels[iCluster1] < 0 || labels[iCluster1] == firstLabel) continue; | |
320 | ||
321 | // is there any chance the hypothetical decay chain can be tagged reconstructible? | |
322 | Int_t stationId = chId[iCluster1]/2; | |
323 | Int_t stationMask = 1 << stationId; | |
324 | Int_t requestedStations = fRequestedStationMask >> stationId; | |
325 | Bool_t isValid = ((1 & requestedStations) == requestedStations); | |
326 | ||
327 | // if not: check whether we can find a better chain than already found | |
328 | if (!isValid && chId[iCluster1] <= lastCh) break; | |
329 | ||
330 | // count the number of fired chambers on stations 4 & 5 | |
331 | Int_t nChHitInSt45[2] = {0, 0}; | |
332 | nChHitInSt45[stationId-3] = 1; | |
333 | Int_t currentCh = chId[iCluster1]; | |
334 | ||
335 | // get the ancestors | |
336 | TArrayI chainLabels(100); | |
337 | Int_t nParticles = 0; | |
338 | Int_t currentLabel = labels[iCluster1]; | |
339 | do { | |
340 | chainLabels[nParticles++] = currentLabel; | |
341 | if (nParticles >= chainLabels.GetSize()) chainLabels.Set(2*chainLabels.GetSize()); | |
342 | AliMCParticle* currentParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(currentLabel)); | |
343 | currentLabel = (currentParticle) ? currentParticle->GetMother() : -1; | |
344 | } while (currentLabel >= 0); | |
345 | ||
346 | // Loop over prior clusters | |
347 | firstLabel = labels[iCluster1]; | |
348 | Int_t nCompatibleLabel = 1; | |
349 | Int_t currentParticle = 0; | |
350 | for (Int_t iCluster2 = iCluster1-1; iCluster2 >= 0; iCluster2--) { | |
351 | ||
352 | // if the number of clusters left is not enough the conditions cannot be fulfilled | |
353 | if (iCluster2 < halfCluster-nCompatibleLabel) break; | |
354 | ||
355 | if (labels[iCluster2] < 0) continue; | |
356 | ||
357 | // check if the cluster belong to the same particle or one of its ancestors | |
358 | Bool_t matchFound = kFALSE; | |
359 | for (Int_t iParticle = currentParticle; iParticle < nParticles; iParticle++) { | |
360 | if (labels[iCluster2] == chainLabels[iParticle]) { | |
361 | currentParticle = iParticle; | |
362 | matchFound = kTRUE; | |
363 | break; | |
364 | } | |
365 | } | |
366 | if (matchFound) nCompatibleLabel++; | |
367 | else continue; | |
368 | ||
369 | // add this station to the mask | |
370 | stationId = chId[iCluster2]/2; | |
371 | stationMask |= 1 << stationId; | |
372 | ||
373 | // count the number of fired chamber on stations 4 & 5 | |
374 | if (stationId > 2 && chId[iCluster2] < currentCh) { | |
375 | nChHitInSt45[stationId-3]++; | |
376 | currentCh = chId[iCluster2]; | |
377 | } | |
378 | ||
379 | // check if we matched enough clusters to tag the track as a decay | |
380 | if (nCompatibleLabel <= halfCluster || chId[iCluster2] > 3 || chainLabels[currentParticle] == firstLabel) continue; | |
381 | ||
382 | // check if this chain is better than already found | |
383 | if (chId[iCluster1] > lastCh) { | |
384 | decayLabel = firstLabel; | |
385 | lastCh = chId[iCluster1]; | |
386 | } | |
387 | ||
388 | // is there enough matched clusters on station 4 & 5 to make the track reconstructible? | |
389 | Bool_t isEnoughClOnSt45 = fRequest2ChInSameSt45 ? (nChHitInSt45[0] == 2 || nChHitInSt45[1] == 2) | |
390 | : (nChHitInSt45[0]+nChHitInSt45[1] >= 2); | |
391 | ||
392 | // is there any chance the current decay chain can still be tagged reconstructible? | |
393 | requestedStations = fRequestedStationMask >> stationId; | |
394 | isValid = (((stationMask >> stationId) & requestedStations) == requestedStations && | |
395 | (chId[iCluster2] > 5 || isEnoughClOnSt45)); | |
396 | ||
397 | // if not then we cannot do better with this trial | |
398 | if (!isValid) break; | |
399 | ||
400 | // take in priority the decay chain that can be tagged reconstructible | |
401 | if (((stationMask & fRequestedStationMask) == fRequestedStationMask) && isEnoughClOnSt45) { | |
402 | lastCh = chId[iCluster1]; | |
403 | isReconstructible = kTRUE; | |
404 | return firstLabel; | |
405 | } | |
406 | ||
407 | } | |
408 | ||
409 | } | |
410 | ||
411 | return decayLabel; | |
412 | } | |
413 | ||
414 | //---------------------------------------------------------------------- | |
415 | void AliAnalysisTaskESDMCLabelAddition::AddCompatibleClusters(const AliMUONTrack &track, const AliMUONTrack &trackRef, | |
416 | TArrayI *labels, Int_t *nLabels) const | |
417 | { | |
418 | /// Try to match clusters between track and trackRef and add the corresponding MC labels to the arrays | |
419 | ||
420 | Double_t chi2Max = 2. * fSigmaCut * fSigmaCut; // 2 because 2 quantities in chi2 | |
421 | ||
422 | // Loop over clusters of first track | |
423 | Int_t nCl1 = track.GetNClusters(); | |
424 | for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) { | |
425 | AliMUONVCluster *cluster1 = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCl1))->GetClusterPtr(); | |
426 | ||
427 | // Loop over clusters of second track | |
428 | Int_t nCl2 = trackRef.GetNClusters(); | |
429 | for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) { | |
430 | AliMUONVCluster *cluster2 = static_cast<AliMUONTrackParam*>(trackRef.GetTrackParamAtCluster()->UncheckedAt(iCl2))->GetClusterPtr(); | |
431 | ||
432 | // check DE Id | |
433 | if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue; | |
434 | ||
435 | // check local chi2 | |
436 | Double_t dX = cluster1->GetX() - cluster2->GetX(); | |
437 | Double_t dY = cluster1->GetY() - cluster2->GetY(); | |
438 | Double_t chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2()); | |
439 | if (chi2 > chi2Max) continue; | |
440 | ||
441 | // expand array if needed | |
442 | if (nLabels[iCl1] >= labels[iCl1].GetSize()) labels[iCl1].Set(2*labels[iCl1].GetSize()); | |
443 | ||
444 | // save label | |
445 | labels[iCl1][nLabels[iCl1]] = static_cast<Int_t>(trackRef.GetUniqueID()); | |
446 | nLabels[iCl1]++; | |
447 | break; | |
448 | ||
449 | } | |
450 | ||
451 | } | |
452 | ||
453 | } | |
454 | ||
455 | //---------------------------------------------------------------------- | |
456 | Int_t AliAnalysisTaskESDMCLabelAddition::IsDecayByLabel(const AliMUONTrack &track, Bool_t &isReconstructible, | |
457 | Int_t &lastCh) const | |
458 | { | |
459 | /// Check whether this track correspond to a decaying particle by using cluster MC labels. | |
460 | /// "lastCh" contains the chamber Id of the most downstream chamber hit by the decay chain | |
461 | ||
462 | Int_t nClusters = track.GetNClusters(); | |
04c32689 | 463 | if (nClusters <= 0) return -1; |
e5b66c63 | 464 | Int_t *chId = new Int_t[nClusters]; |
465 | Int_t *labels = new Int_t[nClusters]; | |
466 | ||
467 | // copy labels and chamber Ids | |
468 | for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) { | |
469 | AliMUONVCluster* cluster = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr(); | |
470 | chId[iCluster] = cluster->GetChamberId(); | |
471 | labels[iCluster] = cluster->GetMCLabel(); | |
472 | } | |
473 | ||
474 | // look for decay | |
475 | lastCh = 0; | |
476 | Int_t decayLabel = IsDecay(nClusters, chId, labels, isReconstructible, lastCh); | |
477 | ||
478 | delete[] chId; | |
479 | delete[] labels; | |
480 | ||
481 | return decayLabel; | |
482 | } | |
483 | ||
484 | //---------------------------------------------------------------------- | |
485 | Int_t AliAnalysisTaskESDMCLabelAddition::IsDecayByPosition(const AliMUONTrack &track, const AliMUONVTrackStore &trackRefStore, | |
486 | Bool_t &isReconstructible, Int_t &lastCh) const | |
487 | { | |
488 | /// Check whether this track correspond to a decaying particle by comparing clusters position | |
489 | /// All possible combinations of compatible clusters from every trackRefs are considered | |
490 | ||
491 | Int_t nClusters = track.GetNClusters(); | |
04c32689 | 492 | if (nClusters <= 0) return -1; |
e5b66c63 | 493 | Int_t *chId = new Int_t[nClusters]; |
494 | Int_t *nLabels = new Int_t[nClusters]; | |
495 | TArrayI *labels = new TArrayI[nClusters]; | |
496 | ||
497 | // copy chamber Ids | |
498 | for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) { | |
499 | AliMUONVCluster* cluster = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr(); | |
500 | chId[iCluster] = cluster->GetChamberId(); | |
501 | nLabels[iCluster] = 0; | |
502 | labels[iCluster].Set(100); | |
503 | } | |
504 | ||
505 | // loop over trackRef store and add label of compatible clusters | |
506 | TIter next1(trackRefStore.CreateIterator()); | |
507 | AliMUONTrack* trackRef; | |
508 | while ( ( trackRef = static_cast<AliMUONTrack*>(next1()) ) ) | |
509 | AddCompatibleClusters(track, *trackRef, labels, nLabels); | |
510 | ||
511 | // complete the arrays of labels with "-1" if no label was found for a given cluster | |
512 | for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) { | |
513 | if (nLabels[iCluster] == 0) { | |
514 | labels[iCluster][0] = -1; | |
515 | nLabels[iCluster]++; | |
516 | } | |
517 | } | |
518 | ||
519 | // loop over all possible combinations | |
520 | Int_t *iLabel = new Int_t[nClusters]; | |
521 | memset(iLabel,0,nClusters*sizeof(Int_t)); | |
522 | iLabel[nClusters-1] = -1; | |
523 | Int_t *currentLabels = new Int_t[nClusters]; | |
524 | Int_t decayLabel = -1; | |
525 | lastCh = 0; | |
526 | isReconstructible = kFALSE; | |
527 | while (kTRUE) { | |
528 | ||
529 | // go to the next combination | |
530 | Int_t iCl = nClusters-1; | |
531 | while (++iLabel[iCl] >= nLabels[iCl] && iCl > 0) iLabel[iCl--] = 0; | |
532 | if (iLabel[iCl] >= nLabels[iCl]) break; // no more combination | |
533 | ||
534 | // copy labels | |
535 | for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) | |
536 | currentLabels[iCluster] = labels[iCluster][iLabel[iCluster]]; | |
537 | ||
538 | // look for decay | |
539 | Int_t currentDecayLabel = IsDecay(nClusters, chId, currentLabels, isReconstructible, lastCh); | |
540 | if (currentDecayLabel >= 0) { | |
541 | decayLabel = currentDecayLabel; | |
542 | if (isReconstructible) break; | |
543 | } | |
544 | ||
545 | } | |
546 | ||
547 | delete[] chId; | |
548 | delete[] nLabels; | |
549 | delete[] labels; | |
550 | delete[] iLabel; | |
551 | delete[] currentLabels; | |
552 | ||
553 | return decayLabel; | |
554 | } | |
555 |