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 
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
double getQuality(const double chi2, const int NDF)
Get quality of fit.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
static struct JRECONSTRUCTION::JMuonGandalf::@56 compare
Auxiliary data structure for sorting of hits.
Energy loss of muon.
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
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:36
JRegressor< JLine3Z, JGandalf > JRegressor_t
Definition: JMuonGandalf.hh:76
static parameter_type pT()
Definition: JLine1Z.hh:182
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:225
JDirection3D getDirection(const Vec &dir)
Get direction.
static const int JMUONGANDALF
do set_variable OUTPUT_DIRECTORY $WORKDIR T
double E_GeV
Energy of muon at vertex [GeV].
JPosition3D getPosition(const Vec &pos)
Get position.
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v2.0.0 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.
static parameter_type pX()
Definition: JLine1Z.hh:180
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
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.
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
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:248
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
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 source JAcoustics sh $DETECTOR_ID typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:36
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
JTOOLS::JRange< double > JZRange
Definition: JFit/JModel.hh:21