Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JShowerFit.hh
Go to the documentation of this file.
1 #ifndef JSHOWERFIT_INCLUDE
2 #define JSHOWERFIT_INCLUDE
3 
4 #include <string>
5 #include <iostream>
6 #include <set>
7 #include <vector>
8 #include <algorithm>
9 #include <memory>
10 
13 
14 #include "JTrigger/JHit.hh"
15 #include "JTrigger/JTimeslice.hh"
16 #include "JTrigger/JHitL0.hh"
17 #include "JTrigger/JHitL1.hh"
18 #include "JTrigger/JHitR1.hh"
19 #include "JTrigger/JBuildL0.hh"
20 #include "JTrigger/JBuildL2.hh"
21 #include "JTrigger/JAlgorithm.hh"
22 #include "JTrigger/JMatch3G.hh"
23 #include "JTrigger/JBind2nd.hh"
24 
26 
28 #include "JFit/JFitToolkit.hh"
29 #include "JFit/JPoint4D.hh"
30 #include "JFit/JModel.hh"
31 #include "JFit/JGandalf.hh"
32 
33 #include "JAAnet/JAAnetToolkit.hh"
34 
36 #include "JReconstruction/JEvt.hh"
40 
41 #include "JDetector/JDetector.hh"
43 
44 #include "JGeometry3D/JVector3D.hh"
46 #include "JGeometry3D/JOmega3D.hh"
47 #include "JGeometry3D/JTrack3D.hh"
48 #include "JGeometry3D/JTrack3E.hh"
49 #include "JGeometry3D/JShower3E.hh"
50 
51 #include "JLang/JSharedPointer.hh"
52 
53 
54 /**
55  * \author adomi
56  */
57 namespace JRECONSTRUCTION {}
58 namespace JPP { using namespace JRECONSTRUCTION; }
59 
60 namespace JRECONSTRUCTION {
61 
65  using JFIT::JRegressor;
66  using JFIT::JEnergy;
67  using JFIT::JShower3EZ;
68  using JFIT::JGandalf;
69  /**
70  * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
71  */
72  class JShowerFit :
74  public JRegressor<JShower3EZ, JGandalf>
75  {
76 
78  using JRegressor_t::operator();
79 
80  public:
81 
82  /**
83  * Parameterized constructor
84  *
85  * \param parameters struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA.
86  * \param router module router, this is built via detector file.
87  * \param summary summary router
88  * \param pdfFile PDF file
89  * \param correct energy correction
90  * \param debug debug
91  */
93  const JModuleRouter& router,
94  const JSummaryRouter& summary,
95  const std::string pdfFile,
97  const int debug = 0):
98  JShowerFitParameters_t(parameters),
99  JRegressor_t(pdfFile),
100  router(router),
101  summary(summary),
102  correct(correct)
103  {
104  using namespace JPP;
105 
107  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
108  JRegressor_t::Vmax_npe = parameters.Vmax_npe;
110  JRegressor_t::EPSILON = parameters.epsilon;
111  this->parameters.resize(3);
112 
113  this->parameters[0] = JShower3EZ::pDX();
114  this->parameters[1] = JShower3EZ::pDY();
115  this->parameters[2] = JShower3EZ::pE();
116 
117  this->estimator.reset(getMEstimator(parameters.mestimator));
118  }
119 
120  /**
121  * Declaration of the member function that actually performs the reconstruction
122  *
123  * \param event = JDAQEvent
124  * \param in = input fits
125  */
127  {
128  using namespace std;
129  using namespace JPP;
130 
131  typedef vector<JHitL0> JDataL0_t;
132  JBuildL0<JHitL0> buildL0;
133 
134  JEvt out;
135 
137 
138  JDataL0_t dataL0;
139 
140  buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
141 
142  for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
143 
144  const JVertex3D Vertex(getPosition(*shower), shower->getT());
145 
147 
148  const double distance = DMax_m + DStep_m * log10(shower->getE());
149 
150  const JFIT::JModel<JPoint4D> match(Vertex, distance, JRegressor_t::T_ns);
151 
152  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
153  if (match(*i)) {
154  top.insert(i->getPMTIdentifier());
155  }
156  }
157 
158  const JDetectorSubset_t subdetector(detector, Vertex.getPosition(), distance);
159 
160  const JRotation3D R(getDirection(*shower));
161 
162  vector<JPMTW0> buffer;
163 
164  for (JDetectorSubset_t::const_iterator module = subdetector.begin();
165  module != subdetector.end(); ++module) {
166 
167  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
168 
169  JModule dom(*module);
170 
171  dom.rotate(R);
172 
173  for (size_t i = 0; i != dom.size(); ++i) {
174 
175  if (getDAQStatus(frame, *module, i) &&
176  getPMTStatus(frame, *module, i) &&
177  frame[i].is_valid() &&
178  !module->getPMT(i).has(PMT_DISABLE)) {
179  const JDAQPMTIdentifier id(module->getID(), i);
180 
181  const double rate_Hz = summary.getRate(id);
182  const size_t count = top.count(id);
183 
184  buffer.push_back(JPMTW0(dom.getPMT(i), rate_Hz, count));
185  }
186  }
187  }
188 
189  double chi2 = (*this)(JShower3EZ(JVertex3D(JVector3D(0,0,0), shower->getT()), JVersor3Z(),
190  shower->getE()), buffer.begin(), buffer.end());
191 
192  double NDF = getCount(buffer.begin(), buffer.end()) - this->parameters.size();
193 
194  JShower3E sh_fit(this->value.getPosition(), this->value.getDirection(),
195  this->value.getT(), correct(this->value.getE()));
196 
197  sh_fit.rotate_back(R);
198 
199  sh_fit.add(Vertex.getPosition());
200 
201  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), sh_fit, getQuality(chi2),
202  NDF, sh_fit.getE()));
203 
204  out.rbegin()->setW(JSHOWERFIT_ENERGY, this->value.getE()); // Uncorrected Energy
205  }
206 
207  return out;
208  }
209 
213  };
214 }
215 
216 #endif
217 
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
double TMax_ns
maximum time for local coincidences [ns]
static parameter_type pDX()
Definition: JShower3Z.hh:171
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:323
JEvt operator()(const KM3NETDAQ::JDAQEvent &event, const JFIT::JEvt &in)
Declaration of the member function that actually performs the reconstruction.
Definition: JShowerFit.hh:126
JShowerFit(const JShowerFitParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string pdfFile, const JShowerEnergyCorrection &correct, const int debug=0)
Parameterized constructor.
Definition: JShowerFit.hh:92
Data structure for a composite optical module.
Definition: JModule.hh:67
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
double TMin_ns
minimum time for local coincidences [ns]
Auxiliary class for handling PMT geometry, rate and response.
Definition: JPMTW0.hh:22
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class for correction of energy determined by JShowerEnergy.cc.
double DMax_m
maximal distance to optical module [m]
Router for direct addressing of module data in detector data structure.
Rotation matrix.
Definition: JRotation3D.hh:111
double getRate() const
Get default rate.
then usage $script< input file >[option] nPossible options count
Definition: JVolume1D.sh:31
*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
const JModuleRouter & router
Definition: JShowerFit.hh:210
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
Basic data structure for time and time over threshold information of hit.
3D track with energy.
Definition: JTrack3E.hh:30
Data structure for detector geometry and calibration.
JTime & add(const JTime &value)
Addition operator.
const JShowerEnergyCorrection & correct
Definition: JShowerFit.hh:212
Basic data structure for L0 hit.
Auxiliary class to extract a subset of optical modules from a detector.
const JSummaryRouter & summary
Definition: JShowerFit.hh:211
static const int JSHOWERFIT_ENERGY
uncorrected energy [GeV] from JShowerFit.cc
Data structure for fit of straight line in positive z-direction with energy.
Definition: JShower3EZ.hh:27
Detector file.
Definition: JHead.hh:226
JDirection3D getDirection(const Vec &dir)
Get direction.
double getE() const
Get energy.
Definition: JTrack3E.hh:93
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
set_variable E_E log10(E_{fit}/E_{#mu})"
Data storage class for rate measurements of all PMTs in one module.
static const int PMT_DISABLE
KM3NeT Data Definitions v3.3.0-2-g5cc95cf https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Data time slice.
JAxis3D & rotate_back(const JRotation3D &R)
Rotate back axis.
Definition: JAxis3D.hh:240
double epsilon
tolerance for minimisation
static parameter_type pE()
Definition: JShower3EZ.hh:205
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
void rotate(const JRotation3D &R)
Rotate module.
Definition: JModule.hh:314
JRegressor< JShower3EZ, JGandalf > JRegressor_t
Definition: JShowerFit.hh:77
Detector subset without binary search functionality.
static const int JSHOWERCOMPLETEFIT
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:84
Reduced data structure for L1 hit.
Regressor function object for JShower3EZ fit using JGandalf minimiser.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
bool getPMTStatus(const JStatus &status)
Test status of PMT.
JFIT::JHistory JHistory
Definition: JHistory.hh:354
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
Data structure for set of track fit results.
JFit getFit(const int id, const JMODEL::JString &string)
Get fit parameters of string.
Data structure for fit parameters.
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 pDY()
Definition: JShower3Z.hh:172
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
int getCount(const T &hit)
Get hit count.
Data regression method for JFIT::JShower3EZ.
static double Vmax_npe
Maximal integral of PDF [npe].
Data structure for fit of energy.
Definition: JEnergy.hh:28
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:322
double Vmax_npe
maximum number of of photo-electrons
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203
class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA ...
Definition: JShowerFit.hh:72
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
Template specialisation of class JModel to match hit with bright point.
Definition: JFit/JModel.hh:121
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
Match operator for Cherenkov light from shower in any direction.
Basic data structure for L1 hit.
int debug
debug level
double DStep_m
step increase for the distance to optical module [m]