df8047056ea31f32670eb6eab24a250172bcd42f
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Rec Point in the PHOS EM calorimeter 
18 //*-- Author : Dmitri Peressounko RRC KI
19 //////////////////////////////////////////////////////////////////////////////
20
21 // --- ROOT system ---
22 #include "TPad.h"
23 #include "TH2.h"
24 #include "TMath.h" 
25 #include "TCanvas.h" 
26
27 // --- Standard library ---
28
29 #include <iostream>
30
31 // --- AliRoot header files ---
32
33 #include "AliPHOSGeometry.h"
34 #include "AliPHOSEmcRecPoint.h"
35 #include "AliRun.h"
36
37 ClassImp(AliPHOSEmcRecPoint)
38
39 //____________________________________________________________________________
40 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
41   : AliPHOSRecPoint()
42 {
43   // ctor
44
45   fMulDigit   = 0 ;  
46   fAmp   = 0. ;   
47   fEnergyList = new Float_t[fMaxDigit]; 
48   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
49   fDelta     =  phosgeom->GetCrystalSize(0) ; 
50   fW0        = W0 ;          
51   fLocMaxCut = LocMaxCut ; 
52   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
53 }
54
55 //____________________________________________________________________________
56 AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() 
57 {
58   // dtor 
59 }
60
61 //____________________________________________________________________________
62 void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
63 {
64   // adds a digit to the digits list
65   // and accumulates the total amplitude and the multiplicity 
66   
67   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
68     fMaxDigit*=2 ; 
69     int * tempo = new ( int[fMaxDigit] ) ; 
70     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
71
72     Int_t index ;     
73     for ( index = 0 ; index < fMulDigit ; index++ ){
74       tempo[index] = fDigitsList[index] ;
75       tempoE[index] = fEnergyList[index] ; 
76     }
77     
78     delete [] fDigitsList ; 
79     fDigitsList =  new ( int[fMaxDigit] ) ;
80  
81     delete [] fEnergyList ;
82     fEnergyList =  new ( Float_t[fMaxDigit] ) ;
83
84     for ( index = 0 ; index < fMulDigit ; index++ ){
85       fDigitsList[index] = tempo[index] ;
86       fEnergyList[index] = tempoE[index] ; 
87     }
88  
89     delete [] tempo ;
90     delete [] tempoE ; 
91   } // if
92   
93   fDigitsList[fMulDigit]   =  (int) &digit  ; 
94   fEnergyList[fMulDigit++] = Energy ;
95   fAmp += Energy ; 
96 }
97
98 //____________________________________________________________________________
99 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
100 {
101   
102   Bool_t aren = kFALSE ;
103   
104   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
105   Int_t relid1[4] ; 
106   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
107
108   Int_t relid2[4] ; 
109   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
110   
111   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
112   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
113
114   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
115     aren = kTRUE ;
116   
117   return aren ;
118 }
119
120 //____________________________________________________________________________
121 Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
122 {
123   Int_t rv ; 
124
125   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; 
126
127  
128   Int_t phosmod1 = this->GetPHOSMod() ;
129   Int_t phosmod2 = clu->GetPHOSMod() ;
130
131   TVector3 locpos1; 
132   this->GetLocalPosition(locpos1) ;
133   TVector3 locpos2;  
134   clu->GetLocalPosition(locpos2) ;  
135
136   if(phosmod1 == phosmod2 ) {
137     Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
138     if (rowdif> 0) 
139       rv = -1 ;
140      else if(rowdif < 0) 
141        rv = 1 ;
142     else if(locpos1.Z()>locpos2.Z()) 
143       rv = -1 ;
144     else 
145       rv = 1 ; 
146      }
147
148   else {
149     if(phosmod1 < phosmod2 ) 
150       rv = -1 ;
151     else 
152       rv = 1 ;
153   }
154
155   return rv ; 
156 }
157
158 //______________________________________________________________________________
159 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
160 {
161   //Execute action corresponding to one event
162   //  This member function is called when a AliPHOSRecPoint is clicked with the locator
163   //
164   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
165   //  and switched off when the mouse button is released.
166   //
167
168   //   static Int_t pxold, pyold;
169
170    static TGraph *  digitgraph = 0 ;
171
172    if (!gPad->IsEditable()) return;
173
174    TH2F * histo = 0 ;
175    TCanvas * histocanvas ; 
176    
177    switch (event) {
178    
179    case kButton1Down: {
180      AliPHOSDigit * digit ;
181      AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
182      Int_t iDigit;
183      Int_t relid[4] ;
184      Float_t xi[fMulDigit] ;
185      Float_t zi[fMulDigit] ;
186
187      // create the histogram for the single cluster 
188      // 1. gets histogram boundaries
189      Float_t ximax = -999. ; 
190      Float_t zimax = -999. ; 
191      Float_t ximin = 999. ; 
192      Float_t zimin = 999. ;
193  
194      for(iDigit=0; iDigit<fMulDigit; iDigit++) {
195        digit = (AliPHOSDigit *) fDigitsList[iDigit];
196        phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
197        phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
198        if ( xi[iDigit] > ximax )
199          ximax = xi[iDigit] ; 
200        if ( xi[iDigit] < ximin )
201          ximin = xi[iDigit] ; 
202        if ( zi[iDigit] > zimax )
203          zimax = zi[iDigit] ; 
204        if ( zi[iDigit] < zimin )
205          zimin = zi[iDigit] ;     
206      }
207      ximax += phosgeom->GetCrystalSize(0) / 2. ;
208      zimax += phosgeom->GetCrystalSize(2) / 2. ;
209      ximin -= phosgeom->GetCrystalSize(0) / 2. ;
210      zimin -= phosgeom->GetCrystalSize(2) / 2. ;
211      Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5  ) ; 
212      Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
213  
214      // 2. gets the histogram title
215
216      Text_t title[100] ; 
217      sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
218   
219      if (!histo) {
220        delete histo ; 
221        histo = 0 ; 
222      }
223      histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
224
225      Float_t x, z ; 
226      for(iDigit=0; iDigit<fMulDigit; iDigit++) {
227        digit = (AliPHOSDigit *) fDigitsList[iDigit];
228        phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
229        phosgeom->RelPosInModule(relid, x, z);
230        histo->Fill(x, z, fEnergyList[iDigit] ) ;
231      }
232
233      if (!digitgraph) {
234        digitgraph = new TGraph(fMulDigit,xi,zi);
235        digitgraph-> SetMarkerStyle(5) ; 
236        digitgraph-> SetMarkerSize(1.) ;
237        digitgraph-> SetMarkerColor(1) ;
238        digitgraph-> Paint("P") ;
239      }
240
241      Print() ;
242      histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
243      histocanvas->Draw() ; 
244      histo->Draw("lego1") ; 
245
246      break;
247    }
248
249    case kButton1Up: 
250      if (digitgraph) {
251        delete digitgraph  ;
252        digitgraph = 0 ;
253      }
254      break;
255   
256    }
257 }
258
259 //____________________________________________________________________________
260 Float_t  AliPHOSEmcRecPoint::GetDispersion() 
261 {
262   Float_t d    = 0 ;
263   Float_t wtot = 0 ;
264
265   TVector3 locpos;
266   GetLocalPosition(locpos);
267   Float_t x = locpos.X() ;
268   Float_t z = locpos.Z() ;
269   //  Int_t i = GetPHOSMod() ;
270
271   AliPHOSDigit * digit ;
272   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
273   
274   Int_t iDigit;
275   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
276     digit = (AliPHOSDigit *) fDigitsList[iDigit];
277     Int_t relid[4] ;
278     Float_t xi ;
279     Float_t zi ;
280     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
281     phosgeom->RelPosInModule(relid, xi, zi);
282     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
283     d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
284     wtot+=w ;
285   }
286
287   d /= wtot ;
288
289   return TMath::Sqrt(d) ;
290 }
291
292 //____________________________________________________________________________
293 void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
294 {
295   Float_t wtot = 0. ;
296   Float_t x    = 0.;
297   Float_t z    = 0.;
298   Float_t dxx  = 0.;
299   Float_t dzz  = 0.;
300   Float_t dxz  = 0.;
301
302   AliPHOSDigit * digit ;
303   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
304   Int_t iDigit;
305
306   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
307     digit = (AliPHOSDigit *) fDigitsList[iDigit];
308     Int_t relid[4] ;
309     Float_t xi ;
310     Float_t zi ;
311     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
312     phosgeom->RelPosInModule(relid, xi, zi);
313     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
314     dxx  += w * xi * xi ;
315     x    += w * xi ;
316     dzz  += w * zi * zi ;
317     z    += w * zi ; 
318     dxz  += w * xi * zi ; 
319     wtot += w ;
320   }
321
322   dxx /= wtot ;
323   x   /= wtot ;
324   dxx -= x * x ;
325   dzz /= wtot ;
326   z   /= wtot ;
327   dzz -= z * z ;
328   dxz /= wtot ;
329   dxz -= x * z ;
330
331   lambda[0] = TMath::Sqrt( 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
332   lambda[1] = TMath::Sqrt( 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
333 }
334
335 //____________________________________________________________________________
336 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
337 {
338   Float_t menergy = 0. ;
339
340   Int_t iDigit;
341
342   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
343  
344     if(fEnergyList[iDigit] > menergy) 
345       menergy = fEnergyList[iDigit] ;
346   }
347   return menergy ;
348 }
349
350 //____________________________________________________________________________
351 Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) 
352 {
353   Int_t multipl   = 0 ;
354   Int_t iDigit ;
355   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
356
357     if(fEnergyList[iDigit] > H * fAmp) 
358       multipl++ ;
359   }
360   return multipl ;
361 }
362
363 //____________________________________________________________________________
364 Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
365
366   AliPHOSDigit * digit ;
367   AliPHOSDigit * digitN ;
368   
369
370   Int_t iDigitN ;
371   Int_t iDigit ;
372
373   for(iDigit = 0; iDigit < fMulDigit; iDigit++){
374     maxAt[iDigit] = fDigitsList[iDigit] ;
375   }
376   
377   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
378     if(maxAt[iDigit] != -1) {
379       digit = (AliPHOSDigit *) maxAt[iDigit] ;
380          
381       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {        
382         digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ; 
383         
384         if ( AreNeighbours(digit, digitN) ) {
385           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
386             maxAt[iDigitN] = -1 ;
387             // but may be digit too is not local max ?
388             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) 
389               maxAt[iDigit] = -1 ;
390           }
391           else {
392             maxAt[iDigit] = -1 ;
393             // but may be digitN too is not local max ?
394             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) 
395               maxAt[iDigitN] = -1 ; 
396           } 
397         } // if Areneighbours
398       } // while digitN
399     } // slot not empty
400   } // while digit
401   
402   iDigitN = 0 ;
403   for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
404     if(maxAt[iDigit] != -1){
405       maxAt[iDigitN] = maxAt[iDigit] ;
406       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
407       iDigitN++ ; 
408     }
409   }
410   return iDigitN ;
411 }
412
413 //____________________________________________________________________________
414 void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
415 {
416   if( fLocPos.X() < 1000000.) { // already evaluated
417    LPos = fLocPos ;
418    return ;
419   }
420
421   Float_t wtot = 0. ;
422  
423   Int_t relid[4] ;
424
425   Float_t x = 0. ;
426   Float_t z = 0. ;
427   
428   AliPHOSDigit * digit ;
429
430   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
431
432   Int_t iDigit;
433
434
435   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
436     digit = (AliPHOSDigit *) fDigitsList[iDigit];
437
438     Float_t xi ;
439     Float_t zi ;
440     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
441     phosgeom->RelPosInModule(relid, xi, zi);
442     Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
443     x    += xi * w ;
444     z    += zi * w ;
445     wtot += w ;
446
447   }
448
449   x /= wtot ;
450   z /= wtot ;
451   fLocPos.SetX(x)  ;
452   fLocPos.SetY(0.) ;
453   fLocPos.SetZ(z)  ;
454
455   LPos = fLocPos ;
456 }
457
458 // //____________________________________________________________________________
459 // AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) 
460 // {
461 //   int * dl = Clu.GetDigitsList() ; 
462  
463 //  if(fDigitsList) 
464 //     delete fDigitsList  ;
465
466 //   AliPHOSDigit * digit ;
467  
468 //   Int_t iDigit;
469
470 //   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
471 //     digit = (AliPHOSDigit *) dl[iDigit];
472 //     AddDigit(*digit) ;
473 //   }
474
475 //   fAmp       = Clu.GetTotalEnergy() ;
476 //   fGeom      = Clu.GetGeom() ;
477 //   TVector3 locpos;
478 //   Clu.GetLocalPosition(locpos) ;
479 //   fLocPos    = locpos;
480 //   fMulDigit       = Clu.GetMultiplicity() ;
481 //   fMaxDigit  = Clu.GetMaximumMultiplicity() ;
482 //   fPHOSMod   = Clu.GetPHOSMod() ;
483 //   fW0        = Clu.GetLogWeightCut() ; 
484 //   fDelta     = Clu.GetDelta() ;
485 //   fLocMaxCut = Clu.GetLocMaxCut() ;
486   
487 //   delete dl ; 
488  
489 //   return *this ;
490 // }
491
492 //____________________________________________________________________________
493 void AliPHOSEmcRecPoint::Print(Option_t * option) 
494 {
495   cout << "AliPHOSEmcRecPoint: " << endl ;
496
497   AliPHOSDigit * digit ; 
498   Int_t iDigit;
499   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
500
501   Float_t xi ;
502   Float_t zi ;
503   Int_t relid[4] ; 
504  
505   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
506     digit = (AliPHOSDigit *) fDigitsList[iDigit];
507     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
508     phosgeom->RelPosInModule(relid, xi, zi);
509     cout << " Id = " << digit->GetId() ;  
510     cout << "   module  = " << relid[0] ;  
511     cout << "   x  = " << xi ;  
512     cout << "   z  = " << zi ;  
513     cout << "   Energy = " << fEnergyList[iDigit] << endl ;
514   }
515   cout << "       Multiplicity    = " << fMulDigit  << endl ;
516   cout << "       Cluster Energy  = " << fAmp << endl ;
517   
518 }
519
520 //______________________________________________________________________________
521 void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
522 {
523    // Stream an object of class AliPHOSEmcRecPoint.
524
525    if (R__b.IsReading()) {
526       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
527       AliPHOSRecPoint::Streamer(R__b);
528       R__b >> fDelta;
529       R__b >> fLocMaxCut;
530       R__b.ReadArray(fEnergyList);
531       R__b >> fW0;
532    } else {
533       R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
534       AliPHOSRecPoint::Streamer(R__b);
535       R__b << fDelta;
536       R__b << fLocMaxCut;
537       R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
538       R__b << fW0;
539    }
540 }
541
542