HDShower Class Reference

#include <HDShower.h>

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

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

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

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

Detailed Description

Definition at line 22 of file HDShower.h.

Member Typedef Documentation

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

Definition at line 21 of file HDShower.cc.

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

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

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

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

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

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

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

  // Get the Famos Histos pointer
  //  myHistos = FamosHistos::instance();
  //  std::cout << " Hello FamosShower " << std::endl;
  
  theECALproperties = theParam->ecalProperties();
  theHCALproperties = theParam->hcalProperties();


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

  if( e < emid ) {
    theParam->setCase(1);
    // avoid "underflow" below Emin (for parameters calculation only)
    if(e < emin) effective = emin; 
  } 
  else 
    theParam->setCase(2);
 
  // Avoid "overflow" beyond Emax (for parameters)
  if(effective > 0.5 * emax) {
    eSpotSize *= 2.5;
    if(effective > emax) {
      effective = emax; 
      eSpotSize *= 4.; 
      depthStep *= 2.;
      if(effective > 1000.)
      eSpotSize *= 2.; 
    }
  }

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


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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

  if(e < criticalEnergy ) nmoresteps = 1;  
  else                    nmoresteps = nDepthSteps;
  
  depthECAL  = 0.;
  depthGAP   = 0.;
  depthGAPx0 = 0.; 
  if(onECAL ) {
    depthECAL  = theGrid->ecalTotalL0();         // ECAL depth segment
    //depthECAL  = theGrid->ecalTotalL0() + theGrid->ps1TotalL0() + theGrid->ps2TotalL0() + theGrid->ps2eeTotalL0(); //TEST: include preshower
    depthGAP   = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
    depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
  }
  
  depthHCAL   = theGrid->hcalTotalL0();    // HCAL depth segment
  depthToHCAL = depthECAL + depthGAP;    

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


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

  if(depthHCAL < depthStep) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " << depthHCAL << " depthStart -> forced to 0 !!!" << std::endl;
    depthStart = 0.;    
    nmoresteps = 0;
    
    if(depthECAL < depthStep) {
      nsteps = -1;
      LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - " << " particle is lost !!! " << std::endl; 
    }
  }

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

  if( onEcal ) {
    if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
    if(depthStart < depthECAL) {
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
      if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.1) {
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step" << std::endl;
        
        //  ECAL - one step
        nsteps++; 
        sum1   += depthECAL;             // at the end of step
        sum2   += depthECAL-depthStart; 
        sum3   += sum2 * lambdaEM / x0EM;
        lamtotal.push_back(sum1);
        lamdepth.push_back(sum2);
        lamcurr.push_back(lambdaEM);
        lamstep.push_back(depthECAL-depthStart);
        x0depth.push_back(sum3);
        x0curr.push_back(x0EM);
        detector.push_back(1);
        mip = 0;    
        
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : " << " in ECAL sum1, sum2 " << sum1 << " " << sum2 << std::endl;
        
        // Gap - no additional step after ECAL
        // just move further to HCAL over the gap
        sum1   += depthGAP;          
        sum2   += depthGAP; 
        sum3   += depthGAPx0;
      }
      else { // Just shift starting point to HCAL
        //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
        if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;

        depthStart = depthToHCAL;
        sum1 += depthStart;     
      }
    }
    else { // GAP or HCAL         
      if(depthStart >= depthECAL &&  depthStart < depthToHCAL ) {
        depthStart = depthToHCAL; // just a shift to HCAL for simplicity
      }
      sum1 += depthStart;
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
    }
  }
  else {   // Forward 
    if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
    sum1 += depthStart;
    transFactor = 0.5;   // makes narower tresverse size of shower     
  }
 
  for (int i = 0; i < nmoresteps ; i++) {
    sum1 += depthStep;
    if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; 
    sum2 += depthStep;
    sum3 += sum2 * lambdaHD / x0HD;
    lamtotal.push_back(sum1);
    lamdepth.push_back(sum2);
    lamcurr.push_back(lambdaHD);
    lamstep.push_back(depthStep);
    x0depth.push_back(sum3);
    x0curr.push_back(x0HD);
    detector.push_back(3);
    nsteps++;
  }

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

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

}

virtual HDShower::~HDShower ( )

inlinevirtual

Definition at line 44 of file HDShower.h.

References compute().

{;}

Member Function Documentation

bool HDShower::compute

(

)

Compute the shower longitudinal and lateral development.

Definition at line 447 of file HDShower.cc.

References HcalHitMaker::addHit(), EcalHitMaker::addHit(), KineDebug3::count(), gather_cfg::cout, debug, detector, digitizers_cfi::ecal, eStep, RandomEngineAndDistribution::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, mps_fire::i, diffTreeTool::index, indexFinder(), spu::inf(), infinity, lamstep, lamtotal, lossesOpt, M_PI, maxTRfactor, nspots, nTRsteps, phi, TtFullHadEvtBuilder_cfi::prob, TCMET_cfi::radius, random, mps_fire::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), mps_update::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation(), and ~HDShower().

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

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

  if(numLongit > 0) {

    status = true;    
    // Prepare the trsanverse probability function
    std::vector<double> Fhist;
    std::vector<double> rhist; 
    for (int j = 0; j < nTRsteps + 1; j++) {
      rhist.push_back(maxTRfactor * j / nTRsteps );  
      Fhist.push_back(transProb(maxTRfactor,1.,rhist[j]));
      if(debug == 3) LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
    }
    
    // Longitudinal steps
    for (int i = 0; i < numLongit ; i++) {
      
      double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
      // vary the longitudinal profile if needed
      if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor;                     
      if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i] << "  currentDepthL0 = " << currentDepthL0 << std::endl;
      
      double maxTRsize   = maxTRfactor * rlamStep[i];     // in lambda units
      double rbinsize    = maxTRsize / nTRsteps; 
      double espot       = eStep[i] / (double)nspots[i];  // re-adjust espot

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

      int ecal = 0;
      if(detector[i] != 1) { 
        bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
    
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of " << " theHcalHitMaker->setDepth(currentDepthL0) is " 
                    << setHDdepth << std::endl;
    
        if(!setHDdepth) {
          currentDepthL0 -= lamstep[i];
          setHDdepth =  theHcalHitMaker->setDepth(currentDepthL0);
        }
        
        if(!setHDdepth) continue;

        theHcalHitMaker->setSpotEnergy(espot);
        
        //fill hcal longitudinal distribution histogram
      }
      else {
        ecal = 1;
        bool status = theGrid->getPads(currentDepthL0);   
    
        if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;
    
        if(!status) continue; 

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

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

      for (int j = 0; j < inf; j++) {
        if(nok == nspots[i]) break;
        count ++;
        
        double prob   = random->flatShoot();
        int index     = indexFinder(prob,Fhist);
        double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize; // in-bin  
        double phi = 2.*M_PI*random->flatShoot();
    
        if(debug == 2)  LogInfo("FastCalorimetry") << std::endl << " FamosHDShower::compute " << " r = " << radius 
                        << "    phi = " << phi << std::endl;
    
        bool result;
        if(ecal) {
          result = theGrid->addHit(radius,phi,0);
      
          if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theGrid->addHit result = " << result << std::endl;
          
          //fill ecal transverse distribution histogram
        }
        else {
          result = theHcalHitMaker->addHit(radius,phi,0); 
 
          if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theHcalHitMaker->addHit result = " << result << std::endl;
          
          //fill hcal transverse distribution histogram
        }    
        
        if(result) nok ++; 



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

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

  return status;

}

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

inlineprivate

Definition at line 52 of file HDShower.h.

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

Referenced by makeSteps().

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

int HDShower::getmip ( )

inline

Definition at line 42 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 582 of file HDShower.cc.

References debug, and findQualityFiles::size.

Referenced by compute(), and transProb().

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

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

  for (iter = 0; iter < size ; iter++) {

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

    if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break;
    prevdir = actudir;
    if(x > Fhist[curr]) {actudir =  1;}
    else                {actudir = -1;}
    if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
    curr += actudir * step;
    if(curr > size) curr = size;
    else { if(curr < 1) {curr = 1;}}
 
    if(debug == 3)
      LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter 
       << " curr, F[curr-1], F[curr] = "
       << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;
    
  }

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


  return curr-1;
}

void HDShower::makeSteps ( int  nsteps )

private

Definition at line 339 of file HDShower.cc.

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

Referenced by HDShower(), and transProb().

                                   {

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

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

  int count = 0;
  for (int i = 0; i < nsteps; i++) {    
    
    double deplam = lamdepth[i] - 0.5 * lamstep[i];
    double depx0  = x0depth[i]  - 0.5 * lamstep[i] / x0curr[i]; 
    double     x = betEM * depx0;
    double     y = betHD * deplam;
    
    if(debug == 2)
      LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps " 
                                  << " - step " << i
                  << "   depx0, x = " << depx0 << ", " << x 
                  << "   deplam, y = " << deplam << ", "
                  << y << std::endl;
    
    double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
          (x0curr[i] * tgamEM) + 
          (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
    
    // protection ...
    if(est < 0.) {
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " - negative step energy !!!" 
       << std::endl;
      est = 0.;
      break ; 
    }

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

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

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

    // good step - to proceed

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

  // fluctuations in ECAL and re-distribution of remaining energy in HCAL
  if(detector[0] == 1 && count > 1) {
    double oldECALenergy = temp[0];
    double oldHCALenergy = sumes - oldECALenergy ;
    double newECALenergy = 2. * sumes;
    for (int i = 0; newECALenergy > sumes && i < infinity; i++)
      newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy; 
     
    if(debug == 2)
      LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " ECAL fraction : old/new - "
       << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl;

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

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

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

  }

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

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

}

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

inlineprivate

Definition at line 57 of file HDShower.h.

References indexFinder(), makeSteps(), dttmaxenums::R, and x.

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 80 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpEM

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpHD

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::balanceEH

private

Definition at line 123 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betEM

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betHD

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::criticalEnergy

private

Definition at line 119 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthECAL

private

Definition at line 131 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthGAP

private

Definition at line 131 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthGAPx0

private

Definition at line 131 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthHCAL

private

Definition at line 131 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStart

private

Definition at line 79 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStep

private

Definition at line 117 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthToHCAL

private

Definition at line 131 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::detector

private

Definition at line 82 of file HDShower.h.

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

double HDShower::e

private

Definition at line 102 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::eSpotSize

private

Definition at line 115 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::eStep

private

Definition at line 83 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::hcalDepthFactor

private

Definition at line 125 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::infinity

private

Definition at line 87 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::lambdaEM

private

Definition at line 78 of file HDShower.h.

Referenced by HDShower().

double HDShower::lambdaHD

private

Definition at line 78 of file HDShower.h.

Referenced by HDShower().

std::vector<double> HDShower::lamcurr

private

Definition at line 85 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamdepth

private

Definition at line 85 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamstep

private

Definition at line 85 of file HDShower.h.

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

std::vector<double> HDShower::lamtotal

private

Definition at line 85 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::lossesOpt

private

Definition at line 105 of file HDShower.h.

Referenced by compute(), and HDShower().

double HDShower::maxTRfactor

private

Definition at line 121 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::mip

private

Definition at line 99 of file HDShower.h.

Referenced by getmip(), and HDShower().

int HDShower::nDepthSteps

private

Definition at line 107 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::nspots

private

Definition at line 82 of file HDShower.h.

Referenced by compute(), and makeSteps().

int HDShower::nTRsteps

private

Definition at line 109 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::onEcal

private

Definition at line 96 of file HDShower.h.

Referenced by HDShower().

double HDShower::part

private

Definition at line 75 of file HDShower.h.

const RandomEngineAndDistribution* HDShower::random

private

Definition at line 128 of file HDShower.h.

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

std::vector<double> HDShower::rlamStep

private

Definition at line 83 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::tgamEM

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::tgamHD

private

Definition at line 75 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

const ECALProperties* HDShower::theECALproperties

private

Definition at line 70 of file HDShower.h.

Referenced by HDShower().

EcalHitMaker* HDShower::theGrid

private

Definition at line 90 of file HDShower.h.

Referenced by compute(), and HDShower().

HcalHitMaker* HDShower::theHcalHitMaker

private

Definition at line 93 of file HDShower.h.

Referenced by compute().

const HCALProperties* HDShower::theHCALproperties

private

Definition at line 71 of file HDShower.h.

Referenced by HDShower().

HDShowerParametrization* HDShower::theParam

private

Definition at line 67 of file HDShower.h.

Referenced by HDShower().

double HDShower::theR1

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR2

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR3

private

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transFactor

private

Definition at line 113 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transParam

private

Definition at line 111 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0curr

private

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0depth

private

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::x0EM

private

Definition at line 78 of file HDShower.h.

Referenced by HDShower().

double HDShower::x0HD

private

Definition at line 78 of file HDShower.h.

Referenced by HDShower().