Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:26:46 

0001 //FastSimulation Headers
0002 #include "FastSimulation/ShowerDevelopment/interface/HFShower.h"
0003 //#include "FastSimulation/Utilities/interface/Histos.h"
0004 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0005 
0006 //////////////////////////////////////////////////////////////////////
0007 // What's this?
0008 //#include "FastSimulation/FamosCalorimeters/interface/FASTCalorimeter.h"
0009 
0010 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0011 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0012 
0013 // CMSSW headers
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 ///////////////////////////////////////////////////////////////
0017 // And This???? Doesn't seem to be needed
0018 // #include "Calorimetry/CaloDetector/interface/CellGeometry.h"
0019 
0020 #include <cmath>
0021 
0022 // number attempts for transverse distribution if exit on a spec. condition
0023 #define infinity 5000
0024 // debugging flag ( 0, 1, 2, 3)
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   // To get an access to constants read in FASTCalorimeter
0037   //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();
0038 
0039   // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
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   // Special tr.size fluctuations
0052   transParam *= (1. + random->flatShoot());
0053 
0054   // Special ad hoc long. extension + some fluctuations
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   // normally 1, in HF - might be smaller to take into account
0067   // a narrowness of the HF shower (Cherenkov light)
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   // simple protection ...
0074   if (e < 0)
0075     e = 0.;
0076 
0077   // Get the Famos Histos pointer
0078   //  myHistos = FamosHistos::instance();
0079   //  std::cout << " Hello FamosShower " << std::endl;
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     // avoid "underflow" below Emin (for parameters calculation only)
0092     if (e < emin)
0093       effective = emin;
0094   } else
0095     theParam->setCase(2);
0096 
0097   // A bit coarse espot size for HF...
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   // private members to set
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.;  // protection - just in case of
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   // Starting point of the shower
0198   // try first with ECAL lambda
0199 
0200   double sum1 = 0.;  // lambda path from the ECAL/HF entrance;
0201   double sum2 = 0.;  // lambda path from the interaction point;
0202   double sum3 = 0.;  // x0     path from the interaction point;
0203   int nsteps = 0;    // full number of longitudinal steps (counter);
0204 
0205   int nmoresteps;  // how many longitudinal steps in addition to
0206                    // one (if interaction happens there) in ECAL
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();          // ECAL depth segment
0218     depthGAP = theGrid->ecalHcalGapTotalL0();    // GAP  depth segment
0219     depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
0220   }
0221 
0222   double depthHCAL = theGrid->hcalTotalL0();  // HCAL depth segment
0223   double depthToHCAL = depthECAL + depthGAP;
0224 
0225   //---------------------------------------------------------------------------
0226   // Depth simulation & various protections, among them
0227   // if too deep - get flat random in the allowed region
0228   // if no HCAL material behind - force to deposit in ECAL
0229   double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
0230 
0231   double depthStart = std::log(1. / random->flatShoot());  // starting point lambda unts
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         //  ECAL - one step
0290         nsteps++;
0291         sum1 += depthECAL;  // at the end of step
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         //                           // Gap - no additional step after ECAL
0307         //                           // just move further to HCAL over the gap
0308         sum1 += depthGAP;
0309         sum2 += depthGAP;
0310         sum3 += depthGAPx0;
0311       }
0312       // Just shift starting point to HCAL
0313       else {
0314         //  cout << " FamosHFShower : not enough space to make ECAL step" << std::endl;
0315         if (debug)
0316           LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
0317 
0318         depthStart = depthToHCAL;
0319         sum1 += depthStart;
0320       }
0321     } else {  // GAP or HCAL
0322 
0323       if (depthStart >= depthECAL && depthStart < depthToHCAL) {
0324         depthStart = depthToHCAL;  // just a shift to HCAL for simplicity
0325       }
0326       sum1 += depthStart;
0327 
0328       if (debug)
0329         LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
0330     }
0331   } else {  // Forward
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   // Make fractions of energy and transverse radii at each step
0354 
0355   // PV
0356   //  std::cout << "HFShower::HFShower() : Nsteps = " << nsteps << std::endl;
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     // protection ...
0389     if (est < 0.) {
0390       LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps "
0391                                   << " - negative step energy !!!" << std::endl;
0392       est = 0.;
0393       break;
0394     }
0395 
0396     // for estimates only
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     // good step - to proceed
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   // fluctuations in ECAL and re-distribution of remaining energy in HCAL
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   // final re-normalization of the energy fractions
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   // The only step is in ECAL - let's make the size bigger ...
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   //  TimeMe theT("FamosHFShower::compute");
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     // Prepare the trsanverse probability function
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     // Longitudinal steps
0484     //================================================================
0485     for (int i = 0; i < numLongit; i++) {
0486       double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
0487       // vary the longitudinal profile if needed
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];  // in lambda units
0495       double rbinsize = maxTRsize / nTRsteps;
0496       double espot = eStep[i] / (double)nspots[i];  // re-adjust espot
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       // Transverse distribution
0532       //------------------------------------------------------------
0533       int nok = 0;  // counter of OK
0534       int count = 0;
0535       int inf = infinity;
0536       if (lossesOpt)
0537         inf = nspots[i];  // losses are enabled, otherwise
0538       // only OK points are counted ...
0539 
0540       // Total energy in this layer
0541       double eremaining = eStep[i];
0542       bool converged = false;
0543 
0544       while (eremaining > 0. && !converged && count < inf) {
0545         ++count;
0546 
0547         // energy spot  (HFL)
0548         double newespot = espot;
0549 
0550         // We need to know a priori if this energy spot if for
0551         // a long (1) or short (2) fiber
0552 
0553         unsigned layer = 1;
0554         if (currentDepthL0 < 1.3)  // first 22 cm = 1.3 lambda - only HFL
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         // process transverse distribution
0566 
0567         double prob = random->flatShoot();
0568         int index = indexFinder(prob, Fhist);
0569         double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;  // in-bin
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         // add hit
0578 
0579         theHcalHitMaker->setSpotEnergy(newespot);
0580         theGrid->setSpotEnergy(newespot);
0581 
0582         bool result;
0583         if (ecal) {
0584           result = theGrid->addHit(radius, phi, 0);  // shouldn't get here !
0585 
0586           if (debug == 2)
0587             LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0588                                         << " theGrid->addHit result = " << result << std::endl;
0589 
0590         } else {
0591           // PV assign espot to long/short fibers
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           //      std::cout << "Transverse : " << nok << " "
0605           //            << " , E= "     << newespot
0606           //            << " , Erem = " << eremaining
0607           //            << std::endl;
0608         } else {
0609           //      std::cout << "WARNING : hit not added" << std::endl;
0610         }
0611       }
0612 
0613       // end of tranverse simulation
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     }  // end of longitudinal steps
0630   }    // end of no steps
0631   return status;
0632 }
0633 
0634 int HFShower::indexFinder(double x, const std::vector<double>& Fhist) {
0635   // binary search in the vector of doubles
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 }
[ 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!