PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx

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