PWGGA/EMCALTasks/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
  16 //_________________________________________________________________________
  17 // This analysis provides a new list of clusters to be used in other analysis
  18 //
  19 // Author: Gustavo Conesa Balbastre,
  20 //         Adapted from analysis class from Deepa Thomas
  21 //
  22 // $Id$
  23 //_________________________________________________________________________
  24
  25 // --- Root ---
  26 #include <TString.h>
  27 #include <TRefArray.h>
  28 #include <TClonesArray.h>
  29 #include <TTree.h>
  30 #include <TGeoManager.h>
  31 #include <TROOT.h>
  32 #include <TInterpreter.h>
  33 #include <TFile.h>
  34
  35 // --- AliRoot Analysis Steering
  36 #include "AliAnalysisTask.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliESDEvent.h"
  39 #include "AliGeomManager.h"
  40 #include "AliVCaloCells.h"
  41 #include "AliAODCaloCluster.h"
  42 #include "AliCDBManager.h"
  43 #include "AliCDBStorage.h"
  44 #include "AliCDBEntry.h"
  45 #include "AliLog.h"
  46 #include "AliVEventHandler.h"
  47 #include "AliAODInputHandler.h"
  48 #include "AliOADBContainer.h"
  49 #include "AliAODMCParticle.h"
  50
  51 // --- EMCAL
  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
  62 #include "AliAnalysisTaskEMCALClusterize.h"
  63
  64 ClassImp(AliAnalysisTaskEMCALClusterize)
  65
  66 //______________________________________________________________________________
  67 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  68 : AliAnalysisTaskSE(name)
  69 , fEvent(0)
  70 , fGeom(0),               fGeomName("")
  71 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  72 , fCalibData(0),          fPedestalData(0)
  73 , fOCDBpath(""),          fAccessOCDB(kFALSE)
  74 , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
  75 , fRecParam(0),           fClusterizer(0)
  76 , fUnfolder(0),           fJustUnfold(kFALSE)
  77 , fOutputAODBranch(0),    fOutputAODBranchName(""),   fOutputAODBranchSet(0)
  78 , fFillAODFile(kFALSE),   fFillAODHeader(0)
  79 , fFillAODCaloCells(0),   fRun(-1)
  80 , fRecoUtils(0),          fConfigName("")
  81 , fOrgClusterCellId()
  82 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
  83 , fCellMatchdEta(),       fCellMatchdPhi()
  84 , fRecalibrateWithClusterTime(0)
  85 , fMaxEvent(0),           fDoTrackMatching(kFALSE)
  86 , fSelectCell(kFALSE),    fSelectCellMinE(0),         fSelectCellMinFrac(0)
  87 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
  88 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
  89 , fOADBSet(kFALSE),       fAccessOADB(kTRUE),         fOADBFilePath("")
  90 , fCentralityClass(""),   fSelectEMCALEvent(0)
  91 , fEMCALEnergyCut(0.),    fEMCALNcellsCut (0)
  92 , fSetCellMCLabelFromCluster(0)
  93 , fRemapMCLabelForAODs(0)
  94 {
  95   // Constructor
  96
  97   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
  98
  99   ResetArrays();
 100
 101   fCentralityBin[0] = fCentralityBin[1]=-1;
 102
 103 }
 104
 105 //______________________________________________________________
 106 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
 107 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
 108 , fEvent(0)
 109 , fGeom(0),                 fGeomName("")
 110 , fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
 111 , fCalibData(0),            fPedestalData(0)
 112 , fOCDBpath(""),            fAccessOCDB(kFALSE)
 113 , fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
 114 , fRecParam(0),             fClusterizer(0)
 115 , fUnfolder(0),             fJustUnfold(kFALSE)
 116 , fOutputAODBranch(0),      fOutputAODBranchName(""),   fOutputAODBranchSet(0)
 117 , fFillAODFile(kFALSE),     fFillAODHeader(0)
 118 , fFillAODCaloCells(0),     fRun(-1)
 119 , fRecoUtils(0),            fConfigName("")
 120 , fOrgClusterCellId()
 121 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
 122 , fCellMatchdEta(),         fCellMatchdPhi()
 123 , fRecalibrateWithClusterTime(0)
 124 , fMaxEvent(0),             fDoTrackMatching(kFALSE)
 125 , fSelectCell(kFALSE),      fSelectCellMinE(0),         fSelectCellMinFrac(0)
 126 , fRemoveLEDEvents(kTRUE),  fRemoveExoticEvents(kFALSE)
 127 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
 128 , fOADBSet(kFALSE),         fAccessOADB(kTRUE),        fOADBFilePath("")
 129 , fCentralityClass(""),     fSelectEMCALEvent(0)
 130 , fEMCALEnergyCut(0.),      fEMCALNcellsCut (0)
 131 , fSetCellMCLabelFromCluster(0)
 132 , fRemapMCLabelForAODs(0)
 133 {
 134   // Constructor
 135
 136   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 137
 138   ResetArrays();
 139
 140   fCentralityBin[0] = fCentralityBin[1]=-1;
 141
 142 }
 143
 144
 145 //_______________________________________________________________
 146 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 147 {
 148   //dtor
 149
 150   if (fDigitsArr)
 151   {
 152     fDigitsArr->Clear("C");
 153     delete fDigitsArr;
 154   }
 155
 156   if (fClusterArr)
 157   {
 158     fClusterArr->Delete();
 159     delete fClusterArr;
 160   }
 161
 162   if (fCaloClusterArr)
 163   {
 164     fCaloClusterArr->Delete();
 165     delete fCaloClusterArr;
 166   }
 167
 168   if(fClusterizer) delete fClusterizer;
 169   if(fUnfolder)    delete fUnfolder;
 170   if(fRecoUtils)   delete fRecoUtils;
 171
 172 }
 173
 174 //_______________________________________________________
 175 Bool_t AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL()
 176 {
 177   // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
 178
 179   if(!fSelectEMCALEvent)   return kTRUE; // accept
 180
 181   if(fEMCALEnergyCut <= 0) return kTRUE; // accept
 182
 183   Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
 184   AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
 185   Int_t           bc       = InputEvent() -> GetBunchCrossNumber();
 186
 187   for(Int_t icalo = 0; icalo < nCluster; icalo++)
 188   {
 189     AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
 190
 191     if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
 192        fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
 193     {
 194
 195       if (fDebug > 0)
 196         printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n",
 197                              clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
 198
 199       return kTRUE;
 200     }
 201
 202   }// loop
 203
 204   if (fDebug > 0)
 205     printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Reject \n");
 206
 207   return kFALSE;
 208
 209 }
 210
 211 //_______________________________________________
 212 void AliAnalysisTaskEMCALClusterize::AccessOADB()
 213 {
 214   // Set the AODB calibration, bad channels etc. parameters at least once
 215   // alignment matrices from OADB done in SetGeometryMatrices
 216
 217   //Set it only once
 218   if(fOADBSet) return ;
 219
 220   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 221   TString pass      = GetPass();
 222
 223   printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
 224
 225   Int_t nSM = fGeom->GetNumberOfSuperModules();
 226
 227   // Bad map
 228   if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
 229   {
 230     AliOADBContainer *contBC=new AliOADBContainer("");
 231     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
 232
 233     TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
 234
 235     if(arrayBC)
 236     {
 237         fRecoUtils->SwitchOnDistToBadChannelRecalculation();
 238         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
 239
 240         for (Int_t i=0; i<nSM; ++i)
 241         {
 242           TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
 243
 244           if (hbm)
 245             delete hbm;
 246
 247           hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
 248
 249           if (!hbm)
 250           {
 251             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
 252             continue;
 253           }
 254
 255           hbm->SetDirectory(0);
 256           fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
 257
 258         } // loop
 259     } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
 260   }  // Remove bad
 261
 262   // Energy Recalibration
 263   if(fRecoUtils->IsRecalibrationOn())
 264   {
 265     AliOADBContainer *contRF=new AliOADBContainer("");
 266
 267     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
 268
 269     TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
 270
 271     if(recal)
 272     {
 273       TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
 274
 275       if(recalpass)
 276       {
 277         TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
 278
 279         if(recalib)
 280         {
 281           printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
 282           for (Int_t i=0; i<nSM; ++i)
 283           {
 284             TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
 285
 286             if (h)
 287               delete h;
 288
 289             h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
 290
 291             if (!h)
 292             {
 293               AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
 294               continue;
 295             }
 296
 297             h->SetDirectory(0);
 298
 299             fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
 300           } // SM loop
 301         }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok
 302       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
 303     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
 304
 305   } // Recalibration on
 306
 307   // Energy Recalibration, apply on top of previous calibration factors
 308   if(fRecoUtils->IsRunDepRecalibrationOn())
 309   {
 310     AliOADBContainer *contRFTD=new AliOADBContainer("");
 311
 312     contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
 313
 314     TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
 315
 316     //If it did not exist for this run, get closes one
 317     if (!htd)
 318     {
 319       AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
 320       // let's get the closest runnumber instead then..
 321       Int_t lower = 0;
 322       Int_t ic = 0;
 323       Int_t maxEntry = contRFTD->GetNumberOfEntries();
 324
 325       while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
 326         lower = ic;
 327         ic++;
 328       }
 329
 330       Int_t closest = lower;
 331       if ( (ic<maxEntry) &&
 332           (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
 333         closest = ic;
 334       }
 335
 336       AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
 337       htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
 338     }
 339
 340     if(htd)
 341     {
 342       printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
 343
 344       for (Int_t ism=0; ism<nSM; ++ism)
 345       {
 346         for (Int_t icol=0; icol<48; ++icol)
 347         {
 348           for (Int_t irow=0; irow<24; ++irow)
 349           {
 350             Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
 351
 352             Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
 353             factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
 354             //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
 355             //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
 356             fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
 357           } // columns
 358         } // rows
 359       } // SM loop
 360     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
 361   } // Run by Run T calibration
 362
 363   // Time Recalibration
 364   if(fRecoUtils->IsTimeRecalibrationOn())
 365   {
 366     AliOADBContainer *contTRF=new AliOADBContainer("");
 367
 368     contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
 369
 370     TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
 371
 372     if(trecal)
 373     {
 374       TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
 375
 376       if(trecalpass)
 377       {
 378         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n");
 379         for (Int_t ibc = 0; ibc < 4; ++ibc)
 380         {
 381           TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
 382
 383           if (h)
 384             delete h;
 385
 386           h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
 387
 388           if (!h)
 389           {
 390             AliError(Form("Could not load hAllTimeAvBC%d",ibc));
 391             continue;
 392           }
 393
 394           h->SetDirectory(0);
 395
 396           fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
 397         } // bunch crossing loop
 398       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
 399     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
 400
 401   } // Time recalibration on
 402
 403   // Parameters already set once, so do not it again
 404   fOADBSet = kTRUE;
 405
 406 }
 407
 408 //_________________________________________________
 409 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 410 {
 411   //Access to OCDB stuff
 412
 413   fEvent = InputEvent();
 414   if (!fEvent)
 415   {
 416     Warning("AccessOCDB","Event not available!!!");
 417     return kFALSE;
 418   }
 419
 420   if (fEvent->GetRunNumber()==fRun)
 421     return kTRUE;
 422   fRun = fEvent->GetRunNumber();
 423
 424   if(DebugLevel() > 1 )
 425     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Begin");
 426
 427   AliCDBManager *cdb = AliCDBManager::Instance();
 428
 429
 430   if (fOCDBpath.Length())
 431   {
 432     cdb->SetDefaultStorage(fOCDBpath.Data());
 433     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
 434   }
 435
 436   cdb->SetRun(fEvent->GetRunNumber());
 437
 438   //
 439   // EMCAL from RAW OCDB
 440   if (fOCDBpath.Contains("alien:"))
 441   {
 442     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 443     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 444   }
 445
 446   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 447
 448   // init parameters:
 449
 450   //Get calibration parameters
 451   if(!fCalibData)
 452   {
 453     AliCDBEntry *entry = (AliCDBEntry*)
 454     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 455
 456     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 457   }
 458
 459   if(!fCalibData)
 460     AliFatal("Calibration parameters not found in CDB!");
 461
 462   //Get calibration parameters
 463   if(!fPedestalData)
 464   {
 465     AliCDBEntry *entry = (AliCDBEntry*)
 466     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 467
 468     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 469   }
 470
 471   if(!fPedestalData)
 472     AliFatal("Dead map not found in CDB!");
 473
 474   return kTRUE;
 475 }
 476
 477 //_____________________________________________________
 478 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
 479 {
 480   // Get the input event, it can depend in embedded events what you want to get
 481   // Also check if the quality of the event is good if not reject it
 482
 483   fEvent = 0x0;
 484
 485   //Process events if there is a high energy cluster
 486   if(!AcceptEventEMCAL())  return ;
 487
 488   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 489   Int_t eventN = Entry();
 490   if(aodIH) eventN = aodIH->GetReadEntry();
 491
 492   if (eventN > fMaxEvent)
 493     return ;
 494
 495   //printf("Clusterizer --- Event %d-- \n",eventN);
 496
 497   //Check if input event are embedded events
 498   //If so, take output event
 499   if (aodIH && aodIH->GetMergeEvents())
 500   {
 501     fEvent  = AODEvent();
 502
 503     if(!aodIH->GetMergeEMCALCells())
 504       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
 505
 506     if(DebugLevel() > 1)
 507     {
 508       printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
 509
 510       printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
 511              InputEvent()->GetEMCALCells()->GetNumberOfCells());
 512
 513       printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
 514              aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
 515
 516       for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
 517       {
 518         AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
 519         if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
 520       }
 521
 522       printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
 523              AODEvent()->GetEMCALCells()->GetNumberOfCells());
 524     }
 525   }
 526   else
 527   {
 528     fEvent =  InputEvent();
 529     if(fFillAODCaloCells) FillAODCaloCells();
 530     if(fFillAODHeader)    FillAODHeader();
 531   }
 532
 533   if (!fEvent)
 534   {
 535     Error("UserExec","Event not available");
 536     return ;
 537   }
 538
 539   //-------------------------------------------------------------------------------------
 540   // Reject events if LED was firing, use only for LHC11a data
 541   // Reject event if triggered by exotic cell and remove exotic cells if not triggered
 542   //-------------------------------------------------------------------------------------
 543
 544   if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
 545
 546   if( IsExoticEvent() )                            { fEvent = 0x0 ; return ; }
 547
 548   //-------------------------------------------------------------------------------------
 549   // Set the cluster array in the event (output or input)
 550   //-------------------------------------------------------------------------------------
 551
 552   if     ( fFillAODFile )
 553   {
 554     //Magic line to write events to AOD filem put after event rejection
 555     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 556   }
 557   else if( !fOutputAODBranchSet )
 558   {
 559     // Create array and put it in the input event, if output AOD not selected, only once
 560     InputEvent()->AddObject(fOutputAODBranch);
 561     fOutputAODBranchSet = kTRUE;
 562     printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data());
 563   }
 564
 565 }
 566
 567 //____________________________________________________
 568 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
 569 {
 570   // Recluster calocells, transform them into digits,
 571   // feed the clusterizer with them and get new list of clusters
 572
 573   //In case of MC, first loop on the clusters and fill MC label to array
 574   Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
 575   Int_t nClustersOrg  = 0;
 576
 577   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 578   if(aodIH && aodIH->GetEventToMerge())  //Embedding
 579     nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
 580
 581   ResetArrays();
 582
 583   for (Int_t i = 0; i < nClusters; i++)
 584   {
 585     AliVCluster *clus = 0;
 586     if(aodIH && aodIH->GetEventToMerge()) //Embedding
 587       clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
 588     else
 589       clus = fEvent->GetCaloCluster(i);
 590
 591     if(!clus) return;
 592
 593     if(clus->IsEMCAL())
 594     {
 595       Int_t label = clus->GetLabel();
 596       Int_t label2 = -1 ;
 597       //printf("Org cluster E %f, Time  %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label );
 598       //printf("Original list of labels from old cluster : \n");
 599       //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
 600
 601       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 602       UShort_t * index    = clus->GetCellsAbsId() ;
 603       for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
 604       {
 605         //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell]));
 606         fOrgClusterCellId[index[icell]] = i;
 607         fCellLabels[index[icell]]       = label;
 608         fCellSecondLabels[index[icell]] = label2;
 609         fCellTime[index[icell]]         = clus->GetTOF();
 610         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
 611         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
 612       }
 613       nClustersOrg++;
 614     }
 615     // printf("\n");
 616   }
 617
 618   // Transform CaloCells into Digits
 619
 620   Int_t    idigit =  0;
 621   Int_t    id     = -1;
 622   Float_t  amp    = -1;
 623   Double_t time   = -1;
 624
 625   AliVCaloCells *cells = fEvent->GetEMCALCells();
 626
 627   Int_t bc = InputEvent()->GetBunchCrossNumber();
 628
 629   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 630   {
 631     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
 632     id = cells->GetCellNumber(icell);
 633     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
 634
 635     // Do not include cells with too low energy, nor exotic cell
 636     if( amp  < fRecParam->GetMinECut() ||
 637         time > fRecParam->GetTimeMax() ||
 638         time < fRecParam->GetTimeMin()    ) accept = kFALSE;
 639
 640     // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
 641     if (fRecalibrateWithClusterTime)
 642     {
 643       time = fCellTime[id];
 644       //printf("cell %d time org %f - ",id, time*1.e9);
 645       fRecoUtils->RecalibrateCellTime(id,bc,time);
 646       //printf("recal %f\n",time*1.e9);
 647     }
 648
 649     //Exotic?
 650     if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
 651         accept = kFALSE;
 652
 653     if( !accept )
 654     {
 655       if( DebugLevel() > 2 )
 656         printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
 657                id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
 658       continue;
 659     }
 660
 661     Int_t mcLabel = cells->GetMCLabel(icell);
 662     //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel);
 663
 664     //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]);
 665     if     ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
 666     else if( fSetCellMCLabelFromCluster == 0 && fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel);
 667     else mcLabel = -1; // found later
 668
 669     //printf("\t new label %d\n",mcLabel);
 670
 671     // Create the digit, put a fake primary deposited energy to trick the clusterizer
 672     // when checking the most likely primary
 673
 674     Float_t efrac = cells->GetEFraction(icell);
 675
 676     //When checking the MC of digits, give weight to cells with embedded signal
 677     if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
 678
 679     //printf("******* Cell %d, id %d, e %f,  fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel);
 680
 681     new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
 682     // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that
 683     // we give more weight to the MC label of the cell with highest energy in the cluster
 684
 685     idigit++;
 686   }
 687
 688   fDigitsArr->Sort();
 689
 690   //-------------------------------------------------------------------------------------
 691   //Do the clusterization
 692   //-------------------------------------------------------------------------------------
 693
 694   fClusterizer->Digits2Clusters("");
 695
 696   //-------------------------------------------------------------------------------------
 697   //Transform the recpoints into AliVClusters
 698   //-------------------------------------------------------------------------------------
 699
 700   RecPoints2Clusters();
 701
 702   if(!fCaloClusterArr)
 703   {
 704     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 705     return;
 706   }
 707
 708   if( DebugLevel() > 0 )
 709   {
 710     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
 711
 712     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
 713     {
 714       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
 715              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
 716              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 717     }
 718   }
 719 }
 720
 721 //_____________________________________________________
 722 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
 723 {
 724   // Take the event clusters and unfold them
 725
 726   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 727   Double_t cellAmplitude = 0;
 728   Double_t cellTime      = 0;
 729   Short_t  cellNumber    = 0;
 730   Int_t    cellMCLabel   = 0;
 731   Double_t cellEFrac     = 0;
 732   Int_t    nClustersOrg  = 0;
 733
 734   // Fill the array with the EMCAL clusters, copy them
 735   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
 736   {
 737     AliVCluster *clus = fEvent->GetCaloCluster(i);
 738     if(clus->IsEMCAL())
 739     {
 740       //recalibrate/remove bad channels/etc if requested
 741       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
 742       {
 743         continue;
 744       }
 745
 746       if(fRecoUtils->IsRecalibrationOn())
 747       {
 748         //Calibrate cluster
 749         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 750
 751         //CalibrateCells
 752         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 753         {
 754           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
 755             break;
 756
 757           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 758           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 759           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 760
 761           //Do not include bad channels found in analysis?
 762           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 763               fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
 764             continue;
 765
 766           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 767
 768         }// cells loop
 769       }// recalibrate
 770
 771       //Cast to ESD or AOD, needed to create the cluster array
 772       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 773       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 774
 775       if     (esdCluster)
 776       {
 777         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 778       }//ESD
 779       else if(aodCluster)
 780       {
 781         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 782       }//AOD
 783       else
 784         Warning("UserExec()"," - Wrong CaloCluster type?");
 785
 786       nClustersOrg++;
 787     }
 788   }
 789
 790   //Do the unfolding
 791   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 792
 793   //CLEAN-UP
 794   fUnfolder->Clear();
 795
 796 }
 797
 798 //_____________________________________________________
 799 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
 800 {
 801   // Put calo cells in standard branch
 802   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
 803   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 804
 805   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 806   aodEMcells.CreateContainer(nEMcell);
 807   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 808   Double_t calibFactor = 1.;
 809   for (Int_t iCell = 0; iCell < nEMcell; iCell++)
 810   {
 811     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 812     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 813     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 814
 815     if(fRecoUtils->IsRecalibrationOn())
 816     {
 817       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 818     }
 819
 820     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
 821     { //Channel is not declared as bad
 822       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
 823                          eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
 824     }
 825     else
 826     {
 827       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
 828     }
 829   }
 830   aodEMcells.Sort();
 831
 832 }
 833
 834 //__________________________________________________
 835 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
 836 {
 837   //Put event header information in standard AOD branch
 838
 839   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
 840   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
 841
 842   Double_t pos[3]   ;
 843   Double_t covVtx[6];
 844   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 845
 846   AliAODHeader* header = AODEvent()->GetHeader();
 847   header->SetRunNumber(fEvent->GetRunNumber());
 848
 849   if(esdevent)
 850   {
 851     TTree* tree = fInputHandler->GetTree();
 852     if (tree)
 853     {
 854       TFile* file = tree->GetCurrentFile();
 855       if (file) header->SetESDFileName(file->GetName());
 856     }
 857   }
 858   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 859
 860   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
 861   header->SetOrbitNumber(fEvent->GetOrbitNumber());
 862   header->SetPeriodNumber(fEvent->GetPeriodNumber());
 863   header->SetEventType(fEvent->GetEventType());
 864
 865   //Centrality
 866   if(fEvent->GetCentrality())
 867     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
 868   else
 869     header->SetCentrality(0);
 870
 871   //Trigger
 872   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 873   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 874   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 875
 876   header->SetTriggerMask(fEvent->GetTriggerMask());
 877   header->SetTriggerCluster(fEvent->GetTriggerCluster());
 878
 879   if      (esdevent)
 880   {
 881     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 882     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 883     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 884   }
 885   else if (aodevent)
 886   {
 887     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 888     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 889     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 890   }
 891
 892   header->SetMagneticField(fEvent->GetMagneticField());
 893   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 894
 895   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 896   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 897   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 898   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 899   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
 900
 901   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
 902   Float_t diamcov[3];
 903   fEvent->GetDiamondCovXY(diamcov);
 904   header->SetDiamond(diamxy,diamcov);
 905   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 906   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 907
 908   //
 909   Int_t nVertices = 1 ;/* = prim. vtx*/;
 910   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
 911
 912   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 913
 914   // Access to the AOD container of vertices
 915   TClonesArray &vertices = *(AODEvent()->GetVertices());
 916   Int_t jVertices=0;
 917
 918   // Add primary vertex. The primary tracks will be defined
 919   // after the loops on the composite objects (V0, cascades, kinks)
 920   fEvent->GetPrimaryVertex()->GetXYZ(pos);
 921   Float_t chi = 0;
 922   if      (esdevent)
 923   {
 924     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 925     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 926   }
 927   else if (aodevent)
 928   {
 929     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 930     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 931   }
 932
 933   AliAODVertex * primary = new(vertices[jVertices++])
 934   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 935   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
 936   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
 937
 938 }
 939
 940 //___________________________________________________________
 941 void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent()
 942 {
 943   // Get the CaloClusters array, do some final calculations
 944   // and put the clusters in the output or input event
 945   // as a separate branch
 946
 947   //First recalculate track-matching for the new clusters
 948   if(fDoTrackMatching)
 949   {
 950     fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
 951   }
 952   //Put the new clusters in the AOD list
 953
 954   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
 955
 956   for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
 957   {
 958     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 959
 960     newCluster->SetID(i);
 961
 962     // Correct cluster energy non linearity
 963     newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
 964
 965     //Add matched track
 966     if(fDoTrackMatching)
 967     {
 968       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 969       if(trackIndex >= 0)
 970       {
 971         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
 972         if(DebugLevel() > 1)
 973           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
 974       }
 975
 976       Float_t dR = 999., dZ = 999.;
 977       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
 978       newCluster->SetTrackDistance(dR,dZ);
 979
 980     }
 981     else
 982     {// Assign previously assigned matched track in reco, very very rough
 983       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
 984       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
 985     }
 986
 987     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
 988     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
 989     //printf("\n");
 990
 991     // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
 992     fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
 993
 994     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 995
 996     if(DebugLevel() > 1 )
 997       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f, mc label %d \n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel());
 998
 999   } // cluster loop
1000
1001   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1002
1003
1004 }
1005
1006 //_______________________________________________
1007 TString AliAnalysisTaskEMCALClusterize::GetPass()
1008 {
1009   // Get passx from filename.
1010
1011   if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1012   {
1013     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
1014     return TString("");
1015   }
1016
1017   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1018   {
1019     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
1020     return TString("");
1021   }
1022
1023   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1024   if      (pass.Contains("ass1")) return TString("pass1");
1025   else if (pass.Contains("ass2")) return TString("pass2");
1026   else if (pass.Contains("ass3")) return TString("pass3");
1027   else if (pass.Contains("ass4")) return TString("pass4");
1028   else if (pass.Contains("ass5")) return TString("pass5");
1029
1030   // No condition fullfilled, give a default value
1031   printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
1032   return TString("");
1033
1034 }
1035
1036 //_________________________________________
1037 void AliAnalysisTaskEMCALClusterize::Init()
1038 {
1039   //Init analysis with configuration macro if available
1040
1041   fOADBSet           = kFALSE;
1042   if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
1043
1044   fBranchNames       = "ESD:AliESDHeader.,EMCALCells.";
1045
1046   if(!fRecParam)     fRecParam  = new AliEMCALRecParam;
1047   if(!fRecoUtils)    fRecoUtils = new AliEMCALRecoUtils();
1048
1049   if(fMaxEvent          <= 0) fMaxEvent          = 1000000000;
1050   if(fSelectCellMinE    <= 0) fSelectCellMinE    = 0.005;
1051   if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1052
1053   //Centrality
1054   if(fCentralityClass  == "") fCentralityClass  = "V0M";
1055
1056   if (fOCDBpath            == "") fOCDBpath            = "raw://" ;
1057   if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1058
1059   if(gROOT->LoadMacro(fConfigName) >=0)
1060   {
1061     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
1062     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1063     fGeomName         = clus->fGeomName;
1064     fLoadGeomMatrices = clus->fLoadGeomMatrices;
1065     fOCDBpath         = clus->fOCDBpath;
1066     fAccessOCDB       = clus->fAccessOCDB;
1067     fRecParam         = clus->fRecParam;
1068     fJustUnfold       = clus->fJustUnfold;
1069     fFillAODFile      = clus->fFillAODFile;
1070     fRecoUtils        = clus->fRecoUtils;
1071     fConfigName       = clus->fConfigName;
1072     fMaxEvent         = clus->fMaxEvent;
1073     fDoTrackMatching  = clus->fDoTrackMatching;
1074     fOutputAODBranchName = clus->fOutputAODBranchName;
1075     for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1076     fCentralityClass  = clus->fCentralityClass;
1077     fCentralityBin[0] = clus->fCentralityBin[0];
1078     fCentralityBin[1] = clus->fCentralityBin[1];
1079   }
1080
1081   // Init geometry, I do not like much to do it like this ...
1082   if(fImportGeometryFromFile && !gGeoManager)
1083   {
1084     if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1085     printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
1086     TGeoManager::Import(fImportGeometryFilePath) ;
1087   }
1088
1089 }
1090
1091 //_______________________________________________________
1092 void AliAnalysisTaskEMCALClusterize::InitClusterization()
1093 {
1094   //Select clusterization/unfolding algorithm and set all the needed parameters
1095
1096   if (fJustUnfold)
1097   {
1098     // init the unfolding afterburner
1099     delete fUnfolder;
1100     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1101     return;
1102   }
1103
1104   //First init the clusterizer
1105   delete fClusterizer;
1106   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1107     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
1108   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1109     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
1110   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1111   {
1112     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
1113     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1114     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1115   }
1116   else
1117   {
1118     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1119   }
1120
1121   //Now set the parameters
1122   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1123   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1124   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
1125   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1126   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1127   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1128   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1129   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1130   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1131   fClusterizer->SetJustClusters          ( kTRUE                               );
1132
1133   // Initialize the cluster rec points and digits arrays and get them.
1134   fClusterizer->SetOutput(0);
1135   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1136   fDigitsArr  = fClusterizer->GetDigits();
1137
1138   //In case of unfolding after clusterization is requested, set the corresponding parameters
1139   if(fRecParam->GetUnfold())
1140   {
1141     Int_t i=0;
1142     for (i = 0; i < 8; i++)
1143     {
1144       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1145     }//end of loop over parameters
1146
1147     for (i = 0; i < 3; i++)
1148     {
1149       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1150       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1151     }//end of loop over parameters
1152
1153     fClusterizer->InitClusterUnfolding();
1154
1155   }// to unfold
1156 }
1157
1158 //_________________________________________________
1159 void AliAnalysisTaskEMCALClusterize::InitGeometry()
1160 {
1161   // Init geometry and set the geometry matrix, for the first event, skip the rest
1162   // Also set once the run dependent calibrations
1163
1164   if(fGeomMatrixSet) return;
1165
1166   Int_t runnumber = InputEvent()->GetRunNumber() ;
1167   if (!fGeom)
1168   {
1169     if(fGeomName=="")
1170     {
1171       if     (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
1172       else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
1173       else                        fGeomName = "EMCAL_COMPLETE12SMV1";
1174       printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
1175     }
1176
1177                 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1178
1179                 if(fDebug > 0)
1180     {
1181                         printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
1182                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1183                         printf("\n");
1184                 }
1185         } // geometry pointer did not exist before
1186
1187   if(fLoadGeomMatrices)
1188   {
1189     printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
1190
1191     // OADB if available
1192     AliOADBContainer emcGeoMat("AliEMCALgeo");
1193     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1194     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1195
1196
1197     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1198     {
1199
1200       if (!fGeomMatrix[mod]) // Get it from OADB
1201       {
1202         if(fDebug > 1 )
1203           printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1204                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1205         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1206
1207         fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1208       }
1209
1210       if(fGeomMatrix[mod])
1211       {
1212         if(DebugLevel() > 1)
1213           fGeomMatrix[mod]->Print();
1214
1215         fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1216       }
1217
1218       fGeomMatrixSet=kTRUE;
1219
1220     }//SM loop
1221   }//Load matrices
1222   else if(!gGeoManager)
1223   {
1224     printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1225     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1226     if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1227     {
1228       if(DebugLevel() > 1)
1229         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1230     }//AOD
1231     else
1232     {
1233       if(DebugLevel() > 1)
1234         printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1235
1236       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1237
1238       if(!esd)
1239       {
1240         Error("InitGeometry"," - This event does not contain ESDs?");
1241         if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1242         return;
1243       }
1244
1245       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1246       {
1247         if(DebugLevel() > 1)
1248           esd->GetEMCALMatrix(mod)->Print();
1249
1250         if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1251
1252       }
1253
1254       fGeomMatrixSet=kTRUE;
1255
1256     }//ESD
1257   }//Load matrices from Data
1258
1259 }
1260
1261 //____________________________________________________
1262 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1263 {
1264
1265   // Check if event is exotic, get an exotic cell and compare with triggered patch
1266   // If there is a match remove event ... to be completed, filled with something provisional
1267
1268   if(!fRemoveExoticEvents) return kFALSE;
1269
1270   // Loop on cells
1271   AliVCaloCells * cells = fEvent->GetEMCALCells();
1272   Float_t totCellE = 0;
1273   Int_t bc = InputEvent()->GetBunchCrossNumber();
1274   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1275   {
1276
1277     Float_t  ecell = 0 ;
1278     Double_t tcell = 0 ;
1279
1280     Int_t absID   = cells->GetCellNumber(icell);
1281     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1282     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1283   }
1284
1285   //  TString triggerclasses = "";
1286   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
1287   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
1288   //    //
1289   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1290   //    if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1291   //    return;
1292   //
1293
1294   //printf("TotE cell %f\n",totCellE);
1295   if(totCellE < 1) return kTRUE;
1296
1297   return kFALSE;
1298
1299 }
1300
1301 //________________________________________________________________
1302 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1303 {
1304   //Check if event is LED
1305
1306   if(!fRemoveLEDEvents) return kFALSE;
1307
1308   //check events of LHC11a period
1309   if(run < 146858 || run > 146860) return kFALSE ;
1310
1311   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1312   Int_t ncellsSM3 = 0;
1313   AliVCaloCells * cells = fEvent->GetEMCALCells();
1314   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1315   {
1316     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1317   }
1318
1319   TString triggerclasses = "";
1320
1321   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
1322   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
1323   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
1324
1325   Int_t ncellcut = 21;
1326   if(triggerclasses.Contains("EMC")) ncellcut = 35;
1327
1328   if( ncellsSM3 >= ncellcut)
1329   {
1330     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d\n",(Int_t)Entry(),ncellsSM3);
1331     if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1332     return kTRUE;
1333   }
1334
1335   return kFALSE;
1336
1337 }
1338
1339 //_______________________________________________________
1340 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
1341 {
1342   // Restore clusters from recPoints
1343   // Cluster energy, global position, cells and their amplitude
1344   // fractions are restored
1345
1346   Int_t j = 0;
1347   for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1348   {
1349     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1350
1351     const Int_t ncells = recPoint->GetMultiplicity();
1352     Int_t ncellsTrue = 0;
1353
1354     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1355
1356     // cells and their amplitude fractions
1357     UShort_t   absIds[ncells];
1358     Double32_t ratios[ncells];
1359
1360     //For later check embedding
1361     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1362
1363     Float_t clusterE = 0;
1364     for (Int_t c = 0; c < ncells; c++)
1365     {
1366       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1367
1368       absIds[ncellsTrue] = digit->GetId();
1369       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1370
1371       // In case of unfolding, remove digits with unfolded energy too low
1372       if(fSelectCell)
1373       {
1374         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1375         {
1376           if(DebugLevel() > 1)
1377           {
1378             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1379                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1380           }
1381
1382           continue;
1383
1384         } // if cuts
1385       }// Select cells
1386
1387       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1388       clusterE  +=recPoint->GetEnergiesList()[c];
1389
1390       // In case of embedding, fill ratio with amount of signal,
1391       if (aodIH && aodIH->GetMergeEvents())
1392       {
1393         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1394         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1395         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1396
1397         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1398         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1399         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1400         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1401
1402         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1403         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1404
1405       }//Embedding
1406
1407       ncellsTrue++;
1408
1409     }// cluster cell loop
1410
1411     if (ncellsTrue < 1)
1412     {
1413       if (DebugLevel() > 1)
1414         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1415                recPoint->GetEnergy(), ncells);
1416       continue;
1417     }
1418
1419     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1420
1421     if(clusterE <  fRecParam->GetClusteringThreshold())
1422     {
1423       if (DebugLevel()>1)
1424         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1425       continue;
1426     }
1427
1428     TVector3 gpos;
1429     Float_t g[3];
1430
1431     // calculate new cluster position
1432
1433     recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1434     recPoint->GetGlobalPosition(gpos);
1435     gpos.GetXYZ(g);
1436
1437     // create a new cluster
1438
1439     (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1440     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1441     j++;
1442     clus->SetType(AliVCluster::kEMCALClusterv1);
1443     clus->SetE(clusterE);
1444     clus->SetPosition(g);
1445     clus->SetNCells(ncellsTrue);
1446     clus->SetCellsAbsId(absIds);
1447     clus->SetCellsAmplitudeFraction(ratios);
1448     clus->SetChi2(-1); //not yet implemented
1449     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1450     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1451     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1452
1453     if(ncells == ncellsTrue)
1454     {
1455       Float_t elipAxis[2];
1456       recPoint->GetElipsAxis(elipAxis);
1457       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1458       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1459       clus->SetDispersion(recPoint->GetDispersion());
1460     }
1461     else if(fSelectCell)
1462     {
1463       // In case some cells rejected, in unfolding case, recalculate
1464       // shower shape parameters and position
1465       if(DebugLevel() > 1)
1466         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1467
1468       AliVCaloCells* cells = 0x0;
1469       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
1470       else                                  cells = InputEvent()->GetEMCALCells();
1471
1472       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1473       fRecoUtils->RecalculateClusterPID(clus);
1474       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1475
1476     }
1477
1478     // MC
1479
1480     if     ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ;
1481     else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ;
1482     else
1483     {
1484       // Normal case, trust what the clusterizer has found
1485       Int_t  parentMult = 0;
1486       Int_t *parentList = recPoint->GetParents(parentMult);
1487       clus->SetLabel(parentList, parentMult);
1488 //      printf("Label list : ");
1489 //      for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]);
1490 //      printf("\n");
1491     }
1492
1493   } // recPoints loop
1494
1495 }
1496
1497 //___________________________________________________________________________
1498 void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label)
1499 {
1500   // MC label for Cells not remapped after ESD filtering, do it here.
1501
1502   if(label < 0) return ;
1503
1504   AliAODEvent  * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1505   if(!evt) return ;
1506
1507   TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1508   if(!arr) return ;
1509
1510   if(label < arr->GetEntriesFast())
1511   {
1512     AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1513     if(!particle) return ;
1514
1515     if(label == particle->Label()) return ; // label already OK
1516     //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label  %d - AOD stack %d \n",label, particle->Label());
1517   }
1518   //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label  %d > AOD labels %d \n",label, arr->GetEntriesFast());
1519
1520   // loop on the particles list and check if there is one with the same label
1521   for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1522   {
1523     AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1524     if(!particle) continue ;
1525
1526     if(label == particle->Label())
1527     {
1528       label = ind;
1529       //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index  %d \n",label);
1530       return;
1531     }
1532   }
1533
1534   label = -1;
1535
1536   //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1537
1538 }
1539
1540 //________________________________________________
1541 void AliAnalysisTaskEMCALClusterize::ResetArrays()
1542 {
1543   // Reset arrays containing information for all possible cells
1544   for(Int_t j = 0; j < 12672; j++)
1545   {
1546     fOrgClusterCellId[j] = -1;
1547     fCellLabels[j]       = -1;
1548     fCellSecondLabels[j] = -1;
1549     fCellTime[j]         =  0.;
1550     fCellMatchdEta[j]    = -999;
1551     fCellMatchdPhi[j]    = -999;
1552   }
1553 }
1554
1555 //_____________________________________________________________________________________________________
1556 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint  * recPoint,
1557                                                                            AliAODCaloCluster * clus)
1558 {
1559   // Set the cluster MC label, the digizer was filled with most likely MC label for all cells in original cluster
1560   // Now check the second most likely MC label and add it to the new cluster
1561
1562   Int_t  parentMult = 0;
1563   Int_t *parentList = recPoint->GetParents(parentMult);
1564   clus->SetLabel(parentList, parentMult);
1565
1566   //Write the second major contributor to each MC cluster.
1567   Int_t iNewLabel ;
1568   for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1569   {
1570
1571     Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1572     if(idCell>=0)
1573     {
1574       iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1575       for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1576       {
1577         if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1578           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1579             }
1580       if (iNewLabel == 1)
1581       {
1582         Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1583         for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1584         {
1585           newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1586         }
1587
1588         newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1589         clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1590         delete [] newLabelArray;
1591       }
1592     }//positive cell number
1593   }
1594 }
1595
1596 //___________________________________________________________________________________________________
1597 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus)
1598 {
1599   // Get the original clusters that contribute to the new cluster, assign the labels of such clusters
1600   // to the new cluster.
1601   // Only approximatedly valid  when output are V1 clusters, handle with care
1602
1603   TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1604   clArray.Reset();
1605   Int_t nClu = 0;
1606   Int_t nLabTotOrg = 0;
1607   Float_t emax = -1;
1608   Int_t idMax = -1;
1609
1610   AliVEvent * event = fEvent;
1611
1612   //In case of embedding the MC signal is in the added event
1613   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1614   if(aodIH && aodIH->GetEventToMerge())  //Embedding
1615     event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1616
1617
1618   //Find the clusters that originally had the cells
1619   for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1620   {
1621     Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1622
1623     if(idCell>=0)
1624     {
1625       Int_t idCluster = fOrgClusterCellId[idCell];
1626
1627       Bool_t set = kTRUE;
1628       for(Int_t icl =0; icl < nClu; icl++)
1629       {
1630         if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1631         if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1632         //   printf("\t \t icell %d  Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1633         //          iLoopCell,                 icl,  idCluster,((Int_t)clArray.GetAt(icl)),set);
1634       }
1635       if( set && idCluster >= 0)
1636       {
1637         clArray.SetAt(idCluster,nClu++);
1638         //printf("******** idCluster %d \n",idCluster);
1639         nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1640
1641         //printf("Cluster in array %d, IdCluster %d\n",nClu-1,  idCluster);
1642
1643         //Search highest E cluster
1644         AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1645         //printf("\t E %f\n",clOrg->E());
1646         if(emax < clOrg->E())
1647         {
1648           emax  = clOrg->E();
1649           idMax = idCluster;
1650         }
1651       }
1652     }
1653   }// cell loop
1654
1655   if(nClu==0 || nLabTotOrg == 0)
1656   {
1657     //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f,  n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
1658     //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1659     //printf("\n");
1660   }
1661
1662   // Put the first in the list the cluster with highest energy
1663   if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1664   {
1665     Int_t maxIndex = -1;
1666     Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1667     for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1668     {
1669       if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1670     }
1671
1672     if(firstCluster >=0 && idMax >=0)
1673     {
1674       clArray.SetAt(idMax,0);
1675       clArray.SetAt(firstCluster,maxIndex);
1676     }
1677   }
1678
1679   // Get the labels list in the original clusters, assign all to the new cluster
1680   TArrayI clMCArray(nLabTotOrg) ;
1681   clMCArray.Reset();
1682
1683   Int_t nLabTot = 0;
1684   for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1685   {
1686     Int_t idCluster = clArray.GetAt(iLoopCluster);
1687     //printf("New Cluster in Array %d,  idCluster %d \n",iLoopCluster,idCluster);
1688     AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1689     Int_t nLab = clOrg->GetNLabels();
1690
1691     for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1692     {
1693       Int_t lab = clOrg->GetLabelAt(iLab) ;
1694       if(lab>=0)
1695       {
1696         Bool_t set = kTRUE;
1697         //printf("\t \t Set Label %d \n", lab);
1698         for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1699         {
1700           if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1701           //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d  Label ID in array %d, label in array %d, set %d\n",
1702           //       iLoopCluster,                           iLab,     lab,           iLabTot,  ((Int_t)clMCArray.GetAt(iLabTot)),set);
1703         }
1704         if( set ) clMCArray.SetAt(lab,nLabTot++);
1705       }
1706     }
1707   }// cluster loop
1708
1709   // Set the final list of labels
1710
1711   clus->SetLabel(clMCArray.GetArray(), nLabTot);
1712
1713 //  printf("Final list of labels for new cluster : \n");
1714 //  for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1715 //  {
1716 //    printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1717 //    Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1718 //    printf(" org %d ",label);
1719 //    RemapMCLabelForAODs(label);
1720 //    printf(" new %d \n",label);
1721 //  }
1722 //  for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1723 }
1724
1725
1726 //____________________________________________________________
1727 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1728 {
1729   // Init geometry, create list of output clusters
1730
1731
1732   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1733
1734   if(fOutputAODBranchName.Length()==0)
1735   {
1736     fOutputAODBranchName = "newEMCALClustersArray";
1737     printf("Cluster branch name not set, set it to newEMCALClustersArray \n");
1738   }
1739
1740   fOutputAODBranch->SetName(fOutputAODBranchName);
1741
1742   if( fFillAODFile )
1743   {
1744     //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1745
1746
1747     //fOutputAODBranch->SetOwner(kFALSE);
1748
1749     AddAODBranch("TClonesArray", &fOutputAODBranch);
1750   }
1751
1752 }
1753
1754 //_______________________________________________________
1755 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1756 {
1757   // Do clusterization event by event, execute different steps
1758   // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed, initializations
1759   //    load things and clear arrays
1760   // 2) Clusterize a) just unfolding existing clusters (fJustUnfold)
1761   //               b) recluster cells
1762   //                   +  convert cells into digits (calibrating them)
1763   //                   +  recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN)
1764   //                      with methods in AliEMCALClusterizer
1765   //                   +  transform recPoints into CaloClusters
1766   // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch
1767
1768   //-------
1769   // Step 1
1770
1771
1772   //Remove the contents of AOD branch output list set in the previous event
1773   fOutputAODBranch->Clear("C");
1774
1775   LoadBranches();
1776
1777   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
1778   //If we need a centrality bin, we select only those events in the corresponding bin.
1779   if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1780   {
1781     Float_t cen = GetEventCentrality();
1782     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
1783   }
1784
1785   // intermediate array with new clusters : init the array only once or clear from previous event
1786   if(!fCaloClusterArr) fCaloClusterArr    = new TObjArray(10000);
1787   else                 fCaloClusterArr->Delete();//Clear("C"); it leaks?
1788
1789   InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1790
1791   //Get the event, do some checks and settings
1792   CheckAndGetEvent() ;
1793
1794   if (!fEvent)
1795   {
1796     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1797     return ;
1798   }
1799
1800   //Init pointers, geometry, clusterizer, ocdb, aodb
1801
1802
1803   if(fAccessOCDB) AccessOCDB();
1804   if(fAccessOADB) AccessOADB(); // only once
1805
1806   InitClusterization();
1807
1808   //-------
1809   // Step 2
1810
1811   // Make clusters
1812   if (fJustUnfold) ClusterUnfolding();
1813   else             ClusterizeCells() ;
1814
1815   //-------
1816   // Step 3
1817
1818   FillCaloClusterInEvent();
1819
1820 }
1821