#include "TVector3.h"
#include "TRotation.h"
#include "TParticle.h"
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include "AliLog.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSEMCAGeometry.h"
#include "AliPHOSRecPoint.h"
ClassImp(AliPHOSGeometry)
AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ;
Bool_t AliPHOSGeometry::fgInit = kFALSE ;
AliPHOSGeometry::AliPHOSGeometry() :
AliPHOSGeoUtils(),
fAngle(0.f),
fPHOSAngle(0),
fIPtoUpperCPVsurface(0),
fCrystalShift(0),
fCryCellShift(0),
fRotMatrixArray(0)
{
fgGeom = 0 ;
fPHOSParams[0] = 0.;
fPHOSParams[1] = 0.;
fPHOSParams[2] = 0.;
fPHOSParams[3] = 0.;
}
AliPHOSGeometry::AliPHOSGeometry(const AliPHOSGeometry & rhs)
: AliPHOSGeoUtils(rhs),
fAngle(rhs.fAngle),
fPHOSAngle(0),
fIPtoUpperCPVsurface(rhs.fIPtoUpperCPVsurface),
fCrystalShift(rhs.fCrystalShift),
fCryCellShift(rhs.fCryCellShift),
fRotMatrixArray(0)
{
Fatal("cpy ctor", "not implemented") ;
}
AliPHOSGeometry::AliPHOSGeometry(const Text_t* name, const Text_t* title)
: AliPHOSGeoUtils(name, title),
fAngle(0.f),
fPHOSAngle(0),
fIPtoUpperCPVsurface(0),
fCrystalShift(0),
fCryCellShift(0),
fRotMatrixArray(0)
{
Init() ;
fgGeom = this;
}
AliPHOSGeometry::~AliPHOSGeometry(void)
{
if (fRotMatrixArray) fRotMatrixArray->Delete() ;
if (fRotMatrixArray) delete fRotMatrixArray ;
if (fPHOSAngle ) delete[] fPHOSAngle ;
}
void AliPHOSGeometry::Init(void)
{
fgInit = kTRUE ;
fAngle = 20;
fPHOSAngle = new Float_t[fNModules] ;
const Float_t * emcParams = fGeometryEMCA->GetEMCParams() ;
fPHOSParams[0] = TMath::Max((Double_t)fGeometryCPV->GetCPVBoxSize(0)/2.,
(Double_t)(emcParams[0] - (emcParams[1]-emcParams[0])*
fGeometryCPV->GetCPVBoxSize(1)/2/emcParams[3]));
fPHOSParams[1] = emcParams[1] ;
fPHOSParams[2] = TMath::Max((Double_t)emcParams[2], (Double_t)fGeometryCPV->GetCPVBoxSize(2)/2.);
fPHOSParams[3] = emcParams[3] + fGeometryCPV->GetCPVBoxSize(1)/2. ;
fIPtoUpperCPVsurface = fGeometryEMCA->GetIPtoOuterCoverDistance() - fGeometryCPV->GetCPVBoxSize(1) ;
const Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ;
const Float_t * strip = fGeometryEMCA->GetStripHalfSize() ;
const Float_t * splate = fGeometryEMCA->GetSupportPlateHalfSize();
const Float_t * crystal = fGeometryEMCA->GetCrystalHalfSize() ;
const Float_t * pin = fGeometryEMCA->GetAPDHalfSize() ;
const Float_t * preamp = fGeometryEMCA->GetPreampHalfSize() ;
fCrystalShift=-inthermo[1]+strip[1]+splate[1]+crystal[1]-fGeometryEMCA->GetAirGapLed()/2.+pin[1]+preamp[1] ;
fCryCellShift=crystal[1]-(fGeometryEMCA->GetAirGapLed()-2*pin[1]-2*preamp[1])/2;
Int_t index ;
for ( index = 0; index < fNModules; index++ )
fPHOSAngle[index] = 0.0 ;
fRotMatrixArray = new TObjArray(fNModules) ;
SetPHOSAngles();
Double_t const kRADDEG = 180.0 / TMath::Pi() ;
Float_t r = GetIPtoOuterCoverDistance() + fPHOSParams[3] - GetCPVBoxSize(1) ;
for (Int_t iModule=0; iModule<fNModules; iModule++) {
fModuleCenter[iModule][0] = r * TMath::Sin(fPHOSAngle[iModule] / kRADDEG );
fModuleCenter[iModule][1] =-r * TMath::Cos(fPHOSAngle[iModule] / kRADDEG );
fModuleCenter[iModule][2] = 0.;
fModuleAngle[iModule][0][0] = 90;
fModuleAngle[iModule][0][1] = fPHOSAngle[iModule];
fModuleAngle[iModule][1][0] = 0;
fModuleAngle[iModule][1][1] = 0;
fModuleAngle[iModule][2][0] = 90;
fModuleAngle[iModule][2][1] = 270 + fPHOSAngle[iModule];
}
}
AliPHOSGeometry * AliPHOSGeometry::GetInstance()
{
return static_cast<AliPHOSGeometry *>( fgGeom ) ;
}
AliPHOSGeometry * AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t* title)
{
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 )
::Error("GetInstance", "Current geometry is %s. You cannot call %s",
fgGeom->GetName(), name) ;
else
rv = (AliPHOSGeometry *) fgGeom ;
}
return rv ;
}
void AliPHOSGeometry::SetPHOSAngles()
{
Double_t const kRADDEG = 180.0 / TMath::Pi() ;
Float_t pphi = 2 * TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ;
pphi *= kRADDEG ;
if (pphi > fAngle){
AliError(Form("PHOS modules overlap!\n pphi = %f fAngle = %f",
pphi, fAngle));
}
pphi = fAngle;
for( Int_t i = 1; i <= fNModules ; i++ ) {
Float_t angle = pphi * ( i - fNModules / 2.0 - 0.5 ) ;
fPHOSAngle[i-1] = - angle ;
}
}
void AliPHOSGeometry::GetGlobal(const AliRecPoint* , TVector3 & ) const
{
AliFatal(Form("Please use GetGlobalPHOS(recPoint,gpos) instead of GetGlobal!"));
}
void AliPHOSGeometry::GetGlobalPHOS(const AliPHOSRecPoint* recPoint, TVector3 & gpos) const
{
const AliPHOSRecPoint * tmpPHOS = recPoint ;
TVector3 localposition ;
tmpPHOS->GetLocalPosition(gpos) ;
if (!gGeoManager){
AliFatal("Geo manager not initialized\n");
}
char path[100] ;
Double_t dy ;
if(tmpPHOS->IsEmc()){
TString spath="/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1";
snprintf(path,spath.Length(),spath.Data(),tmpPHOS->GetPHOSMod()) ;
dy=fCrystalShift ;
}
else{
TString spath="/ALIC_1/PHOS_%d/PCPV_1";
snprintf(path,spath.Length(),spath.Data(),tmpPHOS->GetPHOSMod()) ;
dy= GetCPVBoxSize(1)/2. ;
}
Double_t pos[3]={gpos.X(),gpos.Y()-dy,gpos.Z()} ;
if(tmpPHOS->IsEmc())
pos[2]=-pos[2] ;
Double_t posC[3] = {};
if (!gGeoManager->cd(path)){
AliFatal("Geo manager can not find path \n");
}
TGeoHMatrix *m = gGeoManager->GetCurrentMatrix();
if (m){
m->LocalToMaster(pos,posC);
}
else{
AliFatal("Geo matrixes are not loaded \n") ;
}
gpos.SetXYZ(posC[0],posC[1],posC[2]) ;
}
void AliPHOSGeometry::GetModuleCenter(TVector3& center,
const char *det,
Int_t module) const
{
Float_t rDet = 0.;
if (strcmp(det,"CPV") == 0) rDet = GetIPtoCPVDistance ();
else if (strcmp(det,"EMC") == 0) rDet = GetIPtoCrystalSurface();
else
AliFatal(Form("Wrong detector name %s",det));
Float_t angle = GetPHOSAngle(module);
angle *= TMath::Pi()/180;
angle += 3*TMath::Pi()/2.;
center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.);
}