diff --git a/LArCalorimeter/LArG4/LArG4EC/src/BarretteEnergyCalculator.cc b/LArCalorimeter/LArG4/LArG4EC/src/BarretteEnergyCalculator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..06181a9adb2b151d1a5044f199121b45df91e809
--- /dev/null
+++ b/LArCalorimeter/LArG4/LArG4EC/src/BarretteEnergyCalculator.cc
@@ -0,0 +1,827 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <cmath>
+#include <cassert>
+#include <string>
+
+#include "G4AffineTransform.hh"
+#include "G4NavigationHistory.hh"
+#include "G4VTouchable.hh"
+#include "G4TouchableHistory.hh"
+#include "G4ThreeVector.hh"
+#include "G4Step.hh"
+#include "globals.hh"
+
+#include "LArG4Code/LArG4BirksLaw.h"
+#include "LArG4Code/ILArCalibCalculatorSvc.h"
+
+#include "EnergyCalculator.h"
+#include "AthenaKernel/Units.h"
+#include "HVHelper.h"
+
+using namespace LArG4::EC;
+namespace Units = Athena::Units;
+
+const G4double EnergyCalculator::s_ColdCorrection = 1.0036256;
+const G4double EnergyCalculator::s_LongBarThickness = 20.*Units::mm;
+const G4double EnergyCalculator::s_StripWidth = 3.*Units::mm/s_ColdCorrection;
+const G4double EnergyCalculator::s_KapGap = 1.*Units::mm/s_ColdCorrection;
+const G4double EnergyCalculator::s_EdgeWidth = 1.*Units::mm;
+const G4double EnergyCalculator::s_DistOfEndofCuFromBack = 22.77*Units::mm/s_ColdCorrection;
+const G4double EnergyCalculator::s_DistOfStartofCuFromBack= 31.*Units::mm; // frontface of the barrette
+const G4double EnergyCalculator::s_ZmaxOfSignal = s_DistOfStartofCuFromBack
+ - s_DistOfEndofCuFromBack + s_EdgeWidth;
+G4double EnergyCalculator::s_RefzDist = 0.;
+G4bool EnergyCalculator::s_SetConstOuterBarrett = false;
+G4bool EnergyCalculator::s_SetConstInnerBarrett = false;
+
+const G4double EnergyCalculator::s_S3_Etalim[21]={
+ 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95,
+ 2.00, 2.05, 2.10, 2.15, 2.20, 2.25, 2.30, 2.35, 2.40, 2.45, 2.5
+};
+const G4double EnergyCalculator::s_Rmeas_outer[50]={
+ 11.59, 25.22, 57.28, 71.30, 85.90, 98.94, 103.09, 116.68, 130.42, 146.27,
+ 147.19, 11.59, 15., 56.91, 44.37, 15.13, 14.93, 45.87, 35.03, 15.40,
+ 14.04, 39.67, 26.83, 15.64, 14.90, 30.26, 14.70, 29.09, 43.12, 34.51,
+ 25.08, 11.88, 14.39, 19.54, 17.80, 12.70, 15.31, 13.96, 11.79, -99.,
+ 23.57, 34.64, 55.32, 65.39, 76.34, 10.83, 94.84, 98.00, -99., -99.
+};
+const G4double EnergyCalculator::s_Zmeas_outer[2]={3.81, 7.81};
+
+G4double EnergyCalculator::s_S3_Rlim[21];
+G4double EnergyCalculator::s_rlim[50];
+G4double EnergyCalculator::s_zlim[4];
+
+
+G4double EnergyCalculator::getPhiStartOfPhiDiv(const G4Step* step) const {
+ // --- moved from FindIdentifier_Barrett --- //
+ // get Module and Phidiv number (result is put into VolumeNumber and PhiDivNumber)
+ G4int ModuleNumber = -999;
+ G4int PhiDivNumber = -999;
+
+ G4bool validhit = GetVolumeIndex(step, ModuleNumber, PhiDivNumber);
+ if (!validhit) {
+ ATH_MSG_FATAL(" ERROR ::Process_Barrett::Module, Phidiv is not found"
+ <<" ModuleNumber= "<<ModuleNumber<<" PhiDivNumber= "<<PhiDivNumber);
+ }
+
+ if(lwc()->type() == LArG4::OuterAbsorberWheel){ // for wheel calculation
+ //PhiStartOfPhiDiv =
+ return lwc()->GetFanStepOnPhi()/2. + ModuleNumber * CLHEP::twopi/8.
+ + PhiDivNumber * lwc()->GetFanStepOnPhi();
+ //phi_inb = PhiStartOfPhiDiv + pinLocal.phi(); //in ::BackOuterBarrettes
+ //if(phi_inb < 0.) phi_inb = phi_inb + CLHEP::twopi;
+ //if(phi_inb > CLHEP::twopi) phi_inb = phi_inb - CLHEP::twopi;
+ //phi_inb = CLHEP::twopi - phi_inb; // phi in ::EmecMother system;
+ }
+
+ if(lwc()->type() == LArG4::OuterAbsorberModule){ // for TB modul calculation
+
+ // M_PI_2 - M_PI/8. is from EMECSupportConstruction
+ // PhiStartOfPhiDiv =
+ return M_PI_2 - M_PI/8. + lwc()->GetFanStepOnPhi()/2 + PhiDivNumber * lwc()->GetFanStepOnPhi();
+
+ //phi_inb = PhiStartOfPhiDiv + pinLocal.phi(); //in BackOuterBarrettes;
+ //phi_inb = M_PI - phi_inb; // phi in ::EmecMother system;
+ }
+ // --- moved from FindIdentifier_Barrett --- //
+ ATH_MSG_ERROR(" ERROR ::Process_Barrett::getPhiStartOfPhiDiv, LArWheelType is unknown ("<<lwc()->type()<<"), PhiStartOfPhiDiv is set zero!"
+ <<" ModuleNumber= "<<ModuleNumber<<" PhiDivNumber= "<<PhiDivNumber);
+ return 0.0;
+}
+
+// ****************************************************************************
+G4bool EnergyCalculator::Process_Barrett(const G4Step* step, std::vector<LArHitData>& hdata) const{
+ // ****************************************************************************
+ const G4double PhiStartOfPhiDiv = getPhiStartOfPhiDiv(step);
+ G4ThreeVector startPointinLocal, endPointinLocal;
+ if(!FindIdentifier_Barrett(step, PhiStartOfPhiDiv, hdata, startPointinLocal, endPointinLocal)){
+ //cell id is not found
+ if(!lwc()->GetisBarretteCalib()) return false; // ==> normal calculator return false
+ else{ // if Calibration Calculator for Barrett:
+ // compute DM identifier;
+ // In calibration calculator mode this process is used as "geometry calculator".
+ // We have to be prepared to return
+ // either the emec cell id or , if the hit does not belong to any cell,
+ // the DeadMatter (DM) id.
+ // DM identifier for Barrett is defined
+ // by the EmecSupportCalibrationCalculator.( In general it is
+ // activated by the atlas_calo.py)
+ const G4bool status=FindDMIdentifier_Barrett(step, hdata);
+ if(!status)
+ ATH_MSG_WARNING(" WARNING!!EnergyCalculator::Process_Barrett:"
+ << " DM identifier not found????");
+ return status;
+ } // end of calibr. calculator
+ } // end if cell id not found
+ // compute signal in 'normal' calculator mode, if cellid is found
+ G4double E = step->GetTotalEnergyDeposit() * step->GetTrack()->GetWeight();
+ if (m_birksLaw){
+ const G4ThreeVector midpoint = (startPointinLocal + endPointinLocal) * 0.5;
+ const G4double wholeStepLengthCm = step->GetStepLength() / CLHEP::cm;
+ const G4double gap = GetGapSize_Barrett(midpoint); // const method
+ G4double phi;
+ G4int compartment, eta_bin;
+ if(!GetBarrettePCE(midpoint, PhiStartOfPhiDiv, phi, compartment, eta_bin)){
+ hdata[0].energy = 0.;
+ return true;
+ }
+ const G4double efield =
+ 0.01 * m_HVHelper->GetVoltageBarrett(phi, compartment, eta_bin) / gap; // estimate Efield[KV/cm]
+ E = (*m_birksLaw)(E, wholeStepLengthCm, efield);
+ }
+
+ hdata[0].energy = (this->*m_ecorr_method)(E, startPointinLocal, endPointinLocal, PhiStartOfPhiDiv);
+ return true;
+}
+
+// ****************************************************************************
+void EnergyCalculator::SetConst_InnerBarrett(void){
+ // ****************************************************************************
+ if(s_SetConstInnerBarrett) return;
+ s_SetConstInnerBarrett=true;
+
+ std::cout <<" ===>>> ERROR!! SetConst_InnerBarrett is called!!!" <<std::endl;
+ exit(99);
+}
+
+// ****************************************************************************
+void EnergyCalculator::SetConst_OuterBarrett(void){
+ // ****************************************************************************
+
+ if(s_SetConstOuterBarrett) return;
+ s_SetConstOuterBarrett=true;
+
+ for(G4int i=0;i<=20;++i){
+ const G4double teta = 2.*atan( exp(-s_S3_Etalim[i]));
+ s_S3_Rlim[i] = s_RefzDist*tan(teta);
+ }
+
+ const G4double inv_ColdCorrection = 1. / s_ColdCorrection;
+ s_rlim[0] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[0] /*11.59 */ * inv_ColdCorrection;
+ s_rlim[1] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[1] /*25.22 */ * inv_ColdCorrection;
+ s_rlim[2] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[2] /*57.28 */ * inv_ColdCorrection;
+ s_rlim[3] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[3] /*71.30 */ * inv_ColdCorrection;
+ s_rlim[4] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[4] /*85.90 */ * inv_ColdCorrection;
+ s_rlim[5] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[5] /*98.94 */ * inv_ColdCorrection;
+ s_rlim[6] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[6] /*103.09 */ * inv_ColdCorrection;
+ s_rlim[7] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[7] /*116.68 */ * inv_ColdCorrection;
+ s_rlim[8] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[8] /*130.42 */ * inv_ColdCorrection;
+ s_rlim[9] = s_S3_Rlim[3] + s_KapGap/2. + s_Rmeas_outer[9] /*146.27 */ * inv_ColdCorrection + s_EdgeWidth;
+ s_rlim[10]= s_rlim[8] + s_Rmeas_outer[10]/*147.19 */ * inv_ColdCorrection;
+
+ s_rlim[11]= s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[11]/*11.59 */ * inv_ColdCorrection; //eta=1.65
+ s_rlim[12]= s_S3_Rlim[3] - s_KapGap - s_Rmeas_outer[12]/*15. */ * inv_ColdCorrection;
+
+ s_rlim[13]= s_S3_Rlim[4] + s_KapGap + s_Rmeas_outer[13]/*56.91 */ * inv_ColdCorrection; //eta=1.7
+ s_rlim[14]= s_S3_Rlim[4] - s_KapGap - s_Rmeas_outer[14]/*44.37 */ * inv_ColdCorrection;
+
+ s_rlim[15]= s_S3_Rlim[5] + s_KapGap + s_Rmeas_outer[15]/*15.13 */ * inv_ColdCorrection; //eta=1.75
+ s_rlim[16]= s_S3_Rlim[5] - s_KapGap - s_Rmeas_outer[16]/*14.93 */ * inv_ColdCorrection;
+
+ s_rlim[17]= s_S3_Rlim[6] + s_KapGap + s_Rmeas_outer[17]/*45.87 */ * inv_ColdCorrection; //eta=1.80
+ s_rlim[18]= s_S3_Rlim[6] - s_KapGap - s_Rmeas_outer[18]/*35.03 */ * inv_ColdCorrection;
+
+ s_rlim[19]= s_S3_Rlim[7] + s_KapGap + s_Rmeas_outer[19]/*15.40 */ * inv_ColdCorrection; //eta=1.85
+ s_rlim[20]= s_S3_Rlim[7] - s_KapGap - s_Rmeas_outer[20]/*14.04 */ * inv_ColdCorrection;
+
+ s_rlim[21]= s_S3_Rlim[8] + s_KapGap + s_Rmeas_outer[21]/*39.67 */ * inv_ColdCorrection; //eta=1.90
+ s_rlim[22]= s_S3_Rlim[8] - s_KapGap - s_Rmeas_outer[22]/*26.83 */ * inv_ColdCorrection;
+
+ s_rlim[23]= s_S3_Rlim[9] + s_KapGap + s_Rmeas_outer[23]/*15.64 */ * inv_ColdCorrection; //eta=1.95
+ s_rlim[24]= s_S3_Rlim[9] - s_KapGap - s_Rmeas_outer[24]/*14.90 */ * inv_ColdCorrection;
+
+ s_rlim[25]= s_S3_Rlim[10] + s_KapGap + s_Rmeas_outer[25]/*30.26 */ * inv_ColdCorrection; //eta=2.00
+ s_rlim[26]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[26]/*14.70 */ * inv_ColdCorrection;
+
+ s_rlim[27]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[27]/*29.09 */ * inv_ColdCorrection; //eta=2.05
+ s_rlim[28]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[28]/*43.12 */ * inv_ColdCorrection; //SHAPE CHANGE!!ZZ
+
+ s_rlim[29]= s_S3_Rlim[12] + s_KapGap + s_Rmeas_outer[29]/*34.51 */ * inv_ColdCorrection; //eta=2.10
+ s_rlim[30]= s_S3_Rlim[12] - s_KapGap - s_Rmeas_outer[30]/*25.08 */ * inv_ColdCorrection;
+
+ s_rlim[31]= s_S3_Rlim[13] + s_KapGap + s_Rmeas_outer[31]/*11.88 */ * inv_ColdCorrection; //eta=2.15
+ s_rlim[32]= s_S3_Rlim[13] - s_KapGap - s_Rmeas_outer[32]/*14.39 */ * inv_ColdCorrection;
+
+ s_rlim[33]= s_S3_Rlim[14] + s_KapGap + s_Rmeas_outer[33]/*19.54 */ * inv_ColdCorrection; //eta=2.20
+ s_rlim[34]= s_S3_Rlim[14] - s_KapGap - s_Rmeas_outer[34]/*17.80 */ * inv_ColdCorrection; // !!ZZ
+
+ s_rlim[35]= s_S3_Rlim[15] + s_KapGap + s_Rmeas_outer[35]/*12.70 */ * inv_ColdCorrection; //eta=2.25
+ s_rlim[36]= s_S3_Rlim[15] - s_KapGap - s_Rmeas_outer[36]/*15.31 */ * inv_ColdCorrection;
+
+ s_rlim[37]= s_S3_Rlim[16] + s_KapGap + s_Rmeas_outer[37]/*13.96 */ * inv_ColdCorrection; //eta=2.30
+ s_rlim[38]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[38]/*11.79 */ * inv_ColdCorrection; // !!ZZ!!
+
+ s_rlim[40]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[40]/*23.57 */ * inv_ColdCorrection;
+ s_rlim[41]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[41]/*34.64 */ * inv_ColdCorrection;
+ s_rlim[42]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[42]/*55.32 */ * inv_ColdCorrection;
+ s_rlim[43]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[43]/*65.39 */ * inv_ColdCorrection;
+ s_rlim[44]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[44]/*76.34 */ * inv_ColdCorrection;
+ s_rlim[45]= s_rlim[44] - s_Rmeas_outer[45]/*10.83 */ * inv_ColdCorrection;
+ s_rlim[46]= s_S3_Rlim[16] - s_KapGap/2.- s_Rmeas_outer[46]/*94.84 */ * inv_ColdCorrection - s_EdgeWidth;
+ s_rlim[47]= s_S3_Rlim[16] - s_KapGap/2.- s_Rmeas_outer[47]/*98.00 */ * inv_ColdCorrection - s_EdgeWidth;
+
+ s_zlim[0] = - s_EdgeWidth + s_Zmeas_outer[0]/*3.81*/ * inv_ColdCorrection; //rel. to the end_of_HV_Cu
+ s_zlim[1] = - s_KapGap/2. + s_Zmeas_outer[1]/*7.81*/ * inv_ColdCorrection;
+ s_zlim[2] = s_StripWidth + 1./2. * s_KapGap;
+ s_zlim[3] =2*s_StripWidth + 3./2. * s_KapGap;
+
+ for (G4int i=0; i<=3; ++i){
+ s_zlim[i]= (s_ZmaxOfSignal-s_EdgeWidth) - s_zlim[i]; //rel to start of the Barrette
+ if(s_zlim[i] < 0.) s_zlim[i]=0.;
+ }
+
+ return;
+}
+// ****************************************************************************
+G4bool EnergyCalculator::GetCompartment_Barrett(
+ G4ThreeVector pforcell, G4double r_inb, G4double z_inb, G4double eta_inb,
+ G4int &b_compartment, G4int &etabin) const {
+ // ****************************************************************************
+
+ G4double d,d1,d2,rlim1,rlim2,rlim3,zlim1,zlim2,eta1,eta2;
+ G4int i;
+ G4int i0 = 3;
+
+ G4bool validhit=true;
+ b_compartment=-99;
+ etabin=-99;
+
+ if(r_inb > s_rlim[10] || r_inb < s_rlim[47] )
+ {validhit=false; goto label99;}
+ if(z_inb > s_ZmaxOfSignal ){validhit=false; goto label99;}
+
+ if(r_inb > s_rlim[0]) { // Upper corner
+
+ if(r_inb > s_rlim[9]) {
+ if(z_inb > s_zlim[0]) {validhit=false; goto label99;}
+ label1:
+ if(z_inb > s_zlim[1]) {
+ b_compartment = 8;
+ etabin = 4;
+ goto label99;
+ }
+ label2:
+ b_compartment = 9;
+ etabin = 0;
+ goto label99;
+ }
+ if(r_inb > s_rlim[8]) goto label1;
+ if(r_inb > s_rlim[7]) goto label2;
+ if(r_inb > s_rlim[6]) {
+ b_compartment = 8;
+ etabin = 5;
+ goto label99;
+ }
+ if(r_inb > s_rlim[5]) {
+ label3:
+ b_compartment = 8;
+ etabin = 6;
+ goto label99;
+ }
+ if(r_inb > s_rlim[4]) {
+ if(z_inb > s_zlim[1]) goto label3;
+
+ label4:
+ b_compartment = 9;
+ etabin = 1;
+ goto label99;
+ }
+
+ if(r_inb > s_rlim[3]) goto label4;
+
+ if(r_inb > s_rlim[2]) {
+ b_compartment = 8;
+ etabin = 7;
+ goto label99;
+ }
+
+ if(r_inb > s_rlim[1]) {
+ b_compartment = 8;
+ etabin = 8;
+ goto label99;
+ }
+ if(r_inb > s_rlim[0]) {
+ b_compartment = 9;
+ etabin = 2;
+ goto label99;
+ }
+ }
+
+ if(r_inb < s_rlim[38]){ // lower corner
+
+ if( r_inb > s_rlim[40] ) {
+ b_compartment = 9;
+ etabin = 16;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[41] ) {
+ label5:
+ b_compartment = 8;
+ etabin = 37;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[42] ) {
+
+ d=DistanceToEtaLine( pforcell,2.35);
+ if(fabs(d) < s_StripWidth+s_KapGap) {
+ if( d < 0.) {
+
+ label6:
+ b_compartment = 8;
+ etabin = 38;
+ goto label99;
+ }
+ label7:
+ if( z_inb < s_zlim[3] ) goto label5;
+ goto label6;
+ }
+
+ if( d > 0. ) goto label7;
+
+ b_compartment = 9;
+ etabin = 17;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[43] ) {
+ b_compartment = 8;
+ etabin = 39;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[44] ) {
+ label8:
+ b_compartment = 8;
+ etabin = 40;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[45] ) {
+ if( z_inb < s_zlim[3] ) goto label8;
+
+ label9:
+ b_compartment = 9;
+ etabin = 18;
+ goto label99;
+ }
+
+ if( r_inb > s_rlim[46] ) goto label9;
+
+ if( z_inb < s_ZmaxOfSignal/(s_rlim[46]-s_rlim[47])*(r_inb-s_rlim[47])) goto label9;
+
+ validhit=false;
+ goto label99;
+ }
+
+ // medium r region: s_rlim[0] > r > s_rlim[38];
+ // from middle of cellno 2 to middle of cellno. 16
+ //
+
+ for( i=3; i<=17; ++i){ // eta= 1.65 - 2.35
+ if( eta_inb < s_S3_Etalim[i] ) {
+ i0=i;
+ break;
+ }
+ }
+
+ i=i0;
+
+ eta1 = s_S3_Etalim[i-1];
+ eta2 = s_S3_Etalim[i];
+ rlim1 = s_rlim[2*i+5 - 1];
+ rlim2 = s_rlim[2*i+5];
+ zlim1 = s_zlim[2];
+ zlim2 = s_zlim[3];
+
+ if( i == 15 || i == 17) {
+ zlim1 = s_zlim[3];
+ zlim2 = s_zlim[2];
+ }
+
+ switch(i){
+
+ case 3:
+
+ if( fabs( DistanceToEtaLine(pforcell, eta2) ) < s_StripWidth+s_KapGap
+ || z_inb > zlim1 ){
+
+ b_compartment = 8;
+ etabin = 2*i+3;;
+ break;
+ }
+
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+
+ case 4:
+ case 5:
+ case 6:
+ case 7:
+ case 8:
+ case 9:
+ case 10:
+ case 13:
+ case 16:
+
+ d1=fabs( DistanceToEtaLine( pforcell, eta1) );
+ d2=fabs( DistanceToEtaLine( pforcell, eta2) );
+
+ if( d1 < s_StripWidth+s_KapGap ){
+ label11:
+ b_compartment = 8;
+ etabin = 2*i+2;
+ break;
+ }
+ if( d2 < s_StripWidth+s_KapGap ){
+ label12:
+ b_compartment = 8;
+ etabin = 2*i+3;
+ break;
+ }
+
+ if( z_inb < zlim1 || (r_inb > rlim2 && r_inb < rlim1) ) {
+
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+ }
+ if( r_inb > rlim1 ) goto label11;
+ if( r_inb < rlim2 ) goto label12;
+ validhit=false;
+ break;
+ //========================================================
+ case 11:
+
+ rlim3 = s_rlim[28];
+
+ d1=fabs( DistanceToEtaLine( pforcell, eta1) );
+ d2=fabs( DistanceToEtaLine( pforcell, eta2) );
+
+ if( d1 < s_StripWidth+s_KapGap ){
+ label13:
+ b_compartment = 8;
+ etabin = 2*i+2;
+ break;
+ }
+
+ if( r_inb > rlim1 && z_inb > zlim1 ) goto label13;
+ if( r_inb > rlim2 ){
+ label14:
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+ }
+
+ if( d2 < s_StripWidth+s_KapGap ){
+ if( z_inb > zlim1 ) {
+ label15:
+ b_compartment = 8;
+ etabin = 2*i+4;
+ break;
+ }
+
+ label16:
+ b_compartment = 8;
+ etabin = 2*i+3;
+ break;
+ }
+
+ if( z_inb < zlim2) goto label14;
+ if( r_inb > rlim3) goto label16;
+ if( z_inb > zlim1) goto label15;
+ goto label16;
+ //======================================================
+ case 12:
+
+ d1=fabs( DistanceToEtaLine( pforcell, eta1) );
+ d2=fabs( DistanceToEtaLine( pforcell, eta2) );
+
+ if( d1 < s_StripWidth+s_KapGap ){
+ b_compartment = 8;
+ etabin = 2*i+2;
+ break;
+ }
+ if( d2 < s_StripWidth+s_KapGap ){
+ label17:
+ b_compartment = 8;
+ etabin = 2*i+3;
+ break;
+ }
+
+ if( r_inb > rlim2 || z_inb < zlim1 ){
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+ }
+ goto label17;
+
+ //========================================================
+
+ case 14:
+ case 15:
+
+ d1=fabs( DistanceToEtaLine( pforcell, eta1) );
+ d2=fabs( DistanceToEtaLine( pforcell, eta2) );
+
+ if( d1 < s_StripWidth+s_KapGap ){
+ label18:
+ b_compartment = 8;
+ etabin = 2*i+2;
+ break;
+ }
+ if( d2 < s_StripWidth+s_KapGap ){
+ label19:
+ b_compartment = 8;
+ etabin = 2*i+3;
+ break;
+ }
+ if( r_inb > rlim1 && z_inb > zlim1 ) goto label18;
+ if( r_inb < rlim2 && z_inb > zlim2 ) goto label19;
+
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+ //======================================================
+
+ case 17:
+
+ d1=fabs( DistanceToEtaLine( pforcell, eta1) );
+
+ if( d1 < s_StripWidth+s_KapGap || z_inb > zlim1 ){
+ b_compartment = 8;
+ etabin = 2*i+2;
+ break;
+ }
+ b_compartment = 9;
+ etabin = i-1;
+ break;
+ }
+
+ //======================================================
+
+
+ // end of search for compartment and etabin
+
+ label99:
+ return validhit;
+}
+// ****************************************************************************
+G4bool EnergyCalculator::GetVolumeIndex(const G4Step *step, G4int & ModuleNumber, G4int & PhiDivNumber) const{
+ // ****************************************************************************
+
+ const G4StepPoint* pre_step_point = step->GetPreStepPoint();
+ const G4int ndepth = pre_step_point->GetTouchable()->GetHistoryDepth();
+
+ for(G4int i=0;i<=ndepth;++i){
+ G4String ivolname=pre_step_point->GetTouchable()->GetVolume(i)->GetName();
+ if( ivolname.find("BackOuterBarrette::Module") != std::string::npos ){
+ ModuleNumber = pre_step_point->GetTouchable()->GetVolume(i)->GetCopyNo();
+ }
+ if( ivolname.find("BackOuterBarrette::Module::Phidiv") != std::string::npos ){
+ PhiDivNumber = pre_step_point->GetTouchable()->GetVolume(i)->GetCopyNo();
+ }
+ }
+
+ if(!lwc()->GetisModule()){
+ if(ModuleNumber < 0 || PhiDivNumber < 0) {return false;}
+ }
+ else if(PhiDivNumber < 0 ) {return false;}
+
+ return true;
+}
+// ****************************************************************************
+G4bool EnergyCalculator::FindIdentifier_Barrett(
+ const G4Step* step,
+ G4double PhiStartOfPhiDiv,
+ std::vector<LArHitData>& hdata,
+ G4ThreeVector &startPointLocal,
+ G4ThreeVector &endPointLocal
+ ) const {
+ // ****************************************************************************
+ // works only for outer part of the full wheel or of the module in the Barrett
+ // at the back side of EMEC
+
+ // check whether we are in the outer wheel
+ if(lwc()->type() != LArG4::OuterAbsorberWheel && lwc()->type() != LArG4::OuterAbsorberModule) {
+ ATH_MSG_FATAL(" ERROR ::FindIdentifier_Barrett, not yet prepared for solidtype="
+ << LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type()));
+ }
+
+ G4bool validhit=true;
+
+ // Get point coordinates in the Atlas coord. system
+
+ const G4StepPoint* pre_step_point = step->GetPreStepPoint();
+ const G4StepPoint* post_step_point = step->GetPostStepPoint();
+
+ const G4ThreeVector startPoint = pre_step_point->GetPosition();
+ const G4ThreeVector endPoint = post_step_point->GetPosition();
+ // G4ThreeVector p = 0.5 *(startPoint+endPoint);
+ const G4ThreeVector p = startPoint; // bec. middle point maybe out of volume
+
+ // transform point to the coord system of Barrett::Module::Phidiv (alias local)
+
+ const G4AffineTransform transformation =
+ pre_step_point->GetTouchable()->GetHistory()->GetTopTransform();
+ startPointLocal = transformation.TransformPoint(startPoint);
+ endPointLocal = transformation.TransformPoint(endPoint);
+ // G4ThreeVector pinLocal = 0.5 * (startPointLocal + endPointLocal);
+ const G4ThreeVector pinLocal = startPointLocal;
+
+ //------ code transfered to getPhiStartOfPhiDiv ------//
+
+ const G4double z_inb= lwc()->GetdWRPtoFrontFace()/2.-pinLocal.z(); //dist. from front end of the Back Barrettes
+ const G4double r_inb= pinLocal.perp(); //dist from the z axis
+ const G4ThreeVector pforcell=G4ThreeVector( pinLocal.x(), pinLocal.y(),
+ lwc()->GetElecFocaltoWRP()+lwc()->GetdWRPtoFrontFace()+lwc()->GetWheelThickness()+z_inb );
+ const G4double eta_inb=pforcell.pseudoRapidity(); //eta in the system where electrodes were designed
+
+ G4int compartment,etabin;
+
+ // get compartment and etaBin
+
+ validhit=GetCompartment_Barrett(pforcell, r_inb, z_inb, eta_inb, compartment, etabin);
+
+ G4int etaBin = etabin;
+ if (!validhit) {
+ compartment = 9; // to have some 'reasonable' number
+ etaBin = 0;
+ }
+
+ const G4int c = compartment-1;
+
+ G4int sampling = s_geometry[c].sampling;
+ G4int region = s_geometry[c].region;
+ const G4int atlasside = lwc()->GetAtlasZside() * s_geometry[c].zSide;
+
+ if(lwc()->GetisModule() && atlasside < 0 ) {
+ ATH_MSG_FATAL("EnergyCalculator: TB modul should be at pos z");
+ }
+
+ // get phiBin
+ G4int phigap = GetPhiGap_Barrett(pinLocal, PhiStartOfPhiDiv); // in wheel numbering scheme
+ //int phigapwheel=phigap; //for check
+
+ if(lwc()->GetisModule()) {
+ phigap = phigap - lwc()->GetStartGapNumber() + lwc()->GetLastFan()/2; // in module numbering scheme
+
+ if(phigap < lwc()->GetFirstFan() || phigap >= lwc()->GetLastFan()){
+ if (phigap<lwc()->GetFirstFan()) phigap=lwc()->GetFirstFan();
+ if (phigap>=lwc()->GetLastFan()) phigap=lwc()->GetLastFan()-1;
+ validhit=false;
+ }
+ }
+
+ G4int phiBin = phigap / s_geometry[c].gapsPerBin;
+
+ if(atlasside < 0){
+
+ // The following formula assumes that the z<0 endcap was derived
+ // from the positive endcap by rotateY(180.*deg)
+ // 29-March-2004 ML
+
+ phiBin = (s_geometry[c].maxPhi - 1)/2 - phiBin;
+ if(phiBin < 0) phiBin += s_geometry[c].maxPhi + 1;
+ }
+
+ // checks for phiBin and etaBin
+
+ if(phiBin<0) {
+ ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, phiBin < 0");
+ phiBin=0;
+ validhit=false;
+ }
+ if(phiBin>s_geometry[c].maxPhi) {
+ ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, phiBin = " << phiBin
+ << " > geometry[" << c << "].maxPhi="<< s_geometry[c].maxPhi);
+ phiBin=s_geometry[c].maxPhi;
+ validhit=false;
+ }
+ if(etaBin < 0){
+ ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, etaBin < 0");
+ etaBin=0;
+ validhit=false;
+ }
+ if(etaBin > s_geometry[c].maxEta){
+ ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, etaBin = "
+ << etaBin << " > geometry[" << c << "].maxEta="
+ << s_geometry[c].maxEta);
+ etaBin=s_geometry[c].maxEta;
+ validhit=false;
+ }
+
+ hdata[0].id.clear();
+ hdata[0].id << 4 // LArCalorimeter
+ << 1 // LArEM
+ << atlasside
+ << sampling
+ << region
+ << etaBin
+ << phiBin;
+
+ G4double timeOfFlight = 0.5 *
+ ( pre_step_point->GetGlobalTime() +
+ post_step_point->GetGlobalTime() );
+ hdata[0].time = timeOfFlight/Units::ns - p.mag()/CLHEP::c_light/Units::ns;
+
+ return validhit;
+}
+// ****************************************************************************
+G4bool EnergyCalculator::FindDMIdentifier_Barrett(const G4Step* step, std::vector<LArHitData>& hdata) const {
+ // ****************************************************************************
+
+ // G4bool validid = false;
+ // hdata[0].id = LArG4Identifier();
+ // G4bool validid = m_supportCalculator->Process(step, LArG4::VCalibrationCalculator::kOnlyID);
+ // hdata[0].id = m_supportCalculator->identifier();
+ // return validid;
+ hdata[0].id = LArG4Identifier();
+ std::vector<G4double> tmpv;
+ return m_supportCalculator->Process(step, hdata[0].id, tmpv, LArG4::kOnlyID );
+}
+
+G4double EnergyCalculator::_AdjustedPhiOfPoint_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const {
+ // G4double phi=p.phi(); // in Module::Phidiv
+ // phi = PhiStartOfPhiDiv + phi; // in Barrette
+ // if(phi < 0. ) phi=phi+CLHEP::twopi;
+ // if(phi >= CLHEP::twopi ) phi=phi-CLHEP::twopi;
+
+ // if(lwc()->GetisModule()) phi = M_PI - phi; // in EMECMother; TB modul
+ // else phi = CLHEP::twopi - phi; // in EMECMother; full wheel
+
+ return (lwc()->GetisModule()) ?
+ CLHEP::pi - _normalizeAngle2Pi(PhiStartOfPhiDiv + p.phi()) // in EMECMother; TB modul
+ :
+ CLHEP::twopi - _normalizeAngle2Pi(PhiStartOfPhiDiv + p.phi()) // in EMECMother; full wheel
+ ;
+}
+
+//****************************************************************************
+G4int EnergyCalculator::GetPhiGap_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const {
+ // ****************************************************************************
+ const G4double phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
+ if(phi > CLHEP::twopi-lwc()->GetFanStepOnPhi()*0.5) {return 0;}
+ return (G4int( (phi+lwc()->GetFanStepOnPhi()*0.5)/lwc()->GetFanStepOnPhi()) );
+}
+// ****************************************************************************
+G4double EnergyCalculator::GetGapSize_Barrett(const G4ThreeVector& p) const {
+ // ****************************************************************************
+ const G4double r=p.perp();
+ const G4double ta1=1./ sqrt(4.*r/FanAbsThickness()*r/FanAbsThickness() - 1.);
+ const G4double ta2=1./ sqrt(4.*r/FanEleThickness()*r/FanEleThickness() - 1.);
+ const G4double phieff= lwc()->GetFanStepOnPhi()*0.5-atan(ta1)-atan(ta2);
+ return (r*phieff);
+}
+// ****************************************************************************
+G4double EnergyCalculator::distance_to_the_nearest_electrode_Barrett(
+ const G4ThreeVector &p,
+ G4double PhiStartOfPhiDiv
+ ) const {
+ // ****************************************************************************
+ const G4double phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
+
+ G4double r=p.perp();
+ G4int igap;
+ G4double elephi,dphi;
+ if (phi > CLHEP::twopi-lwc()->GetFanStepOnPhi()*0.5) {
+ dphi=phi-CLHEP::twopi;
+ } else {
+ igap=G4int( (phi+lwc()->GetFanStepOnPhi()*0.5)/lwc()->GetFanStepOnPhi() );
+ elephi=igap*lwc()->GetFanStepOnPhi();
+ dphi=phi-elephi;
+ }
+ G4double dist=r*sin(fabs(dphi))-FanEleThickness()*0.5;
+
+ return fabs(dist);
+}
+
+G4bool EnergyCalculator::GetBarrettePCE(
+ const G4ThreeVector& p, G4double PhiStartOfPhiDiv,
+ G4double &phi, G4int &compartment, G4int &eta_bin
+) const
+{
+ // get etasection
+ const G4double z_inb = lwc()->GetdWRPtoFrontFace() * 0.5 - p.z(); //dist. from front end of the Back Barrettes
+ const G4double r_inb = p.perp(); //dist from the z axis
+ const G4ThreeVector pforcell {
+ p.x(), p.y(),
+ lwc()->GetElecFocaltoWRP() + lwc()->GetdWRPtoFrontFace()
+ + lwc()->GetWheelThickness() + z_inb
+ };
+ const G4double eta_inb = pforcell.pseudoRapidity(); //eta in the system where electrodes were designed
+ G4bool validhit = GetCompartment_Barrett(
+ pforcell, r_inb, z_inb, eta_inb, compartment, eta_bin
+ );
+
+ if(!validhit) return false; // p is not in the sens. region
+
+ // get electrode number and side
+ phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
+ return true;
+}
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.cc b/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.cc
index 9630e2e196d1ea77fa48603444085be6807f56e2..56fa2f73cc816cc99c89e3fde44471c27d947f39 100644
--- a/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.cc
+++ b/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.cc
@@ -61,16 +61,11 @@
#include "LArG4Code/LArG4BirksLaw.h"
#include "LArG4Code/ILArCalibCalculatorSvc.h"
-#include "LArHV/LArHVManager.h"
-#include "LArHV/EMECHVManager.h"
-#include "LArHV/EMECHVModule.h"
-#include "LArHV/EMECHVElectrode.h"
-#include "LArHV/EMECHVDescriptor.h"
-
-//#include "EMECSupportCalibrationCalculator.h"
#include "EnergyCalculator.h"
#include "AthenaKernel/Units.h"
+#include "HVHelper.h"
+
#define MSG_VECTOR(v) "(" << v.x() << ", " << v.y() << ", " << v.z() << ")"
using namespace LArG4::EC;
@@ -85,16 +80,10 @@ const G4double EnergyCalculator::s_AverageGap = 1.3*CLHEP::mm;
G4bool EnergyCalculator::s_FieldMapsRead =false;
G4String EnergyCalculator::s_FieldMapVersion = "";
// ****************************************************************************
-const G4String EnergyCalculator::s_HVEMECMapFileName="/afs/cern.ch/atlas/offline/data/lar/emec/efield/HVEMECMap.dat";
-G4bool EnergyCalculator::s_HVMapRead=false;
-const G4double EnergyCalculator::s_HV_Etalim[s_NofEtaSection+1]=
- {1.375,1.5,1.6,1.8,2.,2.1,2.3,2.5,2.8,3.2};
const G4double EnergyCalculator::s_LArTemperature_ECC0=88.15; //K
const G4double EnergyCalculator::s_LArTemperature_ECC1=88.37;
const G4double EnergyCalculator::s_LArTemperature_ECC5=87.97;
const G4double EnergyCalculator::s_LArTemperature_av =88.16;
-G4int EnergyCalculator::s_HV_Start_phi[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide];
-G4double EnergyCalculator::s_HV_Values[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide][s_NofElectrodesOut];
const G4double EnergyCalculator::s_AverageHV=1250.;//[V]
const G4double EnergyCalculator::s_AverageEfield=0.01*s_AverageHV/s_AverageGap;//[kv/cm]
const G4double EnergyCalculator::s_AverageCurrent=1./s_AverageGap*IonReco(s_AverageEfield)
@@ -109,38 +98,6 @@ G4double EnergyCalculator::s_CHCStotal=0.;
#endif
// ****************************************************************************
-const G4double EnergyCalculator::s_ColdCorrection = 1.0036256;
-const G4double EnergyCalculator::s_LongBarThickness = 20.*CLHEP::mm;
-const G4double EnergyCalculator::s_StripWidth = 3.*CLHEP::mm/s_ColdCorrection;
-const G4double EnergyCalculator::s_KapGap = 1.*CLHEP::mm/s_ColdCorrection;
-const G4double EnergyCalculator::s_EdgeWidth = 1.*CLHEP::mm;
-const G4double EnergyCalculator::s_DistOfEndofCuFromBack = 22.77*CLHEP::mm/s_ColdCorrection;
-const G4double EnergyCalculator::s_DistOfStartofCuFromBack= 31.*CLHEP::mm; // frontface of the barrette
-const G4double EnergyCalculator::s_ZmaxOfSignal = s_DistOfStartofCuFromBack
- - s_DistOfEndofCuFromBack + s_EdgeWidth;
-G4double EnergyCalculator::s_RefzDist = 0.;
-G4bool EnergyCalculator::s_SetConstOuterBarrett = false;
-G4bool EnergyCalculator::s_SetConstInnerBarrett = false;
-
-const G4double EnergyCalculator::s_S3_Etalim[21]={
- 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95,
- 2.00, 2.05, 2.10, 2.15, 2.20, 2.25, 2.30, 2.35, 2.40, 2.45, 2.5
-};
-const G4double EnergyCalculator::s_Rmeas_outer[50]={
- 11.59, 25.22, 57.28, 71.30, 85.90, 98.94, 103.09, 116.68, 130.42, 146.27,
- 147.19, 11.59, 15., 56.91, 44.37, 15.13, 14.93, 45.87, 35.03, 15.40,
- 14.04, 39.67, 26.83, 15.64, 14.90, 30.26, 14.70, 29.09, 43.12, 34.51,
- 25.08, 11.88, 14.39, 19.54, 17.80, 12.70, 15.31, 13.96, 11.79, -99.,
- 23.57, 34.64, 55.32, 65.39, 76.34, 10.83, 94.84, 98.00, -99., -99.
-};
-const G4double EnergyCalculator::s_Zmeas_outer[2]={3.81, 7.81};
-
-G4double EnergyCalculator::s_S3_Rlim[21];
-G4double EnergyCalculator::s_rlim[50];
-G4double EnergyCalculator::s_zlim[4];
-
-// ****************************************************************************
-
static const G4String
ECorr_t_option = "type",
CHC_Map_option = "chcmap",
@@ -152,37 +109,13 @@ static G4double ziw[7]; // used as const after initialization
static G4double zsep23[22]; // used as const after initialization
// ****************************************************************************
-#if 0
-s tatic G4String extract_option(G4String &options, G4String option) // Depracated? DM, 29 Jul 2015
-// ****************************************************************************
-{
- size_t l = options.length();
- size_t p1 = options.find(option);
- if(p1 < l){
- size_t p2 = options.find("=", p1);
- p2 ++;
- if(p2 < l){
- size_t p3 = options.find(" ", p2);
- if(p3 > l) p3 = l;
- return options.substr(p2, p3 - p2);
- }
- }
- return "";
-}
-#endif
// ****************************************************************************
-static inline G4double DistanceToEtaLine(const G4ThreeVector &p, G4double eta) {
+inline G4double EnergyCalculator::DistanceToEtaLine(const G4ThreeVector &p, G4double eta) const {
G4double sinTheta = 2. * exp(-eta) / (1. + exp(2. * -eta));
G4double cosTheta = (1. - exp(2. * -eta)) / (1. + exp(2. * -eta));
return p.perp() * cosTheta - p.z() * sinTheta;
}
-
-template<typename T, typename R, typename ...Args, typename ReturnType = R (T::*)(Args...)>
-ReturnType remove_const(R (T::*Func)(Args...) const) {
- return reinterpret_cast<ReturnType>(Func);
-}
-
// ****************************************************************************
G4bool EnergyCalculator::Process(const G4Step* step, std::vector<LArHitData>& hdata) const {
// ****************************************************************************
@@ -200,6 +133,11 @@ StatusCode EnergyCalculator::finalize() {
}
delete m_lwc;
m_lwc = nullptr;
+
+ if(m_HVHelper){
+ delete m_HVHelper;
+ m_HVHelper = nullptr;
+ }
return StatusCode::SUCCESS;
}
@@ -215,10 +153,9 @@ G4bool EnergyCalculator::Process_Default(const G4Step* step, std::vector<LArHitD
G4double E = step->GetTotalEnergyDeposit() * step->GetTrack()->GetWeight();
if (m_birksLaw) {
const G4ThreeVector midpoint = (startPointinLocal + endPointinLocal) * 0.5;
- const G4double wholeStepLengthCm = step->GetStepLength() / CLHEP::cm;
- const G4double gap = (this->*m_GetGapSize_type)(midpoint); // const method
- const G4double efield = 0.01 * (this->*m_GetHV_Value_type)(midpoint, 0.123456) / gap; // estimate Efield[KV/cm]
- // 0.123456 is arbitrary number. This value is not used for this type of process
+ const G4double wholeStepLengthCm = step->GetStepLength() / Units::cm;
+ const G4double gap = GetGapSize(midpoint);
+ const G4double efield = 0.01 * m_HVHelper->GetVoltage(midpoint) / gap; // estimate Efield[KV/cm]
E = (*m_birksLaw)(E, wholeStepLengthCm, efield);
}
@@ -227,80 +164,6 @@ G4bool EnergyCalculator::Process_Default(const G4Step* step, std::vector<LArHitD
return true;
}
-G4double EnergyCalculator::getPhiStartOfPhiDiv(const G4Step* step) const {
- // --- moved from FindIdentifier_Barrett --- //
- // get Module and Phidiv number (result is put into VolumeNumber and PhiDivNumber)
- G4int ModuleNumber = -999;
- G4int PhiDivNumber = -999;
-
- G4bool validhit = GetVolumeIndex(step, ModuleNumber, PhiDivNumber);
- if (!validhit) {
- ATH_MSG_FATAL(" ERROR ::Process_Barrett::Module, Phidiv is not found"
- <<" ModuleNumber= "<<ModuleNumber<<" PhiDivNumber= "<<PhiDivNumber);
- }
-
- if(lwc()->type() == LArG4::OuterAbsorberWheel){ // for wheel calculation
- //PhiStartOfPhiDiv =
- return lwc()->GetFanStepOnPhi()/2. + ModuleNumber * CLHEP::twopi/8.
- + PhiDivNumber * lwc()->GetFanStepOnPhi();
- //phi_inb = PhiStartOfPhiDiv + pinLocal.phi(); //in ::BackOuterBarrettes
- //if(phi_inb < 0.) phi_inb = phi_inb + CLHEP::twopi;
- //if(phi_inb > CLHEP::twopi) phi_inb = phi_inb - CLHEP::twopi;
- //phi_inb = CLHEP::twopi - phi_inb; // phi in ::EmecMother system;
- }
-
- if(lwc()->type() == LArG4::OuterAbsorberModule){ // for TB modul calculation
-
- // M_PI_2 - M_PI/8. is from EMECSupportConstruction
- // PhiStartOfPhiDiv =
- return M_PI_2 - M_PI/8. + lwc()->GetFanStepOnPhi()/2 + PhiDivNumber * lwc()->GetFanStepOnPhi();
-
- //phi_inb = PhiStartOfPhiDiv + pinLocal.phi(); //in BackOuterBarrettes;
- //phi_inb = M_PI - phi_inb; // phi in ::EmecMother system;
- }
- // --- moved from FindIdentifier_Barrett --- //
- ATH_MSG_ERROR(" ERROR ::Process_Barrett::getPhiStartOfPhiDiv, LArWheelType is unknown ("<<lwc()->type()<<"), PhiStartOfPhiDiv is set zero!"
- <<" ModuleNumber= "<<ModuleNumber<<" PhiDivNumber= "<<PhiDivNumber);
- return 0.0;
-}
-
-// ****************************************************************************
-G4bool EnergyCalculator::Process_Barrett(const G4Step* step, std::vector<LArHitData>& hdata) const{
- // ****************************************************************************
- const G4double PhiStartOfPhiDiv = getPhiStartOfPhiDiv(step);
- G4ThreeVector startPointinLocal, endPointinLocal;
- if(!FindIdentifier_Barrett(step, PhiStartOfPhiDiv, hdata, startPointinLocal, endPointinLocal)){
- //cell id is not found
- if(!lwc()->GetisBarretteCalib()) return false; // ==> normal calculator return false
- else{ // if Calibration Calculator for Barrett:
- // compute DM identifier;
- // In calibration calculator mode this process is used as "geometry calculator".
- // We have to be prepared to return
- // either the emec cell id or , if the hit does not belong to any cell,
- // the DeadMatter (DM) id.
- // DM identifier for Barrett is defined
- // by the EmecSupportCalibrationCalculator.( In general it is
- // activated by the atlas_calo.py)
- const G4bool status=FindDMIdentifier_Barrett(step, hdata);
- if(!status)
- ATH_MSG_WARNING(" WARNING!!EnergyCalculator::Process_Barrett:"
- << " DM identifier not found????");
- return status;
- } // end of calibr. calculator
- } // end if cell id not found
- // compute signal in 'normal' calculator mode, if cellid is found
- G4double E = step->GetTotalEnergyDeposit() * step->GetTrack()->GetWeight();
- if (m_birksLaw){
- const G4ThreeVector midpoint = (startPointinLocal + endPointinLocal) * 0.5;
- const G4double wholeStepLengthCm = step->GetStepLength() / CLHEP::cm;
- const G4double gap = (this->*m_GetGapSize_type)(midpoint); // const method
- const G4double efield = 0.01 * (this->*m_GetHV_Value_type)(midpoint, PhiStartOfPhiDiv) / gap; // estimate Efield[KV/cm]
- E = (*m_birksLaw)(E, wholeStepLengthCm, efield);
- }
-
- hdata[0].energy = (this->*m_ecorr_method)(E, startPointinLocal, endPointinLocal, PhiStartOfPhiDiv);
- return true;
-}
// ****************************************************************************
EnergyCalculator::EnergyCalculator(const std::string& name, ISvcLocator *pSvcLocator)
: LArCalculatorSvcImp(name, pSvcLocator)
@@ -326,6 +189,7 @@ EnergyCalculator::EnergyCalculator(const std::string& name, ISvcLocator *pSvcLoc
, m_zside(1)
, m_birksLaw(nullptr)
, m_lwc(nullptr)
+ , m_HVHelper(nullptr)
// ****************************************************************************
{
ATH_MSG_DEBUG("constructor started");
@@ -493,8 +357,8 @@ StatusCode EnergyCalculator::initialize()
<< "Barrett section is not (yet) prepared for "
<< LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type()));
G4ExceptionDescription description;
- description << G4String("Process_Barrett: ") + "Barrett section is not (yet) prepared for solidtype=" +
- LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type());
+ description << "EnergyCalculator: Barrett section is not (yet) prepared for solidtype="
+ << LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type());
G4Exception("EnergyCalculator", "BarretteSectionNotPrepared", FatalException, description);
}
else{
@@ -502,13 +366,12 @@ StatusCode EnergyCalculator::initialize()
<< "Unknown solidtype:"
<< LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type()));
G4ExceptionDescription description;
- description << G4String("Process_Barrett: ") + "Unknown solidtype=" +
- LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type());
+ description << "EnergyCalculator: Unknown solidtype="
+ << LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type());
G4Exception("EnergyCalculator", "UnknownSolidType", FatalException, description);
}
}
m_Process_type = &EnergyCalculator::Process_Barrett;
- m_GetHV_Value_type = &EnergyCalculator::GetHV_Value_Barrett;
m_GetGapSize_type = &EnergyCalculator::GetGapSize_Barrett;
m_distance_to_the_nearest_electrode_type =
&EnergyCalculator::distance_to_the_nearest_electrode_Barrett;
@@ -519,7 +382,6 @@ StatusCode EnergyCalculator::initialize()
} else {
m_Process_type = &EnergyCalculator::Process_Default;
- m_GetHV_Value_type = &EnergyCalculator::GetHV_Value_Default;
m_GetGapSize_type = &EnergyCalculator::GetGapSize_Default;
m_distance_to_the_nearest_electrode_type =
&EnergyCalculator::distance_to_the_nearest_electrode_Default;
@@ -563,34 +425,10 @@ StatusCode EnergyCalculator::initialize()
ATH_MSG_FATAL("EnergyCalculator: unknown correction type "
<< m_correction_type);
G4Exception("EnergyCalculator", "UnknownCorrectionType", FatalException,
- "Process_Barrett: unknown correction type");
+ "EnergyCalculator: unknown correction type");
}
- if(m_HVMapVersion == "v00" || m_HVMapVersion == "v01") {
- m_NofPhiSections=32;
- } else {
- m_NofPhiSections=lwc()->GetNumberOfFans();
- for(G4int i=0;i<s_NofAtlasSide;++i) {
- for(G4int j=0;j<s_NofEtaSection;++j) {
- for(G4int k=0;k<s_NofElectrodeSide;++k) {
- s_HV_Start_phi[i][j][k]=0;
- }
- }
- }
- }
- m_NumberOfElectrodesInPhiSection = lwc()->GetNumberOfFans() / NofPhiSections();
-
- if(s_HVMapRead == false){
- ATH_MSG_DEBUG("EnergyCalculator: getting EMECHVMap version = "
- << m_HVMapVersion << " from jO");
- const G4String HVMapPath = "LArG4EC";
- const G4String HVMapFileName = "HVEMECMap_"+ m_HVMapVersion + ".dat";
- const G4String HVpartialPath = HVMapPath + "/" + HVMapFileName;
- const G4String HVMapLocation = PathResolver::find_file(HVpartialPath, "ATLASCALDATA");
- GetHVMap(HVMapLocation.c_str());
- }
-
- if(m_DB_HV) get_HV_map_from_DB();
+ m_HVHelper = HVHelper::CreateHelper(lwc(), m_HVMapVersion, m_DB_HV);
// if charge collection is required
if(m_correction_type == EMEC_ECOR_CHCL || m_correction_type == EMEC_ECOR_CHCL1) {
@@ -641,7 +479,7 @@ StatusCode EnergyCalculator::initialize()
<< "Charge Collection cannot be prepared for "
<< LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type()));
G4ExceptionDescription description;
- description << G4String("Process_Barrett: ") +
+ description << G4String("EnergyCalculator: ") +
"Charge Collection cannot be prepared for solidtype=" +
LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type());
G4Exception("EnergyCalculator", "IncorrectSolidType", FatalException, description);
@@ -667,7 +505,7 @@ StatusCode EnergyCalculator::initialize()
if(m_electrode_calculator == 0){
ATH_MSG_FATAL("cannot create helper electrode calculator");
G4Exception("EnergyCalculator", "NoElectrodeCalculator", FatalException,
- "Process_Barrett: cannot create helper electrode calculator");
+ "cannot create helper electrode calculator");
}
ATH_MSG_DEBUG("helper electrode calculator created ("
<< LArWheelCalculator::LArWheelCalculatorTypeString(t) << ")");
@@ -698,74 +536,6 @@ G4double EnergyCalculator::CalculateChargeCollection1(
return a_energy;
}
-
-
-// ****************************************************************************
-G4double EnergyCalculator::GapAdjustment_old(G4double a_energy,
- const G4ThreeVector& p1, const G4ThreeVector &p2, G4double /*Barret_PhiStart*/ ) const
-// ****************************************************************************
-{
- // std::cout<<"*** GapAdjustment_old is called, a_energy="<<a_energy
- // <<std::endl;
-
- const G4ThreeVector p = (p1 + p2) * .5;
- return(a_energy / pow( ((this->*m_GetGapSize_type)(p)) / s_AverageGap, 1.3)); // const method
-}
-
-// ****************************************************************************
-G4double EnergyCalculator::GapAdjustment(G4double a_energy,
- const G4ThreeVector& a_startPoint,
- const G4ThreeVector& a_endPoint,
- G4double /*Barret_PhiStart*/) const
-// ****************************************************************************
-{
- static const G4double substpsize = 0.1*CLHEP::mm;
- // std::cout<<"*** GapAdjustment is called, a_energy="<<a_energy
- // <<std::endl;
-
- const G4ThreeVector step( a_endPoint - a_startPoint );
- const G4int nofstep= int(step.mag()/substpsize)+1;
- const G4double s_energy= a_energy/nofstep;
- G4double corr_energy= 0;
- G4ThreeVector vstep;
-
- for(G4int i = 0; i < nofstep; ++ i){ // loop for substeps
- const G4double ds = (i + 0.5) / nofstep;
- vstep = a_startPoint * (1. - ds) + a_endPoint * ds;
- const G4double gap = (this->*m_GetGapSize_type)(vstep); // const method
- corr_energy += s_energy / pow((gap / s_AverageGap), GApower());
- }
- return corr_energy;
-}
-
-// ****************************************************************************
-G4double EnergyCalculator::GapAdjustment_E(G4double a_energy,
- const G4ThreeVector& a_startPoint,
- const G4ThreeVector& a_endPoint,
- G4double Barret_PhiStart) const
-// ****************************************************************************
-{
- static const G4double substpsize = 0.1*CLHEP::mm;
- //std::cout<<"*** GapAdjustment_E is called, a_energy="<<a_energy
- // <<std::endl;
-
- const G4ThreeVector step( a_endPoint - a_startPoint );
- const G4int nofstep = G4int(step.mag()/substpsize) + 1;
- const G4double s_energy = a_energy / nofstep / s_AverageCurrent;
- G4double corr_energy = 0;
-
- for(G4int i = 0; i < nofstep; i ++){
- const G4double ds = (i + 0.5) / nofstep;
- const G4ThreeVector vstep = a_startPoint * (1. - ds) + a_endPoint * ds;
- const G4double gap = (this->*m_GetGapSize_type)(vstep); // const method
- const G4double HV_value = (this->*m_GetHV_Value_type)(vstep, Barret_PhiStart);
- const G4double efield = (HV_value * CLHEP::volt) / (gap * CLHEP::mm) / (CLHEP::kilovolt / CLHEP::cm); // estimate Efield[KV/cm]
- corr_energy += s_energy / /* gap * gap / */ ( gap - CHC_Esr() )
- * IonReco(efield) * DriftVelo(s_LArTemperature_av, efield);
- }
- return corr_energy;
-}
-
// ****************************************************************************
G4double EnergyCalculator::CalculateChargeCollection(
G4double a_energy,
@@ -833,315 +603,6 @@ G4double EnergyCalculator::CalculateChargeCollection(
return current;
}
-// ****************************************************************************
-void EnergyCalculator::GetHVMap(const G4String fname){
- // ****************************************************************************
-
- ATH_MSG_INFO("HVEMECMap is to be read from file: " << fname);
-
- FILE *lun = fopen(fname, "r");
- if(lun == 0){
- ATH_MSG_ERROR("GetHVMap - file '" << fname << "' not found!"
- << std::endl
- << " Trying to read HVEMECMap from local file:"
- << "HVEMECMap.dat");
- lun = fopen("HVEMECMap.dat", "r");
- }
- if(lun == 0){
- ATH_MSG_FATAL("GetHVMap - file "<< "./HVEMECMap.dat not found!"
- << "Cannot obtain HV map.");
- G4Exception("EnergyCalculator", "NoHVMap", FatalException,
- "GetHVMap: could not read map file");
- }
-
- char ch[80], ch1[80];
-
- for(G4int i=0;i<80;i++) {
- ch[i]=0;
- ch1[i]=' ';
- }
-
- fgets(ch,80,lun);
- ATH_MSG_INFO("actual HVMapVersion = " << ch);
- ch1[0]=ch[9];
- ch1[1]=ch[10];
- if(m_HVMapVersion == "v02" || m_HVMapVersion == "v99") {
- ch1[0]=ch[10];
- ch1[1]=ch[11];
- }
- G4int iv=atoi(ch1);
- G4bool ok=false;
- if(m_HVMapVersion == "v00" && iv == 0) ok=true;
- if(m_HVMapVersion == "v01" && iv == 1) ok=true;
- if(m_HVMapVersion == "v02" && iv == 2) ok=true;
- if(m_HVMapVersion == "v99" && iv == 99) ok=true; // this is a test file
-
- if(!ok){
- ATH_MSG_FATAL("GetHVMap - file does not contain the map requested");
- G4Exception("EnergyCalculator", "IncorrectVMap", FatalException,
- "GetHVMap: incorrect map file");
- }
-
- if(m_HVMapVersion == "v00" || m_HVMapVersion == "v01" ) {
- for(G4int i = 0; i < s_NofAtlasSide; ++ i) {
- if (fscanf(lun, "%79s", ch) < 1) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- const G4String AtlasSide = ch;
- ATH_MSG_DEBUG("AtlasSide = " << AtlasSide);
- for(G4int j = 0; j < s_NofEtaSection; ++ j) {
- for(G4int k = 0; k < s_NofElectrodeSide; ++ k) {
- if (fscanf(lun, "%79s%79s", ch, ch1) < 2) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- const G4String EtaSection = ch;
- const G4String ElectrodeSide=ch1;
- ATH_MSG_DEBUG("EtaSection = " << EtaSection
- << " ElectrodeSide = " << ElectrodeSide);
- if (fscanf(lun, "%i", &s_HV_Start_phi[i][j][k]) < 1) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- ATH_MSG_DEBUG("i, j, k = " << i << ", " << j << ", " << k
- << " " <<" HV_Start_phi = " << s_HV_Start_phi[i][j][k]);
- for(G4int l = 0; l < NofPhiSections(); ++ l){
- if (fscanf(lun, "%lg", &s_HV_Values[i][j][k][l]) < 1) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- if(l == 0){
- ATH_MSG_DEBUG(" HV_Values = " << s_HV_Values[i][j][k][l]);
- }
- }
- }
- }
- } // end for(G4int i = 0; i < NofAtlasSide; ++ i)
- } // end if (HVMapVersion == "v00" || HVMapVersion == "v01" )
- if(m_HVMapVersion == "v02" || m_HVMapVersion == "v99") {
- char s[200];
- const G4int iprmx=3;
- //iprmx=3*256;
- for(G4int i=0; i<s_NofAtlasSide; ++i) { //loop for Atlas side
- for(G4int j=0; j<3; ++j) { // read header lines
- fgets(s,200,lun);
- printf("%s",s);
- }
- for(G4int j=0; j<s_NofElectrodesOut; ++j) { //loop for electrodes
- G4int electrodenumber=0;
- if (fscanf(lun,"%i",&electrodenumber) < 1) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- if (j<iprmx || j==s_NofElectrodesOut-1) {
- printf("%3i",electrodenumber);
- }
- if (j==iprmx) {
- printf("...\n");
- }
- for(G4int k=0; k<s_NofEtaSection; ++k) { //loop for etasection
- for(G4int l=0; l<s_NofElectrodeSide; ++l) { //loop for side
- if (fscanf(lun,"%lg", &s_HV_Values[i][k][l][j]) < 1) {
- ATH_MSG_ERROR("GetHVMap: Error reading map file");
- }
- if(j<iprmx || j==s_NofElectrodesOut-1) {
- printf("%8.2f",s_HV_Values[i][k][l][j]);
- }
- }
- }
- if(j<iprmx || j==s_NofElectrodesOut-1) {
- printf("\n");
- }
- }
- fgets(s,200,lun);
- }
- }
- fclose(lun);
- s_HVMapRead = true;
-}
-
-void EnergyCalculator::get_HV_map_from_DB(void) {
- ISvcLocator* svcLocator = Gaudi::svcLocator();
- StoreGateSvc* pDetStore;
-
- StatusCode status = svcLocator->service("DetectorStore", pDetStore);
- if(status.isFailure()){
- ATH_MSG_WARNING("unable to get Detector Store! Use default HV values");
- return;
- }
-
- // get EMECHV Manager
- const LArHVManager *manager = 0;
- if(pDetStore->retrieve(manager) == StatusCode::SUCCESS){
- const EMECHVManager& hvManager =
- manager->getEMECHVManager(lwc()->GetisInner()? EMECHVModule::INNER: EMECHVModule::OUTER);
- ATH_MSG_INFO("got HV Manager for " << (lwc()->GetisInner()? "inner": "outer") << " wheel");
- const EMECHVDescriptor& dsc = hvManager.getDescriptor();
- unsigned int counter = 0;
- // loop over HV modules
- for(unsigned int iSide = hvManager.beginSideIndex();
- iSide < hvManager.endSideIndex(); ++ iSide
- ){
- unsigned short jSide = 1 - iSide; // local numbering is inverse
- for(unsigned int iEta = hvManager.beginEtaIndex();
- iEta < hvManager.endEtaIndex(); ++ iEta
- ){
- unsigned int jEta = iEta;
- if(lwc()->GetisInner()) jEta += 7;
- for(unsigned int iPhi = hvManager.beginPhiIndex();
- iPhi < hvManager.endPhiIndex(); ++ iPhi
- ){
- for(unsigned int iSector = hvManager.beginSectorIndex();
- iSector < hvManager.endSectorIndex(); ++ iSector
- ){
- const EMECHVModule& hvMod = hvManager.getHVModule(iSide, iEta, iPhi, iSector);
- unsigned int nElec = hvMod.getNumElectrodes();
- for(unsigned int iElec = 0; iElec < nElec; ++ iElec){
- const EMECHVElectrode& electrode = hvMod.getElectrode(iElec);
- unsigned int jElec = iElec;
- jElec += iSector*nElec;
- jElec += iPhi*nElec*dsc.getSectorBinning().getNumDivisions();
- if(jSide == 1){
- jElec = lwc()->GetNumberOfFans() + lwc()->GetNumberOfFans() / 2 - jElec;
- if(jElec >= (unsigned int)lwc()->GetNumberOfFans()) jElec -= lwc()->GetNumberOfFans();
- }
- for(unsigned int iGap = 0; iGap < 2; ++ iGap){
- double hv = electrode.voltage(iGap);
- ATH_MSG_DEBUG("Side, Eta, Elec, Gap, hv "
- << jSide << " " << jEta << " "
- << jElec << " " << iGap << " "
- << s_HV_Values[jSide][jEta][iGap][jElec]
- << " -> " << hv);
- if(fabs((s_HV_Values[jSide][jEta][iGap][jElec] - hv)/s_HV_Values[jSide][jEta][iGap][jElec]) > 0.05){
- ATH_MSG_INFO("eta: " << dsc.getEtaBinning().binCenter(iEta) * (jSide == 0? 1: -1) << " "
- << "phi: " << dsc.getPhiBinning().binCenter(iPhi) << " "
- << "ele phi: " << electrode.getPhi()
- << " side " << iGap
- << " change HV from "
- << s_HV_Values[jSide][jEta][iGap][jElec]
- << " to " << hv);
- }
- if(hv > -999.){
- s_HV_Values[jSide][jEta][iGap][jElec] = hv;
- ++ counter;
- }
- }
- }
- }
- }
- }
- }
- ATH_MSG_INFO(counter << " HV values updated from DB");
- } else {
- ATH_MSG_WARNING("Unable to find HV Manager");
- }
-}
-// ****************************************************************************
-G4double EnergyCalculator::GetHV_Value(const G4ThreeVector& p) const
-// ****************************************************************************
-{
- // pickup HV value from the data of power supplies;
- // everything positioned in the Wheel's coord.system;
- // if it is not the same as the Atlas's one, adjustment is needed
- // either in this code or in the data file
-
- const G4int atlasside = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
-
- const G4ThreeVector pforeta(p.x(), p.y(), p.z() + lwc()->GetElecFocaltoWRP() + lwc()->GetdWRPtoFrontFace());
- const G4double mideta = pforeta.pseudoRapidity();
- G4int etasection = -1;
- for(G4int i = 1; i <= s_NofEtaSection; ++ i){
- if(mideta <= s_HV_Etalim[i]){
- etasection = i - 1;
- break;
- }
- }
- if(!(etasection>=0 && etasection <=s_NofEtaSection-1)) throw std::runtime_error("Index out of range");
-
- //assert(etasection >= 0 && etasection <= s_NofEtaSection - 1);
-
- const std::pair<G4int, G4int> gap = lwc()->GetPhiGapAndSide(p);
- G4int electrodeside = 0; //left side of electrode(small phi)
- if(gap.second > 0) electrodeside = 1;
-
- const G4int electrodenumber = lwc()->PhiGapNumberForWheel(gap.first);
- const G4int firstelectrode = s_HV_Start_phi[atlasside][etasection][electrodeside];
- G4int electrodeindex = electrodenumber-firstelectrode;
- if(electrodeindex < 0) electrodeindex += lwc()->GetNumberOfFans();
-
- const G4int phisection = electrodeindex / NumberOfElectrodesInPhiSection();
- assert(phisection>=0 && phisection<=NofPhiSections()-1);
-
- G4double HV_value= s_HV_Values[atlasside][etasection][electrodeside][phisection];
-
- ATH_MSG_DEBUG("***GetHV::="<<HV_value<<" Asde="<<atlasside<<" eta="<<etasection
- <<" Esde="<<electrodeside
- <<" fi="<<phisection<<" xyz="<<p.x()<<" "<<p.y()<<" "<<p.z()
- <<" igap.first="<<lwc()->PhiGapNumberForWheel(gap.first)
- <<" gap.second="<<gap.second);
-
- return HV_value;
-}
-
-/*
- In this type of correction energy deposited close than s_CHC_Esr to an electrod
- is suppressed.
-*/
-// ****************************************************************************
-G4double EnergyCalculator::GapAdjustment_s(G4double a_energy,
- const G4ThreeVector& a_startPoint,
- const G4ThreeVector& a_endPoint,
- G4double Barret_PhiStart) const
-// ****************************************************************************
-
-{
- static const G4double substpsize = 0.1*CLHEP::mm;
- // std::cout<<"*** GapAdjustment_s is called, a_energy="<<a_energy
- // <<std::endl;
-
- const G4ThreeVector step( a_endPoint - a_startPoint );
- const G4int nofstep= int(step.mag()/substpsize)+1;
- const G4double s_energy= a_energy/nofstep;
- G4double corr_energy= 0.;
-
- for(G4int i = 0; i < nofstep; ++ i){ // loop for substeps
- const G4double ds = (i + 0.5) / nofstep;
- const G4ThreeVector vstep = a_startPoint * (1. - ds) + a_endPoint * ds;
- G4ThreeVector tmp = vstep;
- const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(tmp, Barret_PhiStart);
- if( fabs(dte) < CHC_Esr() ) continue;
- const G4double gap = (this->*m_GetGapSize_type)(vstep); // const method
- corr_energy += s_energy / pow((gap / s_AverageGap), GApower())
- * gap / ( gap - CHC_Esr() );
- }
- return corr_energy;
-}
-// ****************************************************************************
-G4double EnergyCalculator::GapAdjustment__sE(G4double a_energy,
- const G4ThreeVector& a_startPoint,
- const G4ThreeVector& a_endPoint,
- G4double Barret_PhiStart) const
-// ****************************************************************************
-{
- static const G4double substpsize = 0.1*CLHEP::mm;
- // std::cout<<"*** GapAdjustment__sE is called, a_energy="<<a_energy
- // <<std::endl;
-
- const G4ThreeVector step( a_endPoint - a_startPoint );
- const G4int nofstep = G4int(step.mag()/substpsize) + 1;
- const G4double step_energy = a_energy / nofstep;
- G4double corr_energy = 0;
-
- for(G4int i = 0; i < nofstep; i ++){
- const G4double ds = (i + 0.5) / nofstep;
- const G4ThreeVector vstep = a_startPoint * (1. - ds) + a_endPoint * ds;
- G4ThreeVector tmp = vstep;
- const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(tmp, Barret_PhiStart);
- if( fabs(dte) < CHC_Esr() ) continue;
- const G4double gap = (this->*m_GetGapSize_type)(vstep); // const method
- const G4double efield = 0.01 * (this->*m_GetHV_Value_type)(vstep, Barret_PhiStart) / gap; // estimate Efield[KV/cm]
- corr_energy += step_energy / s_AverageCurrent / gap *
- IonReco(efield) * DriftVelo(s_LArTemperature_av, efield)
- * gap / ( gap - CHC_Esr() );
- }
- return corr_energy;
-}
// ****************************************************************************
// The static arrays that describe the endcap geometry. Mostly these
// were copied from the ATLAS G3 code.
@@ -1153,18 +614,7 @@ G4double EnergyCalculator::GapAdjustment__sE(G4double a_energy,
// Someday, the following table should be determined from the detector
// database and not hard-coded.
-typedef struct {
- G4int zSide; // +- 3 for inner wheel, +- 2 for outer wheel, z determines sign
- G4int sampling;
- G4int region;
- G4double etaScale; // 1/deta
- G4double etaOffset; // set so that the range of etaBin starts at zero for each compartment
- G4int maxEta; // the maximum value of etaBin in this compartment
- G4int gapsPerBin; // number of phi gaps (in LArWheelSolid) for each cell bin.
- G4int maxPhi; // the maximum value of phiBin in this compartment
-} geometry_t;
-
-static const geometry_t geometry[] =
+const EnergyCalculator::geometry_t EnergyCalculator::s_geometry[] =
// zSide sampling region etaScale etaOffset maxEta gapsPerBin maxPhi
{ { 3, 1, 0, 10, 25 , 6, 4, 63 },
{ 3, 2, 0, 10, 25 , 6, 4, 63 },
@@ -1271,20 +721,20 @@ G4bool EnergyCalculator::FindIdentifier_Default(
// 19-04-2007 AMS use constant m_AtlasZside obtained in constructor
// zSide is negative if z<0.
- // m_AtlasZside = geometry[c].zSide;
+ // m_AtlasZside = s_geometry[c].zSide;
// if(p.z() < 0.) m_AtlasZside = -m_AtlasZside;
// G4int atlasside = m_AtlasZside;
- const G4int atlasside = lwc()->GetAtlasZside() * geometry[c].zSide;
+ const G4int atlasside = lwc()->GetAtlasZside() * s_geometry[c].zSide;
- G4int sampling = geometry[c].sampling;
- G4int region = geometry[c].region;
+ G4int sampling = s_geometry[c].sampling;
+ G4int region = s_geometry[c].region;
- G4int etaBin = G4int(eta * geometry[c].etaScale - geometry[c].etaOffset);
+ G4int etaBin = G4int(eta * s_geometry[c].etaScale - s_geometry[c].etaOffset);
- if(etaBin < 0 || etaBin > geometry[c].maxEta) {
+ if(etaBin < 0 || etaBin > s_geometry[c].maxEta) {
validhit=false;
if (etaBin<0) etaBin=0;
- if (etaBin>geometry[c].maxEta) etaBin=geometry[c].maxEta;
+ if (etaBin>s_geometry[c].maxEta) etaBin=s_geometry[c].maxEta;
}
// ===============
@@ -1500,7 +950,7 @@ G4bool EnergyCalculator::FindIdentifier_Default(
if(compartment == 10 || compartment == 7 ){
- eta_max = (geometry[c].etaOffset+etaBin+1.)/ geometry[c].etaScale;
+ eta_max = (s_geometry[c].etaOffset+etaBin+1.)/ s_geometry[c].etaScale;
dist=DistanceToEtaLine(pforcell,eta_max);
// assert(dist >= 0.);
if (dist>=0. && dist <= DistMaxS1) cnew = c+1;
@@ -1510,8 +960,8 @@ G4bool EnergyCalculator::FindIdentifier_Default(
else if(compartment == 9 ) {
- eta_min=(geometry[c].etaOffset+etaBin) / geometry[c].etaScale;
- eta_max=(geometry[c].etaOffset+etaBin+1.)/ geometry[c].etaScale;
+ eta_min=(s_geometry[c].etaOffset+etaBin) / s_geometry[c].etaScale;
+ eta_max=(s_geometry[c].etaOffset+etaBin+1.)/ s_geometry[c].etaScale;
dist_min =DistanceToEtaLine(pforcell,eta_min);
dist_max =DistanceToEtaLine(pforcell,eta_max);
@@ -1520,7 +970,7 @@ G4bool EnergyCalculator::FindIdentifier_Default(
if(dist_min<=0. && -dist_min <= DistMax) cnew = c-1;
else {
- if(etaBin != geometry[c].maxEta) {
+ if(etaBin != s_geometry[c].maxEta) {
if(dist_max>=0. && dist_max <= DistMax) cnew = c-1;
}
else{ // close to the crack region between wheels, etaBin=maxEta;
@@ -1534,7 +984,7 @@ G4bool EnergyCalculator::FindIdentifier_Default(
else if(compartment == 2) {
- eta_min=(geometry[c].etaOffset+etaBin)/geometry[c].etaScale;
+ eta_min=(s_geometry[c].etaOffset+etaBin)/s_geometry[c].etaScale;
dist = DistanceToEtaLine(pforcell,eta_min);
// assert(dist<=0.);
@@ -1554,9 +1004,9 @@ G4bool EnergyCalculator::FindIdentifier_Default(
c=cnew;
compartment = c + 1;
- sampling = geometry[c].sampling;
- region = geometry[c].region;
- etaBin = G4int(eta * geometry[c].etaScale - geometry[c].etaOffset);
+ sampling = s_geometry[c].sampling;
+ region = s_geometry[c].region;
+ etaBin = G4int(eta * s_geometry[c].etaScale - s_geometry[c].etaOffset);
/* G4cout<<" edep in HV bus: new:comp="<<c+1<<" sampl="<<sampling<<
" eta="<<etaBin<<" reg="<<region
<<G4endl;*/
@@ -1574,28 +1024,28 @@ G4bool EnergyCalculator::FindIdentifier_Default(
}
}
- G4int phiBin = phiGap / geometry[c].gapsPerBin;
+ G4int phiBin = phiGap / s_geometry[c].gapsPerBin;
if(atlasside < 0){
// The following formula assumes that the z<0 endcap was derived
// from the positive endcap by rotateY(180.*deg)
// 29-March-2004 ML
- phiBin = (geometry[c].maxPhi - 1)/2 - phiBin;
- if(phiBin < 0) phiBin += geometry[c].maxPhi + 1;
+ phiBin = (s_geometry[c].maxPhi - 1)/2 - phiBin;
+ if(phiBin < 0) phiBin += s_geometry[c].maxPhi + 1;
}
assert(phiBin >= 0);
- assert(phiBin <= geometry[c].maxPhi);
+ assert(phiBin <= s_geometry[c].maxPhi);
if(phiBin<0) {phiBin=0;validhit=false;}
- if(phiBin>geometry[c].maxPhi) {
- phiBin=geometry[c].maxPhi;
+ if(phiBin>s_geometry[c].maxPhi) {
+ phiBin=s_geometry[c].maxPhi;
validhit=false;
}
- if(etaBin > geometry[c].maxEta){
+ if(etaBin > s_geometry[c].maxEta){
ATH_MSG_WARNING("FindIdentifier: invalid hit, etaBin = "
<< etaBin << " > geometry[" << c << "].maxEta="
- << geometry[c].maxEta);
- etaBin=geometry[c].maxEta;
+ << s_geometry[c].maxEta);
+ etaBin=s_geometry[c].maxEta;
validhit=false;
}
if(etaBin < 0){
@@ -1682,755 +1132,3 @@ double EnergyCalculator::GetGapSize(const G4ThreeVector& p) const
- lwc()->GetFanHalfThickness() + fabs(d3) - ElectrodeFanHalfThickness());
}
}
-
-// ****************************************************************************
-void EnergyCalculator::SetConst_InnerBarrett(void){
- // ****************************************************************************
- if(s_SetConstInnerBarrett) return;
- s_SetConstInnerBarrett=true;
-
- std::cout <<" ===>>> ERROR!! SetConst_InnerBarrett is called!!!" <<std::endl;
- exit(99);
-}
-
-// ****************************************************************************
-void EnergyCalculator::SetConst_OuterBarrett(void){
- // ****************************************************************************
-
- if(s_SetConstOuterBarrett) return;
- s_SetConstOuterBarrett=true;
-
- for(G4int i=0;i<=20;++i){
- const G4double teta = 2.*atan( exp(-s_S3_Etalim[i]));
- s_S3_Rlim[i] = s_RefzDist*tan(teta);
- }
-
- const G4double inv_ColdCorrection = 1. / s_ColdCorrection;
- s_rlim[0] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[0] /*11.59 */ * inv_ColdCorrection;
- s_rlim[1] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[1] /*25.22 */ * inv_ColdCorrection;
- s_rlim[2] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[2] /*57.28 */ * inv_ColdCorrection;
- s_rlim[3] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[3] /*71.30 */ * inv_ColdCorrection;
- s_rlim[4] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[4] /*85.90 */ * inv_ColdCorrection;
- s_rlim[5] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[5] /*98.94 */ * inv_ColdCorrection;
- s_rlim[6] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[6] /*103.09 */ * inv_ColdCorrection;
- s_rlim[7] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[7] /*116.68 */ * inv_ColdCorrection;
- s_rlim[8] = s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[8] /*130.42 */ * inv_ColdCorrection;
- s_rlim[9] = s_S3_Rlim[3] + s_KapGap/2. + s_Rmeas_outer[9] /*146.27 */ * inv_ColdCorrection + s_EdgeWidth;
- s_rlim[10]= s_rlim[8] + s_Rmeas_outer[10]/*147.19 */ * inv_ColdCorrection;
-
- s_rlim[11]= s_S3_Rlim[3] + s_KapGap + s_Rmeas_outer[11]/*11.59 */ * inv_ColdCorrection; //eta=1.65
- s_rlim[12]= s_S3_Rlim[3] - s_KapGap - s_Rmeas_outer[12]/*15. */ * inv_ColdCorrection;
-
- s_rlim[13]= s_S3_Rlim[4] + s_KapGap + s_Rmeas_outer[13]/*56.91 */ * inv_ColdCorrection; //eta=1.7
- s_rlim[14]= s_S3_Rlim[4] - s_KapGap - s_Rmeas_outer[14]/*44.37 */ * inv_ColdCorrection;
-
- s_rlim[15]= s_S3_Rlim[5] + s_KapGap + s_Rmeas_outer[15]/*15.13 */ * inv_ColdCorrection; //eta=1.75
- s_rlim[16]= s_S3_Rlim[5] - s_KapGap - s_Rmeas_outer[16]/*14.93 */ * inv_ColdCorrection;
-
- s_rlim[17]= s_S3_Rlim[6] + s_KapGap + s_Rmeas_outer[17]/*45.87 */ * inv_ColdCorrection; //eta=1.80
- s_rlim[18]= s_S3_Rlim[6] - s_KapGap - s_Rmeas_outer[18]/*35.03 */ * inv_ColdCorrection;
-
- s_rlim[19]= s_S3_Rlim[7] + s_KapGap + s_Rmeas_outer[19]/*15.40 */ * inv_ColdCorrection; //eta=1.85
- s_rlim[20]= s_S3_Rlim[7] - s_KapGap - s_Rmeas_outer[20]/*14.04 */ * inv_ColdCorrection;
-
- s_rlim[21]= s_S3_Rlim[8] + s_KapGap + s_Rmeas_outer[21]/*39.67 */ * inv_ColdCorrection; //eta=1.90
- s_rlim[22]= s_S3_Rlim[8] - s_KapGap - s_Rmeas_outer[22]/*26.83 */ * inv_ColdCorrection;
-
- s_rlim[23]= s_S3_Rlim[9] + s_KapGap + s_Rmeas_outer[23]/*15.64 */ * inv_ColdCorrection; //eta=1.95
- s_rlim[24]= s_S3_Rlim[9] - s_KapGap - s_Rmeas_outer[24]/*14.90 */ * inv_ColdCorrection;
-
- s_rlim[25]= s_S3_Rlim[10] + s_KapGap + s_Rmeas_outer[25]/*30.26 */ * inv_ColdCorrection; //eta=2.00
- s_rlim[26]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[26]/*14.70 */ * inv_ColdCorrection;
-
- s_rlim[27]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[27]/*29.09 */ * inv_ColdCorrection; //eta=2.05
- s_rlim[28]= s_S3_Rlim[10] - s_KapGap - s_Rmeas_outer[28]/*43.12 */ * inv_ColdCorrection; //SHAPE CHANGE!!ZZ
-
- s_rlim[29]= s_S3_Rlim[12] + s_KapGap + s_Rmeas_outer[29]/*34.51 */ * inv_ColdCorrection; //eta=2.10
- s_rlim[30]= s_S3_Rlim[12] - s_KapGap - s_Rmeas_outer[30]/*25.08 */ * inv_ColdCorrection;
-
- s_rlim[31]= s_S3_Rlim[13] + s_KapGap + s_Rmeas_outer[31]/*11.88 */ * inv_ColdCorrection; //eta=2.15
- s_rlim[32]= s_S3_Rlim[13] - s_KapGap - s_Rmeas_outer[32]/*14.39 */ * inv_ColdCorrection;
-
- s_rlim[33]= s_S3_Rlim[14] + s_KapGap + s_Rmeas_outer[33]/*19.54 */ * inv_ColdCorrection; //eta=2.20
- s_rlim[34]= s_S3_Rlim[14] - s_KapGap - s_Rmeas_outer[34]/*17.80 */ * inv_ColdCorrection; // !!ZZ
-
- s_rlim[35]= s_S3_Rlim[15] + s_KapGap + s_Rmeas_outer[35]/*12.70 */ * inv_ColdCorrection; //eta=2.25
- s_rlim[36]= s_S3_Rlim[15] - s_KapGap - s_Rmeas_outer[36]/*15.31 */ * inv_ColdCorrection;
-
- s_rlim[37]= s_S3_Rlim[16] + s_KapGap + s_Rmeas_outer[37]/*13.96 */ * inv_ColdCorrection; //eta=2.30
- s_rlim[38]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[38]/*11.79 */ * inv_ColdCorrection; // !!ZZ!!
-
- s_rlim[40]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[40]/*23.57 */ * inv_ColdCorrection;
- s_rlim[41]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[41]/*34.64 */ * inv_ColdCorrection;
- s_rlim[42]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[42]/*55.32 */ * inv_ColdCorrection;
- s_rlim[43]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[43]/*65.39 */ * inv_ColdCorrection;
- s_rlim[44]= s_S3_Rlim[16] - s_KapGap - s_Rmeas_outer[44]/*76.34 */ * inv_ColdCorrection;
- s_rlim[45]= s_rlim[44] - s_Rmeas_outer[45]/*10.83 */ * inv_ColdCorrection;
- s_rlim[46]= s_S3_Rlim[16] - s_KapGap/2.- s_Rmeas_outer[46]/*94.84 */ * inv_ColdCorrection - s_EdgeWidth;
- s_rlim[47]= s_S3_Rlim[16] - s_KapGap/2.- s_Rmeas_outer[47]/*98.00 */ * inv_ColdCorrection - s_EdgeWidth;
-
- s_zlim[0] = - s_EdgeWidth + s_Zmeas_outer[0]/*3.81*/ * inv_ColdCorrection; //rel. to the end_of_HV_Cu
- s_zlim[1] = - s_KapGap/2. + s_Zmeas_outer[1]/*7.81*/ * inv_ColdCorrection;
- s_zlim[2] = s_StripWidth + 1./2. * s_KapGap;
- s_zlim[3] =2*s_StripWidth + 3./2. * s_KapGap;
-
- for (G4int i=0; i<=3; ++i){
- s_zlim[i]= (s_ZmaxOfSignal-s_EdgeWidth) - s_zlim[i]; //rel to start of the Barrette
- if(s_zlim[i] < 0.) s_zlim[i]=0.;
- }
-
- return;
-}
-// ****************************************************************************
-G4bool EnergyCalculator::GetCompartment_Barrett(
- G4ThreeVector pforcell, G4double r_inb, G4double z_inb, G4double eta_inb,
- G4int &b_compartment, G4int &etabin) const {
- // ****************************************************************************
-
- G4double d,d1,d2,rlim1,rlim2,rlim3,zlim1,zlim2,eta1,eta2;
- G4int i;
- G4int i0 = 3;
-
- G4bool validhit=true;
- b_compartment=-99;
- etabin=-99;
-
- if(r_inb > s_rlim[10] || r_inb < s_rlim[47] )
- {validhit=false; goto label99;}
- if(z_inb > s_ZmaxOfSignal ){validhit=false; goto label99;}
-
- if(r_inb > s_rlim[0]) { // Upper corner
-
- if(r_inb > s_rlim[9]) {
- if(z_inb > s_zlim[0]) {validhit=false; goto label99;}
- label1:
- if(z_inb > s_zlim[1]) {
- b_compartment = 8;
- etabin = 4;
- goto label99;
- }
- label2:
- b_compartment = 9;
- etabin = 0;
- goto label99;
- }
- if(r_inb > s_rlim[8]) goto label1;
- if(r_inb > s_rlim[7]) goto label2;
- if(r_inb > s_rlim[6]) {
- b_compartment = 8;
- etabin = 5;
- goto label99;
- }
- if(r_inb > s_rlim[5]) {
- label3:
- b_compartment = 8;
- etabin = 6;
- goto label99;
- }
- if(r_inb > s_rlim[4]) {
- if(z_inb > s_zlim[1]) goto label3;
-
- label4:
- b_compartment = 9;
- etabin = 1;
- goto label99;
- }
-
- if(r_inb > s_rlim[3]) goto label4;
-
- if(r_inb > s_rlim[2]) {
- b_compartment = 8;
- etabin = 7;
- goto label99;
- }
-
- if(r_inb > s_rlim[1]) {
- b_compartment = 8;
- etabin = 8;
- goto label99;
- }
- if(r_inb > s_rlim[0]) {
- b_compartment = 9;
- etabin = 2;
- goto label99;
- }
- }
-
- if(r_inb < s_rlim[38]){ // lower corner
-
- if( r_inb > s_rlim[40] ) {
- b_compartment = 9;
- etabin = 16;
- goto label99;
- }
-
- if( r_inb > s_rlim[41] ) {
- label5:
- b_compartment = 8;
- etabin = 37;
- goto label99;
- }
-
- if( r_inb > s_rlim[42] ) {
-
- d=DistanceToEtaLine( pforcell,2.35);
- if(fabs(d) < s_StripWidth+s_KapGap) {
- if( d < 0.) {
-
- label6:
- b_compartment = 8;
- etabin = 38;
- goto label99;
- }
- label7:
- if( z_inb < s_zlim[3] ) goto label5;
- goto label6;
- }
-
- if( d > 0. ) goto label7;
-
- b_compartment = 9;
- etabin = 17;
- goto label99;
- }
-
- if( r_inb > s_rlim[43] ) {
- b_compartment = 8;
- etabin = 39;
- goto label99;
- }
-
- if( r_inb > s_rlim[44] ) {
- label8:
- b_compartment = 8;
- etabin = 40;
- goto label99;
- }
-
- if( r_inb > s_rlim[45] ) {
- if( z_inb < s_zlim[3] ) goto label8;
-
- label9:
- b_compartment = 9;
- etabin = 18;
- goto label99;
- }
-
- if( r_inb > s_rlim[46] ) goto label9;
-
- if( z_inb < s_ZmaxOfSignal/(s_rlim[46]-s_rlim[47])*(r_inb-s_rlim[47])) goto label9;
-
- validhit=false;
- goto label99;
- }
-
- // medium r region: s_rlim[0] > r > s_rlim[38];
- // from middle of cellno 2 to middle of cellno. 16
- //
-
- for( i=3; i<=17; ++i){ // eta= 1.65 - 2.35
- if( eta_inb < s_S3_Etalim[i] ) {
- i0=i;
- break;
- }
- }
-
- i=i0;
-
- eta1 = s_S3_Etalim[i-1];
- eta2 = s_S3_Etalim[i];
- rlim1 = s_rlim[2*i+5 - 1];
- rlim2 = s_rlim[2*i+5];
- zlim1 = s_zlim[2];
- zlim2 = s_zlim[3];
-
- if( i == 15 || i == 17) {
- zlim1 = s_zlim[3];
- zlim2 = s_zlim[2];
- }
-
- switch(i){
-
- case 3:
-
- if( fabs( DistanceToEtaLine(pforcell, eta2) ) < s_StripWidth+s_KapGap
- || z_inb > zlim1 ){
-
- b_compartment = 8;
- etabin = 2*i+3;;
- break;
- }
-
- b_compartment = 9;
- etabin = i-1;
- break;
-
- case 4:
- case 5:
- case 6:
- case 7:
- case 8:
- case 9:
- case 10:
- case 13:
- case 16:
-
- d1=fabs( DistanceToEtaLine( pforcell, eta1) );
- d2=fabs( DistanceToEtaLine( pforcell, eta2) );
-
- if( d1 < s_StripWidth+s_KapGap ){
- label11:
- b_compartment = 8;
- etabin = 2*i+2;
- break;
- }
- if( d2 < s_StripWidth+s_KapGap ){
- label12:
- b_compartment = 8;
- etabin = 2*i+3;
- break;
- }
-
- if( z_inb < zlim1 || (r_inb > rlim2 && r_inb < rlim1) ) {
-
- b_compartment = 9;
- etabin = i-1;
- break;
- }
- if( r_inb > rlim1 ) goto label11;
- if( r_inb < rlim2 ) goto label12;
- validhit=false;
- break;
- //========================================================
- case 11:
-
- rlim3 = s_rlim[28];
-
- d1=fabs( DistanceToEtaLine( pforcell, eta1) );
- d2=fabs( DistanceToEtaLine( pforcell, eta2) );
-
- if( d1 < s_StripWidth+s_KapGap ){
- label13:
- b_compartment = 8;
- etabin = 2*i+2;
- break;
- }
-
- if( r_inb > rlim1 && z_inb > zlim1 ) goto label13;
- if( r_inb > rlim2 ){
- label14:
- b_compartment = 9;
- etabin = i-1;
- break;
- }
-
- if( d2 < s_StripWidth+s_KapGap ){
- if( z_inb > zlim1 ) {
- label15:
- b_compartment = 8;
- etabin = 2*i+4;
- break;
- }
-
- label16:
- b_compartment = 8;
- etabin = 2*i+3;
- break;
- }
-
- if( z_inb < zlim2) goto label14;
- if( r_inb > rlim3) goto label16;
- if( z_inb > zlim1) goto label15;
- goto label16;
- //======================================================
- case 12:
-
- d1=fabs( DistanceToEtaLine( pforcell, eta1) );
- d2=fabs( DistanceToEtaLine( pforcell, eta2) );
-
- if( d1 < s_StripWidth+s_KapGap ){
- b_compartment = 8;
- etabin = 2*i+2;
- break;
- }
- if( d2 < s_StripWidth+s_KapGap ){
- label17:
- b_compartment = 8;
- etabin = 2*i+3;
- break;
- }
-
- if( r_inb > rlim2 || z_inb < zlim1 ){
- b_compartment = 9;
- etabin = i-1;
- break;
- }
- goto label17;
-
- //========================================================
-
- case 14:
- case 15:
-
- d1=fabs( DistanceToEtaLine( pforcell, eta1) );
- d2=fabs( DistanceToEtaLine( pforcell, eta2) );
-
- if( d1 < s_StripWidth+s_KapGap ){
- label18:
- b_compartment = 8;
- etabin = 2*i+2;
- break;
- }
- if( d2 < s_StripWidth+s_KapGap ){
- label19:
- b_compartment = 8;
- etabin = 2*i+3;
- break;
- }
- if( r_inb > rlim1 && z_inb > zlim1 ) goto label18;
- if( r_inb < rlim2 && z_inb > zlim2 ) goto label19;
-
- b_compartment = 9;
- etabin = i-1;
- break;
- //======================================================
-
- case 17:
-
- d1=fabs( DistanceToEtaLine( pforcell, eta1) );
-
- if( d1 < s_StripWidth+s_KapGap || z_inb > zlim1 ){
- b_compartment = 8;
- etabin = 2*i+2;
- break;
- }
- b_compartment = 9;
- etabin = i-1;
- break;
- }
-
- //======================================================
-
-
- // end of search for compartment and etabin
-
- label99:
- return validhit;
-}
-// ****************************************************************************
-G4bool EnergyCalculator::GetVolumeIndex(const G4Step *step, G4int & ModuleNumber, G4int & PhiDivNumber) const{
- // ****************************************************************************
-
- const G4StepPoint* pre_step_point = step->GetPreStepPoint();
- const G4int ndepth = pre_step_point->GetTouchable()->GetHistoryDepth();
-
- for(G4int i=0;i<=ndepth;++i){
- G4String ivolname=pre_step_point->GetTouchable()->GetVolume(i)->GetName();
- if( ivolname.find("BackOuterBarrette::Module") != std::string::npos ){
- ModuleNumber = pre_step_point->GetTouchable()->GetVolume(i)->GetCopyNo();
- }
- if( ivolname.find("BackOuterBarrette::Module::Phidiv") != std::string::npos ){
- PhiDivNumber = pre_step_point->GetTouchable()->GetVolume(i)->GetCopyNo();
- }
- }
-
- if(!lwc()->GetisModule()){
- if(ModuleNumber < 0 || PhiDivNumber < 0) {return false;}
- }
- else if(PhiDivNumber < 0 ) {return false;}
-
- return true;
-}
-// ****************************************************************************
-G4bool EnergyCalculator::FindIdentifier_Barrett(
- const G4Step* step,
- G4double PhiStartOfPhiDiv,
- std::vector<LArHitData>& hdata,
- G4ThreeVector &startPointLocal,
- G4ThreeVector &endPointLocal
- ) const {
- // ****************************************************************************
- // works only for outer part of the full wheel or of the module in the Barrett
- // at the back side of EMEC
-
- // check whether we are in the outer wheel
- if(lwc()->type() != LArG4::OuterAbsorberWheel && lwc()->type() != LArG4::OuterAbsorberModule) {
- ATH_MSG_FATAL(" ERROR ::FindIdentifier_Barrett, not yet prepared for solidtype="
- << LArWheelCalculator::LArWheelCalculatorTypeString(lwc()->type()));
- }
-
- G4bool validhit=true;
-
- // Get point coordinates in the Atlas coord. system
-
- const G4StepPoint* pre_step_point = step->GetPreStepPoint();
- const G4StepPoint* post_step_point = step->GetPostStepPoint();
-
- const G4ThreeVector startPoint = pre_step_point->GetPosition();
- const G4ThreeVector endPoint = post_step_point->GetPosition();
- // G4ThreeVector p = 0.5 *(startPoint+endPoint);
- const G4ThreeVector p = startPoint; // bec. middle point maybe out of volume
-
- // transform point to the coord system of Barrett::Module::Phidiv (alias local)
-
- const G4AffineTransform transformation =
- pre_step_point->GetTouchable()->GetHistory()->GetTopTransform();
- startPointLocal = transformation.TransformPoint(startPoint);
- endPointLocal = transformation.TransformPoint(endPoint);
- // G4ThreeVector pinLocal = 0.5 * (startPointLocal + endPointLocal);
- const G4ThreeVector pinLocal = startPointLocal;
-
- //------ code transfered to getPhiStartOfPhiDiv ------//
-
- const G4double z_inb= lwc()->GetdWRPtoFrontFace()/2.-pinLocal.z(); //dist. from front end of the Back Barrettes
- const G4double r_inb= pinLocal.perp(); //dist from the z axis
- const G4ThreeVector pforcell=G4ThreeVector( pinLocal.x(), pinLocal.y(),
- lwc()->GetElecFocaltoWRP()+lwc()->GetdWRPtoFrontFace()+lwc()->GetWheelThickness()+z_inb );
- const G4double eta_inb=pforcell.pseudoRapidity(); //eta in the system where electrodes were designed
-
- G4int compartment,etabin;
-
- // get compartment and etaBin
-
- validhit=GetCompartment_Barrett(pforcell, r_inb, z_inb, eta_inb, compartment, etabin);
-
- G4int etaBin = etabin;
- if (!validhit) {
- compartment = 9; // to have some 'reasonable' number
- etaBin = 0;
- }
-
- const G4int c = compartment-1;
-
- G4int sampling = geometry[c].sampling;
- G4int region = geometry[c].region;
- const G4int atlasside = lwc()->GetAtlasZside() * geometry[c].zSide;
-
- if(lwc()->GetisModule() && atlasside < 0 ) {
- ATH_MSG_FATAL("EnergyCalculator: TB modul should be at pos z");
- }
-
- // get phiBin
- G4int phigap = GetPhiGap_Barrett(pinLocal, PhiStartOfPhiDiv); // in wheel numbering scheme
- //int phigapwheel=phigap; //for check
-
- if(lwc()->GetisModule()) {
- phigap = phigap - lwc()->GetStartGapNumber() + lwc()->GetLastFan()/2; // in module numbering scheme
-
- if(phigap < lwc()->GetFirstFan() || phigap >= lwc()->GetLastFan()){
- if (phigap<lwc()->GetFirstFan()) phigap=lwc()->GetFirstFan();
- if (phigap>=lwc()->GetLastFan()) phigap=lwc()->GetLastFan()-1;
- validhit=false;
- }
- }
-
- G4int phiBin = phigap / geometry[c].gapsPerBin;
-
- if(atlasside < 0){
-
- // The following formula assumes that the z<0 endcap was derived
- // from the positive endcap by rotateY(180.*deg)
- // 29-March-2004 ML
-
- phiBin = (geometry[c].maxPhi - 1)/2 - phiBin;
- if(phiBin < 0) phiBin += geometry[c].maxPhi + 1;
- }
-
- // checks for phiBin and etaBin
-
- if(phiBin<0) {
- ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, phiBin < 0");
- phiBin=0;
- validhit=false;
- }
- if(phiBin>geometry[c].maxPhi) {
- ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, phiBin = " << phiBin
- << " > geometry[" << c << "].maxPhi="<< geometry[c].maxPhi);
- phiBin=geometry[c].maxPhi;
- validhit=false;
- }
- if(etaBin < 0){
- ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, etaBin < 0");
- etaBin=0;
- validhit=false;
- }
- if(etaBin > geometry[c].maxEta){
- ATH_MSG_WARNING("FindIdentifier_Barrett: invalid hit, etaBin = "
- << etaBin << " > geometry[" << c << "].maxEta="
- << geometry[c].maxEta);
- etaBin=geometry[c].maxEta;
- validhit=false;
- }
-
- hdata[0].id.clear();
- hdata[0].id << 4 // LArCalorimeter
- << 1 // LArEM
- << atlasside
- << sampling
- << region
- << etaBin
- << phiBin;
-
- G4double timeOfFlight = 0.5 *
- ( pre_step_point->GetGlobalTime() +
- post_step_point->GetGlobalTime() );
- hdata[0].time = timeOfFlight/Units::ns - p.mag()/CLHEP::c_light/Units::ns;
-
- return validhit;
-}
-// ****************************************************************************
-G4bool EnergyCalculator::FindDMIdentifier_Barrett(const G4Step* step, std::vector<LArHitData>& hdata) const {
- // ****************************************************************************
-
- // G4bool validid = false;
- // hdata[0].id = LArG4Identifier();
- // G4bool validid = m_supportCalculator->Process(step, LArG4::VCalibrationCalculator::kOnlyID);
- // hdata[0].id = m_supportCalculator->identifier();
- // return validid;
- hdata[0].id = LArG4Identifier();
- std::vector<G4double> tmpv;
- return m_supportCalculator->Process(step, hdata[0].id, tmpv, LArG4::kOnlyID );
-}
-
-G4double EnergyCalculator::_AdjustedPhiOfPoint_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const {
- // G4double phi=p.phi(); // in Module::Phidiv
- // phi = PhiStartOfPhiDiv + phi; // in Barrette
- // if(phi < 0. ) phi=phi+CLHEP::twopi;
- // if(phi >= CLHEP::twopi ) phi=phi-CLHEP::twopi;
-
- // if(lwc()->GetisModule()) phi = M_PI - phi; // in EMECMother; TB modul
- // else phi = CLHEP::twopi - phi; // in EMECMother; full wheel
-
- return (lwc()->GetisModule()) ?
- CLHEP::pi - _normalizeAngle2Pi(PhiStartOfPhiDiv + p.phi()) // in EMECMother; TB modul
- :
- CLHEP::twopi - _normalizeAngle2Pi(PhiStartOfPhiDiv + p.phi()) // in EMECMother; full wheel
- ;
-}
-
-//****************************************************************************
-G4int EnergyCalculator::GetPhiGap_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const {
- // ****************************************************************************
- const G4double phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
- if(phi > CLHEP::twopi-lwc()->GetFanStepOnPhi()*0.5) {return 0;}
- return (G4int( (phi+lwc()->GetFanStepOnPhi()*0.5)/lwc()->GetFanStepOnPhi()) );
-}
-// ****************************************************************************
-G4double EnergyCalculator::GetGapSize_Barrett(const G4ThreeVector& p) const {
- // ****************************************************************************
- const G4double r=p.perp();
- const G4double ta1=1./ sqrt(4.*r/FanAbsThickness()*r/FanAbsThickness() - 1.);
- const G4double ta2=1./ sqrt(4.*r/FanEleThickness()*r/FanEleThickness() - 1.);
- const G4double phieff= lwc()->GetFanStepOnPhi()*0.5-atan(ta1)-atan(ta2);
- return (r*phieff);
-}
-// ****************************************************************************
-G4double EnergyCalculator::distance_to_the_nearest_electrode_Barrett(
- const G4ThreeVector &p,
- G4double PhiStartOfPhiDiv
- ) const {
- // ****************************************************************************
- const G4double phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
-
- G4double r=p.perp();
- G4int igap;
- G4double elephi,dphi;
- if (phi > CLHEP::twopi-lwc()->GetFanStepOnPhi()*0.5) {
- dphi=phi-CLHEP::twopi;
- } else {
- igap=G4int( (phi+lwc()->GetFanStepOnPhi()*0.5)/lwc()->GetFanStepOnPhi() );
- elephi=igap*lwc()->GetFanStepOnPhi();
- dphi=phi-elephi;
- }
- G4double dist=r*sin(fabs(dphi))-FanEleThickness()*0.5;
-
- return fabs(dist);
-}
-
-// ****************************************************************************
-G4double EnergyCalculator::GetHV_Value_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const {
- // ****************************************************************************
-
- // get atlasside
- const G4int atlasside = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
-
- // get etasection
- const G4double z_inb = lwc()->GetdWRPtoFrontFace() * 0.5-p.z(); //dist. from front end of the Back Barrettes
- const G4double r_inb= p.perp(); //dist from the z axis
- const G4ThreeVector pforcell=G4ThreeVector( p.x(), p.y(),
- lwc()->GetElecFocaltoWRP()+lwc()->GetdWRPtoFrontFace()+lwc()->GetWheelThickness()+z_inb );
- const G4double eta_inb=pforcell.pseudoRapidity(); //eta in the system where electrodes were designed
-
- G4int compartment=-99; G4int etabin=-99; G4int etasection=-99;
-
- G4bool validhit=GetCompartment_Barrett(pforcell, r_inb, z_inb, eta_inb, compartment, etabin);
-
- if(!validhit) { return 0.;} // p is not in the sens. region
-
- switch (compartment) {
- case 8:
- if( etabin < 0) { validhit=false; break;}
- if( etabin <= 3) { etasection = 0; break;}
- if( etabin <= 7) { etasection = 1; break;}
- if( etabin <= 15) { etasection = 2; break;}
- if( etabin <= 23) { etasection = 3; break;}
- if( etabin <= 27) { etasection = 4; break;}
- if( etabin <= 35) { etasection = 5; break;}
- if( etabin <= 43) { etasection = 6; break;}
- validhit=false;
- break;
-
- case 9:
- if( etabin < 0) { validhit=false; break;}
- if( etabin <= 1) { etasection = 1; break;}
- if( etabin <= 5) { etasection = 2; break;}
- if( etabin <= 9) { etasection = 3; break;}
- if( etabin <= 11) { etasection = 4; break;}
- if( etabin <= 15) { etasection = 5; break;}
- if( etabin <= 19) { etasection = 6; break;}
- validhit=false;
- break;
-
- default:
- validhit=false;
- break;
- }
-
- if(!validhit) { return 0.;} // p is not in the sens. region
-
- // get electrode number and side
-
- const G4double phi = _AdjustedPhiOfPoint_Barrett(p, PhiStartOfPhiDiv);
-
- G4int igap,iside;
- if (phi > CLHEP::twopi-lwc()->GetFanStepOnPhi()*0.5) {
- igap=0;
- iside=0;
- } else {
- igap=G4int( (phi+lwc()->GetFanStepOnPhi()*0.5)/lwc()->GetFanStepOnPhi() );
- const G4double elephi=igap*lwc()->GetFanStepOnPhi();
- const G4double dphi=phi-elephi;
- if (dphi <=0.) iside=0;
- else iside=1;
- }
-
- const G4int electrodenumber = igap;
- const G4int electrodeside = iside;
-
- // prepare indices for different versions of HV Map
-
- const G4int firstelectrode = s_HV_Start_phi[atlasside][etasection][electrodeside];
-
- G4int electrodeindex = electrodenumber-firstelectrode;
- if(electrodeindex < 0) electrodeindex += lwc()->GetNumberOfFans();
-
- const G4int phisection = electrodeindex / NumberOfElectrodesInPhiSection();
-
- assert(phisection>=0 && phisection<=NofPhiSections()-1);
-
- // pick up HV value from the array
-
- G4double HV_value= s_HV_Values[atlasside][etasection][electrodeside][phisection];
-
- return HV_value;
-}
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.h b/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.h
index 2bf547b73f4e31912ee3829e3c918a0412d51b5e..0e29d9ae0f5607ef9b5b0b671bf7b65511c98231 100644
--- a/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.h
+++ b/LArCalorimeter/LArG4/LArG4EC/src/EnergyCalculator.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
// EnergyCalculator.h
@@ -42,11 +42,11 @@
class ILArCalibCalculatorSvc;
class LArG4BirksLaw;
-
namespace LArG4 {
namespace EC {
+ class HVHelper;
class EnergyCalculator : public LArCalculatorSvcImp
{
@@ -73,10 +73,7 @@ namespace LArG4 {
private:
- //G4int m_compartment; -> made local for FindIdentifier functions
-
G4bool (EnergyCalculator::*m_Process_type) (const G4Step*, std::vector<LArHitData>&) const;
- G4double (EnergyCalculator::*m_GetHV_Value_type) (const G4ThreeVector &p, G4double /*Barret_PhiStart*/) const;
G4double (EnergyCalculator::*m_GetGapSize_type) (const G4ThreeVector &p) const;
G4double (EnergyCalculator::*m_distance_to_the_nearest_electrode_type) (const G4ThreeVector &p, G4double /*Barret_PhiStart*/) const;
@@ -87,10 +84,8 @@ namespace LArG4 {
G4bool FindDMIdentifier_Barrett(const G4Step* step, std::vector<LArHitData>&) const;
G4bool GetCompartment_Barrett(G4ThreeVector,G4double,G4double,G4double,
G4int &, G4int &) const;
- G4double GetHV_Value_Default(const G4ThreeVector& p, G4double /*Barret_PhiStart*/) const {
- return GetHV_Value(p);
- }
- G4double GetHV_Value_Barrett(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const;
+
+ G4double GetHV_Value(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const;
G4double GetGapSize_Default(const G4ThreeVector &p) const {
return GetGapSize(p);
}
@@ -179,9 +174,6 @@ namespace LArG4 {
G4double m_FanEleThickness; // used as const after init
G4double m_WaveLength; // used as const after init
- G4int m_NofPhiSections; // used as const after init
- G4int m_NumberOfElectrodesInPhiSection; // used as const after init
-
inline G4double ElectrodeFanHalfThickness() const { return m_ElectrodeFanHalfThickness; };
inline G4double FanEleThicknessOld() const { return m_FanEleThicknessOld; };
inline G4double FanEleFoldRadiusOld() const { return m_FanEleFoldRadiusOld; };
@@ -189,9 +181,6 @@ namespace LArG4 {
inline G4double FanEleThickness() const { return m_FanEleThickness; };
inline G4double WaveLength() const { return m_WaveLength; };
- inline G4int NofPhiSections() const { return m_NofPhiSections; };
- inline G4int NumberOfElectrodesInPhiSection() const { return m_NumberOfElectrodesInPhiSection; };
-
//variables specific for Efield calculation
static G4bool s_FieldMapsRead; // used as const after init
@@ -239,40 +228,24 @@ namespace LArG4 {
{}
};
- //G4double PhiStartOfPhiDiv;
-
void CreateArrays(Wheel_Efield_Map &, G4int);
inline G4int Index(const Fold_Efield_Map* foldmap, G4int i, G4int j, G4int k ) const {
return foldmap->pLayer[i]+j*foldmap->NofPointsinLayer[i]+k;
};
void SetFoldArea(G4double, FoldArea & ) const;
- //HV for current calculation
-
- static G4bool s_HVMapRead; // used only for initialization
std::string m_HVMapVersion; // used only for initialization
+ G4bool m_DB_HV;
+
static const G4double s_AverageHV;
static const G4double s_AverageEfield;
static const G4double s_AverageCurrent;
- static const G4String s_HVEMECMapFileName; //{"HVEMECMap.dat"};
-
- static const G4int s_NofAtlasSide = 2;
- static const G4int s_NofEtaSection = 9;
- static const G4int s_NofElectrodeSide = 2;
- static const G4int s_NofElectrodesOut = 768;
- static const G4double s_HV_Etalim[s_NofEtaSection+1]; // = {1.375,1.5,1.6,1.8,2.,2.1,2.3,2.5,2.8,3.2};
- static G4int s_HV_Start_phi[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide]; // used as const after init
- static G4double s_HV_Values[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide][s_NofElectrodesOut]; // used as const after init
static const G4double s_LArTemperature_ECC0;//={88.15}; //K
static const G4double s_LArTemperature_ECC1;//={88.37};
static const G4double s_LArTemperature_ECC5;//={87.97};
static const G4double s_LArTemperature_av ;// ={88.16};
- void GetHVMap(const G4String); // used only for initialization
- G4double GetHV_Value(const G4ThreeVector& p) const;
- G4double GetHV_Value_ChColl_Wheel( const G4ThreeVector& , G4int , G4int) const;
-
//Efield in [kv/cm], driftvelo in [mm/microsec], Temperature in [K]
inline static G4double IonReco(const G4double Efield) {
@@ -326,8 +299,6 @@ namespace LArG4 {
#endif
private:
- // inline void vector_to_msg(G4ThreeVector &v) const;
-
/* to be developed...
std::pair<double, double>DxToFans(Hep3Vector &p);
double XDistanceToTheNeutralFibre(const Hep3Vector& P) const;
@@ -344,9 +315,6 @@ namespace LArG4 {
LArWheelCalculator *m_lwc;
const LArWheelCalculator * lwc() const { return m_lwc; }
- void get_HV_map_from_DB(void);
- G4bool m_DB_HV;
-
std::string m_suffix;
// Aug 2007 AMS, lost Aug 2008, restored May 2009
@@ -355,8 +323,6 @@ namespace LArG4 {
G4double GetCurrent1(const G4ThreeVector &, const G4ThreeVector &, G4double) const;
- G4double get_HV_value(const G4ThreeVector&, const std::pair<G4int, G4int> &) const;
-
EnergyCalculator (const EnergyCalculator&);
EnergyCalculator& operator= (const EnergyCalculator&);
@@ -382,8 +348,30 @@ namespace LArG4 {
}
G4double getPhiStartOfPhiDiv(const G4Step* step) const;
+
+ private:
+ HVHelper *m_HVHelper;
+ const static G4double s_GA_SubstepSize;
+ G4double DistanceToEtaLine(const G4ThreeVector &p, G4double eta) const;
+
+ struct geometry_t {
+ G4int zSide; // +- 3 for inner wheel, +- 2 for outer wheel, z determines sign
+ G4int sampling;
+ G4int region;
+ G4double etaScale; // 1/deta
+ G4double etaOffset; // set so that the range of etaBin starts at zero for each compartment
+ G4int maxEta; // the maximum value of etaBin in this compartment
+ G4int gapsPerBin; // number of phi gaps (in LArWheelSolid) for each cell bin.
+ G4int maxPhi; // the maximum value of phiBin in this compartment
};
+ static const geometry_t s_geometry[];
+ G4bool GetBarrettePCE(
+ const G4ThreeVector& p, G4double PhiStartOfPhiDiv,
+ G4double &phi, G4int &compartment, G4int &eta_bin
+ ) const;
+
+ };
} // namespace EC
} // namespace LArG4
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/GapAdjustments.cc b/LArCalorimeter/LArG4/LArG4EC/src/GapAdjustments.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d6fa41657a87a3f69ba5282b023aef3dd4a54c62
--- /dev/null
+++ b/LArCalorimeter/LArG4/LArG4EC/src/GapAdjustments.cc
@@ -0,0 +1,159 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+/*
+ * Here is test energy correction options implementation.
+ * We do not need to optimize it for performance.
+ */
+
+#include <cmath>
+#include <cassert>
+#include <string>
+
+#include "AthenaKernel/Units.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "G4ThreeVector.hh"
+#include "globals.hh"
+
+#include "EnergyCalculator.h"
+#include "HVHelper.h"
+
+using namespace LArG4::EC;
+namespace Units = Athena::Units;
+
+const G4double EnergyCalculator::s_GA_SubstepSize = 0.1*Units::mm;
+
+G4double EnergyCalculator::GetHV_Value(const G4ThreeVector& p, G4double PhiStartOfPhiDiv) const
+{
+ if(lwc()->GetisBarrette()){
+ G4double phi;
+ G4int compartment, eta_bin;
+ if(!GetBarrettePCE(p, PhiStartOfPhiDiv, phi, compartment, eta_bin)) return 0.;
+ return m_HVHelper->GetVoltageBarrett(phi, compartment, eta_bin);
+ } else {
+ return m_HVHelper->GetVoltage(p);
+ }
+}
+
+G4double EnergyCalculator::GapAdjustment_old(
+ G4double energy,
+ const G4ThreeVector& startPoint,
+ const G4ThreeVector& endPoint,
+ G4double /*Barret_PhiStart*/
+) const
+{
+ const G4ThreeVector p = (startPoint + endPoint) * .5;
+ return(energy / pow(((this->*m_GetGapSize_type)(p)) / s_AverageGap, 1.3));
+}
+
+G4double EnergyCalculator::GapAdjustment(
+ G4double energy,
+ const G4ThreeVector& startPoint,
+ const G4ThreeVector& endPoint,
+ G4double /*Barret_PhiStart*/
+) const
+{
+ const G4int nofstep =
+ int((endPoint - startPoint).mag() / s_GA_SubstepSize) + 1;
+ const G4double step_energy = energy / nofstep;
+
+ G4double corrected_energy = 0;
+ for(G4int i = 0; i < nofstep; ++ i){
+ const G4double ds = (i + 0.5) / nofstep;
+ const G4ThreeVector vstep = startPoint * (1. - ds) + endPoint * ds;
+ const G4double gap = (this->*m_GetGapSize_type)(vstep);
+ corrected_energy +=
+ step_energy / pow((gap / s_AverageGap), GApower());
+ }
+
+ return corrected_energy;
+}
+
+G4double EnergyCalculator::GapAdjustment_E(
+ G4double energy,
+ const G4ThreeVector& startPoint,
+ const G4ThreeVector& endPoint,
+ G4double Barret_PhiStart
+) const
+{
+ const G4int nofstep =
+ int((endPoint - startPoint).mag() / s_GA_SubstepSize) + 1;
+ const G4double step_energy = energy / nofstep / s_AverageCurrent;
+
+ G4double corrected_energy = 0;
+ for(G4int i = 0; i < nofstep; ++ i){
+ const G4double ds = (i + 0.5) / nofstep;
+ const G4ThreeVector vstep = startPoint * (1. - ds) + endPoint * ds;
+ const G4double gap = (this->*m_GetGapSize_type)(vstep);
+ const G4double HV_value = GetHV_Value(vstep, Barret_PhiStart);
+ const G4double efield = (HV_value * CLHEP::volt)
+ / (gap*Units::mm)
+ / (CLHEP::kilovolt / Units::cm); // estimate Efield[KV/cm]
+ corrected_energy += step_energy / (gap - CHC_Esr())
+ * IonReco(efield)
+ * DriftVelo(s_LArTemperature_av, efield);
+ }
+
+ return corrected_energy;
+}
+
+/*
+ In this type of correction energy deposited close than s_CHC_Esr to an electrod
+ is suppressed.
+*/
+G4double EnergyCalculator::GapAdjustment_s(
+ G4double energy,
+ const G4ThreeVector& startPoint,
+ const G4ThreeVector& endPoint,
+ G4double Barret_PhiStart
+) const
+{
+ const G4int nofstep =
+ int((endPoint - startPoint).mag() / s_GA_SubstepSize) + 1;
+ const G4double step_energy = energy / nofstep;
+
+ G4double corrected_energy = 0.;
+ for(G4int i = 0; i < nofstep; ++ i){
+ const G4double ds = (i + 0.5) / nofstep;
+ const G4ThreeVector vstep = startPoint * (1. - ds) + endPoint * ds;
+ G4ThreeVector tmp{vstep};
+ const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(tmp, Barret_PhiStart);
+ if(fabs(dte) < CHC_Esr()) continue;
+ const G4double gap = (this->*m_GetGapSize_type)(vstep);
+ corrected_energy += step_energy
+ / pow((gap / s_AverageGap), GApower())
+ * gap / (gap - CHC_Esr());
+ }
+
+ return corrected_energy;
+}
+
+G4double EnergyCalculator::GapAdjustment__sE(
+ G4double energy,
+ const G4ThreeVector& startPoint,
+ const G4ThreeVector& endPoint,
+ G4double Barret_PhiStart
+) const
+{
+ const G4int nofstep =
+ int((endPoint - startPoint).mag() / s_GA_SubstepSize) + 1;
+ const G4double step_energy = energy / nofstep;
+
+ G4double corrected_energy = 0;
+ for(G4int i = 0; i < nofstep; ++ i){
+ const G4double ds = (i + 0.5) / nofstep;
+ const G4ThreeVector vstep = startPoint * (1. - ds) + endPoint * ds;
+ G4ThreeVector tmp{vstep};
+ const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(tmp, Barret_PhiStart);
+ if(fabs(dte) < CHC_Esr()) continue;
+ const G4double gap = (this->*m_GetGapSize_type)(vstep); // const method
+ const G4double efield = 0.01 * GetHV_Value(vstep, Barret_PhiStart) / gap; // estimate Efield[KV/cm]
+ corrected_energy += step_energy / s_AverageCurrent / gap
+ * IonReco(efield)
+ * DriftVelo(s_LArTemperature_av, efield)
+ * gap / (gap - CHC_Esr());
+ }
+
+ return corrected_energy;
+}
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/GetCurrent1.cc b/LArCalorimeter/LArG4/LArG4EC/src/GetCurrent1.cc
index 819b52018bd91dc9f86e87522ae594175ba1cc58..f668a8e15b57ed66523bb0a1f62ed31f36c40b43 100644
--- a/LArCalorimeter/LArG4/LArG4EC/src/GetCurrent1.cc
+++ b/LArCalorimeter/LArG4/LArG4EC/src/GetCurrent1.cc
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
#include <cassert>
@@ -8,6 +8,7 @@
#include "globals.hh"
#include "EnergyCalculator.h"
+#include "HVHelper.h"
#include "CLHEP/Units/SystemOfUnits.h"
@@ -215,7 +216,7 @@ G4double LArG4::EC::EnergyCalculator::GetCurrent1(const G4ThreeVector &P1, cons
//std::cout << "\tvmap: (" << vmap[0] << ", " << vmap[1] << ", " << vmap[2] << ")" << std::endl;
- const G4double HV_value = get_HV_value(Pe, gap1);
+ const G4double HV_value = m_HVHelper->GetVoltage(Pe, gap1);
int Pe_fan = 0;
const G4double dte = elc()->DistanceToTheNearestFan(Pe, Pe_fan);
@@ -285,43 +286,6 @@ G4double LArG4::EC::EnergyCalculator::GetCurrent1(const G4ThreeVector &P1, cons
return current;
}
-
- G4double LArG4::EC::EnergyCalculator::get_HV_value(
- const G4ThreeVector& p, const std::pair<G4int, G4int> &gap) const
- {
- const G4int atlas_side = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
-
- G4ThreeVector p1 ( p );
- p1[2] += lwc()->GetElecFocaltoWRP() + lwc()->GetdWRPtoFrontFace();
- const G4double eta = p1.pseudoRapidity();
- G4int eta_section = -1;
- for(G4int i = 1; i <= s_NofEtaSection; ++ i){
- if(eta <= s_HV_Etalim[i]){
- eta_section = i - 1;
- break;
- }
- }
- if(!(eta_section>=0 && eta_section <=s_NofEtaSection-1)) throw std::runtime_error("Index out of range");
-
- //assert(eta_section >= 0 && eta_section < s_NofEtaSection);
-
- /*(right side of e large phi)*/ /*left side of electrode(small phi)*/
- const G4int e_side = (gap.second > 0) ? 1 : 0;
-
- const G4int first_electrode = s_HV_Start_phi[atlas_side][eta_section][e_side];
-
- if(first_electrode < 0 || first_electrode >= lwc()->GetNumberOfFans()){
- ATH_MSG_FATAL(" get_HV_value: first_electrode number is out of range");
- G4Exception("EnergyCalculator", "ElectrodeOutOfRange", FatalException,
- "get_HV_value: first_electrode number is out of range");
- }
-
- G4int phi_section = lwc()->PhiGapNumberForWheel(gap.first) - first_electrode;
- if(phi_section < 0) phi_section += lwc()->GetNumberOfFans();
-
- return s_HV_Values[atlas_side][eta_section][e_side][phi_section];
- }
-
/*
G4double LArG4::EC::EnergyCalculator::transform_z_to_fieldmap(const G4ThreeVector &p)
{
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.cc b/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9f7306435cfacfde3fa19d4ba282955f67a740cc
--- /dev/null
+++ b/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.cc
@@ -0,0 +1,451 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "LArHV/LArHVManager.h"
+#include "LArHV/EMECHVManager.h"
+#include "LArHV/EMECHVModule.h"
+#include "LArHV/EMECHVElectrode.h"
+#include "LArHV/EMECHVDescriptor.h"
+
+#include "GaudiKernel/ISvcLocator.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "PathResolver/PathResolver.h"
+
+#include "HVHelper.h"
+
+#include <iostream>
+
+using namespace LArG4::EC;
+
+const G4double HVHelper::s_EtaLimit[s_NofEtaSection + 1] = {
+ 1.375, 1.5, 1.6, 1.8, 2., 2.1, 2.3, 2.5, 2.8, 3.2
+};
+
+G4bool HVHelper::s_NeedToReadMaps = true;
+G4double HVHelper::s_Values[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide][s_NofElectrodesOut];
+
+G4int *HVHelperV00::s_StartPhi = nullptr;
+
+HVHelper::HVHelper(
+ const LArWheelCalculator *calc,
+ const G4String &version
+) : AthMessaging(Gaudi::svcLocator()->service< IMessageSvc >("MessageSvc"),
+ "EMECHVHelper")
+ , m_calculator(calc)
+ , m_WheelShift(calc->GetElecFocaltoWRP() + calc->GetdWRPtoFrontFace())
+ , m_NofPhiSections(
+ (version == "v00" || version == "v01")? 32: calc->GetNumberOfFans()
+ )
+ , m_EtaMin(calc->GetisInner()? 8: 1)
+ , m_EtaMax(calc->GetisInner()? 9: 7)
+{
+ ATH_MSG_DEBUG("asked for map version = " << version);
+}
+
+void HVHelper::AcquireMaps(const G4String &version, G4bool from_DB)
+{
+ if(s_NeedToReadMaps){
+ ReadMapFromFile(version); // always have values from file as defaults
+ s_NeedToReadMaps = false;
+ }
+ if(from_DB) GetMapFromDB(); // update them from DB if necessary
+}
+
+FILE * HVHelper::OpenFileAndCheckVersion(const G4String &version)
+{
+ const G4String mapPath = "LArG4EC";
+ const G4String mapFileName = "HVEMECMap_"+ version + ".dat";
+ const G4String partialPath = mapPath + "/" + mapFileName;
+ const G4String mapLocation =
+ PathResolver::find_file(partialPath, "ATLASCALDATA");
+ const G4String localFile = "HVEMECMap.dat";
+
+ ATH_MSG_INFO("reading maps from file: " << mapLocation);
+
+ FILE *F = fopen(mapLocation.c_str(), "r");
+ if(F == 0){
+ ATH_MSG_ERROR("cannot open EMEC HV Map file " << mapLocation);
+ ATH_MSG_ERROR("trying to read the map from local file"
+ << localFile);
+ F = fopen(localFile.c_str(), "r");
+ }
+ if(F == 0){
+ ATH_MSG_ERROR("Cannot open map file " << localFile);
+ ATH_MSG_FATAL("Cannot obtain HV maps");
+ G4Exception("LArG4::EC::HVHelper", "NoHVMap", FatalException,
+ "ReadMapFromFile: cannot open file");
+ }
+
+ const size_t buf_size = 80;
+ char buf[buf_size] = { 0 };
+ fgets(buf, buf_size, F);
+ char *v = buf + 9;
+ if(version == "v02" || version == "v99") v ++;
+ v[2] = 0;
+ ATH_MSG_INFO("actual HV Map Version = " << buf);
+ G4int iv = atoi(v);
+ G4bool version_ok = false;
+ if(version == "v00" && iv == 0) version_ok = true;
+ if(version == "v01" && iv == 1) version_ok = true;
+ if(version == "v02" && iv == 2) version_ok = true;
+ if(version == "v99" && iv == 99) version_ok = true; // this is a test file
+ if(!version_ok){
+ ATH_MSG_FATAL("field map version in the file (" << iv
+ << ") is different from expected ("
+ << version << ")");
+ G4Exception("LArG4::EC::HVHelper", "MapVersionMismatch",
+ FatalException,
+ "ReadMapFromFile: map version mismatch");
+ }
+ return F;
+}
+
+void HVHelperV00::ReadMapFromFile(const G4String &version)
+{
+ if(version != "v00" && version != "v01" ){
+ ATH_MSG_ERROR("wrong HV helper version");
+ return;
+ }
+ FILE *F = OpenFileAndCheckVersion(version);
+ for(G4int i = 0; i < s_NofAtlasSide; ++ i){
+ const size_t buf_size = 80;
+ char buf[buf_size] = { 0 };
+ if(fscanf(F, "%79s", buf) < 1){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+ ATH_MSG_DEBUG("AtlasSide = " << buf);
+ for(G4int j = 0; j < s_NofEtaSection; ++ j){
+ for(G4int k = 0; k < s_NofElectrodeSide; ++ k){
+ char buf1[buf_size] = { 0 };
+ if(fscanf(F, "%79s%79s", buf, buf1) < 2){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+ const G4String EtaSection = buf;
+ const G4String ElectrodeSide = buf1;
+ ATH_MSG_DEBUG("EtaSection = " << EtaSection
+ << " ElectrodeSide = " << ElectrodeSide);
+ if(fscanf(F, "%i", &(StartPhi(i, j, k))) < 1){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+ ATH_MSG_DEBUG("i, j, k = " << i << ", " << j << ", " << k
+ << " " <<" HV_Start_phi = " << StartPhi(i, j, k));
+ for(G4int l = 0; l < m_NofPhiSections; ++ l){
+ if(fscanf(F, "%lg", &s_Values[i][j][k][l]) < 1){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+ if(l == 0){
+ ATH_MSG_DEBUG("HV_Values = " << s_Values[i][j][k][l]);
+ }
+ } // phi section
+ } // electrode side
+ } // eta section
+ } // atlas side
+
+ fclose(F);
+}
+
+void HVHelperV02::ReadMapFromFile(const G4String &version)
+{
+ if(version != "v02" && version != "v99"){
+ ATH_MSG_ERROR("wrong HV helper version");
+ return;
+ }
+ FILE *F = OpenFileAndCheckVersion(version);
+
+ const size_t buf_size = 200;
+ char buffer[buf_size];
+
+// TODO: regularize printout when it will be necessary
+ for(G4int i = 0; i < s_NofAtlasSide; ++ i){
+ for(G4int j = 0; j < 3; ++ j){ // read header lines
+ fgets(buffer, buf_size, F);
+// printf("%s", buffer);
+ }
+ for(G4int j = 0; j < s_NofElectrodesOut; ++ j){
+ G4int electrodenumber = 0;
+ if(fscanf(F, "%i", &electrodenumber) < 1){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+/* if (j<iprmx || j==s_NofElectrodesOut-1) {
+ printf("%3i",electrodenumber);
+ }
+ if (j==iprmx) {
+ printf("...\n");
+ }*/
+ for(G4int k = 0; k < s_NofEtaSection; ++ k){
+ for(G4int l = 0; l < s_NofElectrodeSide; ++ l){
+ if(fscanf(F, "%lg", &s_Values[i][k][l][j]) < 1){
+ ATH_MSG_ERROR("ReadMapFromFile: Error reading map file");
+ }
+/* if(j <iprmx || j==s_NofElectrodesOut-1) {
+ printf("%8.2f", s_Values[i][k][l][j]);
+ }*/
+ }
+ }
+/* if(j<iprmx || j==s_NofElectrodesOut-1) {
+ printf("\n");
+ }*/
+ }
+ fgets(buffer, buf_size, F); // read leftover end of line
+ }
+
+ fclose(F);
+}
+
+void HVHelper::GetMapFromDB(void)
+{
+ ISvcLocator* svcLocator = Gaudi::svcLocator();
+ StoreGateSvc* pDetStore;
+
+ StatusCode status = svcLocator->service("DetectorStore", pDetStore);
+ if(status.isFailure()){
+ ATH_MSG_WARNING("unable to get Detector Store! Use default HV values");
+ return;
+ }
+
+ const G4bool isInner = lwc()->GetisInner();
+ // get EMECHV Manager
+ const LArHVManager *manager = 0;
+ if(pDetStore->retrieve(manager) == StatusCode::SUCCESS){
+ const EMECHVManager& hvManager = manager->getEMECHVManager(
+ isInner? EMECHVModule::INNER: EMECHVModule::OUTER
+ );
+ ATH_MSG_INFO("got LAr HV Manager for "
+ << (isInner? "inner": "outer") << " wheel");
+ const EMECHVDescriptor& dsc = hvManager.getDescriptor();
+ const unsigned int nFans = lwc()->GetNumberOfFans();
+ unsigned int counter = 0;
+ // loop over HV modules
+ for(unsigned int iSide = hvManager.beginSideIndex();
+ iSide < hvManager.endSideIndex(); ++ iSide
+ ){
+ unsigned short jSide = 1 - iSide; // local numbering is inverse
+ for(unsigned int iEta = hvManager.beginEtaIndex();
+ iEta < hvManager.endEtaIndex(); ++ iEta
+ ){
+ unsigned int jEta = iEta;
+ if(isInner) jEta += 7;
+ for(unsigned int iPhi = hvManager.beginPhiIndex();
+ iPhi < hvManager.endPhiIndex(); ++ iPhi
+ ){
+ for(unsigned int iSector = hvManager.beginSectorIndex();
+ iSector < hvManager.endSectorIndex(); ++ iSector
+ ){
+ const EMECHVModule& hvMod = hvManager.getHVModule(iSide, iEta, iPhi, iSector);
+ unsigned int nElec = hvMod.getNumElectrodes();
+ for(unsigned int iElec = 0; iElec < nElec; ++ iElec){
+ const EMECHVElectrode& electrode = hvMod.getElectrode(iElec);
+ unsigned int jElec = iElec;
+ jElec += iSector*nElec;
+ jElec += iPhi*nElec*dsc.getSectorBinning().getNumDivisions();
+ if(jSide == 1){
+ jElec = nFans + nFans / 2 - jElec;
+ if(jElec >= nFans) jElec -= nFans;
+ }
+ for(unsigned int iGap = 0; iGap < 2; ++ iGap){
+ double hv = electrode.voltage(iGap);
+ ATH_MSG_DEBUG("Side, Eta, Elec, Gap, hv "
+ << jSide << " " << jEta << " "
+ << jElec << " " << iGap << " "
+ << s_Values[jSide][jEta][iGap][jElec]
+ << " -> " << hv);
+ if(fabs((s_Values[jSide][jEta][iGap][jElec] - hv)/s_Values[jSide][jEta][iGap][jElec]) > 0.05){
+ ATH_MSG_INFO("eta: " << dsc.getEtaBinning().binCenter(iEta) * (jSide == 0? 1: -1) << " "
+ << "phi: " << dsc.getPhiBinning().binCenter(iPhi) << " "
+ << "ele phi: " << electrode.getPhi()
+ << " side " << iGap
+ << " change HV from "
+ << s_Values[jSide][jEta][iGap][jElec]
+ << " to " << hv);
+ }
+ if(hv > -999.){
+ s_Values[jSide][jEta][iGap][jElec] = hv;
+ ++ counter;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ ATH_MSG_INFO(counter << "HV values updated from DB");
+ } else {
+ ATH_MSG_WARNING("Unable to find LAr HV Manager");
+ }
+}
+
+G4double HVHelper::GetVoltage(const G4ThreeVector& p) const
+{
+ // pickup HV value from the data of power supplies;
+ // everything positioned in the Wheel's coord.system;
+ // if it is not the same as the Atlas's one, adjustment is needed
+ // either in this code or in the data file
+
+ const std::pair<G4int, G4int> gap = lwc()->GetPhiGapAndSide(p);
+ return GetVoltage(p, gap);
+}
+
+G4double HVHelper::GetVoltageBarrett(
+ G4double phi, // phi should be adjused in the calling function
+ G4int compartment,
+ G4int eta_bin
+) const
+{
+ const auto atlas_side = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
+
+ G4int eta_section = -1;
+ // just return 0. if the point is not in the sens. region
+ switch(compartment){
+ case 8:
+ if(eta_bin < 0){ return 0.; }
+ if(eta_bin <= 3){ eta_section = 0; break; }
+ if(eta_bin <= 7){ eta_section = 1; break; }
+ if(eta_bin <= 15){ eta_section = 2; break; }
+ if(eta_bin <= 23){ eta_section = 3; break; }
+ if(eta_bin <= 27){ eta_section = 4; break; }
+ if(eta_bin <= 35){ eta_section = 5; break; }
+ if(eta_bin <= 43){ eta_section = 6; break; }
+ return 0.;
+ break;
+ case 9:
+ if(eta_bin < 0){ return 0.; }
+ if(eta_bin <= 1){ eta_section = 1; break; }
+ if(eta_bin <= 5){ eta_section = 2; break; }
+ if(eta_bin <= 9){ eta_section = 3; break; }
+ if(eta_bin <= 11){ eta_section = 4; break; }
+ if(eta_bin <= 15){ eta_section = 5; break; }
+ if(eta_bin <= 19){ eta_section = 6; break; }
+ return 0.;
+ break;
+ default:
+ return 0.;
+ break;
+ }
+
+ G4int igap, electrode_side;
+ if(phi > Gaudi::Units::twopi - lwc()->GetFanStepOnPhi()*0.5){
+ igap = 0;
+ electrode_side = 0;
+ } else {
+ igap = G4int((phi + lwc()->GetFanStepOnPhi()*0.5) / lwc()->GetFanStepOnPhi());
+ const G4double elephi = igap*lwc()->GetFanStepOnPhi();
+ const G4double dphi = phi - elephi;
+ electrode_side = dphi <= 0.? 0: 1;
+ }
+
+ const G4int phi_section = GetPhiSection(
+ igap, atlas_side, eta_section, electrode_side
+ );
+
+ return s_Values[atlas_side][eta_section][electrode_side][phi_section];
+}
+
+G4double HVHelper::GetVoltage(
+ const G4ThreeVector& p, G4int phigap, G4int phihalfgap
+) const
+{
+ const auto atlas_side = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
+
+ const auto eta_section = GetEtaSection(p);
+
+ /*(right side of e large phi)*/ /*left side of electrode(small phi)*/
+ const G4int electrode_side = (phihalfgap%2 == 0 ) ? 1 : 0 ;
+
+ const G4int phi_section = GetPhiSection(
+ phigap - 1, atlas_side, eta_section, electrode_side
+ );
+
+ return s_Values[atlas_side][eta_section][electrode_side][phi_section];
+}
+
+G4int HVHelperV00::GetPhiSection(
+ G4int gap,
+ G4int atlas_side,
+ G4int eta_section,
+ G4int electrode_side
+) const
+{
+ const G4int &first_electrode =
+ StartPhi(atlas_side, eta_section, electrode_side);
+/* TODO: if necessary make this check at the initialization stage
+ if(first_electrode < 0 || first_electrode >= lwc()->GetNumberOfFans()){
+ ATH_MSG_FATAL(" get_HV_value: first_electrode number is out of range");
+ G4Exception("EnergyCalculator", "ElectrodeOutOfRange", FatalException,
+ "get_HV_value: first_electrode number is out of range");
+ }
+*/
+ G4int phi_section = gap - first_electrode;
+ if(phi_section < 0) phi_section += lwc()->GetNumberOfFans();
+ phi_section /= m_NumberOfElectrodesInPhiSection;
+
+ if(phi_section < 0 || phi_section >= m_NofPhiSections){
+ ATH_MSG_FATAL("phi section number is out of range");
+ G4Exception("LArG4::EC::HVHelper", "PhiSectionOutOfRange",
+ FatalException,
+ "Phi Section number is out of range"
+ );
+ }
+
+ return phi_section;
+}
+
+G4int HVHelperV02::GetPhiSection(
+ G4int gap, G4int, G4int, G4int
+) const
+{
+ G4int phi_section = gap;
+// TODO: check if this check is necessary
+// doesn't have it for barrettes
+ if(phi_section < 0) phi_section += lwc()->GetNumberOfFans();
+ return phi_section;
+}
+
+G4double HVHelper::GetVoltage(
+ const G4ThreeVector& p, const std::pair<G4int, G4int> &gap
+) const
+{
+ const auto atlas_side = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
+
+ const auto eta_section = GetEtaSection(p);
+
+ /*(right side of e large phi)*/ /*left side of electrode(small phi)*/
+ const G4int electrode_side = gap.second > 0? 1: 0;
+
+ const auto phi_section = GetPhiSection(
+ lwc()->PhiGapNumberForWheel(gap.first),
+ atlas_side, eta_section, electrode_side
+ );
+
+ return s_Values[atlas_side][eta_section][electrode_side][phi_section];
+}
+
+G4int HVHelper::GetEtaSection(const G4ThreeVector &p) const
+{
+ const G4ThreeVector p1(p.x(), p.y(), p.z() + m_WheelShift);
+ const G4double eta = p1.pseudoRapidity();
+
+ G4int eta_section = -1;
+ for(G4int i = m_EtaMin; i <= m_EtaMax; ++ i){
+ if(eta <= s_EtaLimit[i]){
+ eta_section = i - 1;
+ break;
+ }
+ }
+
+ if(eta_section < 0) throw std::runtime_error("Index out of range");
+
+ return eta_section;
+}
+
+HVHelper *HVHelper::CreateHelper(
+ const LArWheelCalculator *calc,
+ const G4String &version,
+ G4bool fromDB
+) {
+ if(version == "v00" || version == "v01"){
+ return new HVHelperV00(calc, version, fromDB);
+ } else {
+ return new HVHelperV02(calc, version, fromDB);
+ }
+}
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.h b/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.h
new file mode 100644
index 0000000000000000000000000000000000000000..88e493975d3d3dedf22c94766326b98cacf7a2b6
--- /dev/null
+++ b/LArCalorimeter/LArG4/LArG4EC/src/HVHelper.h
@@ -0,0 +1,131 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef LARG4_EC_HVHELPER_H
+#define LARG4_EC_HVHELPER_H
+
+#include <stdio.h>
+#include <stdexcept>
+
+#include "AthenaBaseComps/AthMessaging.h"
+#include "AthenaKernel/Units.h"
+
+#include "GeoSpecialShapes/LArWheelCalculator.h"
+
+#include "G4ThreeVector.hh"
+#include "globals.hh"
+
+namespace LArG4 { namespace EC {
+
+class HVHelper : public AthMessaging
+{
+ protected:
+ static constexpr G4int s_NofAtlasSide = 2;
+ static constexpr G4int s_NofEtaSection = 9;
+ static constexpr G4int s_NofElectrodeSide = 2;
+ static constexpr G4int s_NofElectrodesOut = 768;
+
+ static const G4double s_EtaLimit[s_NofEtaSection + 1];
+
+ static G4double s_Values[s_NofAtlasSide][s_NofEtaSection][s_NofElectrodeSide][s_NofElectrodesOut];
+ static G4bool s_NeedToReadMaps;
+
+ const LArWheelCalculator *m_calculator;
+ const LArWheelCalculator *lwc(void) const { return m_calculator; }
+
+ const G4double m_WheelShift;
+ const G4int m_NofPhiSections;
+
+ const G4int m_EtaMin, m_EtaMax;
+
+ HVHelper(
+ const LArWheelCalculator *calc,
+ const G4String &version
+ );
+ void AcquireMaps(const G4String &version, G4bool from_DB);
+ void GetMapFromDB(void);
+ virtual void ReadMapFromFile(const G4String &version) = 0;
+ FILE * OpenFileAndCheckVersion(const G4String &version);
+
+ G4int GetEtaSection(const G4ThreeVector &p) const;
+
+ virtual G4int GetPhiSection(
+ G4int, G4int, G4int, G4int
+ ) const = 0;
+
+ public:
+ G4double GetVoltageBarrett(
+ G4double phi, G4int compartment, G4int eta_bin
+ ) const;
+
+ G4double GetVoltage(const G4ThreeVector& p) const;
+
+ virtual G4double GetVoltage(const G4ThreeVector&, G4int, G4int) const;
+ virtual G4double GetVoltage(
+ const G4ThreeVector&, const std::pair<G4int, G4int> &
+ ) const;
+
+ static HVHelper *CreateHelper(
+ const LArWheelCalculator *calc,
+ const G4String &version,
+ G4bool fromDB
+ );
+
+ virtual ~HVHelper() {}
+};
+
+class HVHelperV00 : public HVHelper // This serves map versions 00 and 01
+{
+ private:
+ static G4int *s_StartPhi;
+ static G4int &StartPhi(G4int side, G4int eta, G4int ele)
+ {
+ return s_StartPhi[(side*s_NofEtaSection + eta)*s_NofElectrodeSide + ele];
+ }
+
+ const G4int m_NumberOfElectrodesInPhiSection;
+
+ void ReadMapFromFile(const G4String &version) override final;
+
+ virtual G4int GetPhiSection(
+ G4int, G4int, G4int, G4int
+ ) const override final;
+
+ public:
+ HVHelperV00(
+ const LArWheelCalculator *calc,
+ const G4String &version,
+ G4bool fromDB
+ ) : HVHelper(calc, version)
+ , m_NumberOfElectrodesInPhiSection(lwc()->GetNumberOfFans() / m_NofPhiSections)
+ {
+ if(s_StartPhi == nullptr){
+ s_StartPhi = new G4int [
+ s_NofAtlasSide * s_NofEtaSection * s_NofElectrodeSide
+ ];
+ }
+ AcquireMaps(version, fromDB);
+ }
+};
+
+class HVHelperV02 : public HVHelper // This serves map versions 02 and above
+{
+ private:
+ void ReadMapFromFile(const G4String &version) override final;
+
+ virtual G4int GetPhiSection(
+ G4int, G4int, G4int, G4int
+ ) const override final;
+
+ public:
+ HVHelperV02(
+ const LArWheelCalculator *calc,
+ const G4String &version,
+ G4bool fromDB
+ ) : HVHelper(calc, version) { AcquireMaps(version, fromDB); }
+};
+
+}} // namespaces
+
+#endif // LARG4_EC_HVHELPER_H
diff --git a/LArCalorimeter/LArG4/LArG4EC/src/LArEMECChargeCollection.cc b/LArCalorimeter/LArG4/LArG4EC/src/LArEMECChargeCollection.cc
index 729bfb163229d8a4288589eb3ebfee0668edec77..b327e8dcdd11dc8fc12ba80c344d5348cd276500 100644
--- a/LArCalorimeter/LArG4/LArG4EC/src/LArEMECChargeCollection.cc
+++ b/LArCalorimeter/LArG4/LArG4EC/src/LArEMECChargeCollection.cc
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
// LArEMECChargeCollection.cc
@@ -86,6 +86,7 @@
#include "globals.hh"
#include "EnergyCalculator.h"
+#include "HVHelper.h"
#include "CLHEP/Units/PhysicalConstants.h"
@@ -647,45 +648,6 @@ void LArG4::EC::EnergyCalculator::PrepareFieldMap(Wheel_Efield_Map* ChCollWheelT
} //end for foldtype
} //end of function
-// ****************************************************************************
-G4double LArG4::EC::EnergyCalculator::GetHV_Value_ChColl_Wheel(
- const G4ThreeVector& p, G4int phigap, G4int phihalfgap) const{
- // ****************************************************************************
- const G4int atlasside = (lwc()->GetAtlasZside() > 0) ? 0 : 1;
-
- const G4ThreeVector pforeta= G4ThreeVector(p.x(),p.y(),p.z()+lwc()->GetElecFocaltoWRP()+lwc()->GetdWRPtoFrontFace());
- const G4double eta=pforeta.pseudoRapidity();
- G4int etasection=-1;
- for(G4int i=1;i<=s_NofEtaSection;++i){
- if(eta<=s_HV_Etalim[i]) {etasection=i-1;break;}
- }
-
- if(!(etasection>=0 && etasection <=s_NofEtaSection-1)) throw std::runtime_error("Index out of range");
- //assert(etasection>=0 && etasection <=s_NofEtaSection-1);
- /*(right side of e large phi)*/ /*left side of electrode(small phi)*/
- const G4int electrodeside = (phihalfgap%2 == 0 ) ? 1 : 0 ;
-
- const G4int firstelectrode=s_HV_Start_phi[atlasside][etasection][electrodeside];
-
- if(!( firstelectrode>=0 && firstelectrode<= lwc()->GetNumberOfFans()-1)){
- ATH_MSG_FATAL(" GetCurrent:Electrode number is out of range");
- G4Exception("EnergyCalculator", "ElectrodeOutOfRange", FatalException, "GetCurrent: Electrode number is out of range");
- }
- G4int electrodeindex=(phigap-1)-firstelectrode;
- if(electrodeindex<0) electrodeindex=electrodeindex+lwc()->GetNumberOfFans();
- const G4int phisection=electrodeindex/NumberOfElectrodesInPhiSection();//24(8) for outer(inner) wheel
-
- if(!(phisection>=0 && phisection<=NofPhiSections()-1)){
- ATH_MSG_FATAL(" GetCurrent::Electrode number is out of range");
- G4Exception("EnergyCalculator", "ElectrodeOutOfRange", FatalException,"GetCurrent: Electrode number is out of range");
- }
-
- const G4double HV_value = s_HV_Values[atlasside][etasection][electrodeside][phisection];
-
- return HV_value;
-}
-
-
// **********************************************************************
void LArG4::EC::EnergyCalculator::TransFromBarrtoWheel(const G4double vb[], G4double PhiStartOfPhiDiv, G4double v[]) const {
// **********************************************************************
@@ -929,20 +891,25 @@ G4double LArG4::EC::EnergyCalculator::GetCurrent(
}
G4double HV_value;
+ G4ThreeVector p(hvpoint[0], hvpoint[1], hvpoint[2]);
if(lwc()->GetisBarrette()) {
- HV_value=GetHV_Value_Barrett(
- G4ThreeVector(hvpoint[0],hvpoint[1],hvpoint[2]), Barret_PhiStart);
+ G4double phi;
+ G4int compartment, eta_bin;
+ HV_value = GetBarrettePCE(p, Barret_PhiStart, phi, compartment, eta_bin)?
+ m_HVHelper->GetVoltageBarrett(phi, compartment, eta_bin):
+ 0.
+ ;
} else {
- HV_value=GetHV_Value_ChColl_Wheel(
- G4ThreeVector(hvpoint[0],hvpoint[1],hvpoint[2]),igap1,ihalfgap1);
+ HV_value = m_HVHelper->GetVoltage(p, igap1, ihalfgap1);
}
- const G4ThreeVector tmp = G4ThreeVector( hvpoint[0],hvpoint[1],hvpoint[2]);
- const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(tmp, Barret_PhiStart);
+// const G4ThreeVector tmp = G4ThreeVector( hvpoint[0],hvpoint[1],hvpoint[2]);
+// vector p should be const above!
+ const G4double dte = (this->*m_distance_to_the_nearest_electrode_type)(p, Barret_PhiStart);
if( fabs(dte) < CHC_Esr() ) continue; //skip point if too close to the electrode
- const G4double agap=(this->*m_GetGapSize_type)(tmp); //correction to electrode suppression not to
+ const G4double agap=(this->*m_GetGapSize_type)(p); //correction to electrode suppression not to
G4double x = agap/( agap - CHC_Esr() ); // change av. signal in the gap
if(x >=0.) supcorr=x;