Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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
12namespace JSIRENE {}
13namespace JPP { using namespace JSIRENE; }
14
15namespace 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
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Detector simulations.
double weight_neutron(const double ekin)
Definition pythia.hh:149
double weight_klong(const double ekin)
Definition pythia.hh:119
double ngamma_elec(const double ekin)
Definition pythia.hh:58
double weight_kshort(const double ekin)
Definition pythia.hh:108
double opa_efrac(const int ipart, const double ekin)
Definition pythia.hh:173
double opa_weight_high_e(const double ekin)
Definition pythia.hh:20
static const JPythia pythia
Function object for relative light yield as a function of GEANT particle code.
Definition JPythia.hh:96
double weight_proton(const double ekin)
Definition pythia.hh:135
double weight_kaon(const double ekin)
Definition pythia.hh:92
double weight_pion(const double ekin)
Definition pythia.hh:69
Collection of particles.
Definition JPDB.hh:117
static const JPDB & getInstance()
Get particle data book.
Definition JPDB.hh:131
const JParticle & getPDG(const int pdg) const
Get particle corresponding to given PDG code.
Definition JPDB.hh:240
int geant
GEANT code of particle.
Definition JPDB.hh:96