PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // This analysis provides a new list of clusters to be used in other analysis
  19 //
  20 // Author: Gustavo Conesa Balbastre,
  21 //         Adapted from analysis class from Deepa Thomas
  22 //
  23 //
  24 //_________________________________________________________________________
  25
  26 // --- Root ---
  27 #include "TString.h"
  28 #include "TRefArray.h"
  29 #include "TClonesArray.h"
  30 #include "TTree.h"
  31 #include "TGeoManager.h"
  32 #include "TROOT.h"
  33 #include "TInterpreter.h"
  34
  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 "AliEMCALRecPoint.h"
  56 #include "AliEMCALDigit.h"
  57 #include "AliCaloCalibPedestal.h"
  58 #include "AliEMCALCalibData.h"
  59 #include "AliEMCALRecoUtils.h"
  60
  61 #include "AliAnalysisTaskEMCALClusterize.h"
  62
  63 ClassImp(AliAnalysisTaskEMCALClusterize)
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  67   : AliAnalysisTaskSE(name)
  68   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  69   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  70   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  71   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  72   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kTRUE)
  73   , fRun(-1),            fRecoUtils(0),        fConfigName(""), fCellLabels()
  74
  75   {
  76   //ctor
  77   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
  78   for(Int_t j = 0; j < 12672; j++)  fCellLabels[j] = -1;
  79   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  80   fClusterArr      = new TObjArray(100);
  81   fCaloClusterArr  = new TObjArray(100);
  82   fRecParam        = new AliEMCALRecParam;
  83   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  84   fRecoUtils       = new AliEMCALRecoUtils();
  85 }
  86
  87 //________________________________________________________________________
  88 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
  89   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
  90   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  91   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  92   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  93   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  94   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kFALSE)
  95   , fRun(-1),            fRecoUtils(0),        fConfigName(""), fCellLabels()
  96 {
  97   // Constructor
  98   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
  99   for(Int_t j = 0; j < 24*48*11; j++)  fCellLabels[j] = -1;
 100   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
 101   fClusterArr      = new TObjArray(100);
 102   fCaloClusterArr  = new TObjArray(100);
 103   fRecParam        = new AliEMCALRecParam;
 104   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
 105   fRecoUtils       = new AliEMCALRecoUtils();
 106 }
 107
 108
 109 //________________________________________________________________________
 110 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 111 {
 112   //dtor
 113
 114   if (fDigitsArr){
 115     fDigitsArr->Clear("C");
 116     delete fDigitsArr;
 117   }
 118
 119   if (fClusterArr){
 120     fClusterArr->Delete();
 121     delete fClusterArr;
 122   }
 123
 124   if (fCaloClusterArr){
 125     fCaloClusterArr->Delete();
 126     delete fCaloClusterArr;
 127   }
 128
 129   if(fClusterizer) delete fClusterizer;
 130   if(fUnfolder)    delete fUnfolder;
 131   if(fRecoUtils)   delete fRecoUtils;
 132
 133 }
 134
 135 //-------------------------------------------------------------------
 136 void AliAnalysisTaskEMCALClusterize::Init()
 137 {
 138   //Init analysis with configuration macro if available
 139
 140   if(gROOT->LoadMacro(fConfigName) >=0){
 141     printf("Configure analysis with %s\n",fConfigName.Data());
 142     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 143     fGeomName         = clus->fGeomName;
 144     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 145     fOCDBpath         = clus->fOCDBpath;
 146     fRecParam         = clus->fRecParam;
 147     fJustUnfold       = clus->fJustUnfold;
 148     fFillAODFile      = clus->fFillAODFile;
 149     fRecoUtils        = clus->fRecoUtils;
 150     fConfigName       = clus->fConfigName;
 151     fOutputAODBranchName = clus->fOutputAODBranchName;
 152     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
 153
 154   }
 155
 156 }
 157
 158 //-------------------------------------------------------------------
 159 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
 160 {
 161   // Init geometry, create list of output clusters
 162
 163   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;
 164
 165   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 166   fOutputAODBranch->SetName(fOutputAODBranchName);
 167   AddAODBranch("TClonesArray", &fOutputAODBranch);
 168
 169 }
 170
 171 //________________________________________________________________________
 172 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
 173 {
 174   // Main loop
 175   // Called for each event
 176
 177   //printf("--- Event %d -- \n",(Int_t)Entry());
 178   //Remove the contents of output list set in the previous event
 179   fOutputAODBranch->Clear("C");
 180
 181   //Get the event
 182   AliVEvent   * event    = 0;
 183   AliESDEvent * esdevent = 0;
 184
 185   //Check if input event are embedded events
 186   //If so, take output event
 187   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 188   if (aodIH && aodIH->GetMergeEvents()) {
 189     //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
 190     event  = AODEvent();
 191 //    printf("InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
 192 //           InputEvent()->GetEMCALCells()->GetNumberOfCells());
 193 //    printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
 194 //           AODEvent()->GetEMCALCells()->GetNumberOfCells());
 195   }
 196   else {
 197     event =  InputEvent();
 198     esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 199   }
 200
 201   if (!event) {
 202     Error("UserExec","Event not available");
 203     return;
 204   }
 205
 206   //Magic line to write events to AOD file
 207   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
 208   LoadBranches();
 209
 210   AccessOCDB();
 211
 212   //-------------------------------------------------------------------------------------
 213   //Set the geometry matrix, for the first event, skip the rest
 214   //-------------------------------------------------------------------------------------
 215   if(!fGeomMatrixSet){
 216     if(fLoadGeomMatrices){
 217       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 218         if(fGeomMatrix[mod]){
 219           if(DebugLevel() > 1)
 220             fGeomMatrix[mod]->Print();
 221           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
 222         }
 223         fGeomMatrixSet=kTRUE;
 224       }//SM loop
 225     }//Load matrices
 226     else if(!gGeoManager){
 227       Info("UserExec","Get geo matrices from data");
 228       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 229       if(!strcmp(event->GetName(),"AliAODEvent")) {
 230         if(DebugLevel() > 1)
 231           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 232       }//AOD
 233       else {
 234         if(DebugLevel() > 1)
 235           Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 236         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 237         if(!esd) {
 238           Error("UserExec","This event does not contain ESDs?");
 239           return;
 240         }
 241         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 242           if(DebugLevel() > 1)
 243             esd->GetEMCALMatrix(mod)->Print();
 244           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 245         }
 246         fGeomMatrixSet=kTRUE;
 247       }//ESD
 248     }//Load matrices from Data
 249
 250     //Recover time dependent corrections, put then in recalibration histograms. Do it once
 251     fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
 252
 253   }//first event
 254
 255   //Get the list of cells needed for unfolding and reclustering
 256   AliVCaloCells *cells= event->GetEMCALCells();
 257   Int_t nClustersOrg = 0;
 258   //-------------------------------------------
 259   //---------Unfolding clusters----------------
 260   //-------------------------------------------
 261   if (fJustUnfold) {
 262
 263     //Fill the array with the EMCAL clusters, copy them
 264     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
 265     {
 266       AliVCluster *clus = event->GetCaloCluster(i);
 267       if(clus->IsEMCAL()){
 268         //printf("Org Cluster %d, E %f\n",i, clus->E());
 269
 270         //recalibrate/remove bad channels/etc if requested
 271         if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
 272           //printf("Remove cluster\n");
 273           continue;
 274         }
 275
 276         if(fRecoUtils->IsRecalibrationOn()){
 277           //printf("Energy before %f ",clus->E());
 278           fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 279           //printf("; after %f\n",clus->E());
 280         }
 281         //Cast to ESD or AOD, needed to create the cluster array
 282         AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 283         AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 284         if     (esdCluster){
 285           fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 286         }//ESD
 287         else if(aodCluster){
 288           fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 289         }//AOD
 290         else
 291           Warning("UserExec()"," - Wrong CaloCluster type?");
 292         nClustersOrg++;
 293       }
 294     }
 295
 296     //Do the unfolding
 297     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 298
 299     //CLEAN-UP
 300     fUnfolder->Clear();
 301
 302
 303     //printf("nClustersOrg %d\n",nClustersOrg);
 304   }
 305   //-------------------------------------------
 306   //---------- Recluster cells ----------------
 307   //-------------------------------------------
 308
 309   else{
 310     //-------------------------------------------------------------------------------------
 311     //Transform CaloCells into Digits
 312     //-------------------------------------------------------------------------------------
 313
 314     //In case of MC, first loop on the clusters and fill MC label to array
 315     //.....................................................................
 316
 317 //    for(Int_t j = 0; j < 24*48*11; j++)  {
 318 //      if(fCellLabels[j]!=-1) printf("Not well initialized cell %d, label %d\n",j,fCellLabels[j]) ;
 319 //    }
 320
 321     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
 322     {
 323       AliVCluster *clus = event->GetCaloCluster(i);
 324       if(clus->IsEMCAL()){
 325
 326         Int_t label = clus->GetLabel();
 327         UShort_t * index    = clus->GetCellsAbsId() ;
 328         for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
 329           fCellLabels[index[icell]]=label;
 330           //printf("1) absID %d, label %d\n",index[icell], fCellLabels[index[icell]]);
 331         }
 332         nClustersOrg++;
 333       }
 334     }
 335
 336     // Create digits
 337     //.................
 338     Int_t idigit =  0;
 339     Int_t id     = -1;
 340     Float_t amp  = -1;
 341     Float_t time = -1;
 342
 343     Double_t cellAmplitude = 0;
 344     Double_t cellTime      = 0;
 345     Short_t  cellNumber    = 0;
 346
 347     TTree *digitsTree = new TTree("digitstree","digitstree");
 348     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 349
 350
 351     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 352     {
 353       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 354         break;
 355
 356       time = cellTime;
 357       amp  = cellAmplitude;
 358       id   = cellNumber;
 359
 360       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 361       fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
 362       fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 363
 364       //Do not include bad channels found in analysis?
 365       if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 366           fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 367           fCellLabels[id]=-1; //reset the entry in the array for next event
 368         //printf("Remove channel %d\n",id);
 369         continue;
 370       }
 371
 372       //Recalibrate?
 373       if(fRecoUtils->IsRecalibrationOn()){
 374         //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
 375         amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 376       }
 377
 378       //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 379       new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 380       //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
 381       //else                  printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
 382       fCellLabels[id]=-1; //reset the entry in the array for next event
 383
 384       //AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
 385       //digit->SetId(id);
 386       //digit->SetAmplitude(amp);
 387       //digit->SetTime(time);
 388       //digit->SetTimeR(time);
 389       //digit->SetIndexInList(idigit);
 390       //digit->SetType(AliEMCALDigit::kHG);
 391
 392       //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
 393       idigit++;
 394     }
 395
 396     //Fill the tree with digits
 397     digitsTree->Fill();
 398
 399     //-------------------------------------------------------------------------------------
 400     //Do the clusterization
 401     //-------------------------------------------------------------------------------------
 402     TTree *clustersTree = new TTree("clustertree","clustertree");
 403
 404     fClusterizer->SetInput(digitsTree);
 405     fClusterizer->SetOutput(clustersTree);
 406     fClusterizer->Digits2Clusters("");
 407
 408     //-------------------------------------------------------------------------------------
 409     //Transform the recpoints into AliVClusters
 410     //-------------------------------------------------------------------------------------
 411
 412     clustersTree->SetBranchStatus("*",0); //disable all branches
 413     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 414
 415     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 416     branch->SetAddress(&fClusterArr);
 417     branch->GetEntry(0);
 418
 419     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 420
 421     //---CLEAN UP-----
 422     fClusterizer->Clear();
 423     fDigitsArr  ->Clear("C");
 424     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 425
 426     clustersTree->Delete("all");
 427     digitsTree  ->Delete("all");
 428   }
 429
 430   //Recalculate track-matching for the new clusters, only with ESDs
 431   if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
 432
 433
 434   //-------------------------------------------------------------------------------------
 435   //Put the new clusters in the AOD list
 436   //-------------------------------------------------------------------------------------
 437
 438   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
 439   //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
 440   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
 441     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 442     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
 443
 444     //Add matched track, if any, only with ESDs
 445     if(esdevent){
 446       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 447       if(trackIndex >= 0){
 448         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
 449         if(DebugLevel() > 1)
 450           Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
 451       }
 452     }
 453
 454     //In case of new bad channels, recalculate distance to bad channels
 455     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
 456       //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
 457       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
 458       //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
 459     }
 460
 461 //    if(newCluster->GetNLabels()>0){
 462 //      printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
 463 //      UShort_t * newindex    = newCluster->GetCellsAbsId() ;
 464 //      for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
 465 //       printf("\t absID %d\n",newindex[icell]);
 466 //     }
 467 //    }
 468 //
 469 //    printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
 470
 471     newCluster->SetID(i);
 472     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 473   }
 474
 475   //---CLEAN UP-----
 476   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
 477 }
 478
 479 //_____________________________________________________________________
 480 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 481 {
 482   //Access to OCDB stuff
 483
 484   AliVEvent * event = InputEvent();
 485   if (!event)
 486   {
 487     Warning("AccessOCDB","Event not available!!!");
 488     return kFALSE;
 489   }
 490
 491   if (event->GetRunNumber()==fRun)
 492     return kTRUE;
 493   fRun = event->GetRunNumber();
 494
 495   if(DebugLevel() > 1 )
 496     Info("AccessODCD()"," Begin");
 497
 498   fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 499
 500
 501   AliCDBManager *cdb = AliCDBManager::Instance();
 502
 503
 504   if (fOCDBpath.Length()){
 505     cdb->SetDefaultStorage(fOCDBpath.Data());
 506     Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
 507   }
 508
 509   cdb->SetRun(event->GetRunNumber());
 510
 511   //
 512   // EMCAL from RAW OCDB
 513   if (fOCDBpath.Contains("alien:"))
 514   {
 515     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 516     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 517   }
 518
 519   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 520
 521   // init parameters:
 522
 523   //Get calibration parameters
 524   if(!fCalibData)
 525   {
 526     AliCDBEntry *entry = (AliCDBEntry*)
 527       AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 528
 529     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 530   }
 531
 532   if(!fCalibData)
 533     AliFatal("Calibration parameters not found in CDB!");
 534
 535   //Get calibration parameters
 536   if(!fPedestalData)
 537   {
 538     AliCDBEntry *entry = (AliCDBEntry*)
 539       AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 540
 541     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 542   }
 543
 544   if(!fPedestalData)
 545     AliFatal("Dead map not found in CDB!");
 546
 547   InitClusterization();
 548
 549   return kTRUE;
 550 }
 551
 552 //________________________________________________________________________________________
 553 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 554 {
 555   //Select clusterization/unfolding algorithm and set all the needed parameters
 556
 557   if (fJustUnfold){
 558     // init the unfolding afterburner
 559     delete fUnfolder;
 560     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 561     return;
 562  }
 563
 564   //First init the clusterizer
 565   delete fClusterizer;
 566   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 567     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 568   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 569     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 570   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
 571    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 572    clusterizer->SetNRowDiff(2);
 573    clusterizer->SetNColDiff(2);
 574    fClusterizer = clusterizer;
 575   } else{
 576     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 577   }
 578
 579   //Now set the parameters
 580   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 581   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 582   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 583   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 584   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 585   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 586   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 587   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 588   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 589   fClusterizer->SetJustClusters          ( kTRUE                               );
 590   //In case of unfolding after clusterization is requested, set the corresponding parameters
 591   if(fRecParam->GetUnfold()){
 592     Int_t i=0;
 593     for (i = 0; i < 8; i++) {
 594       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 595     }//end of loop over parameters
 596     for (i = 0; i < 3; i++) {
 597       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 598       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 599     }//end of loop over parameters
 600
 601     fClusterizer->InitClusterUnfolding();
 602   }// to unfold
 603 }
 604
 605 //________________________________________________________________________________________
 606 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
 607 {
 608   // Restore clusters from recPoints
 609   // Cluster energy, global position, cells and their amplitude fractions are restored
 610
 611   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
 612   {
 613     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
 614
 615     const Int_t ncells = recPoint->GetMultiplicity();
 616     Int_t ncells_true = 0;
 617
 618     // cells and their amplitude fractions
 619     UShort_t   absIds[ncells];
 620     Double32_t ratios[ncells];
 621
 622     for (Int_t c = 0; c < ncells; c++) {
 623       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
 624
 625       absIds[ncells_true] = digit->GetId();
 626       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 627
 628       if (ratios[ncells_true] > 0.001) ncells_true++;
 629     }
 630
 631     if (ncells_true < 1) {
 632       Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
 633       continue;
 634     }
 635
 636     TVector3 gpos;
 637     Float_t g[3];
 638
 639     // calculate new cluster position
 640     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
 641     recPoint->GetGlobalPosition(gpos);
 642     gpos.GetXYZ(g);
 643
 644     // create a new cluster
 645     AliAODCaloCluster *clus = new AliAODCaloCluster();
 646     clus->SetType(AliVCluster::kEMCALClusterv1);
 647     clus->SetE(recPoint->GetEnergy());
 648     clus->SetPosition(g);
 649     clus->SetNCells(ncells_true);
 650     clus->SetCellsAbsId(absIds);
 651     clus->SetCellsAmplitudeFraction(ratios);
 652     clus->SetDispersion(recPoint->GetDispersion());
 653     clus->SetChi2(-1); //not yet implemented
 654     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
 655     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
 656     Float_t elipAxis[2];
 657     recPoint->GetElipsAxis(elipAxis);
 658     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 659     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 660     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
 661
 662     //MC
 663     Int_t  parentMult  = 0;
 664     Int_t *parentList = recPoint->GetParents(parentMult);
 665     clus->SetLabel(parentList, parentMult);
 666
 667     clusArray->Add(clus);
 668   } // recPoints loop
 669 }