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