---

C0EmEcalEc1.cxx

Go to the documentation of this file.

// ================================================
// Implementation of C0EmEcalEc1
// ================================================

#include "FastShowerUtils/Samplers/C0EmEcalEc1.h"

#include "FastShowerUtils/ParticleParameters.h"
#include "FastShowerUtils/PolyArgs.h"
#include "FastShowerUtils/IFn.h"
#include "FastShowerUtils/CoreSamples.h"
#include "FastShowerUtils/LinearProcessor.h"
#include "FastShowerUtils/UpdatingGaussian.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 C0EmEcalEc1::s_nSigma1(3.0);
  const double C0EmEcalEc1::s_nSigma2(5.0);
  //const double C0EmEcalEc1::s_tailFluctScale(0.025);
  const double C0EmEcalEc1::s_tailFluctScale(0.125);

  C0EmEcalEc1::C0EmEcalEc1(IUpdatingGaussian* g, 
                           IProcessedDist* t,
                           IInTail* inTail,
                           IFn* ts):
    ISampler(), ICell0(), DebugBase("C0EmEcalEc1"), 
    m_peak(g), m_tail(t), m_inTail(inTail), m_tailLowerBound(ts){}
  C0EmEcalEc1::C0EmEcalEc1(const IConfigurer* configurer,
                           const std::string&
                             ):
    DebugBase("C0"), m_peak(0), m_tail(0), m_inTail(0), m_tailLowerBound(0){
    m_peak           = configurer->findIUG( text()+"Peak" ) ;
    m_tail           = configurer->makeProcessedFlat( text() );
    m_inTail         = configurer->makeIInTail( text() );
    m_tailLowerBound = configurer->findFn( text()+"TailLowerBound");
    cout<<text()<<" finishing construction"<<endl;
  }
  ISampler* C0EmEcalEc1::clone() const {return new C0EmEcalEc1(*this);}

  void C0EmEcalEc1::sample(const PolyArgs& pa, CoreSamples& cs) const{

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

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

    double fracBelow01 = 0.0;

    // handle very low fraction deposits for low energies
    double energy = pa.pp()->energy();
    double aDelPhi = pa.pp()->delPhi();

    fracBelow01 = (0.7-0.03*energy)+(0.3*energy-3.0)*aDelPhi;
    fracBelow01 = max(fracBelow01, 0.0);

    SplitDecision below01(fracBelow01);
    if (below01.lower()) {
      double base = 0.1+RandFlat::shoot()*0.05;
      return (RandFlat::shoot())*base;
    } else {
      //calculate the upper and lower limits of the tails
      //
      // upper tail
      //double uLimit2 = peakMean + s_nSigma2*peakSigma;
      double uLimit2 = (aDelPhi<0.037)? 
        min(peakMean+s_nSigma2*peakSigma, 0.99) : 0.99;
      uLimit2 = min(uLimit2, 1.0);
      //
      double lLimit2 = peakMean + s_nSigma1*peakSigma;
      if(lLimit2>uLimit2){lLimit2 = uLimit2-peakSigma;}
      //
      // lower tail
      double uLimit1 = peakMean - s_nSigma1*peakSigma;
      //
      double lLimit1 = uLimit1 - ( m_tailLowerBound->value(pa) );
      //Add a little fluctuation on the lower tail
      lLimit1 -= s_tailFluctScale*( RandFlat::shoot() );
      if (fracBelow01>0.0) { lLimit1 = max(lLimit1, 0.1); } 
      
      lLimit1 = max(lLimit1,0.0);
      //
      if(uLimit1<lLimit1){ uLimit1=peakMean;}
      
      // Sample from the flat tail
      // decide which tail we are in?
      double base1 = uLimit1-lLimit1;
      double base2 = uLimit2-lLimit2;

      SplitDecision sp(0.422-0.002*energy);
      if ( sp.lower() ){
        LinearProcessor lp(base2,lLimit2);
        return m_tail->sample(&lp);
      }else{
        LinearProcessor lp(base1,lLimit1);
        return m_tail->sample(&lp);
      }
    }
  }
  double C0EmEcalEc1::evalPeak(const PolyArgs& pa)const{
    return m_peak->sample(0.0,s_nSigma1,1.0,s_nSigma1,pa);
  }
  double C0EmEcalEc1::lastValue(const CoreSamples& cs) const {
    return cs.give(this);
  }
  void C0EmEcalEc1::components(IDebug::Cpts& v) const{
    v.push_back(m_peak);
    v.push_back(m_tail);
    v.push_back(m_inTail);
    v.push_back(m_tailLowerBound);
  }
 
}//namespace

---