PWGGA/EMCALTasks/AliAnalysisTaskEMCALClusterize.cxx

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