HLT/EMCAL/AliHLTEMCALUtils.cxx

   1 #include "AliHLTEMCALUtils.h"
   2
   3 #include "TClass.h"
   4 #include "TFile.h"
   5 #include "TGeoManager.h"
   6 #include "TUnixSystem.h"
   7 #include "TTree.h"
   8 #include "TClonesArray.h"
   9 #include "TBranch.h"
  10 #include "TArrayS.h"
  11
  12 #include "AliCDBManager.h"
  13 #include "AliCDBEntry.h"
  14 #include "AliRawReaderMemory.h"
  15
  16 #include "AliEMCALRecParam.h"
  17 #include "AliEMCALReconstructor.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliEMCALClusterizerv1.h"
  20 #include "AliEMCALRawUtils.h"
  21
  22 #include "AliESDCaloCells.h"
  23 #include "AliESDCaloCluster.h"
  24 #include "AliEMCALDigit.h"
  25 #include "AliESDEvent.h"
  26 #include "AliEMCALRecPoint.h"
  27 #include "AliEMCALPID.h"
  28
  29 #include <fstream>
  30 #include <iostream>
  31 using namespace std;
  32
  33 ClassImp(AliHLTEMCALUtils);
  34
  35 Int_t                  AliHLTEMCALUtils::fgDebug        = 0;
  36 AliEMCALGeometry*      AliHLTEMCALUtils::fgGeom         = NULL;
  37 AliEMCALClusterizerv1* AliHLTEMCALUtils::fgClusterizer  = NULL;
  38 AliEMCALRecParam*      AliHLTEMCALUtils::fgRecParam     = NULL;
  39 AliEMCALRawUtils*      AliHLTEMCALUtils::fgRawUtils     = NULL;
  40 TFile*                 AliHLTEMCALUtils::fgGeometryFile = NULL;
  41 TGeoManager*           AliHLTEMCALUtils::fgGeoManager   = NULL;
  42
  43 TTree*                 AliHLTEMCALUtils::fgClustersTreeTmp = NULL;
  44 TTree*                 AliHLTEMCALUtils::fgDigitsTreeTmp   = NULL;
  45 TClonesArray*          AliHLTEMCALUtils::fgDigitsArrTmp    = NULL;
  46 TBranch*               AliHLTEMCALUtils::fgDigitsBranchTmp = NULL;
  47
  48
  49 //____________________________________________________________________________
  50 AliHLTEMCALUtils::AliHLTEMCALUtils()
  51   : TObject()
  52 {
  53   //
  54   // Default constructor
  55   //
  56   ;
  57 }
  58
  59 //____________________________________________________________________________
  60 AliHLTEMCALUtils::~AliHLTEMCALUtils()
  61 {
  62   //
  63   // Destructor
  64   //
  65   ;
  66 }
  67
  68 //____________________________________________________________________________
  69 void AliHLTEMCALUtils::Cleanup()
  70 {
  71   DeleteStaticMembers();
  72   DeleteReconstructionTrees();
  73 }
  74
  75 //____________________________________________________________________________
  76 void AliHLTEMCALUtils::DeleteStaticMembers()
  77 {
  78   if (fgRecParam != AliEMCALReconstructor::GetRecParam())
  79     {
  80       delete fgRecParam;
  81       fgRecParam = NULL;
  82     }
  83
  84   delete fgRawUtils;
  85   fgRawUtils = NULL;
  86
  87   DeleteGeoManager();
  88
  89   delete fgGeom;
  90   fgGeom = NULL;
  91
  92   delete fgClusterizer;
  93   fgClusterizer = NULL;
  94 }
  95
  96 //____________________________________________________________________________
  97 AliHLTEMCALUtils::AliHLTEMCALUtils(const AliHLTEMCALUtils & /*t*/)
  98   : TObject()
  99 {
 100   //
 101   // copy ctor not to be used
 102   //
 103   AliFatal("May not use.");
 104 }
 105
 106 //____________________________________________________________________________
 107 AliHLTEMCALUtils& AliHLTEMCALUtils::operator = (const AliHLTEMCALUtils & /*t*/)
 108 {
 109   //
 110   // assignement operator not to be used
 111   //
 112   AliFatal("May not use.") ;
 113   return *this ;
 114 }
 115
 116 //____________________________________________________________________________
 117 void AliHLTEMCALUtils::InitRecParam()
 118 {
 119   //
 120   // Please check the AliEMCALReconstructor for comparison
 121   // Check if the instance of AliEMCALRecParam exists,
 122   // if not, get it from OCDB if available, otherwise create a default one
 123   //
 124
 125   fgRecParam = (AliEMCALRecParam*) AliEMCALReconstructor::GetRecParam();
 126   if (fgRecParam)
 127     return;
 128
 129  if (!fgRecParam  && (AliCDBManager::Instance()->IsDefaultStorageSet()))
 130    {
 131      AliCDBEntry *entry = (AliCDBEntry*)
 132        AliCDBManager::Instance()->Get("EMCAL/Config/RecParam");
 133      if (entry) fgRecParam = (AliEMCALRecParam*) entry->GetObject();
 134    }
 135
 136   if(!fgRecParam)
 137     {
 138       AliWarningClass("The Reconstruction parameters initialized to default.");
 139       fgRecParam = new AliEMCALRecParam;
 140     }
 141
 142   if (!fgRecParam)
 143     {
 144       AliErrorClass("Unable to init the reco params. Something is really wrong. Memory?");
 145     }
 146
 147   AliEMCALReconstructor::SetRecParam(fgRecParam);
 148 }
 149
 150 //____________________________________________________________________________
 151 AliEMCALRecParam* AliHLTEMCALUtils::GetRecParam()
 152 {
 153   //
 154   // Init the parameters and reuse the Reconstructor
 155   //
 156
 157   AliHLTEMCALUtils::InitRecParam();
 158   return fgRecParam;
 159 }
 160
 161 //____________________________________________________________________________
 162 AliEMCALRawUtils*      AliHLTEMCALUtils::GetRawUtils(AliEMCALGeometry *pGeometry)
 163 {
 164   //
 165   // Init EMCAL raw utils
 166   // Use the geometry or try to figure out the default
 167   // The OCDB must be initialized.
 168   // Otherwise the default contructor of the raw utils will complain.
 169   //
 170
 171   if (fgRawUtils == NULL)
 172     {
 173       if (AliCDBManager::Instance()->IsDefaultStorageSet())
 174         {
 175           if (pGeometry)
 176             {
 177               fgRawUtils = new AliEMCALRawUtils(pGeometry);
 178             }
 179           else
 180             {
 181               fgRawUtils = new AliEMCALRawUtils(AliHLTEMCALUtils::GetGeometry());
 182             }
 183
 184           AliInfoClass("Raw Utils initialized.");
 185         }
 186       else
 187         {
 188           AliErrorClass("OCDB not initialized. Unable to init raw utils.");
 189         }
 190     }
 191
 192   return fgRawUtils;
 193 }
 194
 195 //____________________________________________________________________________
 196 AliEMCALClusterizerv1* AliHLTEMCALUtils::GetClusterizer()
 197 {
 198   //
 199   // Init EMCAL clusterizer
 200   //
 201
 202   if (fgClusterizer == NULL)
 203     {
 204       AliEMCALGeometry *geom = GetGeometry();
 205       fgClusterizer = new AliEMCALClusterizerv1(geom);
 206       AliInfoClass("ClusterizerV1 initialized.");
 207    }
 208
 209   return fgClusterizer;
 210 }
 211
 212 //____________________________________________________________________________
 213 AliEMCALGeometry*      AliHLTEMCALUtils::GetGeometry()
 214 {
 215   //
 216   // Init EMCAL geometry
 217   //
 218
 219   if (fgGeom == NULL)
 220     {
 221       AliInfoClass(Form("Using default geometry"));
 222       fgGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
 223     }
 224   return fgGeom;
 225 }
 226
 227 //____________________________________________________________________________
 228 void AliHLTEMCALUtils::DeleteGeoManager()
 229 {
 230   //
 231   // Delete the geom manager and the geom TFile
 232   //
 233
 234   if (fgGeometryFile)
 235     {
 236       fgGeometryFile->Close();
 237       delete fgGeometryFile;
 238       fgGeometryFile = NULL;
 239     }
 240
 241   if (fgGeoManager)
 242     {
 243       delete fgGeoManager;
 244       fgGeoManager = NULL;
 245     }
 246 }
 247
 248 //____________________________________________________________________________
 249 Bool_t AliHLTEMCALUtils::LoadGeoManagerFromFile(const char *fname)
 250 {
 251   //
 252   // Open the geom file and get the geom manager
 253   //
 254
 255   Bool_t bReturn = kFALSE;
 256
 257   DeleteGeoManager();
 258
 259   fgGeometryFile = TFile::Open(fname);
 260   if (fgGeometryFile)
 261     {
 262       fgGeoManager = (TGeoManager *)fgGeometryFile->Get("Geometry");
 263     }
 264
 265   if (!fgGeoManager)
 266     {
 267       DeleteGeoManager();
 268       AliErrorClass(Form("Unable to retrieve TGeoManager <Geometry> from %s", fname));
 269     }
 270
 271   bReturn = kTRUE;
 272
 273   return bReturn;
 274 }
 275
 276 //____________________________________________________________________________
 277 Bool_t AliHLTEMCALUtils::InitFakeCDB(const char *cdbpath, Int_t runid)
 278 {
 279   //
 280   // Init fake OCDB - use with care.
 281   // For dummy usage only!
 282   // Make sure there is no conflict with the existing ocdb instance.
 283   //
 284
 285   Bool_t bReturn = kFALSE;
 286
 287   TString sPath = cdbpath;
 288
 289   if (sPath.Length() <= 0)
 290     {
 291       sPath  = "local://";
 292       sPath += gSystem->Getenv("ALICE_ROOT");
 293       AliInfoClass(Form("Setting cdbpath to %s", sPath.Data()));
 294     }
 295
 296   AliCDBManager* pCDB = AliCDBManager::Instance();
 297   if (!pCDB)
 298     {
 299       AliErrorClass("Unable to get the CDBManager instance.");
 300     }
 301   else
 302     {
 303       pCDB->SetRun(runid); // THIS HAS TO BE RETRIEVED !!!
 304       pCDB->SetDefaultStorage(sPath.Data());
 305       AliInfoClass(Form("Fake CDB initialized with path %s", sPath.Data()));
 306       bReturn = kTRUE;
 307     }
 308
 309   return bReturn;
 310 }
 311
 312 //____________________________________________________________________________
 313 Bool_t AliHLTEMCALUtils::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree, Option_t* sOption)
 314 {
 315   // This method is a reusage the corresponding one of AliEMCALReconstructor
 316   // Conversion from raw data to EMCAL digits.
 317   // Works on a single-event basis
 318
 319   // Make sure raw utils are initialized (they depend on geometry and OCDB) !
 320
 321   // We Assume the rawReader, digitsTree and sOption are valid!
 322
 323   // Here we assume the raw reader is setup properly, so we do not need extra reset!
 324   // rawReader->Reset() ;
 325
 326   Bool_t bReturn = kFALSE;
 327   if (!fgRawUtils)
 328     {
 329       // Try to initialize - this is potentially dangerous
 330       // We will use default geometry - most probably
 331       // The OCDB must exist! otherwise we will fail.
 332
 333       if (GetRawUtils() == NULL)
 334         {
 335           // This is no go.
 336           AliErrorClass("Failed. No RawUtils.");
 337           return bReturn;
 338         }
 339     }
 340
 341   if (!fgRecParam)
 342     {
 343       // Again we try to init the RecParams to default
 344       // One should better do it properly before!
 345       // The ocdb has to be there! for the recparams init
 346
 347       if (GetRecParam() == NULL)
 348         {
 349           // This is no go.
 350           AliErrorClass("Failed. No RecParams.");
 351           return bReturn;
 352         }
 353     }
 354
 355   TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",200);
 356   Int_t bufsize = 32000;
 357   digitsTree->Branch("EMCAL", &digitsArr, bufsize);
 358
 359   fgRawUtils->SetOption(sOption);
 360
 361   fgRawUtils->SetRawFormatHighLowGainFactor(fgRecParam->GetHighLowGainFactor());
 362   fgRawUtils->SetRawFormatOrder(fgRecParam->GetOrderParameter());
 363   fgRawUtils->SetRawFormatTau(fgRecParam->GetTau());
 364   fgRawUtils->SetNoiseThreshold(fgRecParam->GetNoiseThreshold());
 365   fgRawUtils->SetNPedSamples(fgRecParam->GetNPedSamples());
 366
 367   fgRawUtils->Raw2Digits(rawReader,digitsArr);
 368
 369   digitsTree->Fill();
 370   digitsArr->Delete();
 371   delete digitsArr;
 372
 373   bReturn = kTRUE;
 374   return bReturn;
 375 }
 376
 377 //____________________________________________________________________________
 378 TTree *AliHLTEMCALUtils::GetDigitsTree()
 379 {
 380   //
 381   // Create the digits tree and the digits arrays
 382   //
 383
 384   if (fgDigitsTreeTmp == NULL)
 385     {
 386       fgDigitsTreeTmp = new TTree("EMCALdigits", "EMCAL digits default");
 387     }
 388
 389   if (fgDigitsArrTmp == NULL)
 390     {
 391       fgDigitsArrTmp = new TClonesArray("AliEMCALDigit",200);
 392     }
 393
 394   if (fgDigitsBranchTmp == NULL)
 395     {
 396       Int_t bufsize = 32000;
 397       fgDigitsBranchTmp = fgDigitsTreeTmp->Branch("EMCAL", &fgDigitsArrTmp, bufsize);
 398     }
 399
 400   fgDigitsBranchTmp->SetAddress(&fgDigitsArrTmp);
 401
 402   return fgDigitsTreeTmp;
 403 }
 404
 405 //____________________________________________________________________________
 406 TTree *AliHLTEMCALUtils::GetClustersTree()
 407 {
 408   //
 409   // Create the clusters tree if necessary
 410   //
 411
 412   if (fgClustersTreeTmp == NULL)
 413     {
 414       fgClustersTreeTmp = new TTree("EMCALclusters", "EMCAL clusters");
 415     }
 416
 417   return fgClustersTreeTmp;
 418 }
 419
 420 //____________________________________________________________________________
 421 void AliHLTEMCALUtils::ResetReconstructionTrees()
 422 {
 423   //
 424   // Reset the content of the reconstruction trees (digits and clusters)
 425   //
 426
 427   TTree *tmp = 0;
 428   tmp = GetDigitsTree();
 429   if (tmp)
 430     {
 431       fgDigitsArrTmp->Delete();
 432       tmp->Reset();
 433     }
 434
 435   tmp = GetClustersTree();
 436   if (tmp)
 437     {
 438       tmp->Reset();
 439     }
 440 }
 441
 442 //____________________________________________________________________________
 443 void AliHLTEMCALUtils::DeleteReconstructionTrees()
 444 {
 445   //
 446   // Delete the reconstruction trees (digits and clusters)
 447   //
 448
 449   TTree *tmp = 0;
 450   tmp = GetDigitsTree();
 451   if (tmp)
 452     {
 453       fgDigitsArrTmp->Delete();
 454       delete fgClustersTreeTmp;
 455       fgClustersTreeTmp = NULL;
 456     }
 457
 458   tmp = GetClustersTree();
 459   if (tmp)
 460     {
 461       delete fgDigitsTreeTmp;
 462       fgDigitsTreeTmp = 0;
 463     }
 464 }
 465
 466 //____________________________________________________________________________
 467 Bool_t AliHLTEMCALUtils::Raw2Clusters(AliRawReader* rawReader, TTree* clustersTree, Option_t* sDigitsOption)
 468 {
 469   //
 470   // Here we go from raw to clusters (not realeasing digits on the way)
 471   // This method combines ConvertDigits and Reconstruct of AliEMCALReconstructor
 472   //
 473
 474   // Here we assume the raw reader is setup properly, so we do not need extra reset!
 475   // rawReader->Reset() ;
 476
 477   Bool_t bReturn = kFALSE;
 478   if (!fgRawUtils)
 479     {
 480       // Try to initialize - this is potentially dangerous
 481       // We will use default geometry - most probably
 482       // The OCDB must exist! otherwise we will fail.
 483
 484       if (GetRawUtils() == NULL)
 485         {
 486           // This is no go.
 487           AliErrorClass("Failed. No RawUtils.");
 488           return bReturn;
 489         }
 490     }
 491
 492   if (!fgRecParam)
 493     {
 494       // Again we try to init the RecParams to default
 495       // One should better do it properly before!
 496       // The ocdb has to be there! for the recparams init
 497
 498       if (GetRecParam() == NULL)
 499         {
 500           // This is no go.
 501           AliErrorClass("Failed. No RecParams.");
 502           return bReturn;
 503         }
 504     }
 505
 506   if (GetDigitsTree() == NULL)
 507     {
 508       // This is no go.
 509       AliErrorClass("Failed. No Digits Tree.");
 510       return bReturn;
 511     }
 512
 513   if (GetClusterizer() == NULL)
 514     {
 515       // This is no go.
 516       AliErrorClass("Failed. No clusterizer.");
 517       return bReturn;
 518     }
 519
 520   // raw->digits
 521   fgRawUtils->SetOption(sDigitsOption);
 522   fgRawUtils->SetRawFormatHighLowGainFactor(fgRecParam->GetHighLowGainFactor());
 523   fgRawUtils->SetRawFormatOrder(fgRecParam->GetOrderParameter());
 524   fgRawUtils->SetRawFormatTau(fgRecParam->GetTau());
 525   fgRawUtils->SetNoiseThreshold(fgRecParam->GetNoiseThreshold());
 526   fgRawUtils->SetNPedSamples(fgRecParam->GetNPedSamples());
 527
 528   fgDigitsArrTmp->Delete();
 529   fgRawUtils->Raw2Digits(rawReader, fgDigitsArrTmp);
 530   fgDigitsTreeTmp->Fill();
 531
 532   // digits->clusters
 533   fgClusterizer->SetOutput(clustersTree);
 534
 535   if (fgDigitsArrTmp->GetEntries() > 0)
 536     {
 537       fgClusterizer->SetInput(fgDigitsTreeTmp);
 538
 539       if (fgDebug > 0)
 540         {
 541           fgClusterizer->Digits2Clusters("deb all") ;
 542         }
 543       else
 544         {
 545           fgClusterizer->Digits2Clusters("");
 546         }
 547
 548       fgClusterizer->Clear();
 549     }
 550
 551   bReturn = kTRUE;
 552   return bReturn;
 553 }
 554
 555 //____________________________________________________________________________
 556 Bool_t AliHLTEMCALUtils::RawBuffer2Clusters(UChar_t *buffer, ULong_t buffSize,
 557                                             Int_t eqID,
 558                                             TTree* clustersTree, Option_t* sDigitsOption)
 559 {
 560   //
 561   // Method provided for convenience
 562   // Take the raw (DDL) mem buffer and put the clusters into the tree
 563   //
 564
 565   Bool_t bReturn = kFALSE;
 566
 567   AliRawReaderMemory rawReader;
 568   rawReader.Reset();
 569   rawReader.SetEquipmentID(eqID);
 570   rawReader.SetMemory(buffer, buffSize);
 571
 572   bReturn = AliHLTEMCALUtils::Raw2Clusters(&rawReader, clustersTree, sDigitsOption);
 573
 574   return bReturn;
 575 }
 576
 577 //____________________________________________________________________________
 578 TTree* AliHLTEMCALUtils::RawBuffer2Clusters(UChar_t *buffer, ULong_t buffSize,
 579                                             Int_t eqID,
 580                                             Option_t* sDigitsOption)
 581 {
 582   //
 583   // Method provided for convenience
 584   // Take the raw (DDL) mem buffer and put the clusters into the tree
 585   //
 586
 587   Bool_t bReturn = kFALSE;
 588
 589   AliRawReaderMemory rawReader;
 590   rawReader.Reset();
 591   rawReader.SetMemory(buffer, buffSize);
 592   rawReader.SetEquipmentID(eqID);
 593
 594   TTree *clustersTree = GetClustersTree();
 595
 596   bReturn = AliHLTEMCALUtils::Raw2Clusters(&rawReader, clustersTree, sDigitsOption);
 597   if (bReturn == kFALSE)
 598     {
 599       clustersTree->Delete();
 600       delete clustersTree;
 601       clustersTree = 0;
 602     }
 603
 604   return clustersTree;
 605 }
 606
 607 //____________________________________________________________________________
 608 UChar_t *AliHLTEMCALUtils::FileToMemory(const char *fname, UChar_t *outpbuffer, ULong_t &buffSize)
 609 {
 610   //
 611   // Method provided for convenience
 612   // Take the raw (DDL) file and copy it to a mem buffer
 613   //
 614
 615   buffSize = 0;
 616   outpbuffer = NULL;
 617
 618   ifstream inputFile(fname, ifstream::in);
 619   if (!inputFile.good())
 620     {
 621       AliErrorClass(Form("Unable to open file", fname));
 622       return outpbuffer;
 623     }
 624
 625   inputFile.seekg (0, ios::end);
 626   buffSize = inputFile.tellg();
 627   inputFile.seekg (0, ios::beg);
 628
 629   void *buff = calloc(buffSize, sizeof(char));
 630
 631   if (!buff)
 632     {
 633       AliErrorClass(Form("Unable to allocate memory: %d bytes", buffSize));
 634       return outpbuffer;
 635     }
 636
 637   inputFile.read((char*)buff, buffSize);
 638   outpbuffer = (UChar_t*)buff;
 639
 640   AliInfoClass(Form("File read: %s. Allocated %d bytes at 0x%x", fname, buffSize, outpbuffer));
 641   return outpbuffer;
 642 }
 643
 644 //____________________________________________________________________________
 645 Bool_t AliHLTEMCALUtils::FillESD(TTree* digitsTree, TTree* clustersTree, AliESDEvent* esd)
 646 {
 647   //
 648   // Fill the ESD just like the AliEMCALReconstructor
 649   //
 650   // This is offline code taken directly from the AliEMCALReconstructor
 651   // Factorization of the function in the AliEMCALReconstructor should allow calling the code directly
 652
 653   // Trigger part on raw data has been removed.
 654
 655   //########################################
 656   //##############Fill CaloCells############
 657   //########################################
 658
 659   if (digitsTree == 0)
 660     {
 661       AliErrorClass("Digits tree is NULL!");
 662       return kFALSE;
 663     }
 664
 665   if (clustersTree == 0)
 666     {
 667       AliErrorClass("Clusters tree is NULL!");
 668       return kFALSE;
 669     }
 670
 671   if (esd == 0)
 672     {
 673       AliErrorClass("Pointer to ESD not set!");
 674       return kFALSE;
 675     }
 676
 677   TClonesArray *digits = new TClonesArray("AliEMCALDigit",1000);
 678   TBranch *branchdig = digitsTree->GetBranch("EMCAL");
 679   if (!branchdig)
 680     {
 681       AliErrorClass("Can't get the branch with the PHOS digits !");
 682       return kFALSE;
 683     }
 684
 685   branchdig->SetAddress(&digits);
 686   digitsTree->GetEvent(0);
 687   Int_t nDigits = digits->GetEntries(), idignew = 0 ;
 688
 689   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
 690   emcCells.CreateContainer(nDigits);
 691   emcCells.SetType(AliESDCaloCells::kEMCALCell);
 692   for (Int_t idig = 0 ; idig < nDigits ; idig++)
 693     {
 694       const AliEMCALDigit * dig = (const AliEMCALDigit*)digits->At(idig);
 695       if(dig->GetAmp() > 0 )
 696         {
 697           emcCells.SetCell(idignew,dig->GetId(),dig->GetAmp(), dig->GetTime());
 698           idignew++;
 699         }
 700     }
 701   emcCells.SetNumberOfCells(idignew);
 702   emcCells.Sort();
 703
 704   //------------------------------------------------------------
 705   //-----------------CLUSTERS-----------------------------
 706   //------------------------------------------------------------
 707
 708   TObjArray *clusters = new TObjArray(100);
 709   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 710   branch->SetAddress(&clusters);
 711   clustersTree->GetEvent(0);
 712
 713   Int_t nClusters = clusters->GetEntries(),  nClustersNew=0;
 714
 715   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters
 716
 717   //######################################################
 718   //#######################TRACK MATCHING###############
 719   //######################################################
 720   //Fill list of integers, each one is index of track to which the cluster belongs.
 721
 722   // step 1 - initialize array of matched track indexes
 723   Int_t *matchedTrack = new Int_t[nClusters];
 724   for (Int_t iclus = 0; iclus < nClusters; iclus++)
 725     matchedTrack[iclus] = -1;  // neg. index --> no matched track
 726
 727   // step 2, change the flag for all matched clusters found in tracks
 728   Int_t iemcalMatch = -1;
 729   Int_t endtpc = esd->GetNumberOfTracks();
 730   for (Int_t itrack = 0; itrack < endtpc; itrack++)
 731     {
 732       AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
 733       iemcalMatch = track->GetEMCALcluster();
 734       if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack;
 735     }
 736
 737   //########################################
 738   //##############Fill CaloClusters#############
 739   //########################################
 740   esd->SetNumberOfEMCALClusters(nClusters);
 741   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++)
 742     {
 743       const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)clusters->At(iClust);
 744       //if(clust->GetClusterType()== AliESDCaloCluster::kEMCALClusterv1) nRP++; else nPC++;
 745
 746       // Get information from EMCAL reconstruction points
 747       Float_t xyz[3];
 748       TVector3 gpos;
 749       clust->GetGlobalPosition(gpos);
 750       for (Int_t ixyz=0; ixyz<3; ixyz++)
 751         xyz[ixyz] = gpos[ixyz];
 752       Float_t elipAxis[2];
 753       clust->GetElipsAxis(elipAxis);
 754       //Create digits lists
 755       Int_t cellMult = clust->GetMultiplicity();
 756       //TArrayS digiList(digitMult);
 757       Float_t *amplFloat = clust->GetEnergiesList();
 758       Int_t   *digitInts = clust->GetAbsId();
 759       TArrayS absIdList(cellMult);
 760       //Uncomment when unfolding is done
 761       //TArrayD fracList(cellMult);
 762
 763       Int_t newCellMult = 0;
 764       for (Int_t iCell=0; iCell<cellMult; iCell++)
 765         {
 766           if (amplFloat[iCell] > 0) {
 767             absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
 768             //Uncomment when unfolding is done
 769             //fracList[newCellMult] = amplFloat[iCell]/emcCells.GetCellAmplitude(digitInts[iCell]);
 770             newCellMult++;
 771           }
 772         }
 773
 774       absIdList.Set(newCellMult);
 775       //Uncomment when unfolding is done
 776       //fracList.Set(newCellMult);
 777
 778       if(newCellMult > 0) { // accept cluster if it has some digit
 779         nClustersNew++;
 780         //Primaries
 781         Int_t  parentMult  = 0;
 782         Int_t *parentList =  clust->GetParents(parentMult);
 783         // fills the ESDCaloCluster
 784         AliESDCaloCluster * ec = new AliESDCaloCluster() ;
 785         ec->SetClusterType(AliESDCaloCluster::kEMCALClusterv1);
 786         ec->SetPosition(xyz);
 787         ec->SetE(clust->GetEnergy());
 788         ec->SetNCells(newCellMult);
 789         //Change type of list from short to ushort
 790         UShort_t *newAbsIdList  = new UShort_t[newCellMult];
 791         //Uncomment when unfolding is done
 792         //Double_t *newFracList  = new Double_t[newCellMult];
 793         for(Int_t i = 0; i < newCellMult ; i++)
 794           {
 795             newAbsIdList[i]=absIdList[i];
 796             //Uncomment when unfolding is done
 797             //newFracList[i]=fracList[i];
 798           }
 799
 800         ec->SetCellsAbsId(newAbsIdList);
 801         //Uncomment when unfolding is done
 802         //ec->SetCellsAmplitudeFraction(newFracList);
 803         ec->SetClusterDisp(clust->GetDispersion());
 804         ec->SetClusterChi2(-1); //not yet implemented
 805         ec->SetM02(elipAxis[0]*elipAxis[0]) ;
 806         ec->SetM20(elipAxis[1]*elipAxis[1]) ;
 807         ec->SetTOF(clust->GetTime()) ; //time-of-fligh
 808         ec->SetNExMax(clust->GetNExMax());          //number of local maxima
 809         TArrayI arrayTrackMatched(1);// Only one track, temporal solution.
 810         arrayTrackMatched[0]= matchedTrack[iClust];
 811         ec->AddTracksMatched(arrayTrackMatched);
 812
 813         TArrayI arrayParents(parentMult,parentList);
 814         ec->AddLabels(arrayParents);
 815
 816         // add the cluster to the esd object
 817         esd->AddCaloCluster(ec);
 818         delete ec;
 819         delete [] newAbsIdList ;
 820         //delete [] newFracList ;
 821       }
 822     } // cycle on clusters
 823
 824   delete [] matchedTrack;
 825
 826   esd->SetNumberOfEMCALClusters(nClustersNew);
 827   //if(nClustersNew != nClusters)
 828   //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
 829
 830   //Fill ESDCaloCluster with PID weights
 831   AliEMCALPID *pid = new AliEMCALPID;
 832   //pid->SetPrintInfo(kTRUE);
 833   pid->SetReconstructor(kTRUE);
 834   pid->RunPID(esd);
 835   delete pid;
 836
 837   delete digits;
 838   delete clusters;
 839
 840   // printf(" ## AliEMCALReconstructor::FillESD() is ended : ncl %i -> %i ### \n ",nClusters, nClustersNew);
 841   return kTRUE;
 842 }