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