Jpp  master_rocky
the software that should make you happy
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 
7 #include "JDetector/JDetector.hh"
9 
10 #include "JAcoustics/JToA.hh"
11 #include "JAcoustics/JSupport.hh"
12 #include "JAcoustics/JReceiver.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 
27 namespace JROOT {
28  template<>
29  class JRootCreateFlatTree<JACOUSTICS::JToA> : public std::true_type
30  {};
31 }
32 
33 namespace {
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 = 0.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 pow(10.0, 0.5 * (xmin + xmax));
68  }
69 }
70 
71 /**
72  * \file
73  * Auxiliary program to write signal acoustic data.
74  * \author mdejong
75  */
76 int main(int argc, char **argv)
77 {
78  using namespace std;
79  using namespace JPP;
80 
82 
83  typedef JRange<double> JRange_t;
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  UInt_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 
159  h1.SetBinContent(i, getRange(E, parameters.l_abs));
160  }
161 
162  outputFile.open();
163 
164  outputFile.put(JMeta(argc, argv));
165  outputFile.put(h1);
166 
167  JToA toa;
168 
169  toa.DETID = detector.getID();
170  toa.RUN = run;
171  toa.WAVEFORMID = id;
172  toa.QUALITYFACTOR = 2.0e3;
173  toa.QUALITYNORMALISATION = 0.0e3;
174 
175  for (size_t i = 0; i != numberOfEvents; ++i) {
176 
177  JPosition3D pos(gRandom->Uniform(X.getLowerLimit(), X.getUpperLimit()),
178  gRandom->Uniform(Y.getLowerLimit(), Y.getUpperLimit()),
179  gRandom->Uniform(Z.getLowerLimit(), Z.getUpperLimit()));
180 
181  pos.add(center);
182 
183  Double_t x, y, z;
184 
185  gRandom->Sphere(x, y, z, 1.0);
186 
187  const JVersor3D dir(x,y,z);
188 
189  const JRotation3D R(dir);
190 
191  const double t0 = T_s * i; // event time
192 
193  size_t number_of_hits = 0;
194 
195  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
196 
197  if (module->getFloor() != 0) {
198 
199  JPosition3D P(module->getPosition());
200 
201  P.sub(pos);
202  P.rotate(R);
203 
204  if (sqrt(P.getX()*P.getX() + P.getY()*P.getY()) <= RMax_m && fabs(P.getZ()) <= parameters.ZMax_m) {
205 
206  const JReceiver receiver(module->getID(), module->getPosition(), module->getT0() * 1.0e-9);
207 
208  const double t1 = V.getTime(pos.getDistance(receiver), pos.getZ(), receiver.getZ());
209 
210  toa.DOMID = receiver.getID();
211  toa.TOA_NS = llrint(1E9*receiver.putT(t0 + t1));
212 
213  outputFile.put(toa);
214 
215  number_of_hits += 1;
216  }
217  }
218  }
219 
220  STATUS("event: " << setw(6) << i << ' '
221  << FIXED(12,2) << t0 << ' '
222  << FIXED(9,2) << (pos.getX() - center.getX()) << ' '
223  << FIXED(9,2) << (pos.getY() - center.getY()) << ' '
224  << FIXED(9,2) << (pos.getZ() - center.getZ()) << ' '
225  << FIXED(6,3) << dir.getDX() << ' '
226  << FIXED(6,3) << dir.getDY() << ' '
227  << FIXED(6,3) << dir.getDZ() << ' '
228  << setw(4) << number_of_hits << endl);
229  }
230  outputFile.close();
231 
232  return 0;
233 }
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:69
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.
Definition: JProperties.hh:501
double getY() const
Get y position.
Definition: JVector2D.hh:74
double getX() const
Get x position.
Definition: JVector2D.hh:63
Cylinder object.
Definition: JCylinder3D.hh:41
double getZmin() const
Get minimal z position.
Definition: JCylinder3D.hh:111
double getZmax() const
Get maximal z position.
Definition: JCylinder3D.hh:122
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
Rotation matrix.
Definition: JRotation3D.hh:114
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
double getY() const
Get y position.
Definition: JVector3D.hh:104
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.
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
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
Definition: JSTDTypes.hh:14
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:26
uint32_t DOMID
DAQ run number.
Definition: JToA.hh:32
int64_t TOA_NS
Unique ID of the waveform that best described the signal around TOA_NS.
Definition: JToA.hh:34
uint32_t QUALITYFACTOR
The ticks (16ns) part of the DAQ frame timestamp.
Definition: JToA.hh:37
uint32_t QUALITYNORMALISATION
A measure of how good the waveform match was to the signal.
Definition: JToA.hh:38
int32_t WAVEFORMID
DOM unique identifeir.
Definition: JToA.hh:33
int32_t DETID
Definition: JToA.hh:30
int32_t RUN
detector identifier
Definition: JToA.hh:31
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