]> 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
Moved coordinates
[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
126 Int_t relid1[4] ;
127 phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
128
129 Int_t relid2[4] ;
130 phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
131
132 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
133 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
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;
210// Int_t relid[4] ;
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 ;
314
315 if (IsInPRE())
316 cyl_radius = emcalgeom->GetIP2PRESection() ;
88cb7938 317 else if (IsInECA())
318 cyl_radius = emcalgeom->GetIP2ECASection() ;
319 else if (IsInHCA())
320 cyl_radius = emcalgeom->GetIP2HCASection() ;
12c037a8 321 else
322 Fatal("EvalDispersion", "Unexpected tower section!") ;
ab48128d 323
12c037a8 324 Float_t x = fLocPos.X() ;
325 Float_t y = fLocPos.Y() ;
326 Float_t z = fLocPos.Z() ;
ab48128d 327
12c037a8 328 if (gDebug == 2)
329 Info("EvalDispersion", "x,y,z = %f,%f,%f", x, y, z) ;
330
ab48128d 331// Calculates the dispersion in coordinates
332 wtot = 0.;
333 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
334 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
335 Float_t thetai = 0. ;
336 Float_t phii = 0.;
12c037a8 337 emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
ab48128d 338 Float_t xi = cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
339 Float_t yi = cyl_radius * TMath::Sin(phii * kDeg2Rad ) ;
91ca893f 340 Float_t zi = cyl_radius / TMath::Tan(thetai * kDeg2Rad ) ;
ab48128d 341
12c037a8 342 if (gDebug == 2)
343 Info("EvalDispersion", "id = %d, xi,yi,zi = %f,%f,%f", digit->GetId(), xi, yi, zi) ;
344
ab48128d 345 Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
12c037a8 346 d += w * ( (xi-x)*(xi-x) + (zi-z)*(zi-z) ) ;
ab48128d 347 wtot+=w ;
348 }
349
e7f14e3c 350 if ( wtot > 0 )
351 d /= wtot ;
352 else
353 d = 0. ;
ab48128d 354
355 fDispersion = TMath::Sqrt(d) ;
356
357}
358//______________________________________________________________________________
359void AliEMCALTowerRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
360{
361 // This function calculates energy in the core,
362 // i.e. within a radius rad = 3cm around the center. Beyond this radius
363 // in accordance with shower profile the energy deposition
364 // should be less than 2%
365
366 Float_t coreRadius = 10. ;
367
368 AliEMCALDigit * digit ;
ab48128d 369 Float_t wtot = 0. ;
370
88cb7938 371 AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry();
ab48128d 372 Int_t iDigit;
373
374 if (!fTheta || !fPhi ) {
375 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
376 digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
377
378 Float_t thetai ;
379 Float_t phii ;
12c037a8 380 emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
ab48128d 381 Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
382 fTheta = fTheta + thetai * w ;
383 fPhi += (phii * w );
384 wtot += w ;
385 }
386
e7f14e3c 387 if (wtot > 0 ) {
388 fTheta /= wtot ;
389 fPhi /= wtot ;
390 } else {
391 fTheta = -1 ;
392 fPhi = -1 ;
393 }
ab48128d 394 }
12c037a8 395
396 const Float_t kDeg2Rad = TMath::DegToRad() ;
397
88cb7938 398 Float_t cyl_radius = emcalgeom->GetIP2ECASection();
ab48128d 399 Float_t x = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
400 Float_t y = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
401 Float_t z = cyl_radius * TMath::Tan(fTheta * kDeg2Rad ) ;
402
403 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
404 digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
ab48128d 405 Float_t thetai = 0. ;
406 Float_t phii = 0. ;
12c037a8 407 emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
ab48128d 408
409 Float_t xi = cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
410 Float_t yi = cyl_radius * TMath::Sin(phii * kDeg2Rad ) ;
411 Float_t zi = cyl_radius * TMath::Tan(thetai * kDeg2Rad ) ;
412
413 Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(yi-y)*(yi-y)+(zi-z)*(zi-z)) ;
414 if(distance < coreRadius)
415 fCoreEnergy += fEnergyList[iDigit] ;
416 }
12c037a8 417
ab48128d 418}
419
420//____________________________________________________________________________
421void AliEMCALTowerRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
422{
423 // Calculates the axis of the shower ellipsoid
424
425 Double_t wtot = 0. ;
426 Double_t x = 0.;
427 Double_t z = 0.;
428 Double_t dxx = 0.;
429 Double_t dzz = 0.;
430 Double_t dxz = 0.;
431
432 AliEMCALDigit * digit ;
433
88cb7938 434 AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry();
ab48128d 435
436 Int_t iDigit;
12c037a8 437 const Float_t kDeg2Rad = TMath::DegToRad() ;
438
439 Float_t cyl_radius = 0 ;
ab48128d 440
12c037a8 441 if (IsInPRE())
442 cyl_radius = emcalgeom->GetIP2PRESection() ;
88cb7938 443 else if (IsInECA())
444 cyl_radius = emcalgeom->GetIP2ECASection() ;
445 else if (IsInHCA())
446 cyl_radius = emcalgeom->GetIP2HCASection() ;
12c037a8 447 else
448 Fatal("EvalDispersion", "Unexpected tower section!") ;
ab48128d 449
450 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
451 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
ab48128d 452 Float_t thetai = 0. ;
453 Float_t phii = 0. ;
12c037a8 454 emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
ab48128d 455 Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
456 Float_t xi = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
12c037a8 457 Float_t zi = cyl_radius / TMath::Tan(fTheta * kDeg2Rad ) ;
ab48128d 458 dxx += w * xi * xi ;
459 x += w * xi ;
460 dzz += w * zi * zi ;
461 z += w * zi ;
462 dxz += w * xi * zi ;
463 wtot += w ;
464 }
e7f14e3c 465 if ( wtot > 0 ) {
466 dxx /= wtot ;
467 x /= wtot ;
468 dxx -= x * x ;
469 dzz /= wtot ;
470 z /= wtot ;
471 dzz -= z * z ;
472 dxz /= wtot ;
473 dxz -= x * z ;
474
475
476 // //Apply correction due to non-perpendicular incidence
ab48128d 477// Double_t CosX ;
478// Double_t CosZ ;
88cb7938 479// AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry();
ab48128d 480 // Double_t DistanceToIP= (Double_t ) emcalgeom->GetIPDistance() ;
481
482// CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
483// CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
484
485// dxx = dxx/(CosX*CosX) ;
486// dzz = dzz/(CosZ*CosZ) ;
487// dxz = dxz/(CosX*CosZ) ;
488
489
e7f14e3c 490 fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
491 if(fLambda[0] > 0)
492 fLambda[0] = TMath::Sqrt(fLambda[0]) ;
493
494 fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
495 if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
496 fLambda[1] = TMath::Sqrt(fLambda[1]) ;
497 else
498 fLambda[1]= 0. ;
499 } else {
500 fLambda[0]= 0. ;
ab48128d 501 fLambda[1]= 0. ;
e7f14e3c 502 }
ab48128d 503}
504
505//____________________________________________________________________________
506void AliEMCALTowerRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits )
507{
508 // Evaluates all shower parameters
509
510 AliEMCALRecPoint::EvalAll(logWeight,digits) ;
511 EvalGlobalPosition(logWeight, digits) ;
512 EvalElipsAxis(logWeight, digits) ;
513 EvalDispersion(logWeight, digits) ;
514 EvalCoreEnergy(logWeight, digits);
515 EvalTime(digits) ;
516}
517
518//____________________________________________________________________________
519void AliEMCALTowerRecPoint::EvalGlobalPosition(Float_t logWeight, TClonesArray * digits)
520{
521 // Calculates the center of gravity in the local EMCAL-module coordinates
522 Float_t wtot = 0. ;
523
12c037a8 524 // Int_t relid[4] ;
ab48128d 525
526 AliEMCALDigit * digit ;
88cb7938 527 AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry();
ab48128d 528 Int_t iDigit;
529
530 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
531 digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
532
533 Float_t thetai ;
534 Float_t phii ;
12c037a8 535 emcalgeom->PosInAlice(digit->GetId(), thetai, phii);
ab48128d 536 Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
537 fTheta = fTheta + thetai * w ;
538 fPhi += (phii * w );
539 wtot += w ;
540 }
541
e7f14e3c 542 if ( wtot > 0 ) {
543 fTheta /= wtot ;
544 fPhi /= wtot ;
545 } else {
546 fTheta = -1 ;
547 fPhi = -1.;
548 }
549
12c037a8 550
551 const Float_t kDeg2Rad = TMath::DegToRad() ;
552
553 Float_t cyl_radius = 0 ;
554
555 if (IsInPRE())
556 cyl_radius = emcalgeom->GetIP2PRESection() ;
88cb7938 557 else if (IsInECA())
558 cyl_radius = emcalgeom->GetIP2ECASection() ;
559 else if (IsInHCA())
560 cyl_radius = emcalgeom->GetIP2HCASection() ;
12c037a8 561 else
562 Fatal("EvalGlobalPosition", "Unexpected tower section!") ;
563
564 Float_t x = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
565 Float_t y = cyl_radius * TMath::Sin(fPhi * kDeg2Rad ) ;
566 Float_t z = cyl_radius / TMath::Tan(fTheta * kDeg2Rad ) ;
567
568 fLocPos.SetX(x) ;
569 fLocPos.SetY(y) ;
570 fLocPos.SetZ(z) ;
571
572 if (gDebug==2)
573 Info("EvalGlobalPosition", "x,y,z = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ;
574
ab48128d 575
576 fLocPosM = 0 ;
577}
578
579//____________________________________________________________________________
580Float_t AliEMCALTowerRecPoint::GetMaximalEnergy(void) const
581{
582 // Finds the maximum energy in the cluster
583
584 Float_t menergy = 0. ;
585
586 Int_t iDigit;
587
588 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
589
590 if(fEnergyList[iDigit] > menergy)
591 menergy = fEnergyList[iDigit] ;
592 }
593 return menergy ;
594}
595
596//____________________________________________________________________________
597Int_t AliEMCALTowerRecPoint::GetMultiplicityAtLevel(const Float_t H) const
598{
599 // Calculates the multiplicity of digits with energy larger than H*energy
600
601 Int_t multipl = 0 ;
602 Int_t iDigit ;
603 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
604
605 if(fEnergyList[iDigit] > H * fAmp)
606 multipl++ ;
607 }
608 return multipl ;
609}
610
611//____________________________________________________________________________
a0636361 612Int_t AliEMCALTowerRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
ab48128d 613 Float_t locMaxCut,TClonesArray * digits) const
614{
615 // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
616 // energy difference between two local maxima
617
618 AliEMCALDigit * digit ;
619 AliEMCALDigit * digitN ;
620
621
622 Int_t iDigitN ;
623 Int_t iDigit ;
624
625 for(iDigit = 0; iDigit < fMulDigit; iDigit++)
a0636361 626 maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit]) ;
ab48128d 627
628
629 for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
a0636361 630 if(maxAt[iDigit]) {
631 digit = maxAt[iDigit] ;
ab48128d 632
633 for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
634 digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ;
635
636 if ( AreNeighbours(digit, digitN) ) {
637 if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
a0636361 638 maxAt[iDigitN] = 0 ;
ab48128d 639 // but may be digit too is not local max ?
640 if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
a0636361 641 maxAt[iDigit] = 0 ;
ab48128d 642 }
643 else {
a0636361 644 maxAt[iDigit] = 0 ;
ab48128d 645 // but may be digitN too is not local max ?
646 if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
a0636361 647 maxAt[iDigitN] = 0 ;
ab48128d 648 }
649 } // if Areneighbours
650 } // while digitN
651 } // slot not empty
652 } // while digit
653
654 iDigitN = 0 ;
655 for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
a0636361 656 if(maxAt[iDigit] ){
ab48128d 657 maxAt[iDigitN] = maxAt[iDigit] ;
658 maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
659 iDigitN++ ;
660 }
661 }
662 return iDigitN ;
663}
664//____________________________________________________________________________
665void AliEMCALTowerRecPoint::EvalTime(TClonesArray * digits){
666
667 Float_t maxE = 0;
668 Int_t maxAt = 0;
669 for(Int_t idig=0; idig < fMulDigit; idig++){
670 if(fEnergyList[idig] > maxE){
671 maxE = fEnergyList[idig] ;
672 maxAt = idig;
673 }
674 }
675 fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
676
677}
678//____________________________________________________________________________
9e5d2067 679void AliEMCALTowerRecPoint::Print(Option_t *)
ab48128d 680{
681 // Print the list of digits belonging to the cluster
682
9859bfc0 683 TString message("\n") ;
ab48128d 684
685 Int_t iDigit;
9859bfc0 686 message += "digits # = " ;
687 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
688 message += fDigitsList[iDigit] ;
689 message += " " ;
690 }
691
692 message += "\nEnergies = " ;
693 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
694 message += fEnergyList[iDigit] ;
695 message += " " ;
696 }
697
698 message += "\nPrimaries " ;
699 for(iDigit = 0;iDigit < fMulTrack; iDigit++) {
700 message += fTracksList[iDigit] ;
701 message += " " ;
702 }
703 message += "\n Multiplicity = " ;
704 message += fMulDigit ;
705 message += "\n Cluster Energy = " ;
706 message += fAmp ;
707 message += "\n Number of primaries " ;
708 message += fMulTrack ;
709 message += "\n Stored at position " ;
710 message += GetIndexInList() ;
711
712 Info("Print", message.Data() ) ;
ab48128d 713}
714
12c037a8 715//____________________________________________________________________________
716const TVector3 AliEMCALTowerRecPoint::XYZInAlice(Float_t r, Float_t theta, Float_t phi) const
717{
718 // spherical coordinates of recpoint in Alice reference frame
719
720 if (gDebug == 2)
721 Info("XYZInAlice", "this= %d , r = %f, theta = %f, phi = %f", this, r, theta, phi) ;
722
723 if (theta == 9999. || phi == 9999. || r == 9999.) {
724 TVector3 globalpos;
725 GetGlobalPosition(globalpos);
726 phi = globalpos.X() * TMath::DegToRad() ;
727 r = globalpos.Y() ;
728 theta = globalpos.Z() * TMath::DegToRad() ;
729 }
730 else {
731 theta *= TMath::DegToRad() ;
732 phi *= TMath::DegToRad() ;
733 }
734
735 Float_t y = r * TMath::Cos(phi) ;
736 Float_t x = r * TMath::Sin(phi) * TMath::Sin(theta) ;
737 Float_t z = r * TMath::Sin(phi) * TMath::Cos(theta) ;
738
739 TVector3 vec(z, x, y) ;
740 return vec ;
741}