Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-29 06:08:20 

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