File indexing completed on 2021-02-14 13:26:46
0001
0002 #include "FastSimulation/ShowerDevelopment/interface/HFShower.h"
0003
0004 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0005
0006
0007
0008
0009
0010 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0011 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0012
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016
0017
0018
0019
0020 #include <cmath>
0021
0022
0023 #define infinity 5000
0024
0025 #define debug 0
0026
0027 using namespace edm;
0028
0029 HFShower::HFShower(const RandomEngineAndDistribution* engine,
0030 HDShowerParametrization* myParam,
0031 EcalHitMaker* myGrid,
0032 HcalHitMaker* myHcalHitMaker,
0033 int onECAL,
0034 double epart)
0035 : theParam(myParam), theGrid(myGrid), theHcalHitMaker(myHcalHitMaker), onEcal(onECAL), e(epart), random(engine) {
0036
0037
0038
0039
0040 lossesOpt = myParam->hsParameters()->getHDlossesOpt();
0041 nDepthSteps = myParam->hsParameters()->getHDnDepthSteps();
0042 nTRsteps = myParam->hsParameters()->getHDnTRsteps();
0043 transParam = myParam->hsParameters()->getHDtransParam();
0044 eSpotSize = myParam->hsParameters()->getHDeSpotSize();
0045 depthStep = myParam->hsParameters()->getHDdepthStep();
0046 criticalEnergy = myParam->hsParameters()->getHDcriticalEnergy();
0047 maxTRfactor = myParam->hsParameters()->getHDmaxTRfactor();
0048 balanceEH = myParam->hsParameters()->getHDbalanceEH();
0049 hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();
0050
0051
0052 transParam *= (1. + random->flatShoot());
0053
0054
0055 double depthExt;
0056 if (e < 50.)
0057 depthExt = 0.8 * (50. - e) / 50. + 0.3;
0058 else {
0059 if (e < 500.)
0060 depthExt = (500. - e) / 500. * 0.4 - 0.1;
0061 else
0062 depthExt = -0.1;
0063 }
0064 hcalDepthFactor += depthExt + 0.05 * (2. * random->flatShoot() - 1.);
0065
0066
0067
0068 if (e < 50.)
0069 transFactor = 0.5 - (50. - e) / 50. * 0.2;
0070 else
0071 transFactor = 0.7 - (1000. - e) / 1000. * 0.2;
0072
0073
0074 if (e < 0)
0075 e = 0.;
0076
0077
0078
0079
0080
0081 theECALproperties = theParam->ecalProperties();
0082 theHCALproperties = theParam->hcalProperties();
0083
0084 double emax = theParam->emax();
0085 double emid = theParam->emid();
0086 double emin = theParam->emin();
0087 double effective = e;
0088
0089 if (e < emid) {
0090 theParam->setCase(1);
0091
0092 if (e < emin)
0093 effective = emin;
0094 } else
0095 theParam->setCase(2);
0096
0097
0098 eSpotSize *= 2.5;
0099 if (effective > 0.5 * emax) {
0100 eSpotSize *= 2.;
0101 if (effective > emax) {
0102 effective = emax;
0103 eSpotSize *= 3.;
0104 depthStep *= 2.;
0105 }
0106 }
0107
0108 if (debug == 2)
0109 LogDebug("FastCalorimetry") << " HFShower : " << std::endl
0110 << " Energy " << e << std::endl
0111 << " lossesOpt " << lossesOpt << std::endl
0112 << " nDepthSteps " << nDepthSteps << std::endl
0113 << " nTRsteps " << nTRsteps << std::endl
0114 << " transParam " << transParam << std::endl
0115 << " eSpotSize " << eSpotSize << std::endl
0116 << " criticalEnergy " << criticalEnergy << std::endl
0117 << " maxTRfactor " << maxTRfactor << std::endl
0118 << " balanceEH " << balanceEH << std::endl
0119 << "hcalDepthFactor " << hcalDepthFactor << std::endl;
0120
0121 double alpEM1 = theParam->alpe1();
0122 double alpEM2 = theParam->alpe2();
0123
0124 double betEM1 = theParam->bete1();
0125 double betEM2 = theParam->bete2();
0126
0127 double alpHD1 = theParam->alph1();
0128 double alpHD2 = theParam->alph2();
0129
0130 double betHD1 = theParam->beth1();
0131 double betHD2 = theParam->beth2();
0132
0133 double part1 = theParam->part1();
0134 double part2 = theParam->part2();
0135
0136 aloge = std::log(effective);
0137
0138 double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
0139 double aedep = std::log(edpar);
0140
0141 if (debug == 2)
0142 LogDebug("FastCalorimetry") << " HFShower : " << std::endl
0143 << " edpar " << edpar << " aedep " << aedep << std::endl
0144 << " alpEM1 " << alpEM1 << std::endl
0145 << " alpEM2 " << alpEM2 << std::endl
0146 << " betEM1 " << betEM1 << std::endl
0147 << " betEM2 " << betEM2 << std::endl
0148 << " alpHD1 " << alpHD1 << std::endl
0149 << " alpHD2 " << alpHD2 << std::endl
0150 << " betHD1 " << betHD1 << std::endl
0151 << " betHD2 " << betHD2 << std::endl
0152 << " part1 " << part1 << std::endl
0153 << " part2 " << part2 << std::endl;
0154
0155
0156 theR1 = theParam->r1();
0157 theR2 = theParam->r2();
0158 theR3 = theParam->r3();
0159
0160 alpEM = alpEM1 + alpEM2 * aedep;
0161 tgamEM = tgamma(alpEM);
0162 betEM = betEM1 - betEM2 * aedep;
0163 alpHD = alpHD1 + alpHD2 * aedep;
0164 tgamHD = tgamma(alpHD);
0165 betHD = betHD1 - betHD2 * aedep;
0166 part = part1 - part2 * aedep;
0167 if (part > 1.)
0168 part = 1.;
0169
0170 if (debug == 2)
0171 LogDebug("FastCalorimetry") << " HFShower : " << std::endl
0172 << " alpEM " << alpEM << std::endl
0173 << " tgamEM " << tgamEM << std::endl
0174 << " betEM " << betEM << std::endl
0175 << " alpHD " << alpHD << std::endl
0176 << " tgamHD " << tgamHD << std::endl
0177 << " betHD " << betHD << std::endl
0178 << " part " << part << std::endl;
0179
0180 if (onECAL) {
0181 lambdaEM = theParam->ecalProperties()->interactionLength();
0182 x0EM = theParam->ecalProperties()->radLenIncm();
0183 } else {
0184 lambdaEM = 0.;
0185 x0EM = 0.;
0186 }
0187 lambdaHD = theParam->hcalProperties()->interactionLength();
0188 x0HD = theParam->hcalProperties()->radLenIncm();
0189
0190 if (debug == 2)
0191 LogDebug("FastCalorimetry") << " HFShower e " << e << std::endl
0192 << " x0EM = " << x0EM << std::endl
0193 << " x0HD = " << x0HD << std::endl
0194 << " lamEM = " << lambdaEM << std::endl
0195 << " lamHD = " << lambdaHD << std::endl;
0196
0197
0198
0199
0200 double sum1 = 0.;
0201 double sum2 = 0.;
0202 double sum3 = 0.;
0203 int nsteps = 0;
0204
0205 int nmoresteps;
0206
0207
0208 if (e < criticalEnergy)
0209 nmoresteps = 1;
0210 else
0211 nmoresteps = nDepthSteps;
0212
0213 double depthECAL = 0.;
0214 double depthGAP = 0.;
0215 double depthGAPx0 = 0.;
0216 if (onECAL) {
0217 depthECAL = theGrid->ecalTotalL0();
0218 depthGAP = theGrid->ecalHcalGapTotalL0();
0219 depthGAPx0 = theGrid->ecalHcalGapTotalX0();
0220 }
0221
0222 double depthHCAL = theGrid->hcalTotalL0();
0223 double depthToHCAL = depthECAL + depthGAP;
0224
0225
0226
0227
0228
0229 double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
0230
0231 double depthStart = std::log(1. / random->flatShoot());
0232
0233 if (e < emin) {
0234 if (debug)
0235 LogDebug("FastCalorimetry") << " FamosHFShower : e <emin -> depthStart = 0" << std::endl;
0236 depthStart = 0.;
0237 }
0238
0239 if (depthStart > maxDepth) {
0240 if (debug)
0241 LogDebug("FastCalorimetry") << " FamosHFShower : depthStart too big ... = " << depthStart << std::endl;
0242
0243 depthStart = maxDepth * random->flatShoot();
0244 if (depthStart < 0.)
0245 depthStart = 0.;
0246 if (debug)
0247 LogDebug("FastCalorimetry") << " FamosHFShower : depthStart re-calculated = " << depthStart << std::endl;
0248 }
0249
0250 if (onECAL && e < emid) {
0251 if ((depthECAL - depthStart) / depthECAL > 0.2 && depthECAL > depthStep) {
0252 depthStart = 0.5 * depthECAL * random->flatShoot();
0253 if (debug)
0254 LogDebug("FastCalorimetry") << " FamosHFShower : small energy, "
0255 << " depthStart reduced to = " << depthStart << std::endl;
0256 }
0257 }
0258
0259 if (depthHCAL < depthStep) {
0260 if (debug)
0261 LogDebug("FastCalorimetry") << " FamosHFShower : depthHCAL too small ... = " << depthHCAL
0262 << " depthStart -> forced to 0 !!!" << std::endl;
0263 depthStart = 0.;
0264 nmoresteps = 0;
0265
0266 if (depthECAL < depthStep) {
0267 nsteps = -1;
0268 LogInfo("FastCalorimetry") << " FamosHFShower : too small ECAL and HCAL depths - "
0269 << " particle is lost !!! " << std::endl;
0270 }
0271 }
0272
0273 if (debug)
0274 LogDebug("FastCalorimetry") << " FamosHFShower depths(lam) - " << std::endl
0275 << " ECAL = " << depthECAL << std::endl
0276 << " GAP = " << depthGAP << std::endl
0277 << " HCAL = " << depthHCAL << std::endl
0278 << " starting point = " << depthStart << std::endl;
0279
0280 if (onEcal) {
0281 if (debug)
0282 LogDebug("FastCalorimetry") << " FamosHFShower : onECAL" << std::endl;
0283 if (depthStart < depthECAL) {
0284 if (debug)
0285 LogDebug("FastCalorimetry") << " FamosHFShower : depthStart < depthECAL" << std::endl;
0286 if ((depthECAL - depthStart) / depthECAL > 0.25 && depthECAL > depthStep) {
0287 if (debug)
0288 LogDebug("FastCalorimetry") << " FamosHFShower : enough space to make ECAL step" << std::endl;
0289
0290 nsteps++;
0291 sum1 += depthECAL;
0292 sum2 += depthECAL - depthStart;
0293 sum3 += sum2 * lambdaEM / x0EM;
0294 lamtotal.push_back(sum1);
0295 lamdepth.push_back(sum2);
0296 lamcurr.push_back(lambdaEM);
0297 lamstep.push_back(depthECAL - depthStart);
0298 x0depth.push_back(sum3);
0299 x0curr.push_back(x0EM);
0300 detector.push_back(1);
0301
0302 if (debug)
0303 LogDebug("FastCalorimetry") << " FamosHFShower : "
0304 << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;
0305
0306
0307
0308 sum1 += depthGAP;
0309 sum2 += depthGAP;
0310 sum3 += depthGAPx0;
0311 }
0312
0313 else {
0314
0315 if (debug)
0316 LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
0317
0318 depthStart = depthToHCAL;
0319 sum1 += depthStart;
0320 }
0321 } else {
0322
0323 if (depthStart >= depthECAL && depthStart < depthToHCAL) {
0324 depthStart = depthToHCAL;
0325 }
0326 sum1 += depthStart;
0327
0328 if (debug)
0329 LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
0330 }
0331 } else {
0332 if (debug)
0333 LogDebug("FastCalorimetry") << " FamosHFShower : forward" << std::endl;
0334 sum1 += depthStart;
0335 }
0336
0337 for (int i = 0; i < nmoresteps; i++) {
0338 sum1 += depthStep;
0339 if (sum1 > (depthECAL + depthGAP + depthHCAL))
0340 break;
0341 sum2 += depthStep;
0342 sum3 += sum2 * lambdaHD / x0HD;
0343 lamtotal.push_back(sum1);
0344 lamdepth.push_back(sum2);
0345 lamcurr.push_back(lambdaHD);
0346 lamstep.push_back(depthStep);
0347 x0depth.push_back(sum3);
0348 x0curr.push_back(x0HD);
0349 detector.push_back(3);
0350 nsteps++;
0351 }
0352
0353
0354
0355
0356
0357
0358 if (nsteps > 0) {
0359 makeSteps(nsteps);
0360 }
0361 }
0362
0363 void HFShower::makeSteps(int nsteps) {
0364 double sumes = 0.;
0365 double sum = 0.;
0366 std::vector<double> temp;
0367
0368 if (debug)
0369 LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
0370 << " nsteps required : " << nsteps << std::endl;
0371
0372 int count = 0;
0373 for (int i = 0; i < nsteps; i++) {
0374 double deplam = lamdepth[i] - 0.5 * lamstep[i];
0375 double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i];
0376 double x = betEM * depx0;
0377 double y = betHD * deplam;
0378
0379 if (debug == 2)
0380 LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps "
0381 << " - step " << i << " depx0, x = " << depx0 << ", " << x
0382 << " deplam, y = " << deplam << ", " << y << std::endl;
0383
0384 double est = (part * betEM * gam(x, alpEM) * lamcurr[i] / (x0curr[i] * tgamEM) +
0385 (1. - part) * betHD * gam(y, alpHD) / tgamHD) *
0386 lamstep[i];
0387
0388
0389 if (est < 0.) {
0390 LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps "
0391 << " - negative step energy !!!" << std::endl;
0392 est = 0.;
0393 break;
0394 }
0395
0396
0397 sum += est;
0398 int nPest = (int)(est * e / sum / eSpotSize);
0399
0400 if (debug == 2)
0401 LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - nPoints estimate = " << nPest << std::endl;
0402
0403 if (nPest <= 1 && count != 0)
0404 break;
0405
0406
0407
0408 temp.push_back(est);
0409 sumes += est;
0410
0411 rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor);
0412 count++;
0413 }
0414
0415
0416 if (detector[0] == 1 && count > 1) {
0417 double oldECALenergy = temp[0];
0418 double oldHCALenergy = sumes - oldECALenergy;
0419 double newECALenergy = 2. * sumes;
0420 for (int i = 0; newECALenergy > sumes && i < infinity; i++)
0421 newECALenergy = 2. * balanceEH * random->flatShoot() * oldECALenergy;
0422
0423 if (debug == 2)
0424 LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps "
0425 << " ECAL fraction : old/new - " << oldECALenergy / sumes << "/"
0426 << newECALenergy / sumes << std::endl;
0427
0428 temp[0] = newECALenergy;
0429 double newHCALenergy = sumes - newECALenergy;
0430 double newHCALreweight = newHCALenergy / oldHCALenergy;
0431
0432 for (int i = 1; i < count; i++) {
0433 temp[i] *= newHCALreweight;
0434 }
0435 }
0436
0437
0438 double etot = 0.;
0439 for (int i = 0; i < count; i++) {
0440 eStep.push_back(temp[i] * e / sumes);
0441 nspots.push_back((int)(eStep[i] / eSpotSize) + 1);
0442 etot += eStep[i];
0443
0444 if (debug)
0445 LogDebug("FastCalorimetry") << i << " xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
0446 << " Estep func = " << eStep[i] << " Rstep = " << rlamStep[i]
0447 << " Nspots = " << nspots[i] << std::endl;
0448 }
0449
0450
0451 if (count == 1 and detector[0] == 1)
0452 rlamStep[0] *= 2.;
0453
0454 if (debug) {
0455 if (eStep[0] > 0.95 * e && detector[0] == 1)
0456 LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
0457 << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
0458 }
0459 }
0460
0461 bool HFShower::compute() {
0462
0463
0464 bool status = false;
0465 int numLongit = eStep.size();
0466 if (debug)
0467 LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0468 << " N_long.steps required : " << numLongit << std::endl;
0469
0470 if (numLongit > 0) {
0471 status = true;
0472
0473 std::vector<double> Fhist;
0474 std::vector<double> rhist;
0475 for (int j = 0; j < nTRsteps + 1; j++) {
0476 rhist.push_back(maxTRfactor * j / nTRsteps);
0477 Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
0478 if (debug == 3)
0479 LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
0480 }
0481
0482
0483
0484
0485 for (int i = 0; i < numLongit; i++) {
0486 double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
0487
0488 if (detector[i] != 1)
0489 currentDepthL0 *= hcalDepthFactor;
0490 if (debug)
0491 LogDebug("FastCalorimetry") << " FamosHFShower::compute - detector = " << detector[i]
0492 << " currentDepthL0 = " << currentDepthL0 << std::endl;
0493
0494 double maxTRsize = maxTRfactor * rlamStep[i];
0495 double rbinsize = maxTRsize / nTRsteps;
0496 double espot = eStep[i] / (double)nspots[i];
0497
0498 if (espot > 4. || espot < 0.)
0499 LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = " << espot << std::endl;
0500
0501 int ecal = 0;
0502 if (detector[i] != 1) {
0503 bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0504
0505 if (debug)
0506 LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of "
0507 << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;
0508
0509 if (!setHDdepth) {
0510 currentDepthL0 -= lamstep[i];
0511 setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0512 }
0513 if (!setHDdepth)
0514 continue;
0515
0516 theHcalHitMaker->setSpotEnergy(espot);
0517 } else {
0518 ecal = 1;
0519 bool status = theGrid->getPads(currentDepthL0);
0520
0521 if (debug)
0522 LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = " << status << std::endl;
0523
0524 if (!status)
0525 continue;
0526
0527 theGrid->setSpotEnergy(espot);
0528 }
0529
0530
0531
0532
0533 int nok = 0;
0534 int count = 0;
0535 int inf = infinity;
0536 if (lossesOpt)
0537 inf = nspots[i];
0538
0539
0540
0541 double eremaining = eStep[i];
0542 bool converged = false;
0543
0544 while (eremaining > 0. && !converged && count < inf) {
0545 ++count;
0546
0547
0548 double newespot = espot;
0549
0550
0551
0552
0553 unsigned layer = 1;
0554 if (currentDepthL0 < 1.3)
0555 layer = 1;
0556 else
0557 layer = random->flatShoot() < 0.5 ? 1 : 2;
0558
0559 if (layer == 2)
0560 newespot = 2. * espot;
0561
0562 if (eremaining - newespot < 0.)
0563 newespot = eremaining;
0564
0565
0566
0567 double prob = random->flatShoot();
0568 int index = indexFinder(prob, Fhist);
0569 double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;
0570 double phi = 2. * M_PI * random->flatShoot();
0571
0572 if (debug == 2)
0573 LogDebug("FastCalorimetry") << std::endl
0574 << " FamosHFShower::compute "
0575 << " r = " << radius << " phi = " << phi << std::endl;
0576
0577
0578
0579 theHcalHitMaker->setSpotEnergy(newespot);
0580 theGrid->setSpotEnergy(newespot);
0581
0582 bool result;
0583 if (ecal) {
0584 result = theGrid->addHit(radius, phi, 0);
0585
0586 if (debug == 2)
0587 LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0588 << " theGrid->addHit result = " << result << std::endl;
0589
0590 } else {
0591
0592 result = theHcalHitMaker->addHit(radius, phi, layer);
0593
0594 if (debug == 2)
0595 LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0596 << " theHcalHitMaker->addHit result = " << result << std::endl;
0597 }
0598
0599 if (result) {
0600 ++nok;
0601 eremaining -= newespot;
0602 if (eremaining <= 0.)
0603 converged = true;
0604
0605
0606
0607
0608 } else {
0609
0610 }
0611 }
0612
0613
0614
0615
0616 if (count == infinity) {
0617 status = false;
0618 if (debug)
0619 LogDebug("FastCalorimetry") << "*** FamosHFShower::compute "
0620 << " maximum number of"
0621 << " transverse points " << count << " is used !!!" << std::endl;
0622 break;
0623 }
0624
0625 if (debug)
0626 LogDebug("FastCalorimetry") << " FamosHFShower::compute "
0627 << " long.step No." << i << " Ntry, Nok = " << count << " " << nok << std::endl;
0628
0629 }
0630 }
0631 return status;
0632 }
0633
0634 int HFShower::indexFinder(double x, const std::vector<double>& Fhist) {
0635
0636 int size = Fhist.size();
0637
0638 int curr = size / 2;
0639 int step = size / 4;
0640 int iter;
0641 int prevdir = 0;
0642 int actudir = 0;
0643
0644 for (iter = 0; iter < size; iter++) {
0645 if (curr >= size || curr < 1)
0646 LogWarning("FastCalorimetry") << " FamosHFShower::indexFinder - wrong current index = " << curr << " !!!"
0647 << std::endl;
0648
0649 if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
0650 break;
0651 prevdir = actudir;
0652 if (x > Fhist[curr]) {
0653 actudir = 1;
0654 } else {
0655 actudir = -1;
0656 }
0657 if (prevdir * actudir < 0) {
0658 if (step > 1)
0659 step /= 2;
0660 }
0661 curr += actudir * step;
0662 if (curr > size)
0663 curr = size;
0664 else {
0665 if (curr < 1) {
0666 curr = 1;
0667 }
0668 }
0669
0670 if (debug == 3)
0671 LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
0672 << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
0673 }
0674
0675 if (debug == 3)
0676 LogDebug("FastCalorimetry") << " indexFinder x = " << x << " found index = " << curr - 1 << std::endl;
0677
0678 return curr - 1;
0679 }