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