Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
software/JReconstruction/JMuonCompass.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <vector>
5#include <set>
6#include <algorithm>
7#include <memory>
8
9#include "TROOT.h"
10#include "TFile.h"
11#include "TH2D.h"
12
15
16#include "JDAQ/JDAQEventIO.hh"
18
22
24
25#include "JTrigger/JHit.hh"
26#include "JTrigger/JFrame.hh"
28#include "JTrigger/JHitL0.hh"
29#include "JTrigger/JBuildL0.hh"
31
32#include "JSupport/JSupport.hh"
35#include "JSupport/JMeta.hh"
36
41
42#include "JROOT/JRootToolkit.hh"
43
44#include "Jeep/JParser.hh"
45#include "Jeep/JMessage.hh"
46
47/**
48 * \file
49 *
50 * Program to perform compass calibration using reconstructed muon trajectories.
51 *
52 * \author mdejong
53 */
54int main(int argc, char **argv)
55{
56 using namespace std;
57 using namespace JPP;
58 using namespace KM3NETDAQ;
59
60 typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
61 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
63 JACOUSTICS::JEvt>::typelist calibration_types;
64 typedef JMultipleFileScanner<calibration_types> JCalibration_t;
65
66 JParallelFileScanner_t inputFile;
67 string outputFile;
68 JLimit_t& numberOfEvents = inputFile.getLimit();
69 string detectorFile;
70 JCalibration_t calibrationFile;
71 double Tmax_s;
72 string pdfFile;
73 JMuonGandalfParameters_t parameters;
74 int debug;
75
76 try {
77
78 parameters.numberOfPrefits = 1;
79
80 JParser<> zap("Program to perform compass calibration using reconstructed muon trajectories.");
81
82 zap['f'] = make_field(inputFile);
83 zap['o'] = make_field(outputFile) = "compass.root";
84 zap['a'] = make_field(detectorFile);
85 zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
86 zap['T'] = make_field(Tmax_s) = 100.0;
87 zap['n'] = make_field(numberOfEvents) = JLimit::max();
88 zap['P'] = make_field(pdfFile);
89 zap['@'] = make_field(parameters) = JPARSER::initialised();
90 zap['d'] = make_field(debug) = 1;
91
92 zap(argc, argv);
93 }
94 catch(const exception& error) {
95 FATAL(error.what() << endl);
96 }
97
98
99 if (parameters.numberOfPrefits != 1) {
100 WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
101 }
102
104
105 try {
106 load(detectorFile, detector);
107 }
108 catch(const JException& error) {
109 FATAL(error);
110 }
111
112 unique_ptr<JDynamics> dynamics;
113
114 if (!calibrationFile.empty()) {
115
116 try {
117
118 dynamics.reset(new JDynamics(detector, Tmax_s));
119
120 dynamics->load(calibrationFile);
121 }
122 catch(const exception& error) {
123 FATAL(error.what());
124 }
125 }
126
127 const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
128
129 JSummaryFileRouter summary(inputFile, parameters.R_Hz);
130
131 const JMuonGandalf::storage_type storage(pdfFile, parameters.TTS_ns);
132
133 JMuonGandalf fit(parameters, storage, debug);
134
135 typedef vector<JHitL0> JDataL0_t;
136 typedef vector<JHitW0> JDataW0_t;
137
138 const JBuildL0<JHitL0> buildL0;
139
140 const int nx = getNumberOfModules(detector, true);
141 const int ny = 101;
142
143 TH2D h2(h2_t, NULL, nx, -0.5, nx - 0.5, ny, -PI, +PI);
144
145 while (inputFile.hasNext()) {
146
147 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
148
149 multi_pointer_type ps = inputFile.next();
150
151 JDAQEvent* tev = ps;
152 JFIT::JEvt* in = ps;
153
154 summary.update(*tev);
155
156 if (dynamics) {
157 dynamics->update(*tev);
158 }
159
160 in->select(parameters.numberOfPrefits, qualitySorter);
161
162 JDataL0_t dataL0;
163
164 buildL0(*tev, router, true, back_inserter(dataL0));
165
166 for (JFIT::JEvt::const_iterator track = in->begin(); track != in->end(); ++track) {
167
168 const JRotation3D R (getDirection(*track));
169 const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
170 JRange<double> Z_m;
171
172 if (track->hasW(JSTART_LENGTH_METRES)) {
173 Z_m = JZRange(0.0, track->getW(JSTART_LENGTH_METRES)) + JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
174 }
175
176 const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRange<double>(parameters.TMin_ns, parameters.TMax_ns), Z_m);
177
178 // hit selection based on start value
179
181
182 for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
183
184 double rate_Hz = summary.getRate(*i);
185
186 if (rate_Hz <= 0.0) {
187 rate_Hz = summary.getRate();
188 }
189
190 JHitW0 hit(*i, rate_Hz);
191
192 hit.rotate(R);
193
194 if (match(hit)) {
195 data[hit.getModuleID()].push_back(hit);
196 }
197 }
198
199 const JLine3Z ta(tz);
200
201 for (auto& i : data) {
202
203 const JModule& module = router.getModule(i.first);
204
205 // select first hit
206
207 sort(i.second.begin(), i.second.end(), JHitL0::compare);
208
209 JDataW0_t::iterator __end = unique(i.second.begin(), i.second.end(), equal_to<JDAQPMTIdentifier>());
210
211 for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
212
213 const double x = router.getIndex(i.first);
214 const double y = h2.GetYaxis()->GetBinCenter(iy);
215
216 const JRotation3Z Rz(y);
217
218 double chi2 = 0.0;
219
220 for (JDataW0_t::const_iterator p = i.second.begin(); p != __end; ++p) {
221
222 JHitW0 hit(*p);
223
224 hit.rotate_back(R);
225
226 hit.sub(module.getPosition());
227 hit.rotate(Rz);
228 hit.add(module.getPosition());
229
230 hit.rotate(R);
231
232 chi2 += fit(ta, hit);
233 }
234
235 h2.Fill(x, y, chi2);
236 }
237 }
238 }
239 }
240 STATUS(endl);
241
242 TFile out(outputFile.c_str(), "recreate");
243
244 putObject(out, JMeta(argc, argv));
245
246 out << h2;
247
248 out.Write();
249 out.Close();
250}
string outputFile
Data structure for detector geometry and calibration.
Dynamic detector calibration.
Basic data structure for L0 hit.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define WARNING(A)
Definition JMessage.hh:65
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
const int getIndex(const JObjectID &id) const
Get index of module.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Data structure for a composite optical module.
Definition JModule.hh:75
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
Data structure for fit of straight line in positive z-direction.
Definition JLine3Z.hh:40
JAxis3D & rotate_back(const JRotation3D &R)
Rotate back axis.
Definition JAxis3D.hh:240
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition JAxis3D.hh:225
JPosition3D & rotate(const JRotation3D &R)
Rotate.
const JPosition3D & getPosition() const
Get position.
Rotation around Z-axis.
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition JVector3D.hh:158
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class for a hit with background rate value.
Definition JHitW0.hh:23
General purpose class for object reading from a list of file names.
General purpose class for parallel reading of objects from a single file or multiple files.
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
double getRate() const
Get default rate.
Range of values.
Definition JRange.hh:42
Template L0 hit builder.
Definition JBuildL0.hh:38
int getModuleID() const
Get module identifier.
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
JDirection3D getDirection(const Vec &dir)
Get direction.
JPosition3D getPosition(const Vec &pos)
Get position.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JTOOLS::JRange< double > JZRange
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
const char *const h2_t
Name of histogram with results from JMuonCompass.cc.
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
int main(int argc, char **argv)
Detector file.
Definition JHead.hh:227
Acoustic event fit.
Orientation of module.
Dynamic detector calibration.
Definition JDynamics.hh:86
Auxiliary class to match data points with given model.
Auxiliary class for recursive type list generation.
Definition JTypeList.hh:351
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Wrapper class to make final fit of muon trajectory.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
const JLimit & getLimit() const
Get limit.
Definition JLimit.hh:84
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
Auxiliary data structure for sorting of hits.
Definition JHitL0.hh:85