1e0390ee3efba3579b1bcffc3b9f4447a87ba88b
[u/mrichter/AliRoot.git] / 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 }