File indexing completed on 2024-10-04 22:55:04
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 <CLHEP/Units/SystemOfUnits.h>
0020 using CLHEP::degree;
0021 using CLHEP::m;
0022
0023
0024
0025
0026
0027
0028 HcalTB02HcalNumberingScheme::HcalTB02HcalNumberingScheme()
0029 : HcalTB02NumberingScheme(), phiScale(1000000), etaScale(10000) {
0030 edm::LogVerbatim("HcalTBSim") << "Creating HcalTB02HcalNumberingScheme";
0031 }
0032
0033 HcalTB02HcalNumberingScheme::~HcalTB02HcalNumberingScheme() {
0034 #ifdef EDM_ML_DEBUG
0035 edm::LogVerbatim("HcalTBSim") << "Deleting HcalTB02HcalNumberingScheme";
0036 #endif
0037 }
0038
0039
0040
0041
0042
0043 int HcalTB02HcalNumberingScheme::getUnitID(const G4Step* aStep) const {
0044 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0045 const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0046 float hx = hitPoint.x();
0047 float hy = hitPoint.y();
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 17;
0054 if (hr > 3.830 * m)
0055 return 18;
0056
0057
0058 int scintID = 0;
0059 float hz = hitPoint.z();
0060 float hR = hitPoint.mag();
0061 float htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz / hR, float(1.)), float(-1.))));
0062 float hsintheta = sin(htheta);
0063 float hphi = (hR * hsintheta == 0. ? 0. : acos(std::max(std::min(hx / (hR * hsintheta), float(1.)), float(-1.))));
0064 float heta = (std::fabs(hsintheta) == 1. ? 0. : -std::log(std::fabs(tan(htheta / 2.))));
0065 int eta = int(heta / 0.087);
0066 int phi = int(hphi / (5. * degree));
0067
0068 G4VPhysicalVolume* thePV = preStepPoint->GetPhysicalVolume();
0069 int ilayer = ((thePV->GetCopyNo()) / 10) % 100;
0070 #ifdef EDM_ML_DEBUG
0071 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: Layer " << thePV->GetName()
0072 << " found at phi = " << phi << " eta = " << eta << " lay = " << thePV->GetCopyNo()
0073 << " " << ilayer;
0074 #endif
0075 scintID = phiScale * phi + etaScale * eta + ilayer;
0076 if (hy < 0.)
0077 scintID = -scintID;
0078
0079 return scintID;
0080 }
0081
0082 int HcalTB02HcalNumberingScheme::getlayerID(int sID) const {
0083 sID = std::abs(sID);
0084 int layerID = sID;
0085 if ((layerID != 17) && (layerID != 18))
0086 layerID = sID - int(float(sID) / float(etaScale)) * etaScale;
0087
0088 #ifdef EDM_ML_DEBUG
0089 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " layer = " << layerID;
0090 #endif
0091 return layerID;
0092 }
0093
0094 int HcalTB02HcalNumberingScheme::getphiID(int sID) const {
0095 float IDsign = 1.;
0096 if (sID < 0)
0097 IDsign = -1;
0098 sID = std::abs(sID);
0099 int phiID = int(float(sID) / float(phiScale));
0100 #ifdef EDM_ML_DEBUG
0101 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " phi = " << phiID;
0102 #endif
0103 if (IDsign > 0) {
0104 phiID += 4;
0105 } else {
0106 phiID = std::abs(phiID - 3);
0107 }
0108 return phiID;
0109 }
0110
0111 int HcalTB02HcalNumberingScheme::getetaID(int sID) const {
0112 sID = std::abs(sID);
0113 int aux = sID - int(float(sID) / float(phiScale)) * phiScale;
0114 int etaID = int(float(aux) / float(etaScale));
0115
0116 #ifdef EDM_ML_DEBUG
0117 edm::LogVerbatim("HcalTBSim") << "HcalTB02HcalNumberingScheme:: scintID " << sID << " eta = " << etaID;
0118 #endif
0119 return etaID;
0120 }