Jpp  15.0.0-rc.2
the software that should make you happy
 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 JFIT;
102 
103  if (this->getRmax() < roadWidth_m) {
104  roadWidth_m = this->getRmax();
105  }
106 
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 JFIT;
133  using namespace JTRIGGER;
134 
135  const JBuildL0<hit_type> buildL0;
136 
137  buffer_type dataL0;
138 
139  buildL0(event, router, true, back_inserter(dataL0));
140 
141  return (*this)(dataL0, in);
142  }
143 
144 
145  /**
146  * Fit function.
147  *
148  * \param data hit data
149  * \param in start values
150  * \return fit results
151  */
152  JEvt operator()(const buffer_type& data, const JEvt& in)
153  {
154  using namespace std;
155  using namespace JFIT;
156  using namespace JGEOMETRY3D;
157 
158  JEvt out;
159 
160  for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
161 
162  const JRotation3D R (getDirection(*track));
163  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
164  JRange<double> Z_m;
165 
166  if (track->hasW(JSTART_LENGTH_METRES) &&
167  track->getW(JSTART_LENGTH_METRES) > 0.0) {
168  Z_m = JZRange(ZMin_m, ZMax_m + track->getW(JSTART_LENGTH_METRES));
169  }
170 
171  const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns, Z_m);
172 
173  // hit selection based on start value
174 
175  vector<JHitW0> buffer;
176 
177  for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
178 
179  JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));
180 
181  hit.rotate(R);
182 
183  if (match(hit)) {
184  buffer.push_back(hit);
185  }
186  }
187 
188  // select first hit
189 
190  sort(buffer.begin(), buffer.end(), JHitL0::compare);
191 
192  vector<JHitW0>::iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
193 
194 
195  const int NDF = distance(buffer.begin(), __end) - this->parameters.size();
196 
197  if (NDF > 0) {
198 
199  // set fit parameters
200 
201  if (track->getE() > 0.1)
202  JRegressor_t::E_GeV = track->getE();
203  else
205 
206  const double chi2 = (*this)(JLine3Z(tz), buffer.begin(), __end);
207 
208  JTrack3D tb(this->value);
209 
210  tb.rotate_back(R);
211 
212  out.push_back(getFit(JHistory(track->getHistory()).add(JMUONGANDALF), tb, getQuality(chi2), NDF));
213 
214  // set additional values
215 
216  out.rbegin()->setV(this->V.size(), this->V);
217 
218  out.rbegin()->setW(track->getW());
219  out.rbegin()->setW(JGANDALF_BETA0_RAD, sqrt(this->error.getDX() * this->error.getDX() +
220  this->error.getDY() * this->error.getDY()));
221  out.rbegin()->setW(JGANDALF_BETA1_RAD, sqrt(this->error.getDX() * this->error.getDY()));
222  out.rbegin()->setW(JGANDALF_CHI2, chi2);
223  out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS, distance(buffer.begin(), __end));
224  out.rbegin()->setW(JGANDALF_LAMBDA, this->lambda);
225  out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, this->numberOfIterations);
226  }
227  }
228 
229  return out;
230  }
231 
232 
235  int debug;
236  };
237 }
238 
239 #endif
240 
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
double TTS_ns
transition-time spread [ns]
Auxiliary methods for PDF calculations.
const JModuleRouter & router
double getQuality(const double chi2, const int NDF)
Get quality of fit.
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
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
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
double VMax_npe
maximum number of of photo-electrons
Wrapper class to make final fit of muon trajectory.
Definition: JMuonGandalf.hh:72
std::vector< hit_type > buffer_type
Definition: JMuonGandalf.hh:78
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.
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
static const int JMUONGANDALF
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.1.0-7-g97845ea 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:43
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:280
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 CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
JTOOLS::JRange< double > JZRange
Definition: JFit/JModel.hh:21
static struct JTRIGGER::JHitL0::@77 compare
Auxiliary data structure for sorting of hits.