Jpp
JSimplex.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <vector>
6 #include <map>
7 #include <algorithm>
8 
11 #include "JDAQ/JDAQEventIO.hh"
12 #include "JDAQ/JDAQTimesliceIO.hh"
14 
15 #include "JDetector/JDetector.hh"
18 
19 #include "JTrigger/JHit.hh"
20 #include "JTrigger/JFrame.hh"
21 #include "JTrigger/JTimeslice.hh"
23 #include "JTrigger/JHitL0.hh"
24 #include "JTrigger/JHitL1.hh"
25 #include "JTrigger/JHitR1.hh"
26 #include "JTrigger/JBuildL0.hh"
27 #include "JTrigger/JBuildL1.hh"
28 #include "JTrigger/JBuildL2.hh"
30 
34 #include "JSupport/JSupport.hh"
35 #include "JSupport/JMeta.hh"
36 
37 #include "JTools/JConstants.hh"
38 
39 #include "JFit/JLine3Z.hh"
40 #include "JFit/JSimplex.hh"
41 #include "JFit/JLine3ZRegressor.hh"
42 #include "JFit/JFitToolkit.hh"
43 #include "JFit/JEvt.hh"
44 #include "JFit/JEvtToolkit.hh"
45 #include "JFit/JModel.hh"
46 
47 #include "Jeep/JParser.hh"
48 #include "Jeep/JMessage.hh"
49 
50 
51 /**
52  * \file
53  *
54  * Program to perform JFIT::JRegressor<JLine3Z,JSimplex> fit with I/O of JFIT::JEvt data.
55  * \author mdejong
56  */
57 int main(int argc, char **argv)
58 {
59  using namespace std;
60  using namespace JPP;
61  using namespace KM3NETDAQ;
62 
63  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
64  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
66 
67  JParallelFileScanner_t inputFile;
69  JLimit_t& numberOfEvents = inputFile.getLimit();
70  string detectorFile;
71  double Tmax_ns;
72  double roadWidth_m;
73  double ctMin;
74  size_t numberOfPrefits;
75  double sigma_ns;
76  bool useL0;
77  bool reprocess;
78  int debug;
79 
80  try {
81 
82  JParser<> zap("Program to perform intermediate fit of muon trajectory to data.");
83 
84  zap['f'] = make_field(inputFile);
85  zap['o'] = make_field(outputFile) = "simplex.root";
86  zap['a'] = make_field(detectorFile);
87  zap['n'] = make_field(numberOfEvents) = JLimit::max();
88  zap['T'] = make_field(Tmax_ns) = 20.0; // [ns]
89  zap['R'] = make_field(roadWidth_m) = 200.0; // [m]
90  zap['C'] = make_field(ctMin) = 0.0; //
91  zap['N'] = make_field(numberOfPrefits) = 0;
92  zap['S'] = make_field(sigma_ns);
93  zap['U'] = make_field(useL0);
94  zap['r'] = make_field(reprocess);
95  zap['d'] = make_field(debug) = 1;
96 
97  zap(argc, argv);
98  }
99  catch(const exception& error) {
100  FATAL(error.what() << endl);
101  }
102 
103 
104 
105  cout.tie(&cerr);
106 
107 
108  const double TMAX_NS = 50.0; // [ns]
109 
111 
112  try {
113  load(detectorFile, detector);
114  }
115  catch(const JException& error) {
116  FATAL(error);
117  }
118 
119  const JModuleRouter moduleRouter(detector);
120 
121 
122  typedef vector<JHitL0> JDataL0_t;
123  typedef vector<JHitL1> JDataL1_t;
124  typedef vector<JHitR1> JDataR1_t;
125 
126  JSuperFrame2D<JHit> buffer;
127  const JBuildL0<JHitL0> buildL0;
128  const JBuildL2<JHitL1> buildL2(JL2Parameters(2, Tmax_ns, ctMin));
129 
130 
131  typedef JRegressor<JLine3Z, JSimplex> JRegressor_t;
132 
133  JRegressor_t fit(sigma_ns);
134 
135  fit.estimator.reset(new JMEstimatorLorentzian());
136 
138  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
139 
140 
141  outputFile.open();
142 
143  outputFile.put(JMeta(argc, argv));
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  JEvt* evt = ps;
153  JEvt out;
154 
155  JEvt::iterator __end = evt->end();
156 
157  if (reprocess) {
158  __end = partition(evt->begin(), __end, JHistory::is_not_event(JMUONSIMPLEX));
159  }
160 
161  if (evt->begin() != __end) {
162 
163  copy(evt->begin(), __end, back_inserter(out));
164 
165  if (numberOfPrefits > 0) {
166  advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
167  }
168 
169  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
170 
171  __end = partition(evt->begin(), __end, JHistory::is_event(evt->begin()->getHistory()));
172 
173 
174  JDataL0_t dataL0;
175  JDataL1_t dataL1;
176 
177  JDAQTimeslice timeslice(*tev, true);
178 
179  for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
180 
181  if (moduleRouter.hasModule(i->getModuleID())) {
182 
183  buffer(*i, moduleRouter.getModule(i->getModuleID()));
184 
185  buildL0(buffer, back_inserter(dataL0));
186  buildL2(buffer, back_inserter(dataL1));
187  }
188  }
189 
190 
191  JDataR1_t data;
192 
193  for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
194 
195  const JRotation3D R (getDirection(*track));
196  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
197  const JModel<JLine1Z> match(tz, roadWidth_m, JTimeRange(-TMAX_NS, +TMAX_NS));
198 
199  data.clear();
200 
201  for (JDataL1_t::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {
202 
203  JHitR1 hit(*i);
204 
205  hit.rotate(R);
206 
207  if (match(hit)) {
208  data.push_back(hit);
209  }
210  }
211 
212  if (useL0) {
213 
214  data.reserve(data.size() + dataL0.size());
215 
216  JDataR1_t::iterator __end1 = data.end();
217 
218  for (JDataL0_t::iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
219 
220  if (find_if(data.begin(), __end1, bind2nd(equal_to<JDAQModuleIdentifier>(), i->getModuleID())) == __end1) {
221 
222  JHitR1 hit(*i);
223 
224  hit.rotate(R);
225 
226  if (match(hit)) {
227  data.push_back(hit);
228  }
229  }
230  }
231  }
232 
233 
234  fit.step.resize(5);
235 
236  fit.step[0] = JLine3Z(JLine1Z(JVector3D(0.5, 0.0, 0.0), 0.0));
237  fit.step[1] = JLine3Z(JLine1Z(JVector3D(0.0, 0.5, 0.0), 0.0));
238  fit.step[2] = JLine3Z(JLine1Z(JVector3D(0.0, 0.0, 0.0), 1.0));
239  fit.step[3] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.005, 0.0));
240  fit.step[4] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.0, 0.005));
241 
242  const int NDF = getCount(data.begin(), data.end()) - fit.step.size();
243 
244  if (NDF >= 0) {
245 
246  const double chi2 = fit(JLine3Z(tz), data.begin(), data.end());
247 
248  JTrack3D tb(fit.value);
249 
250  tb.rotate_back(R);
251 
252  out.push_back(getFit(JHistory(track->getHistory()).add(JMUONSIMPLEX), tb, getQuality(chi2, NDF), NDF));
253 
254  out.rbegin()->setW(track->getW());
255  }
256  }
257 
258  // apply default sorter
259 
260  sort(out.begin(), out.end(), qualitySorter);
261  }
262 
263  outputFile.put(out);
264  }
265  STATUS(endl);
266 
268 
269  io >> outputFile;
270 
271  outputFile.close();
272 }
JFIT::JMEstimatorLorentzian
Lorentzian M-estimator.
Definition: JMEstimator.hh:79
JMeta.hh
Head.hh
main
int main(int argc, char **argv)
Definition: JSimplex.cc:57
KM3NETDAQ::JDAQEvent
DAQ Event.
Definition: JDAQEvent.hh:30
JLine3Z.hh
JTriggerParameters.hh
JFileRecorder.hh
JSuperFrame2D.hh
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JFIT::JLine3Z
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:35
JMessage.hh
JTRIGGER::JL2Parameters
Data structure for L2 parameters.
Definition: JTriggerParameters.hh:33
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
JFitToolkit.hh
JLine3ZRegressor.hh
JGEOMETRY3D::JVersor3Z
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:37
qualitySorter
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
Definition: quality_sorter.cc:10
JDETECTOR::JModuleRouter::getModule
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Definition: JModuleRouter.hh:89
JTimeslice.hh
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:476
JSimplex.hh
std::vector< JHitL0 >
JEvtToolkit.hh
Evt.hh
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:30
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
JFIT::JEvt
Data structure for set of track fit results.
Definition: JEvt.hh:293
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
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
JSupport.hh
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JBuildL1.hh
JEvt.hh
JBuildL0.hh
JDAQEventIO.hh
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition: JVector3D.hh:33
JDAQTimesliceIO.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
JFIT::JModel< JLine1Z >
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JModel.hh:34
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
JMUONSIMPLEX
static const int JMUONSIMPLEX
Definition: reconstruction.hh:15
JHit.hh
JModuleRouter.hh
JFIT::JRegressor< JLine3Z, JSimplex >
Regressor function object for JLine3Z fit using JSimplex minimiser.
Definition: JLine3ZRegressor.hh:52
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:101
JParallelFileScanner.hh
JHitR1.hh
JMultipleFileScanner.hh
JLANG::JTypeList
Type list.
Definition: JTypeList.hh:22
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JBuildL2.hh
JParser.hh
JDetectorToolkit.hh
JFIT::JLine1Z
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JFIT::JHistory::is_not_event
Auxiliary class to test history.
Definition: JHistory.hh:149
JDETECTOR::JModuleRouter
Router for direct addressing of module data in detector data structure.
Definition: JModuleRouter.hh:34
JHitL1.hh
JModel.hh
JTRIGGER::JHitR1
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
JAANET::getDirection
JDirection3D getDirection(const Vec &v)
Get direction.
Definition: JAAnetToolkit.hh:221
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JSUPPORT::JParallelFileScanner
General purpose class for parallel reading of objects from a single file or multiple files.
Definition: JParallelFileScanner.hh:31
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JAANET::getPosition
JPosition3D getPosition(const Vec &v)
Get position.
Definition: JAAnetToolkit.hh:197
JAANET::detector
Detector file.
Definition: JHead.hh:130
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
JDETECTOR::JModuleRouter::hasModule
bool hasModule(const JObjectID &id) const
Has module.
Definition: JModuleRouter.hh:101
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JGEOMETRY3D::JTrack3D
3D track.
Definition: JTrack3D.hh:30
JDAQSummarysliceIO.hh
JFIT::getFit
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=0)
Get fit.
Definition: JEvtToolkit.hh:116
JTRIGGER::JBuildL0< JHitL0 >
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:101
JFIT::JHistory
Container for historical events.
Definition: JHistory.hh:95
JDetector.hh
JHitL0.hh
JTRIGGER::JSuperFrame2D
2-dimensional frame with time calibrated data from one optical module.
Definition: JSuperFrame2D.hh:41
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JSUPPORT::JFileRecorder
Object writing to file.
Definition: JFileRecorder.hh:41
JFIT::getQuality
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
JLANG::JException
General exception.
Definition: JException.hh:23
JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition: JRotation3D.hh:111
JSUPPORT::JSingleFileScanner
Object reading from a list of files.
Definition: JSingleFileScanner.hh:75
JFrame.hh
JFIT::getCount
int getCount(const JHitL0 &hit)
Get hit count.
Definition: JEvtToolkit.hh:728
JTRIGGER::JBuildL2
Template L2 builder.
Definition: JBuildL2.hh:45