Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JPlotPDF_L1.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <fstream>
5#include <iomanip>
6#include <vector>
7
8#include "TROOT.h"
9#include "TFile.h"
10#include "TH1D.h"
11
16#include "JPhysics/JPDFTable.hh"
17#include "JPhysics/JPDFTypes.hh"
21#include "JDetector/JModule.hh"
24#include "Jeep/JPrint.hh"
25#include "Jeep/JParser.hh"
26#include "Jeep/JMessage.hh"
27
28
29/**
30 * Auxiliary program to determine PDF of L1 hit.
31 *
32 * \author mdejong
33 */
34int main(int argc, char **argv)
35{
36 using namespace std;
37 using namespace JPP;
38
40
41 string inputFile;
42 string outputFile;
43 double E;
44 double R;
45 double R_Hz;
46 JAngle3D dir;
47 double T_ns;
48 JHistogram_t histogram;
49 int debug;
50
51 try {
52
53 JParser<> zap("Auxiliary program to determine PDF of L1 hit.");
54
55 zap['f'] = make_field(inputFile);
56 zap['o'] = make_field(outputFile) = "pdf.root";
57 zap['E'] = make_field(E, "muon energy [GeV]") = 1.0;
58 zap['R'] = make_field(R, "distance of approach [m]");
59 zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
60 zap['T'] = make_field(T_ns, "time window [ns]") = 10.0; // [ns]
61 zap['B'] = make_field(R_Hz, "background rate [Hz]") = 0.0;
62 zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t();
63 zap['d'] = make_field(debug) = 0;
64
65 zap(argc, argv);
66 }
67 catch(const exception &error) {
68 FATAL(error.what() << endl);
69 }
70
71 typedef JSplineFunction1S_t JFunction1D_t;
74 JPolint1FunctionalGridMap>::maplist JMapList_t;
76
77 JFunction1D_t::JSupervisor supervisor(new JFunction1D_t::JDefaultResult(zero));
78
79 const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
80 SCATTERED_LIGHT_FROM_MUON,
81 DIRECT_LIGHT_FROM_DELTARAYS,
82 SCATTERED_LIGHT_FROM_DELTARAYS,
83 DIRECT_LIGHT_FROM_EMSHOWERS,
84 SCATTERED_LIGHT_FROM_EMSHOWERS };
85
86 const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
87
88 JPDF_t pdf[N];
89
90 for (int i = 0; i != N; ++i) {
91
92 try {
93
94 const string file_name = getFilename(inputFile, pdf_t[i]);
95
96 NOTICE("loading input from file " << file_name << "... " << flush);
97
98 pdf[i].load(file_name.c_str());
99
100 NOTICE("OK" << endl);
101 }
102 catch(const JException& error) {
103 FATAL(error.what() << endl);
104 }
105
106 pdf[i].setExceptionHandler(supervisor);
107 }
108
109 if (!histogram.is_valid()) {
110
111 histogram = JHistogram_t(-50.0, +50.0);
112
113 histogram.setBinWidth(0.1);
114 }
115
116
117 JModule module = getModule<JKM3NeT_t>(1);
118
119 module.rotate(JRotation3D(dir));
120
121
122 TFile out(outputFile.c_str(), "recreate");
123
124 TH1D h0("h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
125
127
128 for (double t1 = histogram.getLowerLimit(); t1 <= histogram.getUpperLimit(); t1 += 0.1) {
129
130 const double t2 = t1 + T_ns;
131
132 JFunction1D_t::result_type y1;
133 JFunction1D_t::result_type y2;
134
135 for (int i = 0; i != N; ++i) {
136
137 JFunction1D_t::result_type p1;
138 JFunction1D_t::result_type p2;
139
140 // add PMTs
141
142 for (const auto& pmt : module) {
143 p1 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t1);
144 p2 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t2);
145 }
146
147 if (is_deltarays(pdf_t[i])) {
148 p1 *= getDeltaRaysFromMuon(E);
149 p2 *= getDeltaRaysFromMuon(E);
150 } else if (is_bremsstrahlung(pdf_t[i])) {
151 p1 *= E;
152 p2 *= E;
153 }
154
155 y1 += p1;
156 y2 += p2;
157 }
158
159 // add background
160
161 y1.f += R_Hz * 1e-9; // function value
162 y1.v += R_Hz * 1e-9 * (t1 - histogram.getLowerLimit()); // integral [Tmin,t1]
163 y2.v += R_Hz * 1e-9 * (t2 - histogram.getLowerLimit()); // integral [Tmin,t2]
164
165 zsp[t1] = y1.f * (1.0 - exp(y1.v - y2.v)); // 1 - P(0)
166 }
167
168 zsp.compile();
169
170 for (int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
171
172 const double t1 = h0.GetBinCenter(ix);
173
174 try {
175
176 const JFunction1D_t::result_type result = zsp(t1);
177
178 const double W = exp(-result.v) * result.f / (1.0 - exp(-result.V));
179
180 h0.SetBinContent(ix, W);
181 }
182 catch(const exception&) {}
183 }
184
185 out.Write();
186 out.Close();
187}
string outputFile
Detector support kit.
TPaveText * p1
Various implementations of functional maps.
General purpose messaging.
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Data structure for optical module.
Auxiliary methods for PDF calculations.
Numbering scheme for PDF types.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
int main(int argc, char **argv)
Auxiliary program to determine PDF of L1 hit.
I/O formatting auxiliaries.
Data structure for a composite optical module.
Definition JModule.hh:75
Data structure for angles in three dimensions.
Definition JAngle3D.hh:35
General exception.
Definition JException.hh:24
virtual const char * what() const override
Get error message.
Definition JException.hh:64
Utility class to parse command line options.
Definition JParser.hh:1698
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
JPDFType_t
PDF types.
Definition JPDFTypes.hh:24
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Simple data structure for histogram binning.
void setBinWidth(const abscissa_type dx, int option=0)
Set bin width.
bool is_valid() const
Check validity of histogram binning.
int getNumberOfBins() const
Get number of bins.
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.