38 theHcalHitMaker(myHcalHitMaker),
64 if (
e < 50.) depthExt = 0.8 * (50. -
e) / 50. + 0.3;
66 if (
e < 500.) depthExt = (500. -
e) / 500. * 0.4 - 0.1;
95 if(
e < emin) effective = emin;
102 if(effective > 0.5 * emax) {
104 if(effective > emax) {
112 LogDebug(
"FastCalorimetry") <<
" HFShower : " << std::endl
113 <<
" Energy " <<
e << std::endl
114 <<
" lossesOpt " <<
lossesOpt << std::endl
116 <<
" nTRsteps " <<
nTRsteps << std::endl
118 <<
" eSpotSize " <<
eSpotSize << std::endl
121 <<
" balanceEH " <<
balanceEH << std::endl
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;
163 alpEM = alpEM1 + alpEM2 * aedep;
165 betEM = betEM1 - betEM2 * aedep;
166 alpHD = alpHD1 + alpHD2 * aedep;
168 betHD = betHD1 - betHD2 * aedep;
169 part = part1 - part2 * aedep;
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;
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;
217 double depthECAL = 0.;
218 double depthGAP = 0.;
219 double depthGAPx0 = 0.;
227 double depthToHCAL = depthECAL + depthGAP;
239 LogDebug(
"FastCalorimetry") <<
" FamosHFShower : e <emin -> depthStart = 0" << std::endl;
243 if(depthStart > maxDepth) {
244 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : depthStart too big ... = " 245 << depthStart << std::endl;
248 if( depthStart < 0.) depthStart = 0.;
249 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : depthStart re-calculated = " 250 << depthStart << std::endl;
253 if( onECAL &&
e < emid ) {
254 if((depthECAL - depthStart)/depthECAL > 0.2 && depthECAL >
depthStep ) {
258 LogDebug(
"FastCalorimetry") <<
" FamosHFShower : small energy, " 259 <<
" depthStart reduced to = " << depthStart << std::endl;
265 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : depthHCAL too small ... = " 266 << depthHCAL <<
" depthStart -> forced to 0 !!!" 273 LogInfo(
"FastCalorimetry") <<
" FamosHFShower : too small ECAL and HCAL depths - " 274 <<
" particle is lost !!! " << std::endl;
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;
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" 303 lamstep.push_back(depthECAL-depthStart);
308 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : " <<
" in ECAL sum1, sum2 " 309 << sum1 <<
" " << sum2 << std::endl;
320 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : goto HCAL" << std::endl;
322 depthStart = depthToHCAL;
328 if(depthStart >= depthECAL && depthStart < depthToHCAL ) {
329 depthStart = depthToHCAL;
333 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : goto HCAL" << std::endl;
337 if(
debug)
LogDebug(
"FastCalorimetry") <<
" FamosHFShower : forward" << std::endl;
341 for (
int i = 0;
i < nmoresteps ;
i++) {
343 if (sum1 > (depthECAL + depthGAP + depthHCAL))
break;
370 std::vector<double>
temp;
374 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::makeSteps - " 375 <<
" nsteps required : " << nsteps << std::endl;
378 for (
int i = 0;
i < nsteps;
i++) {
383 double y =
betHD * deplam;
386 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::makeSteps " 388 <<
" depx0, x = " << depx0 <<
", " << x
389 <<
" deplam, y = " << deplam <<
", " 398 LogDebug(
"FastCalorimetry") <<
"*** FamosHFShower::makeSteps " <<
" - negative step energy !!!" 409 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::makeSteps - nPoints estimate = " 410 << nPest << std::endl;
412 if(nPest <= 1 && count !=0 )
break;
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++)
433 LogDebug(
"FastCalorimetry") <<
"*** FamosHFShower::makeSteps " <<
" ECAL fraction : old/new - " 434 << oldECALenergy/sumes <<
"/" << newECALenergy/sumes << std::endl;
436 temp[0] = newECALenergy;
437 double newHCALenergy = sumes - newECALenergy;
438 double newHCALreweight = newHCALenergy / oldHCALenergy;
441 temp[
i] *= newHCALreweight;
450 eStep.push_back(temp[
i] *
e / sumes );
455 LogDebug(
"FastCalorimetry") <<
i <<
" xO and lamdepth at the end of step = " 469 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::makeSteps - " <<
"ECAL energy = " <<
eStep[0]
470 <<
" out of total = " <<
e << std::endl;
480 int numLongit =
eStep.size();
482 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - " 483 <<
" N_long.steps required : " << numLongit << std::endl;
489 std::vector<double> Fhist;
490 std::vector<double> rhist;
491 for (
int j = 0; j <
nTRsteps + 1; j++) {
495 LogDebug(
"FastCalorimetry") <<
"indexFinder - i, Fhist[i] = " << j <<
" " << Fhist[j] << std::endl;
501 for (
int i = 0;
i < numLongit ;
i++) {
507 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - detector = " 509 <<
" currentDepthL0 = " 510 << currentDepthL0 << std::endl;
513 double rbinsize = maxTRsize /
nTRsteps;
516 if( espot > 4. || espot < 0. )
517 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - unphysical espot = " 518 << espot << std::endl;
525 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - status of " 526 <<
" theHcalHitMaker->setDepth(currentDepthL0) is " 527 << setHDdepth << std::endl;
534 if(!setHDdepth)
continue;
544 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - status of Grid = " 545 << status << std::endl;
547 if(!status)
continue;
563 double eremaining =
eStep[
i];
564 bool converged =
false;
566 while (eremaining > 0. && !converged && count<inf ) {
571 double newespot = espot;
577 if( currentDepthL0 < 1.3 )
582 if ( layer == 2 ) newespot = 2. * espot;
584 if ( eremaining - newespot < 0. ) newespot = eremaining;
595 LogDebug(
"FastCalorimetry") << std::endl <<
" FamosHFShower::compute " <<
" r = " << radius
596 <<
" phi = " << phi << std::endl;
608 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - " 609 <<
" theGrid->addHit result = " 610 << result << std::endl;
618 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute - " 619 <<
" theHcalHitMaker->addHit result = " 620 << result << std::endl;
626 eremaining -= newespot;
627 if ( eremaining <= 0. ) converged =
true;
643 LogDebug(
"FastCalorimetry") <<
"*** FamosHFShower::compute " <<
" maximum number of" 644 <<
" transverse points " << count <<
" is used !!!" << std::endl;
649 LogDebug(
"FastCalorimetry") <<
" FamosHFShower::compute " <<
" long.step No." <<
i 650 <<
" Ntry, Nok = " << count
651 <<
" " << nok << std::endl;
661 int size = Fhist.size();
669 for (iter = 0; iter <
size ; iter++) {
671 if( curr >= size || curr < 1 )
672 LogWarning(
"FastCalorimetry") <<
" FamosHFShower::indexFinder - wrong current index = " 673 << curr <<
" !!!" << std::endl;
675 if ((x <= Fhist[curr]) && (x > Fhist[curr-1]))
break;
677 if(x > Fhist[curr]) {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;}}
685 LogDebug(
"FastCalorimetry") <<
" indexFinder - end of iter." << iter
686 <<
" curr, F[curr-1], F[curr] = " 687 << curr <<
" " << Fhist[curr-1] <<
" " << Fhist[curr] << std::endl;
692 LogDebug(
"FastCalorimetry") <<
" indexFinder x = " << x <<
" found index = " << curr-1
std::vector< double > lamdepth
std::vector< double > eStep
const RandomEngineAndDistribution * random
std::vector< double > rlamStep
bool addHit(double r, double phi, unsigned layer=0) override
std::vector< int > nspots
double flatShoot(double xmin=0.0, double xmax=1.0) const
int getHDnTRsteps() const
double getHDeSpotSize() const
std::vector< double > lamstep
int indexFinder(double x, const std::vector< double > &Fhist)
double radLenIncm() const override
Radiation length in cm.
double getHDhcalDepthFactor() const
std::vector< int > detector
double interactionLength() const override
Interaction length in cm.
double radLenIncm() const override
Radiation length in cm.
const ECALProperties * theECALproperties
int getHDlossesOpt() const
double interactionLength() const override
Interaction length in cm: 18.5 for Standard ECAL.
bool compute()
Compute the shower longitudinal and lateral development.
bool addHit(double r, double phi, unsigned layer=0) override
add the hit in the HCAL in local coordinates
const HSParameters * hsParameters() const
int getHDnDepthSteps() const
double gam(double x, double a) const
HcalHitMaker * theHcalHitMaker
std::vector< double > x0curr
double transProb(double factor, double R, double r)
void setSpotEnergy(double e) override
Set the spot energy.
std::vector< double > lamcurr
double getHDdepthStep() const
const HCALProperties * theHCALproperties
HFShower(const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
double getHDbalanceEH() const
double ecalHcalGapTotalX0() const
ECAL-HCAL transition.
void makeSteps(int nsteps)
std::vector< double > lamtotal
bool getPads(double depth, bool inCm=false)
const HCALProperties * hcalProperties() const
double getHDcriticalEnergy() const
const ECALProperties * ecalProperties() const
double getHDmaxTRfactor() const
double getHDtransParam() const
double ecalTotalL0() const
in the ECAL
HDShowerParametrization * theParam
bool setDepth(double, bool inCm=false)
set the depth in X0 or Lambda0 units depending on showerType
void setSpotEnergy(double e) override
double ecalHcalGapTotalL0() const
ECAL-HCAL transition.
std::vector< double > x0depth
double hcalTotalL0() const
in the HCAL