]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PHOS/AliPHOSCpvRecPoint.cxx
Bag in ::Compare() fixed
[u/mrichter/AliRoot.git] / PHOS / AliPHOSCpvRecPoint.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 /* $Id$ */
17
18 //_________________________________________________________________________
19 //  RecPoint implementation for PHOS-CPV
20 //  An CpvRecPoint is a cluster of digits   
21 //*-- Author: Yuri Kharlov
22 //  (after Dmitri Peressounko (RRC KI & SUBATECH))
23 //  30 October 2000 
24
25 // --- ROOT system ---
26 #include "TPad.h"
27 #include "TH2.h"
28 #include "TMath.h" 
29 #include "TCanvas.h" 
30
31 // --- Standard library ---
32
33 #include <iostream.h> 
34
35 // --- AliRoot header files ---
36
37 #include "AliPHOSGeometry.h"
38 #include "AliPHOSCpvRecPoint.h"
39 #include "AliPHOSPpsdRecPoint.h"
40 #include "AliRun.h"
41 #include "AliPHOSIndexToObject.h"
42
43 ClassImp(AliPHOSCpvRecPoint)
44
45 //____________________________________________________________________________
46 AliPHOSCpvRecPoint::AliPHOSCpvRecPoint(Float_t W0, Float_t LocMaxCut)
47   : AliPHOSRecPoint()
48 {
49   // ctor
50
51   fMulDigit   = 0 ;  
52   fAmp   = 0. ;   
53   fEnergyList = new Float_t[fMaxDigit]; 
54   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
55   fDelta     =  phosgeom->GetCrystalSize(0) ; 
56   fW0        = W0 ;          
57   fLocMaxCut = LocMaxCut ; 
58   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
59   fLengX = -1;
60   fLengZ = -1;
61 }
62
63 //____________________________________________________________________________
64 AliPHOSCpvRecPoint::~AliPHOSCpvRecPoint()
65 {
66   // dtor
67   if ( fEnergyList ) delete[] fEnergyList ; 
68 }
69
70 //____________________________________________________________________________
71 void AliPHOSCpvRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy)
72 {
73   // Adds a digit to the RecPoint
74   //  and accumulates the total amplitude and the multiplicity 
75   
76   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
77     fMaxDigit*=2 ; 
78     Int_t * tempo = new ( Int_t[fMaxDigit] ) ; 
79     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
80
81     Int_t index ;     
82     for ( index = 0 ; index < fMulDigit ; index++ ){
83       tempo[index]  = fDigitsList[index] ;
84       tempoE[index] = fEnergyList[index] ; 
85     }
86     
87     delete [] fDigitsList ; 
88     fDigitsList =  new ( Int_t[fMaxDigit] ) ;
89  
90     delete [] fEnergyList ;
91     fEnergyList =  new ( Float_t[fMaxDigit] ) ;
92
93     for ( index = 0 ; index < fMulDigit ; index++ ){
94       fDigitsList[index] = tempo[index] ;
95       fEnergyList[index] = tempoE[index] ; 
96     }
97  
98     delete [] tempo ;
99     delete [] tempoE ; 
100   } // if
101   
102   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ; 
103   fEnergyList[fMulDigit]   = Energy ;
104   fMulDigit++ ; 
105   fAmp += Energy ; 
106 }
107
108 //____________________________________________________________________________
109 Bool_t AliPHOSCpvRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
110 {
111   // Tells if (true) or not (false) two digits are neighbors)
112   
113   Bool_t aren = kFALSE ;
114   
115   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
116   Int_t relid1[4] ; 
117   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
118
119   Int_t relid2[4] ; 
120   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
121   
122   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
123   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
124
125   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
126     aren = kTRUE ;
127   
128   return aren ;
129 }
130
131 //____________________________________________________________________________
132 Int_t AliPHOSCpvRecPoint::Compare(TObject * obj)
133 {
134   // Compares two RecPoints according to their position in the PHOS modules
135
136   Int_t rv ; 
137
138   if( (strcmp(obj->ClassName() , "AliPHOSPpsdRecPoint" )) == 0)  // PPSD Rec Point
139     {
140       AliPHOSPpsdRecPoint * clu = (AliPHOSPpsdRecPoint *)obj ; 
141       if(this->GetPHOSMod()  < clu->GetPHOSMod() ) 
142         rv = -1 ;
143       else 
144         rv = 1 ;
145       return rv ;
146     }
147   else
148     {
149       AliPHOSCpvRecPoint * clu  = (AliPHOSCpvRecPoint *) obj ; 
150       
151       Int_t phosmod1 = this->GetPHOSMod() ;
152       Int_t phosmod2 = clu->GetPHOSMod() ;
153       
154       TVector3 locpos1; 
155       this->GetLocalPosition(locpos1) ;
156       TVector3 locpos2;  
157       clu->GetLocalPosition(locpos2) ;  
158       
159       if(phosmod1 == phosmod2 ) {
160         Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
161         if (rowdif> 0) 
162           rv = -1 ;
163         else if(rowdif < 0) 
164           rv = 1 ;
165         else if(locpos1.Z()>locpos2.Z()) 
166           rv = -1 ;
167         else 
168           rv = 1 ; 
169       }
170       
171       else {
172         if(phosmod1 < phosmod2 ) 
173           rv = -1 ;
174         else 
175           rv = 1 ;
176       }
177       
178       return rv ; 
179     }
180 }
181
182 //______________________________________________________________________________
183 void AliPHOSCpvRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
184 {
185   // Execute action corresponding to one event
186   //  This member function is called when a AliPHOSRecPoint is clicked with the locator
187   //
188   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
189   //  and switched off when the mouse button is released.
190   //
191
192   //   static Int_t pxold, pyold;
193
194   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
195   
196   static TGraph *  digitgraph = 0 ;
197   
198   if (!gPad->IsEditable()) return;
199   
200   TH2F * histo = 0 ;
201   TCanvas * histocanvas ; 
202   
203   switch (event) {
204     
205   case kButton1Down: {
206     AliPHOSDigit * digit ;
207     AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
208     Int_t iDigit;
209     Int_t relid[4] ;
210     
211     const Int_t kMulDigit = AliPHOSCpvRecPoint::GetDigitsMultiplicity() ; 
212     Float_t * xi = new Float_t[kMulDigit] ; 
213     Float_t * zi = new Float_t[kMulDigit] ; 
214     
215     // create the histogram for the single cluster 
216     // 1. gets histogram boundaries
217     Float_t ximax = -999. ; 
218     Float_t zimax = -999. ; 
219     Float_t ximin = 999. ; 
220     Float_t zimin = 999. ;
221     
222     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
223       digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
224       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
225       phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
226       if ( xi[iDigit] > ximax )
227         ximax = xi[iDigit] ; 
228       if ( xi[iDigit] < ximin )
229         ximin = xi[iDigit] ; 
230       if ( zi[iDigit] > zimax )
231         zimax = zi[iDigit] ; 
232       if ( zi[iDigit] < zimin )
233         zimin = zi[iDigit] ;     
234     }
235     ximax += phosgeom->GetCrystalSize(0) / 2. ;
236     zimax += phosgeom->GetCrystalSize(2) / 2. ;
237     ximin -= phosgeom->GetCrystalSize(0) / 2. ;
238     zimin -= phosgeom->GetCrystalSize(2) / 2. ;
239     Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5  ) ; 
240     Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
241     
242     // 2. gets the histogram title
243     
244     Text_t title[100] ; 
245     sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
246     
247     if (!histo) {
248       delete histo ; 
249       histo = 0 ; 
250     }
251     histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
252     
253     Float_t x, z ; 
254     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
255       digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
256       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
257       phosgeom->RelPosInModule(relid, x, z);
258       histo->Fill(x, z, fEnergyList[iDigit] ) ;
259     }
260     
261     if (!digitgraph) {
262       digitgraph = new TGraph(kMulDigit,xi,zi);
263       digitgraph-> SetMarkerStyle(5) ; 
264       digitgraph-> SetMarkerSize(1.) ;
265       digitgraph-> SetMarkerColor(1) ;
266       digitgraph-> Paint("P") ;
267     }
268     
269     Print() ;
270     histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
271     histocanvas->Draw() ; 
272     histo->Draw("lego1") ; 
273     
274     delete[] xi ; 
275     delete[] zi ; 
276     
277     break;
278   }
279   
280   case kButton1Up: 
281     if (digitgraph) {
282       delete digitgraph  ;
283       digitgraph = 0 ;
284     }
285     break;
286   
287    }
288 }
289
290 //____________________________________________________________________________
291 Float_t  AliPHOSCpvRecPoint::GetDispersion() 
292 {
293   // Calculates the dispersion of the shower at the origine of the RecPoint
294
295   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
296
297   Float_t d    = 0 ;
298   Float_t wtot = 0 ;
299
300   TVector3 locpos;
301   GetLocalPosition(locpos);
302   Float_t x = locpos.X() ;
303   Float_t z = locpos.Z() ;
304
305   AliPHOSDigit * digit ;
306   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
307   
308   Int_t iDigit;
309   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
310     digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
311     Int_t relid[4] ;
312     Float_t xi ;
313     Float_t zi ;
314     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
315     phosgeom->RelPosInModule(relid, xi, zi);
316     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
317     d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
318     wtot+=w ;
319   }
320
321   d /= wtot ;
322
323   return TMath::Sqrt(d) ;
324 }
325
326 //____________________________________________________________________________
327 void  AliPHOSCpvRecPoint::GetElipsAxis(Float_t * lambda)
328 {
329   // Calculates the axis of the shower ellipsoid
330
331   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
332
333   Float_t wtot = 0. ;
334   Float_t x    = 0.;
335   Float_t z    = 0.;
336   Float_t dxx  = 0.;
337   Float_t dzz  = 0.;
338   Float_t dxz  = 0.;
339
340   AliPHOSDigit * digit ;
341   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
342   Int_t iDigit;
343
344   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
345     digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
346     Int_t relid[4] ;
347     Float_t xi ;
348     Float_t zi ;
349     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
350     phosgeom->RelPosInModule(relid, xi, zi);
351     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
352     dxx  += w * xi * xi ;
353     x    += w * xi ;
354     dzz  += w * zi * zi ;
355     z    += w * zi ; 
356     dxz  += w * xi * zi ; 
357     wtot += w ;
358   }
359
360   dxx /= wtot ;
361   x   /= wtot ;
362   dxx -= x * x ;
363   dzz /= wtot ;
364   z   /= wtot ;
365   dzz -= z * z ;
366   dxz /= wtot ;
367   dxz -= x * z ;
368
369   Float_t tmp0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz );
370   Float_t tmp1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz );
371   if (tmp0>=0) lambda[0] = TMath::Sqrt(tmp0);
372   else lambda[0] = -TMath::Sqrt(-tmp0);
373   if (tmp1>=0) lambda[1] = TMath::Sqrt(tmp1);
374   else lambda[1] = -TMath::Sqrt(-tmp1);
375 //    lambda[0] = TMath::Sqrt( 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
376 //    lambda[1] = TMath::Sqrt( 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
377 }
378
379 //____________________________________________________________________________
380 void AliPHOSCpvRecPoint::GetLocalPosition(TVector3 &LPos)
381 {
382   // Calculates the center of gravity in the local PHOS-module coordinates 
383
384   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
385
386   if( fLocPos.X() < 1000000.) { // already evaluated
387    LPos = fLocPos ;
388    return ;
389   }
390
391   Float_t wtot = 0. ;
392  
393   Int_t relid[4] ;
394
395   Float_t x = 0. ;
396   Float_t z = 0. ;
397   
398   AliPHOSDigit * digit ;
399
400   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
401
402   Int_t iDigit;
403
404   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
405     digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) );
406
407     Float_t xi ;
408     Float_t zi ;
409     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
410     phosgeom->RelPosInModule(relid, xi, zi);
411     Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
412     x    += xi * w ;
413     z    += zi * w ;
414     wtot += w ;
415   }
416
417   if (wtot != 0) {
418     x /= wtot ;
419     z /= wtot ;
420   } else {
421     x = -1e6 ;
422     z = -1e6 ;
423     if (fMulDigit != 0) cout << "AliPHOSCpvRecPoint: too low log weight factor "
424                              << "to evaluate cluster's center\n";
425   }
426   fLocPos.SetX(x)  ;
427   fLocPos.SetY(0.) ;
428   fLocPos.SetZ(z)  ;
429
430   LPos = fLocPos ;
431 }
432
433 //____________________________________________________________________________
434 void AliPHOSCpvRecPoint::GetClusterLengths(Int_t &lengX, Int_t &lengZ)
435 {
436   // Calculates the cluster lengths along X and Z axes
437   // These characteristics are needed for CPV to tune
438   // digitization+reconstruction to experimental data
439   // Yuri Kharlov. 24 October 2000
440
441   if( fLengX != -1 && fLengZ != -1 ) { // already evaluated
442     lengX = fLengX ;
443     lengZ = fLengZ ;
444     return ;
445   }
446
447   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
448
449   Int_t relid[4] ;
450
451   AliPHOSDigit * digit ;
452
453   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
454
455   const Int_t kMaxLeng=20;
456   Int_t idX[kMaxLeng], idZ[kMaxLeng];
457   lengX = 0;
458   lengZ = 0;
459   Bool_t dejavu;
460
461   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
462     digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) );
463     Int_t absId = digit->GetId();
464     phosgeom->AbsToRelNumbering(absId, relid) ;
465
466     Int_t i;
467     dejavu=kFALSE;
468     for (i=0; i<lengX; i++) if (relid[2]==idX[i]) { dejavu=kTRUE; break; }
469     if (!dejavu) {
470       idX[lengX]=relid[2];
471       lengX++;
472       lengX = TMath::Min(lengX,kMaxLeng);
473     }
474
475     dejavu=kFALSE;
476     for (i=0; i<lengZ; i++) if (relid[3]==idZ[i]) { dejavu=kTRUE; break; }
477     if (!dejavu) {
478       idZ[lengZ]=relid[3];
479       lengZ++;
480       lengZ = TMath::Min(lengZ,kMaxLeng);
481     }
482   }
483   fLengX = lengX;
484   fLengZ = lengZ;
485 }
486
487 //____________________________________________________________________________
488 Float_t AliPHOSCpvRecPoint::GetMaximalEnergy(void)
489 {
490   // Finds the maximum energy in the cluster
491   
492   Float_t menergy = 0. ;
493
494   Int_t iDigit;
495
496   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
497  
498     if(fEnergyList[iDigit] > menergy) 
499       menergy = fEnergyList[iDigit] ;
500   }
501   return menergy ;
502 }
503
504 //____________________________________________________________________________
505 Int_t AliPHOSCpvRecPoint::GetMultiplicityAtLevel(const Float_t H) 
506 {
507   // Calculates the multiplicity of digits with energy larger than H*energy 
508   
509   Int_t multipl   = 0 ;
510   Int_t iDigit ;
511   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
512
513     if(fEnergyList[iDigit] > H * fAmp) 
514       multipl++ ;
515   }
516   return multipl ;
517 }
518
519 //____________________________________________________________________________
520 Int_t  AliPHOSCpvRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
521
522   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
523   //  energy difference between two local maxima
524
525   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
526
527   AliPHOSDigit * digit ;
528   AliPHOSDigit * digitN ;
529   
530
531   Int_t iDigitN ;
532   Int_t iDigit ;
533
534   for(iDigit = 0; iDigit < fMulDigit; iDigit++){
535     maxAt[iDigit] = (Int_t) ( please->GimeDigit(fDigitsList[iDigit]) ) ;
536   }
537   
538   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
539     if(maxAt[iDigit] != -1) {
540       digit = (AliPHOSDigit *) maxAt[iDigit] ;
541           
542       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {        
543         digitN = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigitN]) ) ; 
544         
545         if ( AreNeighbours(digit, digitN) ) {
546           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
547             maxAt[iDigitN] = -1 ;
548             // but may be digit too is not local max ?
549             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) 
550               maxAt[iDigit] = -1 ;
551           }
552           else {
553             maxAt[iDigit] = -1 ;
554             // but may be digitN too is not local max ?
555             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) 
556               maxAt[iDigitN] = -1 ; 
557           } 
558         } // if Areneighbours
559       } // while digitN
560     } // slot not empty
561   } // while digit
562   
563   iDigitN = 0 ;
564   for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
565     if(maxAt[iDigit] != -1){
566       maxAt[iDigitN] = maxAt[iDigit] ;
567       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
568       iDigitN++ ; 
569     }
570   }
571   return iDigitN ;
572 }
573
574
575 //____________________________________________________________________________
576 void AliPHOSCpvRecPoint::Print(Option_t * option) 
577 {
578   // Print the list of digits belonging to the cluster
579   
580   cout << "AliPHOSCpvRecPoint: " << endl ;
581
582   AliPHOSDigit * digit ; 
583   Int_t iDigit;
584   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
585
586   Float_t xi ;
587   Float_t zi ;
588   Int_t relid[4] ; 
589   
590   AliPHOSIndexToObject * please =  AliPHOSIndexToObject::GetInstance() ; 
591
592   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
593     digit = please->GimeDigit( fDigitsList[iDigit] ) ; 
594     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
595     phosgeom->RelPosInModule(relid, xi, zi);
596     cout << " Id = " << digit->GetId() ;  
597     cout << "   module  = " << relid[0] ;  
598     cout << "   x  = " << xi ;  
599     cout << "   z  = " << zi ;  
600     cout << "   Energy = " << fEnergyList[iDigit] << endl ;
601   }
602   cout << "       Multiplicity    = " << fMulDigit  << endl ;
603   cout << "       Cluster Energy  = " << fAmp << endl ;
604   cout << "       Stored at position " << GetIndexInList() << endl ; 
605  
606 }