Initialization of some data members (Christian)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPoint.cxx
CommitLineData
ab48128d 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//_________________________________________________________________________
70a93198 17// Reconstructed Points for the EMCAL
18// A RecPoint is a cluster of digits
d64c959b 19//*-- Author: Yves Schutz (SUBATECH)
70a93198 20//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
21//*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
ab48128d 22
23// --- ROOT system ---
1d59832c 24class Riostream;
e52475ed 25#include <TPad.h>
1d59832c 26class TGraph;
27class TPaveText;
e52475ed 28#include <TClonesArray.h>
29#include <TMath.h>
ab48128d 30
31// --- Standard library ---
ab48128d 32
33// --- AliRoot header files ---
1d59832c 34//#include "AliGenerator.h"
35class AliGenerator;
4635df1f 36#include "AliRunLoader.h"
37#include "AliRun.h"
1d59832c 38class AliEMCAL;
4635df1f 39#include "AliEMCALLoader.h"
ab48128d 40#include "AliEMCALGeometry.h"
4635df1f 41#include "AliEMCALHit.h"
ab48128d 42#include "AliEMCALDigit.h"
43#include "AliEMCALRecPoint.h"
ab48128d 44
45ClassImp(AliEMCALRecPoint)
46
ab48128d 47//____________________________________________________________________________
48AliEMCALRecPoint::AliEMCALRecPoint()
18a21c7c 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
af5bdd85 59 fDETracksList(0),
18a21c7c 60 fMulParent(0),
61 fMaxParent(0),
62 fParentsList(0),
af5bdd85 63 fDEParentsList(0),
1ae500a2 64 fSuperModuleNumber(0),
65 fDigitIndMax(-1)
ab48128d 66{
67 // ctor
4635df1f 68 AliRunLoader *rl = AliRunLoader::GetRunLoader();
e1a51e6e 69 if (rl && rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
af5bdd85 70 fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
71 else
e1a51e6e 72 fGeomPtr = AliEMCALGeometry::GetInstance();
73 //fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
706863b6 74 fLambda[0] = 0;
75 fLambda[1] = 0;
e1a51e6e 76 // fGeomPtr->GetTransformationForSM(); // Global <-> Local
ab48128d 77}
78
79//____________________________________________________________________________
18a21c7c 80AliEMCALRecPoint::AliEMCALRecPoint(const char * opt)
81 : AliRecPoint(opt),
82 fGeomPtr(0),
83 fClusterType(-1),
84 fCoreEnergy(0),
85 fDispersion(0),
86 fEnergyList(0),
87 fTimeList(0),
88 fAbsIdList(0),
89 fTime(-1.),
90 fCoreRadius(10), //HG check this
af5bdd85 91 fDETracksList(0),
18a21c7c 92 fMulParent(0),
93 fMaxParent(1000),
94 fParentsList(0),
af5bdd85 95 fDEParentsList(0),
1ae500a2 96 fSuperModuleNumber(0),
97 fDigitIndMax(-1)
ab48128d 98{
99 // ctor
94478418 100 // Increase fMaxTrack for EMCAL.
101 delete [] fTracksList;
ff1e7e2f 102 fMaxTrack = 1000 ;
94478418 103 fTracksList = new Int_t[fMaxTrack];
af5bdd85 104 fDETracksList = new Float_t[fMaxTrack];
94478418 105
87cdc3be 106 fParentsList = new Int_t[fMaxParent];
af5bdd85 107 fDEParentsList = new Float_t[fMaxParent];
108 for (Int_t i = 0; i < fMaxTrack; i++)
109 fDETracksList[i] = 0;
94478418 110 for (Int_t i = 0; i < fMaxParent; i++) {
111 fParentsList[i] = -1;
af5bdd85 112 fDEParentsList[i] = 0;
94478418 113 }
18a21c7c 114
4635df1f 115 AliRunLoader *rl = AliRunLoader::GetRunLoader();
e1a51e6e 116 if (rl && rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
af5bdd85 117 fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
118 else
e1a51e6e 119 fGeomPtr = AliEMCALGeometry::GetInstance();
706863b6 120 fLambda[0] = 0;
121 fLambda[1] = 0;
e1a51e6e 122 // fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
123 // fGeomPtr->GetTransformationForSM(); // Global <-> Local
70a93198 124}
18a21c7c 125
70a93198 126//____________________________________________________________________________
18a21c7c 127AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp)
128 : AliRecPoint(rp),
129 fGeomPtr(rp.fGeomPtr),
130 fClusterType(rp.fClusterType),
131 fCoreEnergy(rp.fCoreEnergy),
132 fDispersion(rp.fDispersion),
133 fEnergyList(0),
134 fTimeList(0),
135 fAbsIdList(0),
136 fTime(rp.fTime),
137 fCoreRadius(rp.fCoreRadius),
af5bdd85 138 fDETracksList(0),
18a21c7c 139 fMulParent(rp.fMulParent),
140 fMaxParent(rp.fMaxParent),
141 fParentsList(0),
af5bdd85 142 fDEParentsList(0),
1ae500a2 143 fSuperModuleNumber(rp.fSuperModuleNumber),
144 fDigitIndMax(rp.fDigitIndMax)
0a4cb131 145{
146 //copy ctor
0a4cb131 147 fLambda[0] = rp.fLambda[0];
148 fLambda[1] = rp.fLambda[1];
18a21c7c 149
af5bdd85 150 fEnergyList = new Float_t[rp.fMaxDigit];
151 fTimeList = new Float_t[rp.fMaxDigit];
152 fAbsIdList = new Int_t[rp.fMaxDigit];
0a4cb131 153 for(Int_t i = 0; i < rp.fMulDigit; i++) {
154 fEnergyList[i] = rp.fEnergyList[i];
155 fTimeList[i] = rp.fTimeList[i];
156 fAbsIdList[i] = rp.fAbsIdList[i];
157 }
af5bdd85 158 fDETracksList = new Float_t[rp.fMaxTrack];
159 for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i];
160 fParentsList = new Int_t[rp.fMaxParent];
0a4cb131 161 for(Int_t i = 0; i < rp.fMulParent; i++) fParentsList[i] = rp.fParentsList[i];
af5bdd85 162 fDEParentsList = new Float_t[rp.fMaxParent];
163 for(Int_t i = 0; i < rp.fMulParent; i++) fDEParentsList[i] = rp.fDEParentsList[i];
0a4cb131 164
165}
166//____________________________________________________________________________
70a93198 167AliEMCALRecPoint::~AliEMCALRecPoint()
168{
169 // dtor
170 if ( fEnergyList )
171 delete[] fEnergyList ;
85c60a8e 172 if ( fTimeList )
173 delete[] fTimeList ;
e52475ed 174 if ( fAbsIdList )
175 delete[] fAbsIdList ;
af5bdd85 176 if ( fDETracksList)
177 delete[] fDETracksList;
87cdc3be 178 if ( fParentsList)
179 delete[] fParentsList;
af5bdd85 180 if ( fDEParentsList)
181 delete[] fDEParentsList;
70a93198 182}
183
184//____________________________________________________________________________
185void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
186{
187 // Adds a digit to the RecPoint
188 // and accumulates the total amplitude and the multiplicity
189
190 if(fEnergyList == 0)
191 fEnergyList = new Float_t[fMaxDigit];
85c60a8e 192 if(fTimeList == 0)
193 fTimeList = new Float_t[fMaxDigit];
e52475ed 194 if(fAbsIdList == 0) {
195 fAbsIdList = new Int_t[fMaxDigit];
14ce0a6e 196 fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
e52475ed 197 }
70a93198 198
199 if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
200 fMaxDigit*=2 ;
a64a06d6 201 Int_t * tempo = new Int_t[fMaxDigit];
202 Float_t * tempoE = new Float_t[fMaxDigit];
85c60a8e 203 Float_t * tempoT = new Float_t[fMaxDigit];
e52475ed 204 Int_t * tempoId = new Int_t[fMaxDigit];
70a93198 205
206 Int_t index ;
207 for ( index = 0 ; index < fMulDigit ; index++ ){
e52475ed 208 tempo[index] = fDigitsList[index] ;
209 tempoE[index] = fEnergyList[index] ;
85c60a8e 210 tempoT[index] = fTimeList[index] ;
e52475ed 211 tempoId[index] = fAbsIdList[index] ;
70a93198 212 }
213
94478418 214 delete [] fDigitsList ;
70a93198 215 delete [] fEnergyList ;
85c60a8e 216 delete [] fTimeList ;
e52475ed 217 delete [] fAbsIdList ;
e52475ed 218
94478418 219 fDigitsList = tempo;
220 fEnergyList = tempoE;
221 fTimeList = tempoT;
222 fAbsIdList = tempoId;
70a93198 223 } // if
224
225 fDigitsList[fMulDigit] = digit.GetIndexInList() ;
226 fEnergyList[fMulDigit] = Energy ;
94478418 227 fTimeList[fMulDigit] = digit.GetTimeR() ;
e52475ed 228 fAbsIdList[fMulDigit] = digit.GetId();
70a93198 229 fMulDigit++ ;
230 fAmp += Energy ;
231
232}
233//____________________________________________________________________________
234Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
235{
236 // Tells if (true) or not (false) two digits are neighbours
237 // A neighbour is defined as being two digits which share a corner
238
e52475ed 239 static Bool_t areNeighbours = kFALSE ;
2bb3725c 240 static Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
241 static int nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
e52475ed 242 static Int_t relid1[2] , relid2[2] ; // ieta, iphi
243 static Int_t rowdiff=0, coldiff=0;
70a93198 244
e52475ed 245 areNeighbours = kFALSE ;
246
2bb3725c 247 fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nModule,nIphi,nIeta);
248 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
e52475ed 249
2bb3725c 250 fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
251 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
70a93198 252
e52475ed 253 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
254 coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
70a93198 255
256 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
e52475ed 257 areNeighbours = kTRUE ;
ab48128d 258
70a93198 259 return areNeighbours;
260}
261
262//____________________________________________________________________________
263Int_t AliEMCALRecPoint::Compare(const TObject * obj) const
264{
265 // Compares two RecPoints according to their position in the EMCAL modules
266
267 Float_t delta = 1 ; //Width of "Sorting row". If you change this
268 //value (what is senseless) change as well delta in
269 //AliEMCALTrackSegmentMakerv* and other RecPoints...
270 Int_t rv ;
271
272 AliEMCALRecPoint * clu = (AliEMCALRecPoint *)obj ;
273
274 TVector3 locpos1;
275 GetLocalPosition(locpos1);
276 TVector3 locpos2;
277 clu->GetLocalPosition(locpos2);
278
9848d950 279 Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ;
70a93198 280 if (rowdif> 0)
281 rv = 1 ;
282 else if(rowdif < 0)
283 rv = -1 ;
284 else if(locpos1.Y()>locpos2.Y())
285 rv = -1 ;
286 else
287 rv = 1 ;
288
289 return rv ;
ab48128d 290}
291
292//____________________________________________________________________________
293Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py)
294{
295 // Compute distance from point px,py to a AliEMCALRecPoint considered as a Tmarker
296 // Compute the closest distance of approach from point px,py to this marker.
297 // The distance is computed in pixels units.
70a93198 298 // HG Still need to update -> Not sure what this should achieve
ab48128d 299
300 TVector3 pos(0.,0.,0.) ;
70a93198 301 GetLocalPosition(pos) ;
ab48128d 302 Float_t x = pos.X() ;
70a93198 303 Float_t y = pos.Y() ;
ab48128d 304 const Int_t kMaxDiff = 10;
305 Int_t pxm = gPad->XtoAbsPixel(x);
306 Int_t pym = gPad->YtoAbsPixel(y);
307 Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym);
308
309 if (dist > kMaxDiff) return 9999;
310 return dist;
311}
312
313//___________________________________________________________________________
314 void AliEMCALRecPoint::Draw(Option_t *option)
315 {
316 // Draw this AliEMCALRecPoint with its current attributes
317
318 AppendPad(option);
319 }
320
321//______________________________________________________________________________
70a93198 322void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t)
ab48128d 323{
324 // Execute action corresponding to one event
325 // This member function is called when a AliEMCALRecPoint is clicked with the locator
326 //
327 // If Left button is clicked on AliEMCALRecPoint, the digits are switched on
328 // and switched off when the mouse button is released.
329
330 // static Int_t pxold, pyold;
331
70a93198 332 /* static TGraph * digitgraph = 0 ;
ab48128d 333 static TPaveText* clustertext = 0 ;
334
335 if (!gPad->IsEditable()) return;
336
337 switch (event) {
338
339
340 case kButton1Down:{
341 AliEMCALDigit * digit ;
88cb7938 342 AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry() ;
ab48128d 343
344 Int_t iDigit;
70a93198 345 Int_t relid[2] ;
ab48128d 346
347 const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ;
348 Float_t * xi = new Float_t [kMulDigit] ;
349 Float_t * zi = new Float_t [kMulDigit] ;
350
351 for(iDigit = 0; iDigit < kMulDigit; iDigit++) {
352 Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code");
353 digit = 0 ; //dynamic_cast<AliEMCALDigit *>((fDigitsList)[iDigit]);
354 emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
355 emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ;
356 }
357
358 if (!digitgraph) {
359 digitgraph = new TGraph(fMulDigit,xi,zi);
360 digitgraph-> SetMarkerStyle(5) ;
361 digitgraph-> SetMarkerSize(1.) ;
362 digitgraph-> SetMarkerColor(1) ;
363 digitgraph-> Draw("P") ;
364 }
365 if (!clustertext) {
366
367 TVector3 pos(0.,0.,0.) ;
368 GetLocalPosition(pos) ;
369 clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ;
370 Text_t line1[40] ;
371 Text_t line2[40] ;
372 sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ;
373 sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ;
374 clustertext ->AddText(line1) ;
375 clustertext ->AddText(line2) ;
376 clustertext ->Draw("");
377 }
378 gPad->Update() ;
9e5d2067 379 Print("") ;
ab48128d 380 delete[] xi ;
381 delete[] zi ;
382 }
383
384break;
385
386 case kButton1Up:
387 if (digitgraph) {
388 delete digitgraph ;
389 digitgraph = 0 ;
390 }
391 if (clustertext) {
392 delete clustertext ;
393 clustertext = 0 ;
394 }
395
396 break;
397
70a93198 398 }*/
399}
400//____________________________________________________________________________
401void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
402{
403 // Evaluates all shower parameters
70a93198 404 EvalLocalPosition(logWeight, digits) ;
405 EvalElipsAxis(logWeight, digits) ;
406 EvalDispersion(logWeight, digits) ;
4635df1f 407 //EvalCoreEnergy(logWeight, digits);
70a93198 408 EvalTime(digits) ;
87cdc3be 409 EvalPrimaries(digits) ;
410 EvalParents(digits);
70a93198 411}
412
413//____________________________________________________________________________
414void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
415{
416 // Calculates the dispersion of the shower at the origin of the RecPoint
1d46d1f6 417 // in cell units - Nov 16,2006
70a93198 418
1d46d1f6 419 Double_t d = 0., wtot = 0., w = 0.;
af5bdd85 420 Int_t iDigit=0, nstat=0;
70a93198 421 AliEMCALDigit * digit ;
70a93198 422
1d46d1f6 423 // Calculates the dispersion in cell units
424 Double_t etai, phii, etaMean=0.0, phiMean=0.0;
2bb3725c 425 int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
1d46d1f6 426 int iphi=0, ieta=0;
427 // Calculate mean values
70a93198 428 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
429 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
1963b290 430
1d46d1f6 431 if (fAmp>0 && fEnergyList[iDigit]>0) {
2bb3725c 432 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
433 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 434 etai=(Double_t)ieta;
435 phii=(Double_t)iphi;
436 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
437
438 if(w>0.0) {
439 phiMean += phii*w;
440 etaMean += etai*w;
441 wtot += w;
442 }
443 }
444 }
445 if (wtot>0) {
446 phiMean /= wtot ;
447 etaMean /= wtot ;
448 } else AliError(Form("Wrong weight %f\n", wtot));
70a93198 449
1d46d1f6 450 // Calculate dispersion
451 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
452 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
453
454 if (fAmp>0 && fEnergyList[iDigit]>0) {
2bb3725c 455 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
456 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 457 etai=(Double_t)ieta;
458 phii=(Double_t)iphi;
459 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
460
461 if(w>0.0) {
462 nstat++;
463 d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
e52475ed 464 }
465 }
ab48128d 466 }
70a93198 467
e52475ed 468 if ( wtot > 0 && nstat>1) d /= wtot ;
469 else d = 0. ;
70a93198 470
471 fDispersion = TMath::Sqrt(d) ;
ab48128d 472}
70a93198 473
ab48128d 474//____________________________________________________________________________
70a93198 475void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
88cb7938 476{
70a93198 477 // Calculates the center of gravity in the local EMCAL-module coordinates
e52475ed 478 // Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing
70a93198 479
1ae500a2 480 static Double_t dist;
481
e52475ed 482 AliEMCALDigit * digit;
1ae500a2 483 Int_t i=0, nstat=0, idMax=-1;
e52475ed 484 Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
70a93198 485
1ae500a2 486 //printf(" dist : %f e : %f \n", dist, fAmp);
e52475ed 487 for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
70a93198 488 digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
1ae500a2 489 if(iDigit==0) {
490 idMax = digit->GetId(); // is it correct
491 dist = TmaxInCm(Double_t(fAmp));
492 }
493 fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]);
494 //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), dist, xyzi[0], xyzi[1], xyzi[2]);
70a93198 495
1ae500a2 496 //fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
497 //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), 0.0, xyzi[0], xyzi[1], xyzi[2]);
498 // if(fAmp>102.) assert(0);
e52475ed 499
af5bdd85 500 if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
501 else w = fEnergyList[iDigit]; // just energy
ab48128d 502
e52475ed 503 if(w>0.0) {
504 wtot += w ;
505 nstat++;
506 for(i=0; i<3; i++ ) {
507 clXYZ[i] += (w*xyzi[i]);
508 clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
509 }
510 }
511 }
512 // cout << " wtot " << wtot << endl;
70a93198 513 if ( wtot > 0 ) {
e52475ed 514 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
515 for(i=0; i<3; i++ ) {
516 clXYZ[i] /= wtot;
517 if(nstat>1) {
518 clRmsXYZ[i] /= (wtot*wtot);
519 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
520 if(clRmsXYZ[i] > 0.0) {
521 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
522 } else clRmsXYZ[i] = 0;
523 } else clRmsXYZ[i] = 0;
524 }
70a93198 525 } else {
e52475ed 526 for(i=0; i<3; i++ ) {
527 clXYZ[i] = clRmsXYZ[i] = -1.;
528 }
70a93198 529 }
e52475ed 530 // clRmsXYZ[i] ??
531 fLocPos.SetX(clXYZ[0]);
532 fLocPos.SetY(clXYZ[1]);
533 fLocPos.SetZ(clXYZ[2]);
70a93198 534
1963b290 535// if (gDebug==2)
500aeccc 536// printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ;
e52475ed 537 fLocPosM = 0 ; // covariance matrix
ab48128d 538}
539
e52475ed 540//void AliEMCALRecPoint::EvalLocalPositionSimple()
541//{ // Weight is proportional of cell energy
542//}
1ae500a2 543//____________________________________________________________________________
544void AliEMCALRecPoint::EvalLocalPositionFit(Double_t deff, Double_t logWeight,
545Double_t phiSlope, TClonesArray * digits)
546{
547 // Aug 14-16, 2007 - for fit
548 // Aug 31 - should be static ??
549 static Double_t dist, ycorr;
550 static AliEMCALDigit *digit;
551
552 Int_t i=0, nstat=0, idMax=-1;
553 Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
554
555 for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
556 digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
557 if(iDigit==0) {
558 idMax = digit->GetId(); // is it correct
559 dist = TmaxInCm(Double_t(fAmp));
560 }
561
562 dist = deff;
563 fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]);
564
565 if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
566 else w = fEnergyList[iDigit]; // just energy
567
568 if(w>0.0) {
569 wtot += w ;
570 nstat++;
571 for(i=0; i<3; i++ ) {
572 clXYZ[i] += (w*xyzi[i]);
573 clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
574 }
575 }
576 }
577 // cout << " wtot " << wtot << endl;
578 if ( wtot > 0 ) {
579 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
580 for(i=0; i<3; i++ ) {
581 clXYZ[i] /= wtot;
582 if(nstat>1) {
583 clRmsXYZ[i] /= (wtot*wtot);
584 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
585 if(clRmsXYZ[i] > 0.0) {
586 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
587 } else clRmsXYZ[i] = 0;
588 } else clRmsXYZ[i] = 0;
589 }
590 } else {
591 for(i=0; i<3; i++ ) {
592 clXYZ[i] = clRmsXYZ[i] = -1.;
593 }
594 }
595 // clRmsXYZ[i] ??
596 if(phiSlope != 0.0 && logWeight > 0.0 && wtot) {
597 // Correction in phi direction (y - coords here); Aug 16;
598 // May be put to global level or seperate method
599 ycorr = clXYZ[1] * (1. + phiSlope);
600 //printf(" y %f : ycorr %f : slope %f \n", clXYZ[1], ycorr, phiSlope);
601 clXYZ[1] = ycorr;
602 }
603 fLocPos.SetX(clXYZ[0]);
604 fLocPos.SetY(clXYZ[1]);
605 fLocPos.SetZ(clXYZ[2]);
606
607// if (gDebug==2)
608// printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ;
609 fLocPosM = 0 ; // covariance matrix
610}
611
612Bool_t AliEMCALRecPoint::EvalLocalPosition2(TClonesArray * digits, TArrayD &ed)
613{
614 // Evaluated local position of rec.point using digits
615 // and parametrisation of w0 and deff
616 //printf(" <I> AliEMCALRecPoint::EvalLocalPosition2() \n");
617 return AliEMCALRecPoint::EvalLocalPositionFromDigits(digits, ed, fLocPos);
618}
619
620Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(TClonesArray *digits, TArrayD &ed, TVector3 &locPos)
621{
622 // Used when digits should be recalibrated
623 static Double_t deff, w0, esum;
624 static Int_t iDigit;
625 static AliEMCALDigit *digit;
626
627 if(ed.GetSize() && (digits->GetEntries()!=ed.GetSize())) return kFALSE;
628
629 // Calculate sum energy of digits
630 esum = 0.0;
631 for(iDigit=0; iDigit<ed.GetSize(); iDigit++) esum += ed[iDigit];
632
633 GetDeffW0(esum, deff, w0);
634
635 return EvalLocalPositionFromDigits(esum, deff, w0, digits, ed, locPos);
636}
637
638Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(const Double_t esum, const Double_t deff, const Double_t w0,
639 TClonesArray *digits, TArrayD &ed, TVector3 &locPos)
640{
641 static AliEMCALDigit *digit;
642
643 Int_t i=0, nstat=0, idMax=-1;
644 Double_t clXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
645
646 AliEMCALGeometry* geo = AliEMCALGeometry::GetInstance(); // Get pointer to EMCAL geometry
647
648 for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
649 digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit));
650
651 geo->RelPosCellInSModule(digit->GetId(), idMax, deff, xyzi[0], xyzi[1], xyzi[2]);
652
653 if(w0 > 0.0) w = TMath::Max( 0., w0 + TMath::Log(ed[iDigit] / esum));
654 else w = ed[iDigit]; // just energy
655
656 if(w>0.0) {
657 wtot += w ;
658 nstat++;
659 for(i=0; i<3; i++ ) {
660 clXYZ[i] += (w*xyzi[i]);
661 }
662 }
663 }
664 // cout << " wtot " << wtot << endl;
665 if (wtot > 0) {
666 for(i=0; i<3; i++ ) {
667 clXYZ[i] /= wtot;
668 }
669 locPos.SetX(clXYZ[0]);
670 locPos.SetY(clXYZ[1]);
671 locPos.SetZ(clXYZ[2]);
672 return kTRUE;
673 } else {
674 return kFALSE;
675 }
676
677}
678
679void AliEMCALRecPoint::GetDeffW0(const Double_t esum , Double_t &deff, Double_t &w0)
680{
681 //
682 // Aug 31, 2001
683 // Applied for simulation data with threshold 3 adc
684 // Calculate efective distance (deff) and weigh parameter (w0)
685 // for coordinate calculation; 0.5 GeV < esum <100 GeV.
686 // Look to: http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/CALIB/GEOMCORR/deffandW0VaEgamma_2.gif
687 //
688 static Double_t e=0.0;
689 const Double_t dp0=9.25147, dp1=1.16700; // Hard coded now
690 const Double_t wp0=4.83713, wp1=-2.77970e-01, wp2 = 4.41116;
691
692 // No extrapolation here
693 e = esum<0.5?0.5:esum;
694 e = e>100.?100.:e;
695
696 deff = dp0 + dp1*TMath::Log(e);
697 w0 = wp0 / (1. + TMath::Exp(wp1*(e+wp2)));
698 //printf("<I> AliEMCALRecPoint::GetDeffW0 esum %5.2f : deff %5.2f : w0 %5.2f \n", esum, deff, w0);
699}
e52475ed 700
70a93198 701//______________________________________________________________________________
702void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
703{
704 // This function calculates energy in the core,
4635df1f 705 // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
70a93198 706 // in accordance with shower profile the energy deposition
707 // should be less than 2%
1d46d1f6 708 // Unfinished - Nov 15,2006
709 // Distance is calculate in (phi,eta) units
70a93198 710
711 AliEMCALDigit * digit ;
5dee926e 712
70a93198 713 Int_t iDigit;
714
e52475ed 715 if (!fLocPos.Mag()) {
70a93198 716 EvalLocalPosition(logWeight, digits);
717 }
718
1d46d1f6 719 Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta();
720 Double_t eta, phi, distance;
70a93198 721 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
722 digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
4635df1f 723
1d46d1f6 724 eta = phi = 0.0;
725 fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ;
726 phi = phi * TMath::DegToRad();
70a93198 727
1d46d1f6 728 distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint));
70a93198 729 if(distance < fCoreRadius)
730 fCoreEnergy += fEnergyList[iDigit] ;
731 }
732
733}
ab48128d 734//____________________________________________________________________________
70a93198 735void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
ab48128d 736{
70a93198 737 // Calculates the axis of the shower ellipsoid in eta and phi
1d46d1f6 738 // in cell units
ab48128d 739
1d46d1f6 740 static TString gn(fGeomPtr->GetName());
741
742 Double_t wtot = 0.;
70a93198 743 Double_t x = 0.;
744 Double_t z = 0.;
745 Double_t dxx = 0.;
746 Double_t dzz = 0.;
747 Double_t dxz = 0.;
ab48128d 748
1d46d1f6 749 AliEMCALDigit * digit = 0;
70a93198 750
1d46d1f6 751 Double_t etai , phii, w;
2bb3725c 752 int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
1d46d1f6 753 int iphi=0, ieta=0;
754 for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
70a93198 755 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
1d46d1f6 756 etai = phii = 0.;
757 if(gn.Contains("SHISH")) {
758 // Nov 15,2006 - use cell numbers as coordinates
759 // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi
760 // We can use the eta,phi(or coordinates) of cell
2bb3725c 761 nSupMod = nModule = nIphi = nIeta = iphi = ieta = 0;
1d46d1f6 762
2bb3725c 763 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
764 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 765 etai=(Double_t)ieta;
766 phii=(Double_t)iphi;
767 } else { //
14ce0a6e 768 fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
1d46d1f6 769 phii = phii * TMath::DegToRad();
1963b290 770 }
771
1d46d1f6 772 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
773 // fAmp summed amplitude of digits, i.e. energy of recpoint
774 // Gives smaller value of lambda than log weight
775 // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy
ff1e7e2f 776
70a93198 777 dxx += w * etai * etai ;
778 x += w * etai ;
779 dzz += w * phii * phii ;
780 z += w * phii ;
1963b290 781
ff1e7e2f 782 dxz += w * etai * phii ;
1963b290 783
70a93198 784 wtot += w ;
785 }
ff1e7e2f 786
70a93198 787 if ( wtot > 0 ) {
788 dxx /= wtot ;
789 x /= wtot ;
790 dxx -= x * x ;
791 dzz /= wtot ;
792 z /= wtot ;
793 dzz -= z * z ;
794 dxz /= wtot ;
795 dxz -= x * z ;
ab48128d 796
70a93198 797 fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
798 if(fLambda[0] > 0)
799 fLambda[0] = TMath::Sqrt(fLambda[0]) ;
800 else
801 fLambda[0] = 0;
802
803 fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
ff1e7e2f 804
70a93198 805 if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
806 fLambda[1] = TMath::Sqrt(fLambda[1]) ;
807 else
808 fLambda[1]= 0. ;
809 } else {
810 fLambda[0]= 0. ;
811 fLambda[1]= 0. ;
ab48128d 812 }
ff1e7e2f 813
1963b290 814 // printf("Evalaxis: lambdas = %f,%f", fLambda[0],fLambda[1]) ;
ff1e7e2f 815
ab48128d 816}
817
818//______________________________________________________________________________
819void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
820{
af5bdd85 821 // Constructs the list of primary particles (tracks) which
822 // have contributed to this RecPoint and calculate deposited energy
823 // for each track
ab48128d 824
825 AliEMCALDigit * digit ;
af5bdd85 826 Int_t * primArray = new Int_t[fMaxTrack] ;
827 Float_t * dEPrimArray = new Float_t[fMaxTrack] ;
ab48128d 828
829 Int_t index ;
830 for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
831 digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
832 Int_t nprimaries = digit->GetNprimary() ;
5c0368b8 833 if ( nprimaries == 0 ) continue ;
ab48128d 834 Int_t jndex ;
835 for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
836 if ( fMulTrack > fMaxTrack ) {
f792c312 837 fMulTrack = fMaxTrack ;
94478418 838 Error("EvalPrimaries", "increase fMaxTrack ") ;
ab48128d 839 break ;
840 }
af5bdd85 841 Int_t newPrimary = digit->GetPrimary(jndex+1);
842 Float_t dEPrimary = digit->GetDEPrimary(jndex+1);
ab48128d 843 Int_t kndex ;
844 Bool_t already = kFALSE ;
845 for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
af5bdd85 846 if ( newPrimary == primArray[kndex] ){
ab48128d 847 already = kTRUE ;
af5bdd85 848 dEPrimArray[kndex] += dEPrimary;
ab48128d 849 break ;
850 }
851 } // end of check
5c0368b8 852 if ( !already && (fMulTrack < fMaxTrack)) { // store it
af5bdd85 853 primArray[fMulTrack] = newPrimary ;
854 dEPrimArray[fMulTrack] = dEPrimary ;
ab48128d 855 fMulTrack++ ;
856 } // store it
857 } // all primaries in digit
ab48128d 858 } // all digits
859
af5bdd85 860 Int_t *sortIdx = new Int_t[fMulTrack];
861 TMath::Sort(fMulTrack,dEPrimArray,sortIdx);
862 for(index = 0; index < fMulTrack; index++) {
863 fTracksList[index] = primArray[sortIdx[index]] ;
864 fDETracksList[index] = dEPrimArray[sortIdx[index]] ;
865 }
866 delete [] sortIdx;
867 delete [] primArray ;
868 delete [] dEPrimArray ;
ab48128d 869
870}
7ee5c5be 871
87cdc3be 872//______________________________________________________________________________
873void AliEMCALRecPoint::EvalParents(TClonesArray * digits)
874{
875 // Constructs the list of parent particles (tracks) which have contributed to this RecPoint
876
877 AliEMCALDigit * digit ;
af5bdd85 878 Int_t * parentArray = new Int_t[fMaxTrack] ;
879 Float_t * dEParentArray = new Float_t[fMaxTrack] ;
87cdc3be 880
881 Int_t index ;
882 for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
94478418 883 if (fDigitsList[index] >= digits->GetEntries() || fDigitsList[index] < 0)
884 AliError(Form("Trying to get invalid digit %d (idx in WriteRecPoint %d)",fDigitsList[index],index));
87cdc3be 885 digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
886 Int_t nparents = digit->GetNiparent() ;
5c0368b8 887 if ( nparents == 0 ) continue ;
87cdc3be 888
889 Int_t jndex ;
890 for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit
891 if ( fMulParent > fMaxParent ) {
892 fMulTrack = - 1 ;
94478418 893 Error("EvalParents", "increase fMaxParent") ;
87cdc3be 894 break ;
895 }
af5bdd85 896 Int_t newParent = digit->GetIparent(jndex+1) ;
897 Float_t newdEParent = digit->GetDEParent(jndex+1) ;
87cdc3be 898 Int_t kndex ;
899 Bool_t already = kFALSE ;
f1d429fd 900 for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored
af5bdd85 901 if ( newParent == parentArray[kndex] ){
902 dEParentArray[kndex] += newdEParent;
87cdc3be 903 already = kTRUE ;
904 break ;
905 }
906 } // end of check
94478418 907 if ( !already && (fMulParent < fMaxParent)) { // store it
af5bdd85 908 parentArray[fMulParent] = newParent ;
909 dEParentArray[fMulParent] = newdEParent ;
87cdc3be 910 fMulParent++ ;
911 } // store it
912 } // all parents in digit
87cdc3be 913 } // all digits
914
27e2a47c 915 if (fMulParent>0) {
af5bdd85 916 Int_t *sortIdx = new Int_t[fMulParent];
917 TMath::Sort(fMulParent,dEParentArray,sortIdx);
918 for(index = 0; index < fMulParent; index++) {
919 fParentsList[index] = parentArray[sortIdx[index]] ;
920 fDEParentsList[index] = dEParentArray[sortIdx[index]] ;
921 }
922 delete [] sortIdx;
27e2a47c 923 }
87cdc3be 924
af5bdd85 925 delete [] parentArray;
926 delete [] dEParentArray;
87cdc3be 927}
928
70a93198 929//____________________________________________________________________________
930void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const
931{
932 // returns the position of the cluster in the local reference system of ALICE
70a93198 933
934 lpos.SetX(fLocPos.X()) ;
935 lpos.SetY(fLocPos.Y()) ;
936 lpos.SetZ(fLocPos.Z()) ;
937}
938
ab48128d 939//____________________________________________________________________________
940void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
941{
942 // returns the position of the cluster in the global reference system of ALICE
70a93198 943 // These are now the Cartesian X, Y and Z
e52475ed 944 // cout<<" geom "<<geom<<endl;
14ce0a6e 945 fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
70a93198 946}
947
948//____________________________________________________________________________
949Float_t AliEMCALRecPoint::GetMaximalEnergy(void) const
950{
951 // Finds the maximum energy in the cluster
ab48128d 952
70a93198 953 Float_t menergy = 0. ;
954
955 Int_t iDigit;
956
957 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
958
959 if(fEnergyList[iDigit] > menergy)
960 menergy = fEnergyList[iDigit] ;
961 }
962 return menergy ;
ab48128d 963}
964
aad8e277 965//____________________________________________________________________________
70a93198 966Int_t AliEMCALRecPoint::GetMultiplicityAtLevel(Float_t H) const
aad8e277 967{
70a93198 968 // Calculates the multiplicity of digits with energy larger than H*energy
969
970 Int_t multipl = 0 ;
971 Int_t iDigit ;
972 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
973
974 if(fEnergyList[iDigit] > H * fAmp)
975 multipl++ ;
976 }
977 return multipl ;
978}
979
980//____________________________________________________________________________
981Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
982 Float_t locMaxCut,TClonesArray * digits) const
983{
984 // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
985 // energy difference between two local maxima
986
987 AliEMCALDigit * digit ;
988 AliEMCALDigit * digitN ;
989
990 Int_t iDigitN ;
991 Int_t iDigit ;
992
993 for(iDigit = 0; iDigit < fMulDigit; iDigit++)
994 maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit]) ;
995
996 for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
997 if(maxAt[iDigit]) {
998 digit = maxAt[iDigit] ;
999
1000 for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
1001 digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ;
1002
1003 if ( AreNeighbours(digit, digitN) ) {
1004 if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
1005 maxAt[iDigitN] = 0 ;
1006 // but may be digit too is not local max ?
1007 if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
1008 maxAt[iDigit] = 0 ;
1009 }
1010 else {
1011 maxAt[iDigit] = 0 ;
1012 // but may be digitN too is not local max ?
1013 if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
1014 maxAt[iDigitN] = 0 ;
1015 }
1016 } // if Areneighbours
1017 } // while digitN
1018 } // slot not empty
1019 } // while digit
1020
1021 iDigitN = 0 ;
1022 for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
1023 if(maxAt[iDigit] ){
1024 maxAt[iDigitN] = maxAt[iDigit] ;
1025 maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
1026 iDigitN++ ;
1027 }
1028 }
1029 return iDigitN ;
1030}
4635df1f 1031
1032//____________________________________________________________________________
1033Int_t AliEMCALRecPoint::GetPrimaryIndex() const
1034{
1035 // Get the primary track index in TreeK which deposits the most energy
af5bdd85 1036 // in Digits which forms RecPoint.
4635df1f 1037
af5bdd85 1038 if (fMulTrack)
1039 return fTracksList[0];
1040 return -12345;
4635df1f 1041}
1042
70a93198 1043//____________________________________________________________________________
1044void AliEMCALRecPoint::EvalTime(TClonesArray * digits){
1045 // time is set to the time of the digit with the maximum energy
1046
1047 Float_t maxE = 0;
1048 Int_t maxAt = 0;
1049 for(Int_t idig=0; idig < fMulDigit; idig++){
1050 if(fEnergyList[idig] > maxE){
1051 maxE = fEnergyList[idig] ;
1052 maxAt = idig;
1053 }
1054 }
1055 fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
aad8e277 1056
aad8e277 1057}
ab48128d 1058
1059//______________________________________________________________________________
1060void AliEMCALRecPoint::Paint(Option_t *)
1061{
1062 // Paint this ALiRecPoint as a TMarker with its current attributes
1063
1064 TVector3 pos(0.,0.,0.) ;
1065 GetLocalPosition(pos) ;
1066 Coord_t x = pos.X() ;
1067 Coord_t y = pos.Z() ;
1068 Color_t markercolor = 1 ;
1069 Size_t markersize = 1. ;
1070 Style_t markerstyle = 5 ;
1071
1072 if (!gPad->IsBatch()) {
1073 gVirtualX->SetMarkerColor(markercolor) ;
1074 gVirtualX->SetMarkerSize (markersize) ;
1075 gVirtualX->SetMarkerStyle(markerstyle) ;
1076 }
1077 gPad->SetAttMarkerPS(markercolor,markerstyle,markersize) ;
1078 gPad->PaintPolyMarker(1,&x,&y,"") ;
1079}
70a93198 1080
1ae500a2 1081Double_t AliEMCALRecPoint::TmaxInCm(const Double_t e , const Int_t key)
1082{
1083 // e energ in in GeV)
1084 // key = 0(gamma, default)
1085 // != 0(electron)
1086 static Double_t ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
1087 static Double_t X0 = 1.23; // radiation lenght (cm)
1088 static Double_t tmax = 0.; // position of electromagnetic shower max in cm
1089
1090 tmax = 0.0;
1091 if(e>0.1) {
1092 tmax = TMath::Log(e) + ca;
1093 if (key==0) tmax += 0.5;
1094 else tmax -= 0.5;
1095 tmax *= X0; // convert to cm
1096 }
1097 return tmax;
1098}
1099
70a93198 1100//______________________________________________________________________________
1101Float_t AliEMCALRecPoint::EtaToTheta(Float_t arg) const
1102{
1103 //Converts Theta (Radians) to Eta(Radians)
1104 return (2.*TMath::ATan(TMath::Exp(-arg)));
1105}
1106
1107//______________________________________________________________________________
1108Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const
1109{
1110 //Converts Eta (Radians) to Theta(Radians)
1111 return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
1112}
261b1065 1113
1114//____________________________________________________________________________
e1a51e6e 1115void AliEMCALRecPoint::Print(Option_t *opt) const
261b1065 1116{
1117 // Print the list of digits belonging to the cluster
e1a51e6e 1118 if(strlen(opt)==0) return;
261b1065 1119 TString message ;
4800667c 1120 message = "AliEMCALRecPoint:\n" ;
261b1065 1121 message += " digits # = " ;
1122 Info("Print", message.Data()) ;
1123
1124 Int_t iDigit;
1125 for(iDigit=0; iDigit<fMulDigit; iDigit++)
1126 printf(" %d ", fDigitsList[iDigit] ) ;
e52475ed 1127 printf("\n");
1128
261b1065 1129 Info("Print", " Energies = ") ;
1130 for(iDigit=0; iDigit<fMulDigit; iDigit++)
1131 printf(" %f ", fEnergyList[iDigit] ) ;
e52475ed 1132 printf("\n");
1133
1134 Info("Print", "\n Abs Ids = ") ;
1135 for(iDigit=0; iDigit<fMulDigit; iDigit++)
1136 printf(" %i ", fAbsIdList[iDigit] ) ;
1137 printf("\n");
1138
1139 Info("Print", " Primaries ") ;
261b1065 1140 for(iDigit = 0;iDigit < fMulTrack; iDigit++)
1141 printf(" %d ", fTracksList[iDigit]) ;
e52475ed 1142
1143 printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]);
1144
85c60a8e 1145 message = " ClusterType = %d" ;
1146 message += " Multiplicity = %d" ;
261b1065 1147 message += " Cluster Energy = %f" ;
1148 message += " Core energy = %f" ;
1149 message += " Core radius = %f" ;
1150 message += " Number of primaries %d" ;
1151 message += " Stored at position %d" ;
85c60a8e 1152 Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;
261b1065 1153}
1d46d1f6 1154
1155Double_t AliEMCALRecPoint::GetPointEnergy() const
1156{
1157 static double e;
1158 e=0.0;
1159 for(int ic=0; ic<GetMultiplicity(); ic++) e += double(fEnergyList[ic]);
1160 return e;
1161}