Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonStart.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 #include <set>
7 #include <algorithm>
8 #include <memory>
9 
13 #include "JDAQ/JDAQEventIO.hh"
14 #include "JDAQ/JDAQTimesliceIO.hh"
16 
17 #include "JDetector/JDetector.hh"
21 
22 #include "JDynamics/JDynamics.hh"
23 
24 #include "JPhysics/JK40Rates.hh"
25 
26 #include "JTrigger/JHit.hh"
27 #include "JTrigger/JHitR1.hh"
28 #include "JTrigger/JBuildL0.hh"
30 
34 #include "JSupport/JSupport.hh"
35 #include "JSupport/JMeta.hh"
36 
37 #include "JFit/JLine1Z.hh"
38 #include "JFit/JEnergyRegressor.hh"
39 #include "JFit/JFitToolkit.hh"
40 #include "JFit/JNPEHit.hh"
41 #include "JFit/JMEstimator.hh"
42 #include "JFit/JModel.hh"
43 
44 #include "JReconstruction/JEvt.hh"
48 
49 #include "JPhysics/JConstants.hh"
50 #include "JTools/JFunction1D_t.hh"
52 
53 #include "JPhysics/JPDFTable.hh"
54 #include "JPhysics/JPDFToolkit.hh"
55 #include "JPhysics/JNPETable.hh"
56 
57 #include "JLang/JComparator.hh"
58 
59 #include "Jeep/JPrint.hh"
60 #include "Jeep/JParser.hh"
61 #include "Jeep/JMessage.hh"
62 
63 
64 /**
65  * \file
66  *
67  * Program to determine start position of muon trajectory.
68  * \author mdejong
69  */
70 int main(int argc, char **argv)
71 {
72  using namespace std;
73  using namespace JPP;
74  using namespace KM3NETDAQ;
75 
77  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
78  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
80  JACOUSTICS::JEvt>::typelist calibration_types;
81  typedef JMultipleFileScanner<calibration_types> JCalibration_t;
82 
83  JParallelFileScanner_t inputFile;
85  JLimit_t& numberOfEvents = inputFile.getLimit();
86  string detectorFile;
87  JCalibration_t calibrationFile;
88  double Tmax_s;
89  string pdfFile;
91  int debug;
92 
93  try {
94 
95  JParser<> zap("Program to perform PDF fit of muon trajectory to data.");
96 
97  zap['f'] = make_field(inputFile);
98  zap['o'] = make_field(outputFile) = "start.root";
99  zap['a'] = make_field(detectorFile);
100  zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
101  zap['T'] = make_field(Tmax_s) = 100.0;
102  zap['n'] = make_field(numberOfEvents) = JLimit::max();
103  zap['P'] = make_field(pdfFile);
105  zap['d'] = make_field(debug) = 1;
106 
107  zap(argc, argv);
108  }
109  catch(const exception& error) {
110  FATAL(error.what() << endl);
111  }
112 
113 
114  cout.tie(&cerr);
115 
116 
118 
119  try {
120  load(detectorFile, detector);
121  }
122  catch(const JException& error) {
123  FATAL(error);
124  }
125 
126 
127  getMechanics.load(detector.getID());
128 
129  unique_ptr<JDynamics> dynamics;
130 
131  try {
132 
133  dynamics.reset(new JDynamics(detector, Tmax_s));
134 
135  dynamics->load(calibrationFile);
136  }
137  catch(const exception& error) {
138  if (!calibrationFile.empty()) {
139  FATAL(error.what());
140  }
141  }
142 
143  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
144 
145 
146  typedef JHitR1 hit_type;
147 
148  typedef vector<hit_type> buffer_type;
149  const JBuildL0<hit_type> buildL0;
150 
151 
152  typedef JRegressor<JEnergy> JRegressor_t;
153 
155  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
156 
157  const JStart start(parameters.Pmin1, parameters.Pmin2);
158 
159  JRegressor_t fit(pdfFile);
160 
161  if (fit.getRmax() < parameters.roadWidth_m) {
162  parameters.roadWidth_m = fit.getRmax();
163  }
164 
165 
166  outputFile.open();
167 
168  outputFile.put(JMeta(argc, argv));
169 
170  while (inputFile.hasNext()) {
171 
172  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
173 
174  multi_pointer_type ps = inputFile.next();
175 
176  const JDAQEvent* tev = ps;
177  const JFIT::JEvt* in = ps;
178 
179  if (dynamics) {
180  dynamics->update(*tev);
181  }
182 
183  // select start values
184 
185  JFIT::JEvt cp = *in;
186  JFIT::JEvt out;
187 
188  if (parameters.reprocess) {
190  }
191 
192  cp.select(parameters.numberOfPrefits, qualitySorter);
193 
194  if (!cp.empty()) {
195  cp.select(JHistory::is_event(cp.begin()->getHistory()));
196  }
197 
198 
199  buffer_type data;
200 
201  buildL0(*tev, router, true, back_inserter(data));
202 
203 
204  for (JFIT::JEvt::iterator track = cp.begin(); track != cp.end(); ++track) {
205 
206  const JRotation3D R (getDirection(*track));
207  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
208  const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRegressor_t::T_ns);
209 
210 
212 
213  for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
214 
215  hit_type hit(*i);
216 
217  hit.rotate(R);
218 
219  if (match(hit)) {
220  top.insert(hit.getModuleIdentifier());
221  }
222  }
223 
224  vector<JNPEHit> data;
225 
226  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
227 
228  const size_t count = top.count(i->getID());
229 
230  if (count != 0) {
231 
232  JPosition3D pos(i->getPosition());
233 
234  pos.transform(R, tz.getPosition());
235 
236  if (pos.getX() <= parameters.roadWidth_m) {
237 
238  JModule module = *i;
239 
240  module.transform(R, tz.getPosition());
241 
242  data.push_back(JNPEHit(fit.getNPE(module.begin(), module.end(), parameters.R_Hz), count));
243  }
244  }
245  }
246 
247 
248  double Zmin = 0.0; // relative start position [m]
249  double Zmax = 0.0; // relative end position [m]
250  double npe = 0.0; // npe due to mip up to barycenter of hits
251  double npe_total = 0.0; // npe due to mip up to end of track
252 
253 
254  if (!data.empty()) {
255 
256  sort(data.begin(), data.end(), make_comparator(&JNPEHit::getZ));
257 
258  vector<JNPEHit>::const_iterator track_start = start.find(data. begin(), data. end());
259  vector<JNPEHit>::const_reverse_iterator track_end = start.find(data.rbegin(), data.rend());
260 
261  if (track_start == data. end()) { track_start = data. begin(); }
262  if (track_end == data.rend()) { track_end = data.rbegin(); }
263 
264  Zmin = track_start->getZ();
265  Zmax = track_end ->getZ();
266 
267  for (vector<JNPEHit>::const_iterator i = track_start; i != track_end.base(); ++i) {
268 
269  if (i->getZ() <= 0.0) {
270  npe += i->getYA();
271  }
272 
273  npe_total += i->getYA();
274  }
275 
276  DEBUG("vertex " << FIXED(6,1) << track->getW(JCOPY_Z_M, 0.0) << endl);
277 
278  for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
279  DEBUG(FIXED(6,1) << i->getZ() << ' ' <<
280  FIXED(6,4) << i->getY0() << ' ' <<
281  FIXED(6,4) << i->getYA() << ' ' <<
282  setw(1) << i->getN() << ' ' <<
283  FIXED(6,4) << i->getP() << endl);
284  }
285 
286  DEBUG(" --> " << FIXED(6,1) << Zmin << ' ' << FIXED(6,1) << Zmax << endl);
287  }
288 
289  // move track
290 
291  track->move(Zmin, getSpeedOfLight());
292 
293  out.push_back(JFIT::JFit(*track).add(JMUONSTART));
294 
295  out.rbegin()->setW(track->getW());
296  out.rbegin()->setW(JSTART_NPE_MIP, npe);
297  out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, npe_total);
298  out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
299  }
300 
301  if (dynamics) {
302 
303  const JDynamics::coverage_type coverage = dynamics->getCoverage();
304 
305  for (JFIT::JEvt::iterator i = out.begin(); i != out.end(); ++i) {
306  i->setW(JPP_COVERAGE_ORIENTATION, coverage.orientation);
307  i->setW(JPP_COVERAGE_POSITION, coverage.position);
308  }
309  }
310 
311  // apply default sorter
312 
313  sort(out.begin(), out.end(), qualitySorter);
314 
315  copy(in->begin(), in->end(), back_inserter(out));
316 
317  outputFile.put(out);
318  }
319  STATUS(endl);
320 
322 
323  io >> outputFile;
324 
325  outputFile.close();
326 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
static const int JMUONSTART
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
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Auxiliary methods for PDF calculations.
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
Data structure for a composite optical module.
Definition: JModule.hh:57
Auxiliary method to locate start and end point of muon trajectory.
double position
coverage of detector by available position calibration [0,1]
Definition: JDynamics.hh:536
static JDetectorMechanics getMechanics
Function object to get string mechanics.
Definition: JMechanics.hh:243
Regressor function object for JEnergy fit.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Recording of objects on file according a format that follows from the file name extension.
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
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
T find(T __begin, T __end) const
Get start or end point of muon trajectory.
Definition: JStart.hh:108
Basic data structure for time and time over threshold information of hit.
JFit & add(const int type)
Add event to history.
string outputFile
Data structure for detector geometry and calibration.
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
Various implementations of functional maps.
static const int JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Data structure for track fit results.
double orientation
coverage of detector by available orientation calibration [0,1]
Definition: JDynamics.hh:535
Type list.
Definition: JTypeList.hh:22
Scanning of objects from a single file according a format that follows from the extension of each fil...
Auxiliary class to extract a subset of optical modules from a detector.
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
static const int JCOPY_Z_M
true vertex position along track [m] from JCopy.cc
Data structure for fit parameters.
I/O formatting auxiliaries.
Detector file.
Definition: JHead.hh:196
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
Auxiliary class to test history.
Definition: JHistory.hh:104
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
ROOT I/O of application specific meta data.
JPosition3D getPosition(const Vec &pos)
Get position.
Physics constants.
Dynamic detector calibration.
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
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
General purpose messaging.
Data structure for coverage of dynamic calibrations.
Definition: JDynamics.hh:534
#define FATAL(A)
Definition: JMessage.hh:67
Direct access to module in detector data structure.
static const int JPP_COVERAGE_ORIENTATION
coverage of dynamic orientation calibration from any Jpp application
Reduced data structure for L1 hit.
Dynamic detector calibration.
Definition: JDynamics.hh:61
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
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.
std::vector< int > count
Definition: JAlgorithm.hh:180
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
void transform(const JRotation3D &R, const JVector3D &pos)
Transformation of geometry (see method JGEOMETRY3D::JPosition3D::transform(const JRotation3D&amp;, const JVector3D&amp;)).
Definition: JModule.hh:382
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
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
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
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
Data regression method for JFIT::JEnergy.
Template L0 hit builder.
Definition: JBuildL0.hh:35
do set_variable DETECTOR_TXT $WORKDIR detector
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:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
void select(const JSelector_t &selector)
Select fits.
Maximum likelihood estimator (M-estimators).
bool qualitySorter(const JRECONSTRUCTION::JFit &first, const JRECONSTRUCTION::JFit &second)
Comparison of fit results.
static const int JPP_COVERAGE_POSITION
coverage of dynamic position calibration from any Jpp application