CMS 3D CMS Logo
HFShower.cc
Go to the documentation of this file.
1 // File: HFShower.cc
3 // Description: Sensitive Detector class for calorimeters
5 
6 #include "SimG4CMS/Calo/interface/HFShower.h"
7 #include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
8 #include "DetectorDescription/Core/interface/DDFilter.h"
9 #include "DetectorDescription/Core/interface/DDFilteredView.h"
10 #include "DetectorDescription/Core/interface/DDSolid.h"
11 #include "DetectorDescription/Core/interface/DDSplit.h"
12 #include "DetectorDescription/Core/interface/DDValue.h"
13 
14 #include "G4NavigationHistory.hh"
15 #include "G4Step.hh"
16 #include "G4Track.hh"
17 #include "G4VSolid.hh"
18 #include "Randomize.hh"
19 #include "CLHEP/Units/GlobalPhysicalConstants.h"
20 #include "CLHEP/Units/GlobalSystemOfUnits.h"
21 
22 //#define EDM_ML_DEBUG
23 
24 #include <iostream>
25 
26 HFShower::HFShower(const std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk)
27  : cherenkov(nullptr), fibre(nullptr), chkFibre(chk) {
28  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
29  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
30  probMax = m_HF.getParameter<double>("ProbMax");
31 
32  edm::LogVerbatim("HFShower") << "HFShower:: Maximum probability cut off " << probMax << " Check flag " << chkFibre;
33 
34  cherenkov = new HFCherenkov(m_HF);
35  fibre = new HFFibre(name, cpv, p);
36 }
37 
38 HFShower::~HFShower() {
39  if (cherenkov)
40  delete cherenkov;
41  if (fibre)
42  delete fibre;
43 }
44 
45 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, double weight) {
46  std::vector<HFShower::Hit> hits;
47  int nHit = 0;
48 
49  double edep = weight * (aStep->GetTotalEnergyDeposit());
50 #ifdef EDM_ML_DEBUG
51  edm::LogVerbatim("HFShower") << "HFShower::getHits: energy " << aStep->GetTotalEnergyDeposit() << " weight " << weight
52  << " edep " << edep;
53 #endif
54  double stepl = 0.;
55 
56  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
57  stepl = aStep->GetStepLength();
58  if ((edep == 0.) || (stepl == 0.)) {
59 #ifdef EDM_ML_DEBUG
60  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
61 #endif
62  return hits;
63  }
64  const G4Track *aTrack = aStep->GetTrack();
65  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
66 
67  HFShower::Hit hit;
68  double energy = aParticle->GetTotalEnergy();
69  double momentum = aParticle->GetTotalMomentum();
70  double pBeta = momentum / energy;
71  double dose = 0.;
72  int npeDose = 0;
73 
74  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
75  const G4ParticleDefinition *particleDef = aTrack->GetDefinition();
76 
77  const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
78  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
79  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
80  //double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
81  double zv = std::abs(globalPos.z()) - gpar[4];
82  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
83  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
84  // @@ Here the depth should be changed (Fibers all long in Geometry!)
85  int depth = 1;
86  int npmt = 0;
87  bool ok = true;
88  if (zv < 0. || zv > gpar[1]) {
89  ok = false; // beyond the fiber in Z
90  }
91  if (ok && applyFidCut) {
92  npmt = HFFibreFiducial::PMTNumber(globalPos);
93  if (npmt <= 0) {
94  ok = false;
95  } else if (npmt > 24) { // a short fibre
96  if (zv > gpar[0]) {
97  depth = 2;
98  } else {
99  ok = false;
100  }
101  }
102 #ifdef EDM_ML_DEBUG
103  edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar[4]
104  << ":" << zv << ":" << gpar[0] << " ok " << ok << " depth " << depth;
105 #endif
106  } else {
107  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
108  }
109  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
110 
111  double u = localMom.x();
112  double v = localMom.y();
113  double w = localMom.z();
114  double zCoor = localPos.z();
115  double zFibre = (0.5 * gpar[1] - zCoor - translation.z());
116  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
117  double time = fibre->tShift(localPos, depth, chkFibre);
118 
119 #ifdef EDM_ML_DEBUG
120  edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
121  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
122  << "\n Direction " << momentumDir << " Local " << localMom;
123 #endif
124  // Here npe should be 0 if there is no fiber (@@ M.K.)
125  int npe = 1;
126  std::vector<double> wavelength;
127  std::vector<double> momz;
128  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
129  if (ok)
130  npe = cherenkov->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
131  wavelength = cherenkov->getWL();
132  momz = cherenkov->getMom();
133  } // ^^^^^ End of Tmp close of the cherenkov function
134  if (ok && npe > 0) {
135  for (int i = 0; i < npe; ++i) {
136  double p = 1.;
137  if (!applyFidCut)
138  p = fibre->attLength(wavelength[i]);
139  double r1 = G4UniformRand();
140  double r2 = G4UniformRand();
141 #ifdef EDM_ML_DEBUG
142  edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
143  << " r2 " << r2 << ":" << probMax
144  << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax);
145 #endif
146  if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax)) {
147  hit.depth = depth;
148  hit.time = tSlice + time;
149  if (!applyFidCut) {
150  hit.wavelength = wavelength[i]; // Tmp
151  hit.momentum = momz[i]; // Tmp
152  } else {
153  hit.wavelength = 300.; // Tmp
154  hit.momentum = 1.; // Tmp
155  }
156  hit.position = globalPos;
157  hits.push_back(hit);
158  nHit++;
159  }
160  }
161  }
162 
163 #ifdef EDM_ML_DEBUG
164  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
165  for (int i = 0; i < nHit; ++i)
166  edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
167  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
168  << hits[i].position;
169 #endif
170  return hits;
171 }
172 
173 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibraryProducer, double zoffset) {
174  std::vector<HFShower::Hit> hits;
175  int nHit = 0;
176 
177  double edep = aStep->GetTotalEnergyDeposit();
178  double stepl = 0.;
179 
180  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
181  stepl = aStep->GetStepLength();
182  if ((edep == 0.) || (stepl == 0.)) {
183 #ifdef EDM_ML_DEBUG
184  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
185 #endif
186  return hits;
187  }
188  const G4Track *aTrack = aStep->GetTrack();
189  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
190 
191  HFShower::Hit hit;
192  double energy = aParticle->GetTotalEnergy();
193  double momentum = aParticle->GetTotalMomentum();
194  double pBeta = momentum / energy;
195  double dose = 0.;
196  int npeDose = 0;
197 
198  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
199  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
200 
201  G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
202  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
203  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
204  //double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
205  //double zv = std::abs(globalPos.z()) - gpar[4];
206  double zv = gpar[1] - (std::abs(globalPos.z()) - zoffset);
207  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
208  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
209  // @@ Here the depth should be changed (Fibers all long in Geometry!)
210  int depth = 1;
211  int npmt = 0;
212  bool ok = true;
213  if (zv < 0. || zv > gpar[1]) {
214  ok = false; // beyond the fiber in Z
215  }
216  if (ok && applyFidCut) {
217  npmt = HFFibreFiducial::PMTNumber(globalPos);
218  if (npmt <= 0) {
219  ok = false;
220  } else if (npmt > 24) { // a short fibre
221  if (zv > gpar[0]) {
222  depth = 2;
223  } else {
224  ok = false;
225  }
226  }
227 #ifdef EDM_ML_DEBUG
228  edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar[4]
229  << ":" << zv << ":" << gpar[0] << " ok " << ok << " depth " << depth;
230 #endif
231  } else {
232  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
233  }
234  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
235 
236  double u = localMom.x();
237  double v = localMom.y();
238  double w = localMom.z();
239  double zCoor = localPos.z();
240  double zFibre = (0.5 * gpar[1] - zCoor - translation.z());
241  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
242  double time = fibre->tShift(localPos, depth, chkFibre);
243 
244 #ifdef EDM_ML_DEBUG
245  edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
246  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
247  << "\n Direction " << momentumDir << " Local " << localMom;
248 #endif
249  // Here npe should be 0 if there is no fiber (@@ M.K.)
250  int npe = 1;
251  std::vector<double> wavelength;
252  std::vector<double> momz;
253  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
254  if (ok)
255  npe = cherenkov->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
256  wavelength = cherenkov->getWL();
257  momz = cherenkov->getMom();
258  } // ^^^^^ End of Tmp close of the cherenkov function
259  if (ok && npe > 0) {
260  for (int i = 0; i < npe; ++i) {
261  double p = 1.;
262  if (!applyFidCut)
263  p = fibre->attLength(wavelength[i]);
264  double r1 = G4UniformRand();
265  double r2 = G4UniformRand();
266 #ifdef EDM_ML_DEBUG
267  edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
268  << " r2 " << r2 << ":" << probMax
269  << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax);
270 #endif
271  if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax)) {
272  hit.depth = depth;
273  hit.time = tSlice + time;
274  if (!applyFidCut) {
275  hit.wavelength = wavelength[i]; // Tmp
276  hit.momentum = momz[i]; // Tmp
277  } else {
278  hit.wavelength = 300.; // Tmp
279  hit.momentum = 1.; // Tmp
280  }
281  hit.position = globalPos;
282  hits.push_back(hit);
283  nHit++;
284  }
285  }
286  }
287 
288 #ifdef EDM_ML_DEBUG
289  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
290  for (int i = 0; i < nHit; ++i)
291  edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
292  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
293  << hits[i].position;
294 #endif
295  return hits;
296 }
297 
298 std::vector<HFShower::Hit> HFShower::getHits(const G4Step *aStep, bool forLibrary) {
299  std::vector<HFShower::Hit> hits;
300  int nHit = 0;
301 
302  double edep = aStep->GetTotalEnergyDeposit();
303  double stepl = 0.;
304 
305  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
306  stepl = aStep->GetStepLength();
307  if ((edep == 0.) || (stepl == 0.)) {
308 #ifdef EDM_ML_DEBUG
309  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
310 #endif
311  return hits;
312  }
313 
314  const G4Track *aTrack = aStep->GetTrack();
315  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
316 
317  HFShower::Hit hit;
318  double energy = aParticle->GetTotalEnergy();
319  double momentum = aParticle->GetTotalMomentum();
320  double pBeta = momentum / energy;
321  double dose = 0.;
322  int npeDose = 0;
323 
324  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
325  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
326 
327  const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
328  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
329  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
330  double zv = std::abs(globalPos.z()) - gpar[4] - 0.5 * gpar[1];
331  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
332  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
333  // @@ Here the depth should be changed (Fibers are all long in Geometry!)
334  int depth = 1;
335  int npmt = 0;
336  bool ok = true;
337  if (zv < 0 || zv > gpar[1]) {
338  ok = false; // beyond the fiber in Z
339  }
340  if (ok && applyFidCut) {
341  npmt = HFFibreFiducial::PMTNumber(globalPos);
342  if (npmt <= 0) {
343  ok = false;
344  } else if (npmt > 24) { // a short fibre
345  if (zv > gpar[0]) {
346  depth = 2;
347  } else {
348  ok = false;
349  }
350  }
351 #ifdef EDM_ML_DEBUG
352  edm::LogVerbatim("HFShower") << "HFShower:getHits(SL): npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":"
353  << gpar[4] << ":" << zv << ":" << gpar[0] << " ok " << ok << " depth " << depth;
354 #endif
355  } else {
356  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
357  }
358  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
359 
360  double u = localMom.x();
361  double v = localMom.y();
362  double w = localMom.z();
363  double zCoor = localPos.z();
364  double zFibre = (0.5 * gpar[1] - zCoor - translation.z());
365  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
366  double time = fibre->tShift(localPos, depth, chkFibre);
367 
368 #ifdef EDM_ML_DEBUG
369  edm::LogVerbatim("HFShower") << "HFShower::getHits(SL): in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
370  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
371  << "\n Direction " << momentumDir << " Local " << localMom;
372 #endif
373  // npe should be 0
374  int npe = 0;
375  if (ok)
376  npe = cherenkov->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
377  std::vector<double> wavelength = cherenkov->getWL();
378  std::vector<double> momz = cherenkov->getMom();
379 
380  for (int i = 0; i < npe; ++i) {
381  hit.time = tSlice + time;
382  hit.wavelength = wavelength[i];
383  hit.momentum = momz[i];
384  hit.position = globalPos;
385  hits.push_back(hit);
386  nHit++;
387  }
388 
389  return hits;
390 }
391 
392 std::vector<double> HFShower::getDDDArray(const std::string &str, const DDsvalues_type &sv, int &nmin) {
393 #ifdef EDM_ML_DEBUG
394  edm::LogVerbatim("HFShower") << "HFShower:getDDDArray called for " << str << " with nMin " << nmin;
395 #endif
396  DDValue value(str);
397  if (DDfetch(&sv, value)) {
398 #ifdef EDM_ML_DEBUG
399  edm::LogVerbatim("HFShower") << value;
400 #endif
401  const std::vector<double> &fvec = value.doubles();
402  int nval = fvec.size();
403  if (nmin > 0) {
404  if (nval < nmin) {
405  edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
406  throw cms::Exception("Unknown", "HFShower") << "nval < nmin for array " << str << "\n";
407  }
408  } else {
409  if (nval < 2) {
410  edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval << " < 2 ==> illegal"
411  << " (nmin=" << nmin << ")";
412  throw cms::Exception("Unknown", "HFShower") << "nval < 2 for array " << str;
413  }
414  }
415  nmin = nval;
416 
417  return fvec;
418  } else {
419  edm::LogError("HFShower") << "HFShower : cannot get array " << str;
420  throw cms::Exception("Unknown", "HFShower") << "cannot get array " << str << "\n";
421  }
422 }
423 
424 void HFShower::initRun(const HcalDDDSimConstants *hcons) {
425  //Special Geometry parameters
426  gpar = hcons->getGparHF();
427  if (fibre) {
428  fibre->initRun(hcons);
429  }
430 }
HFFibre::tShift
double tShift(const G4ThreeVector &point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:113
HFShower::Hit::wavelength
double wavelength
Definition: HFShower.h:32
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HFShower::fibre
HFFibre * fibre
Definition: HFShower.h:48
DDValue::doubles
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:140
mps_fire.i
i
Definition: mps_fire.py:338
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HFFibre::initRun
void initRun(const HcalDDDSimConstants *)
Definition: HFFibre.cc:81
w
const double w
Definition: UKUtility.cc:23
HFShower::Hit::momentum
double momentum
Definition: HFShower.h:33
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
HFFibre::attLength
double attLength(double lambda)
Definition: HFFibre.cc:97
nullptr
#define nullptr
Definition: GCC11Compatibility.h:37
HFShower::applyFidCut
bool applyFidCut
Definition: HFShower.h:44
HFShower::Hit::position
G4ThreeVector position
Definition: HFShower.h:34
weight
Definition: weight.py:1
pfDeepBoostedJetPreprocessParams_cfi.sv
sv
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:352
DDCompactView
Compact representation of the geometrical detector hierarchy.
Definition: DDCompactView.h:80
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:81
HFShower::getHits
std::vector< Hit > getHits(const G4Step *aStep, double weight)
Definition: HFShower.cc:45
HFShower::cherenkov
HFCherenkov * cherenkov
Definition: HFShower.h:47
HFShower::chkFibre
int chkFibre
Definition: HFShower.h:50
HFShower::~HFShower
virtual ~HFShower()
Definition: HFShower.h:41
HFShower::probMax
double probMax
Definition: HFShower.h:51
HFCherenkov::getWL
std::vector< double > getWL()
Definition: HFCherenkov.cc:322
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HFCherenkov::getMom
std::vector< double > getMom()
Definition: HFCherenkov.cc:327
HFShower::gpar
std::vector< double > gpar
Definition: HFShower.h:52
HFShower::Hit::time
double time
Definition: HFShower.h:31
HcalDDDSimConstants::getGparHF
const std::vector< double > & getGparHF() const
Definition: HcalDDDSimConstants.h:51
HFShower::HFShower
HFShower(const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
Definition: HFShower.cc:30
HFShower::initRun
void initRun(const HcalDDDSimConstants *)
Definition: HFShower.cc:424
HFFibreFiducial::PMTNumber
int PMTNumber(const G4ThreeVector &pe_effect)
Definition: HFFibreFiducial.cc:9
protons_cff.time
time
Definition: protons_cff.py:35
HFShower::getDDDArray
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
Definition: HFShower.cc:392
DDsvalues_type
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
Definition: DDsvalues.h:12
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
HFCherenkov::computeNPE
int computeNPE(const G4Step *step, const G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:77
str
#define str(s)
Definition: TestProcessor.cc:48
egammaForCoreTracking_cff.depth
depth
Definition: egammaForCoreTracking_cff.py:29
Generated for CMSSW Reference Manual by doxygen 1.8.11