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 fAlignOCDBpath(""),
  57 fRecoParamOCDBpath(""),
  58 fRequestedStationMask(0),
  59 fRequest2ChInSameSt45(kFALSE),
  60 fExternalTrkSigmaCut(-1.),
  61 fSigmaCut(-1.),
  62 fExternalTrgSigmaCut(-1.),
  63 fSigmaCutTrig(-1.),
  64 fDecayAsFake(kFALSE)
  65 {
  66   /// Default constructor
  67 }
  68
  69
  70 //----------------------------------------------------------------------
  71 AliAnalysisTaskESDMCLabelAddition::AliAnalysisTaskESDMCLabelAddition(const char* name):
  72 AliAnalysisTaskSE(name),
  73 fDefaultStorage("raw://"),
  74 fAlignOCDBpath(""),
  75 fRecoParamOCDBpath(""),
  76 fRequestedStationMask(0),
  77 fRequest2ChInSameSt45(kFALSE),
  78 fExternalTrkSigmaCut(-1.),
  79 fSigmaCut(-1.),
  80 fExternalTrgSigmaCut(-1.),
  81 fSigmaCutTrig(-1.),
  82 fDecayAsFake(kFALSE)
  83 {
  84   /// Constructor
  85 }
  86
  87
  88 //----------------------------------------------------------------------
  89 void AliAnalysisTaskESDMCLabelAddition::UserCreateOutputObjects()
  90 {
  91   /// Create output objects
  92 }
  93
  94
  95 //----------------------------------------------------------------------
  96 void AliAnalysisTaskESDMCLabelAddition::NotifyRun()
  97 {
  98   /// Load OCDB inputs
  99
 100   // load OCDB objects only once
 101   if (fSigmaCut > 0) return;
 102
 103   // set OCDB location
 104   AliCDBManager* cdbm = AliCDBManager::Instance();
 105   if (cdbm->IsDefaultStorageSet()) printf("MCLabelAddition: CDB default storage already set!\n");
 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   }
 111   if (cdbm->GetRun() > -1) printf("MCLabelAddition: run number already set!\n");
 112   else cdbm->SetRun(fCurrentRunNumber);
 113
 114   // load recoParam
 115   const AliMUONRecoParam* recoParam = (AliMUONESDInterface::GetTracker())
 116   ? AliMUONESDInterface::GetTracker()->GetRecoParam()
 117   : AliMUONCDB::LoadRecoParam();
 118
 119   if (!recoParam) {
 120     fRequestedStationMask = 0;
 121     fRequest2ChInSameSt45 = kFALSE;
 122     fSigmaCut = -1.;
 123     fSigmaCutTrig = -1.;
 124     return;
 125   }
 126
 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();
 133
 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();
 142
 143 }
 144
 145
 146 //----------------------------------------------------------------------
 147 void AliAnalysisTaskESDMCLabelAddition::UserExec(Option_t */*option*/)
 148 {
 149   /// Execute analysis for current event
 150
 151   AliDebug(1, Form("MCLabel Addition: Analysing event # %5d\n",(Int_t) Entry()));
 152
 153   // make sure necessary information from OCDB have been loaded
 154   if (fSigmaCut < 0) return;
 155
 156   /// Load ESD event
 157   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 158   if (!esd) {
 159     AliError("Cannot get input event");
 160     return;
 161   }
 162
 163   // Load MC event
 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
 170
 171   // Get reference tracks
 172   AliMUONRecoCheck rc(esd,mcH);
 173   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
 174   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
 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
 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
 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);
 202
 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
 263       } else if (decayLabel >= 0 && !fDecayAsFake) {
 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       }
 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     }
 293
 294   }
 295
 296 }
 297
 298
 299 //----------------------------------------------------------------------
 300 void AliAnalysisTaskESDMCLabelAddition::Terminate(Option_t */*option*/)
 301 {
 302   /// Terminate analysis
 303   AliDebug(2, "Terminate()");
 304 }
 305
 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();
 473   if (nClusters <= 0) return -1;
 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();
 502   if (nClusters <= 0) return -1;
 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