c0b8207f94f50ce542e3d3c3d3bdbe92cc42f772
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPoint.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 //  Reconstructed Points for the EMCAL
18 //  A RecPoint is a cluster of digits
19 //*-- Author: Yves Schutz (SUBATECH)
20 //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
21 //*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
22
23 // --- ROOT system ---
24 class Riostream;
25 #include <TPad.h>
26 class TGraph;
27 class TPaveText;
28 #include <TClonesArray.h>
29 #include <TMath.h>
30
31 // --- Standard library ---
32
33 // --- AliRoot header files ---
34 //#include "AliGenerator.h"
35 class AliGenerator;
36 #include "AliRunLoader.h"
37 #include "AliRun.h"
38 class AliEMCAL;
39 #include "AliEMCALLoader.h"
40 #include "AliEMCALGeometry.h"
41 #include "AliEMCALHit.h"
42 #include "AliEMCALDigit.h"
43 #include "AliEMCALRecPoint.h"
44
45 ClassImp(AliEMCALRecPoint)
46
47 //____________________________________________________________________________
48 AliEMCALRecPoint::AliEMCALRecPoint()
49   : AliRecPoint(),
50     fGeomPtr(0),
51     fClusterType(-1),
52     fCoreEnergy(0),
53     fDispersion(0),
54     fEnergyList(0),
55     fTimeList(0),
56     fAbsIdList(0),
57     fTime(0.),
58     fCoreRadius(10),  //HG check this
59     fDETracksList(0),
60     fMulParent(0),
61     fMaxParent(0),
62     fParentsList(0),
63     fDEParentsList(0),
64     fSuperModuleNumber(0),
65     fDigitIndMax(-1)
66 {
67   // ctor
68   AliRunLoader *rl = AliRunLoader::GetRunLoader();
69   if (rl && rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
70     fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
71   else
72     fGeomPtr = AliEMCALGeometry::GetInstance();
73   //fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
74     //  fGeomPtr->GetTransformationForSM(); // Global <-> Local
75 }
76
77 //____________________________________________________________________________
78 AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) 
79   : AliRecPoint(opt),
80     fGeomPtr(0),
81     fClusterType(-1),
82     fCoreEnergy(0),
83     fDispersion(0),
84     fEnergyList(0),
85     fTimeList(0),
86     fAbsIdList(0),
87     fTime(-1.),
88     fCoreRadius(10),  //HG check this
89     fDETracksList(0),
90     fMulParent(0),
91     fMaxParent(1000),
92     fParentsList(0),
93     fDEParentsList(0),
94     fSuperModuleNumber(0),
95     fDigitIndMax(-1)
96 {
97   // ctor
98   // Increase fMaxTrack for EMCAL.
99   delete [] fTracksList;
100   fMaxTrack = 1000 ;
101   fTracksList = new Int_t[fMaxTrack];
102   fDETracksList = new Float_t[fMaxTrack];
103
104   fParentsList = new Int_t[fMaxParent];
105   fDEParentsList = new Float_t[fMaxParent];
106   for (Int_t i = 0; i < fMaxTrack; i++)
107     fDETracksList[i] = 0;
108   for (Int_t i = 0; i < fMaxParent; i++) {
109     fParentsList[i] = -1;
110     fDEParentsList[i] = 0;
111   }
112
113   AliRunLoader *rl = AliRunLoader::GetRunLoader();
114   if (rl && rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
115     fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
116   else
117     fGeomPtr = AliEMCALGeometry::GetInstance();
118   //    fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
119   //  fGeomPtr->GetTransformationForSM(); // Global <-> Local
120 }
121
122 //____________________________________________________________________________
123 AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp) 
124   : AliRecPoint(rp),
125     fGeomPtr(rp.fGeomPtr),
126     fClusterType(rp.fClusterType),
127     fCoreEnergy(rp.fCoreEnergy),
128     fDispersion(rp.fDispersion),
129     fEnergyList(0),
130     fTimeList(0),
131     fAbsIdList(0),
132     fTime(rp.fTime),
133     fCoreRadius(rp.fCoreRadius),
134     fDETracksList(0),
135     fMulParent(rp.fMulParent),
136     fMaxParent(rp.fMaxParent),
137     fParentsList(0),
138     fDEParentsList(0),
139     fSuperModuleNumber(rp.fSuperModuleNumber),
140     fDigitIndMax(rp.fDigitIndMax)
141 {
142   //copy ctor
143   fLambda[0] = rp.fLambda[0];
144   fLambda[1] = rp.fLambda[1];
145
146   fEnergyList = new Float_t[rp.fMaxDigit];
147   fTimeList = new Float_t[rp.fMaxDigit];
148   fAbsIdList = new Int_t[rp.fMaxDigit];
149   for(Int_t i = 0; i < rp.fMulDigit; i++) {
150     fEnergyList[i] = rp.fEnergyList[i];
151     fTimeList[i] = rp.fTimeList[i];
152     fAbsIdList[i] = rp.fAbsIdList[i];
153   }
154   fDETracksList = new Float_t[rp.fMaxTrack];
155   for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i];
156   fParentsList = new Int_t[rp.fMaxParent];
157   for(Int_t i = 0; i < rp.fMulParent; i++) fParentsList[i] = rp.fParentsList[i];
158   fDEParentsList = new Float_t[rp.fMaxParent];
159   for(Int_t i = 0; i < rp.fMulParent; i++) fDEParentsList[i] = rp.fDEParentsList[i];
160
161 }
162 //____________________________________________________________________________
163 AliEMCALRecPoint::~AliEMCALRecPoint()
164 {
165   // dtor
166   if ( fEnergyList )
167     delete[] fEnergyList ; 
168   if ( fTimeList )
169     delete[] fTimeList ; 
170   if ( fAbsIdList )
171     delete[] fAbsIdList ; 
172    if ( fDETracksList)
173     delete[] fDETracksList;
174    if ( fParentsList)
175     delete[] fParentsList;
176    if ( fDEParentsList)
177     delete[] fDEParentsList;
178 }
179
180 //____________________________________________________________________________
181 void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
182 {
183   // Adds a digit to the RecPoint
184   // and accumulates the total amplitude and the multiplicity 
185   
186   if(fEnergyList == 0)
187     fEnergyList =  new Float_t[fMaxDigit]; 
188   if(fTimeList == 0)
189     fTimeList =  new Float_t[fMaxDigit]; 
190   if(fAbsIdList == 0) {
191     fAbsIdList =  new Int_t[fMaxDigit];
192     fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
193   }
194
195   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
196     fMaxDigit*=2 ; 
197     Int_t * tempo = new Int_t[fMaxDigit]; 
198     Float_t * tempoE =  new Float_t[fMaxDigit];
199     Float_t * tempoT =  new Float_t[fMaxDigit];
200     Int_t * tempoId = new Int_t[fMaxDigit]; 
201
202     Int_t index ;     
203     for ( index = 0 ; index < fMulDigit ; index++ ){
204       tempo[index]   = fDigitsList[index] ;
205       tempoE[index]  = fEnergyList[index] ; 
206       tempoT[index]  = fTimeList[index] ; 
207       tempoId[index] = fAbsIdList[index] ; 
208     }
209     
210     delete [] fDigitsList ;
211     delete [] fEnergyList ;
212     delete [] fTimeList ;
213     delete [] fAbsIdList ;
214
215     fDigitsList = tempo;
216     fEnergyList = tempoE; 
217     fTimeList = tempoT;
218     fAbsIdList = tempoId;
219   } // if
220   
221   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ; 
222   fEnergyList[fMulDigit]   = Energy ;
223   fTimeList[fMulDigit]     = digit.GetTimeR() ;
224   fAbsIdList[fMulDigit]    = digit.GetId();
225   fMulDigit++ ; 
226   fAmp += Energy ; 
227
228 }
229 //____________________________________________________________________________
230 Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
231 {
232   // Tells if (true) or not (false) two digits are neighbours
233   // A neighbour is defined as being two digits which share a corner
234   
235   static Bool_t areNeighbours = kFALSE ;
236   static Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
237   static int nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
238   static Int_t relid1[2] , relid2[2] ; // ieta, iphi
239   static Int_t rowdiff=0, coldiff=0;
240
241   areNeighbours = kFALSE ;
242
243   fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nModule,nIphi,nIeta);
244   fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
245
246   fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
247   fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
248   
249   rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
250   coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
251
252   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
253   areNeighbours = kTRUE ;
254   
255   return areNeighbours;
256 }
257
258 //____________________________________________________________________________
259 Int_t AliEMCALRecPoint::Compare(const TObject * obj) const
260 {
261   // Compares two RecPoints according to their position in the EMCAL modules
262
263   Float_t delta = 1 ; //Width of "Sorting row". If you change this 
264                       //value (what is senseless) change as well delta in
265                       //AliEMCALTrackSegmentMakerv* and other RecPoints...
266   Int_t rv ; 
267
268   AliEMCALRecPoint * clu = (AliEMCALRecPoint *)obj ; 
269
270   TVector3 locpos1; 
271   GetLocalPosition(locpos1);
272   TVector3 locpos2;  
273   clu->GetLocalPosition(locpos2);  
274
275   Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ;
276   if (rowdif> 0) 
277     rv = 1 ;
278   else if(rowdif < 0) 
279     rv = -1 ;
280   else if(locpos1.Y()>locpos2.Y()) 
281     rv = -1 ;
282   else 
283     rv = 1 ; 
284
285   return rv ; 
286 }
287
288 //____________________________________________________________________________
289 Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py)
290 {
291   // Compute distance from point px,py to  a AliEMCALRecPoint considered as a Tmarker
292   // Compute the closest distance of approach from point px,py to this marker.
293   // The distance is computed in pixels units.
294   // HG Still need to update -> Not sure what this should achieve
295
296   TVector3 pos(0.,0.,0.) ;
297   GetLocalPosition(pos) ;
298   Float_t x =  pos.X() ;
299   Float_t y =  pos.Y() ;
300   const Int_t kMaxDiff = 10;
301   Int_t pxm  = gPad->XtoAbsPixel(x);
302   Int_t pym  = gPad->YtoAbsPixel(y);
303   Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym);
304   
305   if (dist > kMaxDiff) return 9999;
306   return dist;
307 }
308
309 //___________________________________________________________________________
310  void AliEMCALRecPoint::Draw(Option_t *option)
311  {
312    // Draw this AliEMCALRecPoint with its current attributes
313    
314    AppendPad(option);
315  }
316
317 //______________________________________________________________________________
318 void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t)
319 {
320   // Execute action corresponding to one event
321   // This member function is called when a AliEMCALRecPoint is clicked with the locator
322   //
323   // If Left button is clicked on AliEMCALRecPoint, the digits are switched on    
324   // and switched off when the mouse button is released.
325
326   //  static Int_t pxold, pyold;
327
328   /* static TGraph *  digitgraph = 0 ;
329   static TPaveText* clustertext = 0 ;
330   
331   if (!gPad->IsEditable()) return;
332   
333   switch (event) {
334     
335     
336   case kButton1Down:{
337     AliEMCALDigit * digit ;
338     AliEMCALGeometry * emcalgeom =  (AliEMCALGetter::Instance())->EMCALGeometry() ;
339
340     Int_t iDigit;
341     Int_t relid[2] ;
342   
343     const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ;
344     Float_t * xi = new Float_t [kMulDigit] ; 
345     Float_t * zi = new Float_t [kMulDigit] ;
346     
347     for(iDigit = 0; iDigit < kMulDigit; iDigit++) {
348       Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code"); 
349       digit = 0 ; //dynamic_cast<AliEMCALDigit *>((fDigitsList)[iDigit]);
350       emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
351       emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ;
352     }
353     
354     if (!digitgraph) {
355       digitgraph = new TGraph(fMulDigit,xi,zi);
356       digitgraph-> SetMarkerStyle(5) ; 
357       digitgraph-> SetMarkerSize(1.) ;
358       digitgraph-> SetMarkerColor(1) ;
359       digitgraph-> Draw("P") ;
360     }
361     if (!clustertext) {
362       
363       TVector3 pos(0.,0.,0.) ;
364       GetLocalPosition(pos) ;
365       clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ;
366       Text_t  line1[40] ;
367       Text_t  line2[40] ;
368       sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ;
369       sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ;
370       clustertext ->AddText(line1) ;
371       clustertext ->AddText(line2) ;
372       clustertext ->Draw("");
373     }
374     gPad->Update() ; 
375     Print("") ;
376     delete[] xi ; 
377     delete[] zi ; 
378    }
379   
380 break;
381   
382   case kButton1Up:
383     if (digitgraph) {
384       delete digitgraph  ;
385       digitgraph = 0 ;
386     }
387     if (clustertext) {
388       delete clustertext ;
389       clustertext = 0 ;
390     }
391     
392     break;
393     
394     }*/
395 }
396 //____________________________________________________________________________
397 void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) 
398 {
399   // Evaluates all shower parameters
400   EvalLocalPosition(logWeight, digits) ;
401   EvalElipsAxis(logWeight, digits) ;
402   EvalDispersion(logWeight, digits) ;
403   //EvalCoreEnergy(logWeight, digits);
404   EvalTime(digits) ;
405   EvalPrimaries(digits) ;
406   EvalParents(digits);
407 }
408
409 //____________________________________________________________________________
410 void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
411 {
412   // Calculates the dispersion of the shower at the origin of the RecPoint
413   // in cell units - Nov 16,2006
414
415   Double_t d = 0., wtot = 0., w = 0.;
416   Int_t iDigit=0, nstat=0;
417   AliEMCALDigit * digit ;
418   
419   // Calculates the dispersion in cell units 
420   Double_t etai, phii, etaMean=0.0, phiMean=0.0; 
421   int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
422   int iphi=0, ieta=0;
423   // Calculate mean values
424   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
425     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
426
427     if (fAmp>0 && fEnergyList[iDigit]>0) {
428       fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
429       fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
430       etai=(Double_t)ieta;
431       phii=(Double_t)iphi;
432       w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
433
434       if(w>0.0) {
435         phiMean += phii*w;
436         etaMean += etai*w;
437         wtot    += w;
438       }
439     }
440   }
441   if (wtot>0) {
442     phiMean /= wtot ;
443     etaMean /= wtot ;
444   } else AliError(Form("Wrong weight %f\n", wtot));
445
446   // Calculate dispersion
447   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
448     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
449
450     if (fAmp>0 && fEnergyList[iDigit]>0) {
451       fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
452       fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
453       etai=(Double_t)ieta;
454       phii=(Double_t)iphi;
455       w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
456
457       if(w>0.0) {
458         nstat++;
459         d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
460       }
461     }
462   }
463   
464   if ( wtot > 0 && nstat>1) d /= wtot ;
465   else                      d = 0. ; 
466
467   fDispersion = TMath::Sqrt(d) ;
468 }
469
470 //____________________________________________________________________________
471 void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
472 {
473   // Calculates the center of gravity in the local EMCAL-module coordinates 
474   //  Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing
475   
476   static Double_t dist;
477
478   AliEMCALDigit * digit;
479   Int_t i=0, nstat=0, idMax=-1;
480   Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
481
482   //printf(" dist : %f  e : %f \n", dist, fAmp);
483   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
484     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
485     if(iDigit==0) {
486       idMax = digit->GetId(); // is it correct
487       dist  = TmaxInCm(Double_t(fAmp));
488     }
489     fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]);
490     //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), dist, xyzi[0], xyzi[1], xyzi[2]);
491
492     //fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
493     //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), 0.0, xyzi[0], xyzi[1], xyzi[2]);
494     //    if(fAmp>102.) assert(0);
495
496     if(logWeight > 0.0)  w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
497     else  w = fEnergyList[iDigit]; // just energy
498
499     if(w>0.0) {
500       wtot += w ;
501       nstat++;
502       for(i=0; i<3; i++ ) {
503         clXYZ[i]    += (w*xyzi[i]);
504         clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
505       }
506     }
507   }
508   //  cout << " wtot " << wtot << endl;
509   if ( wtot > 0 ) { 
510     //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
511     for(i=0; i<3; i++ ) {
512       clXYZ[i] /= wtot;
513       if(nstat>1) {
514         clRmsXYZ[i] /= (wtot*wtot);  
515         clRmsXYZ[i]  =  clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
516         if(clRmsXYZ[i] > 0.0) {
517           clRmsXYZ[i] =  TMath::Sqrt(clRmsXYZ[i]);
518         } else      clRmsXYZ[i] = 0;
519       } else        clRmsXYZ[i] = 0;
520     }    
521   } else {
522     for(i=0; i<3; i++ ) {
523       clXYZ[i] = clRmsXYZ[i] = -1.;
524     }
525   }
526   // clRmsXYZ[i] ??
527   fLocPos.SetX(clXYZ[0]);
528   fLocPos.SetY(clXYZ[1]);
529   fLocPos.SetZ(clXYZ[2]);
530     
531 //  if (gDebug==2)
532 //    printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; 
533   fLocPosM = 0 ; // covariance matrix
534 }
535
536 //void AliEMCALRecPoint::EvalLocalPositionSimple()
537 //{ // Weight is proportional of cell energy 
538 //}
539 //____________________________________________________________________________
540 void AliEMCALRecPoint::EvalLocalPositionFit(Double_t deff, Double_t logWeight, 
541 Double_t phiSlope, TClonesArray * digits)
542 {
543   // Aug 14-16, 2007 - for fit 
544   // Aug 31 - should be static ??
545   static Double_t dist, ycorr;
546   static AliEMCALDigit *digit;
547
548   Int_t i=0, nstat=0, idMax=-1;
549   Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.; 
550
551   for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
552     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
553     if(iDigit==0) {
554       idMax = digit->GetId(); // is it correct
555       dist  = TmaxInCm(Double_t(fAmp));
556     }
557
558     dist = deff;
559     fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]);
560
561     if(logWeight > 0.0)  w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
562     else                 w = fEnergyList[iDigit]; // just energy
563
564     if(w>0.0) {
565       wtot += w ;
566       nstat++;
567       for(i=0; i<3; i++ ) {
568         clXYZ[i]    += (w*xyzi[i]);
569         clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
570       }
571     }
572   }
573   //  cout << " wtot " << wtot << endl;
574   if ( wtot > 0 ) { 
575     //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
576     for(i=0; i<3; i++ ) {
577       clXYZ[i] /= wtot;
578       if(nstat>1) {
579         clRmsXYZ[i] /= (wtot*wtot);  
580         clRmsXYZ[i]  =  clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
581         if(clRmsXYZ[i] > 0.0) {
582           clRmsXYZ[i] =  TMath::Sqrt(clRmsXYZ[i]);
583         } else      clRmsXYZ[i] = 0;
584       } else        clRmsXYZ[i] = 0;
585     }    
586   } else {
587     for(i=0; i<3; i++ ) {
588       clXYZ[i] = clRmsXYZ[i] = -1.;
589     }
590   }
591   // clRmsXYZ[i] ??
592   if(phiSlope != 0.0 && logWeight > 0.0 && wtot) { 
593     // Correction in phi direction (y - coords here); Aug 16;
594     // May be put to global level or seperate method
595     ycorr = clXYZ[1] * (1. + phiSlope);
596     //printf(" y %f : ycorr %f : slope %f \n", clXYZ[1], ycorr, phiSlope); 
597     clXYZ[1] = ycorr;
598   }
599   fLocPos.SetX(clXYZ[0]);
600   fLocPos.SetY(clXYZ[1]);
601   fLocPos.SetZ(clXYZ[2]);
602     
603 //  if (gDebug==2)
604 //    printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; 
605   fLocPosM = 0 ; // covariance matrix
606 }
607
608 Bool_t AliEMCALRecPoint::EvalLocalPosition2(TClonesArray * digits, TArrayD &ed)
609 {
610   // Evaluated local position of rec.point using digits 
611   // and parametrisation of w0 and deff
612   //printf(" <I> AliEMCALRecPoint::EvalLocalPosition2() \n"); 
613   return AliEMCALRecPoint::EvalLocalPositionFromDigits(digits, ed, fLocPos);
614 }
615
616 Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(TClonesArray *digits, TArrayD &ed, TVector3 &locPos)
617 {
618   // Used when digits should be recalibrated
619   static Double_t deff, w0, esum;
620   static Int_t iDigit;
621   static AliEMCALDigit *digit;
622
623   if(ed.GetSize() && (digits->GetEntries()!=ed.GetSize())) return kFALSE;
624
625   // Calculate sum energy of digits
626   esum = 0.0;
627   for(iDigit=0; iDigit<ed.GetSize(); iDigit++) esum += ed[iDigit];
628
629   GetDeffW0(esum, deff, w0);
630   
631   return EvalLocalPositionFromDigits(esum, deff, w0, digits, ed, locPos); 
632 }
633
634 Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(const Double_t esum, const Double_t deff, const Double_t w0, 
635                                                      TClonesArray *digits, TArrayD &ed, TVector3 &locPos)
636 {
637   static AliEMCALDigit *digit;
638
639   Int_t i=0, nstat=0, idMax=-1;
640   Double_t clXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.; 
641
642   AliEMCALGeometry* geo = AliEMCALGeometry::GetInstance(); // Get pointer to EMCAL geometry
643    
644   for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
645     digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit));
646
647     geo->RelPosCellInSModule(digit->GetId(), idMax, deff, xyzi[0], xyzi[1], xyzi[2]);
648
649     if(w0 > 0.0)  w = TMath::Max( 0., w0 + TMath::Log(ed[iDigit] / esum));
650     else          w = ed[iDigit]; // just energy
651
652     if(w>0.0) {
653       wtot += w ;
654       nstat++;
655       for(i=0; i<3; i++ ) {
656         clXYZ[i] += (w*xyzi[i]);
657       }
658     }
659   }
660   //  cout << " wtot " << wtot << endl;
661   if (wtot > 0) { 
662     for(i=0; i<3; i++ ) {
663       clXYZ[i] /= wtot;
664     }
665     locPos.SetX(clXYZ[0]);
666     locPos.SetY(clXYZ[1]);
667     locPos.SetZ(clXYZ[2]);
668     return kTRUE;
669   } else {
670     return kFALSE;
671   }
672
673 }
674
675 void AliEMCALRecPoint::GetDeffW0(const Double_t esum , Double_t &deff,  Double_t &w0)
676 {
677   //
678   // Aug 31, 2001 
679   // Applied for simulation data with threshold 3 adc
680   // Calculate efective distance (deff) and weigh parameter (w0) 
681   // for coordinate calculation; 0.5 GeV < esum <100 GeV.
682   // Look to:  http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/CALIB/GEOMCORR/deffandW0VaEgamma_2.gif
683   //
684   static Double_t e=0.0;
685   const  Double_t dp0=9.25147, dp1=1.16700; // Hard coded now
686   const  Double_t wp0=4.83713, wp1=-2.77970e-01, wp2 = 4.41116;
687
688   // No extrapolation here
689   e = esum<0.5?0.5:esum;
690   e = e>100.?100.:e;
691
692   deff = dp0 + dp1*TMath::Log(e);
693   w0   = wp0 / (1. + TMath::Exp(wp1*(e+wp2)));
694   //printf("<I> AliEMCALRecPoint::GetDeffW0 esum %5.2f : deff %5.2f : w0 %5.2f \n", esum, deff, w0); 
695 }
696
697 //______________________________________________________________________________
698 void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
699 {
700   // This function calculates energy in the core, 
701   // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
702   // in accordance with shower profile the energy deposition 
703   // should be less than 2%
704   // Unfinished - Nov 15,2006
705   // Distance is calculate in (phi,eta) units
706
707   AliEMCALDigit * digit ;
708
709   Int_t iDigit;
710
711   if (!fLocPos.Mag()) {
712     EvalLocalPosition(logWeight, digits);
713   }
714   
715   Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta();
716   Double_t eta, phi, distance;
717   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
718     digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
719     
720     eta = phi = 0.0;
721     fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ;
722     phi = phi * TMath::DegToRad();
723   
724     distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint));
725     if(distance < fCoreRadius)
726       fCoreEnergy += fEnergyList[iDigit] ;
727   }
728   
729 }
730 //____________________________________________________________________________
731 void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
732 {
733   // Calculates the axis of the shower ellipsoid in eta and phi
734   // in cell units
735
736   static TString gn(fGeomPtr->GetName());
737
738   Double_t wtot = 0.;
739   Double_t x    = 0.;
740   Double_t z    = 0.;
741   Double_t dxx  = 0.;
742   Double_t dzz  = 0.;
743   Double_t dxz  = 0.;
744
745   AliEMCALDigit * digit = 0;
746
747   Double_t etai , phii, w; 
748   int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
749   int iphi=0, ieta=0;
750   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
751     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
752     etai = phii = 0.; 
753     if(gn.Contains("SHISH")) { 
754     // Nov 15,2006 - use cell numbers as coordinates
755     // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi
756     // We can use the eta,phi(or coordinates) of cell
757       nSupMod = nModule = nIphi = nIeta = iphi = ieta = 0;
758
759       fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
760       fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
761       etai=(Double_t)ieta;
762       phii=(Double_t)iphi;
763     } else { // 
764       fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
765       phii = phii * TMath::DegToRad();
766     }
767
768     w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
769     // fAmp summed amplitude of digits, i.e. energy of recpoint 
770     // Gives smaller value of lambda than log weight  
771     // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy
772
773     dxx  += w * etai * etai ;
774     x    += w * etai ;
775     dzz  += w * phii * phii ;
776     z    += w * phii ; 
777
778     dxz  += w * etai * phii ; 
779
780     wtot += w ;
781   }
782
783   if ( wtot > 0 ) { 
784     dxx /= wtot ;
785     x   /= wtot ;
786     dxx -= x * x ;
787     dzz /= wtot ;
788     z   /= wtot ;
789     dzz -= z * z ;
790     dxz /= wtot ;
791     dxz -= x * z ;
792
793     fLambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
794     if(fLambda[0] > 0)
795       fLambda[0] = TMath::Sqrt(fLambda[0]) ;
796     else
797       fLambda[0] = 0;
798     
799     fLambda[1] =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
800
801     if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
802       fLambda[1] = TMath::Sqrt(fLambda[1]) ;
803     else
804       fLambda[1]= 0. ;
805   } else { 
806     fLambda[0]= 0. ;
807     fLambda[1]= 0. ;
808   }
809
810   //  printf("Evalaxis: lambdas  = %f,%f", fLambda[0],fLambda[1]) ; 
811
812 }
813
814 //______________________________________________________________________________
815 void  AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
816 {
817   // Constructs the list of primary particles (tracks) which 
818   // have contributed to this RecPoint and calculate deposited energy 
819   // for each track
820   
821   AliEMCALDigit * digit ;
822   Int_t * primArray = new Int_t[fMaxTrack] ;
823   Float_t * dEPrimArray = new Float_t[fMaxTrack] ;
824
825   Int_t index ;  
826   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
827     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ; 
828     Int_t nprimaries = digit->GetNprimary() ;
829     if ( nprimaries == 0 ) continue ;
830     Int_t jndex ;
831     for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
832       if ( fMulTrack > fMaxTrack ) {
833         fMulTrack = fMaxTrack ;
834         Error("EvalPrimaries", "increase fMaxTrack ")  ;
835         break ;
836       }
837       Int_t newPrimary = digit->GetPrimary(jndex+1);
838       Float_t dEPrimary = digit->GetDEPrimary(jndex+1);
839       Int_t kndex ;
840       Bool_t already = kFALSE ;
841       for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
842         if ( newPrimary == primArray[kndex] ){
843           already = kTRUE ;
844           dEPrimArray[kndex] += dEPrimary; 
845           break ;
846         }
847       } // end of check
848       if ( !already && (fMulTrack < fMaxTrack)) { // store it
849         primArray[fMulTrack] = newPrimary ; 
850         dEPrimArray[fMulTrack] = dEPrimary ; 
851         fMulTrack++ ;
852       } // store it
853     } // all primaries in digit
854   } // all digits
855
856   Int_t *sortIdx = new Int_t[fMulTrack];
857   TMath::Sort(fMulTrack,dEPrimArray,sortIdx); 
858   for(index = 0; index < fMulTrack; index++) {
859     fTracksList[index] = primArray[sortIdx[index]] ;    
860     fDETracksList[index] = dEPrimArray[sortIdx[index]] ;
861   }
862   delete [] sortIdx;
863   delete [] primArray ;
864   delete [] dEPrimArray ;
865
866 }
867
868 //______________________________________________________________________________
869 void  AliEMCALRecPoint::EvalParents(TClonesArray * digits)
870 {
871   // Constructs the list of parent particles (tracks) which have contributed to this RecPoint
872  
873   AliEMCALDigit * digit ;
874   Int_t * parentArray = new Int_t[fMaxTrack] ;
875   Float_t * dEParentArray = new Float_t[fMaxTrack] ;
876
877   Int_t index ;  
878   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
879     if (fDigitsList[index] >= digits->GetEntries() || fDigitsList[index] < 0)
880        AliError(Form("Trying to get invalid digit %d (idx in WriteRecPoint %d)",fDigitsList[index],index));
881     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ; 
882     Int_t nparents = digit->GetNiparent() ;
883     if ( nparents == 0 ) continue ;
884
885     Int_t jndex ;
886     for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit
887       if ( fMulParent > fMaxParent ) {
888         fMulTrack = - 1 ;
889         Error("EvalParents", "increase fMaxParent")  ;
890         break ;
891       }
892       Int_t newParent = digit->GetIparent(jndex+1) ;
893       Float_t newdEParent = digit->GetDEParent(jndex+1) ;
894       Int_t kndex ;
895       Bool_t already = kFALSE ;
896       for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored
897         if ( newParent == parentArray[kndex] ){
898           dEParentArray[kndex] += newdEParent;
899           already = kTRUE ;
900           break ;
901         }
902       } // end of check
903       if ( !already && (fMulParent < fMaxParent)) { // store it
904         parentArray[fMulParent] = newParent ; 
905         dEParentArray[fMulParent] = newdEParent ; 
906         fMulParent++ ;
907       } // store it
908     } // all parents in digit
909   } // all digits
910
911   if (fMulParent>0) {
912     Int_t *sortIdx = new Int_t[fMulParent];
913     TMath::Sort(fMulParent,dEParentArray,sortIdx); 
914     for(index = 0; index < fMulParent; index++) {
915       fParentsList[index] = parentArray[sortIdx[index]] ;      
916       fDEParentsList[index] = dEParentArray[sortIdx[index]] ;
917     }
918     delete [] sortIdx;
919   }
920  
921   delete [] parentArray;
922   delete [] dEParentArray;
923 }
924
925 //____________________________________________________________________________
926 void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const
927 {
928   // returns the position of the cluster in the local reference system of ALICE
929   
930   lpos.SetX(fLocPos.X()) ;
931   lpos.SetY(fLocPos.Y()) ;
932   lpos.SetZ(fLocPos.Z()) ;
933 }
934
935 //____________________________________________________________________________
936 void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
937 {
938   // returns the position of the cluster in the global reference system of ALICE
939   // These are now the Cartesian X, Y and Z
940   //  cout<<" geom "<<geom<<endl;
941   fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
942 }
943
944 //____________________________________________________________________________
945 Float_t AliEMCALRecPoint::GetMaximalEnergy(void) const
946 {
947   // Finds the maximum energy in the cluster
948   
949   Float_t menergy = 0. ;
950
951   Int_t iDigit;
952
953   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
954  
955     if(fEnergyList[iDigit] > menergy) 
956       menergy = fEnergyList[iDigit] ;
957   }
958   return menergy ;
959 }
960
961 //____________________________________________________________________________
962 Int_t AliEMCALRecPoint::GetMultiplicityAtLevel(Float_t H) const
963 {
964   // Calculates the multiplicity of digits with energy larger than H*energy 
965   
966   Int_t multipl   = 0 ;
967   Int_t iDigit ;
968   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
969
970     if(fEnergyList[iDigit] > H * fAmp) 
971       multipl++ ;
972   }
973   return multipl ;
974 }
975
976 //____________________________________________________________________________
977 Int_t  AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t * maxAtEnergy,
978                                                Float_t locMaxCut,TClonesArray * digits) const
979
980   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
981   // energy difference between two local maxima
982
983   AliEMCALDigit * digit ;
984   AliEMCALDigit * digitN ;
985   
986   Int_t iDigitN ;
987   Int_t iDigit ;
988
989   for(iDigit = 0; iDigit < fMulDigit; iDigit++)
990     maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit])  ;
991   
992   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
993     if(maxAt[iDigit]) {
994       digit = maxAt[iDigit] ;
995           
996       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {        
997         digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ; 
998         
999         if ( AreNeighbours(digit, digitN) ) {
1000           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
1001             maxAt[iDigitN] = 0 ;
1002             // but may be digit too is not local max ?
1003             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut) 
1004               maxAt[iDigit] = 0 ;
1005           }
1006           else {
1007             maxAt[iDigit] = 0 ;
1008             // but may be digitN too is not local max ?
1009             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) 
1010               maxAt[iDigitN] = 0 ; 
1011           } 
1012         } // if Areneighbours
1013       } // while digitN
1014     } // slot not empty
1015   } // while digit
1016   
1017   iDigitN = 0 ;
1018   for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
1019     if(maxAt[iDigit] ){
1020       maxAt[iDigitN] = maxAt[iDigit] ;
1021       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
1022       iDigitN++ ; 
1023     }
1024   }
1025   return iDigitN ;
1026 }
1027
1028 //____________________________________________________________________________
1029 Int_t AliEMCALRecPoint::GetPrimaryIndex() const  
1030 {
1031   // Get the primary track index in TreeK which deposits the most energy 
1032   // in Digits which forms RecPoint. 
1033
1034   if (fMulTrack)
1035     return fTracksList[0];
1036   return -12345;
1037 }
1038
1039 //____________________________________________________________________________
1040 void AliEMCALRecPoint::EvalTime(TClonesArray * digits){
1041   // time is set to the time of the digit with the maximum energy
1042
1043   Float_t maxE = 0;
1044   Int_t maxAt = 0;
1045   for(Int_t idig=0; idig < fMulDigit; idig++){
1046     if(fEnergyList[idig] > maxE){
1047       maxE = fEnergyList[idig] ;
1048       maxAt = idig;
1049     }
1050   }
1051   fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
1052   
1053 }
1054
1055 //______________________________________________________________________________
1056 void AliEMCALRecPoint::Paint(Option_t *)
1057 {
1058   // Paint this ALiRecPoint as a TMarker  with its current attributes
1059   
1060   TVector3 pos(0.,0.,0.)  ;
1061   GetLocalPosition(pos)   ;
1062   Coord_t x = pos.X()     ;
1063   Coord_t y = pos.Z()     ;
1064   Color_t markercolor = 1 ;
1065   Size_t  markersize = 1. ;
1066   Style_t markerstyle = 5 ;
1067   
1068   if (!gPad->IsBatch()) {
1069     gVirtualX->SetMarkerColor(markercolor) ;
1070     gVirtualX->SetMarkerSize (markersize)  ;
1071     gVirtualX->SetMarkerStyle(markerstyle) ;
1072   }
1073   gPad->SetAttMarkerPS(markercolor,markerstyle,markersize) ;
1074   gPad->PaintPolyMarker(1,&x,&y,"") ;
1075 }
1076
1077 Double_t AliEMCALRecPoint::TmaxInCm(const Double_t e , const Int_t key)
1078
1079   // e energ in in GeV)
1080   // key  =  0(gamma, default)
1081   //     !=  0(electron)
1082   static Double_t ca = 4.82;  // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
1083   static Double_t X0 = 1.23;  // radiation lenght (cm)
1084   static Double_t tmax = 0.;   // position of electromagnetic shower max in cm
1085
1086   tmax = 0.0;
1087   if(e>0.1) {
1088     tmax = TMath::Log(e) + ca;
1089     if      (key==0) tmax += 0.5; 
1090     else             tmax -= 0.5;
1091     tmax *= X0; // convert to cm
1092   }
1093   return tmax;
1094 }
1095
1096 //______________________________________________________________________________
1097 Float_t AliEMCALRecPoint::EtaToTheta(Float_t arg) const
1098 {
1099   //Converts Theta (Radians) to Eta(Radians)
1100   return (2.*TMath::ATan(TMath::Exp(-arg)));
1101 }
1102
1103 //______________________________________________________________________________
1104 Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const
1105 {
1106   //Converts Eta (Radians) to Theta(Radians)
1107   return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
1108 }
1109
1110 //____________________________________________________________________________
1111 void AliEMCALRecPoint::Print(Option_t *opt) const
1112 {
1113   // Print the list of digits belonging to the cluster
1114   if(strlen(opt)==0) return;
1115   TString message ; 
1116   message  = "AliEMCALRecPoint:\n" ;
1117   message +=  " digits # = " ; 
1118   Info("Print", message.Data()) ; 
1119
1120   Int_t iDigit;
1121   for(iDigit=0; iDigit<fMulDigit; iDigit++)
1122     printf(" %d ", fDigitsList[iDigit] ) ;  
1123   printf("\n");
1124
1125   Info("Print", " Energies = ") ;
1126   for(iDigit=0; iDigit<fMulDigit; iDigit++) 
1127     printf(" %f ", fEnergyList[iDigit] ) ;
1128   printf("\n");
1129
1130   Info("Print", "\n Abs Ids = ") ;
1131   for(iDigit=0; iDigit<fMulDigit; iDigit++) 
1132     printf(" %i ", fAbsIdList[iDigit] ) ;
1133   printf("\n");
1134
1135   Info("Print", " Primaries  ") ;
1136   for(iDigit = 0;iDigit < fMulTrack; iDigit++)
1137     printf(" %d ", fTracksList[iDigit]) ;
1138
1139   printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]);
1140
1141   message  = "       ClusterType     = %d" ;
1142   message += "       Multiplicity    = %d" ;
1143   message += "       Cluster Energy  = %f" ; 
1144   message += "       Core energy     = %f" ; 
1145   message += "       Core radius     = %f" ; 
1146   message += "       Number of primaries %d" ; 
1147   message += "       Stored at position %d" ; 
1148   Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;  
1149 }
1150
1151 Double_t  AliEMCALRecPoint::GetPointEnergy() const
1152 {
1153   static double e;
1154   e=0.0;
1155   for(int ic=0; ic<GetMultiplicity(); ic++) e += double(fEnergyList[ic]);
1156   return e;
1157 }