FitUtilities.h
Go to the documentation of this file.
1 /// FitUtilities.h
3 
4 enum SystType {
10 };
11 
12 enum SystConfig {
17 };
18 
25 };
26 
27 enum FitType {
31 };
32 
33 // Struct to hold fit configuration options
34 struct FitConfig {
35 
36  unsigned int jobNumber; // Job number
37  unsigned int nJobs; // Number of jobs to parallelise over
38  SystType systType; // Type of systematic fit to run
39  SystConfig systConfig; // What configuration of systs to run with
40  std::string systName; // Name of single systematic, if applicable
41 
42  double epsilon; // epsilon value to add to matrix diagonal
43 
44  std::vector<covmx::Sample> samples; // Samples to fit with
45 
46  FakeDataType fakeDataType; // Fake data type
47  unsigned int fakeDataSet; // Which data set (ie. for fake signal, we may have different definitions of signal)
48  unsigned int fakeDataUniverse; // For datasets with randomly fluctuated universes, which universe to use
49 
50  const IFitVar* xVar; // Variable to fit in X
51  Binning* xBins; // Binning in X
52  const IFitVar* yVar; // Variable to fit in Y
53  Binning* yBins; // Binning in Y
54 
55  std::vector<const IFitVar*> fitVars; // Other variables to fit over
56  std::map<const IFitVar*, std::vector<double>> fitSeeds; // Seed values to use for fit
57 
58  FitType fitType; // Which fitting method to use
59 };
60 
61 /// Parse an option string into a fit configuration
63 
64  std::string optstring(opt.Data()); // transfer TString to std::string so we can use string::find
65 
66  FitConfig cfg; // Empty config to populate and return
67 
68  // Parallelisation options
69  cfg.jobNumber = 999;
70  size_t pos = optstring.find("jobnumber=");
71  if (pos != std::string::npos) {
72  size_t start = pos + 10;
73  size_t end = optstring.find("_", pos);
74  cfg.jobNumber = std::stoi(optstring.substr(start, end-start));
75  } else if (RunningOnGrid()) {
76  cfg.jobNumber = JobNumber();
77  } else if (optstring.find("njobs=")) {
78  cfg.jobNumber = 0; // for testing purposes
79  }
80 
81  cfg.nJobs = 999;
82  pos = optstring.find("njobs=");
83  if (pos != std::string::npos) {
84  size_t start = pos + 6;
85  size_t end = optstring.find("_", pos);
86  cfg.nJobs = std::stoi(optstring.substr(start, end-start));
87  std::cout << "Running surface job " << cfg.jobNumber << " / " << cfg.nJobs << std::endl;
88  }
89 
90  // Get pull and covariance matrix syst options
91  if (CheckOption(opt, "statsonly")) cfg.systType = kStatsOnly;
92  else if (CheckOption(opt, "covmxsysts")) cfg.systType = kCovMxSysts;
93  else if (optstring.find("covmxsyst=") != std::string::npos) {
94  cfg.systType = kCovMxSingleSyst;
95  pos = optstring.find("covmxsyst=");
96  size_t start = pos + 10;
97  size_t end = optstring.find("_", pos);
98  cfg.systName = optstring.substr(start, end-start);
99  std::cout << "Single CovMx Systematic Fit: " << cfg.systName << std::endl;
100  }
101  else if (optstring.find("pullsyst=") != std::string::npos) {
102  pos = optstring.find("pullsyst=");
103  size_t start = pos + 9;
104  size_t end = optstring.find("_", pos);
105  cfg.systName = optstring.substr(start, end-start);
106  if (CheckOption(opt, "covmxsysts")) {
107  cfg.systType = kCovMxWithSinglePullSyst;
108  std::cout << "Fitting with covariance matrix but with "
109  << cfg.systName << " as a pull term." << std::endl;
110  } else {
111  cfg.systType = kPullSingleSyst;
112  std::cout << "Fitting with " << cfg.systName << " as a pull term." << std::endl;
113  }
114  }
115 
116  // Systematic configuration
117  pos = optstring.find("systconfig=");
118  if (pos != std::string::npos) {
119  size_t start = pos + 11;
120  size_t end = optstring.find("_", start);
121  std::string config = optstring.substr(start, end-start);
122  if (config == "detector") cfg.systConfig = kDetectorSysts;
123  else if (config == "beam") cfg.systConfig = kBeamSysts;
124  else if (config == "xsec") cfg.systConfig = kXsecSysts;
125  else if (config == "all") cfg.systConfig = kAllSysts;
126  else {
127  std::ostringstream err;
128  err << "Systematic config \"" << config << "\" not recognised!";
129  throw std::runtime_error(err.str());
130  }
131  } else {
132  cfg.systConfig = kAllSysts; // Default to all systematics
133  }
134 
135  // Get sample configuration
136  cfg.samples = GetSamplesFromOptString(opt);
137 
138  // Get fake data configuration
139  pos = optstring.find("fakedata=");
140  if (pos != std::string::npos) {
141  size_t start = pos + 9;
142  if (optstring.substr(start, 6) == "asimov") {
143  cfg.fakeDataType = kAsimov;
144  } else if (optstring.substr(start, 5) == "stats") {
145  cfg.fakeDataType = kStatFluctuated;
146  size_t end = optstring.find("_", pos);
147  cfg.fakeDataUniverse = std::stoul(optstring.substr(start+5, end-(start+5)));
148  } else if (optstring.substr(start, 5) == "systs") {
149  cfg.fakeDataType = kStatAndSystFluctuated;
150  size_t end = optstring.find("_", pos);
151  cfg.fakeDataUniverse = std::stoul(optstring.substr(start+5, end-(start+5)));
152  } else if (optstring.substr(start, 10) == "fakesignal") {
153  std::string err = "Fake signal data must take form \"fakesignal-i-j\", where i is the fake signal set (1-3) and j is the universe (1-1000)!";
154  size_t end = optstring.find("-", start);
155  if (end-start != 10) throw std::runtime_error(err);
156  start += 11;
157  // Look for another dash, which is an indicator we want to load a pseudoexperiment
158  size_t nextDash = optstring.find("-", start);
159  size_t nextUnderscore = optstring.find("_", start);
160  if (nextDash < nextUnderscore) { // Loading a pseudoexperiment
161  cfg.fakeDataType = kFakeSignalFluctuated;
162  end = nextDash;
163  cfg.fakeDataSet = std::stoul(optstring.substr(start, end-start));
164  start = end+1;
165  end = optstring.find("_", start);
166  cfg.fakeDataUniverse = std::stoul(optstring.substr(start, end-start));
167  } else { // NOT loading a pseudoexperiment
168  cfg.fakeDataType = kFakeSignal;
169  end = nextUnderscore;
170  cfg.fakeDataSet = std::stoul(optstring.substr(start, end-start));
171  }
172  } // if fake signal
173  else throw std::runtime_error("Fake data type not recognised!");
174  }
175  else throw std::runtime_error("Fake data type not specified!");
176 
177  // Configure oscillation parameters
178  if (CheckOption(opt, "th24vsth34")) { // Th24 vs Th34
179  cfg.xVar = &kFitTheta34InDegreesSterile;
180  cfg.xBins = new Binning(Binning::Simple(50, 0, 50));
181  cfg.yVar = &kFitTheta24InDegreesSterile;
182  cfg.yBins = new Binning(Binning::Simple(50, 0, 45));
183  // cfg.fitSeeds[&kFitDmSq41Sterile] = { 0.5 };
184  } else if (CheckOption(opt, "th24vsdm41")) { // Th24 vs Dm41
185  cfg.xVar = &kFitSinSqTheta24Sterile;
186  cfg.xBins = new Binning(Binning::LogUniform(50, 1e-4, 1));
187  cfg.yVar = &kFitDmSq41Sterile;
188  cfg.yBins = new Binning(Binning::LogUniform(50, 1e-2, 100));
189  if (!CheckOption(opt, "noprofssth34")) {
190  cfg.fitVars.push_back(&kFitSinSqTheta34Sterile);
191  cfg.fitSeeds[&kFitSinSqTheta34Sterile] = { 0., 0.5, 1. };
192  }
193  } else if (CheckOption(opt, "th34vsdm41")) { // Th34 vs Dm41
194  cfg.xVar = &kFitSinSqTheta34Sterile;
195  cfg.xBins = new Binning(Binning::Simple(50, 0, 1));
196  cfg.yVar = &kFitDmSq41Sterile;
197  cfg.yBins = new Binning(Binning::LogUniform(50, 1e-2, 100));
198  if (!CheckOption(opt, "noprofssth24")) {
199  cfg.fitVars.push_back(&kFitSinSqTheta24Sterile);
200  cfg.fitSeeds[&kFitSinSqTheta24Sterile] = { 0., 1. };
201  }
202  } else {
203  throw std::runtime_error("You must specify a type of fit! (\"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\")");
204  }
205 
206  // Common fit vars
207  if (!CheckOption(opt, "noprofth23")) {
208  cfg.fitVars.push_back(&kFitTheta23Sterile);
209  }
210  if (!CheckOption(opt, "noprofdm32")) {
211  cfg.fitVars.push_back(&kFitDmSq32Sterile);
212  }
213  if (!CheckOption(opt, "noprofdelta24")) {
214  cfg.fitVars.push_back(&kFitDelta24InPiUnitsSterile);
215  }
216  // Set seeds regardless of whether they're being fit
217  cfg.fitSeeds[&kFitTheta23Sterile] = { kNux20Th23 };
218  cfg.fitSeeds[&kFitDmSq32Sterile] = { kNux20Dm32 };
219  cfg.fitSeeds[&kFitDelta24InPiUnitsSterile] = { 0.5 };
220 
221  // Set fit type
222  if (CheckOption(opt, "loglikelihood")) {
223  cfg.fitType = kPoissonLL;
224  } else if (CheckOption(opt, "neymanpearson")) {
225  cfg.fitType = kCNP;
226  } else {
227  cfg.fitType = kDynamicScaling;
228  }
229 
230  // Add in matrix epsilon, if necessary
231  cfg.epsilon = 5e-6;
232  pos = optstring.find("epsilon=");
233  if (pos != std::string::npos) {
234  if (cfg.fitType != kPoissonLL) {
235  throw std::runtime_error("Epsilon value should only be provided for Poisson likelihood matrix method!");
236  }
237  size_t start = pos + 8;
238  size_t end = optstring.find("_", start);
239  cfg.epsilon = std::stod(optstring.substr(start, end-start));
240  } // if epsilon provided
241 
242  return cfg;
243 }
const IFitVar * yVar
Definition: FitUtilities.h:52
const double kNux20Th23
Definition: Utilities.h:40
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:366
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
std::vector< covmx::Sample > samples
Definition: FitUtilities.h:44
double stod(const std::string &val)
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:359
const FitDmSq32Sterile kFitDmSq32Sterile
unsigned int fakeDataUniverse
Definition: FitUtilities.h:48
SystType
FitUtilities.h.
Definition: FitUtilities.h:4
Definition: config.py:1
FitType
Definition: FitUtilities.h:27
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
FakeDataType fakeDataType
Definition: FitUtilities.h:46
std::map< const IFitVar *, std::vector< double > > fitSeeds
Definition: FitUtilities.h:56
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
std::string systName
Definition: FitUtilities.h:40
if(dump)
SystConfig systConfig
Definition: FitUtilities.h:39
unsigned int nJobs
Definition: FitUtilities.h:37
unsigned int fakeDataSet
Definition: FitUtilities.h:47
FitType fitType
Definition: FitUtilities.h:58
bool CheckOption(TString opts, TString opt)
unsigned int jobNumber
Definition: FitUtilities.h:36
OStream cout
Definition: OStream.cxx:6
double epsilon
Definition: FitUtilities.h:42
SystType systType
Definition: FitUtilities.h:38
FakeDataType
Definition: FitUtilities.h:19
const IFitVar * xVar
Definition: FitUtilities.h:50
std::vector< covmx::Sample > GetSamplesFromOptString(TString optString)
Function to take an option TString and return a vector of associated covmx::Samples.
Definition: Utilities.h:384
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
SystConfig
Definition: FitUtilities.h:12
Binning * xBins
Definition: FitUtilities.h:51
Binning * yBins
Definition: FitUtilities.h:53
const double kNux20Dm32
Definition: Utilities.h:35
FitConfig GetConfig(TString opt)
Parse an option string into a fit configuration.
Definition: FitUtilities.h:62
Float_t e
Definition: plot.C:35
std::vector< const IFitVar * > fitVars
Definition: FitUtilities.h:55
const FitDmSq41Sterile kFitDmSq41Sterile
const FitTheta23Sterile kFitTheta23Sterile
enum BeamMode string