git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
----------------------------------------------------------------------
[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
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
35ClassImp(AliPHOSEmcRecPoint)
36
37//____________________________________________________________________________
38AliPHOSEmcRecPoint::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//____________________________________________________________________________
60void 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//____________________________________________________________________________
87Bool_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//____________________________________________________________________________
109Int_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//____________________________________________________________________________
147Float_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//____________________________________________________________________________
180void 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//____________________________________________________________________________
223Float_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//____________________________________________________________________________
238Int_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//____________________________________________________________________________
251Int_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//____________________________________________________________________________
300void 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//____________________________________________________________________________
379void 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//______________________________________________________________________________
397void 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