Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-23 07:40:34 

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   for (int i = 0; i < count; i++) {
0439     eStep.push_back(temp[i] * e / sumes);
0440     nspots.push_back((int)(eStep[i] / eSpotSize) + 1);
0441 
0442     if (debug)
0443       LogDebug("FastCalorimetry") << i << "  xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
0444                                   << "   Estep func = " << eStep[i] << "   Rstep = " << rlamStep[i]
0445                                   << "  Nspots = " << nspots[i] << std::endl;
0446   }
0447 
0448   // The only step is in ECAL - let's make the size bigger ...
0449   if (count == 1 and detector[0] == 1)
0450     rlamStep[0] *= 2.;
0451 
0452   if (debug) {
0453     if (eStep[0] > 0.95 * e && detector[0] == 1)
0454       LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - "
0455                                   << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
0456   }
0457 }
0458 
0459 bool HFShower::compute() {
0460   //  TimeMe theT("FamosHFShower::compute");
0461 
0462   bool status = false;
0463   int numLongit = eStep.size();
0464   if (debug)
0465     LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0466                                 << " N_long.steps required : " << numLongit << std::endl;
0467 
0468   if (numLongit > 0) {
0469     status = true;
0470     // Prepare the trsanverse probability function
0471     std::vector<double> Fhist;
0472     std::vector<double> rhist;
0473     for (int j = 0; j < nTRsteps + 1; j++) {
0474       rhist.push_back(maxTRfactor * j / nTRsteps);
0475       Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
0476       if (debug == 3)
0477         LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
0478     }
0479 
0480     //================================================================
0481     // Longitudinal steps
0482     //================================================================
0483     for (int i = 0; i < numLongit; i++) {
0484       double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
0485       // vary the longitudinal profile if needed
0486       if (detector[i] != 1)
0487         currentDepthL0 *= hcalDepthFactor;
0488       if (debug)
0489         LogDebug("FastCalorimetry") << " FamosHFShower::compute - detector = " << detector[i]
0490                                     << "  currentDepthL0 = " << currentDepthL0 << std::endl;
0491 
0492       double maxTRsize = maxTRfactor * rlamStep[i];  // in lambda units
0493       double rbinsize = maxTRsize / nTRsteps;
0494       double espot = eStep[i] / (double)nspots[i];  // re-adjust espot
0495 
0496       if (espot > 4. || espot < 0.)
0497         LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = " << espot << std::endl;
0498 
0499       int ecal = 0;
0500       if (detector[i] != 1) {
0501         bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0502 
0503         if (debug)
0504           LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of "
0505                                       << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;
0506 
0507         if (!setHDdepth) {
0508           currentDepthL0 -= lamstep[i];
0509           setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
0510         }
0511         if (!setHDdepth)
0512           continue;
0513 
0514         theHcalHitMaker->setSpotEnergy(espot);
0515       } else {
0516         ecal = 1;
0517         bool status = theGrid->getPads(currentDepthL0);
0518 
0519         if (debug)
0520           LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = " << status << std::endl;
0521 
0522         if (!status)
0523           continue;
0524 
0525         theGrid->setSpotEnergy(espot);
0526       }
0527 
0528       //------------------------------------------------------------
0529       // Transverse distribution
0530       //------------------------------------------------------------
0531       int nok = 0;  // counter of OK
0532       int count = 0;
0533       int inf = infinity;
0534       if (lossesOpt)
0535         inf = nspots[i];  // losses are enabled, otherwise
0536       // only OK points are counted ...
0537 
0538       // Total energy in this layer
0539       double eremaining = eStep[i];
0540       bool converged = false;
0541 
0542       while (eremaining > 0. && !converged && count < inf) {
0543         ++count;
0544 
0545         // energy spot  (HFL)
0546         double newespot = espot;
0547 
0548         // We need to know a priori if this energy spot if for
0549         // a long (1) or short (2) fiber
0550 
0551         unsigned layer = 1;
0552         if (currentDepthL0 < 1.3)  // first 22 cm = 1.3 lambda - only HFL
0553           layer = 1;
0554         else
0555           layer = random->flatShoot() < 0.5 ? 1 : 2;
0556 
0557         if (layer == 2)
0558           newespot = 2. * espot;
0559 
0560         if (eremaining - newespot < 0.)
0561           newespot = eremaining;
0562 
0563         // process transverse distribution
0564 
0565         double prob = random->flatShoot();
0566         int index = indexFinder(prob, Fhist);
0567         double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;  // in-bin
0568         double phi = 2. * M_PI * random->flatShoot();
0569 
0570         if (debug == 2)
0571           LogDebug("FastCalorimetry") << std::endl
0572                                       << " FamosHFShower::compute "
0573                                       << " r = " << radius << "    phi = " << phi << std::endl;
0574 
0575         // add hit
0576 
0577         theHcalHitMaker->setSpotEnergy(newespot);
0578         theGrid->setSpotEnergy(newespot);
0579 
0580         bool result;
0581         if (ecal) {
0582           result = theGrid->addHit(radius, phi, 0);  // shouldn't get here !
0583 
0584           if (debug == 2)
0585             LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0586                                         << " theGrid->addHit result = " << result << std::endl;
0587 
0588         } else {
0589           // PV assign espot to long/short fibers
0590           result = theHcalHitMaker->addHit(radius, phi, layer);
0591 
0592           if (debug == 2)
0593             LogDebug("FastCalorimetry") << " FamosHFShower::compute - "
0594                                         << " theHcalHitMaker->addHit result = " << result << std::endl;
0595         }
0596 
0597         if (result) {
0598           ++nok;
0599           eremaining -= newespot;
0600           if (eremaining <= 0.)
0601             converged = true;
0602           //      std::cout << "Transverse : " << nok << " "
0603           //            << " , E= "     << newespot
0604           //            << " , Erem = " << eremaining
0605           //            << std::endl;
0606         } else {
0607           //      std::cout << "WARNING : hit not added" << std::endl;
0608         }
0609       }
0610 
0611       // end of tranverse simulation
0612       //-----------------------------------------------------
0613 
0614       if (count == infinity) {
0615         status = false;
0616         if (debug)
0617           LogDebug("FastCalorimetry") << "*** FamosHFShower::compute "
0618                                       << " maximum number of"
0619                                       << " transverse points " << count << " is used !!!" << std::endl;
0620         break;
0621       }
0622 
0623       if (debug)
0624         LogDebug("FastCalorimetry") << " FamosHFShower::compute "
0625                                     << " long.step No." << i << "   Ntry, Nok = " << count << " " << nok << std::endl;
0626 
0627     }  // end of longitudinal steps
0628   }    // end of no steps
0629   return status;
0630 }
0631 
0632 int HFShower::indexFinder(double x, const std::vector<double>& Fhist) {
0633   // binary search in the vector of doubles
0634   int size = Fhist.size();
0635 
0636   int curr = size / 2;
0637   int step = size / 4;
0638   int iter;
0639   int prevdir = 0;
0640   int actudir = 0;
0641 
0642   for (iter = 0; iter < size; iter++) {
0643     if (curr >= size || curr < 1)
0644       LogWarning("FastCalorimetry") << " FamosHFShower::indexFinder - wrong current index = " << curr << " !!!"
0645                                     << std::endl;
0646 
0647     if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
0648       break;
0649     prevdir = actudir;
0650     if (x > Fhist[curr]) {
0651       actudir = 1;
0652     } else {
0653       actudir = -1;
0654     }
0655     if (prevdir * actudir < 0) {
0656       if (step > 1)
0657         step /= 2;
0658     }
0659     curr += actudir * step;
0660     if (curr > size)
0661       curr = size;
0662     else {
0663       if (curr < 1) {
0664         curr = 1;
0665       }
0666     }
0667 
0668     if (debug == 3)
0669       LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
0670                                   << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
0671   }
0672 
0673   if (debug == 3)
0674     LogDebug("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr - 1 << std::endl;
0675 
0676   return curr - 1;
0677 }
[ 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!