Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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