Jpp  19.1.0
the software that should make you happy
JShowerPositionFit.hh
Go to the documentation of this file.
1 #ifndef JSHOWERPOSITIONFIT_INCLUDE
2 #define JSHOWERPOSITIONFIT_INCLUDE
3 
4 #include <string>
5 #include <iostream>
6 #include <set>
7 #include <vector>
8 #include <algorithm>
9 #include <memory>
10 #include <math.h>
11 
14 
15 #include "JTrigger/JHit.hh"
16 #include "JTrigger/JTimeslice.hh"
17 #include "JTrigger/JHitL0.hh"
18 #include "JTrigger/JHitL1.hh"
19 #include "JTrigger/JHitR1.hh"
20 #include "JTrigger/JBuildL0.hh"
21 #include "JTrigger/JBuildL2.hh"
22 #include "JTrigger/JAlgorithm.hh"
23 #include "JTrigger/JMatch3G.hh"
24 #include "JTrigger/JBind2nd.hh"
25 
27 
28 #include "JFit/JFitToolkit.hh"
29 #include "JFit/JEstimator.hh"
30 #include "JFit/JPoint4E.hh"
31 #include "JFit/JModel.hh"
32 #include "JFit/JSimplex.hh"
34 
36 #include "JReconstruction/JEvt.hh"
39 
40 #include "JGeometry3D/JVector3D.hh"
41 #include "JGeometry3D/JShower3D.hh"
42 #include "JGeometry3D/JShower3E.hh"
43 #include "JGeometry3D/JTrack3D.hh"
44 #include "JGeometry3D/JTrack3E.hh"
46 
48 #include "JDetector/JDetector.hh"
50 
51 #include "JTools/JPermutation.hh"
52 #include "JTools/JRange.hh"
53 #include "JLang/JSharedPointer.hh"
54 
55 /**
56  * \author adomi, vcarretero
57  */
58 namespace JRECONSTRUCTION {}
59 namespace JPP { using namespace JRECONSTRUCTION; }
60 
61 namespace JRECONSTRUCTION {
62 
65  using JFIT::JPoint4D;
66  using JFIT::JPoint4E;
67  using JFIT::JGandalf;
68  using JFIT::JRegressor;
69 
70  /**
71  * class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA
72  */
74 
76  public JRegressor<JPoint4E, JGandalf>
77 
78  {
79 
81  using JRegressor_t::operator();
83  public:
84 
85  /**
86  * Parameterized constructor
87  *
88  * \param parameters struct that holds default-optimized parameters for the reconstruction
89  * \param router module router, this is built via detector file.
90  * \param summary summary router
91  * \param pdfFile file name with PDF
92  * \param debug debug
93  */
95  const JModuleRouter& router,
96  const JSummaryRouter& summary,
97  const std::string pdfFile,
98  const int debug = 0):
100  JRegressor_t(pdfFile, parameters.TTS_ns),
101  router(router),
103  {
104  using namespace JPP;
105 
107  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
108  JRegressor_t::Vmax_npe = VMax_npe;
109  JRegressor_t::MAXIMUM_ITERATIONS = NMax;
110  JRegressor_t::EPSILON = 1e-3;
111 
112  if (Emin_GeV > Emax_GeV || En <= 1) {
113  THROW(JException, "Invalid energy input " << Emin_GeV << ' ' << Emax_GeV << ' ' << En);
114  }
115 
116  const double base = std::pow((Emax_GeV / Emin_GeV), 1.0 / (En - 1));
117 
118  for (int i = 0; i != En; ++i) {
119  Ev.push_back(Emin_GeV * std::pow(base, i));
120  }
121 
122  this->parameters.resize(5);
123 
124  this->parameters[0] = JPoint4E::pX();
125  this->parameters[1] = JPoint4E::pY();
126  this->parameters[2] = JPoint4E::pZ();
127  this->parameters[3] = JPoint4E::pT();
128  this->parameters[4] = JPoint4E::pE();
129  }
130 
131  /**
132  * Declaration of the member function that actually performs the reconstruction
133  *
134  * \param event which is a JDAQEvent
135  * \param in input fits
136  */
138  {
139  using namespace std;
140  using namespace JPP;
141 
142  typedef vector<JHitL0> JDataL0_t;
143  const JBuildL0<JHitL0> buildL0;
144 
145  JEvt out;
146 
147  JDataL0_t dataL0;
148  buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
149 
150  for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
151 
152  JPoint4D vx(getPosition(*shower), shower->getT());
153 
155 
156  const JFIT::JModel<JPoint4D> match(vx, DMax_m, JRegressor_t::T_ns);
157 
158  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
159 
160  const double rate_Hz = summary.getRate(*i);
161 
162  if (match(*i)) {
163  data.push_back(JHitW0(*i, rate_Hz));
164  }
165  }
166 
167  // select first hit
168 
169  sort(data.begin(), data.end(), JHitL0::compare);
170 
171  vector<JHitW0>::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
172 
173  const int NDF = distance(data.begin(), __end) - this->parameters.size();
174 
175  if (NDF > 0) {
176 
177  // set fit parameters
178  for (vector<double>::iterator e = Ev.begin(); e != Ev.end(); ++e) {
179  JPoint4E sh(vx, *e);
180  double chi2 = (*this)(sh, data.begin(), __end);
181 
182  JShower3E sh_fit(this->value.getPosition(), JDirection3D(),
183  this->value.getT(), this->value.getE());
184 
185  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), sh_fit,
186  getQuality(chi2), NDF, sh_fit.getE()));
187 
188  }
189  }
190  }
191 
192  return out;
193  }
194 
195 
198  };
199 }
200 
201 #endif
202 
Algorithms for hit clustering and sorting.
Auxiliary class to extract a subset of optical modules from a detector.
Data structure for detector geometry and calibration.
Linear fit methods.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Auxiliary methods to evaluate Poisson probabilities and chi2.
Basic data structure for L0 hit.
Basic data structure for L1 hit.
Reduced data structure for L1 hit.
Match operator for Cherenkov light from shower in any direction.
int debug
debug level
Definition: JSirene.cc:69
Direct access to module in detector data structure.
Auxiliary class to define a range between two values.
Data regression method for JFIT::JPoint4E from a bright point isoptropic emission PDF.
Basic data structure for time and time over threshold information of hit.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Router for direct addressing of module data in detector data structure.
Data structure for set of track fit results.
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:86
Data structure for vertex fit.
Definition: JPoint4D.hh:24
Data structure for vertex fit.
Definition: JPoint4E.hh:24
static parameter_type pZ()
Definition: JPoint4E.hh:71
static parameter_type pX()
Definition: JPoint4E.hh:69
static parameter_type pY()
Definition: JPoint4E.hh:70
static parameter_type pE()
Definition: JPoint4E.hh:73
static parameter_type pT()
Definition: JPoint4E.hh:72
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:35
3D track with energy.
Definition: JTrack3E.hh:32
General exception.
Definition: JException.hh:24
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:23
class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA
JShowerPositionFit(const JShowerPositionFitParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string pdfFile, const int debug=0)
Parameterized constructor.
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JFIT::JEvt &in)
Declaration of the member function that actually performs the reconstruction.
JRegressor< JPoint4E, JGandalf > JRegressor_t
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
double getRate() const
Get default rate.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:105
static const int JSHOWERPOSITIONFIT
JPosition3D getPosition(const Vec &pos)
Get position.
double getQuality(const double chi2, const int NDF)
Get quality of fit.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
JFIT::JHistory JHistory
Definition: JHistory.hh:354
Definition: JSTDTypes.hh:14
JHistory & add(const int type)
Add event to history.
Definition: JHistory.hh:295
Template specialisation of class JModel to match hit with bright point.
Definition: JFit/JModel.hh:123
Regressor function object for JPoint4E fit using JGandalf minimiser.
Template definition of a data regressor of given model.
Definition: JRegressor.hh:70
double VMax_npe
maximum number of of photo-electrons
double TMin_ns
minimum time for local coincidences [ns]
double TMax_ns
maximum time for local coincidences [ns]
int En
number of points to scan in energy range
double DMax_m
maximal distance to optical module [m]
Auxiliary data structure for sorting of hits.
Definition: JHitL0.hh:85