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 #include <TGeoManager.h>
  35 #include <TGeoMatrix.h>
  36
  37 // --- Standard library ---
  38
  39 // --- AliRoot header files ---
  40 #include "AliLog.h"
  41 #include "AliPHOSGeometry.h"
  42 #include "AliPHOSEMCAGeometry.h"
  43 #include "AliPHOSRecPoint.h"
  44
  45 ClassImp(AliPHOSGeometry)
  46
  47 // these initialisations are needed for a singleton
  48 AliPHOSGeometry  * AliPHOSGeometry::fgGeom = 0 ;
  49 Bool_t             AliPHOSGeometry::fgInit = kFALSE ;
  50
  51 //____________________________________________________________________________
  52 AliPHOSGeometry::AliPHOSGeometry() :
  53                     fNModules(0),
  54                     fAngle(0.f),
  55                     fPHOSAngle(0),
  56                     fIPtoUpperCPVsurface(0),
  57                     fCrystalShift(0),
  58                     fCryCellShift(0),
  59                     fRotMatrixArray(0),
  60                     fGeometryEMCA(0),
  61                     fGeometryCPV(0),
  62                     fGeometrySUPP(0)
  63 {
  64     // default ctor
  65     // must be kept public for root persistency purposes, but should never be called by the outside world
  66     fgGeom          = 0 ;
  67 }
  68
  69 //____________________________________________________________________________
  70 AliPHOSGeometry::AliPHOSGeometry(const AliPHOSGeometry & rhs)
  71                     : AliGeometry(rhs),
  72                       fNModules(rhs.fNModules),
  73                       fAngle(rhs.fAngle),
  74                       fPHOSAngle(0),
  75                       fIPtoUpperCPVsurface(rhs.fIPtoUpperCPVsurface),
  76                       fCrystalShift(rhs.fCrystalShift),
  77                       fCryCellShift(rhs.fCryCellShift),
  78                       fRotMatrixArray(0),
  79                       fGeometryEMCA(0),
  80                       fGeometryCPV(0),
  81                       fGeometrySUPP(0)
  82 {
  83   Fatal("cpy ctor", "not implemented") ;
  84 }
  85
  86 //____________________________________________________________________________
  87 AliPHOSGeometry::AliPHOSGeometry(const Text_t* name, const Text_t* title)
  88                   : AliGeometry(name, title),
  89                     fNModules(0),
  90                     fAngle(0.f),
  91                     fPHOSAngle(0),
  92                     fIPtoUpperCPVsurface(0),
  93                     fCrystalShift(0),
  94                     fCryCellShift(0),
  95                     fRotMatrixArray(0),
  96                     fGeometryEMCA(0),
  97                     fGeometryCPV(0),
  98                     fGeometrySUPP(0)
  99 {
 100   // ctor only for internal usage (singleton)
 101   Init() ;
 102   fgGeom = this;
 103 }
 104
 105 //____________________________________________________________________________
 106 AliPHOSGeometry::~AliPHOSGeometry(void)
 107 {
 108   // dtor
 109
 110   if (fRotMatrixArray) fRotMatrixArray->Delete() ;
 111   if (fRotMatrixArray) delete fRotMatrixArray ;
 112   if (fPHOSAngle     ) delete[] fPHOSAngle ;
 113 }
 114
 115 //____________________________________________________________________________
 116 void AliPHOSGeometry::Init(void)
 117 {
 118   // Initializes the PHOS parameters :
 119   //  IHEP is the Protvino CPV (cathode pad chambers)
 120
 121   TString test(GetName()) ;
 122   if (test != "IHEP" && test != "noCPV") {
 123     AliFatal(Form("%s is not a known geometry (choose among IHEP)",
 124                   test.Data() )) ;
 125   }
 126
 127   fgInit     = kTRUE ;
 128
 129   fNModules     = 5;
 130   fAngle        = 20;
 131
 132   fGeometryEMCA = new AliPHOSEMCAGeometry();
 133
 134   fGeometryCPV  = new AliPHOSCPVGeometry ();
 135
 136   fGeometrySUPP = new AliPHOSSupportGeometry();
 137
 138   fPHOSAngle = new Float_t[fNModules] ;
 139
 140   Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
 141
 142   fPHOSParams[0] =  TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2.,
 143                                (Double_t)(emcParams[0] - (emcParams[1]-emcParams[0])*
 144                                           fGeometryCPV->GetCPVBoxSize(1)/2/emcParams[3]));
 145   fPHOSParams[1] = emcParams[1] ;
 146   fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
 147   fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
 148
 149   fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
 150
 151   //calculate offset to crystal surface
 152   Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
 153   Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
 154   Float_t* splate = fGeometryEMCA->GetSupportPlateHalfSize();
 155   Float_t * crystal = fGeometryEMCA->GetCrystalHalfSize() ;
 156   Float_t * pin = fGeometryEMCA->GetAPDHalfSize() ;
 157   Float_t * preamp = fGeometryEMCA->GetPreampHalfSize() ;
 158   fCrystalShift=-inthermo[1]+strip[1]+splate[1]+crystal[1]-fGeometryEMCA->GetAirGapLed()/2.+pin[1]+preamp[1] ;
 159   fCryCellShift=crystal[1]-(fGeometryEMCA->GetAirGapLed()-2*pin[1]-2*preamp[1])/2;
 160
 161   Int_t index ;
 162   for ( index = 0; index < fNModules; index++ )
 163     fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()
 164
 165   fRotMatrixArray = new TObjArray(fNModules) ;
 166
 167   // Geometry parameters are calculated
 168
 169   SetPHOSAngles();
 170   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
 171   Float_t r = GetIPtoOuterCoverDistance() + fPHOSParams[3] - GetCPVBoxSize(1) ;
 172   for (Int_t iModule=0; iModule<fNModules; iModule++) {
 173     fModuleCenter[iModule][0] = r * TMath::Sin(fPHOSAngle[iModule] / kRADDEG );
 174     fModuleCenter[iModule][1] =-r * TMath::Cos(fPHOSAngle[iModule] / kRADDEG );
 175     fModuleCenter[iModule][2] = 0.;
 176
 177     fModuleAngle[iModule][0][0] =  90;
 178     fModuleAngle[iModule][0][1] =   fPHOSAngle[iModule];
 179     fModuleAngle[iModule][1][0] =   0;
 180     fModuleAngle[iModule][1][1] =   0;
 181     fModuleAngle[iModule][2][0] =  90;
 182     fModuleAngle[iModule][2][1] = 270 + fPHOSAngle[iModule];
 183   }
 184
 185 }
 186
 187 //____________________________________________________________________________
 188 AliPHOSGeometry *  AliPHOSGeometry::GetInstance()
 189 {
 190   // Returns the pointer of the unique instance; singleton specific
 191
 192   return static_cast<AliPHOSGeometry *>( fgGeom ) ;
 193 }
 194
 195 //____________________________________________________________________________
 196 AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title)
 197 {
 198   // Returns the pointer of the unique instance
 199   // Creates it with the specified options (name, title) if it does not exist yet
 200
 201   AliPHOSGeometry * rv = 0  ;
 202   if ( fgGeom == 0 ) {
 203     if ( strcmp(name,"") == 0 )
 204       rv = 0 ;
 205     else {
 206       fgGeom = new AliPHOSGeometry(name, title) ;
 207       if ( fgInit )
 208         rv = (AliPHOSGeometry * ) fgGeom ;
 209       else {
 210         rv = 0 ;
 211         delete fgGeom ;
 212         fgGeom = 0 ;
 213       }
 214     }
 215   }
 216   else {
 217     if ( strcmp(fgGeom->GetName(), name) != 0 )
 218       ::Error("GetInstance", "Current geometry is %s. You cannot call %s",
 219                       fgGeom->GetName(), name) ;
 220     else
 221       rv = (AliPHOSGeometry *) fgGeom ;
 222   }
 223   return rv ;
 224 }
 225
 226 //____________________________________________________________________________
 227 void AliPHOSGeometry::SetPHOSAngles()
 228 {
 229   // Calculates the position of the PHOS modules in ALICE global coordinate system
 230   // in ideal geometry
 231
 232   Double_t const kRADDEG = 180.0 / TMath::Pi() ;
 233   Float_t pphi =  2 * TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
 234   pphi *= kRADDEG ;
 235   if (pphi > fAngle){
 236     AliError(Form("PHOS modules overlap!\n pphi = %f fAngle = %f",
 237                   pphi, fAngle));
 238
 239   }
 240   pphi = fAngle;
 241
 242   for( Int_t i = 1; i <= fNModules ; i++ ) {
 243     Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
 244     fPHOSAngle[i-1] = -  angle ;
 245   }
 246 }
 247
 248 //____________________________________________________________________________
 249 Bool_t AliPHOSGeometry::AbsToRelNumbering(Int_t absId, Int_t * relid) const
 250 {
 251   // Converts the absolute numbering into the following array
 252   //  relid[0] = PHOS Module number 1:fNModules
 253   //  relid[1] = 0 if PbW04
 254   //           = -1 if CPV
 255   //  relid[2] = Row number inside a PHOS module
 256   //  relid[3] = Column number inside a PHOS module
 257
 258   Bool_t rv  = kTRUE ;
 259   Float_t id = absId ;
 260
 261   Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ;
 262
 263   if ( phosmodulenumber >  GetNModules() ) { // it is a CPV pad
 264
 265     id -=  GetNPhi() * GetNZ() *  GetNModules() ;
 266     Float_t nCPV  = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ;
 267     relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ;
 268     relid[1] = -1 ;
 269     id -= ( relid[0] - 1 ) * nCPV ;
 270     relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ;
 271     relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ;
 272   }
 273   else { // it is a PW04 crystal
 274
 275     relid[0] = phosmodulenumber ;
 276     relid[1] = 0 ;
 277     id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ;
 278     relid[2] = (Int_t)TMath::Ceil( id / GetNZ() )  ;
 279     relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ;
 280   }
 281   return rv ;
 282 }
 283 //____________________________________________________________________________
 284 void AliPHOSGeometry::GetGlobal(const AliRecPoint* , TVector3 & ) const
 285 {
 286   AliFatal(Form("Please use GetGlobalPHOS(recPoint,gpos) instead of GetGlobal!"));
 287 }
 288
 289 //____________________________________________________________________________
 290 void AliPHOSGeometry::GetGlobalPHOS(const AliPHOSRecPoint* recPoint, TVector3 & gpos) const
 291 {
 292   // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system
 293
 294   const AliPHOSRecPoint * tmpPHOS = recPoint ;
 295   TVector3 localposition ;
 296
 297   tmpPHOS->GetLocalPosition(gpos) ;
 298
 299   if (!gGeoManager){
 300     AliFatal("Geo manager not initialized\n");
 301   }
 302   //construct module name
 303   char path[100] ;
 304   Double_t dy ;
 305   if(tmpPHOS->IsEmc()){
 306     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",tmpPHOS->GetPHOSMod()) ;
 307 //    sprintf(path,"/ALIC_1/PHOS_%d",tmpPHOS->GetPHOSMod()) ;
 308     dy=fCrystalShift ;
 309   }
 310   else{
 311     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",tmpPHOS->GetPHOSMod()) ;
 312     dy= GetCPVBoxSize(1)/2. ; //center of CPV module
 313   }
 314   Double_t pos[3]={gpos.X(),gpos.Y()-dy,gpos.Z()} ;
 315   if(tmpPHOS->IsEmc())
 316     pos[2]=-pos[2] ; //Opposite z directions in EMC matrix and local frame!!!
 317   Double_t posC[3];
 318   //now apply possible shifts and rotations
 319   if (!gGeoManager->cd(path)){
 320     AliFatal("Geo manager can not find path \n");
 321   }
 322   TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 323   if (m){
 324      m->LocalToMaster(pos,posC);
 325   }
 326   else{
 327     AliFatal("Geo matrixes are not loaded \n") ;
 328   }
 329   gpos.SetXYZ(posC[0],posC[1],posC[2]) ;
 330
 331 }
 332 //____________________________________________________________________________
 333 void AliPHOSGeometry::ImpactOnEmc(Double_t * vtx, Double_t theta, Double_t phi,
 334                                   Int_t & moduleNumber, Double_t & z, Double_t & x) const
 335 {
 336   // calculates the impact coordinates on PHOS of a neutral particle
 337   // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
 338   TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
 339   TVector3 v(vtx[0],vtx[1],vtx[2]) ;
 340
 341   if (!gGeoManager){
 342     AliFatal("Geo manager not initialized\n");
 343   }
 344
 345   for(Int_t imod=1; imod<=GetNModules() ; imod++){
 346     //create vector from (0,0,0) to center of crystal surface of imod module
 347     Double_t tmp[3]={0.,-fCrystalShift,0.} ;
 348
 349     char path[100] ;
 350     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",imod) ;
 351     if (!gGeoManager->cd(path)){
 352       AliFatal("Geo manager can not find path \n");
 353     }
 354     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 355     Double_t posG[3]={0.,0.,0.} ;
 356     if (m) m->LocalToMaster(tmp,posG);
 357     TVector3 n(posG[0],posG[1],posG[2]) ;
 358     Double_t direction=n.Dot(p) ;
 359     if(direction<=0.)
 360       continue ; //momentum directed FROM module
 361     Double_t fr = (n.Mag2()-n.Dot(v))/direction ;
 362     //Calculate direction in module plain
 363     n-=v+fr*p ;
 364     n*=-1. ;
 365     Float_t * sz = fGeometryEMCA->GetInnerThermoHalfSize() ; //Wery close to the zise of the Xtl set
 366     if(TMath::Abs(TMath::Abs(n.Z())<sz[2]) && n.Pt()<sz[0]){
 367       moduleNumber = imod ;
 368       z=n.Z() ;
 369       x=TMath::Sign(n.Pt(),n.X()) ;
 370       //no need to return to local system since we calcilated distance from module center
 371       //and tilts can not be significant.
 372       return ;
 373     }
 374   }
 375   //Not in acceptance
 376   x=0; z=0 ;
 377   moduleNumber=0 ;
 378
 379 }
 380
 381 //____________________________________________________________________________
 382 Bool_t  AliPHOSGeometry::Impact(const TParticle * particle) const
 383 {
 384   // Tells if a particle enters PHOS
 385   Bool_t in=kFALSE;
 386   Int_t moduleNumber=0;
 387   Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ;
 388   Double_t z,x;
 389   ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x);
 390   if(moduleNumber!=0)
 391     in=kTRUE;
 392   return in;
 393 }
 394
 395 //____________________________________________________________________________
 396 Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  absId) const
 397 {
 398   // Converts the relative numbering into the absolute numbering
 399   // EMCA crystals:
 400   //  absId = from 1 to fNModules * fNPhi * fNZ
 401   // CPV pad:
 402   //  absId = from N(total PHOS crystals) + 1
 403   //          to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ
 404
 405   Bool_t rv = kTRUE ;
 406
 407   if ( relid[1] ==  0 ) {                            // it is a Phos crystal
 408     absId =
 409       ( relid[0] - 1 ) * GetNPhi() * GetNZ()         // the offset of PHOS modules
 410       + ( relid[2] - 1 ) * GetNZ()                   // the offset along phi
 411       +   relid[3] ;                                 // the offset along z
 412   }
 413   else { // it is a CPV pad
 414     absId =    GetNPhi() * GetNZ() *  GetNModules()         // the offset to separate EMCA crystals from CPV pads
 415       + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ()   // the pads offset of PHOS modules
 416       + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ()                             // the pads offset of a CPV row
 417       +   relid[3] ;                                                         // the column number
 418   }
 419
 420   return rv ;
 421 }
 422
 423 //____________________________________________________________________________
 424 void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const
 425 {
 426   // Converts the absolute numbering into the global ALICE coordinate system
 427
 428   if (!gGeoManager){
 429     AliFatal("Geo manager not initialized\n");
 430   }
 431
 432   Int_t relid[4] ;
 433
 434   AbsToRelNumbering(id , relid) ;
 435
 436   //construct module name
 437   char path[100] ;
 438   if(relid[1]==0){ //this is EMC
 439
 440     Double_t ps[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal
 441     Double_t psC[3]={0.0,0.0,0.}; //Global position
 442
 443     //Shift and possibly apply misalignment corrections
 444     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
 445     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
 446     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
 447                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
 448     Int_t cellraw= relid[3]%nCellsZInStrip ;
 449     if(cellraw==0)cellraw=nCellsZInStrip ;
 450     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ;
 451     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
 452             relid[0],strip,cell) ;
 453     if (!gGeoManager->cd(path)){
 454       AliFatal("Geo manager can not find path \n");
 455     }
 456     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 457     if (m) m->LocalToMaster(ps,psC);
 458     else{
 459       AliFatal("Geo matrixes are not loaded \n") ;
 460     }
 461     pos.SetXYZ(psC[0],psC[1],psC[2]) ;
 462   }
 463   else{
 464     //first calculate position with respect to CPV plain
 465     Int_t row        = relid[2] ; //offset along x axis
 466     Int_t column     = relid[3] ; //offset along z axis
 467     Double_t ps[3]= {0.0,GetCPVBoxSize(1)/2.,0.}; //Position on top of CPV
 468     Double_t psC[3]={0.0,0.0,0.}; //Global position
 469     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
 470     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
 471
 472     //now apply possible shifts and rotations
 473     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
 474     if (!gGeoManager->cd(path)){
 475       AliFatal("Geo manager can not find path \n");
 476     }
 477     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 478     if (m) m->LocalToMaster(ps,psC);
 479     else{
 480       AliFatal("Geo matrixes are not loaded \n") ;
 481     }
 482     pos.SetXYZ(psC[0],psC[1],-psC[2]) ;
 483   }
 484 }
 485
 486 //____________________________________________________________________________
 487 void AliPHOSGeometry::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const
 488 {
 489   // converts local PHOS-module (x, z) coordinates to absId
 490
 491   //find Global position
 492   if (!gGeoManager){
 493     AliFatal("Geo manager not initialized\n");
 494   }
 495   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
 496   Double_t posG[3] ;
 497   char path[100] ;
 498   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
 499   if (!gGeoManager->cd(path)){
 500     AliFatal("Geo manager can not find path \n");
 501   }
 502   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
 503   if (mPHOS){
 504      mPHOS->LocalToMaster(posL,posG);
 505   }
 506   else{
 507     AliFatal("Geo matrixes are not loaded \n") ;
 508   }
 509
 510   Int_t relid[4] ;
 511   gGeoManager->FindNode(posG[0],posG[1],posG[2]) ;
 512   //Check that path contains PSTR and extract strip number
 513   TString cpath(gGeoManager->GetPath()) ;
 514   Int_t indx = cpath.Index("PCEL") ;
 515   if(indx==-1){ //for the few events when particle hits between srips use ideal geometry
 516     relid[0] = module ;
 517     relid[1] = 0 ;
 518     relid[2] = static_cast<Int_t>(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) );
 519     relid[3] = static_cast<Int_t>(TMath::Ceil(-z/ GetCellStep() + GetNZ()   / 2.) ) ;
 520     if(relid[2]<1)relid[2]=1 ;
 521     if(relid[3]<1)relid[3]=1 ;
 522     if(relid[2]>GetNPhi())relid[2]=GetNPhi() ;
 523     if(relid[3]>GetNZ())relid[3]=GetNZ() ;
 524     RelToAbsNumbering(relid,absId) ;
 525   }
 526   else{
 527     Int_t indx2 = cpath.Index("/",indx) ;
 528     if(indx2==-1)
 529       indx2=cpath.Length() ;
 530     TString cell=cpath(indx+5,indx2-indx-5) ;
 531     Int_t icell=cell.Atoi() ;
 532     indx = cpath.Index("PSTR") ;
 533     indx2 = cpath.Index("/",indx) ;
 534     TString strip=cpath(indx+5,indx2-indx-5) ;
 535     Int_t iStrip = strip.Atoi() ;
 536
 537     Int_t row = fGeometryEMCA->GetNStripZ() - (iStrip - 1) % (fGeometryEMCA->GetNStripZ()) ;
 538     Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fGeometryEMCA->GetNStripZ())) -1 ;
 539
 540     // Absid for 8x2-strips. Looks nice :)
 541     absId = (module-1)*GetNCristalsInModule() +
 542                   row * 2 + (col*fGeometryEMCA->GetNCellsXInStrip() + (icell - 1) / 2)*GetNZ() - (icell & 1 ? 1 : 0);
 543
 544   }
 545
 546 }
 547
 548 //____________________________________________________________________________
 549 void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const
 550 {
 551   // Converts the relative numbering into the local PHOS-module (x, z) coordinates
 552   // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
 553
 554
 555   if (!gGeoManager){
 556     AliFatal("Geo manager not initialized\n");
 557   }
 558   //construct module name
 559   char path[100] ;
 560   if(relid[1]==0){ //this is PHOS
 561
 562 //   Calculations using ideal geometry (obsolete)
 563 //    x = - ( GetNPhi()/2. - relid[2]    + 0.5 ) *  GetCellStep() ; // position of Xtal with respect
 564 //    z = - ( GetNZ()  /2. - relid[3] + 0.5 ) *  GetCellStep() ; // of center of PHOS module
 565
 566     Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal
 567     Double_t posC[3]={0.0,0.0,0.}; //Global position
 568
 569     //Shift and possibly apply misalignment corrections
 570     Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ;
 571     Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ;
 572 //    Int_t strip=1+(relid[3]-1)/fGeometryEMCA->GetNCellsZInStrip()+((relid[2]-1)/nCellsInStrip)*fGeometryEMCA->GetNStripZ() ;
 573     Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()-
 574                 (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ;
 575     Int_t cellraw= relid[3]%nCellsZInStrip ;
 576     if(cellraw==0)cellraw=nCellsZInStrip ;
 577     Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ;
 578     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d",
 579             relid[0],strip,cell) ;
 580     if (!gGeoManager->cd(path)){
 581       AliFatal("Geo manager can not find path \n");
 582     }
 583     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 584     if (m) m->LocalToMaster(pos,posC);
 585     else{
 586       AliFatal("Geo matrixes are not loaded \n") ;
 587     }
 588     //    printf("Local: x=%f, y=%f, z=%f \n",pos[0],pos[1],pos[2]) ;
 589     //    printf("   gl: x=%f, y=%f, z=%f \n",posC[0],posC[1],posC[2]) ;
 590     //Return to PHOS local system
 591     Double_t posL[3]={posC[0],posC[1],posC[2]};
 592     sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",relid[0]) ;
 593     //    sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
 594     if (!gGeoManager->cd(path)){
 595       AliFatal("Geo manager can not find path \n");
 596     }
 597     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
 598     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
 599     else{
 600       AliFatal("Geo matrixes are not loaded \n") ;
 601     }
 602 //printf("RelPosInMod: posL=[%f,%f,%f]\n",posL[0],posL[1],posL[2]) ;
 603 //printf("old: x=%f, z=%f \n",x,z);
 604     x=posL[0] ;
 605     z=-posL[2];
 606     return ;
 607   }
 608   else{//CPV
 609     //first calculate position with respect to CPV plain
 610     Int_t row        = relid[2] ; //offset along x axis
 611     Int_t column     = relid[3] ; //offset along z axis
 612     Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit
 613     Double_t posC[3]={0.0,0.0,0.}; //Global position
 614     //    x = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
 615     //    z = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
 616     pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row    - 0.5 ) * GetPadSizePhi()  ; // position of pad  with respect
 617     pos[2] = - ( GetNumberOfCPVPadsZ()  /2. - column - 0.5 ) * GetPadSizeZ()  ; // of center of PHOS module
 618
 619     //now apply possible shifts and rotations
 620     sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ;
 621     if (!gGeoManager->cd(path)){
 622       AliFatal("Geo manager can not find path \n");
 623     }
 624     TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
 625     if (m) m->LocalToMaster(pos,posC);
 626     else{
 627       AliFatal("Geo matrixes are not loaded \n") ;
 628     }
 629     //Return to PHOS local system
 630     Double_t posL[3]={0.,0.,0.,} ;
 631     sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ;
 632     if (!gGeoManager->cd(path)){
 633       AliFatal("Geo manager can not find path \n");
 634     }
 635     TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
 636     if (mPHOS) mPHOS->MasterToLocal(posC,posL);
 637     else{
 638       AliFatal("Geo matrixes are not loaded \n") ;
 639     }
 640     x=posL[0] ;
 641     z=posL[1];
 642     return ;
 643
 644   }
 645
 646 }
 647
 648 //____________________________________________________________________________
 649
 650 void AliPHOSGeometry::GetModuleCenter(TVector3& center,
 651                                       const char *det,
 652                                       Int_t module) const
 653 {
 654   // Returns a position of the center of the CPV or EMC module
 655   // in ideal (not misaligned) geometry
 656   Float_t rDet = 0.;
 657   if      (strcmp(det,"CPV") == 0) rDet  = GetIPtoCPVDistance   ();
 658   else if (strcmp(det,"EMC") == 0) rDet  = GetIPtoCrystalSurface();
 659   else
 660     AliFatal(Form("Wrong detector name %s",det));
 661
 662   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
 663   angle *= TMath::Pi()/180;
 664   angle += 3*TMath::Pi()/2.;
 665   center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
 666 }
 667
 668 //____________________________________________________________________________
 669
 670 void AliPHOSGeometry::Global2Local(TVector3& localPosition,
 671                                    const TVector3& globalPosition,
 672                                    Int_t module) const
 673 {
 674   // Transforms a global position of the rec.point to the local coordinate system
 675   //Return to PHOS local system
 676   Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ;
 677   Double_t posL[3]={0.,0.,0.} ;
 678   char path[100] ;
 679   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ;
 680 //  sprintf(path,"/ALIC_1/PHOS_%d",module) ;
 681   if (!gGeoManager->cd(path)){
 682     AliFatal("Geo manager can not find path \n");
 683   }
 684   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
 685   if (mPHOS) mPHOS->MasterToLocal(posG,posL);
 686   else{
 687     AliFatal("Geo matrixes are not loaded \n") ;
 688   }
 689   localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ;
 690
 691 /*
 692   Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees
 693   angle *= TMath::Pi()/180;
 694   angle += 3*TMath::Pi()/2.;
 695   localPosition = globalPosition;
 696   localPosition.RotateZ(-angle);
 697 */
 698 }
 699 //____________________________________________________________________________
 700 void AliPHOSGeometry::Local2Global(Int_t mod, Float_t x, Float_t z,
 701                                    TVector3& globalPosition) const
 702 {
 703   char path[100] ;
 704   sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ;
 705 //  sprintf(path,"/ALIC_1/PHOS_%d",mod) ;
 706   if (!gGeoManager->cd(path)){
 707     AliFatal("Geo manager can not find path \n");
 708   }
 709   Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!!
 710   Double_t posG[3] ;
 711   TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix();
 712   if (mPHOS){
 713      mPHOS->LocalToMaster(posL,posG);
 714   }
 715   else{
 716     AliFatal("Geo matrixes are not loaded \n") ;
 717   }
 718   globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ;
 719 }
 720 //____________________________________________________________________________
 721 void AliPHOSGeometry::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const {
 722   //Calculates vector pointing from vertex to current poisition in module local frame
 723   //Note that PHOS local system and ALICE global have opposite z directions
 724
 725   Global2Local(vInc,vtx,module) ;
 726   vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ;
 727 }