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 }