Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
JRECONSTRUCTION::JMuonEnergy Class Reference

Auxiliary class to to determine muon energy. More...

#include <JMuonEnergy.hh>

Inheritance diagram for JRECONSTRUCTION::JMuonEnergy:
JRECONSTRUCTION::JMuonEnergyParameters_t JFIT::JRegressor< JModel_t, JMinimiser_t > TObject

Classes

struct  JResult
 Auxiliary class for energy estimation. More...
 

Public Types

typedef JRegressor< JEnergyJRegressor_t
 
typedef JTRIGGER::JHitR1 hit_type
 
typedef std::vector< hit_typebuffer_type
 
typedef JTOOLS::JRange< JEnergyJEnergyRange
 

Public Member Functions

 JMuonEnergy (const JMuonEnergyParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string &pdfFile, const int debug=0)
 Constructor. More...
 
JEvt operator() (const KM3NETDAQ::JDAQEvent &event, const JEvt &in)
 Fit function. More...
 
void reset ()
 Reset fit parameters. More...
 
bool equals (const JMuonEnergyParameters_t &parameters) const
 Equality. More...
 
 ClassDef (JMuonEnergyParameters_t, 1)
 

Public Attributes

double roadWidth_m
 road width [m] More...
 
double R_Hz
 default rate [Hz] More...
 
size_t numberOfPrefits
 number of prefits More...
 
double EMin_log
 minimal energy [log10(GeV)] More...
 
double EMax_log
 maximal energy [log10(GeV)] More...
 
double TMin_ns
 minimal time w.r.t. Cherenkov hypothesis [ns] More...
 
double TMax_ns
 maximal time w.r.t. Cherenkov hypothesis [ns] More...
 
double ZMin_m
 minimal z-position [m] More...
 
double resolution
 energy resolution [log10(GeV)] More...
 
int mestimator
 M-estimator. More...
 
bool reprocess
 reprocess More...
 

Private Attributes

const JModuleRouterrouter
 
const JSummaryRoutersummary
 
int debug
 

Detailed Description

Auxiliary class to to determine muon energy.

Definition at line 58 of file JMuonEnergy.hh.

Member Typedef Documentation

Definition at line 63 of file JMuonEnergy.hh.

Definition at line 64 of file JMuonEnergy.hh.

Definition at line 65 of file JMuonEnergy.hh.

Definition at line 66 of file JMuonEnergy.hh.

Constructor & Destructor Documentation

JRECONSTRUCTION::JMuonEnergy::JMuonEnergy ( const JMuonEnergyParameters_t parameters,
const JModuleRouter router,
const JSummaryRouter summary,
const std::string &  pdfFile,
const int  debug = 0 
)
inline

Constructor.

Parameters
parametersparameters
routermodule router
summarysummary router
pdfFilePDF file
debugdebug

Definition at line 79 of file JMuonEnergy.hh.

83  :
84  JMuonEnergyParameters_t(parameters),
85  JRegressor_t(pdfFile),
86  router(router),
87  summary(summary),
88  debug(debug)
89  {
90  using namespace JPP;
91 
94 
95  if (this->getRmax() < roadWidth_m) {
96  roadWidth_m = this->getRmax();
97  }
98 
99  this->estimator.reset(getMEstimator(mestimator));
100 
101  if (!this->estimator.is_valid()) {
102  FATAL("Invalid M-Estimator." << endl);
103  }
104  }
static int debug
debug level (default is off).
Definition: JMessage.hh:45
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
const JSummaryRouter & summary
Definition: JMuonEnergy.hh:346
JRegressor< JEnergy > JRegressor_t
Definition: JMuonEnergy.hh:63
const JModuleRouter & router
Definition: JMuonEnergy.hh:345
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
#define FATAL(A)
Definition: JMessage.hh:67
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
JMEstimator * getMEstimator(const int type)
Get M-Estimator.
Definition: JMEstimator.hh:203

Member Function Documentation

JEvt JRECONSTRUCTION::JMuonEnergy::operator() ( const KM3NETDAQ::JDAQEvent event,
const JEvt in 
)
inline

Fit function.

Parameters
eventevent
instart values
Returns
fit results

Definition at line 146 of file JMuonEnergy.hh.

147  {
148  using namespace std;
149  using namespace JPP;
150 
151  JEvt out;
152 
153  typedef vector<JHitL0> JDataL0_t;
154  const JBuildL0<JHitL0> buildL0;
155 
157 
158  JDataL0_t dataL0;
159 
160  buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
161 
162  for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
163 
164  const JRotation3D R (getDirection(*track));
165  const JLine1Z tz(getPosition (*track).rotate(R), track->getT());
167 
169 
171 
172  for (JDataL0_t::const_iterator i = dataL0.begin(); i !=dataL0.end(); ++i) {
173 
174  JHitR1 hit(*i);
175 
176  hit.rotate(R);
177 
178  if (match(hit)) {
179  top.insert(i->getPMTIdentifier());
180  }
181  }
182 
183  JRange<double> Z_m;
184 
185  if ( track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0 ) {
186  Z_m.setLowerLimit(this->ZMin_m);
187  }
188 
189  const JDetectorSubset_t subdetector(detector, getAxis(*track), roadWidth_m, Z_m);
190 
191  for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
192 
193  if (summary.hasSummaryFrame(module->getID())) {
194 
195  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
196 
197  for (size_t i = 0; i != module->size(); ++i) {
198 
199  if (getDAQStatus(frame, *module, i) &&
200  getPMTStatus(frame, *module, i) &&
201  frame[i].is_valid() &&
202  !module->getPMT(i).has(PMT_DISABLE)) {
203 
204  const JDAQPMTIdentifier id(module->getID(), i);
205 
206  const double rate_Hz = summary.getRate(id);
207  const size_t count = top.count(id);
208 
209  data.push_back(JNPEHit(this->getNPE(module->getPMT(i), rate_Hz), count));
210  }
211  }
212  }
213  }
214 
215  const int NDF = distance(data.begin(), data.end()) - 1;
216 
217  if( NDF >= 0 ) {
218 
219  // 5-point search between given limits
220 
221  const int N = 5;
222 
223  JResult result[N];
224 
225  for (int i = 0; i != N; ++i) {
226  result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
227  }
228 
229  map<double, double> buffer;
230 
231  do {
232 
233  int j = 0;
234 
235  for (int i = 0; i != N; ++i) {
236 
237  if (!result[i]) {
238 
239  const JEnergy x = result[i].x;
240  const double chi2 = (*this)(x, data.begin(), data.end());
241 
242  result[i].chi2 = chi2;
243  buffer[chi2] = x.getE();
244  }
245 
246  if (result[i].chi2 < result[j].chi2) {
247  j = i;
248  }
249  }
250 
251 
252  for (int i = 0; i != N; ++i) {
253  DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
254  }
255  DEBUG(endl);
256 
257  // squeeze range
258 
259  switch (j) {
260 
261  case 0:
262  result[4] = result[1];
263  result[2] = result[0];
264  result[0] = JResult(2*result[2].x - result[4].x);
265  break;
266 
267  case 1:
268  result[4] = result[2];
269  result[2] = result[1];
270  break;
271 
272  case 2:
273  result[0] = result[1];
274  result[4] = result[3];
275  break;
276 
277  case 3:
278  result[0] = result[2];
279  result[2] = result[3];
280  break;
281 
282  case 4:
283  result[0] = result[3];
284  result[2] = result[4];
285  result[4] = JResult(2*result[2].x - result[0].x);
286  break;
287  }
288 
289  result[1] = JResult(0.5 * (result[0].x + result[2].x));
290  result[3] = JResult(0.5 * (result[2].x + result[4].x));
291 
292  } while (result[4].x - result[0].x > resolution);
293 
294 
295  if (result[1].chi2 != result[3].chi2) {
296 
297  result[2].x += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);
298  result[2].chi2 = (*this)(result[2].x, data.begin(), data.end());
299 
300  }
301 
302  const double chi2 = result[2].chi2;
303  const double E = result[2].x.getE();
304 
305  // calculate additional variables
306 
307  double Emin = numeric_limits<double>::max();
308  double Emax = numeric_limits<double>::lowest();
309 
310  for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
311  if (i->second < Emin) { Emin = i->second; }
312  if (i->second > Emax) { Emax = i->second; }
313  }
314 
315  const double mu_range = gWater(E); // range of a muon with the reconstructed energy
316 
317  double noise_likelihood = 0.0; // log-likelihood of every hit being from K40
318  int number_of_hits = 0; // number of hits selected for JEnergy
319 
320  for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
321  noise_likelihood += log10(getP(i->getY0(), i->getN())); // probability of getting the observed multiplicity with K40
322  number_of_hits += i->getN(); // hit multiplicity
323  }
324 
325  out.push_back(JFit(*track).add(JMUONENERGY));
326 
327  out.rbegin()->setE(E);
328 
329  out.rbegin()->setW(track->getW());
330  out.rbegin()->setW(JENERGY_ENERGY, E);
331  out.rbegin()->setW(JENERGY_CHI2, chi2);
332  out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
333  out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD, noise_likelihood);
334  out.rbegin()->setW(JENERGY_NDF, NDF);
335  out.rbegin()->setW(JENERGY_NUMBER_OF_HITS, number_of_hits);
336  out.rbegin()->setW(JENERGY_MINIMAL_ENERGY, Emin);
337  out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY, Emax);
338  }
339  }
340 
341  return out;
342  }
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
static const int JENERGY_ENERGY
uncorrected energy [GeV] from JEnergy.cc
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Template specialisation of L0 builder for JHitL0 data type.
Definition: JBuildL0.hh:102
static const int JENERGY_MAXIMAL_ENERGY
maximal energy [GeV] from JEnergy.cc
Detector data structure.
Definition: JDetector.hh:89
const JSummaryRouter & summary
Definition: JMuonEnergy.hh:346
Rotation matrix.
Definition: JRotation3D.hh:111
double getRate() const
Get default rate.
then usage $script< input file >[option] nPossible options count
Definition: JVolume1D.sh:31
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
Definition: JFit/JModel.hh:34
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
const JDAQSummaryFrame & getSummaryFrame() const
Get default summary frame.
JFit & add(const int type)
Add event to history.
static const int JENERGY_CHI2
chi2 from JEnergy.cc
double resolution
energy resolution [log10(GeV)]
static const int JENERGY_NOISE_LIKELIHOOD
log likelihood of every hit being K40 from JEnergy.cc
static const int JENERGY_NDF
number of degrees of freedom from JEnergy.cc
static const int JENERGY_MINIMAL_ENERGY
minimal energy [GeV] from JEnergy.cc
const JModuleRouter & router
Definition: JMuonEnergy.hh:345
Detector file.
Definition: JHead.hh:226
Acoustic event fit.
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
set_variable E_E log10(E_{fit}/E_{#mu})"
Data storage class for rate measurements of all PMTs in one module.
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
static const int PMT_DISABLE
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
Definition: pmt_status.hh:12
JDirection3D getDirection(const JFit &fit)
Get direction.
Data time slice.
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
Detector subset without binary search functionality.
double getP(const double expval, bool hit)
Get Poisson probability to observe a hit or not for given expectation value for the number of hits...
Definition: JFitToolkit.hh:41
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
bool getPMTStatus(const JStatus &status)
Test status of PMT.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
const JClass_t & getReference() const
Get reference to object.
Definition: JReference.hh:38
static const int JSTART_LENGTH_METRES
distance between first and last hits in metres from JStart.cc
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
double EMin_log
minimal energy [log10(GeV)]
then if[[!-f $DETECTOR]] then JDetector sh $DETECTOR fi cat $WORKDIR trigger_parameters txt<< EOFtrigger3DMuon.enabled=1;trigger3DMuon.numberOfHits=5;trigger3DMuon.gridAngle_deg=1;ctMin=0.0;TMaxLocal_ns=15.0;EOF set_variable TRIGGEREFFICIENCY_TRIGGERED_EVENTS_ONLY INPUT_FILES=() for((i=1;$i<=$NUMBER_OF_RUNS;++i));do JSirene.sh $DETECTOR $JPP_DATA/genhen.km3net_wpd_V2_0.evt.gz $WORKDIR/sirene_ ${i}.root JTriggerEfficiency.sh $DETECTOR $DETECTOR $WORKDIR/sirene_ ${i}.root $WORKDIR/trigger_efficiency_ ${i}.root $WORKDIR/trigger_parameters.txt $JPP_DATA/PMT_parameters.txt INPUT_FILES+=($WORKDIR/trigger_efficiency_ ${i}.root) done for ANGLE_DEG in $ANGLES_DEG[*];do set_variable SIGMA_NS 3.0 set_variable OUTLIERS 3 set_variable OUTPUT_FILE $WORKDIR/matrix\[${ANGLE_DEG}\deg\].root $JPP_DIR/examples/JReconstruction-f"$INPUT_FILES[*]"-o $OUTPUT_FILE-S ${SIGMA_NS}-A ${ANGLE_DEG}-O ${OUTLIERS}-d ${DEBUG}--!fiif[[$OPTION=="plot"]];then if((0));then for H1 in h0 h1;do JPlot1D-f"$WORKDIR/matrix["${^ANGLES_DEG}" deg].root:${H1}"-y"1 2e3"-Y-L TR-T""-\^"number of events [a.u.]"-> o chi2
Definition: JMatrixNZ.sh:106
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
Reduced data structure for L1 hit.
Definition: JHitR1.hh:31
Data structure for fit of energy.
Definition: JEnergy.hh:28
JPosition3D getPosition(const JFit &fit)
Get position.
double getNPE(const Hit &hit)
Get true charge of hit.
int j
Definition: JPolint.hh:792
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
Definition: JPosition3D.hh:186
JAxis3D getAxis(const JFit &fit)
Get axis.
static const int JENERGY_MUON_RANGE_METRES
range of a muon with the reconstructed energy [m] from JEnergy.cc
void setLowerLimit(argument_type x)
Set lower limit.
Definition: JRange.hh:224
static const int JMUONENERGY
double EMax_log
maximal energy [log10(GeV)]
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
void JRECONSTRUCTION::JMuonEnergyParameters_t::reset ( )
inlineinherited

Reset fit parameters.

Definition at line 44 of file JMuonEnergyParameters_t.hh.

45  {
46  roadWidth_m = std::numeric_limits<double>::max();
47  R_Hz = 6.0e3;
48  numberOfPrefits = 0;
49  EMin_log = 1.0;
50  EMax_log = 8.0;
51  TMin_ns = -50.0;
52  TMax_ns = +450.0;
53  ZMin_m = std::numeric_limits<double>::lowest();
54  resolution = 0.01;
56  reprocess = false;
57  }
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double resolution
energy resolution [log10(GeV)]
double EMin_log
minimal energy [log10(GeV)]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double EMax_log
maximal energy [log10(GeV)]
bool JRECONSTRUCTION::JMuonEnergyParameters_t::equals ( const JMuonEnergyParameters_t parameters) const
inlineinherited

Equality.

Parameters
parametersfit parameters
Returns
true if equals; else false

Definition at line 65 of file JMuonEnergyParameters_t.hh.

66  {
67  return (this->roadWidth_m == parameters.roadWidth_m &&
68  this->R_Hz == parameters.R_Hz &&
69  this->numberOfPrefits == parameters.numberOfPrefits &&
70  this->EMin_log == parameters.EMin_log &&
71  this->EMax_log == parameters.EMax_log &&
72  this->TMin_ns == parameters.TMin_ns &&
73  this->TMax_ns == parameters.TMax_ns &&
74  this->ZMin_m == parameters.ZMin_m &&
75  this->resolution == parameters.resolution &&
76  this->mestimator == parameters.mestimator &&
77  this->reprocess == parameters.reprocess);
78  }
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
then usage $script< detector file >< detectorfile > nIf the range of floors is the first detector file is aligned to the second before the comparison nIn this
double resolution
energy resolution [log10(GeV)]
double EMin_log
minimal energy [log10(GeV)]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double EMax_log
maximal energy [log10(GeV)]
JRECONSTRUCTION::JMuonEnergyParameters_t::ClassDef ( JMuonEnergyParameters_t  ,
 
)
inherited

Member Data Documentation

const JModuleRouter& JRECONSTRUCTION::JMuonEnergy::router
private

Definition at line 345 of file JMuonEnergy.hh.

const JSummaryRouter& JRECONSTRUCTION::JMuonEnergy::summary
private

Definition at line 346 of file JMuonEnergy.hh.

int JRECONSTRUCTION::JMuonEnergy::debug
private

Definition at line 347 of file JMuonEnergy.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::roadWidth_m
inherited

road width [m]

Definition at line 82 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::R_Hz
inherited

default rate [Hz]

Definition at line 83 of file JMuonEnergyParameters_t.hh.

size_t JRECONSTRUCTION::JMuonEnergyParameters_t::numberOfPrefits
inherited

number of prefits

Definition at line 84 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::EMin_log
inherited

minimal energy [log10(GeV)]

Definition at line 85 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::EMax_log
inherited

maximal energy [log10(GeV)]

Definition at line 86 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::TMin_ns
inherited

minimal time w.r.t. Cherenkov hypothesis [ns]

Definition at line 87 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::TMax_ns
inherited

maximal time w.r.t. Cherenkov hypothesis [ns]

Definition at line 88 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::ZMin_m
inherited

minimal z-position [m]

Definition at line 89 of file JMuonEnergyParameters_t.hh.

double JRECONSTRUCTION::JMuonEnergyParameters_t::resolution
inherited

energy resolution [log10(GeV)]

Definition at line 90 of file JMuonEnergyParameters_t.hh.

int JRECONSTRUCTION::JMuonEnergyParameters_t::mestimator
inherited

M-estimator.

Definition at line 91 of file JMuonEnergyParameters_t.hh.

bool JRECONSTRUCTION::JMuonEnergyParameters_t::reprocess
inherited

reprocess

Definition at line 92 of file JMuonEnergyParameters_t.hh.


The documentation for this class was generated from the following file: