Jpp  18.0.0-rc.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonSimplex.hh
Go to the documentation of this file.
1 #ifndef __JRECONSTRUCTION__JMUONSIMPLEX__
2 #define __JRECONSTRUCTION__JMUONSIMPLEX__
3 
4 #include <iostream>
5 #include <iomanip>
6 #include <vector>
7 #include <algorithm>
8 
12 
13 #include "JTrigger/JHitR1.hh"
15 #include "JTrigger/JBuildL0.hh"
16 #include "JTrigger/JBuildL2.hh"
17 
18 #include "JFit/JFitToolkit.hh"
19 #include "JFit/JLine1Z.hh"
20 #include "JFit/JLine3Z.hh"
21 #include "JFit/JModel.hh"
22 #include "JFit/JSimplex.hh"
23 #include "JFit/JLine3ZRegressor.hh"
24 #include "JFit/JMEstimator.hh"
25 
26 #include "JReconstruction/JEvt.hh"
29 
30 #include "JLang/JPredicate.hh"
31 
33 
35 
36 #include "Jeep/JMessage.hh"
37 
38 
39 /**
40  * \author mdejong, gmaggi
41  */
42 
43 namespace JRECONSTRUCTION {}
44 namespace JPP { using namespace JRECONSTRUCTION; }
45 
46 namespace JRECONSTRUCTION {
47 
49  using JFIT::JRegressor;
50  using JFIT::JLine3Z;
51  using JFIT::JSimplex;
52 
53 
54  /**
55  * Wrapper class to make intermediate fit of muon trajectory.
56  *
57  * The JMuonSimplex fit uses one or more start values (usually taken from the output of JMuonPrefit) and
58  * produces new start values for subsequent fits (usually JMuonGandalf).\n
59  * All hits of which the PMT position lies within a set road width (JMuonSimplexParameters_t::roadWidth_m) and
60  * time is within a set window (JMuonSimplexParameters_t::TMin_ns, JMuonSimplexParameters_t::TMax_ns) around the Cherenkov hypothesis are taken.\n
61  * In case there are multiple hits from the same PMT is the specified window,
62  * the first hit is taken and the other hits are discarded.\n
63  * The chi-squared is based on an M-estimator of the time residuals.\n
64  */
65  struct JMuonSimplex :
67  public JRegressor<JLine3Z, JSimplex>
68  {
72 
73  using JRegressor_t::operator();
74 
75  /**
76  * Constructor
77  *
78  * \param parameters parameters
79  * \param router module router
80  * \param debug debug
81  */
83  const JModuleRouter& router,
84  const int debug = 0) :
85  JMuonSimplexParameters_t(parameters),
86  JRegressor_t(parameters.sigma_ns),
87  router(router)
88  {
89  using namespace JFIT;
90 
91  this->estimator.reset(new JMEstimatorLorentzian());
92 
95  }
96 
97 
98  /**
99  * Fit function.
100  *
101  * \param event event
102  * \param in start values
103  * \return fit results
104  */
106  {
107  using namespace std;
108  using namespace JFIT;
109  using namespace JTRIGGER;
110 
111  const JBuildL0<hit_type> buildL0;
113 
114  buffer_type dataL0;
115  buffer_type dataL1;
116 
117  const KM3NETDAQ::JDAQTimeslice timeslice(event, true);
118  JSuperFrame2D<JHit> buffer;
119 
120  for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
121 
122  if (router.hasModule(i->getModuleID())) {
123 
124  buffer(*i, router.getModule(i->getModuleID()));
125 
126  if (useL0) {
127  buildL0(buffer, back_inserter(dataL0));
128  }
129  {
130  buildL2(buffer, back_inserter(dataL1));
131  }
132  }
133  }
134 
135  return (*this)(dataL0, dataL1, in);
136  }
137 
138 
139  /**
140  * Fit function.
141  *
142  * \param dataL0 L0 hit data
143  * \param dataL1 L1 hit data
144  * \param in start values
145  * \return fit results
146  */
147  JEvt operator()(const buffer_type& dataL0,
148  const buffer_type& dataL1,
149  const JEvt& in)
150  {
151  using namespace std;
152  using namespace JFIT;
153  using namespace JGEOMETRY3D;
154  using namespace JLANG;
155 
156  JEvt out;
157 
159 
160  data.reserve(dataL0.size() +
161  dataL1.size());
162 
163 
164  for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
165 
166  const JRotation3D R (getDirection(*track));
167  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
169 
170  data.clear();
171 
172  for (buffer_type::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {
173 
174  hit_type hit(*i);
175 
176  hit.rotate(R);
177 
178  if (match(hit)) {
179  data.push_back(hit);
180  }
181  }
182 
183  if (useL0) {
184 
185  buffer_type::iterator __end = data.end();
186 
187  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
188 
189  if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
190 
191  hit_type hit(*i);
192 
193  hit.rotate(R);
194 
195  if (match(hit)) {
196  data.push_back(hit);
197  }
198  }
199  }
200  }
201 
202 
203  this->step.resize(5);
204 
205  this->step[0] = JLine3Z(JLine1Z(JVector3D(0.5, 0.0, 0.0), 0.0));
206  this->step[1] = JLine3Z(JLine1Z(JVector3D(0.0, 0.5, 0.0), 0.0));
207  this->step[2] = JLine3Z(JLine1Z(JVector3D(0.0, 0.0, 0.0), 1.0));
208  this->step[3] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.005, 0.0));
209  this->step[4] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.0, 0.005));
210 
211  const int NDF = getCount(data.begin(), data.end()) - this->step.size();
212 
213  if (NDF > 0) {
214 
215  const double chi2 = (*this)(JLine3Z(tz), data.begin(), data.end());
216 
217  JTrack3D tb(this->value);
218 
219  tb.rotate_back(R);
220 
221  out.push_back(getFit(JHistory(track->getHistory()).add(JMUONSIMPLEX), tb, getQuality(chi2, NDF), NDF));
222 
223  out.rbegin()->setW(track->getW());
224  }
225  }
226 
227  return out;
228  }
229 
230 
232  };
233 }
234 
235 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
static int debug
debug level (default is off).
Definition: JMessage.hh:45
Data regression method for JFIT::JLine3Z.
Template definition of a data regressor of given model.
Definition: JRegressor.hh:68
JEvt operator()(const buffer_type &dataL0, const buffer_type &dataL1, const JEvt &in)
Fit function.
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
int getModuleID() const
Get module identifier.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
const JModuleRouter & router
double getQuality(const double chi2, const int NDF)
Get quality of fit.
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
double TMaxLocal_ns
time window for local coincidences [ns]
*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
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:36
Template L2 builder.
Definition: JBuildL2.hh:45
JDirection3D getDirection(const Vec &dir)
Get direction.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
Wrapper class to make intermediate fit of muon trajectory.
Definition: JMuonSimplex.hh:65
Regressor function object for JLine3Z fit using JSimplex minimiser.
JRegressor< JLine3Z, JSimplex > JRegressor_t
Definition: JMuonSimplex.hh:69
JPosition3D getPosition(const Vec &pos)
Get position.
Data time slice.
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JEvt &in)
Fit function.
General purpose messaging.
Direct access to module in detector data structure.
std::vector< hit_type > buffer_type
Definition: JMuonSimplex.hh:71
Data structure for L2 parameters.
Reduced data structure for L1 hit.
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
JFIT::JHistory JHistory
Definition: JHistory.hh:354
Data structure for set of track fit results.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
static const int JMUONSIMPLEX
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
bool hasModule(const JObjectID &id) const
Has module.
int getCount(const T &hit)
Get hit count.
JMuonSimplex(const JMuonSimplexParameters_t &parameters, const JModuleRouter &router, const int debug=0)
Constructor.
Definition: JMuonSimplex.hh:82
2-dimensional frame with time calibrated data from one optical module.
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JSimplex.hh:237
double ctMin
minimal cosine space angle between PMT axes
Template L0 hit builder.
Definition: JBuildL0.hh:35
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] 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:46
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
Maximum likelihood estimator (M-estimators).
Lorentzian M-estimator.
Definition: JMEstimator.hh:65
int debug
debug level
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).