File indexing completed on 2023-03-17 11:24:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include "SimG4CMS/HcalTestBeam/interface/HcalTB02HcalNumberingScheme.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "G4SystemOfUnits.hh"
0020
0021
0022
0023
0024
0025 HcalTB02HcalNumberingScheme::HcalTB02HcalNumberingScheme()
0026 : HcalTB02NumberingScheme(), phiScale(1000000), etaScale(10000) {
0027 edm::LogVerbatim("HcalTBSim") << "Creating HcalTB02HcalNumberingScheme";
0028 }
0029
0030 HcalTB02HcalNumberingScheme::~HcalTB02HcalNumberingScheme() {
0031 #ifdef EDM_ML_DEBUG
0032 edm::LogVerbatim("HcalTBSim") << "Deleting HcalTB02HcalNumberingScheme";
0033 #endif
0034 }
0035
0036
0037
0038
0039
0040 int HcalTB02HcalNumberingScheme::getUnitID(const G4Step* aStep) const {
0041 int scintID = 0;
0042
0043 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0044 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0045 float hx = hitPoint.x();
0046 float hy = hitPoint.y();
0047 float hz = hitPoint.z();
0048 float hr = std::sqrt(pow(hx, 2) + pow(hy, 2));
0049
0050
0051
0052 if ((hr > 3. * m) && (hr < 3.830 * m))
0053 return scintID = 17;
0054 if (hr > 3.830 * m)
0055 return scintID = 18;
0056
0057
0058
0059 float hR = hitPoint.mag();
0060 float htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz / hR, float(1.)), float(-1.))));
0061 float hsintheta = sin(htheta);
0062 float hphi = (hR * hsintheta == 0. ? 0. : acos(std::max(std::min(hx / (hR * hsintheta), float(1.)), float(-1.))));
0063 float heta = (std::fabs(hsintheta) == 1. ? 0. : -std::log(std::fabs(tan(htheta / 2.))));
0064 int eta = int(heta / 0.087);
0065 int phi = int(hphi / (5. * degree));
0066
0067 G4VPhysicalVolume* thePV = preStepPoint->GetPhysicalVolume();
0068 int ilayer = ((thePV->GetCopyNo()) / 10) % 100;
0069 #ifdef EDM_ML_DEBUG
0070 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: Layer " << thePV->GetName()
0071 << " found at phi = " << phi << " eta = " << eta << " lay = " << thePV->GetCopyNo()
0072 << " " << ilayer;
0073 #endif
0074 scintID = phiScale * phi + etaScale * eta + ilayer;
0075 if (hy < 0.)
0076 scintID = -scintID;
0077
0078 return scintID;
0079 }
0080
0081 int HcalTB02HcalNumberingScheme::getlayerID(int sID) const {
0082 sID = std::abs(sID);
0083 int layerID = sID;
0084 if ((layerID != 17) && (layerID != 18))
0085 layerID = sID - int(float(sID) / float(etaScale)) * etaScale;
0086
0087 #ifdef EDM_ML_DEBUG
0088 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " layer = " << layerID;
0089 #endif
0090 return layerID;
0091 }
0092
0093 int HcalTB02HcalNumberingScheme::getphiID(int sID) const {
0094 float IDsign = 1.;
0095 if (sID < 0)
0096 IDsign = -1;
0097 sID = std::abs(sID);
0098 int phiID = int(float(sID) / float(phiScale));
0099 #ifdef EDM_ML_DEBUG
0100 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " phi = " << phiID;
0101 #endif
0102 if (IDsign > 0) {
0103 phiID += 4;
0104 } else {
0105 phiID = std::abs(phiID - 3);
0106 }
0107 return phiID;
0108 }
0109
0110 int HcalTB02HcalNumberingScheme::getetaID(int sID) const {
0111 sID = std::abs(sID);
0112 int aux = sID - int(float(sID) / float(phiScale)) * phiScale;
0113 int etaID = int(float(aux) / float(etaScale));
0114
0115 #ifdef EDM_ML_DEBUG
0116 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " eta = " << etaID;
0117 #endif
0118 return etaID;
0119 }