Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JMuonCompass.cc File Reference

Program to perform compass calibration using reconstructed muon trajectories. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include <memory>
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDynamics/JDynamics.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JMuonGandalf.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"
#include "JReconstruction/JMuonCompass.hh"
#include "JROOT/JRootToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program to perform compass calibration using reconstructed muon trajectories.

Author
mdejong

Definition in file software/JReconstruction/JMuonCompass.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 54 of file software/JReconstruction/JMuonCompass.cc.

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);
130
131 const JMuonGandalf::storage_type storage(pdfFile, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), 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 JHitW0 hit(*i, summary.getRate(*i, parameters.R_Hz));
185
186 hit.rotate(R);
187
188 if (match(hit)) {
189 data[hit.getModuleID()].push_back(hit);
190 }
191 }
192
193 const JLine3Z ta(tz);
194
195 for (auto& i : data) {
196
197 const JModule& module = router.getModule(i.first);
198
199 // select first hit
200
201 sort(i.second.begin(), i.second.end(), JHitL0::compare);
202
203 JDataW0_t::iterator __end = unique(i.second.begin(), i.second.end(), equal_to<JDAQPMTIdentifier>());
204
205 for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
206
207 const double x = router.getIndex(i.first);
208 const double y = h2.GetYaxis()->GetBinCenter(iy);
209
210 const JRotation3Z Rz(y);
211
212 double chi2 = 0.0;
213
214 for (JDataW0_t::const_iterator p = i.second.begin(); p != __end; ++p) {
215
216 JHitW0 hit(*p);
217
218 hit.rotate_back(R);
219
220 hit.sub(module.getPosition());
221 hit.rotate(Rz);
222 hit.add(module.getPosition());
223
224 hit.rotate(R);
225
226 chi2 += fit(ta, hit);
227 }
228
229 h2.Fill(x, y, chi2);
230 }
231 }
232 }
233 }
234 STATUS(endl);
235
236 TFile out(outputFile.c_str(), "recreate");
237
238 putObject(out, JMeta(argc, argv));
239
240 out << h2;
241
242 out.Write();
243 out.Close();
244}
string outputFile
#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
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
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
JPosition3D & rotate(const JRotation3D &R)
Rotate.
const JPosition3D & getPosition() const
Get position.
Rotation around Z-axis.
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.
Range of values.
Definition JRange.hh:42
Template L0 hit builder.
Definition JBuildL0.hh:38
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.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
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
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