4 #include "KVLightEnergyCsIFull.h"
17 if (beta > beta_delta) rap = TMath::Power(beta / beta_delta, 2.);
19 double deno = log(1.022e6 / (16.*TMath::Power(fZmed, 0.9)) * TMath::Power(beta_delta, 2)) + log(rap);
20 frac = 0.5 * log(rap) / deno;
38 double e_delta = par[3];
41 if (zfit == 2. && afit == 4.) m = 3727.4;
42 if ((zfit == 1.) && (afit == 1.)) m = 938.;
50 beta = sqrt(1. - m * m / ((e + m) * (e + m)));
52 fraction = GetDeltaFraction(beta, beta_delta);
53 if (fraction > 1.) fraction = 1.;
54 if (fraction < 0.) fraction = 0.;
57 double se = sp_e(zfit, afit, e);
58 double sn = sp_n(zfit, afit, e);
59 double deno = 1. + ar * se + an * sn;
60 double deno2 = 1. + (1. - fraction) * ar * se + an * sn;
61 double raps = 1. + sn / se;
62 if (e <= e_delta * afit) {
63 if (se > 0. && raps > 0. && deno > 0.) dlum = (1. + an * sn) / (deno * raps);
66 if (se > 0. && raps > 0. && deno > 0.) {
67 dlum = (1. + an * sn) * (1. - fraction) / (deno2 * raps) + fraction / raps;
86 if (ar < 1.e-8) ar = 0.;
88 double e_delta = par[3];
90 if (zfit == 1. && afit == 1.) m = 938.2;
91 if (zfit == 1. && afit == 2.) m = 1876.;
92 if (zfit == 1. && afit == 3.) m = 2809.4;
93 if (zfit == 2. && afit == 3.) m = 2809.4;
94 if (zfit == 2. && afit == 4.) m = 3727.4;
95 if (zfit == 2. && afit == 6.) m = 5606.6;
96 if (zfit == 2. && afit == 8.) m = 7483.7;
101 beta = sqrt(1. - m * m / ((e + m) * (e + m)));
102 fraction = GetDeltaFraction(beta, beta_delta);
103 if (fraction > 1.) fraction = 1.;
104 if (fraction < 0.) fraction = 0.;
106 double se = sp_e(zfit, afit, e);
107 double sn = sp_n(zfit, afit, e);
109 if (se > 1.e-6) raps = 1. + sn / se;
110 double arg = 1. - ar * se / (1. + an * sn + ar * se);
112 double deno = 1. + an * sn + (1. - fraction) * ar * se;
114 if (e <= e_delta * afit) {
115 if (arg > 0. && raps * ar * se > 0.) dlum = -fact * log(arg) / (raps * ar * se);
120 if (deno > 0.) arg = 1. - (1. - fraction) * ar * se / deno;
121 if (arg > 0. && ar * se > 0. && raps > 0.) dlum = -fact * log(arg) / (raps * ar * se) + fraction / raps;
155 if (e < 1.e-4)
return 0.;
158 double a2_i = 0.6519;
163 double a7_i = 0.1043;
164 double a8_i = -0.1315;
165 double a9_i = -0.009175;
168 double a2_cs = 0.3923;
169 double a3_cs = 152.5;
170 double a4_cs = 8.354;
171 double a5_cs = 2.597;
172 double a6_cs = 4.671;
173 double a7_cs = 0.1136;
174 double a8_cs = -0.1298;
175 double a9_cs = -0.009078;
179 double ei = float(aref) * e / a * 1000.;
181 double shigh_cs = 0.;
187 slow_cs = a1_cs * pow(ei, a2_cs);
188 shigh_cs = 1000.*a3_cs * log(1. + 1000.*a4_cs / ei + a5_cs * ei / 1000.) / ei;
189 slow_i = a1_i * pow(ei, a2_i);
190 shigh_i = 1000.*a3_i * log(1. + 1000.*a4_i / ei + a5_i * ei / 1000.) / ei;
193 se_cs = fact * slow_cs * shigh_cs / (slow_cs + shigh_cs);
194 se_i = fact * slow_i * shigh_i / (slow_i + shigh_i);
197 se_cs = exp(a6_cs + a7_cs * log(1000. / ei) + a8_cs * pow(log(1000. / ei), 2.) + a9_cs * pow(log(1000. / ei), 3.));
198 se_i = exp(a6_i + a7_i * log(1000. / ei) + a8_i * pow(log(1000. / ei), 2.) + a9_i * pow(log(1000. / ei), 3.));
201 double seref = (127.*se_i + 133.*se_cs) / (127. + 133.);
203 se = pow(gamma(z, a, e) *
float(z) / (gamma(zref, aref, ei) *
float(zref)), 2.) * seref;
205 if (z ==
int(zref)) se = seref;
206 if (z == 1.) se = seref / (pow(
float(zref), 2.) * gamma(zref, aref, ei));
208 se *= 6.022 / (10.*fAmed);
227 g0 = (gamma_hbg(z, a, e) + gamma_ziegler(z, a, e)) / 2.;
231 if (e < 0.4 * a) g0 = gamma_ziegler(z, a, 0.4 * a);
243 if (z == 2. && a == 4.) m = 3727.4;
244 double v0 = 7.2983824e-3;
245 double v = sqrt(2.*ei / m);
246 double v1 = 0.886 * (v / v0) * pow(z, -2. / 3.);
247 double v2 = v1 + 0.038 * sin(3.141592654 * v1 / 2.);
248 double xgamma = 1. - (1.034 - 0.17777 * exp(-0.08114 * z)) * exp(-v2);
261 double d = 1.164 + 0.2319 * exp(-0.0004302 * z2);
262 double x1 = d + b * exp(-c * z);
263 double x2 = 8.144 + 0.09876 * log(z2);
264 double x3 = 0.314 + 0.01072 * log(z2);
265 double x4 = 0.5218 + 0.02521 * log(z2);
271 x1 = d + b * exp(-c * z);
281 x1 = d + b * exp(-c * z);
286 double xgamma = 1. - x1 * exp(-x2 * pow(e / a, x3) / pow(z, x4));
287 if (xgamma > 1.) xgamma = 1.;
288 if (xgamma < 0.) xgamma = 0.;
300 if (e < 1.e-4)
return 0.;
301 double ered = 32.53 * fAmed * e * 1000. / (z * fZmed * (a + fAmed) * sqrt(TMath::Power(z, 2. / 3.) + TMath::Power(fZmed, 2. / 3.)));
302 double sn = 1.593 * sqrt(ered);
303 if (ered >= 0.01 && ered <= 10.) sn = 1.7 * sqrt(ered) * (log(ered + exp(1.)) / (1. + 6.8 * ered + 3.4 * TMath::Power(ered, 1.5)));
304 if (ered > 10.) sn = log(0.47 * ered) / (2.*ered);
308 sn = sn * 8.462 * z * fZmed * a / ((a + fAmed) * sqrt(TMath::Power(z, 2. / 3.) + TMath::Power(fZmed, 2. / 3.)));
312 sn = sn * 6.022 / (10.*fAmed);
346 #if ROOT_VERSION_CODE > ROOT_VERSION(5,99,01)
348 fDlight->SetParameters(par);
349 return par[4] + fDlight->Integral(emin, emax, epsilon);
351 return par[4] + fDlight->Integral(emin, emax, par, epsilon);
379 Double_t energie = x[0];
381 Double_t c1 = par[1];
382 Double_t c2 = Z * Z * A * par[2];
383 Double_t c3 = A * par[3];
384 Double_t c4 = par[4];
386 Double_t c4_new = c4 / (1. + TMath::Exp((c3 - energie) / T));
387 Double_t c0 = c4 / (1. + TMath::Exp(c3 / T));
389 Double_t lumcalc = c1 * energie;
391 Double_t xm = -c1 * c0 * c2 * TMath::Log(c2 / (c2 + c3));
392 lumcalc = lumcalc - c1 * c2 * TMath::Log(1. + energie / c2) + c1 * c2 * c4_new * TMath::Log((energie + c2) / (c3 + c2)) + xm;
395 return par[0] + lumcalc;
423 Double_t pied = par[0];
424 Double_t a0 = par[1];
425 Double_t a1 = par[2];
426 Double_t a2 = par[3];
427 Double_t a3 = par[4];
429 Double_t raz2 = A * Z * Z * a1;
430 Double_t si = a0 * (E - raz2 * log(1. + E / raz2) + raz2 * log((E + raz2) /
431 (a3 * A + raz2)) - (1. - a2) * (1. - a2) * raz2 * log((E + (1. - a2) * raz2) /
432 (a3 * A + (1. - a2) * raz2)));
524 TString formula_types[] = {
"Exact",
"ApproxIntegral",
"Approx",
"ApproxSilicon"};
525 for (
int i = 0; i < 4; ++i) {
526 if (formula_types[i] == which) {
531 Error(
"SetLightFormula",
"Unkown light formula type : %s", which.Data());
virtual const Char_t * GetType() const
virtual void SetType(const Char_t *str)
void SetCalibFunction(TF1 *f)
virtual void Print(Option_t *opt="") const
Print a description of the calibration object, including a list of its parameters.
Exact calibration formula for CsI detectors.
TF1 * fDlight
function to integrate to get fLight
Double_t gamma_hbg(double z, double a, double e)
Double_t gamma(double z, double a, double e)
Hubert, Bimbot and Gauvin.
Double_t fZmed
Z of detector material.
void SetOptions(const KVNameValueList &opt)
Double_t dLightIntegralApprox(double *x, double *par)
Double_t GetLight(double *x, double *par)
Double_t dLightIntegral(double *x, double *par)
Double_t GetLightApproxSilicon(double *x, double *par)
Double_t GetLightApprox(double *x, double *par)
Double_t sp_n(double z, double a, double e)
return 0.;
Int_t fLightFormula
light formula (see NIMa of Marian)
Double_t fAmed
A of detector material (CsI)
Double_t sp_e(double z, double a, double e)
if energy E<0.1 keV, return 0
virtual void Print(Option_t *opt="") const
void SetLightFormula(Int_t form)
void init()
default initialisations
Double_t gamma_ziegler(double z, double a, double e)
Light-energy calibration function for CsI detectors using a Fermi-function dependence on energy for d...
Handles lists of named parameters with different types, a list of KVNamedParameter objects.
const Char_t * GetStringValue(const Char_t *name) const