HDShower Class Reference

#include <HDShower.h>

Public 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

Public Member Functions
bool	compute ()
	Compute the shower longitudinal and lateral development. More...
int	getmip ()
	HDShower (const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart, double pmip)
virtual	~HDShower ()

Private Member Functions
double	gam (double x, double a) const
int	indexFinder (double x, const std::vector< double > &Fhist)
void	makeSteps (int nsteps)
double	transProb (double factor, double R, double r)

Private Attributes
double	aloge
double	alpEM
double	alpHD
double	balanceEH
double	betEM
double	betHD
double	criticalEnergy
double	depthECAL
double	depthGAP
double	depthGAPx0
double	depthHCAL
double	depthStart
double	depthStep
double	depthToHCAL
std::vector< int >	detector
double	e
double	eSpotSize
std::vector< double >	eStep
double	hcalDepthFactor
int	infinity
double	lambdaEM
double	lambdaHD
std::vector< double >	lamcurr
std::vector< double >	lamdepth
std::vector< double >	lamstep
std::vector< double >	lamtotal
int	lossesOpt
double	maxTRfactor
int	mip
int	nDepthSteps
std::vector< int >	nspots
int	nTRsteps
int	onEcal
double	part
const RandomEngineAndDistribution *	random
std::vector< double >	rlamStep
double	tgamEM
double	tgamHD
const ECALProperties *	theECALproperties
EcalHitMaker *	theGrid
HcalHitMaker *	theHcalHitMaker
const HCALProperties *	theHCALproperties
HDShowerParametrization *	theParam
double	theR1
double	theR2
double	theR3
double	transFactor
double	transParam
std::vector< double >	x0curr
std::vector< double >	x0depth
double	x0EM
double	x0HD

Detailed Description

Definition at line 22 of file HDShower.h.

Member Typedef Documentation

Spot

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

Definition at line 26 of file HDShower.h.

Step

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

Definition at line 27 of file HDShower.h.

step_iterator

typedef Steps::const_iterator HDShower::step_iterator

Definition at line 29 of file HDShower.h.

Steps

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

Definition at line 28 of file HDShower.h.

XYZPoint

typedef math::XYZVector HDShower::XYZPoint

Definition at line 24 of file HDShower.h.

Constructor & Destructor Documentation

HDShower()

HDShower::HDShower	(	const RandomEngineAndDistribution *	engine,
		HDShowerParametrization *	myParam,
		EcalHitMaker *	myGrid,
		HcalHitMaker *	myHcalHitMaker,
		int	onECAL,
		double	epart,
		double	pmip
	)

Definition at line 21 of file HDShower.cc.

References aloge, HDShowerParametrization::alpe1(), HDShowerParametrization::alpe2(), alpEM, HDShowerParametrization::alph1(), HDShowerParametrization::alph2(), alpHD, balanceEH, HDShowerParametrization::bete1(), HDShowerParametrization::bete2(), betEM, HDShowerParametrization::beth1(), HDShowerParametrization::beth2(), betHD, criticalEnergy, debug, depthECAL, depthGAP, depthGAPx0, depthHCAL, depthStart, depthStep, depthToHCAL, detector, e, HDShowerParametrization::e1(), HDShowerParametrization::e2(), EcalHitMaker::ecalHcalGapTotalL0(), EcalHitMaker::ecalHcalGapTotalX0(), HDShowerParametrization::ecalProperties(), EcalHitMaker::ecalTotalL0(), HDShowerParametrization::emax(), HDShowerParametrization::emid(), HDShowerParametrization::emin(), eSpotSize, RandomEngineAndDistribution::flatShoot(), HSParameters::getHDbalanceEH(), HSParameters::getHDcriticalEnergy(), HSParameters::getHDdepthStep(), HSParameters::getHDeSpotSize(), HSParameters::getHDhcalDepthFactor(), HSParameters::getHDlossesOpt(), HSParameters::getHDmaxTRfactor(), HSParameters::getHDnDepthSteps(), HSParameters::getHDnTRsteps(), HSParameters::getHDtransParam(), hcalDepthFactor, HDShowerParametrization::hcalProperties(), EcalHitMaker::hcalTotalL0(), HDShowerParametrization::hsParameters(), mps_fire::i, ECALProperties::interactionLength(), HCALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, dqm-mbProfile::log, lossesOpt, makeSteps(), HLT_2022v15_cff::maxDepth, maxTRfactor, mip, nDepthSteps, nTRsteps, onEcal, HDShowerParametrization::part1(), HDShowerParametrization::part2(), HDShowerParametrization::r1(), HDShowerParametrization::r2(), HDShowerParametrization::r3(), ECALProperties::radLenIncm(), HCALProperties::radLenIncm(), random, HDShowerParametrization::setCase(), tgamEM, tgamHD, theECALproperties, theGrid, theHCALproperties, theParam, theR1, theR2, theR3, transFactor, transParam, x0curr, x0depth, x0EM, and x0HD.

    : theParam(myParam),
      theGrid(myGrid),
      theHcalHitMaker(myHcalHitMaker),
      onEcal(onECAL),
      e(epart),
      //    pmip(pmip),
      random(engine) {
  // To get an access to constants read in FASTCalorimeter
  //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();

  // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
  lossesOpt = myParam->hsParameters()->getHDlossesOpt();
  nDepthSteps = myParam->hsParameters()->getHDnDepthSteps();
  nTRsteps = myParam->hsParameters()->getHDnTRsteps();
  transParam = myParam->hsParameters()->getHDtransParam();
  eSpotSize = myParam->hsParameters()->getHDeSpotSize();
  depthStep = myParam->hsParameters()->getHDdepthStep();
  criticalEnergy = myParam->hsParameters()->getHDcriticalEnergy();
  maxTRfactor = myParam->hsParameters()->getHDmaxTRfactor();
  balanceEH = myParam->hsParameters()->getHDbalanceEH();
  hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();

  // Special tr.size fluctuations
  transParam *= (1. + random->flatShoot());

  // Special long. fluctuations
  hcalDepthFactor += 0.05 * (2. * random->flatShoot() - 1.);

  transFactor = 1.;  // normally 1, in HF - might be smaller
                     // to take into account
                     // a narrowness of the HF shower (Cherenkov light)

  // simple protection ...
  if (e < 0)
    e = 0.;

  // Get the Famos Histos pointer
  //  myHistos = FamosHistos::instance();
  //  std::cout << " Hello FamosShower " << std::endl;

  theECALproperties = theParam->ecalProperties();
  theHCALproperties = theParam->hcalProperties();

  double emax = theParam->emax();
  double emid = theParam->emid();
  double emin = theParam->emin();
  double effective = e;

  if (e < emid) {
    theParam->setCase(1);
    // avoid "underflow" below Emin (for parameters calculation only)
    if (e < emin)
      effective = emin;
  } else
    theParam->setCase(2);

  // Avoid "overflow" beyond Emax (for parameters)
  if (effective > 0.5 * emax) {
    eSpotSize *= 2.5;
    if (effective > emax) {
      effective = emax;
      eSpotSize *= 4.;
      depthStep *= 2.;
      if (effective > 1000.)
        eSpotSize *= 2.;
    }
  }

  if (debug == 2)
    LogInfo("FastCalorimetry") << " HDShower : " << std::endl
                               << "       Energy   " << e << std::endl
                               << "      lossesOpt " << lossesOpt << std::endl
                               << "    nDepthSteps " << nDepthSteps << std::endl
                               << "       nTRsteps " << nTRsteps << std::endl
                               << "     transParam " << transParam << std::endl
                               << "      eSpotSize " << eSpotSize << std::endl
                               << " criticalEnergy " << criticalEnergy << std::endl
                               << "    maxTRfactor " << maxTRfactor << std::endl
                               << "      balanceEH " << balanceEH << std::endl
                               << "hcalDepthFactor " << hcalDepthFactor << std::endl;

  double alpEM1 = theParam->alpe1();
  double alpEM2 = theParam->alpe2();

  double betEM1 = theParam->bete1();
  double betEM2 = theParam->bete2();

  double alpHD1 = theParam->alph1();
  double alpHD2 = theParam->alph2();

  double betHD1 = theParam->beth1();
  double betHD2 = theParam->beth2();

  double part1 = theParam->part1();
  double part2 = theParam->part2();

  aloge = std::log(effective);

  double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
  double aedep = std::log(edpar);

  if (debug == 2)
    LogInfo("FastCalorimetry") << " HDShower : " << std::endl
                               << "     edpar " << edpar << "   aedep " << aedep << std::endl
                               << "    alpEM1 " << alpEM1 << std::endl
                               << "    alpEM2 " << alpEM2 << std::endl
                               << "    betEM1 " << betEM1 << std::endl
                               << "    betEM2 " << betEM2 << std::endl
                               << "    alpHD1 " << alpHD1 << std::endl
                               << "    alpHD2 " << alpHD2 << std::endl
                               << "    betHD1 " << betHD1 << std::endl
                               << "    betHD2 " << betHD2 << std::endl
                               << "     part1 " << part1 << std::endl
                               << "     part2 " << part2 << std::endl;

  // private members to set
  theR1 = theParam->r1();
  theR2 = theParam->r2();
  theR3 = theParam->r3();

  alpEM = alpEM1 + alpEM2 * aedep;
  tgamEM = tgamma(alpEM);
  betEM = betEM1 - betEM2 * aedep;
  alpHD = alpHD1 + alpHD2 * aedep;
  tgamHD = tgamma(alpHD);
  betHD = betHD1 - betHD2 * aedep;
  part = part1 - part2 * aedep;
  if (part > 1.)
    part = 1.;  // protection - just in case of

  if (debug == 2)
    LogInfo("FastCalorimetry") << " HDShower : " << std::endl
                               << "    alpEM " << alpEM << std::endl
                               << "   tgamEM " << tgamEM << std::endl
                               << "    betEM " << betEM << std::endl
                               << "    alpHD " << alpHD << std::endl
                               << "   tgamHD " << tgamHD << std::endl
                               << "    betHD " << betHD << std::endl
                               << "     part " << part << std::endl;

  if (onECAL) {
    lambdaEM = theParam->ecalProperties()->interactionLength();
    x0EM = theParam->ecalProperties()->radLenIncm();
  } else {
    lambdaEM = 0.;
    x0EM = 0.;
  }
  lambdaHD = theParam->hcalProperties()->interactionLength();
  x0HD = theParam->hcalProperties()->radLenIncm();

  if (debug == 2)
    LogInfo("FastCalorimetry") << " HDShower e " << e << std::endl
                               << "          x0EM = " << x0EM << std::endl
                               << "          x0HD = " << x0HD << std::endl
                               << "         lamEM = " << lambdaEM << std::endl
                               << "         lamHD = " << lambdaHD << std::endl;

  // Starting point of the shower
  // try first with ECAL lambda

  double sum1 = 0.;  // lambda path from the ECAL/HF entrance;
  double sum2 = 0.;  // lambda path from the interaction point;
  double sum3 = 0.;  // x0     path from the interaction point;
  int nsteps = 0;    // full number of longitudinal steps (counter);

  int nmoresteps;  // how many longitudinal steps in addition to
                   // one (if interaction happens there) in ECAL

  mip = 1;  // just to initiate particle as MIP in ECAL

  if (e < criticalEnergy)
    nmoresteps = 1;
  else
    nmoresteps = nDepthSteps;

  depthECAL = 0.;
  depthGAP = 0.;
  depthGAPx0 = 0.;
  if (onECAL) {
    depthECAL = theGrid->ecalTotalL0();  // ECAL depth segment
        //depthECAL  = theGrid->ecalTotalL0() + theGrid->ps1TotalL0() + theGrid->ps2TotalL0() + theGrid->ps2eeTotalL0(); //TEST: include preshower
    depthGAP = theGrid->ecalHcalGapTotalL0();    // GAP  depth segment
    depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
  }

  depthHCAL = theGrid->hcalTotalL0();  // HCAL depth segment
  depthToHCAL = depthECAL + depthGAP;

  //---------------------------------------------------------------------------
  // Depth simulation & various protections, among them
  // if too deep - get flat random in the allowed region
  // if no HCAL material behind - force to deposit in ECAL
  double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
  double depthStart = std::log(1. / random->flatShoot());  // starting point lambda unts
  //std::cout<<"generated depth             "<<depthStart<<std::endl;
  if (pmip == 1) {
    depthStart = depthToHCAL;
  } else {
    depthStart = depthStart * 0.9 * depthECAL / std::log(1. / pmip);
  }
  // std::cout<<"modified  depth             "<<depthStart<<std::endl;

  if (e < emin) {
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl;
    depthStart = 0.;
  }

  if (depthStart > maxDepth) {
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " << depthStart << std::endl;
    depthStart = maxDepth * random->flatShoot();
    if (depthStart < 0.)
      depthStart = 0.;
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : depthStart re-calculated = " << depthStart << std::endl;
  }

  if (onECAL && e < emid) {
    if (depthECAL > depthStep && (depthECAL - depthStart) / depthECAL > 0.2) {
      depthStart = 0.5 * depthECAL * random->flatShoot();
      if (debug)
        LogInfo("FastCalorimetry") << " FamosHDShower : small energy, "
                                   << " depthStart reduced to = " << depthStart << std::endl;
    }
  }

  if (depthHCAL < depthStep) {
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " << depthHCAL
                                 << " depthStart -> forced to 0 !!!" << std::endl;
    depthStart = 0.;
    nmoresteps = 0;

    if (depthECAL < depthStep) {
      nsteps = -1;
      LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - "
                                 << " particle is lost !!! " << std::endl;
    }
  }

  if (debug)
    LogInfo("FastCalorimetry") << " FamosHDShower  depths(lam) - " << std::endl
                               << "          ECAL = " << depthECAL << std::endl
                               << "           GAP = " << depthGAP << std::endl
                               << "          HCAL = " << depthHCAL << std::endl
                               << "  starting point = " << depthStart << std::endl;

  if (onEcal) {
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
    if (depthStart < depthECAL) {
      if (debug)
        LogInfo("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
      if (depthECAL > depthStep && (depthECAL - depthStart) / depthECAL > 0.1) {
        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step" << std::endl;

        //  ECAL - one step
        nsteps++;
        sum1 += depthECAL;  // at the end of step
        sum2 += depthECAL - depthStart;
        sum3 += sum2 * lambdaEM / x0EM;
        lamtotal.push_back(sum1);
        lamdepth.push_back(sum2);
        lamcurr.push_back(lambdaEM);
        lamstep.push_back(depthECAL - depthStart);
        x0depth.push_back(sum3);
        x0curr.push_back(x0EM);
        detector.push_back(1);
        mip = 0;

        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower : "
                                     << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;

        // Gap - no additional step after ECAL
        // just move further to HCAL over the gap
        sum1 += depthGAP;
        sum2 += depthGAP;
        sum3 += depthGAPx0;
      } else {  // Just shift starting point to HCAL
        //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;

        depthStart = depthToHCAL;
        sum1 += depthStart;
      }
    } else {  // GAP or HCAL
      if (depthStart >= depthECAL && depthStart < depthToHCAL) {
        depthStart = depthToHCAL;  // just a shift to HCAL for simplicity
      }
      sum1 += depthStart;
      if (debug)
        LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
    }
  } else {  // Forward
    if (debug)
      LogInfo("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
    sum1 += depthStart;
    transFactor = 0.5;  // makes narower tresverse size of shower
  }

  for (int i = 0; i < nmoresteps; i++) {
    sum1 += depthStep;
    if (sum1 > (depthECAL + depthGAP + depthHCAL))
      break;
    sum2 += depthStep;
    sum3 += sum2 * lambdaHD / x0HD;
    lamtotal.push_back(sum1);
    lamdepth.push_back(sum2);
    lamcurr.push_back(lambdaHD);
    lamstep.push_back(depthStep);
    x0depth.push_back(sum3);
    x0curr.push_back(x0HD);
    detector.push_back(3);
    nsteps++;
  }

  // Make fractions of energy and transverse radii at each step

  if (nsteps > 0) {
    makeSteps(nsteps);
  }
}

~HDShower()

virtual HDShower::~HDShower ( )

inlinevirtual

Definition at line 41 of file HDShower.h.

{ ; }

Member Function Documentation

compute()

bool HDShower::compute

(

)

Compute the shower longitudinal and lateral development.

Definition at line 453 of file HDShower.cc.

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), submitPVResolutionJobs::count, gather_cfg::cout, debug, detector, eStep, RandomEngineAndDistribution::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, mps_fire::i, indexFinder(), dqmiodatasetharvest::inf, infinity, dqmiolumiharvest::j, lamstep, lamtotal, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, TtFullHadEvtBuilder_cfi::prob, CosmicsPD_Skims::radius, random, mps_fire::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), mps_update::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

                       {
  //  TimeMe theT("FamosHDShower::compute");

  bool status = false;
  int numLongit = eStep.size();
  if (debug)
    LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
                               << " N_long.steps required : " << numLongit << std::endl;

  if (numLongit > 0) {
    status = true;
    // Prepare the trsanverse probability function
    std::vector<double> Fhist;
    std::vector<double> rhist;
    for (int j = 0; j < nTRsteps + 1; j++) {
      rhist.push_back(maxTRfactor * j / nTRsteps);
      Fhist.push_back(transProb(maxTRfactor, 1., rhist[j]));
      if (debug == 3)
        LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
    }

    // Longitudinal steps
    for (int i = 0; i < numLongit; i++) {
      double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
      // vary the longitudinal profile if needed
      if (detector[i] != 1)
        currentDepthL0 *= hcalDepthFactor;
      if (debug)
        LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i]
                                   << "  currentDepthL0 = " << currentDepthL0 << std::endl;

      double maxTRsize = maxTRfactor * rlamStep[i];  // in lambda units
      double rbinsize = maxTRsize / nTRsteps;
      double espot = eStep[i] / (double)nspots[i];  // re-adjust espot

      if (espot > 2. || espot < 0.)
        LogInfo("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " << espot << std::endl;

      int ecal = 0;
      if (detector[i] != 1) {
        bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);

        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of "
                                     << " theHcalHitMaker->setDepth(currentDepthL0) is " << setHDdepth << std::endl;

        if (!setHDdepth) {
          currentDepthL0 -= lamstep[i];
          setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
        }

        if (!setHDdepth)
          continue;

        theHcalHitMaker->setSpotEnergy(espot);

        //fill hcal longitudinal distribution histogram
      } else {
        ecal = 1;
        bool status = theGrid->getPads(currentDepthL0);

        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;

        if (!status)
          continue;

        int ntry = nspots[i] * 10;
        if (ntry >= infinity) {  // use max allowed in case of too many spots
          nspots[i] = 0.5 * infinity;
          espot *= 0.1 * (double)ntry / double(nspots[i]);
        } else {
          espot *= 0.1;  // fine-grain energy spots in ECAL
                         // to avoid false ECAL clustering
          nspots[i] = ntry;
        }

        theGrid->setSpotEnergy(espot);

        //fill ecal longitudinal distribution histogram
      }

      // Transverse distribition
      int nok = 0;  // counter of OK
      int count = 0;
      int inf = infinity;
      if (lossesOpt)
        inf = nspots[i];  // if losses are enabled, otherwise
      // only OK points are counted ...
      if (nspots[i] > inf)
        std::cout << " FamosHDShower::compute - at long.step " << i << "  too many spots required : " << nspots[i]
                  << " !!! " << std::endl;

      for (int j = 0; j < inf; j++) {
        if (nok == nspots[i])
          break;
        count++;

        double prob = random->flatShoot();
        int index = indexFinder(prob, Fhist);
        double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize;  // in-bin
        double phi = 2. * M_PI * random->flatShoot();

        if (debug == 2)
          LogInfo("FastCalorimetry") << std::endl
                                     << " FamosHDShower::compute "
                                     << " r = " << radius << "    phi = " << phi << std::endl;

        bool result;
        if (ecal) {
          result = theGrid->addHit(radius, phi, 0);

          if (debug == 2)
            LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
                                       << " theGrid->addHit result = " << result << std::endl;

          //fill ecal transverse distribution histogram
        } else {
          result = theHcalHitMaker->addHit(radius, phi, 0);

          if (debug == 2)
            LogInfo("FastCalorimetry") << " FamosHDShower::compute - "
                                       << " theHcalHitMaker->addHit result = " << result << std::endl;

          //fill hcal transverse distribution histogram
        }

        if (result)
          nok++;

      }  // end of tranverse simulation

      if (count == infinity) {
        if (debug)
          LogInfo("FastCalorimetry") << " FamosHDShower::compute "
                                     << " maximum number of"
                                     << " transverse points " << count << " is used !!!" << std::endl;
      }

      if (debug)
        LogInfo("FastCalorimetry") << " FamosHDShower::compute "
                                   << " long.step No." << i << "   Ntry, Nok = " << count << " " << nok << std::endl;
    }  // end of longitudinal steps
  }    // end of no steps

  return status;
}

gam()

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

inlineprivate

Definition at line 48 of file HDShower.h.

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

Referenced by makeSteps().

{ return pow(x, a - 1.) * exp(-x); }

getmip()

int HDShower::getmip ( )

inline

Definition at line 39 of file HDShower.h.

References mip.

Referenced by CalorimetryManager::HDShowerSimulation().

{ return mip; }

indexFinder()

int HDShower::indexFinder ( double  x, const std::vector< double > &  Fhist  )

private

Definition at line 601 of file HDShower.cc.

References debug, findQualityFiles::size, and x.

Referenced by compute().

                                                                  {
  // binary search in the vector of doubles
  int size = Fhist.size();

  int curr = size / 2;
  int step = size / 4;
  int iter;
  int prevdir = 0;
  int actudir = 0;

  for (iter = 0; iter < size; iter++) {
    if (curr >= size || curr < 1)
      LogWarning("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " << curr << " !!!"
                                    << std::endl;

    if ((x <= Fhist[curr]) && (x > Fhist[curr - 1]))
      break;
    prevdir = actudir;
    if (x > Fhist[curr]) {
      actudir = 1;
    } else {
      actudir = -1;
    }
    if (prevdir * actudir < 0) {
      if (step > 1)
        step /= 2;
    }
    curr += actudir * step;
    if (curr > size)
      curr = size;
    else {
      if (curr < 1) {
        curr = 1;
      }
    }

    if (debug == 3)
      LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter << " curr, F[curr-1], F[curr] = " << curr
                                 << " " << Fhist[curr - 1] << " " << Fhist[curr] << std::endl;
  }

  if (debug == 3)
    LogInfo("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr - 1 << std::endl;

  return curr - 1;
}

makeSteps()

void HDShower::makeSteps ( int  nsteps )

private

Definition at line 355 of file HDShower.cc.

References aloge, alpEM, alpHD, balanceEH, betEM, betHD, submitPVResolutionJobs::count, debug, detector, e, eSpotSize, eStep, RandomEngineAndDistribution::flatShoot(), gam(), mps_fire::i, infinity, createfilelist::int, lamcurr, lamdepth, lamstep, nspots, random, rlamStep, groupFilesInBlocks::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, x, x0curr, x0depth, and y.

Referenced by HDShower().

                                   {
  double sumes = 0.;
  double sum = 0.;
  std::vector<double> temp;

  if (debug)
    LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - "
                               << " nsteps required : " << nsteps << std::endl;

  int count = 0;
  for (int i = 0; i < nsteps; i++) {
    double deplam = lamdepth[i] - 0.5 * lamstep[i];
    double depx0 = x0depth[i] - 0.5 * lamstep[i] / x0curr[i];
    double x = betEM * depx0;
    double y = betHD * deplam;

    if (debug == 2)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps "
                                 << " - step " << i << "   depx0, x = " << depx0 << ", " << x
                                 << "   deplam, y = " << deplam << ", " << y << std::endl;

    double est = (part * betEM * gam(x, alpEM) * lamcurr[i] / (x0curr[i] * tgamEM) +
                  (1. - part) * betHD * gam(y, alpHD) / tgamHD) *
                 lamstep[i];

    // protection ...
    if (est < 0.) {
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps "
                                 << " - negative step energy !!!" << std::endl;
      est = 0.;
      break;
    }

    // for estimates only
    sum += est;
    int nPest = (int)(est * e / sum / eSpotSize);

    if (debug == 2)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " << nPest << std::endl;

    if (nPest <= 1 && count != 0)
      break;

    // good step - to proceed

    temp.push_back(est);
    sumes += est;

    rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge)) * deplam * transFactor);
    count++;
  }

  // fluctuations in ECAL and re-distribution of remaining energy in HCAL
  if (detector[0] == 1 && count > 1) {
    double oldECALenergy = temp[0];
    double oldHCALenergy = sumes - oldECALenergy;
    double newECALenergy = 2. * sumes;
    for (int i = 0; newECALenergy > sumes && i < infinity; i++)
      newECALenergy = 2. * balanceEH * random->flatShoot() * oldECALenergy;

    if (debug == 2)
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps "
                                 << " ECAL fraction : old/new - " << oldECALenergy / sumes << "/"
                                 << newECALenergy / sumes << std::endl;

    temp[0] = newECALenergy;
    double newHCALenergy = sumes - newECALenergy;
    double newHCALreweight = newHCALenergy / oldHCALenergy;

    for (int i = 1; i < count; i++) {
      temp[i] *= newHCALreweight;
    }
  }

  // final re-normalization of the energy fractions
  for (int i = 0; i < count; i++) {
    eStep.push_back(temp[i] * e / sumes);
    nspots.push_back((int)(eStep[i] / eSpotSize) + 1);

    if (debug)
      LogInfo("FastCalorimetry") << " step " << i << "  det: " << detector[i]
                                 << "  xO and lamdepth at the end of step = " << x0depth[i] << " " << lamdepth[i]
                                 << "   Estep func = " << eStep[i] << "   Rstep = " << rlamStep[i]
                                 << "  Nspots = " << nspots[i] << "  espot = " << eStep[i] / (double)nspots[i]
                                 << std::endl;
  }

  // The only step is in ECAL - let's make the size bigger ...
  if (count == 1 and detector[0] == 1)
    rlamStep[0] *= 2.;

  if (debug) {
    if (eStep[0] > 0.95 * e && detector[0] == 1)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - "
                                 << "ECAL energy = " << eStep[0] << " out of total = " << e << std::endl;
  }
}

transProb()

double HDShower::transProb ( double  factor, double  R, double  r  )

inlineprivate

Definition at line 53 of file HDShower.h.

References dttmaxenums::R, and alignCSCRings::r.

Referenced by compute().

                                                      {
    double fsq = factor * factor;
    return ((fsq + 1.) / fsq) * r * r / (r * r + R * R);
  }

Member Data Documentation

aloge

double HDShower::aloge

private

Definition at line 76 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

alpEM

double HDShower::alpEM

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

alpHD

double HDShower::alpHD

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

balanceEH

double HDShower::balanceEH

private

Definition at line 119 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

betEM

double HDShower::betEM

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

betHD

double HDShower::betHD

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

criticalEnergy

double HDShower::criticalEnergy

private

Definition at line 115 of file HDShower.h.

Referenced by HDShower().

depthECAL

double HDShower::depthECAL

private

Definition at line 127 of file HDShower.h.

Referenced by HDShower().

depthGAP

double HDShower::depthGAP

private

Definition at line 127 of file HDShower.h.

Referenced by HDShower().

depthGAPx0

double HDShower::depthGAPx0

private

Definition at line 127 of file HDShower.h.

Referenced by HDShower().

depthHCAL

double HDShower::depthHCAL

private

Definition at line 127 of file HDShower.h.

Referenced by HDShower().

depthStart

double HDShower::depthStart

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower().

depthStep

double HDShower::depthStep

private

Definition at line 113 of file HDShower.h.

Referenced by HDShower().

depthToHCAL

double HDShower::depthToHCAL

private

Definition at line 127 of file HDShower.h.

Referenced by HDShower().

detector

std::vector<int> HDShower::detector

private

Definition at line 78 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

e

double HDShower::e

private

Definition at line 98 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

eSpotSize

double HDShower::eSpotSize

private

Definition at line 111 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

eStep

std::vector<double> HDShower::eStep

private

Definition at line 79 of file HDShower.h.

Referenced by compute(), and makeSteps().

hcalDepthFactor

double HDShower::hcalDepthFactor

private

Definition at line 121 of file HDShower.h.

Referenced by compute(), and HDShower().

infinity

int HDShower::infinity

private

Definition at line 83 of file HDShower.h.

Referenced by compute(), and makeSteps().

lambdaEM

double HDShower::lambdaEM

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower().

lambdaHD

double HDShower::lambdaHD

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower().

lamcurr

std::vector<double> HDShower::lamcurr

private

Definition at line 81 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

lamdepth

std::vector<double> HDShower::lamdepth

private

Definition at line 81 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

lamstep

std::vector<double> HDShower::lamstep

private

Definition at line 81 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

lamtotal

std::vector<double> HDShower::lamtotal

private

Definition at line 81 of file HDShower.h.

Referenced by compute(), and HDShower().

lossesOpt

int HDShower::lossesOpt

private

Definition at line 101 of file HDShower.h.

Referenced by compute(), and HDShower().

maxTRfactor

double HDShower::maxTRfactor

private

Definition at line 117 of file HDShower.h.

Referenced by compute(), and HDShower().

mip

int HDShower::mip

private

Definition at line 95 of file HDShower.h.

Referenced by getmip(), and HDShower().

nDepthSteps

int HDShower::nDepthSteps

private

Definition at line 103 of file HDShower.h.

Referenced by HDShower().

nspots

std::vector<int> HDShower::nspots

private

Definition at line 78 of file HDShower.h.

Referenced by compute(), and makeSteps().

nTRsteps

int HDShower::nTRsteps

private

Definition at line 105 of file HDShower.h.

Referenced by compute(), and HDShower().

onEcal

int HDShower::onEcal

private

Definition at line 92 of file HDShower.h.

Referenced by HDShower().

part

double HDShower::part

private

Definition at line 71 of file HDShower.h.

random

const RandomEngineAndDistribution* HDShower::random

private

Definition at line 124 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

rlamStep

std::vector<double> HDShower::rlamStep

private

Definition at line 79 of file HDShower.h.

Referenced by compute(), and makeSteps().

tgamEM

double HDShower::tgamEM

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

tgamHD

double HDShower::tgamHD

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

theECALproperties

const ECALProperties* HDShower::theECALproperties

private

Definition at line 66 of file HDShower.h.

Referenced by HDShower().

theGrid

EcalHitMaker* HDShower::theGrid

private

Definition at line 86 of file HDShower.h.

Referenced by compute(), and HDShower().

theHcalHitMaker

HcalHitMaker* HDShower::theHcalHitMaker

private

Definition at line 89 of file HDShower.h.

Referenced by compute().

theHCALproperties

const HCALProperties* HDShower::theHCALproperties

private

Definition at line 67 of file HDShower.h.

Referenced by HDShower().

theParam

HDShowerParametrization* HDShower::theParam

private

Definition at line 63 of file HDShower.h.

Referenced by HDShower().

theR1

double HDShower::theR1

private

Definition at line 70 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

theR2

double HDShower::theR2

private

Definition at line 70 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

theR3

double HDShower::theR3

private

Definition at line 70 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

transFactor

double HDShower::transFactor

private

Definition at line 109 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

transParam

double HDShower::transParam

private

Definition at line 107 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

x0curr

std::vector<double> HDShower::x0curr

private

Definition at line 80 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

x0depth

std::vector<double> HDShower::x0depth

private

Definition at line 80 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

x0EM

double HDShower::x0EM

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower().

x0HD

double HDShower::x0HD

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower().