X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSGeometry.cxx;h=b78404bd9fae34d99018c02d95cbacae422c7f48;hb=da2ee9fc96ded6eef1542fbe2f1fe1fbb7cca5d7;hp=bfd942387e140fe445d9af0ee1e5f9a81fbb83fa;hpb=7b51037fc6b1e59c5845434bb31d585c00712919;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSGeometry.cxx b/PHOS/AliPHOSGeometry.cxx index bfd942387e1..b78404bd9fa 100644 --- a/PHOS/AliPHOSGeometry.cxx +++ b/PHOS/AliPHOSGeometry.cxx @@ -24,13 +24,14 @@ // can be easily implemented // The title is used to identify the version of CPV used. // -//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH) +// -- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC "KI" & SUBATECH) // --- ROOT system --- #include "TVector3.h" #include "TRotation.h" #include "TParticle.h" +#include // --- Standard library --- @@ -43,21 +44,57 @@ 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() { +AliPHOSGeometry::AliPHOSGeometry() : + fNModules(0), + fAngle(0.f), + fPHOSAngle(0), + fIPtoUpperCPVsurface(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 - fPHOSAngle = 0 ; - fGeometryEMCA = 0 ; - fGeometrySUPP = 0 ; - fGeometryCPV = 0 ; fgGeom = 0 ; - fRotMatrixArray = 0 ; } +//____________________________________________________________________________ +AliPHOSGeometry::AliPHOSGeometry(const AliPHOSGeometry & rhs) + : AliGeometry(rhs), + fNModules(rhs.fNModules), + fAngle(rhs.fAngle), + fPHOSAngle(0), + fIPtoUpperCPVsurface(rhs.fIPtoUpperCPVsurface), + 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), + fRotMatrixArray(0), + fGeometryEMCA(0), + fGeometryCPV(0), + fGeometrySUPP(0) +{ + // ctor only for internal usage (singleton) + Init() ; + fgGeom = this; +} + //____________________________________________________________________________ AliPHOSGeometry::~AliPHOSGeometry(void) { @@ -67,24 +104,24 @@ AliPHOSGeometry::~AliPHOSGeometry(void) if (fRotMatrixArray) delete fRotMatrixArray ; if (fPHOSAngle ) delete[] fPHOSAngle ; } -//____________________________________________________________________________ +//____________________________________________________________________________ void AliPHOSGeometry::Init(void) { // Initializes the PHOS parameters : // IHEP is the Protvino CPV (cathode pad chambers) TString test(GetName()) ; - if (test != "IHEP" ) { + if (test != "IHEP" && test != "noCPV") { AliFatal(Form("%s is not a known geometry (choose among IHEP)", test.Data() )) ; } fgInit = kTRUE ; - + fNModules = 5; fAngle = 20; - + fGeometryEMCA = new AliPHOSEMCAGeometry(); fGeometryCPV = new AliPHOSCPVGeometry (); @@ -103,14 +140,41 @@ void AliPHOSGeometry::Init(void) 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() - this->SetPHOSAngles() ; fRotMatrixArray = new TObjArray(fNModules) ; - + + // 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; iModuleGetStripHalfSize() ; - Float_t x0 = fGeometryEMCA->GetNStripX()*strip[0] ; - Float_t z0 = fGeometryEMCA->GetNStripZ()*strip[2] ; - Double_t angle = TMath::ATan( x0 / y0 ) ; - phi = phi + 1.5 * TMath::Pi() ; // to follow the convention of the particle generator(PHOS is between 220 and 320 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( z0 / 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) const +//____________________________________________________________________________ +void AliPHOSGeometry::GetGlobal(const AliRecPoint* recPoint, TVector3 & gpos) 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 { - AliWarning(Form("%s unknown option; result in radian", opt)) ; - conv = 1. ; - } - - Float_t y0 = GetIPtoCrystalSurface() ; - 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( - GetIPtoCrystalSurface()) ; - - } - else - { // it is a CPV - gpos.SetY(- GetIPtoUpperCPVsurface() ) ; - } - - Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; - Double_t const kRADDEG = 180.0 / TMath::Pi() ; - 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 coordinates of a RecPoint in the ALICE global coordinate system - - AliPHOSRecPoint * tmpPHOS = (AliPHOSRecPoint *) RecPoint ; - TVector3 localposition ; - tmpPHOS->GetLocalPosition(gpos) ; - - - if ( tmpPHOS->IsEmc() ) // it is a EMC crystal - { gpos.SetY( - GetIPtoCrystalSurface() ) ; - } - else - { // it is a CPV - gpos.SetY(- GetIPtoUpperCPVsurface() ) ; - } - - Float_t phi = GetPHOSAngle( tmpPHOS->GetPHOSMod()) ; - Double_t const kRADDEG = 180.0 / TMath::Pi() ; - 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(Double_t theta, Double_t phi, Int_t & moduleNumber, Double_t & z, Double_t & x) const -{ - // calculates the impact coordinates on PHOS of a neutral particle - // emitted in the direction theta and phi in the ALICE global coordinate system - - // 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 (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); } - if ( moduleNumber != 0 ) { - Float_t phi0 = GetPHOSAngle(moduleNumber) * (TMath::Pi() / 180.) + 1.5 * TMath::Pi() ; - Float_t y0 = GetIPtoCrystalSurface() ; - Double_t angle = phi - phi0; - x = y0 * TMath::Tan(angle) ; - angle = theta - TMath::Pi() / 2 ; - z = y0 * TMath::Tan(angle) ; + //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 ; } -} - -//____________________________________________________________________________ -void AliPHOSGeometry::ImpactOnEmc(const TVector3& vec, Int_t & moduleNumber, Double_t & z, Double_t & x) const -{ - // calculates the impact coordinates on PHOS of a neutral particle - // emitted in the direction theta and phi in the ALICE global coordinate system - // searches for the PHOS EMC module - - Double_t theta = vec.Theta() ; - Double_t phi = vec.Phi() ; + 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]) ; - ImpactOnEmc(theta, phi, moduleNumber, z, x) ; } - //____________________________________________________________________________ -void AliPHOSGeometry::ImpactOnEmc(const TParticle& p, Int_t & moduleNumber, Double_t & z, Double_t & x) 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 impact coordinates on PHOS of a neutral particle - // emitted in the direction theta and phi in the ALICE global coordinate system + // 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]) ; - // searches for the PHOS EMC module - Double_t theta = p.Theta() ; - Double_t phi = p.Phi() ; + 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"); + } + 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())Vx(),particle->Vy(),particle->Vz()} ; Double_t z,x; - ImpactOnEmc(particle->Theta(),particle->Phi(),moduleNumber,z,x); - if(moduleNumber) + ImpactOnEmc(vtx,particle->Theta(),particle->Phi(),moduleNumber,z,x); + if(moduleNumber!=0) in=kTRUE; - else - in=kFALSE; return in; } //____________________________________________________________________________ -Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId) const +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 + // 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 ) { // it is a Phos crystal - AbsId = + 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 { // it is a CPV pad - AbsId = GetNPhi() * GetNZ() * GetNModules() // the offset to separate EMCA crystals from CPV pads + 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 @@ -423,55 +414,128 @@ Bool_t AliPHOSGeometry::RelToAbsNumbering(const Int_t * relid, Int_t & AbsId) c } //____________________________________________________________________________ - void AliPHOSGeometry::RelPosInAlice(Int_t id, TVector3 & pos ) const { // Converts the absolute numbering into the global ALICE coordinate system + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } - 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 = - GetIPtoCrystalSurface() ; - else - y0 = - GetIPtoUpperCPVsurface() ; - - Float_t x, z ; - RelPosInModule(relid, x, z) ; - - pos.SetX(x) ; - pos.SetZ(z) ; - pos.SetY(y0) ; - - Float_t phi = GetPHOSAngle( phosmodule) ; - Double_t const kRADDEG = 180.0 / TMath::Pi() ; - Float_t rphi = phi / kRADDEG ; - - TRotation rot ; - rot.RotateZ(-rphi) ; // a rotation around Z by angle + Int_t relid[4] ; - TRotation dummy = rot.Invert() ; // to transform from original frame to rotate frame + AbsToRelNumbering(id , relid) ; - pos.Transform(rot) ; // rotate the baby + //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"); + } + 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]) ; + } + 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::RelPosToAbsId(Int_t module, Double_t x, Double_t z, Int_t & AbsId) const +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] ; - 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.) ) ; - - RelToAbsNumbering(relid,AbsId) ; + 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); + + } + } //____________________________________________________________________________ @@ -480,18 +544,98 @@ void AliPHOSGeometry::RelPosInModule(const Int_t * relid, Float_t & x, Float_t & // 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 row = relid[2] ; //offset along x axis - Int_t column = relid[3] ; //offset along z axis - - if ( relid[1] == 0 ) { // its a PbW04 crystal - x = - ( GetNPhi()/2. - row + 0.5 ) * GetCellStep() ; // position of Xtal with respect - z = - ( GetNZ() /2. - column + 0.5 ) * GetCellStep() ; // of center of PHOS module - } - else { - x = - ( GetNumberOfCPVPadsPhi()/2. - row - 0.5 ) * GetPadSizePhi() ; // position of pad with respect - z = - ( GetNumberOfCPVPadsZ() /2. - column - 0.5 ) * GetPadSizeZ() ; // of center of PHOS module + if (!gGeoManager){ + AliFatal("Geo manager not initialized\n"); + } + //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 ; } + 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"); + } + 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 ; + + } + } //____________________________________________________________________________ @@ -501,6 +645,7 @@ void AliPHOSGeometry::GetModuleCenter(TVector3& center, 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(); @@ -520,9 +665,56 @@ void AliPHOSGeometry::Global2Local(TVector3& localPosition, 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) ; }