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