Jpp  16.0.3
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
software/JReconstruction/JMuonStart.cc File Reference

Program to determine start position of muon trajectory. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <memory>
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDynamics/JDynamics.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JEnergyRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JMEstimator.hh"
#include "JFit/JModel.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonStartParameters_t.hh"
#include "JReconstruction/JStart.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JK40Rates.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JNPETable.hh"
#include "JLang/JComparator.hh"
#include "JLang/JPredicate.hh"
#include "Jeep/JPrint.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 determine start position of muon trajectory.

Author
mdejong

Definition in file software/JReconstruction/JMuonStart.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 121 of file software/JReconstruction/JMuonStart.cc.

122 {
123  using namespace std;
124  using namespace JPP;
125  using namespace KM3NETDAQ;
126 
128  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
129  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
131  JACOUSTICS::JEvt>::typelist calibration_types;
132  typedef JMultipleFileScanner<calibration_types> JCalibration_t;
133 
134  JParallelFileScanner_t inputFile;
136  JLimit_t& numberOfEvents = inputFile.getLimit();
137  string detectorFile;
138  JCalibration_t calibrationFile;
139  double Tmax_s;
140  string pdfFile;
142  JK40Rates rates_Hz;
143  int debug;
144 
145  try {
146 
147  JParser<> zap("Program to perform PDF fit of muon trajectory to data.");
148 
149  zap['f'] = make_field(inputFile);
150  zap['o'] = make_field(outputFile) = "start.root";
151  zap['a'] = make_field(detectorFile);
152  zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
153  zap['T'] = make_field(Tmax_s) = 100.0;
154  zap['n'] = make_field(numberOfEvents) = JLimit::max();
155  zap['P'] = make_field(pdfFile);
157  zap['B'] = make_field(rates_Hz) = KM3NET::getK40Rates();
158  zap['d'] = make_field(debug) = 1;
159 
160  zap(argc, argv);
161  }
162  catch(const exception& error) {
163  FATAL(error.what() << endl);
164  }
165 
166 
167  cout.tie(&cerr);
168 
169 
171 
172  try {
173  load(detectorFile, detector);
174  }
175  catch(const JException& error) {
176  FATAL(error);
177  }
178 
179 
180  getMechanics.load(detector.getID());
181 
182  unique_ptr<JDynamics> dynamics;
183 
184  try {
185 
186  dynamics.reset(new JDynamics(detector, Tmax_s));
187 
188  dynamics->load(calibrationFile);
189  }
190  catch(const exception& error) {
191  if (!calibrationFile.empty()) {
192  FATAL(error.what());
193  }
194  }
195 
196  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
197 
198  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
199 
200 
201  typedef JHitL0 hit_type;
202 
203  typedef vector<hit_type> buffer_type;
204  const JBuildL0<hit_type> buildL0;
205 
206 
207  typedef JRegressor<JEnergy> JRegressor_t;
208 
210  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
211 
212  const JStart start(parameters.Pmin1, parameters.Pmin2, parameters.Nmax2);
213 
214  JRegressor_t fit(pdfFile);
215 
216  if (fit.getRmax() < parameters.roadWidth_m) {
217  parameters.roadWidth_m = fit.getRmax();
218  }
219 
220 
221  outputFile.open();
222 
223  outputFile.put(JMeta(argc, argv));
224 
225  while (inputFile.hasNext()) {
226 
227  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
228 
229  multi_pointer_type ps = inputFile.next();
230 
231  const JDAQEvent* tev = ps;
232  const JFIT::JEvt* in = ps;
233 
234  summary.update(*tev);
235 
236  if (dynamics) {
237  dynamics->update(*tev);
238  }
239 
240  // select start values
241 
242  JFIT::JEvt cp = *in;
243  JFIT::JEvt out;
244 
245  if (parameters.reprocess) {
247  }
248 
249  cp.select(parameters.numberOfPrefits, qualitySorter);
250 
251  if (!cp.empty()) {
252  cp.select(JHistory::is_event(cp.begin()->getHistory()));
253  }
254 
255 
256  buffer_type dataL0;
257 
258  buildL0(*tev, router, true, back_inserter(dataL0));
259 
260 
261  for (JFIT::JEvt::iterator track = cp.begin(); track != cp.end(); ++track) {
262 
263  const JRotation3D R (getDirection(*track));
264  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
265  const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRegressor_t::T_ns);
266 
268 
269  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
270 
271  hit_type hit(*i);
272 
273  hit.rotate(R);
274 
275  if (match(hit)) {
276  top[hit.getModuleID()].insert(hit.getPMTAddress());
277  }
278  }
279 
280  struct JHit_t {
281 
282  double getZ() const { return z; }
283  double getP() const { return p; }
284 
285  double z;
286  double p;
287  };
288 
289  vector<JHit_t> data;
290 
291  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
292 
293  if (summary.hasSummaryFrame(module->getID())) {
294 
295  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
296 
297  JPosition3D pos(module->getPosition());
298 
299  pos.transform(R, tz.getPosition());
300 
301  if (pos.getX() <= parameters.roadWidth_m) {
302 
303  const double z = pos.getZ() - pos.getX() / getTanThetaC();
304  const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), JRegressor_t::T_ns.getLength(), top[module->getID()]);
305 
306  data.push_back({ z, p });
307  }
308  }
309  }
310 
311 
312  double Zmin = 0.0; // relative start position [m]
313  double Zmax = 0.0; // relative end position [m]
314  double npe = 0.0; // npe due to MIP up to barycenter of hits
315  double npe_total = 0.0; // npe due to MIP up to end of track
316 
317  if (!data.empty()) {
318 
319  sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
320 
321  vector<JHit_t>::const_iterator track_start = start.find(data. begin(), data. end());
322  vector<JHit_t>::const_reverse_iterator track_end = start.find(data.rbegin(), data.rend());
323 
324  if (track_start != data. end()) { Zmin = track_start->getZ(); }
325  if (track_end != data.rend()) { Zmax = track_end ->getZ(); }
326 
327  // additional features
328 
329  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
330 
331  JPosition3D pos(module->getPosition());
332 
333  pos.transform(R, tz.getPosition());
334 
335  if (pos.getX() <= parameters.roadWidth_m) {
336 
337  const double z = pos.getZ() - pos.getX() / getTanThetaC();
338 
339  if (z >= Zmin && z <= Zmax) {
340 
341  for (size_t i = 0; i != module->size(); ++i) {
342 
343  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
344 
345  if (getDAQStatus(frame, *module, i) &&
346  getPMTStatus(frame, *module, i) &&
347  !module->getPMT(i).has(PMT_DISABLE)) {
348 
349  JPMT pmt = module->getPMT(i);
350 
351  pmt.transform(R, tz.getPosition());
352 
353  const double ya = fit.getNPE(pmt, 0.0).getYA();
354 
355  npe_total += ya;
356 
357  if (pmt.getZ() <= 0.0) {
358  npe += ya;
359  }
360  }
361  }
362  }
363  }
364  }
365  }
366 
367  // move track
368 
369  track->move(Zmin, getSpeedOfLight());
370 
371  out.push_back(JFIT::JFit(*track).add(JMUONSTART));
372 
373  out.rbegin()->setW(track->getW());
374  out.rbegin()->setW(JSTART_NPE_MIP, npe);
375  out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, npe_total);
376  out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
377  }
378 
379  if (dynamics) {
380 
381  const JDynamics::coverage_type coverage = dynamics->getCoverage();
382 
383  for (JFIT::JEvt::iterator i = out.begin(); i != out.end(); ++i) {
384  i->setW(JPP_COVERAGE_ORIENTATION, coverage.orientation);
385  i->setW(JPP_COVERAGE_POSITION, coverage.position);
386  }
387  }
388 
389  // apply default sorter
390 
391  sort(out.begin(), out.end(), qualitySorter);
392 
393  copy(in->begin(), in->end(), back_inserter(out));
394 
395  outputFile.put(out);
396  }
397  STATUS(endl);
398 
400 
401  io >> outputFile;
402 
403  outputFile.close();
404 }
static const int JMUONSTART
Auxiliary class to set-up Hit.
Definition: JHit_t.hh:25
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
static const int JSTART_NPE_MIP_TOTAL
number of photo-electrons along the whole track from JStart.cc
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition: JComparator.hh:185
double position
coverage of detector by available position calibration [0,1]
Definition: JDynamics.hh:543
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
Regressor function object for fit of muon energy.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
General purpose class for parallel reading of objects from a single file or multiple files...
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
Auxiliary class to test history.
Definition: JHistory.hh:152
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
const JK40Rates & getK40Rates()
Get K40 rates.
Definition: KM3NeT.hh:36
JFit & add(const int type)
Add event to history.
string outputFile
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition: JAxis3D.hh:359
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
static const int JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Data structure for track fit results with history and optional associated values. ...
double orientation
coverage of detector by available orientation calibration [0,1]
Definition: JDynamics.hh:542
Type list.
Definition: JTypeList.hh:22
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Data structure for fit parameters.
Detector file.
Definition: JHead.hh:224
JDirection3D getDirection(const Vec &dir)
Get direction.
Acoustic event fit.
Auxiliary class for recursive type list generation.
Definition: JTypeList.hh:351
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Auxiliary class for start or end point evaluation.
Definition: JStart.hh:21
Data storage class for rate measurements of all PMTs in one module.
Auxiliary class to test history.
Definition: JHistory.hh:104
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
static const int PMT_DISABLE
KM3NeT Data Definitions v2.1.0-12-g9520e6e https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
JPosition3D getPosition(const Vec &pos)
Get position.
File router for fast addressing of summary data.
void load(const std::string &file_name)
Load mechanical model parameters from file.
Definition: JMechanics.hh:142
int debug
debug level
Definition: JSirene.cc:63
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
Data structure for coverage of dynamic calibrations.
Definition: JDynamics.hh:541
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits...
Definition: JFitToolkit.hh:41
#define FATAL(A)
Definition: JMessage.hh:67
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
Dynamic detector calibration.
Definition: JDynamics.hh:62
bool getPMTStatus(const JStatus &status)
Test status of PMT.
const double getSpeedOfLight()
Get speed of light.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Data structure for set of track fit results.
General purpose class for object reading from a list of file names.
Data structure for L0 hit.
Definition: JHitL0.hh:27
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double getProbability(const size_t N, const size_t M, const JK40Rates &R_Hz, const double T_ns)
Get probability due to random background.
Definition: JFitToolkit.hh:248
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
then cp
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Orientation of module.
Object reading from a list of files.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
Template L0 hit builder.
Definition: JBuildL0.hh:35
do set_variable DETECTOR_TXT $WORKDIR detector
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:42
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
void select(const JSelector_t &selector)
Select fits.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
then usage $script< detector file >< tripodfile >< stage > input file nInput files correspond to the output of JAcousticsEventBuilder[.sh] nFirst stage eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY eval JPrintDetector a $DETECTOR O CAN source JAcoustics sh $DETECTOR_ID typeset A CONFIGURATION for key in Tmax_s
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application