Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMuonGandalf.hh
Go to the documentation of this file.
1 #ifndef __JRECONSTRUCTION__JMUONGANDALF__
2 #define __JRECONSTRUCTION__JMUONGANDALF__
3 
4 #include <string>
5 #include <iostream>
6 #include <iomanip>
7 #include <vector>
8 #include <algorithm>
9 
13 
14 #include "JTrigger/JHitL0.hh"
15 #include "JTrigger/JBuildL0.hh"
16 
17 #include "JFit/JFitToolkit.hh"
18 #include "JFit/JLine1Z.hh"
19 #include "JFit/JLine3Z.hh"
20 #include "JFit/JModel.hh"
21 #include "JFit/JGandalf.hh"
22 #include "JFit/JLine3ZRegressor.hh"
23 #include "JFit/JHitW0.hh"
24 
25 #include "JReconstruction/JEvt.hh"
28 
29 #include "JPhysics/JPDFTable.hh"
30 #include "JPhysics/JPDFToolkit.hh"
31 #include "JPhysics/JGeane.hh"
32 
34 
36 
38 
39 #include "JTools/JRange.hh"
40 
41 #include "Jeep/JMessage.hh"
42 
43 
44 /**
45  * \author mdejong, gmaggi
46  */
47 
48 namespace JRECONSTRUCTION {}
49 namespace JPP { using namespace JRECONSTRUCTION; }
50 
51 namespace JRECONSTRUCTION {
52 
55  using JFIT::JRegressor;
56  using JFIT::JLine3Z;
57  using JFIT::JGandalf;
58 
59 
60  /**
61  * Wrapper class to make final fit of muon trajectory.
62  *
63  * The JMuonGandalf fit uses one or more start values (usually taken from the output of JMuonSimplex).\n
64  * All hits of which the PMT position lies within a set road width (JMuonGandalfParameters_t::roadWidth_m) and
65  * time is within a set window (JMuonGandalfParameters_t::TMin_ns, JMuonGandalfParameters_t::TMax_ns) around the Cherenkov hypothesis are taken.\n
66  * In case there are multiple hits from the same PMT is the specified window,
67  * the first hit is taken and the other hits are discarded.\n
68  * The PDF is accordingly evaluated, i.e.\ the normalised probability for a first hit at the given time of the hit is taken.
69  * The normalisation is consistently based on the specified time window.\n
70  * Note that this hit selection is unbiased with respect to the PDF of a single PMT.
71  */
72  struct JMuonGandalf :
74  public JRegressor<JLine3Z, JGandalf>
75  {
79 
80  using JRegressor_t::operator();
81 
82  /**
83  * Constructor
84  *
85  * \param parameters parameters
86  * \param router module router
87  * \param summary summary file router
88  * \param pdf_file PDF file
89  * \param debug debug
90  */
92  const JModuleRouter& router,
93  const JSummaryRouter& summary,
94  const std::string& pdf_file,
95  const int debug = 0) :
96  JMuonGandalfParameters_t(parameters),
97  JRegressor_t(pdf_file, parameters.TTS_ns),
98  router (router),
99  summary(summary)
100  {
101  using namespace JPP;
102 
103  if (this->getRmax() < roadWidth_m) {
104  roadWidth_m = this->getRmax();
105  }
106 
109  JRegressor_t::Vmax_npe = 10.0;
111 
112  this->parameters.resize(5);
113 
114  this->parameters[0] = JLine3Z::pX();
115  this->parameters[1] = JLine3Z::pY();
116  this->parameters[2] = JLine3Z::pT();
117  this->parameters[3] = JLine3Z::pDX();
118  this->parameters[4] = JLine3Z::pDY();
119  }
120 
121 
122  /**
123  * Fit function.
124  *
125  * \param event event
126  * \param in start values
127  * \return fit results
128  */
130  {
131  using namespace std;
132  using namespace JPP;
133 
134  const JBuildL0<hit_type> buildL0;
135 
136  buffer_type dataL0;
137 
138  buildL0(event, router, true, back_inserter(dataL0));
139 
140  return (*this)(dataL0, in);
141  }
142 
143 
144  /**
145  * Fit function.
146  *
147  * \param data hit data
148  * \param in start values
149  * \return fit results
150  */
151  JEvt operator()(const buffer_type& data, const JEvt& in)
152  {
153  using namespace std;
154  using namespace JPP;
155 
156  JEvt out;
157 
158  for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
159 
160  const JRotation3D R (getDirection(*track));
161  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
162  JRange<double> Z_m;
163 
164  if (track->hasW(JSTART_LENGTH_METRES)) {
165  Z_m = JZRange(ZMin_m, ZMax_m + track->getW(JSTART_LENGTH_METRES));
166  }
167 
168  const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns, Z_m);
169 
170  // hit selection based on start value
171 
172  vector<JHitW0> buffer;
173 
174  for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
175 
176  const double rate_Hz = summary.getRate(*i);
177 
178  JHitW0 hit(*i, rate_Hz);
179 
180  hit.rotate(R);
181 
182  if (match(hit)) {
183  buffer.push_back(hit);
184  }
185  }
186 
187  // select first hit
188 
189  sort(buffer.begin(), buffer.end(), compare);
190 
191  vector<JHitW0>::iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
192 
193 
194  const int NDF = distance(buffer.begin(), __end) - this->parameters.size();
195 
196  if (NDF > 0) {
197 
198  // set fit parameters
199 
200  if (track->getE() > 0.1)
201  JRegressor_t::E_GeV = track->getE();
202  else
204 
205  const double chi2 = (*this)(JLine3Z(tz), buffer.begin(), __end);
206 
207  JTrack3D tb(this->value);
208 
209  tb.rotate_back(R);
210 
211  out.push_back(getFit(JHistory(track->getHistory()).add(JMUONGANDALF), tb, getQuality(chi2), NDF));
212 
213  // set additional values
214 
215  out.rbegin()->setW(track->getW());
216  out.rbegin()->setW(JGANDALF_BETA0_RAD, sqrt(this->error.getDX() * this->error.getDX() +
217  this->error.getDY() * this->error.getDY()));
218  out.rbegin()->setW(JGANDALF_BETA1_RAD, sqrt(this->error.getDX() * this->error.getDY()));
219  out.rbegin()->setW(JGANDALF_CHI2, chi2);
220  out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS, distance(buffer.begin(), __end));
221  out.rbegin()->setW(JGANDALF_LAMBDA, this->lambda);
222  out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, this->numberOfIterations);
223  }
224  }
225 
226  return out;
227  }
228 
229 
230  /**
231  * Auxiliary data structure for sorting of hits.
232  */
233  static const struct {
234  /**
235  * Compare hits by PMT identifier and time.
236  *
237  * \param first first hit
238  * \param second second hit
239  * \return true if first before second; else false
240  */
241  template<class T>
242  bool operator()(const T& first, const T& second) const
243  {
244  using namespace std;
245  using namespace JPP;
246  using namespace KM3NETDAQ;
247 
248  if (equal_to<JDAQPMTIdentifier>()(first, second))
249  return less<JHit>()(first, second);
250  else
251  return less<JDAQPMTIdentifier>()(first, second);
252  }
253  } compare;
254 
257  int debug;
258  };
259 }
260 
261 #endif
262 
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
Auxiliary methods for PDF calculations.
const JModuleRouter & router
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Energy loss of muon.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
double getRate() const
Get default rate.
static const int JGANDALF_BETA1_RAD
angular resolution [rad] from JGandalf.cc
JMuonGandalf(const JMuonGandalfParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string &pdf_file, const int debug=0)
Constructor.
Definition: JMuonGandalf.hh:91
*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:35
JRegressor< JLine3Z, JGandalf > JRegressor_t
Definition: JMuonGandalf.hh:76
static parameter_type pT()
Definition: JLine1Z.hh:182
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:206
Wrapper class to make final fit of muon trajectory.
Definition: JMuonGandalf.hh:72
std::vector< hit_type > buffer_type
Definition: JMuonGandalf.hh:78
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
const JSummaryRouter & summary
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Basic data structure for L0 hit.
static const int JGANDALF_NUMBER_OF_HITS
number of hits from JGandalf.cc
static parameter_type pDX()
Definition: JLine3Z.hh:319
static const int JGANDALF_LAMBDA
control parameter from JGandalf.cc
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:236
static const int JMUONGANDALF
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double E_GeV
Energy of muon at vertex [GeV].
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
static struct JRECONSTRUCTION::JMuonGandalf::@51 compare
Auxiliary data structure for sorting of hits.
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v1.2.3 https://git.km3net.de/common/km3net-dataformat.
Router for fast addressing of summary data in JDAQSummaryslice data structure as a function of the op...
static parameter_type pY()
Definition: JLine1Z.hh:181
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
General purpose messaging.
static parameter_type pDY()
Definition: JLine3Z.hh:320
Direct access to module in detector data structure.
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:52
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:126
static parameter_type pX()
Definition: JLine1Z.hh:180
JFIT::JHistory JHistory
Definition: JHistory.hh:283
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.
Definition: JEvt.hh:294
Regressor function object for JLine3Z fit using JGandalf minimiser.
Auxiliary class to define a range between two values.
static double Vmax_npe
Maximal integral of PDF [npe].
Data structure for L0 hit.
Definition: JHitL0.hh:27
JEvt operator()(const buffer_type &data, const JEvt &in)
Fit function.
static const int JGANDALF_NUMBER_OF_ITERATIONS
number of iterations from JGandalf.cc
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
JDirection3D getDirection(const Vec &v)
Get direction.
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:247
static const int JGANDALF_CHI2
chi2 from JGandalf.cc
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JEvt &in)
Fit function.
Template L0 hit builder.
Definition: JBuildL0.hh:35
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:185
JTOOLS::JRange< double > JZRange
Definition: JFit/JModel.hh:21
JPosition3D getPosition(const Vec &v)
Get position.