LArBarrelCalculator Class Reference

#include <LArBarrelCalculator.h>

Inheritance diagram for LArBarrelCalculator:

Collaboration diagram for LArBarrelCalculator:

Public Member Functions
	LArBarrelCalculator (const std::string &name, ISvcLocator *pSvcLocator)
virtual StatusCode	initialize () override final
virtual StatusCode	finalize () override final
	LArBarrelCalculator (const LArBarrelCalculator &)=delete
LArBarrelCalculator &	operator= (const LArBarrelCalculator &)=delete
void	initializeForSDCreation () override final
virtual G4bool	Process (const G4Step *a_step, std::vector< LArHitData > &hdata) const override final
virtual G4bool	isInTime (G4double hitTime) const override final
virtual G4float	OOTcut () const override final

Protected Attributes
Gaudi::Property< bool >	m_BirksLaw {this, "BirksLaw", true}
Gaudi::Property< double >	m_Birksk {this, "Birksk", 0.05832}
Gaudi::Property< double >	m_OOTcut {this, "OOTcut", 300*CLHEP::ns}

Private Member Functions
G4bool	FiducialCuts (G4double, G4double, G4double) const
void	InitHV ()
double	ScaleHV (double, double, double, double) const

Private Attributes
ServiceHandle< ILArBarrelGeometry >	m_geometry {this, "GeometryCalculator", "LArBarrelGeometry"}
const AccMap *	m_accmap {nullptr}
std::unique_ptr< MapEta >	m_etamap1 {}
std::unique_ptr< MapEta >	m_etamap2 {}
std::unique_ptr< MapEta >	m_etamap3 {}
Gaudi::Property< bool >	m_IflCur {this, "EMBCurr", true}
Gaudi::Property< bool >	m_IflMapTrans {this, "EMBEtaTrans", true}
Gaudi::Property< bool >	m_IflXtalk {this, "EMBXtalk", true}
Gaudi::Property< double >	m_dstep {this, "EMBdstep", .2*CLHEP::mm}
const LArG4BirksLaw *	m_birksLaw {nullptr}
Gaudi::Property< bool >	m_doHV {this, "EMBHVEnable", false}
double	m_etaMaxBarrel {0.}
double	m_zMinBarrel {0.}
double	m_zMaxBarrel {0.}
double	m_rMinAccordion {0.}
double	m_rMaxAccordion {0.}
double	m_ThickAbs {0.}
double	m_ThickEle {0.}
int	m_NCellTot {0}
int	m_NCellMax {0}
bool	m_testbeam {false}
double	m_hv [2][1024][7][2] {}

Detailed Description

Definition at line 32 of file LArBarrelCalculator.h.

Constructor & Destructor Documentation

LArBarrelCalculator() [1/2]

LArBarrelCalculator::LArBarrelCalculator	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 47 of file LArBarrelCalculator.cxx.

  : LArCalculatorSvcImp(name, pSvcLocator)
{
}

LArBarrelCalculator() [2/2]

LArBarrelCalculator::LArBarrelCalculator ( const LArBarrelCalculator &  )

delete

Member Function Documentation

FiducialCuts()

G4bool LArBarrelCalculator::FiducialCuts ( G4double  radloc, G4double  zloc, G4double  etaloc  ) const

private

Definition at line 616 of file LArBarrelCalculator.cxx.

{
  if (radloc < m_rMinAccordion ||
      radloc > m_rMaxAccordion ||
      etaloc > m_etaMaxBarrel ||
      zloc   < m_zMinBarrel ||
      zloc   > m_zMaxBarrel)      return false;
  else                           return true;
}

finalize()

StatusCode LArBarrelCalculator::finalize ( )

finaloverridevirtual

Definition at line 146 of file LArBarrelCalculator.cxx.

{
  if (m_BirksLaw) delete m_birksLaw;
  return StatusCode::SUCCESS;
}

InitHV()

void LArBarrelCalculator::InitHV ( )

private

Definition at line 630 of file LArBarrelCalculator.cxx.

{
  ATH_MSG_DEBUG(" **** in LArBarrelCalculator::InitHV() ");
 
 
 
 
  float defaultHvVal=2000.;
  ATH_MSG_DEBUG(" defaultHvVal " << defaultHvVal);
  for (int ipm=0;ipm<2;ipm++) {
    for (int ielec=0;ielec<1024;ielec++) {
      for (int ieta=0;ieta<7;ieta++) {
        for (int iside=0;iside<2;iside++) {
          m_hv[ipm][ielec][ieta][iside] = defaultHvVal;
        }
      }
    }
  }
 
  if (m_doHV) {
    // initialize services
    SmartIF<StoreGateSvc> pDetStore{Gaudi::svcLocator()->service("DetectorStore")};
    if (!pDetStore) {
      ATH_MSG_WARNING("LArBarrelCalculator::InitHV() unable to get Detector Store! Use default HV values.");
      return;
    }
 
    // get EMBHV Manager
    const LArHVManager *manager = nullptr;
    if (pDetStore->retrieve(manager)==StatusCode::SUCCESS) {
      const EMBHVManager& hvManager=manager->getEMBHVManager();
      const EMBHVManager::EMBHVData hvdata = hvManager.getDataSim();
      ATH_MSG_DEBUG(" got HV Manager ");
      // loop over HV modules
      for (unsigned int iSide=0;iSide<2;iSide++) {
        for (unsigned int iPhi=0;iPhi<16;iPhi++) {
          for (unsigned int iSector=0;iSector<2;iSector++) {
            for (unsigned int iEta=0;iEta<7;iEta++) {
              const EMBHVModule& hvMod = hvManager.getHVModule(iSide,iEta+1,iPhi,iSector);
              for (unsigned int ielec=0;ielec<32;ielec++) {
                const EMBHVElectrode& electrode = hvMod.getElectrode(ielec);
                unsigned jElec = ielec+32*iSector+64*iPhi;
                for (unsigned int iGap=0;iGap<2;iGap++) {
                  double hv = hvdata.voltage (electrode, iGap);
                  ATH_MSG_DEBUG(" iSide,jElec,iEta,iGap,hv " << iSide << " " << jElec << " " << iEta << " " << iGap << " " << hv);
                  if (hv>-999.) m_hv[iSide][jElec][iEta][iGap] = hv;
                }
              }
            }
          }
        }
      }
    }
    else {
      ATH_MSG_WARNING(" Unable to find HV Manager ");
    }
  }
 
  return;
}

initialize()

StatusCode LArBarrelCalculator::initialize ( )

finaloverridevirtual

Definition at line 52 of file LArBarrelCalculator.cxx.

{
  ATH_MSG_DEBUG("LArBarrelCalculator: Beginning initialization ");
  if (m_BirksLaw) {
    const double Birks_LAr_density = 1.396;
    ATH_MSG_INFO("Use Birks_LAr_density="<<Birks_LAr_density<<", Birksk="<<(double)m_Birksk);
    m_birksLaw = new LArG4BirksLaw(Birks_LAr_density,m_Birksk);
  }
 
  // Access source of detector parameters.
  const LArVG4DetectorParameters* parameters = LArVG4DetectorParameters::GetInstance();
 
  // Get the out-of-time cut from the detector parameters routine.
  m_OOTcut = parameters->GetValue("LArExpHallOutOfTimeCut"); //FIXME should be done via configurables
  ATH_MSG_DEBUG("**** OutOfTime cut " << m_OOTcut);
 
  // Main Barrel parameters
  // All the UNITS are implicitly the GEANT4 ONES e.g. mm, rad, etc ...
  m_etaMaxBarrel = parameters->GetValue("LArEMBMaxEtaAcceptance");
  m_zMinBarrel = parameters->GetValue("LArEMBfiducialMothZmin");
  m_zMaxBarrel = parameters->GetValue("LArEMBfiducialMothZmax");
  m_NCellMax   = (int) (parameters->GetValue("LArEMBnoOFPhysPhiCell"));
  m_rMinAccordion  =  parameters->GetValue("LArEMBRadiusInnerAccordion");
  m_rMaxAccordion  = parameters->GetValue("LArEMBFiducialRmax");
 
  // absorbers and electrodes thickness
  m_ThickAbs = 0.5*( parameters->GetValue("LArEMBThickAbsGlue")
                     +parameters->GetValue("LArEMBThickAbsIron")
                     +parameters->GetValue("LArEMBThickAbsLead"));
 
  G4double check = 0.5*( parameters->GetValue("LArEMBThinAbsGlue")
                         +parameters->GetValue("LArEMBThinAbsIron")
                         +parameters->GetValue("LArEMBThinAbsLead"));
 
  if (fabs(check-m_ThickAbs)>0.001)
    {
      ATH_MSG_WARNING("Thin and Thick Abs have difference thickness!");
    }
 
  m_ThickEle= 0.5*( parameters->GetValue("LArEMBThickElecCopper")
                    +parameters->GetValue("LArEMBThickElecKapton"));
 
  // === GU 11/06/2003   total number of cells in phi
  // to distinguish 1 module (testbeam case) from full Atlas
  m_NCellTot = (int) (parameters->GetValue("LArEMBnoOFPhysPhiCell"));
  if (m_NCellTot != 1024)
    {
      m_NCellMax=1024;
      m_testbeam=true;
    }
  // ===
  //  access geometry computation class
  ATH_CHECK(m_geometry.retrieve());
 
  // access current maps if required
  if (m_IflCur)
    {
      ATH_MSG_DEBUG(" LArBarrelCalculator: start reading of current maps");
      m_accmap = AccMap::GetAccMap();
      m_etamap1 = std::make_unique<MapEta>(1);
      m_etamap2 = std::make_unique<MapEta>(2);
      if (m_IflMapTrans) m_etamap3 = std::make_unique<MapEta>(3);
      ATH_MSG_DEBUG(" LArBarrelCalculator: end of reading current maps");
    }
  // Initialize HV values
  this->InitHV();
 
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_NCellMax       " << m_NCellMax);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_NCellTot       " << m_NCellTot);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_rMinAccordion  " << m_rMinAccordion);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_rMaxAccordion  " << m_rMaxAccordion);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_zMinBarrel     " << m_zMinBarrel);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_zMaxBarrel     " << m_zMaxBarrel);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_etaMaxBarrel   " << m_etaMaxBarrel);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_ThickAbs       " << m_ThickAbs);
  ATH_MSG_DEBUG(" LArBarrelCalculator: s_ThickEle       " << m_ThickEle);
  if(m_IflCur) ATH_MSG_DEBUG(" LArBarrelCalculator: Deposited Energy  dE/dX  Corrected ==> CURRENT Option ON");
  else ATH_MSG_DEBUG(" LArBarrelCalculator: Crude Deposited Energy  dE/dX  NO CURRENT option");
  if (m_IflCur && m_IflMapTrans) ATH_MSG_DEBUG(" LArBarrelCalculator: Compute effect of E field around eta=0.8 ");
  else ATH_MSG_DEBUG(" LArBarrelCalculator: Ignore effect of E field around eta=0.8 ");
  if(m_BirksLaw)
    {
      ATH_MSG_DEBUG(" LArBarrelCalculator: Birks' law ON ");
      ATH_MSG_DEBUG(" LArBarrelCalculator:   parameter k    " << m_birksLaw->k());
    }
  else
    {
      ATH_MSG_DEBUG(" LArBarrelCalculator: Birks' law OFF");
    }
 
  return StatusCode::SUCCESS;
}

initializeForSDCreation()

void LArBarrelCalculator::initializeForSDCreation ( )

inlinefinaloverridevirtual

Reimplemented from LArCalculatorSvcImp.

Definition at line 43 of file LArBarrelCalculator.h.

{ m_geometry->initializeForSDCreation(); };

isInTime()

virtual G4bool LArBarrelCalculator::isInTime ( G4double  hitTime ) const

inlinefinaloverridevirtual

Definition at line 48 of file LArBarrelCalculator.h.

  {
    return !(std::fabs(hitTime) > m_OOTcut);
  }

OOTcut()

virtual G4float LArBarrelCalculator::OOTcut ( ) const

inlinefinaloverridevirtual

Definition at line 56 of file LArBarrelCalculator.h.

{ return m_OOTcut; }

operator=()

LArBarrelCalculator& LArBarrelCalculator::operator= ( const LArBarrelCalculator &  )

delete

Process()

G4bool LArBarrelCalculator::Process ( const G4Step *  a_step, std::vector< LArHitData > &  hdata  ) const

finaloverridevirtual

Definition at line 153 of file LArBarrelCalculator.cxx.

{
 
  hdata.clear();
 
  LArG4Identifier identifier2;
  LArG4Identifier identifier_xt1;
  LArG4Identifier identifier_xt2;
 
  //  check the Step content is non trivial
  G4double thisStepEnergyDeposit = step->GetTotalEnergyDeposit() * step->GetTrack()->GetWeight();
  G4double thisStepLength = step->GetStepLength() / Units::mm;
 
#ifdef  DEBUGSTEP
  ATH_MSG_DEBUG("******  LArBarrelCalculator:  Step energy,length "
                << thisStepEnergyDeposit << " " << thisStepLength);
#endif
 
  if(thisStepEnergyDeposit <= 0. || thisStepLength <= 0.)
    {
#ifdef DEBUGSTEP2
      ATH_MSG_DEBUG("   Invalid hit  trivial content");
#endif
      return false;
    }
 
  // Get Step geometrical parameters (first and end)
  G4StepPoint *thisStepPoint = step->GetPreStepPoint();
  G4StepPoint *thisStepBackPoint = step->GetPostStepPoint();
  G4ThreeVector startPoint = thisStepPoint->GetPosition();
  G4ThreeVector endPoint = thisStepBackPoint->GetPosition();
 
#ifdef  DEBUGSTEP
  ATH_MSG_DEBUG("   Beginning step position "
                << startPoint.x() << " " << startPoint.y() << " " << startPoint.z());
#endif
 
  // find zside from volume name
  G4int zSide = 1;
  if (!m_testbeam) {
    const G4NavigationHistory* g4navigation = thisStepPoint->GetTouchable()->GetHistory();
    G4int ndep = g4navigation->GetDepth();
    for (G4int ii=0;ii<=ndep;ii++) {
      G4VPhysicalVolume* v1 = g4navigation->GetVolume(ii);
      G4String vname = v1->GetName();
      if ( vname.find("NegPhysical") != std::string::npos) zSide=-1;
    }
  }
 
 
  // BACK directly into the LOCAL half_Barrel Z > 0. PART (mother "stac_phys1")
 
  const G4AffineTransform transformation =
    thisStepPoint->GetTouchable()->GetHistory()->GetTopTransform();
 
  G4ThreeVector startPointinLocal =
    transformation.TransformPoint(startPoint);
  G4ThreeVector   endPointinLocal =
    transformation.TransformPoint  (endPoint);
 
#ifdef  DEBUGSTEP
  ATH_MSG_DEBUG("   Local beginning step position "
                << startPointinLocal.x() << " " << startPointinLocal.y() << " "
                << startPointinLocal.z());
  ATH_MSG_DEBUG("   Local end       step position "
                << endPointinLocal.x() << " " << endPointinLocal.y() << " "
                << endPointinLocal.z());
#endif
 
  G4double energy = step->GetTotalEnergyDeposit() * step->GetTrack()->GetWeight(); // Despite the name, this is only ionization.
 
  if (m_BirksLaw) {
    const double EField = 10.;         // kV/cm, assume constant for now
    const double wholeStepLengthCm = step->GetStepLength() / CLHEP::cm;
    energy = (*m_birksLaw)(energy, wholeStepLengthCm, EField);
  }
 
 
  // compute number of sub steps
  //   = 1 if no charge collection
  //   otherwise 200 microns (dstep) division
 
  G4int nsub_step=1;
  if (m_IflCur) nsub_step=(int) (thisStepLength/m_dstep) + 1;
  // delta is fraction of step between two sub steps
  G4double delta=1./((double) nsub_step);
#ifdef DEBUGSTEP
  ATH_MSG_DEBUG("   nsub_step,delta " << nsub_step << " " << delta);
#endif
 
  energy /= ((float) (nsub_step)); // Divide into substeps.
 
  // loop over sub steps
 
  for (G4int i=0;i<nsub_step;i++) {
 
    bool Xtalk = false;
 
    G4double fraction=(((G4double) i)+0.5)*delta;
    G4ThreeVector subinLocal=startPointinLocal*(1-fraction) + endPointinLocal*fraction;
    G4double xloc   = subinLocal.x();
    G4double yloc   = subinLocal.y();
    G4double zloc   = subinLocal.z();
    G4double etaloc = subinLocal.pseudoRapidity();
    G4double philoc = subinLocal.phi();
    if(philoc<0.) philoc = philoc + 2.*M_PI;
    G4double radloc = sqrt( xloc*xloc + yloc*yloc );
#ifdef DEBUGSTEP
    ATH_MSG_DEBUG("   local position sub_step "
                  << xloc << " " << yloc << " " << zloc);
#endif
 
    //  apply fiducial cuts
    if (!FiducialCuts(radloc,zloc,etaloc)) {
#ifdef DEBUGSTEP
      ATH_MSG_DEBUG("LArBarrelCalculator: outside fiducial cuts");
#endif
      continue;
    }
    LArG4::Barrel::CalcData currentCellData;
    m_geometry->findCell(currentCellData,xloc,yloc,zloc,radloc,etaloc,philoc,m_IflCur);
 
    if (currentCellData.cellID == 0)
      {
#ifdef DEBUGSTEP
        ATH_MSG_DEBUG("LArBarrelCalculator: Invalid hit CELLID == 0 ");
        ATH_MSG_DEBUG("x,y,z local " << xloc << " " << yloc << " " << zloc);
        ATH_MSG_DEBUG("radius " << radloc <<" etaloc "<< etaloc << " philoc " << philoc);
#endif
        continue;
      }
    G4int region = currentCellData.region;
    G4int etaBin = currentCellData.etaBin;
    G4int phiBin = currentCellData.phiBin;
    G4int sampling = currentCellData.sampling;
 
    if( zSide == -1 )
      {
        // following code for an Y-axis rotation to define Z < 0. Barrel part
        if( sampling == 1 && region ==0 )
          {
            phiBin = 31 - phiBin;
            if(phiBin < 0 ) phiBin += 64;
          }
        if( sampling == 1 && region ==1 )
          {
            phiBin = 127 - phiBin;
            if(phiBin < 0 ) phiBin += 256;
          }
        if( sampling >= 2 )
          {
            phiBin = 127 - phiBin;
            if(phiBin < 0 ) phiBin += 256;
          }
      }
#ifdef DEBUGSTEP
    ATH_MSG_DEBUG("   region,side,sampling,eta,phi " << region << " " << zSide << " "
                  << sampling << " " << etaBin << " " << phiBin);
    ATH_MSG_DEBUG("   distance to electrode,abs " << currentCellData.distElec << " "
                  << currentCellData.distAbs);
    ATH_MSG_DEBUG("    local coordinates " << currentCellData.x0 << " " << currentCellData.y0);
#endif
 
    if (std::fabs(currentCellData.distElec)< m_ThickEle ||
        std::fabs(currentCellData.distAbs) < m_ThickAbs) {
#ifdef DEBUGSTEP
      ATH_MSG_DEBUG("   hit in absorber or electrode radius:" << radloc);
#endif
      continue;
    }
 
    identifier2.clear();
    identifier2 << 4          // LArCalorimeter
                << 1          // LArEM
                << zSide
                << sampling
                << region
                << etaBin
                << phiBin;
 
    // time computation is not necessarily correct for test beam
    //    G4double time;
    //    if (m_testbeam)
    //      {
    //        time=0.;
    //      }
    //    else
    //      {
    //        const G4double tof = thisStepPoint->GetGlobalTime() / Units::ns;
    //        time  = tof - ( thisStepPoint->GetPosition().mag()/CLHEP::c_light ) / Units::ns;
    //      }
 
    const G4double time = (m_testbeam)? 0.0 :
      ( thisStepPoint->GetGlobalTime() - ( thisStepPoint->GetPosition().mag()/CLHEP::c_light ) ) / Units::ns;
 
#ifdef DEBUGSTEP
    ATH_MSG_DEBUG("   Energy for sub step " << energy);
#endif
 
    G4double Current;
    G4double Current_xt1,Current_xt2;
    Current_xt1=0;
    Current_xt2=0;
    if (!m_IflCur)  {
      // no charge collection   Current=E from Geant
      Current=energy;
    }
    else  {
      // full charge collection
      G4double xmap,ymap;
      G4int nfold = currentCellData.nfold;
      G4double x0=currentCellData.x0;
      G4double y0=currentCellData.y0;
      if (x0<1500 || x0>1960 || y0>30 || y0<-30) {
        ATH_MSG_WARNING("weird x0,y0 " << x0 << " " << y0);
      }
#ifdef DEBUGSTEP
      G4double rr = sqrt(x0*x0+y0*y0);
      ATH_MSG_DEBUG("   radius,rad0 " << radloc << " " << rr);
#endif
      const CurrMap* mapOfCurrent(nullptr);
      bool UseFold=false;
      // Are we close to a fold ?   (fold 0 has some pathology)
      if ((x0 > m_accmap->GetXmin(nfold) || nfold==0) &&
          x0 < m_accmap->GetXmax(nfold)  &&
          y0 > m_accmap->GetYmin(nfold) &&
          y0 < m_accmap->GetYmax(nfold) &&
          (nfold != 13 || x0 < 1957.5) &&
          (nfold>0 || x0 < 1504.6) ) {
        xmap=x0;
        ymap=y0;
#ifdef DEBUGSTEP
        ATH_MSG_DEBUG("   Map for fold xmap,ymap " << nfold << " " << xmap << " " << ymap);
#endif
        UseFold=true;
        G4int sampMap = currentCellData.sampMap;
        G4int etaMap  = currentCellData.etaMap;
        mapOfCurrent = m_accmap->GetMap(nfold,region,sampMap,etaMap);
        // catch problem to find map
        if (!mapOfCurrent) {
          ATH_MSG_WARNING(" Problem to access map fold = " << nfold);
          ATH_MSG_WARNING(" region,sampling,eta,fold " << region << " " << sampMap << " "
                          << etaMap << " " << nfold);
          return false;
        }
      }
      else {
        G4int n;
        G4int nstraight = currentCellData.nstraight;
        if (nstraight%2==0) n=22;
        else                n=21;
        mapOfCurrent = m_accmap->GetMap(n,region,sampling,etaBin);
        // catch problem to find map
        if (!mapOfCurrent) {
          ATH_MSG_WARNING(" Problem to access map straight = " << nstraight);
          return false;
        }
        xmap = currentCellData.xl;
        ymap = currentCellData.distElec;
        // special case for first straight section, which is shorter
        if (nstraight==0) xmap = 0.5*(xmap+1.);
#ifdef DEBUGSTEP
        ATH_MSG_DEBUG("   Map for straight xl,delec " << xmap << " " << ymap);
#endif
      }   // fold or straight
      double gap;
      double current0,current1,current2;
      //  get current for elementary charge
      mapOfCurrent->GetAll(xmap,ymap,&gap,&current0,&current1,&current2);
      G4double gap2=std::fabs(currentCellData.distElec)+std::fabs(currentCellData.distAbs)
        -m_ThickEle-m_ThickAbs;
 
      // in which HV cell are we ?
      int ipm,ielec,ieta,iside;
      if (zSide==1) ipm=1;    // A side
      else          ipm=0;    // C side
      ielec=currentCellData.phiGap;
      if (zSide==-1) {
        ielec = 511 - ielec;
        if(ielec < 0 ) ielec += 1024;
      }
      ieta=((int) (etaloc*(1./0.2)));
      if (ieta>6) ieta=6;  //part 1.4 to 1.475 is same HV as 1.2 to 1.4
      iside=0;    // phi lower than electrode 0, 1 for phi higher than electrode
      if ((currentCellData.distElec>0 && zSide==1)
          || (currentCellData.distElec<0 && zSide==-1) ) iside=1;
 
      // HV extrapolation
      double current;
      double hv=m_hv[ipm][ielec][ieta][iside];
      if (hv>1995.) current=current0;
      else if (hv>5.) current=ScaleHV(hv,current0,current1,current2);
      else current=0.;
 
      // extrapolation for non nominal gap (allows to include sagging effect)
      //  i ~ (gap/gap2)**1.3
      //  gap = nominal gap from current map
      //  gap2 = real gap from geometry
      //  dgap = gap/gap2 -1
      //   at first order  i ~ (1+1.3*dgap)
      double dgap=0;
      if (gap>1e-3 && gap2 >1e-3) dgap=gap/gap2-1;
      current = current * (1. + 1.3*dgap);
 
#ifdef DEBUGSTEP
      ATH_MSG_DEBUG("   elementary current " << current0);
      ATH_MSG_DEBUG("   Gap from Map/calculator " << gap << " " << gap2);
#endif
 
      Current = energy*current;
 
      // check if pathology...
      if (!UseFold && std::fabs(currentCellData.distElec)>2.1 && current0 < 0.1) {
        ATH_MSG_WARNING(" xl,distEle " << currentCellData.xl << " "
                        << currentCellData.distElec
                        << " str number " << currentCellData.nstraight
                        << " sampling,eta " << sampling << " " << etaBin << " "
                        << " current/E " << current0);
      }
 
      // compute effect around transition at eta=0.8 from electric field
      // read from map weighting factor ( E**1.3 = E*E**0.3)
      // in x  +-9mm around transition (assume symmetry around 0)
      // in y  distance from electrode in gap
      bool InTrans=false;
      //double Current_test=Current;
      if (m_IflMapTrans) {
        float etaTrans=0.8;
        if (fabs(etaloc-etaTrans) < 0.025) {
          double x=std::fabs(zloc-radloc*sinh(etaTrans))/cosh(etaTrans);
          double y=std::fabs(currentCellData.distElec);
          if ( x < m_etamap3->Xmax() ) {
            double resp;
            m_etamap3->GetData0(x,y,&resp);
            Current = Current*resp;
            InTrans=true;
          }
        }     // eta = 0.8 +- 0.025
      }       // if m_IflMapTrans
 
      // simulate cross-talk effects and transitions in eta between cells
      // (only in region 0, samplings 1 and 2)
 
      if (region==0) {
 
        double resp,xt0,xt1,xt2,deta;
        if (sampling==1) {
          deta=etaloc-0.003125*((double)etaBin+0.5);
          m_etamap1->GetData(std::fabs(deta),std::fabs(currentCellData.distElec),
                             &resp,&xt0,&xt1,&xt2);
#ifdef DEBUGSTEP
          ATH_MSG_DEBUG("hit in strip etaloc,etaBin,deta,delec,resp,xt0,xt1,xt2 "
                        << etaloc << " " << etaBin << " " << deta*1000 << " " << currentCellData.distElec
                        << " " << resp << " " << xt0 << " " << xt1 << " " << xt2);
#endif
          if (!InTrans) Current = Current*resp;
          if (m_IflXtalk && etaBin > 1 && etaBin < 446) {
            Xtalk=true;
            if (deta>0) {
              identifier_xt1.clear();
              identifier_xt1 << 4 << 1 << zSide << sampling << region << (etaBin+1) << phiBin;
              Current_xt1 = Current*xt1;
              identifier_xt2.clear();
              identifier_xt2 << 4 << 1 << zSide << sampling << region << (etaBin-1) << phiBin;
              Current_xt2 = Current*xt2;
            }
            else {
              identifier_xt1.clear();
              identifier_xt1 << 4 << 1 << zSide << sampling << region << (etaBin-1) << phiBin;
              Current_xt1 = Current*xt1;
              identifier_xt2.clear();
              identifier_xt2 << 4 << 1 << zSide << sampling << region << (etaBin+1) << phiBin;
              Current_xt2 = Current*xt2;
            }
            Current = Current*xt0;
          }    // m_IflXtalk = true
        }
        else if (sampling==2 && !InTrans) {
          deta=etaloc-0.025*((double)etaBin+0.5);
          m_etamap2->GetData0(std::fabs(deta),std::fabs(currentCellData.distElec),
                              &resp);
#ifdef DEBUGSTEP
          ATH_MSG_DEBUG("hit in middle etaloc,etaBin,deta,delec,resp,xt0,xt1,xt2 "
                        << etaloc << " " << etaBin << " " << deta*1000 << " " << currentCellData.distElec
                        << " " << resp);
#endif
          Current = Current*resp;
        }   // sampling =1 or 2
 
      }      // region=0
 
    }    // switch for current simulation
 
 
    // check if we have a new hit in a different cell, or if we add this substep
    //  to an already existing hit
    G4bool found=false;
    for (unsigned int i=0; i<hdata.size(); i++) {
      if (hdata[i].id==identifier2) {
        hdata[i].energy += Current;
        hdata[i].time += time*Current;
        found=true;
        break;
      }
    }    // loop over hits
    if (!found) {
      hdata.emplace_back(LArHitData{identifier2, time*Current, Current});
    }    // hit was not existing before
 
    if (Xtalk) {
      for (unsigned int i=0; i<hdata.size(); i++) {
        if (hdata[i].id==identifier_xt1) {
          hdata[i].energy += Current_xt1;
          hdata[i].time += time*Current_xt1;
          found=true;
          break;
        }
      }    // loop over hits
      if (!found) {
        hdata.emplace_back(LArHitData{identifier_xt1, time*Current_xt1, Current_xt1});
      }
      found=false;
      for (unsigned int i=0; i<hdata.size(); i++) {
        if (hdata[i].id==identifier_xt2) {
          hdata[i].energy += Current_xt2;
          hdata[i].time += time*Current_xt2;
          found=true;
          break;
        }
      }    // loop over hits
      if (!found) {
        hdata.emplace_back(LArHitData{identifier_xt2, time*Current_xt2, Current_xt2});
      }
    }    // Xtalk true
 
 
  }    // *** End of loop over sub steps
 
#ifdef DEBUGSTEP
  ATH_MSG_DEBUG("Number of hits for this step " << m_nhits << " "
                << hdata.size());
#endif
 
  // finalise time computations
  for (unsigned int i=0;i<hdata.size();i++) {
    if (std::fabs(hdata[i].energy)>1e-6) hdata[i].time=hdata[i].time/hdata[i].energy;
    else hdata[i].time=0.;
#ifdef DEBUGSTEP
    ATH_MSG_DEBUG("Hit Energy/Time "
                  << hdata[i].energy << " " << hdata[i].time);
#endif
  }
 
  if (hdata.size()>0) return true;
  return false;
}

ScaleHV()

double LArBarrelCalculator::ScaleHV ( double  hv, double  curr0, double  curr1, double  curr2  ) const

private

Definition at line 698 of file LArBarrelCalculator.cxx.

{
  double b;
  double resp=0.;
  if (hv>1000.) {
    if (std::fabs(curr1)>1e-6) {
      double x=curr0/curr1;
      if (x>1e-6) {
        b=log(x)*(1./(log(2000.)-log(1000.)));
        resp = curr0*exp(b*(log(hv)-log(2000.)));
      }
    }
  }
  else if (hv>400.) {
    if (std::fabs(curr2)>1e-6) {
      double x=curr1/curr2;
      if (x>1e-6) {
        b=log(x)*(1./(log(1000.)-log(400.)));
        resp = curr1*exp(b*(log(hv)-log(1000.)));
      }
    }
  }
  else {
    if (std::fabs(curr2)>1e-6) {
      b=0.57;
      resp = curr2*exp(b*(log(hv)-log(400.)));
    }
  }
  return resp;
}

Member Data Documentation

m_accmap

const AccMap* LArBarrelCalculator::m_accmap {nullptr}

private

Definition at line 61 of file LArBarrelCalculator.h.

m_Birksk

Gaudi::Property<double> LArCalculatorSvcImp::m_Birksk {this, "Birksk", 0.05832}

protectedinherited

Definition at line 27 of file LArCalculatorSvcImp.h.

m_BirksLaw

Gaudi::Property<bool> LArCalculatorSvcImp::m_BirksLaw {this, "BirksLaw", true}

protectedinherited

Definition at line 23 of file LArCalculatorSvcImp.h.

m_birksLaw

const LArG4BirksLaw* LArBarrelCalculator::m_birksLaw {nullptr}

private

Definition at line 73 of file LArBarrelCalculator.h.

m_doHV

Gaudi::Property<bool> LArBarrelCalculator::m_doHV {this, "EMBHVEnable", false}

private

Definition at line 74 of file LArBarrelCalculator.h.

m_dstep

Gaudi::Property<double> LArBarrelCalculator::m_dstep {this, "EMBdstep", .2*CLHEP::mm}

private

Definition at line 71 of file LArBarrelCalculator.h.

m_etamap1

std::unique_ptr<MapEta> LArBarrelCalculator::m_etamap1 {}

private

Definition at line 62 of file LArBarrelCalculator.h.

m_etamap2

std::unique_ptr<MapEta> LArBarrelCalculator::m_etamap2 {}

private

Definition at line 63 of file LArBarrelCalculator.h.

m_etamap3

std::unique_ptr<MapEta> LArBarrelCalculator::m_etamap3 {}

private

Definition at line 64 of file LArBarrelCalculator.h.

m_etaMaxBarrel

double LArBarrelCalculator::m_etaMaxBarrel {0.}

private

Definition at line 77 of file LArBarrelCalculator.h.

m_geometry

ServiceHandle<ILArBarrelGeometry> LArBarrelCalculator::m_geometry {this, "GeometryCalculator", "LArBarrelGeometry"}

private

Definition at line 60 of file LArBarrelCalculator.h.

m_hv

double LArBarrelCalculator::m_hv[2][1024][7][2] {}

private

Definition at line 100 of file LArBarrelCalculator.h.

m_IflCur

Gaudi::Property<bool> LArBarrelCalculator::m_IflCur {this, "EMBCurr", true}

private

Definition at line 67 of file LArBarrelCalculator.h.

m_IflMapTrans

Gaudi::Property<bool> LArBarrelCalculator::m_IflMapTrans {this, "EMBEtaTrans", true}

private

Definition at line 68 of file LArBarrelCalculator.h.

m_IflXtalk

Gaudi::Property<bool> LArBarrelCalculator::m_IflXtalk {this, "EMBXtalk", true}

private

Definition at line 69 of file LArBarrelCalculator.h.

m_NCellMax

int LArBarrelCalculator::m_NCellMax {0}

private

Definition at line 88 of file LArBarrelCalculator.h.

m_NCellTot

int LArBarrelCalculator::m_NCellTot {0}

private

Definition at line 87 of file LArBarrelCalculator.h.

m_OOTcut

Gaudi::Property<double> LArCalculatorSvcImp::m_OOTcut {this, "OOTcut", 300*CLHEP::ns}

protectedinherited

Definition at line 30 of file LArCalculatorSvcImp.h.

m_rMaxAccordion

double LArBarrelCalculator::m_rMaxAccordion {0.}

private

Definition at line 82 of file LArBarrelCalculator.h.

m_rMinAccordion

double LArBarrelCalculator::m_rMinAccordion {0.}

private

Definition at line 81 of file LArBarrelCalculator.h.

m_testbeam

bool LArBarrelCalculator::m_testbeam {false}

private

Definition at line 91 of file LArBarrelCalculator.h.

m_ThickAbs

double LArBarrelCalculator::m_ThickAbs {0.}

private

Definition at line 84 of file LArBarrelCalculator.h.

m_ThickEle

double LArBarrelCalculator::m_ThickEle {0.}

private

Definition at line 85 of file LArBarrelCalculator.h.

m_zMaxBarrel

double LArBarrelCalculator::m_zMaxBarrel {0.}

private

Definition at line 79 of file LArBarrelCalculator.h.

m_zMinBarrel

double LArBarrelCalculator::m_zMinBarrel {0.}

private

Definition at line 78 of file LArBarrelCalculator.h.

---

The documentation for this class was generated from the following files:

• LArBarrelCalculator.h
• LArBarrelCalculator.cxx