Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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),
102  summary(summary)
103  {
104  using namespace JPP;
105 
107  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
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(4);
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  }
129 
130  /**
131  * Declaration of the member function that actually performs the reconstruction
132  *
133  * \param event which is a JDAQEvent
134  * \param in input fits
135  */
137  {
138  using namespace std;
139  using namespace JPP;
140 
141  typedef vector<JHitL0> JDataL0_t;
142  const JBuildL0<JHitL0> buildL0;
143 
144  JEvt out;
145 
146  JDataL0_t dataL0;
147  buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
148 
149  for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
150 
151  JPoint4D vx(getPosition(*shower), shower->getT());
152 
154 
156 
157  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
158 
159  const double rate_Hz = summary.getRate(*i);
160 
161  if (match(*i)) {
162  data.push_back(JHitW0(*i, rate_Hz));
163  }
164  }
165 
166  // select first hit
167 
168  sort(data.begin(), data.end(), JHitL0::compare);
169 
170  vector<JHitW0>::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
171 
172  const int NDF = distance(data.begin(), __end) - this->parameters.size();
173 
174  if (NDF > 0) {
175 
176  // set fit parameters
177  for (vector<double>::iterator i = Ev.begin(); i != Ev.end(); ++i) {
178  JPoint4E sh(vx, *i);
179  double chi2 = (*this)(sh, data.begin(), __end);
180 
181  JShower3E sh_fit(this->value.getPosition(), JDirection3D(),
182  this->value.getT(), this->value.getE());
183 
184  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), sh_fit,
185  getQuality(chi2), NDF, sh_fit.getE()));
186 
187  }
188  }
189  }
190 
191  return out;
192  }
193 
194 
197  };
198 }
199 
200 #endif
201 
Auxiliary methods to evaluate Poisson probabilities and chi2.
static int debug
debug level (default is off).
Definition: JMessage.hh:45
Template definition of a data regressor of given model.
Definition: JRegressor.hh:70
Linear fit methods.
General exception.
Definition: JException.hh:24
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:336
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Algorithms for hit clustering and sorting.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
Data structure for vertex fit.
Definition: JPoint4D.hh:22
JRegressor< JPoint4E, JGandalf > JRegressor_t
Router for direct addressing of module data in detector data structure.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
double getRate() const
Get default rate.
*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
class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA ...
static double Vmax_npe
Maximal integral of PDF [npe].
Basic data structure for time and time over threshold information of hit.
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JFIT::JEvt &in)
Declaration of the member function that actually performs the reconstruction.
3D track with energy.
Definition: JTrack3E.hh:30
Data structure for detector geometry and calibration.
double TMin_ns
minimum time for local coincidences [ns]
static struct JTRIGGER::JHitL0::compare compare
static const int JSHOWERPOSITIONFIT
double VMax_npe
maximum number of of photo-electrons
Basic data structure for L0 hit.
Auxiliary class to extract a subset of optical modules from a detector.
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
double TMax_ns
maximum time for local coincidences [ns]
int En
number of points to scan in energy range
Data structure for vertex fit.
Definition: JPoint4E.hh:22
JPosition3D getPosition(const Vec &pos)
Get position.
double DMax_m
maximal distance to optical module [m]
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
Data time slice.
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
Direct access to module in detector data structure.
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:84
Reduced data structure for L1 hit.
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
Data regression method for JFIT::JPoint4E from a bright point isoptropic emission PDF...
static parameter_type pT()
Definition: JPoint4E.hh:72
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.
Auxiliary class to define a range between two values.
then fatal The output file must have the wildcard in the e g root fi 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:48
static parameter_type pX()
Definition: JPoint4E.hh:69
static parameter_type pY()
Definition: JPoint4E.hh:70
JShowerPositionFit(const JShowerPositionFitParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string pdfFile, const int debug=0)
Parameterized constructor.
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
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:335
Template specialisation of class JModel to match hit with bright point.
Definition: JFit/JModel.hh:121
Match operator for Cherenkov light from shower in any direction.
Basic data structure for L1 hit.
int debug
debug level
static parameter_type pZ()
Definition: JPoint4E.hh:71
Regressor function object for JPoint4E fit using JGandalf minimiser.