PWGGA/EMCALTasks/AliAnalysisTaskEMCALClusterize.cxx

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