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)
 
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
 
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, funct::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
Log< T >::type log(const T &t)
Definition: Log.h:22
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.
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 22 of file HFShower.cc.

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

23  : cherenkov(0),
24  fibre(0),
25  chkFibre(chk) {
26 
27  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
28  probMax = m_HF.getParameter<double>("ProbMax");
29 
30  edm::LogInfo("HFShower") << "HFShower:: Maximum probability cut off "
31  << probMax << " Check flag " << chkFibre;
32 
33 
34  G4String attribute = "ReadOutName";
35  G4String value = name;
37  DDValue ddv(attribute,value,0);
39  DDFilteredView fv(cpv);
40  fv.addFilter(filter);
41  bool dodet = fv.firstChild();
42  if (dodet) {
43  DDsvalues_type sv(fv.mergedSpecifics());
44 
45  //Special Geometry parameters
46  int ngpar = 7;
47  gpar = getDDDArray("gparHF",sv,ngpar);
48  edm::LogInfo("HFShower") << "HFShower: " << ngpar << " gpar (cm)";
49  for (int ig=0; ig<ngpar; ig++)
50  edm::LogInfo("HFShower") << "HFShower: gpar[" << ig << "] = "
51  << gpar[ig]/cm << " cm";
52  } else {
53  edm::LogError("HFShower") << "HFShower: cannot get filtered "
54  << " view for " << attribute << " matching "
55  << name;
56  throw cms::Exception("Unknown", "HFShower")
57  << "cannot match " << attribute << " to " << name <<"\n";
58  }
59 
60  cherenkov = new HFCherenkov(m_HF);
61  fibre = new HFFibre(name, cpv, p);
62 }
T getParameter(std::string const &) const
HFFibre * fibre
Definition: HFShower.h:50
HFCherenkov * cherenkov
Definition: HFShower.h:49
int chkFibre
Definition: HFShower.h:52
double probMax
Definition: HFShower.h:53
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:54
tuple filter
USE THIS FOR SKIMMED TRACKS process.p = cms.Path(process.hltLevel1GTSeed*process.skimming*process.offlineBeamSpot*process.TrackRefitter2) OTHERWISE USE THIS.
Definition: align_tpl.py:86
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
Definition: HFShower.cc:238
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<infinity ) {
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 funct::exp(), and funct::pow().

Referenced by makeSteps().

49 { return pow(x,a-1.)*exp(-x); }
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
double a
Definition: hdecay.h:121
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 238 of file HFShower.cc.

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

Referenced by HFShower().

240  {
241 
242 #ifdef DebugLog
243  LogDebug("HFShower") << "HFShower:getDDDArray called for " << str
244  << " with nMin " << nmin;
245 #endif
246  DDValue value(str);
247  if (DDfetch(&sv,value)) {
248 #ifdef DebugLog
249  LogDebug("HFShower") << value;
250 #endif
251  const std::vector<double> & fvec = value.doubles();
252  int nval = fvec.size();
253  if (nmin > 0) {
254  if (nval < nmin) {
255  edm::LogError("HFShower") << "HFShower : # of " << str
256  << " bins " << nval << " < " << nmin
257  << " ==> illegal";
258  throw cms::Exception("Unknown", "HFShower")
259  << "nval < nmin for array " << str << "\n";
260  }
261  } else {
262  if (nval < 2) {
263  edm::LogError("HFShower") << "HFShower : # of " << str
264  << " bins " << nval << " < 2 ==> illegal"
265  << " (nmin=" << nmin << ")";
266  throw cms::Exception("Unknown", "HFShower")
267  << "nval < 2 for array " << str << "\n";
268  }
269  }
270  nmin = nval;
271 
272  return fvec;
273  } else {
274  edm::LogError("HFShower") << "HFShower : cannot get array " << str;
275  throw cms::Exception("Unknown", "HFShower")
276  << "cannot get array " << str << "\n";
277  }
278 }
#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 69 of file HFShower.cc.

References abs, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, funct::exp(), fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, AlCaRecoCosmics_cfg::name, L1TEmulatorMonitor_cff::p, HFShower::Hit::position, probMax, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, and HFShower::Hit::wavelength.

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

69  {
70 
71  std::vector<HFShower::Hit> hits;
72  int nHit = 0;
73 
74  double edep = aStep->GetTotalEnergyDeposit();
75  double stepl = 0.;
76 
77  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
78  stepl = aStep->GetStepLength();
79  if ((edep == 0.) || (stepl == 0.)) {
80 #ifdef DebugLog
81  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
82 #endif
83  return hits;
84  }
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 =
102  preStepPoint->GetTouchable()->GetSolid(0)->GetName();
103  double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
104  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(),
105  zv);
106  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
107  GetTopTransform().TransformAxis(momentumDir);
108  int depth =
109  (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10;
110  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
111 
112  double u = localMom.x();
113  double v = localMom.y();
114  double w = localMom.z();
115  double zCoor = localPos.z();
116  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
117  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
118  double time = fibre->tShift(localPos, depth, chkFibre);
119 
120 #ifdef DebugLog
121  LogDebug("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor
122  << "(" << globalPos.z() << ") " << zFibre << " Trans "
123  << translation << " Time " << tSlice << " " << time
124  << "\n Direction " << momentumDir
125  << " Local " << localMom;
126 #endif
127  int npe = cherenkov->computeNPE(particleDef, pBeta, u, v, w, stepl, zFibre,
128  dose, npeDose);
129  std::vector<double> wavelength = cherenkov->getWL();
130  std::vector<double> momz = cherenkov->getMom();
131 
132  for (int i = 0; i<npe; ++i) {
133  double p = fibre->attLength(wavelength[i]);
134  double r1 = G4UniformRand();
135  double r2 = G4UniformRand();
136 #ifdef DebugLog
137  LogDebug("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1
138  <<":" << exp(-p*zFibre) << " r2 " << r2 << ":"
139  << probMax << " Survive: "
140  << (r1 <= exp(-p*zFibre) && r2 <= probMax);
141 #endif
142  if (chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
143  hit.time = tSlice+time;
144  hit.wavelength = wavelength[i];
145  hit.momentum = momz[i];
146  hit.position = globalPos;
147  hits.push_back(hit);
148  nHit++;
149  }
150  }
151 
152 #ifdef DebugLog
153  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
154  for (int i=0; i<nHit; i++)
155  LogDebug("HFShower") << "HFShower::Hit " << i << " WaveLength "
156  << hits[i].wavelength << " Time " << hits[i].time
157  << " Momentum " << hits[i].momentum << " Position "
158  << hits[i].position;
159 #endif
160  return hits;
161 
162 }
#define LogDebug(id)
double wavelength
Definition: HFShower.h:34
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:50
double momentum
Definition: HFShower.h:35
double attLength(double lambda)
Definition: HFFibre.cc:113
#define abs(x)
Definition: mlp_lapack.h:159
G4ThreeVector position
Definition: HFShower.h:36
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
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:85
HFCherenkov * cherenkov
Definition: HFShower.h:49
int chkFibre
Definition: HFShower.h:52
double probMax
Definition: HFShower.h:53
std::vector< double > getWL()
Definition: HFCherenkov.cc:306
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
std::vector< double > getMom()
Definition: HFCherenkov.cc:311
std::vector< double > gpar
Definition: HFShower.h:54
double time
Definition: HFShower.h:33
mathSSE::Vec4< T > v
std::vector< HFShower::Hit > HFShower::getHits ( G4Step *  aStep,
bool  forLibrary 
)

Definition at line 164 of file HFShower.cc.

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

164  {
165 
166  std::vector<HFShower::Hit> hits;
167  int nHit = 0;
168 
169  double edep = aStep->GetTotalEnergyDeposit();
170  double stepl = 0.;
171 
172  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
173  stepl = aStep->GetStepLength();
174  if ((edep == 0.) || (stepl == 0.)) {
175 #ifdef DebugLog
176  LogDebug("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
177 #endif
178  return hits;
179  }
180 
181  G4Track *aTrack = aStep->GetTrack();
182  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
183 
185  double energy = aParticle->GetTotalEnergy();
186  double momentum = aParticle->GetTotalMomentum();
187  double pBeta = momentum / energy;
188  double dose = 0.;
189  int npeDose = 0;
190 
191  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
192  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
193 
194  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
195  G4ThreeVector globalPos = preStepPoint->GetPosition();
196  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
197  double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
198  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(),
199  zv);
200  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
201  GetTopTransform().TransformAxis(momentumDir);
202  int depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10;
203  G4ThreeVector translation = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
204 
205  double u = localMom.x();
206  double v = localMom.y();
207  double w = localMom.z();
208  double zCoor = localPos.z();
209  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
210  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
211  double time = fibre->tShift(localPos, depth, chkFibre);
212 
213 #ifdef DebugLog
214  LogDebug("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor
215  << "(" << globalPos.z() << ") " << zFibre << " Trans "
216  << translation << " Time " << tSlice << " " << time
217  << "\n Direction " << momentumDir
218  << " Local " << localMom;
219 #endif
220  int npe = cherenkov->computeNPE(particleDef, pBeta, u, v, w, stepl, zFibre,
221  dose, npeDose);
222  std::vector<double> wavelength = cherenkov->getWL();
223  std::vector<double> momz = cherenkov->getMom();
224 
225  for (int i = 0; i<npe; ++i) {
226  hit.time = tSlice+time;
227  hit.wavelength = wavelength[i];
228  hit.momentum = momz[i];
229  hit.position = globalPos;
230  hits.push_back(hit);
231  nHit++;
232  }
233 
234  return hits;
235 
236 }
#define LogDebug(id)
double wavelength
Definition: HFShower.h:34
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:50
double momentum
Definition: HFShower.h:35
#define abs(x)
Definition: mlp_lapack.h:159
G4ThreeVector position
Definition: HFShower.h:36
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:85
HFCherenkov * cherenkov
Definition: HFShower.h:49
int chkFibre
Definition: HFShower.h:52
std::vector< double > getWL()
Definition: HFCherenkov.cc:306
double tShift(G4ThreeVector point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:131
std::vector< double > getMom()
Definition: HFCherenkov.cc:311
std::vector< double > gpar
Definition: HFShower.h:54
double time
Definition: HFShower.h:33
mathSSE::Vec4< T > v
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 ExpressReco_HICollisions_FallBack::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)
#define debug
Definition: MEtoEDMFormat.h:34
tuple size
Write out results.
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, cond::rpcobtemp::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, ExpressReco_HICollisions_FallBack::x, x0curr, x0depth, and ExpressReco_HICollisions_FallBack::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
#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().

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 49 of file HFShower.h.

Referenced by getHits(), and HFShower().

int HFShower::chkFibre
private

Definition at line 52 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 50 of file HFShower.h.

Referenced by getHits(), and HFShower().

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

Definition at line 54 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 53 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().