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 RandomEngine *engine, HDShowerParametrization *myParam, EcalHitMaker *myGrid, HcalHitMaker *myHcalHitMaker, int onECAL, double epart)
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	depthStart
double	depthStep
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 RandomEngine *	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

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

Definition at line 29 of file HDShower.h.

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

Definition at line 30 of file HDShower.h.

typedef Steps::const_iterator HDShower::step_iterator

Definition at line 32 of file HDShower.h.

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

Definition at line 31 of file HDShower.h.

typedef math::XYZVector HDShower::XYZPoint

Definition at line 27 of file HDShower.h.

Constructor & Destructor Documentation

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

Definition at line 30 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, depthStart, depthStep, detector, e, HDShowerParametrization::e1(), HDShowerParametrization::e2(), EcalHitMaker::ecalHcalGapTotalL0(), EcalHitMaker::ecalHcalGapTotalX0(), HDShowerParametrization::ecalProperties(), EcalHitMaker::ecalTotalL0(), HDShowerParametrization::emax(), HDShowerParametrization::emid(), HDShowerParametrization::emin(), eSpotSize, RandomEngine::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, 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),
    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 )
    LogDebug("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)
    LogDebug("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 )
    LogDebug("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)
    LogDebug("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;
  
  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") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl; 
    depthStart = 0.;
  }
 
  if(depthStart > maxDepth) {
    if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " 
           << depthStart << std::endl; 
    
    depthStart = maxDepth *  random->flatShoot();
    if( depthStart < 0.) depthStart = 0.;
    if(debug) LogDebug("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) 
    LogDebug("FastCalorimetry") << " FamosHDShower : small energy, "
         << " depthStart reduced to = " << depthStart << std::endl; 
      
    }
  }

  if( depthHCAL < depthStep) {
    if(debug) LogDebug("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)
    LogDebug("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) LogDebug("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
    if(depthStart < depthECAL) {
      if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
      if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.1) {
    if(debug) LogDebug("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) LogDebug("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;
      }
      // Just shift starting point to HCAL
      else { 
    //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
    if(debug)  LogDebug("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) LogDebug("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
    }
  }
  else {   // Forward 
    if(debug)  LogDebug("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); }

}

virtual HDShower::~HDShower ( )

inlinevirtual

Definition at line 43 of file HDShower.h.

{;}

Member Function Documentation

bool HDShower::compute

(

)

Compute the shower longitudinal and lateral development.

Definition at line 471 of file HDShower.cc.

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), prof2calltree::count, gather_cfg::cout, debug, detector, patCandidatesForDimuonsSequences_cff::ecal, eStep, RandomEngine::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, getHLTprescales::index, indexFinder(), EcalCondDB::inf, infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, 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("FamosHDShower::compute");

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

      int ecal = 0;
      if(detector[i] != 1) { 
    bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
    
    if(debug)
      LogDebug("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);        
      }
      else {
    ecal = 1;
    bool status = theGrid->getPads(currentDepthL0);   
    
    if(debug)
      LogDebug("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);
      }  

      
      // 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)
      LogDebug("FastCalorimetry") << std::endl << " FamosHDShower::compute " << " r = " << radius 
           << "    phi = " << phi << std::endl;
    
    bool result;
    if(ecal) {
      result = theGrid->addHit(radius,phi,0);
      
      if(debug == 2)
        LogDebug("FastCalorimetry") << " FamosHDShower::compute - " 
         << " theGrid->addHit result = " 
         << result << std::endl;
    }
    else {
      result = theHcalHitMaker->addHit(radius,phi,0); 
      
      if(debug == 2)
        LogDebug("FastCalorimetry") << " FamosHDShower::compute - " 
         << " theHcalHitMaker->addHit result = " 
         << result << std::endl;
    }    
    if(result) nok ++; 
    
      } // end of tranverse simulation
      if(count == infinity) { 
    if(debug)
      LogDebug("FastCalorimetry") 
        << " FamosHDShower::compute " << " maximum number of" 
        << " transverse points " << count << " is used !!!" << std::endl; 
      }

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

    } // end of longitudinal steps
  } // end of no steps

  return status;

}

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

inlineprivate

Definition at line 51 of file HDShower.h.

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

Referenced by makeSteps().

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

int HDShower::getmip ( )

inline

Definition at line 41 of file HDShower.h.

References mip.

Referenced by CalorimetryManager::HDShowerSimulation().

{return mip;}

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

private

Definition at line 624 of file HDShower.cc.

References debug, LogDebug, findQualityFiles::size, and launcher::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") << " 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)
      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 HDShower::makeSteps ( int  nsteps )

private

Definition at line 359 of file HDShower.cc.

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

Referenced by HDShower().

                                   {

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


  if(debug)
    LogDebug("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)
      LogDebug("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.) {
      LogDebug("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)
      LogDebug("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)
      LogDebug("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)
     LogDebug("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) 
      LogDebug("FastCalorimetry") << " FamosHDShower::makeSteps - " << "ECAL energy = " << eStep[0]
       << " out of total = " << e << std::endl;  
  }

}

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

inlineprivate

Definition at line 56 of file HDShower.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 HDShower::aloge

private

Definition at line 79 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpEM

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpHD

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::balanceEH

private

Definition at line 122 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betEM

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betHD

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::criticalEnergy

private

Definition at line 118 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStart

private

Definition at line 78 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStep

private

Definition at line 116 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::detector

private

Definition at line 81 of file HDShower.h.

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

double HDShower::e

private

Definition at line 101 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::eSpotSize

private

Definition at line 114 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::eStep

private

Definition at line 82 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::hcalDepthFactor

private

Definition at line 124 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::infinity

private

Definition at line 86 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::lambdaEM

private

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

double HDShower::lambdaHD

private

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

std::vector<double> HDShower::lamcurr

private

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamdepth

private

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamstep

private

Definition at line 84 of file HDShower.h.

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

std::vector<double> HDShower::lamtotal

private

Definition at line 84 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::lossesOpt

private

Definition at line 104 of file HDShower.h.

Referenced by compute(), and HDShower().

double HDShower::maxTRfactor

private

Definition at line 120 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::mip

private

Definition at line 98 of file HDShower.h.

Referenced by getmip(), and HDShower().

int HDShower::nDepthSteps

private

Definition at line 106 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::nspots

private

Definition at line 81 of file HDShower.h.

Referenced by compute(), and makeSteps().

int HDShower::nTRsteps

private

Definition at line 108 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::onEcal

private

Definition at line 95 of file HDShower.h.

Referenced by HDShower().

double HDShower::part

private

Definition at line 74 of file HDShower.h.

const RandomEngine* HDShower::random

private

Definition at line 127 of file HDShower.h.

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

std::vector<double> HDShower::rlamStep

private

Definition at line 82 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::tgamEM

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::tgamHD

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

const ECALProperties* HDShower::theECALproperties

private

Definition at line 69 of file HDShower.h.

Referenced by HDShower().

EcalHitMaker* HDShower::theGrid

private

Definition at line 89 of file HDShower.h.

Referenced by compute(), and HDShower().

HcalHitMaker* HDShower::theHcalHitMaker

private

Definition at line 92 of file HDShower.h.

Referenced by compute().

const HCALProperties* HDShower::theHCALproperties

private

Definition at line 70 of file HDShower.h.

Referenced by HDShower().

HDShowerParametrization* HDShower::theParam

private

Definition at line 66 of file HDShower.h.

Referenced by HDShower().

double HDShower::theR1

private

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR2

private

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR3

private

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transFactor

private

Definition at line 112 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transParam

private

Definition at line 110 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0curr

private

Definition at line 83 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0depth

private

Definition at line 83 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::x0EM

private

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

double HDShower::x0HD

private

Definition at line 77 of file HDShower.h.

Referenced by HDShower().