CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 (G4Step *aStep)
 
std::vector< HitgetHits (G4Step *aStep, bool forLibrary)
 
 HFShower (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk=0)
 
 HFShower (const RandomEngine *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
 
void initRun (G4ParticleTable *)
 
virtual ~HFShower ()
 
virtual ~HFShower ()
 

Private Member Functions

double gam (double x, double a) const
 
std::vector< double > getDDDArray (const std::string &, const DDsvalues_type &, int &)
 
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
 
HFCherenkovcherenkov
 
int chkFibre
 
double criticalEnergy
 
double depthStart
 
double depthStep
 
std::vector< int > detector
 
double e
 
double eSpotSize
 
std::vector< double > eStep
 
HFFibrefibre
 
std::vector< double > gpar
 
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 RandomEnginerandom
 
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

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

Definition at line 29 of file HFShower.h.

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

Definition at line 30 of file HFShower.h.

typedef Steps::const_iterator HFShower::step_iterator

Definition at line 32 of file HFShower.h.

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

Definition at line 31 of file HFShower.h.

Definition at line 27 of file HFShower.h.

Constructor & Destructor Documentation

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

Definition at line 30 of file HFShower.cc.

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, RandomEngine::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(), i, ECALProperties::interactionLength(), HCALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, create_public_lumi_plots::log, LogDebug, lossesOpt, makeSteps(), maxTRfactor, nDepthSteps, nTRsteps, onEcal, HDShowerParametrization::part1(), HDShowerParametrization::part2(), HDShowerParametrization::r1(), HDShowerParametrization::r2(), HDShowerParametrization::r3(), ECALProperties::radLenIncm(), HCALProperties::radLenIncm(), random, HDShowerParametrization::setCase(), tgamEM, tgamHD, theECALproperties, theGrid, theHCALproperties, theParam, theR1, theR2, theR3, transFactor, transParam, x0curr, x0depth, x0EM, and x0HD.

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

Definition at line 41 of file HFShower.h.

41 {;}
HFShower::HFShower ( std::string &  name,
const DDCompactView cpv,
edm::ParameterSet const &  p,
int  chk = 0 
)

Definition at line 25 of file HFShower.cc.

References DDFilteredView::addFilter(), applyFidCut, cherenkov, chkFibre, DDSpecificsFilter::equals, edm::hlt::Exception, fibre, alcazmumu_cfi::filter, DDFilteredView::firstChild(), getDDDArray(), edm::ParameterSet::getParameter(), gpar, DDFilteredView::mergedSpecifics(), mergeVDriftHistosByStation::name, probMax, DDSpecificsFilter::setCriteria(), and relativeConstraints::value.

26  : cherenkov(0),
27  fibre(0),
28  chkFibre(chk) {
29 
30  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
31  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
32  probMax = m_HF.getParameter<double>("ProbMax");
33 
34  edm::LogInfo("HFShower") << "HFShower:: Maximum probability cut off "
35  << probMax << " Check flag " << chkFibre;
36 
37  G4String attribute = "ReadOutName";
38  G4String value = name;
40  DDValue ddv(attribute,value,0);
42  DDFilteredView fv(cpv);
43  fv.addFilter(filter);
44  bool dodet = fv.firstChild();
45  if ( dodet ) {
46  DDsvalues_type sv(fv.mergedSpecifics());
47 
48  //Special Geometry parameters
49  int ngpar = 7;
50  gpar = getDDDArray("gparHF", sv,ngpar);
51  edm::LogInfo("HFShower") << "HFShower: " << ngpar << " gpar (cm)";
52  for (int ig=0; ig<ngpar; ig++)
53  edm::LogInfo("HFShower") << "HFShower: gpar[" << ig << "] = "
54  << gpar[ig]/cm << " cm";
55  } else {
56  edm::LogError("HFShower") << "HFShower: cannot get filtered "
57  << " view for " << attribute << " matching " << name;
58  throw cms::Exception("Unknown", "HFShower")
59  << "cannot match " << attribute << " to " << name <<"\n";
60  }
61  cherenkov = new HFCherenkov(m_HF);
62  fibre = new HFFibre(name, cpv, p);
63 }
T getParameter(std::string const &) const
HFFibre * fibre
Definition: HFShower.h:52
bool applyFidCut
Definition: HFShower.h:47
HFCherenkov * cherenkov
Definition: HFShower.h:51
int chkFibre
Definition: HFShower.h:54
double probMax
Definition: HFShower.h:55
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
std::vector< double > gpar
Definition: HFShower.h:56
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
Definition: HFShower.cc:301
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37
virtual HFShower::~HFShower ( )
virtual

Member Function Documentation

bool HFShower::compute ( )

Compute the shower longitudinal and lateral development.

Definition at line 475 of file HFShower.cc.

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), prof2calltree::count, debug, detector, patCandidatesForDimuonsSequences_cff::ecal, eStep, RandomEngine::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, getHLTprescales::index, indexFinder(), EcalCondDB::inf, infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, CosmicsPD_Skims::radius, random, query::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), ntuplemaker::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

475  {
476 
477  // TimeMe theT("FamosHFShower::compute");
478 
479  bool status = false;
480  int numLongit = eStep.size();
481  if(debug)
482  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
483  << " N_long.steps required : " << numLongit << std::endl;
484 
485  if(numLongit > 0) {
486 
487  status = true;
488  // Prepare the trsanverse probability function
489  std::vector<double> Fhist;
490  std::vector<double> rhist;
491  for (int j = 0; j < nTRsteps + 1; j++) {
492  rhist.push_back(maxTRfactor * j / nTRsteps );
493  Fhist.push_back(transProb(maxTRfactor,1.,rhist[j]));
494  if(debug == 3)
495  LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
496  }
497 
498  //================================================================
499  // Longitudinal steps
500  //================================================================
501  for (int i = 0; i < numLongit ; i++) {
502 
503  double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
504  // vary the longitudinal profile if needed
505  if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor;
506  if(debug)
507  LogDebug("FastCalorimetry") << " FamosHFShower::compute - detector = "
508  << detector[i]
509  << " currentDepthL0 = "
510  << currentDepthL0 << std::endl;
511 
512  double maxTRsize = maxTRfactor * rlamStep[i]; // in lambda units
513  double rbinsize = maxTRsize / nTRsteps;
514  double espot = eStep[i] / (double)nspots[i]; // re-adjust espot
515 
516  if( espot > 4. || espot < 0. )
517  LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = "
518  << espot << std::endl;
519 
520  int ecal = 0;
521  if(detector[i] != 1) {
522  bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
523 
524  if(debug)
525  LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of "
526  << " theHcalHitMaker->setDepth(currentDepthL0) is "
527  << setHDdepth << std::endl;
528 
529  if(!setHDdepth)
530  {
531  currentDepthL0 -= lamstep[i];
532  setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
533  }
534  if(!setHDdepth) continue;
535 
537  }
538  else {
539 
540  ecal = 1;
541  bool status = theGrid->getPads(currentDepthL0);
542 
543  if(debug)
544  LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = "
545  << status << std::endl;
546 
547  if(!status) continue;
548 
549  theGrid->setSpotEnergy(espot);
550  }
551 
552 
553  //------------------------------------------------------------
554  // Transverse distribution
555  //------------------------------------------------------------
556  int nok = 0; // counter of OK
557  int count = 0;
558  int inf = infinity;
559  if(lossesOpt) inf = nspots[i]; // losses are enabled, otherwise
560  // only OK points are counted ...
561 
562  // Total energy in this layer
563  double eremaining = eStep[i];
564  bool converged = false;
565 
566  while (eremaining > 0. && !converged && count<inf ) {
567 
568  ++count;
569 
570  // energy spot (HFL)
571  double newespot = espot;
572 
573  // We need to know a priori if this energy spot if for
574  // a long (1) or short (2) fiber
575 
576  unsigned layer = 1;
577  if( currentDepthL0 < 1.3 ) // first 22 cm = 1.3 lambda - only HFL
578  layer = 1;
579  else
580  layer = random->flatShoot() < 0.5 ? 1 : 2;
581 
582  if ( layer == 2 ) newespot = 2. * espot;
583 
584  if ( eremaining - newespot < 0. ) newespot = eremaining;
585 
586  // process transverse distribution
587 
588  double prob = random->flatShoot();
589  int index = indexFinder(prob,Fhist);
590  double radius = rlamStep[i] * rhist[index] +
591  random->flatShoot() * rbinsize; // in-bin
592  double phi = 2.*M_PI*random->flatShoot();
593 
594  if(debug == 2)
595  LogDebug("FastCalorimetry") << std::endl << " FamosHFShower::compute " << " r = " << radius
596  << " phi = " << phi << std::endl;
597 
598  // add hit
599 
600  theHcalHitMaker->setSpotEnergy(newespot);
601  theGrid->setSpotEnergy(newespot);
602 
603  bool result;
604  if(ecal) {
605  result = theGrid->addHit(radius,phi,0); // shouldn't get here !
606 
607  if(debug == 2)
608  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
609  << " theGrid->addHit result = "
610  << result << std::endl;
611 
612  }
613  else {
614  // PV assign espot to long/short fibers
615  result = theHcalHitMaker->addHit(radius,phi,layer);
616 
617  if(debug == 2)
618  LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
619  << " theHcalHitMaker->addHit result = "
620  << result << std::endl;
621 
622  }
623 
624  if (result) {
625  ++nok;
626  eremaining -= newespot;
627  if ( eremaining <= 0. ) converged = true;
628  // std::cout << "Transverse : " << nok << " "
629  // << " , E= " << newespot
630  // << " , Erem = " << eremaining
631  // << std::endl;
632  } else {
633  // std::cout << "WARNING : hit not added" << std::endl;
634  }
635  }
636 
637  // end of tranverse simulation
638  //-----------------------------------------------------
639 
640  if(count == infinity) {
641  status = false;
642  if(debug)
643  LogDebug("FastCalorimetry") << "*** FamosHFShower::compute " << " maximum number of"
644  << " transverse points " << count << " is used !!!" << std::endl;
645  break;
646  }
647 
648  if(debug)
649  LogDebug("FastCalorimetry") << " FamosHFShower::compute " << " long.step No." << i
650  << " Ntry, Nok = " << count
651  << " " << nok << std::endl;
652 
653  } // end of longitudinal steps
654  } // end of no steps
655  return status;
656 
657 }
#define LogDebug(id)
void setSpotEnergy(double e)
Set the spot energy.
Definition: HcalHitMaker.h:28
std::vector< double > eStep
Definition: HFShower.h:80
int i
Definition: DBlmapReader.cc:9
std::vector< double > rlamStep
Definition: HFShower.h:80
std::vector< int > nspots
Definition: HFShower.h:79
bool addHit(double r, double phi, unsigned layer=0)
void setSpotEnergy(double e)
Definition: EcalHitMaker.h:117
int infinity
Definition: HFShower.h:84
std::vector< double > lamstep
Definition: HFShower.h:82
int indexFinder(double x, const std::vector< double > &Fhist)
Definition: HFShower.cc:659
std::vector< int > detector
Definition: HFShower.h:79
HcalHitMaker * theHcalHitMaker
Definition: HFShower.h:90
double transProb(double factor, double R, double r)
Definition: HFShower.h:54
tuple result
Definition: query.py:137
string inf
Definition: EcalCondDB.py:94
double maxTRfactor
Definition: HFShower.h:115
int j
Definition: DBlmapReader.cc:9
double hcalDepthFactor
Definition: HFShower.h:119
EcalHitMaker * theGrid
Definition: HFShower.h:87
#define M_PI
Definition: BFit3D.cc:3
int lossesOpt
Definition: HFShower.h:99
std::vector< double > lamtotal
Definition: HFShower.h:82
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
bool getPads(double depth, bool inCm=false)
int nTRsteps
Definition: HFShower.h:103
const RandomEngine * random
Definition: HFShower.h:122
tuple status
Definition: ntuplemaker.py:245
#define debug
Definition: MEtoEDMFormat.h:34
bool addHit(double r, double phi, unsigned layer=0)
add the hit in the HCAL in local coordinates
Definition: HcalHitMaker.cc:28
bool setDepth(double, bool inCm=false)
set the depth in X0 or Lambda0 units depending on showerType
Definition: DDAxes.h:10
double HFShower::gam ( double  x,
double  a 
) const
inlineprivate

Definition at line 49 of file HFShower.h.

References create_public_lumi_plots::exp, and funct::pow().

Referenced by makeSteps().

49 { return pow(x,a-1.)*exp(-x); }
double a
Definition: hdecay.h:121
x
Definition: VDTMath.h:216
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
std::vector< double > HFShower::getDDDArray ( const std::string &  str,
const DDsvalues_type sv,
int &  nmin 
)
private

Definition at line 301 of file HFShower.cc.

References DDfetch(), DDValue::doubles(), edm::hlt::Exception, LogDebug, and relativeConstraints::value.

Referenced by HFShower().

304 {
305 #ifdef DebugLog
306  LogDebug("HFShower") << "HFShower:getDDDArray called for " << str
307  << " with nMin " << nmin;
308 #endif
309  DDValue value(str);
310  if (DDfetch(&sv, value)) {
311 #ifdef DebugLog
312  LogDebug("HFShower") << value;
313 #endif
314  const std::vector<double> & fvec = value.doubles();
315  int nval = fvec.size();
316  if (nmin > 0) {
317  if (nval < nmin) {
318  edm::LogError("HFShower") << "HFShower : # of " << str << " bins "
319  << nval << " < " << nmin << " ==> illegal";
320  throw cms::Exception("Unknown", "HFShower")
321  << "nval < nmin for array " << str << "\n";
322  }
323  } else {
324  if (nval < 2) {
325  edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval
326  << " < 2 ==> illegal" << " (nmin=" << nmin << ")";
327  throw cms::Exception("Unknown", "HFShower")
328  << "nval < 2 for array " << str << "\n";
329  }
330  }
331  nmin = nval;
332 
333  return fvec;
334  } else {
335  edm::LogError("HFShower") << "HFShower : cannot get array " << str;
336  throw cms::Exception("Unknown", "HFShower") << "cannot get array " << str << "\n";
337  }
338 }
#define LogDebug(id)
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
std::vector< HFShower::Hit > HFShower::getHits ( G4Step *  aStep)

Definition at line 70 of file HFShower.cc.

References abs, applyFidCut, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), HFShower::Hit::depth, relval_parameters_module::energy, create_public_lumi_plots::exp, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, probMax, diffTwoXMLs::r1, diffTwoXMLs::r2, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, w(), and HFShower::Hit::wavelength.

Referenced by HCalSD::hitForFibre(), and FiberSD::ProcessHits().

70  {
71 
72  std::vector<HFShower::Hit> hits;
73  int nHit = 0;
74 
75  double edep = aStep->GetTotalEnergyDeposit();
76  double stepl = 0.;
77 
78  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
79  stepl = aStep->GetStepLength();
80  if ((edep == 0.) || (stepl == 0.)) {
81 #ifdef DebugLog
82  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
83 #endif
84  return hits;
85  }
86  G4Track *aTrack = aStep->GetTrack();
87  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
88 
90  double energy = aParticle->GetTotalEnergy();
91  double momentum = aParticle->GetTotalMomentum();
92  double pBeta = momentum / energy;
93  double dose = 0.;
94  int npeDose = 0;
95 
96  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
97  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
98 
99  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
100  G4ThreeVector globalPos = preStepPoint->GetPosition();
101  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
102  //double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
103  double zv = std::abs(globalPos.z()) - gpar[4];
104  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
105  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
106  GetTopTransform().TransformAxis(momentumDir);
107  // @@ Here the depth should be changed (Fibers all long in Geometry!)
108  int depth = 1;
109  int npmt = 0;
110  bool ok = true;
111  if (zv < 0. || zv > gpar[1]) {
112  ok = false; // beyond the fiber in Z
113  }
114  if (ok && applyFidCut) {
115  npmt = HFFibreFiducial::PMTNumber(globalPos);
116  if (npmt <= 0) {
117  ok = false;
118  } else if (npmt > 24) { // a short fibre
119  if (zv > gpar[0]) {
120  depth = 2;
121  } else {
122  ok = false;
123  }
124  }
125 #ifdef DebugLog
126  edm::LogInfo("HFShower") << "HFShower: npmt " << npmt
127  << " zv " << std::abs(pe_effect.z())
128  << ":" << gpar[4] << ":" << zv << ":"
129  << gpar[0] << " ok " << ok << " depth " << depth;
130 #endif
131  } else {
132  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
133  }
134  G4ThreeVector translation =
135  preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
136 
137  double u = localMom.x();
138  double v = localMom.y();
139  double w = localMom.z();
140  double zCoor = localPos.z();
141  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
142  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
143  double time = fibre->tShift(localPos, depth, chkFibre);
144 
145 #ifdef DebugLog
146  LogDebug("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor
147  << "(" << globalPos.z() << ") " << zFibre << " Trans "
148  << translation << " Time " << tSlice << " " << time
149  << "\n Direction " << momentumDir
150  << " Local " << localMom;
151 #endif
152  // Here npe should be 0 if there is no fiber (@@ M.K.)
153  int npe = 1;
154  std::vector<double> wavelength;
155  std::vector<double> momz;
156  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
157  if (ok) npe = cherenkov->computeNPE(particleDef,pBeta,u,v,w,stepl,zFibre,dose, npeDose);
158  wavelength = cherenkov->getWL();
159  momz = cherenkov->getMom();
160  } // ^^^^^ End of Tmp close of the cherenkov function
161  if(ok && npe>0) {
162  for (int i = 0; i<npe; ++i) {
163  double p=1.;
164  if (!applyFidCut) p = fibre->attLength(wavelength[i]);
165  double r1 = G4UniformRand();
166  double r2 = G4UniformRand();
167 #ifdef DebugLog
168  LogDebug("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1
169  <<":" << exp(-p*zFibre) << " r2 " << r2 << ":"
170  << probMax << " Survive: "
171  << (r1 <= exp(-p*zFibre) && r2 <= probMax);
172 #endif
173  if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
174  hit.depth = depth;
175  hit.time = tSlice+time;
176  if (!applyFidCut) {
177  hit.wavelength = wavelength[i]; // Tmp
178  hit.momentum = momz[i]; // Tmp
179  } else {
180  hit.wavelength = 300.; // Tmp
181  hit.momentum = 1.; // Tmp
182  }
183  hit.position = globalPos;
184  hits.push_back(hit);
185  nHit++;
186  }
187  }
188  }
189 
190 #ifdef DebugLog
191  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
192  for (int i=0; i<nHit; ++i)
193  LogDebug("HFShower") << "HFShower::Hit " << i << " WaveLength "
194  << hits[i].wavelength << " Time " << hits[i].time
195  << " Momentum " << hits[i].momentum << " Position "
196  << hits[i].position;
197 #endif
198  return hits;
199 
200 }
#define LogDebug(id)
double wavelength
Definition: HFShower.h:35
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:52
static int PMTNumber(G4ThreeVector pe_effect)
double momentum
Definition: HFShower.h:36
double attLength(double lambda)
Definition: HFFibre.cc:113
bool applyFidCut
Definition: HFShower.h:47
#define abs(x)
Definition: mlp_lapack.h:159
G4ThreeVector position
Definition: HFShower.h:37
int computeNPE(G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:80
HFCherenkov * cherenkov
Definition: HFShower.h:51
int chkFibre
Definition: HFShower.h:54
double probMax
Definition: HFShower.h:55
std::vector< double > getWL()
Definition: HFCherenkov.cc:291
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
std::vector< double > getMom()
Definition: HFCherenkov.cc:296
std::vector< double > gpar
Definition: HFShower.h:56
double time
Definition: HFShower.h:34
mathSSE::Vec4< T > v
T w() const
std::vector< HFShower::Hit > HFShower::getHits ( G4Step *  aStep,
bool  forLibrary 
)

Definition at line 202 of file HFShower.cc.

References abs, applyFidCut, cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, w(), and HFShower::Hit::wavelength.

202  {
203  std::vector<HFShower::Hit> hits;
204  int nHit = 0;
205 
206  double edep = aStep->GetTotalEnergyDeposit();
207  double stepl = 0.;
208 
209  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
210  stepl = aStep->GetStepLength();
211  if ((edep == 0.) || (stepl == 0.)) {
212 #ifdef DebugLog
213  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
214 #endif
215  return hits;
216  }
217 
218  G4Track *aTrack = aStep->GetTrack();
219  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
220 
222  double energy = aParticle->GetTotalEnergy();
223  double momentum = aParticle->GetTotalMomentum();
224  double pBeta = momentum / energy;
225  double dose = 0.;
226  int npeDose = 0;
227 
228  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
229  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
230 
231  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
232  G4ThreeVector globalPos = preStepPoint->GetPosition();
233  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
234  double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
235  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
236  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
237  GetTopTransform().TransformAxis(momentumDir);
238  // @@ Here the depth should be changed (Fibers are all long in Geometry!)
239  int depth = 1;
240  int npmt = 0;
241  bool ok = true;
242  if (zv < 0 || zv > gpar[1]) {
243  ok = false; // beyond the fiber in Z
244  }
245  if (ok && applyFidCut) {
246  npmt = HFFibreFiducial:: PMTNumber(globalPos);
247  if (npmt <= 0) {
248  ok = false;
249  } else if (npmt > 24) { // a short fibre
250  if (zv > gpar[0]) {
251  depth = 2;
252  } else {
253  ok = false;
254  }
255  }
256 #ifdef DebugLog
257  edm::LogInfo("HFShower") << "HFShower:getHits(SL): npmt " << npmt
258  << " zv " << std::abs(pe_effect.z())
259  << ":" << gpar[4] << ":" << zv << ":"
260  << gpar[0] << " ok " << ok << " depth " << depth;
261 #endif
262  } else {
263  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
264  }
265  G4ThreeVector translation =
266  preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
267 
268  double u = localMom.x();
269  double v = localMom.y();
270  double w = localMom.z();
271  double zCoor = localPos.z();
272  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
273  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
274  double time = fibre->tShift(localPos, depth, chkFibre);
275 
276 #ifdef DebugLog
277  LogDebug("HFShower") << "HFShower::getHits(SL): in " << name << " Z " <<zCoor
278  << "(" << globalPos.z() << ") " << zFibre << " Trans "
279  << translation << " Time " << tSlice << " " << time
280  << "\n Direction " << momentumDir
281  << " Local " << localMom;
282 #endif
283  // npe should be 0
284  int npe = 0;
285  if(ok) npe = cherenkov->computeNPE(particleDef,pBeta,u,v,w, stepl,zFibre, dose, npeDose);
286  std::vector<double> wavelength = cherenkov->getWL();
287  std::vector<double> momz = cherenkov->getMom();
288 
289  for (int i = 0; i<npe; ++i) {
290  hit.time = tSlice+time;
291  hit.wavelength = wavelength[i];
292  hit.momentum = momz[i];
293  hit.position = globalPos;
294  hits.push_back(hit);
295  nHit++;
296  }
297 
298  return hits;
299 }
#define LogDebug(id)
double wavelength
Definition: HFShower.h:35
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:52
static int PMTNumber(G4ThreeVector pe_effect)
double momentum
Definition: HFShower.h:36
bool applyFidCut
Definition: HFShower.h:47
#define abs(x)
Definition: mlp_lapack.h:159
G4ThreeVector position
Definition: HFShower.h:37
int computeNPE(G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:80
HFCherenkov * cherenkov
Definition: HFShower.h:51
int chkFibre
Definition: HFShower.h:54
std::vector< double > getWL()
Definition: HFCherenkov.cc:291
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
std::vector< double > getMom()
Definition: HFCherenkov.cc:296
std::vector< double > gpar
Definition: HFShower.h:56
double time
Definition: HFShower.h:34
mathSSE::Vec4< T > v
T w() const
int HFShower::indexFinder ( double  x,
const std::vector< double > &  Fhist 
)
private

Definition at line 659 of file HFShower.cc.

References debug, LogDebug, findQualityFiles::size, and launcher::step.

Referenced by compute().

659  {
660  // binary search in the vector of doubles
661  int size = Fhist.size();
662 
663  int curr = size / 2;
664  int step = size / 4;
665  int iter;
666  int prevdir = 0;
667  int actudir = 0;
668 
669  for (iter = 0; iter < size ; iter++) {
670 
671  if( curr >= size || curr < 1 )
672  LogWarning("FastCalorimetry") << " FamosHFShower::indexFinder - wrong current index = "
673  << curr << " !!!" << std::endl;
674 
675  if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break;
676  prevdir = actudir;
677  if(x > Fhist[curr]) {actudir = 1;}
678  else {actudir = -1;}
679  if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
680  curr += actudir * step;
681  if(curr > size) curr = size;
682  else { if(curr < 1) {curr = 1;}}
683 
684  if(debug == 3)
685  LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter
686  << " curr, F[curr-1], F[curr] = "
687  << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;
688 
689  }
690 
691  if(debug == 3)
692  LogDebug("FastCalorimetry") << " indexFinder x = " << x << " found index = " << curr-1
693  << std::endl;
694 
695 
696  return curr-1;
697 }
#define LogDebug(id)
list step
Definition: launcher.py:15
x
Definition: VDTMath.h:216
#define debug
Definition: MEtoEDMFormat.h:34
tuple size
Write out results.
void HFShower::initRun ( G4ParticleTable *  theParticleTable)

Definition at line 340 of file HFShower.cc.

Referenced by HCalSD::initRun().

340  {// Define PDG codes
341 
342  /*
343  emPDG = theParticleTable->FindParticle("e-")->GetPDGEncoding();
344  epPDG = theParticleTable->FindParticle("e+")->GetPDGEncoding();
345  gammaPDG = theParticleTable->FindParticle("gamma")->GetPDGEncoding();
346 #ifdef DebugLog
347  edm::LogInfo("HFShower") << "HFShower: Particle codes"
348  << " for e- = " << emPDG
349  << " for e+ = " << epPDG
350  << " for gamma = " << gammaPDG;
351 #endif
352  */
353 }
void HFShower::makeSteps ( int  nsteps)
private

Definition at line 366 of file HFShower.cc.

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

Referenced by HFShower().

366  {
367 
368  double sumes = 0.;
369  double sum = 0.;
370  std::vector<double> temp;
371 
372 
373  if(debug)
374  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
375  << " nsteps required : " << nsteps << std::endl;
376 
377  int count = 0;
378  for (int i = 0; i < nsteps; i++) {
379 
380  double deplam = lamdepth[i] - 0.5 * lamstep[i];
381  double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i];
382  double x = betEM * depx0;
383  double y = betHD * deplam;
384 
385  if(debug == 2)
386  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps "
387  << " - step " << i
388  << " depx0, x = " << depx0 << ", " << x
389  << " deplam, y = " << deplam << ", "
390  << y << std::endl;
391 
392  double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
393  (x0curr[i] * tgamEM) +
394  (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
395 
396  // protection ...
397  if(est < 0.) {
398  LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " - negative step energy !!!"
399  << std::endl;
400  est = 0.;
401  break ;
402  }
403 
404  // for estimates only
405  sum += est;
406  int nPest = (int) (est * e / sum / eSpotSize) ;
407 
408  if(debug == 2)
409  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - nPoints estimate = "
410  << nPest << std::endl;
411 
412  if(nPest <= 1 && count !=0 ) break;
413 
414  // good step - to proceed
415 
416  temp.push_back(est);
417  sumes += est;
418 
419  rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge))
420  * deplam * transFactor);
421  count ++;
422  }
423 
424  // fluctuations in ECAL and re-distribution of remaining energy in HCAL
425  if(detector[0] == 1 && count > 1) {
426  double oldECALenergy = temp[0];
427  double oldHCALenergy = sumes - oldECALenergy ;
428  double newECALenergy = 2. * sumes;
429  for (int i = 0; newECALenergy > sumes && i < infinity; i++)
430  newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy;
431 
432  if(debug == 2)
433  LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " ECAL fraction : old/new - "
434  << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl;
435 
436  temp[0] = newECALenergy;
437  double newHCALenergy = sumes - newECALenergy;
438  double newHCALreweight = newHCALenergy / oldHCALenergy;
439 
440  for (int i = 1; i < count; i++) {
441  temp[i] *= newHCALreweight;
442  }
443 
444  }
445 
446 
447  // final re-normalization of the energy fractions
448  double etot = 0.;
449  for (int i = 0; i < count ; i++) {
450  eStep.push_back(temp[i] * e / sumes );
451  nspots.push_back((int)(eStep[i]/eSpotSize)+1);
452  etot += eStep[i];
453 
454  if(debug)
455  LogDebug("FastCalorimetry") << i << " xO and lamdepth at the end of step = "
456  << x0depth[i] << " "
457  << lamdepth[i] << " Estep func = " << eStep[i]
458  << " Rstep = " << rlamStep[i] << " Nspots = " << nspots[i]
459  << std::endl;
460 
461  }
462 
463  // The only step is in ECAL - let's make the size bigger ...
464  if(count == 1 and detector[0] == 1) rlamStep[0] *= 2.;
465 
466 
467  if(debug) {
468  if(eStep[0] > 0.95 * e && detector[0] == 1)
469  LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - " << "ECAL energy = " << eStep[0]
470  << " out of total = " << e << std::endl;
471  }
472 
473 }
#define LogDebug(id)
std::vector< double > lamdepth
Definition: HFShower.h:82
std::vector< double > eStep
Definition: HFShower.h:80
double balanceEH
Definition: HFShower.h:117
int i
Definition: DBlmapReader.cc:9
double alpHD
Definition: HFShower.h:72
std::vector< double > rlamStep
Definition: HFShower.h:80
std::vector< int > nspots
Definition: HFShower.h:79
double theR2
Definition: HFShower.h:71
double transParam
Definition: HFShower.h:105
int infinity
Definition: HFShower.h:84
double betHD
Definition: HFShower.h:72
std::vector< double > lamstep
Definition: HFShower.h:82
double theR1
Definition: HFShower.h:71
double eSpotSize
Definition: HFShower.h:109
double transFactor
Definition: HFShower.h:107
std::vector< int > detector
Definition: HFShower.h:79
double tgamEM
Definition: HFShower.h:72
double aloge
Definition: HFShower.h:77
double theR3
Definition: HFShower.h:71
double gam(double x, double a) const
Definition: HFShower.h:49
std::vector< double > x0curr
Definition: HFShower.h:81
double e
Definition: HFShower.h:96
std::vector< double > lamcurr
Definition: HFShower.h:82
part
Definition: HCALResponse.h:21
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
const RandomEngine * random
Definition: HFShower.h:122
x
Definition: VDTMath.h:216
#define debug
Definition: MEtoEDMFormat.h:34
double tgamHD
Definition: HFShower.h:72
double betEM
Definition: HFShower.h:72
double alpEM
Definition: HFShower.h:72
std::vector< double > x0depth
Definition: HFShower.h:81
double HFShower::transProb ( double  factor,
double  R,
double  r 
)
inlineprivate

Definition at line 54 of file HFShower.h.

References dttmaxenums::R.

Referenced by compute().

54  {
55  double fsq = factor * factor;
56  return ((fsq + 1.)/fsq) * r * r / (r*r + R*R) ;
57  }

Member Data Documentation

double HFShower::aloge
private

Definition at line 77 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpEM
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpHD
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

bool HFShower::applyFidCut
private

Definition at line 47 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::balanceEH
private

Definition at line 117 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betEM
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betHD
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

HFCherenkov* HFShower::cherenkov
private

Definition at line 51 of file HFShower.h.

Referenced by getHits(), and HFShower().

int HFShower::chkFibre
private

Definition at line 54 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::criticalEnergy
private

Definition at line 113 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStart
private

Definition at line 76 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStep
private

Definition at line 111 of file HFShower.h.

Referenced by HFShower().

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

Definition at line 79 of file HFShower.h.

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

double HFShower::e
private

Definition at line 96 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::eSpotSize
private

Definition at line 109 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

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

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

HFFibre* HFShower::fibre
private

Definition at line 52 of file HFShower.h.

Referenced by getHits(), and HFShower().

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

Definition at line 56 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::hcalDepthFactor
private

Definition at line 119 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::infinity
private

Definition at line 84 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::lambdaEM
private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::lambdaHD
private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

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

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

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

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

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

Definition at line 82 of file HFShower.h.

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

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

Definition at line 82 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::lossesOpt
private

Definition at line 99 of file HFShower.h.

Referenced by compute(), and HFShower().

double HFShower::maxTRfactor
private

Definition at line 115 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::nDepthSteps
private

Definition at line 101 of file HFShower.h.

Referenced by HFShower().

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

Definition at line 79 of file HFShower.h.

Referenced by compute(), and makeSteps().

int HFShower::nTRsteps
private

Definition at line 103 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::onEcal
private

Definition at line 93 of file HFShower.h.

Referenced by HFShower().

double HFShower::part
private

Definition at line 72 of file HFShower.h.

double HFShower::probMax
private

Definition at line 55 of file HFShower.h.

Referenced by getHits(), and HFShower().

const RandomEngine* HFShower::random
private

Definition at line 122 of file HFShower.h.

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

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

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::tgamEM
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::tgamHD
private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

const ECALProperties* HFShower::theECALproperties
private

Definition at line 67 of file HFShower.h.

Referenced by HFShower().

EcalHitMaker* HFShower::theGrid
private

Definition at line 87 of file HFShower.h.

Referenced by compute(), and HFShower().

HcalHitMaker* HFShower::theHcalHitMaker
private

Definition at line 90 of file HFShower.h.

Referenced by compute().

const HCALProperties* HFShower::theHCALproperties
private

Definition at line 68 of file HFShower.h.

Referenced by HFShower().

HDShowerParametrization* HFShower::theParam
private

Definition at line 64 of file HFShower.h.

Referenced by HFShower().

double HFShower::theR1
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR2
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR3
private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transFactor
private

Definition at line 107 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transParam
private

Definition at line 105 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

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

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

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

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::x0EM
private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::x0HD
private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().