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