00001 #include "FastShowerUtils/Normalisers/HadEarlyEcalCoreFracEc.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 HadEarlyEcalCoreFracEc::value(const ParticleParameters& pp) const{
00014
00015 double energy = pp.energy();
00016
00017 double p1 = 0.1561 - 0.001613*energy;
00018 double p2 = 0.1257*std::exp(-0.08885*energy) - 0.001882;
00019 double p3 = -0.8500*std::exp(-0.02366*energy) + 1.286;
00020 double p4 = 1.381 - 0.01140*energy;
00021 double p5 = 0.1839*std::exp(-0.04761*energy) + 0.01200;
00022 double p6 = 5.636;
00023 double p7 = 0.08336*std::exp(-0.12240*energy) + 1.0340;
00024 double p8 = 0.09318*std::exp(-0.09204*energy) + 1.0260;
00025 double p9 = 0.15480*std::exp(-0.06403*energy) + 0.1302;
00026
00027 double yref = RandFlat::shoot();
00028
00029 double x1o3 = 0.4;
00030 double yref1 = p4*std::pow(x1o3,p6) + p5;
00031 double yref3 = p1*std::pow(x1o3,p3) + p2;
00032
00033 double x2 = 0.8;
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
00057 result = max(result, 0.0);
00058 result = min(result, 1.0);
00059
00060 return result;
00061 }
00062
00063 IFnOfParticleParameters* HadEarlyEcalCoreFracEc::clone() const {
00064 return new HadEarlyEcalCoreFracEc(*this);
00065 }
00066 }
