---

HadEcalHcalSharesBase.cxx

Go to the documentation of this file.

#include "FastShowerUtils/Normalisers/HadEcalHcalSharesBase.h"

#include "FastShowerUtils/Normalisers/DistRandomiser2D.h"
#include "FastShowerUtils/ParticleParameters.h"
//#include "FastShowerUtils/IDoubleShowererSelectorConfig.h"
#include "FastShowerUtils/ITripleShowererSelectorConfig.h"

#include <float.h>
#include <cmath>
#include <iostream>
namespace FastShower{
  //
  //HadEcalHcalSharesBase::HadEcalHcalSharesBase(const IDoubleShowererSelectorConfig* c, 
  //                                             const std::string& s):
  HadEcalHcalSharesBase::HadEcalHcalSharesBase(const ITripleShowererSelectorConfig* c, 
                                               const std::string& s):
    IFnOfParticleParameters2(), DebugBase(s){

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

    // loop through files, and setup 2d-dist-randomiser
    cout<<"HadEcalHcalSharesBase: 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<<"HadEcalHcalSharesBase:  hist = \""
          <<(*itr).second<<"\",  energy = "<<(*itr).first<<"GeV"<<endl;
      // pass histogram-file-name to 2d-dist-randomiser
      DistRandomiser2D* randomiser = new DistRandomiser2D((*itr).second);
      // set the corresponding energy (energy point)
      m_energyDists[(*itr).first] = randomiser;
    }
    m_energyDists[FLT_MAX] = 0;  // null pointers for end-points
    cout<<"HadEcalHcalSharesBase: distributions done"<<endl;
  }
  //
  std::pair<double,double> 
  HadEcalHcalSharesBase::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);

    DistRandomiser2D* randomiser = (*mapEntryPoint).second;
    double energyPoint = (*mapEntryPoint).first;
    std::pair<double,double> ehFracs = randomiser->sample();

    nudge(energyPoint,energy,ehFracs);

    return ehFracs;
  }
  //
  void HadEcalHcalSharesBase::nudge(const double energyPoint, 
                                    const double energy,
                                    std::pair<double,double>& ehFracs) const {

    double energyPointFracSum = ehFracs.first + ehFracs.second;
    double energyFracSum = energyPointFracSum;
    
    if (energy<=0.15) {
    // no interpolation at very small Energies
    } else if (energy<=1.0) {
    // interpolate mean but not sigma
      energyFracSum = energyPointFracSum - mean(energyPoint) + mean(energy);
    } else {
    // interpolate both mean and sigma
      double sigmaScale = sigma(energy)/sigma(energyPoint);
      energyFracSum = (energyPointFracSum - mean(energyPoint))*sigmaScale 
                     + mean(energy);
    }
    ehFracs.first *= (energyFracSum/energyPointFracSum);
    ehFracs.second *= (energyFracSum/energyPointFracSum);
  }
  //
  HadEcalHcalSharesBase::~HadEcalHcalSharesBase(){
    std::map<double, DistRandomiser2D*>::iterator itr = m_energyDists.begin();
    std::map<double, DistRandomiser2D*>::iterator end = m_energyDists.end();
    for(; itr!= end; ++itr){
      delete (*itr).second;
    }
    m_energyDists.clear();
  }
}//namespace

---