#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.
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.

  : 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 ( ) [inline, virtual]

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 EcalHitMaker::addHit(), HcalHitMaker::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 `[inline, private]`

Definition at line 51 of file HDShower.h.

References funct::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 ExpressReco_HICollisions_FallBack::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, cond::rpcobtemp::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, ExpressReco_HICollisions_FallBack::x, x0curr, x0depth, and ExpressReco_HICollisions_FallBack::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
	)		`[inline, private]`

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().

HDShower Class Reference

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation