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