gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Flav_reader.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
26 
30 
31 
32 namespace Gambit
33 {
34 
35  namespace FlavBit
36  {
37 
38  using namespace std;
39 
41  void operator >> (const YAML::Node& node, Correlation& c)
42  {
43  // safety
44  string name = node["name"].as<std::string>();
45  if (name == "NONE") return;
46  // reading correlation
47  c.corr_name= name;
48  c.corr_val = node["value"].as<double>();
49  }
50 
52  void operator >> (const YAML::Node& node, Measurement& v)
53  {
54  v.name=node["name"].as<std::string>();
55  v.is_limit = node["islimit"].as<bool>();
56  v.exp_value= node["exp_value"].as<double>();
57  if(node["exp_90_CL"]) v.exp_90_CL = node["exp_90_CL"].as<double>();
58  if(node["exp_95_CL"]) v.exp_95_CL = node["exp_95_CL"].as<double>();
59  if(node["one_sided"]) v.exp_one_sided = node["one_sided"].as<bool>();
60  v.exp_source = node["exp_source"].as<std::string>();
61  v.exp_stat_error = node["exp_stat_error"].as<double>();
62  v.exp_sys_error = node["exp_sys_error"].as<double>();
63  if(v.is_limit and !v.exp_stat_error and (v.exp_90_CL or v.exp_95_CL))
64  {
65  v.exp_stat_error =
67  }
69  v.th_error=node["th_error"].as<double>();
70  v.th_error_type=node["th_error_type"].as<std::string>();
71  v.th_error_source = node["th_error_source"].as<std::string>();
72  const YAML::Node& correlations = node["correlation"];
73  v.corr.clear();
74  for (unsigned i=0; i<correlations.size(); ++i)
75  {
76  Correlation corr_tmp;
77  correlations[i] >> corr_tmp;
78  if (corr_tmp.corr_name != "") v.corr.push_back(corr_tmp);
79  }
80  }
81 
84  {
85  measurements=vector< Measurement >(0);
86  measurement_location=loc;
87  debug=false;
88  use_S=true;
89  use_P=false;
90  number_measurements=0;
91  }
92 
94  void Flav_reader::read_yaml(string name)
95  {
96  string path=measurement_location+"/"+name;
97  std::ifstream fin(path.c_str());
98  YAML::Node doc = YAML::Load(fin);
99  number_measurements=0;
100  for(unsigned i=0;i<doc.size();++i)
101  {
102  Measurement mes_tmp;
103  doc[i] >> mes_tmp;
104  if(debug) print(mes_tmp);
105  measurements.push_back(mes_tmp);
106  number_measurements++;
107  }
108  if (debug) cout<<"Number of measurements: "<<number_measurements<<endl;
109  }
110 
112  void Flav_reader::read_yaml_measurement(string name, string measurement_name)
113  {
114  string path=measurement_location+"/"+name;
115  std::ifstream fin(path.c_str());
116  YAML::Node doc = YAML::Load(fin);
117  Measurement mes_tmp;
118 
119  if (debug) cout << measurement_name.c_str() << endl;
120 
121  for (unsigned i=0; i<doc.size(); ++i)
122  {
123  Measurement mes_tmp;
124  doc[i] >> mes_tmp;
125  if(mes_tmp.name!=measurement_name.c_str()) continue;
126  measurements.push_back(mes_tmp);
127  number_measurements++;
128  if (debug) cout << "Read in " << measurement_name << " " << mes_tmp.name << " " << mes_tmp.exp_value << endl;
129  }
130  if(debug) cout << "Number of measurements: " << number_measurements << endl;
131  }
132 
135  {
136  cout<<"################### Mesurement"<<endl;
137  cout<<"Name: "<<mes.name<<endl;
138  cout<<(mes.is_limit ? "Limit" : "Value")<<": "<<mes.exp_value<<endl;
139  cout<<"Stat/sys errror: "<< mes.exp_stat_error<<"/"<<mes.exp_sys_error<<endl;
140  cout<<"Error plus/minus: "<<mes.exp_error<<endl;
141  cout<<"Experimental source: "<<mes.exp_source<<endl;
142  cout<<"Correlations:"<<endl;
143  for(unsigned i=0;i<mes.corr.size();++i)
144  {
145  cout<<mes.corr[i].corr_name<<" "<<mes.corr[i].corr_val<<endl;
146  }
147  cout<<"########## END"<<endl;
148  }
149 
152  {
153  // Create the correlation matrix
154  M_cor_cov = boost::numeric::ublas::identity_matrix<double>(number_measurements);
155  for(int i=0; i<number_measurements; ++i)
156  {
157  #ifdef FLAVBIT_DEBUG
158  cout<<"Correlation size: "<< measurements[i].corr.size()<<endl;
159  #endif
160  for ( unsigned icorr=0; icorr< measurements[i].corr.size(); ++icorr)
161  {
162  #ifdef FLAVBIT_DEBUG
163  cout<<"Searching for correlation: "<< measurements[i].corr[icorr].corr_name <<endl;
164  #endif
165  int i_corr_index=get_measurement_for_corr(measurements[i].corr[icorr].corr_name );
166  M_cor_cov(i_corr_index,i)=measurements[i].corr[icorr].corr_val;
167  }
168  }
169  if (debug) print_matrix(M_cor_cov,"Correlation Matrix:");
170 
171  // Make sure the correlation matrix is symmetric and has 1s on the diagonal
172  check_corr_matrix(M_cor_cov);
173 
174  // Convert the correlation matrix to a covariance matrix
175  for(int i=0; i<number_measurements;++i)
176  {
177  for(int j=0;j<number_measurements;++j)
178  {
179  M_cor_cov(i,j)=M_cor_cov(i,j)*(measurements[i].exp_error)*(measurements[j].exp_error);
180  }
181  }
182  if (debug) print_matrix(M_cor_cov,"Covariance Matrix:");
183 
184  // Construct the the measurement vector
185  M_measurements= boost::numeric::ublas::matrix<double> (number_measurements,1);
186  for(int i=0; i<number_measurements;++i)
187  {
188  M_measurements(i,0)=measurements[i].exp_value;
189  }
190  if (debug) print_matrix(M_measurements, "Measurements:", false);
191 
192  // Construct the theory error vector
193  M_th_err = boost::numeric::ublas::matrix< std::pair<double,bool> > (number_measurements,1);
194  for(int i=0; i<number_measurements;++i)
195  {
196  M_th_err(i,0).first = measurements[i].th_error;
197  str err_type = measurements[i].th_error_type;
198  if (err_type == "A") M_th_err(i,0).second = true;
199  else if (err_type == "M") M_th_err(i,0).second = false;
200  else FlavBit_error().raise(LOCAL_INFO, "Unrecognised theory error type in database: "+err_type);
201  }
202  if (debug) print_matrix(M_th_err, "Theory errors:", false);
203 
204  }
205 
208  {
209  for(int i=0;i<number_measurements;++i)
210  {
211  if(measurements[i].name == ss) return i;
212  }
213  cout<<"Error!, didn't find measuremnet: "<<ss<<endl;
214  return -1;
215  }
216 
218  void Flav_reader::print_matrix(boost::numeric::ublas::matrix<double>& M, str name, bool is_true_matrix)
219  {
220  int jmax = is_true_matrix ? number_measurements : 1;
221  cout<<name<<endl;
222  for(int i=0; i < number_measurements; ++i)
223  {
224  for(int j=0 ; j< jmax; ++j) cout<<M(i,j)<<"\t";
225  cout<<endl;
226  }
227  }
228 
230  void Flav_reader::print_matrix(boost::numeric::ublas::matrix< std::pair<double, bool> >& M, str name, bool is_true_matrix)
231  {
232  int jmax = is_true_matrix ? number_measurements : 1;
233  cout<<name<<endl;
234  for(int i=0; i < number_measurements; ++i)
235  {
236  for(int j=0 ; j< jmax; ++j) cout<<M(i,j).first<<":"<<M(i,j).second<<"\t";
237  cout<<endl;
238  }
239  }
240 
242  void Flav_reader::check_corr_matrix(boost::numeric::ublas::matrix<double>& M)
243  {
244  for( int i=0; i < number_measurements; ++i)
245  {
246  for( int j=0 ; j< number_measurements; ++j)
247  {
248  if (i==j && M(i,i)!=1)
249  {
250  FlavBit_error().raise(LOCAL_INFO, "Correlation matrix diagonal element differs from 1!");
251  }
252  if (M(i,j) != M(j,i))
253  {
254  cout << "Correlation matrix " << i << "," << j << " = " << M(i,j) << endl;
255  cout << "Correlation matrix " << i << "," << j << " = " << M(j,i) << endl;
256  FlavBit_error().raise(LOCAL_INFO, "Correlation matrix not symmetric");
257  }
258  }
259  }
260  }
261 
262  double Flav_reader::get_error_from_confidence_levels(double exp_value, double CL_90, double CL_95, bool one_sided)
263  {
264  double error = 0;
265  if(one_sided)
266  {
267  if(CL_90 != 0)
268  error = fabs(exp_value - CL_90)/1.28;
269  else if (CL_95 != 0)
270  error =fabs(exp_value - CL_95)/1.64;
271  }
272  else
273  {
274  if(CL_90 != 0)
275  error = fabs(exp_value - CL_90)/1.64;
276  else if(CL_95 != 0)
277  error = fabs(exp_value - CL_95)/1.96;
278  }
279  return error;
280 
281  }
282  }
283 }
void print(MixMatrix A)
Helper function to print a matrix.
void print_matrix(boost::numeric::ublas::matrix< double > &, str, bool is_true_matrix=true)
Print a boost ublas matrix.
static double get_error_from_confidence_levels(double, double, double, bool)
Calculates the experimental statistical error from confidence levels.
#define LOCAL_INFO
Definition: local_info.hpp:34
STL namespace.
T * matrix(const int xN)
Definition: random_tools.hpp:9
void print(Measurement)
Print a measurement previously read in from the database.
Flav_reader(str loc)
Constructor that takes the location of the database as an argument.
Definition: Flav_reader.cpp:83
GAMBIT error class.
Definition: exceptions.hpp:136
Representation of a single entry in the FlavBit YAML database.
dictionary fin
Rollcall header for module FlavBit.
void initialise_matrices()
Compute the covariance matrix and populate the measurement and theory error vectors.
void operator>>(const YAML::Node &node, Correlation &c)
Extraction operator for correlation.
Definition: Flav_reader.cpp:41
void read_yaml(str name)
Read the entire database into memory.
Definition: Flav_reader.cpp:94
Header file that includes all GAMBIT headers required for a module source file.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
START_MODEL M
void read_yaml_measurement(str name, str measurement_name)
Read a single measurement from the database into memory.
int get_measurement_for_corr(str)
Find the second measurement that corresponds to a given correlation.
void check_corr_matrix(boost::numeric::ublas::matrix< double > &)
Check that a correlation matrix has 1s on the diagonal and is symmetric.
Declaration of reader class for FlavBit YAML database.
Simple structure for holding a correlation value and name of the correlated observable.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
std::vector< Correlation > corr