Using the recommended way of forward declarations for TVector and TMatrix (see v5...
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeometry.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 // Geometry class  for PHOS : singleton  
20 // PHOS consists of the electromagnetic calorimeter (EMCA)
21 // and a charged particle veto either in the Subatech's version (PPSD)
22 // or in the IHEP's one (CPV).
23 // The EMCA/PPSD/CPV modules are parametrized so that any configuration
24 // can be easily implemented 
25 // The title is used to identify the version of CPV used.
26 //                  
27 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH)
28
29 // --- ROOT system ---
30
31 #include "TVector3.h"
32 #include "TRotation.h" 
33 #include "TParticle.h"
34
35 // --- Standard library ---
36
37 // --- AliRoot header files ---
38 #include "AliLog.h"
39 #include "AliPHOSGeometry.h"
40 #include "AliPHOSEMCAGeometry.h" 
41 #include "AliPHOSRecPoint.h"
42
43 ClassImp(AliPHOSGeometry)
44
45 // these initialisations are needed for a singleton
46 AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
47 Bool_t            AliPHOSGeometry::fgInit = kFALSE ;
48
49 //____________________________________________________________________________
50 AliPHOSGeometry::AliPHOSGeometry() {
51     // default ctor 
52     // must be kept public for root persistency purposes, but should never be called by the outside world
53     fPHOSAngle      = 0 ;
54     fGeometryEMCA   = 0 ;
55     fGeometrySUPP   = 0 ;
56     fGeometryCPV    = 0 ;
57     fgGeom          = 0 ;
58     fRotMatrixArray = 0 ;  
59 }  
60
61 //____________________________________________________________________________
62 AliPHOSGeometry::~AliPHOSGeometry(void)
63 {
64   // dtor
65
66   if (fRotMatrixArray) fRotMatrixArray->Delete() ; 
67   if (fRotMatrixArray) delete fRotMatrixArray ; 
68   if (fPHOSAngle     ) delete[] fPHOSAngle ; 
69 }
70 //____________________________________________________________________________
71
72 void AliPHOSGeometry::Init(void)
73 {
74   // Initializes the PHOS parameters :
75   //  IHEP is the Protvino CPV (cathode pad chambers)
76   
77   TString test(GetName()) ; 
78   if (test != "IHEP" ) {
79     AliFatal(Form("%s is not a known geometry (choose among IHEP)", 
80                   test.Data() )) ; 
81   }
82
83   fgInit     = kTRUE ; 
84   
85   fNModules     = 5;
86   fAngle        = 20;
87   
88   fGeometryEMCA = new AliPHOSEMCAGeometry();
89   
90   fGeometryCPV  = new AliPHOSCPVGeometry ();
91   
92   fGeometrySUPP = new AliPHOSSupportGeometry();
93   
94   fPHOSAngle = new Float_t[fNModules] ;
95   
96   Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
97   
98   fPHOSParams[0] =  TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2., 
99                                (Double_t)(emcParams[0] - (emcParams[1]-emcParams[0])*
100                                           fGeometryCPV->GetCPVBoxSize(1)/2/emcParams[3]));
101   fPHOSParams[1] = emcParams[1] ;
102   fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
103   fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
104   
105   fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
106   
107   Int_t index ;
108   for ( index = 0; index < fNModules; index++ )
109     fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
110   
111   this->SetPHOSAngles() ; 
112   fRotMatrixArray = new TObjArray(fNModules) ; 
113   
114 }
115
116 //____________________________________________________________________________
117 AliPHOSGeometry *  AliPHOSGeometry::GetInstance() 
118
119   // Returns the pointer of the unique instance; singleton specific
120   
121   return static_cast<AliPHOSGeometry *>( fgGeom ) ; 
122 }
123
124 //____________________________________________________________________________
125 AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title) 
126 {
127   // Returns the pointer of the unique instance
128   // Creates it with the specified options (name, title) if it does not exist yet
129
130   AliPHOSGeometry * rv = 0  ; 
131   if ( fgGeom == 0 ) {
132     if ( strcmp(name,"") == 0 ) 
133       rv = 0 ;
134     else {    
135       fgGeom = new AliPHOSGeometry(name, title) ;
136       if ( fgInit )
137         rv = (AliPHOSGeometry * ) fgGeom ;
138       else {
139         rv = 0 ; 
140         delete fgGeom ; 
141         fgGeom = 0 ; 
142       }
143     }
144   }
145   else {
146     if ( strcmp(fgGeom->GetName(), name) != 0 ) 
147       ::Error("GetInstance", "Current geometry is %s. You cannot call %s", 
148                       fgGeom->GetName(), name) ; 
149     else
150       rv = (AliPHOSGeometry *) fgGeom ; 
151   } 
152   return rv ; 
153 }
154
155 //____________________________________________________________________________
156 void AliPHOSGeometry::SetPHOSAngles() 
157
158   // Calculates the position of the PHOS modules in ALICE global coordinate system
159   
160   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
161   Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
162   pphi *= kRADDEG ;
163   if (pphi > fAngle){ 
164     AliError(Form("PHOS modules overlap!\n pphi = %f fAngle = %f", 
165                   pphi, fAngle));
166
167   }
168   pphi = fAngle;
169   
170   for( Int_t i = 1; i <= fNModules ; i++ ) {
171     Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
172     fPHOSAngle[i-1] = -  angle ;
173   } 
174 }
175
176 //____________________________________________________________________________
177 Bool_t AliPHOSGeometry::AbsToRelNumbering(Int_t AbsId, Int_t * relid) const
178 {
179   // Converts the absolute numbering into the following array/
180   //  relid[0] = PHOS Module number 1:fNModules 
181   //  relid[1] = 0 if PbW04
182   //           = -1 if CPV
183   //  relid[2] = Row number inside a PHOS module
184   //  relid[3] = Column number inside a PHOS module
185
186   Bool_t rv  = kTRUE ; 
187   Float_t id = AbsId ;
188
189   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ; 
190   
191   if ( phosmodulenumber >  GetNModules() ) { // it is a CPV pad
192     
193     id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
194     Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
195     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
196     relid[1] = -1 ;
197     id -= ( relid[0] - 1 ) * nCPV ; 
198     relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
199     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; 
200   } 
201   else { // it is a PW04 crystal
202
203     relid[0] = phosmodulenumber ;
204     relid[1] = 0 ;
205     id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ; 
206     relid[2] = (Int_t)TMath::Ceil( id / GetNZ() )  ;
207     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ; 
208   } 
209   return rv ; 
210 }
211
212 //____________________________________________________________________________  
213 void AliPHOSGeometry::EmcModuleCoverage(Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) const 
214 {
215   // calculates the angular coverage in theta and phi of one EMC (=PHOS) module
216
217  Double_t conv ; 
218   if ( opt == Radian() ) 
219     conv = 1. ; 
220   else if ( opt == Degre() )
221     conv = 180. / TMath::Pi() ; 
222   else {
223     AliWarning(Form("%s unknown option; result in radian", opt)) ; 
224     conv = 1. ;
225       }
226
227   Float_t phi = GetPHOSAngle(mod) *  (TMath::Pi() / 180.)  ;  
228   Float_t y0  = GetIPtoCrystalSurface() ; 
229   Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
230   Float_t x0  = fGeometryEMCA->GetNStripX()*strip[0] ;
231   Float_t z0  = fGeometryEMCA->GetNStripZ()*strip[2] ;
232   Double_t angle = TMath::ATan( x0 / y0 ) ;
233   phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 220 and 320 deg.)
234   Double_t max  = phi - angle ;
235   Double_t min   = phi + angle ;
236   pM = TMath::Max(max, min) * conv ;
237   pm = TMath::Min(max, min) * conv ; 
238   
239   angle =  TMath::ATan( z0 /  y0 ) ;
240   max  = TMath::Pi() / 2.  + angle ; // to follow the convention of the particle generator(PHOS is at 90 deg.)
241   min  = TMath::Pi() / 2.  - angle ;
242   tM = TMath::Max(max, min) * conv ;
243   tm = TMath::Min(max, min) * conv ; 
244  
245 }
246
247 //____________________________________________________________________________  
248 void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t * opt) const
249 {
250   // calculates the angular coverage in theta and phi of a single crystal in a EMC(=PHOS) module
251
252   Double_t conv ; 
253   if ( opt == Radian() ) 
254     conv = 1. ; 
255   else if ( opt == Degre() )
256     conv = 180. / TMath::Pi() ; 
257   else {
258     AliWarning(Form("%s unknown option; result in radian", opt)) ;  
259     conv = 1. ;
260       }
261
262   Float_t y0  = GetIPtoCrystalSurface() ; 
263   theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
264   phi   = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
265 }
266  
267
268 //____________________________________________________________________________
269 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrixF & /*gmat*/) const
270 {
271   // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
272  
273   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
274   TVector3 localposition ;
275
276   tmpPHOS->GetLocalPosition(gpos) ;
277
278
279   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
280     {  gpos.SetY( - GetIPtoCrystalSurface()) ;  
281
282     }
283   else
284     { // it is a CPV
285       gpos.SetY(- GetIPtoUpperCPVsurface()  ) ; 
286     }  
287
288   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
289   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
290   Float_t rphi          = phi / kRADDEG ; 
291   
292   TRotation rot ;
293   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
294   
295   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
296   gpos.Transform(rot) ; // rotate the baby 
297
298 }
299
300 //____________________________________________________________________________
301 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const 
302 {
303   // Calculates the coordinates of a RecPoint in the ALICE global coordinate system 
304
305   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
306   TVector3 localposition ;
307   tmpPHOS->GetLocalPosition(gpos) ;
308
309
310   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
311     {  gpos.SetY( - GetIPtoCrystalSurface()  ) ;  
312     }
313   else
314     { // it is a CPV
315           gpos.SetY(- GetIPtoUpperCPVsurface()  ) ; 
316     }  
317
318   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
319   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
320   Float_t rphi          = phi / kRADDEG ; 
321   
322   TRotation rot ;
323   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
324   
325   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
326   gpos.Transform(rot) ; // rotate the baby 
327 }
328
329 //____________________________________________________________________________
330 void AliPHOSGeometry::ImpactOnEmc(Double_t theta, Double_t phi, Int_t & moduleNumber, Double_t & z, Double_t & x) const
331 {
332   // calculates the impact coordinates on PHOS of a neutral particle  
333   // emitted in the direction theta and phi in the ALICE global coordinate system
334
335   // searches for the PHOS EMC module
336
337   moduleNumber = 0 ; 
338   Double_t tm, tM, pm, pM ; 
339   Int_t index = 1 ; 
340   while ( moduleNumber == 0 && index <= GetNModules() ) { 
341     EmcModuleCoverage(index, tm, tM, pm, pM) ; 
342     if ( (theta >= tm && theta <= tM) && (phi >= pm && phi <= pM ) ) 
343       moduleNumber = index ; 
344     index++ ;    
345   }
346   if ( moduleNumber != 0 ) {
347     Float_t phi0 =  GetPHOSAngle(moduleNumber) *  (TMath::Pi() / 180.) + 1.5 * TMath::Pi()  ;  
348     Float_t y0  =  GetIPtoCrystalSurface()  ;   
349     Double_t angle = phi - phi0; 
350     x = y0 * TMath::Tan(angle) ; 
351     angle = theta - TMath::Pi() / 2 ; 
352     z = y0 * TMath::Tan(angle) ; 
353   }
354 }
355
356 //____________________________________________________________________________
357 void AliPHOSGeometry::ImpactOnEmc(const TVector3& vec, Int_t & moduleNumber, Double_t & z, Double_t & x) const
358 {
359   // calculates the impact coordinates on PHOS of a neutral particle  
360   // emitted in the direction theta and phi in the ALICE global coordinate system
361   // searches for the PHOS EMC module
362
363   Double_t theta = vec.Theta() ; 
364   Double_t phi   = vec.Phi() ; 
365
366   ImpactOnEmc(theta, phi, moduleNumber, z, x) ;
367 }
368
369 //____________________________________________________________________________
370 void AliPHOSGeometry::ImpactOnEmc(const TParticle& p, Int_t & moduleNumber, Double_t & z, Double_t & x) const
371 {
372   // calculates the impact coordinates on PHOS of a neutral particle  
373   // emitted in the direction theta and phi in the ALICE global coordinate system
374
375   // searches for the PHOS EMC module
376   Double_t theta = p.Theta() ; 
377   Double_t phi   = p.Phi() ; 
378
379   ImpactOnEmc(theta, phi, moduleNumber, z, x) ;
380 }
381
382 //____________________________________________________________________________
383 Bool_t  AliPHOSGeometry::Impact(const TParticle * particle) const 
384 {
385   // Tells if a particle enters PHOS
386   Bool_t in=kFALSE;
387   Int_t moduleNumber=0;
388   Double_t z,x;
389   ImpactOnEmc(particle->Theta(),particle->Phi(),moduleNumber,z,x);
390   if(moduleNumber) 
391     in=kTRUE;
392   else 
393     in=kFALSE;
394   return in;
395 }
396
397 //____________________________________________________________________________
398 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId) const
399 {
400   // Converts the relative numbering into the absolute numbering
401   // EMCA crystals:
402   //  AbsId = from 1 to fNModules * fNPhi * fNZ
403   // CPV pad:
404   //  AbsId = from N(total PHOS crystals) + 1
405   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
406
407   Bool_t rv = kTRUE ; 
408   
409   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
410     AbsId =
411       ( relid[0] - 1 ) * GetNPhi() * GetNZ()         // the offset of PHOS modules
412       + ( relid[2] - 1 ) * GetNZ()                   // the offset along phi
413       +   relid[3] ;                                 // the offset along z
414   }
415   else { // it is a CPV pad
416     AbsId =    GetNPhi() * GetNZ() *  GetNModules()         // the offset to separate EMCA crystals from CPV pads
417       + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()   // the pads offset of PHOS modules 
418       + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                             // the pads offset of a CPV row
419       +   relid[3] ;                                                         // the column number
420   }
421   
422   return rv ; 
423 }
424
425 //____________________________________________________________________________
426
427 void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const
428 {
429   // Converts the absolute numbering into the global ALICE coordinate system
430   
431     
432     Int_t relid[4] ;
433     
434     AbsToRelNumbering(id , relid) ;
435     
436     Int_t phosmodule = relid[0] ; 
437     
438     Float_t y0 = 0 ; 
439     
440     if ( relid[1] == 0 )  // it is a PbW04 crystal 
441       y0 =  - GetIPtoCrystalSurface() ;  
442     else
443       y0 =  - GetIPtoUpperCPVsurface() ; 
444
445     Float_t x, z ; 
446     RelPosInModule(relid, x, z) ; 
447     
448     pos.SetX(x) ;
449     pos.SetZ(z) ;
450     pos.SetY(y0) ;
451     
452     Float_t phi           = GetPHOSAngle( phosmodule) ; 
453     Double_t const kRADDEG = 180.0 / TMath::Pi() ;
454     Float_t rphi          = phi / kRADDEG ; 
455     
456     TRotation rot ;
457     rot.RotateZ(-rphi) ; // a rotation around Z by angle  
458     
459     TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
460     
461     pos.Transform(rot) ; // rotate the baby 
462
463
464 //____________________________________________________________________________
465 void AliPHOSGeometry::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & AbsId) const
466 {
467   // converts local PHOS-module (x, z) coordinates to absId 
468   Int_t relid[4] ;
469   relid[0] = module ;
470   relid[1] = 0 ;
471   relid[2] = static_cast<Int_t>(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) );
472   relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ GetCellStep() + GetNZ()   / 2.) ) ;
473   
474   RelToAbsNumbering(relid,AbsId) ;
475 }
476
477 //____________________________________________________________________________
478 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const 
479 {
480   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
481   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
482   
483   Int_t row        = relid[2] ; //offset along x axis
484   Int_t column     = relid[3] ; //offset along z axis
485
486   
487   if ( relid[1] == 0 ) { // its a PbW04 crystal
488     x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCellStep() ; // position of Xtal with respect
489     z = - ( GetNZ()  /2. - column + 0.5 ) *  GetCellStep() ; // of center of PHOS module  
490   }  
491   else  {    
492     x = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
493     z = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module  
494   }
495 }
496
497 //____________________________________________________________________________
498
499 void AliPHOSGeometry::GetModuleCenter(TVector3& center, 
500                                       const char *det,
501                                       Int_t module) const
502 {
503   // Returns a position of the center of the CPV or EMC module
504   Float_t rDet = 0.;
505   if      (strcmp(det,"CPV") == 0) rDet  = GetIPtoCPVDistance   ();
506   else if (strcmp(det,"EMC") == 0) rDet  = GetIPtoCrystalSurface();
507   else 
508     AliFatal(Form("Wrong detector name %s",det));
509
510   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
511   angle *= TMath::Pi()/180;
512   angle += 3*TMath::Pi()/2.;
513   center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
514 }
515
516 //____________________________________________________________________________
517
518 void AliPHOSGeometry::Global2Local(TVector3& localPosition,
519                                    const TVector3& globalPosition,
520                                    Int_t module) const
521 {
522   // Transforms a global position of the rec.point to the local coordinate system
523   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
524   angle *= TMath::Pi()/180;
525   angle += 3*TMath::Pi()/2.;
526   localPosition = globalPosition;
527   localPosition.RotateZ(-angle);
528 }