Jpp  18.2.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pythia.hh
Go to the documentation of this file.
1 #ifndef __pythia__
2 #define __pythia__
3 
4 #include "JAAnet/JPDB.hh"
5 
6 /**
7  * \file
8  * This file contains converted Fortran code from km3.
9  * \author mdejong
10  */
11 
12 namespace JSIRENE {}
13 namespace JPP { using namespace JSIRENE; }
14 
15 namespace JSIRENE {
16 
17  using JAANET::JPDB;
18 
19  // *** high-energy fit to weights for a high-energy pion (same as v4r5 version)
20  inline double opa_weight_high_e(const double ekin)
21  {
22  static const double Ma = 72.425;
23  static const double Mb = -49.417;
24  static const double Mc = 5.858;
25  static const double Md = 207.252;
26  static const double Me = 132.784;
27  static const double Mf = -10.277;
28  static const double Mg = -19.441;
29  static const double Mh = 58.598;
30  static const double Mi = 53.161;
31  static const double Mkref = 2.698;
32  // leading_coef Mlc=(Ma-Mf)/Mkref
33  static const double Mlc = (Ma-Mf)/Mkref;
34 
35  double num, denom, lognp, Ee, logE;
36 
37  // implement the energy-fractions according to M. Dentler: ANTARES-SOFT-2012-009. E must be the kinetic energy in GeV. Here, 'everything else' is treated as a pion. Other than charged pions, mostly kaons of all varieties and protons fall into this category. These fits were made mostly at high energies - at low energies, other fits are used.
38 
39  if (ekin > 0.2)
40  logE=log10(ekin);
41  else
42  // nphotons constant below 0.2 GeV (lowest fit point)
43  return 0.292;
44 
45  num = Me + Md*logE + Mc*pow(logE,2) + Mb*pow(logE,3) + Ma *pow(logE,4) + Mlc*pow(logE,5);
46  denom = Mi + Mh*logE + Mg*pow(logE,2) + Mf*pow(logE,3) + Mlc*pow(logE,4);
47 
48  if (denom <= 0) {
49  return 0.;
50  }
51 
52  lognp = num/denom;
53  Ee = pow(10.0, lognp-Mkref);
54 
55  return Ee/ekin;
56  }
57 
58  inline double ngamma_elec(const double ekin)
59  {
60  static const double eleca = 1.33356e5;
61  static const double elecb = 1.66113e2;
62  static const double elecc = 16.4949;
63  static const double elecd = 1.5385e5;
64  static const double elece = 6.04871e5;
65 
66  return (eleca*ekin+elecb) * exp(-ekin/elecc) + (elecd*ekin+elece)* (1.-exp(-ekin/elecc));
67  }
68 
69  inline double weight_pion(const double ekin)
70  {
71  static const double pia = 0.538346;
72  static const double pib = 1.32041;
73  static const double pic = 0.737415;
74  static const double pid = -0.813861;
75  static const double pie = -2.22444;
76  static const double pif = -2.15795;
77  static const double pig = -6.47242;
78  static const double pih = -2.7567;
79  static const double pix = 8.83498;
80 
81  if (ekin < 0.06)
82  return 1e4*pia/ngamma_elec(ekin);
83  else if (ekin < 0.15)
84  return pix*1e4*ekin/ngamma_elec(ekin);
85  else {
86  const double lek=log(ekin);
87  return (pib*1e5*ekin + (pow(ekin,(1.-1./pic))) * (pid*1e4 + 1e4*pie*lek + 1e4*pif*pow(lek,2) + 1e3*pig*pow(lek,3) + 1e2*pih*pow(lek,4)))/ngamma_elec(ekin);
88  }
89  }
90 
91  // **** calculates the expected number of photons for pions as a function of energy
92  inline double weight_kaon(const double ekin)
93  {
94  static const double ka = 12.7537;
95  static const double kb = -1.052;
96  static const double kc = 3.49559;
97  static const double kd = 16.7914;
98  static const double ke = -0.355066;
99  static const double kf = 2.77116;
100 
101  if (ekin > 0.26)
102  return (1e4*ka*(ekin+kb)*(1.-exp(-ekin/kc)) + 1e4*(kd*ekin+ke)*exp(-ekin/kc))/ngamma_elec(ekin);
103  else
104  return kf*1e4/ngamma_elec(ekin);
105  }
106 
107  // **** calculates the expected number of photons for K0-short as a function of energy
108  inline double weight_kshort(const double ekin)
109  {
110  static const double k0sa = 0.351242;
111  static const double k0sb = 0.613076;
112  static const double k0sc = 6.24682;
113  static const double k0sd = 2.8858;
114 
115  return (k0sa*1e5 + k0sb*1e5*ekin + ekin*k0sc*1e4*log(k0sd + 1./ekin))/ngamma_elec(ekin);
116  }
117 
118  // **** calculates the expected number of photons for K0-short as a function of energy
119  inline double weight_klong(const double ekin)
120  {
121  static const double k0la = 2.18152;
122  static const double k0lb = -0.632798;
123  static const double k0lc = 0.999514;
124  static const double k0ld = 1.36052;
125  static const double k0le = 4.22484;
126  static const double k0lf = -0.307859;
127 
128  if (ekin < 1.5)
129  return (1e4*k0la + ekin*1e5*k0lb*log(ekin) + 1e5*k0lc*pow(ekin,1.5))/ngamma_elec(ekin);
130  else
131  return (ekin*k0ld*1e5 + pow(ekin,1.-1./k0le)*k0lf*1e5)/ngamma_elec(ekin);
132  }
133 
134  // **** calculates the expected number of photons for protons as a function of energy
135  inline double weight_proton(const double ekin)
136  {
137  static const double pa = 12.1281;
138  static const double pb = -17.1528;
139  static const double pc = 0.573158;
140  static const double pd = 34.1436;
141  static const double pe = -0.28944;
142  static const double pf = -13.2619;
143  static const double pg = 24.1357;
144 
145  return (1e4*(pa*ekin+pb)*(1.-exp(-ekin/pc)) + 1e4*(pd*ekin +pe+ pf*pow(ekin,2) + pg*pow(ekin,3))*exp(-ekin/pc)) / ngamma_elec(ekin);
146  }
147 
148  // **** calculates the expected number of photons for neutrons as a function of energy
149  inline double weight_neutron(const double ekin)
150  {
151  static const double na = 1.24605;
152  static const double nb = 0.63819;
153  static const double nc = -0.802822;
154  static const double nd = -0.935327;
155  static const double ne = -6.1126;
156  static const double nf = -1.96894;
157  static const double ng = 0.716954;
158  static const double nh = 2.68246;
159  static const double ni = -9.39464;
160  static const double nj = 15.0737;
161 
162  if (ekin > 0.5) {
163  const double lek=log(ekin);
164  return (na*1e5*ekin + 1e3*pow(ekin,1.-1./nb) * (100.*nc+100.*nd*lek + 10.*ne*pow(lek,2) + 11.*nf*pow(lek,3)))/ngamma_elec(ekin);
165  } else
166  return (1e3*ng + 1e4*nh*ekin + 1e4*ni*pow(ekin,2) + 1e4*nj*pow(ekin,3))/ngamma_elec(ekin);
167  }
168 
169 
170  // Routine to calculate the effective electron-shower energy using the
171  // one-particle approximation.
172  // includes low-energy fits per particle, and high-energy fits for pions
173  inline double opa_efrac(const int ipart, const double ekin)
174  {
175  using namespace JPP;
176 
177  double etemp, weight_part;
178 
179  // neutrino or muon
180  if (ipart == GEANT4_TYPE_NEUTRINO ||
181  ipart == GEANT4_TYPE_ANTIMUON ||
182  ipart == GEANT4_TYPE_MUON ||
183  ipart == GEANT4_TYPE_ANTITAU ||
184  ipart == GEANT4_TYPE_TAU) {
185  return 0.0;
186  }
187 
188  // 100% energy to EM showers if electron, positron, pi0, or gamma
189  if (ipart == GEANT4_TYPE_PHOTON ||
190  ipart == GEANT4_TYPE_ANTIELECTRON ||
191  ipart == GEANT4_TYPE_ELECTRON ||
192  ipart == GEANT4_TYPE_NEUTRAL_PION ||
193  ipart == GEANT4_TYPE_ETA) {
194  return 1.0;
195  }
196 
197  // low-energy fits went up to 40 GeV
198  if (ekin <= 40)
199  etemp=ekin;
200  else
201  etemp=40.;
202 
203  // otherwise, calculate fractions the long way. First we get the
204  // expected number of photons from the particle:
205  if (ipart == GEANT4_TYPE_PION_PLUS ||
206  ipart == GEANT4_TYPE_PION_MINUS)
207  weight_part=weight_pion(etemp);
208  else if (ipart == GEANT4_TYPE_KAON_LONG)
209  weight_part=weight_klong(etemp);
210  else if (ipart == GEANT4_TYPE_KAON_SHORT)
211  weight_part=weight_kshort(etemp);
212  else if (ipart == GEANT4_TYPE_KAON_PLUS ||
213  ipart == GEANT4_TYPE_KAON_MINUS)
214  weight_part=weight_kaon(etemp);
215  else if (ipart == GEANT4_TYPE_NEUTRON)
216  weight_part=weight_neutron(etemp);
217  else if (ipart == GEANT4_TYPE_PROTON ||
218  ipart == GEANT4_TYPE_ANTIPROTON)
219  weight_part=weight_proton(etemp);
220  else
221  // when in doubt, treat it as a proton
222  weight_part=weight_proton(etemp);
223 
224  if (ekin <= 40)
225  return weight_part;
226  else if (ekin >= 1e7)
227  return opa_weight_high_e(ekin);
228  else {
229  // smooth transition between low-energy weights and high-energy pion weight
230  const double wp40 = weight_part;
231  const double whe40 = opa_weight_high_e(40.);
232  const double whex = opa_weight_high_e(ekin);
233  return whex - (whe40-wp40)*(7.-log10(ekin))/5.398;
234  }
235  }
236 
237 
238  /**
239  * Get equivalent EM-energy for given pion energy.
240  *
241  * \param type particle type [PDG]
242  * \param E particle energy [GeV]
243  * \return EM-equivalent energy [GeV]
244  */
245  inline double pythia(const int type, const double E)
246  {
247  int ipart = JPDB::getInstance().getPDG(type).geant;
248  float ekin = (float) E;
249 
250  return opa_efrac(ipart, ekin) * E;
251  }
252 }
253 
254 #endif
JSIRENE::weight_klong
double weight_klong(const double ekin)
Definition: pythia.hh:119
JSIRENE::opa_weight_high_e
double opa_weight_high_e(const double ekin)
Definition: pythia.hh:20
E
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
JSIRENE::weight_proton
double weight_proton(const double ekin)
Definition: pythia.hh:135
JAANET::JPDB::getPDG
const JParticle & getPDG(const int pdg) const
Get particle corresponding to given PDG code.
Definition: JPDB.hh:234
JSIRENE::pythia
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition: JPythia.hh:96
JAANET::JPDB
Collection of particles.
Definition: JPDB.hh:109
JSIRENE::ngamma_elec
double ngamma_elec(const double ekin)
Definition: pythia.hh:58
JSIRENE::weight_pion
double weight_pion(const double ekin)
Definition: pythia.hh:69
JAANET::JPDB::getInstance
static const JPDB & getInstance()
Get particle data book.
Definition: JPDB.hh:125
JSIRENE::weight_kaon
double weight_kaon(const double ekin)
Definition: pythia.hh:92
log10
set_variable E_E log10(E_{fit}/E_{#mu})"
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
log
then set_variable DIR else fatal Wrong number of arguments fi for INPUT_FILE in ls rt $DIR stage * log
Definition: plot-acoustics-fit.sh:30
JAANET::JParticle::geant
int geant
GEANT code of particle.
Definition: JPDB.hh:91
JSIRENE::opa_efrac
double opa_efrac(const int ipart, const double ekin)
Definition: pythia.hh:173
JSIRENE::weight_kshort
double weight_kshort(const double ekin)
Definition: pythia.hh:108
JSIRENE::weight_neutron
double weight_neutron(const double ekin)
Definition: pythia.hh:149
exp
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
Definition: toashort_monitor.sh:46
Generated by   doxygen 1.8.5