---

SNEmEcalEc1.cxx

Go to the documentation of this file.

// ================================================
// Implementation of SNEmEcalEc1
// ================================================

#include "FastShowerUtils/Samplers/SNEmEcalEc1.h"

#include "FastShowerUtils/PolyArgs.h"
#include "FastShowerUtils/IFn.h"
#include "FastShowerUtils/CoreSamples.h"
#include "FastShowerUtils/LinearProcessor.h"
#include "FastShowerUtils/IProcessedDist.h"
#include "FastShowerUtils/ProcessedNormal.h"
#include "FastShowerUtils/IInTail.h"
#include "FastShowerUtils/SplitDecision.h"
#include "FastShowerUtils/IConfigurer.h"
#include "CLHEP/Random/RandFlat.h"
#include <cmath>
#include <algorithm>
#include <iostream>

namespace FastShower{

  using std::pair;
  
  const double SNEmEcalEc1::s_nSigma(3.0);
  SNEmEcalEc1::SNEmEcalEc1(IUpdatingGaussian* g,
                           IProcessedDist* f,
                           IInTail* it):
    ISampler(), ICellSN(), DebugBase("SNEmEcalEc1"), 
    m_peak(g), m_tail(f), m_inTail(it){}
  SNEmEcalEc1::SNEmEcalEc1(const IConfigurer* configurer,
                           const std::string&):
    DebugBase("SN1"), m_peak(0), m_tail(0), m_inTail(0){
    m_peak           = configurer->findIUG( text()+"Peak" ) ;
    m_tail           = configurer->makeProcessedFlat( text() );
    m_inTail         = configurer->makeIInTail( text() );
    cout<<text()<<" finishing construction"<<endl;
  }
  ISampler* SNEmEcalEc1::clone() const {return new SNEmEcalEc1(*this);}

  void SNEmEcalEc1::sample(const PolyArgs& pa, CoreSamples& cs) const{
    
    if( m_inTail->operator()(pa) ){
      cs.fill(this, evalTail(pa) );
    }else{
      cs.fill(this, evalPeak(pa) );
    }
  }
  double SNEmEcalEc1::evalTail(const PolyArgs& pa) const{

    pair<double, double>  peakParams = m_peak->parameters(pa);
    double peakMean  = peakParams.first;
    double peakSigma = peakParams.second;

    //calculate the upper and lower limits of the tails
    double uLimit = peakMean - s_nSigma*peakSigma;
    double lLimit = 0.0;

    double energy = pa.pp()->energy();
    double cell0  = pa.cs()->cell0();

    //find a better name for "limit" (see TailFracSNEmEcalEc1.cxx)
    double limit = 0.0; 
    if (energy<15.0) {limit = 0.45;}
    else if (energy<25.0) {limit = 0.40;}

    if (cell0 < limit) {
      double fracBin1 = 0.0;
      fracBin1 = (energy<10.0)?  0.35/(0.35+0.25) :  0.0;
      fracBin1 = (energy>=10.0 && energy<25.0)?  1.09 - 0.0507*energy :  0.0;
      fracBin1 = max(fracBin1, 0.0);
      fracBin1 = min(fracBin1, 1.0);
      SplitDecision inBin1(fracBin1);
      if (inBin1.lower()) {
        lLimit = 0.0;
        uLimit = 0.05;
      } else {
        lLimit = 0.05;
      }
    } else {
      if (energy<15.0) {
        lLimit = 0.50*peakMean;
      } else if (energy<30.0) {
        lLimit = 0.75*peakMean;
      } else {
        lLimit = 0.95*peakMean;
      }
    }
    //
    if(lLimit<0.){lLimit=0.;}
    //IMPROVEME - the logic of setting limits to be made clearer
    if(uLimit<0.){ uLimit=peakMean;}
    if(uLimit<=lLimit){ uLimit=peakMean;}
    //
    LinearProcessor lp( (uLimit-lLimit), lLimit );
    return m_tail->sample(&lp);
  }
  double SNEmEcalEc1::evalPeak(const PolyArgs& pa) const{
    pair<double, double>  peakParams = m_peak->parameters(pa);
    double peakMean  = peakParams.first;
    if (peakMean<=0.0) return 0.0;
    return m_peak->sample(0.0, s_nSigma, peakMean, s_nSigma, pa);
  }
  double SNEmEcalEc1::lastValue(const CoreSamples& cs) const {
    return cs.give(this);
  }
  void SNEmEcalEc1::components(IDebug::Cpts& v) const{
    v.push_back(m_peak);
    v.push_back(m_tail);
    v.push_back(m_inTail);
  }
}//namespace

---