Jpp  18.5.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JDrawPDF.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 
10 #include "JPhysics/JPDF.hh"
11 #include "JPhysics/Antares.hh"
12 #include "JPhysics/KM3NeT.hh"
13 #include "JTools/JElement.hh"
14 #include "JTools/JFunction1D_t.hh"
15 #include "JTools/JQuantiles.hh"
17 #include "JGeometry3D/JAngle3D.hh"
18 #include "Jeep/JProperties.hh"
19 #include "Jeep/JPrint.hh"
20 #include "Jeep/JParser.hh"
21 #include "Jeep/JMessage.hh"
22 
23 
24 /**
25  * Scaling of absorption and scattering length.
26  */
29 
30 inline double getAbsorptionLength(const double lambda)
31 {
33 }
34 
35 
36 inline double getScatteringLength(const double lambda)
37 {
39 }
40 
41 
42 /**
43  * \file
44  *
45  * Auxiliary program to draw PDF of Cherenkov light from muon.
46  * \author mdejong
47  */
48 int main(int argc, char **argv)
49 {
50  using namespace std;
51  using namespace JPP;
52 
53  typedef JAbstractHistogram<double> JHistogram_t;
54 
55  string outputFile;
56  int numberOfPoints;
57  double epsilon;
58  double E;
59  double R;
60  JAngle3D dir;
61  vector<int> function;
62  JHistogram_t histogram;
63  int debug;
64 
65  try {
66 
67  JProperties properties;
68 
69  properties.insert(gmake_property(JPHYSICS::MODULE_RADIUS_M));
70 
71  JParser<> zap("Auxiliary program to draw PDF of Cherenkov light from muon.");
72 
73  zap['@'] = make_field(properties) = JPARSER::initialised();
74  zap['o'] = make_field(outputFile) = "";
75  zap['n'] = make_field(numberOfPoints, "points for integration") = 25;
76  zap['e'] = make_field(epsilon, "precision for integration") = 1.0e-10;
77  zap['A'] = make_field(absorptionLengthFactor, "scaling factor") = 1.0;
78  zap['S'] = make_field(scatteringLengthFactor, "scaling factor") = 1.0;
79  zap['E'] = make_field(E, "muon energy [GeV]") = 1.0;
80  zap['R'] = make_field(R, "distance of approach [m]");
81  zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
82  zap['F'] = make_field(function, "PDF type");
83  zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t();
84  zap['d'] = make_field(debug) = 0;
85 
86  zap(argc, argv);
87  }
88  catch(const exception &error) {
89  FATAL(error.what() << endl);
90  }
91 
92 
93  const JPDF_C
104  epsilon);
105 
106 
107  if (outputFile == "") {
108 
109  for (double dt; cin >> dt; ) {
110 
111  for (vector<int>::const_iterator F = function.begin(); F != function.end(); ++F) {
112 
113  cout << setw(2) << *F << ' '
114  << SCIENTIFIC(7,1) << E << ' '
115  << FIXED(5,1) << R << ' '
116  << FIXED(5,2) << dir.getTheta() << ' '
117  << FIXED(5,2) << dir.getPhi() << ' '
118  << FIXED(5,1) << dt << ' '
119  << SCIENTIFIC(9,3) << pdf.getLightFromMuon(*F, E, R, dir.getTheta(), dir.getPhi(), dt) * E << endl;
120  }
121  }
122 
123  return 0;
124  }
125 
126 
127  TFile out(outputFile.c_str(), "recreate");
128 
129  //const double z0 = -100.0; // [m]
130  //const double t0 = (-z0 + R * getTanThetaC()) / C; // time [ns]
131  const double t0 = 0.0; // time [ns]
132 
133  if (!histogram.is_valid()) {
134 
135  if (function.size() == 1 && function[0] == DIRECT_LIGHT_FROM_MUON) {
136 
137  histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);
138 
139  histogram.setBinWidth(0.1);
140 
141  } else {
142 
143  histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);
144 
145  histogram.setBinWidth(0.5);
146  }
147  }
148 
149  TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
150 
151  JSplineFunction1S_t f1;
152 
153  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
154 
155  const double dt = h0.GetBinCenter(i) - t0;
156 
157  double value = 0.0;
158 
159  for (vector<int>::const_iterator F = function.begin(); F != function.end(); ++F) {
160  value += pdf.getLightFromMuon(*F, E, R, dir.getTheta(), dir.getPhi(), dt);
161  }
162 
163  h0.SetBinContent(i, value);
164 
165  f1[dt] = value;
166  }
167 
168  f1.setExceptionHandler(new JSplineFunction1S_t::JDefaultResult(JMATH::zero));
169  f1.compile();
170 
171  const double T = 5; // [ns]
172 
173  JQuantiles quantiles(f1);
174 
175  const double t1 = quantiles.getX();
176  const double y = f1(t1 + T).v - f1(t1 - T).v;
177 
178  DEBUG("E " << E << endl);
179  DEBUG("R " << R << endl);
180  DEBUG("theta " << dir.getTheta() << endl);
181  DEBUG("phi " << dir.getPhi() << endl);
182  DEBUG("int " << quantiles.getIntegral() << endl);
183  DEBUG("t1 " << t1 << endl);
184  DEBUG("max " << quantiles.getY() << endl);
185  DEBUG("FWHM " << quantiles.getFWHM() << endl);
186  DEBUG("int[] " << y << endl);
187 
188 
189  out.Write();
190  out.Close();
191 }
Utility class to parse command line options.
Definition: JParser.hh:1514
int main(int argc, char *argv[])
Definition: Main.cc:15
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
The elements in a collection are sorted according to their abscissa values and a given distance opera...
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Utility class to parse parameter values.
Definition: JProperties.hh:497
Properties of Antares PMT and deep-sea water.
double absorptionLengthFactor
Scaling of absorption and scattering length.
Definition: JDrawPD0.cc:24
double getScatteringProbability(const double x)
Function to describe light scattering in water.
Definition: Antares.hh:210
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
double getMinimalWavelength()
Get minimal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:26
then warning Cannot perform comparison test for histogram
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
string outputFile
Utility class to parse parameter values.
direct light from muon
Definition: JPDFTypes.hh:26
const JPolynome f1(1.0, 2.0, 3.0)
Function.
static double MODULE_RADIUS_M
Radius of optical module [m].
Definition: JPDF.hh:40
Properties of KM3NeT PMT and deep-sea water.
double getMaximalWavelength()
Get maximal wavelength for PDF evaluations.
Definition: JPDFToolkit.hh:37
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
double getQE(const double R, const double mu)
Get QE for given ratio of hit probabilities and expectation value of the number of photo-electrons...
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double getAmbientPressure()
Get ambient pressure.
Definition: Antares.hh:40
double getPhotocathodeArea()
Get photo-cathode area of PMT.
Definition: Antares.hh:51
then awk F
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
double getAbsorptionLength(const double lambda)
Get absorption length.
Definition: Antares.hh:63
Utility class to parse command line options.
double getScatteringLength(const double lambda)
Get scattering length.
Definition: Antares.hh:148
int numberOfPoints
Definition: JResultPDF.cc:22
double scatteringLengthFactor
Definition: JDrawPD0.cc:25
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
const double epsilon
Definition: JQuadrature.cc:21
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:68