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