Calo/src/HFShowerFibreBundle.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFShowerFibreBundle.cc
0003 // Description: Hits in the fibre bundles
0004 ///////////////////////////////////////////////////////////////////////////////
0005
0006 #include "SimG4CMS/Calo/interface/HFShowerFibreBundle.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "G4NavigationHistory.hh"
0010 #include "G4VPhysicalVolume.hh"
0011 #include "G4Step.hh"
0012 #include "G4Track.hh"
0013 #include <CLHEP/Units/PhysicalConstants.h>
0014 #include <CLHEP/Units/SystemOfUnits.h>
0015 #include <sstream>
0016
0017 //#define EDM_ML_DEBUG
0018
0019 HFShowerFibreBundle::HFShowerFibreBundle(const std::string& name,
0020                                          const HcalDDDSimConstants* hcons,
0021                                          const HcalSimulationParameters* hps,
0022                                          edm::ParameterSet const& p)
0023     : hcalConstant_(hcons), hcalsimpar_(hps) {
0024   edm::ParameterSet m_HF1 = p.getParameter<edm::ParameterSet>("HFShowerStraightBundle");
0025   facTube = m_HF1.getParameter<double>("FactorBundle");
0026   cherenkov1_ = std::make_unique<HFCherenkov>(m_HF1);
0027   edm::ParameterSet m_HF2 = p.getParameter<edm::ParameterSet>("HFShowerConicalBundle");
0028   facCone = m_HF2.getParameter<double>("FactorBundle");
0029   cherenkov2_ = std::make_unique<HFCherenkov>(m_HF2);
0030   edm::LogVerbatim("HFShower") << "HFShowerFibreBundle intialized with factors: " << facTube
0031                                << " for the straight portion and " << facCone << " for the curved portion";
0032
0033   //Special Geometry parameters
0034   pmtR1 = hcalsimpar_->pmtRight_;
0035   pmtFib1 = hcalsimpar_->pmtFiberRight_;
0036   pmtR2 = hcalsimpar_->pmtLeft_;
0037   pmtFib2 = hcalsimpar_->pmtFiberLeft_;
0038 #ifdef EDM_ML_DEBUG
0039   edm::LogVerbatim("HFShower") << "HFShowerPMT: gets the Index matches for " << pmtR1.size() << " PMTs";
0040   for (unsigned int ii = 0; ii < pmtR1.size(); ii++) {
0041     edm::LogVerbatim("HFShower") << "HFShowerPMT: rIndexR[" << ii << "] = " << pmtR1[ii] << " fibreR[" << ii
0042                                  << "] = " << pmtFib1[ii] << " rIndexL[" << ii << "] = " << pmtR2[ii] << " fibreL["
0043                                  << ii << "] = " << pmtFib2[ii];
0044   }
0045 #endif
0046
0047   // Other Geometry parameters
0048   rTable = hcalConstant_->getRTableHF();
0049 #ifdef EDM_ML_DEBUG
0050   std::stringstream sss;
0051   for (unsigned int ig = 0; ig < rTable.size(); ig++) {
0052     if (ig / 10 * 10 == ig) {
0053       sss << "\n";
0054     }
0055     sss << "  " << rTable[ig] / CLHEP::cm;
0056   }
0057   edm::LogVerbatim("HFShower") << "HFShowerFibreBundle: " << rTable.size() << " rTable(cm):" << sss.str();
0058 #endif
0059 }
0060
0061 HFShowerFibreBundle::~HFShowerFibreBundle() {}
0062
0063 double HFShowerFibreBundle::getHits(const G4Step* aStep, bool type) {
0064   indexR = indexF = -1;
0065
0066   const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0067   const G4VTouchable* touch = preStepPoint->GetTouchable();
0068   int boxNo = touch->GetReplicaNumber(1);
0069   int pmtNo = touch->GetReplicaNumber(0);
0070   if (boxNo <= 1) {
0071     indexR = pmtR1[pmtNo - 1];
0072     indexF = pmtFib1[pmtNo - 1];
0073   } else {
0074     indexR = pmtR2[pmtNo - 1];
0075     indexF = pmtFib2[pmtNo - 1];
0076   }
0077
0078 #ifdef EDM_ML_DEBUG
0079   double edep = aStep->GetTotalEnergyDeposit();
0080   edm::LogVerbatim("HFShower") << "HFShowerFibreBundle: Box " << boxNo << " PMT " << pmtNo << " Mapped Indices "
0081                                << indexR << ", " << indexF << " Edeposit " << edep / CLHEP::MeV << " MeV";
0082 #endif
0083
0084   double photons = 0;
0085   if (indexR >= 0 && indexF > 0) {
0086     const G4Track* aTrack = aStep->GetTrack();
0087     const G4ParticleDefinition* particleDef = aTrack->GetDefinition();
0088     double stepl = aStep->GetStepLength();
0089     double beta = preStepPoint->GetBeta();
0090     G4ThreeVector pDir = aTrack->GetDynamicParticle()->GetMomentumDirection();
0091     G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(pDir);
0092     if (type) {
0093       photons =
0094           facCone * cherenkov2_->computeNPEinPMT(particleDef, beta, localMom.x(), localMom.y(), localMom.z(), stepl);
0095     } else {
0096       photons =
0097           facTube * cherenkov1_->computeNPEinPMT(particleDef, beta, localMom.x(), localMom.y(), localMom.z(), stepl);
0098     }
0099 #ifdef EDM_ML_DEBUG
0100     edm::LogVerbatim("HFShower") << "HFShowerFibreBundle::getHits: for particle " << particleDef->GetParticleName()
0101                                  << " Step " << stepl << " Beta " << beta << " Direction " << pDir << " Local "
0102                                  << localMom << " p.e. " << photons;
0103 #endif
0104   }
0105   return photons;
0106 }
0107
0108 double HFShowerFibreBundle::getRadius() {
0109   double r = 0.;
0110   if (indexR >= 0 && indexR + 1 < (int)(rTable.size()))
0111     r = 0.5 * (rTable[indexR] + rTable[indexR + 1]);
0112 #ifdef EDM_ML_DEBUG
0113   else
0114     edm::LogVerbatim("HFShower") << "HFShowerFibreBundle::getRadius: R " << indexR << " F " << indexF;
0115 #endif
0116   if (indexF == 2)
0117     r = -r;
0118 #ifdef EDM_ML_DEBUG
0119   edm::LogVerbatim("HFShower") << "HFShowerFibreBundle: Radius (" << indexR << "/" << indexF << ") " << r;
0120 #endif
0121   return r;
0122 }