PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // This analysis provides a new list of clusters to be used in other analysis
  19 //
  20 // Author: Gustavo Conesa Balbastre,
  21 //         Adapted from analysis class from Deepa Thomas
  22 //
  23 //
  24 //_________________________________________________________________________
  25
  26 // --- Root ---
  27 #include "TString.h"
  28 #include "TRefArray.h"
  29 #include "TClonesArray.h"
  30 #include "TTree.h"
  31 #include "TGeoManager.h"
  32 #include "TROOT.h"
  33 #include "TInterpreter.h"
  34 #include "TFile.h"
  35 //#include "string.h"
  36
  37 // --- AliRoot Analysis Steering
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDEvent.h"
  41 #include "AliGeomManager.h"
  42 #include "AliVCaloCells.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliCDBManager.h"
  45 #include "AliCDBStorage.h"
  46 #include "AliCDBEntry.h"
  47 #include "AliLog.h"
  48 #include "AliVEventHandler.h"
  49 #include "AliAODInputHandler.h"
  50
  51 // --- EMCAL
  52 #include "AliEMCALRecParam.h"
  53 #include "AliEMCALAfterBurnerUF.h"
  54 #include "AliEMCALGeometry.h"
  55 #include "AliEMCALClusterizerNxN.h"
  56 #include "AliEMCALClusterizerv1.h"
  57 #include "AliEMCALClusterizerv2.h"
  58 #include "AliEMCALRecPoint.h"
  59 #include "AliEMCALDigit.h"
  60 #include "AliCaloCalibPedestal.h"
  61 #include "AliEMCALCalibData.h"
  62 #include "AliEMCALRecoUtils.h"
  63
  64 #include "AliAnalysisTaskEMCALClusterize.h"
  65
  66 ClassImp(AliAnalysisTaskEMCALClusterize)
  67
  68 //________________________________________________________________________
  69 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  70   : AliAnalysisTaskSE(name)
  71   , fGeom(0), fGeomName("EMCAL_COMPLETEV1"),   fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  72   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
  73   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  74   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  75   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
  76   , fFillAODFile(kTRUE), fFillAODHeader(0),    fFillAODCaloCells(0)
  77   , fRun(-1),            fRecoUtils(0),        fConfigName("")
  78   , fCellLabels(),       fCellSecondLabels(),  fCellTime()
  79   , fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
  80   , fSelectCell(kFALSE), fSelectCellMinE(0.005), fSelectCellMinFrac(0.001)
  81 {
  82   //ctor
  83   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
  84   for(Int_t j = 0; j < 24*48*11; j++)  {
  85     fCellLabels[j]       = -1;
  86     fCellSecondLabels[j] = -1;
  87     fCellTime[j]         =  0.;
  88   }
  89
  90   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  91   fClusterArr      = new TObjArray(100);
  92   fCaloClusterArr  = new TObjArray(1000);
  93   fRecParam        = new AliEMCALRecParam;
  94   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  95   fRecoUtils       = new AliEMCALRecoUtils();
  96
  97 }
  98
  99 //________________________________________________________________________
 100 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
 101   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
 102   , fGeom(0), fGeomName("EMCAL_COMPLETEV1"),    fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
 103   , fCalibData(0),        fPedestalData(0),     fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
 104   , fDigitsArr(0),        fClusterArr(0),       fCaloClusterArr(0)
 105   , fRecParam(0),         fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
 106   , fOutputAODBranch(0),  fOutputAODBranchName("newEMCALClusters")
 107   , fFillAODFile(kFALSE), fFillAODHeader(0),    fFillAODCaloCells(0)
 108   , fRun(-1),             fRecoUtils(0),        fConfigName("")
 109   , fCellLabels(),        fCellSecondLabels(),  fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
 110   , fSelectCell(kFALSE),  fSelectCellMinE(0.005), fSelectCellMinFrac(0.001)
 111 {
 112   // Constructor
 113   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
 114   for(Int_t j = 0; j < 24*48*11; j++)  {
 115     fCellLabels[j]       = -1;
 116     fCellSecondLabels[j] = -1;
 117     fCellTime[j]         =  0.;
 118   }
 119   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
 120   fClusterArr      = new TObjArray(100);
 121   fCaloClusterArr  = new TObjArray(100);
 122   fRecParam        = new AliEMCALRecParam;
 123   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
 124   fRecoUtils       = new AliEMCALRecoUtils();
 125 }
 126
 127
 128 //________________________________________________________________________
 129 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 130 {
 131   //dtor
 132
 133   if (fDigitsArr){
 134     fDigitsArr->Clear("C");
 135     delete fDigitsArr;
 136   }
 137
 138   if (fClusterArr){
 139     fClusterArr->Delete();
 140     delete fClusterArr;
 141   }
 142
 143   if (fCaloClusterArr){
 144     fCaloClusterArr->Delete();
 145     delete fCaloClusterArr;
 146   }
 147
 148   if(fClusterizer) delete fClusterizer;
 149   if(fUnfolder)    delete fUnfolder;
 150   if(fRecoUtils)   delete fRecoUtils;
 151
 152 }
 153
 154 //-------------------------------------------------------------------
 155 void AliAnalysisTaskEMCALClusterize::Init()
 156 {
 157   //Init analysis with configuration macro if available
 158
 159   if(gROOT->LoadMacro(fConfigName) >=0){
 160     printf("Configure analysis with %s\n",fConfigName.Data());
 161     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 162     fGeomName         = clus->fGeomName;
 163     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 164     fOCDBpath         = clus->fOCDBpath;
 165     fAccessOCDB       = clus->fAccessOCDB;
 166     fRecParam         = clus->fRecParam;
 167     fJustUnfold       = clus->fJustUnfold;
 168     fFillAODFile      = clus->fFillAODFile;
 169     fRecoUtils        = clus->fRecoUtils;
 170     fConfigName       = clus->fConfigName;
 171     fMaxEvent         = clus->fMaxEvent;
 172     fDoTrackMatching  = clus->fDoTrackMatching;
 173     fOutputAODBranchName = clus->fOutputAODBranchName;
 174     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
 175
 176   }
 177
 178 }
 179
 180 //-------------------------------------------------------------------
 181 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
 182 {
 183   // Init geometry, create list of output clusters
 184
 185   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;
 186   if(fOutputAODBranchName.Length()!=0){
 187     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 188     fOutputAODBranch->SetName(fOutputAODBranchName);
 189     AddAODBranch("TClonesArray", &fOutputAODBranch);
 190   }
 191   else {
 192     AliFatal("fOutputAODBranchName not set\n");
 193   }
 194 }
 195
 196 //________________________________________________________________________
 197 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
 198   // Put calo cells in standard branch
 199   AliVEvent * event = InputEvent();
 200   AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
 201   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 202
 203   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 204   aodEMcells.CreateContainer(nEMcell);
 205   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 206   Double_t calibFactor = 1.;
 207   for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 208     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 209     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 210     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 211
 212     if(fRecoUtils->IsRecalibrationOn()){
 213       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 214     }
 215
 216     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 217       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 218     }
 219     else {
 220       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 221     }
 222   }
 223   aodEMcells.Sort();
 224
 225 }
 226
 227 //________________________________________________________________________
 228 void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
 229   //Put event header information in standard AOD branch
 230
 231   AliVEvent* event      = InputEvent();
 232   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
 233   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
 234
 235   Double_t pos[3]   ;
 236   Double_t covVtx[6];
 237   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 238
 239   AliAODHeader* header = AODEvent()->GetHeader();
 240   header->SetRunNumber(event->GetRunNumber());
 241
 242   if(esdevent){
 243     TTree* tree = fInputHandler->GetTree();
 244     if (tree) {
 245       TFile* file = tree->GetCurrentFile();
 246       if (file) header->SetESDFileName(file->GetName());
 247     }
 248   }
 249   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 250
 251   header->SetBunchCrossNumber(event->GetBunchCrossNumber());
 252   header->SetOrbitNumber(event->GetOrbitNumber());
 253   header->SetPeriodNumber(event->GetPeriodNumber());
 254   header->SetEventType(event->GetEventType());
 255
 256   //Centrality
 257   if(event->GetCentrality()){
 258     header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
 259   }
 260   else{
 261     header->SetCentrality(0);
 262   }
 263
 264   //Trigger
 265   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 266   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 267   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 268   header->SetTriggerMask(event->GetTriggerMask());
 269   header->SetTriggerCluster(event->GetTriggerCluster());
 270   if(esdevent){
 271     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 272     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 273     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 274   }
 275   else if (aodevent){
 276     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 277     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 278     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 279   }
 280
 281   header->SetMagneticField(event->GetMagneticField());
 282   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 283
 284   header->SetZDCN1Energy(event->GetZDCN1Energy());
 285   header->SetZDCP1Energy(event->GetZDCP1Energy());
 286   header->SetZDCN2Energy(event->GetZDCN2Energy());
 287   header->SetZDCP2Energy(event->GetZDCP2Energy());
 288   header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
 289
 290   Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
 291   Float_t diamcov[3];
 292   event->GetDiamondCovXY(diamcov);
 293   header->SetDiamond(diamxy,diamcov);
 294   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 295   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 296
 297   //
 298   //
 299   Int_t nVertices = 1 ;/* = prim. vtx*/;
 300   Int_t nCaloClus = event->GetNumberOfCaloClusters();
 301
 302   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 303
 304   // Access to the AOD container of vertices
 305   TClonesArray &vertices = *(AODEvent()->GetVertices());
 306   Int_t jVertices=0;
 307
 308   // Add primary vertex. The primary tracks will be defined
 309   // after the loops on the composite objects (V0, cascades, kinks)
 310   event->GetPrimaryVertex()->GetXYZ(pos);
 311   Float_t chi = 0;
 312   if      (esdevent){
 313     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 314     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 315   }
 316   else if (aodevent){
 317     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 318     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 319   }
 320
 321   AliAODVertex * primary = new(vertices[jVertices++])
 322   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 323   primary->SetName(event->GetPrimaryVertex()->GetName());
 324   primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
 325
 326 }
 327
 328 //________________________________________________________________________
 329 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
 330 {
 331   // Main loop
 332   // Called for each event
 333
 334   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 335   Int_t eventN = Entry();
 336   if(aodIH) eventN = aodIH->GetReadEntry();
 337
 338   if (eventN > fMaxEvent) return;
 339   //printf("Clusterizer --- Event %d-- \n",eventN);
 340
 341   //Remove the contents of output list set in the previous event
 342   fOutputAODBranch->Clear("C");
 343
 344   //Magic line to write events to AOD file
 345   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
 346   LoadBranches();
 347
 348   //Init pointers, clusterizer, ocdb
 349   if(fAccessOCDB) AccessOCDB();
 350   InitClusterization();
 351
 352   //Get the event
 353   AliVEvent   * event    = 0;
 354   AliESDEvent * esdevent = 0;
 355
 356   //Fill output event with header
 357
 358   //Check if input event are embedded events
 359   //If so, take output event
 360   if (aodIH && aodIH->GetMergeEvents()) {
 361     //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
 362     event  = AODEvent();
 363
 364     if(!aodIH->GetMergeEMCALCells())
 365       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
 366
 367     //    printf("InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
 368     //           InputEvent()->GetEMCALCells()->GetNumberOfCells());
 369     //    printf("MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
 370     //           aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
 371     //    for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
 372     //        AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
 373     //        if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
 374     //    }
 375     //    printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
 376     //           AODEvent()->GetEMCALCells()->GetNumberOfCells());
 377
 378   }
 379   else {
 380     event =  InputEvent();
 381     esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 382     if(fFillAODCaloCells) FillAODCaloCells();
 383     if(fFillAODHeader)    FillAODHeader();
 384   }
 385
 386   if (!event) {
 387     Error("UserExec","Event not available");
 388     return;
 389   }
 390
 391   //-------------------------------------------------------------------------------------
 392   //Set the geometry matrix, for the first event, skip the rest
 393   //-------------------------------------------------------------------------------------
 394   if(!fGeomMatrixSet){
 395     if(fLoadGeomMatrices){
 396       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 397         if(fGeomMatrix[mod]){
 398           if(DebugLevel() > 1)
 399             fGeomMatrix[mod]->Print();
 400           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
 401         }
 402         fGeomMatrixSet=kTRUE;
 403       }//SM loop
 404     }//Load matrices
 405     else if(!gGeoManager){
 406       Info("UserExec","Get geo matrices from data");
 407       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 408       if(!strcmp(event->GetName(),"AliAODEvent")) {
 409         if(DebugLevel() > 1)
 410           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 411       }//AOD
 412       else {
 413         if(DebugLevel() > 1)
 414           Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 415         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 416         if(!esd) {
 417           Error("UserExec","This event does not contain ESDs?");
 418           return;
 419         }
 420         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 421           if(DebugLevel() > 1)
 422             esd->GetEMCALMatrix(mod)->Print();
 423           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 424         }
 425         fGeomMatrixSet=kTRUE;
 426       }//ESD
 427     }//Load matrices from Data
 428
 429     //Recover time dependent corrections, put then in recalibration histograms. Do it once
 430     fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
 431
 432   }//first event
 433
 434   //Get the list of cells needed for unfolding and reclustering
 435   AliVCaloCells *cells   = event->GetEMCALCells();
 436   Int_t    nClustersOrg  = 0;
 437   Double_t cellAmplitude = 0;
 438   Double_t cellTime      = 0;
 439   Short_t  cellNumber    = 0;
 440
 441   //-------------------------------------------
 442   //---------Unfolding clusters----------------
 443   //-------------------------------------------
 444   if (fJustUnfold) {
 445
 446     //Fill the array with the EMCAL clusters, copy them
 447     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
 448     {
 449       AliVCluster *clus = event->GetCaloCluster(i);
 450       if(clus->IsEMCAL()){
 451
 452         //recalibrate/remove bad channels/etc if requested
 453         if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
 454           continue;
 455         }
 456
 457         if(fRecoUtils->IsRecalibrationOn()){
 458
 459           //Calibrate cluster
 460           //printf("Energy before %f ",clus->E());
 461           fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 462           //printf("; after %f\n",clus->E());
 463
 464           //CalibrateCells
 465           for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 466           {
 467             if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 468               break;
 469
 470
 471             Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 472             fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 473             fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 474
 475             //Do not include bad channels found in analysis?
 476             if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 477                 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 478               fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
 479               fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
 480               fCellTime[cellNumber]        = 0.;
 481               continue;
 482             }
 483
 484             cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 485
 486           }// cells loop
 487         }// recalibrate
 488
 489         //Cast to ESD or AOD, needed to create the cluster array
 490         AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 491         AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 492         if     (esdCluster){
 493           fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 494         }//ESD
 495         else if(aodCluster){
 496           fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 497         }//AOD
 498         else
 499           Warning("UserExec()"," - Wrong CaloCluster type?");
 500         nClustersOrg++;
 501       }
 502     }
 503
 504     //Do the unfolding
 505     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 506
 507     //CLEAN-UP
 508     fUnfolder->Clear();
 509
 510   } // just unfold ESD/AOD cluster
 511
 512   //-------------------------------------------
 513   //---------- Recluster cells ----------------
 514   //-------------------------------------------
 515
 516   else{
 517     //-------------------------------------------------------------------------------------
 518     //Transform CaloCells into Digits
 519     //-------------------------------------------------------------------------------------
 520
 521     //In case of MC, first loop on the clusters and fill MC label to array
 522     //.....................................................................
 523
 524     //    for(Int_t j = 0; j < 24*48*11; j++)  {
 525     //      if(fCellLabels[j]      !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j]      ) ;
 526     //      if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
 527     //    }
 528
 529     Int_t nClusters = event->GetNumberOfCaloClusters();
 530     if(aodIH && aodIH->GetEventToMerge())  //Embedding
 531       nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
 532     for (Int_t i = 0; i < nClusters; i++)
 533     {
 534       AliVCluster *clus = 0;
 535       if(aodIH && aodIH->GetEventToMerge()) //Embedding
 536         clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
 537       else
 538         clus = event->GetCaloCluster(i);
 539
 540       if(!clus) {
 541         printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
 542         return;
 543       }
 544
 545       if(clus->IsEMCAL()){
 546         //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
 547         Int_t label = clus->GetLabel();
 548         Int_t label2 = -1 ;
 549         if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 550         UShort_t * index    = clus->GetCellsAbsId() ;
 551         for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
 552           fCellLabels[index[icell]]       = label;
 553           fCellSecondLabels[index[icell]] = label2;
 554           fCellTime[icell]                = clus->GetTOF();
 555
 556           //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
 557         }
 558         nClustersOrg++;
 559       }
 560     }
 561
 562     // Create digits
 563     //.................
 564     Int_t idigit =  0;
 565     Int_t id     = -1;
 566     Float_t amp  = -1;
 567     Float_t time = -1;
 568
 569     TTree *digitsTree = new TTree("digitstree","digitstree");
 570     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 571
 572     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 573     {
 574       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 575         break;
 576
 577       time = cellTime;
 578       amp  = cellAmplitude;
 579       id   = cellNumber;
 580
 581       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 582       fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
 583       fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 584
 585       //Do not include bad channels found in analysis?
 586       if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 587          fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 588         fCellLabels[id]      =-1; //reset the entry in the array for next event
 589         fCellSecondLabels[id]=-1; //reset the entry in the array for next event
 590         fCellTime[id]        = 0.;
 591         //printf("Remove channel %d\n",id);
 592         continue;
 593       }
 594
 595       //Recalibrate?
 596       if(fRecoUtils->IsRecalibrationOn()){
 597         //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
 598         amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 599       }
 600
 601       // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
 602       if (time*1e9 < 1.) time = fCellTime[id];
 603
 604       //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 605       new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 606       //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
 607       //else                  printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
 608       fCellLabels[id]      =-1; //reset the entry in the array for next event
 609
 610       //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
 611       idigit++;
 612     }
 613
 614     //Fill the tree with digits
 615     digitsTree->Fill();
 616
 617     //-------------------------------------------------------------------------------------
 618     //Do the clusterization
 619     //-------------------------------------------------------------------------------------
 620     TTree *clustersTree = new TTree("clustertree","clustertree");
 621
 622     fClusterizer->SetInput(digitsTree);
 623     fClusterizer->SetOutput(clustersTree);
 624     fClusterizer->Digits2Clusters("");
 625
 626     //-------------------------------------------------------------------------------------
 627     //Transform the recpoints into AliVClusters
 628     //-------------------------------------------------------------------------------------
 629
 630     clustersTree->SetBranchStatus("*",0); //disable all branches
 631     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 632
 633     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 634     branch->SetAddress(&fClusterArr);
 635     branch->GetEntry(0);
 636
 637     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 638
 639     if(!fCaloClusterArr){
 640       printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 641       return;
 642     }
 643
 644     if( DebugLevel() > 0 && fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
 645       printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
 646              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
 647              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 648     }
 649
 650     //Reset the array with second labels for this event
 651     memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
 652
 653     //---CLEAN UP-----
 654     fClusterizer->Clear();
 655     fDigitsArr  ->Clear("C");
 656     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 657
 658     clustersTree->Delete("all");
 659     digitsTree  ->Delete("all");
 660   }
 661
 662   //Recalculate track-matching for the new clusters, only with ESDs
 663   if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
 664
 665
 666   //-------------------------------------------------------------------------------------
 667   //Put the new clusters in the AOD list
 668   //-------------------------------------------------------------------------------------
 669
 670   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
 671   //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
 672   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
 673     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 674     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
 675
 676     //Add matched track, if any, only with ESDs
 677     if(esdevent && fDoTrackMatching){
 678       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 679       if(trackIndex >= 0){
 680         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
 681         if(DebugLevel() > 1)
 682           Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
 683       }
 684     }
 685
 686     //In case of new bad channels, recalculate distance to bad channels
 687     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
 688       //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
 689       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
 690       //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
 691     }
 692
 693     //    if(newCluster->GetNLabels()>0){
 694     //      printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
 695     //      UShort_t * newindex    = newCluster->GetCellsAbsId() ;
 696     //      for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
 697     //       printf("\t absID %d\n",newindex[icell]);
 698     //     }
 699     //    }
 700
 701     //    printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
 702     //    printf("labels: ");
 703     //    for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
 704     //      printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
 705     //     printf("\n ");
 706
 707     newCluster->SetID(i);
 708     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 709   }
 710
 711   //if(fOutputAODBranchName.Length()!=0)
 712   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
 713
 714   //---CLEAN UP-----
 715   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
 716 }
 717
 718 //_____________________________________________________________________
 719 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 720 {
 721   //Access to OCDB stuff
 722
 723   AliVEvent * event = InputEvent();
 724   if (!event)
 725   {
 726     Warning("AccessOCDB","Event not available!!!");
 727     return kFALSE;
 728   }
 729
 730   if (event->GetRunNumber()==fRun)
 731     return kTRUE;
 732   fRun = event->GetRunNumber();
 733
 734   if(DebugLevel() > 1 )
 735     Info("AccessODCD()"," Begin");
 736
 737   //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 738
 739   AliCDBManager *cdb = AliCDBManager::Instance();
 740
 741
 742   if (fOCDBpath.Length()){
 743     cdb->SetDefaultStorage(fOCDBpath.Data());
 744     Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
 745   }
 746
 747   cdb->SetRun(event->GetRunNumber());
 748
 749   //
 750   // EMCAL from RAW OCDB
 751   if (fOCDBpath.Contains("alien:"))
 752   {
 753     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 754     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 755   }
 756
 757   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 758
 759   // init parameters:
 760
 761   //Get calibration parameters
 762   if(!fCalibData)
 763   {
 764     AliCDBEntry *entry = (AliCDBEntry*)
 765     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 766
 767     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 768   }
 769
 770   if(!fCalibData)
 771     AliFatal("Calibration parameters not found in CDB!");
 772
 773   //Get calibration parameters
 774   if(!fPedestalData)
 775   {
 776     AliCDBEntry *entry = (AliCDBEntry*)
 777     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 778
 779     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 780   }
 781
 782   if(!fPedestalData)
 783     AliFatal("Dead map not found in CDB!");
 784
 785   return kTRUE;
 786 }
 787
 788 //________________________________________________________________________________________
 789 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 790 {
 791   //Select clusterization/unfolding algorithm and set all the needed parameters
 792
 793   if (fJustUnfold){
 794     // init the unfolding afterburner
 795     delete fUnfolder;
 796     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 797     return;
 798   }
 799
 800   //First init the clusterizer
 801   delete fClusterizer;
 802   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 803     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 804   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 805     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
 806   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){
 807     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 808     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
 809     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
 810   } else {
 811     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 812   }
 813
 814   //Now set the parameters
 815   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 816   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 817   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 818   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 819   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 820   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 821   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 822   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 823   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 824   fClusterizer->SetJustClusters          ( kTRUE                               );
 825   //In case of unfolding after clusterization is requested, set the corresponding parameters
 826   if(fRecParam->GetUnfold()){
 827     Int_t i=0;
 828     for (i = 0; i < 8; i++) {
 829       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 830     }//end of loop over parameters
 831     for (i = 0; i < 3; i++) {
 832       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 833       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 834     }//end of loop over parameters
 835
 836     fClusterizer->InitClusterUnfolding();
 837
 838   }// to unfold
 839 }
 840
 841 //________________________________________________________________________________________
 842 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
 843 {
 844   // Restore clusters from recPoints
 845   // Cluster energy, global position, cells and their amplitude fractions are restored
 846   Int_t j = 0;
 847   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
 848   {
 849     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
 850
 851     const Int_t ncells = recPoint->GetMultiplicity();
 852     Int_t ncellsTrue = 0;
 853
 854     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
 855
 856     // cells and their amplitude fractions
 857     UShort_t   absIds[ncells];
 858     Double32_t ratios[ncells];
 859
 860     //For later check embedding
 861     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 862
 863     Float_t clusterE = 0;
 864     for (Int_t c = 0; c < ncells; c++) {
 865       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
 866
 867       absIds[ncellsTrue] = digit->GetId();
 868       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 869
 870       // In case of unfolding, remove digits with unfolded energy too low
 871       if(fSelectCell){
 872         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  {
 873
 874           if(DebugLevel() > 1)  {
 875             printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
 876                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
 877           }
 878
 879           continue;
 880
 881         } // if cuts
 882       }// Select cells
 883
 884       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
 885       ncellsTrue++;
 886       clusterE  +=recPoint->GetEnergiesList()[c];
 887
 888       // In case of embedding, fill ratio with amount of signal,
 889       if (aodIH && aodIH->GetMergeEvents()) {
 890
 891         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
 892         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
 893         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
 894
 895         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 896         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 897         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 898         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
 899
 900         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
 901         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
 902
 903       }//Embedding
 904
 905     }// cluster cell loop
 906
 907     if (ncellsTrue < 1) {
 908       if (DebugLevel() > 1)
 909         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
 910                recPoint->GetEnergy(), ncells);
 911       continue;
 912     }
 913
 914     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
 915
 916     if(clusterE <  fRecParam->GetClusteringThreshold()) {
 917       if (DebugLevel()>1)
 918         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
 919       continue;
 920     }
 921
 922     TVector3 gpos;
 923     Float_t g[3];
 924
 925     // calculate new cluster position
 926     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
 927     recPoint->GetGlobalPosition(gpos);
 928     gpos.GetXYZ(g);
 929
 930     // create a new cluster
 931     //AliAODCaloCluster *clus = new AliAODCaloCluster();
 932     (*clusArray)[j] = new AliAODCaloCluster() ;
 933     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
 934     j++;
 935     clus->SetType(AliVCluster::kEMCALClusterv1);
 936     clus->SetE(clusterE);//recPoint->GetEnergy());
 937     clus->SetPosition(g);
 938     clus->SetNCells(ncellsTrue);
 939     clus->SetCellsAbsId(absIds);
 940     clus->SetCellsAmplitudeFraction(ratios);
 941     clus->SetChi2(-1); //not yet implemented
 942     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
 943     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
 944     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
 945
 946     if(ncells == ncellsTrue){
 947       Float_t elipAxis[2];
 948       recPoint->GetElipsAxis(elipAxis);
 949       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 950       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 951       clus->SetDispersion(recPoint->GetDispersion());
 952     }
 953     else if(fSelectCell){
 954       // In case some cells rejected, in unfolding case, recalculate
 955       // shower shape parameters and position
 956       AliVCaloCells* cells = 0x0;
 957       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
 958       else                                  cells = InputEvent()->GetEMCALCells();
 959       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
 960       fRecoUtils->RecalculateClusterPID(clus);
 961       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
 962
 963 //      Float_t elipAxis[2];
 964 //      recPoint->GetElipsAxis(elipAxis);
 965 //      if(ncellsTrue > 2)
 966 //      printf("SS, old: l0 %f, l1 %f, D %f; New l0 %f, l1 %f, D %f\n",
 967 //             elipAxis[0]*elipAxis[0],elipAxis[1]*elipAxis[1], recPoint->GetDispersion(),
 968 //             clus->GetM02(),clus->GetM20(),clus->GetDispersion());
 969     }
 970
 971     //MC
 972     Int_t  parentMult  = 0;
 973     Int_t *parentList = recPoint->GetParents(parentMult);
 974     clus->SetLabel(parentList, parentMult);
 975
 976     //    if(parentMult)printf("RecToESD: Labels %d",parentMult);
 977     //    for (Int_t iii = 0; iii < parentMult; iii++) {
 978     //      printf("\t label %d\n",parentList[iii]);
 979     //    }
 980
 981     //Write the second major contributor to each MC cluster.
 982     Int_t iNewLabel ;
 983     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
 984
 985       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
 986       if(idCell>=0){
 987         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
 988         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
 989         {
 990           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
 991           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
 992         }
 993         if (iNewLabel == 1)
 994         {
 995           //Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
 996           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
 997           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
 998             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
 999           }
1000
1001           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1002           //if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
1003           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1004           delete [] newLabelArray;
1005         }
1006       }//positive cell number
1007     }
1008
1009   } // recPoints loop
1010
1011 }