---

HadHcalEtaResponseBar.cxx

Go to the documentation of this file.

#include "FastShowerUtils/Normalisers/HadHcalEtaResponseBar.h"

#include "FastShowerUtils/ParticleParameters.h"

#include "CLHEP/Random/RandFlat.h"

#include <cmath>

namespace FastShower{

  double HadHcalEtaResponseBar::value(const ParticleParameters& pp) const{
  
    // response symmetric about eta=0.0
    double aEta = fabs(pp.eta());
    
    // spline knots
    double x[14]={0.05, 0.15, 0.35, 0.55, 0.75, 0.95, 1.25, 
                  1.45, 1.75, 2.05, 2.25, 2.65, 2.95, 3.15 };

    // 2nd index is flag: (0)->e.m.+had.  &   (1)->e.m.
    double a[2][14] = {{.000, .839, .837, .860, .847, .846, .793,
                        .853, .815, .839, .839, .857, .823, .775},
                       {.000, .341, .337, .390, .442, .444, .515,
                        .543, .440, .489, .431, .433, .368, .412}};
    double b[2][14] = {{.000,-.0508, .0520, .0288,-.0113,-.184,-.0645,
                        .00211,-.140, .0839, .0523, .0354,-.0434,-.545},
                       {.000,-.0950, .0604, .359, .0759, .151, .413,
                        -.413,-.258, .0672,-.105,-.0408, .0417,-.0468}};
    double c[2][14] = {{.000, .000, 1.03,-1.14, .944,-1.81, 2.40,
                        -2.18, 1.47,-.722, .616,-.700, .503,-2.18},
                       {.000, .000, 1.55,-.0633,-1.35, 1.73, -.420,
                        -2.33, 3.11,-2.03, 1.45,-1.13, 1.34,-1.63}};
    double d[2][14] = {{.000, 3.43,-3.62, 3.48,-4.59, 7.02,-5.09,
                        6.08,-2.43, 1.49,-2.19, 1.00,-2.98, 3.63},
                       {.000, 5.18,-2.70,-2.15, 5.13,-3.58,-2.13,
                        9.07,-5.71, 3.87,-4.31, 2.06,-3.30, 2.72}};


    // iEta according to the abscissa of the spline function
    int iEta = (int)(aEta*10) + 1;
    double xval = (double)(iEta)*0.1 - 0.05;

    // default indices
    int i = 10;
    int j = 11;

    if (iEta>25) {
      // e.m.+had. response saturates at |eta|=2.5
      aEta=2.5;
    } else {
      // everything ends-up here if not saturated
      i=-1;
      do {i+=1;}
      while (i<13 && xval>x[i+1]);
    }
    
    // response at |eta|
    // e.m.+had. response at |eta|
    double emPLUShad = a[0][j] + b[0][j] * (aEta-x[i])
                               + c[0][j] * std::pow(aEta-x[i],2)
                               + d[0][j] * std::pow(aEta-x[i],3);
    // e.m. response at |eta|
    double em = a[1][j] + b[1][j] * (aEta-x[i])
                        + c[1][j] * std::pow(aEta-x[i],2)
                        + d[1][j] * std::pow(aEta-x[i],3);
    // had. response at |eta|
    double result = emPLUShad - em;

    // normalisation
    // e.m. normalisation at |eta|=0.2
    double emNorm = a[1][13] + b[1][13] * (0.2-x[1])
                             + c[1][13] * std::pow(0.2-x[1],2)
                             + d[1][13] * std::pow(0.2-x[1],3);

    // e.m.+had. normalisation at |eta|=0.2
    double emhadNorm = a[0][j] + b[0][j] * (0.2-x[1])
                               + c[0][j] * std::pow(0.2-x[1],2)
                               + d[0][j] * std::pow(0.2-x[1],3);
    // had. normalisation at |eta|=0.2
    double hadNorm = emhadNorm - emNorm;


    // relative response factors wrt had. response at |eta|=0.2
    return result/hadNorm;
  }
  //
  IFnOfParticleParameters* HadHcalEtaResponseBar::clone() const {
    return new HadHcalEtaResponseBar(*this);
  }
}//namespace

---