00001 #include "FastShowerUtils/Normalisers/HadEarlyHcalCoreFracEc.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 HadEarlyHcalCoreFracEc::value(const ParticleParameters& pp) const{
00014
00015 double energy = pp.energy();
00016
00017 double p1 = -0.02474*std::exp(0.05275*energy) + 0.5719;
00018 double p2 = 0.5922*std::exp(-0.1043*energy) + 0.00576;
00019 double p3 = -0.6738*std::exp(-0.03873*energy) + 1.083;
00020 double p4 = -0.9529*std::exp(-0.04834*energy) + 1.240;
00021 double p5 = 0.5962*std::exp(-0.06673*energy) + 0.01706;
00022 double p6 = -3.278*std::exp(-0.01745*energy) + 4.014;
00023 double p7 = 0.2431*std::exp(-0.1905*energy) + 0.9941;
00024 double p8 = 0.6580*std::exp(-0.1478*energy) + 0.9878;
00025 double p9 = 0.8339*std::exp(-0.1028*energy) + 0.2240;
00026
00027 double x1o3 = 0.25;
00028 double x2 = 0.75;
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* HadEarlyHcalCoreFracEc::clone() const {
00063 return new HadEarlyHcalCoreFracEc(*this);
00064 }
00065 }
00066
00067
00068
00069