LArBarrelPresamplerGeometry.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
/********************************************************************
 
NAME:     LArBarrelPresamplerGeometry.cxx
PACKAGE:  offline/LArCalorimeter/LArG4/LArG4Barrel
 
AUTHORS:  G. Unal, L. Carminati
CREATED:  September, 2004
 
PURPOSE:  'geometrical' methods used by the LArBarrelPresamplerCalculator.
          These original methods (previously in LArBarrelPresampler Calculator) were
          written by Bill Seligman
 
UPDATES:  - Calculate identifier method used by PresamplerCalibrationCalculator.
          (sept 2004)
          - Extended (added findCell method) to compute distance to electrode
          (GU oct 2005). Also improved computation of etaBin taking into account
          more correctly effect of tilted electrodes
 
********************************************************************/
 
//#undef DEBUGHITS
 
#include "LArBarrelPresamplerGeometry.h"
#include "LArG4Code/LArVG4DetectorParameters.h"
#include "LArG4Code/LArG4Identifier.h"
 
#include "G4AffineTransform.hh"
#include "G4NavigationHistory.hh"
 
#include "G4ThreeVector.hh"
#include "G4StepPoint.hh"
#include "G4Step.hh"
 
#include "G4ios.hh"
#include "globals.hh"
 
#include <cmath>
#include <limits.h>
#include <iostream>
 
namespace LArG4 {
 
  namespace BarrelPresampler {
 
    Geometry::Geometry(const std::string& name, ISvcLocator *pSvcLocator)
      : base_class(name, pSvcLocator)
    {
    }
 
    // ==================================================================
    StatusCode Geometry::initialize()
    {
      const double eps=0.007*Athena::Units::mm;
      m_end_module[0]=(m_mod[0][0]*m_cmm+2*eps)+m_zminPS+eps;
      for (int i=1;i<8;i++) m_end_module[i]=m_end_module[i-1]+(m_mod[i][0]*m_cmm+2*eps)+eps;
#ifdef DEBUGHITS
      for (int i=0;i<8;i++) {
        ATH_MSG_VERBOSE("Module length " << m_mod[0][0]*m_cmm+2*eps);
        ATH_MSG_VERBOSE("End of Module " << m_end_module[i]);
      }
#endif
 
      m_cat_th=m_cathode_th*m_cmm;
      m_first_cathod[0]=m_zminPS+m_mod[0][5]*m_cmm+m_cat_th/2.+2*eps;
      for (int i=1;i<8;i++) m_first_cathod[i]=m_end_module[i-1]+m_mod[i][5]*m_cmm+m_cat_th/2.+2*eps;
 
#ifdef DEBUGHITS
      for (int i=0;i<8;i++) ATH_MSG_VERBOSE("position of first cathode " << m_first_cathod[i]);
#endif
 
      // number of cells in eta per module
      for (int i=0;i<7;i++) m_ncell_module[i]=8;
      m_ncell_module[7]=5;
 
      // electrode tild in rad
      for (int i=0;i<8;i++) m_tilt[i]=m_mod[i][3]*CLHEP::deg;
 
      // number of gaps per cell    module 7 is somewhat pathological (last cell is shorter)
      for (int i=0;i<7;i++) m_ngap_cell[i]=(int)((m_mod[i][1]+0.1)/m_ncell_module[i]);
      m_ngap_cell[7]=18;
#ifdef DEBUGHITS
      for (int i=0;i<8;i++) ATH_MSG_VERBOSE("ngap per cell " << m_ngap_cell[i]);
#endif
 
      // pitch in z of gaps
      for (int i=0;i<8;i++) m_pitch[i]=m_mod[i][4]*m_cmm;
 
      return StatusCode::SUCCESS;
    }
 
    //======================================================================================
    //======================================================================================
 
    LArG4Identifier Geometry::CalculateIdentifier(const G4Step* a_step) const
    {
      const static G4String fullPSName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler" : G4String(m_detectorName+"::LAr::Barrel::Presampler");
      const static G4String fullCryoName = (m_detectorName.empty()) ? "LAr::TBBarrel::Cryostat::LAr" : G4String(m_detectorName+"::LAr::TBBarrel::Cryostat::LAr");
      const static G4String fullModuleName = (m_detectorName.empty()) ? "LAr::Barrel::Presampler::Module" : G4String(m_detectorName+"::LAr::Barrel::Presampler::Module");

      const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
      const G4int ndep = g4navigation->GetDepth();
      bool iactive(false);
      G4int idep(-999);

      for (G4int ii=0;ii<=ndep;ii++) {
        // FIXME Need to find a way to avoid these string-comparisons
        const G4String& vname = g4navigation->GetVolume(ii)->GetName();
        if (idep<0 && vname==fullPSName) idep=ii;    // half barrel
        else if (!iactive && vname==fullModuleName) iactive=true;
      }
 
      if (idep < 0) std::abort();
 
      if ( iactive ) {
        return this->CalculatePSActiveIdentifier( a_step , idep );
      }
      return this->CalculatePS_DMIdentifier( a_step , idep );
    }
 
    // ==========================================================================================
 
    LArG4Identifier Geometry::CalculatePSActiveIdentifier(const G4Step* a_step, const G4int ind) const
    {
      LArG4Identifier PSActiveID = LArG4Identifier();
 
      const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
 
#ifdef  DEBUGHITS
      ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
      ATH_MSG_VERBOSE("Eta and Phi in the atlas frame                  --> " << p.eta() << " " << p.phi());
#endif
 
      /* to get coordinates in the Presampler frame */
 
      const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
      const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
 
      const G4int zSide(this->determineZSide(p.z()));
 
      CalcData currentCellData;
      (void)this->findCell(currentCellData,ploc.x(),ploc.y(),ploc.z());
 
      if( zSide == -1 )
        {
          currentCellData.phiBin = 31 - currentCellData.phiBin;
          if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
        }

      PSActiveID   << 4          // LArCalorimeter
                   << 1          // LArEM
                   << zSide
                   << currentCellData.sampling
                   << currentCellData.region
                   << currentCellData.etaBin
                   << currentCellData.phiBin;
 
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE("Here the identifier for the presampler ACTIVE ----> ");
      ATH_MSG_VERBOSE("m_zSide  ----> " << zSide);
      ATH_MSG_VERBOSE("m_sampling  ----> " << currentCellData.sampling);
      ATH_MSG_VERBOSE("m_region  ----> " <<  currentCellData.region);
      ATH_MSG_VERBOSE("currentCellData.etaBin  ----> " << currentCellData.etaBin);
      ATH_MSG_VERBOSE("currentCellData.phiBin  ----> " << currentCellData.phiBin);
#endif
 
      return PSActiveID ;
    }
 
    //==========================================================================================
 
    LArG4Identifier Geometry::CalculatePS_DMIdentifier(const G4Step* a_step, const G4int ind) const
    {
 
      /******************************************************************************
  CaloDM_ID identifier:
 
  detector_system/subdet/type/sampling/region/eta/phi
  detector system = 10   -> Calorimeters
  subdet          = +/-4 -> LAr dead materials
  type            = 1    -> dead materials outside accordion and active presampler layers
  sampling        = 1    -> dead materials in front and in active LAr calorimeters
                              (starting from front warm cryostat walls)
  regions:        = 0 barrel warm wall and solenoid in front of the barrel presampler, 0 < |eta| < 1.5
                  = 1 barrel cryostat cold wall in front of the barrel presampler, 0 < |eta| < 1.5
                  = 2 all materials in front of the barrel presampler at radius larger than cold wall
                      outer radius, 0 < |eta| < 1.5
                  = 3 all materials from the active layer of the barrel presampler to the active layer
                      of accordion, 0 < |eta| < 1.5
 
   ---> Granularity : deta       0.1          granularity within region
                      dphi       pi/32 ~ 0.1  granularity within region
 
      ***********************************************************************************/
 
      const G4ThreeVector p = (a_step->GetPostStepPoint()->GetPosition() + a_step->GetPreStepPoint()->GetPosition()) * 0.5;
 
#ifdef  DEBUGHITS
      ATH_MSG_VERBOSE("Position of the step in the ATLAS frame (x,y,z) --> " << p.x() << " " << p.y() << " " << p.z());
      ATH_MSG_VERBOSE("Eta and Phi in the atlas frame                  --> " << p.eta() << " " << p.phi());
#endif

 
      const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
      const G4ThreeVector ploc = g4navigation->GetTransform(ind).TransformPoint(p);
      const G4double radius=sqrt(ploc.x()*ploc.x() + ploc.y()*ploc.y());

      const G4ThreeVector ploc2(ploc.x(),ploc.y(),ploc.z()+m_zpres+m_zminPS);
 
#ifdef  DEBUGHITS
      ATH_MSG_VERBOSE("Position of the step after traslation (x,y,z) --> " << ploc2.x() << " " << ploc2.y() << " " << ploc2.z());
      ATH_MSG_VERBOSE("Eta and Phi after translation                 --> " << ploc2.eta() << " " << ploc2.phi());
#endif
 
      // 01-Feb-2001 WGS: Add zSide calculation.
      const G4int zSide(this->determineZSide(p.z()));

      const G4double phi = (ploc2.phi() < 0.) ? ploc2.phi()+2.*M_PI : ploc2.phi();
      const G4double eta = ploc2.eta();
      //G4double z2=fabs(ploc2.z());

 
      const G4int numberPhiMod = 32;
      const G4double dphi = ( 2.*M_PI ) / numberPhiMod;
      const G4double inv_dphi = 1. / dphi;
      const G4double PSModuleRmean = 1420 ;
      const G4double phicheck = phi - int(phi * inv_dphi) * dphi  - (dphi /2.);
      const G4double Rcheck =  PSModuleRmean  / cos(phicheck);
      CalcData currentCellData;
      if (radius > Rcheck) {
        currentCellData.region = 3;
      } else {
        currentCellData.region = 2;
      }
 
      const G4double detaDM = 0.1 ;
      const G4double dphiDM = ( 2 * M_PI ) / 64. ;
 
      currentCellData.phiBin = G4int( phi * (1. / dphiDM) );
      currentCellData.etaBin = G4int( eta * (1. / detaDM) );
 
      if( zSide == -1 )
        {
          currentCellData.phiBin = 31 - currentCellData.phiBin;
          if(currentCellData.phiBin < 0 ) currentCellData.phiBin += 64;
        }
 
      // 07-Jul-2005 WGS: Handle an extremely rare rounding problem.
      if ( currentCellData.phiBin == 64 ) currentCellData.phiBin = 0;

      LArG4Identifier PS_DM_ID = LArG4Identifier();
      PS_DM_ID     << 10               // ATLAS
                   << zSide*4          // LArEM
                   << 1
                   << 1
                   << currentCellData.region
                   << currentCellData.etaBin
                   << currentCellData.phiBin;
 
#ifdef  DEBUGHITS
      ATH_MSG_VERBOSE("Here the identifier for the presampler DEAD materials ---->");
      ATH_MSG_VERBOSE("m_zSide  ----> " << zSide*4);
      ATH_MSG_VERBOSE("region   ----> " << currentCellData.region);
      ATH_MSG_VERBOSE("etaBin   ----> " << currentCellData.etaBin);
      ATH_MSG_VERBOSE("phiBin   ----> " << currentCellData.phiBin);
#endif
 
      return PS_DM_ID ;
    }

 
    bool  Geometry::findCell(CalcData & currentCellData, G4double xloc,G4double yloc,G4double zloc) const
    {
 
      currentCellData.sampling = 0;
      currentCellData.region   = 0;

      G4double phi = atan2(yloc,xloc);
      if ( phi < 0. ) phi += 2.*M_PI;
      const G4double z2=fabs(zloc+m_zpres+m_zminPS);

      const G4int numberPhiBins = 64;
      const G4double dphi = ( 2.*M_PI ) / numberPhiBins;
      const G4double inv_dphi = 1. / dphi;
      currentCellData.phiBin = G4int( phi * inv_dphi );
      if (currentCellData.phiBin >63) currentCellData.phiBin=63;
      if (currentCellData.phiBin <0)  currentCellData.phiBin=0;

      if (z2 < m_zminPS ) {
        currentCellData.etaBin=0;
        return false;
      }
      if (z2 > m_end_module[7]) {
        currentCellData.etaBin=60;
        return false;
      }

      currentCellData.module=0;
      for (int i=1;i<8;i++) {
        if (m_first_cathod[i]>=z2) break;
        currentCellData.module++;
      }
      if (currentCellData.module <0 || currentCellData.module > 7)
        {
          G4cerr << "LArBarrelPresampler: invalid module/hit " << currentCellData.module << " " << z2 << G4endl;
          if (currentCellData.module<0) currentCellData.etaBin=0;
          if (currentCellData.module>7) currentCellData.etaBin=60;
          return false;
        }

      const G4int isect=G4int(phi*m_nsectors/(2.*M_PI));
      const G4double phi0= ((double)isect+0.5)*2.*M_PI/((double) m_nsectors);
      static const G4double r0=1420.4*CLHEP::mm;   // FIXME should be recomputed from database
      currentCellData.dist=(xloc*cos(phi0)+yloc*sin(phi0)-r0);
 
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE("sector number, dist along height " << isect << " " << currentCellData.dist);
      ATH_MSG_VERBOSE("z2,module number,m_first_cathod " << z2 << " " << currentCellData.module << " "
                      << m_first_cathod[currentCellData.module]);
#endif

      G4double deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
      if (deltaz<0 ) {
        if (currentCellData.module>0) {
          currentCellData.module=currentCellData.module-1;
          deltaz=z2-(m_first_cathod[currentCellData.module]+currentCellData.dist*tan(m_tilt[currentCellData.module]));
        }
        else deltaz=0;
      }

      currentCellData.gap = ((int)(deltaz/m_pitch[currentCellData.module]));
 
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE("deltaz from first cathode,gap number " << deltaz << " " << currentCellData.gap);
#endif

      currentCellData.etaBin= currentCellData.gap/m_ngap_cell[currentCellData.module];
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE("etaBin inside module " << currentCellData.etaBin);
#endif
      if (currentCellData.etaBin >= m_ncell_module[currentCellData.module]) currentCellData.etaBin=m_ncell_module[currentCellData.module]-1;
 
      for (int i=0;i<currentCellData.module;i++) currentCellData.etaBin=currentCellData.etaBin+m_ncell_module[i];
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE(" final etaBin " << currentCellData.etaBin);
#endif
 
      if (currentCellData.etaBin < 0 || currentCellData.etaBin > 60) {
        ATH_MSG_WARNING("LArBarrelPresamplerGeometry::findCell  etaBin outside range " << currentCellData.etaBin);
      }

      const G4double zmiddle=m_first_cathod[currentCellData.module]+((double)(currentCellData.gap+0.5))*m_pitch[currentCellData.module];

      currentCellData.xElec=currentCellData.dist*cos(m_tilt[currentCellData.module])+(z2-zmiddle)*sin(m_tilt[currentCellData.module]);
      currentCellData.distElec=(z2-zmiddle)*cos(m_tilt[currentCellData.module]) - currentCellData.dist*sin(m_tilt[currentCellData.module]);
#ifdef DEBUGHITS
      ATH_MSG_VERBOSE("zmiddle,xloc,yloc " << zmiddle << " " << currentCellData.distElec << " " << currentCellData.xElec);
#endif
 
      bool status=true;
      if (fabs(currentCellData.dist)>m_halfThickLAr) {
#ifdef DEBUGHITS
        ATH_MSG_VERBOSE("Outside normal LAr 13mm gap "),
#endif
          status=false;
      }
 
      return status;
    }
 
  } 
 
}