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 | |
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 | |