PWG/muondep/AliAnalysisTaskESDMCLabelAddition.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // ROOT includes
  19 #include <TArrayI.h>
  20 #include <TClonesArray.h>
  21
  22 // STEER includes
  23 #include "AliESDEvent.h"
  24 #include "AliESDMuonTrack.h"
  25 #include "AliLog.h"
  26 #include "AliMCEventHandler.h"
  27 #include "AliCDBManager.h"
  28 #include "AliMCParticle.h"
  29 #include "AliMCEvent.h"
  30
  31 // ANALYSIS includes
  32 #include "AliAnalysisManager.h"
  33
  34 // MUON includes
  35 #include "AliMUONCDB.h"
  36 #include "AliMUONRecoParam.h"
  37 #include "AliMUONRecoCheck.h"
  38 #include "AliMUONESDInterface.h"
  39 #include "AliMUONVTrackReconstructor.h"
  40 #include "AliMUONTrack.h"
  41 #include "AliMUONTrackParam.h"
  42 #include "AliMUONVTrackStore.h"
  43 #include "AliMUONTriggerTrack.h"
  44 #include "AliMUONVTriggerTrackStore.h"
  45 #include "AliMUONLocalTrigger.h"
  46 #include "AliMUONVCluster.h"
  47
  48 #include "AliAnalysisTaskESDMCLabelAddition.h"
  49
  50 ClassImp(AliAnalysisTaskESDMCLabelAddition)
  51
  52 //----------------------------------------------------------------------
  53 AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition():
  54 AliAnalysisTaskSE(),
  55 fDefaultStorage(""),
  56 fRequestedStationMask(0),
  57 fRequest2ChInSameSt45(kFALSE),
  58 fExternalTrkSigmaCut(-1.),
  59 fSigmaCut(-1.),
  60 fExternalTrgSigmaCut(-1.),
  61 fSigmaCutTrig(-1.)
  62 {
  63   /// Default constructor
  64 }
  65
  66
  67 //----------------------------------------------------------------------
  68 AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition(const char* name):
  69 AliAnalysisTaskSE(name),
  70 fDefaultStorage("raw://"),
  71 fRequestedStationMask(0),
  72 fRequest2ChInSameSt45(kFALSE),
  73 fExternalTrkSigmaCut(-1.),
  74 fSigmaCut(-1.),
  75 fExternalTrgSigmaCut(-1.),
  76 fSigmaCutTrig(-1.)
  77 {
  78   /// Constructor
  79 }
  80
  81
  82 //----------------------------------------------------------------------
  83 void AliAnalysisTaskESDMCLabelAddition::UserCreateOutputObjects()
  84 {
  85   /// Create output objects
  86 }
  87
  88
  89 //----------------------------------------------------------------------
  90 void AliAnalysisTaskESDMCLabelAddition::NotifyRun()
  91 {
  92   /// Load OCDB inputs
  93
  94   // load OCDB objects only once
  95   if (fSigmaCut > 0) return;
  96
  97   // set OCDB location
  98   AliCDBManager* cdbm = AliCDBManager::Instance();
  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);
 103
 104   // load recoParam
 105   const AliMUONRecoParam* recoParam = (AliMUONESDInterface::GetTracker())
 106   ? AliMUONESDInterface::GetTracker()->GetRecoParam()
 107   : AliMUONCDB::LoadRecoParam();
 108
 109   if (!recoParam) {
 110     fRequestedStationMask = 0;
 111     fRequest2ChInSameSt45 = kFALSE;
 112     fSigmaCut = -1.;
 113     fSigmaCutTrig = -1.;
 114     return;
 115   }
 116
 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();
 123
 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();
 132
 133 }
 134
 135
 136 //----------------------------------------------------------------------
 137 void AliAnalysisTaskESDMCLabelAddition::UserExec(Option_t */*option*/)
 138 {
 139   /// Execute analysis for current event
 140
 141   AliDebug(1, Form("MCLabel Addition: Analysing event # %5d\n",(Int_t) Entry()));
 142
 143   // make sure necessary information from OCDB have been loaded
 144   if (fSigmaCut < 0) return;
 145
 146   /// Load ESD event
 147   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 148   if (!esd) {
 149     AliError("Cannot get input event");
 150     return;
 151   }
 152
 153   // Load MC event
 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
 160
 161   // Get reference tracks
 162   AliMUONRecoCheck rc(esd,mcH);
 163   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
 164   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
 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
 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
 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);
 192
 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       }
 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     }
 283
 284   }
 285
 286 }
 287
 288
 289 //----------------------------------------------------------------------
 290 void AliAnalysisTaskESDMCLabelAddition::Terminate(Option_t */*option*/)
 291 {
 292   /// Terminate analysis
 293   AliDebug(2, "Terminate()");
 294 }
 295
 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();
 463   if (nClusters <= 0) return -1;
 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();
 492   if (nClusters <= 0) return -1;
 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