Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JSimplex.cc File Reference

Program to perform JFIT::JRegressor<JLine3Z,JSimplex> fit with I/O of JFIT::JEvt data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JTools/JConstants.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JModel.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 perform JFIT::JRegressor<JLine3Z,JSimplex> fit with I/O of JFIT::JEvt data.

Author
mdejong

Definition in file JSimplex.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 57 of file JSimplex.cc.

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;
68  JFileRecorder<typelist> outputFile;
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 
110  JDetector detector;
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(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  }
255 
256  // apply default sorter
257 
258  sort(out.begin(), out.end(), qualitySorter);
259  }
260 
261  outputFile.put(out);
262  }
263  STATUS(endl);
264 
265  JSingleFileScanner<JRemove<typelist, JEvt>::typelist> io(inputFile);
266 
267  io >> outputFile;
268 
269  outputFile.close();
270 }
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
double getCount(TH1D *hptr, int muon_threshold)
#define STATUS(A)
Definition: JMessage.hh:61
string outputFile
JRange< double > JTimeRange
Type definition for time range.
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
JLimit JLimit_t
Type definition of limit.
Definition: JLimit.hh:214
Type list.
Definition: JTypeList.hh:22
const_iterator< T > end() const
Get end of data.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Data time slice.
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
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
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
bool qualitySorter(const JFIT::JFit &first, const JFIT::JFit &second)
Comparison of fit results.
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
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JPosition3D getPosition(const Vec &v)
Get position.