Jpp
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"
17 #include "JDetector/JDetectorSubset.hh"
18 #include "JDetector/JDetectorToolkit.hh"
19 #include "JDetector/JModuleRouter.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"
28 #include "JTrigger/JTriggerParameters.hh"
29 
30 #include "JSupport/JMultipleFileScanner.hh"
31 #include "JSupport/JParallelFileScanner.hh"
32 #include "JSupport/JFileRecorder.hh"
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"
49 #include "JTools/JFunctionalMap_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;
76  typedef JTYPELIST<JAllTypes_t, JEvt>::typelist typelist;
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) = 0;
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 
136  JRegressor_t::debug = debug;
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 
192  multiset<JDAQModuleIdentifier> top;
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(track->getW());
268  out.rbegin()->setW(JSTART_NPE_MIP, npe);
269  out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, npe_total);
270  out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
271  }
272 
273  // apply default sorter
274 
275  sort(out.begin(), out.end(), qualitySorter);
276  }
277 
278  outputFile.put(out);
279  }
280  STATUS(endl);
281 
282  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
283 
284  io >> outputFile;
285 
286  outputFile.close();
287 }
JEnergyRegressor.hh
JMeta.hh
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:34
JTriggerParameters.hh
JFileRecorder.hh
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JTOOLS::getSpeedOfLight
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
JMessage.hh
JROOT::advance
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
JPrint.hh
JAANET::getAxis
JAxis3D getAxis(const Trk &track)
Get axis.
Definition: JAAnetToolkit.hh:244
JLANG::getInstance
T & getInstance(const T &object)
Get static instance from temporary object.
Definition: JObject.hh:73
JFitToolkit.hh
qualitySorter
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
Definition: quality_sorter.cc:10
JTimeslice.hh
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JComparator.hh
JFIT::JMUONSTART
JStart.cc.
Definition: JFitApplications.hh:27
JPDFToolkit.hh
std::vector< JHitL0 >
JEvtToolkit.hh
JSUPPORT::JLimit_t
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:215
JK40Rates.hh
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
JDAQTimeslice.hh
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JNPEHit.hh
main
int main(int argc, char **argv)
Definition: JStart.cc:68
JDAQSummaryslice.hh
JGEOMETRY3D::JPosition3D::rotate
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
KM3NETDAQ::JDAQEvent::end
const_iterator< T > end() const
Get end of data.
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JEvt.hh
JFunction1D_t.hh
JStart.hh
JBuildL0.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
JFunctionalMap_t.hh
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
std::multiset
Definition: JSTDTypes.hh:14
JHit.hh
JModuleRouter.hh
JNPETable.hh
JParallelFileScanner.hh
JHitR1.hh
JMultipleFileScanner.hh
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JParser.hh
JDetectorToolkit.hh
JFIT::JSTART_NPE_MIP
number of photo-electrons up to the barycentre from JStart.cc
Definition: JFitParameters.hh:25
JModel.hh
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
JFIT::JCOPY_Z_M
true vertex position along track [m] from JCopy.cc
Definition: JFitParameters.hh:34
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JFitParameters.hh
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
JFIT::JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
Definition: JFitParameters.hh:27
JMEstimator.hh
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JDAQEvent.hh
JDetector.hh
JLANG::make_comparator
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
JDetectorSubset.hh
JHitL0.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JFIT::JSTART_NPE_MIP_TOTAL
number of photo-electrons along the whole track from JStart.cc
Definition: JFitParameters.hh:26
JPDFTable.hh
JFrame.hh
JLine1Z.hh
Generated by   doxygen 1.8.16