214 JParser<> zap(
"Program to determine the energy loss due to visible delta-rays.");
218 zap[
'N'] =
make_field(use_numerical,
"perform numeric integration");
220 zap[
'L'] =
make_field(lepton) = muon, tau, positron, electron;
221 zap[
'A'] =
make_field(A,
"atomic mass") = 18.0;
222 zap[
'Z'] =
make_field(Z,
"atomic number") = 10.0;
229 catch(
const exception &error) {
230 FATAL(error.what() << endl);
233 if (option.find(
'R') == string::npos) { option +=
'R'; }
234 if (option.find(
'S') == string::npos) { option +=
'S'; }
241 else if (lepton == tau)
243 else if (lepton == positron)
245 else if (lepton == electron)
248 FATAL(
"Invalid lepton " << lepton << endl);
250 const double Tmin = (T_GeV.
is_valid() ?
255 const double error = 1.0e-12;
257 const double xmin = log10(MASS_LEPTON);
258 const double xmax = 9.25;
260 TH1D h0(
"h0", NULL, 320, xmin, xmax);
261 TH1D h1(
"h1", NULL, 320, xmin, xmax);
263 for (
int i = 1; i <= h0.GetNbinsX(); ++i) {
265 const double x = h0.GetBinCenter(i);
266 const double E =
pow(10.0,x);
268 double Tmax = (lepton != electron ?
279 const double gamma = E / MASS_LEPTON;
280 const double beta = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
282 const JFormFactor F(0.5/(E*E), -beta*beta/Tmax, 1.0);
296 DEBUG(
"dE/dx [MeV g^-1 cm^2] " <<
FIXED(5,4) << dEdx << endl);
297 DEBUG(
"dE/dx [g^-1 cm^2] " <<
FIXED(5,2) << dNdx << endl);
299 h0.SetBinContent(i, dEdx);
300 h0.SetBinError (i, error);
301 h1.SetBinContent(i, dNdx);
306 TF1 f1(
"f1",
"[0] + [1]*x + [2]*x*x + [3]*x*x*x");
308 if (option.find(
'W') == string::npos) {
312 f1.SetParameter(0, h0.GetMinimum());
313 f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
314 f1.SetParameter(2, 0.0);
315 f1.SetParameter(3, 0.0);
317 f1.SetLineColor(kRed);
321 const TFitResultPtr
result = h0.Fit(&f1, option.c_str(),
"same", xmin, xmax);
323 cout <<
"chi2/NDF " <<
result->Chi2() <<
'/' <<
result->Ndf() << endl;
325 cout <<
"\t// " << lepton << endl;
326 cout <<
"\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
328 for (
int i = 0; i != 4; ++i) {
329 cout <<
"\tstatic const double " << (char) (
'a' + i) <<
" = " <<
SCIENTIFIC(10,3) << f1.GetParameter(i) <<
";" << endl;
332 double Emin = 1 * MASS_LEPTON;
333 double Emax = 5 * MASS_LEPTON;
335 for (
double E = 0.5 * (Emin + Emax); ; E = 0.5 * (Emin + Emax)) {
339 if (fabs(Tmax - Tmin) < precision) {
340 cout <<
"\tstatic const double Emin = " <<
FIXED(7,5) << E <<
"; // [GeV]" << endl;
double getDeltaRays(const double E, const double M, const double Tmin, const double Tmax, const double Z, const double A)
Get equivalent EM-shower energy due to delta-rays per unit track length for an ionising particle with...