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