----------------------------------------------------------------------
[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 //_________________________________________________________________________
17 // Rec Point in the PHOS EM calorimeter 
18 //*-- Author : Dmitri Peressounko RRC KI
19 //////////////////////////////////////////////////////////////////////////////
20
21 // --- ROOT system ---
22
23 #include "TMath.h" 
24
25 // --- Standard library ---
26
27 #include "iostream.h"
28
29 // --- AliRoot header files ---
30
31 #include "AliPHOSGeometry.h"
32 #include "AliPHOSEmcRecPoint.h"
33 #include "AliRun.h"
34
35 ClassImp(AliPHOSEmcRecPoint)
36
37 //____________________________________________________________________________
38 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
39   : AliPHOSRecPoint()
40 {
41   // ctor
42
43   fMulDigit   = 0 ;  
44   fAmp   = 0. ;   
45   fEnergyList = new Float_t[fMaxDigit]; 
46   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
47   fDelta     =  PHOSGeom->GetCrystalSize(0) ; 
48   fW0        = W0 ;          
49   fLocMaxCut = LocMaxCut ; 
50   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
51 }
52
53 // //____________________________________________________________________________
54 // AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() 
55 // {
56 //   // dtor 
57 // }
58
59 //____________________________________________________________________________
60 void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
61 {
62   // adds a digit to the digits list
63   // and accumulates the total amplitude and the multiplicity 
64   
65   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
66     int * tempo = new ( int[fMaxDigit*=2] ) ; 
67     Float_t * tempoE =  new ( Float_t[fMaxDigit*=2] ) ;
68     Int_t index ; 
69     
70     for ( index = 0 ; index < fMulDigit ; index++ ){
71       tempo[index] = fDigitsList[index] ;
72       tempoE[index] = fEnergyList[index] ; 
73     }
74     
75     delete fDigitsList ; 
76     delete fEnergyList ;
77     fDigitsList = tempo ;
78     fEnergyList = tempoE ; 
79   }
80   
81   fDigitsList[fMulDigit]  =  (int) &digit  ; 
82   fEnergyList[fMulDigit++]= Energy ;
83   fAmp += Energy ; 
84 }
85
86 //____________________________________________________________________________
87 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
88 {
89   
90   Bool_t aren = kFALSE ;
91   
92   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
93   Int_t relid1[4] ; 
94   PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
95
96   Int_t relid2[4] ; 
97   PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
98   
99   Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
100   Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
101
102   if (( ColDiff<=1 )  && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0)) 
103     aren = kTRUE ;
104   
105   return aren ;
106 }
107
108 //____________________________________________________________________________
109 Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
110 {
111   Int_t rv ; 
112
113   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; 
114
115  
116   Int_t PHOSMod1 = this->GetPHOSMod() ;
117   Int_t PHOSMod2 = clu->GetPHOSMod() ;
118
119   TVector3 LocPos1; 
120   this->GetLocalPosition(LocPos1) ;
121   TVector3 LocPos2;  
122   clu->GetLocalPosition(LocPos2) ;  
123
124   if(PHOSMod1 == PHOSMod2 ) {
125     Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ;
126     if (rowdif> 0) 
127       rv = -1 ;
128      else if(rowdif < 0) 
129        rv = 1 ;
130     else if(LocPos1.Z()>LocPos2.Z()) 
131       rv = -1 ;
132     else 
133       rv = 1 ; 
134      }
135
136   else {
137     if(PHOSMod1 < PHOSMod2 ) 
138       rv = -1 ;
139     else 
140       rv = 1 ;
141   }
142
143   return rv ; 
144 }
145
146 //____________________________________________________________________________
147 Float_t  AliPHOSEmcRecPoint::GetDispersion() 
148 {
149   Float_t D    = 0 ;
150   Float_t wtot = 0 ;
151
152   TVector3 LocPos;
153   GetLocalPosition(LocPos);
154   Float_t x = LocPos.X() ;
155   Float_t z = LocPos.Z() ;
156   //  Int_t i = GetPHOSMod() ;
157
158   AliPHOSDigit * digit ;
159   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
160   
161   Int_t iDigit;
162   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
163     digit = (AliPHOSDigit *) fDigitsList[iDigit];
164     Int_t relid[4] ;
165     Float_t xi ;
166     Float_t zi ;
167     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
168     PHOSGeom->RelPosInModule(relid, xi, zi);
169     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
170     D += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
171     wtot+=w ;
172   }
173
174   D /= wtot ;
175
176   return TMath::Sqrt(D) ;
177 }
178
179 //____________________________________________________________________________
180 void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
181 {
182   Float_t wtot = 0. ;
183   Float_t x    = 0.;
184   Float_t z    = 0.;
185   Float_t Dxx  = 0.;
186   Float_t Dzz  = 0.;
187   Float_t Dxz  = 0.;
188
189   AliPHOSDigit * digit ;
190   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
191   Int_t iDigit;
192
193   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
194     digit = (AliPHOSDigit *) fDigitsList[iDigit];
195     Int_t relid[4] ;
196     Float_t xi ;
197     Float_t zi ;
198     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
199     PHOSGeom->RelPosInModule(relid, xi, zi);
200     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
201     Dxx  += w * xi * xi ;
202     x    += w * xi ;
203     Dzz  += w * zi * zi ;
204     z    += w * zi ; 
205     Dxz  += w * xi * zi ; 
206     wtot += w ;
207   }
208
209   Dxx /= wtot ;
210   x   /= wtot ;
211   Dxx -= x * x ;
212   Dzz /= wtot ;
213   z   /= wtot ;
214   Dzz -= z * z ;
215   Dxz /= wtot ;
216   Dxz -= x * z ;
217
218   lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
219   lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
220 }
221
222 //____________________________________________________________________________
223 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
224 {
225   Float_t menergy = 0. ;
226
227   Int_t iDigit;
228
229   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
230  
231     if(fEnergyList[iDigit] > menergy) 
232       menergy = fEnergyList[iDigit] ;
233   }
234   return menergy ;
235 }
236
237 //____________________________________________________________________________
238 Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) 
239 {
240   Int_t multipl   = 0 ;
241   Int_t iDigit ;
242   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
243
244     if(fEnergyList[iDigit] > H * fAmp) 
245       multipl++ ;
246   }
247   return multipl ;
248 }
249
250 //____________________________________________________________________________
251 Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(int *  maxAt, Float_t * maxAtEnergy) 
252
253   AliPHOSDigit * digit ;
254   AliPHOSDigit * digitN ;
255   
256
257   Int_t iDigitN ;
258   Int_t iDigit ;
259
260   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
261     maxAt[iDigit] = fDigitsList[iDigit] ;
262   }
263   
264   for(iDigit=0 ; iDigit<fMulDigit; iDigit++) {   
265     if(maxAt[iDigit]!= -1) {
266       digit = (AliPHOSDigit *) maxAt[iDigit];
267          
268       for(iDigitN=0; iDigitN<fMulDigit; iDigitN++) {    
269         digitN = (AliPHOSDigit *) fDigitsList[iDigitN]; 
270         
271         if ( AreNeighbours(digit,digitN) ) {
272           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
273             maxAt[iDigitN] = -1 ;
274             //but may be digit is not local max too ?
275             if(fEnergyList[iDigit]<fEnergyList[iDigitN]+fLocMaxCut) 
276               maxAt[iDigit] = -1 ;
277           }
278           else {
279             maxAt[iDigit] = -1 ;
280             //but may be digitN is not local max too ?
281             if(fEnergyList[iDigit] >fEnergyList[iDigitN]-fLocMaxCut) 
282               maxAt[iDigitN] = -1 ; 
283           } 
284         } // if Areneighbours
285       } // while digitN
286     } // slot not empty
287   } // while digit
288   
289   iDigitN = 0 ;
290   for(iDigit=0; iDigit<fMulDigit; iDigit++) { 
291     if(maxAt[iDigit] != -1){
292       maxAt[iDigitN] = maxAt[iDigit] ;
293       maxAtEnergy[iDigitN++] = fEnergyList[iDigit] ;
294     }
295   }
296   return iDigitN ;
297 }
298
299 //____________________________________________________________________________
300 void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
301 {
302   if( fLocPos.X() < 1000000.) { //allready evaluated
303    LPos = fLocPos ;
304    return ;
305   }
306
307   Float_t wtot = 0. ;
308  
309   Int_t relid[4] ;
310
311   Float_t x = 0. ;
312   Float_t z = 0. ;
313   
314   AliPHOSDigit * digit ;
315
316   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
317
318   Int_t iDigit;
319
320
321   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
322     digit = (AliPHOSDigit *) fDigitsList[iDigit];
323
324     Float_t xi ;
325     Float_t zi ;
326     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
327     PHOSGeom->RelPosInModule(relid, xi, zi);
328     Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
329     x    += xi * w ;
330     z    += zi * w ;
331     wtot += w ;
332
333   }
334
335   x /= wtot ;
336   z /= wtot ;
337   fLocPos.SetX(x)  ;
338   fLocPos.SetY(0.) ;
339   fLocPos.SetZ(z)  ;
340
341   LPos = fLocPos ;
342 }
343
344 // //____________________________________________________________________________
345 // AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) 
346 // {
347 //   int * DL = Clu.GetDigitsList() ; 
348  
349 //  if(fDigitsList) 
350 //     delete fDigitsList  ;
351
352 //   AliPHOSDigit * digit ;
353  
354 //   Int_t iDigit;
355
356 //   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
357 //     digit = (AliPHOSDigit *) DL[iDigit];
358 //     AddDigit(*digit) ;
359 //   }
360
361 //   fAmp       = Clu.GetTotalEnergy() ;
362 //   fGeom      = Clu.GetGeom() ;
363 //   TVector3 LocPos;
364 //   Clu.GetLocalPosition(LocPos) ;
365 //   fLocPos    = LocPos;
366 //   fMulDigit       = Clu.GetMultiplicity() ;
367 //   fMaxDigit  = Clu.GetMaximumMultiplicity() ;
368 //   fPHOSMod   = Clu.GetPHOSMod() ;
369 //   fW0        = Clu.GetLogWeightCut() ; 
370 //   fDelta     = Clu.GetDelta() ;
371 //   fLocMaxCut = Clu.GetLocMaxCut() ;
372   
373 //   delete DL ; 
374  
375 //   return *this ;
376 // }
377
378 //____________________________________________________________________________
379 void AliPHOSEmcRecPoint::Print(Option_t * option) 
380 {
381   cout << "AliPHOSEmcRecPoint: " << endl ;
382
383   AliPHOSDigit * digit ; 
384   Int_t iDigit;
385
386   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
387     digit = (AliPHOSDigit *) fDigitsList[iDigit];
388     cout << "     digit Id          = " << digit->GetId()  
389          << "     digit Energy      = " << fEnergyList[iDigit] << endl ;
390   }
391   cout << "       Multiplicity    = " << fMulDigit  << endl ;
392   cout << "       Cluster Energy  = " << fAmp << endl ;
393   
394 }
395
396 //______________________________________________________________________________
397 void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
398 {
399    // Stream an object of class AliPHOSEmcRecPoint.
400
401    if (R__b.IsReading()) {
402       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
403       AliPHOSRecPoint::Streamer(R__b);
404       R__b >> fDelta;
405       R__b >> fLocMaxCut;
406       R__b.ReadArray(fEnergyList);
407       R__b >> fW0;
408    } else {
409       R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
410       AliPHOSRecPoint::Streamer(R__b);
411       R__b << fDelta;
412       R__b << fLocMaxCut;
413       R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
414       R__b << fW0;
415    }
416 }
417
418