Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonStart.hh
Go to the documentation of this file.
1 #ifndef __JRECONSTRUCTION__JMUONSTART__
2 #define __JRECONSTRUCTION__JMUONSTART__
3 
4 #include <string>
5 #include <istream>
6 #include <ostream>
7 #include <iomanip>
8 #include <vector>
9 #include <set>
10 #include <map>
11 #include <algorithm>
12 #include <cmath>
13 
17 
18 #include "JDetector/JDetector.hh"
21 
22 #include "JTrigger/JHitL0.hh"
23 #include "JTrigger/JHitR1.hh"
24 #include "JTrigger/JBuildL0.hh"
25 
27 
28 #include "JFit/JLine1Z.hh"
29 #include "JFit/JEnergyRegressor.hh"
30 #include "JFit/JFitToolkit.hh"
31 #include "JFit/JModel.hh"
32 #include "JFit/JNPEHit.hh"
33 #include "JFit/JEnergy.hh"
34 #include "JFit/JMEstimator.hh"
35 
36 #include "JReconstruction/JEvt.hh"
40 
41 #include "JPhysics/JConstants.hh"
42 #include "JPhysics/JK40Rates.hh"
43 
44 #include "JLang/JComparator.hh"
45 #include "JLang/JPredicate.hh"
46 
47 
48 /**
49  * \author mdejong, azegarelli, scelli
50  */
51 namespace JRECONSTRUCTION {}
52 namespace JPP { using namespace JRECONSTRUCTION; }
53 
54 namespace JRECONSTRUCTION {
55 
57  using JDETECTOR::JModule;
59  using JFIT::JRegressor;
60  using JFIT::JEnergy;
61  using JPHYSICS::JK40Rates;
62  using JPHYSICS::JRateL1_t;
64 
65  /**
66  * Get probability of given response in optical module due to random background.
67  *
68  * \param module module
69  * \param frame summary frame
70  * \param R_Hz multiples rates [Hz]
71  * \param T_ns time window [ns]
72  * \param top list of indentifiers of PMTs with hit in given module
73  * \return probability
74  */
75  inline double getProbability(const JModule& module,
76  const JDAQSummaryFrame& frame,
77  const JRateL1_t& R_Hz,
78  const double T_ns,
79  const std::multiset<int>& top)
80  {
81  using namespace std;
82  using namespace JPP;
83  using namespace KM3NETDAQ;
84 
85  size_t N = 0; // number of active PMTs
86  size_t M = 0; // multiplicity
87  double R = 0.0; // total rate
88 
89  for (size_t i = 0; i != module.size(); ++i) {
90 
91  if (getDAQStatus(frame, module, i) &&
92  getPMTStatus(frame, module, i) &&
93  frame[i].is_valid() &&
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  * Auxiliary class to determine start and end position of muon trajectory.
111  */
112  class JMuonStart :
113  public JMuonStartParameters_t,
114  public JRegressor<JEnergy>
115  {
116  public:
120 
121  using JRegressor_t::operator();
122 
123  /**
124  * Constructor.
125  *
126  * \param parameters parameters
127  * \param router module router
128  * \param summary summary router
129  * \param pdfFile PDF file
130  * \param rates_Hz K40 rates [Hz]
131  * \param debug debug
132  */
134  const JModuleRouter& router,
135  const JSummaryRouter& summary,
136  const std::string& pdfFile,
137  const JK40Rates& rates_Hz,
138  const int debug = 0):
139  JMuonStartParameters_t(parameters),
140  JRegressor_t(pdfFile),
141  router(router),
142  summary(summary),
143  rates_Hz(rates_Hz),
144  debug(debug)
145  {
147  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
148 
149  if (this->getRmax() < parameters.roadWidth_m) {
150  roadWidth_m = this->getRmax();
151  }
152  }
153 
154 
155  /**
156  * Fit function.
157  *
158  * \param event event
159  * \param in start values
160  * \return fit results
161  */
163  {
164  using namespace std;
165  using namespace JPP;
166 
167  const JBuildL0<hit_type> buildL0;
168 
169  buffer_type dataL0;
170 
171  buildL0(event, router, true, back_inserter(dataL0));
172 
173  return (*this)(dataL0, in);
174  }
175 
176 
177  /**
178  * Fit function.
179  *
180  * \param data dataL0
181  * \param in start values
182  * \return fit results
183  */
185  {
186  using namespace std;
187  using namespace JPP;
188 
189  JEvt out;
190 
191  const JDetector& detector = router.getReference();
192 
193  const JStart start(Pmin1, Pmin2, Nmax2);
194 
195  for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
196 
197  const JRotation3D R (getDirection(*track));
198  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
199  const JFIT::JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
200 
202 
203  for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
204 
205  hit_type hit(*i);
206 
207  hit.rotate(R);
208 
209  if (match(hit)) {
210  top[hit.getModuleID()].insert(hit.getPMTAddress());
211  }
212  }
213 
214  struct JHit_t {
215 
216  double getZ() const { return z; }
217  double getP() const { return p; }
218 
219  double z;
220  double p;
221  };
222 
224 
225  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
226 
227  if (summary.hasSummaryFrame(module->getID())) {
228 
229  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
230 
231  JPosition3D pos(module->getPosition());
232 
233  pos.transform(R, tz.getPosition());
234 
235  if (pos.getX() <= roadWidth_m) {
236 
237  const double z = pos.getZ() - pos.getX() / getTanThetaC();
238 
239  const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), JRegressor_t::T_ns.getLength(), top[module->getID()]);
240 
241  data.push_back({ z, p });
242  }
243  }
244  }
245 
246  double Zmin = 0.0; // relative start position [m]
247  double Zmax = 0.0; // relative end position [m]
248  double npe_total = 0.0; // npe due to mip
249  double npe_missed = 0.0; // npe due to mip not detected
250 
251  if (!data.empty()) {
252 
253  sort(data.begin(), data.end(), JLANG::make_comparator(&JHit_t::getZ));
254 
255  vector<JHit_t>::const_iterator track_start = start.find(data. begin(), data. end());
256  vector<JHit_t>::const_reverse_iterator track_end = start.find(data.rbegin(), data.rend());
257 
258  if (track_start != data. end()) { Zmin = track_start->getZ(); }
259  if (track_end != data.rend()) { Zmax = track_end ->getZ(); }
260 
261  // additional features
262 
263  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
264 
265  JPosition3D pos(module->getPosition());
266 
267  pos.transform(R, tz.getPosition());
268 
269  if (pos.getX() <= roadWidth_m) {
270 
271  const double z = pos.getZ() - pos.getX() / getTanThetaC();
272 
273  if (z >= Zmin && z <= Zmax) {
274 
275  for (size_t i = 0; i != module->size(); ++i) {
276 
277  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
278 
279  if (getDAQStatus(frame, *module, i) &&
280  getPMTStatus(frame, *module, i) &&
281  frame[i].is_valid() &&
282  !module->getPMT(i).has(PMT_DISABLE)) {
283 
284  JPMT pmt = module->getPMT(i);
285 
286  pmt.transform(R, tz.getPosition());
287 
288  const double ya = this->getNPE(pmt, 0.0).getYA();
289 
290  npe_total += ya;
291 
292  if (top[module->getID()].count(i) == 0) {
293  npe_missed += ya;
294  }
295  }
296  }
297  }
298  }
299  }
300  }
301 
302  JFit fit = *track;
303 
304  // move track
305 
306  fit.move(Zmin, getSpeedOfLight());
307 
308  out.push_back(fit.add(JMUONSTART));
309 
310  out.rbegin()->setW(track->getW());
311  out.rbegin()->setW(JSTART_NPE_MIP_TOTAL, npe_total);
312  out.rbegin()->setW(JSTART_NPE_MIP_MISSED, npe_missed);
313  out.rbegin()->setW(JSTART_LENGTH_METRES, Zmax - Zmin);
314  }
315 
316  return out;
317  }
318 
319  private:
323  int debug;
324  };
325 }
326 
327 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
static const int JMUONSTART
JRegressor< JEnergy > JRegressor_t
Definition: JMuonStart.hh:117
Template definition of a data regressor of given model.
Definition: JRegressor.hh:70
static const int JSTART_NPE_MIP_TOTAL
number of photo-electrons along the whole track from JStart.cc
int getModuleID() const
Get module identifier.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Data structure for a composite optical module.
Definition: JModule.hh:67
Auxiliary method to locate start and end point of muon trajectory.
Regressor function object for fit of muon energy.
Detector data structure.
Definition: JDetector.hh:89
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
T find(T __begin, T __end) const
Get start point of muon trajectory.
Definition: JStart.hh:77
JFit & add(const int type)
Add event to history.
Data structure for detector geometry and calibration.
void transform(const JRotation3D &R, const JVector3D &pos)
Transform position.
Definition: JPosition3D.hh:331
void transform(const JAxis3D &axis)
Transform axis to reference frame of given axis.
Definition: JAxis3D.hh:359
int getPMTAddress() const
Get PMT identifier.
Basic data structure for L0 hit.
Data structure for track fit results with history and optional associated values. ...
Auxiliary class to extract a subset of optical modules from a detector.
const JK40Rates rates_Hz
Definition: JMuonStart.hh:322
Data structure for fit parameters.
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Auxiliary class to determine start and end position of muon trajectory.
Definition: JMuonStart.hh:112
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:225
Detector file.
Definition: JHead.hh:226
JDirection3D getDirection(const Vec &dir)
Get direction.
JMuonStart(const JMuonStartParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string &pdfFile, const JK40Rates &rates_Hz, const int debug=0)
Constructor.
Definition: JMuonStart.hh:133
static const int JSTART_NPE_MIP_MISSED
number of photo-electrons missed from JStart.cc
bool has(const int bit) const
Test PMT status.
Definition: JStatus.hh:120
Auxiliary class for start or end point evaluation.
Definition: JStart.hh:21
Data storage class for rate measurements of all PMTs in one module.
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
static const int PMT_DISABLE
KM3NeT Data Definitions v3.3.0-2-g5cc95cf https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
JPosition3D getPosition(const Vec &pos)
Get position.
Physics constants.
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JEvt &in)
Fit function.
Definition: JMuonStart.hh:162
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
Direct access to module in detector data structure.
const JModuleRouter & router
Definition: JMuonStart.hh:320
Reduced data structure for L1 hit.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
bool getPMTStatus(const JStatus &status)
Test status of PMT.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
JTRIGGER::JHitL0 hit_type
Definition: JMuonStart.hh:118
const double getSpeedOfLight()
Get speed of light.
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
Data structure for set of track fit results.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
Data structure for L0 hit.
Definition: JHitL0.hh:27
Auxiliary class to set-up Hit.
Definition: JSirene.hh:57
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
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
const JSummaryRouter & summary
Definition: JMuonStart.hh:321
Data structure for fit of energy.
Definition: JEnergy.hh:28
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
double getNPE(const Hit &hit)
Get true charge of hit.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
Data regression method for JFIT::JEnergy.
void move(const double step, const double velocity)
Move vertex along this track with given velocity.
Template L0 hit builder.
Definition: JBuildL0.hh:35
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
JEvt operator()(const buffer_type &data, const JEvt &in)
Fit function.
Definition: JMuonStart.hh:184
std::vector< double > JRateL1_t
Type definition of count rate as a function of multiplicty [Hz] The multiples rate start counting at ...
Definition: JK40Rates.hh:27
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
std::vector< hit_type > buffer_type
Definition: JMuonStart.hh:119
Maximum likelihood estimator (M-estimators).
int debug
debug level
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41