HcalTestBeam/src/HcalTB02HcalNumberingScheme.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB02HcalNumberingScheme
0005 //
0006 // Implementation:
0007 //     Numbering scheme for hadron calorimeter in 2002 test beam
0008 //
0009 // Original Author:
0010 //         Created:  Sun 21 10:14:34 CEST 2006
0011 //
0012
0013 // system include files
0014
0015 // user include files
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 //#define EDM_ML_DEBUG
0024 //
0025 // constructors and destructor
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 // member functions
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   // Check if hit happened in first HO layer or second.
0051
0052   if ((hr > 3. * m) && (hr < 3.830 * m))
0053     return 17;
0054   if (hr > 3.830 * m)
0055     return 18;
0056
0057   // Compute the scintID in the HB.
0058   int scintID = 0;
0059   float hz = hitPoint.z();
0060   float hR = hitPoint.mag();  //sqrt( pow(hx,2)+pow(hy,2)+pow(hz,2) );
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 }