X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSGeometry.cxx;h=ec31b1677a888ce508d85a4b8b30187c0ae3c551;hb=ed82d35ec5e88ef7df1d7a5baa882467571fe3c7;hp=ff4fb6f3eae35df23c5dbac56205102995bf8424;hpb=a4e98857bb4086fdb904ae5fbc807ae6065d5080;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSGeometry.cxx b/PHOS/AliPHOSGeometry.cxx index ff4fb6f3eae..ec31b1677a8 100644 --- a/PHOS/AliPHOSGeometry.cxx +++ b/PHOS/AliPHOSGeometry.cxx @@ -24,29 +24,83 @@ // can be easily implemented // The title is used to identify the version of CPV used. // -//*-- Author: Yves Schutz (SUBATECH) +// -- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH) // --- ROOT system --- #include "TVector3.h" #include "TRotation.h" +#include "TParticle.h" +#include +#include // --- Standard library --- -#include - // --- AliRoot header files --- - +#include "AliLog.h" #include "AliPHOSGeometry.h" #include "AliPHOSEMCAGeometry.h" -#include "AliPHOSPpsdRecPoint.h" -#include "AliConst.h" +#include "AliPHOSRecPoint.h" -ClassImp(AliPHOSGeometry) ; +ClassImp(AliPHOSGeometry) // these initialisations are needed for a singleton -AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ; -Bool_t AliPHOSGeometry::fgInit = kFALSE ; +AliPHOSGeometry * AliPHOSGeometry::fgGeom = 0 ; +Bool_t AliPHOSGeometry::fgInit = kFALSE ; + +//____________________________________________________________________________ +AliPHOSGeometry::AliPHOSGeometry() : + fNModules(0), + fAngle(0.f), + fPHOSAngle(0), + fIPtoUpperCPVsurface(0), + fCrystalShift(0), + fCryCellShift(0), + fRotMatrixArray(0), + fGeometryEMCA(0), + fGeometryCPV(0), + fGeometrySUPP(0) +{ + // default ctor + // must be kept public for root persistency purposes, but should never be called by the outside world + fgGeom = 0 ; +} + +//____________________________________________________________________________ +AliPHOSGeometry::AliPHOSGeometry(const AliPHOSGeometry & rhs) + : AliGeometry(rhs), + fNModules(rhs.fNModules), + fAngle(rhs.fAngle), + fPHOSAngle(0), + fIPtoUpperCPVsurface(rhs.fIPtoUpperCPVsurface), + fCrystalShift(rhs.fCrystalShift), + fCryCellShift(rhs.fCryCellShift), + fRotMatrixArray(0), + fGeometryEMCA(0), + fGeometryCPV(0), + fGeometrySUPP(0) +{ + Fatal("cpy ctor", "not implemented") ; +} + +//____________________________________________________________________________ +AliPHOSGeometry::AliPHOSGeometry(const Text_t* name, const Text_t* title) + : AliGeometry(name, title), + fNModules(0), + fAngle(0.f), + fPHOSAngle(0), + fIPtoUpperCPVsurface(0), + fCrystalShift(0), + fCryCellShift(0), + fRotMatrixArray(0), + fGeometryEMCA(0), + fGeometryCPV(0), + fGeometrySUPP(0) +{ + // ctor only for internal usage (singleton) + Init() ; + fgGeom = this; +} //____________________________________________________________________________ AliPHOSGeometry::~AliPHOSGeometry(void) @@ -55,80 +109,89 @@ AliPHOSGeometry::~AliPHOSGeometry(void) if (fRotMatrixArray) fRotMatrixArray->Delete() ; if (fRotMatrixArray) delete fRotMatrixArray ; - if (fPHOSAngle ) delete fPHOSAngle ; + if (fPHOSAngle ) delete[] fPHOSAngle ; } //____________________________________________________________________________ - void AliPHOSGeometry::Init(void) { // Initializes the PHOS parameters : // IHEP is the Protvino CPV (cathode pad chambers) - // GPS2 is the Subatech Pre-Shower (two micromegas sandwiching a passive lead converter) - // MIXT 4 PHOS modules withe the IHEP CPV qnd one PHOS module with the Subatche Pre-Shower - - if ( ((strcmp( fName, "GPS2" )) == 0) || - ((strcmp( fName, "IHEP" )) == 0) || - ((strcmp( fName, "MIXT" )) == 0) ) { - fgInit = kTRUE ; - - fNModules = 5; - fNPPSDModules = 0; - fAngle = 20; - - fGeometryEMCA = new AliPHOSEMCAGeometry(); - if ( ((strcmp( fName, "GPS2" )) == 0) ) { - fGeometryPPSD = new AliPHOSPPSDGeometry(); - fGeometryCPV = 0; - fNPPSDModules = fNModules; - } - else if ( ((strcmp( fName, "IHEP" )) == 0) ) { - fGeometryCPV = new AliPHOSCPVGeometry (); - fGeometryPPSD = 0; - fNPPSDModules = 0; - } - else if ( ((strcmp( fName, "MIXT" )) == 0) ) { - fGeometryCPV = new AliPHOSCPVGeometry (); - fGeometryPPSD = new AliPHOSPPSDGeometry(); - fNPPSDModules = 1; - } - fGeometrySUPP = new AliPHOSSupportGeometry(); + +/* + TString test(GetName()) ; + if (test != "IHEP" && test != "noCPV") { + AliFatal(Form("%s is not a known geometry (choose among IHEP)", + test.Data() )) ; + } +*/ + + fgInit = kTRUE ; + + fNModules = 5; + fAngle = 20; - fPHOSAngle = new Float_t[fNModules] ; - Int_t index ; - for ( index = 0; index < fNModules; index++ ) + fGeometryEMCA = new AliPHOSEMCAGeometry(); + + fGeometryCPV = new AliPHOSCPVGeometry (); + + fGeometrySUPP = new AliPHOSSupportGeometry(); + + fPHOSAngle = new Float_t[fNModules] ; + + 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) ; + + //calculate offset to crystal surface + Float_t * inthermo = fGeometryEMCA->GetInnerThermoHalfSize() ; + Float_t * strip = fGeometryEMCA->GetStripHalfSize() ; + Float_t* splate = fGeometryEMCA->GetSupportPlateHalfSize(); + Float_t * crystal = fGeometryEMCA->GetCrystalHalfSize() ; + Float_t * pin = fGeometryEMCA->GetAPDHalfSize() ; + 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 ; // Module position angles are set in CreateGeometry() + + fRotMatrixArray = new TObjArray(fNModules) ; - this->SetPHOSAngles() ; - fRotMatrixArray = new TObjArray(fNModules) ; - } - else { - fgInit = kFALSE ; - cout << "PHOS Geometry setup: option not defined " << fName << endl ; + // Geometry parameters are calculated + + SetPHOSAngles(); + Double_t const kRADDEG = 180.0 / TMath::Pi() ; + Float_t r = GetIPtoOuterCoverDistance() + fPHOSParams[3] - GetCPVBoxSize(1) ; + for (Int_t iModule=0; iModuleGetCPVBoxSize(index); - else if (strcmp(fName,"IHEP")==0) - return fGeometryCPV ->GetCPVBoxSize(index); - else if (strcmp(fName,"MIXT")==0) - return TMath::Max(fGeometryCPV ->GetCPVBoxSize(index), fGeometryPPSD->GetCPVBoxSize(index)); - else - return 0; -} +} //____________________________________________________________________________ AliPHOSGeometry * AliPHOSGeometry::GetInstance() { // Returns the pointer of the unique instance; singleton specific - return (AliPHOSGeometry *) fgGeom ; + return static_cast( fgGeom ) ; } //____________________________________________________________________________ @@ -153,10 +216,9 @@ AliPHOSGeometry * AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t } } else { - if ( strcmp(fgGeom->GetName(), name) != 0 ) { - cout << "AliPHOSGeometry : current geometry is " << fgGeom->GetName() << endl - << " you cannot call " << name << endl ; - } + if ( strcmp(fgGeom->GetName(), name) != 0 ) + ::Error("GetInstance", "Current geometry is %s. You cannot call %s", + fgGeom->GetName(), name) ; else rv = (AliPHOSGeometry *) fgGeom ; } @@ -167,11 +229,16 @@ AliPHOSGeometry * AliPHOSGeometry::GetInstance(const Text_t* name, const Text_t void AliPHOSGeometry::SetPHOSAngles() { // Calculates the position of the PHOS modules in ALICE global coordinate system + // in ideal geometry - Double_t const kRADDEG = 180.0 / kPI ; - Float_t pphi = 2 * TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoOuterCoverDistance() ) ) ; + Double_t const kRADDEG = 180.0 / TMath::Pi() ; + Float_t pphi = 2 * TMath::ATan( GetOuterBoxSize(0) / ( 2.0 * GetIPtoUpperCPVsurface() ) ) ; pphi *= kRADDEG ; - if (pphi > fAngle) cout << "AliPHOSGeometry: PHOS modules overlap!\n"; + 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++ ) { @@ -181,284 +248,174 @@ void AliPHOSGeometry::SetPHOSAngles() } //____________________________________________________________________________ -Bool_t AliPHOSGeometry::AbsToRelNumbering(const Int_t AbsId, Int_t * relid) +Bool_t AliPHOSGeometry::AbsToRelNumbering(Int_t absId, Int_t * relid) const { - // Converts the absolute numbering into the following array/ + // 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 + // = -1 if CPV + // relid[2] = Row number inside a PHOS module + // relid[3] = Column number inside a PHOS module Bool_t rv = kTRUE ; - Float_t id = AbsId ; + Float_t id = absId ; - Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / ( GetNPhi() * GetNZ() ) ) ; + Int_t phosmodulenumber = (Int_t)TMath:: Ceil( id / GetNCristalsInModule() ) ; - 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() ; - Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ; - relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ; - relid[1] = 1 ; - id -= ( relid[0] - 1 ) * nCPV ; - relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ; - relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; - } - else if ( strcmp(fName,"MIXT") == 0 ) { - id -= GetNPhi() * GetNZ() * GetNModules() ; - Float_t nPPSD = 2 * GetNumberOfModulesPhi() * GetNumberOfModulesZ() * GetNumberOfPadsPhi() * GetNumberOfPadsZ() ; - Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ; - if (id <= nCPV*GetNCPVModules()) { // this pad belons to CPV - relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ; - relid[1] = 1 ; - id -= ( relid[0] - 1 ) * nCPV ; - relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ; - relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; - } - else { // this pad belons to PPSD - id -= nCPV*GetNCPVModules(); - relid[0] = (Int_t)TMath::Ceil( id / nPPSD ); - id -= ( relid[0] - 1 ) * nPPSD ; - relid[0] += GetNCPVModules(); - 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() ) ; - } - } + if ( phosmodulenumber > GetNModules() ) { // it is a CPV pad + + id -= GetNPhi() * GetNZ() * GetNModules() ; + Float_t nCPV = GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() ; + relid[0] = (Int_t) TMath::Ceil( id / nCPV ) ; + relid[1] = -1 ; + id -= ( relid[0] - 1 ) * nCPV ; + relid[2] = (Int_t) TMath::Ceil( id / GetNumberOfCPVPadsZ() ) ; + relid[3] = (Int_t) ( id - ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() ) ; } - else { // its a PW04 crystal + else { // it is 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() ) ; + relid[2] = (Int_t)TMath::Ceil( id / GetNZ() ) ; + relid[3] = (Int_t)( id - ( relid[2] - 1 ) * GetNZ() ) ; } 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 one EMC (=PHOS) module - - Double_t conv ; - if ( opt == Radian() ) - conv = 1. ; - else if ( opt == Degre() ) - conv = 180. / TMath::Pi() ; - else { - cout << " 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) +//____________________________________________________________________________ +void AliPHOSGeometry::GetGlobal(const AliRecPoint* , TVector3 & ) const { - // calculates the angular coverage in theta and phi of a single crystal in a EMC(=PHOS) module - - Double_t conv ; - if ( opt == Radian() ) - conv = 1. ; - else if ( opt == Degre() ) - conv = 180. / TMath::Pi() ; - else { - cout << " 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 ; + AliFatal(Form("Please use GetGlobalPHOS(recPoint,gpos) instead of GetGlobal!")); } - //____________________________________________________________________________ -void AliPHOSGeometry::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos, TMatrix & gmat) const +void AliPHOSGeometry::GetGlobalPHOS(const AliPHOSRecPoint* recPoint, TVector3 & gpos) const { // Calculates the coordinates of a RecPoint and the error matrix in the ALICE global coordinate system - AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ; + const AliPHOSRecPoint * tmpPHOS = 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 + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } + //construct module name + char path[100] ; + Double_t dy ; + if(tmpPHOS->IsEmc()){ + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",tmpPHOS->GetPHOSMod()) ; +// sprintf(path,"/ALIC_1/PHOS_%d",tmpPHOS->GetPHOSMod()) ; + dy=fCrystalShift ; + } + else{ + sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",tmpPHOS->GetPHOSMod()) ; + dy= GetCPVBoxSize(1)/2. ; //center of CPV module + } + Double_t pos[3]={gpos.X(),gpos.Y()-dy,gpos.Z()} ; + if(tmpPHOS->IsEmc()) + pos[2]=-pos[2] ; //Opposite z directions in EMC matrix and local frame!!! + Double_t posC[3]; + //now apply possible shifts and rotations + 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::GetGlobal(const AliRecPoint* RecPoint, TVector3 & gpos) const +void AliPHOSGeometry::ImpactOnEmc(Double_t * vtx, Double_t theta, Double_t phi, + Int_t & moduleNumber, Double_t & z, Double_t & x) const { - // Calculates the coordinates of a RecPoint in the ALICE global coordinate system - - AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ; - TVector3 localposition ; - tmpPHOS->GetLocalPosition(gpos) ; - + // calculates the impact coordinates on PHOS of a neutral particle + // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system + TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ; + TVector3 v(vtx[0],vtx[1],vtx[2]) ; - if ( tmpPHOS->IsEmc() ) // it is a EMC crystal - { gpos.SetY( -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() + - GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ) ; + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } + + for(Int_t imod=1; imod<=GetNModules() ; imod++){ + //create vector from (0,0,0) to center of crystal surface of imod module + Double_t tmp[3]={0.,-fCrystalShift,0.} ; + + char path[100] ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",imod) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); } - 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 + TGeoHMatrix *m = gGeoManager->GetCurrentMatrix(); + Double_t posG[3]={0.,0.,0.} ; + if (m) m->LocalToMaster(tmp,posG); + TVector3 n(posG[0],posG[1],posG[2]) ; + Double_t direction=n.Dot(p) ; + if(direction<=0.) + continue ; //momentum directed FROM module + Double_t fr = (n.Mag2()-n.Dot(v))/direction ; + //Calculate direction in module plain + n-=v+fr*p ; + n*=-1. ; + Float_t * sz = fGeometryEMCA->GetInnerThermoHalfSize() ; //Wery close to the zise of the Xtl set + if(TMath::Abs(TMath::Abs(n.Z())= 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) ; - } + // Tells if a particle enters PHOS + Bool_t in=kFALSE; + Int_t moduleNumber=0; + Double_t vtx[3]={particle->Vx(),particle->Vy(),particle->Vz()} ; + Double_t z,x; + ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x); + if(moduleNumber!=0) + in=kTRUE; + return in; } //____________________________________________________________________________ -Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId) +Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & absId) const { // Converts the relative numbering into the absolute numbering // EMCA crystals: - // AbsId = from 1 to fNModules * fNPhi * fNZ - // PPSD gas cell: - // AbsId = from N(total EMCA crystals) + 1 - // to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ + - // fNModules * 2 * (fNumberOfModulesPhi * fNumberOfModulesZ) * fNumberOfPadsPhi * fNumberOfPadsZ + // absId = from 1 to fNModules * fNPhi * fNZ // CPV pad: - // AbsId = from N(total PHOS crystals) + 1 + // absId = from N(total PHOS crystals) + 1 // to NCPVModules * fNumberOfCPVPadsPhi * fNumberOfCPVPadsZ Bool_t rv = kTRUE ; - - if ( relid[1] > 0 && strcmp(fName,"GPS2")==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 PPSD 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 && strcmp(fName,"MIXT")==0) { // it is a PPSD pad - AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from PPSD pads - + GetNCPVModules() * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of CPV modules if any - + ( relid[0] - 1 - GetNCPVModules()) - * GetNumberOfModulesPhi() * GetNumberOfModulesZ() // the pads offset of PPSD 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 + + if ( relid[1] == 0 ) { // it is a Phos crystal + absId = + ( relid[0] - 1 ) * GetNPhi() * GetNZ() // the offset of PHOS modules + + ( relid[2] - 1 ) * GetNZ() // the offset along phi + + relid[3] ; // the offset along z } - - 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 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of PHOS modules - + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() // the pads offset of a CPV row + else { // it is a CPV pad + absId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from CPV pads + + ( relid[0] - 1 ) * GetNumberOfCPVPadsPhi() * GetNumberOfCPVPadsZ() // the pads offset of PHOS modules + + ( relid[2] - 1 ) * GetNumberOfCPVPadsZ() // the pads offset of a CPV row + relid[3] ; // the column number } @@ -466,114 +423,307 @@ Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId) } //____________________________________________________________________________ - -void AliPHOSGeometry::RelPosInAlice(const Int_t id, TVector3 & pos ) +void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const { // Converts the absolute numbering into the global ALICE coordinate system - // It works only for the GPS2 geometry - if (id > 0 && strcmp(fName,"GPS2")==0) { - - Int_t relid[4] ; - - AbsToRelNumbering(id , relid) ; + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } - Int_t phosmodule = relid[0] ; + Int_t relid[4] ; - Float_t y0 = 0 ; + AbsToRelNumbering(id , relid) ; - if ( relid[1] == 0 ) { // it is a PbW04 crystal - y0 = -(GetIPtoOuterCoverDistance() + GetUpperPlateThickness() - + GetSecondUpperPlateThickness() + GetUpperCoolingPlateThickness()) ; + //construct module name + char path[100] ; + if(relid[1]==0){ //this is EMC + + Double_t ps[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal + Double_t psC[3]={0.0,0.0,0.}; //Global position + + //Shift and possibly apply misalignment corrections + Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ; + Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ; + Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()- + (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ; + Int_t cellraw= relid[3]%nCellsZInStrip ; + if(cellraw==0)cellraw=nCellsZInStrip ; + Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d", + relid[0],strip,cell) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); } - 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) ; + TGeoHMatrix *m = gGeoManager->GetCurrentMatrix(); + if (m) m->LocalToMaster(ps,psC); + else{ + AliFatal("Geo matrixes are not loaded \n") ; } - - 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 + pos.SetXYZ(psC[0],psC[1],psC[2]) ; } - else { - pos.SetX(0.); - pos.SetY(0.); - pos.SetZ(0.); + else{ + //first calculate position with respect to CPV plain + Int_t row = relid[2] ; //offset along x axis + Int_t column = relid[3] ; //offset along z axis + Double_t ps[3]= {0.0,GetCPVBoxSize(1)/2.,0.}; //Position on top of CPV + Double_t psC[3]={0.0,0.0,0.}; //Global position + pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row - 0.5 ) * GetPadSizePhi() ; // position of pad with respect + pos[2] = - ( GetNumberOfCPVPadsZ() /2. - column - 0.5 ) * GetPadSizeZ() ; // of center of PHOS module + + //now apply possible shifts and rotations + sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); + } + TGeoHMatrix *m = gGeoManager->GetCurrentMatrix(); + if (m) m->LocalToMaster(ps,psC); + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + pos.SetXYZ(psC[0],psC[1],-psC[2]) ; } } //____________________________________________________________________________ -void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) +void AliPHOSGeometry::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & absId) const +{ + // converts local PHOS-module (x, z) coordinates to absId + + //find Global position + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } + Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!! + Double_t posG[3] ; + char path[100] ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); + } + TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix(); + if (mPHOS){ + mPHOS->LocalToMaster(posL,posG); + } + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + + Int_t relid[4] ; + gGeoManager->FindNode(posG[0],posG[1],posG[2]) ; + //Check that path contains PSTR and extract strip number + TString cpath(gGeoManager->GetPath()) ; + Int_t indx = cpath.Index("PCEL") ; + if(indx==-1){ //for the few events when particle hits between srips use ideal geometry + relid[0] = module ; + relid[1] = 0 ; + relid[2] = static_cast(TMath::Ceil( x/ GetCellStep() + GetNPhi() / 2.) ); + relid[3] = static_cast(TMath::Ceil(-z/ GetCellStep() + GetNZ() / 2.) ) ; + if(relid[2]<1)relid[2]=1 ; + if(relid[3]<1)relid[3]=1 ; + if(relid[2]>GetNPhi())relid[2]=GetNPhi() ; + if(relid[3]>GetNZ())relid[3]=GetNZ() ; + RelToAbsNumbering(relid,absId) ; + } + else{ + Int_t indx2 = cpath.Index("/",indx) ; + if(indx2==-1) + indx2=cpath.Length() ; + TString cell=cpath(indx+5,indx2-indx-5) ; + Int_t icell=cell.Atoi() ; + indx = cpath.Index("PSTR") ; + indx2 = cpath.Index("/",indx) ; + TString strip=cpath(indx+5,indx2-indx-5) ; + Int_t iStrip = strip.Atoi() ; + + Int_t row = fGeometryEMCA->GetNStripZ() - (iStrip - 1) % (fGeometryEMCA->GetNStripZ()) ; + Int_t col = (Int_t) TMath::Ceil((Double_t) iStrip/(fGeometryEMCA->GetNStripZ())) -1 ; + + // Absid for 8x2-strips. Looks nice :) + absId = (module-1)*GetNCristalsInModule() + + row * 2 + (col*fGeometryEMCA->GetNCellsXInStrip() + (icell - 1) / 2)*GetNZ() - (icell & 1 ? 1 : 0); + + } + +} + +//____________________________________________________________________________ +void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & z) const { // 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) - - Bool_t padOfCPV = (strcmp(fName,"IHEP")==0) || - ((strcmp(fName,"MIXT")==0) && relid[0]<=GetNCPVModules()) ; - Bool_t padOfPPSD = (strcmp(fName,"GPS2")==0) || - ((strcmp(fName,"MIXT")==0) && relid[0]> GetNCPVModules()) ; - 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 = 0; - Float_t padsizeX = 0; - Int_t nOfPadsPhi = 0; - Int_t nOfPadsZ = 0; - if ( padOfPPSD ) { - padsizeZ = GetPPSDModuleSize(2) / GetNumberOfPadsZ(); - padsizeX = GetPPSDModuleSize(0) / GetNumberOfPadsPhi(); - nOfPadsPhi = GetNumberOfPadsPhi(); - nOfPadsZ = GetNumberOfPadsZ(); + + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); } - else if ( padOfCPV ) { - padsizeZ = GetPadSizeZ(); - padsizeX = GetPadSizePhi(); - nOfPadsPhi = GetNumberOfCPVPadsPhi(); - nOfPadsZ = GetNumberOfCPVPadsZ(); + //construct module name + char path[100] ; + if(relid[1]==0){ //this is PHOS + +// Calculations using ideal geometry (obsolete) +// x = - ( GetNPhi()/2. - relid[2] + 0.5 ) * GetCellStep() ; // position of Xtal with respect +// z = - ( GetNZ() /2. - relid[3] + 0.5 ) * GetCellStep() ; // of center of PHOS module + + Double_t pos[3]= {0.0,-fCryCellShift,0.}; //Position incide the crystal + Double_t posC[3]={0.0,0.0,0.}; //Global position + + //Shift and possibly apply misalignment corrections + Int_t nCellsXInStrip=fGeometryEMCA->GetNCellsXInStrip() ; + Int_t nCellsZInStrip=fGeometryEMCA->GetNCellsZInStrip() ; +// Int_t strip=1+(relid[3]-1)/fGeometryEMCA->GetNCellsZInStrip()+((relid[2]-1)/nCellsInStrip)*fGeometryEMCA->GetNStripZ() ; + Int_t strip=1+((Int_t) TMath::Ceil((Double_t)relid[2]/nCellsXInStrip))*fGeometryEMCA->GetNStripZ()- + (Int_t) TMath::Ceil((Double_t)relid[3]/nCellsZInStrip) ; + Int_t cellraw= relid[3]%nCellsZInStrip ; + if(cellraw==0)cellraw=nCellsZInStrip ; + Int_t cell= ((relid[2]-1)%nCellsXInStrip)*nCellsZInStrip + cellraw ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1/PSTR_%d/PCEL_%d", + relid[0],strip,cell) ; + 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") ; + } + // printf("Local: x=%f, y=%f, z=%f \n",pos[0],pos[1],pos[2]) ; + // printf(" gl: x=%f, y=%f, z=%f \n",posC[0],posC[1],posC[2]) ; + //Return to PHOS local system + Double_t posL[3]={posC[0],posC[1],posC[2]}; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",relid[0]) ; + // sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); + } + TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix(); + if (mPHOS) mPHOS->MasterToLocal(posC,posL); + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } +//printf("RelPosInMod: posL=[%f,%f,%f]\n",posL[0],posL[1],posL[2]) ; +//printf("old: x=%f, z=%f \n",x,z); + x=posL[0] ; + z=-posL[2]; + return ; } - - 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 ( padOfPPSD ) { - 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() ; - x0 = ( GetNumberOfModulesPhi() / 2. - modrow + 0.5 ) * GetPPSDModuleSize(0) ; - z0 = ( GetNumberOfModulesZ() / 2. - modcol + 0.5 ) * GetPPSDModuleSize(2) ; - } else { - x0 = 0; - z0 = 0; + else{//CPV + //first calculate position with respect to CPV plain + Int_t row = relid[2] ; //offset along x axis + Int_t column = relid[3] ; //offset along z axis + Double_t pos[3]= {0.0,0.0,0.}; //Position incide the CPV printed circuit + Double_t posC[3]={0.0,0.0,0.}; //Global position + // x = - ( GetNumberOfCPVPadsPhi()/2. - row - 0.5 ) * GetPadSizePhi() ; // position of pad with respect + // z = - ( GetNumberOfCPVPadsZ() /2. - column - 0.5 ) * GetPadSizeZ() ; // of center of PHOS module + pos[0] = - ( GetNumberOfCPVPadsPhi()/2. - row - 0.5 ) * GetPadSizePhi() ; // position of pad with respect + pos[2] = - ( GetNumberOfCPVPadsZ() /2. - column - 0.5 ) * GetPadSizeZ() ; // of center of PHOS module + + //now apply possible shifts and rotations + sprintf(path,"/ALIC_1/PHOS_%d/PCPV_1",relid[0]) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); } - x = - ( nOfPadsPhi/2. - row - 0.5 ) * padsizeX + x0 ; // position of pad with respect - z = ( nOfPadsZ /2. - column - 0.5 ) * padsizeZ - z0 ; // of center of PHOS module + TGeoHMatrix *m = gGeoManager->GetCurrentMatrix(); + if (m) m->LocalToMaster(pos,posC); + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + //Return to PHOS local system + Double_t posL[3]={0.,0.,0.,} ; + sprintf(path,"/ALIC_1/PHOS_%d",relid[0]) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); + } + TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix(); + if (mPHOS) mPHOS->MasterToLocal(posC,posL); + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + x=posL[0] ; + z=posL[1]; + return ; + + } + +} + +//____________________________________________________________________________ + +void AliPHOSGeometry::GetModuleCenter(TVector3& center, + const char *det, + Int_t module) const +{ + // Returns a position of the center of the CPV or EMC module + // in ideal (not misaligned) geometry + 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); // (40,20,0,-20,-40) degrees + angle *= TMath::Pi()/180; + angle += 3*TMath::Pi()/2.; + center.SetXYZ(rDet*TMath::Cos(angle), rDet*TMath::Sin(angle), 0.); +} + +//____________________________________________________________________________ + +void AliPHOSGeometry::Global2Local(TVector3& localPosition, + const TVector3& globalPosition, + Int_t module) const +{ + // Transforms a global position of the rec.point to the local coordinate system + //Return to PHOS local system + Double_t posG[3]={globalPosition.X(),globalPosition.Y(),globalPosition.Z()} ; + Double_t posL[3]={0.,0.,0.} ; + char path[100] ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",module) ; +// sprintf(path,"/ALIC_1/PHOS_%d",module) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); } + TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix(); + if (mPHOS) mPHOS->MasterToLocal(posG,posL); + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + localPosition.SetXYZ(posL[0],posL[1]+fCrystalShift,-posL[2]) ; + +/* + Float_t angle = GetPHOSAngle(module); // (40,20,0,-20,-40) degrees + angle *= TMath::Pi()/180; + angle += 3*TMath::Pi()/2.; + localPosition = globalPosition; + localPosition.RotateZ(-angle); +*/ +} +//____________________________________________________________________________ +void AliPHOSGeometry::Local2Global(Int_t mod, Float_t x, Float_t z, + TVector3& globalPosition) const +{ + char path[100] ; + sprintf(path,"/ALIC_1/PHOS_%d/PEMC_1/PCOL_1/PTIO_1/PCOR_1/PAGA_1/PTII_1",mod) ; +// sprintf(path,"/ALIC_1/PHOS_%d",mod) ; + if (!gGeoManager->cd(path)){ + AliFatal("Geo manager can not find path \n"); + } + Double_t posL[3]={x,-fCrystalShift,-z} ; //Only for EMC!!! + Double_t posG[3] ; + TGeoHMatrix *mPHOS = gGeoManager->GetCurrentMatrix(); + if (mPHOS){ + mPHOS->LocalToMaster(posL,posG); + } + else{ + AliFatal("Geo matrixes are not loaded \n") ; + } + globalPosition.SetXYZ(posG[0],posG[1],posG[2]) ; +} +//____________________________________________________________________________ +void AliPHOSGeometry::GetIncidentVector(const TVector3 &vtx, Int_t module, Float_t x,Float_t z, TVector3 &vInc) const { + //Calculates vector pointing from vertex to current poisition in module local frame + //Note that PHOS local system and ALICE global have opposite z directions + + Global2Local(vInc,vtx,module) ; + vInc.SetXYZ(vInc.X()+x,vInc.Y(),vInc.Z()+z) ; }