PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.cxx

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