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