Fit function.
162 for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
172 for (JDataL0_t::const_iterator
i = dataL0.begin();
i !=dataL0.end(); ++
i) {
179 top.insert(
i->getPMTIdentifier());
191 for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
197 for (
size_t i = 0;
i != module->size(); ++
i) {
206 const size_t count = top.count(
id);
208 data.push_back(
JNPEHit(this->
getNPE(module->getPMT(
i), rate_Hz), count));
214 const int NDF =
distance(data.begin(), data.end()) - 1;
224 for (
int i = 0;
i !=
N; ++
i) {
234 for (
int i = 0;
i !=
N; ++
i) {
239 const double chi2 = (*this)(
x, data.begin(), data.end());
241 result[
i].chi2 =
chi2;
242 buffer[
chi2] = x.getE();
245 if (result[i].chi2 < result[j].chi2) {
251 for (
int i = 0; i !=
N; ++
i) {
252 DEBUG(
' ' <<
FIXED(5,2) << result[i].x <<
' ' <<
FIXED(9,3) << result[i].chi2);
261 result[4] = result[1];
262 result[2] = result[0];
263 result[0] = JResult(2*result[2].x - result[4].x);
267 result[4] = result[2];
268 result[2] = result[1];
272 result[0] = result[1];
273 result[4] = result[3];
277 result[0] = result[2];
278 result[2] = result[3];
282 result[0] = result[3];
283 result[2] = result[4];
284 result[4] = JResult(2*result[2].x - result[0].x);
288 result[1] = JResult(0.5 * (result[0].x + result[2].x));
289 result[3] = JResult(0.5 * (result[2].x + result[4].x));
291 }
while (result[4].x - result[0].x >
resolution);
294 if (result[1].chi2 != result[3].chi2) {
296 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);
297 result[2].chi2 = (*this)(result[2].x, data.begin(), data.end());
301 const double chi2 = result[2].chi2;
302 const double E = result[2].x.getE();
306 double Emin = numeric_limits<double>::max();
307 double Emax = numeric_limits<double>::lowest();
310 if (i->second < Emin) { Emin = i->second; }
311 if (i->second > Emax) { Emax = i->second; }
314 const double mu_range =
gWater(E);
316 double noise_likelihood = 0.0;
317 int number_of_hits = 0;
320 noise_likelihood +=
log10(
getP(i->getY0(), i->getN()));
321 number_of_hits += i->getN();
326 out.rbegin()->setE(E);
328 out.rbegin()->setW(track->getW());
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
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.
static const int JENERGY_MAXIMAL_ENERGY
maximal energy [GeV] from JEnergy.cc
double roadWidth_m
road width [m]
const JSummaryRouter & summary
double ZMin_m
minimal z-position [m]
double getRate() const
Get default rate.
Template specialisation of class JModel to match hit with muon trajectory along z-axis.
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Auxiliary data structure for floating point format specification.
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
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.1.0 https://git.km3net.de/common/km3net-dataformat.
JDirection3D getDirection(const JFit &fit)
Get direction.
Auxiliary class for simultaneously handling light yields and response of PMT.
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...
then JCookie sh JDataQuality D $DETECTOR_ID R
bool getPMTStatus(const JStatus &status)
Test status of PMT.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
const JClass_t & getReference() const
Get reference to object.
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
Data structure for fit of straight line paralel to z-axis.
Reduced data structure for L1 hit.
Data structure for fit of energy.
JPosition3D getPosition(const JFit &fit)
Get position.
double getNPE(const Hit &hit)
Get true charge of hit.
bool getDAQStatus(const JDAQFrameStatus &frame, const JStatus &status)
Test status of DAQ.
JPosition3D & rotate(const JRotation3D &R)
Rotate.
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
static const int JMUONENERGY
double EMax_log
maximal energy [log10(GeV)]
#define DEBUG(A)
Message macros.
1.8.5