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