CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
FastHFShowerLibrary.cc
Go to the documentation of this file.
1 // File: FastHFShowerLibrary.cc
3 // Description: Shower library for Very forward hadron calorimeter
5 
20 
21 #include "Randomize.hh"
22 #include "CLHEP/Units/GlobalSystemOfUnits.h"
23 #include "CLHEP/Units/GlobalPhysicalConstants.h"
24 
25 // Geant4 headers
26 #include "G4ParticleDefinition.hh"
27 #include "G4DynamicParticle.hh"
28 #include "G4DecayPhysics.hh"
29 #include "G4ParticleTable.hh"
30 #include "G4ParticleTypes.hh"
31 
32 // STL headers
33 #include <iostream>
34 #include <mutex>
35 #include <vector>
36 
37 //#define DebugLog
38 
39 static std::once_flag initializeOnce;
40 
42  : fast(p), hcalDDDSimConstantsESToken_(iC.esConsumes()), hcalSimulationConstantsESToken_(iC.esConsumes()) {
43  edm::ParameterSet m_HS = p.getParameter<edm::ParameterSet>("HFShowerLibrary");
44  applyFidCut = m_HS.getParameter<bool>("ApplyFiducialCut");
45 }
46 
48  edm::LogInfo("FastCalorimetry") << "initHFShowerLibrary::initialization";
49 
52 
53  std::string name = "HcalHits";
54  numberingFromDDD = std::make_unique<HcalNumberingFromDDD>(hcalConstants);
55  hfshower = std::make_unique<HFShowerLibrary>(name, hcalConstants, hsps->hcalsimpar(), fast);
56 
57  //only one thread can be allowed to setup the G4 physics table.
58  std::call_once(initializeOnce, []() {
59  // Geant4 particles
60  G4DecayPhysics decays;
61  decays.ConstructParticle();
62  G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
63  partTable->SetReadiness();
64  });
65  //G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
66  //hfshower->initRun(partTable, hcalConstants); // init particle code
67 }
68 
70  // define Geant4 engine per thread
71  G4Random::setTheEngine(&(rnd->theEngine()));
72  LogDebug("FastHFShowerLibrary::recoHFShowerLibrary")
73  << "Begin of event " << G4UniformRand() << " " << rnd->theEngine().name() << " " << rnd->theEngine();
74 }
75 
77 #ifdef DebugLog
78  edm::LogInfo("FastCalorimetry") << "FastHFShowerLibrary: recoHFShowerLibrary ";
79 #endif
80 
81  if (!myTrack.onVFcal()) {
82 #ifdef DebugLog
83  edm::LogInfo("FastCalorimetry") << "FastHFShowerLibrary: we should not be here ";
84 #endif
85  }
86 
87  hitMap.clear();
88  double eGen = 1000. * myTrack.vfcalEntrance().e(); // energy in [MeV]
89  double delZv = (myTrack.vfcalEntrance().vertex().Z() > 0.0) ? 50.0 : -50.0;
90  G4ThreeVector vertex(10. * myTrack.vfcalEntrance().vertex().X(),
91  10. * myTrack.vfcalEntrance().vertex().Y(),
92  10. * myTrack.vfcalEntrance().vertex().Z() + delZv); // in [mm]
93 
94  G4ThreeVector direction(
95  myTrack.vfcalEntrance().Vect().X(), myTrack.vfcalEntrance().Vect().Y(), myTrack.vfcalEntrance().Vect().Z());
96 
97  bool ok;
98  double weight = 1.0; // rad. damage
99  int parCode = myTrack.type();
100  double tSlice = 0.1 * vertex.mag() / 29.98;
101 
102  std::vector<HFShowerLibrary::Hit> hits =
103  hfshower->fillHits(vertex, direction, parCode, eGen, ok, weight, tSlice, false);
104 
105  for (unsigned int i = 0; i < hits.size(); ++i) {
106  G4ThreeVector pos = hits[i].position;
107  int depth = hits[i].depth;
108  double time = hits[i].time;
109  if (!applyFidCut || (HFFibreFiducial::PMTNumber(pos) > 0)) {
110  // if (!applyFidCut || (applyFidCut && HFFibreFiducial::PMTNumber(pos)>0)) {
111  int det = 5;
112  int lay = 1;
113  uint32_t id = 0;
115  numberingFromDDD->unitID(det, math::XYZVectorD(pos.x(), pos.y(), pos.z()), depth, lay);
116  modifyDepth(tmp);
117  id = numberingScheme.getUnitID(tmp);
118 
119  CaloHitID current_id(id, time, myTrack.id());
120  std::map<CaloHitID, float>::iterator cellitr;
121  cellitr = hitMap.find(current_id);
122  if (cellitr == hitMap.end()) {
123  hitMap.insert(std::pair<CaloHitID, float>(current_id, 1.0));
124  } else {
125  cellitr->second += 1.0;
126  }
127  } // end of isItinFidVolume check
128  } // end loop over hits
129 }
130 
132  if (id.subdet == HcalForward) {
133  int ieta = (id.zside == 0) ? -id.etaR : id.etaR;
134  if (hcalConstants->maxHFDepth(ieta, id.phis) > 2) {
135  if (id.depth <= 2) {
136  if (G4UniformRand() > 0.5)
137  id.depth += 2;
138  }
139  }
140  }
141 }
const edm::ESGetToken< HcalSimulationConstants, HcalSimNumberingRecord > hcalSimulationConstantsESToken_
std::map< CaloHitID, float > hitMap
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
const RawParticle & vfcalEntrance() const
The particle at VFCAL entrance.
Definition: FSimTrack.h:149
const edm::ParameterSet fast
void recoHFShowerLibrary(const FSimTrack &myTrack)
virtual uint32_t getUnitID(const HcalNumberingFromDDD::HcalID &id)
const HcalSimulationParameters * hcalsimpar() const
bool getData(T &iHolder) const
Definition: EventSetup.h:122
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
HcalNumberingScheme numberingScheme
double e() const
energy of the momentum
Definition: RawParticle.h:305
const edm::ESGetToken< HcalDDDSimConstants, HcalSimNumberingRecord > hcalDDDSimConstantsESToken_
FastHFShowerLibrary(edm::ParameterSet const &, edm::ConsumesCollector &&)
static std::once_flag initializeOnce
int onVFcal() const
Definition: FSimTrack.h:121
CLHEP::HepRandomEngine & theEngine() const
std::unique_ptr< HcalNumberingFromDDD > numberingFromDDD
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
Log< level::Info, false > LogInfo
int PMTNumber(const G4ThreeVector &pe_effect)
void modifyDepth(HcalNumberingFromDDD::HcalID &id)
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
void SetRandom(const RandomEngineAndDistribution *)
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
std::unique_ptr< HFShowerLibrary > hfshower
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:96
const HcalDDDSimConstants * hcalConstants
int weight
Definition: histoStyle.py:51
void const initHFShowerLibrary(const edm::EventSetup &)
tmp
align.sh
Definition: createJobs.py:716
int maxHFDepth(const int &ieta, const int &iphi) const
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
#define LogDebug(id)