Jpp  debug
the software that should make you happy
examples/JReconstruction/JMuonStart.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 #include <set>
6 #include <map>
7 #include <algorithm>
8 
9 #include "TROOT.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 
19 
20 #include "JDAQ/JDAQEventIO.hh"
22 
23 #include "JDetector/JDetector.hh"
26 #include "JDetector/JPMTRouter.hh"
27 
28 #include "JTrigger/JHitL0.hh"
29 #include "JTrigger/JBuildL0.hh"
31 
32 #include "JSupport/JSupport.hh"
36 #include "JAAnet/JAAnetToolkit.hh"
37 
38 #include "JFit/JLine1Z.hh"
39 #include "JFit/JNPEHit.hh"
40 #include "JFit/JModel.hh"
41 
42 #include "JPhysics/KM3NeT.hh"
43 #include "JPhysics/JGeane.hh"
44 
46 #include "JReconstruction/JEvt.hh"
50 
51 #include "JROOT/JManager.hh"
53 
54 #include "JLang/JPredicate.hh"
55 #include "JLang/JComparator.hh"
56 
57 #include "Jeep/JPrint.hh"
58 #include "Jeep/JParser.hh"
59 #include "Jeep/JMessage.hh"
60 
61 namespace {
62 
63  using namespace JPP;
64  using namespace KM3NETDAQ;
65 
66  /**
67  * Get probability of given response in optical module due to random background.
68  *
69  * \param module module
70  * \param frame summary frame
71  * \param R_Hz multiples rates [Hz]
72  * \param T_ns time window [ns]
73  * \param top list of indentifiers of PMTs with hit in given module
74  * \return probability
75  */
76  inline double getProbability(const JModule& module,
77  const JDAQSummaryFrame& frame,
78  const JRateL1_t& R_Hz,
79  const double T_ns,
80  const std::multiset<int>& top)
81  {
82  using namespace std;
83  using namespace JPP;
84  using namespace KM3NETDAQ;
85 
86  size_t N = 0; // number of active PMTs
87  size_t M = 0; // multiplicity
88  double R = 0.0; // total rate
89 
90  for (size_t i = 0; i != module.size(); ++i) {
91 
92  if (getDAQStatus(frame, module, i) &&
93  getPMTStatus(frame, module, i) &&
94  !module.getPMT(i).has(PMT_DISABLE)) {
95 
96  N += 1;
97  M += top.count(i);
98  R += frame.getRate(i);
99  }
100  }
101 
102  if (N > 0)
103  return JFIT::getProbability(N, M, JK40Rates(R/N, R_Hz), T_ns);
104  else
105  return (M == 0 ? 1.0 : 0.0);
106  }
107 }
108 
109 
110 /**
111  * \file
112  *
113  * Program to test JMuonStart.cc.
114  *
115  * \author mdejong
116  */
117 int main(int argc, char **argv)
118 {
119  using namespace std;
120  using namespace JPP;
121  using namespace KM3NETDAQ;
122 
123  typedef JTriggeredFileScanner<JEvt, JSingleFileScanner> JTriggeredFileScanner_t;
124  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
125 
127 
128  JTriggeredFileScanner_t inputFile;
129  JLimit_t& numberOfEvents = inputFile.getLimit();
130  string detectorFile;
131  string outputFile;
132  JMuonStartParameters_t parameters;
133  JK40Rates rates_Hz;
134  JHistogram_t histogram;
135  int debug;
136 
137  try {
138 
139  parameters.numberOfPrefits = 1;
140 
141  JParser<> zap("Program to test JMuonStart.");
142 
143  zap['f'] = make_field(inputFile, "input file (output of JXXXMuonReconstruction.sh)");
144  zap['a'] = make_field(detectorFile);
145  zap['n'] = make_field(numberOfEvents) = JLimit::max();
146  zap['o'] = make_field(outputFile) = "start.root";
147  zap['@'] = make_field(parameters) = JPARSER::initialised();
148  zap['B'] = make_field(rates_Hz) = KM3NET::getK40Rates();
149  zap['H'] = make_field(histogram, "histogram binning") = JHistogram_t(200, -250.0, +250.0);
150  zap['d'] = make_field(debug) = 1;
151 
152  zap(argc, argv);
153  }
154  catch(const exception& error) {
155  FATAL(error.what() << endl);
156  }
157 
158 
160 
161  try {
162  load(detectorFile, detector);
163  }
164  catch(const JException& error) {
165  FATAL(error);
166  }
167 
168  const JModuleRouter moduleRouter(detector);
169  const JPMTRouter pmtRouter (detector);
170 
171  JSummaryFileRouter summary(inputFile, parameters.R_Hz);
172 
173  const JTimeRange T_ns (parameters.TMin_ns, parameters.TMax_ns);
174  const JStart start(parameters.Pmin1, parameters.Pmin2, parameters.Nmax2);
175 
176  typedef JHitL0 hit_type;
177 
179  const JBuildL0<hit_type> buildL0;
180 
181 
182  JHead head;
183 
184  try {
185  head = getHeader(inputFile);
186  }
187  catch(const JException& error) {
188  FATAL(error);
189  }
190 
191  const Vec offset = getOffset(head);
192 
193  JManager<string, TH1D> H1(new TH1D("[%].1L", NULL,
194  histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit()));
195  JManager<string, TH2D> H2(new TH2D("[%].2D", NULL,
196  histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit(),
197  histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit()));
198 
199  while (inputFile.hasNext()) {
200 
201  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
202 
203  multi_pointer_type ps = inputFile.next();
204 
205  JDAQEvent* tev = ps;
206  JEvt* evt = ps;
207  Evt* event = ps;
208 
209  summary.update(*tev);
210 
211  vector<Trk>::iterator muon = event->mc_trks.end();
212 
213  for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
214 
215  if (is_muon(*track)) {
216  if (muon == event->mc_trks.end() || track->E > muon->E) {
217  muon = track;
218  }
219  }
220  }
221 
222  if (muon == event->mc_trks.end()) {
223  continue;
224  }
225 
226  if (muon->len == 0.0) {
227  muon->len = gWater(muon->E);
228  }
229 
230 
231  typedef JRange<double> JRange_t;
232 
233  JRange_t za(JRange_t::DEFAULT_RANGE()); // Monte Carlo
234  JRange_t zb(JRange_t::DEFAULT_RANGE()); // original result
235  JRange_t zd(JRange_t::DEFAULT_RANGE()); // new
236 
237  {
238  const JRotation3D R(getDirection(*muon));
239 
240  JTrack3D ta = getTrack(*muon);
241 
242  ta.add(getPosition(offset));
243  ta.rotate(R);
244 
245  za.setLowerLimit(ta.getZ());
246  za.setLength(fabs(muon->len));
247  }
248 
249 
250  buffer_type dataL0;
251 
252  buildL0(*tev, moduleRouter, true, back_inserter(dataL0));
253 
254 
255  if (!evt->empty()) {
256 
257  sort(evt->begin(), evt->end(), qualitySorter);
258 
259  JEvt::const_iterator track = evt->begin();
260 
261  const JRotation3D R (getDirection(*track));
262  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
263  const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns);
264 
265  if (track->has(JMUONSTART)) {
266  if (track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0) {
267  zb = JRange_t(tz.getZ(), tz.getZ() + track->getW(JSTART_LENGTH_METRES));
268  }
269  }
270 
271  {
273 
274  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
275 
276  hit_type hit(*i);
277 
278  hit.rotate(R);
279 
280  if (match(hit)) {
281  top[hit.getModuleID()].insert(hit.getPMTAddress());
282  }
283  }
284 
285  struct JHit_t {
286 
287  double getZ() const { return z; }
288  double getP() const { return p; }
289 
290  double z;
291  double p;
292  };
293 
295 
296  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
297 
298  if (summary.hasSummaryFrame(module->getID())) {
299 
300  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
301 
302  JPosition3D pos(module->getPosition());
303 
304  pos.transform(R, tz.getPosition());
305 
306  if (pos.getX() <= parameters.roadWidth_m) {
307 
308  const double z = pos.getZ() - pos.getX() / getTanThetaC();
309  const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), T_ns.getLength(), top[module->getID()]);
310 
311  data.push_back({ z, p });
312  }
313  }
314  }
315 
316  if (!data.empty()) {
317 
318  sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
319 
320  vector<JHit_t>::const_iterator track_start = start.find(data. begin(), data. end());
321  vector<JHit_t>::const_reverse_iterator track_end = start.find(data.rbegin(), data.rend());
322 
323  if (track_start == data. end()) { track_start = data. begin(); }
324  if (track_end == data.rend()) { track_end = data.rbegin(); }
325 
326  zd = JRange_t(track_start->getZ(), track_end->getZ());
327 
328  zd.add(tz.getZ());
329  }
330  }
331  }
332 
333  H1["B"]->Fill(zb.getLength() - za.getLength());
334  H1["D"]->Fill(zd.getLength() - za.getLength());
335 
336  H2["B"]->Fill(zb.getLowerLimit() - za.getLowerLimit(), zb.getUpperLimit() - za.getUpperLimit());
337  H2["D"]->Fill(zd.getLowerLimit() - za.getLowerLimit(), zd.getUpperLimit() - za.getUpperLimit());
338  }
339  STATUS(endl);
340 
341 
342  TFile out(outputFile.c_str(), "recreate");
343 
344  out << H1 << H2;
345 
346  out.Write();
347  out.Close();
348 }
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
KM3NeT DAQ constants, bit handling, etc.
Data structure for detector geometry and calibration.
Energy loss of muon.
Basic data structure for L0 hit.
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define STATUS(A)
Definition: JMessage.hh:63
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
Direct access to module in detector data structure.
Direct access to PMT in detector data structure.
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
I/O formatting auxiliaries.
Auxiliary method to locate start and end point of muon trajectory.
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
Properties of KM3NeT PMT and deep-sea water.
Monte Carlo run header.
Definition: JHead.hh:1236
Detector data structure.
Definition: JDetector.hh:96
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition: JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:37
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:29
double getZ(const JPosition3D &pos) const
Get point of emission of Cherenkov light along muon path.
Definition: JLine1Z.hh:134
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:225
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
Rotation matrix.
Definition: JRotation3D.hh:114
JTime & add(const JTime &value)
Addition operator.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getX() const
Get x position.
Definition: JVector3D.hh:94
General exception.
Definition: JException.hh:24
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Utility class to parse command line options.
Definition: JParser.hh:1714
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:304
File router for fast addressing of summary data.
void update(const JDAQHeader &header)
Update router.
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
void setLength(argument_type length)
Set length (difference between upper and lower limit).
Definition: JRange.hh:300
T getLength() const
Get length (difference between upper and lower limit).
Definition: JRange.hh:289
range_type & add(argument_type x)
Add offset.
Definition: JRange.hh:446
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
void setLowerLimit(argument_type x)
Set lower limit.
Definition: JRange.hh:224
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Template L0 hit builder.
Definition: JBuildL0.hh:38
Data structure for L0 hit.
Definition: JHitL0.hh:31
Data storage class for rate measurements of all PMTs in one module.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
static const int JMUONSTART
int main(int argc, char **argv)
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
JDirection3D getDirection(const Vec &dir)
Get direction.
JTrack3E getTrack(const Trk &track)
Get track.
JPosition3D getPosition(const Vec &pos)
Get position.
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Vec getOffset(const JHead &header)
Get offset.
JAbstractHistogram< double > JHistogram_t
Type definition for scan along axis.
Definition: JBillabong.cc:61
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double getProbability(const size_t N, const size_t M, const JK40Rates &R_Hz, const double T_ns)
Get probability due to random background.
Definition: JFitToolkit.hh:248
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits.
Definition: JFitToolkit.hh:41
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
std::vector< JHitW0 > buffer_type
hits
Definition: JPerth.cc:70
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
bool getPMTStatus(const JStatus &status)
Test status of PMT.
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
const JK40Rates & getK40Rates()
Get K40 rates.
Definition: KM3NeT.hh:36
Definition: JSTDTypes.hh:14
static const int PMT_DISABLE
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:21
Type definition of range.
Definition: JHead.hh:43
Detector file.
Definition: JHead.hh:227
Acoustic event fit.
bool has(const int bit) const
Test PMT status.
Definition: JStatus.hh:120
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:36
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
const JRateL1_t & getMultiplesRates() const
Get multiples rate.
Definition: JK40Rates.hh:82
Data structure for fit parameters.
int Nmax2
maximal number for twofold observations
double Pmin1
minimal probability single observation
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double Pmin2
minimal probability for twofold observations
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Auxiliary class for start or end point evaluation.
Definition: JStart.hh:21
T find(T __begin, T __end) const
Get start point of muon trajectory.
Definition: JStart.hh:77
Auxiliary class to set-up Hit.
Definition: JSirene.hh:58
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition: JLimit.hh:128
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Simple data structure for histogram binning.
int getNumberOfBins() const
Get number of bins.
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:13