Jpp  18.2.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 "JFit/JFitToolkit.hh"
22 #include "JFit/JEstimator.hh"
23 #include "JFit/JPoint4D.hh"
24 
26 
27 #include "JReconstruction/JEvt.hh"
30 
31 #include "JLang/JComparator.hh"
32 
33 #include "JGeometry3D/JTrack3D.hh"
34 
35 /**
36  * \author adomi, gmaggi
37  */
38 
39 namespace JRECONSTRUCTION
40 {
41 
44 
45  /**
46  * class to handle first step of the shower reconstruction in ORCA:
47  * it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits
48  */
49 
51  {
52 
53  public:
54 
57 
58  /**
59  * Parameterized constructor
60  *
61  * \param parameters struct that holds default-optimized parameters for the reconstruction.
62  * \param router module router, this is built via detector file.
63  * \param debug debug
64  */
65 
67  const JModuleRouter& router,
68  const int debug = 0):
69  JShowerPrefitParameters_t(parameters),
70  router(router)
71  {}
72 
73  /**
74  * Declaration of the operator that performs the reconstruction
75  *
76  * \param event which is a JDAQEvent
77  */
78 
79  JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
80  {
81 
82  using namespace std;
83  using namespace JPP;
84 
85  const double STANDARD_DEVIATIONS = 3.0;
86 
87  typedef JEstimator<JPoint4D> JEstimator_t;
88 
89  const JBuildL0<hit_type> buildL0;
91  const JMatch3G<hit_type> match3G(roadWidth_m, TMaxLocal_ns); // causality relation for showers
92 
93  //define empty vectors to be filled with JHitL0 and JHitL1
94  buffer_type dataL0;
95  buffer_type dataL1;
96 
97  JEvt out;
98 
99  buffer_type buffer;
100 
101  buildL0(JDAQTimeslice(event, true), router, back_inserter(buffer)); // true => snapshot
102  buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1)); // false => triggered
103 
104  copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));
105 
106  for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
107 
108  if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxLocal_ns)) == dataL1.end()) {
109  dataL0.push_back(*i);
110  }
111  }
112 
113  for (buffer_type::const_iterator root = dataL1.begin(); root != dataL1.end(); ++root) {
114 
115  buffer_type data(1, *root);
116 
117  JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
118 
119  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
120 
121  if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
122  data.push_back(*i);
123  }
124  }
125 
126  buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
127 
128  // 4D fit
129  JEstimator_t fit;
130  JPoint4D vx;
131  double chi2 = numeric_limits<double>::max();
132  int NDF = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
133  int N = getCount(data.begin(), __end1);
134 
135  if(distance(data.begin(), __end1) <= factoryLimit){
136 
137  double ymin = numeric_limits<double>::max();
138 
139  buffer_type::iterator __end2 = __end1;
140 
141  for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) >=
142  JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
143 
144  sort(data.begin() + 1, __end1, hit_type::compare);
145 
146  do {
147  try {
148 
149  fit(data.begin(), __end2);
150 
151  double y = getChi2(fit, data.begin(), __end2, sigma_ns);
152 
153  if (y < ymin) {
154  ymin = y;
155  vx = fit;
156  chi2 = ymin;
157  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
158  N = getCount(data.begin(), __end2);
159  }
160  }
161  catch(JException& error) { }
162 
163  } while (next_permutation(data.begin() + 1, __end2, __end1, hit_type::compare));
164 
165  ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
166  }
167 
168  } else {
169 
170  const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;
171 
172  buffer_type::iterator __end2 = __end1;
173 
174  for (int n = 0; n <= number_of_outliers; ++n) {
175 
176  try{
177 
178  fit(data.begin(), __end2);
179  vx = fit;
180  chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
181  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
182  N = getCount(data.begin(), __end2);
183 
184  }
185  catch(const JException& error){ }
186 
187  double ymax = 0;
188  buffer_type::iterator imax = __end2;
189 
190  for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {
191 
192  double y = getChi2(fit, *i, sigma_ns);
193 
194  if (y > ymax) {
195  ymax = y;
196  imax = i;
197  }
198  }
199 
200  if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
201  --__end2;
202  swap(*imax, *__end2);
203  } else {
204  break;
205  }
206  }
207  }
208 
209  if (NDF >= 0 && chi2 > numeric_limits<double>::lowest()) {
210 
211  const JGEOMETRY3D::JTrack3D fitResult(vx, JGEOMETRY3D::JVersor3Z());
212 
213  out.push_back(getFit(JHistory(JSHOWERPREFIT), fitResult, getQuality(chi2, N), NDF));
214 
215  out.rbegin()->setW(13, chi2);
216  out.rbegin()->setW(14, N);
217  }
218  }
219 
220  return out;
221 
222  }
223 
225 
226  /**
227  * Auxiliary class to match hit to root hit.
228  */
229  struct match_t {
230  /**
231  * Constructor.
232  *
233  * \param root root hit
234  * \param TMax_ns maximal time [ns]
235  */
236  match_t(const hit_type& root, const double TMax_ns) :
237  root(root),
238  TMax_ns(TMax_ns)
239  {}
240 
241  /**
242  * Test match.
243  *
244  * \param hit hit
245  * \return true if mach; else false
246  */
247  bool operator()(const hit_type& hit) const
248  {
249  return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
250  }
251 
252  const hit_type& root;
253  double TMax_ns;
254  };
255  };
256 }
257 
258 #endif
259 
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.
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:697
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
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.
do JPlot2D f $WORKDIR canberra[${EMITTER}] root
Definition: JCanberra.sh:132
double TMaxLocal_ns
time window for local coincidences [ns]
int debug
debug level