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