Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-23 23:48:32 

0001 //updated by Reza Goldouzian
0002 //FastSimulation Headers
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 // CMSSW headers
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 
0012 #include <cmath>
0013 
0014 // number attempts for transverse distribution if exit on a spec. condition
0015 #define infinity 10000
0016 // debugging flag ( 0, 1, 2, 3)
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       //    pmip(pmip),
0034       random(engine) {
0035   init(engine, myParam, myGrid, myHcalHitMaker, onECAL, epart, pmip);
0036 }
0037 
0038 void HDShower::init(const RandomEngineAndDistribution* engine,
0039                     HDShowerParametrization* myParam,
0040                     EcalHitMaker* myGrid,
0041                     HcalHitMaker* myHcalHitMaker,
0042                     int onECAL,
0043                     double epart,
0044                     double pmip) {
0045   // To get an access to constants read in FASTCalorimeter
0046   //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();
0047 
0048   // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
0049   lossesOpt = myParam->hsParameters()->getHDlossesOpt();
0050   nDepthSteps = myParam->hsParameters()->getHDnDepthSteps();
0051   nTRsteps = myParam->hsParameters()->getHDnTRsteps();
0052   transParam = myParam->hsParameters()->getHDtransParam();
0053   eSpotSize = myParam->hsParameters()->getHDeSpotSize();
0054   depthStep = myParam->hsParameters()->getHDdepthStep();
0055   criticalEnergy = myParam->hsParameters()->getHDcriticalEnergy();
0056   maxTRfactor = myParam->hsParameters()->getHDmaxTRfactor();
0057   balanceEH = myParam->hsParameters()->getHDbalanceEH();
0058   hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();
0059 
0060   // Special tr.size fluctuations
0061   transParam *= (1. + random->flatShoot());
0062 
0063   // Special long. fluctuations
0064   hcalDepthFactor += 0.05 * (2. * random->flatShoot() - 1.);
0065 
0066   transFactor = 1.;  // normally 1, in HF - might be smaller
0067                      // to take into account
0068                      // a narrowness of the HF shower (Cherenkov light)
0069 
0070   // simple protection ...
0071   if (e < 0)
0072     e = 0.;
0073 
0074   // Get the Famos Histos pointer
0075   //  myHistos = FamosHistos::instance();
0076   //  std::cout << " Hello FamosShower " << std::endl;
0077 
0078   theECALproperties = theParam->ecalProperties();
0079   theHCALproperties = theParam->hcalProperties();
0080 
0081   double emax = theParam->emax();
0082   double emid = theParam->emid();
0083   double emin = theParam->emin();
0084   double effective = e;
0085 
0086   if (e < emid) {
0087     theParam->setCase(1);
0088     // avoid "underflow" below Emin (for parameters calculation only)
0089     if (e < emin)
0090       effective = emin;
0091   } else
0092     theParam->setCase(2);
0093 
0094   // Avoid "overflow" beyond Emax (for parameters)
0095   if (effective > 0.5 * emax) {
0096     eSpotSize *= 2.5;
0097     if (effective > emax) {
0098       effective = emax;
0099       eSpotSize *= 4.;
0100       depthStep *= 2.;
0101       if (effective > 1000.)
0102         eSpotSize *= 2.;
0103     }
0104   }
0105 
0106   if (debug == 2)
0107     LogInfo("FastCalorimetry") << " HDShower : " << std::endl
0108                                << "       Energy   " << e << std::endl
0109                                << "      lossesOpt " << lossesOpt << std::endl
0110                                << "    nDepthSteps " << nDepthSteps << std::endl
0111                                << "       nTRsteps " << nTRsteps << std::endl
0112                                << "     transParam " << transParam << std::endl
0113                                << "      eSpotSize " << eSpotSize << std::endl
0114                                << " criticalEnergy " << criticalEnergy << std::endl
0115                                << "    maxTRfactor " << maxTRfactor << std::endl
0116                                << "      balanceEH " << balanceEH << std::endl
0117                                << "hcalDepthFactor " << hcalDepthFactor << std::endl;
0118 
0119   double alpEM1 = theParam->alpe1();
0120   double alpEM2 = theParam->alpe2();
0121 
0122   double betEM1 = theParam->bete1();
0123   double betEM2 = theParam->bete2();
0124 
0125   double alpHD1 = theParam->alph1();
0126   double alpHD2 = theParam->alph2();
0127 
0128   double betHD1 = theParam->beth1();
0129   double betHD2 = theParam->beth2();
0130 
0131   double part1 = theParam->part1();
0132   double part2 = theParam->part2();
0133 
0134   aloge = std::log(effective);
0135 
0136   double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
0137   double aedep = std::log(edpar);
0138 
0139   if (debug == 2)
0140     LogInfo("FastCalorimetry") << " HDShower : " << std::endl
0141                                << "     edpar " << edpar << "   aedep " << aedep << std::endl
0142                                << "    alpEM1 " << alpEM1 << std::endl
0143                                << "    alpEM2 " << alpEM2 << std::endl
0144                                << "    betEM1 " << betEM1 << std::endl
0145                                << "    betEM2 " << betEM2 << std::endl
0146                                << "    alpHD1 " << alpHD1 << std::endl
0147                                << "    alpHD2 " << alpHD2 << std::endl
0148                                << "    betHD1 " << betHD1 << std::endl
0149                                << "    betHD2 " << betHD2 << std::endl
0150                                << "     part1 " << part1 << std::endl
0151                                << "     part2 " << part2 << std::endl;
0152 
0153   // private members to set
0154   theR1 = theParam->r1();
0155   theR2 = theParam->r2();
0156   theR3 = theParam->r3();
0157 
0158   alpEM = alpEM1 + alpEM2 * aedep;
0159   tgamEM = tgamma(alpEM);
0160   betEM = betEM1 - betEM2 * aedep;
0161   alpHD = alpHD1 + alpHD2 * aedep;
0162   tgamHD = tgamma(alpHD);
0163   betHD = betHD1 - betHD2 * aedep;
0164   part = part1 - part2 * aedep;
0165   if (part > 1.)
0166     part = 1.;  // protection - just in case of
0167 
0168   if (debug == 2)
0169     LogInfo("FastCalorimetry") << " HDShower : " << std::endl
0170                                << "    alpEM " << alpEM << std::endl
0171                                << "   tgamEM " << tgamEM << std::endl
0172                                << "    betEM " << betEM << std::endl
0173                                << "    alpHD " << alpHD << std::endl
0174                                << "   tgamHD " << tgamHD << std::endl
0175                                << "    betHD " << betHD << std::endl
0176                                << "     part " << part << std::endl;
0177 
0178   if (onECAL) {
0179     lambdaEM = theParam->ecalProperties()->interactionLength();
0180     x0EM = theParam->ecalProperties()->radLenIncm();
0181   } else {
0182     lambdaEM = 0.;
0183     x0EM = 0.;
0184   }
0185   lambdaHD = theParam->hcalProperties()->interactionLength();
0186   x0HD = theParam->hcalProperties()->radLenIncm();
0187 
0188   if (debug == 2)
0189     LogInfo("FastCalorimetry") << " HDShower e " << e << std::endl
0190                                << "          x0EM = " << x0EM << std::endl
0191                                << "          x0HD = " << x0HD << std::endl
0192                                << "         lamEM = " << lambdaEM << std::endl
0193                                << "         lamHD = " << lambdaHD << std::endl;
0194 
0195   // Starting point of the shower
0196   // try first with ECAL lambda
0197 
0198   double sum1 = 0.;  // lambda path from the ECAL/HF entrance;
0199   double sum2 = 0.;  // lambda path from the interaction point;
0200   double sum3 = 0.;  // x0     path from the interaction point;
0201   int nsteps = 0;    // full number of longitudinal steps (counter);
0202 
0203   int nmoresteps;  // how many longitudinal steps in addition to
0204                    // one (if interaction happens there) in ECAL
0205 
0206   mip = 1;  // just to initiate particle as MIP in ECAL
0207 
0208   if (e < criticalEnergy)
0209     nmoresteps = 1;
0210   else
0211     nmoresteps = nDepthSteps;
0212 
0213   depthECAL = 0.;
0214   depthGAP = 0.;
0215   depthGAPx0 = 0.;
0216   if (onECAL) {
0217     depthECAL = theGrid->ecalTotalL0();  // ECAL depth segment
0218         //depthECAL  = theGrid->ecalTotalL0() + theGrid->ps1TotalL0() + theGrid->ps2TotalL0() + theGrid->ps2eeTotalL0(); //TEST: include preshower
0219     depthGAP = theGrid->ecalHcalGapTotalL0();    // GAP  depth segment
0220     depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
0221   }
0222 
0223   depthHCAL = theGrid->hcalTotalL0();  // HCAL depth segment
0224   depthToHCAL = depthECAL + depthGAP;
0225 
0226   //---------------------------------------------------------------------------
0227   // Depth simulation & various protections, among them
0228   // if too deep - get flat random in the allowed region
0229   // if no HCAL material behind - force to deposit in ECAL
0230   double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
0231   double depthStart = std::log(1. / random->flatShoot());  // starting point lambda unts
0232   //std::cout<<"generated depth             "<<depthStart<<std::endl;
0233   if (pmip == 1) {
0234     depthStart = depthToHCAL;
0235   } else {
0236     depthStart = depthStart * 0.9 * depthECAL / std::log(1. / pmip);
0237   }
0238   // std::cout<<"modified  depth             "<<depthStart<<std::endl;
0239 
0240   if (e < emin) {
0241     if (debug)
0242       LogInfo("FastCalorimetry") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl;
0243     depthStart = 0.;
0244   }
0245 
0246   if (depthStart > maxDepth) {
0247     if (debug)
0248       LogInfo("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " << depthStart << std::endl;
0249     depthStart = maxDepth * random->flatShoot();
0250     if (depthStart < 0.)
0251       depthStart = 0.;
0252     if (debug)
0253       LogInfo("FastCalorimetry") << " FamosHDShower : depthStart re-calculated = " << depthStart << std::endl;
0254   }
0255 
0256   if (onECAL && e < emid) {
0257     if (depthECAL > depthStep && (depthECAL - depthStart) / depthECAL > 0.2) {
0258       depthStart = 0.5 * depthECAL * random->flatShoot();
0259       if (debug)
0260         LogInfo("FastCalorimetry") << " FamosHDShower : small energy, "
0261                                    << " depthStart reduced to = " << depthStart << std::endl;
0262     }
0263   }
0264 
0265   if (depthHCAL < depthStep) {
0266     if (debug)
0267       LogInfo("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " << depthHCAL
0268                                  << " depthStart -> forced to 0 !!!" << std::endl;
0269     depthStart = 0.;
0270     nmoresteps = 0;
0271 
0272     if (depthECAL < depthStep) {
0273       nsteps = -1;
0274       LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - "
0275                                  << " particle is lost !!! " << std::endl;
0276     }
0277   }
0278 
0279   if (debug)
0280     LogInfo("FastCalorimetry") << " FamosHDShower  depths(lam) - " << std::endl
0281                                << "          ECAL = " << depthECAL << std::endl
0282                                << "           GAP = " << depthGAP << std::endl
0283                                << "          HCAL = " << depthHCAL << std::endl
0284                                << "  starting point = " << depthStart << std::endl;
0285 
0286   if (onEcal) {
0287     if (debug)
0288       LogInfo("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
0289     if (depthStart < depthECAL) {
0290       if (debug)
0291         LogInfo("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
0292       if (depthECAL > depthStep && (depthECAL - depthStart) / depthECAL > 0.1) {
0293         if (debug)
0294           LogInfo("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step" << std::endl;
0295 
0296         //  ECAL - one step
0297         nsteps++;
0298         sum1 += depthECAL;  // at the end of step
0299         sum2 += depthECAL - depthStart;
0300         sum3 += sum2 * lambdaEM / x0EM;
0301         lamtotal.push_back(sum1);
0302         lamdepth.push_back(sum2);
0303         lamcurr.push_back(lambdaEM);
0304         lamstep.push_back(depthECAL - depthStart);
0305         x0depth.push_back(sum3);
0306         x0curr.push_back(x0EM);
0307         detector.push_back(1);
0308         mip = 0;
0309 
0310         if (debug)
0311           LogInfo("FastCalorimetry") << " FamosHDShower : "
0312                                      << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;
0313 
0314         // Gap - no additional step after ECAL
0315         // just move further to HCAL over the gap
0316         sum1 += depthGAP;
0317         sum2 += depthGAP;
0318         sum3 += depthGAPx0;
0319       } else {  // Just shift starting point to HCAL
0320         //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
0321         if (debug)
0322           LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
0323 
0324         depthStart = depthToHCAL;
0325         sum1 += depthStart;
0326       }
0327     } else {  // GAP or HCAL
0328       if (depthStart >= depthECAL && depthStart < depthToHCAL) {
0329         depthStart = depthToHCAL;  // just a shift to HCAL for simplicity
0330       }
0331       sum1 += depthStart;
0332       if (debug)
0333         LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
0334     }
0335   } else {  // Forward
0336     if (debug)
0337       LogInfo("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
0338     sum1 += depthStart;
0339     transFactor = 0.5;  // makes narower tresverse size of shower
0340   }
0341 
0342   for (int i = 0; i < nmoresteps; i++) {
0343     sum1 += depthStep;
0344     if (sum1 > (depthECAL + depthGAP + depthHCAL))
0345       break;
0346     sum2 += depthStep;
0347     sum3 += sum2 * lambdaHD / x0HD;
0348     lamtotal.push_back(sum1);
0349     lamdepth.push_back(sum2);
0350     lamcurr.push_back(lambdaHD);
0351     lamstep.push_back(depthStep);
0352     x0depth.push_back(sum3);
0353     x0curr.push_back(x0HD);
0354     detector.push_back(3);
0355     nsteps++;
0356   }
0357 
0358   // Make fractions of energy and transverse radii at each step
0359 
0360   if (nsteps > 0) {
0361     makeSteps(nsteps);
0362   }
0363 }
0364 
0365 void HDShower::makeSteps(int nsteps) {
0366   double sumes = 0.;
0367   double sum = 0.;
0368   std::vector<double> temp;
0369 
0370   if (debug)
0371     LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - "
0372                                << " nsteps required : " << nsteps << std::endl;
0373 
0374   int count = 0;
0375   for (int i = 0; i < nsteps; i++) {
0376     double deplam = lamdepth[i] - 0.5 * lamstep[i];
0377     double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i];
0378     double x = betEM * depx0;
0379     double y = betHD * deplam;
0380 
0381     if (debug == 2)
0382       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps "
0383                                  << " - step " << i << "   depx0, x = " << depx0 << ", " << x
0384                                  << "   deplam, y = " << deplam << ", " << y << std::endl;
0385 
0386     double est = (part * betEM * gam(x, alpEM) * lamcurr[i] / (x0curr[i] * tgamEM) +
0387                   (1. - part) * betHD * gam(y, alpHD) / tgamHD) *
0388                  lamstep[i];
0389 
0390     // protection ...
0391     if (est < 0.) {
0392       LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps "
0393                                  << " - negative step energy !!!" << std::endl;
0394       break;
0395     }
0396 
0397     // for estimates only
0398     sum += est;
0399     int nPest = (int)(est * e / sum / eSpotSize);
0400 
0401     if (debug == 2)
0402       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " << nPest << std::endl;
0403 
0404     if (nPest <= 1 && count != 0)
0405       break;
0406 
0407     // good step - to proceed
0408 
0409     temp.push_back(est);
0410     sumes += est;
0411 
0412     rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor);
0413     count++;
0414   }
0415 
0416   // fluctuations in ECAL and re-distribution of remaining energy in HCAL
0417   if (detector[0] == 1 && count > 1) {
0418     double oldECALenergy = temp[0];
0419     double oldHCALenergy = sumes - oldECALenergy;
0420     double newECALenergy = 2. * sumes;
0421     for (int i = 0; newECALenergy > sumes && i < infinity; i++)
0422       newECALenergy = 2. * balanceEH * random->flatShoot() * oldECALenergy;
0423 
0424     if (debug == 2)
0425       LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps "
0426                                  << " ECAL fraction : old/new - " << oldECALenergy / sumes << "/"
0427                                  << newECALenergy / sumes << std::endl;
0428 
0429     temp[0] = newECALenergy;
0430     double newHCALenergy = sumes - newECALenergy;
0431     double newHCALreweight = newHCALenergy / oldHCALenergy;
0432 
0433     for (int i = 1; i < count; i++) {
0434       temp[i] *= newHCALreweight;
0435     }
0436   }
0437 
0438   // final re-normalization of the energy fractions
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 
0443     if (debug)
0444       LogInfo("FastCalorimetry") << " step " << i << "  det: " << detector[i]
0445                                  << "  xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
0446                                  << "   Estep func = " << eStep[i] << "   Rstep = " << rlamStep[i]
0447                                  << "  Nspots = " << nspots[i] << "  espot = " << eStep[i] / (double)nspots[i]
0448                                  << std::endl;
0449   }
0450 
0451   // The only step is in ECAL - let's make the size bigger ...
0452   if (count == 1 and detector[0] == 1)
0453     rlamStep[0] *= 2.;
0454 
0455   if (debug) {
0456     if (eStep[0] > 0.95 * e && detector[0] == 1)
0457       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - "
0458                                  << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
0459   }
0460 }
0461 
0462 bool HDShower::compute() {
0463   //  TimeMe theT("FamosHDShower::compute");
0464 
0465   bool status = false;
0466   int numLongit = eStep.size();
0467   if (debug)
0468     LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0469                                << " N_long.steps required : " << numLongit << std::endl;
0470 
0471   if (numLongit > 0) {
0472     status = true;
0473     // Prepare the trsanverse probability function
0474     std::vector<double> Fhist;
0475     std::vector<double> rhist;
0476     for (int j = 0; j < nTRsteps + 1; j++) {
0477       rhist.push_back(maxTRfactor * j / nTRsteps);
0478       Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
0479       if (debug == 3)
0480         LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
0481     }
0482 
0483     // Longitudinal steps
0484     for (int i = 0; i < numLongit; i++) {
0485       double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
0486       // vary the longitudinal profile if needed
0487       if (detector[i] != 1)
0488         currentDepthL0 *= hcalDepthFactor;
0489       if (debug)
0490         LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i]
0491                                    << "  currentDepthL0 = " << currentDepthL0 << std::endl;
0492 
0493       double maxTRsize = maxTRfactor * rlamStep[i];  // in lambda units
0494       double rbinsize = maxTRsize / nTRsteps;
0495       double espot = eStep[i] / (double)nspots[i];  // re-adjust espot
0496 
0497       if (espot > 2. || espot < 0.)
0498         LogInfo("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " << espot << std::endl;
0499 
0500       int ecal = 0;
0501       if (detector[i] != 1) {
0502         bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0503 
0504         if (debug)
0505           LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of "
0506                                      << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;
0507 
0508         if (!setHDdepth) {
0509           currentDepthL0 -= lamstep[i];
0510           setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0511         }
0512 
0513         if (!setHDdepth)
0514           continue;
0515 
0516         theHcalHitMaker->setSpotEnergy(espot);
0517 
0518         //fill hcal longitudinal distribution histogram
0519       } else {
0520         ecal = 1;
0521         bool status = theGrid->getPads(currentDepthL0);
0522 
0523         if (debug)
0524           LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;
0525 
0526         if (!status)
0527           continue;
0528 
0529         int ntry = nspots[i] * 10;
0530         if (ntry >= infinity) {  // use max allowed in case of too many spots
0531           nspots[i] = 0.5 * infinity;
0532           espot *= 0.1 * (double)ntry / double(nspots[i]);
0533         } else {
0534           espot *= 0.1;  // fine-grain energy spots in ECAL
0535                          // to avoid false ECAL clustering
0536           nspots[i] = ntry;
0537         }
0538 
0539         theGrid->setSpotEnergy(espot);
0540 
0541         //fill ecal longitudinal distribution histogram
0542       }
0543 
0544       // Transverse distribition
0545       int nok = 0;  // counter of OK
0546       int count = 0;
0547       int inf = infinity;
0548       if (lossesOpt)
0549         inf = nspots[i];  // if losses are enabled, otherwise
0550       // only OK points are counted ...
0551       if (nspots[i] > inf)
0552         std::cout << " FamosHDShower::compute - at long.step " << i << "  too many spots required : " << nspots[i]
0553                   << " !!! " << std::endl;
0554 
0555       for (int j = 0; j < inf; j++) {
0556         if (nok == nspots[i])
0557           break;
0558         count++;
0559 
0560         double prob = random->flatShoot();
0561         int index = indexFinder(prob, Fhist);
0562         double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;  // in-bin
0563         double phi = 2. * M_PI * random->flatShoot();
0564 
0565         if (debug == 2)
0566           LogInfo("FastCalorimetry") << std::endl
0567                                      << " FamosHDShower::compute "
0568                                      << " r = " << radius << "    phi = " << phi << std::endl;
0569 
0570         bool result;
0571         if (ecal) {
0572           result = theGrid->addHit(radius, phi, 0);
0573 
0574           if (debug == 2)
0575             LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0576                                        << " theGrid->addHit result = " << result << std::endl;
0577 
0578           //fill ecal transverse distribution histogram
0579         } else {
0580           result = theHcalHitMaker->addHit(radius, phi, 0);
0581 
0582           if (debug == 2)
0583             LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
0584                                        << " theHcalHitMaker->addHit result = " << result << std::endl;
0585 
0586           //fill hcal transverse distribution histogram
0587         }
0588 
0589         if (result)
0590           nok++;
0591 
0592       }  // end of tranverse simulation
0593 
0594       if (count == infinity) {
0595         if (debug)
0596           LogInfo("FastCalorimetry") << " FamosHDShower::compute "
0597                                      << " maximum number of"
0598                                      << " transverse points " << count << " is used !!!" << std::endl;
0599       }
0600 
0601       if (debug)
0602         LogInfo("FastCalorimetry") << " FamosHDShower::compute "
0603                                    << " long.step No." << i << "   Ntry, Nok = " << count << " " << nok << std::endl;
0604     }  // end of longitudinal steps
0605   }  // end of no steps
0606 
0607   return status;
0608 }
0609 
0610 int HDShower::indexFinder(double x, const std::vector<double>& Fhist) {
0611   // binary search in the vector of doubles
0612   int size = Fhist.size();
0613 
0614   int curr = size / 2;
0615   int step = size / 4;
0616   int iter;
0617   int prevdir = 0;
0618   int actudir = 0;
0619 
0620   for (iter = 0; iter < size; iter++) {
0621     if (curr >= size || curr < 1)
0622       LogWarning("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " << curr << " !!!"
0623                                     << std::endl;
0624 
0625     if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
0626       break;
0627     prevdir = actudir;
0628     if (x > Fhist[curr]) {
0629       actudir = 1;
0630     } else {
0631       actudir = -1;
0632     }
0633     if (prevdir * actudir < 0) {
0634       if (step > 1)
0635         step /= 2;
0636     }
0637     curr += actudir * step;
0638     if (curr > size)
0639       curr = size;
0640     else {
0641       if (curr < 1) {
0642         curr = 1;
0643       }
0644     }
0645 
0646     if (debug == 3)
0647       LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
0648                                  << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
0649   }
0650 
0651   if (debug == 3)
0652     LogInfo("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr - 1 << std::endl;
0653 
0654   return curr - 1;
0655 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!