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