HFShower Class Reference

#include <HFShower.h>

Classes
struct	Hit

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...
std::vector< Hit >	getHits (G4Step *aStep, double weight)
std::vector< Hit >	getHits (G4Step *aStep, bool forLibrary)
std::vector< Hit >	getHits (G4Step *aStep, bool forLibraryProducer, double zoffset)
	HFShower (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk=0)
	HFShower (const RandomEngineAndDistribution *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
void	initRun (G4ParticleTable *)
virtual	~HFShower ()
virtual	~HFShower ()

Private Member Functions
double	gam (double x, double a) const
std::vector< double >	getDDDArray (const std::string &, const DDsvalues_type &, int &)
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
bool	applyFidCut
double	balanceEH
double	betEM
double	betHD
HFCherenkov *	cherenkov
int	chkFibre
double	criticalEnergy
double	depthStart
double	depthStep
std::vector< int >	detector
double	e
double	eSpotSize
std::vector< double >	eStep
HFFibre *	fibre
std::vector< double >	gpar
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	nDepthSteps
std::vector< int >	nspots
int	nTRsteps
int	onEcal
double	part
double	probMax
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 HFShower.h.

Member Typedef Documentation

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

Definition at line 29 of file HFShower.h.

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

Definition at line 30 of file HFShower.h.

typedef Steps::const_iterator HFShower::step_iterator

Definition at line 32 of file HFShower.h.

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

Definition at line 31 of file HFShower.h.

typedef math::XYZVector HFShower::XYZPoint

Definition at line 27 of file HFShower.h.

Constructor & Destructor Documentation

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

Definition at line 30 of file HFShower.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, depthStart, depthStep, 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(), i, ECALProperties::interactionLength(), HCALProperties::interactionLength(), lambdaEM, lambdaHD, lamcurr, lamdepth, lamstep, lamtotal, create_public_lumi_plots::log, LogDebug, lossesOpt, makeSteps(), maxTRfactor, 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),
    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 ad hoc long. extension + some fluctuations
  double   depthExt;
  if (e < 50.)   depthExt =  0.8 * (50.  - e) / 50.  + 0.3;
  else {
    if (e < 500.)  depthExt =  (500. - e) / 500. * 0.4 - 0.1;
    else   depthExt =  -0.1;
  }
  hcalDepthFactor += depthExt + 0.05 * (2.* random->flatShoot() - 1.);

  // normally 1, in HF - might be smaller to take into account
  // a narrowness of the HF shower (Cherenkov light)
  if( e < 50. )  transFactor = 0.5 - (50.   - e ) / 50.  * 0.2;
  else           transFactor = 0.7 - (1000. - e ) /1000. * 0.2;

  // 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);

  // A bit coarse espot size for HF...
  eSpotSize *= 2.5;
  if(effective > 0.5 * emax) {
    eSpotSize *= 2.;
    if(effective > emax) {
      effective = emax;
      eSpotSize *= 3.;
      depthStep *= 2.;
     }
  }

  if(debug == 2 )
    LogDebug("FastCalorimetry") <<  " HFShower : " << 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)
    LogDebug("FastCalorimetry") << " HFShower : " << 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 )
    LogDebug("FastCalorimetry") << " HFShower : " << 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)
    LogDebug("FastCalorimetry") << " HFShower 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


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

  double depthECAL  = 0.;
  double depthGAP   = 0.;
  double depthGAPx0 = 0.;
  if(onECAL ) {
    depthECAL  = theGrid->ecalTotalL0();         // ECAL depth segment
    depthGAP   = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
    depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
  }

  double depthHCAL   = theGrid->hcalTotalL0();    // HCAL depth segment
  double 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

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

  if(depthStart > maxDepth) {
    if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart too big ...   = "
           << depthStart << std::endl;

    depthStart = maxDepth *  random->flatShoot();
    if( depthStart < 0.) depthStart = 0.;
    if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart re-calculated = "
           << depthStart << std::endl;
  }

  if( onECAL && e < emid ) {
    if((depthECAL - depthStart)/depthECAL > 0.2 && depthECAL > depthStep ) {

      depthStart = 0.5 * depthECAL * random->flatShoot();
      if(debug)
        LogDebug("FastCalorimetry") << " FamosHFShower : small energy, "
         << " depthStart reduced to = " << depthStart << std::endl;

    }
  }

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

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




  if(debug)
    LogDebug("FastCalorimetry") << " FamosHFShower  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) LogDebug("FastCalorimetry") << " FamosHFShower : onECAL" << std::endl;
    if(depthStart < depthECAL) {
      if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : depthStart < depthECAL" << std::endl;
      if((depthECAL - depthStart)/depthECAL > 0.25 && depthECAL > depthStep) {
    if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : 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);

    if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : " << " 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;
      }
      // Just shift starting point to HCAL
      else {
    //  cout << " FamosHFShower : not enough space to make ECAL step" << std::endl;
    if(debug)  LogDebug("FastCalorimetry") << " FamosHFShower : 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) LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
    }
  }
  else {   // Forward
    if(debug)  LogDebug("FastCalorimetry") << " FamosHFShower : forward" << std::endl;
    sum1 += depthStart;
  }

  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

  // PV
  //  std::cout << "HFShower::HFShower() : Nsteps = " << nsteps << std::endl;

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

}

HFShower::~HFShower ( )

inlinevirtual

Definition at line 41 of file HFShower.h.

{;}

HFShower::HFShower	(	std::string &	name,
		const DDCompactView &	cpv,
		edm::ParameterSet const &	p,
		int	chk = 0
	)

Definition at line 26 of file HFShower.cc.

References DDFilteredView::addFilter(), applyFidCut, cherenkov, chkFibre, DDSpecificsFilter::equals, edm::hlt::Exception, fibre, alcazmumu_cfi::filter, DDFilteredView::firstChild(), getDDDArray(), edm::ParameterSet::getParameter(), gpar, DDFilteredView::mergedSpecifics(), mergeVDriftHistosByStation::name, probMax, DDSpecificsFilter::setCriteria(), and relativeConstraints::value.

                                                       : cherenkov(0),
                                                           fibre(0),
                                                           chkFibre(chk) {

  edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
  applyFidCut            = m_HF.getParameter<bool>("ApplyFiducialCut");
  probMax                = m_HF.getParameter<double>("ProbMax");

  edm::LogInfo("HFShower") << "HFShower:: Maximum probability cut off " 
                           << probMax << " Check flag " << chkFibre;

  G4String attribute = "ReadOutName";
  G4String value     = name;
  DDSpecificsFilter filter;
  DDValue           ddv(attribute,value,0);
  filter.setCriteria(ddv,DDSpecificsFilter::equals);
  DDFilteredView fv(cpv);
  fv.addFilter(filter);
  bool dodet = fv.firstChild();
  if ( dodet ) {
    DDsvalues_type sv(fv.mergedSpecifics());

    //Special Geometry parameters
    int ngpar = 7;
    gpar      = getDDDArray("gparHF", sv,ngpar);
    edm::LogInfo("HFShower") << "HFShower: " << ngpar << " gpar (cm)";
    for (int ig=0; ig<ngpar; ig++)
      edm::LogInfo("HFShower") << "HFShower: gpar[" << ig << "] = "
                               << gpar[ig]/cm << " cm";
  } else {
    edm::LogError("HFShower") << "HFShower: cannot get filtered "
                              << " view for " << attribute << " matching " << name;
    throw cms::Exception("Unknown", "HFShower")
      << "cannot match " << attribute << " to " << name <<"\n";
  }
  cherenkov = new HFCherenkov(m_HF);
  fibre     = new HFFibre(name, cpv, p);
}

virtual HFShower::~HFShower ( )

virtual

Member Function Documentation

bool HFShower::compute

(

)

Compute the shower longitudinal and lateral development.

Definition at line 475 of file HFShower.cc.

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), prof2calltree::count, debug, detector, patCandidatesForDimuonsSequences_cff::ecal, eStep, RandomEngineAndDistribution::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, getHLTprescales::index, indexFinder(), EcalCondDB::inf, infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, mix_2012_Summer_inTimeOnly_cff::prob, CosmicsPD_Skims::radius, random, query::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), ntuplemaker::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

                       {

  //  TimeMe theT("FamosHFShower::compute");

  bool status = false;
  int numLongit = eStep.size();
  if(debug)
    LogDebug("FastCalorimetry") << " FamosHFShower::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)
    LogDebug("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)
    LogDebug("FastCalorimetry") << " FamosHFShower::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 > 4. || espot < 0. )
    LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = "
         << espot << std::endl;

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

    if(debug)
      LogDebug("FastCalorimetry") << " FamosHFShower::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);
      }
      else {

    ecal = 1;
    bool status = theGrid->getPads(currentDepthL0);

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

    if(!status) continue;

    theGrid->setSpotEnergy(espot);
      }


      //------------------------------------------------------------
      // Transverse distribution
      //------------------------------------------------------------
      int nok   = 0;                          // counter of OK
      int count = 0;
      int inf   = infinity;
      if(lossesOpt) inf = nspots[i];          // losses are enabled, otherwise
      // only OK points are counted ...

      // Total energy in this layer
      double eremaining = eStep[i];
      bool converged    = false;

      while (eremaining > 0. && !converged && count<inf ) {

    ++count;

    // energy spot  (HFL)
        double newespot =  espot;

    // We need to know a priori if this energy spot if for
        // a long (1) or short (2) fiber

        unsigned layer = 1;
        if( currentDepthL0 < 1.3 ) // first 22 cm = 1.3 lambda - only HFL
      layer = 1;
        else
      layer = random->flatShoot() < 0.5 ? 1 : 2;

    if ( layer == 2 ) newespot = 2. * espot;

    if ( eremaining - newespot < 0. ) newespot = eremaining;

    // process transverse distribution

    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)
      LogDebug("FastCalorimetry") << std::endl << " FamosHFShower::compute " << " r = " << radius
           << "    phi = " << phi << std::endl;

    // add hit

    theHcalHitMaker->setSpotEnergy(newespot);
    theGrid->setSpotEnergy(newespot);

    bool result;
    if(ecal) {
      result = theGrid->addHit(radius,phi,0);  // shouldn't get here !

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

    }
    else {
      // PV assign espot to long/short fibers
          result = theHcalHitMaker->addHit(radius,phi,layer);

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

    }

    if (result) {
      ++nok;
      eremaining -= newespot;
      if ( eremaining <= 0. ) converged = true;
      //      std::cout << "Transverse : " << nok << " "
      //            << " , E= "     << newespot
      //            << " , Erem = " << eremaining
      //            << std::endl;
    } else {
      //      std::cout << "WARNING : hit not added" << std::endl;
    }
      }

    // end of tranverse simulation
    //-----------------------------------------------------

      if(count == infinity) {
        status = false;
    if(debug)
      LogDebug("FastCalorimetry") << "*** FamosHFShower::compute " << " maximum number of"
           << " transverse points " << count << " is used !!!" << std::endl;
        break;
      }

      if(debug)
    LogDebug("FastCalorimetry") << " FamosHFShower::compute " << " long.step No." << i
         << "   Ntry, Nok = " << count
         << " " << nok << std::endl;

    } // end of longitudinal steps
  } // end of no steps
  return status;

}

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

inlineprivate

Definition at line 49 of file HFShower.h.

References create_public_lumi_plots::exp, and funct::pow().

Referenced by makeSteps().

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

std::vector< double > HFShower::getDDDArray ( const std::string &  str, const DDsvalues_type &  sv, int &  nmin  )

private

Definition at line 440 of file HFShower.cc.

References DDfetch(), DDValue::doubles(), edm::hlt::Exception, LogDebug, and relativeConstraints::value.

Referenced by HFShower().

                                  {
#ifdef DebugLog
  LogDebug("HFShower") << "HFShower:getDDDArray called for " << str 
                       << " with nMin " << nmin;
#endif
  DDValue value(str);
  if (DDfetch(&sv, value)) {
#ifdef DebugLog
    LogDebug("HFShower") << value;
#endif
    const std::vector<double> & fvec = value.doubles();
    int nval = fvec.size();
    if (nmin > 0) {
      if (nval < nmin) {
        edm::LogError("HFShower") << "HFShower : # of " << str << " bins "
                                  << nval << " < " << nmin  << " ==> illegal";
        throw cms::Exception("Unknown", "HFShower")
                                  << "nval < nmin for array " << str << "\n";
      }
    } else {
      if (nval < 2) {
        edm::LogError("HFShower") << "HFShower : # of " << str << " bins " << nval
                                  << " < 2 ==> illegal" << " (nmin=" << nmin << ")";
        throw cms::Exception("Unknown", "HFShower")
      << "nval < 2 for array " << str << "\n";
      }
    }
    nmin = nval;

    return fvec;
  } else {
    edm::LogError("HFShower") << "HFShower : cannot get array " << str;
    throw cms::Exception("Unknown", "HFShower") << "cannot get array " << str << "\n";
  }
}

std::vector< HFShower::Hit > HFShower::getHits	(	G4Step *	aStep,
		double	weight
	)

Definition at line 71 of file HFShower.cc.

References funct::abs(), applyFidCut, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), HFShower::Hit::depth, relval_parameters_module::energy, create_public_lumi_plots::exp, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, probMax, diffTwoXMLs::r1, diffTwoXMLs::r2, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), findQualityFiles::v, w(), and HFShower::Hit::wavelength.

Referenced by HCalSD::hitForFibre(), and FiberSD::ProcessHits().

                                                                      {

  std::vector<HFShower::Hit> hits;
  int    nHit    = 0;

  double edep    = weight*(aStep->GetTotalEnergyDeposit());
#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits: energy " << aStep->GetTotalEnergyDeposit() << " weight " << weight << " edep " << edep;
#endif
  double stepl   = 0.;

  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
    stepl = aStep->GetStepLength();
  if ((edep == 0.) || (stepl == 0.)) {
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
#endif
    return hits;
  }
  G4Track *aTrack = aStep->GetTrack();
  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
 
  HFShower::Hit hit;
  double energy   = aParticle->GetTotalEnergy();
  double momentum = aParticle->GetTotalMomentum();
  double pBeta    = momentum / energy;
  double dose     = 0.;
  int    npeDose  = 0;

  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
     
  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector globalPos    = preStepPoint->GetPosition();
  G4String      name         = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
  //double        zv           = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
  double        zv           = std::abs(globalPos.z()) - gpar[4];
  G4ThreeVector localPos     = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
  G4ThreeVector localMom     = preStepPoint->GetTouchable()->GetHistory()->
                               GetTopTransform().TransformAxis(momentumDir);
  // @@ Here the depth should be changed  (Fibers all long in Geometry!)
  int  depth = 1;
  int  npmt  = 0;
  bool ok    = true;
  if (zv < 0. || zv > gpar[1]) {
    ok = false; // beyond the fiber in Z
  }
  if (ok && applyFidCut) {
    npmt = HFFibreFiducial::PMTNumber(globalPos);
    if (npmt <= 0) {
      ok = false;
    } else if (npmt > 24) { // a short fibre
      if (zv > gpar[0]) {
    depth = 2; 
      } else {
        ok = false;
      }
    }
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower: npmt " << npmt 
                             << " zv " << std::abs(globalPos.z()) 
                             << ":" << gpar[4] << ":" << zv << ":" 
                             << gpar[0] << " ok " << ok << " depth " << depth;
#endif
  } else {
    depth      = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
  }
  G4ThreeVector translation =
                        preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();

  double u        = localMom.x();
  double v        = localMom.y();
  double w        = localMom.z();
  double zCoor    = localPos.z();
  double zFibre   = (0.5*gpar[1]-zCoor-translation.z());
  double tSlice   = (aStep->GetPostStepPoint()->GetGlobalTime());
  double time     = fibre->tShift(localPos, depth, chkFibre);

#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits: in " << name << " Z " 
               << zCoor << "(" << globalPos.z() << ") " << zFibre 
               << " Trans " << translation << " Time " << tSlice  
               << " " << time << "\n                  Direction " 
               << momentumDir << " Local " << localMom;
#endif 
  // Here npe should be 0 if there is no fiber (@@ M.K.)
  int npe = 1;
  std::vector<double> wavelength;
  std::vector<double> momz;
  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
    if (ok) npe = cherenkov->computeNPE(aStep,particleDef,pBeta,u,v,w,stepl,zFibre,dose, npeDose);
    wavelength = cherenkov->getWL();
    momz       = cherenkov->getMom();
  } // ^^^^^ End of Tmp close of the cherenkov function
  if(ok && npe>0) {
    for (int i = 0; i<npe; ++i) {
      double p=1.;
      if (!applyFidCut) p   = fibre->attLength(wavelength[i]); 
      double r1  = G4UniformRand();
      double r2  = G4UniformRand();
#ifdef DebugLog
      edm::LogInfo("HFShower") << "HFShower::getHits: " << i << " attenuation "
                   << r1 <<":" << exp(-p*zFibre) << " r2 " << r2 
                   << ":" << probMax << " Survive: " 
                   << (r1 <= exp(-p*zFibre) && r2 <= probMax);
#endif
      if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
    hit.depth      = depth;
    hit.time       = tSlice+time;
    if (!applyFidCut) {
      hit.wavelength = wavelength[i]; // Tmp
      hit.momentum   = momz[i]; // Tmp
    } else {
      hit.wavelength = 300.; // Tmp
      hit.momentum   = 1.; // Tmp
    }
    hit.position   = globalPos;
    hits.push_back(hit);
    nHit++;
      }
    }
  }

#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
  for (int i=0; i<nHit; ++i)
    edm::LogInfo("HFShower") << "HFShower::Hit " << i << " WaveLength " 
                 << hits[i].wavelength << " Time " << hits[i].time
                 << " Momentum " << hits[i].momentum <<" Position "
                 << hits[i].position;
#endif
  return hits;

} 

std::vector< HFShower::Hit > HFShower::getHits	(	G4Step *	aStep,
		bool	forLibrary
	)

Definition at line 341 of file HFShower.cc.

References funct::abs(), applyFidCut, cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), findQualityFiles::v, w(), and HFShower::Hit::wavelength.

                                                                        {
  std::vector<HFShower::Hit> hits;
  int    nHit    = 0;

  double edep    = aStep->GetTotalEnergyDeposit();
  double stepl   = 0.;

  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
    stepl = aStep->GetStepLength();
  if ((edep == 0.) || (stepl == 0.)) {
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
#endif
    return hits;
  }

  G4Track *aTrack = aStep->GetTrack();
  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();

  HFShower::Hit hit;
  double energy   = aParticle->GetTotalEnergy();
  double momentum = aParticle->GetTotalMomentum();
  double pBeta    = momentum / energy;
  double dose     = 0.;
  int    npeDose  = 0;

  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
  G4ParticleDefinition *particleDef = aTrack->GetDefinition();

  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector globalPos    = preStepPoint->GetPosition();
  G4String      name         = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
  double        zv           = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
  G4ThreeVector localPos     = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
  G4ThreeVector localMom     = preStepPoint->GetTouchable()->GetHistory()->
    GetTopTransform().TransformAxis(momentumDir);
  // @@ Here the depth should be changed (Fibers are all long in Geometry!)
  int depth  = 1;
  int npmt   = 0;
  bool ok    = true;
  if (zv < 0 || zv > gpar[1]) {
    ok = false; // beyond the fiber in Z
  }
  if (ok && applyFidCut) {
    npmt = HFFibreFiducial:: PMTNumber(globalPos);
    if (npmt <= 0) {
      ok = false;
    } else if (npmt > 24) { // a short fibre
      if (zv > gpar[0]) {
    depth = 2; 
      } else {
        ok = false;
      }
    }
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower:getHits(SL): npmt " << npmt 
                             << " zv " << std::abs(globalPos.z()) 
                             << ":" << gpar[4] << ":" << zv << ":" 
                             << gpar[0] << " ok " << ok << " depth " << depth;
#endif
  } else {
    depth      = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
  }
  G4ThreeVector translation  =
                        preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();

  double u        = localMom.x();
  double v        = localMom.y();
  double w        = localMom.z();
  double zCoor    = localPos.z();
  double zFibre   = (0.5*gpar[1]-zCoor-translation.z());
  double tSlice   = (aStep->GetPostStepPoint()->GetGlobalTime());
  double time     = fibre->tShift(localPos, depth, chkFibre);

#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits(SL): in " << name << " Z " 
               << zCoor << "(" << globalPos.z() << ") " << zFibre 
               << " Trans " << translation << " Time " << tSlice  
               << " " << time << "\n                  Direction " 
               << momentumDir << " Local " << localMom;
#endif
  // npe should be 0
  int npe = 0;
  if(ok) npe = cherenkov->computeNPE(aStep,particleDef,pBeta,u,v,w, stepl,zFibre, dose, npeDose);
  std::vector<double> wavelength = cherenkov->getWL();
  std::vector<double> momz       = cherenkov->getMom();

  for (int i = 0; i<npe; ++i) {
    hit.time = tSlice+time;
    hit.wavelength = wavelength[i];
    hit.momentum   = momz[i];
    hit.position   = globalPos;
    hits.push_back(hit);
    nHit++;
  }

  return hits;
}

std::vector< HFShower::Hit > HFShower::getHits	(	G4Step *	aStep,
		bool	forLibraryProducer,
		double	zoffset
	)

Definition at line 206 of file HFShower.cc.

References funct::abs(), applyFidCut, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), HFShower::Hit::depth, relval_parameters_module::energy, create_public_lumi_plots::exp, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, HFShower::Hit::momentum, mergeVDriftHistosByStation::name, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, HFFibreFiducial::PMTNumber(), HFShower::Hit::position, probMax, diffTwoXMLs::r1, diffTwoXMLs::r2, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), findQualityFiles::v, w(), and HFShower::Hit::wavelength.

                                         {

  std::vector<HFShower::Hit> hits;
  int    nHit    = 0;

  double edep    = aStep->GetTotalEnergyDeposit();
  double stepl   = 0.;

  if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.)
    stepl = aStep->GetStepLength();
  if ((edep == 0.) || (stepl == 0.)) {
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
#endif
    return hits;
  }
  G4Track *aTrack = aStep->GetTrack();
  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
 
  HFShower::Hit hit;
  double energy   = aParticle->GetTotalEnergy();
  double momentum = aParticle->GetTotalMomentum();
  double pBeta    = momentum / energy;
  double dose     = 0.;
  int    npeDose  = 0;

  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
     
  G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector globalPos    = preStepPoint->GetPosition();
  G4String      name         = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
  //double        zv           = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
  //double        zv           = std::abs(globalPos.z()) - gpar[4];
  double        zv           = gpar[1]-(std::abs(globalPos.z()) - zoffset);
  G4ThreeVector localPos     = G4ThreeVector(globalPos.x(),globalPos.y(), zv);
  G4ThreeVector localMom     = preStepPoint->GetTouchable()->GetHistory()->
                               GetTopTransform().TransformAxis(momentumDir);
  // @@ Here the depth should be changed  (Fibers all long in Geometry!)
  int  depth = 1;
  int  npmt  = 0;
  bool ok    = true;
  if (zv < 0. || zv > gpar[1]) {
    ok = false; // beyond the fiber in Z
  }
  if (ok && applyFidCut) {
    npmt = HFFibreFiducial::PMTNumber(globalPos);
    if (npmt <= 0) {
      ok = false;
    } else if (npmt > 24) { // a short fibre
      if (zv > gpar[0]) {
    depth = 2; 
      } else {
        ok = false;
      }
    }
#ifdef DebugLog
    edm::LogInfo("HFShower") << "HFShower: npmt " << npmt 
                             << " zv " << std::abs(globalPos.z()) 
                             << ":" << gpar[4] << ":" << zv << ":" 
                             << gpar[0] << " ok " << ok << " depth " << depth;
#endif
  } else {
    depth      = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10; // All LONG!
  }
  G4ThreeVector translation =
                        preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();

  double u        = localMom.x();
  double v        = localMom.y();
  double w        = localMom.z();
  double zCoor    = localPos.z();
  double zFibre   = (0.5*gpar[1]-zCoor-translation.z());
  double tSlice   = (aStep->GetPostStepPoint()->GetGlobalTime());
  double time     = fibre->tShift(localPos, depth, chkFibre);

#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits: in " << name << " Z " <<zCoor
               << "(" << globalPos.z() <<") " << zFibre <<" Trans "
               << translation << " Time " << tSlice  << " " << time
               << "\n                  Direction " << momentumDir 
               << " Local " << localMom;
#endif 
  // Here npe should be 0 if there is no fiber (@@ M.K.)
  int npe = 1;
  std::vector<double> wavelength;
  std::vector<double> momz;
  if (!applyFidCut) { // _____ Tmp close of the cherenkov function
    if (ok) npe = cherenkov->computeNPE(aStep,particleDef,pBeta,u,v,w,stepl,zFibre,dose, npeDose);
    wavelength = cherenkov->getWL();
    momz       = cherenkov->getMom();
  } // ^^^^^ End of Tmp close of the cherenkov function
  if (ok && npe>0) {
    for (int i = 0; i<npe; ++i) {
      double p=1.;
      if (!applyFidCut) p   = fibre->attLength(wavelength[i]); 
      double r1  = G4UniformRand();
      double r2  = G4UniformRand();
#ifdef DebugLog
      edm::LogInfo("HFShower") << "HFShower::getHits: " << i << " attenuation "
                   << r1 << ":" << exp(-p*zFibre) << " r2 " << r2 
                   << ":" << probMax << " Survive: " 
                   << (r1 <= exp(-p*zFibre) && r2 <= probMax);
#endif
      if (applyFidCut || chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
    hit.depth      = depth;
    hit.time       = tSlice+time;
    if (!applyFidCut) {
      hit.wavelength = wavelength[i]; // Tmp
      hit.momentum   = momz[i]; // Tmp
    } else {
      hit.wavelength = 300.; // Tmp
      hit.momentum   = 1.; // Tmp
    }
    hit.position   = globalPos;
    hits.push_back(hit);
    nHit++;
      }
    }
  }

#ifdef DebugLog
  edm::LogInfo("HFShower") << "HFShower::getHits: Number of Hits " << nHit;
  for (int i=0; i<nHit; ++i)
    edm::LogInfo("HFShower") << "HFShower::Hit " << i << " WaveLength " 
                 << hits[i].wavelength << " Time " << hits[i].time
                 << " Momentum " << hits[i].momentum <<" Position "
                 << hits[i].position;
#endif
  return hits;

} 

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

private

Definition at line 659 of file HFShower.cc.

References debug, LogDebug, findQualityFiles::size, and relval_parameters_module::step.

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") << " FamosHFShower::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)
      LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter
       << " curr, F[curr-1], F[curr] = "
       << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;

  }

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


  return curr-1;
}

void HFShower::initRun

(

G4ParticleTable *

)

Definition at line 478 of file HFShower.cc.

Referenced by HCalSD::initRun().

                                        {// Define PDG codes
}

void HFShower::makeSteps ( int  nsteps )

private

Definition at line 366 of file HFShower.cc.

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

Referenced by HFShower().

                                   {

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


  if(debug)
    LogDebug("FastCalorimetry") << " FamosHFShower::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)
      LogDebug("FastCalorimetry") << " FamosHFShower::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.) {
      LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " - negative step energy !!!"
       << std::endl;
      est = 0.;
      break ;
    }

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

    if(debug == 2)
      LogDebug("FastCalorimetry") << " FamosHFShower::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)
      LogDebug("FastCalorimetry") << "*** FamosHFShower::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
  double etot = 0.;
  for (int i = 0; i < count ; i++) {
    eStep.push_back(temp[i] * e / sumes );
    nspots.push_back((int)(eStep[i]/eSpotSize)+1);
    etot += eStep[i];

   if(debug)
     LogDebug("FastCalorimetry") << i << "  xO and lamdepth at the end of step = "
      << x0depth[i] << " "
      << lamdepth[i] << "   Estep func = " <<  eStep[i]
      << "   Rstep = " << rlamStep[i] << "  Nspots = " <<  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)
      LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - " << "ECAL energy = " << eStep[0]
       << " out of total = " << e << std::endl;
  }

}

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

inlineprivate

Definition at line 54 of file HFShower.h.

References dttmaxenums::R.

Referenced by compute().

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

Member Data Documentation

double HFShower::aloge

private

Definition at line 77 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpEM

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpHD

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

bool HFShower::applyFidCut

private

Definition at line 49 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::balanceEH

private

Definition at line 117 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betEM

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betHD

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

HFCherenkov* HFShower::cherenkov

private

Definition at line 53 of file HFShower.h.

Referenced by getHits(), and HFShower().

int HFShower::chkFibre

private

Definition at line 56 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::criticalEnergy

private

Definition at line 113 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStart

private

Definition at line 76 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStep

private

Definition at line 111 of file HFShower.h.

Referenced by HFShower().

std::vector<int> HFShower::detector

private

Definition at line 79 of file HFShower.h.

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

double HFShower::e

private

Definition at line 96 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::eSpotSize

private

Definition at line 109 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::eStep

private

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

HFFibre* HFShower::fibre

private

Definition at line 54 of file HFShower.h.

Referenced by getHits(), and HFShower().

std::vector<double> HFShower::gpar

private

Definition at line 58 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::hcalDepthFactor

private

Definition at line 119 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::infinity

private

Definition at line 84 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::lambdaEM

private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::lambdaHD

private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

std::vector<double> HFShower::lamcurr

private

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::lamdepth

private

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::lamstep

private

Definition at line 82 of file HFShower.h.

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

std::vector<double> HFShower::lamtotal

private

Definition at line 82 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::lossesOpt

private

Definition at line 99 of file HFShower.h.

Referenced by compute(), and HFShower().

double HFShower::maxTRfactor

private

Definition at line 115 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::nDepthSteps

private

Definition at line 101 of file HFShower.h.

Referenced by HFShower().

std::vector<int> HFShower::nspots

private

Definition at line 79 of file HFShower.h.

Referenced by compute(), and makeSteps().

int HFShower::nTRsteps

private

Definition at line 103 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::onEcal

private

Definition at line 93 of file HFShower.h.

Referenced by HFShower().

double HFShower::part

private

Definition at line 72 of file HFShower.h.

double HFShower::probMax

private

Definition at line 57 of file HFShower.h.

Referenced by getHits(), and HFShower().

const RandomEngineAndDistribution* HFShower::random

private

Definition at line 122 of file HFShower.h.

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

std::vector<double> HFShower::rlamStep

private

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::tgamEM

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::tgamHD

private

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

const ECALProperties* HFShower::theECALproperties

private

Definition at line 67 of file HFShower.h.

Referenced by HFShower().

EcalHitMaker* HFShower::theGrid

private

Definition at line 87 of file HFShower.h.

Referenced by compute(), and HFShower().

HcalHitMaker* HFShower::theHcalHitMaker

private

Definition at line 90 of file HFShower.h.

Referenced by compute().

const HCALProperties* HFShower::theHCALproperties

private

Definition at line 68 of file HFShower.h.

Referenced by HFShower().

HDShowerParametrization* HFShower::theParam

private

Definition at line 64 of file HFShower.h.

Referenced by HFShower().

double HFShower::theR1

private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR2

private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR3

private

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transFactor

private

Definition at line 107 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transParam

private

Definition at line 105 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::x0curr

private

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::x0depth

private

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::x0EM

private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::x0HD

private

Definition at line 75 of file HFShower.h.

Referenced by HFShower().