]>
Commit | Line | Data |
---|---|---|
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 | ||
b2a60966 | 16 | /* $Id$ */ |
17 | ||
d15a28e7 | 18 | //_________________________________________________________________________ |
b2a60966 | 19 | // RecPoint implementation for PHOS-EMC |
20 | // An EmcRecPoint is a cluster of digits | |
baef0810 | 21 | //*-- |
b2a60966 | 22 | //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH) |
23 | ||
d15a28e7 | 24 | |
25 | // --- ROOT system --- | |
9f616d61 | 26 | #include "TPad.h" |
27 | #include "TH2.h" | |
d15a28e7 | 28 | #include "TMath.h" |
9f616d61 | 29 | #include "TCanvas.h" |
55fe9d13 | 30 | #include "TGraph.h" |
d15a28e7 | 31 | |
32 | // --- Standard library --- | |
33 | ||
d15a28e7 | 34 | // --- AliRoot header files --- |
35 | ||
a6eedfad | 36 | #include "AliGenerator.h" |
d15a28e7 | 37 | #include "AliPHOSGeometry.h" |
38 | #include "AliPHOSEmcRecPoint.h" | |
39 | #include "AliRun.h" | |
7b7c1533 | 40 | #include "AliPHOSGetter.h" |
d15a28e7 | 41 | |
42 | ClassImp(AliPHOSEmcRecPoint) | |
43 | ||
44 | //____________________________________________________________________________ | |
7932f811 | 45 | AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint() |
d15a28e7 | 46 | { |
47 | // ctor | |
48 | ||
49 | fMulDigit = 0 ; | |
50 | fAmp = 0. ; | |
7932f811 | 51 | fCoreEnergy = 0 ; |
52 | fEnergyList = 0 ; | |
94256c23 | 53 | fNExMax = 0 ; //Not unfolded yet |
9688c1dd | 54 | fTime = -1. ; |
a6eedfad | 55 | fLocPos.SetX(1000000.) ; //Local position should be evaluated |
ce2a9a95 | 56 | fDebug=0; |
7b7c1533 | 57 | |
83974468 | 58 | } |
59 | ||
73a68ccb | 60 | //____________________________________________________________________________ |
61 | AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : AliPHOSRecPoint(opt) | |
62 | { | |
63 | // ctor | |
64 | ||
65 | fMulDigit = 0 ; | |
66 | fAmp = 0. ; | |
94256c23 | 67 | fNExMax = 0 ; //Not unfolded yet |
73a68ccb | 68 | fCoreEnergy = 0 ; |
69 | fEnergyList = 0 ; | |
70 | fTime = -1. ; | |
71 | fLocPos.SetX(1000000.) ; //Local position should be evaluated | |
ce2a9a95 | 72 | fDebug=0; |
73a68ccb | 73 | |
74 | } | |
75 | ||
55fe9d13 | 76 | //____________________________________________________________________________ |
77 | AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const AliPHOSEmcRecPoint & rp) : AliPHOSRecPoint(rp) | |
78 | { | |
79 | // cpy ctor | |
80 | ||
81 | fMulDigit = rp.fMulDigit ; | |
82 | fAmp = rp.fAmp ; | |
83 | fCoreEnergy = rp.fCoreEnergy ; | |
84 | fEnergyList = new Float_t[rp.fMulDigit] ; | |
85 | Int_t index ; | |
86 | for(index = 0 ; index < fMulDigit ; index++) | |
87 | fEnergyList[index] = rp.fEnergyList[index] ; | |
88 | fNExMax = rp.fNExMax ; | |
89 | fTime = rp.fTime ; | |
90 | } | |
91 | ||
83974468 | 92 | //____________________________________________________________________________ |
93 | AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() | |
94 | { | |
88714635 | 95 | // dtor |
a4e98857 | 96 | |
83974468 | 97 | if ( fEnergyList ) |
98 | delete[] fEnergyList ; | |
d15a28e7 | 99 | } |
100 | ||
d15a28e7 | 101 | //____________________________________________________________________________ |
83974468 | 102 | void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) |
d15a28e7 | 103 | { |
b2a60966 | 104 | // Adds a digit to the RecPoint |
a4e98857 | 105 | // and accumulates the total amplitude and the multiplicity |
d15a28e7 | 106 | |
7932f811 | 107 | if(fEnergyList == 0) |
108 | fEnergyList = new Float_t[fMaxDigit]; | |
109 | ||
d15a28e7 | 110 | if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists |
9f616d61 | 111 | fMaxDigit*=2 ; |
83974468 | 112 | Int_t * tempo = new ( Int_t[fMaxDigit] ) ; |
9f616d61 | 113 | Float_t * tempoE = new ( Float_t[fMaxDigit] ) ; |
114 | ||
115 | Int_t index ; | |
d15a28e7 | 116 | for ( index = 0 ; index < fMulDigit ; index++ ){ |
83974468 | 117 | tempo[index] = fDigitsList[index] ; |
d15a28e7 | 118 | tempoE[index] = fEnergyList[index] ; |
119 | } | |
120 | ||
9f616d61 | 121 | delete [] fDigitsList ; |
83974468 | 122 | fDigitsList = new ( Int_t[fMaxDigit] ) ; |
9f616d61 | 123 | |
124 | delete [] fEnergyList ; | |
125 | fEnergyList = new ( Float_t[fMaxDigit] ) ; | |
126 | ||
127 | for ( index = 0 ; index < fMulDigit ; index++ ){ | |
128 | fDigitsList[index] = tempo[index] ; | |
129 | fEnergyList[index] = tempoE[index] ; | |
130 | } | |
131 | ||
132 | delete [] tempo ; | |
133 | delete [] tempoE ; | |
134 | } // if | |
d15a28e7 | 135 | |
83974468 | 136 | fDigitsList[fMulDigit] = digit.GetIndexInList() ; |
137 | fEnergyList[fMulDigit] = Energy ; | |
138 | fMulDigit++ ; | |
d15a28e7 | 139 | fAmp += Energy ; |
7932f811 | 140 | |
141 | EvalPHOSMod(&digit) ; | |
d15a28e7 | 142 | } |
143 | ||
144 | //____________________________________________________________________________ | |
ad8cfaf4 | 145 | Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const |
d15a28e7 | 146 | { |
a4e98857 | 147 | // Tells if (true) or not (false) two digits are neighbors |
d15a28e7 | 148 | |
149 | Bool_t aren = kFALSE ; | |
150 | ||
7b7c1533 | 151 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
152 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
153 | ||
d15a28e7 | 154 | Int_t relid1[4] ; |
92862013 | 155 | phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; |
d15a28e7 | 156 | |
157 | Int_t relid2[4] ; | |
92862013 | 158 | phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; |
d15a28e7 | 159 | |
92862013 | 160 | Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; |
161 | Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; | |
d15a28e7 | 162 | |
92862013 | 163 | if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) |
d15a28e7 | 164 | aren = kTRUE ; |
165 | ||
166 | return aren ; | |
167 | } | |
168 | ||
169 | //____________________________________________________________________________ | |
2a941f4e | 170 | Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const |
d15a28e7 | 171 | { |
b2a60966 | 172 | // Compares two RecPoints according to their position in the PHOS modules |
173 | ||
7932f811 | 174 | Float_t delta = 1 ; //Width of "Sorting row". If you changibg this |
175 | //value (what is senseless) change as vell delta in | |
176 | //AliPHOSTrackSegmentMakerv* and other RecPoints... | |
d15a28e7 | 177 | Int_t rv ; |
178 | ||
179 | AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; | |
180 | ||
181 | ||
ad8cfaf4 | 182 | Int_t phosmod1 = GetPHOSMod() ; |
92862013 | 183 | Int_t phosmod2 = clu->GetPHOSMod() ; |
d15a28e7 | 184 | |
92862013 | 185 | TVector3 locpos1; |
7932f811 | 186 | GetLocalPosition(locpos1) ; |
92862013 | 187 | TVector3 locpos2; |
188 | clu->GetLocalPosition(locpos2) ; | |
d15a28e7 | 189 | |
92862013 | 190 | if(phosmod1 == phosmod2 ) { |
7932f811 | 191 | Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ; |
d15a28e7 | 192 | if (rowdif> 0) |
7932f811 | 193 | rv = 1 ; |
d15a28e7 | 194 | else if(rowdif < 0) |
7932f811 | 195 | rv = -1 ; |
92862013 | 196 | else if(locpos1.Z()>locpos2.Z()) |
d15a28e7 | 197 | rv = -1 ; |
198 | else | |
199 | rv = 1 ; | |
200 | } | |
201 | ||
202 | else { | |
92862013 | 203 | if(phosmod1 < phosmod2 ) |
d15a28e7 | 204 | rv = -1 ; |
205 | else | |
206 | rv = 1 ; | |
207 | } | |
208 | ||
209 | return rv ; | |
210 | } | |
9f616d61 | 211 | //______________________________________________________________________________ |
baef0810 | 212 | void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const |
9f616d61 | 213 | { |
9f616d61 | 214 | |
9688c1dd | 215 | // Execute action corresponding to one event |
216 | // This member function is called when a AliPHOSRecPoint is clicked with the locator | |
217 | // | |
218 | // If Left button is clicked on AliPHOSRecPoint, the digits are switched on | |
219 | // and switched off when the mouse button is released. | |
220 | ||
221 | ||
222 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; | |
223 | if(!gime) return ; | |
224 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
83974468 | 225 | |
9688c1dd | 226 | static TGraph * digitgraph = 0 ; |
83974468 | 227 | |
9688c1dd | 228 | if (!gPad->IsEditable()) return; |
83974468 | 229 | |
9688c1dd | 230 | TH2F * histo = 0 ; |
231 | TCanvas * histocanvas ; | |
232 | ||
afa51c4e | 233 | const TClonesArray * digits = gime->Digits() ; |
9688c1dd | 234 | |
235 | switch (event) { | |
83974468 | 236 | |
9688c1dd | 237 | case kButton1Down: { |
238 | AliPHOSDigit * digit ; | |
239 | Int_t iDigit; | |
240 | Int_t relid[4] ; | |
83974468 | 241 | |
9688c1dd | 242 | const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ; |
243 | Float_t * xi = new Float_t[kMulDigit] ; | |
244 | Float_t * zi = new Float_t[kMulDigit] ; | |
83974468 | 245 | |
9688c1dd | 246 | // create the histogram for the single cluster |
247 | // 1. gets histogram boundaries | |
248 | Float_t ximax = -999. ; | |
249 | Float_t zimax = -999. ; | |
250 | Float_t ximin = 999. ; | |
251 | Float_t zimin = 999. ; | |
83974468 | 252 | |
9688c1dd | 253 | for(iDigit=0; iDigit<kMulDigit; iDigit++) { |
254 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
255 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
256 | phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]); | |
257 | if ( xi[iDigit] > ximax ) | |
258 | ximax = xi[iDigit] ; | |
259 | if ( xi[iDigit] < ximin ) | |
260 | ximin = xi[iDigit] ; | |
261 | if ( zi[iDigit] > zimax ) | |
262 | zimax = zi[iDigit] ; | |
263 | if ( zi[iDigit] < zimin ) | |
264 | zimin = zi[iDigit] ; | |
265 | } | |
266 | ximax += phosgeom->GetCrystalSize(0) / 2. ; | |
267 | zimax += phosgeom->GetCrystalSize(2) / 2. ; | |
268 | ximin -= phosgeom->GetCrystalSize(0) / 2. ; | |
269 | zimin -= phosgeom->GetCrystalSize(2) / 2. ; | |
270 | Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5 ) ; | |
271 | Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ; | |
83974468 | 272 | |
9688c1dd | 273 | // 2. gets the histogram title |
83974468 | 274 | |
9688c1dd | 275 | Text_t title[100] ; |
276 | sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ; | |
83974468 | 277 | |
9688c1dd | 278 | if (!histo) { |
279 | delete histo ; | |
280 | histo = 0 ; | |
281 | } | |
282 | histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ; | |
83974468 | 283 | |
9688c1dd | 284 | Float_t x, z ; |
285 | for(iDigit=0; iDigit<kMulDigit; iDigit++) { | |
286 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
287 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
288 | phosgeom->RelPosInModule(relid, x, z); | |
289 | histo->Fill(x, z, fEnergyList[iDigit] ) ; | |
290 | } | |
83974468 | 291 | |
9688c1dd | 292 | if (!digitgraph) { |
293 | digitgraph = new TGraph(kMulDigit,xi,zi); | |
294 | digitgraph-> SetMarkerStyle(5) ; | |
295 | digitgraph-> SetMarkerSize(1.) ; | |
296 | digitgraph-> SetMarkerColor(1) ; | |
297 | digitgraph-> Paint("P") ; | |
298 | } | |
83974468 | 299 | |
9688c1dd | 300 | // Print() ; |
301 | histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ; | |
302 | histocanvas->Draw() ; | |
303 | histo->Draw("lego1") ; | |
83974468 | 304 | |
9688c1dd | 305 | delete[] xi ; |
306 | delete[] zi ; | |
83974468 | 307 | |
9688c1dd | 308 | break; |
309 | } | |
83974468 | 310 | |
9688c1dd | 311 | case kButton1Up: |
312 | if (digitgraph) { | |
313 | delete digitgraph ; | |
314 | digitgraph = 0 ; | |
315 | } | |
316 | break; | |
9f616d61 | 317 | |
9688c1dd | 318 | } |
9f616d61 | 319 | } |
320 | ||
d15a28e7 | 321 | //____________________________________________________________________________ |
7932f811 | 322 | void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits) |
d15a28e7 | 323 | { |
b2a60966 | 324 | // Calculates the dispersion of the shower at the origine of the RecPoint |
325 | ||
e5b16749 | 326 | Float_t d = 0. ; |
327 | Float_t wtot = 0. ; | |
d15a28e7 | 328 | |
d084d50d | 329 | Float_t x = 0.; |
330 | Float_t z = 0.; | |
d15a28e7 | 331 | |
332 | AliPHOSDigit * digit ; | |
7b7c1533 | 333 | |
334 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; | |
335 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
d15a28e7 | 336 | |
e5b16749 | 337 | |
338 | // Calculates the center of gravity in the local PHOS-module coordinates | |
339 | ||
d15a28e7 | 340 | Int_t iDigit; |
e5b16749 | 341 | |
342 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
343 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
344 | Int_t relid[4] ; | |
345 | Float_t xi ; | |
346 | Float_t zi ; | |
347 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
348 | phosgeom->RelPosInModule(relid, xi, zi); | |
349 | Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; | |
350 | x += xi * w ; | |
351 | z += zi * w ; | |
352 | wtot += w ; | |
353 | } | |
354 | x /= wtot ; | |
355 | z /= wtot ; | |
356 | ||
357 | ||
358 | // Calculates the dispersion in coordinates | |
359 | wtot = 0.; | |
88714635 | 360 | for(iDigit=0; iDigit < fMulDigit; iDigit++) { |
7932f811 | 361 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; |
d15a28e7 | 362 | Int_t relid[4] ; |
363 | Float_t xi ; | |
364 | Float_t zi ; | |
92862013 | 365 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; |
366 | phosgeom->RelPosInModule(relid, xi, zi); | |
7932f811 | 367 | Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; |
92862013 | 368 | d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; |
d15a28e7 | 369 | wtot+=w ; |
d084d50d | 370 | |
371 | ||
d15a28e7 | 372 | } |
e5b16749 | 373 | |
7932f811 | 374 | |
92862013 | 375 | d /= wtot ; |
d15a28e7 | 376 | |
7932f811 | 377 | fDispersion = TMath::Sqrt(d) ; |
c63c49e9 | 378 | |
d15a28e7 | 379 | } |
fad3e5b9 | 380 | //______________________________________________________________________________ |
e5b16749 | 381 | void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) |
fad3e5b9 | 382 | { |
a4e98857 | 383 | // This function calculates energy in the core, |
384 | // i.e. within a radius rad = 3cm around the center. Beyond this radius | |
385 | // in accordance with shower profile the energy deposition | |
fad3e5b9 | 386 | // should be less than 2% |
387 | ||
fad3e5b9 | 388 | Float_t coreRadius = 3 ; |
389 | ||
e5b16749 | 390 | Float_t x = 0 ; |
391 | Float_t z = 0 ; | |
fad3e5b9 | 392 | |
393 | AliPHOSDigit * digit ; | |
7b7c1533 | 394 | |
395 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; | |
396 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
397 | ||
fad3e5b9 | 398 | Int_t iDigit; |
e5b16749 | 399 | |
400 | // Calculates the center of gravity in the local PHOS-module coordinates | |
401 | Float_t wtot = 0; | |
402 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
403 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
404 | Int_t relid[4] ; | |
405 | Float_t xi ; | |
406 | Float_t zi ; | |
407 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
408 | phosgeom->RelPosInModule(relid, xi, zi); | |
409 | Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; | |
410 | x += xi * w ; | |
411 | z += zi * w ; | |
412 | wtot += w ; | |
413 | } | |
414 | x /= wtot ; | |
415 | z /= wtot ; | |
416 | ||
417 | ||
fad3e5b9 | 418 | for(iDigit=0; iDigit < fMulDigit; iDigit++) { |
7932f811 | 419 | digit = (AliPHOSDigit *) ( digits->At(fDigitsList[iDigit]) ) ; |
fad3e5b9 | 420 | Int_t relid[4] ; |
421 | Float_t xi ; | |
422 | Float_t zi ; | |
423 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
424 | phosgeom->RelPosInModule(relid, xi, zi); | |
425 | Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(zi-z)*(zi-z)) ; | |
426 | if(distance < coreRadius) | |
7932f811 | 427 | fCoreEnergy += fEnergyList[iDigit] ; |
fad3e5b9 | 428 | } |
429 | ||
fad3e5b9 | 430 | } |
d15a28e7 | 431 | |
432 | //____________________________________________________________________________ | |
7932f811 | 433 | void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) |
d15a28e7 | 434 | { |
b2a60966 | 435 | // Calculates the axis of the shower ellipsoid |
83974468 | 436 | |
e8dbb96e | 437 | Double_t wtot = 0. ; |
438 | Double_t x = 0.; | |
439 | Double_t z = 0.; | |
440 | Double_t dxx = 0.; | |
441 | Double_t dzz = 0.; | |
442 | Double_t dxz = 0.; | |
d15a28e7 | 443 | |
444 | AliPHOSDigit * digit ; | |
7b7c1533 | 445 | |
446 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; | |
447 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
448 | ||
d15a28e7 | 449 | Int_t iDigit; |
450 | ||
e5b16749 | 451 | |
d15a28e7 | 452 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { |
7932f811 | 453 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; |
d15a28e7 | 454 | Int_t relid[4] ; |
455 | Float_t xi ; | |
456 | Float_t zi ; | |
92862013 | 457 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; |
458 | phosgeom->RelPosInModule(relid, xi, zi); | |
7932f811 | 459 | Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; |
92862013 | 460 | dxx += w * xi * xi ; |
d15a28e7 | 461 | x += w * xi ; |
92862013 | 462 | dzz += w * zi * zi ; |
d15a28e7 | 463 | z += w * zi ; |
92862013 | 464 | dxz += w * xi * zi ; |
d15a28e7 | 465 | wtot += w ; |
466 | } | |
92862013 | 467 | dxx /= wtot ; |
d15a28e7 | 468 | x /= wtot ; |
92862013 | 469 | dxx -= x * x ; |
470 | dzz /= wtot ; | |
d15a28e7 | 471 | z /= wtot ; |
92862013 | 472 | dzz -= z * z ; |
473 | dxz /= wtot ; | |
474 | dxz -= x * z ; | |
d15a28e7 | 475 | |
69183710 | 476 | // //Apply correction due to non-perpendicular incidence |
477 | // Double_t CosX ; | |
478 | // Double_t CosZ ; | |
7b7c1533 | 479 | // AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
480 | // AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
481 | // Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ; | |
69183710 | 482 | |
483 | // CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ; | |
484 | // CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ; | |
485 | ||
486 | // dxx = dxx/(CosX*CosX) ; | |
487 | // dzz = dzz/(CosZ*CosZ) ; | |
488 | // dxz = dxz/(CosX*CosZ) ; | |
489 | ||
490 | ||
7932f811 | 491 | fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; |
492 | if(fLambda[0] > 0) | |
493 | fLambda[0] = TMath::Sqrt(fLambda[0]) ; | |
e8dbb96e | 494 | |
7932f811 | 495 | fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; |
496 | if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda. | |
497 | fLambda[1] = TMath::Sqrt(fLambda[1]) ; | |
e8dbb96e | 498 | else |
7932f811 | 499 | fLambda[1]= 0. ; |
d15a28e7 | 500 | } |
501 | ||
ce2a9a95 | 502 | //____________________________________________________________________________ |
503 | void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) | |
504 | { | |
505 | // Calculate the shower moments in the eigen reference system | |
506 | // M2x, M2z, M3x, M4z | |
507 | // Calculate the angle between the shower position vector and the eigen vector | |
508 | ||
509 | Double_t wtot = 0. ; | |
510 | Double_t x = 0.; | |
511 | Double_t z = 0.; | |
512 | Double_t dxx = 0.; | |
513 | Double_t dzz = 0.; | |
514 | Double_t dxz = 0.; | |
515 | Double_t lambda0=0, lambda1=0; | |
516 | ||
517 | AliPHOSDigit * digit ; | |
518 | ||
519 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; | |
520 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
521 | ||
522 | Int_t iDigit; | |
523 | ||
524 | // 1) Find covariance matrix elements: | |
525 | // || dxx dxz || | |
526 | // || dxz dzz || | |
527 | ||
528 | Float_t xi ; | |
529 | Float_t zi ; | |
530 | Int_t relid[4] ; | |
531 | Double_t w; | |
532 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
533 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
534 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
535 | phosgeom->RelPosInModule(relid, xi, zi); | |
536 | w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; | |
537 | x += w * xi ; | |
538 | z += w * zi ; | |
539 | dxx += w * xi * xi ; | |
540 | dzz += w * zi * zi ; | |
541 | dxz += w * xi * zi ; | |
542 | wtot += w ; | |
543 | } | |
544 | x /= wtot ; | |
545 | z /= wtot ; | |
546 | dxx /= wtot ; | |
547 | dzz /= wtot ; | |
548 | dxz /= wtot ; | |
549 | dxx -= x * x ; | |
550 | dzz -= z * z ; | |
551 | dxz -= x * z ; | |
552 | ||
553 | // 2) Find covariance matrix eigen values lambda0 and lambda1 | |
554 | ||
555 | lambda0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; | |
556 | lambda1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; | |
557 | ||
558 | // 3) Find covariance matrix eigen vectors e0 and e1 | |
559 | ||
560 | TVector2 e0,e1; | |
561 | if (dxz != 0) | |
562 | e0.Set(1.,(lambda0-dxx)/dxz); | |
563 | else | |
564 | e0.Set(0.,1.); | |
565 | ||
566 | e0 = e0.Unit(); | |
567 | e1.Set(-e0.Y(),e0.X()); | |
568 | ||
569 | // 4) Rotate cluster tensor from (x,z) to (e0,e1) system | |
570 | // and calculate moments M3x and M4z | |
571 | ||
572 | Float_t cosPhi = e0.X(); | |
573 | Float_t sinPhi = e0.Y(); | |
574 | ||
575 | Float_t xiPHOS ; | |
576 | Float_t ziPHOS ; | |
577 | Double_t dx3, dz3, dz4; | |
578 | wtot = 0.; | |
579 | x = 0.; | |
580 | z = 0.; | |
581 | dxx = 0.; | |
582 | dzz = 0.; | |
583 | dxz = 0.; | |
584 | dx3 = 0.; | |
585 | dz3 = 0.; | |
586 | dz4 = 0.; | |
587 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
588 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; | |
589 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
590 | phosgeom->RelPosInModule(relid, xiPHOS, ziPHOS); | |
591 | xi = xiPHOS*cosPhi + ziPHOS*sinPhi; | |
592 | zi = ziPHOS*cosPhi - xiPHOS*sinPhi; | |
593 | w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; | |
594 | x += w * xi ; | |
595 | z += w * zi ; | |
596 | dxx += w * xi * xi ; | |
597 | dzz += w * zi * zi ; | |
598 | dxz += w * xi * zi ; | |
599 | dx3 += w * xi * xi * xi; | |
600 | dz3 += w * zi * zi * zi ; | |
601 | dz4 += w * zi * zi * zi * zi ; | |
602 | wtot += w ; | |
603 | } | |
604 | x /= wtot ; | |
605 | z /= wtot ; | |
606 | dxx /= wtot ; | |
607 | dzz /= wtot ; | |
608 | dxz /= wtot ; | |
609 | dx3 /= wtot ; | |
610 | dz3 /= wtot ; | |
611 | dz4 /= wtot ; | |
612 | dx3 += -3*dxx*x + 2*x*x*x; | |
613 | dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z; | |
614 | dxx -= x * x ; | |
615 | dzz -= z * z ; | |
616 | dxz -= x * z ; | |
617 | ||
618 | // 5) Find an angle between cluster center vector and eigen vector e0 | |
619 | ||
620 | Float_t phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z)); | |
621 | ||
622 | fM2x = lambda0; | |
623 | fM2z = lambda1; | |
624 | fM3x = dx3; | |
625 | fM4z = dz4; | |
626 | fPhixe = phi; | |
627 | ||
628 | } | |
b2a60966 | 629 | //____________________________________________________________________________ |
7932f811 | 630 | void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits ) |
ad8cfaf4 | 631 | { |
baef0810 | 632 | // Evaluates all shower parameters |
633 | ||
7932f811 | 634 | EvalLocalPosition(logWeight, digits) ; |
635 | EvalElipsAxis(logWeight, digits) ; | |
ce2a9a95 | 636 | EvalMoments(logWeight, digits) ; |
7932f811 | 637 | EvalDispersion(logWeight, digits) ; |
e5b16749 | 638 | EvalCoreEnergy(logWeight, digits); |
a6eedfad | 639 | EvalTime(digits) ; |
bb53c80c | 640 | AliPHOSRecPoint::EvalAll(logWeight,digits) ; |
ad8cfaf4 | 641 | } |
642 | //____________________________________________________________________________ | |
7932f811 | 643 | void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits) |
b2a60966 | 644 | { |
645 | // Calculates the center of gravity in the local PHOS-module coordinates | |
b2a60966 | 646 | Float_t wtot = 0. ; |
647 | ||
648 | Int_t relid[4] ; | |
649 | ||
650 | Float_t x = 0. ; | |
651 | Float_t z = 0. ; | |
652 | ||
653 | AliPHOSDigit * digit ; | |
654 | ||
7b7c1533 | 655 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
656 | AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); | |
b2a60966 | 657 | |
658 | Int_t iDigit; | |
659 | ||
b2a60966 | 660 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { |
7932f811 | 661 | digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; |
b2a60966 | 662 | |
663 | Float_t xi ; | |
664 | Float_t zi ; | |
665 | phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
666 | phosgeom->RelPosInModule(relid, xi, zi); | |
7932f811 | 667 | Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; |
b2a60966 | 668 | x += xi * w ; |
669 | z += zi * w ; | |
670 | wtot += w ; | |
b2a60966 | 671 | |
9ce1a8d9 | 672 | } |
69183710 | 673 | |
674 | x /= wtot ; | |
675 | z /= wtot ; | |
ad8cfaf4 | 676 | |
677 | // Correction for the depth of the shower starting point (TDR p 127) | |
678 | Float_t para = 0.925 ; | |
679 | Float_t parb = 6.52 ; | |
680 | ||
1b799736 | 681 | Float_t xo,yo,zo ; //Coordinates of the origin |
682 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
683 | ||
5830e1d9 | 684 | Float_t phi = phosgeom->GetPHOSAngle(relid[0]) ; |
1b799736 | 685 | |
686 | //Transform to the local ref.frame | |
687 | Float_t xoL,yoL ; | |
688 | xoL = xo*TMath::Cos(phi)-yo*TMath::Sin(phi) ; | |
689 | yoL = xo*TMath::Sin(phi)+yo*TMath::Cos(phi) ; | |
690 | ||
fd84cb73 | 691 | Float_t radius = phosgeom->GetIPtoCrystalSurface()-yoL; |
1b799736 | 692 | |
693 | Float_t incidencephi = TMath::ATan((x-xoL ) / radius) ; | |
694 | Float_t incidencetheta = TMath::ATan((z-zo) / radius) ; | |
ad8cfaf4 | 695 | |
696 | Float_t depthx = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencephi) ; | |
697 | Float_t depthz = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencetheta) ; | |
ad8cfaf4 | 698 | |
699 | fLocPos.SetX(x - depthx) ; | |
b2a60966 | 700 | fLocPos.SetY(0.) ; |
ad8cfaf4 | 701 | fLocPos.SetZ(z - depthz) ; |
b2a60966 | 702 | |
a6eedfad | 703 | fLocPosM = 0 ; |
b2a60966 | 704 | } |
705 | ||
d15a28e7 | 706 | //____________________________________________________________________________ |
ad8cfaf4 | 707 | Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) const |
d15a28e7 | 708 | { |
b2a60966 | 709 | // Finds the maximum energy in the cluster |
710 | ||
d15a28e7 | 711 | Float_t menergy = 0. ; |
712 | ||
713 | Int_t iDigit; | |
714 | ||
715 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
716 | ||
717 | if(fEnergyList[iDigit] > menergy) | |
718 | menergy = fEnergyList[iDigit] ; | |
719 | } | |
720 | return menergy ; | |
721 | } | |
722 | ||
723 | //____________________________________________________________________________ | |
ad8cfaf4 | 724 | Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const |
d15a28e7 | 725 | { |
b2a60966 | 726 | // Calculates the multiplicity of digits with energy larger than H*energy |
727 | ||
d15a28e7 | 728 | Int_t multipl = 0 ; |
729 | Int_t iDigit ; | |
730 | for(iDigit=0; iDigit<fMulDigit; iDigit++) { | |
731 | ||
732 | if(fEnergyList[iDigit] > H * fAmp) | |
733 | multipl++ ; | |
734 | } | |
735 | return multipl ; | |
736 | } | |
737 | ||
738 | //____________________________________________________________________________ | |
a0636361 | 739 | Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax( AliPHOSDigit ** maxAt, Float_t * maxAtEnergy, |
7932f811 | 740 | Float_t locMaxCut,TClonesArray * digits) const |
d15a28e7 | 741 | { |
b2a60966 | 742 | // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum |
a4e98857 | 743 | // energy difference between two local maxima |
b2a60966 | 744 | |
d15a28e7 | 745 | AliPHOSDigit * digit ; |
746 | AliPHOSDigit * digitN ; | |
747 | ||
748 | ||
749 | Int_t iDigitN ; | |
750 | Int_t iDigit ; | |
751 | ||
7932f811 | 752 | for(iDigit = 0; iDigit < fMulDigit; iDigit++) |
a0636361 | 753 | maxAt[iDigit] = (AliPHOSDigit*) digits->At(fDigitsList[iDigit]) ; |
7932f811 | 754 | |
d15a28e7 | 755 | |
6ad0bfa0 | 756 | for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) { |
a0636361 | 757 | if(maxAt[iDigit]) { |
758 | digit = maxAt[iDigit] ; | |
83974468 | 759 | |
6ad0bfa0 | 760 | for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) { |
7932f811 | 761 | digitN = (AliPHOSDigit *) digits->At(fDigitsList[iDigitN]) ; |
d15a28e7 | 762 | |
9f616d61 | 763 | if ( AreNeighbours(digit, digitN) ) { |
d15a28e7 | 764 | if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) { |
a0636361 | 765 | maxAt[iDigitN] = 0 ; |
6ad0bfa0 | 766 | // but may be digit too is not local max ? |
7932f811 | 767 | if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut) |
a0636361 | 768 | maxAt[iDigit] = 0 ; |
d15a28e7 | 769 | } |
770 | else { | |
a0636361 | 771 | maxAt[iDigit] = 0 ; |
6ad0bfa0 | 772 | // but may be digitN too is not local max ? |
7932f811 | 773 | if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) |
a0636361 | 774 | maxAt[iDigitN] = 0 ; |
d15a28e7 | 775 | } |
776 | } // if Areneighbours | |
777 | } // while digitN | |
778 | } // slot not empty | |
779 | } // while digit | |
780 | ||
781 | iDigitN = 0 ; | |
6ad0bfa0 | 782 | for(iDigit = 0; iDigit < fMulDigit; iDigit++) { |
a0636361 | 783 | if(maxAt[iDigit]){ |
d15a28e7 | 784 | maxAt[iDigitN] = maxAt[iDigit] ; |
9f616d61 | 785 | maxAtEnergy[iDigitN] = fEnergyList[iDigit] ; |
786 | iDigitN++ ; | |
d15a28e7 | 787 | } |
788 | } | |
789 | return iDigitN ; | |
790 | } | |
9688c1dd | 791 | //____________________________________________________________________________ |
0bc3b8ed | 792 | void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits) |
793 | { | |
794 | // Define a rec.point time as a time in the cell with the maximum energy | |
795 | ||
9688c1dd | 796 | Float_t maxE = 0; |
797 | Int_t maxAt = 0; | |
798 | for(Int_t idig=0; idig < fMulDigit; idig++){ | |
799 | if(fEnergyList[idig] > maxE){ | |
800 | maxE = fEnergyList[idig] ; | |
801 | maxAt = idig; | |
802 | } | |
803 | } | |
804 | fTime = ((AliPHOSDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ; | |
805 | ||
806 | } | |
d15a28e7 | 807 | //____________________________________________________________________________ |
092b50ba | 808 | void AliPHOSEmcRecPoint::Purify(Float_t threshold){ |
809 | //Removes digits below threshold | |
810 | ||
811 | Int_t * tempo = new ( Int_t[fMaxDigit] ) ; | |
812 | Float_t * tempoE = new ( Float_t[fMaxDigit] ) ; | |
813 | ||
814 | Int_t mult = 0 ; | |
815 | for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ | |
816 | if(fEnergyList[iDigit] > threshold){ | |
817 | tempo[mult] = fDigitsList[iDigit] ; | |
818 | tempoE[mult] = fEnergyList[iDigit] ; | |
819 | mult++ ; | |
820 | } | |
821 | } | |
822 | ||
823 | fMulDigit = mult ; | |
824 | delete [] fDigitsList ; | |
825 | delete [] fEnergyList ; | |
826 | fDigitsList = new (Int_t[fMulDigit]) ; | |
827 | fEnergyList = new ( Float_t[fMulDigit]) ; | |
828 | ||
829 | for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ | |
830 | fDigitsList[iDigit] = tempo[iDigit]; | |
831 | fEnergyList[iDigit] = tempoE[iDigit] ; | |
832 | } | |
833 | ||
834 | delete [] tempo ; | |
835 | delete [] tempoE ; | |
836 | ||
837 | } | |
838 | //____________________________________________________________________________ | |
839 | void AliPHOSEmcRecPoint::Print(Option_t * option) const | |
d15a28e7 | 840 | { |
b2a60966 | 841 | // Print the list of digits belonging to the cluster |
842 | ||
21cd0c07 | 843 | TString message ; |
844 | message = "AliPHOSEmcRecPoint:\n" ; | |
845 | message += " digits # = " ; | |
846 | Info("Print", message.Data()) ; | |
d15a28e7 | 847 | |
d15a28e7 | 848 | Int_t iDigit; |
7932f811 | 849 | for(iDigit=0; iDigit<fMulDigit; iDigit++) |
53ad38d8 | 850 | printf(" %d ", fDigitsList[iDigit] ) ; |
7932f811 | 851 | |
21cd0c07 | 852 | Info("Print", " Energies = ") ; |
7932f811 | 853 | for(iDigit=0; iDigit<fMulDigit; iDigit++) |
53ad38d8 | 854 | printf(" %f ", fEnergyList[iDigit] ) ; |
855 | printf("\n") ; | |
21cd0c07 | 856 | Info("Print", " Primaries ") ; |
bc68d12c | 857 | for(iDigit = 0;iDigit < fMulTrack; iDigit++) |
53ad38d8 | 858 | printf(" %d ", fTracksList[iDigit]) ; |
859 | printf("\n") ; | |
21cd0c07 | 860 | message = " Multiplicity = %d" ; |
861 | message += " Cluster Energy = %f" ; | |
862 | message += " Number of primaries %d" ; | |
863 | message += " Stored at position %d" ; | |
83974468 | 864 | |
21cd0c07 | 865 | Info("Print", message.Data(), fMulDigit, fAmp, fMulTrack,GetIndexInList() ) ; |
d15a28e7 | 866 | } |
88714635 | 867 | |
7932f811 | 868 |