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 "G4SystemOfUnits.hh"
0020 //#define EDM_ML_DEBUG
0021 //
0022 // constructors and destructor
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 // member functions
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   // Check if hit happened in first HO layer or second.
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   // Compute the scintID in the HB.
0058
0059   float hR = hitPoint.mag();  //sqrt( pow(hx,2)+pow(hy,2)+pow(hz,2) );
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 }