Inserting TMath.h where required by the new version of ROOT
[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 /* History of cvs commits:
19  *
20  * $Log$
21  * Revision 1.23  2005/12/20 14:28:47  hristov
22  * Additional protection
23  *
24  * Revision 1.22  2005/05/28 14:19:04  schutz
25  * Compilation warnings fixed by T.P.
26  *
27  */
28
29 //_________________________________________________________________________
30 //  RecPoint implementation for PHOS-CPV
31 //  An CpvRecPoint is a cluster of digits   
32 //*-- Author: Yuri Kharlov
33 //  (after Dmitri Peressounko (RRC KI & SUBATECH))
34 //  30 October 2000 
35
36 // --- ROOT system ---
37
38 #include "TMath.h" 
39 #include "TClonesArray.h" 
40
41 // --- Standard library ---
42
43 // --- AliRoot header files ---
44 #include "AliLog.h"
45 #include "AliPHOSGeometry.h" 
46 #include "AliPHOSDigit.h"
47 #include "AliPHOSCpvRecPoint.h"
48 #include "AliPHOSLoader.h"
49
50 ClassImp(AliPHOSCpvRecPoint)
51
52 //____________________________________________________________________________
53 AliPHOSCpvRecPoint::AliPHOSCpvRecPoint() : 
54   AliPHOSEmcRecPoint(),
55   fLengX(-1),
56   fLengZ(-1)
57 {
58   // ctor
59 }
60
61 //____________________________________________________________________________
62 AliPHOSCpvRecPoint::AliPHOSCpvRecPoint(const char * opt) : 
63   AliPHOSEmcRecPoint(opt),
64   fLengX(-1),
65   fLengZ(-1)
66 {
67    // ctor
68 }
69
70 //____________________________________________________________________________
71 AliPHOSCpvRecPoint::~AliPHOSCpvRecPoint()
72 {
73   // dtor
74 }
75
76 //____________________________________________________________________________
77 Bool_t AliPHOSCpvRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const
78 {
79   // Tells if (true) or not (false) two digits are neighbors)
80   
81   Bool_t aren = kFALSE ;
82   
83   AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
84
85   Int_t relid1[4] ; 
86   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
87
88   Int_t relid2[4] ; 
89   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
90   
91   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
92   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
93
94   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
95     aren = kTRUE ;
96   
97   return aren ;
98 }
99
100 //____________________________________________________________________________
101 Int_t AliPHOSCpvRecPoint::Compare(const TObject * obj) const
102 {
103   // Compares two RecPoints according to their position in the PHOS modules
104
105   Float_t delta = 1 ; //Width of "Sorting row". If you changibg this 
106                       //value (what is senseless) change as vell delta in
107                       //AliPHOSTrackSegmentMakerv* and other RecPoints...
108
109   Int_t rv ; 
110
111   AliPHOSCpvRecPoint * clu  = (AliPHOSCpvRecPoint *) obj ; 
112   
113   Int_t phosmod1 = GetPHOSMod() ;
114   Int_t phosmod2 = clu->GetPHOSMod() ;
115   
116   TVector3 locpos1; 
117   GetLocalPosition(locpos1) ;
118   TVector3 locpos2;  
119   clu->GetLocalPosition(locpos2) ;  
120   
121   if(phosmod1 == phosmod2 ) {
122     Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
123     if (rowdif> 0) 
124       rv = 1 ;
125     else if(rowdif < 0) 
126       rv = -1 ;
127     else if(locpos1.Z()>locpos2.Z()) 
128       rv = -1 ;
129     else 
130       rv = 1 ; 
131   }
132   
133   else {
134     if(phosmod1 < phosmod2 ) 
135       rv = -1 ;
136     else 
137       rv = 1 ;
138   }
139   
140   return rv ; 
141
142 }
143
144 //______________________________________________________________________________
145 void AliPHOSCpvRecPoint::ExecuteEvent(Int_t, Int_t, Int_t ) /*const*/
146 {
147 //   // Execute action corresponding to one event
148 //   //  This member function is called when a AliPHOSRecPoint is clicked with the locator
149 //   //
150 //   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
151 //   //  and switched off when the mouse button is released.
152 //   //
153
154 //   //   static Int_t pxold, pyold;
155
156 //   AliPHOSLoader * gime =  AliPHOSLoader::GetInstance() ; 
157   
158 //   static TGraph *  digitgraph = 0 ;
159   
160 //   if (!gPad->IsEditable()) return;
161   
162 //   TH2F * histo = 0 ;
163 //   TCanvas * histocanvas ; 
164   
165 //   switch (event) {
166     
167 //   case kButton1Down: {
168 //     AliPHOSDigit * digit ;
169 //     AliPHOSGeometry * phosgeom =  (AliPHOSGeometry *) fGeom ;
170 //     Int_t iDigit;
171 //     Int_t relid[4] ;
172     
173 //     const Int_t kMulDigit = AliPHOSCpvRecPoint::GetDigitsMultiplicity() ; 
174 //     Float_t * xi = new Float_t[kMulDigit] ; 
175 //     Float_t * zi = new Float_t[kMulDigit] ; 
176     
177 //     // create the histogram for the single cluster 
178 //     // 1. gets histogram boundaries
179 //     Float_t ximax = -999. ; 
180 //     Float_t zimax = -999. ; 
181 //     Float_t ximin = 999. ; 
182 //     Float_t zimin = 999. ;
183     
184 //     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
185 //       digit = (AliPHOSDigit *) ( gime->Digit(fDigitsList[iDigit]) ) ;
186 //       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
187 //       phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
188 //       if ( xi[iDigit] > ximax )
189 //      ximax = xi[iDigit] ; 
190 //       if ( xi[iDigit] < ximin )
191 //      ximin = xi[iDigit] ; 
192 //       if ( zi[iDigit] > zimax )
193 //      zimax = zi[iDigit] ; 
194 //       if ( zi[iDigit] < zimin )
195 //      zimin = zi[iDigit] ;     
196 //     }
197 //     ximax += phosgeom->GetCrystalSize(0) / 2. ;
198 //     zimax += phosgeom->GetCrystalSize(2) / 2. ;
199 //     ximin -= phosgeom->GetCrystalSize(0) / 2. ;
200 //     zimin -= phosgeom->GetCrystalSize(2) / 2. ;
201 //     Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5  ) ; 
202 //     Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
203     
204 //     // 2. gets the histogram title
205     
206 //     Text_t title[100] ; 
207 //     sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
208     
209 //     if (!histo) {
210 //       delete histo ; 
211 //       histo = 0 ; 
212 //     }
213 //     histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
214     
215 //     Float_t x, z ; 
216 //     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
217 //       digit = (AliPHOSDigit *) ( gime->Digit(fDigitsList[iDigit]) ) ;
218 //       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
219 //       phosgeom->RelPosInModule(relid, x, z);
220 //       histo->Fill(x, z, fEnergyList[iDigit] ) ;
221 //     }
222     
223 //     if (!digitgraph) {
224 //       digitgraph = new TGraph(kMulDigit,xi,zi);
225 //       digitgraph-> SetMarkerStyle(5) ; 
226 //       digitgraph-> SetMarkerSize(1.) ;
227 //       digitgraph-> SetMarkerColor(1) ;
228 //       digitgraph-> Paint("P") ;
229 //     }
230     
231 //     Print() ;
232 //     histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
233 //     histocanvas->Draw() ; 
234 //     histo->Draw("lego1") ; 
235     
236 //     delete[] xi ; 
237 //     delete[] zi ; 
238     
239 //     break;
240 //   }
241   
242 //   case kButton1Up: 
243 //     if (digitgraph) {
244 //       delete digitgraph  ;
245 //       digitgraph = 0 ;
246 //     }
247 //     break;
248   
249 //    }
250 }
251
252 //____________________________________________________________________________
253 void AliPHOSCpvRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
254 {
255   // wraps other methods
256   AliPHOSEmcRecPoint::EvalAll(logWeight,digits) ;
257   EvalClusterLengths(digits) ;
258 }
259 //____________________________________________________________________________
260 void AliPHOSCpvRecPoint::EvalLocalPosition(Float_t logWeight,TClonesArray * digits)
261 {
262   // Calculates the center of gravity in the local PHOS-module coordinates 
263
264   Float_t wtot = 0. ;
265  
266   Int_t relid[4] ;
267
268   Float_t x = 0. ;
269   Float_t z = 0. ;
270   
271   AliPHOSDigit * digit ;
272
273   AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
274   
275   Int_t iDigit;
276
277   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
278     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]);
279
280     Float_t xi ;
281     Float_t zi ;
282     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
283     phosgeom->RelPosInModule(relid, xi, zi);
284     if (fAmp>0 && fEnergyList[iDigit]>0) {
285       Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
286       x    += xi * w ;
287       z    += zi * w ;
288       wtot += w ;
289     }
290     else
291       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
292   }
293
294   if (wtot != 0) {
295     x /= wtot ;
296     z /= wtot ;
297   } else {
298     x = -1e6 ;
299     z = -1e6 ;
300     if (fMulDigit != 0) 
301       AliWarning(Form("Too low log weight factor to evaluate cluster's center" )) ;
302   }
303   fLocPos.SetX(x)  ;
304   fLocPos.SetY(0.) ;
305   fLocPos.SetZ(z)  ;
306   fLocPosM = 0 ;
307
308 }
309
310 //____________________________________________________________________________
311 void AliPHOSCpvRecPoint::EvalClusterLengths(TClonesArray * digits)
312 {
313   //Modified 15.03.2001 by Dmitri Peressounko
314   
315   // Calculates the cluster lengths along X and Z axes
316   // These characteristics are needed for CPV to tune
317   // digitization+reconstruction to experimental data
318   // Yuri Kharlov. 24 October 2000
319
320   Int_t relid[4] ;
321
322   AliPHOSDigit * digit ;
323
324   AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
325
326   const Int_t kMaxLeng=20;
327   Int_t idX[kMaxLeng], idZ[kMaxLeng];
328   fLengX = 0;
329   fLengZ = 0;
330   Bool_t dejavu;
331
332   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
333     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
334     Int_t absId = digit->GetId();
335     phosgeom->AbsToRelNumbering(absId, relid) ;
336
337     Int_t i;
338     dejavu=kFALSE;
339     for (i=0; i<fLengX; i++) if (relid[2]==idX[i]) { dejavu=kTRUE; break; }
340     if (!dejavu) {
341       idX[fLengX]=relid[2];
342       fLengX++;
343       fLengX = TMath::Min(fLengX,kMaxLeng);
344     }
345
346     dejavu=kFALSE;
347     for (i=0; i<fLengZ; i++) if (relid[3]==idZ[i]) { dejavu=kTRUE; break; }
348     if (!dejavu) {
349       idZ[fLengZ]=relid[3];
350       fLengZ++;
351       fLengZ = TMath::Min(fLengZ,kMaxLeng);
352     }
353   }
354 }
355
356
357
358 //____________________________________________________________________________
359 void AliPHOSCpvRecPoint::Print(const Option_t *) const
360 {
361   // Print the list of digits belonging to the cluster
362   
363   TString message ; 
364   message  =  "AliPHOSCpvRecPoint: " ;
365   message +=  "Digits #   " ;
366   AliInfo(Form(message.Data())) ; 
367   
368   Int_t iDigit;
369
370   for(iDigit=0; iDigit<fMulDigit; iDigit++) 
371     printf(" %d \n", fDigitsList[iDigit]) ; 
372
373   printf("Energies: \n")  ;
374   for(iDigit=0; iDigit<fMulDigit; iDigit++) 
375     printf(" %f ", fEnergyList[iDigit]) ; 
376   
377   message  = "       Multiplicity    = %d\n" ;
378   message += "       Cluster Energy  = %f\n" ;
379   message += "       Stored at position %d\n" ; 
380  
381   printf(message.Data(), fMulDigit, fAmp, GetIndexInList() ) ; 
382
383 }