00001 #include "FastShowerUtils/Normalisers/HadEarlyHcalCoreFracBar.h"
00002
00003 #include "FastShowerUtils/ParticleParameters.h"
00004
00005 #include "CLHEP/Random/RandFlat.h"
00006
00007 #include <float.h>
00008 #include <algorithm>
00009 #include <cmath>
00010
00011 namespace FastShower{
00012
00013 double HadEarlyHcalCoreFracBar::value(const ParticleParameters& pp) const{
00014
00015 double energy = pp.energy();
00016
00017 double p1 = 0.6066*std::exp(-0.0433*energy) + 0.09131;
00018 double p2 = 0.2930*std::exp(-0.1202*energy) + 0.00100;
00019 double p3 = -2.505*std::exp(-0.01723*energy) + 2.8010;
00020 double p4 = -0.9982*std::exp(-0.037595*energy) + 1.3277;
00021 double p5 = 0.7788*std::exp(-0.1131*energy) + 0.01989;
00022 double p6 = -7.2420*std::exp(-0.01934*energy) + 8.045;
00023 double p7 = -0.002219*energy + 1.248;
00024 double p8 = 0.3260*std::exp(-0.02749*energy) + 0.9911;
00025 double p9 = 0.3969*std::exp(-0.03985*energy) + 0.1169;
00026
00027 double x1o3 = 0.4;
00028 double x2 = 0.8;
00029
00030 double yref = RandFlat::shoot();
00031
00032 double yref1 = p4*std::pow(x1o3,p6) + p5;
00033 double yref3 = p1*std::pow(x1o3,p3) + p2;
00034 double yref2 = p4*std::pow(x2,p6) + p5;
00035
00036 double result;
00037 if (yref>yref2) {
00038 result = p8 - p9*std::sqrt(fabs(-2.*std::log(yref/p7)));
00039 if ((yref >= 1.) && (result < 1.)) {
00040 result += RandFlat::shoot()*(1.-result);
00041 }
00042 }else if (yref>yref1) {
00043 if (yref<yref3) {
00044 if (RandFlat::shoot()<0.5) {
00045 result = std::pow(max((yref-p2),DBL_EPSILON)/p1,1/p3);
00046 }else{
00047 result = std::pow(max((yref-p5),DBL_EPSILON)/p4,1/p6);
00048 }
00049 }else{
00050 result = std::pow(max((yref-p5),DBL_EPSILON)/p4,1/p6);
00051 }
00052 }else{
00053 result = std::pow(max((yref-p2),DBL_EPSILON)/p1,1/p3);
00054 }
00055
00056 result = max(result, 0.0);
00057 result = min(result, 1.0);
00058
00059 return result;
00060 }
00061
00062 IFnOfParticleParameters* HadEarlyHcalCoreFracBar::clone() const {
00063 return new HadEarlyHcalCoreFracBar(*this);
00064 }
00065 }
