Changes imposed by HP-UX which ignores the new C++ standart
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
CommitLineData
d15a28e7 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
b2a60966 16/* $Id$ */
17
d15a28e7 18//_________________________________________________________________________
b2a60966 19// RecPoint implementation for PHOS-EMC
20// An EmcRecPoint is a cluster of digits
21//
22//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
23
d15a28e7 24
25// --- ROOT system ---
9f616d61 26#include "TPad.h"
27#include "TH2.h"
d15a28e7 28#include "TMath.h"
9f616d61 29#include "TCanvas.h"
d15a28e7 30
31// --- Standard library ---
32
de9ec31b 33#include <iostream.h>
d15a28e7 34
35// --- AliRoot header files ---
36
37#include "AliPHOSGeometry.h"
38#include "AliPHOSEmcRecPoint.h"
39#include "AliRun.h"
40
41ClassImp(AliPHOSEmcRecPoint)
42
43//____________________________________________________________________________
44AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
45 : AliPHOSRecPoint()
46{
47 // ctor
48
49 fMulDigit = 0 ;
50 fAmp = 0. ;
51 fEnergyList = new Float_t[fMaxDigit];
92862013 52 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
53 fDelta = phosgeom->GetCrystalSize(0) ;
d15a28e7 54 fW0 = W0 ;
55 fLocMaxCut = LocMaxCut ;
56 fLocPos.SetX(1000000.) ; //Local position should be evaluated
57}
58
9f616d61 59//____________________________________________________________________________
d15a28e7 60void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
61{
b2a60966 62 // Adds a digit to the RecPoint
63 // and accumulates the total amplitude and the multiplicity
d15a28e7 64
65 if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
9f616d61 66 fMaxDigit*=2 ;
67 int * tempo = new ( int[fMaxDigit] ) ;
68 Float_t * tempoE = new ( Float_t[fMaxDigit] ) ;
69
70 Int_t index ;
d15a28e7 71 for ( index = 0 ; index < fMulDigit ; index++ ){
72 tempo[index] = fDigitsList[index] ;
73 tempoE[index] = fEnergyList[index] ;
74 }
75
9f616d61 76 delete [] fDigitsList ;
77 fDigitsList = new ( int[fMaxDigit] ) ;
78
79 delete [] fEnergyList ;
80 fEnergyList = new ( Float_t[fMaxDigit] ) ;
81
82 for ( index = 0 ; index < fMulDigit ; index++ ){
83 fDigitsList[index] = tempo[index] ;
84 fEnergyList[index] = tempoE[index] ;
85 }
86
87 delete [] tempo ;
88 delete [] tempoE ;
89 } // if
d15a28e7 90
9f616d61 91 fDigitsList[fMulDigit] = (int) &digit ;
92 fEnergyList[fMulDigit++] = Energy ;
d15a28e7 93 fAmp += Energy ;
94}
95
96//____________________________________________________________________________
97Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 )
98{
b2a60966 99 // Tells if (true) or not (false) two digits are neighbors)
d15a28e7 100
101 Bool_t aren = kFALSE ;
102
92862013 103 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7 104 Int_t relid1[4] ;
92862013 105 phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
d15a28e7 106
107 Int_t relid2[4] ;
92862013 108 phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
d15a28e7 109
92862013 110 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
111 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7 112
92862013 113 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
d15a28e7 114 aren = kTRUE ;
115
116 return aren ;
117}
118
119//____________________________________________________________________________
120Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
121{
b2a60966 122 // Compares two RecPoints according to their position in the PHOS modules
123
d15a28e7 124 Int_t rv ;
125
126 AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ;
127
128
92862013 129 Int_t phosmod1 = this->GetPHOSMod() ;
130 Int_t phosmod2 = clu->GetPHOSMod() ;
d15a28e7 131
92862013 132 TVector3 locpos1;
133 this->GetLocalPosition(locpos1) ;
134 TVector3 locpos2;
135 clu->GetLocalPosition(locpos2) ;
d15a28e7 136
92862013 137 if(phosmod1 == phosmod2 ) {
138 Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
d15a28e7 139 if (rowdif> 0)
140 rv = -1 ;
141 else if(rowdif < 0)
142 rv = 1 ;
92862013 143 else if(locpos1.Z()>locpos2.Z())
d15a28e7 144 rv = -1 ;
145 else
146 rv = 1 ;
147 }
148
149 else {
92862013 150 if(phosmod1 < phosmod2 )
d15a28e7 151 rv = -1 ;
152 else
153 rv = 1 ;
154 }
155
156 return rv ;
157}
158
9f616d61 159//______________________________________________________________________________
160void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
161{
b2a60966 162 // Execute action corresponding to one event
6ad0bfa0 163 // This member function is called when a AliPHOSRecPoint is clicked with the locator
164 //
165 // If Left button is clicked on AliPHOSRecPoint, the digits are switched on
166 // and switched off when the mouse button is released.
167 //
9f616d61 168
169 // static Int_t pxold, pyold;
170
92862013 171 static TGraph * digitgraph = 0 ;
9f616d61 172
173 if (!gPad->IsEditable()) return;
174
92862013 175 TH2F * histo = 0 ;
176 TCanvas * histocanvas ;
9f616d61 177
178 switch (event) {
179
180 case kButton1Down: {
181 AliPHOSDigit * digit ;
92862013 182 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
9f616d61 183 Int_t iDigit;
184 Int_t relid[4] ;
31aa6d6c 185
186 const Int_t fMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ;
187 Float_t * xi = new Float_t[fMulDigit] ;
188 Float_t * zi = new Float_t[fMulDigit] ;
9f616d61 189
190 // create the histogram for the single cluster
191 // 1. gets histogram boundaries
192 Float_t ximax = -999. ;
193 Float_t zimax = -999. ;
194 Float_t ximin = 999. ;
195 Float_t zimin = 999. ;
196
197 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
198 digit = (AliPHOSDigit *) fDigitsList[iDigit];
92862013 199 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
200 phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
9f616d61 201 if ( xi[iDigit] > ximax )
202 ximax = xi[iDigit] ;
203 if ( xi[iDigit] < ximin )
204 ximin = xi[iDigit] ;
205 if ( zi[iDigit] > zimax )
206 zimax = zi[iDigit] ;
207 if ( zi[iDigit] < zimin )
208 zimin = zi[iDigit] ;
209 }
92862013 210 ximax += phosgeom->GetCrystalSize(0) / 2. ;
211 zimax += phosgeom->GetCrystalSize(2) / 2. ;
212 ximin -= phosgeom->GetCrystalSize(0) / 2. ;
213 zimin -= phosgeom->GetCrystalSize(2) / 2. ;
214 Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5 ) ;
215 Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
9f616d61 216
217 // 2. gets the histogram title
218
219 Text_t title[100] ;
220 sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
221
92862013 222 if (!histo) {
223 delete histo ;
224 histo = 0 ;
9f616d61 225 }
92862013 226 histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
9f616d61 227
228 Float_t x, z ;
229 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
230 digit = (AliPHOSDigit *) fDigitsList[iDigit];
92862013 231 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
232 phosgeom->RelPosInModule(relid, x, z);
233 histo->Fill(x, z, fEnergyList[iDigit] ) ;
9f616d61 234 }
235
92862013 236 if (!digitgraph) {
237 digitgraph = new TGraph(fMulDigit,xi,zi);
238 digitgraph-> SetMarkerStyle(5) ;
239 digitgraph-> SetMarkerSize(1.) ;
240 digitgraph-> SetMarkerColor(1) ;
241 digitgraph-> Paint("P") ;
9f616d61 242 }
243
244 Print() ;
92862013 245 histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ;
246 histocanvas->Draw() ;
247 histo->Draw("lego1") ;
9f616d61 248
31aa6d6c 249 delete[] xi ;
250 delete[] zi ;
251
9f616d61 252 break;
253 }
254
255 case kButton1Up:
92862013 256 if (digitgraph) {
257 delete digitgraph ;
258 digitgraph = 0 ;
9f616d61 259 }
260 break;
261
262 }
263}
264
d15a28e7 265//____________________________________________________________________________
266Float_t AliPHOSEmcRecPoint::GetDispersion()
267{
b2a60966 268 // Calculates the dispersion of the shower at the origine of the RecPoint
269
92862013 270 Float_t d = 0 ;
d15a28e7 271 Float_t wtot = 0 ;
272
92862013 273 TVector3 locpos;
274 GetLocalPosition(locpos);
275 Float_t x = locpos.X() ;
276 Float_t z = locpos.Z() ;
d15a28e7 277
278 AliPHOSDigit * digit ;
92862013 279 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7 280
281 Int_t iDigit;
282 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
283 digit = (AliPHOSDigit *) fDigitsList[iDigit];
284 Int_t relid[4] ;
285 Float_t xi ;
286 Float_t zi ;
92862013 287 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
288 phosgeom->RelPosInModule(relid, xi, zi);
d15a28e7 289 Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
92862013 290 d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ;
d15a28e7 291 wtot+=w ;
292 }
293
92862013 294 d /= wtot ;
d15a28e7 295
92862013 296 return TMath::Sqrt(d) ;
d15a28e7 297}
298
299//____________________________________________________________________________
300void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
301{
b2a60966 302 // Calculates the axis of the shower ellipsoid
303
d15a28e7 304 Float_t wtot = 0. ;
305 Float_t x = 0.;
306 Float_t z = 0.;
92862013 307 Float_t dxx = 0.;
308 Float_t dzz = 0.;
309 Float_t dxz = 0.;
d15a28e7 310
311 AliPHOSDigit * digit ;
92862013 312 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7 313 Int_t iDigit;
314
315 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
316 digit = (AliPHOSDigit *) fDigitsList[iDigit];
317 Int_t relid[4] ;
318 Float_t xi ;
319 Float_t zi ;
92862013 320 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
321 phosgeom->RelPosInModule(relid, xi, zi);
d15a28e7 322 Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
92862013 323 dxx += w * xi * xi ;
d15a28e7 324 x += w * xi ;
92862013 325 dzz += w * zi * zi ;
d15a28e7 326 z += w * zi ;
92862013 327 dxz += w * xi * zi ;
d15a28e7 328 wtot += w ;
329 }
330
92862013 331 dxx /= wtot ;
d15a28e7 332 x /= wtot ;
92862013 333 dxx -= x * x ;
334 dzz /= wtot ;
d15a28e7 335 z /= wtot ;
92862013 336 dzz -= z * z ;
337 dxz /= wtot ;
338 dxz -= x * z ;
d15a28e7 339
92862013 340 lambda[0] = TMath::Sqrt( 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
341 lambda[1] = TMath::Sqrt( 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
d15a28e7 342}
343
344//____________________________________________________________________________
b2a60966 345void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
346{
347 // Calculates the center of gravity in the local PHOS-module coordinates
348
349 if( fLocPos.X() < 1000000.) { // already evaluated
350 LPos = fLocPos ;
351 return ;
352 }
353
354 Float_t wtot = 0. ;
355
356 Int_t relid[4] ;
357
358 Float_t x = 0. ;
359 Float_t z = 0. ;
360
361 AliPHOSDigit * digit ;
362
363 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
364
365 Int_t iDigit;
366
367
368 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
369 digit = (AliPHOSDigit *) fDigitsList[iDigit];
370
371 Float_t xi ;
372 Float_t zi ;
373 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
374 phosgeom->RelPosInModule(relid, xi, zi);
375 Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
376 x += xi * w ;
377 z += zi * w ;
378 wtot += w ;
379
380 }
381
382 x /= wtot ;
383 z /= wtot ;
384 fLocPos.SetX(x) ;
385 fLocPos.SetY(0.) ;
386 fLocPos.SetZ(z) ;
387
388 LPos = fLocPos ;
389}
390
391//____________________________________________________________________________
d15a28e7 392Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
393{
b2a60966 394 // Finds the maximum energy in the cluster
395
d15a28e7 396 Float_t menergy = 0. ;
397
398 Int_t iDigit;
399
400 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
401
402 if(fEnergyList[iDigit] > menergy)
403 menergy = fEnergyList[iDigit] ;
404 }
405 return menergy ;
406}
407
408//____________________________________________________________________________
31aa6d6c 409Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H)
d15a28e7 410{
b2a60966 411 // Calculates the multiplicity of digits with energy larger than H*energy
412
d15a28e7 413 Int_t multipl = 0 ;
414 Int_t iDigit ;
415 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
416
417 if(fEnergyList[iDigit] > H * fAmp)
418 multipl++ ;
419 }
420 return multipl ;
421}
422
423//____________________________________________________________________________
6ad0bfa0 424Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy)
d15a28e7 425{
b2a60966 426 // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
427 // energy difference between two local maxima
428
d15a28e7 429 AliPHOSDigit * digit ;
430 AliPHOSDigit * digitN ;
431
432
433 Int_t iDigitN ;
434 Int_t iDigit ;
435
6ad0bfa0 436 for(iDigit = 0; iDigit < fMulDigit; iDigit++){
d15a28e7 437 maxAt[iDigit] = fDigitsList[iDigit] ;
438 }
439
6ad0bfa0 440 for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
9f616d61 441 if(maxAt[iDigit] != -1) {
442 digit = (AliPHOSDigit *) maxAt[iDigit] ;
d15a28e7 443
6ad0bfa0 444 for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
9f616d61 445 digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ;
d15a28e7 446
9f616d61 447 if ( AreNeighbours(digit, digitN) ) {
d15a28e7 448 if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
449 maxAt[iDigitN] = -1 ;
6ad0bfa0 450 // but may be digit too is not local max ?
451 if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut)
d15a28e7 452 maxAt[iDigit] = -1 ;
453 }
454 else {
455 maxAt[iDigit] = -1 ;
6ad0bfa0 456 // but may be digitN too is not local max ?
457 if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut)
d15a28e7 458 maxAt[iDigitN] = -1 ;
459 }
460 } // if Areneighbours
461 } // while digitN
462 } // slot not empty
463 } // while digit
464
465 iDigitN = 0 ;
6ad0bfa0 466 for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
d15a28e7 467 if(maxAt[iDigit] != -1){
468 maxAt[iDigitN] = maxAt[iDigit] ;
9f616d61 469 maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
470 iDigitN++ ;
d15a28e7 471 }
472 }
473 return iDigitN ;
474}
475
d15a28e7 476
477// //____________________________________________________________________________
478// AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu)
479// {
92862013 480// int * dl = Clu.GetDigitsList() ;
d15a28e7 481
482// if(fDigitsList)
483// delete fDigitsList ;
484
485// AliPHOSDigit * digit ;
486
487// Int_t iDigit;
488
489// for(iDigit=0; iDigit<fMulDigit; iDigit++) {
92862013 490// digit = (AliPHOSDigit *) dl[iDigit];
d15a28e7 491// AddDigit(*digit) ;
492// }
493
494// fAmp = Clu.GetTotalEnergy() ;
495// fGeom = Clu.GetGeom() ;
92862013 496// TVector3 locpos;
497// Clu.GetLocalPosition(locpos) ;
498// fLocPos = locpos;
d15a28e7 499// fMulDigit = Clu.GetMultiplicity() ;
500// fMaxDigit = Clu.GetMaximumMultiplicity() ;
501// fPHOSMod = Clu.GetPHOSMod() ;
502// fW0 = Clu.GetLogWeightCut() ;
503// fDelta = Clu.GetDelta() ;
504// fLocMaxCut = Clu.GetLocMaxCut() ;
505
92862013 506// delete dl ;
d15a28e7 507
508// return *this ;
509// }
510
511//____________________________________________________________________________
512void AliPHOSEmcRecPoint::Print(Option_t * option)
513{
b2a60966 514 // Print the list of digits belonging to the cluster
515
d15a28e7 516 cout << "AliPHOSEmcRecPoint: " << endl ;
517
518 AliPHOSDigit * digit ;
519 Int_t iDigit;
92862013 520 AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7 521
9f616d61 522 Float_t xi ;
523 Float_t zi ;
524 Int_t relid[4] ;
525
d15a28e7 526 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
527 digit = (AliPHOSDigit *) fDigitsList[iDigit];
92862013 528 phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
529 phosgeom->RelPosInModule(relid, xi, zi);
9f616d61 530 cout << " Id = " << digit->GetId() ;
531 cout << " module = " << relid[0] ;
532 cout << " x = " << xi ;
533 cout << " z = " << zi ;
534 cout << " Energy = " << fEnergyList[iDigit] << endl ;
d15a28e7 535 }
536 cout << " Multiplicity = " << fMulDigit << endl ;
537 cout << " Cluster Energy = " << fAmp << endl ;
538
539}
540
541//______________________________________________________________________________
542void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
543{
b2a60966 544 // Stream an object of class AliPHOSEmcRecPoint.
545 // Needed because of the array fEnergyList
546
d15a28e7 547 if (R__b.IsReading()) {
548 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
549 AliPHOSRecPoint::Streamer(R__b);
550 R__b >> fDelta;
551 R__b >> fLocMaxCut;
552 R__b.ReadArray(fEnergyList);
553 R__b >> fW0;
554 } else {
555 R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
556 AliPHOSRecPoint::Streamer(R__b);
557 R__b << fDelta;
558 R__b << fLocMaxCut;
559 R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
560 R__b << fW0;
561 }
562}
563
564