#include <HFShower.h>

Classes
struct	Hit
Public Types
typedef std::pair< XYZPoint, double >	Spot
typedef std::pair< unsigned int, double >	Step
typedef Steps::const_iterator	step_iterator
typedef std::vector< Step >	Steps
typedef math::XYZVector	XYZPoint
Public Member Functions
bool	compute ()
	Compute the shower longitudinal and lateral development.
std::vector< Hit >	getHits (G4Step *aStep, bool forLibrary)
std::vector< Hit >	getHits (G4Step *aStep)
	HFShower (std::string &name, const DDCompactView &cpv, edm::ParameterSet const &p, int chk=0)
	HFShower (const RandomEngine engine, HDShowerParametrization myParam, EcalHitMaker myGrid, HcalHitMaker myHcalHitMaker, int onECAL, double epart)
virtual	~HFShower ()
virtual	~HFShower ()
Private Member Functions
double	gam (double x, double a) const
std::vector< double >	getDDDArray (const std::string &, const DDsvalues_type &, int &)
int	indexFinder (double x, const std::vector< double > &Fhist)
void	makeSteps (int nsteps)
double	transProb (double factor, double R, double r)
Private Attributes
double	aloge
double	alpEM
double	alpHD
double	balanceEH
double	betEM
double	betHD
HFCherenkov *	cherenkov
int	chkFibre
double	criticalEnergy
double	depthStart
double	depthStep
std::vector< int >	detector
double	e
double	eSpotSize
std::vector< double >	eStep
HFFibre *	fibre
std::vector< double >	gpar
double	hcalDepthFactor
int	infinity
double	lambdaEM
double	lambdaHD
std::vector< double >	lamcurr
std::vector< double >	lamdepth
std::vector< double >	lamstep
std::vector< double >	lamtotal
int	lossesOpt
double	maxTRfactor
int	nDepthSteps
std::vector< int >	nspots
int	nTRsteps
int	onEcal
double	part
double	probMax
const 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 HFShower.h.

Member Typedef Documentation

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

Definition at line 29 of file HFShower.h.

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

Definition at line 30 of file HFShower.h.

typedef Steps::const_iterator HFShower::step_iterator

Definition at line 32 of file HFShower.h.

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

Definition at line 31 of file HFShower.h.

typedef math::XYZVector HFShower::XYZPoint

Definition at line 27 of file HFShower.h.

Constructor & Destructor Documentation

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

Definition at line 30 of file HFShower.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 ad hoc long. extension + some fluctuations
  double   depthExt;
  if (e < 50.)   depthExt =  0.8 * (50.  - e) / 50.  + 0.3; 
  else { 
    if (e < 500.)  depthExt =  (500. - e) / 500. * 0.4 - 0.1; 
    else   depthExt =  -0.1; 
  }
  hcalDepthFactor += depthExt + 0.05 * (2.* random->flatShoot() - 1.);

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

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

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

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

  if( e < emid ) {
    theParam->setCase(1);
    // avoid "underflow" below Emin (for parameters calculation only)
    if(e < emin) effective = emin; 
  } 
  else 
    theParam->setCase(2);
 
  // A bit coarse espot size for HF...
  eSpotSize *= 2.5; 
  if(effective > 0.5 * emax) {
    eSpotSize *= 2.;
    if(effective > emax) {
      effective = emax; 
      eSpotSize *= 3.; 
      depthStep *= 2.;
     }
  }

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

  //---------------------------------------------------------------------------
  // Depth simulation & various protections, among them
  // if too deep - get flat random in the allowed region
  // if no HCAL material behind - force to deposit in ECAL
  double maxDepth    = depthToHCAL + depthHCAL - 1.1 * depthStep;

  double depthStart  = std::log(1./random->flatShoot()); // starting point lambda unts

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

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




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

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

        depthStart = depthToHCAL;
        sum1 += depthStart;     
      }
    }
    else { // GAP or HCAL   
      
      if(depthStart >= depthECAL &&  depthStart < depthToHCAL ) {
      depthStart = depthToHCAL; // just a shift to HCAL for simplicity
      }
      sum1 += depthStart;

      if(debug) LogDebug("FastCalorimetry") << " FamosHFShower : goto HCAL" << std::endl;
    }
  }
  else {   // Forward 
    if(debug)  LogDebug("FastCalorimetry") << " FamosHFShower : forward" << std::endl;
    sum1 += depthStart;
  }
 
  for (int i = 0; i < nmoresteps ; i++) {
    sum1 += depthStep;
    if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; 
    sum2 += depthStep;
    sum3 += sum2 * lambdaHD / x0HD;
    lamtotal.push_back(sum1);
    lamdepth.push_back(sum2);
    lamcurr.push_back(lambdaHD);
    lamstep.push_back(depthStep);
    x0depth.push_back(sum3);
    x0curr.push_back(x0HD);
    detector.push_back(3);
    nsteps++;
  }
  

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

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

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

}

HFShower::~HFShower ( ) [inline, virtual]

Definition at line 41 of file HFShower.h.

{;}

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

Definition at line 22 of file HFShower.cc.

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

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

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

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

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

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

  cherenkov = new HFCherenkov(m_HF);
  fibre     = new HFFibre(name, cpv, p);
}

virtual HFShower::~HFShower ( ) [virtual]

Member Function Documentation

bool HFShower::compute ( )

Compute the shower longitudinal and lateral development.

Definition at line 475 of file HFShower.cc.

References EcalHitMaker::addHit(), HcalHitMaker::addHit(), prof2calltree::count, 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("FamosHFShower::compute");

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

  if(numLongit > 0) {

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

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

      if( espot > 4. || espot < 0. ) 
        LogDebug("FastCalorimetry") << " FamosHFShower::compute - unphysical espot = " 
             << espot << std::endl;

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

        ecal = 1;
        bool status = theGrid->getPads(currentDepthL0);   
        
        if(debug)
          LogDebug("FastCalorimetry") << " FamosHFShower::compute - status of Grid = " 
               << status << std::endl;
        
        if(!status) continue; 

        theGrid->setSpotEnergy(espot);
      }  


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

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

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

        ++count;

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

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

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

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

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

        // process transverse distribution
        
        double prob   = random->flatShoot();
        int index     = indexFinder(prob,Fhist);
        double radius = rlamStep[i] * rhist[index] +
          random->flatShoot() * rbinsize; // in-bin  
        double phi = 2.*M_PI*random->flatShoot();
        
        if(debug == 2)
          LogDebug("FastCalorimetry") << std::endl << " FamosHFShower::compute " << " r = " << radius 
               << "    phi = " << phi << std::endl;

        // add hit
        
        theHcalHitMaker->setSpotEnergy(newespot);
        theGrid->setSpotEnergy(newespot);
        
        bool result;
        if(ecal) {                                 
          result = theGrid->addHit(radius,phi,0);  // shouldn't get here !
          
          if(debug == 2)
            LogDebug("FastCalorimetry") << " FamosHFShower::compute - " 
                 << " theGrid->addHit result = " 
                 << result << std::endl;

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

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

        }    

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

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

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

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

}

double HFShower::gam	(	double	x,
		double	a
	)		const `[inline, private]`

Definition at line 49 of file HFShower.h.

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

Referenced by makeSteps().

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

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

Definition at line 238 of file HFShower.cc.

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

Referenced by HFShower().

                                                      {

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

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

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

Definition at line 164 of file HFShower.cc.

References abs, cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, AlCaRecoCosmics_cfg::name, HFShower::Hit::position, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, and HFShower::Hit::wavelength.

                                                                        {

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

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

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

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

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

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

  G4StepPoint *     preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector     globalPos    = preStepPoint->GetPosition();
  G4String          name         = preStepPoint->GetTouchable()->GetSolid(0)->GetName();
  double            zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
  G4ThreeVector     localPos     = G4ThreeVector(globalPos.x(),globalPos.y(),
                                                 zv);
  G4ThreeVector     localMom     = preStepPoint->GetTouchable()->GetHistory()->
    GetTopTransform().TransformAxis(momentumDir);
  int               depth        = (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10;
  G4ThreeVector     translation  = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();

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

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

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

  return hits;

}

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

Definition at line 69 of file HFShower.cc.

References abs, HFFibre::attLength(), cherenkov, chkFibre, HFCherenkov::computeNPE(), relval_parameters_module::energy, funct::exp(), fibre, HFCherenkov::getMom(), HFCherenkov::getWL(), gpar, i, LogDebug, HFShower::Hit::momentum, AlCaRecoCosmics_cfg::name, L1TEmulatorMonitor_cff::p, HFShower::Hit::position, probMax, HFShower::Hit::time, cond::rpcobgas::time, HFFibre::tShift(), v, and HFShower::Hit::wavelength.

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

                                                       {

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

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

  G4Track *aTrack = aStep->GetTrack();
  const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle();
 
  HFShower::Hit hit;
  double energy   = aParticle->GetTotalEnergy();
  double momentum = aParticle->GetTotalMomentum();
  double pBeta    = momentum / energy;
  double dose     = 0.;
  int    npeDose  = 0;

  G4ThreeVector momentumDir = aParticle->GetMomentumDirection();
  G4ParticleDefinition *particleDef = aTrack->GetDefinition();
     
  G4StepPoint *     preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector     globalPos    = preStepPoint->GetPosition();
  G4String          name         = 
    preStepPoint->GetTouchable()->GetSolid(0)->GetName();
  double            zv = std::abs(globalPos.z()) - gpar[4] - 0.5*gpar[1];
  G4ThreeVector     localPos     = G4ThreeVector(globalPos.x(),globalPos.y(),
                                                 zv);
  G4ThreeVector     localMom     = preStepPoint->GetTouchable()->GetHistory()->
    GetTopTransform().TransformAxis(momentumDir);
  int               depth        = 
    (preStepPoint->GetTouchable()->GetReplicaNumber(0))%10;
  G4ThreeVector     translation  = preStepPoint->GetTouchable()->GetVolume(1)->GetObjectTranslation();

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

#ifdef DebugLog
  LogDebug("HFShower") << "HFShower::getHits: in " << name << " Z " << zCoor 
                       << "(" << globalPos.z() << ") " << zFibre << " Trans "
                       << translation << " Time " << tSlice  << " " << time 
                       << "\n                  Direction " << momentumDir 
                       << " Local " << localMom;
#endif 
  int npe = cherenkov->computeNPE(particleDef, pBeta, u, v, w, stepl, zFibre, 
                                  dose, npeDose);
  std::vector<double> wavelength = cherenkov->getWL();
  std::vector<double> momz       = cherenkov->getMom();
  
  for (int i = 0; i<npe; ++i) {
    double p   = fibre->attLength(wavelength[i]);
    double r1  = G4UniformRand();
    double r2  = G4UniformRand();
#ifdef DebugLog
    LogDebug("HFShower") << "HFShower::getHits: " << i << " attenuation " << r1
                         <<":" << exp(-p*zFibre) << " r2 " << r2 << ":" 
                         << probMax << " Survive: " 
                         << (r1 <= exp(-p*zFibre) && r2 <= probMax);
#endif
    if (chkFibre < 0 || (r1 <= exp(-p*zFibre) && r2 <= probMax)) {
      hit.time = tSlice+time;
      hit.wavelength = wavelength[i];
      hit.momentum   = momz[i];
      hit.position   = globalPos;
      hits.push_back(hit);
      nHit++;
    }
  }

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

}

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

Definition at line 659 of file HFShower.cc.

References debug, LogDebug, findQualityFiles::size, and 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") << " FamosHFShower::indexFinder - wrong current index = " 
           << curr << " !!!" << std::endl;

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

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


  return curr-1;
}

void HFShower::makeSteps ( int nsteps ) [private]

Definition at line 366 of file HFShower.cc.

References aloge, alpEM, alpHD, balanceEH, betEM, betHD, prof2calltree::count, debug, detector, e, eSpotSize, eStep, 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 HFShower().

                                   {

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


  if(debug)
    LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps - " 
       << " nsteps required : " << nsteps << std::endl;

  int count = 0;
  for (int i = 0; i < nsteps; i++) {    
    
    double deplam = lamdepth[i] - 0.5 * lamstep[i];
    double depx0  = x0depth[i]  - 0.5 * lamstep[i] / x0curr[i]; 
    double     x = betEM * depx0;
    double     y = betHD * deplam;
    
    if(debug == 2)
      LogDebug("FastCalorimetry") << " FamosHFShower::makeSteps " 
                                  << " - step " << i
                                  << "   depx0, x = " << depx0 << ", " << x 
                                  << "   deplam, y = " << deplam << ", "
                                  << y << std::endl;

    double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
                  (x0curr[i] * tgamEM) + 
                  (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
    
    // protection ...
    if(est < 0.) {
      LogDebug("FastCalorimetry") << "*** FamosHFShower::makeSteps " << " - negative step energy !!!" 
           << std::endl;
      est = 0.;
      break ; 
    }

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

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

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

    // good step - to proceed

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

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

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

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

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

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

  }

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


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

}

double HFShower::transProb	(	double	factor,
		double	R,
		double	r
	)		`[inline, private]`

Definition at line 54 of file HFShower.h.

References dttmaxenums::R.

Referenced by compute().

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

Member Data Documentation

double HFShower::aloge [private]

Definition at line 77 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpEM [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::alpHD [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::balanceEH [private]

Definition at line 117 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betEM [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::betHD [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

HFCherenkov* HFShower::cherenkov [private]

Definition at line 49 of file HFShower.h.

Referenced by getHits(), and HFShower().

int HFShower::chkFibre [private]

Definition at line 52 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::criticalEnergy [private]

Definition at line 113 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStart [private]

Definition at line 76 of file HFShower.h.

Referenced by HFShower().

double HFShower::depthStep [private]

Definition at line 111 of file HFShower.h.

Referenced by HFShower().

std::vector<int> HFShower::detector [private]

Definition at line 79 of file HFShower.h.

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

double HFShower::e [private]

Definition at line 96 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::eSpotSize [private]

Definition at line 109 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::eStep [private]

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

HFFibre* HFShower::fibre [private]

Definition at line 50 of file HFShower.h.

Referenced by getHits(), and HFShower().

std::vector<double> HFShower::gpar [private]

Definition at line 54 of file HFShower.h.

Referenced by getHits(), and HFShower().

double HFShower::hcalDepthFactor [private]

Definition at line 119 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::infinity [private]

Definition at line 84 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::lambdaEM [private]

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::lambdaHD [private]

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

std::vector<double> HFShower::lamcurr [private]

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::lamdepth [private]

Definition at line 82 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::lamstep [private]

Definition at line 82 of file HFShower.h.

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

std::vector<double> HFShower::lamtotal [private]

Definition at line 82 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::lossesOpt [private]

Definition at line 99 of file HFShower.h.

Referenced by compute(), and HFShower().

double HFShower::maxTRfactor [private]

Definition at line 115 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::nDepthSteps [private]

Definition at line 101 of file HFShower.h.

Referenced by HFShower().

std::vector<int> HFShower::nspots [private]

Definition at line 79 of file HFShower.h.

Referenced by compute(), and makeSteps().

int HFShower::nTRsteps [private]

Definition at line 103 of file HFShower.h.

Referenced by compute(), and HFShower().

int HFShower::onEcal [private]

Definition at line 93 of file HFShower.h.

Referenced by HFShower().

double HFShower::part [private]

Definition at line 72 of file HFShower.h.

double HFShower::probMax [private]

Definition at line 53 of file HFShower.h.

Referenced by getHits(), and HFShower().

const RandomEngine* HFShower::random [private]

Definition at line 122 of file HFShower.h.

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

std::vector<double> HFShower::rlamStep [private]

Definition at line 80 of file HFShower.h.

Referenced by compute(), and makeSteps().

double HFShower::tgamEM [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::tgamHD [private]

Definition at line 72 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

const ECALProperties* HFShower::theECALproperties [private]

Definition at line 67 of file HFShower.h.

Referenced by HFShower().

EcalHitMaker* HFShower::theGrid [private]

Definition at line 87 of file HFShower.h.

Referenced by compute(), and HFShower().

HcalHitMaker* HFShower::theHcalHitMaker [private]

Definition at line 90 of file HFShower.h.

Referenced by compute().

const HCALProperties* HFShower::theHCALproperties [private]

Definition at line 68 of file HFShower.h.

Referenced by HFShower().

HDShowerParametrization* HFShower::theParam [private]

Definition at line 64 of file HFShower.h.

Referenced by HFShower().

double HFShower::theR1 [private]

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR2 [private]

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::theR3 [private]

Definition at line 71 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transFactor [private]

Definition at line 107 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::transParam [private]

Definition at line 105 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::x0curr [private]

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

std::vector<double> HFShower::x0depth [private]

Definition at line 81 of file HFShower.h.

Referenced by HFShower(), and makeSteps().

double HFShower::x0EM [private]

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

double HFShower::x0HD [private]

Definition at line 75 of file HFShower.h.

Referenced by HFShower().

HFShower Class Reference

Classes

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation