Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JStart.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 
10 #include "evt/Head.hh"
11 #include "evt/Evt.hh"
12 #include "JDAQ/JDAQEvent.hh"
13 #include "JDAQ/JDAQTimeslice.hh"
14 #include "JDAQ/JDAQSummaryslice.hh"
15 
16 #include "JDetector/JDetector.hh"
20 #include "JDetector/JK40Rates.hh"
21 
22 #include "JTrigger/JHit.hh"
23 #include "JTrigger/JFrame.hh"
24 #include "JTrigger/JTimeslice.hh"
25 #include "JTrigger/JHitL0.hh"
26 #include "JTrigger/JHitR1.hh"
27 #include "JTrigger/JBuildL0.hh"
29 
33 #include "JSupport/JSupport.hh"
34 #include "JSupport/JMeta.hh"
35 
36 #include "JFit/JLine1Z.hh"
37 #include "JFit/JEnergyRegressor.hh"
38 #include "JFit/JFitToolkit.hh"
39 #include "JFit/JFitParameters.hh"
40 #include "JFit/JEvt.hh"
41 #include "JFit/JEvtToolkit.hh"
42 #include "JFit/JNPEHit.hh"
43 #include "JFit/JStart.hh"
44 #include "JFit/JMEstimator.hh"
45 #include "JFit/JModel.hh"
46 
47 #include "JTools/JConstants.hh"
48 #include "JTools/JFunction1D_t.hh"
50 
51 #include "JPhysics/JPDFTable.hh"
52 #include "JPhysics/JPDFToolkit.hh"
53 #include "JPhysics/JNPETable.hh"
54 
55 #include "JLang/JComparator.hh"
56 
57 #include "Jeep/JPrint.hh"
58 #include "Jeep/JParser.hh"
59 #include "Jeep/JMessage.hh"
60 
61 
62 /**
63  * \file
64  *
65  * Program to determine start position of muon trajectory.
66  * \author mdejong
67  */
68 int main(int argc, char **argv)
69 {
70  using namespace std;
71  using namespace JPP;
72  using namespace KM3NETDAQ;
73 
74  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
75  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
77 
78  JParallelFileScanner_t inputFile;
79  JFileRecorder<typelist> outputFile;
80  JLimit_t& numberOfEvents = inputFile.getLimit();
81  string detectorFile;
82  string pdfFile;
83  double roadWidth_m;
84  JK40Rates rates_Hz;
85  JTimeRange T_ns;
86  size_t numberOfPrefits;
87  JStart start;
88  bool reprocess;
89  int debug;
90 
91  try {
92 
93  JParser<> zap("Program to perform PDF fit of muon trajectory to data.");
94 
95  zap['f'] = make_field(inputFile);
96  zap['o'] = make_field(outputFile) = "start.root";
97  zap['a'] = make_field(detectorFile);
98  zap['n'] = make_field(numberOfEvents) = JLimit::max();
99  zap['P'] = make_field(pdfFile);
100  zap['R'] = make_field(roadWidth_m) = numeric_limits<double>::max();
101  zap['B'] = make_field(rates_Hz) = JK40Rates::getInstance();
102  zap['T'] = make_field(T_ns) = JTimeRange(-25.0, +50.0);
103  zap['N'] = make_field(numberOfPrefits) = 1;
104  zap['s'] = make_field(start) = JStart(1.0e-3, 1.0e-2);
105  zap['r'] = make_field(reprocess);
106  zap['d'] = make_field(debug) = 1;
107 
108  zap(argc, argv);
109  }
110  catch(const exception& error) {
111  FATAL(error.what() << endl);
112  }
113 
114 
115 
116  cout.tie(&cerr);
117 
118 
119  JDetector detector;
120 
121  try {
122  load(detectorFile, detector);
123  }
124  catch(const JException& error) {
125  FATAL(error);
126  }
127 
128  const JModuleRouter moduleRouter(detector);
129 
130  typedef vector<JHitL0> JDataL0_t;
131  const JBuildL0<JHitL0> buildL0;
132 
133 
134  typedef JRegressor<JEnergy> JRegressor_t;
135 
137  JRegressor_t::T_ns = T_ns;
138 
139 
140  JRegressor_t fit(pdfFile);
141 
142  if (fit.getRmax() < roadWidth_m) {
143  roadWidth_m = fit.getRmax();
144  }
145 
146 
147  outputFile.open();
148 
149  outputFile.put(JMeta(argc, argv));
150 
151  while (inputFile.hasNext()) {
152 
153  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
154 
155  multi_pointer_type ps = inputFile.next();
156 
157  JDAQEvent* tev = ps;
158  JEvt* evt = ps;
159  JEvt out;
160 
161  JEvt::iterator __end = evt->end();
162 
163  if (reprocess) {
164  __end = partition(evt->begin(), __end, JHistory::is_not_event(JMUONSTART));
165  }
166 
167  if (evt->begin() != __end) {
168 
169  copy(evt->begin(), __end, back_inserter(out));
170 
171  if (numberOfPrefits > 0) {
172  advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
173  }
174 
175  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
176 
177  __end = partition(evt->begin(), __end, JHistory::is_event(evt->begin()->getHistory()));
178 
179 
180  JDataL0_t dataL0;
181 
182  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
183 
184 
185  for (JEvt::iterator track = evt->begin(); track != __end; ++track) {
186 
187  const JRotation3D R (getDirection(*track));
188  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
189  const JModel<JLine1Z> match(tz, roadWidth_m, T_ns);
190 
191 
193 
194  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
195 
196  JHitR1 hit(*i);
197 
198  hit.rotate(R);
199 
200  if (match(hit)) {
201  top.insert(hit.getModuleIdentifier());
202  }
203  }
204 
205 
206  vector<JNPEHit> data;
207 
208  const JDetectorSubset_t subdetector(detector, getAxis(*track), roadWidth_m);
209 
210  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
211 
212  size_t count = top.count(module->getID());
213 
214  if (count != 0) {
215  data.push_back(JNPEHit(fit.getNPE(*module, rates_Hz), count));
216  }
217  }
218 
219 
220  double Zmin = 0.0; // relative start position [m]
221  double Zmax = 0.0; // relative end position [m]
222  double npe = 0.0; // npe due to mip up to barycenter of hits
223  double npe_total = 0.0; // npe due to mip up to end of track
224 
225 
226  if (!data.empty()) {
227 
228  sort(data.begin(), data.end(), make_comparator(&JNPEHit::getZ));
229 
230  vector<JNPEHit>::const_iterator track_start = start.find(data. begin(), data. end());
231  vector<JNPEHit>::const_reverse_iterator track_end = start.find(data.rbegin(), data.rend());
232 
233  if (track_start == data. end()) { track_start = data. begin(); }
234  if (track_end == data.rend()) { track_end = data.rbegin(); }
235 
236  Zmin = track_start->getZ();
237  Zmax = track_end ->getZ();
238 
239  for (vector<JNPEHit>::const_iterator i = track_start; i != track_end.base(); ++i) {
240 
241  if (i->getZ() <= 0.0) {
242  npe += i->getYA();
243  }
244 
245  npe_total += i->getYA();
246  }
247 
248  DEBUG("vertex " << FIXED(6,1) << track->getW(JCOPY_Z_M, 0.0) << endl);
249 
250  for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
251  DEBUG(FIXED(6,1) << i->getZ() << ' ' <<
252  FIXED(6,4) << i->getY0() << ' ' <<
253  FIXED(6,4) << i->getYA() << ' ' <<
254  setw(1) << i->getN() << ' ' <<
255  FIXED(6,4) << i->getP() << endl);
256  }
257 
258  DEBUG(" --> " << FIXED(6,1) << Zmin << ' ' << FIXED(6,1) << Zmax << endl);
259  }
260 
261  // move track
262 
263  track->move(Zmin, getSpeedOfLight());
264 
265  out.push_back(JFit(*track).add(JMUONSTART));
266 
267  out.rbegin()->setW(JSTART_NPE_MIP, npe);
268  out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, npe_total);
269  out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
270  }
271 
272  // apply default sorter
273 
274  sort(out.begin(), out.end(), qualitySorter);
275  }
276 
277  outputFile.put(out);
278  }
279  STATUS(endl);
280 
281  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
282 
283  io >> outputFile;
284 
285  outputFile.close();
286 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Utility class to parse command line options.
Definition: JParser.hh:1410
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:143
Auxiliary method to locate start and end point of muon trajectory.
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
#define STATUS(A)
Definition: JMessage.hh:61
Recording of objects on file according a format that follows from the file name extension.
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
string outputFile
Data structure for detector geometry and calibration.
Various implementations of functional maps.
JRange< double > JTimeRange
Type definition for time range.
Basic data structure for L0 hit.
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
Definition of fit parameters from various applications.
Basic data structure for time and time over threshold information of hit.
Type list.
Definition: JTypeList.hh:22
Auxiliary class to extract a subset of optical modules from a detector.
const_iterator< T > end() const
Get end of data.
Constants.
I/O formatting auxiliaries.
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:73
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
JAxis3D getAxis(const Trk &track)
Get axis.
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.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
true vertex position along track [m] from JCopy.cc
General purpose messaging.
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JCounter.hh:35
#define FATAL(A)
Definition: JMessage.hh:65
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Direct access to module in detector data structure.
number of photo-electrons up to the barycentre from JStart.cc
Reduced data structure for L1 hit.
Utility class to parse command line options.
ROOT TTree parameter settings.
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
distance between first and last hits in metres from JStart.cc
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:40
JDirection3D getDirection(const Vec &v)
Get direction.
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:72
Data regression method for JFIT::JEnergy.
number of photo-electrons along the whole track from JStart.cc
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
Maximum likelihood estimator (M-estimators).
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JPosition3D getPosition(const Vec &v)
Get position.
int main(int argc, char *argv[])
Definition: Main.cpp:15