Jpp  19.1.0-rc.1
the software that should make you happy
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, vcarretero
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, TMaxExtra_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 
104  buildL0(JDAQTimeslice(event, true), router, back_inserter(buffer)); // true => snapshot
105  buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1)); // false => triggered
106 
107  copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));
108 
109  if(dataL1.size() <= numberOfL1){
110  for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
111  if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxExtra_ns)) == dataL1.end()) {
112  dataL0.push_back(*i);
113  }
114  }
115  }
116  for (buffer_type::const_iterator root = dataL1.begin(); root != dataL1.end(); ++root) {
117 
118  buffer_type data(1, *root);
119 
120  JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
121 
122  for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
123 
124  if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
125  data.push_back(*i);
126  }
127  }
128 
129  buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
130 
131  // 4D fit
132  JEstimator_t fit;
133  JPoint4D vx;
134  double chi2 = numeric_limits<double>::max();
135  int NDF = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
136  int N = getCount(data.begin(), __end1);
137 
138  if(NDF > 0){
139  if(distance(data.begin(), __end1) <= factoryLimit){
140 
141  double ymin = numeric_limits<double>::max();
142 
143  buffer_type::iterator __end2 = __end1;
144 
145  for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) >
146  JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
147 
148  sort(data.begin() + 1, __end1, hit_type::compare);
149 
150  do {
151  try {
152 
153  fit(data.begin(), __end2);
154 
155  double y = getChi2(fit, data.begin(), __end2, sigma_ns);
156 
157  if (y < ymin) {
158  ymin = y;
159  vx = fit;
160  chi2 = ymin;
161  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
162  N = getCount(data.begin(), __end2);
163  }
164  }
165  catch(JException& error) { }
166 
167  } while (next_permutation(data.begin() + 1, __end2, __end1, hit_type::compare));
168 
169  ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
170  }
171 
172  } else {
173 
174  const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;
175 
176  buffer_type::iterator __end2 = __end1;
177 
178  for (int n = 0; n <= number_of_outliers; ++n) {
179 
180  try{
181 
182  fit(data.begin(), __end2);
183  vx = fit;
184  chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
185  NDF = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
186  N = getCount(data.begin(), __end2);
187 
188  }
189  catch(const JException& error){ }
190 
191  double ymax = 0;
192  buffer_type::iterator imax = __end2;
193 
194  for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {
195 
196  double y = getChi2(fit, *i, sigma_ns);
197 
198  if (y > ymax) {
199  ymax = y;
200  imax = i;
201  }
202  }
203 
204  if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
205  --__end2;
206  swap(*imax, *__end2);
207  } else {
208  break;
209  }
210  }
211  }
212 
213  out.push_back(getFit(JHistory(JSHOWERPREFIT), JTrack3D(vx, JGEOMETRY3D::JVersor3Z()), getQuality(chi2, N), NDF));
214 
215  }
216  }
217 
219 
220  size_t solutions = out.size();
221 
222  for(size_t i=0; i < solutions; i++){
223  for(int x = -pos_grid_m; x < pos_grid_m + pos_step_m/2.; x += pos_step_m){
224  for(int y = -pos_grid_m; y < pos_grid_m + pos_step_m/2.; y += pos_step_m){
225  for(int z = -pos_grid_m; z < pos_grid_m + pos_step_m/2.; z += pos_step_m){
226  for(int t = -time_grid_ns; t < time_grid_ns + time_step_ns/2.; t += time_step_ns){
227  if (x != 0 || y != 0 || z != 0 || t != 0) {
228 
229  out.push_back(getFit(JHistory(JSHOWERPREFIT),
230  JTrack3D(JPoint4D(getPosition(out[i]) + JPosition3D(x,y,z), out[i].getT()+t), JVersor3Z()),
231  0,
232  0));
233 
234  }
235  }
236  }
237  }
238  }
239  }
240  return out;
241 
242  }
243 
245 
246  /**
247  * Auxiliary class to match hit to root hit.
248  */
249  struct match_t {
250  /**
251  * Constructor.
252  *
253  * \param root root hit
254  * \param TMax_ns maximal time [ns]
255  */
256  match_t(const hit_type& root, const double TMax_ns) :
257  root(root),
258  TMax_ns(TMax_ns)
259  {}
260 
261  /**
262  * Test match.
263  *
264  * \param hit hit
265  * \return true if mach; else false
266  */
267  bool operator()(const hit_type& hit) const
268  {
269  return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
270  }
271 
272  const hit_type& root;
273  double TMax_ns;
274  };
275  };
276 }
277 
278 #endif
279 
Algorithms for hit clustering and sorting.
Linear fit methods.
Auxiliary methods to evaluate Poisson probabilities and chi2.
Reduced data structure for L1 hit.
Match operator for Cherenkov light from shower in any direction.
int debug
debug level
Definition: JSirene.cc:69
Direct access to module in detector data structure.
Linear fit of JFIT::JPoint4D.
Basic data structure for time and time over threshold information of hit.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Router for direct addressing of module data in detector data structure.
Linear fit of bright point (position and time) between hits (objects with position and time).
Data structure for set of track fit results.
void select(const JSelector_t &selector)
Select fits.
Data structure for vertex fit.
Definition: JPoint4D.hh:24
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:41
General exception.
Definition: JException.hh:24
The template JSharedPointer class can be used to share a pointer to an object.
class to handle first step of the shower reconstruction in ORCA: it reconstructs the shower vertex,...
JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
Declaration of the operator that performs the reconstruction.
const JModuleRouter & router
std::vector< hit_type > buffer_type
JShowerPrefit(const JShowerPrefitParameters_t &parameters, const JModuleRouter &router, const int debug=0)
Parameterized constructor.
Auxiliary class to convert binary JMatch operator and given hit to unary match operator.
Definition: JBind2nd.hh:24
Template L0 hit builder.
Definition: JBuildL0.hh:38
Template L2 builder.
Definition: JBuildL2.hh:49
Reduced data structure for L1 hit.
Definition: JHitR1.hh:35
double getT() const
Get calibrated time of hit.
3G match criterion.
Definition: JMatch3G.hh:31
int getModuleID() const
Get module identifier.
static const int JSHOWERPREFIT
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition: JVectorize.hh:261
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Model fits to data.
double getQuality(const double chi2, const int N, const int NDF)
Get quality of fit.
JPosition3D getPosition(const JFit &fit)
Get position.
JFIT::JHistory JHistory
Definition: JHistory.hh:354
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
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.
const int n
Definition: JPolint.hh:786
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
JBinder2nd< JHit_t > JBind2nd(const JMatch< JHit_t > &match, const JHit_t &second)
Auxiliary method to create JBinder2nd object.
Definition: JBind2nd.hh:66
static const struct JTRIGGER::clusterizeWeight clusterizeWeight
Definition: root.py:1
Definition: JSTDTypes.hh:14
int factoryLimit
factory limit for combinatorics
double DMax_m
maximal distance to optical module [m]
double ctMin
minimal cosine space angle between PMT axes
size_t numberOfGrids
number of prefits to be used to build a grid around
double TMaxLocal_ns
time window for local coincidences [ns]
double TMaxExtra_ns
time window for extra coincidences [ns]
Auxiliary class to match hit to root hit.
bool operator()(const hit_type &hit) const
Test match.
match_t(const hit_type &root, const double TMax_ns)
Constructor.
Auxiliary data structure for sorting of hits.
Definition: JHitR1.hh:203
Data structure for L2 parameters.