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       //printf("Org cluster E %f, Time  %e, Id = ", clus->E(), clus->GetTOF() );
 336       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 337       UShort_t * index    = clus->GetCellsAbsId() ;
 338       for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
 339         fCellLabels[index[icell]]       = label;
 340         fCellSecondLabels[index[icell]] = label2;
 341         fCellTime[index[icell]]         = clus->GetTOF();
 342         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
 343         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
 344         //printf(" %d,", index[icell] );
 345       }
 346       nClustersOrg++;
 347     }
 348      // printf("\n");
 349   }
 350
 351   // Transform CaloCells into Digits
 352
 353   Int_t    idigit =  0;
 354   Int_t    id     = -1;
 355   Float_t  amp    = -1;
 356   Double_t time   = -1;
 357
 358   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 359
 360   TTree *digitsTree = new TTree("digitstree","digitstree");
 361   digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 362   Int_t bc = InputEvent()->GetBunchCrossNumber();
 363   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 364   {
 365
 366     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
 367     id = cells->GetCellNumber(icell);
 368     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
 369
 370     // Do not include cells with too low energy, nor exotic cell
 371     if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
 372
 373     // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
 374     if (time*1e9 < 1.) {
 375       time = fCellTime[id];
 376       //printf("cell %d time cluster %e\n",id, time);
 377       fRecoUtils->RecalibrateCellTime(id,bc,time);
 378     }
 379
 380     if(  accept && fRecoUtils->IsExoticCell(id,cells,bc)){
 381       accept = kFALSE;
 382     }
 383
 384     if( !accept ){
 385       fCellLabels[id]      =-1; //reset the entry in the array for next event
 386       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
 387       fCellTime[id]        = 0.;
 388       fCellMatchdEta[id]   =-999;
 389       fCellMatchdPhi[id]   =-999;
 390       if( DebugLevel() > 2 ) printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
 391                                     id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
 392       continue;
 393     }
 394
 395     //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 396     new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 397
 398     fCellLabels[id]      =-1; //reset the entry in the array for next event
 399
 400     idigit++;
 401   }
 402
 403   //Fill the tree with digits
 404   digitsTree->Fill();
 405
 406   //-------------------------------------------------------------------------------------
 407   //Do the clusterization
 408   //-------------------------------------------------------------------------------------
 409   TTree *clustersTree = new TTree("clustertree","clustertree");
 410
 411   fClusterizer->SetInput(digitsTree);
 412   fClusterizer->SetOutput(clustersTree);
 413   fClusterizer->Digits2Clusters("");
 414
 415   //-------------------------------------------------------------------------------------
 416   //Transform the recpoints into AliVClusters
 417   //-------------------------------------------------------------------------------------
 418
 419   clustersTree->SetBranchStatus("*",0); //disable all branches
 420   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 421
 422   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 423   branch->SetAddress(&fClusterArr);
 424   branch->GetEntry(0);
 425
 426   RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 427
 428   if(!fCaloClusterArr){
 429     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 430     return;
 431   }
 432
 433   if( DebugLevel() > 0 ){
 434
 435     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
 436
 437     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
 438       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
 439              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
 440              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 441     }
 442   }
 443
 444   //Reset the array with second labels for this event
 445   memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
 446
 447   //---CLEAN UP-----
 448   fClusterizer->Clear();
 449   fDigitsArr  ->Clear("C");
 450   fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 451
 452   clustersTree->Delete("all");
 453   digitsTree  ->Delete("all");
 454
 455 }
 456
 457 //_____________________________________________________
 458 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
 459 {
 460   // Take the event clusters and unfold them
 461
 462   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 463   Double_t cellAmplitude = 0;
 464   Double_t cellTime      = 0;
 465   Short_t  cellNumber    = 0;
 466   Int_t    nClustersOrg  = 0;
 467
 468   // Fill the array with the EMCAL clusters, copy them
 469   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
 470   {
 471     AliVCluster *clus = fEvent->GetCaloCluster(i);
 472     if(clus->IsEMCAL()){
 473
 474       //recalibrate/remove bad channels/etc if requested
 475       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
 476         continue;
 477       }
 478
 479       if(fRecoUtils->IsRecalibrationOn()){
 480
 481         //Calibrate cluster
 482         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 483
 484         //CalibrateCells
 485         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 486         {
 487           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 488             break;
 489
 490           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 491           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 492           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 493
 494           //Do not include bad channels found in analysis?
 495           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 496              fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 497             fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
 498             fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
 499             fCellTime[cellNumber]        = 0.;
 500             fCellMatchdEta[cellNumber]   =-999;
 501             fCellMatchdPhi[cellNumber]   =-999;
 502             continue;
 503           }
 504
 505           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 506
 507         }// cells loop
 508       }// recalibrate
 509
 510       //Cast to ESD or AOD, needed to create the cluster array
 511       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 512       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 513       if     (esdCluster){
 514         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 515       }//ESD
 516       else if(aodCluster){
 517         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 518       }//AOD
 519       else
 520         Warning("UserExec()"," - Wrong CaloCluster type?");
 521       nClustersOrg++;
 522     }
 523   }
 524
 525   //Do the unfolding
 526   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 527
 528   //CLEAN-UP
 529   fUnfolder->Clear();
 530
 531 }
 532
 533 //_____________________________________________________
 534 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
 535 {
 536   // Put calo cells in standard branch
 537   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
 538   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 539
 540   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 541   aodEMcells.CreateContainer(nEMcell);
 542   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 543   Double_t calibFactor = 1.;
 544   for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 545     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 546     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 547     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 548
 549     if(fRecoUtils->IsRecalibrationOn()){
 550       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 551     }
 552
 553     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 554       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 555     }
 556     else {
 557       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 558     }
 559   }
 560   aodEMcells.Sort();
 561
 562 }
 563
 564 //__________________________________________________
 565 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
 566 {
 567   //Put event header information in standard AOD branch
 568
 569   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
 570   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
 571
 572   Double_t pos[3]   ;
 573   Double_t covVtx[6];
 574   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 575
 576   AliAODHeader* header = AODEvent()->GetHeader();
 577   header->SetRunNumber(fEvent->GetRunNumber());
 578
 579   if(esdevent){
 580     TTree* tree = fInputHandler->GetTree();
 581     if (tree) {
 582       TFile* file = tree->GetCurrentFile();
 583       if (file) header->SetESDFileName(file->GetName());
 584     }
 585   }
 586   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 587
 588   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
 589   header->SetOrbitNumber(fEvent->GetOrbitNumber());
 590   header->SetPeriodNumber(fEvent->GetPeriodNumber());
 591   header->SetEventType(fEvent->GetEventType());
 592
 593   //Centrality
 594   if(fEvent->GetCentrality()){
 595     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
 596   }
 597   else{
 598     header->SetCentrality(0);
 599   }
 600
 601   //Trigger
 602   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 603   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 604   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 605   header->SetTriggerMask(fEvent->GetTriggerMask());
 606   header->SetTriggerCluster(fEvent->GetTriggerCluster());
 607   if(esdevent){
 608     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 609     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 610     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 611   }
 612   else if (aodevent){
 613     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 614     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 615     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 616   }
 617
 618   header->SetMagneticField(fEvent->GetMagneticField());
 619   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 620
 621   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 622   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 623   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 624   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 625   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
 626
 627   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
 628   Float_t diamcov[3];
 629   fEvent->GetDiamondCovXY(diamcov);
 630   header->SetDiamond(diamxy,diamcov);
 631   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 632   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 633
 634   //
 635   //
 636   Int_t nVertices = 1 ;/* = prim. vtx*/;
 637   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
 638
 639   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 640
 641   // Access to the AOD container of vertices
 642   TClonesArray &vertices = *(AODEvent()->GetVertices());
 643   Int_t jVertices=0;
 644
 645   // Add primary vertex. The primary tracks will be defined
 646   // after the loops on the composite objects (V0, cascades, kinks)
 647   fEvent->GetPrimaryVertex()->GetXYZ(pos);
 648   Float_t chi = 0;
 649   if      (esdevent){
 650     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 651     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 652   }
 653   else if (aodevent){
 654     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 655     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 656   }
 657
 658   AliAODVertex * primary = new(vertices[jVertices++])
 659   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 660   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
 661   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
 662
 663 }
 664
 665 //_________________________________________
 666 void AliAnalysisTaskEMCALClusterize::Init()
 667 {
 668   //Init analysis with configuration macro if available
 669
 670   if(gROOT->LoadMacro(fConfigName) >=0){
 671     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
 672     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 673     fGeomName         = clus->fGeomName;
 674     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 675     fOCDBpath         = clus->fOCDBpath;
 676     fAccessOCDB       = clus->fAccessOCDB;
 677     fRecParam         = clus->fRecParam;
 678     fJustUnfold       = clus->fJustUnfold;
 679     fFillAODFile      = clus->fFillAODFile;
 680     fRecoUtils        = clus->fRecoUtils;
 681     fConfigName       = clus->fConfigName;
 682     fMaxEvent         = clus->fMaxEvent;
 683     fDoTrackMatching  = clus->fDoTrackMatching;
 684     fOutputAODBranchName = clus->fOutputAODBranchName;
 685     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
 686
 687   }
 688
 689   // Init geometry, I do not like much to do it like this ...
 690   if(fImportGeometryFromFile && !gGeoManager) {
 691     printf("AliAnalysisTaskEMCALClusterize::Init() - Import geometry.root file\n");
 692     TGeoManager::Import(Form("%s/geometry.root", fImportGeometryFilePath.Data())) ; // default need file "geometry.root" in local dir!!!!
 693   }
 694
 695 }
 696
 697 //_______________________________________________________
 698 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 699 {
 700   //Select clusterization/unfolding algorithm and set all the needed parameters
 701
 702   if (fJustUnfold){
 703     // init the unfolding afterburner
 704     delete fUnfolder;
 705     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
 706     return;
 707   }
 708
 709   //First init the clusterizer
 710   delete fClusterizer;
 711   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 712     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 713   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 714     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
 715   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){
 716     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 717     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
 718     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
 719   } else {
 720     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 721   }
 722
 723   //Now set the parameters
 724   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 725   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 726   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 727   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 728   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 729   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 730   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 731   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 732   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 733   fClusterizer->SetJustClusters          ( kTRUE                               );
 734   //In case of unfolding after clusterization is requested, set the corresponding parameters
 735   if(fRecParam->GetUnfold()){
 736     Int_t i=0;
 737     for (i = 0; i < 8; i++) {
 738       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 739     }//end of loop over parameters
 740     for (i = 0; i < 3; i++) {
 741       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 742       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 743     }//end of loop over parameters
 744
 745     fClusterizer->InitClusterUnfolding();
 746
 747   }// to unfold
 748 }
 749
 750 //_________________________________________________
 751 void AliAnalysisTaskEMCALClusterize::InitGeometry()
 752 {
 753   // Set the geometry matrix, for the first event, skip the rest
 754   // Also set once the run dependent calibrations
 755
 756   if(!fGeomMatrixSet){
 757     if(fLoadGeomMatrices){
 758       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 759         if(fGeomMatrix[mod]){
 760           if(DebugLevel() > 1)
 761             fGeomMatrix[mod]->Print();
 762           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
 763         }
 764         fGeomMatrixSet=kTRUE;
 765       }//SM loop
 766     }//Load matrices
 767     else if(!gGeoManager){
 768       printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
 769       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 770       if(!strcmp(fEvent->GetName(),"AliAODEvent")) {
 771         if(DebugLevel() > 1)
 772           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 773       }//AOD
 774       else {
 775         if(DebugLevel() > 1)
 776           printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 777         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
 778         if(!esd) {
 779           Error("InitGeometry"," - This event does not contain ESDs?");
 780           AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
 781           return;
 782         }
 783         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 784           if(DebugLevel() > 1)
 785             esd->GetEMCALMatrix(mod)->Print();
 786           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 787         }
 788         fGeomMatrixSet=kTRUE;
 789       }//ESD
 790     }//Load matrices from Data
 791
 792     //Recover time dependent corrections, put then in recalibration histograms. Do it once
 793     fRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
 794
 795   }//first event
 796
 797 }
 798
 799 //____________________________________________________
 800 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
 801 {
 802
 803   // Check if event is exotic, get an exotic cell and compare with triggered patch
 804   // If there is a match remove event ... to be completed, filled with something provisional
 805
 806   if(!fRemoveExoticEvents) return kFALSE;
 807
 808   // Loop on cells
 809   AliVCaloCells * cells = fEvent->GetEMCALCells();
 810   Float_t totCellE = 0;
 811   Int_t bc = InputEvent()->GetBunchCrossNumber();
 812   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
 813
 814     Float_t  ecell = 0 ;
 815     Double_t tcell = 0 ;
 816
 817     Int_t absID   = cells->GetCellNumber(icell);
 818     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
 819     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
 820   }
 821
 822   //  TString triggerclasses = "";
 823   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
 824   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
 825   //    //
 826   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
 827   //    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
 828   //    return;
 829   //
 830
 831   //printf("TotE cell %f\n",totCellE);
 832   if(totCellE < 1) return kTRUE;
 833
 834   return kFALSE;
 835
 836 }
 837
 838 //_________________________________________________
 839 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent()
 840 {
 841   //Check if event is LED
 842
 843   if(!fRemoveLEDEvents) return kFALSE;
 844
 845   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
 846   Int_t ncellsSM3 = 0;
 847   Int_t ncellsSM4 = 0;
 848   AliVCaloCells * cells = fEvent->GetEMCALCells();
 849   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
 850     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
 851     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
 852   }
 853
 854   TString triggerclasses = "";
 855
 856   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
 857   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
 858   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
 859
 860   Int_t ncellcut = 21;
 861   if(triggerclasses.Contains("EMC")) ncellcut = 35;
 862
 863   if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ){
 864     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
 865     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
 866     return kTRUE;
 867   }
 868
 869   return kFALSE;
 870
 871 }
 872
 873 //______________________________________________________________________________
 874 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
 875                                                         TObjArray *recPoints,
 876                                                         TObjArray *clusArray)
 877 {
 878   // Restore clusters from recPoints
 879   // Cluster energy, global position, cells and their amplitude fractions are restored
 880   Int_t j = 0;
 881   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
 882   {
 883     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
 884
 885     const Int_t ncells = recPoint->GetMultiplicity();
 886     Int_t ncellsTrue = 0;
 887
 888     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
 889
 890     // cells and their amplitude fractions
 891     UShort_t   absIds[ncells];
 892     Double32_t ratios[ncells];
 893
 894     //For later check embedding
 895     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 896
 897     Float_t clusterE = 0;
 898     for (Int_t c = 0; c < ncells; c++) {
 899       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
 900
 901       absIds[ncellsTrue] = digit->GetId();
 902       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 903
 904       // In case of unfolding, remove digits with unfolded energy too low
 905       if(fSelectCell){
 906         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  {
 907
 908           if(DebugLevel() > 1)  {
 909             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
 910                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
 911           }
 912
 913           continue;
 914
 915         } // if cuts
 916       }// Select cells
 917
 918       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
 919       ncellsTrue++;
 920       clusterE  +=recPoint->GetEnergiesList()[c];
 921
 922       // In case of embedding, fill ratio with amount of signal,
 923       if (aodIH && aodIH->GetMergeEvents()) {
 924
 925         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
 926         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
 927         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
 928
 929         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 930         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 931         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
 932         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
 933
 934         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
 935         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
 936
 937       }//Embedding
 938
 939     }// cluster cell loop
 940
 941     if (ncellsTrue < 1) {
 942       if (DebugLevel() > 1)
 943         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
 944                recPoint->GetEnergy(), ncells);
 945       continue;
 946     }
 947
 948     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
 949
 950     if(clusterE <  fRecParam->GetClusteringThreshold()) {
 951       if (DebugLevel()>1)
 952         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
 953       continue;
 954     }
 955
 956     TVector3 gpos;
 957     Float_t g[3];
 958
 959     // calculate new cluster position
 960     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
 961     recPoint->GetGlobalPosition(gpos);
 962     gpos.GetXYZ(g);
 963
 964     // create a new cluster
 965     (*clusArray)[j] = new AliAODCaloCluster() ;
 966     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
 967     j++;
 968     clus->SetType(AliVCluster::kEMCALClusterv1);
 969     clus->SetE(clusterE);
 970     clus->SetPosition(g);
 971     clus->SetNCells(ncellsTrue);
 972     clus->SetCellsAbsId(absIds);
 973     clus->SetCellsAmplitudeFraction(ratios);
 974     clus->SetChi2(-1); //not yet implemented
 975     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
 976     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
 977     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
 978
 979     if(ncells == ncellsTrue){
 980       Float_t elipAxis[2];
 981       recPoint->GetElipsAxis(elipAxis);
 982       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 983       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 984       clus->SetDispersion(recPoint->GetDispersion());
 985     }
 986     else if(fSelectCell){
 987       // In case some cells rejected, in unfolding case, recalculate
 988       // shower shape parameters and position
 989       AliVCaloCells* cells = 0x0;
 990       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
 991       else                                  cells = InputEvent()->GetEMCALCells();
 992       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
 993       fRecoUtils->RecalculateClusterPID(clus);
 994       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
 995
 996     }
 997
 998     //MC
 999     Int_t  parentMult = 0;
1000     Int_t *parentList = recPoint->GetParents(parentMult);
1001     clus->SetLabel(parentList, parentMult);
1002
1003     //Write the second major contributor to each MC cluster.
1004     Int_t iNewLabel ;
1005     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
1006
1007       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1008       if(idCell>=0){
1009         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1010         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1011         {
1012           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1013           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1014         }
1015         if (iNewLabel == 1)
1016         {
1017           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1018           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
1019             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1020           }
1021
1022           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1023           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1024           delete [] newLabelArray;
1025         }
1026       }//positive cell number
1027     }
1028
1029   } // recPoints loop
1030
1031 }
1032
1033 //____________________________________________________________
1034 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1035 {
1036   // Init geometry, create list of output clusters
1037
1038   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;
1039   if(fOutputAODBranchName.Length()!=0){
1040     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1041     fOutputAODBranch->SetName(fOutputAODBranchName);
1042     //fOutputAODBranch->SetOwner(kFALSE);
1043     AddAODBranch("TClonesArray", &fOutputAODBranch);
1044   }
1045   else {
1046     AliFatal("fOutputAODBranchName not set\n");
1047   }
1048
1049   //PostData(0,fOutputAODBranch);
1050
1051 }
1052
1053 //_______________________________________________________
1054 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1055 {
1056   // Main loop
1057   // Called for each event
1058
1059   //Remove the contents of output list set in the previous event
1060   fOutputAODBranch->Clear("C");
1061
1062   LoadBranches();
1063
1064   //Get the event, do some checks and settings
1065   CheckAndGetEvent() ;
1066
1067   if (!fEvent) {
1068     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1069     return ;
1070   }
1071
1072   //Init pointers, clusterizer, ocdb
1073   if(fAccessOCDB) AccessOCDB();
1074
1075   InitClusterization();
1076
1077   InitGeometry(); // only once
1078
1079   // Make clusters
1080   if (fJustUnfold) ClusterUnfolding();
1081   else             ClusterizeCells() ;
1082
1083   //Recalculate track-matching for the new clusters
1084   if(fDoTrackMatching) fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
1085
1086   //Put the new clusters in the AOD list
1087
1088   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
1089   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
1090     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1091
1092     //Add matched track
1093     if(fDoTrackMatching){
1094       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1095       if(trackIndex >= 0){
1096         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1097         if(DebugLevel() > 1)
1098           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
1099       }
1100
1101       Float_t dR = 999., dZ = 999.;
1102       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
1103       newCluster->SetTrackDistance(dR,dZ);
1104
1105     }
1106     else {// Assign previously assigned matched track in reco, very very rough
1107       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1108       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1109     }
1110
1111     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1112     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1113     //printf("\n");
1114
1115     //In case of new bad channels, recalculate distance to bad channels
1116     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
1117       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1118
1119     newCluster->SetID(i);
1120     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
1121
1122     if(DebugLevel() > 1 )
1123       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
1124
1125   } // cluster loop
1126
1127   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1128
1129   // Clean up
1130   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
1131
1132   //PostData(0,fOutputAODBranch);
1133
1134 }
1135