Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JLine1Z.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <fstream>
6#include <vector>
7
8#include "TROOT.h"
9#include "TFile.h"
10#include "TH1D.h"
11#include "TRandom3.h"
12
15#include "JFit/JLine1Z.hh"
17#include "JROOT/JRootToolkit.hh"
18#include "JTools/JQuantile.hh"
19
20#include "Jeep/JParser.hh"
21#include "Jeep/JMessage.hh"
22
23
24/**
25 * \file
26 *
27 * Test program for JFIT::JEstimator<JLine1Z> fit.
28 * \author mdejong
29 */
30int main(int argc, char **argv)
31{
32 using namespace std;
33 using namespace JPP;
34
35 string inputFile;
36 string outputFile;
37 int numberOfEvents;
38 int numberOfPMTs;
39 double precision;
40 int debug;
41
42 try {
43
44 JParser<> zap("Test program for vertex fit.");
45
46 zap['f'] = make_field(inputFile) = "";
47 zap['o'] = make_field(outputFile) = "";
48 zap['n'] = make_field(numberOfEvents) = 1000;
49 zap['N'] = make_field(numberOfPMTs) = 10;
50 zap['e'] = make_field(precision) = 1.0e-5;
51 zap['d'] = make_field(debug) = 3;
52
53 zap(argc, argv);
54 }
55 catch(const exception& error) {
56 FATAL(error.what() << endl);
57 }
58
59
60 TH1D hx("hx", NULL, 101, -1.0, +1.0);
61 TH1D hy("hy", NULL, 101, -1.0, +1.0);
62 TH1D ht("ht", NULL, 101, -1.0, +1.0);
63
64 JQuantile Qx("x");
65 JQuantile Qy("y");
66 JQuantile Qt("t");
67
68
69 typedef vector<JVector3D> JDetector_t;
70
71 JDetector_t detector;
72
73 if (inputFile != "") {
74
75 ifstream in(inputFile.c_str());
76
77 for (double x, y, z; in >> x >> y >> z; ) {
78 detector.push_back(JVector3D(x,y,z));
79 }
80
81 in.close();
82
83 } else {
84
85 const JCylinder3D cylinder(JCircle2D(JVector2D(0.0, 0.0), 100.0), -50.0, +50.0);
86
87 for (int i = 0; i != numberOfPMTs; ++i) {
88 detector.push_back(getRandomPosition(cylinder));
89 }
90 }
91
92
93 typedef JVertex3D JHit_t;
94
95 const double xmin = -1.0;
96 const double xmax = +1.0;
97
98 const double tmin = -1.0;
99 const double tmax = +1.0;
100
101 for (int i = 0; i != numberOfEvents; ++i) {
102
103 // generate line
104
105 const double x = 0.0 + gRandom->Uniform(xmin, xmax);
106 const double y = 0.0 + gRandom->Uniform(xmin, xmax);
107 const double z = 0.0;
108 const double t = 0.0 + gRandom->Uniform(tmin, tmax);
109
110 JLine1Z line(JVector3D(x,y,z), t);
111
112 DEBUG("line "
113 << FIXED(12,5) << line.getX() << ' '
114 << FIXED(12,5) << line.getY() << ' '
115 << FIXED(12,5) << line.getZ() << ' '
116 << FIXED(12,5) << line.getT() << endl);
117
118 // generate hits
119
120 vector<JHit_t> data;
121
122 for (JDetector_t::const_iterator pos = detector.begin(); pos != detector.end(); ++pos) {
123 data.push_back(JHit_t(*pos, line.getT(*pos)));
124 }
125
126 for (vector<JHit_t>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
127 DEBUG("hit "
128 << FIXED(12,5) << hit->getX() << ' '
129 << FIXED(12,5) << hit->getY() << ' '
130 << FIXED(12,5) << hit->getZ() << ' '
131 << FIXED(12,5) << hit->getT() << endl);
132 }
133
134 // fit
135
136 JEstimator<JLine1Z> result(data.begin(), data.end());
137
138 result.setZ(line.getZ(), getSpeedOfLight());
139
140 DEBUG("result "
141 << FIXED(12,5) << result.getX() << ' '
142 << FIXED(12,5) << result.getY() << ' '
143 << FIXED(12,5) << result.getZ() << ' '
144 << FIXED(12,5) << result.getT() << endl);
145
146 hx.Fill(line.getX() - result.getX());
147 hy.Fill(line.getY() - result.getY());
148 ht.Fill(line.getT() - result.getT());
149
150 Qx.put(line.getX() - result.getX());
151 Qy.put(line.getY() - result.getY());
152 Qt.put(line.getT() - result.getT());
153 }
154
155 if (debug >= debug_t) {
156 Qx.print(cout);
157 Qy.print(cout);
158 Qt.print(cout);
159 }
160
161 if (outputFile != "") {
162
163 TFile out(outputFile.c_str(), "recreate");
164
165 out << hx << hy << ht;
166
167 out.Write();
168 out.Close();
169 }
170
171 ASSERT(numberOfEvents > 0);
172
173 ASSERT(fabs(Qx.getMean()) <= precision);
174 ASSERT(fabs(Qy.getMean()) <= precision);
175 ASSERT(fabs(Qt.getMean()) <= precision);
176
177 ASSERT(Qx.getSTDev() <= precision);
178 ASSERT(Qy.getSTDev() <= precision);
179 ASSERT(Qt.getSTDev() <= precision);
180
181 return 0;
182}
string outputFile
Linear fit of JFIT::JLine1Z.
int main(int argc, char **argv)
Definition JLine1Z.cc:30
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Template definition of linear fit.
Definition JEstimator.hh:25
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JLine1Z.hh:114
double getZ(const JPosition3D &pos) const
Get point of emission of Cherenkov light along muon path.
Definition JLine1Z.hh:134
Data structure for circle in two dimensions.
Definition JCircle2D.hh:35
Data structure for vector in two dimensions.
Definition JVector2D.hh:34
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
double getY() const
Get y position.
Definition JVector3D.hh:104
double getX() const
Get x position.
Definition JVector3D.hh:94
Utility class to parse command line options.
Definition JParser.hh:1698
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
Detector file.
Definition JHead.hh:227
Auxiliary class to set-up Hit.
Definition JSirene.hh:58
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition JQuantile.hh:382
double getSTDev() const
Get standard deviation.
Definition JQuantile.hh:281
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
double getMean() const
Get mean value.
Definition JQuantile.hh:252