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, double weight)
 
std::vector< HitgetHits (G4Step *aStep, bool forLibrary)
 
std::vector< HitgetHits (G4Step *aStep, bool forLibraryProducer, double zoffset)
 
 HFShower (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk=0)
 
 HFShower (const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
 
void initRun (G4ParticleTable *, HcalDDDSimConstants *)
 
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 RandomEngineAndDistributionrandom
 
std::vector< double > rlamStep
 
double tgamEM
 
double tgamHD
 
const ECALPropertiestheECALproperties
 
EcalHitMakertheGrid
 
HcalHitMakertheHcalHitMaker
 
const HCALPropertiestheHCALproperties
 
HDShowerParametrizationtheParam
 
double theR1
 
double theR2
 
double theR3
 
double transFactor
 
double transParam
 
std::vector< double > x0curr
 
std::vector< double > x0depth
 
double x0EM
 
double x0HD
 

Detailed Description

Definition at line 22 of file HFShower.h.

Member Typedef Documentation

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 RandomEngineAndDistribution 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, RandomEngineAndDistribution::flatShoot(), HSParameters::getHDbalanceEH(), HSParameters::getHDcriticalEnergy(), HSParameters::getHDdepthStep(), HSParameters::getHDeSpotSize(), HSParameters::getHDhcalDepthFactor(), HSParameters::getHDlossesOpt(), HSParameters::getHDmaxTRfactor(), HSParameters::getHDnDepthSteps(), HSParameters::getHDnTRsteps(), HSParameters::getHDtransParam(), hcalDepthFactor, HDShowerParametrization::hcalProperties(), EcalHitMaker::hcalTotalL0(), HDShowerParametrization::hsParameters(), i, ECALProperties::interactionLength(), HCALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, dqm-mbProfile::log, LogDebug, lossesOpt, makeSteps(), HLT_FULL_cff::maxDepth, maxTRfactor, nDepthSteps, nTRsteps, onEcal, HDShowerParametrization::part1(), HDShowerParametrization::part2(), HDShowerParametrization::r1(), HDShowerParametrization::r2(), HDShowerParametrization::r3(), ECALProperties::radLenIncm(), HCALProperties::radLenIncm(), random, HDShowerParametrization::setCase(), combinedConstraintHelpers::sum2(), tgamEM, tgamHD, theECALproperties, theGrid, theHCALproperties, theParam, theR1, theR2, theR3, transFactor, transParam, x0curr, x0depth, x0EM, and x0HD.

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
const RandomEngineAndDistribution * random
Definition: HFShower.h:122
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 flatShoot(double xmin=0.0, double xmax=1.0) const
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
#define debug
Definition: HDRShower.cc:19
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:20
void makeSteps(int nsteps)
Definition: HFShower.cc:366
std::vector< double > lamtotal
Definition: HFShower.h:82
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
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 26 of file HFShower.cc.

References applyFidCut, cherenkov, chkFibre, fibre, edm::ParameterSet::getParameter(), and probMax.

27  : cherenkov(0),
28  fibre(0),
29  chkFibre(chk) {
30 
31  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
32  applyFidCut = m_HF.getParameter<bool>("ApplyFiducialCut");
33  probMax = m_HF.getParameter<double>("ProbMax");
34 
35  edm::LogInfo("HFShower") << "HFShower:: Maximum probability cut off "
36  << probMax << " Check flag " << chkFibre;
37 
38  cherenkov = new HFCherenkov(m_HF);
39  fibre = new HFFibre(name, cpv, p);
40 }
T getParameter(std::string const &) const
HFFibre * fibre
Definition: HFShower.h:55
bool applyFidCut
Definition: HFShower.h:50
HFCherenkov * cherenkov
Definition: HFShower.h:54
int chkFibre
Definition: HFShower.h:57
double probMax
Definition: HFShower.h:58
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(), KineDebug3::count(), debug, detector, patCandidatesForDimuonsSequences_cff::ecal, eStep, RandomEngineAndDistribution::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, cmsHarvester::index, indexFinder(), EcalCondDB::inf, infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, CosmicsPD_Skims::radius, random, mps_fire::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), mps_update::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

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
const RandomEngineAndDistribution * random
Definition: HFShower.h:122
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
double flatShoot(double xmin=0.0, double xmax=1.0) const
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
tuple result
Definition: mps_fire.py:84
HcalHitMaker * theHcalHitMaker
Definition: HFShower.h:90
double transProb(double factor, double R, double r)
Definition: HFShower.h:54
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
#define debug
Definition: HDRShower.cc:19
int lossesOpt
Definition: HFShower.h:99
std::vector< double > lamtotal
Definition: HFShower.h:82
bool getPads(double depth, bool inCm=false)
int nTRsteps
Definition: HFShower.h:103
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
tuple status
Definition: mps_update.py:57
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
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 416 of file HFShower.cc.

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

418  {
419 #ifdef DebugLog
420  LogDebug("HFShower") << "HFShower:getDDDArray called for " << str
421  << " with nMin " << nmin;
422 #endif
423  DDValue value(str);
424  if (DDfetch(&sv, value)) {
425 #ifdef DebugLog
426  LogDebug("HFShower") << value;
427 #endif
428  const std::vector<double> & fvec = value.doubles();
429  int nval = fvec.size();
430  if (nmin > 0) {
431  if (nval < nmin) {
432  edm::LogError("HFShower") << "HFShower : # of " << str << " bins "
433  << nval << " < " << nmin << " ==> illegal";
434  throw cms::Exception("Unknown", "HFShower")
435  << "nval < nmin for array " << str << "\n";
436  }
437  } else {
438  if (nval < 2) {
439  edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval
440  << " < 2 ==> illegal" << " (nmin=" << nmin << ")";
441  throw cms::Exception("Unknown", "HFShower")
442  << "nval < 2 for array " << str << "\n";
443  }
444  }
445  nmin = nval;
446 
447  return fvec;
448  } else {
449  edm::LogError("HFShower") << "HFShower : cannot get array " << str;
450  throw cms::Exception("Unknown", "HFShower") << "cannot get array " << str << "\n";
451  }
452 }
#define LogDebug(id)
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:80
std::vector< HFShower::Hit > HFShower::getHits ( G4Step *  aStep,
double  weight 
)

Definition at line 47 of file HFShower.cc.

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

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

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

Definition at line 317 of file HFShower.cc.

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

317  {
318  std::vector<HFShower::Hit> hits;
319  int nHit = 0;
320 
321  double edep = aStep->GetTotalEnergyDeposit();
322  double stepl = 0.;
323 
324  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
325  stepl = aStep->GetStepLength();
326  if ((edep == 0.) || (stepl == 0.)) {
327 #ifdef DebugLog
328  edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
329 #endif
330  return hits;
331  }
332 
333  G4Track *aTrack = aStep->GetTrack();
334  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
335 
337  double energy = aParticle->GetTotalEnergy();
338  double momentum = aParticle->GetTotalMomentum();
339  double pBeta = momentum / energy;
340  double dose = 0.;
341  int npeDose = 0;
342 
343  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
344  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
345 
346  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
347  G4ThreeVector globalPos = preStepPoint->GetPosition();
348  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
349  double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
350  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
351  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
352  GetTopTransform().TransformAxis(momentumDir);
353  // @@ Here the depth should be changed (Fibers are all long in Geometry!)
354  int depth = 1;
355  int npmt = 0;
356  bool ok = true;
357  if (zv < 0 || zv > gpar[1]) {
358  ok = false; // beyond the fiber in Z
359  }
360  if (ok && applyFidCut) {
361  npmt = HFFibreFiducial:: PMTNumber(globalPos);
362  if (npmt <= 0) {
363  ok = false;
364  } else if (npmt > 24) { // a short fibre
365  if (zv > gpar[0]) {
366  depth = 2;
367  } else {
368  ok = false;
369  }
370  }
371 #ifdef DebugLog
372  edm::LogInfo("HFShower") << "HFShower:getHits(SL): npmt " << npmt
373  << " zv " << std::abs(globalPos.z())
374  << ":" << gpar[4] << ":" << zv << ":"
375  << gpar[0] << " ok " << ok << " depth " << depth;
376 #endif
377  } else {
378  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
379  }
380  G4ThreeVector translation =
381  preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
382 
383  double u = localMom.x();
384  double v = localMom.y();
385  double w = localMom.z();
386  double zCoor = localPos.z();
387  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
388  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
389  double time = fibre->tShift(localPos, depth, chkFibre);
390 
391 #ifdef DebugLog
392  edm::LogInfo("HFShower") << "HFShower::getHits(SL): in " << name << " Z "
393  << zCoor << "(" << globalPos.z() << ") " << zFibre
394  << " Trans " << translation << " Time " << tSlice
395  << " " << time << "\n Direction "
396  << momentumDir << " Local " << localMom;
397 #endif
398  // npe should be 0
399  int npe = 0;
400  if(ok) npe = cherenkov->computeNPE(aStep,particleDef,pBeta,u,v,w, stepl,zFibre, dose, npeDose);
401  std::vector<double> wavelength = cherenkov->getWL();
402  std::vector<double> momz = cherenkov->getMom();
403 
404  for (int i = 0; i<npe; ++i) {
405  hit.time = tSlice+time;
406  hit.wavelength = wavelength[i];
407  hit.momentum = momz[i];
408  hit.position = globalPos;
409  hits.push_back(hit);
410  nHit++;
411  }
412 
413  return hits;
414 }
double tShift(const G4ThreeVector &point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:115
double wavelength
Definition: HFShower.h:36
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:55
const double w
Definition: UKUtility.cc:23
double momentum
Definition: HFShower.h:37
bool applyFidCut
Definition: HFShower.h:50
G4ThreeVector position
Definition: HFShower.h:38
static int PMTNumber(const G4ThreeVector &pe_effect)
HFCherenkov * cherenkov
Definition: HFShower.h:54
int chkFibre
Definition: HFShower.h:57
std::vector< double > getWL()
Definition: HFCherenkov.cc:336
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > getMom()
Definition: HFCherenkov.cc:341
std::vector< double > gpar
Definition: HFShower.h:59
double time
Definition: HFShower.h:35
int computeNPE(G4Step *step, G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:82
std::vector< HFShower::Hit > HFShower::getHits ( G4Step *  aStep,
bool  forLibraryProducer,
double  zoffset 
)

Definition at line 182 of file HFShower.cc.

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

184  {
185 
186  std::vector<HFShower::Hit> hits;
187  int nHit = 0;
188 
189  double edep = aStep->GetTotalEnergyDeposit();
190  double stepl = 0.;
191 
192  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
193  stepl = aStep->GetStepLength();
194  if ((edep == 0.) || (stepl == 0.)) {
195 #ifdef DebugLog
196  edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
197 #endif
198  return hits;
199  }
200  G4Track *aTrack = aStep->GetTrack();
201  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
202 
204  double energy = aParticle->GetTotalEnergy();
205  double momentum = aParticle->GetTotalMomentum();
206  double pBeta = momentum / energy;
207  double dose = 0.;
208  int npeDose = 0;
209 
210  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
211  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
212 
213  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
214  G4ThreeVector globalPos = preStepPoint->GetPosition();
215  G4String name = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
216  //double zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
217  //double zv = std::abs(globalPos.z()) - gpar[4];
218  double zv = gpar[1]-(std::abs(globalPos.z()) - zoffset);
219  G4ThreeVector localPos = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
220  G4ThreeVector localMom = preStepPoint->GetTouchable()->GetHistory()->
221  GetTopTransform().TransformAxis(momentumDir);
222  // @@ Here the depth should be changed (Fibers all long in Geometry!)
223  int depth = 1;
224  int npmt = 0;
225  bool ok = true;
226  if (zv < 0. || zv > gpar[1]) {
227  ok = false; // beyond the fiber in Z
228  }
229  if (ok && applyFidCut) {
230  npmt = HFFibreFiducial::PMTNumber(globalPos);
231  if (npmt <= 0) {
232  ok = false;
233  } else if (npmt > 24) { // a short fibre
234  if (zv > gpar[0]) {
235  depth = 2;
236  } else {
237  ok = false;
238  }
239  }
240 #ifdef DebugLog
241  edm::LogInfo("HFShower") << "HFShower: npmt " << npmt
242  << " zv " << std::abs(globalPos.z())
243  << ":" << gpar[4] << ":" << zv << ":"
244  << gpar[0] << " ok " << ok << " depth " << depth;
245 #endif
246  } else {
247  depth = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
248  }
249  G4ThreeVector translation =
250  preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();
251 
252  double u = localMom.x();
253  double v = localMom.y();
254  double w = localMom.z();
255  double zCoor = localPos.z();
256  double zFibre = (0.5*gpar[1]-zCoor-translation.z());
257  double tSlice = (aStep->GetPostStepPoint()->GetGlobalTime());
258  double time = fibre->tShift(localPos, depth, chkFibre);
259 
260 #ifdef DebugLog
261  edm::LogInfo("HFShower") << "HFShower::getHits: in " << name << " Z " <<zCoor
262  << "(" << globalPos.z() <<") " << zFibre <<" Trans "
263  << translation << " Time " << tSlice << " " << time
264  << "\n Direction " << momentumDir
265  << " Local " << localMom;
266 #endif
267  // Here npe should be 0 if there is no fiber (@@ M.K.)
268  int npe = 1;
269  std::vector<double> wavelength;
270  std::vector<double> momz;
271  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
272  if (ok) npe = cherenkov->computeNPE(aStep,particleDef,pBeta,u,v,w,stepl,zFibre,dose, npeDose);
273  wavelength = cherenkov->getWL();
274  momz = cherenkov->getMom();
275  } // ^^^^^ End of Tmp close of the cherenkov function
276  if (ok && npe>0) {
277  for (int i = 0; i<npe; ++i) {
278  double p=1.;
279  if (!applyFidCut) p = fibre->attLength(wavelength[i]);
280  double r1 = G4UniformRand();
281  double r2 = G4UniformRand();
282 #ifdef DebugLog
283  edm::LogInfo("HFShower") << "HFShower::getHits: " << i << " attenuation "
284  << r1 << ":" << exp(-p*zFibre) << " r2 " << r2
285  << ":" << probMax << " Survive: "
286  << (r1 <= exp(-p*zFibre) && r2 <= probMax);
287 #endif
288  if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
289  hit.depth = depth;
290  hit.time = tSlice+time;
291  if (!applyFidCut) {
292  hit.wavelength = wavelength[i]; // Tmp
293  hit.momentum = momz[i]; // Tmp
294  } else {
295  hit.wavelength = 300.; // Tmp
296  hit.momentum = 1.; // Tmp
297  }
298  hit.position = globalPos;
299  hits.push_back(hit);
300  nHit++;
301  }
302  }
303  }
304 
305 #ifdef DebugLog
306  edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
307  for (int i=0; i<nHit; ++i)
308  edm::LogInfo("HFShower") << "HFShower::Hit " << i << " WaveLength "
309  << hits[i].wavelength << " Time " << hits[i].time
310  << " Momentum " << hits[i].momentum <<" Position "
311  << hits[i].position;
312 #endif
313  return hits;
314 
315 }
double tShift(const G4ThreeVector &point, int depth, int fromEndAbs=0)
Definition: HFFibre.cc:115
double wavelength
Definition: HFShower.h:36
int i
Definition: DBlmapReader.cc:9
HFFibre * fibre
Definition: HFShower.h:55
const double w
Definition: UKUtility.cc:23
double momentum
Definition: HFShower.h:37
double attLength(double lambda)
Definition: HFFibre.cc:97
bool applyFidCut
Definition: HFShower.h:50
G4ThreeVector position
Definition: HFShower.h:38
static int PMTNumber(const G4ThreeVector &pe_effect)
HFCherenkov * cherenkov
Definition: HFShower.h:54
int chkFibre
Definition: HFShower.h:57
double probMax
Definition: HFShower.h:58
std::vector< double > getWL()
Definition: HFCherenkov.cc:336
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > getMom()
Definition: HFCherenkov.cc:341
std::vector< double > gpar
Definition: HFShower.h:59
double time
Definition: HFShower.h:35
int computeNPE(G4Step *step, G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:82
int HFShower::indexFinder ( double  x,
const std::vector< double > &  Fhist 
)
private

Definition at line 659 of file HFShower.cc.

References debug, LogDebug, and findQualityFiles::size.

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: HDRShower.cc:19
step
tuple size
Write out results.
void HFShower::initRun ( G4ParticleTable *  ,
HcalDDDSimConstants hcons 
)

Definition at line 454 of file HFShower.cc.

References fibre, HcalDDDSimConstants::getGparHF(), gpar, and HFFibre::initRun().

Referenced by HCalSD::initRun(), and FiberSD::update().

454  {
455 
456  //Special Geometry parameters
457  gpar = hcons->getGparHF();
458  edm::LogInfo("HFShower") << "HFShower: " << gpar.size() << " gpar (cm)";
459  for (unsigned int ig=0; ig<gpar.size(); ig++)
460  edm::LogInfo("HFShower") << "HFShower: gpar[" << ig << "] = "
461  << gpar[ig]/cm << " cm";
462 
463  if (fibre) fibre->initRun(hcons);
464 }
HFFibre * fibre
Definition: HFShower.h:55
void initRun(HcalDDDSimConstants *)
Definition: HFFibre.cc:80
std::vector< double > gpar
Definition: HFShower.h:59
const std::vector< double > & getGparHF() const
void HFShower::makeSteps ( int  nsteps)
private

Definition at line 366 of file HFShower.cc.

References aloge, alpEM, alpHD, balanceEH, betEM, betHD, KineDebug3::count(), debug, detector, e, eSpotSize, eStep, RandomEngineAndDistribution::flatShoot(), gam(), i, infinity, lamcurr, lamdepth, lamstep, LogDebug, nspots, random, rlamStep, groupFilesInBlocks::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, x, x0curr, x0depth, and 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
const RandomEngineAndDistribution * random
Definition: HFShower.h:122
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
double flatShoot(double xmin=0.0, double xmax=1.0) const
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
#define debug
Definition: HDRShower.cc:19
part
Definition: HCALResponse.h:20
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 V0MonitoringClient_cfi::factor, and 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 50 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 54 of file HFShower.h.

Referenced by getHits(), and HFShower().

int HFShower::chkFibre
private

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

Referenced by getHits(), HFShower(), and initRun().

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

Definition at line 59 of file HFShower.h.

Referenced by getHits(), and initRun().

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

Referenced by getHits(), and HFShower().

const RandomEngineAndDistribution* 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().