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