File indexing completed on 2024-10-29 06:08: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 break;
0385 }
0386
0387
0388 sum += est;
0389 int nPest = (int)(est * e / sum / eSpotSize);
0390
0391 if (debug == 2)
0392 LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " << nPest << std::endl;
0393
0394 if (nPest <= 1 && count != 0)
0395 break;
0396
0397
0398
0399 temp.push_back(est);
0400 sumes += est;
0401
0402 rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor);
0403 count++;
0404 }
0405
0406
0407 if (detector[0] == 1 && count > 1) {
0408 double oldECALenergy = temp[0];
0409 double oldHCALenergy = sumes - oldECALenergy;
0410 double newECALenergy = 2. * sumes;
0411 for (int i = 0; newECALenergy > sumes && i < infinity; i++)
0412 newECALenergy = 2. * balanceEH * random->flatShoot() * oldECALenergy;
0413
0414 if (debug == 2)
0415 LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps "
0416 << " ECAL fraction : old/new - " << oldECALenergy / sumes << "/"
0417 << newECALenergy / sumes << std::endl;
0418
0419 temp[0] = newECALenergy;
0420 double newHCALenergy = sumes - newECALenergy;
0421 double newHCALreweight = newHCALenergy / oldHCALenergy;
0422
0423 for (int i = 1; i < count; i++) {
0424 temp[i] *= newHCALreweight;
0425 }
0426 }
0427
0428
0429 for (int i = 0; i < count; i++) {
0430 eStep.push_back(temp[i] * e / sumes);
0431 nspots.push_back((int)(eStep[i] / eSpotSize) + 1);
0432
0433 if (debug)
0434 LogInfo("FastCalorimetry") << " step " << i << " det: " << detector[i]
0435 << " xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
0436 << " Estep func = " << eStep[i] << " Rstep = " << rlamStep[i]
0437 << " Nspots = " << nspots[i] << " espot = " << eStep[i] / (double)nspots[i]
0438 << std::endl;
0439 }
0440
0441
0442 if (count == 1 and detector[0] == 1)
0443 rlamStep[0] *= 2.;
0444
0445 if (debug) {
0446 if (eStep[0] > 0.95 * e && detector[0] == 1)
0447 LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - "
0448 << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
0449 }
0450 }
0451
0452 bool HDShower::compute() {
0453
0454
0455 bool status = false;
0456 int numLongit = eStep.size();
0457 if (debug)
0458 LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0459 << " N_long.steps required : " << numLongit << std::endl;
0460
0461 if (numLongit > 0) {
0462 status = true;
0463
0464 std::vector<double> Fhist;
0465 std::vector<double> rhist;
0466 for (int j = 0; j < nTRsteps + 1; j++) {
0467 rhist.push_back(maxTRfactor * j / nTRsteps);
0468 Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
0469 if (debug == 3)
0470 LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
0471 }
0472
0473
0474 for (int i = 0; i < numLongit; i++) {
0475 double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
0476
0477 if (detector[i] != 1)
0478 currentDepthL0 *= hcalDepthFactor;
0479 if (debug)
0480 LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i]
0481 << " currentDepthL0 = " << currentDepthL0 << std::endl;
0482
0483 double maxTRsize = maxTRfactor * rlamStep[i];
0484 double rbinsize = maxTRsize / nTRsteps;
0485 double espot = eStep[i] / (double)nspots[i];
0486
0487 if (espot > 2. || espot < 0.)
0488 LogInfo("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " << espot << std::endl;
0489
0490 int ecal = 0;
0491 if (detector[i] != 1) {
0492 bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0493
0494 if (debug)
0495 LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of "
0496 << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;
0497
0498 if (!setHDdepth) {
0499 currentDepthL0 -= lamstep[i];
0500 setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0501 }
0502
0503 if (!setHDdepth)
0504 continue;
0505
0506 theHcalHitMaker->setSpotEnergy(espot);
0507
0508
0509 } else {
0510 ecal = 1;
0511 bool status = theGrid->getPads(currentDepthL0);
0512
0513 if (debug)
0514 LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;
0515
0516 if (!status)
0517 continue;
0518
0519 int ntry = nspots[i] * 10;
0520 if (ntry >= infinity) {
0521 nspots[i] = 0.5 * infinity;
0522 espot *= 0.1 * (double)ntry / double(nspots[i]);
0523 } else {
0524 espot *= 0.1;
0525
0526 nspots[i] = ntry;
0527 }
0528
0529 theGrid->setSpotEnergy(espot);
0530
0531
0532 }
0533
0534
0535 int nok = 0;
0536 int count = 0;
0537 int inf = infinity;
0538 if (lossesOpt)
0539 inf = nspots[i];
0540
0541 if (nspots[i] > inf)
0542 std::cout << " FamosHDShower::compute - at long.step " << i << " too many spots required : " << nspots[i]
0543 << " !!! " << std::endl;
0544
0545 for (int j = 0; j < inf; j++) {
0546 if (nok == nspots[i])
0547 break;
0548 count++;
0549
0550 double prob = random->flatShoot();
0551 int index = indexFinder(prob, Fhist);
0552 double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;
0553 double phi = 2. * M_PI * random->flatShoot();
0554
0555 if (debug == 2)
0556 LogInfo("FastCalorimetry") << std::endl
0557 << " FamosHDShower::compute "
0558 << " r = " << radius << " phi = " << phi << std::endl;
0559
0560 bool result;
0561 if (ecal) {
0562 result = theGrid->addHit(radius, phi, 0);
0563
0564 if (debug == 2)
0565 LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0566 << " theGrid->addHit result = " << result << std::endl;
0567
0568
0569 } else {
0570 result = theHcalHitMaker->addHit(radius, phi, 0);
0571
0572 if (debug == 2)
0573 LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0574 << " theHcalHitMaker->addHit result = " << result << std::endl;
0575
0576
0577 }
0578
0579 if (result)
0580 nok++;
0581
0582 }
0583
0584 if (count == infinity) {
0585 if (debug)
0586 LogInfo("FastCalorimetry") << " FamosHDShower::compute "
0587 << " maximum number of"
0588 << " transverse points " << count << " is used !!!" << std::endl;
0589 }
0590
0591 if (debug)
0592 LogInfo("FastCalorimetry") << " FamosHDShower::compute "
0593 << " long.step No." << i << " Ntry, Nok = " << count << " " << nok << std::endl;
0594 }
0595 }
0596
0597 return status;
0598 }
0599
0600 int HDShower::indexFinder(double x, const std::vector<double>& Fhist) {
0601
0602 int size = Fhist.size();
0603
0604 int curr = size / 2;
0605 int step = size / 4;
0606 int iter;
0607 int prevdir = 0;
0608 int actudir = 0;
0609
0610 for (iter = 0; iter < size; iter++) {
0611 if (curr >= size || curr < 1)
0612 LogWarning("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " << curr << " !!!"
0613 << std::endl;
0614
0615 if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
0616 break;
0617 prevdir = actudir;
0618 if (x > Fhist[curr]) {
0619 actudir = 1;
0620 } else {
0621 actudir = -1;
0622 }
0623 if (prevdir * actudir < 0) {
0624 if (step > 1)
0625 step /= 2;
0626 }
0627 curr += actudir * step;
0628 if (curr > size)
0629 curr = size;
0630 else {
0631 if (curr < 1) {
0632 curr = 1;
0633 }
0634 }
0635
0636 if (debug == 3)
0637 LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
0638 << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
0639 }
0640
0641 if (debug == 3)
0642 LogInfo("FastCalorimetry") << " indexFinder x = " << x << " found index = " << curr - 1 << std::endl;
0643
0644 return curr - 1;
0645 }