[u/mrichter/AliRoot.git] / PHOS / AliPHOSGeometry.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//_________________________________________________________________________
// Geometry class  for PHOS : singleton  
// PHOS consists of the electromagnetic calorimeter (EMCA)
// and a charged particle veto either in the Subatech's version (PPSD)
// or in the IHEP's one (CPV).
// The EMCA/PPSD/CPV modules are parametrized so that any configuration
// can be easily implemented 
// The title is used to identify the version of CPV used.
//                  
//*-- Author: Yves Schutz (SUBATECH)

// --- ROOT system ---

#include "TVector3.h"
#include "TRotation.h" 

// --- Standard library ---

#include <iostream.h>

// --- AliRoot header files ---

#include "AliPHOSGeometry.h"
#include "AliPHOSPpsdRecPoint.h"
#include "AliConst.h"

ClassImp(AliPHOSGeometry) ;

AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
Bool_t            AliPHOSGeometry::fgInit = kFALSE ;

//____________________________________________________________________________
AliPHOSGeometry::~AliPHOSGeometry(void)
{
  // dtor

  if (fRotMatrixArray) fRotMatrixArray->Delete() ; 
  if (fRotMatrixArray) delete fRotMatrixArray ; 
  if (fPHOSAngle     ) delete fPHOSAngle ; 
//    if (fGeometryEMCA  ) detete fGeometryEMCA;
//    if (fGeometryCPV   ) detete fGeometryCPV ;
}

//____________________________________________________________________________

void AliPHOSGeometry::Init(void)
{
  // Initializes the PHOS parameters

  if ( ((strcmp( fName, "default" )) == 0) || 
       ((strcmp( fName, "GPS2" ))    == 0) ||
       ((strcmp( fName, "IHEP" ))    == 0) ) {
    fgInit     = kTRUE ; 
                                             fGeometryEMCA = new AliPHOSEMCAGeometry();
    if ( ((strcmp( fName, "GPS2" ))  == 0) ) fGeometryCPV  = new AliPHOSPPSDGeometry();
    if ( ((strcmp( fName, "IHEP" ))  == 0) ) fGeometryCPV  = new AliPHOSCPVGeometry ();
                                             fGeometrySUPP = new AliPHOSSupportGeometry();
    fNModules = 5;
    fPHOSAngle = new Float_t[fNModules] ;
    Int_t index ;
    for ( index = 0; index < fNModules; index++ )
    fPHOSAngle[index] = 0.0 ; // Module position angles are set in CreateGeometry()

    this->SetPHOSAngles() ; 
    fRotMatrixArray = new TObjArray(fNModules) ; 
  }
 else {
   fgInit = kFALSE ; 
   cout << "PHOS Geometry setup: option not defined " << fName << endl ; 
 }
}

//____________________________________________________________________________
AliPHOSGeometry *  AliPHOSGeometry::GetInstance() 
{ 
  // Returns the pointer of the unique instance
  return (AliPHOSGeometry *) fgGeom ; 
}

//____________________________________________________________________________
AliPHOSGeometry *  AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title) 
{
  // Returns the pointer of the unique instance
  AliPHOSGeometry * rv = 0  ; 
  if ( fgGeom == 0 ) {
    if ( strcmp(name,"") == 0 ) 
      rv = 0 ;
    else {    
      fgGeom = new AliPHOSGeometry(name, title) ;
      if ( fgInit )
	rv = (AliPHOSGeometry * ) fgGeom ;
      else {
	rv = 0 ; 
	delete fgGeom ; 
	fgGeom = 0 ; 
      }
    }
  }
  else {
    if ( strcmp(fgGeom->GetName(), name) != 0 ) {
      cout << "AliPHOSGeometry <E> : current geometry is " << fgGeom->GetName() << endl
	   << "                      you cannot call     " << name << endl ; 
    }
    else
      rv = (AliPHOSGeometry *) fgGeom ; 
  } 
  return rv ; 
}

//____________________________________________________________________________
void AliPHOSGeometry::SetPHOSAngles() 
{ 
  // Calculates the position in ALICE of the PHOS modules
  
  Double_t const kRADDEG = 180.0 / kPI ;
  Float_t pphi =  TMath::ATan( GetOuterBoxSize(0)  / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ;
  pphi *= kRADDEG ;
  
  for( Int_t i = 1; i <= fNModules ; i++ ) {
    Float_t angle = pphi * 2 * ( i - fNModules / 2.0 - 0.5 ) ;
    fPHOSAngle[i-1] = -  angle ;
  } 
}

//____________________________________________________________________________
Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid)
{
  // Converts the absolute numbering into the following array/
  //  relid[0] = PHOS Module number 1:fNModules 
  //  relid[1] = 0 if PbW04
  //           = PPSD Module number 1:fNumberOfModulesPhi*fNumberOfModulesZ*2 (2->up and bottom level)
  //  relid[2] = Row number inside a PHOS or PPSD module
  //  relid[3] = Column number inside a PHOS or PPSD module

  Bool_t rv  = kTRUE ; 
  Float_t id = AbsId ;

  Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / ( GetNPhi() * GetNZ() ) ) ; 
  
  if ( phosmodulenumber >  GetNModules() ) { // it is a PPSD or CPV pad

    if      ( strcmp(fName,"GPS2") == 0 ) {
      id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
      Float_t tempo = 2 *  GetNumberOfModulesPhi() * GetNumberOfModulesZ() *  GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; 
      relid[0] = (Int_t)TMath::Ceil( id / tempo ) ; 
      id -= ( relid[0] - 1 ) * tempo ;
      relid[1] = (Int_t)TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ; 
      id -= ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ;
      relid[2] = (Int_t)TMath::Ceil( id / GetNumberOfPadsPhi() ) ;
      relid[3] = (Int_t) ( id - ( relid[2] - 1 )  * GetNumberOfPadsPhi() ) ; 
    }
    else if ( strcmp(fName,"IHEP") == 0 ) {
      id -=  GetNPhi() * GetNZ() *  GetNModules() ; 
      relid[0] = (Int_t) TMath::Ceil( id / ( GetNumberOfPadsPhi() * GetNumberOfPadsZ() ) ) ;
      relid[1] = 1 ;
      id -= ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()  ; 
      relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfPadsZ() ) ;
      relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfPadsZ() ) ; 
    }
  } 
  else { // its a PW04 crystal

    relid[0] = phosmodulenumber ;
    relid[1] = 0 ;
    id -= ( phosmodulenumber - 1 ) *  GetNPhi() * GetNZ() ; 
    relid[2] = (Int_t)TMath::Ceil( id / GetNPhi() ) ;
    relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNPhi() ) ; 
  } 
  return rv ; 
}

//____________________________________________________________________________  
void AliPHOSGeometry::EmcModuleCoverage(const Int_t mod, Double_t & tm, Double_t & tM, Double_t & pm, Double_t & pM, Option_t * opt) 
{
  // calculates the angular coverage in theta and phi of a EMC module

 Double_t conv ; 
  if ( opt == Radian() ) 
    conv = 1. ; 
  else if ( opt == Degre() )
    conv = 180. / TMath::Pi() ; 
  else {
    cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
    conv = 1. ;
      }

  Float_t phi =  GetPHOSAngle(mod) *  (TMath::Pi() / 180.)  ;  
  Float_t y0  =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
		  + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;  
  
  Double_t angle = TMath::ATan( GetCrystalSize(0)*GetNPhi() / (2 * y0) ) ;
  phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 230 and 310 deg.)
  Double_t max  = phi - angle ;
  Double_t min   = phi + angle ;
  pM = TMath::Max(max, min) * conv ;
  pm = TMath::Min(max, min) * conv ; 
  
  angle =  TMath::ATan( GetCrystalSize(2)*GetNZ() / (2 * y0) ) ;
  max  = TMath::Pi() / 2.  + angle ; // to follow the convention of the particle generator(PHOS is at 90 deg.)
  min  = TMath::Pi() / 2.  - angle ;
  tM = TMath::Max(max, min) * conv ;
  tm = TMath::Min(max, min) * conv ; 
 
}

//____________________________________________________________________________  
void AliPHOSGeometry::EmcXtalCoverage(Double_t & theta, Double_t & phi, Option_t * opt) 
{
  // calculates the angular coverage in theta and phi of a single crystal in a EMC module

  Double_t conv ; 
  if ( opt == Radian() ) 
    conv = 1. ; 
  else if ( opt == Degre() )
    conv = 180. / TMath::Pi() ; 
  else {
    cout << "<I>  AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ; 
    conv = 1. ;
      }

  Float_t y0   =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
    + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;  
  theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
  phi   = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
}
 

//____________________________________________________________________________
void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
{
  // Calculates the ALICE global coordinates of a RecPoint and the error matrix
 
  AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
  TVector3 localposition ;

  tmpPHOS->GetLocalPosition(gpos) ;


  if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
    {  gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() +
		    GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ;  

    }
  else
    { // it is a PPSD pad
      AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ;
      if (tmpPpsd->GetUp() ) // it is an upper module
	{
	  gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - 
		       GetLeadToMicro2Gap() - GetLeadConverterThickness() -  
		       GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 )  ) ; 
	} 
      else // it is a lower module
	gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; 
    }  

  Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
  Double_t const kRADDEG = 180.0 / kPI ;
  Float_t rphi          = phi / kRADDEG ; 
  
  TRotation rot ;
  rot.RotateZ(-rphi) ; // a rotation around Z by angle  
  
  TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
  gpos.Transform(rot) ; // rotate the baby 

}

//____________________________________________________________________________
void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const 
{
  // Calculates the ALICE global coordinates of a RecPoint 

  AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;  
  TVector3 localposition ;
  tmpPHOS->GetLocalPosition(gpos) ;


  if ( tmpPHOS->IsEmc() ) // it is a EMC crystal 
    {  gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() +
		    GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ;  
    }
  else
    { // it is a PPSD pad
      AliPHOSPpsdRecPoint * tmpPpsd = (AliPHOSPpsdRecPoint *) RecPoint ;
      if (tmpPpsd->GetUp() ) // it is an upper module
	{
	  gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - 
		       GetLeadToMicro2Gap() - GetLeadConverterThickness() -  
		       GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0 )  ) ; 
	} 
      else // it is a lower module
	gpos.SetY(-( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ) ; 
    }  

  Float_t phi           = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; 
  Double_t const kRADDEG = 180.0 / kPI ;
  Float_t rphi          = phi / kRADDEG ; 
  
  TRotation rot ;
  rot.RotateZ(-rphi) ; // a rotation around Z by angle  
  
  TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
  gpos.Transform(rot) ; // rotate the baby 
}

//____________________________________________________________________________
void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x) 
{
  // calculates the impact coordinates of a neutral particle  
  // emitted in direction theta and phi in ALICE

  // searches for the PHOS EMC module
  ModuleNumber = 0 ; 
  Double_t tm, tM, pm, pM ; 
  Int_t index = 1 ; 
  while ( ModuleNumber == 0 && index <= GetNModules() ) { 
    EmcModuleCoverage(index, tm, tM, pm, pM) ; 
    if ( (theta >= tm && theta <= tM) && (phi >= pm && phi <= pM ) ) 
      ModuleNumber = index ; 
    index++ ;    
  }
  if ( ModuleNumber != 0 ) {
    Float_t phi0 =  GetPHOSAngle(ModuleNumber) *  (TMath::Pi() / 180.) + 1.5 * TMath::Pi()  ;  
    Float_t y0  =  GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
      + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()  ;   
    Double_t angle = phi - phi0; 
    x = y0 * TMath::Tan(angle) ; 
    angle = theta - TMath::Pi() / 2 ; 
    z = y0 * TMath::Tan(angle) ; 
  }
}

//____________________________________________________________________________
Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t &  AbsId)
{
  // Converts the relative numbering into the absolute numbering
  //  AbsId = 1 to fNModules * fNPhi * fNZ                                                                       -> PbWO4
  //  AbsId = N(total PHOS crystals) +
  //          1 to fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ -> PPSD
  //  AbsId = N(total PHOS crystals) +
  //          1:fNModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ                                              -> CPV

  Bool_t rv = kTRUE ; 
 
  if      ( relid[1]  >  0 ) { // it is a PPSD pad
    AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from PPSD pads
      + ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ()   // the pads offset of PHOS modules 
                         * GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
      + ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PPSD modules 
      + ( relid[2] - 1 ) * GetNumberOfPadsPhi()                              // the pads offset of a PPSD row
      +   relid[3] ;                                                         // the column number
  } 

  else if ( relid[1] ==  0 ) { // it is a Phos crystal
    AbsId =
        ( relid[0] - 1 ) * GetNPhi() * GetNZ()                               // the offset of PHOS modules
      + ( relid[2] - 1 ) * GetNPhi()                                         // the offset of a xtal row
      +   relid[3] ;                                                         // the column number
  }

  else if ( relid[1] == -1 ) { // it is a CPV pad
    AbsId =    GetNPhi() * GetNZ() *  GetNModules()                          // the offset to separate EMCA crystals from CPV pads
      + ( relid[0] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ()         // the pads offset of PHOS modules 
      + ( relid[2] - 1 ) * GetNumberOfPadsZ()                                // the pads offset of a CPV row
      +   relid[3] ;                                                         // the column number
  }
  
  return rv ; 
}

//____________________________________________________________________________

void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) 
{
  // Converts the absolute numbering into the global ALICE coordinates
  
   if (id > 0) { 

  Int_t relid[4] ;
 
  AbsToRelNumbering(id , relid) ;

  Int_t phosmodule = relid[0] ; 

  Float_t y0 = 0 ; 

  if ( relid[1] == 0 ) // it is a PbW04 crystal 
  {  y0 =  -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
      + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness())  ;  
  }
  if ( relid[1] > 0 ) { // its a PPSD pad
    if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() ) // its an bottom module
     {
       y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0)  ;
     } 
    else // its an upper module
      y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
	-  GetLeadConverterThickness() -  GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ; 
  }

  Float_t x, z ; 
  RelPosInModule(relid, x, z) ; 

  pos.SetX(x) ;
  pos.SetZ(z) ;
  pos.SetY( TMath::Sqrt(x*x + z*z + y0*y0) ) ; 


   Float_t phi           = GetPHOSAngle( phosmodule) ; 
   Double_t const kRADDEG = 180.0 / kPI ;
   Float_t rphi          = phi / kRADDEG ; 

   TRotation rot ;
   rot.RotateZ(-rphi) ; // a rotation around Z by angle  
  
   TRotation dummy = rot.Invert() ;  // to transform from original frame to rotate frame
  
   pos.Transform(rot) ; // rotate the baby 
  }
  else {
 pos.SetX(0.);
 pos.SetY(0.);
 pos.SetZ(0.);
       }
} 

//____________________________________________________________________________
void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) 
{
  // Converts the relative numbering into the local PHOS-module (x, z) coordinates
  // Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
  
  Int_t ppsdmodule  ; 
  Float_t x0,z0;
  Int_t row        = relid[2] ; //offset along x axiz
  Int_t column     = relid[3] ; //offset along z axiz

  Float_t padsizeZ = GetPadSizeZ();
  Float_t padsizeX = GetPadSizePhi();

  if ( relid[1] == 0 ) { // its a PbW04 crystal 
    x = - ( GetNPhi()/2. - row    + 0.5 ) *  GetCrystalSize(0) ; // position ox Xtal with respect
    z =   ( GetNZ()  /2. - column + 0.5 ) *  GetCrystalSize(2) ; // of center of PHOS module  
  }  
  else  {    
    if ( relid[1] >  GetNumberOfModulesPhi() *  GetNumberOfModulesZ() )
      ppsdmodule =  relid[1]-GetNumberOfModulesPhi() *  GetNumberOfModulesZ(); 
    else
      ppsdmodule =  relid[1] ;
    Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ; 
    Int_t modcol = ppsdmodule -  ( modrow - 1 ) * GetNumberOfModulesPhi() ;     
    if ( ((strcmp( fName, "GPS2" ))  == 0) ) {
      x0 = (  GetNumberOfModulesPhi() / 2.  - modrow  + 0.5 ) * GetPPSDModuleSize(0) ;
      z0 = (  GetNumberOfModulesZ()   / 2.  - modcol  + 0.5 ) * GetPPSDModuleSize(2)  ;     
    } else {
      x0 = 0;
      z0 = 0;
    }
    x = - ( GetNumberOfPadsPhi()/2. - row    - 0.5 ) * padsizeX + x0 ; // position of pad  with respect
    z =   ( GetNumberOfPadsZ()  /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module  
  }
}
Commit	Line	Data
d15a28e7	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
b2a60966	16	/* $Id$ */
b2a60966	17
d15a28e7	18	//_________________________________________________________________________
b2a60966	19	// Geometry class for PHOS : singleton
a3dfe79c	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.
b2a60966	26	//
b2a60966	27	//*-- Author: Yves Schutz (SUBATECH)
d15a28e7	28
	29	// --- ROOT system ---
	30
	31	#include "TVector3.h"
	32	#include "TRotation.h"
	33
	34	// --- Standard library ---
	35
de9ec31b	36	#include <iostream.h>
d15a28e7	37
	38	// --- AliRoot header files ---
	39
	40	#include "AliPHOSGeometry.h"
	41	#include "AliPHOSPpsdRecPoint.h"
	42	#include "AliConst.h"
	43
9ec91567	44	ClassImp(AliPHOSGeometry) ;
d15a28e7	45
9ec91567	46	AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
282c5906	47	Bool_t AliPHOSGeometry::fgInit = kFALSE ;
9ec91567	48
d15a28e7	49	//____________________________________________________________________________
	50	AliPHOSGeometry::~AliPHOSGeometry(void)
	51	{
b2a60966	52	// dtor
b2a60966	53
52a36ffd	54	if (fRotMatrixArray) fRotMatrixArray->Delete() ;
	55	if (fRotMatrixArray) delete fRotMatrixArray ;
	56	if (fPHOSAngle ) delete fPHOSAngle ;
eb92d866	57	// if (fGeometryEMCA ) detete fGeometryEMCA;
eb92d866	58	// if (fGeometryCPV ) detete fGeometryCPV ;
52a36ffd	59	}
	60
	61	//____________________________________________________________________________
	62
	63	void AliPHOSGeometry::Init(void)
	64	{
	65	// Initializes the PHOS parameters
	66
52a36ffd	67	if ( ((strcmp( fName, "default" )) == 0) \|\|
	68	((strcmp( fName, "GPS2" )) == 0) \|\|
	69	((strcmp( fName, "IHEP" )) == 0) ) {
	70	fgInit = kTRUE ;
eb92d866	71	fGeometryEMCA = new AliPHOSEMCAGeometry();
	72	if ( ((strcmp( fName, "GPS2" )) == 0) ) fGeometryCPV = new AliPHOSPPSDGeometry();
	73	if ( ((strcmp( fName, "IHEP" )) == 0) ) fGeometryCPV = new AliPHOSCPVGeometry ();
ed19b2e1	74	fGeometrySUPP = new AliPHOSSupportGeometry();
52a36ffd	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	}
4697edca	126
52a36ffd	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	}
d15a28e7	140	}
	141
	142	//____________________________________________________________________________
92862013	143	Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid)
d15a28e7	144	{
b2a60966	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:fNumberOfModulesPhifNumberOfModulesZ2 (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
d15a28e7	151
d15a28e7	152	Bool_t rv = kTRUE ;
92862013	153	Float_t id = AbsId ;
d15a28e7	154
92862013	155	Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / ( GetNPhi() * GetNZ() ) ) ;
d15a28e7	156
52a36ffd	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	}
d15a28e7	177	}
	178	else { // its a PW04 crystal
	179
92862013	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() ) ;
d15a28e7	185	}
	186	return rv ;
	187	}
52a36ffd	188
9f616d61	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 ;
cf0c2bc1	195	if ( opt == Radian() )
9f616d61	196	conv = 1. ;
cf0c2bc1	197	else if ( opt == Degre() )
9f616d61	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.) ;
92862013	205	Float_t y0 = GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
9f616d61	206	+ GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness() ;
9f616d61	207
92862013	208	Double_t angle = TMath::ATan( GetCrystalSize(0)GetNPhi() / (2 y0) ) ;
9f616d61	209	phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 230 and 310 deg.)
92862013	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 ;
9f616d61	214
92862013	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 ;
9f616d61	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 ;
cf0c2bc1	229	if ( opt == Radian() )
9f616d61	230	conv = 1. ;
cf0c2bc1	231	else if ( opt == Degre() )
9f616d61	232	conv = 180. / TMath::Pi() ;
	233	else {
	234	cout << "<I> AliPHOSGeometry::EmcXtalCoverage : " << opt << " unknown option; result in radian " << endl ;
	235	conv = 1. ;
	236	}
	237
92862013	238	Float_t y0 = GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
9f616d61	239	+ GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness() ;
92862013	240	theta = 2 * TMath::ATan( GetCrystalSize(2) / (2 * y0) ) * conv ;
92862013	241	phi = 2 * TMath::ATan( GetCrystalSize(0) / (2 * y0) ) * conv ;
9f616d61	242	}
	243
	244
	245	//____________________________________________________________________________
52a36ffd	246	void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const
d15a28e7	247	{
b2a60966	248	// Calculates the ALICE global coordinates of a RecPoint and the error matrix
b2a60966	249
d15a28e7	250	AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
92862013	251	TVector3 localposition ;
d15a28e7	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
92862013	274	Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
	275	Double_t const kRADDEG = 180.0 / kPI ;
	276	Float_t rphi = phi / kRADDEG ;
d15a28e7	277
92862013	278	TRotation rot ;
92862013	279	rot.RotateZ(-rphi) ; // a rotation around Z by angle
d15a28e7	280
92862013	281	TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
92862013	282	gpos.Transform(rot) ; // rotate the baby
6ad0bfa0	283
d15a28e7	284	}
	285
	286	//____________________________________________________________________________
5cda30f6	287	void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const
d15a28e7	288	{
b2a60966	289	// Calculates the ALICE global coordinates of a RecPoint
b2a60966	290
d15a28e7	291	AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ;
92862013	292	TVector3 localposition ;
d15a28e7	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
92862013	313	Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ;
	314	Double_t const kRADDEG = 180.0 / kPI ;
	315	Float_t rphi = phi / kRADDEG ;
d15a28e7	316
92862013	317	TRotation rot ;
92862013	318	rot.RotateZ(-rphi) ; // a rotation around Z by angle
d15a28e7	319
92862013	320	TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
92862013	321	gpos.Transform(rot) ; // rotate the baby
d15a28e7	322	}
	323
	324	//____________________________________________________________________________
52a36ffd	325	void AliPHOSGeometry::ImpactOnEmc(const Double_t theta, const Double_t phi, Int_t & ModuleNumber, Double_t & z, Double_t & x)
d15a28e7	326	{
52a36ffd	327	// calculates the impact coordinates of a neutral particle
52a36ffd	328	// emitted in direction theta and phi in ALICE
d15a28e7	329
52a36ffd	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++ ;
d15a28e7	339	}
52a36ffd	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) ;
d15a28e7	348	}
d15a28e7	349	}
	350
	351	//____________________________________________________________________________
92862013	352	Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId)
d15a28e7	353	{
b2a60966	354	// Converts the relative numbering into the absolute numbering
52a36ffd	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
d15a28e7	360
	361	Bool_t rv = kTRUE ;
	362
52a36ffd	363	if ( relid[1] > 0 ) { // it is a PPSD pad
52a36ffd	364	AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from PPSD pads
92862013	365	+ ( relid[0] - 1 ) * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PHOS modules
d15a28e7	366	* GetNumberOfPadsPhi() * GetNumberOfPadsZ() * 2
92862013	367	+ ( relid[1] - 1 ) * GetNumberOfPadsPhi() * GetNumberOfPadsZ() // the pads offset of PPSD modules
92862013	368	+ ( relid[2] - 1 ) * GetNumberOfPadsPhi() // the pads offset of a PPSD row
52a36ffd	369	+ relid[3] ; // the column number
d15a28e7	370	}
52a36ffd	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
d15a28e7	377	}
d15a28e7	378
52a36ffd	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
d15a28e7	386	return rv ;
	387	}
	388
	389	//____________________________________________________________________________
	390
92862013	391	void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos )
d15a28e7	392	{
b2a60966	393	// Converts the absolute numbering into the global ALICE coordinates
b2a60966	394
92862013	395	if (id > 0) {
d15a28e7	396
92862013	397	Int_t relid[4] ;
d15a28e7	398
92862013	399	AbsToRelNumbering(id , relid) ;
d15a28e7	400
92862013	401	Int_t phosmodule = relid[0] ;
d15a28e7	402
92862013	403	Float_t y0 = 0 ;
9f616d61	404
92862013	405	if ( relid[1] == 0 ) // it is a PbW04 crystal
92862013	406	{ y0 = -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness()
9f616d61	407	+ GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ;
d15a28e7	408	}
92862013	409	if ( relid[1] > 0 ) { // its a PPSD pad
92862013	410	if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() ) // its an bottom module
d15a28e7	411	{
92862013	412	y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() / 2.0) ;
d15a28e7	413	}
d15a28e7	414	else // its an upper module
92862013	415	y0 = -( GetIPtoOuterCoverDistance() - GetMicromegas2Thickness() - GetLeadToMicro2Gap()
9f616d61	416	- GetLeadConverterThickness() - GetMicro1ToLeadGap() - GetMicromegas1Thickness() / 2.0) ;
d15a28e7	417	}
	418
	419	Float_t x, z ;
92862013	420	RelPosInModule(relid, x, z) ;
d15a28e7	421
9f616d61	422	pos.SetX(x) ;
9f616d61	423	pos.SetZ(z) ;
92862013	424	pos.SetY( TMath::Sqrt(xx + zz + y0*y0) ) ;
9f616d61	425
d15a28e7	426
d15a28e7	427
92862013	428	Float_t phi = GetPHOSAngle( phosmodule) ;
	429	Double_t const kRADDEG = 180.0 / kPI ;
	430	Float_t rphi = phi / kRADDEG ;
d15a28e7	431
92862013	432	TRotation rot ;
92862013	433	rot.RotateZ(-rphi) ; // a rotation around Z by angle
d15a28e7	434
92862013	435	TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame
d15a28e7	436
92862013	437	pos.Transform(rot) ; // rotate the baby
d15a28e7	438	}
	439	else {
	440	pos.SetX(0.);
	441	pos.SetY(0.);
	442	pos.SetZ(0.);
	443	}
	444	}
	445
	446	//____________________________________________________________________________
92862013	447	void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z)
d15a28e7	448	{
b2a60966	449	// Converts the relative numbering into the local PHOS-module (x, z) coordinates
52a36ffd	450	// Note: sign of z differs from that in the previous version (Yu.Kharlov, 12 Oct 2000)
b2a60966	451
92862013	452	Int_t ppsdmodule ;
a3dfe79c	453	Float_t x0,z0;
52a36ffd	454	Int_t row = relid[2] ; //offset along x axiz
52a36ffd	455	Int_t column = relid[3] ; //offset along z axiz
d15a28e7	456
52a36ffd	457	Float_t padsizeZ = GetPadSizeZ();
52a36ffd	458	Float_t padsizeX = GetPadSizePhi();
d15a28e7	459
92862013	460	if ( relid[1] == 0 ) { // its a PbW04 crystal
52a36ffd	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 {
92862013	465	if ( relid[1] > GetNumberOfModulesPhi() * GetNumberOfModulesZ() )
52a36ffd	466	ppsdmodule = relid[1]-GetNumberOfModulesPhi() * GetNumberOfModulesZ();
	467	else
	468	ppsdmodule = relid[1] ;
92862013	469	Int_t modrow = 1+(Int_t)TMath::Ceil( (Float_t)ppsdmodule / GetNumberOfModulesPhi()-1. ) ;
92862013	470	Int_t modcol = ppsdmodule - ( modrow - 1 ) * GetNumberOfModulesPhi() ;
a3dfe79c	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	}
52a36ffd	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	}
2f3366b6	481	}