new classes for track segments
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTowerRecPoint.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
16 /* $Id$ */
17
18 //_________________________________________________________________________
19 //  RecPoint implementation for EMCAL-EMC 
20 //  An TowerRecPoint is a cluster of digits   
21 //*--
22 //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
23
24
25 // --- ROOT system ---
26 #include "TPad.h"
27 #include "TH2.h"
28 #include "TMath.h" 
29 #include "TCanvas.h" 
30
31 // --- Standard library ---
32
33 // --- AliRoot header files ---
34
35  #include "AliGenerator.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliEMCALTowerRecPoint.h"
38 #include "AliRun.h"
39 #include "AliEMCALGetter.h"
40
41 ClassImp(AliEMCALTowerRecPoint)
42
43 //____________________________________________________________________________
44 AliEMCALTowerRecPoint::AliEMCALTowerRecPoint() : AliEMCALRecPoint()
45 {
46   // ctor
47
48   fMulDigit   = 0 ;  
49   fAmp   = 0. ;   
50   fCoreEnergy = 0 ; 
51   fEnergyList = 0 ;
52   fTime = 0. ;
53   fLocPos.SetX(0.)  ;      //Local position should be evaluated
54 }
55
56 //____________________________________________________________________________
57 AliEMCALTowerRecPoint::AliEMCALTowerRecPoint(const char * opt) : AliEMCALRecPoint(opt)
58 {
59   // ctor
60   
61   fMulDigit   = 0 ;  
62   fAmp   = 0. ;   
63   fCoreEnergy = 0 ; 
64   fEnergyList = 0 ;
65   fTime = -1. ;
66   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated  
67 }
68
69 //____________________________________________________________________________
70 AliEMCALTowerRecPoint::~AliEMCALTowerRecPoint()
71 {
72   // dtor
73
74   if ( fEnergyList )
75     delete[] fEnergyList ; 
76 }
77
78 //____________________________________________________________________________
79 void AliEMCALTowerRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
80 {
81   // Adds a digit to the RecPoint
82   // and accumulates the total amplitude and the multiplicity 
83   
84   if(fEnergyList == 0)
85     fEnergyList =  new Float_t[fMaxDigit]; 
86
87   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
88     fMaxDigit*=2 ; 
89     Int_t * tempo = new ( Int_t[fMaxDigit] ) ; 
90     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
91
92     Int_t index ;     
93     for ( index = 0 ; index < fMulDigit ; index++ ){
94       tempo[index]  = fDigitsList[index] ;
95       tempoE[index] = fEnergyList[index] ; 
96     }
97     
98     delete [] fDigitsList ; 
99     fDigitsList =  new ( Int_t[fMaxDigit] ) ;
100  
101     delete [] fEnergyList ;
102     fEnergyList =  new ( Float_t[fMaxDigit] ) ;
103
104     for ( index = 0 ; index < fMulDigit ; index++ ){
105       fDigitsList[index] = tempo[index] ;
106       fEnergyList[index] = tempoE[index] ; 
107     }
108  
109     delete [] tempo ;
110     delete [] tempoE ; 
111   } // if
112   
113   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ; 
114   fEnergyList[fMulDigit]   = Energy ;
115   fMulDigit++ ; 
116   fAmp += Energy ; 
117
118   // EvalEMCALMod(&digit) ;
119 }
120
121 //____________________________________________________________________________
122 Bool_t AliEMCALTowerRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
123 {
124   // Tells if (true) or not (false) two digits are neighbors
125   
126   Bool_t aren = kFALSE ;
127   
128   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
129   AliEMCALGeometry * phosgeom =  (AliEMCALGeometry*)gime->EMCALGeometry();
130
131   Int_t relid1[4] ; 
132   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
133
134   Int_t relid2[4] ; 
135   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
136   
137   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
138   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
139
140   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
141     aren = kTRUE ;
142   
143   return aren ;
144 }
145
146 //____________________________________________________________________________
147 Int_t AliEMCALTowerRecPoint::Compare(const TObject * obj) const
148 {
149   // Compares two RecPoints according to their position in the EMCAL modules
150
151   Float_t delta = 1 ; //Width of "Sorting row". If you change this 
152                       //value (what is senseless) change as vell delta in
153                       //AliEMCALTrackSegmentMakerv* and other RecPoints...
154   Int_t rv ; 
155
156   AliEMCALTowerRecPoint * clu = (AliEMCALTowerRecPoint *)obj ; 
157
158  
159   Int_t phosmod1 = GetEMCALArm() ;
160   Int_t phosmod2 = clu->GetEMCALArm() ;
161
162   TVector3 locpos1; 
163   GetLocalPosition(locpos1) ;
164   TVector3 locpos2;  
165   clu->GetLocalPosition(locpos2) ;  
166
167   if(phosmod1 == phosmod2 ) {
168     Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
169     if (rowdif> 0) 
170       rv = 1 ;
171      else if(rowdif < 0) 
172        rv = -1 ;
173     else if(locpos1.Z()>locpos2.Z()) 
174       rv = -1 ;
175     else 
176       rv = 1 ; 
177      }
178
179   else {
180     if(phosmod1 < phosmod2 ) 
181       rv = -1 ;
182     else 
183       rv = 1 ;
184   }
185
186   return rv ; 
187 }
188 //______________________________________________________________________________
189 void AliEMCALTowerRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
190 {
191   
192   // Execute action corresponding to one event
193   //  This member function is called when a AliEMCALRecPoint is clicked with the locator
194   //
195   //  If Left button is clicked on AliEMCALRecPoint, the digits are switched on    
196   //  and switched off when the mouse button is released.
197   
198     
199   // AliEMCALGetter * gime =  AliEMCALGetter::GetInstance() ; 
200 //   if(!gime) return ;
201 //   AliEMCALGeometry * emcalgeom =  (AliEMCALGeometry*)gime->EMCALGeometry();
202   
203 //   static TGraph *  digitgraph = 0 ;
204   
205 //   if (!gPad->IsEditable()) return;
206   
207 //   TH2F * histo = 0 ;
208 //   TCanvas * histocanvas ; 
209
210 //   const TClonesArray * digits = gime->Digits() ;
211   
212 //   switch (event) {
213     
214 //   case kButton1Down: {
215 //     AliEMCALDigit * digit ;
216 //     Int_t iDigit;
217 //     Int_t relid[4] ;
218     
219 //     const Int_t kMulDigit = AliEMCALTowerRecPoint::GetDigitsMultiplicity() ; 
220 //     Float_t * xi = new Float_t[kMulDigit] ; 
221 //     Float_t * zi = new Float_t[kMulDigit] ; 
222     
223 //     // create the histogram for the single cluster 
224 //     // 1. gets histogram boundaries
225 //     Float_t ximax = -999. ; 
226 //     Float_t zimax = -999. ; 
227 //     Float_t ximin = 999. ; 
228 //     Float_t zimin = 999. ;
229     
230 //     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
231 //       digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
232 //       emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
233 //       emcalgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
234 //       if ( xi[iDigit] > ximax )
235 //      ximax = xi[iDigit] ; 
236 //       if ( xi[iDigit] < ximin )
237 //      ximin = xi[iDigit] ; 
238 //       if ( zi[iDigit] > zimax )
239 //      zimax = zi[iDigit] ; 
240 //       if ( zi[iDigit] < zimin )
241 //      zimin = zi[iDigit] ;     
242 //     }
243 //     ximax += emcalgeom->GetCrystalSize(0) / 2. ;
244 //     zimax += emcalgeom->GetCrystalSize(2) / 2. ;
245 //     ximin -= emcalgeom->GetCrystalSize(0) / 2. ;
246 //     zimin -= emcalgeom->GetCrystalSize(2) / 2. ;
247 //     Int_t xdim = (int)( (ximax - ximin ) / emcalgeom->GetCrystalSize(0) + 0.5  ) ; 
248 //     Int_t zdim = (int)( (zimax - zimin ) / emcalgeom->GetCrystalSize(2) + 0.5 ) ;
249     
250 //     // 2. gets the histogram title
251     
252 //     Text_t title[100] ; 
253 //     sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
254     
255 //     if (!histo) {
256 //       delete histo ; 
257 //       histo = 0 ; 
258 //     }
259 //     histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
260     
261 //     Float_t x, z ; 
262 //     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
263 //       digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
264 //       emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
265 //       emcalgeom->RelPosInModule(relid, x, z);
266 //       histo->Fill(x, z, fEnergyList[iDigit] ) ;
267 //     }
268     
269 //     if (!digitgraph) {
270 //       digitgraph = new TGraph(kMulDigit,xi,zi);
271 //       digitgraph-> SetMarkerStyle(5) ; 
272 //       digitgraph-> SetMarkerSize(1.) ;
273 //       digitgraph-> SetMarkerColor(1) ;
274 //       digitgraph-> Paint("P") ;
275 //     }
276     
277 //     //    Print() ;
278 //     histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ; 
279 //     histocanvas->Draw() ; 
280 //     histo->Draw("lego1") ; 
281     
282 //     delete[] xi ; 
283 //     delete[] zi ; 
284     
285 //     break;
286 //   }
287   
288 //   case kButton1Up: 
289 //     if (digitgraph) {
290 //       delete digitgraph  ;
291 //       digitgraph = 0 ;
292 //     }
293 //     break;
294   
295 //    }
296 }
297
298 //____________________________________________________________________________
299 void  AliEMCALTowerRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits)
300 {
301   // Calculates the dispersion of the shower at the origine of the RecPoint
302
303   Float_t d    = 0. ;
304   Float_t wtot = 0. ;
305
306   AliEMCALDigit * digit ;
307  
308   AliEMCALGeometry * emcalgeom = AliEMCALGetter::GetInstance()->EMCALGeometry();
309   
310
311  // Calculates the center of gravity in the local EMCAL-module coordinates 
312   
313   Int_t iDigit;
314
315   if (!fTheta || !fPhi ) 
316     EvalGlobalPosition(logWeight, digits) ;
317   
318   const Float_t kDeg2Rad = TMath::DegToRad() ; 
319     
320   Float_t cyl_radius = 0 ;  
321   
322   if (IsInPRE()) 
323     cyl_radius = emcalgeom->GetIP2PRESection() ;
324   else if (IsInECAL()) 
325     cyl_radius = emcalgeom->GetIP2ECALSection() ;
326   else if (IsInHCAL()) 
327     cyl_radius = emcalgeom->GetIP2HCALSection() ;
328   else 
329     Fatal("EvalDispersion", "Unexpected tower section!") ; 
330   
331   Float_t x =  fLocPos.X() ; 
332   Float_t y =  fLocPos.Y() ; 
333   Float_t z =  fLocPos.Z() ; 
334   
335   if (gDebug == 2) 
336     Info("EvalDispersion", "x,y,z = %f,%f,%f", x, y, z) ;
337
338 // Calculates the dispersion in coordinates 
339   wtot = 0.;
340   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
341     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
342     Float_t thetai = 0. ;
343     Float_t phii = 0.;
344     emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
345     Float_t xi =  cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
346     Float_t yi =  cyl_radius * TMath::Sin(phii * kDeg2Rad ) ; 
347     Float_t zi =  cyl_radius / TMath::Tan(thetai * kDeg2Rad ) ; 
348
349     if (gDebug == 2) 
350       Info("EvalDispersion", "id = %d, xi,yi,zi = %f,%f,%f", digit->GetId(), xi, yi, zi) ;
351
352     Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
353     d += w * ( (xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
354     wtot+=w ;
355   }
356   
357   if ( wtot > 0 ) 
358     d /= wtot ;
359   else 
360     d = 0. ; 
361
362   fDispersion = TMath::Sqrt(d) ;
363  
364 }
365 //______________________________________________________________________________
366 void AliEMCALTowerRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
367 {
368   // This function calculates energy in the core, 
369   // i.e. within a radius rad = 3cm around the center. Beyond this radius
370   // in accordance with shower profile the energy deposition 
371   // should be less than 2%
372
373   Float_t coreRadius = 10. ;
374
375   AliEMCALDigit * digit ;
376   Float_t wtot = 0. ;
377
378   AliEMCALGeometry * emcalgeom = AliEMCALGetter::GetInstance()->EMCALGeometry();    
379   Int_t iDigit;
380
381   if (!fTheta || !fPhi ) {
382     for(iDigit=0; iDigit<fMulDigit; iDigit++) {
383       digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
384       
385       Float_t thetai ;
386       Float_t phii ;
387       emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
388       Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
389       fTheta = fTheta + thetai * w ;
390       fPhi   += (phii * w );
391       wtot  += w ;
392     }
393     
394     if (wtot > 0 ) { 
395       fTheta /= wtot ;
396       fPhi   /= wtot ;
397     } else { 
398       fTheta = -1 ; 
399       fPhi   = -1 ; 
400     }
401   }
402   
403   const Float_t kDeg2Rad = TMath::DegToRad() ; 
404   
405   Float_t cyl_radius = emcalgeom->GetIP2ECALSection();
406   Float_t x =  cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
407   Float_t y =  cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ; 
408   Float_t z =  cyl_radius * TMath::Tan(fTheta * kDeg2Rad ) ; 
409
410   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
411     digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
412     Float_t thetai = 0. ;
413     Float_t phii = 0. ;
414     emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
415     
416     Float_t xi =  cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
417     Float_t yi =  cyl_radius * TMath::Sin(phii * kDeg2Rad ) ; 
418     Float_t zi =  cyl_radius * TMath::Tan(thetai * kDeg2Rad ) ; 
419      
420     Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(yi-y)*(yi-y)+(zi-z)*(zi-z)) ;
421     if(distance < coreRadius)
422       fCoreEnergy += fEnergyList[iDigit] ;
423   }
424   
425 }
426
427 //____________________________________________________________________________
428 void  AliEMCALTowerRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
429 {
430   // Calculates the axis of the shower ellipsoid
431
432   Double_t wtot = 0. ;
433   Double_t x    = 0.;
434   Double_t z    = 0.;
435   Double_t dxx  = 0.;
436   Double_t dzz  = 0.;
437   Double_t dxz  = 0.;
438
439   AliEMCALDigit * digit ;
440
441   AliEMCALGeometry * emcalgeom = AliEMCALGetter::GetInstance()->EMCALGeometry();
442
443   Int_t iDigit;
444   const Float_t kDeg2Rad = TMath::DegToRad() ; 
445   
446    Float_t cyl_radius = 0 ;  
447   
448   if (IsInPRE()) 
449     cyl_radius = emcalgeom->GetIP2PRESection() ;
450   else if (IsInECAL()) 
451     cyl_radius = emcalgeom->GetIP2ECALSection() ;
452   else if (IsInHCAL()) 
453     cyl_radius = emcalgeom->GetIP2HCALSection() ;
454   else 
455     Fatal("EvalDispersion", "Unexpected tower section!") ; 
456
457   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
458     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
459     Float_t thetai = 0. ;
460     Float_t phii = 0. ; 
461     emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
462     Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
463     Float_t xi =  cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
464     Float_t zi =  cyl_radius / TMath::Tan(fTheta * kDeg2Rad ) ; 
465     dxx  += w * xi * xi ;
466     x    += w * xi ;
467     dzz  += w * zi * zi ;
468     z    += w * zi ; 
469     dxz  += w * xi * zi ; 
470     wtot += w ;
471   }
472   if ( wtot > 0 ) { 
473     dxx /= wtot ;
474     x   /= wtot ;
475     dxx -= x * x ;
476     dzz /= wtot ;
477     z   /= wtot ;
478     dzz -= z * z ;
479     dxz /= wtot ;
480     dxz -= x * z ;
481
482     
483     //   //Apply correction due to non-perpendicular incidence
484 //   Double_t CosX ;
485 //   Double_t CosZ ;
486 //   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
487 //   AliEMCALGeometry * emcalgeom =  (AliEMCALGeometry*)gime->EMCALGeometry();
488   //   Double_t DistanceToIP= (Double_t ) emcalgeom->GetIPDistance() ;
489   
490 //   CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
491 //   CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
492
493 //   dxx = dxx/(CosX*CosX) ;
494 //   dzz = dzz/(CosZ*CosZ) ;
495 //   dxz = dxz/(CosX*CosZ) ;
496
497
498     fLambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
499     if(fLambda[0] > 0)
500       fLambda[0] = TMath::Sqrt(fLambda[0]) ;
501     
502     fLambda[1] =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
503     if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
504       fLambda[1] = TMath::Sqrt(fLambda[1]) ;
505     else
506       fLambda[1]= 0. ;
507   } else { 
508     fLambda[0]= 0. ;
509     fLambda[1]= 0. ;
510   }
511 }
512
513 //____________________________________________________________________________
514 void AliEMCALTowerRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits )
515 {
516   // Evaluates all shower parameters
517
518   AliEMCALRecPoint::EvalAll(logWeight,digits) ;
519   EvalGlobalPosition(logWeight, digits) ;
520   EvalElipsAxis(logWeight, digits) ;
521   EvalDispersion(logWeight, digits) ;
522   EvalCoreEnergy(logWeight, digits);
523   EvalTime(digits) ;
524 }
525
526 //____________________________________________________________________________
527 void AliEMCALTowerRecPoint::EvalGlobalPosition(Float_t logWeight, TClonesArray * digits)
528 {
529   // Calculates the center of gravity in the local EMCAL-module coordinates 
530   Float_t wtot = 0. ;
531  
532   //  Int_t relid[4] ;
533   
534   AliEMCALDigit * digit ;
535   AliEMCALGeometry * emcalgeom  =  AliEMCALGetter::GetInstance()->EMCALGeometry();
536   Int_t iDigit;
537
538   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
539     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
540
541     Float_t thetai ;
542     Float_t phii ;
543     emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
544     Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
545     fTheta = fTheta + thetai * w ;
546     fPhi   += (phii * w );
547     wtot  += w ;
548   }
549
550   if ( wtot > 0 ) { 
551     fTheta /= wtot ;
552     fPhi   /= wtot ;
553   } else {
554     fTheta = -1 ; 
555     fPhi   = -1.; 
556   }
557   
558
559   const Float_t kDeg2Rad = TMath::DegToRad() ; 
560   
561   Float_t cyl_radius = 0 ;  
562
563   if (IsInPRE()) 
564     cyl_radius = emcalgeom->GetIP2PRESection() ;
565   else if (IsInECAL()) 
566     cyl_radius = emcalgeom->GetIP2ECALSection() ;
567   else if (IsInHCAL()) 
568     cyl_radius = emcalgeom->GetIP2HCALSection() ;
569   else 
570     Fatal("EvalGlobalPosition", "Unexpected tower section!") ; 
571   
572   Float_t x =  cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
573   Float_t y =  cyl_radius * TMath::Sin(fPhi * kDeg2Rad ) ; 
574   Float_t z =  cyl_radius / TMath::Tan(fTheta * kDeg2Rad ) ; 
575   
576   fLocPos.SetX(x)  ;
577   fLocPos.SetY(y) ;
578   fLocPos.SetZ(z)  ;
579     
580   if (gDebug==2)
581     Info("EvalGlobalPosition", "x,y,z = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; 
582
583
584   fLocPosM = 0 ;
585 }
586
587 //____________________________________________________________________________
588 Float_t AliEMCALTowerRecPoint::GetMaximalEnergy(void) const
589 {
590   // Finds the maximum energy in the cluster
591   
592   Float_t menergy = 0. ;
593
594   Int_t iDigit;
595
596   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
597  
598     if(fEnergyList[iDigit] > menergy) 
599       menergy = fEnergyList[iDigit] ;
600   }
601   return menergy ;
602 }
603
604 //____________________________________________________________________________
605 Int_t AliEMCALTowerRecPoint::GetMultiplicityAtLevel(const Float_t H) const
606 {
607   // Calculates the multiplicity of digits with energy larger than H*energy 
608   
609   Int_t multipl   = 0 ;
610   Int_t iDigit ;
611   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
612
613     if(fEnergyList[iDigit] > H * fAmp) 
614       multipl++ ;
615   }
616   return multipl ;
617 }
618
619 //____________________________________________________________________________
620 Int_t  AliEMCALTowerRecPoint::GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t * maxAtEnergy,
621                                                Float_t locMaxCut,TClonesArray * digits) const
622
623   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
624   // energy difference between two local maxima
625
626   AliEMCALDigit * digit ;
627   AliEMCALDigit * digitN ;
628   
629
630   Int_t iDigitN ;
631   Int_t iDigit ;
632
633   for(iDigit = 0; iDigit < fMulDigit; iDigit++)
634     maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit])  ;
635
636   
637   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
638     if(maxAt[iDigit]) {
639       digit = maxAt[iDigit] ;
640           
641       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {        
642         digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ; 
643         
644         if ( AreNeighbours(digit, digitN) ) {
645           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
646             maxAt[iDigitN] = 0 ;
647             // but may be digit too is not local max ?
648             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut) 
649               maxAt[iDigit] = 0 ;
650           }
651           else {
652             maxAt[iDigit] = 0 ;
653             // but may be digitN too is not local max ?
654             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) 
655               maxAt[iDigitN] = 0 ; 
656           } 
657         } // if Areneighbours
658       } // while digitN
659     } // slot not empty
660   } // while digit
661   
662   iDigitN = 0 ;
663   for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
664     if(maxAt[iDigit] ){
665       maxAt[iDigitN] = maxAt[iDigit] ;
666       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
667       iDigitN++ ; 
668     }
669   }
670   return iDigitN ;
671 }
672 //____________________________________________________________________________
673 void AliEMCALTowerRecPoint::EvalTime(TClonesArray * digits){
674   
675   Float_t maxE = 0;
676   Int_t maxAt = 0;
677   for(Int_t idig=0; idig < fMulDigit; idig++){
678     if(fEnergyList[idig] > maxE){
679       maxE = fEnergyList[idig] ;
680       maxAt = idig;
681     }
682   }
683   fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
684   
685 }
686 //____________________________________________________________________________
687 void AliEMCALTowerRecPoint::Print(Option_t * option) 
688 {
689   // Print the list of digits belonging to the cluster
690   
691   TString message("\n") ; 
692
693   Int_t iDigit;
694   message += "digits # = " ;
695   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
696     message += fDigitsList[iDigit] ; 
697     message += "  " ;
698   } 
699   
700   message += "\nEnergies = " ;
701   for(iDigit=0; iDigit<fMulDigit; iDigit++) { 
702     message += fEnergyList[iDigit] ; 
703     message += "  " ;
704   }
705   
706    message += "\nPrimaries  " ;
707    for(iDigit = 0;iDigit < fMulTrack; iDigit++) {
708      message += fTracksList[iDigit] ;
709      message += " " ;
710    }
711    message += "\n       Multiplicity    = " ; 
712    message += fMulDigit ;
713    message += "\n       Cluster Energy  = " ; 
714    message += fAmp ;
715    message += "\n       Number of primaries " ; 
716    message += fMulTrack ;
717    message += "\n       Stored at position " ;
718    message += GetIndexInList() ; 
719    
720    Info("Print", message.Data() ) ; 
721 }
722  
723 //____________________________________________________________________________
724 const TVector3 AliEMCALTowerRecPoint::XYZInAlice(Float_t r, Float_t theta, Float_t phi) const 
725 {
726   // spherical coordinates of recpoint in Alice reference frame
727
728   if (gDebug == 2) 
729     Info("XYZInAlice", "this= %d , r = %f, theta = %f, phi = %f", this, r, theta, phi) ; 
730
731   if (theta == 9999. || phi == 9999. || r == 9999.) {
732     TVector3  globalpos;  
733     GetGlobalPosition(globalpos);
734     phi   =  globalpos.X() * TMath::DegToRad() ; 
735     r     =  globalpos.Y() ; 
736     theta =  globalpos.Z() * TMath::DegToRad() ; 
737   }
738   else {
739     theta *= TMath::DegToRad() ; 
740     phi   *= TMath::DegToRad() ; 
741   }
742   
743   Float_t y = r * TMath::Cos(phi) ; 
744   Float_t x = r * TMath::Sin(phi) * TMath::Sin(theta) ; 
745   Float_t z = r * TMath::Sin(phi) * TMath::Cos(theta) ; 
746   
747   TVector3 vec(z, x, y) ; 
748   return vec ; 
749 }