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