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 "TFolder.h"
  34 #include "TROOT.h"
  35
  36 // --- Standard library ---
  37
  38 #include <iostream.h>
  39 #include <stdlib.h>
  40
  41 // --- AliRoot header files ---
  42
  43 #include "AliPHOSGeometry.h"
  44 #include "AliPHOSEMCAGeometry.h"
  45 #include "AliPHOSRecPoint.h"
  46 #include "AliConst.h"
  47
  48 ClassImp(AliPHOSGeometry) ;
  49
  50 // these initialisations are needed for a singleton
  51 AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
  52 Bool_t            AliPHOSGeometry::fgInit = kFALSE ;
  53
  54 //____________________________________________________________________________
  55 AliPHOSGeometry::~AliPHOSGeometry(void)
  56 {
  57   // dtor
  58
  59   if (fRotMatrixArray) fRotMatrixArray->Delete() ;
  60   if (fRotMatrixArray) delete fRotMatrixArray ;
  61   if (fPHOSAngle     ) delete fPHOSAngle ;
  62 }
  63 //____________________________________________________________________________
  64
  65 void AliPHOSGeometry::Init(void)
  66 {
  67   // Initializes the PHOS parameters :
  68   //  IHEP is the Protvino CPV (cathode pad chambers)
  69   //  GPS2 is the Subatech Pre-Shower (two micromegas sandwiching a passive lead converter)
  70   //  MIXT 4 PHOS modules withe the IHEP CPV qnd one PHOS module with the Subatche Pre-Shower
  71
  72   fgInit     = kTRUE ;
  73
  74   fNModules     = 5;
  75   fAngle        = 20;
  76
  77   fGeometryEMCA = new AliPHOSEMCAGeometry();
  78
  79   fGeometryCPV  = new AliPHOSCPVGeometry ();
  80
  81   fGeometrySUPP = new AliPHOSSupportGeometry();
  82
  83   fPHOSAngle = new Float_t[fNModules] ;
  84
  85   Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
  86
  87   fPHOSParams[0] =  TMath::Max(fGeometryCPV->GetCPVBoxSize(0)/2.,
  88                                (emcParams[0]*(fGeometryCPV->GetCPVBoxSize(1)+emcParams[3]) -
  89                                 emcParams[1]* fGeometryCPV->GetCPVBoxSize(1))/emcParams[3] ) ;
  90   fPHOSParams[1] = emcParams[1] ;
  91   fPHOSParams[2] = TMath::Max(emcParams[2], fGeometryCPV->GetCPVBoxSize(2)/2.);
  92   fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
  93
  94   fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
  95
  96   Int_t index ;
  97   for ( index = 0; index < fNModules; index++ )
  98     fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
  99
 100   this->SetPHOSAngles() ;
 101   fRotMatrixArray = new TObjArray(fNModules) ;
 102
 103   // post the geometry into the appropriate folder
 104   TFolder * folder = (TFolder*)gROOT->FindObjectAny("YSAlice/WhiteBoard/Geometry/PHOS");
 105   if ( !folder ) {
 106     cerr << "ERROR: AliPHOSGeometry::Init -> No WhiteBoard/Geometry/PHOS found !" << endl ;
 107     abort();
 108   } else {
 109     folder->SetOwner() ;
 110     folder->Add(this) ;
 111   }
 112
 113 }
 114
 115
 116 //____________________________________________________________________________
 117 AliPHOSGeometry *  AliPHOSGeometry::GetInstance()
 118 {
 119   // Returns the pointer of the unique instance; singleton specific
 120
 121   return (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       cout << "AliPHOSGeometry <E> : current geometry is " << fgGeom->GetName() << endl
 148            << "                      you cannot call     " << name << endl ;
 149     }
 150     else
 151       rv = (AliPHOSGeometry *) fgGeom ;
 152   }
 153   return rv ;
 154 }
 155
 156 //____________________________________________________________________________
 157 void AliPHOSGeometry::SetPHOSAngles()
 158 {
 159   // Calculates the position of the PHOS modules in ALICE global coordinate system
 160
 161   Double_t const kRADDEG = 180.0 / kPI ;
 162   Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
 163   pphi *= kRADDEG ;
 164   if (pphi > fAngle){
 165     cout << "AliPHOSGeometry: PHOS modules overlap!\n";
 166     cout <<  "pphi = " << pphi << " fAngle  " << fAngle << endl ;
 167
 168   }
 169   pphi = fAngle;
 170
 171   for( Int_t i = 1; i <= fNModules ; i++ ) {
 172     Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
 173     fPHOSAngle[i-1] = -  angle ;
 174   }
 175 }
 176
 177 //____________________________________________________________________________
 178 Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid) const
 179 {
 180   // Converts the absolute numbering into the following array/
 181   //  relid[0] = PHOS Module number 1:fNModules
 182   //  relid[1] = 0 if PbW04
 183   //           = -1 if CPV
 184   //  relid[2] = Row number inside a PHOS module
 185   //  relid[3] = Column number inside a PHOS module
 186
 187   Bool_t rv  = kTRUE ;
 188   Float_t id = AbsId ;
 189
 190   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ;
 191
 192   if ( phosmodulenumber >  GetNModules() ) { // it is a CPV pad
 193
 194     id -=  GetNPhi() * GetNZ() *  GetNModules() ;
 195     Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
 196     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
 197     relid[1] = -1 ;
 198     id -= ( relid[0] - 1 ) * nCPV ;
 199     relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
 200     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ;
 201   }
 202   else { // it is a PW04 crystal
 203
 204     relid[0] = phosmodulenumber ;
 205     relid[1] = 0 ;
 206     id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ;
 207     relid[2] = (Int_t)TMath::Ceil( id / GetNZ() )  ;
 208     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ;
 209   }
 210   return rv ;
 211 }
 212
 213 //____________________________________________________________________________
 214 void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) const
 215 {
 216   // calculates the angular coverage in theta and phi of one EMC (=PHOS) module
 217
 218  Double_t conv ;
 219   if ( opt == Radian() )
 220     conv = 1. ;
 221   else if ( opt == Degre() )
 222     conv = 180. / TMath::Pi() ;
 223   else {
 224     cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ;
 225     conv = 1. ;
 226       }
 227
 228   Float_t phi = GetPHOSAngle(mod) *  (TMath::Pi() / 180.)  ;
 229   Float_t y0  = GetIPtoCrystalSurface() ;
 230   Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
 231   Float_t x0  = fGeometryEMCA->GetNStripX()*strip[0] ;
 232   Float_t z0  = fGeometryEMCA->GetNStripZ()*strip[2] ;
 233   Double_t angle = TMath::ATan( x0 / y0 ) ;
 234   phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 220 and 320 deg.)
 235   Double_t max  = phi - angle ;
 236   Double_t min   = phi + angle ;
 237   pM = TMath::Max(max, min) * conv ;
 238   pm = TMath::Min(max, min) * conv ;
 239
 240   angle =  TMath::ATan( z0 /  y0 ) ;
 241   max  = TMath::Pi() / 2.  + angle ; // to follow the convention of the particle generator(PHOS is at 90 deg.)
 242   min  = TMath::Pi() / 2.  - angle ;
 243   tM = TMath::Max(max, min) * conv ;
 244   tm = TMath::Min(max, min) * conv ;
 245
 246 }
 247
 248 //____________________________________________________________________________
 249 void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t * opt) const
 250 {
 251   // calculates the angular coverage in theta and phi of a single crystal in a EMC(=PHOS) module
 252
 253   Double_t conv ;
 254   if ( opt == Radian() )
 255     conv = 1. ;
 256   else if ( opt == Degre() )
 257     conv = 180. / TMath::Pi() ;
 258   else {
 259     cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ;
 260     conv = 1. ;
 261       }
 262
 263   Float_t y0  = GetIPtoCrystalSurface() ;
 264   theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
 265   phi   = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
 266 }
 267
 268
 269 //____________________________________________________________________________
 270 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
 271 {
 272   // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
 273
 274   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
 275   TVector3 localposition ;
 276
 277   tmpPHOS->GetLocalPosition(gpos) ;
 278
 279
 280   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal
 281     {  gpos.SetY( - GetIPtoCrystalSurface()) ;
 282
 283     }
 284   else
 285     { // it is a CPV
 286       gpos.SetY(- GetIPtoUpperCPVsurface()  ) ;
 287     }
 288
 289   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
 290   Double_t const kRADDEG = 180.0 / kPI ;
 291   Float_t rphi          = phi / kRADDEG ;
 292
 293   TRotation rot ;
 294   rot.RotateZ(-rphi) ; // a rotation around Z by angle
 295
 296   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
 297   gpos.Transform(rot) ; // rotate the baby
 298
 299 }
 300
 301 //____________________________________________________________________________
 302 void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const
 303 {
 304   // Calculates the coordinates of a RecPoint in the ALICE global coordinate system
 305
 306   AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
 307   TVector3 localposition ;
 308   tmpPHOS->GetLocalPosition(gpos) ;
 309
 310
 311   if ( tmpPHOS->IsEmc() ) // it is a EMC crystal
 312     {  gpos.SetY( - GetIPtoCrystalSurface()  ) ;
 313     }
 314   else
 315     { // it is a CPV
 316           gpos.SetY(- GetIPtoUpperCPVsurface()  ) ;
 317     }
 318
 319   Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
 320   Double_t const kRADDEG = 180.0 / kPI ;
 321   Float_t rphi          = phi / kRADDEG ;
 322
 323   TRotation rot ;
 324   rot.RotateZ(-rphi) ; // a rotation around Z by angle
 325
 326   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
 327   gpos.Transform(rot) ; // rotate the baby
 328 }
 329
 330 //____________________________________________________________________________
 331 void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x) const
 332 {
 333   // calculates the impact coordinates on PHOS of a neutral particle
 334   // emitted in the direction theta and phi in the ALICE global coordinate system
 335
 336   // searches for the PHOS EMC module
 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 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId) const
 358 {
 359   // Converts the relative numbering into the absolute numbering
 360   // EMCA crystals:
 361   //  AbsId = from 1 to fNModules * fNPhi * fNZ
 362   // CPV pad:
 363   //  AbsId = from N(total PHOS crystals) + 1
 364   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
 365
 366   Bool_t rv = kTRUE ;
 367
 368   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 ) * GetNZ()                   // the offset along phi
 372       +   relid[3] ;                                 // the offset along z
 373   }
 374   else { // it is a CPV pad
 375     AbsId =    GetNPhi() * GetNZ() *  GetNModules()         // the offset to separate EMCA crystals from CPV pads
 376       + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()   // the pads offset of PHOS modules
 377       + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                             // the pads offset of a CPV row
 378       +   relid[3] ;                                                         // the column number
 379   }
 380
 381   return rv ;
 382 }
 383
 384 //____________________________________________________________________________
 385
 386 void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) const
 387 {
 388   // Converts the absolute numbering into the global ALICE coordinate system
 389
 390
 391     Int_t relid[4] ;
 392
 393     AbsToRelNumbering(id , relid) ;
 394
 395     Int_t phosmodule = relid[0] ;
 396
 397     Float_t y0 = 0 ;
 398
 399     if ( relid[1] == 0 )  // it is a PbW04 crystal
 400       y0 =  - GetIPtoCrystalSurface() ;
 401     else
 402       y0 =  - GetIPtoUpperCPVsurface() ;
 403
 404     Float_t x, z ;
 405     RelPosInModule(relid, x, z) ;
 406
 407     pos.SetX(x) ;
 408     pos.SetZ(z) ;
 409     pos.SetY(y0) ;
 410
 411     Float_t phi           = GetPHOSAngle( phosmodule) ;
 412     Double_t const kRADDEG = 180.0 / kPI ;
 413     Float_t rphi          = phi / kRADDEG ;
 414
 415     TRotation rot ;
 416     rot.RotateZ(-rphi) ; // a rotation around Z by angle
 417
 418     TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
 419
 420     pos.Transform(rot) ; // rotate the baby
 421 }
 422
 423 //____________________________________________________________________________
 424 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const
 425 {
 426   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
 427   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
 428
 429   Int_t row        = relid[2] ; //offset along x axiz
 430   Int_t column     = relid[3] ; //offset along z axiz
 431
 432
 433   if ( relid[1] == 0 ) { // its a PbW04 crystal
 434     x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCrystalSize(0) ; // position ox Xtal with respect
 435     z =   ( GetNZ()  /2. - column + 0.5 ) *  GetCrystalSize(2) ; // of center of PHOS module
 436   }
 437   else  {
 438     x = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
 439     z =   ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
 440   }
 441 }