22 static const double K = 0.307075;
23 static const double Z = 10.0;
24 static const double A = 18.0;
26 static const char*
const electron =
"electron";
27 static const char*
const positron =
"positron";
28 static const char*
const muon =
"muon";
29 static const char*
const tau =
"tau";
41 int main(
int argc,
char **argv)
57 JParser<> zap(
"Program to determine the energy loss due to visible delta-rays.");
69 catch(
const exception &error) {
70 FATAL(error.what() << endl);
79 else if (lepton == tau)
81 else if (lepton == positron)
86 FATAL(
"Invalid lepton " << lepton << endl);
100 NOTICE(
"Tmin [GeV] " <<
FIXED(8,6) << Tmin <<
' ' <<
FIXED(8,6) << T_GeV.getLowerLimit() << endl);
102 if (T_GeV.is_valid()) {
103 if (Tmin < T_GeV.getLowerLimit()) {
104 Tmin = T_GeV.getLowerLimit();
111 TH1D h0(
"h0", NULL, 80,
log10(MASS_LEPTON), +9.25);
112 TH1D h1(
"h1", NULL, 80,
log10(MASS_LEPTON), +9.25);
114 for (
int i = 1;
i <= h0.GetNbinsX(); ++
i) {
116 const double x = h0.GetBinCenter(
i);
117 const double E =
pow(10.0,x);
119 const double gamma = E / MASS_LEPTON;
120 const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
126 (0.5*sqrt(E*E - MASS_ELECTRON*MASS_ELECTRON)) :
127 (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
129 if (T_GeV.is_valid()) {
130 if (Tmax > T_GeV.getUpperLimit()) {
131 Tmax = T_GeV.getUpperLimit();
141 const double a = 1.0;
142 const double b = -beta*beta/Tmax;
143 const double c = 0.5/(E*
E);
145 const double W = 0.5 *
K * (
Z/
A) * (1.0/(beta*beta));
152 for (
double x = xmin; x <=
xmax; x += dx) {
154 const double T =
exp(x);
156 const double F = a + T*(b + T*(
c));
157 const double y = W * F * dx;
163 const double xmin = 1.0/Tmax;
164 const double xmax = 1.0/Tmin;
167 for (
double x = xmin; x <=
xmax; x += dx) {
169 const double T = 1.0/
x;
171 const double F = a + T*(b + T*(
c));
172 const double y = W * F * dx;
182 DEBUG(
"dE/dx [MeV g^-1 cm^2] " <<
FIXED(5,4) << weight << endl);
183 DEBUG(
"dE/dx [g^-1 cm^2] " <<
FIXED(5,2) << count << endl);
185 h0.SetBinContent(
i, weight);
186 h1.SetBinContent(
i, count);
190 TF1
f1(
"f1",
"[0] + [1]*x + [2]*x*x + [3]*x*x*x");
192 if (option.find(
'W') == string::npos) {
196 f1.SetParameter(0, h0.GetMinimum());
197 f1.SetParameter(1, (h0.GetMaximum() - h0.GetMinimum()) / (h0.GetXaxis()->GetXmax() - h0.GetXaxis()->GetXmin()));
198 f1.SetParameter(2, 0.0);
199 f1.SetParameter(3, 0.0);
201 f1.SetLineColor(kRed);
205 h0.Fit(&f1, option.c_str(),
"same");
208 cout <<
"\t// " << lepton << endl;
209 cout <<
"\t// dE/dX = a + bx + cx^2 + dx^3; // [MeV g^-1 cm^2]; x = log10(E/GeV);" << endl;
211 for (
int i = 0;
i != 4; ++
i) {
212 cout <<
"\tstatic const double " << (char) (
'a' +
i) <<
" = " <<
SCIENTIFIC(10,3) << f1.GetParameter(
i) <<
";" << endl;
216 double Tmin = sqrt(EMIN*EMIN - MASS_ELECTRON*MASS_ELECTRON);
218 double Emin = 1 * MASS_LEPTON;
219 double Emax = 5 * MASS_LEPTON;
221 for (
double E = 0.5 * (Emin + Emax); ;
E = 0.5 * (Emin + Emax)) {
223 const double gamma =
E / MASS_LEPTON;
224 const double beta = sqrt(1.0 - 1.0/(gamma*gamma));
226 const double Tmax = (lepton ==
electron ?
227 (0.5*sqrt(
E*
E - MASS_ELECTRON*MASS_ELECTRON)) :
228 (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) / (1.0 + 2.0*gamma*RATIO + RATIO*RATIO));
230 if (fabs(Tmax - Tmin) < precision) {
231 cout <<
"\tstatic const double Emin = " <<
FIXED(7,5) <<
E <<
"; // [GeV]" << endl;
Utility class to parse command line options.
int main(int argc, char *argv[])
static const uint32_t K[64]
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
then set_variable DIR else fatal Wrong number of arguments fi for INPUT_FILE in ls rt $DIR stage * log
static const double MASS_MUON
muon mass [GeV]
static const double INDEX_OF_REFRACTION_WATER
Average index of refraction of water corresponding to the group velocity.
Auxiliary data structure for floating point format specification.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
Type definition of range.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
set_variable E_E log10(E_{fit}/E_{#mu})"
static const double MASS_TAU
tau mass [GeV]
do set_variable OUTPUT_DIRECTORY $WORKDIR T
T pow(const T &x, const double y)
Power .
static const double MASS_ELECTRON
electron mass [GeV]
General purpose messaging.
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
Auxiliary class to define a range between two values.
Utility class to parse command line options.
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
Auxiliary data structure for floating point format specification.
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
then for LEPTON in muon tau positron electron
#define DEBUG(A)
Message macros.