Project CMSSW displayed by LXR CMSSW_15_0_X_2024-12-13-2300/​FastSimulation/​ShowerDevelopment/​src/​EMShower.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_15_0_X_2024-12-13-2300 CMSSW_15_0_X_2024-12-12-2300 CMSSW_15_0_X_2024-12-11-2300 CMSSW_15_0_X_2024-12-10-2300 CMSSW_15_0_X_2024-12-09-2300 CMSSW_15_0_X_2024-12-08-2300 CMSSW_15_0_X_2024-12-06-2300 Revert   Change

File indexing completed on 2024-04-06 12:11:19

 //FAMOS Headers
 #include "FastSimulation/ShowerDevelopment/interface/EMShower.h"
 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
 
 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
 
 #include <cmath>
 
 //#include "FastSimulation/Utilities/interface/Histos.h"
 
 using std::vector;
 
 EMShower::EMShower(const RandomEngineAndDistribution* engine,
                    GammaFunctionGenerator* gamma,
                    EMECALShowerParametrization* const myParam,
                    vector<const RawParticle*>* const myPart,
                    EcalHitMaker* const myGrid,
                    PreshowerHitMaker* const myPresh,
                    bool bFixedLength)
 
     : theParam(myParam),
       thePart(myPart),
       theGrid(myGrid),
       thePreshower(myPresh),
       random(engine),
       myGammaGenerator(gamma),
       bFixedLength_(bFixedLength) {
   // Get the Famos Histos pointer
   //  myHistos = Histos::instance();
   //  myGammaGenerator = GammaFunctionGenerator::instance();
   stepsCalculated = false;
   hasPreshower = myPresh != nullptr;
   theECAL = myParam->ecalProperties();
   theHCAL = myParam->hcalProperties();
   theLayer1 = myParam->layer1Properties();
   theLayer2 = myParam->layer2Properties();
 
   double fotos = theECAL->photoStatistics() * theECAL->lightCollectionEfficiency();
 
   nPart = thePart->size();
   totalEnergy = 0.;
   globalMaximum = 0.;
   //double meanDepth = 0.;
   // Initialize the shower parameters for each particle
 
   for (unsigned int i = 0; i < nPart; ++i) {
     //    std::cout << " AAA " << *(*thePart)[i] << std::endl;
     // The particle and the shower energy
     Etot.push_back(0.);
     E.push_back(((*thePart)[i])->e());
     totalEnergy += E[i];
 
     double lny = std::log(E[i] / theECAL->criticalEnergy());
 
     // Average and Sigma for T and alpha
     double theMeanT = myParam->meanT(lny);
     double theMeanAlpha = myParam->meanAlpha(lny);
     double theMeanLnT = myParam->meanLnT(lny);
     double theMeanLnAlpha = myParam->meanLnAlpha(lny);
     double theSigmaLnT = myParam->sigmaLnT(lny);
     double theSigmaLnAlpha = myParam->sigmaLnAlpha(lny);
 
     // The correlation matrix
     double theCorrelation = myParam->correlationAlphaT(lny);
     double rhop = std::sqrt((1. + theCorrelation) / 2.);
     double rhom = std::sqrt((1. - theCorrelation) / 2.);
 
     // The number of spots in ECAL / HCAL
     theNumberOfSpots.push_back(myParam->nSpots(E[i]));
     //    theNumberOfSpots.push_back(myParam->nSpots(E[i])*spotFraction);
     //theNumberOfSpots = random->poissonShoot(myParam->nSpots(myPart->e()));
 
     // Photo-statistics
     photos.push_back(E[i] * fotos);
 
     // The longitudinal shower development parameters
     // Fluctuations of alpha, T and beta
     double z1 = 0.;
     double z2 = 0.;
     double aa = 0.;
 
     // Protect against too large fluctuations (a < 1) for small energies
     while (aa <= 1.) {
       z1 = random->gaussShoot(0., 1.);
       z2 = random->gaussShoot(0., 1.);
       aa = std::exp(theMeanLnAlpha + theSigmaLnAlpha * (z1 * rhop - z2 * rhom));
     }
 
     a.push_back(aa);
     T.push_back(std::exp(theMeanLnT + theSigmaLnT * (z1 * rhop + z2 * rhom)));
     b.push_back((a[i] - 1.) / T[i]);
     maximumOfShower.push_back((a[i] - 1.) / b[i]);
     globalMaximum += maximumOfShower[i] * E[i];
     //meanDepth += a[i] / b[i] * E[i];
     //    std::cout << " Adding max " << maximumOfShower[i] << " " << E[i] << " " <<maximumOfShower[i]*E[i] << std::endl;
     //    std::cout << std::setw(8) << std::setprecision(5) << " a /b " << a[i] << " " << b[i] << std::endl;
     Ti.push_back(a[i] / b[i] * (std::exp(theMeanLnAlpha) - 1.) / std::exp(theMeanLnAlpha));
 
     // The parameters for the number of energy spots
     TSpot.push_back(theParam->meanTSpot(theMeanT));
     aSpot.push_back(theParam->meanAlphaSpot(theMeanAlpha));
     bSpot.push_back((aSpot[i] - 1.) / TSpot[i]);
     //    myHistos->fill("h7000",a[i]);
     //    myHistos->fill("h7002",E[i],a[i]);
   }
   //  std::cout << " PS1 : " << myGrid->ps1TotalX0()
   //         << " PS2 : " << myGrid->ps2TotalX0()
   //         << " ECAL : " << myGrid->ecalTotalX0()
   //         << " HCAL : " << myGrid->hcalTotalX0()
   //         << " Offset : " << myGrid->x0DepthOffset()
   //         << std::endl;
 
   globalMaximum /= totalEnergy;
   //meanDepth /= totalEnergy;
   // std::cout << " Total Energy " << totalEnergy << " Global max " << globalMaximum << std::endl;
 }
 
 void EMShower::prepareSteps() {
   //  TimeMe theT("EMShower::compute");
 
   // Determine the longitudinal intervals
   //  std::cout << " EMShower compute" << std::endl;
   double dt;
   double radlen;
   int stps;
   int first_Ecal_step = 0;
   int last_Ecal_step = 0;
 
   // The maximum is in principe 8 (with 5X0 steps in the ECAL)
   steps.reserve(24);
 
   /*  
   std::cout << " PS1 : " << theGrid->ps1TotalX0()
         << " PS2 : " << theGrid->ps2TotalX0()
         << " PS2 and ECAL : " << theGrid->ps2eeTotalX0()
         << " ECAL : " << theGrid->ecalTotalX0()
         << " HCAL : " << theGrid->hcalTotalX0() 
         << " Offset : " << theGrid->x0DepthOffset()
         << std::endl;
   */
 
   radlen = -theGrid->x0DepthOffset();
 
   // Preshower Layer 1
   radlen += theGrid->ps1TotalX0();
   if (radlen > 0.) {
     steps.push_back(Step(0, radlen));
     radlen = 0.;
   }
 
   // Preshower Layer 2
   radlen += theGrid->ps2TotalX0();
   if (radlen > 0.) {
     steps.push_back(Step(1, radlen));
     radlen = 0.;
   }
 
   // add a step between preshower and ee
   radlen += theGrid->ps2eeTotalX0();
   if (radlen > 0.) {
     steps.push_back(Step(5, radlen));
     radlen = 0.;
   }
 
   // ECAL
   radlen += theGrid->ecalTotalX0();
 
   //  std::cout << "theGrid->ecalTotalX0() = " << theGrid->ecalTotalX0() << std::endl;
 
   if (radlen > 0.) {
     if (!bFixedLength_) {
       stps = (int)((radlen + 2.5) / 5.);
       //    stps=(int)((radlen+.5)/1.);
       if (stps == 0)
         stps = 1;
       dt = radlen / (double)stps;
       Step step(2, dt);
       first_Ecal_step = steps.size();
       for (int ist = 0; ist < stps; ++ist)
         steps.push_back(step);
       last_Ecal_step = steps.size() - 1;
       radlen = 0.;
     } else {
       dt = 1.0;
       stps = static_cast<int>(radlen);
       if (stps == 0)
         stps = 1;
       Step step(2, dt);
       first_Ecal_step = steps.size();
       for (int ist = 0; ist < stps; ++ist)
         steps.push_back(step);
       dt = radlen - stps;
       if (dt > 0) {
         Step stepLast(2, dt);
         steps.push_back(stepLast);
       }
       last_Ecal_step = steps.size() - 1;
       //      std::cout << "radlen = "  << radlen << " stps = " << stps << " dt = " << dt << std::endl;
       radlen = 0.;
     }
   }
 
   // I should had a gap here !
 
   // HCAL
   radlen += theGrid->hcalTotalX0();
   if (radlen > 0.) {
     double dtFrontHcal = theGrid->totalX0() - theGrid->hcalTotalX0();
     // One single step for the full HCAL
     if (dtFrontHcal < 30.) {
       dt = 30. - dtFrontHcal;
       Step step(3, dt);
       steps.push_back(step);
     }
   }
 
   nSteps = steps.size();
   if (nSteps == 0)
     return;
   double ESliceTot = 0.;
   double MeanDepth = 0.;
   depositedEnergy.resize(nSteps);
   meanDepth.resize(nSteps);
   double t = 0.;
 
   int offset = 0;
   for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
     ESliceTot = 0.;
     MeanDepth = 0.;
     double realTotalEnergy = 0;
     dt = steps[iStep].second;
     t += dt;
     for (unsigned int i = 0; i < nPart; ++i) {
       depositedEnergy[iStep].push_back(deposit(t, a[i], b[i], dt));
       ESliceTot += depositedEnergy[iStep][i];
       MeanDepth += deposit(t, a[i] + 1., b[i], dt) / b[i] * a[i];
       realTotalEnergy += depositedEnergy[iStep][i] * E[i];
     }
 
     if (ESliceTot > 0.)  // can happen for the shower tails; this depth will be skipped anyway
       MeanDepth /= ESliceTot;
     else
       MeanDepth = t - dt;
 
     meanDepth[iStep] = MeanDepth;
     if (realTotalEnergy < 0.001) {
       offset -= 1;
     }
   }
 
   innerDepth = meanDepth[first_Ecal_step];
   if (last_Ecal_step + offset >= 0)
     outerDepth = meanDepth[last_Ecal_step + offset];
   else
     outerDepth = innerDepth;
 
   stepsCalculated = true;
 }
 
 void EMShower::compute() {
   double t = 0.;
   double dt = 0.;
   if (!stepsCalculated)
     prepareSteps();
 
   // Prepare the grids in EcalHitMaker
   // theGrid->setInnerAndOuterDepth(innerDepth,outerDepth);
   //float pstot = 0.;
   //float ps2tot = 0.;
   //float ps1tot = 0.;
   bool status = false;
   //  double E1 = 0.;  // Energy layer 1
   //  double E2 = 0.;  // Energy layer 2
   //  double n1 = 0.;  // #mips layer 1
   //  double n2 = 0.;  // #mips layer 2
   //  double E9 = 0.;  // Energy ECAL
 
   // Loop over all segments for the longitudinal development
   //double totECalc = 0;
 
   for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
     // The length of the shower in this segment
     dt = steps[iStep].second;
 
     //    std::cout << " Detector " << steps[iStep].first << " t " << t << " " << dt << std::endl;
 
     // The elapsed length
     t += dt;
 
     // In what detector are we ?
     unsigned detector = steps[iStep].first;
 
     bool presh1 = detector == 0;
     bool presh2 = detector == 1;
     bool ecal = detector == 2;
     bool hcal = detector == 3;
     bool vfcal = detector == 4;
     bool gap = detector == 5;
 
     // Temporary. Will be removed
     if (theHCAL == nullptr)
       hcal = false;
 
     // Keep only ECAL for now
     if (vfcal)
       continue;
 
     // Nothing to do in the gap
     if (gap)
       continue;
 
     //    cout << " t = " << t << endl;
     // Build the grid of crystals at this ECAL depth
     // Actually, it might be useful to check if this grid is empty or not.
     // If it is empty (because no crystal at this depth), it is of no use
     // (and time consuming) to generate the spots
 
     // middle of the step
     double tt = t - 0.5 * dt;
 
     double realTotalEnergy = 0.;
     for (unsigned int i = 0; i < nPart; ++i) {
       realTotalEnergy += depositedEnergy[iStep][i] * E[i];
     }
 
     //    std::cout << " Step " << tt << std::endl;
     //    std::cout << "ecal " << ecal << " hcal "  << hcal <<std::endl;
 
     // If the amount of energy is greater than 1 MeV, make a new grid
     // otherwise put in the previous one.
     bool usePreviousGrid = (realTotalEnergy < 0.001);
 
     // If the amount of energy is greater than 1 MeV, make a new grid
     // otherwise put in the previous one.
 
     // If less than 1 kEV. Just skip
     if (iStep > 2 && realTotalEnergy < 0.000001)
       continue;
 
     if (ecal && !usePreviousGrid) {
       status = theGrid->getPads(meanDepth[iStep]);
     }
     if (hcal) {
       status = theHcalHitMaker->setDepth(tt);
     }
     if ((ecal || hcal) && !status)
       continue;
 
     bool detailedShowerTail = false;
     // check if a detailed treatment of the rear leakage should be applied
     if (ecal && !usePreviousGrid) {
       detailedShowerTail = (t - dt > theGrid->getX0back());
     }
 
     // The particles of the shower are processed in parallel
     for (unsigned int i = 0; i < nPart; ++i) {
       //      double Edepo=deposit(t,a[i],b[i],dt);
 
       //  integration of the shower profile between t-dt and t
       double dE = (!hcal) ? depositedEnergy[iStep][i] : 1. - deposit(a[i], b[i], t - dt);
 
       // no need to do the full machinery if there is ~nothing to distribute)
       if (dE * E[i] < 0.000001)
         continue;
 
       if (ecal && !theECAL->isHom()) {
         double mean = dE * E[i];
         double sigma = theECAL->resE() * sqrt(mean);
 
         /*
       double meanLn = log(mean);
       double kLn = sigma/mean+1;
       double sigmaLn = log(kLn);
     */
 
         double dE0 = dE;
 
         //    std::cout << "dE before shoot = " << dE << std::endl;
         dE = random->gaussShoot(mean, sigma) / E[i];
 
         //    myGammaGenerator->setParameters(aSam,bSam,0);
         //    dE = myGammaGenerator->shoot()/E[i];
         //    std::cout << "dE shooted = " << dE << " E[i] = " << E[i] << std::endl;
         if (dE * E[i] < 0.000001)
           continue;
         photos[i] = photos[i] * dE / dE0;
       }
 
       /*
       if (ecal && !theParam->ecalProperties()->isHom()){
 
     double cSquare = TMath::Power(theParam->ecalProperties()->resE(),2);
     double aSam = dE/cSquare;
     double bSam = 1./cSquare;
 
     //  dE = dE*gam(bSam*dE, aSam)/tgamma(aSam);
       }
       */
 
       //totECalc += dE;
 
       // The number of energy spots (or mips)
       double nS = 0;
 
       // ECAL case : Account for photostatistics and long'al non-uniformity
       if (ecal) {
         //  double aSam = E[i]*dE*one_over_resoSquare;
         //  double bSam = one_over_resoSquare;
 
         dE = random->poissonShoot(dE * photos[i]) / photos[i];
         double z0 = random->gaussShoot(0., 1.);
         dE *= 1. + z0 * theECAL->lightCollectionUniformity();
 
         // Expected spot number
         nS = (theNumberOfSpots[i] * gam(bSpot[i] * tt, aSpot[i]) * bSpot[i] * dt / tgamma(aSpot[i]));
 
         // Preshower : Expected number of mips + fluctuation
       } else if (hcal) {
         nS = (theNumberOfSpots[i] * gam(bSpot[i] * tt, aSpot[i]) * bSpot[i] * dt / tgamma(aSpot[i])) *
              theHCAL->spotFraction();
         double nSo = nS;
         nS = random->poissonShoot(nS);
         // 'Quick and dirty' fix (but this line should be better removed):
         if (nSo > 0. && nS / nSo < 10.)
           dE *= nS / nSo;
 
         //  if(true)
         //    {
         //      std::cout << " theHCAL->spotFraction = " <<theHCAL->spotFraction() <<std::endl;
         //      std::cout << " nSpot Ecal : " << nSo/theHCAL->spotFraction() << " Final " << nS << std::endl;
         //    }
       } else if (presh1) {
         nS = random->poissonShoot(dE * E[i] * theLayer1->mipsPerGeV());
         //  std::cout << " dE *E[i] (1)" << dE*E[i] << " " << dE*E[i]*theLayer1->mipsPerGeV() << "  "<< nS << std::endl;
         //pstot += dE * E[i];
         //ps1tot += dE * E[i];
         dE = nS / (E[i] * theLayer1->mipsPerGeV());
 
         //        E1 += dE*E[i];
         //  n1 += nS;
         //  if (presh2) { E2 += SpotEnergy; ++n2; }
 
       } else if (presh2) {
         nS = random->poissonShoot(dE * E[i] * theLayer2->mipsPerGeV());
         //  std::cout << " dE *E[i] (2) " << dE*E[i] << " " << dE*E[i]*theLayer2->mipsPerGeV() << "  "<< nS << std::endl;
         //pstot += dE * E[i];
         //ps2tot += dE * E[i];
         dE = nS / (E[i] * theLayer2->mipsPerGeV());
 
         //        E2 += dE*E[i];
         //  n2 += nS;
       }
 
       if (detailedShowerTail)
         myGammaGenerator->setParameters(floor(a[i] + 0.5), b[i], t - dt);
 
       //    myHistos->fill("h100",t,dE);
 
       // The lateral development parameters
 
       // Energy of the spots
       double eSpot = (nS > 0.) ? dE / nS : 0.;
       double SpotEnergy = eSpot * E[i];
 
       if (hasPreshower && (presh1 || presh2))
         thePreshower->setSpotEnergy(0.00009);
       if (hcal) {
         SpotEnergy *= theHCAL->hOverPi();
         theHcalHitMaker->setSpotEnergy(SpotEnergy);
       }
       // Poissonian fluctuations for the number of spots
       //    int nSpot = random->poissonShoot(nS);
       int nSpot = (int)(nS + 0.5);
 
       // Fig. 11 (right) *** Does not match.
       //    myHistos->fill("h101",t,(double)nSpot/theNumberOfSpots);
 
       //double taui = t/T;
       double taui = tt / Ti[i];
       double proba = theParam->p(taui, E[i]);
       double theRC = theParam->rC(taui, E[i]);
       double theRT = theParam->rT(taui, E[i]);
 
       // Fig. 10
       //    myHistos->fill("h300",taui,theRC);
       //    myHistos->fill("h301",taui,theRT);
       //    myHistos->fill("h302",taui,proba);
 
       double dSpotsCore = random->gaussShoot(proba * nSpot, std::sqrt(proba * (1. - proba) * nSpot));
 
       if (dSpotsCore < 0)
         dSpotsCore = 0;
 
       unsigned nSpots_core = (unsigned)(dSpotsCore + 0.5);
       unsigned nSpots_tail = ((unsigned)nSpot > nSpots_core) ? nSpot - nSpots_core : 0;
 
       for (unsigned icomp = 0; icomp < 2; ++icomp) {
         double theR = (icomp == 0) ? theRC : theRT;
         unsigned ncompspots = (icomp == 0) ? nSpots_core : nSpots_tail;
 
         RadialInterval radInterval(theR, ncompspots, SpotEnergy, random);
         if (ecal) {
           if (icomp == 0) {
             setIntervals(icomp, radInterval);
           } else {
             setIntervals(icomp, radInterval);
           }
         } else {
           radInterval.addInterval(100., 1.);  // 100% of the spots
         }
 
         radInterval.compute();
         // irad = 0 : central circle; irad=1 : outside
 
         unsigned nrad = radInterval.nIntervals();
 
         for (unsigned irad = 0; irad < nrad; ++irad) {
           double spote = radInterval.getSpotEnergy(irad);
           if (ecal)
             theGrid->setSpotEnergy(spote);
           if (hcal)
             theHcalHitMaker->setSpotEnergy(spote);
           unsigned nradspots = radInterval.getNumberOfSpots(irad);
           double umin = radInterval.getUmin(irad);
           double umax = radInterval.getUmax(irad);
           // Go for the lateral development
           //           std::cout << "Couche " << iStep << " irad = " << irad << " Ene = " << E[i] << " eSpot = " << eSpot << " spote = " << spote << " nSpot = " << nS << std::endl;
 
           for (unsigned ispot = 0; ispot < nradspots; ++ispot) {
             double z3 = random->flatShoot(umin, umax);
             double ri = theR * std::sqrt(z3 / (1. - z3));
 
             // Generate phi
             double phi = 2. * M_PI * random->flatShoot();
 
             // Add the hit in the crystal
             //  if( ecal ) theGrid->addHit(ri*theECAL->moliereRadius(),phi);
             // Now the *moliereRadius is done in EcalHitMaker
             if (ecal) {
               if (detailedShowerTail) {
                 //             std::cout << "About to call addHitDepth " << std::endl;
                 double depth;
                 do {
                   depth = myGammaGenerator->shoot(random);
                 } while (depth > t);
                 theGrid->addHitDepth(ri, phi, depth);
                 //             std::cout << " Done " << std::endl;
               } else
                 theGrid->addHit(ri, phi);
             } else if (hasPreshower && presh1)
               thePreshower->addHit(ri, phi, 1);
             else if (hasPreshower && presh2)
               thePreshower->addHit(ri, phi, 2);
             else if (hcal) {
               //               std::cout << " About to add a spot in the HCAL" << status << std::endl;
               theHcalHitMaker->addHit(ri, phi);
               //               std::cout << " Added a spot in the HCAL" << status << std::endl;
             }
             //  if (ecal) E9 += SpotEnergy;
             //  if (presh1) { E1 += SpotEnergy; ++n1; }
             //  if (presh2) { E2 += SpotEnergy; ++n2; }
 
             Etot[i] += spote;
           }
         }
       }
       //      std::cout << " Done with the step " << std::endl;
       // The shower!
       //myHistos->fill("h500",theSpot.z(),theSpot.perp());
     }
     //    std::cout << " nPart " << nPart << std::endl;
   }
   //  std::cout << " Finshed the step loop " << std::endl;
   //  myHistos->fill("h500",E1+0.7*E2,E9);
   //  myHistos->fill("h501",n1+0.7*n2,E9);
   //  myHistos->fill("h400",n1);
   //  myHistos->fill("h401",n2);
   //  myHistos->fill("h402",E9+E1+0.7*E2);
   //  if(!standalone)theGrid->printGrid();
   /*
   double Etotal = 0.;
   for (unsigned i = 0; i < nPart; ++i) {
     //      myHistos->fill("h10",Etot[i]);
     Etotal += Etot[i];
   }
   */
 
   //  std::cout << "Etotal = " << Etotal << " nPart = "<< nPart << std::endl;
   //  std::cout << "totECalc = " << totECalc << std::endl;
 
   //  myHistos->fill("h20",Etotal);
   //  if(thePreshower)
   //    std::cout << " PS " << thePreshower->layer1Calibrated() << " " << thePreshower->layer2Calibrated() << " " << thePreshower->totalCalibrated() << " " << ps1tot << " " <<ps2tot << " " << pstot << std::endl;
 }
 
 double EMShower::gam(double x, double a) const {
   // A stupid gamma function
   return std::pow(x, a - 1.) * std::exp(-x);
 }
 
 //double
 //EMShower::deposit(double t, double a, double b, double dt) {
 //
 //  // The number of integration steps (about 1 / X0)
 //  int numberOfSteps = (int)dt+1;
 //
 //  // The size if the integration step
 //  double integrationstep = dt/(double)numberOfSteps;
 //
 //  // Half this size
 //  double halfis = 0.5*integrationstep;
 //
 //  double dE = 0.;
 //
 //  for(double tt=t-dt+halfis;tt<t;tt+=integrationstep) {
 //
 //    // Simpson integration over each of these steps
 //    dE +=   gam(b*(tt-halfis),a)
 //       + 4.*gam(b* tt        ,a)
 //       +    gam(b*(tt+halfis),a);
 //
 //  }
 //
 //  // Normalization
 //  dE *= b*integrationstep/tgamma(a)/6.;
 //
 //  // There we go.
 //  return dE;
 //}
 
 double EMShower::deposit(double t, double a, double b, double dt) {
   myIncompleteGamma.a().setValue(a);
   double b1 = b * (t - dt);
   double b2 = b * t;
   double result = 0.;
   double rb1 = (b1 != 0.) ? myIncompleteGamma(b1) : 0.;
   double rb2 = (b2 != 0.) ? myIncompleteGamma(b2) : 0.;
   result = (rb2 - rb1);
   return result;
 }
 
 void EMShower::setIntervals(unsigned icomp, RadialInterval& rad) {
   //  std::cout << " Got the pointer " << std::endl;
   const std::vector<double>& myValues((icomp) ? theParam->getTailIntervals() : theParam->getCoreIntervals());
   //  std::cout << " Got the vector " << myValues.size () << std::endl;
   unsigned nvals = myValues.size() / 2;
   for (unsigned iv = 0; iv < nvals; ++iv) {
     //      std::cout << myValues[2*iv] << " " <<  myValues[2*iv+1] <<std::endl;
     rad.addInterval(myValues[2 * iv], myValues[2 * iv + 1]);
   }
 }
 
 void EMShower::setPreshower(PreshowerHitMaker* const myPresh) {
   if (myPresh != nullptr) {
     thePreshower = myPresh;
     hasPreshower = true;
   }
 }
 
 void EMShower::setHcal(HcalHitMaker* const myHcal) { theHcalHitMaker = myHcal; }
 
 double EMShower::deposit(double a, double b, double t) {
   //  std::cout << " Deposit " << std::endl;
   myIncompleteGamma.a().setValue(a);
   double b2 = b * t;
   double result = 0.;
   if (fabs(b2) < 1.e-9)
     b2 = 1.e-9;
   result = myIncompleteGamma(b2);
   //  std::cout << " deposit t = " << t  << " "  << result <<std::endl;
   return result;
 }

This page was automatically generated by the 2.2.1 LXR engine. The LXR team