RefineSeeds.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/IFitVar.h"
5 
6 #include "OscLib/IOscCalc.h"
7 
8 #include "Utilities/func/MathUtil.h"
9 
10 namespace ana
11 {
12  // --------------------------------------------------------------------------
14  const IExperiment* expt,
15  const std::vector<const IFitVar*>& vars,
16  const osc::IOscCalcAdjustable* calc0,
17  double dchisq_max)
18  {
19  struct Result
20  {
21  Result(double c, const Seed& p) : chisq(c), pt(p) {}
22 
23  bool operator<(const Result& r) const {return chisq < r.chisq;}
24 
25  double Distance(const Seed& s) const
26  {
27  auto a = pt.GetVals();
28  auto b = s.GetVals();
29 
30  // Compute Euclidean distance naively. Approximately OK since we
31  // normally scale our parameters to have human-ish scale.
32  double dist2 = 0;
33  for(auto it: a) dist2 += util::sqr(it.second - b[it.first]);
34  return sqrt(dist2);
35  }
36 
37  double chisq;
38  Seed pt;
39  };
40 
41  std::vector<Result> results;
42 
43  // Explore all the seeds to find out where a no-systs minimization of them
44  // ends up
45  for(const Seed& seed: seeds.GetSeeds()){
46  osc::IOscCalcAdjustable* calc = calc0->Copy();
47  seed.ResetCalc(calc);
48 
49  MinuitFitter fit(expt, vars);
50  const double chi = fit.Fit(calc, IFitter::kQuiet)->EvalMetricVal();
51 
52  std::map<const IFitVar*, double> valmap;
53  for(const IFitVar* v: vars) valmap[v] = v->GetValue(calc);
54  results.emplace_back(chi, valmap);
55  }
56 
57  // Lowest chisq first
58  std::sort(results.begin(), results.end());
59 
60  const double bestChisq = results[0].chisq;
61 
62  // Cull the list to a list of unique results not too far above the best
63  // result
64  std::vector<Seed> ret;
65  for(const Result& result: results){
66  // Stop once seed results become too bad
67  if(dchisq_max > 0 && result.chisq > bestChisq + dchisq_max) break;
68 
69  bool ok = true;
70  for(const Seed& already: ret){
71  // Reject any points that are very nearby other previously added
72  // seeds. The threshold distance is arbitrary. We're really trying to
73  // determine if these are the same minimum.
74  if(result.Distance(already) < 1e-3){
75  ok = false;
76  break;
77  }
78  } // end for already
79  if(ok) ret.push_back(result.pt);
80  }
81 
82  return ret;
83  }
84 }
virtual _IOscCalcAdjustable< T > * Copy() const =0
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
const std::vector< Seed > & GetSeeds() const
Definition: SeedList.h:47
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
double Distance(double x1, double y1, double x2, double y2)
const std::map< const IFitVar *, double > GetVals() const
Definition: SeedList.h:25
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
const XML_Char * s
Definition: expat.h:262
expt
Definition: demo5.py:34
unsigned int seed
Definition: runWimpSim.h:102
const double a
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:64
GenericCut< T > operator<(double c, const GenericVar< T > &v)
Definition: Var.h:153
const std::map< std::pair< std::string, std::string >, Variable > vars
Base class defining interface for experiments.
Definition: IExperiment.h:14
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
Interface definition for fittable variables.
Definition: IFitVar.h:16
Float_t e
Definition: plot.C:35
SeedList RefineSeeds(const SeedList &seeds, const IExperiment *expt, const std::vector< const IFitVar * > &vars, const osc::IOscCalcAdjustable *calc0, double dchisq_max)
Refine an initial list of seeds by exploring stats-only fits.
Definition: RefineSeeds.cxx:13
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17