0daf2d5401cb747d2e1869f1df35fd117156d3c0
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
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-EMC 
20 //  An EmcRecPoint 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
33 #include <iostream.h> 
34
35 // --- AliRoot header files ---
36
37 #include "AliPHOSGeometry.h"
38 #include "AliPHOSEmcRecPoint.h"
39 #include "AliRun.h"
40
41 ClassImp(AliPHOSEmcRecPoint)
42
43 //____________________________________________________________________________
44 AliPHOSEmcRecPoint::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]; 
52   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
53   fDelta     =  phosgeom->GetCrystalSize(0) ; 
54   fW0        = W0 ;          
55   fLocMaxCut = LocMaxCut ; 
56   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
57 }
58
59 //____________________________________________________________________________
60 void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
61 {
62   // Adds a digit to the RecPoint
63   //  and accumulates the total amplitude and the multiplicity 
64   
65   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
66     fMaxDigit*=2 ; 
67     int * tempo = new ( int[fMaxDigit] ) ; 
68     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
69
70     Int_t index ;     
71     for ( index = 0 ; index < fMulDigit ; index++ ){
72       tempo[index] = fDigitsList[index] ;
73       tempoE[index] = fEnergyList[index] ; 
74     }
75     
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
90   
91   fDigitsList[fMulDigit]   =  (int) &digit  ; 
92   fEnergyList[fMulDigit++] = Energy ;
93   fAmp += Energy ; 
94 }
95
96 //____________________________________________________________________________
97 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
98 {
99   // Tells if (true) or not (false) two digits are neighbors)
100   
101   Bool_t aren = kFALSE ;
102   
103   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
104   Int_t relid1[4] ; 
105   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
106
107   Int_t relid2[4] ; 
108   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
109   
110   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
111   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
112
113   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
114     aren = kTRUE ;
115   
116   return aren ;
117 }
118
119 //____________________________________________________________________________
120 Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
121 {
122   // Compares two RecPoints according to their position in the PHOS modules
123
124   Int_t rv ; 
125
126   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; 
127
128  
129   Int_t phosmod1 = this->GetPHOSMod() ;
130   Int_t phosmod2 = clu->GetPHOSMod() ;
131
132   TVector3 locpos1; 
133   this->GetLocalPosition(locpos1) ;
134   TVector3 locpos2;  
135   clu->GetLocalPosition(locpos2) ;  
136
137   if(phosmod1 == phosmod2 ) {
138     Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
139     if (rowdif> 0) 
140       rv = -1 ;
141      else if(rowdif < 0) 
142        rv = 1 ;
143     else if(locpos1.Z()>locpos2.Z()) 
144       rv = -1 ;
145     else 
146       rv = 1 ; 
147      }
148
149   else {
150     if(phosmod1 < phosmod2 ) 
151       rv = -1 ;
152     else 
153       rv = 1 ;
154   }
155
156   return rv ; 
157 }
158
159 //______________________________________________________________________________
160 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
161 {
162   // Execute action corresponding to one event
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   //
168
169   //   static Int_t pxold, pyold;
170
171    static TGraph *  digitgraph = 0 ;
172
173    if (!gPad->IsEditable()) return;
174
175    TH2F * histo = 0 ;
176    TCanvas * histocanvas ; 
177    
178    switch (event) {
179    
180    case kButton1Down: {
181      AliPHOSDigit * digit ;
182      AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
183      Int_t iDigit;
184      Int_t relid[4] ;
185      
186      const Int_t fMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ; 
187      Float_t * xi = new Float_t[fMulDigit] ; 
188      Float_t * zi = new Float_t[fMulDigit] ; 
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];
199        phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
200        phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
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      }
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 ) ;
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   
222      if (!histo) {
223        delete histo ; 
224        histo = 0 ; 
225      }
226      histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
227
228      Float_t x, z ; 
229      for(iDigit=0; iDigit<fMulDigit; iDigit++) {
230        digit = (AliPHOSDigit *) fDigitsList[iDigit];
231        phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
232        phosgeom->RelPosInModule(relid, x, z);
233        histo->Fill(x, z, fEnergyList[iDigit] ) ;
234      }
235
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") ;
242      }
243
244      Print() ;
245      histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
246      histocanvas->Draw() ; 
247      histo->Draw("lego1") ; 
248
249      delete[] xi ; 
250      delete[] zi ; 
251      
252      break;
253    }
254
255    case kButton1Up: 
256      if (digitgraph) {
257        delete digitgraph  ;
258        digitgraph = 0 ;
259      }
260      break;
261   
262    }
263 }
264
265 //____________________________________________________________________________
266 Float_t  AliPHOSEmcRecPoint::GetDispersion() 
267 {
268   // Calculates the dispersion of the shower at the origine of the RecPoint
269
270   Float_t d    = 0 ;
271   Float_t wtot = 0 ;
272
273   TVector3 locpos;
274   GetLocalPosition(locpos);
275   Float_t x = locpos.X() ;
276   Float_t z = locpos.Z() ;
277
278   AliPHOSDigit * digit ;
279   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
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 ;
287     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
288     phosgeom->RelPosInModule(relid, xi, zi);
289     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
290     d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
291     wtot+=w ;
292   }
293
294   d /= wtot ;
295
296   return TMath::Sqrt(d) ;
297 }
298
299 //____________________________________________________________________________
300 void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
301 {
302   // Calculates the axis of the shower ellipsoid
303   
304   Float_t wtot = 0. ;
305   Float_t x    = 0.;
306   Float_t z    = 0.;
307   Float_t dxx  = 0.;
308   Float_t dzz  = 0.;
309   Float_t dxz  = 0.;
310
311   AliPHOSDigit * digit ;
312   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
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 ;
320     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
321     phosgeom->RelPosInModule(relid, xi, zi);
322     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
323     dxx  += w * xi * xi ;
324     x    += w * xi ;
325     dzz  += w * zi * zi ;
326     z    += w * zi ; 
327     dxz  += w * xi * zi ; 
328     wtot += w ;
329   }
330
331   dxx /= wtot ;
332   x   /= wtot ;
333   dxx -= x * x ;
334   dzz /= wtot ;
335   z   /= wtot ;
336   dzz -= z * z ;
337   dxz /= wtot ;
338   dxz -= x * z ;
339
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 ) ) ;
342 }
343
344 //____________________________________________________________________________
345 void 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 //____________________________________________________________________________
392 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
393 {
394   // Finds the maximum energy in the cluster
395   
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 //____________________________________________________________________________
409 Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) 
410 {
411   // Calculates the multiplicity of digits with energy larger than H*energy 
412   
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 //____________________________________________________________________________
424 Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
425
426   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
427   //  energy difference between two local maxima
428
429   AliPHOSDigit * digit ;
430   AliPHOSDigit * digitN ;
431   
432
433   Int_t iDigitN ;
434   Int_t iDigit ;
435
436   for(iDigit = 0; iDigit < fMulDigit; iDigit++){
437     maxAt[iDigit] = fDigitsList[iDigit] ;
438   }
439   
440   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
441     if(maxAt[iDigit] != -1) {
442       digit = (AliPHOSDigit *) maxAt[iDigit] ;
443          
444       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {        
445         digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ; 
446         
447         if ( AreNeighbours(digit, digitN) ) {
448           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
449             maxAt[iDigitN] = -1 ;
450             // but may be digit too is not local max ?
451             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) 
452               maxAt[iDigit] = -1 ;
453           }
454           else {
455             maxAt[iDigit] = -1 ;
456             // but may be digitN too is not local max ?
457             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) 
458               maxAt[iDigitN] = -1 ; 
459           } 
460         } // if Areneighbours
461       } // while digitN
462     } // slot not empty
463   } // while digit
464   
465   iDigitN = 0 ;
466   for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
467     if(maxAt[iDigit] != -1){
468       maxAt[iDigitN] = maxAt[iDigit] ;
469       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
470       iDigitN++ ; 
471     }
472   }
473   return iDigitN ;
474 }
475
476
477 // //____________________________________________________________________________
478 // AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) 
479 // {
480 //   int * dl = Clu.GetDigitsList() ; 
481  
482 //  if(fDigitsList) 
483 //     delete fDigitsList  ;
484
485 //   AliPHOSDigit * digit ;
486  
487 //   Int_t iDigit;
488
489 //   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
490 //     digit = (AliPHOSDigit *) dl[iDigit];
491 //     AddDigit(*digit) ;
492 //   }
493
494 //   fAmp       = Clu.GetTotalEnergy() ;
495 //   fGeom      = Clu.GetGeom() ;
496 //   TVector3 locpos;
497 //   Clu.GetLocalPosition(locpos) ;
498 //   fLocPos    = locpos;
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   
506 //   delete dl ; 
507  
508 //   return *this ;
509 // }
510
511 //____________________________________________________________________________
512 void AliPHOSEmcRecPoint::Print(Option_t * option) 
513 {
514   // Print the list of digits belonging to the cluster
515   
516   cout << "AliPHOSEmcRecPoint: " << endl ;
517
518   AliPHOSDigit * digit ; 
519   Int_t iDigit;
520   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
521
522   Float_t xi ;
523   Float_t zi ;
524   Int_t relid[4] ; 
525  
526   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
527     digit = (AliPHOSDigit *) fDigitsList[iDigit];
528     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
529     phosgeom->RelPosInModule(relid, xi, zi);
530     cout << " Id = " << digit->GetId() ;  
531     cout << "   module  = " << relid[0] ;  
532     cout << "   x  = " << xi ;  
533     cout << "   z  = " << zi ;  
534     cout << "   Energy = " << fEnergyList[iDigit] << endl ;
535   }
536   cout << "       Multiplicity    = " << fMulDigit  << endl ;
537   cout << "       Cluster Energy  = " << fAmp << endl ;
538   
539 }
540
541 //______________________________________________________________________________
542 void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
543 {
544   // Stream an object of class AliPHOSEmcRecPoint.
545   // Needed because of the array fEnergyList
546   
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