ShowerDevelopment/src/FastHFShowerLibrary.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: FastHFShowerLibrary.cc
0003 // Description: Shower library for Very forward hadron calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005
0006 #include "FastSimulation/ShowerDevelopment/interface/FastHFShowerLibrary.h"
0007 #include "FastSimulation/Event/interface/FSimEvent.h"
0008 #include "FastSimulation/Event/interface/FSimTrack.h"
0009 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0010 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012 #include "FWCore/Framework/interface/ESTransientHandle.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0016 #include "Geometry/Records/interface/HcalParametersRcd.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0019 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0020
0021 #include "Randomize.hh"
0022 #include <CLHEP/Units/SystemOfUnits.h>
0023 #include <CLHEP/Units/PhysicalConstants.h>
0024
0025 // Geant4 headers
0026 #include "G4ParticleDefinition.hh"
0027 #include "G4DynamicParticle.hh"
0028 #include "G4DecayPhysics.hh"
0029 #include "G4ParticleTable.hh"
0030 #include "G4ParticleTypes.hh"
0031
0032 // STL headers
0033 #include <iostream>
0034 #include <mutex>
0035 #include <vector>
0036
0037 //#define DebugLog
0038
0039 static std::once_flag initializeOnce;
0040
0041 FastHFShowerLibrary::FastHFShowerLibrary(edm::ParameterSet const& p, edm::ConsumesCollector&& iC)
0042     : fast(p), hcalDDDSimConstantsESToken_(iC.esConsumes()), hcalSimulationConstantsESToken_(iC.esConsumes()) {
0043   edm::ParameterSet m_HS = p.getParameter<edm::ParameterSet>("HFShowerLibrary");
0044   applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
0045 }
0046
0047 void const FastHFShowerLibrary::initHFShowerLibrary(const edm::EventSetup& iSetup) {
0048   edm::LogInfo("FastCalorimetry") << "initHFShowerLibrary::initialization";
0049
0050   hcalConstants = &iSetup.getData(hcalDDDSimConstantsESToken_);
0051   const HcalSimulationConstants* hsps = &iSetup.getData(hcalSimulationConstantsESToken_);
0052
0053   std::string name = "HcalHits";
0054   numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcalConstants);
0055   hfshower = std::make_unique<HFShowerLibrary>(name, hcalConstants, hsps->hcalsimpar(), fast);
0056
0057   //only one thread can be allowed to setup the G4 physics table.
0058   std::call_once(initializeOnce, []() {
0059     // Geant4 particles
0060     G4DecayPhysics decays;
0061     decays.ConstructParticle();
0062     G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0063     partTable->SetReadiness();
0064   });
0065   //G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0066   //hfshower->initRun(partTable, hcalConstants);  // init particle code
0067 }
0068
0069 void FastHFShowerLibrary::SetRandom(const RandomEngineAndDistribution* rnd) {
0070   // define Geant4 engine per thread
0071   G4Random::setTheEngine(&(rnd->theEngine()));
0072   LogDebug("FastHFShowerLibrary::recoHFShowerLibrary")
0073       << "Begin of event " << G4UniformRand() << "  " << rnd->theEngine().name() << "  " << rnd->theEngine();
0074 }
0075
0076 void FastHFShowerLibrary::recoHFShowerLibrary(const FSimTrack& myTrack) {
0077 #ifdef DebugLog
0078   edm::LogInfo("FastCalorimetry") << "FastHFShowerLibrary: recoHFShowerLibrary ";
0079 #endif
0080
0081   if (!myTrack.onVFcal()) {
0082 #ifdef DebugLog
0083     edm::LogInfo("FastCalorimetry") << "FastHFShowerLibrary: we should not be here ";
0084 #endif
0085   }
0086
0087   hitMap.clear();
0088   double eGen = 1000. * myTrack.vfcalEntrance().e();  // energy in [MeV]
0089   double delZv = (myTrack.vfcalEntrance().vertex().Z() > 0.0) ? 50.0 : -50.0;
0090   G4ThreeVector vertex(10. * myTrack.vfcalEntrance().vertex().X(),
0091                        10. * myTrack.vfcalEntrance().vertex().Y(),
0092                        10. * myTrack.vfcalEntrance().vertex().Z() + delZv);  // in [mm]
0093
0094   G4ThreeVector direction(
0095       myTrack.vfcalEntrance().Vect().X(), myTrack.vfcalEntrance().Vect().Y(), myTrack.vfcalEntrance().Vect().Z());
0096
0097   bool ok;
0098   double weight = 1.0;  // rad. damage
0099   int parCode = myTrack.type();
0100   double tSlice = 0.1 * vertex.mag() / 29.98;
0101
0102   std::vector<HFShowerLibrary::Hit> hits =
0103       hfshower->fillHits(vertex, direction, parCode, eGen, ok, weight, tSlice, false);
0104
0105   for (unsigned int i = 0; i < hits.size(); ++i) {
0106     G4ThreeVector pos = hits[i].position;
0107     int depth = hits[i].depth;
0108     double time = hits[i].time;
0109     if (!applyFidCut || (HFFibreFiducial::PMTNumber(pos) > 0)) {
0110       //    if (!applyFidCut || (applyFidCut && HFFibreFiducial::PMTNumber(pos)>0)) {
0111       int det = 5;
0112       int lay = 1;
0113       uint32_t id = 0;
0114       HcalNumberingFromDDD::HcalID tmp =
0115           numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
0116       modifyDepth(tmp);
0117       id = numberingScheme.getUnitID(tmp);
0118
0119       CaloHitID current_id(id, time, myTrack.id());
0120       std::map<CaloHitID, float>::iterator cellitr;
0121       cellitr = hitMap.find(current_id);
0122       if (cellitr == hitMap.end()) {
0123         hitMap.insert(std::pair<CaloHitID, float>(current_id, 1.0));
0124       } else {
0125         cellitr->second += 1.0;
0126       }
0127     }  // end of isItinFidVolume check
0128   }  // end loop over hits
0129 }
0130
0131 void FastHFShowerLibrary::modifyDepth(HcalNumberingFromDDD::HcalID& id) {
0132   if (id.subdet == HcalForward) {
0133     int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
0134     if (hcalConstants->maxHFDepth(ieta, id.phis) > 2) {
0135       if (id.depth <= 2) {
0136         if (G4UniformRand() > 0.5)
0137           id.depth += 2;
0138       }
0139     }
0140   }
0141 }