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