ShowerDevelopment/src/RadialInterval.cc

0001 //FAMOS header
0002 #include "FastSimulation/ShowerDevelopment/interface/RadialInterval.h"
0003 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0004
0005 #include <cmath>
0006
0007 RadialInterval::RadialInterval(double RC, unsigned nSpots, double energy, const RandomEngineAndDistribution* engine)
0008     : theR(RC), theNumberOfSpots(nSpots), theSpotEnergy(energy), random(engine) {
0009   currentRad = 0.;
0010   currentEnergyFraction = 0.;
0011   currentUlim = 0.;
0012   nInter = 0;
0013 }
0014
0015 void RadialInterval::addInterval(double radius, double spotf) {
0016   double radiussq = radius * radius;
0017   double enFrac = energyFractionInRadius(radius);
0018   if (radius > 10)
0019     enFrac = 1.;
0020   double energyFrac = enFrac - currentEnergyFraction;
0021   currentEnergyFraction = enFrac;
0022   //  std::cout << " Rad " << nInter << " Energy frac " << energyFrac << std::endl;
0023
0024   // Calculates the number of spots. Add binomial fluctuations
0025   double dspot =
0026       random->gaussShoot(theNumberOfSpots * energyFrac, sqrt(energyFrac * (1. - energyFrac) * theNumberOfSpots));
0027   //  std::cout << " Normal : " << theNumberOfSpots*energyFrac << " "<< dspot << std::endl;
0028   unsigned nspot = (unsigned)(dspot * spotf + 0.5);
0029   //  if(nspot<100)
0030   //    {
0031   //      spotf=1.;
0032   //      nspot=(unsigned)(theNumberOfSpots*energyFrac+0.5);
0033   //    }
0034
0035   dspotsunscaled.push_back(dspot);
0036   spotfraction.push_back(spotf);
0037
0038   double spotEnergy = theSpotEnergy / spotf;
0039   //  std::cout << " The number of spots " << theNumberOfSpots << " " << nspot << std::endl;
0040
0041   // This is not correct for the last interval, but will be overriden later
0042   nspots.push_back(nspot);
0043   spotE.push_back(spotEnergy);
0044   // computes the limits
0045   double umax = radiussq / (radiussq + theR * theR);
0046   if (radius > 10) {
0047     umax = 1.;
0048   }
0049
0050   // Stores the limits
0051   uMax.push_back(umax);
0052   uMin.push_back(currentUlim);
0053   currentUlim = umax;
0054
0055   // Stores the energy
0056   //  std::cout << " SpotE " << theSpotEnergy<< " " << spotf << " " << theSpotEnergy/spotf<< std::endl;
0057
0058   ++nInter;
0059 }
0060
0061 double RadialInterval::energyFractionInRadius(double rm) {
0062   double rm2 = rm * rm;
0063   return (rm2 / (rm2 + theR * theR));
0064 }
0065
0066 void RadialInterval::compute() {
0067   //  std::cout << " The number of Spots " << theNumberOfSpots << std::endl;
0068   //  std::cout << " Size : " << nInter << " " << nspots.size() << " " << dspotsunscaled.size() << std::endl;
0069   double ntotspots = 0.;
0070   for (unsigned irad = 0; irad < nInter - 1; ++irad) {
0071     ntotspots += dspotsunscaled[irad];
0072     //    std::cout << " In the loop " << ntotspots << std::endl;
0073   }
0074
0075   // This can happen (fluctuations)
0076   if (ntotspots > theNumberOfSpots)
0077     ntotspots = (double)theNumberOfSpots;
0078   //  std::cout << " Sous-total " << ntotspots << std::endl;
0079   dspotsunscaled[nInter - 1] = (double)theNumberOfSpots - ntotspots;
0080
0081   nspots[nInter - 1] = (unsigned)(dspotsunscaled[nInter - 1] * spotfraction[nInter - 1] + 0.5);
0082   //    std::cout << " Nlast " << nspots[nInter-1] << std::endl;
0083 }