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
JShowerPrefit.hh
Go to the documentation of this file.
1 #ifndef JSHOWERPREFIT_INCLUDE
2 #define JSHOWERPREFIT_INCLUDE
3 
4 #include <iostream>
5 #include <vector>
6 #include <algorithm>
7 #include <functional>
8 
11 
12 #include "JTrigger/JHit.hh"
13 #include "JTrigger/JHitR1.hh"
14 #include "JTrigger/JBuildL0.hh"
15 #include "JTrigger/JBuildL2.hh"
16 #include "JTrigger/JAlgorithm.hh"
17 #include "JTrigger/JMatch3G.hh"
18 #include "JTrigger/JBind2nd.hh"
19 
20 #include "JTools/JPermutation.hh"
21 
22 #include "JFit/JFitToolkit.hh"
24 #include "JFit/JEstimator.hh"
25 #include "JFit/JPoint4D.hh"
26 
28 
29 #include "JReconstruction/JEvt.hh"
32 
33 #include "JLang/JComparator.hh"
34 
35 #include "JGeometry3D/JTrack3D.hh"
36 
37 /**
38  * \author adomi, gmaggi
39  */
40 
41 namespace JRECONSTRUCTION
42 {
43 
46 
47  /**
48  * class to handle first step of the shower reconstruction in ORCA:
49  * it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits
50  */
51 
53  {
54 
55  public:
56 
59 
60  /**
61  * Parameterized constructor
62  *
63  * \param parameters struct that holds default-optimized parameters for the reconstruction.
64  * \param router module router, this is built via detector file.
65  * \param debug debug
66  */
67 
69  const JModuleRouter& router,
70  const int debug = 0):
71  JShowerPrefitParameters_t(parameters),
72  router(router)
73  {}
74 
75  /**
76  * Declaration of the operator that performs the reconstruction
77  *
78  * \param event which is a JDAQEvent
79  */
80 
81  JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
82  {
83 
84  using namespace std;
85  using namespace JPP;
86 
87  const double STANDARD_DEVIATIONS = 3.0;
88 
89  typedef JEstimator<JPoint4D> JEstimator_t;
90 
91  const JBuildL0<hit_type> buildL0;
93  const JMatch3G<hit_type> match3G(DMax_m, TMaxLocal_ns); // causality relation for showers
94 
95  //define empty vectors to be filled with JHitL0 and JHitL1
96  buffer_type dataL0;
97  buffer_type dataL1;
98 
99  JEvt out;
100 
101  buffer_type buffer;
102 
103  buildL0(JDAQTimeslice(event, true), router, back_inserter(buffer)); // true => snapshot
104  buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1)); // false => triggered
105 
106  copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));
107 
108  for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
109 
110  if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxLocal_ns)) == dataL1.end()) {
111  dataL0.push_back(*i);
112  }
113  }
114 
115  for (buffer_type::const_iterator root = dataL1.begin(); root != dataL1.end(); ++root) {
116 
117  buffer_type data(1, *root);
118 
119  JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
120 
121  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
122 
123  if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
124  data.push_back(*i);
125  }
126  }
127 
128  buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
129 
130  // 4D fit
131  JEstimator_t fit;
132  JPoint4D vx;
133  double chi2 = numeric_limits<double>::max();
134  int NDF = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
135  int N = getCount(data.begin(), __end1);
136 
137  if(distance(data.begin(), __end1) <= factoryLimit){
138 
139  double ymin = numeric_limits<double>::max();
140 
141  buffer_type::iterator __end2 = __end1;
142 
143  for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) >=
144  JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
145 
146  sort(data.begin() + 1, __end1, hit_type::compare);
147 
148  do {
149  try {
150 
151  fit(data.begin(), __end2);
152 
153  double y = getChi2(fit, data.begin(), __end2, sigma_ns);
154 
155  if (y < ymin) {
156  ymin = y;
157  vx = fit;
158  chi2 = ymin;
159  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
160  N = getCount(data.begin(), __end2);
161  }
162  }
163  catch(JException& error) { }
164 
165  } while (next_permutation(data.begin() + 1, __end2, __end1, hit_type::compare));
166 
167  ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
168  }
169 
170  } else {
171 
172  const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;
173 
174  buffer_type::iterator __end2 = __end1;
175 
176  for (int n = 0; n <= number_of_outliers; ++n) {
177 
178  try{
179 
180  fit(data.begin(), __end2);
181  vx = fit;
182  chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
183  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
184  N = getCount(data.begin(), __end2);
185 
186  }
187  catch(const JException& error){ }
188 
189  double ymax = 0;
190  buffer_type::iterator imax = __end2;
191 
192  for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {
193 
194  double y = getChi2(fit, *i, sigma_ns);
195 
196  if (y > ymax) {
197  ymax = y;
198  imax = i;
199  }
200  }
201 
202  if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
203  --__end2;
204  swap(*imax, *__end2);
205  } else {
206  break;
207  }
208  }
209  }
210 
211  if (NDF >= 0 && chi2 > numeric_limits<double>::lowest()) {
212 
213  const JGEOMETRY3D::JTrack3D fitResult(vx, JGEOMETRY3D::JVersor3Z());
214 
215  out.push_back(getFit(JHistory(JSHOWERPREFIT), fitResult, getQuality(chi2, N), NDF));
216 
217  out.rbegin()->setW(13, chi2);
218  out.rbegin()->setW(14, N);
219  }
220  }
221 
222  return out;
223 
224  }
225 
227 
228  /**
229  * Auxiliary class to match hit to root hit.
230  */
231  struct match_t {
232  /**
233  * Constructor.
234  *
235  * \param root root hit
236  * \param TMax_ns maximal time [ns]
237  */
238  match_t(const hit_type& root, const double TMax_ns) :
239  root(root),
240  TMax_ns(TMax_ns)
241  {}
242 
243  /**
244  * Test match.
245  *
246  * \param hit hit
247  * \return true if mach; else false
248  */
249  bool operator()(const hit_type& hit) const
250  {
251  return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
252  }
253 
254  const hit_type& root;
255  double TMax_ns;
256  };
257  };
258 }
259 
260 #endif
261 
Auxiliary methods to evaluate Poisson probabilities and chi2.
int factoryLimit
factory limit for combinatorics
Linear fit methods.
double getT() const
Get calibrated time of hit.
General exception.
Definition: JException.hh:24
3G match criterion.
Definition: JMatch3G.hh:29
static struct JTRIGGER::clusterizeWeight clusterizeWeight
int getModuleID() const
Get module identifier.
JBinder2nd< JHit_t > JBind2nd(const JMatch< JHit_t > &match, const JHit_t &second)
Auxiliary method to create JBinder2nd object.
Definition: JBind2nd.hh:66
Linear fit of bright point (position and time) between hits (objects with position and time)...
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.
Data structure for vertex fit.
Definition: JPoint4D.hh:22
Router for direct addressing of module data in detector data structure.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=SINGLE_STAGE)
Get fit.
*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
Auxiliary class to convert binary JMatch operator and given hit to unary match operator.
Definition: JBind2nd.hh:24
class to handle first step of the shower reconstruction in ORCA: it reconstructs the shower vertex...
Basic data structure for time and time over threshold information of hit.
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
std::vector< hit_type > buffer_type
JShowerPrefit(const JShowerPrefitParameters_t &parameters, const JModuleRouter &router, const int debug=0)
Parameterized constructor.
do JPlot2D f $WORKDIR canberra[${EMITTER}\] root
JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
Declaration of the operator that performs the reconstruction.
static struct JTRIGGER::JHitR1::compare compare
bool operator()(const hit_type &hit) const
Test match.
const int n
Definition: JPolint.hh:786
static const int JSHOWERPREFIT
Template L2 builder.
Definition: JBuildL2.hh:45
The template JSharedPointer class can be used to share a pointer to an object.
const JModuleRouter & router
double ctMin
minimal cosine space angle between PMT axes
Data time slice.
Direct access to module in detector data structure.
Data structure for L2 parameters.
Reduced data structure for L1 hit.
Linear fit of JFIT::JPoint4D.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
JFIT::JHistory JHistory
Definition: JHistory.hh:354
Data structure for set of track fit results.
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.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
double DMax_m
maximal distance to optical module [m]
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
bool next_permutation(T __begin, T __last, T __end, JComparator_t compare, std::bidirectional_iterator_tag)
Implementation of method next_permutation for bidirectional iterators.
Definition: JPermutation.hh:20
Template L0 hit builder.
Definition: JBuildL0.hh:35
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
match_t(const hit_type &root, const double TMax_ns)
Constructor.
Auxiliary class to match hit to root hit.
Match operator for Cherenkov light from shower in any direction.
double TMaxLocal_ns
time window for local coincidences [ns]
int debug
debug level