---

SmearSoftPhotonBase.cxx

Go to the documentation of this file.

#include "FastShowerUtils/Normalisers/SmearSoftPhotonBase.h"

#include "FastShowerUtils/Normalisers/DistRandomiser1D.h"
#include "FastShowerUtils/ParticleParameters.h"
#include "FastShowerUtils/ISingleShowererSelectorConfig.h"

#include <float.h>
#include <cmath>
#include <iostream>
namespace FastShower{
  SmearSoftPhotonBase::SmearSoftPhotonBase(const ISingleShowererSelectorConfig* c, 
                                           const std::string& s): 
    IFnOfParticleParameters(), DebugBase(s){

    // get histograms' file-names and energies
    cout<<"SmearSoftPhotonBase: getting histograms ..."<<endl;
    m_histograms = c->histograms();

    // loop through files, and setup 1d-dist-randomiser
    cout<<"SmearSoftPhotonBase: setting up distributions ..."<<endl;
    MCItr_Type1 itr = m_histograms.begin();
    MCItr_Type1 end = m_histograms.end();
    m_energyDists[0.0] = 0;  // null pointers for end-points
    for (;itr!=end; ++itr){
      cout<<"SmearSoftPhotonBase:  hist = \""
          <<(*itr).second<<"\",  energy = "<<(*itr).first<<"GeV"<<endl;
      // pass histogram-file-name to 1d-dist-randomiser
      DistRandomiser1D* randomiser = new DistRandomiser1D((*itr).second);
      // set the corresponding energy (energy point)
      m_energyDists[(*itr).first] = randomiser;
    }
    m_energyDists[FLT_MAX] = 0;  // null pointers for end-points
    cout<<"SmearSoftPhotonBase: distributions done"<<endl;
  }
  //
  double SmearSoftPhotonBase::value(const ParticleParameters& pp) const{
    double energy = pp.rawEnergy();

    energy = max((float)energy,FLT_EPSILON);

    MCItr_Type2 floor = m_energyDists.begin();
    MCItr_Type2 roof  = m_energyDists.end();
    /* 
     * usual case: energy differs from key values in the map => lower==upper
     * rare case: energy overlaps with a key entry in the map => lower==upper-1
     * energies below/above the floor/roof  =>  valid iterators
     * energies overlapping with an energyPoint (i.e. a key) are special!
     */

    // the following logic works only if the map has at least 3 entries, 
    // i.e. at least 1 valid entry
    std::pair<MCItr_Type2,MCItr_Type2> luItrs = m_energyDists.equal_range(energy);
    MCItr_Type2 lowerBound = luItrs.first;
    MCItr_Type2 upperBound = luItrs.second;

    // force the lower/upperBound iterators to point at 2 consecutive entries
    if (lowerBound==upperBound) {--lowerBound;}  // takes care of the usual case
    // find the closest boundary
    MCItr_Type2 mapEntryPoint = ((energy - (*lowerBound).first)<
                         ((*upperBound).first - energy))? 
      lowerBound : upperBound;

    // move up to the 1st valid point if energy is between 0 and the 1st valid point
    if (mapEntryPoint == floor) {++mapEntryPoint;}

    // check for strange condition that energy is closer to FLT_MAX 
    // than to the last valid point
    MCItr_Type2 ceiling = --roof;
    assert(mapEntryPoint != ceiling);

    DistRandomiser1D* randomiser = (*mapEntryPoint).second;
    double energyPoint = (*mapEntryPoint).first;
    double eFrac = randomiser->sample();
    nudge(energyPoint,energy,eFrac);
    return eFrac;
  }
  //
  void 
  SmearSoftPhotonBase::nudge(const double energyPoint, 
                             const double energy, double& eFrac) const {
    // interpolation
    eFrac += (mean(energy) - mean(energyPoint));
  }
  //
  SmearSoftPhotonBase::~SmearSoftPhotonBase(){
    m_histograms.clear();
    std::map<double,DistRandomiser1D*>::iterator itr2 = m_energyDists.begin();
    std::map<double,DistRandomiser1D*>::iterator end2 = m_energyDists.end();
    for(; itr2!= end2; ++itr2){
      delete (*itr2).second;
    }
    m_energyDists.clear();
  }
}//namespace

---