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