Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JGandalf.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 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH2D.h"
12 
13 #include "evt/Head.hh"
14 #include "evt/Evt.hh"
15 #include "JDAQ/JDAQEvent.hh"
16 #include "JDAQ/JDAQTimeslice.hh"
17 #include "JDAQ/JDAQSummaryslice.hh"
18 
19 #include "JDetector/JDetector.hh"
22 
23 #include "JLang/JPredicate.hh"
24 
25 #include "JTrigger/JHit.hh"
26 #include "JTrigger/JFrame.hh"
27 #include "JTrigger/JTimeslice.hh"
28 #include "JTrigger/JHitL0.hh"
29 #include "JTrigger/JBuildL0.hh"
31 
35 #include "JSupport/JSupport.hh"
36 #include "JSupport/JMeta.hh"
37 
38 #include "JFit/JHitW0.hh"
39 #include "JFit/JPMTW0.hh"
40 #include "JFit/JLine3Z.hh"
41 #include "JFit/JGandalf.hh"
42 #include "JFit/JLine3ZRegressor.hh"
43 #include "JFit/JFitToolkit.hh"
44 #include "JFit/JEvt.hh"
45 #include "JFit/JEvtToolkit.hh"
46 #include "JFit/JRegressor.hh"
47 #include "JFit/JFitParameters.hh"
48 #include "JFit/JModel.hh"
49 
50 #include "JTools/JConstants.hh"
51 #include "JTools/JFunction1D_t.hh"
53 
54 #include "JPhysics/JPDFTable.hh"
55 #include "JPhysics/JPDFToolkit.hh"
56 
57 #include "Jeep/JParser.hh"
58 #include "Jeep/JMessage.hh"
59 
60 
61 namespace {
62 
63  using namespace JPP;
64 
65  /**
66  * Compare hits by PMT identifier and time.
67  *
68  * \param first first hit
69  * \param second second hit
70  * \return true if first before second; else false
71  */
72  inline bool compare(const JHitL0& first, const JHitL0& second)
73  {
74  if (std::equal_to<KM3NETDAQ::JDAQPMTIdentifier>()(first, second))
75  return std::less<JTRIGGER::JHit>()(first, second);
76  else
77  return std::less<KM3NETDAQ::JDAQPMTIdentifier>()(first, second);
78  }
79 }
80 
81 
82 /**
83  * \file
84  *
85  * Program to perform JFIT::JRegressor<JLine3Z,JGandalf> fit with I/O of JFIT::JEvt data.
86  * \author mdejong
87  */
88 int main(int argc, char **argv)
89 {
90  using namespace std;
91  using namespace JPP;
92  using namespace KM3NETDAQ;
93 
94  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> > JParallelFileScanner_t;
95  typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
97 
98  JParallelFileScanner_t inputFile;
100  JLimit_t& numberOfEvents = inputFile.getLimit();
101  string detectorFile;
102  string pdfFile;
103  double roadWidth_m;
104  double R_Hz;
105  size_t numberOfPrefits;
106  double TTS_ns;
107  double E_GeV;
108  bool reprocess;
109  bool calibrate;
110  int debug;
111 
112  try {
113 
114  JParser<> zap("Program to perform PDF fit of muon trajectory to data.");
115 
116  zap['f'] = make_field(inputFile);
117  zap['o'] = make_field(outputFile) = "gandalf.root";
118  zap['a'] = make_field(detectorFile);
119  zap['n'] = make_field(numberOfEvents) = JLimit::max();
120  zap['P'] = make_field(pdfFile);
121  zap['R'] = make_field(roadWidth_m) = numeric_limits<double>::max();
122  zap['B'] = make_field(R_Hz) = 6.0e3;
123  zap['N'] = make_field(numberOfPrefits) = 0;
124  zap['T'] = make_field(TTS_ns);
125  zap['E'] = make_field(E_GeV) = 1.0e3;
126  zap['r'] = make_field(reprocess);
127  zap['c'] = make_field(calibrate, "Enable time calibration per optical module.");
128  zap['d'] = make_field(debug) = 1;
129 
130  zap(argc, argv);
131  }
132  catch(const exception& error) {
133  FATAL(error.what() << endl);
134  }
135 
136 
137 
138  cout.tie(&cerr);
139 
140 
142 
143  try {
144  load(detectorFile, detector);
145  }
146  catch(const JException& error) {
147  FATAL(error);
148  }
149 
150  const JModuleRouter moduleRouter(detector);
151 
152 
153  typedef vector<JHitL0> JDataL0_t;
154  typedef vector<JHitW0> JDataW0_t;
155  const JBuildL0<JHitL0> buildL0;
156 
157 
158  typedef JRegressor<JLine3Z, JGandalf> JRegressor_t;
159 
161  JRegressor_t::T_ns.setRange(-50.0, +450.0); // ns
162  JRegressor_t::Vmax_npe = 10.0; // npe
163  JRegressor_t::MAXIMUM_ITERATIONS = 1000;
164 
165  JRegressor_t fit(pdfFile, TTS_ns);
166 
167  fit.parameters.resize(5);
168 
169  fit.parameters[0] = JLine3Z::pX();
170  fit.parameters[1] = JLine3Z::pY();
171  fit.parameters[2] = JLine3Z::pT();
172  fit.parameters[3] = JLine3Z::pDX();
173  fit.parameters[4] = JLine3Z::pDY();
174 
175  if (fit.getRmax() < roadWidth_m) {
176  roadWidth_m = fit.getRmax();
177  }
178 
179  TH2D h2("h2", NULL, detector.size(), -0.5, detector.size() - 0.5, 1001, -250.0, +250.0);
180 
181  h2.Sumw2();
182 
183  outputFile.open();
184 
185  outputFile.put(JMeta(argc, argv));
186 
187  while (inputFile.hasNext()) {
188 
189  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
190 
191  multi_pointer_type ps = inputFile.next();
192 
193  JDAQEvent* tev = ps;
194  JEvt* evt = ps;
195  JEvt out;
196 
197  JEvt::iterator __end = evt->end();
198 
199  if (reprocess) {
200  __end = partition(evt->begin(), __end, JHistory::is_not_event(JMUONGANDALF));
201  }
202 
203  if (evt->begin() != __end) {
204 
205  copy(evt->begin(), __end, back_inserter(out));
206 
207  if (numberOfPrefits > 0) {
208  advance(__end = evt->begin(), min(numberOfPrefits, out.size()));
209  }
210 
211  partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
212 
213  __end = partition(evt->begin(), __end, JHistory::is_event(evt->begin()->getHistory()));
214 
215 
216  JDataL0_t dataL0;
217 
218  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
219 
220 
221  if (dataL0.size() >= fit.parameters.size()) {
222 
223  for (JEvt::const_iterator track = evt->begin(); track != __end; ++track) {
224 
225  const JRotation3D R (getDirection(*track));
226  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
227  const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
228 
229  // hit selection based on start value
230 
231  JDataW0_t data;
232 
233  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
234 
235  JHitW0 hit(*i, R_Hz);
236 
237  hit.rotate(R);
238 
239  if (match(hit)) {
240  data.push_back(hit);
241  }
242  }
243 
244  // select first hit
245 
246  sort(data.begin(), data.end(), compare);
247 
248  JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
249 
250 
251  const int NDF = distance(data.begin(), __end) - fit.parameters.size();
252 
253  if (NDF >= 0) {
254 
255  // set fit parameters
256 
257  if (track->getE() > 0.1)
258  fit.E_GeV = track->getE();
259  else
260  fit.E_GeV = E_GeV;
261 
262  const double chi2 = fit(JLine3Z(tz), data.begin(), __end);
263 
264  JTrack3D tb(fit.value);
265 
266  tb.rotate_back(R);
267 
268  out.push_back(getFit(JHistory(track->getHistory()).add(JMUONGANDALF), tb, getQuality(chi2), NDF));
269 
270  // set additional values
271 
272  out.rbegin()->setW(JGANDALF_BETA0_RAD, sqrt(fit.error.getDX() * fit.error.getDX() +
273  fit.error.getDY() * fit.error.getDY()));
274  out.rbegin()->setW(JGANDALF_BETA1_RAD, sqrt(fit.error.getDX() * fit.error.getDY()));
275  out.rbegin()->setW(JGANDALF_CHI2, chi2);
276  out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS, distance(data.begin(), __end));
277  out.rbegin()->setW(JGANDALF_LAMBDA, fit.lambda);
278  out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, fit.numberOfIterations);
279 
280 
281  if (calibrate) {
282 
283  set<int> buffer;
284 
285  for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
286  buffer.insert(hit->getModuleID());
287  }
288 
289  const JLine3Z ta = fit.value;
290 
291  for (set<int>::const_iterator id = buffer.begin(); id != buffer.end(); ++id) {
292 
293  JDataW0_t::iterator q = partition(data.begin(), __end, make_predicate(&JHitW0::getModuleID, *id, JComparison::ne()));
294 
295  if (distance(data.begin(), q) - fit.parameters.size() > 0) {
296 
297  fit(ta, data.begin(), q);
298 
299  sort(q, __end, compare); // select first hit in module
300 
301  h2.Fill(moduleRouter.getIndex(*id), q->getT() - fit.value.getT(q->getPosition()));
302  }
303  }
304  }
305  }
306  }
307  }
308 
309  // apply default sorter
310 
311  sort(out.begin(), out.end(), qualitySorter);
312  }
313 
314  outputFile.put(out);
315  }
316  STATUS(endl);
317 
318  if (calibrate) {
319  outputFile.put(h2);
320  }
321 
323 
324  io >> outputFile;
325 
326  outputFile.close();
327 }
Auxiliary methods to evaluate Poisson probabilities and chi2.
Data regression method for JFIT::JLine3Z.
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:71
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1410
General exception.
Definition: JException.hh:40
General purpose data regression method.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
Auxiliary methods for PDF calculations.
number of iterations from JGandalf.cc
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
#define STATUS(A)
Definition: JMessage.hh:61
Detector data structure.
Definition: JDetector.hh:77
const int getIndex(const JObjectID &id) const
Get index of module.
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:108
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: JModel.hh:31
Auxiliary class to test history.
Definition: JHistory.hh:149
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:35
Container for historical events.
Definition: JHistory.hh:95
chi2 from JGandalf.cc
string outputFile
Data structure for detector geometry and calibration.
angular resolution [rad] from JGandalf.cc
Various implementations of functional maps.
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
Basic data structure for L0 hit.
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 for defining the range of iterations of objects.
Definition: JLimit.hh:41
Constants.
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:235
Detector file.
Definition: JHead.hh:126
control parameter from JGandalf.cc
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1836
Auxiliary class to test history.
Definition: JHistory.hh:101
Parallel scanning of objects from a single file or multiple files according a format that follows fro...
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
ROOT I/O of application specific meta data.
angular resolution [rad] from JGandalf.cc
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
int debug
debug level
Definition: JSirene.cc:59
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 hits from JGandalf.cc
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
Data structure for set of track fit results.
Definition: JEvt.hh:312
Regressor function object for JLine3Z fit using JGandalf minimiser.
Utility class to parse command line options.
Data structure for L0 hit.
Definition: JHitL0.hh:27
ROOT TTree parameter settings.
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
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.
Object reading from a list of files.
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.
int main(int argc, char *argv[])
Definition: Main.cpp:15