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