PHOS support added
[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)
28
29 // --- ROOT system ---
30
31 #include "TVector3.h"
32 #include "TRotation.h" 
33
34 // --- Standard library ---
35
36 #include <iostream.h>
37
38 // --- AliRoot header files ---
39
40 #include "AliPHOSGeometry.h"
41 #include "AliPHOSPpsdRecPoint.h"
42 #include "AliConst.h"
43
44 ClassImp(AliPHOSGeometry) ;
45
46 AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
47 Bool_t            AliPHOSGeometry::fgInit = kFALSE ;
48
49 //____________________________________________________________________________
50 AliPHOSGeometry::~AliPHOSGeometry(void)
51 {
52   // dtor
53
54   if (fRotMatrixArray) fRotMatrixArray->Delete() ; 
55   if (fRotMatrixArray) delete fRotMatrixArray ; 
56   if (fPHOSAngle     ) delete fPHOSAngle ; 
57 //    if (fGeometryEMCA  ) detete fGeometryEMCA;
58 //    if (fGeometryCPV   ) detete fGeometryCPV ;
59 }
60
61 //____________________________________________________________________________
62
63 void AliPHOSGeometry::Init(void)
64 {
65   // Initializes the PHOS parameters
66
67   if ( ((strcmp( fName, "default" )) == 0) || 
68        ((strcmp( fName, "GPS2" ))    == 0) ||
69        ((strcmp( fName, "IHEP" ))    == 0) ) {
70     fgInit     = kTRUE ; 
71                                              fGeometryEMCA = new AliPHOSEMCAGeometry();
72     if ( ((strcmp( fName, "GPS2" ))  == 0) ) fGeometryCPV  = new AliPHOSPPSDGeometry();
73     if ( ((strcmp( fName, "IHEP" ))  == 0) ) fGeometryCPV  = new AliPHOSCPVGeometry ();
74                                              fGeometrySUPP = new AliPHOSSupportGeometry();
75     fNModules = 5;
76     fPHOSAngle = new Float_t[fNModules] ;
77     Int_t index ;
78     for ( index = 0; index < fNModules; index++ )
79     fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
80
81     this->SetPHOSAngles() ; 
82     fRotMatrixArray = new TObjArray(fNModules) ; 
83   }
84  else {
85    fgInit = kFALSE ; 
86    cout << "PHOS Geometry setup: option not defined " << fName << endl ; 
87  }
88 }
89
90 //____________________________________________________________________________
91 AliPHOSGeometry *  AliPHOSGeometry::GetInstance() 
92
93   // Returns the pointer of the unique instance
94   return (AliPHOSGeometry *) fgGeom ; 
95 }
96
97 //____________________________________________________________________________
98 AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title) 
99 {
100   // Returns the pointer of the unique instance
101   AliPHOSGeometry * rv = 0  ; 
102   if ( fgGeom == 0 ) {
103     if ( strcmp(name,"") == 0 ) 
104       rv = 0 ;
105     else {    
106       fgGeom = new AliPHOSGeometry(name, title) ;
107       if ( fgInit )
108         rv = (AliPHOSGeometry * ) fgGeom ;
109       else {
110         rv = 0 ; 
111         delete fgGeom ; 
112         fgGeom = 0 ; 
113       }
114     }
115   }
116   else {
117     if ( strcmp(fgGeom->GetName(), name) != 0 ) {
118       cout << "AliPHOSGeometry <E> : current geometry is " << fgGeom->GetName() << endl
119            << "                      you cannot call     " << name << endl ; 
120     }
121     else
122       rv = (AliPHOSGeometry *) fgGeom ; 
123   } 
124   return rv ; 
125 }
126
127 //____________________________________________________________________________
128 void AliPHOSGeometry::SetPHOSAngles() 
129
130   // Calculates the position in ALICE of the PHOS modules
131   
132   Double_t const kRADDEG = 180.0 / kPI ;
133   Float_t pphi =  TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
134   pphi *= kRADDEG ;
135   
136   for( Int_t i = 1; i <= fNModules ; i++ ) {
137     Float_t angle = pphi * 2 * ( i - fNModules / 2.0 - 0.5 ) ;
138     fPHOSAngle[i-1] = -  angle ;
139   } 
140 }
141
142 //____________________________________________________________________________
143 Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid)
144 {
145   // Converts the absolute numbering into the following array/
146   //  relid[0] = PHOS Module number 1:fNModules 
147   //  relid[1] = 0 if PbW04
148   //           = PPSD Module number 1:fNumberOfModulesPhi*fNumberOfModulesZ*2 (2->up and bottom level)
149   //  relid[2] = Row number inside a PHOS or PPSD module
150   //  relid[3] = Column number inside a PHOS or PPSD module
151
152   Bool_t rv  = kTRUE ; 
153   Float_t id = AbsId ;
154
155   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / ( GetNPhi() * GetNZ() ) ) ; 
156   
157   if ( phosmodulenumber >  GetNModules() ) { // it is a PPSD or CPV pad
158
159     if      ( strcmp(fName,"GPS2") == 0 ) {
160       id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
161       Float_t tempo = 2 *  GetNumberOfModulesPhi() * GetNumberOfModulesZ() *  GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; 
162       relid[0] = (Int_t)TMath::Ceil( id / tempo ) ; 
163       id -= ( relid[0] - 1 ) * tempo ;
164       relid[1] = (Int_t)TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ; 
165       id -= ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
166       relid[2] = (Int_t)TMath::Ceil( id / GetNumberOfPadsPhi() ) ;
167       relid[3] = (Int_t) ( id - ( relid[2] - 1 )  * GetNumberOfPadsPhi() ) ; 
168     }
169     else if ( strcmp(fName,"IHEP") == 0 ) {
170       id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
171       relid[0] = (Int_t) TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ;
172       relid[1] = 1 ;
173       id -= ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()  ; 
174       relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfPadsZ() ) ;
175       relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfPadsZ() ) ; 
176     }
177   } 
178   else { // its a PW04 crystal
179
180     relid[0] = phosmodulenumber ;
181     relid[1] = 0 ;
182     id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ; 
183     relid[2] = (Int_t)TMath::Ceil( id / GetNPhi() ) ;
184     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNPhi() ) ; 
185   } 
186   return rv ; 
187 }
188
189 //____________________________________________________________________________  
190 void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) 
191 {
192   // calculates the angular coverage in theta and phi of a EMC module
193
194  Double_t conv ; 
195   if ( opt == Radian() ) 
196     conv = 1. ; 
197   else if ( opt == Degre() )
198     conv = 180. / TMath::Pi() ; 
199   else {
200     cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
201     conv = 1. ;
202       }
203
204   Float_t phi =  GetPHOSAngle(mod) *  (TMath::Pi() / 180.)  ;  
205   Float_t y0  =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
206                   + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;  
207   
208   Double_t angle = TMath::ATan( GetCrystalSize(0)*GetNPhi() / (2 * y0) ) ;
209   phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 230 and 310 deg.)
210   Double_t max  = phi - angle ;
211   Double_t min   = phi + angle ;
212   pM = TMath::Max(max, min) * conv ;
213   pm = TMath::Min(max, min) * conv ; 
214   
215   angle =  TMath::ATan( GetCrystalSize(2)*GetNZ() / (2 * y0) ) ;
216   max  = TMath::Pi() / 2.  + angle ; // to follow the convention of the particle generator(PHOS is at 90 deg.)
217   min  = TMath::Pi() / 2.  - angle ;
218   tM = TMath::Max(max, min) * conv ;
219   tm = TMath::Min(max, min) * conv ; 
220  
221 }
222
223 //____________________________________________________________________________  
224 void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t * opt) 
225 {
226   // calculates the angular coverage in theta and phi of a single crystal in a EMC module
227
228   Double_t conv ; 
229   if ( opt == Radian() ) 
230     conv = 1. ; 
231   else if ( opt == Degre() )
232     conv = 180. / TMath::Pi() ; 
233   else {
234     cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
235     conv = 1. ;
236       }
237
238   Float_t y0   =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
239     + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;  
240   theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
241   phi   = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
242 }
243  
244
245 //____________________________________________________________________________
246 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
247 {
248   // Calculates the ALICE global coordinates of a RecPoint and the error matrix
249  
250   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
251   TVector3 localposition ;
252
253   tmpPHOS->GetLocalPosition(gpos) ;
254
255
256   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
257     {  gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() +
258                     GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ;  
259
260     }
261   else
262     { // it is a PPSD pad
263       AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ;
264       if (tmpPpsd->GetUp() ) // it is an upper module
265         {
266           gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - 
267                        GetLeadToMicro2Gap() - GetLeadConverterThickness() -  
268                        GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 )  ) ; 
269         } 
270       else // it is a lower module
271         gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; 
272     }  
273
274   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
275   Double_t const kRADDEG = 180.0 / kPI ;
276   Float_t rphi          = phi / kRADDEG ; 
277   
278   TRotation rot ;
279   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
280   
281   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
282   gpos.Transform(rot) ; // rotate the baby 
283
284 }
285
286 //____________________________________________________________________________
287 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const 
288 {
289   // Calculates the ALICE global coordinates of a RecPoint 
290
291   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
292   TVector3 localposition ;
293   tmpPHOS->GetLocalPosition(gpos) ;
294
295
296   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
297     {  gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() +
298                     GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ;  
299     }
300   else
301     { // it is a PPSD pad
302       AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ;
303       if (tmpPpsd->GetUp() ) // it is an upper module
304         {
305           gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - 
306                        GetLeadToMicro2Gap() - GetLeadConverterThickness() -  
307                        GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 )  ) ; 
308         } 
309       else // it is a lower module
310         gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; 
311     }  
312
313   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
314   Double_t const kRADDEG = 180.0 / kPI ;
315   Float_t rphi          = phi / kRADDEG ; 
316   
317   TRotation rot ;
318   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
319   
320   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
321   gpos.Transform(rot) ; // rotate the baby 
322 }
323
324 //____________________________________________________________________________
325 void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x) 
326 {
327   // calculates the impact coordinates of a neutral particle  
328   // emitted in direction theta and phi in ALICE
329
330   // searches for the PHOS EMC module
331   ModuleNumber = 0 ; 
332   Double_t tm, tM, pm, pM ; 
333   Int_t index = 1 ; 
334   while ( ModuleNumber == 0 && index <= GetNModules() ) { 
335     EmcModuleCoverage(index, tm, tM, pm, pM) ; 
336     if ( (theta >= tm && theta <= tM) && (phi >= pm && phi <= pM ) ) 
337       ModuleNumber = index ; 
338     index++ ;    
339   }
340   if ( ModuleNumber != 0 ) {
341     Float_t phi0 =  GetPHOSAngle(ModuleNumber) *  (TMath::Pi() / 180.) + 1.5 * TMath::Pi()  ;  
342     Float_t y0  =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
343       + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;   
344     Double_t angle = phi - phi0; 
345     x = y0 * TMath::Tan(angle) ; 
346     angle = theta - TMath::Pi() / 2 ; 
347     z = y0 * TMath::Tan(angle) ; 
348   }
349 }
350
351 //____________________________________________________________________________
352 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId)
353 {
354   // Converts the relative numbering into the absolute numbering
355   //  AbsId = 1 to fNModules * fNPhi * fNZ                                                                       -> PbWO4
356   //  AbsId = N(total PHOS crystals) +
357   //          1 to fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ -> PPSD
358   //  AbsId = N(total PHOS crystals) +
359   //          1:fNModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ                                              -> CPV
360
361   Bool_t rv = kTRUE ; 
362  
363   if      ( relid[1]  >  0 ) { // it is a PPSD pad
364     AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from PPSD pads
365       + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ()   // the pads offset of PHOS modules 
366                          * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
367       + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PPSD modules 
368       + ( relid[2] - 1 ) * GetNumberOfPadsPhi()                              // the pads offset of a PPSD row
369       +   relid[3] ;                                                         // the column number
370   } 
371
372   else if ( relid[1] ==  0 ) { // it is a Phos crystal
373     AbsId =
374         ( relid[0] - 1 ) * GetNPhi() * GetNZ()                               // the offset of PHOS modules
375       + ( relid[2] - 1 ) * GetNPhi()                                         // the offset of a xtal row
376       +   relid[3] ;                                                         // the column number
377   }
378
379   else if ( relid[1] == -1 ) { // it is a CPV pad
380     AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from CPV pads
381       + ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PHOS modules 
382       + ( relid[2] - 1 ) * GetNumberOfPadsZ()                                // the pads offset of a CPV row
383       +   relid[3] ;                                                         // the column number
384   }
385   
386   return rv ; 
387 }
388
389 //____________________________________________________________________________
390
391 void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) 
392 {
393   // Converts the absolute numbering into the global ALICE coordinates
394   
395    if (id > 0) { 
396
397   Int_t relid[4] ;
398  
399   AbsToRelNumbering(id , relid) ;
400
401   Int_t phosmodule = relid[0] ; 
402
403   Float_t y0 = 0 ; 
404
405   if ( relid[1] == 0 ) // it is a PbW04 crystal 
406   {  y0 =  -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
407       + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness())  ;  
408   }
409   if ( relid[1] > 0 ) { // its a PPSD pad
410     if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() ) // its an bottom module
411      {
412        y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0)  ;
413      } 
414     else // its an upper module
415       y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
416         -  GetLeadConverterThickness() -  GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ; 
417   }
418
419   Float_t x, z ; 
420   RelPosInModule(relid, x, z) ; 
421
422   pos.SetX(x) ;
423   pos.SetZ(z) ;
424   pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ; 
425
426
427
428    Float_t phi           = GetPHOSAngle( phosmodule) ; 
429    Double_t const kRADDEG = 180.0 / kPI ;
430    Float_t rphi          = phi / kRADDEG ; 
431
432    TRotation rot ;
433    rot.RotateZ(-rphi) ; // a rotation around Z by angle  
434   
435    TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
436   
437    pos.Transform(rot) ; // rotate the baby 
438   }
439   else {
440  pos.SetX(0.);
441  pos.SetY(0.);
442  pos.SetZ(0.);
443        }
444
445
446 //____________________________________________________________________________
447 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) 
448 {
449   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
450   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
451   
452   Int_t ppsdmodule  ; 
453   Float_t x0,z0;
454   Int_t row        = relid[2] ; //offset along x axiz
455   Int_t column     = relid[3] ; //offset along z axiz
456
457   Float_t padsizeZ = GetPadSizeZ();
458   Float_t padsizeX = GetPadSizePhi();
459
460   if ( relid[1] == 0 ) { // its a PbW04 crystal 
461     x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCrystalSize(0) ; // position ox Xtal with respect
462     z =   ( GetNZ()  /2. - column + 0.5 ) *  GetCrystalSize(2) ; // of center of PHOS module  
463   }  
464   else  {    
465     if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() )
466       ppsdmodule =  relid[1]-GetNumberOfModulesPhi() *  GetNumberOfModulesZ(); 
467     else
468       ppsdmodule =  relid[1] ;
469     Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ; 
470     Int_t modcol = ppsdmodule -  ( modrow - 1 ) * GetNumberOfModulesPhi() ;     
471     if ( ((strcmp( fName, "GPS2" ))  == 0) ) {
472       x0 = (  GetNumberOfModulesPhi() / 2.  - modrow  + 0.5 ) * GetPPSDModuleSize(0) ;
473       z0 = (  GetNumberOfModulesZ()   / 2.  - modcol  + 0.5 ) * GetPPSDModuleSize(2)  ;     
474     } else {
475       x0 = 0;
476       z0 = 0;
477     }
478     x = - ( GetNumberOfPadsPhi()/2. - row    - 0.5 ) * padsizeX + x0 ; // position of pad  with respect
479     z =   ( GetNumberOfPadsZ()  /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module  
480   }
481 }