1 #ifndef __JFIT__JGRADIENT__
2 #define __JFIT__JGRADIENT__
21 namespace JPP {
using namespace JFIT; }
40 virtual void apply(
const double step) = 0;
48 public std::shared_ptr<JParameter_t>
91 const int debug = 3) :
123 return numeric_limits<double>::max();
128 const size_t N = this->size();
133 for (
size_t i = 0;
i !=
N; ++
i) {
152 for (
double ds = 1.0; ds > 1.0e-3; ) {
154 this->
move(+1.0 * ds);
158 DEBUG(
"chi2[3] " << setw(4) << m <<
' ' <<
FIXED(12,5) << chi2[3] <<
' ' <<
FIXED(12,5) << ds << endl);
160 if (chi2[3] < chi2[2]) {
184 for ( ; m != 0; --m) {
185 this->
move(-1.0 * ds);
192 this->
move(-1.0 * ds);
194 if (chi2[2] < chi2[3]) {
202 const double f21 = chi2[2] - chi2[1];
203 const double f23 = chi2[2] - chi2[3];
205 const double xs = 0.5 * (f21 - f23) / (f23 + f21);
207 this->
move(+1.0 * xs * ds);
211 if (chi2[3] < chi2[2]) {
217 this->
move(-1.0 * xs * ds);
232 if (fabs(chi2[2] - chi2[0]) <
epsilon * 0.5 * (fabs(chi2[0]) + fabs(chi2[2]))) {
244 for (
size_t i = 0;
i !=
N; ++
i){
255 for (
size_t i = 0;
i !=
N; ++
i){
257 H[
i] = G[
i] + dgg * H[
i];
287 const size_t N = this->size();
299 for (
size_t i = 0;
i !=
N; ++
i) {
300 if ((*
this)[
i].
name.size() > width) {
301 width = (*this)[
i].name.size();
307 for (
size_t i = 0;
i !=
N; ++
i) {
309 if ((*
this)[
i].value != 0.0) {
311 (*this)[
i]->apply(+0.5 * (*
this)[
i].value);
315 (*this)[
i]->apply(-0.5 * (*
this)[
i].value);
316 (*this)[
i]->apply(-0.5 * (*
this)[
i].value);
322 (*this)[
i]->apply(+0.5 * (*
this)[
i].value);
334 DEBUG(setw(width) << left <<
"|gradient|" << right <<
' ' <<
FIXED(12,5) << sqrt(
V) << endl);
348 for (
size_t i = 0;
i != this->size(); ++
i) {
349 (*this)[
i ]->apply((*
this)[
i ].value *
gradient[
i ] * factor);
351 }
else if (factor < 0.0) {
352 for (
size_t i = this->size();
i != 0; --
i) {
353 (*this)[
i-1]->apply((*
this)[
i-1].value *
gradient[
i-1] * factor);
JModifier_t(const std::string &name, JParameter_t *parameter, const double value)
Constructor.
void move(const double factor)
Move.
virtual ~JParameter_t()
Virtual destructor.
double operator()(const T &getChi2)
Fit.
static const double H
Planck constant [eV s].
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
size_t Nextra
maximum number of extra steps
do set_variable OUTPUT_DIRECTORY $WORKDIR T
JGradient(const size_t Nmax=std::numeric_limits< size_t >::max(), const size_t Nextra=0, const double epsilon=1.0e-4, const int debug=3)
Constructor.
std::vector< double > gradient
virtual void apply(const double step)=0
Apply step.
size_t Nmax
maximum number of iterations
Auxiliary data structure for fit parameter.
General purpose messaging.
then fatal The output file must have the wildcard in the name
Auxiliary data structure for editable parameter.
double evaluate(const T &getChi2)
Evaluate gradient.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
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
size_t numberOfIterations
double getChi2(const double P)
Get chi2 corresponding to given probability.
Auxiliary data structure for floating point format specification.
#define DEBUG(A)
Message macros.