Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:08

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HFShowerPMT.cc
0003 // Description: Parametrized version of HF hits
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "SimG4CMS/Calo/interface/HFShowerPMT.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/GlobalPhysicalConstants.h"
0014 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0015 #include <sstream>
0016 
0017 //#define EDM_ML_DEBUG
0018 
0019 HFShowerPMT::HFShowerPMT(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_HF = p.getParameter<edm::ParameterSet>("HFShowerPMT");
0025   pePerGeV = m_HF.getParameter<double>("PEPerGeVPMT");
0026 
0027   //Special Geometry parameters
0028   pmtR1 = hcalsimpar_->pmtRight_;
0029   pmtFib1 = hcalsimpar_->pmtFiberRight_;
0030   pmtR2 = hcalsimpar_->pmtLeft_;
0031   pmtFib2 = hcalsimpar_->pmtFiberLeft_;
0032 #ifdef EDM_ML_DEBUG
0033   edm::LogVerbatim("HFShower") << "HFShowerPMT: gets the Index matches for " << pmtR1.size() << " PMTs";
0034   for (unsigned int ii = 0; ii < pmtR1.size(); ii++) {
0035     edm::LogVerbatim("HFShower") << "HFShowerPMT: rIndexR[" << ii << "] = " << pmtR1[ii] << " fibreR[" << ii
0036                                  << "] = " << pmtFib1[ii] << " rIndexL[" << ii << "] = " << pmtR2[ii] << " fibreL["
0037                                  << ii << "] = " << pmtFib2[ii];
0038   }
0039 #endif
0040   cherenkov_ = std::make_unique<HFCherenkov>(m_HF);
0041 
0042   // Special Geometry parameters
0043   rTable = hcalConstant_->getRTableHF();
0044 #ifdef EDM_ML_DEBUG
0045   std::stringstream sss;
0046   for (unsigned int ig = 0; ig < rTable.size(); ++ig) {
0047     if (ig / 10 * 10 == ig) {
0048       sss << "\n";
0049     }
0050     sss << "  " << rTable[ig] / cm;
0051   }
0052   edm::LogVerbatim("HFShowerPMT") << "HFShowerPMT: " << rTable.size() << " rTable(cm):" << sss.str();
0053 #endif
0054 }
0055 
0056 HFShowerPMT::~HFShowerPMT() {}
0057 
0058 double HFShowerPMT::getHits(const G4Step* aStep) {
0059   indexR = indexF = -1;
0060 
0061   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0062   const G4VTouchable* touch = preStepPoint->GetTouchable();
0063   int boxNo = touch->GetReplicaNumber(2);
0064   int pmtNo = touch->GetReplicaNumber(1);
0065   if (boxNo <= 1) {
0066     indexR = pmtR1[pmtNo - 1];
0067     indexF = pmtFib1[pmtNo - 1];
0068   } else {
0069     indexR = pmtR2[pmtNo - 1];
0070     indexF = pmtFib2[pmtNo - 1];
0071   }
0072 
0073 #ifdef EDM_ML_DEBUG
0074   double edep = aStep->GetTotalEnergyDeposit();
0075   edm::LogVerbatim("HFShower") << "HFShowerPMT: Box " << boxNo << " PMT " << pmtNo << " Mapped Indices " << indexR
0076                                << ", " << indexF << " Edeposit " << edep / MeV << " MeV; PE " << edep * pePerGeV / GeV;
0077 #endif
0078 
0079   double photons = 0;
0080   if (indexR >= 0 && indexF > 0) {
0081     G4Track* aTrack = aStep->GetTrack();
0082     G4ParticleDefinition* particleDef = aTrack->GetDefinition();
0083     double stepl = aStep->GetStepLength();
0084     double beta = preStepPoint->GetBeta();
0085     G4ThreeVector pDir = aTrack->GetDynamicParticle()->GetMomentumDirection();
0086     G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(pDir);
0087     photons = cherenkov_->computeNPEinPMT(particleDef, beta, localMom.x(), localMom.y(), localMom.z(), stepl);
0088 #ifdef EDM_ML_DEBUG
0089     edm::LogVerbatim("HFShower") << "HFShowerPMT::getHits: for particle " << particleDef->GetParticleName() << " Step "
0090                                  << stepl << " Beta " << beta << " Direction " << pDir << " Local " << localMom
0091                                  << " p.e. " << photons;
0092 #endif
0093   }
0094   return photons;
0095 }
0096 
0097 double HFShowerPMT::getRadius() {
0098   double r = 0.;
0099   if (indexR >= 0 && indexR + 1 < (int)(rTable.size()))
0100     r = 0.5 * (rTable[indexR] + rTable[indexR + 1]);
0101 #ifdef EDM_ML_DEBUG
0102   else
0103     edm::LogVerbatim("HFShower") << "HFShowerPMT::getRadius: R " << indexR << " F " << indexF;
0104 #endif
0105   if (indexF == 2)
0106     r = -r;
0107 #ifdef EDM_ML_DEBUG
0108   edm::LogVerbatim("HFShower") << "HFShowerPMT: Radius (" << indexR << "/" << indexF << ") " << r;
0109 #endif
0110   return r;
0111 }