Jpp
Classes | Public Member Functions | Static Public Member Functions | Private Types | Private Attributes | List of all members
JFIT::JORCAShowerFit Class Reference

#include <JORCAShowerFit.hh>

Classes

struct  JShowerHypothesis
 Container of shower hypothesis: During the fit the algorithm creates a grid in direction and energy space, for each of these points a shower hypothesis is saved: direction, energy and LogLikelihood value. More...
 

Public Member Functions

 JORCAShowerFit ()
 Default constructor. More...
 
 JORCAShowerFit (const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &moduleRouter, const JFIT::JShowerFitParameters_t &parameters, const std::string pdfFile)
 Parameterized constructor. More...
 
 ~JORCAShowerFit ()
 
void getJEvt (const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0, JFIT::JEvt &InPreFits, JFIT::JEvt &OutFits) const
 Declaration of the member function that actually performs the reconstruction. More...
 

Static Public Member Functions

static double estimateInitialEnergy (int nHits)
 Function to estimate the initial em shower energy. More...
 
static bool sortLogLik (JShowerHypothesis A, JShowerHypothesis B)
 Function to sort different shower hypotheses by their likelihood. More...
 

Private Types

typedef JFIT::JRegressor< JFIT::JShower3EZ, JFIT::JSimplexJRegressor_t
 

Private Attributes

JRegressor_t fit_
 
JLANG::JSharedPointer< const JDETECTOR::JModuleRoutermoduleRouter_
 
JFIT::JShowerFitParameters_t parameters_
 

Detailed Description

Definition at line 64 of file JORCAShowerFit.hh.

Member Typedef Documentation

◆ JRegressor_t

Definition at line 68 of file JORCAShowerFit.hh.

Constructor & Destructor Documentation

◆ JORCAShowerFit() [1/2]

JFIT::JORCAShowerFit::JORCAShowerFit ( )
inline

Default constructor.

Definition at line 79 of file JORCAShowerFit.hh.

80  {}

◆ JORCAShowerFit() [2/2]

JFIT::JORCAShowerFit::JORCAShowerFit ( const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &  moduleRouter,
const JFIT::JShowerFitParameters_t parameters,
const std::string  pdfFile 
)
inline

Parameterized constructor.

Parameters
moduleRouterJSharedPointer of JModuleRouter, this is built via detector file.
parametersstruct that holds default-optimized parameters for the reconstruction, which can be found in $JPP_DATA.
pdfFilePDF file

Definition at line 90 of file JORCAShowerFit.hh.

92  :
93  moduleRouter_(moduleRouter),
94  parameters_(parameters)
95  {
97  JRegressor_t::T_ns.setRange(-parameters.Tmax_ns,parameters.Tmax_ns);
98  JRegressor_t::Vmax_npe = 20.0;
99  JRegressor_t::MAXIMUM_ITERATIONS = 1000;
100 
101  fit_ = {pdfFile};
102  fit_.estimator.reset(new JMEstimatorNull());
103  }

◆ ~JORCAShowerFit()

JFIT::JORCAShowerFit::~JORCAShowerFit ( )
inline

Definition at line 105 of file JORCAShowerFit.hh.

105 {}

Member Function Documentation

◆ estimateInitialEnergy()

static double JFIT::JORCAShowerFit::estimateInitialEnergy ( int  nHits)
inlinestatic

Function to estimate the initial em shower energy.

Definition at line 125 of file JORCAShowerFit.hh.

125  {
126 
127  double Energy = exp((nHits + 26)/26);
128  if(Energy > 100) Energy = 100;
129 
130  return Energy;
131  }

◆ sortLogLik()

static bool JFIT::JORCAShowerFit::sortLogLik ( JShowerHypothesis  A,
JShowerHypothesis  B 
)
inlinestatic

Function to sort different shower hypotheses by their likelihood.

Definition at line 138 of file JORCAShowerFit.hh.

139  {
140  return(A.LogLik < B.LogLik);
141  };

◆ getJEvt()

void JFIT::JORCAShowerFit::getJEvt ( const KM3NETDAQ::JDAQTimeslice timeSliceBuildL0,
JFIT::JEvt InPreFits,
JFIT::JEvt OutFits 
) const

Declaration of the member function that actually performs the reconstruction.

member function definition

Parameters
timeSliceBuildL0which is contructed via a JDAQEvent
InPreFitsinput fits, normally produced by JORCAShowerPositionFit
OutFitsoutput fits

Definition at line 165 of file JORCAShowerFit.hh.

168 {
169 
170  const double Dmax_m = parameters_.Dmax_m;
171  const double R_Hz = parameters_.R_Hz;
172  const size_t numberOfPrefits = parameters_.numberOfPrefits;
173 
174  if ( InPreFits.empty() ) return;
175 
176  typedef vector<JTRIGGER::JHitL0> JDataL0_t;
177  typedef vector<JFIT::JPMTW0> JPMTW0_t;
178  JBuildL0<JTRIGGER::JHitL0> buildL0;
179 
180  const JDETECTOR::JDetector &detector = moduleRouter_->getReference();
181 
182  JEvt::iterator __end = InPreFits.end();
183  std::copy(InPreFits.begin(), __end, back_inserter(OutFits));
184 
185  if (numberOfPrefits > 0) {
186  advance(__end = InPreFits.begin(), min(numberOfPrefits, InPreFits.size()));
187  }
188 
189  partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
190 
191  JDataL0_t dataL0;
192  buildL0(timeSliceBuildL0, *moduleRouter_, back_inserter(dataL0));
193 
194  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
195  const JVector3D pos(getPosition(*shower));
196 
197  JDetectorSubset_t subdetector(detector, pos, Dmax_m);
198 
199  JPoint4D pt(pos, shower->getT());
200  JVector3D start_dir(0,0,0);
201 
202  // get the detector center and find the distance between the reco vertex and it.
203  const JCenter3D center(detector.begin(), detector.end());
204  const double radius = pos.getDistance(center);
205 
206  multiset<KM3NETDAQ::JDAQPMTIdentifier> top;
207 
208  double Nhits = 0;
209 
210  for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
211 
212  const JHitW0 hit(*i, R_Hz);
213  const JPosition3D hit_pos(hit.getPosition());
214 
215  const double D = hit_pos.getDistance(pt.getPosition());
216  const double t_res = hit.getT() - pt.getT(hit_pos);
217 
218  const JDirection3D photonDir(hit_pos - pt.getPosition());
219  const double cosT = photonDir.getDot(hit.getDirection());
220 
221  if(D < Dmax_m && (t_res >= -25 && t_res <= 25) && (cosT >= -1 && cosT <= 0.1)){
222  top.insert(hit.getPMTIdentifier());
223  const JVector3D d(photonDir.getDX(), photonDir.getDY(), photonDir.getDZ());
224  start_dir.add(d);
225  }
226  if(D < Dmax_m && (t_res > -30 && t_res < 20)){
227  Nhits++;
228  }
229  }
230 
231  double start_E = estimateInitialEnergy(Nhits);
232 
233  // rotate subdetector wrt start_dir, subdetector is already subtracted by pos
234 
235  const JDirection3D conversion(start_dir.getX(), start_dir.getY(), start_dir.getZ());
236  const JRotation3D R(conversion);
237 
238  vector<JPMTW0> buffer;
239 
240  JFIT::JORCAShowerFit::JShowerHypothesis shower_hypothesis;
241  vector<JFIT::JORCAShowerFit::JShowerHypothesis> shower_hypotheses;
242 
243  double chi2 = 0.0;
244  double max_theta = 0, max_phi = 0, scan_step = 0;
245  double Nmin = 0;
246 
247  if(start_E > 10 && radius < 100){
248  max_theta = 30;
249  max_phi = 30;
250  scan_step = 5;
251  Nmin = 10;
252  } else {
253  max_theta = 35;
254  max_phi = 35;
255  scan_step = 5;
256  Nmin = 20;
257  }
258 
259  const JAngle3D start_dir_angles(start_dir.getX(), start_dir.getY(), start_dir.getZ());
260 
261  for(double E = -15; E <= 5; E += 7.5){ // scan in energy
262  for(double th = -max_theta; th <= max_theta; th += scan_step){
263  for(double ph = -max_phi; ph <= max_phi; ph += scan_step){
264 
265  buffer.clear();
266  chi2 = 0.0;
267 
268  //todo: use precise value of pi
269  const double theta = th * PI / 180;
270  const double phi = ph * PI / 180;
271  const double scan_E = start_E + E; // [GeV]
272 
273  const JAngle3D dir_shifted_angles(start_dir_angles.getTheta() + theta,
274  start_dir_angles.getPhi() + phi);
275  const JVector3D dir_shifted(dir_shifted_angles.getDX(),
276  dir_shifted_angles.getDY(),
277  dir_shifted_angles.getDZ());
278  const JDirection3D conversion(dir_shifted.getX(),
279  dir_shifted.getY(),
280  dir_shifted.getZ());
281  const JRotation3D R(conversion);
282  pt.rotate(R);
283 
284  for (JDetectorSubset_t::iterator module = subdetector.begin();
285  module != subdetector.end(); ++module) {
286 
287  module->rotate(R);
288  JPosition3D module_pos(module->getPosition());
289 
290  for (unsigned int i = 0; i != module->size(); ++i) {
291  const JDAQPMTIdentifier id(module->getID(), i);
292  JPMT pmt(module->getPMT(i));
293  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
294  }
295  }
296 
297  for(JPMTW0_t::iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt ){
298  chi2 += fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), *pmt);
299  }
300 
301  shower_hypothesis.LogLik = chi2;
302  shower_hypothesis.Energy = scan_E;
303  shower_hypothesis.Direction = dir_shifted;
304 
305  shower_hypotheses.push_back(shower_hypothesis);
306 
307  pt.rotate_back(R);
308 
309  }
310  }
311  }
312 
313  vector<JShowerHypothesis>::iterator __end2 = shower_hypotheses.begin() + Nmin;
314  partial_sort(shower_hypotheses.begin(), __end2, shower_hypotheses.end(), sortLogLik);
315 
316  for(vector<JShowerHypothesis>::iterator hp = shower_hypotheses.begin(); hp != __end2; ++hp){
317 
318  fit_.step.resize(3);
319  const double dir_step = 0.05;
320  fit_.step[0] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(dir_step, 0.0), JEnergy());
321  fit_.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JEnergy());
322  fit_.step[2] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.05);
323 
324  const double scan_E = hp->Energy;
325 
326  buffer.clear();
327 
328  const JDirection3D startVersor(hp->Direction);
329 
330  const JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
331  JRotation3D R(conversion);
332  pt.rotate(R);
333 
334  for(JDetector::const_iterator module = detector.begin();
335  module != detector.end(); ++module) {
336  JPosition3D pos(module->getPosition());
337  pos.rotate(R);
338 
339  if (pt.getDistance(pos) < Dmax_m) {
340  for (unsigned int i = 0; i != module->size(); ++i) {
341  const JDAQPMTIdentifier id(module->getID(), i);
342  JPMT pmt(module->getPMT(i));
343  pmt.rotate(R);
344  buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
345  }
346  }
347  }
348 
349  const int NDF = buffer.size() - fit_.step.size();
350 
351  chi2 = fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), buffer.begin(), buffer.end());
352 
353  pt.rotate_back(R);
354 
355  const JVector3D resPos(fit_.value);
356  const JVersor3D resDir(fit_.value.getDX(), fit_.value.getDY(), fit_.value.getDZ());
357 
358  const JTrack3D td(resPos, resDir, shower->getT());
359  double E_reco_corrected = abs(fit_.value.getE() - 4);
360 
361  // Now that JShowerFit uses a Param PDF, the energy has to be corrected with a linear function obtained
362  // from a fit to the reconstruction error: this can be improved.
363 
364  const double p0 = -2.028;
365  const double p1 = 0.7683;
366  if(E_reco_corrected > 6){
367  E_reco_corrected = (E_reco_corrected - p0) / (1 + p1);
368  }
369 
370  JTrack3E tb(td, E_reco_corrected);
371  tb.rotate_back(R);
372 
373  const JFit &outFit = getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT),
374  tb, -chi2, NDF, E_reco_corrected,
376 
377  OutFits.push_back(outFit);
378  }
379  }
380  sort(OutFits.begin(), OutFits.end(), qualitySorter);
381 }

Member Data Documentation

◆ fit_

JRegressor_t JFIT::JORCAShowerFit::fit_
mutableprivate

Definition at line 69 of file JORCAShowerFit.hh.

◆ moduleRouter_

JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> JFIT::JORCAShowerFit::moduleRouter_
private

Definition at line 71 of file JORCAShowerFit.hh.

◆ parameters_

JFIT::JShowerFitParameters_t JFIT::JORCAShowerFit::parameters_
private

Definition at line 72 of file JORCAShowerFit.hh.


The documentation for this class was generated from the following file:
JFIT::JORCAShowerFit::JShowerHypothesis::LogLik
double LogLik
Definition: JORCAShowerFit.hh:118
JFIT::JORCAShowerFit::JShowerHypothesis::Direction
JVector3D Direction
Definition: JORCAShowerFit.hh:116
JFIT::JORCAShowerFit::sortLogLik
static bool sortLogLik(JShowerHypothesis A, JShowerHypothesis B)
Function to sort different shower hypotheses by their likelihood.
Definition: JORCAShowerFit.hh:138
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
JFIT::JORCAShowerFit::JShowerHypothesis
Container of shower hypothesis: During the fit the algorithm creates a grid in direction and energy s...
Definition: JORCAShowerFit.hh:114
JFIT::JORCAShowerFit::estimateInitialEnergy
static double estimateInitialEnergy(int nHits)
Function to estimate the initial em shower energy.
Definition: JORCAShowerFit.hh:125
JFIT::JShowerFitParameters_t::numberOfPrefits
std::size_t numberOfPrefits
Definition: JShowerFitParameters_t.hh:25
JFIT::JORCAShowerFit::fit_
JRegressor_t fit_
Definition: JORCAShowerFit.hh:69
JFIT::JShowerFitParameters_t::Dmax_m
double Dmax_m
Definition: JShowerFitParameters_t.hh:23
JAANET::copy
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:152
JFIT::qualitySorter
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
Definition: JEvtToolkit.hh:223
debug
int debug
debug level
Definition: JSirene.cc:59
JLANG::JReference::getReference
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
JFIT::JShowerFitParameters_t::R_Hz
double R_Hz
Definition: JShowerFitParameters_t.hh:24
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JFIT::JORCAShowerFit::JShowerHypothesis::Energy
double Energy
Definition: JORCAShowerFit.hh:117
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JFIT::JSHOWERCOMPLETEFIT
JShowerFit.cc.
Definition: JFitApplications.hh:34
JFIT::JShowerFitParameters_t::Tmax_ns
double Tmax_ns
Definition: JShowerFitParameters_t.hh:22
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
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
JFIT::JORCAShowerFit::moduleRouter_
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > moduleRouter_
Definition: JORCAShowerFit.hh:71
JFIT::getPosition
JPosition3D getPosition(const JFit &fit)
Get position.
Definition: JEvtToolkit.hh:63
JFIT::JORCAShowerFit::parameters_
JFIT::JShowerFitParameters_t parameters_
Definition: JORCAShowerFit.hh:72
JFIT::OKAY
Definition: JFitStatus.hh:22