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