Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JORCAShowerFit.hh
Go to the documentation of this file.
1 #ifndef JORCASHOWERFIT_INCLUDE
2 #define JORCASHOWERFIT_INCLUDE
3 
4 /*TODO-NOTES :
5  - This algorithm erases element from STL vectors, that is slow!, try to use another STL container as a list.
6  - Prefer to accesses to container elements via at() instead of [], the latter doesn't check the range.
7  - This algorithm uses several magic numbers, that should be avoided and those should clearly be defined.
8  - An exact value of pi should be used.
9  - It would be healthy to check whether some numbers are nan, e.g. , when using getDot.
10 */
11 
12 #include <string>
13 #include <iostream>
14 #include <set>
15 #include <vector>
16 #include <algorithm>
17 #include <memory>
18 
19 #include "JDAQ/JDAQTimeslice.hh"
20 #include "JDAQ/JDAQSummaryslice.hh"
21 
22 #include "JTrigger/JHit.hh"
23 #include "JTrigger/JTimeslice.hh"
24 #include "JTrigger/JHitL0.hh"
25 #include "JTrigger/JHitL1.hh"
26 #include "JTrigger/JHitR1.hh"
27 #include "JTrigger/JBuildL0.hh"
28 #include "JTrigger/JBuildL2.hh"
29 #include "JTrigger/JAlgorithm.hh"
30 #include "JTrigger/JMatch3G.hh"
31 #include "JTrigger/JBind2nd.hh"
32 
33 #include "JFit/JFitStatus.hh"
35 #include "JFit/JHitW0.hh"
36 #include "JFit/JEvt.hh"
37 #include "JFit/JEvtToolkit.hh"
38 #include "JFit/JFitToolkit.hh"
39 #include "JFit/JPoint4D.hh"
40 #include "JFit/JModel.hh"
41 #include "JFit/JSimplex.hh"
44 
45 #include "JDAQ/JDAQTimeslice.hh"
46 #include "JDetector/JDetector.hh"
47 #include "JGeometry3D/JVector3D.hh"
48 #include "JLang/JSharedPointer.hh"
49 
50 
51 namespace JFIT
52 {
53  /**
54  * class to handle the third step of the shower reconstruction, mainly dedicated for ORCA
55  */
56 
58  {
59  private:
60 
62  mutable JRegressor_t fit_;
63 
66 
67  public:
68 
69  /**
70  * Default constructor
71  */
73  {}
74 
75  /**
76  * Parameterized constructor
77  *
78  * \param moduleRouter JSharedPointer of JModuleRouter, this is built via detector file.
79  * \param parameters struct that holds default-optimized parameters for the reconstruction, which can be found in $JPP_DATA.
80  * \param pdfFile PDF file
81  */
82 
84  const JFIT::JShowerFitParameters_t &parameters,
85  const std::string pdfFile):
86  moduleRouter_(moduleRouter),
87  parameters_(parameters)
88  {
90  JRegressor_t::T_ns.setRange(-parameters.Tmax_ns,parameters.Tmax_ns);
91  JRegressor_t::Vmax_npe = 20.0;
92  JRegressor_t::MAXIMUM_ITERATIONS = 1000;
93 
94  fit_ = {pdfFile};
95  }
96 
98 
99 
100  /**
101  * Container of shower hypothesis:
102  * During the fit the algorithm creates a grid in direction and energy space,
103  * for each of these points a shower hypothesis is saved: direction, energy and LogLikelihood value.
104  */
105 
107  {
109  double Energy;
110  double LogLik;
111  };
112 
113  /**
114  * Function to sort different shower hypotheses by their likelihood.
115  */
116 
118  {
119  return(A.LogLik < B.LogLik);
120  };
121 
122  /**
123  * Declaration of the member function that actually performs the reconstruction
124  *
125  * \param timeSliceBuildL0 which is contructed via a JDAQEvent
126  * \param InPreFits input fits, normally produced by JORCAShowerPositionFit
127  * \param OutFits output fits
128  */
129 
130  void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
131  JFIT::JEvt &InPreFits,
132  JFIT::JEvt &OutFits) const;
133 
134  };
135 
136 }
137 
138 
139 /**
140  * member function definition
141  */
142 
143 void
145  JFIT::JEvt &InPreFits,
146  JFIT::JEvt &OutFits) const
147 {
148  using JFIT::JShower3EZ;
149  using JFIT::JPoint4D;
150  using JFIT::JVersor3Z;
152 
153  const double Dmax_m = parameters_.Dmax_m;
154  const double R_Hz = parameters_.R_Hz;
155  const std::size_t numberOfPrefits = parameters_.numberOfPrefits;
156 
157  if ( InPreFits.empty() ) return;
158 
159  typedef std::vector<JTRIGGER::JHitL0> JDataL0_t;
160  typedef std::vector<JFIT::JPMTW0> JPMTW0_t;
162 
163  const JDETECTOR::JDetector &detector = moduleRouter_->getReference();
164 
165  JEvt::iterator __end = InPreFits.end();
166  std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));
167 
168  if (numberOfPrefits > 0) {
169  std::advance(__end = InPreFits.begin(), std::min(numberOfPrefits, InPreFits.size()));
170  }
171 
172  std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
173 
174  JDataL0_t dataL0;
175  buildL0(timeSliceBuildL0, *moduleRouter_, std::back_inserter(dataL0));
176 
177  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
178  const JVector3D pos(getPosition(*shower)); //.rotate(R));
179 
180  JFIT::JPoint4D pt(pos, shower->getT());
181  JVector3D start_dir(0,0,0);
182 
183  const JVector3D DetC(0.0, 0.0, 117);
184  const double radius = pos.getDistance(DetC);
185 
187 
188  double Nhits = 0;
189 
190  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
191 
192  const JFIT::JHitW0 hit(*i, R_Hz);
193  const JGEOMETRY3D::JPosition3D hit_pos(hit.getPosition());
194 
195  const double D = hit_pos.getDistance(pt.getPosition());
196  const double t_res = hit.getT() - pt.getT(hit_pos);
197 
198  const JGEOMETRY3D::JDirection3D photonDir(hit_pos - pt.getPosition());
199  const double cosT = photonDir.getDot(hit.getDirection());
200 
201  if(D < Dmax_m && (t_res >= -25 && t_res <= 25) && (cosT >= -1 && cosT <= 0.1)){
202  top.insert(hit.getPMTIdentifier());
203  const JVector3D d(photonDir.getDX(), photonDir.getDY(), photonDir.getDZ());
204  start_dir.add(d);
205  }
206  if(D < Dmax_m && (t_res > -30 && t_res < 20)){
207  Nhits++;
208  }
209  }
210 
211  double start_E = exp((Nhits + 26)/26);
212  if(start_E > 100) start_E = 100;
213 
214  JFIT::JORCAShowerFit::JShowerHypothesis shower_hypothesis;
216 
218  double chi2 = 0.0;
219  double max_theta = 0, max_phi = 0, scan_step = 0;
220  double Nmin = 0;
221 
222  if(start_E > 10 && radius < 100){
223  max_theta = 30;
224  max_phi = 30;
225  scan_step = 5;
226  Nmin = 10;
227  } else {
228  max_theta = 35;
229  max_phi = 35;
230  scan_step = 5;
231  Nmin = 20;
232  }
233 
234  const JGEOMETRY3D::JAngle3D start_dir_angles(start_dir.getX(), start_dir.getY(), start_dir.getZ());
235 
236  for(double E = -15; E <= 5; E += 7.5){ // scan in energy
237  for(double th = -max_theta; th <= max_theta; th += scan_step){
238  for(double ph = -max_phi; ph <= max_phi; ph += scan_step){
239 
240  buffer.clear();
241  chi2 = 0.0;
242 
243  //todo: use precise value of pi
244  const double theta = th * JTOOLS::PI / 180;
245  const double phi = ph * JTOOLS::PI / 180;
246  const double scan_E = start_E + E; // [GeV]
247 
248  const JGEOMETRY3D::JAngle3D dir_shifted_angles(start_dir_angles.getTheta() + theta,
249  start_dir_angles.getPhi() + phi);
250  const JVector3D dir_shifted(dir_shifted_angles.getDX(),
251  dir_shifted_angles.getDY(),
252  dir_shifted_angles.getDZ());
253  const JGEOMETRY3D::JDirection3D conversion(dir_shifted.getX(),
254  dir_shifted.getY(),
255  dir_shifted.getZ());
256  const JGEOMETRY3D::JRotation3D R(conversion);
257  pt.rotate(R);
258 
259  for (JDETECTOR::JDetector::const_iterator module = detector.begin();
260  module != detector.end(); ++module) {
261 
262  JGEOMETRY3D::JPosition3D pos(module->getPosition());
263  pos.rotate(R);
264 
265  if (pt.getDistance(pos) < Dmax_m) {
266  for (unsigned int i = 0; i != module->size(); ++i) {
267  const KM3NETDAQ::JDAQPMTIdentifier id(module->getID(), i);
268  JDETECTOR::JPMT pmt(module->getPMT(i));
269  pmt.rotate(R);
270  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
271  }
272  }
273  }
274 
275  for(JPMTW0_t::iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt ){
276  chi2 += fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), *pmt);
277  }
278 
279  shower_hypothesis.LogLik = chi2;
280  shower_hypothesis.Energy = scan_E;
281  shower_hypothesis.Direction = dir_shifted;
282 
283  shower_hypotheses.push_back(shower_hypothesis);
284 
285  pt.rotate_back(R);
286 
287  }
288  }
289  }
290 
291  std::vector<JShowerHypothesis>::iterator __end2 = shower_hypotheses.begin() + Nmin;
292  std::partial_sort(shower_hypotheses.begin(), __end2, shower_hypotheses.end(), sortLogLik);
293 
294  for(std::vector<JShowerHypothesis>::iterator hp = shower_hypotheses.begin(); hp != __end2; ++hp){
295 
296  fit_.step.resize(3);
297  const double dir_step = 0.05;
298  fit_.step[0] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(dir_step, 0.0), JFIT::JEnergy());
299  fit_.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JFIT::JEnergy());
300  fit_.step[2] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.05);
301 
302  const double scan_E = hp->Energy;
303 
304  buffer.clear();
305 
306  const JDirection3D startVersor(hp->Direction);
307 
308  const JGEOMETRY3D::JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
309  JGEOMETRY3D::JRotation3D R(conversion);
310  pt.rotate(R);
311 
312  for(JDETECTOR::JDetector::const_iterator module = detector.begin();
313  module != detector.end(); ++module) {
314  JPosition3D pos(module->getPosition());
315  pos.rotate(R);
316 
317  if (pt.getDistance(pos) < Dmax_m) {
318  for (unsigned int i = 0; i != module->size(); ++i) {
319  const KM3NETDAQ::JDAQPMTIdentifier id(module->getID(), i);
320  JDETECTOR::JPMT pmt(module->getPMT(i));
321  pmt.rotate(R);
322  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
323  }
324  }
325  }
326 
327  const int NDF = buffer.size() - fit_.step.size();
328 
329  chi2 = fit_(JShower3EZ(pt, JVersor3Z(), JFIT::JEnergy(log10(scan_E))), buffer.begin(), buffer.end());
330 
331  pt.rotate_back(R);
332 
333  const JVector3D resPos(fit_.value);
334  const JVersor3D resDir(fit_.value.getDX(), fit_.value.getDY(), fit_.value.getDZ());
335 
336  const JGEOMETRY3D::JTrack3D td(resPos, resDir, shower->getT());
337  const double E_reco_corrected = fit_.value.getE() - 4;
338  JGEOMETRY3D::JTrack3E tb(td, E_reco_corrected);
339  tb.rotate_back(R);
340 
341  const JFIT::JFit &outFit = JFIT::getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT),
342  tb, -chi2, NDF,
344 
345  OutFits.push_back(outFit);
346  OutFits.rbegin()->setE(E_reco_corrected);
347  }
348  }
349  std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
350 }
351 
352 #endif
353 
354 
355 
Auxiliary methods to evaluate Poisson probabilities and chi2.
Data structure for angles in three dimensions.
Definition: JAngle3D.hh:30
double getT() const
Get calibrated time of hit.
Definition: JHit.hh:143
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
Definition for fit results, A good fit should hold: OKAY.
Algorithms for hit clustering and sorting.
JFIT::JRegressor< JFIT::JShower3EZ, JFIT::JSimplex > JRegressor_t
const JDAQPMTIdentifier & getPMTIdentifier() const
Get PMT identifier.
Auxiliary class for handling PMT geometry, rate and response.
Definition: JPMTW0.hh:22
Data structure for vertex fit.
Definition: JPoint4D.hh:22
Detector data structure.
Definition: JDetector.hh:77
const JDirection3D & getDirection() const
Get direction.
Rotation matrix.
Definition: JRotation3D.hh:108
static const double PI
Constants.
Definition: JConstants.hh:20
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > moduleRouter_
Container for historical events.
Definition: JHistory.hh:95
class to handle the third step of the shower reconstruction, mainly dedicated for ORCA ...
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:266
3D track with energy.
Definition: JTrack3E.hh:29
Data structure for detector geometry and calibration.
double getDot(const JAngle3D &angle) const
Get dot product.
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
Basic data structure for L0 hit.
Data structure for track fit results.
Definition: JEvt.hh:32
Basic data structure for time and time over threshold information of hit.
Data structure for fit of straight line in positive z-direction with energy.
Definition: JShower3EZ.hh:25
JAxis3D & rotate(const JRotation3D &R)
Rotate axis.
Definition: JAxis3D.hh:235
Data structure for vector in three dimensions.
Definition: JVector3D.hh:32
static bool sortLogLik(JShowerHypothesis A, JShowerHypothesis B)
Function to sort different shower hypotheses by their likelihood.
double getDY() const
Get y direction.
Definition: JVersor3D.hh:101
double getDX() const
Get x direction.
Definition: JVersor3D.hh:90
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:52
Auxiliary class for a hit with background rate value.
Definition: JHitW0.hh:21
JFIT::JShowerFitParameters_t parameters_
void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0, JFIT::JEvt &InPreFits, JFIT::JEvt &OutFits) const
Declaration of the member function that actually performs the reconstruction.
Data time slice.
JAxis3D & rotate_back(const JRotation3D &R)
Rotate back axis.
Definition: JAxis3D.hh:250
int debug
debug level
Definition: JSirene.cc:59
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
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
Reduced data structure for L1 hit.
JPosition3D getPosition(const JFit &fit)
Get position.
Definition: JEvtToolkit.hh:63
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
Data structure for set of track fit results.
Definition: JEvt.hh:312
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:40
Data regression method for JFIT::JShower3Z.
Data structure for fit of energy.
Definition: JEnergy.hh:24
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Definition: JEvtToolkit.hh:223
JORCAShowerFit(const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &moduleRouter, const JFIT::JShowerFitParameters_t &parameters, const std::string pdfFile)
Parameterized constructor.
Template L0 hit builder.
Definition: JBuildL0.hh:35
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:36
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:112
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:140
Match operator for Cherenkov light from shower in any direction.
struct that holds the Parameters used for JORCAShowerFit
Basic data structure for L1 hit.
JORCAShowerFit()
Default constructor.
Container of shower hypothesis: During the fit the algorithm creates a grid in direction and energy s...