ShowerDevelopment/src/EMShower.cc

0001 //FAMOS Headers
0002 #include "FastSimulation/ShowerDevelopment/interface/EMShower.h"
0003 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0004 #include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
0005 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0006
0007 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0008 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
0009
0010 #include <cmath>
0011
0012 //#include "FastSimulation/Utilities/interface/Histos.h"
0013
0014 using std::vector;
0015
0016 EMShower::EMShower(const RandomEngineAndDistribution* engine,
0017                    GammaFunctionGenerator* gamma,
0018                    EMECALShowerParametrization* const myParam,
0019                    vector<const RawParticle*>* const myPart,
0020                    EcalHitMaker* const myGrid,
0021                    PreshowerHitMaker* const myPresh,
0022                    bool bFixedLength)
0023
0024     : theParam(myParam),
0025       thePart(myPart),
0026       theGrid(myGrid),
0027       thePreshower(myPresh),
0028       random(engine),
0029       myGammaGenerator(gamma),
0030       bFixedLength_(bFixedLength) {
0031   // Get the Famos Histos pointer
0032   //  myHistos = Histos::instance();
0033   //  myGammaGenerator = GammaFunctionGenerator::instance();
0034   stepsCalculated = false;
0035   hasPreshower = myPresh != nullptr;
0036   theECAL = myParam->ecalProperties();
0037   theHCAL = myParam->hcalProperties();
0038   theLayer1 = myParam->layer1Properties();
0039   theLayer2 = myParam->layer2Properties();
0040
0041   double fotos = theECAL->photoStatistics() * theECAL->lightCollectionEfficiency();
0042
0043   nPart = thePart->size();
0044   totalEnergy = 0.;
0045   globalMaximum = 0.;
0046   //double meanDepth = 0.;
0047   // Initialize the shower parameters for each particle
0048
0049   for (unsigned int i = 0; i < nPart; ++i) {
0050     //    std::cout << " AAA " << *(*thePart)[i] << std::endl;
0051     // The particle and the shower energy
0052     Etot.push_back(0.);
0053     E.push_back(((*thePart)[i])->e());
0054     totalEnergy += E[i];
0055
0056     double lny = std::log(E[i] / theECAL->criticalEnergy());
0057
0058     // Average and Sigma for T and alpha
0059     double theMeanT = myParam->meanT(lny);
0060     double theMeanAlpha = myParam->meanAlpha(lny);
0061     double theMeanLnT = myParam->meanLnT(lny);
0062     double theMeanLnAlpha = myParam->meanLnAlpha(lny);
0063     double theSigmaLnT = myParam->sigmaLnT(lny);
0064     double theSigmaLnAlpha = myParam->sigmaLnAlpha(lny);
0065
0066     // The correlation matrix
0067     double theCorrelation = myParam->correlationAlphaT(lny);
0068     double rhop = std::sqrt((1. + theCorrelation) / 2.);
0069     double rhom = std::sqrt((1. - theCorrelation) / 2.);
0070
0071     // The number of spots in ECAL / HCAL
0072     theNumberOfSpots.push_back(myParam->nSpots(E[i]));
0073     //    theNumberOfSpots.push_back(myParam->nSpots(E[i])*spotFraction);
0074     //theNumberOfSpots = random->poissonShoot(myParam->nSpots(myPart->e()));
0075
0076     // Photo-statistics
0077     photos.push_back(E[i] * fotos);
0078
0079     // The longitudinal shower development parameters
0080     // Fluctuations of alpha, T and beta
0081     double z1 = 0.;
0082     double z2 = 0.;
0083     double aa = 0.;
0084
0085     // Protect against too large fluctuations (a < 1) for small energies
0086     while (aa <= 1.) {
0087       z1 = random->gaussShoot(0., 1.);
0088       z2 = random->gaussShoot(0., 1.);
0089       aa = std::exp(theMeanLnAlpha + theSigmaLnAlpha * (z1 * rhop - z2 * rhom));
0090     }
0091
0092     a.push_back(aa);
0093     T.push_back(std::exp(theMeanLnT + theSigmaLnT * (z1 * rhop + z2 * rhom)));
0094     b.push_back((a[i] - 1.) / T[i]);
0095     maximumOfShower.push_back((a[i] - 1.) / b[i]);
0096     globalMaximum += maximumOfShower[i] * E[i];
0097     //meanDepth += a[i] / b[i] * E[i];
0098     //    std::cout << " Adding max " << maximumOfShower[i] << " " << E[i] << " " <<maximumOfShower[i]*E[i] << std::endl;
0099     //    std::cout << std::setw(8) << std::setprecision(5) << " a /b " << a[i] << " " << b[i] << std::endl;
0100     Ti.push_back(a[i] / b[i] * (std::exp(theMeanLnAlpha) - 1.) / std::exp(theMeanLnAlpha));
0101
0102     // The parameters for the number of energy spots
0103     TSpot.push_back(theParam->meanTSpot(theMeanT));
0104     aSpot.push_back(theParam->meanAlphaSpot(theMeanAlpha));
0105     bSpot.push_back((aSpot[i] - 1.) / TSpot[i]);
0106     //    myHistos->fill("h7000",a[i]);
0107     //    myHistos->fill("h7002",E[i],a[i]);
0108   }
0109   //  std::cout << " PS1 : " << myGrid->ps1TotalX0()
0110   //         << " PS2 : " << myGrid->ps2TotalX0()
0111   //         << " ECAL : " << myGrid->ecalTotalX0()
0112   //         << " HCAL : " << myGrid->hcalTotalX0()
0113   //         << " Offset : " << myGrid->x0DepthOffset()
0114   //         << std::endl;
0115
0116   globalMaximum /= totalEnergy;
0117   //meanDepth /= totalEnergy;
0118   // std::cout << " Total Energy " << totalEnergy << " Global max " << globalMaximum << std::endl;
0119 }
0120
0121 void EMShower::prepareSteps() {
0122   //  TimeMe theT("EMShower::compute");
0123
0124   // Determine the longitudinal intervals
0125   //  std::cout << " EMShower compute" << std::endl;
0126   double dt;
0127   double radlen;
0128   int stps;
0129   int first_Ecal_step = 0;
0130   int last_Ecal_step = 0;
0131
0132   // The maximum is in principe 8 (with 5X0 steps in the ECAL)
0133   steps.reserve(24);
0134
0135   /*
0136   std::cout << " PS1 : " << theGrid->ps1TotalX0()
0137         << " PS2 : " << theGrid->ps2TotalX0()
0138         << " PS2 and ECAL : " << theGrid->ps2eeTotalX0()
0139         << " ECAL : " << theGrid->ecalTotalX0()
0140         << " HCAL : " << theGrid->hcalTotalX0()
0141         << " Offset : " << theGrid->x0DepthOffset()
0142         << std::endl;
0143   */
0144
0145   radlen = -theGrid->x0DepthOffset();
0146
0147   // Preshower Layer 1
0148   radlen += theGrid->ps1TotalX0();
0149   if (radlen > 0.) {
0150     steps.push_back(Step(0, radlen));
0151     radlen = 0.;
0152   }
0153
0154   // Preshower Layer 2
0155   radlen += theGrid->ps2TotalX0();
0156   if (radlen > 0.) {
0157     steps.push_back(Step(1, radlen));
0158     radlen = 0.;
0159   }
0160
0161   // add a step between preshower and ee
0162   radlen += theGrid->ps2eeTotalX0();
0163   if (radlen > 0.) {
0164     steps.push_back(Step(5, radlen));
0165     radlen = 0.;
0166   }
0167
0168   // ECAL
0169   radlen += theGrid->ecalTotalX0();
0170
0171   //  std::cout << "theGrid->ecalTotalX0() = " << theGrid->ecalTotalX0() << std::endl;
0172
0173   if (radlen > 0.) {
0174     if (!bFixedLength_) {
0175       stps = (int)((radlen + 2.5) / 5.);
0176       //    stps=(int)((radlen+.5)/1.);
0177       if (stps == 0)
0178         stps = 1;
0179       dt = radlen / (double)stps;
0180       Step step(2, dt);
0181       first_Ecal_step = steps.size();
0182       for (int ist = 0; ist < stps; ++ist)
0183         steps.push_back(step);
0184       last_Ecal_step = steps.size() - 1;
0185       radlen = 0.;
0186     } else {
0187       dt = 1.0;
0188       stps = static_cast<int>(radlen);
0189       if (stps == 0)
0190         stps = 1;
0191       Step step(2, dt);
0192       first_Ecal_step = steps.size();
0193       for (int ist = 0; ist < stps; ++ist)
0194         steps.push_back(step);
0195       dt = radlen - stps;
0196       if (dt > 0) {
0197         Step stepLast(2, dt);
0198         steps.push_back(stepLast);
0199       }
0200       last_Ecal_step = steps.size() - 1;
0201       //      std::cout << "radlen = "  << radlen << " stps = " << stps << " dt = " << dt << std::endl;
0202       radlen = 0.;
0203     }
0204   }
0205
0206   // I should had a gap here !
0207
0208   // HCAL
0209   radlen += theGrid->hcalTotalX0();
0210   if (radlen > 0.) {
0211     double dtFrontHcal = theGrid->totalX0() - theGrid->hcalTotalX0();
0212     // One single step for the full HCAL
0213     if (dtFrontHcal < 30.) {
0214       dt = 30. - dtFrontHcal;
0215       Step step(3, dt);
0216       steps.push_back(step);
0217     }
0218   }
0219
0220   nSteps = steps.size();
0221   if (nSteps == 0)
0222     return;
0223   double ESliceTot = 0.;
0224   double MeanDepth = 0.;
0225   depositedEnergy.resize(nSteps);
0226   meanDepth.resize(nSteps);
0227   double t = 0.;
0228
0229   int offset = 0;
0230   for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
0231     ESliceTot = 0.;
0232     MeanDepth = 0.;
0233     double realTotalEnergy = 0;
0234     dt = steps[iStep].second;
0235     t += dt;
0236     for (unsigned int i = 0; i < nPart; ++i) {
0237       depositedEnergy[iStep].push_back(deposit(t, a[i], b[i], dt));
0238       ESliceTot += depositedEnergy[iStep][i];
0239       MeanDepth += deposit(t, a[i] + 1., b[i], dt) / b[i] * a[i];
0240       realTotalEnergy += depositedEnergy[iStep][i] * E[i];
0241     }
0242
0243     if (ESliceTot > 0.)  // can happen for the shower tails; this depth will be skipped anyway
0244       MeanDepth /= ESliceTot;
0245     else
0246       MeanDepth = t - dt;
0247
0248     meanDepth[iStep] = MeanDepth;
0249     if (realTotalEnergy < 0.001) {
0250       offset -= 1;
0251     }
0252   }
0253
0254   innerDepth = meanDepth[first_Ecal_step];
0255   if (last_Ecal_step + offset >= 0)
0256     outerDepth = meanDepth[last_Ecal_step + offset];
0257   else
0258     outerDepth = innerDepth;
0259
0260   stepsCalculated = true;
0261 }
0262
0263 void EMShower::compute() {
0264   double t = 0.;
0265   double dt = 0.;
0266   if (!stepsCalculated)
0267     prepareSteps();
0268
0269   // Prepare the grids in EcalHitMaker
0270   // theGrid->setInnerAndOuterDepth(innerDepth,outerDepth);
0271   //float pstot = 0.;
0272   //float ps2tot = 0.;
0273   //float ps1tot = 0.;
0274   bool status = false;
0275   //  double E1 = 0.;  // Energy layer 1
0276   //  double E2 = 0.;  // Energy layer 2
0277   //  double n1 = 0.;  // #mips layer 1
0278   //  double n2 = 0.;  // #mips layer 2
0279   //  double E9 = 0.;  // Energy ECAL
0280
0281   // Loop over all segments for the longitudinal development
0282   //double totECalc = 0;
0283
0284   for (unsigned iStep = 0; iStep < nSteps; ++iStep) {
0285     // The length of the shower in this segment
0286     dt = steps[iStep].second;
0287
0288     //    std::cout << " Detector " << steps[iStep].first << " t " << t << " " << dt << std::endl;
0289
0290     // The elapsed length
0291     t += dt;
0292
0293     // In what detector are we ?
0294     unsigned detector = steps[iStep].first;
0295
0296     bool presh1 = detector == 0;
0297     bool presh2 = detector == 1;
0298     bool ecal = detector == 2;
0299     bool hcal = detector == 3;
0300     bool vfcal = detector == 4;
0301     bool gap = detector == 5;
0302
0303     // Temporary. Will be removed
0304     if (theHCAL == nullptr)
0305       hcal = false;
0306
0307     // Keep only ECAL for now
0308     if (vfcal)
0309       continue;
0310
0311     // Nothing to do in the gap
0312     if (gap)
0313       continue;
0314
0315     //    cout << " t = " << t << endl;
0316     // Build the grid of crystals at this ECAL depth
0317     // Actually, it might be useful to check if this grid is empty or not.
0318     // If it is empty (because no crystal at this depth), it is of no use
0319     // (and time consuming) to generate the spots
0320
0321     // middle of the step
0322     double tt = t - 0.5 * dt;
0323
0324     double realTotalEnergy = 0.;
0325     for (unsigned int i = 0; i < nPart; ++i) {
0326       realTotalEnergy += depositedEnergy[iStep][i] * E[i];
0327     }
0328
0329     //    std::cout << " Step " << tt << std::endl;
0330     //    std::cout << "ecal " << ecal << " hcal "  << hcal <<std::endl;
0331
0332     // If the amount of energy is greater than 1 MeV, make a new grid
0333     // otherwise put in the previous one.
0334     bool usePreviousGrid = (realTotalEnergy < 0.001);
0335
0336     // If the amount of energy is greater than 1 MeV, make a new grid
0337     // otherwise put in the previous one.
0338
0339     // If less than 1 kEV. Just skip
0340     if (iStep > 2 && realTotalEnergy < 0.000001)
0341       continue;
0342
0343     if (ecal && !usePreviousGrid) {
0344       status = theGrid->getPads(meanDepth[iStep]);
0345     }
0346     if (hcal) {
0347       status = theHcalHitMaker->setDepth(tt);
0348     }
0349     if ((ecal || hcal) && !status)
0350       continue;
0351
0352     bool detailedShowerTail = false;
0353     // check if a detailed treatment of the rear leakage should be applied
0354     if (ecal && !usePreviousGrid) {
0355       detailedShowerTail = (t - dt > theGrid->getX0back());
0356     }
0357
0358     // The particles of the shower are processed in parallel
0359     for (unsigned int i = 0; i < nPart; ++i) {
0360       //      double Edepo=deposit(t,a[i],b[i],dt);
0361
0362       //  integration of the shower profile between t-dt and t
0363       double dE = (!hcal) ? depositedEnergy[iStep][i] : 1. - deposit(a[i], b[i], t - dt);
0364
0365       // no need to do the full machinery if there is ~nothing to distribute)
0366       if (dE * E[i] < 0.000001)
0367         continue;
0368
0369       if (ecal && !theECAL->isHom()) {
0370         double mean = dE * E[i];
0371         double sigma = theECAL->resE() * sqrt(mean);
0372
0373         /*
0374       double meanLn = log(mean);
0375       double kLn = sigma/mean+1;
0376       double sigmaLn = log(kLn);
0377     */
0378
0379         double dE0 = dE;
0380
0381         //    std::cout << "dE before shoot = " << dE << std::endl;
0382         dE = random->gaussShoot(mean, sigma) / E[i];
0383
0384         //    myGammaGenerator->setParameters(aSam,bSam,0);
0385         //    dE = myGammaGenerator->shoot()/E[i];
0386         //    std::cout << "dE shooted = " << dE << " E[i] = " << E[i] << std::endl;
0387         if (dE * E[i] < 0.000001)
0388           continue;
0389         photos[i] = photos[i] * dE / dE0;
0390       }
0391
0392       /*
0393       if (ecal && !theParam->ecalProperties()->isHom()){
0394
0395     double cSquare = TMath::Power(theParam->ecalProperties()->resE(),2);
0396     double aSam = dE/cSquare;
0397     double bSam = 1./cSquare;
0398
0399     //  dE = dE*gam(bSam*dE, aSam)/tgamma(aSam);
0400       }
0401       */
0402
0403       //totECalc += dE;
0404
0405       // The number of energy spots (or mips)
0406       double nS = 0;
0407
0408       // ECAL case : Account for photostatistics and long'al non-uniformity
0409       if (ecal) {
0410         //  double aSam = E[i]*dE*one_over_resoSquare;
0411         //  double bSam = one_over_resoSquare;
0412
0413         dE = random->poissonShoot(dE * photos[i]) / photos[i];
0414         double z0 = random->gaussShoot(0., 1.);
0415         dE *= 1. + z0 * theECAL->lightCollectionUniformity();
0416
0417         // Expected spot number
0418         nS = (theNumberOfSpots[i] * gam(bSpot[i] * tt, aSpot[i]) * bSpot[i] * dt / tgamma(aSpot[i]));
0419
0420         // Preshower : Expected number of mips + fluctuation
0421       } else if (hcal) {
0422         nS = (theNumberOfSpots[i] * gam(bSpot[i] * tt, aSpot[i]) * bSpot[i] * dt / tgamma(aSpot[i])) *
0423              theHCAL->spotFraction();
0424         double nSo = nS;
0425         nS = random->poissonShoot(nS);
0426         // 'Quick and dirty' fix (but this line should be better removed):
0427         if (nSo > 0. && nS / nSo < 10.)
0428           dE *= nS / nSo;
0429
0430         //  if(true)
0431         //    {
0432         //      std::cout << " theHCAL->spotFraction = " <<theHCAL->spotFraction() <<std::endl;
0433         //      std::cout << " nSpot Ecal : " << nSo/theHCAL->spotFraction() << " Final " << nS << std::endl;
0434         //    }
0435       } else if (presh1) {
0436         nS = random->poissonShoot(dE * E[i] * theLayer1->mipsPerGeV());
0437         //  std::cout << " dE *E[i] (1)" << dE*E[i] << " " << dE*E[i]*theLayer1->mipsPerGeV() << "  "<< nS << std::endl;
0438         //pstot += dE * E[i];
0439         //ps1tot += dE * E[i];
0440         dE = nS / (E[i] * theLayer1->mipsPerGeV());
0441
0442         //        E1 += dE*E[i];
0443         //  n1 += nS;
0444         //  if (presh2) { E2 += SpotEnergy; ++n2; }
0445
0446       } else if (presh2) {
0447         nS = random->poissonShoot(dE * E[i] * theLayer2->mipsPerGeV());
0448         //  std::cout << " dE *E[i] (2) " << dE*E[i] << " " << dE*E[i]*theLayer2->mipsPerGeV() << "  "<< nS << std::endl;
0449         //pstot += dE * E[i];
0450         //ps2tot += dE * E[i];
0451         dE = nS / (E[i] * theLayer2->mipsPerGeV());
0452
0453         //        E2 += dE*E[i];
0454         //  n2 += nS;
0455       }
0456
0457       if (detailedShowerTail)
0458         myGammaGenerator->setParameters(floor(a[i] + 0.5), b[i], t - dt);
0459
0460       //    myHistos->fill("h100",t,dE);
0461
0462       // The lateral development parameters
0463
0464       // Energy of the spots
0465       double eSpot = (nS > 0.) ? dE / nS : 0.;
0466       double SpotEnergy = eSpot * E[i];
0467
0468       if (hasPreshower && (presh1 || presh2))
0469         thePreshower->setSpotEnergy(0.00009);
0470       if (hcal) {
0471         SpotEnergy *= theHCAL->hOverPi();
0472         theHcalHitMaker->setSpotEnergy(SpotEnergy);
0473       }
0474       // Poissonian fluctuations for the number of spots
0475       //    int nSpot = random->poissonShoot(nS);
0476       int nSpot = (int)(nS + 0.5);
0477
0478       // Fig. 11 (right) *** Does not match.
0479       //    myHistos->fill("h101",t,(double)nSpot/theNumberOfSpots);
0480
0481       //double taui = t/T;
0482       double taui = tt / Ti[i];
0483       double proba = theParam->p(taui, E[i]);
0484       double theRC = theParam->rC(taui, E[i]);
0485       double theRT = theParam->rT(taui, E[i]);
0486
0487       // Fig. 10
0488       //    myHistos->fill("h300",taui,theRC);
0489       //    myHistos->fill("h301",taui,theRT);
0490       //    myHistos->fill("h302",taui,proba);
0491
0492       double dSpotsCore = random->gaussShoot(proba * nSpot, std::sqrt(proba * (1. - proba) * nSpot));
0493
0494       if (dSpotsCore < 0)
0495         dSpotsCore = 0;
0496
0497       unsigned nSpots_core = (unsigned)(dSpotsCore + 0.5);
0498       unsigned nSpots_tail = ((unsigned)nSpot > nSpots_core) ? nSpot - nSpots_core : 0;
0499
0500       for (unsigned icomp = 0; icomp < 2; ++icomp) {
0501         double theR = (icomp == 0) ? theRC : theRT;
0502         unsigned ncompspots = (icomp == 0) ? nSpots_core : nSpots_tail;
0503
0504         RadialInterval radInterval(theR, ncompspots, SpotEnergy, random);
0505         if (ecal) {
0506           if (icomp == 0) {
0507             setIntervals(icomp, radInterval);
0508           } else {
0509             setIntervals(icomp, radInterval);
0510           }
0511         } else {
0512           radInterval.addInterval(100., 1.);  // 100% of the spots
0513         }
0514
0515         radInterval.compute();
0516         // irad = 0 : central circle; irad=1 : outside
0517
0518         unsigned nrad = radInterval.nIntervals();
0519
0520         for (unsigned irad = 0; irad < nrad; ++irad) {
0521           double spote = radInterval.getSpotEnergy(irad);
0522           if (ecal)
0523             theGrid->setSpotEnergy(spote);
0524           if (hcal)
0525             theHcalHitMaker->setSpotEnergy(spote);
0526           unsigned nradspots = radInterval.getNumberOfSpots(irad);
0527           double umin = radInterval.getUmin(irad);
0528           double umax = radInterval.getUmax(irad);
0529           // Go for the lateral development
0530           //           std::cout << "Couche " << iStep << " irad = " << irad << " Ene = " << E[i] << " eSpot = " << eSpot << " spote = " << spote << " nSpot = " << nS << std::endl;
0531
0532           for (unsigned ispot = 0; ispot < nradspots; ++ispot) {
0533             double z3 = random->flatShoot(umin, umax);
0534             double ri = theR * std::sqrt(z3 / (1. - z3));
0535
0536             // Generate phi
0537             double phi = 2. * M_PI * random->flatShoot();
0538
0539             // Add the hit in the crystal
0540             //  if( ecal ) theGrid->addHit(ri*theECAL->moliereRadius(),phi);
0541             // Now the *moliereRadius is done in EcalHitMaker
0542             if (ecal) {
0543               if (detailedShowerTail) {
0544                 //             std::cout << "About to call addHitDepth " << std::endl;
0545                 double depth;
0546                 do {
0547                   depth = myGammaGenerator->shoot(random);
0548                 } while (depth > t);
0549                 theGrid->addHitDepth(ri, phi, depth);
0550                 //             std::cout << " Done " << std::endl;
0551               } else
0552                 theGrid->addHit(ri, phi);
0553             } else if (hasPreshower && presh1)
0554               thePreshower->addHit(ri, phi, 1);
0555             else if (hasPreshower && presh2)
0556               thePreshower->addHit(ri, phi, 2);
0557             else if (hcal) {
0558               //               std::cout << " About to add a spot in the HCAL" << status << std::endl;
0559               theHcalHitMaker->addHit(ri, phi);
0560               //               std::cout << " Added a spot in the HCAL" << status << std::endl;
0561             }
0562             //  if (ecal) E9 += SpotEnergy;
0563             //  if (presh1) { E1 += SpotEnergy; ++n1; }
0564             //  if (presh2) { E2 += SpotEnergy; ++n2; }
0565
0566             Etot[i] += spote;
0567           }
0568         }
0569       }
0570       //      std::cout << " Done with the step " << std::endl;
0571       // The shower!
0572       //myHistos->fill("h500",theSpot.z(),theSpot.perp());
0573     }
0574     //    std::cout << " nPart " << nPart << std::endl;
0575   }
0576   //  std::cout << " Finshed the step loop " << std::endl;
0577   //  myHistos->fill("h500",E1+0.7*E2,E9);
0578   //  myHistos->fill("h501",n1+0.7*n2,E9);
0579   //  myHistos->fill("h400",n1);
0580   //  myHistos->fill("h401",n2);
0581   //  myHistos->fill("h402",E9+E1+0.7*E2);
0582   //  if(!standalone)theGrid->printGrid();
0583   /*
0584   double Etotal = 0.;
0585   for (unsigned i = 0; i < nPart; ++i) {
0586     //      myHistos->fill("h10",Etot[i]);
0587     Etotal += Etot[i];
0588   }
0589   */
0590
0591   //  std::cout << "Etotal = " << Etotal << " nPart = "<< nPart << std::endl;
0592   //  std::cout << "totECalc = " << totECalc << std::endl;
0593
0594   //  myHistos->fill("h20",Etotal);
0595   //  if(thePreshower)
0596   //    std::cout << " PS " << thePreshower->layer1Calibrated() << " " << thePreshower->layer2Calibrated() << " " << thePreshower->totalCalibrated() << " " << ps1tot << " " <<ps2tot << " " << pstot << std::endl;
0597 }
0598
0599 double EMShower::gam(double x, double a) const {
0600   // A stupid gamma function
0601   return std::pow(x, a - 1.) * std::exp(-x);
0602 }
0603
0604 //double
0605 //EMShower::deposit(double t, double a, double b, double dt) {
0606 //
0607 //  // The number of integration steps (about 1 / X0)
0608 //  int numberOfSteps = (int)dt+1;
0609 //
0610 //  // The size if the integration step
0611 //  double integrationstep = dt/(double)numberOfSteps;
0612 //
0613 //  // Half this size
0614 //  double halfis = 0.5*integrationstep;
0615 //
0616 //  double dE = 0.;
0617 //
0618 //  for(double tt=t-dt+halfis;tt<t;tt+=integrationstep) {
0619 //
0620 //    // Simpson integration over each of these steps
0621 //    dE +=   gam(b*(tt-halfis),a)
0622 //       + 4.*gam(b* tt        ,a)
0623 //       +    gam(b*(tt+halfis),a);
0624 //
0625 //  }
0626 //
0627 //  // Normalization
0628 //  dE *= b*integrationstep/tgamma(a)/6.;
0629 //
0630 //  // There we go.
0631 //  return dE;
0632 //}
0633
0634 double EMShower::deposit(double t, double a, double b, double dt) {
0635   myIncompleteGamma.a().setValue(a);
0636   double b1 = b * (t - dt);
0637   double b2 = b * t;
0638   double result = 0.;
0639   double rb1 = (b1 != 0.) ? myIncompleteGamma(b1) : 0.;
0640   double rb2 = (b2 != 0.) ? myIncompleteGamma(b2) : 0.;
0641   result = (rb2 - rb1);
0642   return result;
0643 }
0644
0645 void EMShower::setIntervals(unsigned icomp, RadialInterval& rad) {
0646   //  std::cout << " Got the pointer " << std::endl;
0647   const std::vector<double>& myValues((icomp) ? theParam->getTailIntervals() : theParam->getCoreIntervals());
0648   //  std::cout << " Got the vector " << myValues.size () << std::endl;
0649   unsigned nvals = myValues.size() / 2;
0650   for (unsigned iv = 0; iv < nvals; ++iv) {
0651     //      std::cout << myValues[2*iv] << " " <<  myValues[2*iv+1] <<std::endl;
0652     rad.addInterval(myValues[2 * iv], myValues[2 * iv + 1]);
0653   }
0654 }
0655
0656 void EMShower::setPreshower(PreshowerHitMaker* const myPresh) {
0657   if (myPresh != nullptr) {
0658     thePreshower = myPresh;
0659     hasPreshower = true;
0660   }
0661 }
0662
0663 void EMShower::setHcal(HcalHitMaker* const myHcal) { theHcalHitMaker = myHcal; }
0664
0665 double EMShower::deposit(double a, double b, double t) {
0666   //  std::cout << " Deposit " << std::endl;
0667   myIncompleteGamma.a().setValue(a);
0668   double b2 = b * t;
0669   double result = 0.;
0670   if (fabs(b2) < 1.e-9)
0671     b2 = 1.e-9;
0672   result = myIncompleteGamma(b2);
0673   //  std::cout << " deposit t = " << t  << " "  << result <<std::endl;
0674   return result;
0675 }