Jpp  pmt_effective_area_update_2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JDrawLED.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 "Jeep/JParser.hh"
11 #include "Jeep/JMessage.hh"
12 #include "JPhysics/JLED.hh"
13 #include "JPhysics/JPDFToolkit.hh"
14 #include "JPhysics/Antares.hh"
15 #include "JPhysics/KM3NeT.hh"
16 
17 #include "JTools/JSpline.hh"
18 #include "JTools/JQuantiles.hh"
19 #include "JTools/JToolsToolkit.hh"
20 
21 
23 
24 std::istream& operator>>(std::istream& in, orientation& x) { return in >> x.first >> x.second; }
25 std::ostream& operator<<(std::ostream& out, const orientation& x) { return out << x.first << ' ' << x.second; }
26 
27 
28 inline double henyey_greenstein(const double ct)
29 {
30  static const double a = 0.924;
31 
32  return ANTARES::henyey_greenstein(a, ct);
33 }
34 
35 
36 inline double rayleigh(const double ct)
37 {
38  static const double a = 0.853;
39 
40  return ANTARES::rayleigh(a, ct);
41 }
42 
43 
44 /**
45  * Light yield from %LED (number of p.e. per unit solid angle per unit time).
46  */
47 class LED : public JPHYSICS::JAbstractLED {
48 public:
49  /**
50  * Constructor.
51  */
52  LED()
53  {}
54 
55  /**
56  * Light yield from %LED (number of p.e. per unit solid angle per unit time).
57  *
58  * \param ct zenith angle of emission
59  * \param phi azimuth angle of emission
60  * \param dt time of emission [ns]
61  * \return d^2P / dOmega dt [npe/ns/sr]
62  */
63  double getLightFromLED(const double ct,
64  const double phi,
65  const double dt) const
66  {
67  using namespace JPP;
68 
69  static const double sigma = 2.0;
70 
71  const double x = dt / sigma;
72 
73  return exp(-0.5*x*x) / (sigma*sqrt(2*PI)) / (4*PI);
74  }
75 };
76 
77 
78 /**
79  * Angular acceptence of PMT
80  *
81  * \param x cosine of angle of incidence
82  * \return probability
83  */
84 inline double getAngularAcceptance(const double x)
85 {
87 }
88 
89 
90 /**
91  * \file
92  *
93  * Example program to draw PDF from %LED beacon.
94  * \author mdejong
95  */
96 int main(int argc, char **argv)
97 {
98  using namespace std;
99 
100  string outputFile;
101  int numberOfPoints;
102  int number_of_points;
103  double epsilon;
104  double D;
105  double ct;
106  orientation dir;
107  double wavelength;
108  double absorptionLength;
109  double scatteringLength;
110  int debug;
111 
112  try {
113 
114  JParser<> zap("Example program to draw PDF from LED beacon.");
115 
116  zap['o'] = make_field(outputFile) = "led.root";
117  zap['n'] = make_field(numberOfPoints) = 25;
118  zap['N'] = make_field(number_of_points) = 25;
119  zap['e'] = make_field(epsilon) = 1.0e-10;
120  zap['A'] = make_field(absorptionLength) = 50.0;
121  zap['S'] = make_field(scatteringLength) = 50.0;
122  zap['R'] = make_field(D);
123  zap['c'] = make_field(ct);
124  zap['D'] = make_field(dir);
125  zap['w'] = make_field(wavelength) = 470;
126  zap['d'] = make_field(debug) = 2;
127 
128  zap(argc, argv);
129  }
130  catch(const exception &error) {
131  FATAL(error.what() << endl);
132  }
133 
134 
135  const double theta = dir.first;
136  const double phi = dir.second;
137 
138  using namespace JPP;
139 
140  // set global parameters
141 
142  const double P_atm = 240.0; // ambient pressure [Atm]
143  const double tmin = -10.0; // minimal time of emission [ns]
144  const double tmax = +10.0; // maximal time of emission [ns]
145  const double A = 440e-4; // ANTARES photo-cathode area [m2]
146  const double R_Hz = 0.0e3; // singles rate [Hz]
147 
148  const JDispersion dispersion(P_atm);
149 
150  //const double ng = dispersion.getIndexOfRefractionGroup(wavelength);
151  const double lm = scatteringLength / 0.83;
152  const double lr = scatteringLength / 0.17;
153 
154  const double cs = 0.83 * 0.92; // average cosine scattering angle
155 
156  const double l_abs = absorptionLength;
157  const double ls = scatteringLength;
158  const double l_att = l_abs * lr / (l_abs + lr);
159 
160 
161  LED* led = new LED();
162 
163  cout << "Rayleigh scattering length " << lr << " m" << endl;
164  cout << "Mie scattering length " << lm << " m" << endl;
165  cout << "Absorption length " << l_abs << " m" << endl;
166 
167  const JLED_C
168  pdfMie(A,
169  led,
172  l_abs,
173  lm,
175  P_atm,
176  wavelength,
177  tmin,
178  tmax,
179  JCotangent(number_of_points),
181  epsilon);
182 
183 
184  const JLED_C
185  pdfRayleigh(A,
186  led,
189  l_att,
190  lr,
191  rayleigh,
192  P_atm,
193  wavelength,
194  tmin,
195  tmax,
196  JCotangent(number_of_points),
198  epsilon);
199 
200 
201 
202  TFile out(outputFile.c_str(), "recreate");
203 
204  TH1D h0("h0", NULL, 430, -15.0, +200.0);
205  TH1D h1("h1", NULL, 430, -15.0, +200.0);
206  TH1D h2("h2", NULL, 430, -15.0, +200.0);
207  TH1D h3("h3", NULL, 430, -15.0, +200.0);
208  TH1D ha("ha", NULL, 430, -15.0, +200.0);
209 
210  JSplineFunction1S_t f[4];
211  JSplineFunction1S_t f1;
212 
213 
214  for (int i = 1; i <= h1.GetNbinsX(); ++i) {
215 
216  const double t1 = h1.GetBinCenter(i);
217 
218  const double F1 = pdfMie .getDirectLightFromLED (D, ct, theta, phi, t1);
219  const double F2 = pdfMie .getScatteredLightFromLED(D, ct, theta, phi, t1);
220  const double F3 = pdfRayleigh.getScatteredLightFromLED(D, ct, theta, phi, t1);
221 
222 
223  h0.SetBinContent(i, F1 + F2 + F3);
224  h1.SetBinContent(i, F1);
225  h2.SetBinContent(i, F2);
226  h3.SetBinContent(i, F3);
227 
228  f[0][t1] = F1;
229  f[1][t1] = F2;
230  f[2][t1] = F3;
231  f[3][t1] = F1 + F2 + F3;
232 
233  f1[t1] = F1 + F2 + F3;
234  }
235 
236  f1.compile();
237 
238  for (int i = 3; i != sizeof(f)/sizeof(f[0]); ++i) {
239 
240  f[i].compile();
241 
242  JQuantiles quantiles(f[i]);
243 
244 
245  const double lz = l_abs * ls / (l_abs*(1.0-cs) + ls);
246 
247  NOTICE("int " << quantiles.getIntegral() * D*D * exp(D/lz) << endl);
248 
249  DEBUG("D " << D << endl);
250  DEBUG("ct " << ct << endl);
251  DEBUG("theta " << theta << endl);
252  DEBUG("phi " << phi << endl);
253  DEBUG("int " << quantiles.getIntegral() << endl);
254  DEBUG("t1 " << quantiles.getX() << endl);
255  DEBUG("max " << quantiles.getY() << endl);
256  DEBUG("FWHM " << quantiles.getFWHM() << endl);
257  }
258 
259 
260  const double Tmin = ha.GetXaxis()->GetXmin();
261  const double Tmax = ha.GetXaxis()->GetXmax();
262 
263  const double V = f1.rbegin()->getIntegral() + R_Hz * 1e-9 * (Tmax - Tmin); // integral [Tmin,Tmax]
264 
265  for (int i = 1; i <= ha.GetNbinsX(); ++i) {
266 
267  const double t1 = ha.GetBinCenter(i);
268 
269  JSplineFunction1S_t::result_type p = f1(t1);
270 
271  double v = p.v + R_Hz * 1e-9 * (t1 - Tmin);
272  double y = p.f + R_Hz * 1e-9; // function value
273 
274  const double W = exp(-v) * y / (1.0 - exp(-V));
275 
276  ha.SetBinContent(i,W);
277  }
278 
279  delete led;
280 
281  out.Write();
282  out.Close();
283 }
Utility class to parse command line options.
Definition: JParser.hh:1500
int main(int argc, char *argv[])
Definition: Main.cc:15
double rayleigh(const double ct)
Definition: JDrawLED.cc:36
Auxiliary methods for PDF calculations.
double getQE(const double lambda, const bool option)
Quantum efficiency of 10-inch Hamamatsu PMT.
Definition: Antares.hh:376
Optical properties of Antares deep-sea site.
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Interface for emission profile from LED.
Definition: JLED.hh:34
string outputFile
Light yield from LED (number of p.e.
Definition: JDrawLED.cc:47
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"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
double getLightFromLED(const double ct, const double phi, const double dt) const
Light yield from LED (number of p.e.
Definition: JDrawLED.cc:63
then JPizza f
Definition: JPizza.sh:46
Optical properties of KM3NeT deep-sea site.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
This include file contains various recursive methods to operate on multi-dimensional collections...
#define NOTICE(A)
Definition: JMessage.hh:64
static const double PI
Mathematical constants.
int debug
debug level
Definition: JSirene.cc:63
double henyey_greenstein(const double g, const double x)
Auxiliary method to describe light scattering in water (Heneyey-Greenstein)
Definition: Antares.hh:201
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
std::istream & operator>>(std::istream &in, JAANET::JHead &header)
Read header from input.
Definition: JHead.hh:1618
then JCalibrateToT a
Definition: JTuneHV.sh:116
Utility class to parse command line options.
double getAngularAcceptance(const double x)
Angular acceptence of Antares PMT.
Definition: Antares.hh:363
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
int numberOfPoints
Definition: JResultPDF.cc:22
double henyey_greenstein(const double ct)
Definition: JDrawLED.cc:28
std::ostream & operator<<(std::ostream &stream, const CLBCommonHeader &header)
std::pair< double, double > orientation
Definition: JDrawLED.cc:22
LED()
Constructor.
Definition: JDrawLED.cc:52
data_type v[N+1][M+1]
Definition: JPolint.hh:740
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:40
double rayleigh(const double a, const double x)
Auxiliary method to describe light scattering in water (Rayleigh)
Definition: Antares.hh:233
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
double getAngularAcceptance(const double x)
Angular acceptence of PMT.
Definition: JDrawLED.cc:84