Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JPlotNPE1D.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <fstream>
4#include <iomanip>
5#include <vector>
6
7#include "TROOT.h"
8#include "TFile.h"
9#include "TH1D.h"
10
15#include "JPhysics/JPDFTable.hh"
16#include "JPhysics/JPDFTypes.hh"
17#include "JPhysics/JNPETable.hh"
18#include "JPhysics/JGeanz.hh"
20#include "Jeep/JPrint.hh"
21#include "Jeep/JParser.hh"
22#include "Jeep/JMessage.hh"
23
24
25/**
26 * \file
27 *
28 * Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.
29 * \author jseneca
30 */
31int main(int argc, char **argv)
32{
33 using namespace std;
34 using namespace JPP;
35
37
38 vector<string> inputFile;
39 string outputFile;
40 JAngle3D dir;
41 double D;
43 double E;
45 int debug;
46
47 try {
48
49 JParser<> zap("Program to plot PDF of Cherenkov light from EM-shower using interpolation tables.");
50
51 zap['f'] = make_field(inputFile);
52 zap['o'] = make_field(outputFile) = "";
53 zap['D'] = make_field(dir, "(theta, phi) of PMT [rad]");
54 zap['R'] = make_field(D, "distance from vertex to PMT [m]");
55 zap['x'] = make_field(x, "histogram x-binning") = JHistogram_t();
56 zap['E'] = make_field(E, "Energy [GeV]") = 0.0;
57 zap['N'] = make_field(numberOfPoints) = 10;
58 zap['d'] = make_field(debug) = 2;
59
60 zap(argc, argv);
61 }
62 catch(const exception &error) {
63 FATAL(error.what() << endl);
64 }
65
66
68
69 if (E > 0.0) {
70 for (int i = 0; i != numberOfPoints; ++i) {
71 Z.push_back(geanz.getLength(E, (i + 0.5) / (double) numberOfPoints));
72 }
73 }
74
75
76 typedef JPolint0Function1D_t JFunction1D_t;
77 typedef
84
86
87 const int N = inputFile.size();
88
89 JPDF_t pdf[N];
90 JNPE_t npe[N];
91
92 try {
93
94 for (int i = 0; i != N; ++i) {
95
96 NOTICE("loading input from file " << inputFile[i] << "... " << flush);
97
98 pdf[i].load(inputFile[i].c_str());
99
100 pdf[i].setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
101
102 npe[i] = JNPE_t(pdf[i]);
103
104 npe[i].setExceptionHandler(new JFunction<double, double>::JDefaultResult(JMATH::zero));
105
106 NOTICE("OK" << endl);
107 }
108 }
109 catch(const JException& error) {
110 FATAL(error.what() << endl);
111 }
112
113
114 if (outputFile != "") {
115
116 TFile out(outputFile.c_str(), "recreate");
117
118 if (!x.is_valid()) { x = JHistogram_t(100000, -1.0, +1.0); }
119
120 TH1D h0("h0", "PDF Projection; D [m]; P [npe]", x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit());
121
122
123 for (int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
124
125 const double cd = h0.GetBinCenter(ix);
126
127 double Y = 0.0;
128
129 if (!Z.empty()) {
130
131 const double W = 1.0 / (double) Z.size();
132
133 for (vector<double>::const_iterator z = Z.begin(); z != Z.end(); ++z) {
134
135 const double __D = sqrt(D*D - 2.0*(D*cd)*(*z) + (*z)*(*z));
136 const double __cd = (D * cd - (*z)) / __D;
137
138 for (int i = 0; i != N; ++i) {
139 try {
140 Y += W * npe[i](__D, __cd, dir.getTheta(), dir.getPhi());
141 }
142 catch(const exception& error) {}
143 }
144 }
145
146 } else {
147
148 for (int i = 0; i != N; ++i) {
149 Y += npe[i](D, cd, dir.getTheta(), dir.getPhi());
150 }
151 }
152
153 h0.SetBinContent(ix, Y);
154 }
155
156 out.Write();
157 out.Close();
158
159 } else {
160
161 for (double cd; ; ) {
162
163 cout << "> " << flush;
164 cin >> D >> cd;
165
166 for (int i = 0; i != N; ++i) {
167 cout << i << ' ' << FIXED(9,5) << npe[i](D, cd, dir.getTheta(), dir.getPhi()) << endl;
168 }
169 }
170 }
171}
string outputFile
Various implementations of functional maps.
Longitudinal emission profile EM-shower.
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
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)
Definition JPlotNPE1D.cc:31
I/O formatting auxiliaries.
int numberOfPoints
Definition JResultPDF.cc:22
Data structure for angles in three dimensions.
Definition JAngle3D.hh:35
double getTheta() const
Get theta angle.
Definition JAngle3D.hh:86
double getPhi() const
Get phi angle.
Definition JAngle3D.hh:97
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Custom class for integrated values of the PDF of the arrival time of Cherenkov light.
Definition JNPETable.hh:46
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Simple data structure for histogram binning.
Template class for distance evaluation.
Definition JDistance.hh:24
Template definition of function object interface in multidimensions.
Map list.
Definition JMapList.hh:25
Type definition of a zero degree polynomial interpolation with result type double.
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.