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