Jpp
JORCAShowerPositionFit.hh
Go to the documentation of this file.
1 #ifndef JORCASHOWERPOSITIONFIT_INCLUDE
2 #define JORCASHOWERPOSITIONFIT_INCLUDE
3 
4 /*TODO-suggested
5  This code erases elements from a std::vector, that is slow!, try to use another STL container
6 */
7 
8 #include <string>
9 #include <iostream>
10 #include <vector>
11 #include <algorithm>
12 #include <memory>
13 
14 #include "TFile.h"
15 
16 #include "JDAQ/JDAQTimeslice.hh"
17 #include "JDAQ/JDAQSummaryslice.hh"
18 
19 #include "JTrigger/JHit.hh"
20 #include "JTrigger/JTimeslice.hh"
21 #include "JTrigger/JHitL0.hh"
22 #include "JTrigger/JHitL1.hh"
23 #include "JTrigger/JHitR1.hh"
24 #include "JTrigger/JBuildL0.hh"
25 #include "JTrigger/JBuildL2.hh"
26 #include "JTrigger/JAlgorithm.hh"
27 #include "JTrigger/JMatch3G.hh"
28 #include "JTrigger/JBind2nd.hh"
29 
30 #include "JFit/JFitStatus.hh"
31 #include "JFit/JEvt.hh"
32 #include "JFit/JEvtToolkit.hh"
33 #include "JFit/JFitToolkit.hh"
36 #include "JFit/JEstimator.hh"
37 #include "JFit/JPoint4D.hh"
38 #include "JFit/JModel.hh"
39 #include "JFit/JSimplex.hh"
42 
44 #include "JGeometry3D/JVector3D.hh"
45 #include "JTools/JPermutation.hh"
46 #include "JTools/JRange.hh"
47 #include "JLang/JSharedPointer.hh"
48 
49 /**
50  * \author adomi, gmaggi
51 */
52 
53 namespace JFIT
54 {
55  /**
56  * class to handle the second step of the shower reconstruction, mainly dedicated for ORCA
57  */
59  {
60  private:
61 
63  std::shared_ptr<JRegressor_t> fit0_;
64  std::shared_ptr<JRegressor_t> fit_;
65 
68 
69  public:
70  /**
71  * Default constructor
72  */
74  {}
75 
76  /**
77  * Parameterized constructor
78  *
79  * \param moduleRouter JSharedPointer of JModuleRouter, this is built via detector file.
80  * \param parameters struct that holds default-optimized parameters for the reconstruction, which can be found in $JPP_DATA.
81  * \param pdfFile file name with PDF
82  */
83 
85  const JShowerPositionFitParameters_t &parameters,
86  const std::string pdfFile):
87  moduleRouter_(moduleRouter),
88  parameters_(parameters)
89  {
90  TFile* pdf_file = new TFile(pdfFile.c_str());
91  TH2D* hpdf = (TH2D*)pdf_file->Get("hPDF2Dist");
92 
94  JRegressor_t::MAXIMUM_ITERATIONS = 10000;
95 
96  fit0_ = std::make_shared<JRegressor_t>(true);
97  fit0_->estimator.reset(new JMEstimatorTukey(10.0));
98 
99  fit_ = std::make_shared<JRegressor_t>(hpdf, false);
100  }
101 
103 
104  /**
105  * Container of vertex hypothesis:
106  * During the fit the algorithm creates a grid in position space (x, y, z, t),
107  * for each of these points a vertex hypothesis is saved: 4D-Point and LogLikelihood value.
108  */
109 
111  {
113  double LogLik;
114  };
115 
116 
117  /**
118  * Function to sort different vertex hypotheses by their likelihood.
119  */
120 
122  {
123  return(A.LogLik < B.LogLik);
124  }
125 
126  /**
127  * Declaration of Member function that actually performs the reconstruction
128  *
129  * \param timeSliceBuildL0 which is contructed via a JDAQEvent
130  * \param timeSliceBuildL2 which is contructed via a JDAQEvent
131  * \param InPreFits input fits (normally made by JORCAShowerPrefit)
132  * \param OutFits output fits
133  */
134 
135  void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
136  const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2,
137  JFIT::JEvt &InPreFits,
138  JFIT::JEvt &OutFits) const
139  {
140 
141  if ( InPreFits.empty() ) return;
142 
143  const int numberOfPrefits = parameters_.numberOfPrefits;
144  const double Tmax_ns = parameters_.Tmax_ns; // [ns]
145  const double ctMin = parameters_.ctMin; // optimal value
146  const double Dmax_m = parameters_.Dmax_m; // [m] optimal value
148  const double pos_step = parameters_.pos_step; // [m] density of the grid
149  const JTOOLS::JRange<double> time_grid = parameters_.time_grid;
150  const double time_step = parameters_.time_step; // [ns] density of the grid
151 
152  typedef std::vector<JTRIGGER::JHitL0> JDataL0_t;
153  typedef std::vector<JTRIGGER::JHitL1> JDataL1_t;
154 
156  const JTRIGGER::JBuildL2<JTRIGGER::JHitL1> buildL2(JTRIGGER::JL2Parameters(3, Tmax_ns, ctMin));
157  const JTRIGGER::JMatch3G<JTRIGGER::JHitL1> match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS
158 
159 
160  JEvt::iterator __end = InPreFits.end();
161  if (numberOfPrefits > 0) {
162  std::advance(__end = InPreFits.begin(), std::min(static_cast<std::size_t>(numberOfPrefits), InPreFits.size()));
163  }
164  std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
165 
166  if(!InPreFits.empty()) {
167 
168  std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));
169 
170  JDataL0_t dataL0;
171  JDataL0_t dataL0_selected;
172  JDataL1_t data;
173 
174  buildL0(timeSliceBuildL0, *moduleRouter_, std::back_inserter(dataL0));
175  buildL2(timeSliceBuildL2,*moduleRouter_, std::back_inserter(data));
176 
177  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
178 
179  const JFIT::JPoint4D vertex(JFIT::getPosition(*shower), shower->getT());
180 
181  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
182 
183  const JTRIGGER::JHitL0 hit(*i);
184  const JGEOMETRY3D::JPosition3D hit_pos(hit.getPosition());
185  const double D = hit_pos.getDistance(vertex);
186  const double t_res = hit.getT() - vertex.getT(hit_pos);
187  const JGEOMETRY3D::JDirection3D photonDir(hit_pos - vertex.getPosition());
188  const double cosT = photonDir.getDot(hit.getDirection());
189 
190  if(data.size()>0){
191  if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
192  dataL0_selected.push_back(hit);
193  }
194  } else {
195  if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
196  dataL0_selected.push_back(hit);
197  }
198  }
199 
200  }
201 
202  if(dataL0_selected.size() != 0){
203 
204  double chi2 = 0.0;
205 
208 
209  for(double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit();
210  x += pos_step){
211  for(double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit();
212  y += pos_step){
213  for(double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit();
214  z += pos_step){
215  for(double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit();
216  t += time_step){
217 
218  chi2 = 0.0;
219 
220  const JGEOMETRY3D::JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
221  const JFIT::JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
222 
223  for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
224  const JHitL0 hit(*i);
225  chi2 += (*fit0_)(point_shifted, hit);
226  }
227 
228  vertex_hypothesis.LogLik = chi2;
229  vertex_hypothesis.Vertex = point_shifted;
230 
231  vertex_hypotheses.push_back(vertex_hypothesis);
232 
233  }
234  }
235  }
236  }
237 
238  std::vector<JVertexHypothesis>::iterator __end2 = vertex_hypotheses.begin() + 35;
239  std::partial_sort(vertex_hypotheses.begin(), __end2, vertex_hypotheses.end(), sortLogLik);
240 
241  for(std::vector<JVertexHypothesis>::iterator hp = vertex_hypotheses.begin(); hp != __end2; ++hp){
242 
243  double simplex_step = 1;
244  fit0_->step.resize(4);
245  fit0_->step[0] = JFIT::JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
246  fit0_->step[1] = JFIT::JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
247  fit0_->step[2] = JFIT::JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
248  fit0_->step[3] = JFIT::JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);
249 
250  chi2 = (*fit0_)(JPoint4D(hp->Vertex), dataL0_selected.begin(), dataL0_selected.end());
251 
252  const JGEOMETRY3D::JVector3D fitPos0(fit0_->value);
253  const JFIT::JPoint4D pt(fitPos0, fit0_->value.getT());
254 
255  simplex_step = 0.3;
256  fit_->step.resize(4);
257  fit_->step[0] = JFIT::JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
258  fit_->step[1] = JFIT::JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
259  fit_->step[2] = JFIT::JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
260  fit_->step[3] = JFIT::JPoint4D(JVector3D(0.0, 0.0, 0.0), 10);
261 
262  chi2 = (*fit_)(pt, dataL0_selected.begin(), dataL0_selected.end());
263 
264  const int NDF = dataL0_selected.size() - fit_->step.size();
265 
266  const JGEOMETRY3D::JVector3D fitPos(fit_->value);
267  const JGEOMETRY3D::JTrack3D fitResult(fitPos, JVersor3Z(0,0), fit_->value.getT());
268 
269  const double energy(0);
270  const JFIT::JFit &outFit = JFIT::getFit(JFIT::JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT),
271  fitResult, getQuality(chi2, NDF), NDF,
272  energy,JFIT::JFitStatus_t::OKAY);
273 
274  OutFits.push_back(outFit);
275  OutFits.rbegin()->setW(13, chi2);
276  OutFits.rbegin()->setW(14, NDF);
277  }
278  } else{
279  std::cout<<"Too few hits " << std::endl;
280  }
281  }
282  std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
283  }
284 
285  }
286 
287  };
288 }
289 
290 #endif
291 
292 
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:215
JFIT::JORCAShowerPositionFit::fit_
std::shared_ptr< JRegressor_t > fit_
Definition: JORCAShowerPositionFit.hh:64
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:226
JVector3D.hh
JFIT::JORCAShowerPositionFit
class to handle the second step of the shower reconstruction, mainly dedicated for ORCA
Definition: JORCAShowerPositionFit.hh:58
JFIT
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
JTRIGGER::JL2Parameters
Data structure for L2 parameters.
Definition: JTriggerParameters.hh:33
JFIT::JORCAShowerPositionFit::fit0_
std::shared_ptr< JRegressor_t > fit0_
Definition: JORCAShowerPositionFit.hh:63
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JCounter.hh:35
JPermutation.hh
JFitToolkit.hh
JGEOMETRY3D::JVersor3Z
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:37
JPoint4D.hh
JFIT::JORCAShowerPositionFit::moduleRouter_
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > moduleRouter_
Definition: JORCAShowerPositionFit.hh:66
JFitStatus.hh
JGEOMETRY3D::JDirection3D::getDirection
const JDirection3D & getDirection() const
Get direction.
Definition: JDirection3D.hh:106
JTimeslice.hh
JFIT::JShowerPositionFitParameters_t::pos_grid
JTOOLS::JRange< double > pos_grid
Definition: JShowerPositionFitParameters_t.hh:29
JFIT::JSHOWERPOSITIONFIT
JShowerPositionFit.cc.
Definition: JFitApplications.hh:33
JGEOMETRY3D::JVector3D::getDistance
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:268
JSharedPointer.hh
JGEOMETRY3D::JDirection3D::getDot
double getDot(const JAngle3D &angle) const
Get dot product.
Definition: JDirection3D.hh:333
JSimplex.hh
std::vector
Definition: JSTDTypes.hh:12
JPoint3DEstimator.hh
JTOOLS::JRange< double >
JPoint4DEstimator.hh
JEvtToolkit.hh
JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
KM3NETDAQ::JDAQTimeslice
Data time slice.
Definition: JDAQTimeslice.hh:36
JFIT::JShowerPositionFitParameters_t::numberOfPrefits
int numberOfPrefits
Definition: JShowerPositionFitParameters_t.hh:25
JPoint4DRegressor.hh
JFIT::JShowerPositionFitParameters_t::Tmax_ns
double Tmax_ns
Definition: JShowerPositionFitParameters_t.hh:26
JDAQTimeslice.hh
JFIT::JEvt
Data structure for set of track fit results.
Definition: JEvt.hh:292
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JDAQSummaryslice.hh
JTRIGGER::JHit::getT
double getT() const
Get calibrated time of hit.
Definition: JHit.hh:143
JShowerParameters.hh
JEvt.hh
JFIT::JORCAShowerPositionFit::JVertexHypothesis::LogLik
double LogLik
Definition: JORCAShowerPositionFit.hh:113
JTRIGGER::JBuildL0
Template L0 hit builder.
Definition: JBuildL0.hh:34
JRange.hh
JBuildL0.hh
JFIT::JShowerPositionFitParameters_t
Definition: JShowerPositionFitParameters_t.hh:22
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition: JVector3D.hh:33
JFIT::qualitySorter
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Definition: JEvtToolkit.hh:223
JMatch3G.hh
JEstimator.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JFIT::JORCAShowerPositionFit::~JORCAShowerPositionFit
~JORCAShowerPositionFit()
Definition: JORCAShowerPositionFit.hh:102
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JFIT::JORCAShowerPositionFit::getJEvt
void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0, const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2, JFIT::JEvt &InPreFits, JFIT::JEvt &OutFits) const
Declaration of Member function that actually performs the reconstruction.
Definition: JORCAShowerPositionFit.hh:135
JFIT::JORCAShowerPositionFit::JVertexHypothesis::Vertex
JPoint4D Vertex
Definition: JORCAShowerPositionFit.hh:112
JFIT::JORCAShowerPositionFit::sortLogLik
static bool sortLogLik(JVertexHypothesis A, JVertexHypothesis B)
Function to sort different vertex hypotheses by their likelihood.
Definition: JORCAShowerPositionFit.hh:121
JFIT::JORCAShowerPositionFit::JORCAShowerPositionFit
JORCAShowerPositionFit()
Default constructor.
Definition: JORCAShowerPositionFit.hh:73
JHit.hh
JFIT::JRegressor
Template definition of a data regressor of given model.
Definition: JRegressor.hh:68
JModuleRouter.hh
JBind2nd.hh
JFIT::JShowerPositionFitParameters_t::ctMin
double ctMin
Definition: JShowerPositionFitParameters_t.hh:27
JHitR1.hh
JTRIGGER::JMatch3G
3G match criterion.
Definition: JMatch3G.hh:29
JBuildL2.hh
JFIT::JPoint4D
Data structure for vertex fit.
Definition: JPoint4D.hh:22
JFIT::JShowerPositionFitParameters_t::Dmax_m
double Dmax_m
Definition: JShowerPositionFitParameters_t.hh:28
JHitL1.hh
JFIT::JORCAShowerPositionFit::JVertexHypothesis
Container of vertex hypothesis: During the fit the algorithm creates a grid in position space (x,...
Definition: JORCAShowerPositionFit.hh:110
JModel.hh
JFIT::JMEstimatorTukey
Tukey's biweight M-estimator.
Definition: JMEstimator.hh:105
JFIT::JShowerPositionFitParameters_t::time_grid
JTOOLS::JRange< double > time_grid
Definition: JShowerPositionFitParameters_t.hh:30
JFIT::JShowerPositionFitParameters_t::time_step
double time_step
Definition: JShowerPositionFitParameters_t.hh:32
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
JTRIGGER::JHitL0
Data structure for L0 hit.
Definition: JHitL0.hh:27
JFIT::JORCAShowerPositionFit::parameters_
JShowerPositionFitParameters_t parameters_
Definition: JORCAShowerPositionFit.hh:67
JGEOMETRY3D::JTrack3D
3D track.
Definition: JTrack3D.hh:30
JFIT::JORCAShowerPositionFit::JORCAShowerPositionFit
JORCAShowerPositionFit(const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &moduleRouter, const JShowerPositionFitParameters_t &parameters, const std::string pdfFile)
Parameterized constructor.
Definition: JORCAShowerPositionFit.hh:84
JFIT::JShowerPositionFitParameters_t::pos_step
double pos_step
Definition: JShowerPositionFitParameters_t.hh:31
JFIT::getFit
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.
Definition: JEvtToolkit.hh:116
JFIT::JHistory
Container for historical events.
Definition: JHistory.hh:95
JHitL0.hh
JFIT::getQuality
double getQuality(const double chi2, const int NDF)
Get quality of fit.
Definition: JEvtToolkit.hh:195
JFIT::getPosition
JPosition3D getPosition(const JFit &fit)
Get position.
Definition: JEvtToolkit.hh:63
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter >
JFIT::JFit
Data structure for track fit results.
Definition: JEvt.hh:31
JFIT::OKAY
Definition: JFitStatus.hh:22
JTRIGGER::JBuildL2
Template L2 builder.
Definition: JBuildL2.hh:45
JFIT::JORCAShowerPositionFit::JRegressor_t
JFIT::JRegressor< JFIT::JPoint4D, JFIT::JSimplex > JRegressor_t
Definition: JORCAShowerPositionFit.hh:62
JAlgorithm.hh