EMShower Class Reference

#include <EMShower.h>

Public Member Functions
void	compute ()
	Compute the shower longitudinal and lateral development. More...
	EMShower (RandomEngineAndDistribution const *engine, GammaFunctionGenerator *gamma, EMECALShowerParametrization *const myParam, std::vector< const RawParticle *> *const myPart, EcalHitMaker *const myGrid=nullptr, PreshowerHitMaker *const myPreshower=nullptr, bool bFixedLength=false)
double	getMaximumOfShower () const
	get the depth of the centre of gravity of the shower(s) More...
void	prepareSteps ()
	Computes the steps before the real compute. More...
void	setGrid (EcalHitMaker *const myGrid)
	set the grid address More...
void	setHcal (HcalHitMaker *const myHcal)
	set the HCAL address More...
void	setPreshower (PreshowerHitMaker *const myPresh)
	set the preshower address More...
virtual	~EMShower ()

Private Types
typedef std::pair< XYZPoint, double >	Spot
typedef std::pair< unsigned int, double >	Step
typedef Steps::const_iterator	step_iterator
typedef std::vector< Step >	Steps
typedef math::XYZVector	XYZPoint

Private Member Functions
double	deposit (double t, double a, double b, double dt)
double	deposit (double a, double b, double t)
double	gam (double x, double a) const
void	setIntervals (unsigned icomp, RadialInterval &rad)

Private Attributes
std::vector< double >	a
std::vector< double >	aSpot
std::vector< double >	b
bool	bFixedLength_
std::vector< double >	bSpot
std::vector< std::vector< double > >	depositedEnergy
std::vector< double >	E
std::vector< double >	Etot
double	globalMaximum
bool	hasPreshower
double	innerDepth
std::vector< double >	maximumOfShower
std::vector< double >	meanDepth
GammaFunctionGenerator *	myGammaGenerator
Genfun::IncompleteGamma	myIncompleteGamma
unsigned int	nPart
unsigned	nSteps
double	outerDepth
std::vector< double >	photos
const RandomEngineAndDistribution *	random
Steps	steps
bool	stepsCalculated
std::vector< double >	T
const ECALProperties *	theECAL
EcalHitMaker *	theGrid
const HCALProperties *	theHCAL
HcalHitMaker *	theHcalHitMaker
const PreshowerLayer1Properties *	theLayer1
const PreshowerLayer2Properties *	theLayer2
std::vector< double >	theNumberOfSpots
EMECALShowerParametrization *const	theParam
std::vector< const RawParticle * > *const	thePart
PreshowerHitMaker *	thePreshower
std::vector< double >	Ti
double	totalEnergy
std::vector< double >	TSpot

Detailed Description

Definition at line 26 of file EMShower.h.

Member Typedef Documentation

Spot

typedef std::pair<XYZPoint, double> EMShower::Spot

private

Definition at line 29 of file EMShower.h.

Step

typedef std::pair<unsigned int, double> EMShower::Step

private

Definition at line 30 of file EMShower.h.

step_iterator

typedef Steps::const_iterator EMShower::step_iterator

private

Definition at line 32 of file EMShower.h.

Steps

typedef std::vector<Step> EMShower::Steps

private

Definition at line 31 of file EMShower.h.

XYZPoint

typedef math::XYZVector EMShower::XYZPoint

private

Definition at line 27 of file EMShower.h.

Constructor & Destructor Documentation

EMShower()

EMShower::EMShower	(	RandomEngineAndDistribution const *	engine,
		GammaFunctionGenerator *	gamma,
		EMECALShowerParametrization *const	myParam,
		std::vector< const RawParticle *> *const	myPart,
		EcalHitMaker *const	myGrid = nullptr,
		PreshowerHitMaker *const	myPreshower = nullptr,
		bool	bFixedLength = false
	)

Definition at line 16 of file EMShower.cc.

References a, aSpot, b, bSpot, EMECALShowerParametrization::correlationAlphaT(), ECALProperties::criticalEnergy(), MillePedeFileConverter_cfg::e, E, EMECALShowerParametrization::ecalProperties(), Etot, JetChargeProducer_cfi::exp, RandomEngineAndDistribution::gaussShoot(), globalMaximum, hasPreshower, EMECALShowerParametrization::hcalProperties(), mps_fire::i, EMECALShowerParametrization::layer1Properties(), EMECALShowerParametrization::layer2Properties(), ECALProperties::lightCollectionEfficiency(), dqm-mbProfile::log, maximumOfShower, EMECALShowerParametrization::meanAlpha(), EMECALShowerParametrization::meanAlphaSpot(), EMECALShowerParametrization::meanLnAlpha(), EMECALShowerParametrization::meanLnT(), EMECALShowerParametrization::meanT(), EMECALShowerParametrization::meanTSpot(), nPart, EMECALShowerParametrization::nSpots(), photos, ECALProperties::photoStatistics(), random, jetAnalyzer_cfi::rhom, EMECALShowerParametrization::sigmaLnAlpha(), EMECALShowerParametrization::sigmaLnT(), mathSSE::sqrt(), stepsCalculated, theECAL, theHCAL, theLayer1, theLayer2, theNumberOfSpots, theParam, thePart, Ti, totalEnergy, TSpot, and testProducerWithPsetDescEmpty_cfi::z2.

    : 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;
}

~EMShower()

virtual EMShower::~EMShower ( )

inlinevirtual

Definition at line 43 of file EMShower.h.

{ ; }

Member Function Documentation

compute()

void EMShower::compute

(

)

Compute the shower longitudinal and lateral development.

Definition at line 263 of file EMShower.cc.

References a, PreshowerHitMaker::addHit(), HcalHitMaker::addHit(), EcalHitMaker::addHit(), EcalHitMaker::addHitDepth(), RadialInterval::addInterval(), aSpot, b, bSpot, RadialInterval::compute(), deposit(), depositedEnergy, LEDCalibrationChannels::depth, hgcalTestNeighbor_cfi::detector, dt, E, Etot, RandomEngineAndDistribution::flatShoot(), gam(), HLT_2023v11_cff::gap, RandomEngineAndDistribution::gaussShoot(), RadialInterval::getNumberOfSpots(), EcalHitMaker::getPads(), RadialInterval::getSpotEnergy(), RadialInterval::getUmax(), RadialInterval::getUmin(), EcalHitMaker::getX0back(), hasPreshower, hltEgammaHLTExtra_cfi::hcal, HCALProperties::hOverPi(), mps_fire::i, createfilelist::int, ECALProperties::isHom(), ECALProperties::lightCollectionUniformity(), M_PI, SiStripPI::mean, meanDepth, PreshowerLayer1Properties::mipsPerGeV(), PreshowerLayer2Properties::mipsPerGeV(), myGammaGenerator, RadialInterval::nIntervals(), nPart, nSteps, EMECALShowerParametrization::p(), phi, photos, RandomEngineAndDistribution::poissonShoot(), prepareSteps(), random, EMECALShowerParametrization::rC(), ECALProperties::resE(), EMECALShowerParametrization::rT(), HcalHitMaker::setDepth(), setIntervals(), GammaFunctionGenerator::setParameters(), HcalHitMaker::setSpotEnergy(), PreshowerHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), GammaFunctionGenerator::shoot(), HCALProperties::spotFraction(), mathSSE::sqrt(), mps_update::status, steps, stepsCalculated, submitPVValidationJobs::t, theECAL, theGrid, theHCAL, theHcalHitMaker, theLayer1, theLayer2, theNumberOfSpots, theParam, thePreshower, Ti, and HLTMuonOfflineAnalyzer_cfi::z0.

Referenced by CalorimetryManager::EMShowerSimulation().

                       {
  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;
}

deposit() [1/2]

double EMShower::deposit ( double  t, double  a, double  b, double  dt  )

private

Definition at line 634 of file EMShower.cc.

References a, b, b1, b2, dt, myIncompleteGamma, mps_fire::result, and submitPVValidationJobs::t.

Referenced by compute(), and prepareSteps().

                                                                {
  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;
}

deposit() [2/2]

double EMShower::deposit ( double  a, double  b, double  t  )

private

Definition at line 665 of file EMShower.cc.

References a, b, b2, MillePedeFileConverter_cfg::e, myIncompleteGamma, mps_fire::result, and submitPVValidationJobs::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;
}

gam()

double EMShower::gam ( double  x, double  a  ) const

private

Definition at line 599 of file EMShower.cc.

References a, JetChargeProducer_cfi::exp, funct::pow(), and x.

Referenced by compute().

                                             {
  // A stupid gamma function
  return std::pow(x, a - 1.) * std::exp(-x);
}

getMaximumOfShower()

double EMShower::getMaximumOfShower ( ) const

inline

get the depth of the centre of gravity of the shower(s)

get the depth of the maximum of the shower

Definition at line 55 of file EMShower.h.

References globalMaximum.

Referenced by CalorimetryManager::EMShowerSimulation().

{ return globalMaximum; }

prepareSteps()

void EMShower::prepareSteps

(

)

Computes the steps before the real compute.

Definition at line 121 of file EMShower.cc.

References a, b, bFixedLength_, deposit(), depositedEnergy, dt, E, EcalHitMaker::ecalTotalX0(), EcalHitMaker::hcalTotalX0(), mps_fire::i, innerDepth, createfilelist::int, meanDepth, nPart, nSteps, hltrates_dqm_sourceclient-live_cfg::offset, outerDepth, EcalHitMaker::ps1TotalX0(), EcalHitMaker::ps2eeTotalX0(), EcalHitMaker::ps2TotalX0(), steps, stepsCalculated, submitPVValidationJobs::t, theGrid, EcalHitMaker::totalX0(), and EcalHitMaker::x0DepthOffset().

Referenced by compute().

                            {
  //  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;
}

setGrid()

void EMShower::setGrid ( EcalHitMaker *const  myGrid )

inline

set the grid address

Definition at line 58 of file EMShower.h.

References theGrid.

Referenced by CalorimetryManager::EMShowerSimulation().

{ theGrid = myGrid; }

setHcal()

void EMShower::setHcal

(

HcalHitMaker *const

myHcal

)

set the HCAL address

Definition at line 663 of file EMShower.cc.

References theHcalHitMaker.

Referenced by CalorimetryManager::EMShowerSimulation().

{ theHcalHitMaker = myHcal; }

setIntervals()

void EMShower::setIntervals ( unsigned  icomp, RadialInterval &  rad  )

private

Definition at line 645 of file EMShower.cc.

References RadialInterval::addInterval(), EMECALShowerParametrization::getCoreIntervals(), EMECALShowerParametrization::getTailIntervals(), gpuVertexFinder::iv, and theParam.

Referenced by compute().

                                                               {
  //  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]);
  }
}

setPreshower()

void EMShower::setPreshower

(

PreshowerHitMaker *const

myPresh

)

set the preshower address

Definition at line 656 of file EMShower.cc.

References hasPreshower, and thePreshower.

Referenced by CalorimetryManager::EMShowerSimulation().

                                                            {
  if (myPresh != nullptr) {
    thePreshower = myPresh;
    hasPreshower = true;
  }
}

Member Data Documentation

a

std::vector<double> EMShower::a

private

Definition at line 98 of file EMShower.h.

Referenced by compute(), deposit(), EMShower(), gam(), and prepareSteps().

aSpot

std::vector<double> EMShower::aSpot

private

Definition at line 102 of file EMShower.h.

Referenced by compute(), and EMShower().

b

std::vector<double> EMShower::b

private

Definition at line 99 of file EMShower.h.

Referenced by compute(), deposit(), EMShower(), and prepareSteps().

bFixedLength_

bool EMShower::bFixedLength_

private

Definition at line 144 of file EMShower.h.

Referenced by prepareSteps().

bSpot

std::vector<double> EMShower::bSpot

private

Definition at line 103 of file EMShower.h.

Referenced by compute(), and EMShower().

depositedEnergy

std::vector<std::vector<double> > EMShower::depositedEnergy

private

Definition at line 109 of file EMShower.h.

Referenced by compute(), and prepareSteps().

E

std::vector<double> EMShower::E

private

Definition at line 95 of file EMShower.h.

Referenced by compute(), EMShower(), and prepareSteps().

Etot

std::vector<double> EMShower::Etot

private

Definition at line 94 of file EMShower.h.

Referenced by compute(), and EMShower().

globalMaximum

double EMShower::globalMaximum

private

Definition at line 113 of file EMShower.h.

Referenced by EMShower(), and getMaximumOfShower().

hasPreshower

bool EMShower::hasPreshower

private

Definition at line 132 of file EMShower.h.

Referenced by compute(), EMShower(), and setPreshower().

innerDepth

double EMShower::innerDepth

private

Definition at line 111 of file EMShower.h.

Referenced by prepareSteps().

maximumOfShower

std::vector<double> EMShower::maximumOfShower

private

Definition at line 108 of file EMShower.h.

Referenced by EMShower().

meanDepth

std::vector<double> EMShower::meanDepth

private

Definition at line 110 of file EMShower.h.

Referenced by compute(), and prepareSteps().

myGammaGenerator

GammaFunctionGenerator* EMShower::myGammaGenerator

private

Definition at line 142 of file EMShower.h.

Referenced by compute().

myIncompleteGamma

Genfun::IncompleteGamma EMShower::myIncompleteGamma

private

Definition at line 136 of file EMShower.h.

Referenced by deposit().

nPart

unsigned int EMShower::nPart

private

Definition at line 90 of file EMShower.h.

Referenced by compute(), EMShower(), and prepareSteps().

nSteps

unsigned EMShower::nSteps

private

Definition at line 119 of file EMShower.h.

Referenced by compute(), and prepareSteps().

outerDepth

double EMShower::outerDepth

private

Definition at line 111 of file EMShower.h.

Referenced by prepareSteps().

photos

std::vector<double> EMShower::photos

private

Definition at line 96 of file EMShower.h.

Referenced by compute(), and EMShower().

random

const RandomEngineAndDistribution* EMShower::random

private

Definition at line 139 of file EMShower.h.

Referenced by compute(), and EMShower().

steps

Steps EMShower::steps

private

Definition at line 118 of file EMShower.h.

Referenced by compute(), and prepareSteps().

stepsCalculated

bool EMShower::stepsCalculated

private

Definition at line 120 of file EMShower.h.

Referenced by compute(), EMShower(), and prepareSteps().

T

std::vector<double> EMShower::T

private

Definition at line 97 of file EMShower.h.

theECAL

const ECALProperties* EMShower::theECAL

private

Definition at line 83 of file EMShower.h.

Referenced by compute(), and EMShower().

theGrid

EcalHitMaker* EMShower::theGrid

private

Definition at line 123 of file EMShower.h.

Referenced by compute(), prepareSteps(), and setGrid().

theHCAL

const HCALProperties* EMShower::theHCAL

private

Definition at line 84 of file EMShower.h.

Referenced by compute(), and EMShower().

theHcalHitMaker

HcalHitMaker* EMShower::theHcalHitMaker

private

Definition at line 129 of file EMShower.h.

Referenced by compute(), and setHcal().

theLayer1

const PreshowerLayer1Properties* EMShower::theLayer1

private

Definition at line 85 of file EMShower.h.

Referenced by compute(), and EMShower().

theLayer2

const PreshowerLayer2Properties* EMShower::theLayer2

private

Definition at line 86 of file EMShower.h.

Referenced by compute(), and EMShower().

theNumberOfSpots

std::vector<double> EMShower::theNumberOfSpots

private

Definition at line 93 of file EMShower.h.

Referenced by compute(), and EMShower().

theParam

EMECALShowerParametrization* const EMShower::theParam

private

Definition at line 80 of file EMShower.h.

Referenced by compute(), EMShower(), and setIntervals().

thePart

std::vector<const RawParticle*>* const EMShower::thePart

private

Definition at line 89 of file EMShower.h.

Referenced by EMShower().

thePreshower

PreshowerHitMaker* EMShower::thePreshower

private

Definition at line 126 of file EMShower.h.

Referenced by compute(), and setPreshower().

Ti

std::vector<double> EMShower::Ti

private

Definition at line 100 of file EMShower.h.

Referenced by compute(), and EMShower().

totalEnergy

double EMShower::totalEnergy

private

Definition at line 115 of file EMShower.h.

Referenced by EMShower().

TSpot

std::vector<double> EMShower::TSpot

private

Definition at line 101 of file EMShower.h.

Referenced by EMShower().