CMS 3D CMS Logo

List of all members | Classes | Public Types | Public Member Functions | Private Member Functions | Private Attributes
HFShower Class Reference

#include <HFShower.h>

Classes

struct  Hit
 

Public Types

typedef std::pair< XYZPoint, double > Spot
 
typedef std::pair< unsigned int, double > Step
 
typedef Steps::const_iterator step_iterator
 
typedef std::vector< StepSteps
 
typedef math::XYZVector XYZPoint
 

Public Member Functions

bool compute ()
 Compute the shower longitudinal and lateral development. More...
 
std::vector< HitgetHits (const G4Step *aStep, bool forLibrary)
 
std::vector< HitgetHits (const G4Step *aStep, bool forLibraryProducer, double zoffset)
 
std::vector< HitgetHits (const G4Step *aStep, double weight)
 
 HFShower (const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
 
 HFShower (const std::string &name, const HcalDDDSimConstants *hcons, const HcalSimulationParameters *hps, edm::ParameterSet const &p, int chk=0)
 
virtual ~HFShower ()
 
virtual ~HFShower ()
 

Private Member Functions

double gam (double x, double a) const
 
int indexFinder (double x, const std::vector< double > &Fhist)
 
void makeSteps (int nsteps)
 
double transProb (double factor, double R, double r)
 

Private Attributes

double aloge
 
double alpEM
 
double alpHD
 
bool applyFidCut_
 
double balanceEH
 
double betEM
 
double betHD
 
std::unique_ptr< HFCherenkovcherenkov_
 
int chkFibre_
 
double criticalEnergy
 
double depthStart
 
double depthStep
 
std::vector< int > detector
 
double e
 
double eSpotSize
 
std::vector< double > eStep
 
std::unique_ptr< HFFibrefibre_
 
std::vector< double > gpar_
 
const HcalDDDSimConstantshcalConstant_
 
double hcalDepthFactor
 
int infinity
 
double lambdaEM
 
double lambdaHD
 
std::vector< double > lamcurr
 
std::vector< double > lamdepth
 
std::vector< double > lamstep
 
std::vector< double > lamtotal
 
int lossesOpt
 
double maxTRfactor
 
int nDepthSteps
 
std::vector< int > nspots
 
int nTRsteps
 
int onEcal
 
double part
 
double probMax_
 
const RandomEngineAndDistributionrandom
 
std::vector< double > rlamStep
 
double tgamEM
 
double tgamHD
 
const ECALPropertiestheECALproperties
 
EcalHitMakertheGrid
 
HcalHitMakertheHcalHitMaker
 
const HCALPropertiestheHCALproperties
 
HDShowerParametrizationtheParam
 
double theR1
 
double theR2
 
double theR3
 
double transFactor
 
double transParam
 
std::vector< double > x0curr
 
std::vector< double > x0depth
 
double x0EM
 
double x0HD
 

Detailed Description

Definition at line 22 of file HFShower.h.

Member Typedef Documentation

◆ Spot

typedef std::pair<XYZPoint, double> HFShower::Spot

Definition at line 26 of file HFShower.h.

◆ Step

typedef std::pair<unsigned int, double> HFShower::Step

Definition at line 27 of file HFShower.h.

◆ step_iterator

typedef Steps::const_iterator HFShower::step_iterator

Definition at line 29 of file HFShower.h.

◆ Steps

typedef std::vector<Step> HFShower::Steps

Definition at line 28 of file HFShower.h.

◆ XYZPoint

Definition at line 24 of file HFShower.h.

Constructor & Destructor Documentation

◆ HFShower() [1/2]

HFShower::HFShower ( const RandomEngineAndDistribution engine,
HDShowerParametrization myParam,
EcalHitMaker myGrid,
HcalHitMaker myHcalHitMaker,
int  onECAL,
double  epart 
)

Definition at line 29 of file HFShower.cc.

35  : theParam(myParam), theGrid(myGrid), theHcalHitMaker(myHcalHitMaker), onEcal(onECAL), e(epart), random(engine) {
36  // To get an access to constants read in FASTCalorimeter
37  // FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();
38 
39  // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
40  lossesOpt = myParam->hsParameters()->getHDlossesOpt();
42  nTRsteps = myParam->hsParameters()->getHDnTRsteps();
43  transParam = myParam->hsParameters()->getHDtransParam();
44  eSpotSize = myParam->hsParameters()->getHDeSpotSize();
45  depthStep = myParam->hsParameters()->getHDdepthStep();
48  balanceEH = myParam->hsParameters()->getHDbalanceEH();
50 
51  // Special tr.size fluctuations
52  transParam *= (1. + random->flatShoot());
53 
54  // Special ad hoc long. extension + some fluctuations
55  double depthExt;
56  if (e < 50.)
57  depthExt = 0.8 * (50. - e) / 50. + 0.3;
58  else {
59  if (e < 500.)
60  depthExt = (500. - e) / 500. * 0.4 - 0.1;
61  else
62  depthExt = -0.1;
63  }
64  hcalDepthFactor += depthExt + 0.05 * (2. * random->flatShoot() - 1.);
65 
66  // normally 1, in HF - might be smaller to take into account
67  // a narrowness of the HF shower (Cherenkov light)
68  if (e < 50.)
69  transFactor = 0.5 - (50. - e) / 50. * 0.2;
70  else
71  transFactor = 0.7 - (1000. - e) / 1000. * 0.2;
72 
73  // simple protection ...
74  if (e < 0)
75  e = 0.;
76 
77  // Get the Famos Histos pointer
78  // myHistos = FamosHistos::instance();
79  // std::cout << " Hello FamosShower " << std::endl;
80 
83 
84  double emax = theParam->emax();
85  double emid = theParam->emid();
86  double emin = theParam->emin();
87  double effective = e;
88 
89  if (e < emid) {
90  theParam->setCase(1);
91  // avoid "underflow" below Emin (for parameters calculation only)
92  if (e < emin)
93  effective = emin;
94  } else
95  theParam->setCase(2);
96 
97  // A bit coarse espot size for HF...
98  eSpotSize *= 2.5;
99  if (effective > 0.5 * emax) {
100  eSpotSize *= 2.;
101  if (effective > emax) {
102  effective = emax;
103  eSpotSize *= 3.;
104  depthStep *= 2.;
105  }
106  }
107 
108  if (debug == 2)
109  LogDebug("FastCalorimetry") << " HFShower : " << std::endl
110  << " Energy " << e << std::endl
111  << " lossesOpt " << lossesOpt << std::endl
112  << " nDepthSteps " << nDepthSteps << std::endl
113  << " nTRsteps " << nTRsteps << std::endl
114  << " transParam " << transParam << std::endl
115  << " eSpotSize " << eSpotSize << std::endl
116  << " criticalEnergy " << criticalEnergy << std::endl
117  << " maxTRfactor " << maxTRfactor << std::endl
118  << " balanceEH " << balanceEH << std::endl
119  << "hcalDepthFactor " << hcalDepthFactor << std::endl;
120 
121  double alpEM1 = theParam->alpe1();
122  double alpEM2 = theParam->alpe2();
123 
124  double betEM1 = theParam->bete1();
125  double betEM2 = theParam->bete2();
126 
127  double alpHD1 = theParam->alph1();
128  double alpHD2 = theParam->alph2();
129 
130  double betHD1 = theParam->beth1();
131  double betHD2 = theParam->beth2();
132 
133  double part1 = theParam->part1();
134  double part2 = theParam->part2();
135 
136  aloge = std::log(effective);
137 
138  double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
139  double aedep = std::log(edpar);
140 
141  if (debug == 2)
142  LogDebug("FastCalorimetry") << " HFShower : " << std::endl
143  << " edpar " << edpar << " aedep " << aedep << std::endl
144  << " alpEM1 " << alpEM1 << std::endl
145  << " alpEM2 " << alpEM2 << std::endl
146  << " betEM1 " << betEM1 << std::endl
147  << " betEM2 " << betEM2 << std::endl
148  << " alpHD1 " << alpHD1 << std::endl
149  << " alpHD2 " << alpHD2 << std::endl
150  << " betHD1 " << betHD1 << std::endl
151  << " betHD2 " << betHD2 << std::endl
152  << " part1 " << part1 << std::endl
153  << " part2 " << part2 << std::endl;
154 
155  // private members to set
156  theR1 = theParam->r1();
157  theR2 = theParam->r2();
158  theR3 = theParam->r3();
159 
160  alpEM = alpEM1 + alpEM2 * aedep;
161  tgamEM = tgamma(alpEM);
162  betEM = betEM1 - betEM2 * aedep;
163  alpHD = alpHD1 + alpHD2 * aedep;
164  tgamHD = tgamma(alpHD);
165  betHD = betHD1 - betHD2 * aedep;
166  part = part1 - part2 * aedep;
167  if (part > 1.)
168  part = 1.; // protection - just in case of
169 
170  if (debug == 2)
171  LogDebug("FastCalorimetry") << " HFShower : " << std::endl
172  << " alpEM " << alpEM << std::endl
173  << " tgamEM " << tgamEM << std::endl
174  << " betEM " << betEM << std::endl
175  << " alpHD " << alpHD << std::endl
176  << " tgamHD " << tgamHD << std::endl
177  << " betHD " << betHD << std::endl
178  << " part " << part << std::endl;
179 
180  if (onECAL) {
183  } else {
184  lambdaEM = 0.;
185  x0EM = 0.;
186  }
189 
190  if (debug == 2)
191  LogDebug("FastCalorimetry") << " HFShower e " << e << std::endl
192  << " x0EM = " << x0EM << std::endl
193  << " x0HD = " << x0HD << std::endl
194  << " lamEM = " << lambdaEM << std::endl
195  << " lamHD = " << lambdaHD << std::endl;
196 
197  // Starting point of the shower
198  // try first with ECAL lambda
199 
200  double sum1 = 0.; // lambda path from the ECAL/HF entrance;
201  double sum2 = 0.; // lambda path from the interaction point;
202  double sum3 = 0.; // x0 path from the interaction point;
203  int nsteps = 0; // full number of longitudinal steps (counter);
204 
205  int nmoresteps; // how many longitudinal steps in addition to
206  // one (if interaction happens there) in ECAL
207 
208  if (e < criticalEnergy)
209  nmoresteps = 1;
210  else
211  nmoresteps = nDepthSteps;
212 
213  double depthECAL = 0.;
214  double depthGAP = 0.;
215  double depthGAPx0 = 0.;
216  if (onECAL) {
217  depthECAL = theGrid->ecalTotalL0(); // ECAL depth segment
218  depthGAP = theGrid->ecalHcalGapTotalL0(); // GAP depth segment
219  depthGAPx0 = theGrid->ecalHcalGapTotalX0(); // GAP depth x0
220  }
221 
222  double depthHCAL = theGrid->hcalTotalL0(); // HCAL depth segment
223  double depthToHCAL = depthECAL + depthGAP;
224 
225  //---------------------------------------------------------------------------
226  // Depth simulation & various protections, among them
227  // if too deep - get flat random in the allowed region
228  // if no HCAL material behind - force to deposit in ECAL
229  double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
230 
231  double depthStart = std::log(1. / random->flatShoot()); // starting point lambda unts
232 
233  if (e < emin) {
234  if (debug)
235  LogDebug("FastCalorimetry") << " FamosHFShower : e <emin -> depthStart = 0" << std::endl;
236  depthStart = 0.;
237  }
238 
239  if (depthStart > maxDepth) {
240  if (debug)
241  LogDebug("FastCalorimetry") << " FamosHFShower : depthStart too big ... = " << depthStart << std::endl;
242 
244  if (depthStart < 0.)
245  depthStart = 0.;
246  if (debug)
247  LogDebug("FastCalorimetry") << " FamosHFShower : depthStart re-calculated = " << depthStart << std::endl;
248  }
249 
250  if (onECAL && e < emid) {
251  if ((depthECAL - depthStart) / depthECAL > 0.2 && depthECAL > depthStep) {
252  depthStart = 0.5 * depthECAL * random->flatShoot();
253  if (debug)
254  LogDebug("FastCalorimetry") << " FamosHFShower : small energy, "
255  << " depthStart reduced to = " << depthStart << std::endl;
256  }
257  }
258 
259  if (depthHCAL < depthStep) {
260  if (debug)
261  LogDebug("FastCalorimetry") << " FamosHFShower : depthHCAL too small ... = " << depthHCAL
262  << " depthStart -> forced to 0 !!!" << std::endl;
263  depthStart = 0.;
264  nmoresteps = 0;
265 
266  if (depthECAL < depthStep) {
267  nsteps = -1;
268  LogInfo("FastCalorimetry") << " FamosHFShower : too small ECAL and HCAL depths - "
269  << " particle is lost !!! " << std::endl;
270  }
271  }
272 
273  if (debug)
274  LogDebug("FastCalorimetry") << " FamosHFShower depths(lam) - " << std::endl
275  << " ECAL = " << depthECAL << std::endl
276  << " GAP = " << depthGAP << std::endl
277  << " HCAL = " << depthHCAL << std::endl
278  << " starting point = " << depthStart << std::endl;
279 
280  if (onEcal) {
281  if (debug)
282  LogDebug("FastCalorimetry") << " FamosHFShower : onECAL" << std::endl;
283  if (depthStart < depthECAL) {
284  if (debug)
285  LogDebug("FastCalorimetry") << " FamosHFShower : depthStart < depthECAL" << std::endl;
286  if ((depthECAL - depthStart) / depthECAL > 0.25 && depthECAL > depthStep) {
287  if (debug)
288  LogDebug("FastCalorimetry") << " FamosHFShower : enough space to make ECAL step" << std::endl;
289  // ECAL - one step
290  nsteps++;
291  sum1 += depthECAL; // at the end of step
292  sum2 += depthECAL - depthStart;
293  sum3 += sum2 * lambdaEM / x0EM;
294  lamtotal.push_back(sum1);
295  lamdepth.push_back(sum2);
296  lamcurr.push_back(lambdaEM);
297  lamstep.push_back(depthECAL - depthStart);
298  x0depth.push_back(sum3);
299  x0curr.push_back(x0EM);
300  detector.push_back(1);
301 
302  if (debug)
303  LogDebug("FastCalorimetry") << " FamosHFShower : "
304  << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;
305 
306  // // Gap - no additional step after ECAL
307  // // just move further to HCAL over the gap
308  sum1 += depthGAP;
309  sum2 += depthGAP;
310  sum3 += depthGAPx0;
311  }
312  // Just shift starting point to HCAL
313  else {
314  // cout << " FamosHFShower : not enough space to make ECAL step" << std::endl;
315  if (debug)
316  LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
317 
318  depthStart = depthToHCAL;
319  sum1 += depthStart;
320  }
321  } else { // GAP or HCAL
322 
323  if (depthStart >= depthECAL && depthStart < depthToHCAL) {
324  depthStart = depthToHCAL; // just a shift to HCAL for simplicity
325  }
326  sum1 += depthStart;
327 
328  if (debug)
329  LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
330  }
331  } else { // Forward
332  if (debug)
333  LogDebug("FastCalorimetry") << " FamosHFShower : forward" << std::endl;
334  sum1 += depthStart;
335  }
336 
337  for (int i = 0; i < nmoresteps; i++) {
338  sum1 += depthStep;
339  if (sum1 > (depthECAL + depthGAP + depthHCAL))
340  break;
341  sum2 += depthStep;
342  sum3 += sum2 * lambdaHD / x0HD;
343  lamtotal.push_back(sum1);
344  lamdepth.push_back(sum2);
345  lamcurr.push_back(lambdaHD);
346  lamstep.push_back(depthStep);
347  x0depth.push_back(sum3);
348  x0curr.push_back(x0HD);
349  detector.push_back(3);
350  nsteps++;
351  }
352 
353  // Make fractions of energy and transverse radii at each step
354 
355  // PV
356  // std::cout << "HFShower::HFShower() : Nsteps = " << nsteps << std::endl;
357 
358  if (nsteps > 0) {
359  makeSteps(nsteps);
360  }
361 }

References aloge, HDShowerParametrization::alpe1(), HDShowerParametrization::alpe2(), alpEM, HDShowerParametrization::alph1(), HDShowerParametrization::alph2(), alpHD, balanceEH, HDShowerParametrization::bete1(), HDShowerParametrization::bete2(), betEM, HDShowerParametrization::beth1(), HDShowerParametrization::beth2(), betHD, criticalEnergy, debug, depthStart, depthStep, detector, e, HDShowerParametrization::e1(), HDShowerParametrization::e2(), EcalHitMaker::ecalHcalGapTotalL0(), EcalHitMaker::ecalHcalGapTotalX0(), HDShowerParametrization::ecalProperties(), EcalHitMaker::ecalTotalL0(), HDShowerParametrization::emax(), HDShowerParametrization::emid(), HDShowerParametrization::emin(), eSpotSize, RandomEngineAndDistribution::flatShoot(), HSParameters::getHDbalanceEH(), HSParameters::getHDcriticalEnergy(), HSParameters::getHDdepthStep(), HSParameters::getHDeSpotSize(), HSParameters::getHDhcalDepthFactor(), HSParameters::getHDlossesOpt(), HSParameters::getHDmaxTRfactor(), HSParameters::getHDnDepthSteps(), HSParameters::getHDnTRsteps(), HSParameters::getHDtransParam(), hcalDepthFactor, HDShowerParametrization::hcalProperties(), EcalHitMaker::hcalTotalL0(), HDShowerParametrization::hsParameters(), mps_fire::i, ECALProperties::interactionLength(), HCALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, dqm-mbProfile::log, LogDebug, lossesOpt, makeSteps(), HLT_2018_cff::maxDepth, maxTRfactor, nDepthSteps, nTRsteps, onEcal, HDShowerParametrization::part1(), HDShowerParametrization::part2(), HDShowerParametrization::r1(), HDShowerParametrization::r2(), HDShowerParametrization::r3(), ECALProperties::radLenIncm(), HCALProperties::radLenIncm(), random, HDShowerParametrization::setCase(), combinedConstraintHelpers::sum2(), tgamEM, tgamHD, theECALproperties, theGrid, theHCALproperties, theParam, theR1, theR2, theR3, transFactor, transParam, x0curr, x0depth, x0EM, and x0HD.

◆ ~HFShower() [1/2]

HFShower::~HFShower ( )
inlinevirtual

Definition at line 38 of file HFShower.h.

38 { ; }

◆ HFShower() [2/2]

HFShower::HFShower ( const std::string &  name,
const HcalDDDSimConstants hcons,
const HcalSimulationParameters hps,
edm::ParameterSet const &  p,
int  chk = 0 
)

Definition at line 21 of file HFShower.cc.

26  : hcalConstant_(hcons), chkFibre_(chk) {
27  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
28  applyFidCut_ = m_HF.getParameter<bool>("ApplyFiducialCut");
29  probMax_ = m_HF.getParameter<double>("ProbMax");
30 
31  edm::LogVerbatim("HFShower") << "HFShower:: Maximum probability cut off " << probMax_ << " Check flag " << chkFibre_;
32 
33  cherenkov_ = std::make_unique<HFCherenkov>(m_HF);
34  fibre_ = std::make_unique<HFFibre>(name, hcalConstant_, hps, p);
35 
36  //Special Geometry parameters
38 }

References applyFidCut_, cherenkov_, chkFibre_, fibre_, HcalDDDSimConstants::getGparHF(), edm::ParameterSet::getParameter(), gpar_, hcalConstant_, Skims_PA_cff::name, AlCaHLTBitMon_ParallelJobs::p, and probMax_.

◆ ~HFShower() [2/2]

virtual HFShower::~HFShower ( )
virtual

Member Function Documentation

◆ compute()

bool HFShower::compute ( )

Compute the shower longitudinal and lateral development.

Definition at line 461 of file HFShower.cc.

461  {
462  // TimeMe theT("FamosHFShower::compute");
463 
464  bool status = false;
465  int numLongit = eStep.size();
466  if (debug)
467  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
468  << " N_long.steps required : " << numLongit << std::endl;
469 
470  if (numLongit > 0) {
471  status = true;
472  // Prepare the trsanverse probability function
473  std::vector<double> Fhist;
474  std::vector<double> rhist;
475  for (int j = 0; j < nTRsteps + 1; j++) {
476  rhist.push_back(maxTRfactor * j / nTRsteps);
477  Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
478  if (debug == 3)
479  LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
480  }
481 
482  //================================================================
483  // Longitudinal steps
484  //================================================================
485  for (int i = 0; i < numLongit; i++) {
486  double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
487  // vary the longitudinal profile if needed
488  if (detector[i] != 1)
489  currentDepthL0 *= hcalDepthFactor;
490  if (debug)
491  LogDebug("FastCalorimetry") << " FamosHFShower::compute - detector = " << detector[i]
492  << " currentDepthL0 = " << currentDepthL0 << std::endl;
493 
494  double maxTRsize = maxTRfactor * rlamStep[i]; // in lambda units
495  double rbinsize = maxTRsize / nTRsteps;
496  double espot = eStep[i] / (double)nspots[i]; // re-adjust espot
497 
498  if (espot > 4. || espot < 0.)
499  LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = " << espot << std::endl;
500 
501  int ecal = 0;
502  if (detector[i] != 1) {
503  bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
504 
505  if (debug)
506  LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of "
507  << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;
508 
509  if (!setHDdepth) {
510  currentDepthL0 -= lamstep[i];
511  setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
512  }
513  if (!setHDdepth)
514  continue;
515 
517  } else {
518  ecal = 1;
519  bool status = theGrid->getPads(currentDepthL0);
520 
521  if (debug)
522  LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = " << status << std::endl;
523 
524  if (!status)
525  continue;
526 
527  theGrid->setSpotEnergy(espot);
528  }
529 
530  //------------------------------------------------------------
531  // Transverse distribution
532  //------------------------------------------------------------
533  int nok = 0; // counter of OK
534  int count = 0;
535  int inf = infinity;
536  if (lossesOpt)
537  inf = nspots[i]; // losses are enabled, otherwise
538  // only OK points are counted ...
539 
540  // Total energy in this layer
541  double eremaining = eStep[i];
542  bool converged = false;
543 
544  while (eremaining > 0. && !converged && count < inf) {
545  ++count;
546 
547  // energy spot (HFL)
548  double newespot = espot;
549 
550  // We need to know a priori if this energy spot if for
551  // a long (1) or short (2) fiber
552 
553  unsigned layer = 1;
554  if (currentDepthL0 < 1.3) // first 22 cm = 1.3 lambda - only HFL
555  layer = 1;
556  else
557  layer = random->flatShoot() < 0.5 ? 1 : 2;
558 
559  if (layer == 2)
560  newespot = 2. * espot;
561 
562  if (eremaining - newespot < 0.)
563  newespot = eremaining;
564 
565  // process transverse distribution
566 
567  double prob = random->flatShoot();
568  int index = indexFinder(prob, Fhist);
569  double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize; // in-bin
570  double phi = 2. * M_PI * random->flatShoot();
571 
572  if (debug == 2)
573  LogDebug("FastCalorimetry") << std::endl
574  << " FamosHFShower::compute "
575  << " r = " << radius << " phi = " << phi << std::endl;
576 
577  // add hit
578 
579  theHcalHitMaker->setSpotEnergy(newespot);
580  theGrid->setSpotEnergy(newespot);
581 
582  bool result;
583  if (ecal) {
584  result = theGrid->addHit(radius, phi, 0); // shouldn't get here !
585 
586  if (debug == 2)
587  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
588  << " theGrid->addHit result = " << result << std::endl;
589 
590  } else {
591  // PV assign espot to long/short fibers
592  result = theHcalHitMaker->addHit(radius, phi, layer);
593 
594  if (debug == 2)
595  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
596  << " theHcalHitMaker->addHit result = " << result << std::endl;
597  }
598 
599  if (result) {
600  ++nok;
601  eremaining -= newespot;
602  if (eremaining <= 0.)
603  converged = true;
604  // std::cout << "Transverse : " << nok << " "
605  // << " , E= " << newespot
606  // << " , Erem = " << eremaining
607  // << std::endl;
608  } else {
609  // std::cout << "WARNING : hit not added" << std::endl;
610  }
611  }
612 
613  // end of tranverse simulation
614  //-----------------------------------------------------
615 
616  if (count == infinity) {
617  status = false;
618  if (debug)
619  LogDebug("FastCalorimetry") << "*** FamosHFShower::compute "
620  << " maximum number of"
621  << " transverse points " << count << " is used !!!" << std::endl;
622  break;
623  }
624 
625  if (debug)
626  LogDebug("FastCalorimetry") << " FamosHFShower::compute "
627  << " long.step No." << i << " Ntry, Nok = " << count << " " << nok << std::endl;
628 
629  } // end of longitudinal steps
630  } // end of no steps
631  return status;
632 }

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), KineDebug3::count(), debug, detector, bsc_activity_cfg::ecal, eStep, RandomEngineAndDistribution::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, mps_fire::i, indexFinder(), dqmiodatasetharvest::inf, infinity, dqmiolumiharvest::j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, TtFullHadEvtBuilder_cfi::prob, CosmicsPD_Skims::radius, random, mps_fire::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), mps_update::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

◆ gam()

double HFShower::gam ( double  x,
double  a 
) const
inlineprivate

Definition at line 45 of file HFShower.h.

45 { return pow(x, a - 1.) * exp(-x); }

References a, JetChargeProducer_cfi::exp, funct::pow(), and x.

Referenced by makeSteps().

◆ getHits() [1/3]

std::vector< HFShower::Hit > HFShower::getHits ( const G4Step *  aStep,
bool  forLibrary 
)

Definition at line 295 of file HFShower.cc.

295  {
296  std::vector<HFShower::Hit> hits;
297  int nHit = 0;
298 
299  double edep = aStep->GetTotalEnergyDeposit();
300  double stepl = 0.;
301 
302  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
303  stepl = aStep->GetStepLength();
304  if ((edep == 0.) || (stepl == 0.)) {
305 #ifdef EDM_ML_DEBUG
306  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
307 #endif
308  return hits;
309  }
310 
311  const G4Track *aTrack = aStep->GetTrack();
312  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
313 
315  double energy = aParticle->GetTotalEnergy();
316  double momentum = aParticle->GetTotalMomentum();
317  double pBeta = momentum / energy;
318  double dose = 0.;
319  int npeDose = 0;
320 
321  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
322  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
323 
324  const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
325  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
326  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
327  double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5 * gpar_[1];
328  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
329  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
330  // @@ Here the depth should be changed (Fibers are all long in Geometry!)
331  int depth = 1;
332  int npmt = 0;
333  bool ok = true;
334  if (zv < 0 || zv > gpar_[1]) {
335  ok = false; // beyond the fiber in Z
336  }
337  if (ok && applyFidCut_) {
338  npmt = HFFibreFiducial::PMTNumber(globalPos);
339  if (npmt <= 0) {
340  ok = false;
341  } else if (npmt > 24) { // a short fibre
342  if (zv > gpar_[0]) {
343  depth = 2;
344  } else {
345  ok = false;
346  }
347  }
348 #ifdef EDM_ML_DEBUG
349  edm::LogVerbatim("HFShower") << "HFShower:getHits(SL): npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":"
350  << gpar_[4] << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
351 #endif
352  } else {
353  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
354  }
355  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
356 
357  double u = localMom.x();
358  double v = localMom.y();
359  double w = localMom.z();
360  double zCoor = localPos.z();
361  double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
362  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
363  double time = fibre_->tShift(localPos, depth, chkFibre_);
364 
365 #ifdef EDM_ML_DEBUG
366  edm::LogVerbatim("HFShower") << "HFShower::getHits(SL): in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
367  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
368  << "\n Direction " << momentumDir << " Local " << localMom;
369 #endif
370  // npe should be 0
371  int npe = 0;
372  if (ok)
373  npe = cherenkov_->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
374  std::vector<double> wavelength = cherenkov_->getWL();
375  std::vector<double> momz = cherenkov_->getMom();
376 
377  for (int i = 0; i < npe; ++i) {
378  hit.time = tSlice + time;
379  hit.wavelength = wavelength[i];
380  hit.momentum = momz[i];
381  hit.position = globalPos;
382  hits.push_back(hit);
383  nHit++;
384  }
385 
386  return hits;
387 }

References funct::abs(), applyFidCut_, cherenkov_, chkFibre_, LEDCalibrationChannels::depth, HCALHighEnergyHPDFilter_cfi::energy, fibre_, gpar_, hfClusterShapes_cfi::hits, mps_fire::i, Skims_PA_cff::name, convertSQLiteXML::ok, HFFibreFiducial::PMTNumber(), ntuplemaker::time, findQualityFiles::v, and w.

◆ getHits() [2/3]

std::vector< HFShower::Hit > HFShower::getHits ( const G4Step *  aStep,
bool  forLibraryProducer,
double  zoffset 
)

Definition at line 170 of file HFShower.cc.

170  {
171  std::vector<HFShower::Hit> hits;
172  int nHit = 0;
173 
174  double edep = aStep->GetTotalEnergyDeposit();
175  double stepl = 0.;
176 
177  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
178  stepl = aStep->GetStepLength();
179  if ((edep == 0.) || (stepl == 0.)) {
180 #ifdef EDM_ML_DEBUG
181  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
182 #endif
183  return hits;
184  }
185  const G4Track *aTrack = aStep->GetTrack();
186  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
187 
189  double energy = aParticle->GetTotalEnergy();
190  double momentum = aParticle->GetTotalMomentum();
191  double pBeta = momentum / energy;
192  double dose = 0.;
193  int npeDose = 0;
194 
195  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
196  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
197 
198  G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
199  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
200  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
201  //double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5*gpar_[1];
202  //double zv = std::abs(globalPos.z()) - gpar_[4];
203  double zv = gpar_[1] - (std::abs(globalPos.z()) - zoffset);
204  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
205  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
206  // @@ Here the depth should be changed (Fibers all long in Geometry!)
207  int depth = 1;
208  int npmt = 0;
209  bool ok = true;
210  if (zv < 0. || zv > gpar_[1]) {
211  ok = false; // beyond the fiber in Z
212  }
213  if (ok && applyFidCut_) {
214  npmt = HFFibreFiducial::PMTNumber(globalPos);
215  if (npmt <= 0) {
216  ok = false;
217  } else if (npmt > 24) { // a short fibre
218  if (zv > gpar_[0]) {
219  depth = 2;
220  } else {
221  ok = false;
222  }
223  }
224 #ifdef EDM_ML_DEBUG
225  edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar_[4]
226  << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
227 #endif
228  } else {
229  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
230  }
231  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
232 
233  double u = localMom.x();
234  double v = localMom.y();
235  double w = localMom.z();
236  double zCoor = localPos.z();
237  double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
238  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
239  double time = fibre_->tShift(localPos, depth, chkFibre_);
240 
241 #ifdef EDM_ML_DEBUG
242  edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
243  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
244  << "\n Direction " << momentumDir << " Local " << localMom;
245 #endif
246  // Here npe should be 0 if there is no fiber (@@ M.K.)
247  int npe = 1;
248  std::vector<double> wavelength;
249  std::vector<double> momz;
250  if (!applyFidCut_) { // _____ Tmp close of the cherenkov function
251  if (ok)
252  npe = cherenkov_->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
253  wavelength = cherenkov_->getWL();
254  momz = cherenkov_->getMom();
255  } // ^^^^^ End of Tmp close of the cherenkov function
256  if (ok && npe > 0) {
257  for (int i = 0; i < npe; ++i) {
258  double p = 1.;
259  if (!applyFidCut_)
260  p = fibre_->attLength(wavelength[i]);
261  double r1 = G4UniformRand();
262  double r2 = G4UniformRand();
263 #ifdef EDM_ML_DEBUG
264  edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
265  << " r2 " << r2 << ":" << probMax_
266  << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax_);
267 #endif
268  if (applyFidCut_ || chkFibre_ < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax_)) {
269  hit.depth = depth;
270  hit.time = tSlice + time;
271  if (!applyFidCut_) {
272  hit.wavelength = wavelength[i]; // Tmp
273  hit.momentum = momz[i]; // Tmp
274  } else {
275  hit.wavelength = 300.; // Tmp
276  hit.momentum = 1.; // Tmp
277  }
278  hit.position = globalPos;
279  hits.push_back(hit);
280  nHit++;
281  }
282  }
283  }
284 
285 #ifdef EDM_ML_DEBUG
286  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
287  for (int i = 0; i < nHit; ++i)
288  edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
289  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
290  << hits[i].position;
291 #endif
292  return hits;
293 }

References funct::abs(), applyFidCut_, cherenkov_, chkFibre_, LEDCalibrationChannels::depth, HCALHighEnergyHPDFilter_cfi::energy, JetChargeProducer_cfi::exp, fibre_, gpar_, hfClusterShapes_cfi::hits, mps_fire::i, Skims_PA_cff::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), probMax_, diffTwoXMLs::r1, diffTwoXMLs::r2, ntuplemaker::time, findQualityFiles::v, and w.

◆ getHits() [3/3]

std::vector< HFShower::Hit > HFShower::getHits ( const G4Step *  aStep,
double  weight 
)

Definition at line 42 of file HFShower.cc.

42  {
43  std::vector<HFShower::Hit> hits;
44  int nHit = 0;
45 
46  double edep = weight * (aStep->GetTotalEnergyDeposit());
47 #ifdef EDM_ML_DEBUG
48  edm::LogVerbatim("HFShower") << "HFShower::getHits: energy " << aStep->GetTotalEnergyDeposit() << " weight " << weight
49  << " edep " << edep;
50 #endif
51  double stepl = 0.;
52 
53  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
54  stepl = aStep->GetStepLength();
55  if ((edep == 0.) || (stepl == 0.)) {
56 #ifdef EDM_ML_DEBUG
57  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
58 #endif
59  return hits;
60  }
61  const G4Track *aTrack = aStep->GetTrack();
62  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
63 
65  double energy = aParticle->GetTotalEnergy();
66  double momentum = aParticle->GetTotalMomentum();
67  double pBeta = momentum / energy;
68  double dose = 0.;
69  int npeDose = 0;
70 
71  const G4ThreeVector &momentumDir = aParticle->GetMomentumDirection();
72  const G4ParticleDefinition *particleDef = aTrack->GetDefinition();
73 
74  const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
75  const G4ThreeVector &globalPos = preStepPoint->GetPosition();
76  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
77  //double zv = std::abs(globalPos.z()) - gpar_[4] - 0.5*gpar_[1];
78  double zv = std::abs(globalPos.z()) - gpar_[4];
79  G4ThreeVector localPos = G4ThreeVector(globalPos.x(), globalPos.y(), zv);
80  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(momentumDir);
81  // @@ Here the depth should be changed (Fibers all long in Geometry!)
82  int depth = 1;
83  int npmt = 0;
84  bool ok = true;
85  if (zv < 0. || zv > gpar_[1]) {
86  ok = false; // beyond the fiber in Z
87  }
88  if (ok && applyFidCut_) {
89  npmt = HFFibreFiducial::PMTNumber(globalPos);
90  if (npmt <= 0) {
91  ok = false;
92  } else if (npmt > 24) { // a short fibre
93  if (zv > gpar_[0]) {
94  depth = 2;
95  } else {
96  ok = false;
97  }
98  }
99 #ifdef EDM_ML_DEBUG
100  edm::LogVerbatim("HFShower") << "HFShower: npmt " << npmt << " zv " << std::abs(globalPos.z()) << ":" << gpar_[4]
101  << ":" << zv << ":" << gpar_[0] << " ok " << ok << " depth " << depth;
102 #endif
103  } else {
104  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0)) % 10; // All LONG!
105  }
106  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
107 
108  double u = localMom.x();
109  double v = localMom.y();
110  double w = localMom.z();
111  double zCoor = localPos.z();
112  double zFibre = (0.5 * gpar_[1] - zCoor - translation.z());
113  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
114  double time = fibre_->tShift(localPos, depth, chkFibre_);
115 
116 #ifdef EDM_ML_DEBUG
117  edm::LogVerbatim("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor << "(" << globalPos.z() << ") "
118  << zFibre << " Trans " << translation << " Time " << tSlice << " " << time
119  << "\n Direction " << momentumDir << " Local " << localMom;
120 #endif
121  // Here npe should be 0 if there is no fiber (@@ M.K.)
122  int npe = 1;
123  std::vector<double> wavelength;
124  std::vector<double> momz;
125  if (!applyFidCut_) { // _____ Tmp close of the cherenkov function
126  if (ok)
127  npe = cherenkov_->computeNPE(aStep, particleDef, pBeta, u, v, w, stepl, zFibre, dose, npeDose);
128  wavelength = cherenkov_->getWL();
129  momz = cherenkov_->getMom();
130  } // ^^^^^ End of Tmp close of the cherenkov function
131  if (ok && npe > 0) {
132  for (int i = 0; i < npe; ++i) {
133  double p = 1.;
134  if (!applyFidCut_)
135  p = fibre_->attLength(wavelength[i]);
136  double r1 = G4UniformRand();
137  double r2 = G4UniformRand();
138 #ifdef EDM_ML_DEBUG
139  edm::LogVerbatim("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1 << ":" << exp(-p * zFibre)
140  << " r2 " << r2 << ":" << probMax_
141  << " Survive: " << (r1 <= exp(-p * zFibre) && r2 <= probMax_);
142 #endif
143  if (applyFidCut_ || chkFibre_ < 0 || (r1 <= exp(-p * zFibre) && r2 <= probMax_)) {
144  hit.depth = depth;
145  hit.time = tSlice + time;
146  if (!applyFidCut_) {
147  hit.wavelength = wavelength[i]; // Tmp
148  hit.momentum = momz[i]; // Tmp
149  } else {
150  hit.wavelength = 300.; // Tmp
151  hit.momentum = 1.; // Tmp
152  }
153  hit.position = globalPos;
154  hits.push_back(hit);
155  nHit++;
156  }
157  }
158  }
159 
160 #ifdef EDM_ML_DEBUG
161  edm::LogVerbatim("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
162  for (int i = 0; i < nHit; ++i)
163  edm::LogVerbatim("HFShower") << "HFShower::Hit " << i << " WaveLength " << hits[i].wavelength << " Time "
164  << hits[i].time << " Momentum " << hits[i].momentum << " Position "
165  << hits[i].position;
166 #endif
167  return hits;
168 }

References funct::abs(), applyFidCut_, cherenkov_, chkFibre_, LEDCalibrationChannels::depth, HCALHighEnergyHPDFilter_cfi::energy, JetChargeProducer_cfi::exp, fibre_, gpar_, hfClusterShapes_cfi::hits, mps_fire::i, Skims_PA_cff::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), probMax_, diffTwoXMLs::r1, diffTwoXMLs::r2, ntuplemaker::time, findQualityFiles::v, and w.

Referenced by FiberSD::ProcessHits().

◆ indexFinder()

int HFShower::indexFinder ( double  x,
const std::vector< double > &  Fhist 
)
private

Definition at line 634 of file HFShower.cc.

634  {
635  // binary search in the vector of doubles
636  int size = Fhist.size();
637 
638  int curr = size / 2;
639  int step = size / 4;
640  int iter;
641  int prevdir = 0;
642  int actudir = 0;
643 
644  for (iter = 0; iter < size; iter++) {
645  if (curr >= size || curr < 1)
646  LogWarning("FastCalorimetry") << " FamosHFShower::indexFinder - wrong current index = " << curr << " !!!"
647  << std::endl;
648 
649  if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
650  break;
651  prevdir = actudir;
652  if (x > Fhist[curr]) {
653  actudir = 1;
654  } else {
655  actudir = -1;
656  }
657  if (prevdir * actudir < 0) {
658  if (step > 1)
659  step /= 2;
660  }
661  curr += actudir * step;
662  if (curr > size)
663  curr = size;
664  else {
665  if (curr < 1) {
666  curr = 1;
667  }
668  }
669 
670  if (debug == 3)
671  LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
672  << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
673  }
674 
675  if (debug == 3)
676  LogDebug("FastCalorimetry") << " indexFinder x = " << x << " found index = " << curr - 1 << std::endl;
677 
678  return curr - 1;
679 }

References debug, LogDebug, findQualityFiles::size, and x.

Referenced by compute().

◆ makeSteps()

void HFShower::makeSteps ( int  nsteps)
private

Definition at line 363 of file HFShower.cc.

363  {
364  double sumes = 0.;
365  double sum = 0.;
366  std::vector<double> temp;
367 
368  if (debug)
369  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
370  << " nsteps required : " << nsteps << std::endl;
371 
372  int count = 0;
373  for (int i = 0; i < nsteps; i++) {
374  double deplam = lamdepth[i] - 0.5 * lamstep[i];
375  double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i];
376  double x = betEM * depx0;
377  double y = betHD * deplam;
378 
379  if (debug == 2)
380  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps "
381  << " - step " << i << " depx0, x = " << depx0 << ", " << x
382  << " deplam, y = " << deplam << ", " << y << std::endl;
383 
384  double est = (part * betEM * gam(x, alpEM) * lamcurr[i] / (x0curr[i] * tgamEM) +
385  (1. - part) * betHD * gam(y, alpHD) / tgamHD) *
386  lamstep[i];
387 
388  // protection ...
389  if (est < 0.) {
390  LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps "
391  << " - negative step energy !!!" << std::endl;
392  est = 0.;
393  break;
394  }
395 
396  // for estimates only
397  sum += est;
398  int nPest = (int)(est * e / sum / eSpotSize);
399 
400  if (debug == 2)
401  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - nPoints estimate = " << nPest << std::endl;
402 
403  if (nPest <= 1 && count != 0)
404  break;
405 
406  // good step - to proceed
407 
408  temp.push_back(est);
409  sumes += est;
410 
411  rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor);
412  count++;
413  }
414 
415  // fluctuations in ECAL and re-distribution of remaining energy in HCAL
416  if (detector[0] == 1 && count > 1) {
417  double oldECALenergy = temp[0];
418  double oldHCALenergy = sumes - oldECALenergy;
419  double newECALenergy = 2. * sumes;
420  for (int i = 0; newECALenergy > sumes && i < infinity; i++)
421  newECALenergy = 2. * balanceEH * random->flatShoot() * oldECALenergy;
422 
423  if (debug == 2)
424  LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps "
425  << " ECAL fraction : old/new - " << oldECALenergy / sumes << "/"
426  << newECALenergy / sumes << std::endl;
427 
428  temp[0] = newECALenergy;
429  double newHCALenergy = sumes - newECALenergy;
430  double newHCALreweight = newHCALenergy / oldHCALenergy;
431 
432  for (int i = 1; i < count; i++) {
433  temp[i] *= newHCALreweight;
434  }
435  }
436 
437  // final re-normalization of the energy fractions
438  double etot = 0.;
439  for (int i = 0; i < count; i++) {
440  eStep.push_back(temp[i] * e / sumes);
441  nspots.push_back((int)(eStep[i] / eSpotSize) + 1);
442  etot += eStep[i];
443 
444  if (debug)
445  LogDebug("FastCalorimetry") << i << " xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
446  << " Estep func = " << eStep[i] << " Rstep = " << rlamStep[i]
447  << " Nspots = " << nspots[i] << std::endl;
448  }
449 
450  // The only step is in ECAL - let's make the size bigger ...
451  if (count == 1 and detector[0] == 1)
452  rlamStep[0] *= 2.;
453 
454  if (debug) {
455  if (eStep[0] > 0.95 * e && detector[0] == 1)
456  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
457  << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
458  }
459 }

References aloge, alpEM, alpHD, balanceEH, betEM, betHD, KineDebug3::count(), debug, detector, e, eSpotSize, eStep, RandomEngineAndDistribution::flatShoot(), gam(), mps_fire::i, infinity, createfilelist::int, lamcurr, lamdepth, lamstep, LogDebug, nspots, random, rlamStep, groupFilesInBlocks::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, x, x0curr, x0depth, and y.

Referenced by HFShower().

◆ transProb()

double HFShower::transProb ( double  factor,
double  R,
double  r 
)
inlineprivate

Definition at line 50 of file HFShower.h.

50  {
51  double fsq = factor * factor;
52  return ((fsq + 1.) / fsq) * r * r / (r * r + R * R);
53  }

References DQMScaleToClient_cfi::factor, dttmaxenums::R, and alignCSCRings::r.

Referenced by compute().

Member Data Documentation

◆ aloge

double HFShower::aloge
private

Definition at line 73 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ alpEM

double HFShower::alpEM
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ alpHD

double HFShower::alpHD
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ applyFidCut_

bool HFShower::applyFidCut_
private

Definition at line 52 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ balanceEH

double HFShower::balanceEH
private

Definition at line 113 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ betEM

double HFShower::betEM
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ betHD

double HFShower::betHD
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ cherenkov_

std::unique_ptr<HFCherenkov> HFShower::cherenkov_
private

Definition at line 48 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ chkFibre_

int HFShower::chkFibre_
private

Definition at line 51 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ criticalEnergy

double HFShower::criticalEnergy
private

Definition at line 109 of file HFShower.h.

Referenced by HFShower().

◆ depthStart

double HFShower::depthStart
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower().

◆ depthStep

double HFShower::depthStep
private

Definition at line 107 of file HFShower.h.

Referenced by HFShower().

◆ detector

std::vector<int> HFShower::detector
private

Definition at line 75 of file HFShower.h.

Referenced by compute(), HFShower(), and makeSteps().

◆ e

double HFShower::e
private

Definition at line 92 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ eSpotSize

double HFShower::eSpotSize
private

Definition at line 105 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ eStep

std::vector<double> HFShower::eStep
private

Definition at line 76 of file HFShower.h.

Referenced by compute(), and makeSteps().

◆ fibre_

std::unique_ptr<HFFibre> HFShower::fibre_
private

Definition at line 49 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ gpar_

std::vector<double> HFShower::gpar_
private

Definition at line 54 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ hcalConstant_

const HcalDDDSimConstants* HFShower::hcalConstant_
private

Definition at line 46 of file HFShower.h.

Referenced by HFShower().

◆ hcalDepthFactor

double HFShower::hcalDepthFactor
private

Definition at line 115 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ infinity

int HFShower::infinity
private

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

◆ lambdaEM

double HFShower::lambdaEM
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower().

◆ lambdaHD

double HFShower::lambdaHD
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower().

◆ lamcurr

std::vector<double> HFShower::lamcurr
private

Definition at line 78 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ lamdepth

std::vector<double> HFShower::lamdepth
private

Definition at line 78 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ lamstep

std::vector<double> HFShower::lamstep
private

Definition at line 78 of file HFShower.h.

Referenced by compute(), HFShower(), and makeSteps().

◆ lamtotal

std::vector<double> HFShower::lamtotal
private

Definition at line 78 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ lossesOpt

int HFShower::lossesOpt
private

Definition at line 95 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ maxTRfactor

double HFShower::maxTRfactor
private

Definition at line 111 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ nDepthSteps

int HFShower::nDepthSteps
private

Definition at line 97 of file HFShower.h.

Referenced by HFShower().

◆ nspots

std::vector<int> HFShower::nspots
private

Definition at line 75 of file HFShower.h.

Referenced by compute(), and makeSteps().

◆ nTRsteps

int HFShower::nTRsteps
private

Definition at line 99 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ onEcal

int HFShower::onEcal
private

Definition at line 89 of file HFShower.h.

Referenced by HFShower().

◆ part

double HFShower::part
private

Definition at line 68 of file HFShower.h.

◆ probMax_

double HFShower::probMax_
private

Definition at line 53 of file HFShower.h.

Referenced by getHits(), and HFShower().

◆ random

const RandomEngineAndDistribution* HFShower::random
private

Definition at line 118 of file HFShower.h.

Referenced by compute(), HFShower(), and makeSteps().

◆ rlamStep

std::vector<double> HFShower::rlamStep
private

Definition at line 76 of file HFShower.h.

Referenced by compute(), and makeSteps().

◆ tgamEM

double HFShower::tgamEM
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ tgamHD

double HFShower::tgamHD
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ theECALproperties

const ECALProperties* HFShower::theECALproperties
private

Definition at line 63 of file HFShower.h.

Referenced by HFShower().

◆ theGrid

EcalHitMaker* HFShower::theGrid
private

Definition at line 83 of file HFShower.h.

Referenced by compute(), and HFShower().

◆ theHcalHitMaker

HcalHitMaker* HFShower::theHcalHitMaker
private

Definition at line 86 of file HFShower.h.

Referenced by compute().

◆ theHCALproperties

const HCALProperties* HFShower::theHCALproperties
private

Definition at line 64 of file HFShower.h.

Referenced by HFShower().

◆ theParam

HDShowerParametrization* HFShower::theParam
private

Definition at line 60 of file HFShower.h.

Referenced by HFShower().

◆ theR1

double HFShower::theR1
private

Definition at line 67 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ theR2

double HFShower::theR2
private

Definition at line 67 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ theR3

double HFShower::theR3
private

Definition at line 67 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ transFactor

double HFShower::transFactor
private

Definition at line 103 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ transParam

double HFShower::transParam
private

Definition at line 101 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ x0curr

std::vector<double> HFShower::x0curr
private

Definition at line 77 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ x0depth

std::vector<double> HFShower::x0depth
private

Definition at line 77 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

◆ x0EM

double HFShower::x0EM
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower().

◆ x0HD

double HFShower::x0HD
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower().

HFShower::theR2
double theR2
Definition: HFShower.h:67
HFShower::depthStep
double depthStep
Definition: HFShower.h:107
HFShower::theR3
double theR3
Definition: HFShower.h:67
HSParameters::getHDnDepthSteps
int getHDnDepthSteps() const
Definition: HSParameters.h:20
DDAxes::y
HFShower::eSpotSize
double eSpotSize
Definition: HFShower.h:105
mps_fire.i
i
Definition: mps_fire.py:355
HFShower::theR1
double theR1
Definition: HFShower.h:67
HFShower::tgamHD
double tgamHD
Definition: HFShower.h:68
HFShower::cherenkov_
std::unique_ptr< HFCherenkov > cherenkov_
Definition: HFShower.h:48
HFShower::detector
std::vector< int > detector
Definition: HFShower.h:75
HFShower::aloge
double aloge
Definition: HFShower.h:73
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
HFShower::applyFidCut_
bool applyFidCut_
Definition: HFShower.h:52
step
step
Definition: StallMonitor.cc:94
HSParameters::getHDmaxTRfactor
double getHDmaxTRfactor() const
Definition: HSParameters.h:26
mps_update.status
status
Definition: mps_update.py:69
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HCALProperties::radLenIncm
double radLenIncm() const override
Radiation length in cm.
Definition: HCALProperties.h:37
HFShower::gpar_
std::vector< double > gpar_
Definition: HFShower.h:54
edm::LogInfo
Definition: MessageLogger.h:254
HDShowerParametrization::beth1
double beth1() const
Definition: HDShowerParametrization.h:81
HcalDDDSimConstants::getGparHF
const std::vector< double > & getGparHF() const
Definition: HcalDDDSimConstants.h:45
HFShower::lambdaHD
double lambdaHD
Definition: HFShower.h:71
HDShowerParametrization::bete2
double bete2() const
Definition: HDShowerParametrization.h:63
HFShower::nTRsteps
int nTRsteps
Definition: HFShower.h:99
DDAxes::x
HFShower::theHcalHitMaker
HcalHitMaker * theHcalHitMaker
Definition: HFShower.h:86
HSParameters::getHDdepthStep
double getHDdepthStep() const
Definition: HSParameters.h:24
findQualityFiles.v
v
Definition: findQualityFiles.py:179
HDShowerParametrization::part1
double part1() const
Definition: HDShowerParametrization.h:93
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
HDShowerParametrization::r2
double r2() const
Definition: HDShowerParametrization.h:106
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
EcalHitMaker::getPads
bool getPads(double depth, bool inCm=false)
Definition: EcalHitMaker.cc:820
HFFibreFiducial::PMTNumber
int PMTNumber(const G4ThreeVector &pe_effect)
Definition: HFFibreFiducial.cc:9
HDShowerParametrization::alph1
double alph1() const
Definition: HDShowerParametrization.h:69
HFShower::maxTRfactor
double maxTRfactor
Definition: HFShower.h:111
HLT_2018_cff.maxDepth
maxDepth
Definition: HLT_2018_cff.py:7356
EcalHitMaker::addHit
bool addHit(double r, double phi, unsigned layer=0) override
Definition: EcalHitMaker.cc:184
HFShower::lambdaEM
double lambdaEM
Definition: HFShower.h:71
HSParameters::getHDtransParam
double getHDtransParam() const
Definition: HSParameters.h:22
HFShower::nDepthSteps
int nDepthSteps
Definition: HFShower.h:97
HDShowerParametrization::emax
double emax() const
Definition: HDShowerParametrization.h:41
part
part
Definition: HCALResponse.h:20
HFShower::criticalEnergy
double criticalEnergy
Definition: HFShower.h:109
HFShower::eStep
std::vector< double > eStep
Definition: HFShower.h:76
HFShower::depthStart
double depthStart
Definition: HFShower.h:72
HDShowerParametrization::r1
double r1() const
Definition: HDShowerParametrization.h:105
w
const double w
Definition: UKUtility.cc:23
HDShowerParametrization::part2
double part2() const
Definition: HDShowerParametrization.h:99
HCALProperties::interactionLength
double interactionLength() const override
Interaction length in cm.
Definition: HCALProperties.h:56
HFShower::chkFibre_
int chkFibre_
Definition: HFShower.h:51
EcalHitMaker::ecalHcalGapTotalL0
double ecalHcalGapTotalL0() const
ECAL-HCAL transition.
Definition: EcalHitMaker.h:91
HFShower::infinity
int infinity
Definition: HFShower.h:80
DQMScaleToClient_cfi.factor
factor
Definition: DQMScaleToClient_cfi.py:8
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
HFShower::hcalConstant_
const HcalDDDSimConstants * hcalConstant_
Definition: HFShower.h:46
HDShowerParametrization::setCase
void setCase(int choice)
Definition: HDShowerParametrization.h:29
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
HFShower::transParam
double transParam
Definition: HFShower.h:101
edm::LogWarning
Definition: MessageLogger.h:141
HFShower::lamstep
std::vector< double > lamstep
Definition: HFShower.h:78
HSParameters::getHDbalanceEH
double getHDbalanceEH() const
Definition: HSParameters.h:27
HFShower::indexFinder
int indexFinder(double x, const std::vector< double > &Fhist)
Definition: HFShower.cc:634
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:670
HSParameters::getHDlossesOpt
int getHDlossesOpt() const
Definition: HSParameters.h:19
edm::ParameterSet
Definition: ParameterSet.h:36
a
double a
Definition: hdecay.h:119
HFShower::balanceEH
double balanceEH
Definition: HFShower.h:113
HFShower::random
const RandomEngineAndDistribution * random
Definition: HFShower.h:118
KineDebug3::count
void count()
Definition: KinematicConstrainedVertexUpdatorT.h:21
HSParameters::getHDhcalDepthFactor
double getHDhcalDepthFactor() const
Definition: HSParameters.h:28
EcalHitMaker::setSpotEnergy
void setSpotEnergy(double e) override
Definition: EcalHitMaker.h:112
diffTwoXMLs.r2
r2
Definition: diffTwoXMLs.py:73
HDShowerParametrization::hsParameters
const HSParameters * hsParameters() const
Definition: HDShowerParametrization.h:26
HFShower::tgamEM
double tgamEM
Definition: HFShower.h:68
createfilelist.int
int
Definition: createfilelist.py:10
HFShower::x0curr
std::vector< double > x0curr
Definition: HFShower.h:77
HDShowerParametrization::e1
double e1() const
Definition: HDShowerParametrization.h:43
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:50
edm::LogVerbatim
Definition: MessageLogger.h:297
HFShower::theGrid
EcalHitMaker * theGrid
Definition: HFShower.h:83
HFShower::transFactor
double transFactor
Definition: HFShower.h:103
HFShower::lamtotal
std::vector< double > lamtotal
Definition: HFShower.h:78
combinedConstraintHelpers::sum2
void sum2(T &x, T y)
Definition: CombinedKinematicConstraintT.h:74
HFShower::gam
double gam(double x, double a) const
Definition: HFShower.h:45
HDShowerParametrization::hcalProperties
const HCALProperties * hcalProperties() const
Definition: HDShowerParametrization.h:24
RandomEngineAndDistribution::flatShoot
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngineAndDistribution.h:27
HcalHitMaker::setDepth
bool setDepth(double, bool inCm=false)
set the depth in X0 or Lambda0 units depending on showerType
Definition: HcalHitMaker.cc:114
HFShower::x0HD
double x0HD
Definition: HFShower.h:71
alignCSCRings.r
r
Definition: alignCSCRings.py:93
HDShowerParametrization::alph2
double alph2() const
Definition: HDShowerParametrization.h:75
HFShower::lossesOpt
int lossesOpt
Definition: HFShower.h:95
HDShowerParametrization::alpe1
double alpe1() const
Definition: HDShowerParametrization.h:45
DDAxes::phi
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HDShowerParametrization::alpe2
double alpe2() const
Definition: HDShowerParametrization.h:51
HFShower::onEcal
int onEcal
Definition: HFShower.h:89
HSParameters::getHDcriticalEnergy
double getHDcriticalEnergy() const
Definition: HSParameters.h:25
ECALProperties::radLenIncm
double radLenIncm() const override
Radiation length in cm.
Definition: ECALProperties.h:32
HDShowerParametrization::ecalProperties
const ECALProperties * ecalProperties() const
Definition: HDShowerParametrization.h:22
HFShower::alpHD
double alpHD
Definition: HFShower.h:68
dqmiodatasetharvest.inf
inf
Definition: dqmiodatasetharvest.py:38
debug
#define debug
Definition: HFShower.cc:25
HFShower::x0depth
std::vector< double > x0depth
Definition: HFShower.h:77
diffTwoXMLs.r1
r1
Definition: diffTwoXMLs.py:53
HFShower::fibre_
std::unique_ptr< HFFibre > fibre_
Definition: HFShower.h:49
HFShower::theECALproperties
const ECALProperties * theECALproperties
Definition: HFShower.h:63
HFShower::alpEM
double alpEM
Definition: HFShower.h:68
HSParameters::getHDnTRsteps
int getHDnTRsteps() const
Definition: HSParameters.h:21
EcalHitMaker::ecalHcalGapTotalX0
double ecalHcalGapTotalX0() const
ECAL-HCAL transition.
Definition: EcalHitMaker.h:70
CosmicsPD_Skims.radius
radius
Definition: CosmicsPD_Skims.py:135
HDShowerParametrization::bete1
double bete1() const
Definition: HDShowerParametrization.h:57
HFShower::lamdepth
std::vector< double > lamdepth
Definition: HFShower.h:78
HFShower::nspots
std::vector< int > nspots
Definition: HFShower.h:75
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
HFShower::probMax_
double probMax_
Definition: HFShower.h:53
HFShower::betHD
double betHD
Definition: HFShower.h:68
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
bsc_activity_cfg.ecal
ecal
Definition: bsc_activity_cfg.py:25
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
ECALProperties::interactionLength
double interactionLength() const override
Interaction length in cm: 18.5 for Standard ECAL.
Definition: ECALProperties.h:49
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
HDShowerParametrization::e2
double e2() const
Definition: HDShowerParametrization.h:44
mps_fire.result
result
Definition: mps_fire.py:303
HcalHitMaker::addHit
bool addHit(double r, double phi, unsigned layer=0) override
add the hit in the HCAL in local coordinates
Definition: HcalHitMaker.cc:30
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HFShower::lamcurr
std::vector< double > lamcurr
Definition: HFShower.h:78
HFShower::betEM
double betEM
Definition: HFShower.h:68
HDShowerParametrization::emin
double emin() const
Definition: HDShowerParametrization.h:37
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
EcalHitMaker::ecalTotalL0
double ecalTotalL0() const
in the ECAL
Definition: EcalHitMaker.h:85
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
HSParameters::getHDeSpotSize
double getHDeSpotSize() const
Definition: HSParameters.h:23
HFShower::makeSteps
void makeSteps(int nsteps)
Definition: HFShower.cc:363
HFShower::x0EM
double x0EM
Definition: HFShower.h:71
ntuplemaker.time
time
Definition: ntuplemaker.py:310
HDShowerParametrization::emid
double emid() const
Definition: HDShowerParametrization.h:39
HFShower::rlamStep
std::vector< double > rlamStep
Definition: HFShower.h:76
EcalHitMaker::hcalTotalL0
double hcalTotalL0() const
in the HCAL
Definition: EcalHitMaker.h:88
HFShower::theParam
HDShowerParametrization * theParam
Definition: HFShower.h:60
HcalHitMaker::setSpotEnergy
void setSpotEnergy(double e) override
Set the spot energy.
Definition: HcalHitMaker.h:26
HFShower::hcalDepthFactor
double hcalDepthFactor
Definition: HFShower.h:115
dttmaxenums::R
Definition: DTTMax.h:29
HDShowerParametrization::r3
double r3() const
Definition: HDShowerParametrization.h:107
HFShower::Hit
Definition: HFShower.h:32
TtFullHadEvtBuilder_cfi.prob
prob
Definition: TtFullHadEvtBuilder_cfi.py:33
HFShower::theHCALproperties
const HCALProperties * theHCALproperties
Definition: HFShower.h:64
weight
Definition: weight.py:1
HFShower::transProb
double transProb(double factor, double R, double r)
Definition: HFShower.h:50
hit
Definition: SiStripHitEffFromCalibTree.cc:88
HFShower::e
double e
Definition: HFShower.h:92
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
HDShowerParametrization::beth2
double beth2() const
Definition: HDShowerParametrization.h:87