Utilities.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
2 #include "CAFAna/Core/Ratio.h"
4 
6 #include "Utilities/func/Stan.h"
7 
8 #include "TArrayD.h"
9 #include "TClass.h"
10 #include "TDirectory.h"
11 #include "TH2.h"
12 #include "TH3.h"
13 #include "TF1.h"
14 #include "TMatrixD.h"
15 #include "TObjString.h"
16 #include "TString.h"
17 #include "TVector3.h"
18 #include "TVectorD.h"
19 
20 #include <cassert>
21 #include <cmath>
22 #include <fstream>
23 #include <unistd.h>
24 #include "sys/stat.h"
25 
26 namespace ana
27 {
28  //----------------------------------------------------------------------
30  {
31  static int N = 0;
32  return TString::Format("cafanauniq%d", N++).Data();
33  }
34 
35  //----------------------------------------------------------------------
37  {
38  fBackup = TH1::AddDirectoryStatus();
39  TH1::AddDirectory(false);
40  }
41 
42  //----------------------------------------------------------------------
44  {
45  TH1::AddDirectory(fBackup);
46  }
47 
48  //----------------------------------------------------------------------
50  {
51  const char* s = getenv("IFDH_SILENT");
52  fSet = (s && s == std::string("0"));
53  if(!fSet) setenv("IFDH_SILENT", "1", 1);
54  if(s) std::cout << "IFDH_SILENT=" << s << std::endl;
55  }
56 
57  //----------------------------------------------------------------------
59  {
60  if(!fSet) unsetenv("IFDH_SILENT");
61  }
62 
63  //----------------------------------------------------------------------
65  {
66  // Don't want any pending FPEs to trigger when we flip exceptions
67  // on. Whoever had them off previously had a reason.
68  feclearexcept(FE_INVALID);
69 
70  fegetexceptflag(&fBackup, FE_INVALID);
71 
72  if(enable)
73  feenableexcept(FE_INVALID);
74  else
75  fedisableexcept(FE_INVALID);
76  }
77 
78  //----------------------------------------------------------------------
80  {
81  fesetexceptflag(&fBackup, FE_INVALID);
82  }
83 
84  //----------------------------------------------------------------------
85  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin, int lastBin)
86  {
87  if (binSets.size() < 2)
88  return std::unique_ptr<TMatrixD>(nullptr);
89 
90  if (lastBin < 0)
91  lastBin = binSets[0]->GetSize() - 1; // indexed from 0
92 
93  int nBins = lastBin - firstBin + 1; // firstBin and lastBin are inclusive
94 
95  std::vector<double> binMeans(nBins);
96  for( const auto & binSet : binSets )
97  {
98  for ( decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++ )
99  binMeans[binIdx] += (*binSet)[binIdx];
100  }
101  for (decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++)
102  binMeans[binIdx] /= binSets.size();
103 
104 
105  auto covmx = std::make_unique<TMatrixD>(nBins, nBins);
106 
107  for( unsigned int hist_idx = 0; hist_idx < binSets.size(); ++hist_idx )
108  {
109  // first calculate the weighted sum of squares of the deviations
110  for( decltype(nBins) i = 0; i < nBins; i++ )
111  {
112  double xi = (*(binSets[hist_idx]))[i];
113  for( decltype(nBins) k = i; k < nBins; k++ )
114  {
115  double xk = (*(binSets[hist_idx]))[k];
116  (*covmx)[i][k] += (xi - binMeans[i]) * (xk - binMeans[k]);
117  if (i != k)
118  (*covmx)[k][i] = (*covmx)[i][k]; // covariance matrices are always symmetric
119  }
120  }
121  } // for (hist_idx)
122 
123  // now divide by N-1 to get sample covariance
124  (*covmx) *= 1./(binSets.size()-1);
125 
126  return covmx;
127  }
128 
129  //----------------------------------------------------------------------
130  double LogLikelihood(const TH1* eh, const TH1* oh, bool useOverflow)
131  {
132  assert(eh->GetNbinsX() == oh->GetNbinsX());
133 
134  double chi = 0;
135 
136  int bufferBins = useOverflow? 2 : 1;
137 
138  for(int i = 0; i < eh->GetNbinsX()+bufferBins; ++i){
139  const double e = eh->GetBinContent(i);
140  const double o = oh->GetBinContent(i);
141 
142  chi += LogLikelihood(e, o);
143  }
144 
145  return chi;
146  }
147 
148  //----------------------------------------------------------------------
150  int nbinsx, double xmin, double xmax, bool xlog,
151  int nbinsy, double ymin, double ymax, bool ylog)
152  {
153  DontAddDirectory guard;
154 
155  if(xlog){xmin = log(xmin); xmax = log(xmax);}
156  if(ylog){ymin = log(ymin); ymax = log(ymax);}
157 
158  // How wide the bins will be once we're done
159  const double xwidth = (xmax-xmin)/(nbinsx-1);
160  const double ywidth = (ymax-ymin)/(nbinsy-1);
161 
162  // Move the bin edges so that the limits occur at the centres
163  xmin -= xwidth/2; ymin -= ywidth/2;
164  xmax += xwidth/2; ymax += ywidth/2;
165 
166  std::vector<double> xedges(nbinsx+1);
167  std::vector<double> yedges(nbinsy+1);
168 
169  for(int i = 0; i <= nbinsx; ++i){
170  xedges[i] = xmin + (xmax-xmin)*i/double(nbinsx);
171  if(xlog) xedges[i] = exp(xedges[i]);
172  }
173  for(int i = 0; i <= nbinsy; ++i){
174  yedges[i] = ymin + (ymax-ymin)*i/double(nbinsy);
175  if(ylog) yedges[i] = exp(yedges[i]);
176  }
177 
178  return new TH2F(UniqueName().c_str(), title.c_str(),
179  nbinsx, &xedges.front(),
180  nbinsy, &yedges.front());
181  }
182 
183  //----------------------------------------------------------------------
185  {
186 #ifdef NOVACMAKE
187  const char* pub = getenv("NOVASOFT_DIR");
188  const char* priv = getenv("NOVASOFT_DIR");
189 #else
190  const char* pub = getenv("FW_BASE");
191  const char* priv = getenv("FW_RELEASE_BASE");
192 #endif
193 
194  if(priv && priv != std::string(".")){
195  const std::string ret = std::string(priv)+"/"+Dir;
196  struct stat junk;
197  if(stat(ret.c_str(), &junk) == 0) return ret;
198  }
199  assert(pub);
200  return std::string(pub)+"/"+Dir;
201  }
202 
203  //----------------------------------------------------------------------
205  {
206  return FindPackageDir("CAFAna");
207  }
208 
209  //----------------------------------------------------------------------
210  std::vector<std::string> LoadFileList(const std::string& listfile)
211  {
212  std::vector<std::string> ret;
213 
214  std::ifstream is(listfile);
215  if(!is.good()){
216  std::cerr << "Can't open file list '" << listfile << "'. Aborting." << std::endl;
217  abort();
218  }
219 
220  while(!is.eof()){
222  is >> fname;
223  if(!fname.empty()) ret.push_back(fname);
224  }
225  return ret;
226  }
227 
228  //----------------------------------------------------------------------
229  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir)
230  {
231  std::map<std::string, std::string> ret;
232 
233  TIter next(dir->GetListOfKeys());
234  while(TObject* key = next()){
235  TObject* obj = dir->Get(key->GetName());
236  assert(obj);
237 
238  const TString className = obj->ClassName();
239  if(className == "TObjString"){
240  ret[key->GetName()] = ((TObjString*)obj)->GetString();
241  }
242  else{
243  std::cerr << "Unexpected object " << key->GetName() << " of type " << className << " while looking for metadata. Ignoring" << std::endl;
244  }
245  }
246 
247  return ret;
248  }
249 
250  //----------------------------------------------------------------------
251  void CombineMetadata(std::map<std::string, std::string>& base,
252  const std::map<std::string, std::string>& add,
253  std::set<std::string>& mask)
254  {
255  for(auto it: add){
256  const std::string& key = it.first;
257 
258  // Needs special handling anyway, leave it blank.
259  if(key == "parents") continue;
260 
261  // Accumulate the runs list
262  if(key == "runs"){
263  const std::string& r1 = base[key];
264  const std::string& r2 = it.second;
265 
266  assert(!r2.empty());
267 
268  // when concatenating MC files,
269  // it's possible to get duplicates
270  // (run+subrun the same, only differing
271  // in cycle number, which isn't part
272  // of this metadata field)
273  // don't repeat, because SAM gets angry.
274  if (r1.find(r2.substr(1, r2.size() - 2)) != std::string::npos)
275  continue;
276 
277  if(r1.empty()){
278  base[key] = r2;
279  continue;
280  }
281 
282  // "[foo]" + "[bar]"
283  std::string sum = r1+&r2[1]; // "[foo]bar]"
284  sum[r1.size()-1] = ','; // "[foo,bar]"
285  base[key] = sum;
286  continue;
287  }
288 
289  if(base.find(key) == base.end()){
290  // If it's new, add it
291  base[key] = it.second;
292  }
293  else{
294  if(key == "simulated.number_of_spills" ||
295  key == "event_count" ||
296  key == "online.totalevents"){
297  // These two fields should be accumulated
298  base[key] = TString::Format("%d",
299  atoi(base[key].c_str()) +
300  atoi(it.second.c_str())).Data();
301  }
302  else{
303  // If it's a clash, record it
304  if(base[key] != it.second) mask.insert(key);
305  }
306  }
307  }
308  }
309 
310 
311  //----------------------------------------------------------------------
312  void WriteCAFMetadata(TDirectory* dir,
313  const std::map<std::string, std::string>& meta)
314  {
315  TDirectory* tmp = gDirectory;
316  dir->cd();
317 
318  for(auto it: meta){
319  TObjString str(it.second.c_str());
320  str.Write(it.first.c_str());
321  }
322 
323  dir->Save();
324 
325  tmp->cd();
326  }
327 
328  //----------------------------------------------------------------------
330  {
331  // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
332  // CAF files and fitting the resulting distribution to a gaussian.
333  // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
334  // and containment cuts.
335  // CAF Files used:
336  // NearDet:
337  // prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
338  // FarDet:
339  // prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08
340  // Unknown for NDOS
341 
342  switch(det){
343  case caf::kFARDET: return TVector3(-6.833e-05,
344  +6.388e-02,
345  +9.980e-1).Unit();
346  case caf::kNDOS: return TVector3(-3.845e-4,
347  +6.517e-2,
348  +9.967e-1).Unit();
349  case caf::kNEARDET: return TVector3(-8.424e-04,
350  -6.174e-02,
351  +9.981e-1).Unit();
352  default:
353  std::cerr<< "Aborting. Bad geometry configuration: " << det << "\n";
354  abort();
355  }
356  }
357 
358  //----------------------------------------------------------------------
360  {
361  static bool cache;
362  static bool cache_set = false;
363  if(!cache_set){
364  cache = (getenv("_CONDOR_SCRATCH_DIR") != 0);
365  cache_set = true;
366  }
367 
368  return cache;
369  }
370 
371  //----------------------------------------------------------------------
372  size_t Stride(bool allow_default)
373  {
374  static int cache = -1;
375 
376  if(cache < 0){
377  char* env = getenv("CAFANA_STRIDE");
378  if(env){
379  cache = std::atoi(env);
380  }
381  else{
382  if(allow_default){
383  cache = 1;
384  }
385  else{
386  std::cout << "Stride() called, but CAFANA_STRIDE is not set (--stride not passed?)" << std::endl;
387  abort();
388  }
389  }
390  }
391 
392  return cache;
393  }
394 
395  //----------------------------------------------------------------------
396  size_t Offset(bool allow_default)
397  {
398  static int cache = -1;
399 
400  if(cache < 0){
401  char* env = getenv("CAFANA_OFFSET");
402  if(env){
403  cache = std::atoi(env);
404  }
405  else{
406  if(allow_default){
407  cache = 0;
408  }
409  else{
410  std::cout << "Offset() called, but CAFANA_OFFSET is not set (--offset not passed?)" << std::endl;
411  abort();
412  }
413  }
414  }
415 
416  return cache;
417  }
418 
419  //----------------------------------------------------------------------
420  int Limit()
421  {
422  static int cache = 0;
423 
424  if(cache == 0){
425  char* env = getenv("CAFANA_LIMIT");
426  if(env){
427  cache = std::atoi(env);
428  }
429  else{
430  cache = -1;
431  }
432  }
433 
434  return cache;
435  }
436 
437  //----------------------------------------------------------------------
438  size_t JobNumber()
439  {
440  if(!RunningOnGrid()){
441  std::cout << "JobNumber() called, but we are not running on the grid" << std::endl;
442  abort();
443  }
444 
445  return Offset(false);
446  }
447 
448  //----------------------------------------------------------------------
449  size_t NumJobs()
450  {
451  if(!RunningOnGrid()){
452  std::cout << "NumJobs() called, but we are not running on the grid" << std::endl;
453  abort();
454  }
455 
456  return Stride(false);
457  }
458 
459  //----------------------------------------------------------------------
460  std::size_t Hash(const std::vector<std::size_t> &vals)
461  {
462  std::size_t seed = 0;
463 
464  for (auto v : vals)
465  seed ^= std::hash<std::size_t>{}(v) + 0x9e3779b9 + (seed<<6) + (seed>>2);
466  return seed;
467  }
468 
469  //----------------------------------------------------------------------
470  std::size_t Hash(const caf::SRSpillProxy *sr)
471  {
472  return Hash({sr->run, sr->subrun, sr->evt});
473  }
474 
475  //----------------------------------------------------------------------
476  std::size_t Hash(const caf::SRProxy *sr)
477  {
478  return Hash({Hash(&(sr->spill)), std::size_t(sr->hdr.cycle), sr->hdr.subevt});
479  }
480 
481  //----------------------------------------------------------------------
482  FitToFourier::FitToFourier(TH1* h, double xlo, double xhi, int NOsc)
483  : fHist(h), fxlo(xlo), fxhi(xhi), fNOsc(NOsc)
484  {
485  }
486 
487  //----------------------------------------------------------------------
489  {
490  }
491 
492  //----------------------------------------------------------------------
493  double FitToFourier::operator()(double *x, double *par) const
494  {
495  double x0 = x[0];
496  double val = par[0];
497  for (int i = 1; i <= fNOsc; i++)
498  val += par[2*i-1]*sin(i*M_PI*x0) + par[2*i]*cos(i*M_PI*x0);
499  return val;
500  }
501 
502  //----------------------------------------------------------------------
503  TF1* FitToFourier::Fit() const
504  {
505  std::vector<double> s(fNOsc, 0.);
506  std::vector<double> c(fNOsc, 0.);
507  int nBins = 0;
508  for(int i = 1; i <= fHist->GetNbinsX(); ++i){
509  const double x = M_PI * fHist->GetXaxis()->GetBinCenter(i);
510  const double y = fHist->GetBinContent(i);
511 
512  if(y == 0) continue;
513  ++nBins;
514 
515  for(int n = 0; n <= fNOsc; ++n){
516  s[n] += y * sin(n*x);
517  c[n] += y * cos(n*x);
518  }
519  }
520 
521  for(int n = 0; n <= fNOsc; ++n){
522  s[n] *= 2./nBins;
523  c[n] *= 2./nBins;
524  }
525 
526  TF1* f = new TF1(UniqueName().c_str(), this, fxlo, fxhi, 2*fNOsc+1);
527 
528  f->SetParameter(0, c[0]/2);
529  for(int n = 1; n <= fNOsc; ++n){
530  f->SetParameter(n*2-1, s[n]);
531  f->SetParameter(n*2, c[n]);
532  }
533 
534  // Because ROOT is having problems drawing f if I don't
535  double min = fHist->GetMinimum();
536  double max = fHist->GetMaximum();
537  f->GetYaxis()->SetRangeUser(0.8*min, 1.2*max);
538  return f;
539  }
540 }
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2143
Near Detector underground.
Definition: SREnums.h:10
const double fxhi
Definition: Utilities.h:221
Det_t
Which NOvA detector?
Definition: SREnums.h:7
int nBins
Definition: plotROC.py:16
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:438
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
set< int >::iterator it
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:359
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Int_t par
Definition: SimpleIterate.C:24
void CombineMetadata(std::map< std::string, std::string > &base, const std::map< std::string, std::string > &add, std::set< std::string > &mask)
base += add
Definition: Utilities.cxx:251
OStream cerr
Definition: OStream.cxx:7
FitToFourier(TH1 *h, double xlo, double xhi, int NOsc)
Definition: Utilities.cxx:482
Float_t tmp
Definition: plot.C:36
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
#define M_PI
Definition: SbMath.h:34
double operator()(double *x, double *par) const
Definition: Utilities.cxx:493
FloatingExceptionOnNaN(bool enable=true)
Definition: Utilities.cxx:64
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
Definition: Utilities.cxx:149
Double_t ymax
Definition: plot.C:25
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
Definition: Utilities.cxx:210
const XML_Char * s
Definition: expat.h:262
int Limit()
Value passed to –limit, or -1 if not specified.
Definition: Utilities.cxx:420
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:250
unsigned int seed
Definition: runWimpSim.h:102
Prototype Near Detector on the Surface.
Definition: SREnums.h:12
std::unique_ptr< TMatrixD > CalcCovMx(const std::vector< TArrayD * > &binSets, int firstBin, int lastBin)
Compute bin-to-bin covariance matrix from a collection of sets of bin contents.
Definition: Utilities.cxx:85
const double fxlo
Definition: Utilities.h:220
std::string getenv(std::string const &name)
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:1367
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Class for parsing *.fcl files for their metadata information.
void WriteCAFMetadata(TDirectory *dir, const std::map< std::string, std::string > &meta)
Write map of metadata parameters into a CAF file.
Definition: Utilities.cxx:312
caf::StandardRecord * sr
size_t Stride(bool allow_default)
Value passed to –stride, or 1 if not specified.
Definition: Utilities.cxx:372
caf::Proxy< int > cycle
Definition: SRProxy.h:230
const TH1 * fHist
Definition: Utilities.h:219
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:184
TF1 * Fit() const
Definition: Utilities.cxx:503
caf::Proxy< unsigned int > subrun
Definition: SRProxy.h:1410
Proxy for caf::SRSpill.
Definition: SRProxy.h:1346
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:460
size_t Offset(bool allow_default)
Value passed to –offset, or 0 if not specified.
Definition: Utilities.cxx:396
T sin(T number)
Definition: d0nt_math.hpp:132
pub
Definition: cvnie.py:6
caf::Proxy< unsigned int > run
Definition: SRProxy.h:1406
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
size_t NumJobs()
Definition: Utilities.cxx:449
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Extract map of metadata parameters from a CAF file.
Definition: Utilities.cxx:229
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
const int fNOsc
Definition: Utilities.h:222
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:329
Double_t ymin
Definition: plot.C:24
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
priv
Definition: cvnie.py:5
Double_t sum
Definition: plot.C:31
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
add("abs", expr_type(int_type()), expr_type(int_type()))
void next()
Definition: show_event.C:84
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string