Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
JSignalToAWriter.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3
4#include "TRandom3.h"
5#include "TH1D.h"
6
9
10#include "JAcoustics/JToA.hh"
14
17
18#include "JTools/JRange.hh"
19
21#include "JSupport/JMeta.hh"
22
23#include "Jeep/JProperties.hh"
24#include "Jeep/JParser.hh"
25#include "Jeep/JMessage.hh"
26
27namespace JROOT {
28 template<>
29 class JRootCreateFlatTree<JACOUSTICS::JToA> : public std::true_type
30 {};
31}
32
33namespace {
34
35 /**
36 * Get range of sound.
37 *
38 * \param E neutrino energy [GeV]
39 * \param l_abs absorption length [m]
40 * \return range [m]
41 */
42 inline double getRange(const double E, const double l_abs)
43 {
44 const double precision = 1.0e-6;
45 const double y0 = 1.0e-3;
46
47 double xmin = -1.0;
48 double xmax = +2.0;
49
50 for (int i = 0; i != 1000; ++i) {
51
52 const double x = 0.5 * (xmin + xmax);
53 const double R = 1.0e3 * pow(10.0, x);
54
55 const double y = E * 2.0e-9 * exp(-R/l_abs) / R;
56
57 if (fabs(y - y0) <= precision) {
58 return R;
59 }
60
61 if (y > y0)
62 xmin = x;
63 else
64 xmax = x;
65 }
66
67 return 0.0;
68 }
69}
70
71/**
72 * \file
73 * Auxiliary program to write signal acoustic data.
74 * \author mdejong
75 */
76int main(int argc, char **argv)
77{
78 using namespace std;
79 using namespace JPP;
80
82
84
85 JFileRecorder_t outputFile;
86 string detectorFile;
87 double T_s;
88 size_t numberOfEvents;
89 double E_GeV;
90 int run;
91 int id;
92 JRange_t X;
93 JRange_t Y;
94 JRange_t Z;
95
96 struct {
97 double l_abs = 5.0e3; // absorption lentgh [m]
98 double ZMax_m = 150.0; // signal length [m]
99 } parameters;
100
101 JSoundVelocity V = getSoundVelocity; // default sound velocity
102 ULong_t seed;
103 int debug;
104
105 try {
106
107 JProperties properties;
108
109 properties.insert(gmake_property(parameters.l_abs));
110 properties.insert(gmake_property(parameters.ZMax_m));
111
112 JParser<> zap(" Auxiliary program to write signal acoustic data.");
113
114 zap['a'] = make_field(detectorFile, "detector.");
115 zap['@'] = make_field(properties) = JPARSER::initialised();
116 zap['o'] = make_field(outputFile, "output file");
117 zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
118 zap['E'] = make_field(E_GeV, "neutrino energy [GeV]");
119 zap['R'] = make_field(run, "run number") = -1;
120 zap['T'] = make_field(T_s, "time interval between events [s]");
121 zap['n'] = make_field(numberOfEvents) = 1;
122 zap['X'] = make_field(X, "x-range") = JRange_t(0.0, 0.0);
123 zap['Y'] = make_field(Y, "y-range") = JRange_t(0.0, 0.0);
124 zap['Z'] = make_field(Z, "z-range") = JRange_t(0.0, 0.0);
125 zap['W'] = make_field(id, "waveform identifier");
126 zap['S'] = make_field(seed, "seed") = 0;
127 zap['d'] = make_field(debug, "debug") = 0;
128
129 zap(argc, argv);
130 }
131 catch(const exception &error) {
132 FATAL(error.what() << endl);
133 }
134
135
136 gRandom->SetSeed(seed);
137
139
140 try {
141 load(detectorFile, detector);
142 }
143 catch(const JException& error) {
144 FATAL(error);
145 }
146
147 const double RMax_m = getRange(E_GeV, parameters.l_abs);
148
149 const JCylinder3D cylinder(detector.begin(), detector.end());
150 const JPosition3D center(cylinder.getX(), cylinder.getY(), 0.5 * (cylinder.getZmin() + cylinder.getZmax()));
151
152 TH1D h1("h1", NULL, 60, 6.0, 12.0);
153
154 for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {
155
156 const double x = h1.GetXaxis()->GetBinCenter(i);
157 const double E = pow(10.0, x);
158 const double R = getRange(E, parameters.l_abs);
159
160 h1.SetBinContent(i, R);
161 }
162
163 outputFile.open();
164
165 outputFile.put(JMeta(argc, argv));
166 outputFile.put(h1);
167
168 JToA toa;
169
170 toa.DETID = detector.getID();
171 toa.RUN = run;
172 toa.WAVEFORMID = id;
173 toa.QUALITYFACTOR = 2.0e3;
174 toa.QUALITYNORMALISATION = 0.0e3;
175
176 for (size_t i = 0; i != numberOfEvents; ++i) {
177
178 JPosition3D pos(gRandom->Uniform(X.getLowerLimit(), X.getUpperLimit()),
179 gRandom->Uniform(Y.getLowerLimit(), Y.getUpperLimit()),
180 gRandom->Uniform(Z.getLowerLimit(), Z.getUpperLimit()));
181
182 pos.add(center);
183
184 Double_t x, y, z;
185
186 gRandom->Sphere(x, y, z, 1.0);
187
188 const JVersor3D dir(x,y,z);
189
190 const JRotation3D R(dir);
191
192 const double t0 = T_s * i; // event time
193
194 size_t number_of_hits = 0;
195
196 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
197
198 if (module->getFloor() != 0) {
199
200 JPosition3D P(module->getPosition());
201
202 P.sub(pos);
203 P.rotate(R);
204
205 if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= RMax_m && fabs(P.getZ()) <= parameters.ZMax_m) {
206
207 const JReceiver receiver(module->getID(), module->getPosition(), module->getT0() * 1.0e-9);
208
209 const double t1 = V.getTime(pos.getDistance(receiver), pos.getZ(), receiver.getZ());
210
211 toa.DOMID = receiver.getID();
212 toa.TOA_NS = llrint(1E9*receiver.putT(t0 + t1));
213
214 outputFile.put(toa);
215
216 number_of_hits += 1;
217 }
218 }
219 }
220
221 STATUS("event: " << setw(6) << i << ' '
222 << FIXED(12,2) << t0 << ' '
223 << FIXED(9,2) << (pos.getX() - center.getX()) << ' '
224 << FIXED(9,2) << (pos.getY() - center.getY()) << ' '
225 << FIXED(9,2) << (pos.getZ() - center.getZ()) << ' '
226 << FIXED(6,3) << dir.getDX() << ' '
227 << FIXED(6,3) << dir.getDY() << ' '
228 << FIXED(6,3) << dir.getDZ() << ' '
229 << setw(4) << number_of_hits << endl);
230 }
231 outputFile.close();
232
233 return 0;
234}
ROOT TTree parameter settings.
string outputFile
Data structure for detector geometry and calibration.
Recording of objects on file according a format that follows from the file name extension.
General purpose messaging.
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
ROOT I/O of application specific meta data.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Acoustic receiver.
int main(int argc, char **argv)
Sound velocity.
Acoustic event.
Detector data structure.
Definition JDetector.hh:96
Utility class to parse parameter values.
double getY() const
Get y position.
Definition JVector2D.hh:74
double getX() const
Get x position.
Definition JVector2D.hh:63
double getZmin() const
Get minimal z position.
double getZmax() const
Get maximal z position.
Data structure for position in three dimensions.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
double getY() const
Get y position.
Definition JVector3D.hh:104
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition JVector3D.hh:270
double getZ() const
Get z position.
Definition JVector3D.hh:115
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition JVector3D.hh:158
double getX() const
Get x position.
Definition JVector3D.hh:94
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
double getDY() const
Get y direction.
Definition JVersor3D.hh:106
double getDX() const
Get x direction.
Definition JVersor3D.hh:95
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
Utility class to parse command line options.
Definition JParser.hh:1698
Object writing to file.
Range of values.
Definition JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
const double xmax
const double xmin
Auxiliary classes and methods for acoustic position calibration.
static const JSoundVelocity getSoundVelocity(1541.0, -17.0e-3, -2000.0)
Function object for velocity of sound.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for ROOT I/O.
range_type getRange(TAxis *pAxis)
Get range of given axis.
Definition JRootfit.hh:52
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Type definition of range.
Definition JHead.hh:43
Detector file.
Definition JHead.hh:227
Acoustic receiver.
Definition JReceiver.hh:30
double putT(const double t_s) const
Get uncorrected time.
Definition JReceiver.hh:84
Implementation for depth dependend velocity of sound.
virtual double getTime(const double D_m, const double z1, const double z2) const override
Get propagation time of sound.
Time-of-arrival data from acoustic piezo sensor or hydrophone.
Definition JToA.hh:28
uint32_t DOMID
DAQ run number.
Definition JToA.hh:34
int64_t TOA_NS
Unique ID of the waveform that best described the signal around TOA_NS.
Definition JToA.hh:36
uint32_t QUALITYFACTOR
The ticks (16ns) part of the DAQ frame timestamp.
Definition JToA.hh:39
uint32_t QUALITYNORMALISATION
A measure of how good the waveform match was to the signal.
Definition JToA.hh:40
int32_t WAVEFORMID
DOM unique identifeir.
Definition JToA.hh:35
int32_t DETID
Definition JToA.hh:32
int32_t RUN
detector identifier
Definition JToA.hh:33
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72