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, TMatrix & /*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 }