Jpp  16.0.0-rc.2
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:23
3G match criterion.
Definition: JMatch3G.hh:29
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:70
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.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
*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:26
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.
bool operator()(const hit_type &hit) const
Test match.
const int n
Definition: JPolint.hh:676
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.
int debug
debug level
Definition: JSirene.cc:63
Direct access to module in detector data structure.
Data structure for L2 parameters.
JFit getFit(const JHistory &history, const JTrack3D &track, const double Q, const int NDF, const double energy=0.0, const int status=0)
Get fit.
Reduced data structure for L1 hit.
Linear fit of JFIT::JPoint4D.
do JPlot2D f $WORKDIR detector root
JFIT::JHistory JHistory
Definition: JHistory.hh:301
Data structure for set of track fit results.
int getCount(const T &hit)
Get hit count.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
static struct JTRIGGER::JHitR1::@82 compare
Auxiliary data structure for sorting of hits.
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
static struct JTRIGGER::@78 clusterizeWeight
Anonymous struct for weighed clustering of hits.
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]