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 #include "wordexp.h"
26 
27 namespace ana
28 {
29  //----------------------------------------------------------------------
31  {
32  static int N = 0;
33  return TString::Format("cafanauniq%d", N++).Data();
34  }
35 
36  //----------------------------------------------------------------------
38  {
39  fBackup = TH1::AddDirectoryStatus();
40  TH1::AddDirectory(false);
41  }
42 
43  //----------------------------------------------------------------------
45  {
46  TH1::AddDirectory(fBackup);
47  }
48 
49  //----------------------------------------------------------------------
51  {
52  const char* s = getenv("IFDH_SILENT");
53  fSet = (s && s == std::string("0"));
54  if(!fSet) setenv("IFDH_SILENT", "1", 1);
55  if(s) std::cout << "IFDH_SILENT=" << s << std::endl;
56  }
57 
58  //----------------------------------------------------------------------
60  {
61  if(!fSet) unsetenv("IFDH_SILENT");
62  }
63 
64  //----------------------------------------------------------------------
66  {
67  // Don't want any pending FPEs to trigger when we flip exceptions
68  // on. Whoever had them off previously had a reason.
69  feclearexcept(FE_INVALID);
70 
71  fegetexceptflag(&fBackup, FE_INVALID);
72 
73 #ifndef DARWINBUILD
74  if(enable)
75  feenableexcept(FE_INVALID);
76  else
77  fedisableexcept(FE_INVALID);
78 #else
79  std::cerr << "WARNING: CAFAna/Core/Utilities.cxx built on OS X, no feenableexcept available" << std::endl;
80 #endif
81  }
82 
83  //----------------------------------------------------------------------
85  {
86  fesetexceptflag(&fBackup, FE_INVALID);
87  }
88 
89  //----------------------------------------------------------------------
90  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin, int lastBin)
91  {
92  if (binSets.size() < 2)
93  return std::unique_ptr<TMatrixD>(nullptr);
94 
95  if (lastBin < 0)
96  lastBin = binSets[0]->GetSize() - 1; // indexed from 0
97 
98  int nBins = lastBin - firstBin + 1; // firstBin and lastBin are inclusive
99 
100  std::vector<double> binMeans(nBins);
101  for( const auto & binSet : binSets )
102  {
103  for ( decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++ )
104  binMeans[binIdx] += (*binSet)[binIdx];
105  }
106  for (decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++)
107  binMeans[binIdx] /= binSets.size();
108 
109 
110  auto covmx = std::make_unique<TMatrixD>(nBins, nBins);
111 
112  for( unsigned int hist_idx = 0; hist_idx < binSets.size(); ++hist_idx )
113  {
114  // first calculate the weighted sum of squares of the deviations
115  for( decltype(nBins) i = 0; i < nBins; i++ )
116  {
117  double xi = (*(binSets[hist_idx]))[i];
118  for( decltype(nBins) k = i; k < nBins; k++ )
119  {
120  double xk = (*(binSets[hist_idx]))[k];
121  (*covmx)[i][k] += (xi - binMeans[i]) * (xk - binMeans[k]);
122  if (i != k)
123  (*covmx)[k][i] = (*covmx)[i][k]; // covariance matrices are always symmetric
124  }
125  }
126  } // for (hist_idx)
127 
128  // now divide by N-1 to get sample covariance
129  (*covmx) *= 1./(binSets.size()-1);
130 
131  return covmx;
132  }
133 
134  //----------------------------------------------------------------------
135  double LogLikelihood(const TH1* eh, const TH1* oh, bool useOverflow)
136  {
137  assert(eh->GetNbinsX() == oh->GetNbinsX());
138 
139  double chi = 0;
140 
141  int bufferBins = useOverflow? 2 : 1;
142 
143  for(int i = 0; i < eh->GetNbinsX()+bufferBins; ++i){
144  const double e = eh->GetBinContent(i);
145  const double o = oh->GetBinContent(i);
146 
147  chi += LogLikelihood(e, o);
148  }
149 
150  return chi;
151  }
152 
153  //----------------------------------------------------------------------
155  int nbinsx, double xmin, double xmax, bool xlog,
156  int nbinsy, double ymin, double ymax, bool ylog)
157  {
158  DontAddDirectory guard;
159 
160  if(xlog){xmin = log(xmin); xmax = log(xmax);}
161  if(ylog){ymin = log(ymin); ymax = log(ymax);}
162 
163  // How wide the bins will be once we're done
164  const double xwidth = (xmax-xmin)/(nbinsx-1);
165  const double ywidth = (ymax-ymin)/(nbinsy-1);
166 
167  // Move the bin edges so that the limits occur at the centres
168  xmin -= xwidth/2; ymin -= ywidth/2;
169  xmax += xwidth/2; ymax += ywidth/2;
170 
171  std::vector<double> xedges(nbinsx+1);
172  std::vector<double> yedges(nbinsy+1);
173 
174  for(int i = 0; i <= nbinsx; ++i){
175  xedges[i] = xmin + (xmax-xmin)*i/double(nbinsx);
176  if(xlog) xedges[i] = exp(xedges[i]);
177  }
178  for(int i = 0; i <= nbinsy; ++i){
179  yedges[i] = ymin + (ymax-ymin)*i/double(nbinsy);
180  if(ylog) yedges[i] = exp(yedges[i]);
181  }
182 
183  return new TH2F(UniqueName().c_str(), title.c_str(),
184  nbinsx, &xedges.front(),
185  nbinsy, &yedges.front());
186  }
187 
188  // Helper functions for MakeTHND().
189  namespace{
190  // Eventually the bin parameters will all be unpacked and we just pass them
191  // on to the regular constructor.
192  template<class T, class... A> T* MakeHist(A... a)
193  {
194  DontAddDirectory guard;
195  return new T(a...);
196  }
197 
198  // This function consumes bins from the start of the argument list and
199  // pushes their translations onto the list of arguments at the end.
200  template<class T, class... A> T* MakeHist(const Binning& firstBin,
201  A... args)
202  {
203  if(firstBin.IsSimple())
204  return MakeHist<T>(args...,
205  firstBin.NBins(), firstBin.Min(), firstBin.Max());
206  else
207  return MakeHist<T>(args...,
208  firstBin.NBins(), &firstBin.Edges().front());
209  }
210  }
211 
212  // Concrete instantiations. MakeHist() requires us to put the bin arguments
213  // first...
214  //----------------------------------------------------------------------
215  TH1D* MakeTH1D(const char* name, const char* title, const Binning& bins)
216  {
217  return MakeHist<TH1D>(bins, name, title);
218  }
219 
220  //----------------------------------------------------------------------
221  TH1F* MakeTH1F(const char* name, const char* title, const Binning& bins)
222  {
223  return MakeHist<TH1F>(bins, name, title);
224  }
225 
226  //----------------------------------------------------------------------
227  TH2D* MakeTH2D(const char* name, const char* title,
228  const Binning& binsx,
229  const Binning& binsy)
230  {
231  return MakeHist<TH2D>(binsx, binsy, name, title);
232  }
233 
234  //----------------------------------------------------------------------
235  TH2F* MakeTH2F(const char* name, const char* title,
236  const Binning& binsx,
237  const Binning& binsy)
238  {
239  return MakeHist<TH2F>(binsx, binsy, name, title);
240  }
241 
242  //----------------------------------------------------------------------
243 
244  TH3D* MakeTH3D(const char* name, const char* title,
245  const Binning& binsx,
246  const Binning& binsy,
247  const Binning& binsz)
248  {
249  return MakeHist<TH3D>(name, title,
250  binsx.NBins(), &binsx.Edges().front(),
251  binsy.NBins(), &binsy.Edges().front(),
252  binsz.NBins(), &binsz.Edges().front());
253  }
254 
255  //----------------------------------------------------------------------
256  TH2* ToTH2(const Spectrum& s, double exposure, ana::EExposureType expotype,
257  const Binning& binsx, const Binning& binsy, ana::EBinType bintype)
258  {
259  DontAddDirectory guard;
260 
261  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
262  return ToTH2Helper(std::move(h1), binsx, binsy, bintype);
263  }
264 
265  //----------------------------------------------------------------------
266  TH2* ToTH2(const Ratio& r,
267  const Binning& binsx, const Binning& binsy)
268  {
269  DontAddDirectory guard;
270 
271  std::unique_ptr<TH1> h1(r.ToTH1());
272  return ToTH2Helper(std::move(h1), binsx, binsy);
273  }
274 
275  //----------------------------------------------------------------------
276  TH2* ToTH2Helper(std::unique_ptr<TH1> h1,
277  const Binning& binsx, const Binning& binsy,
278  ana::EBinType bintype)
279  {
280  // Make sure it's compatible with having been made with this binning
281  assert(h1->GetNbinsX() == binsx.NBins()*binsy.NBins());
282 
283  TH2* ret = MakeTH2D("", UniqueName().c_str(), binsx, binsy);
284 
285  for(int i = 0; i < h1->GetNbinsX(); ++i){
286  const double val = h1->GetBinContent(i+1);
287  const double err = h1->GetBinError(i+1);
288 
289  const int ix = i / binsy.NBins();
290  const int iy = i % binsy.NBins();
291 
292  ret->SetBinContent(ix+1, iy+1, val);
293  ret->SetBinError (ix+1, iy+1, err);
294  }
295 
296  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
297 
298  return ret;
299  }
300 
301  //----------------------------------------------------------------------
302 
303  TH3* ToTH3(const Spectrum& s, double exposure, ana::EExposureType expotype,
304  const Binning& binsx, const Binning& binsy, const Binning& binsz,
305  ana::EBinType bintype)
306  {
307  DontAddDirectory guard;
308 
309  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
310 
311  return ToTH3Helper(std::move(h1), binsx, binsy, binsz, bintype);
312  }
313 
314  //----------------------------------------------------------------------
315 
316  TH3* ToTH3(const Ratio& r,
317  const Binning& binsx, const Binning& binsy, const Binning& binsz)
318  {
319  DontAddDirectory guard;
320 
321  std::unique_ptr<TH1> h1(r.ToTH1());
322 
323  return ToTH3Helper(std::move(h1), binsx, binsy, binsz);
324  }
325 
326  //----------------------------------------------------------------------
327  TH3* ToTH3Helper(std::unique_ptr<TH1> h1,
328  const Binning& binsx,
329  const Binning& binsy,
330  const Binning& binsz,
331  ana::EBinType bintype)
332  {
333 
334  const int nx = binsx.NBins();
335  const int ny = binsy.NBins();
336  const int nz = binsz.NBins();
337 
338  // Make sure it's compatible with having been made with this binning
339  assert(h1->GetNbinsX() == nx*ny*nz);
340 
341  TH3* ret;
342 
343  // If all three axes are simple, we can call a simpler constructor
344  if(binsx.IsSimple() && binsy.IsSimple() && binsz.IsSimple()){
345  ret = new TH3F(UniqueName().c_str(), "",
346  nx, binsx.Min(), binsx.Max(),
347  ny, binsy.Min(), binsy.Max(),
348  nz, binsz.Min(), binsz.Max());
349 
350  }
351  else{
352  // If >= 1 axis has a not simple binning, you have to pass the
353  // edges vector to the TH3 constructor
354  ret = new TH3F(UniqueName().c_str(), "",
355  nx, &binsx.Edges().front(),
356  ny, &binsy.Edges().front(),
357  nz, &binsz.Edges().front());
358  }
359 
360  for(int i = 0; i < h1->GetNbinsX(); ++i){
361  const double val = h1->GetBinContent(i+1);
362  const double err = h1->GetBinError(i+1);
363 
364  const int nynz = ny*nz;
365  const int nmodnynz = i%nynz;
366  const int ix = i/nynz;
367  const int iy = nmodnynz/nz;
368  const int iz = i%nz;
369 
370  ret->SetBinContent(ix+1, iy+1, iz+1, val);
371  ret->SetBinError (ix+1, iy+1, iz+1, err);
372  }
373 
374  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
375 
376  return ret;
377 
378  }
379 
380  //----------------------------------------------------------------------
381  std::vector<std::string> Wildcard(const std::string& wildcardString)
382  {
383  // Expand environment variables and wildcards like the shell would
384  wordexp_t p;
385  const int status = wordexp(wildcardString.c_str(), &p, WRDE_SHOWERR);
386 
387  if(status != 0){
388  std::cerr << "Wildcard string '" << wildcardString
389  << "' returned error " << status << " from wordexp()."
390  << std::endl;
391  return {};
392  }
393 
394  std::vector<std::string> fileList;
395 
396  for(unsigned int i = 0; i < p.we_wordc; ++i){
397  // Check the file exists before adding it
398  struct stat sb;
399  if(stat(p.we_wordv[i], &sb) == 0)
400  fileList.push_back(p.we_wordv[i]);
401  }
402 
403  wordfree(&p);
404 
405  return fileList;
406  }
407 
408  //----------------------------------------------------------------------
410  {
411 #ifdef NOVACMAKE
412  const char* pub = getenv("NOVASOFT_DIR");
413  const char* priv = getenv("NOVASOFT_DIR");
414 #else
415  const char* pub = getenv("FW_BASE");
416  const char* priv = getenv("FW_RELEASE_BASE");
417 #endif
418 
419  if(priv && priv != std::string(".")){
420  const std::string ret = std::string(priv)+"/"+Dir;
421  struct stat junk;
422  if(stat(ret.c_str(), &junk) == 0) return ret;
423  }
424  assert(pub);
425  return std::string(pub)+"/"+Dir;
426  }
427 
428  //----------------------------------------------------------------------
430  {
431  return FindPackageDir("CAFAna");
432  }
433 
434  //----------------------------------------------------------------------
435  std::vector<std::string> LoadFileList(const std::string& listfile)
436  {
437  std::vector<std::string> ret;
438 
439  std::ifstream is(listfile);
440  if(!is.good()){
441  std::cerr << "Can't open file list '" << listfile << "'. Aborting." << std::endl;
442  abort();
443  }
444 
445  while(!is.eof()){
447  is >> fname;
448  if(!fname.empty()) ret.push_back(fname);
449  }
450  return ret;
451  }
452 
453  //----------------------------------------------------------------------
454  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir)
455  {
456  std::map<std::string, std::string> ret;
457 
458  TIter next(dir->GetListOfKeys());
459  while(TObject* key = next()){
460  TObject* obj = dir->Get(key->GetName());
461  assert(obj);
462 
463  const TString className = obj->ClassName();
464  if(className == "TObjString"){
465  ret[key->GetName()] = ((TObjString*)obj)->GetString();
466  }
467  else{
468  std::cerr << "Unexpected object " << key->GetName() << " of type " << className << " while looking for metadata. Ignoring" << std::endl;
469  }
470  }
471 
472  return ret;
473  }
474 
475  //----------------------------------------------------------------------
476  void CombineMetadata(std::map<std::string, std::string>& base,
477  const std::map<std::string, std::string>& add,
478  std::set<std::string>& mask)
479  {
480  for(auto it: add){
481  const std::string& key = it.first;
482 
483  // Needs special handling anyway, leave it blank.
484  if(key == "parents") continue;
485 
486  // Accumulate the runs list
487  if(key == "runs"){
488  const std::string& r1 = base[key];
489  const std::string& r2 = it.second;
490 
491  assert(!r2.empty());
492 
493  // when concatenating MC files,
494  // it's possible to get duplicates
495  // (run+subrun the same, only differing
496  // in cycle number, which isn't part
497  // of this metadata field)
498  // don't repeat, because SAM gets angry.
499  if (r1.find(r2.substr(1, r2.size() - 2)) != std::string::npos)
500  continue;
501 
502  if(r1.empty()){
503  base[key] = r2;
504  continue;
505  }
506 
507  // "[foo]" + "[bar]"
508  std::string sum = r1+&r2[1]; // "[foo]bar]"
509  sum[r1.size()-1] = ','; // "[foo,bar]"
510  base[key] = sum;
511  continue;
512  }
513 
514  if(base.find(key) == base.end()){
515  // If it's new, add it
516  base[key] = it.second;
517  }
518  else{
519  if(key == "simulated.number_of_spills" ||
520  key == "event_count" ||
521  key == "online.totalevents"){
522  // These two fields should be accumulated
523  base[key] = TString::Format("%d",
524  atoi(base[key].c_str()) +
525  atoi(it.second.c_str())).Data();
526  }
527  else{
528  // If it's a clash, record it
529  if(base[key] != it.second) mask.insert(key);
530  }
531  }
532  }
533  }
534 
535 
536  //----------------------------------------------------------------------
537  void WriteCAFMetadata(TDirectory* dir,
538  const std::map<std::string, std::string>& meta)
539  {
540  TDirectory* tmp = gDirectory;
541  dir->cd();
542 
543  for(auto it: meta){
544  TObjString str(it.second.c_str());
545  str.Write(it.first.c_str());
546  }
547 
548  dir->Save();
549 
550  tmp->cd();
551  }
552 
553  //----------------------------------------------------------------------
554 
556  {
557  // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
558  // CAF files and fitting the resulting distribution to a gaussian.
559  // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
560  // and containment cuts.
561  // CAF Files used:
562  // NearDet:
563  // prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
564  // FarDet:
565  // prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08
566  // Unknown for NDOS
567 
568  switch(det){
569  case caf::kFARDET: return TVector3(-6.833e-05,
570  +6.388e-02,
571  +9.980e-1).Unit();
572  case caf::kNDOS: return TVector3(-3.845e-4,
573  +6.517e-2,
574  +9.967e-1).Unit();
575  case caf::kNEARDET: return TVector3(-8.424e-04,
576  -6.174e-02,
577  +9.981e-1).Unit();
578  default:
579  std::cerr<< "Aborting. Bad geometry configuration: " << det << "\n";
580  abort();
581  }
582 
583  }
584 
585 
586  //----------------------------------------------------------------------
588  {
589  return (getenv("_CONDOR_SCRATCH_DIR") != 0);
590  }
591 
592 
593  //----------------------------------------------------------------------
594  size_t JobNumber()
595  {
596  if (!RunningOnGrid())
597  throw std::runtime_error("Not running on grid!");
598 
599  return std::stoi(std::string(getenv("PROCESS")));
600  }
601 
602  //----------------------------------------------------------------------
604  {
605  // I would be much much happier to do this in proper code, but I'm not sure
606  // how, there's no samweb C++ API?
607  return system(TString::Format("samweb list-definitions --defname %s | grep %s", def.c_str(), def.c_str()).Data()) == 0;
608  }
609 
610  //----------------------------------------------------------------------
611  bool AlmostEqual(double a, double b, double eps)
612  {
613  if(a == 0 && b == 0) return true;
614 
615  return fabs(a-b) <= eps * std::max(fabs(a), fabs(b));
616  }
617 
618 
619  //----------------------------------------------------------------------
621  {
622  static bool first = true;
623  static bool onsite = false;
624 
625  if (first && unauth) {
626  first = false;
627  char chostname[255];
628  gethostname(chostname, 255);
629  std::string hostname = chostname;
630 
631  if ( hostname.find("fnal.gov") != std::string::npos ){
632  onsite = true;
633  std::cout << "Using unauthenticated xrootd access (port 1095) while on-site, hostname: " << hostname << std::endl;
634  }
635  else {
636  onsite = false;
637  std::cout << "Using authenticated xrootd access (port 1094) access while off-site, hostname: " << hostname << std::endl;
638  }
639  }
640 
641  if(loc.rfind("/pnfs/", 0) == 0){ // ie begins with
642  if ( onsite && unauth )
643  loc = std::string("root://fndca1.fnal.gov:1095//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
644  else
645  loc = std::string("root://fndca1.fnal.gov:1094//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
646  }
647  return loc;
648  }
649 
650  //----------------------------------------------------------------------
651  std::size_t Hash(const std::vector<std::size_t> &vals)
652  {
653  std::size_t seed = 0;
654 
655  for (auto v : vals)
656  seed ^= std::hash<std::size_t>{}(v) + 0x9e3779b9 + (seed<<6) + (seed>>2);
657  return seed;
658  }
659 
660  //----------------------------------------------------------------------
661  std::size_t Hash(const caf::SRSpillProxy *sr)
662  {
663  return Hash({sr->run, sr->subrun, sr->evt});
664  }
665 
666  //----------------------------------------------------------------------
667  std::size_t Hash(const caf::SRProxy *sr)
668  {
669  return Hash({Hash(&(sr->spill)), std::size_t(sr->hdr.cycle), sr->hdr.subevt});
670  }
671 
672  //----------------------------------------------------------------------
673  FitToFourier::FitToFourier(TH1* h, double xlo, double xhi, int NOsc)
674  : fHist(h), fxlo(xlo), fxhi(xhi), fNOsc(NOsc)
675  {
676  }
677 
678  //----------------------------------------------------------------------
680  {
681  }
682 
683  //----------------------------------------------------------------------
684  double FitToFourier::operator()(double *x, double *par) const
685  {
686  double x0 = x[0];
687  double val = par[0];
688  for (int i = 1; i <= fNOsc; i++)
689  val += par[2*i-1]*sin(i*M_PI*x0) + par[2*i]*cos(i*M_PI*x0);
690  return val;
691  }
692 
693  //----------------------------------------------------------------------
694  TF1* FitToFourier::Fit() const
695  {
696  std::vector<double> s(fNOsc, 0.);
697  std::vector<double> c(fNOsc, 0.);
698  int nBins = 0;
699  for(int i = 1; i <= fHist->GetNbinsX(); ++i){
700  const double x = M_PI * fHist->GetXaxis()->GetBinCenter(i);
701  const double y = fHist->GetBinContent(i);
702 
703  if(y == 0) continue;
704  ++nBins;
705 
706  for(int n = 0; n <= fNOsc; ++n){
707  s[n] += y * sin(n*x);
708  c[n] += y * cos(n*x);
709  }
710  }
711 
712  for(int n = 0; n <= fNOsc; ++n){
713  s[n] *= 2./nBins;
714  c[n] *= 2./nBins;
715  }
716 
717  TF1* f = new TF1(UniqueName().c_str(), this, fxlo, fxhi, 2*fNOsc+1);
718 
719  f->SetParameter(0, c[0]/2);
720  for(int n = 1; n <= fNOsc; ++n){
721  f->SetParameter(n*2-1, s[n]);
722  f->SetParameter(n*2, c[n]);
723  }
724 
725  // Because ROOT is having problems drawing f if I don't
726  double min = fHist->GetMinimum();
727  double max = fHist->GetMaximum();
728  f->GetYaxis()->SetRangeUser(0.8*min, 1.2*max);
729  return f;
730  }
731 }
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2055
Near Detector underground.
Definition: SREnums.h:10
const double fxhi
Definition: Utilities.h:309
T max(const caf::Proxy< T > &a, T b)
Det_t
Which NOvA detector?
Definition: SREnums.h:7
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:594
Far Detector at Ash River.
Definition: SREnums.h:11
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
system("rm -rf microbeam.root")
int status
Definition: fabricate.py:1613
set< int >::iterator it
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:68
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:209
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:587
TH1F * MakeTH1F(const char *name, const char *title, const Binning &bins)
Definition: Utilities.cxx:221
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2049
const char * p
Definition: xmltok.h:285
TH3 * ToTH3Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Helper for ana::ToTH3.
Definition: Utilities.cxx:327
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
Definition: Utilities.cxx:215
TH2F * MakeTH2F(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: Utilities.cxx:235
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:476
Divide bin contents by bin widths.
Definition: Utilities.h:45
OStream cerr
Definition: OStream.cxx:7
FitToFourier(TH1 *h, double xlo, double xhi, int NOsc)
Definition: Utilities.cxx:673
Float_t tmp
Definition: plot.C:36
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:620
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:429
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
void MakeHist(int n)
#define M_PI
Definition: SbMath.h:34
double operator()(double *x, double *par) const
Definition: Utilities.cxx:684
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
FloatingExceptionOnNaN(bool enable=true)
Definition: Utilities.cxx:65
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:154
list fileList
Definition: cafExposure.py:30
bool IsSimple() const
Definition: Binning.h:29
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:435
EBinType
Definition: Utilities.h:42
const XML_Char * s
Definition: expat.h:262
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
Definition: Utilities.cxx:381
double Min() const
Definition: Binning.h:26
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:245
unsigned int seed
Definition: runWimpSim.h:102
bool AlmostEqual(double a, double b, double eps)
Definition: Utilities.cxx:611
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:90
const double fxlo
Definition: Utilities.h:308
const double a
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: Utilities.h:50
std::string getenv(std::string const &name)
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:537
caf::Proxy< int > cycle
Definition: SRProxy.h:225
const TH1 * fHist
Definition: Utilities.h:307
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:409
TF1 * Fit() const
Definition: Utilities.cxx:694
Represent the ratio between two spectra.
Definition: Ratio.h:8
const std::vector< double > & Edges() const
Definition: Binning.h:30
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Utility function to avoid 4-way combinatorial explosion on the bin types.
Definition: Utilities.cxx:227
double Max() const
Definition: Binning.h:27
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:651
static const double A
Definition: Units.h:82
T sin(T number)
Definition: d0nt_math.hpp:132
pub
Definition: cvnie.py:6
bool SAMDefinitionExists(const std::string &def)
Definition: Utilities.cxx:603
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
int NBins() const
Definition: Binning.h:25
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Extract map of metadata parameters from a CAF file.
Definition: Utilities.cxx:454
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
TH2 * ToTH2Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Helper for ana::ToTH2.
Definition: Utilities.cxx:276
const int fNOsc
Definition: Utilities.h:310
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:555
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
Definition: Xdiff_gwt.C:5
TH3D * MakeTH3D(const char *name, const char *title, const Binning &binsx, const Binning &binsy, const Binning &binsz)
Utility function for 3D hist.
Definition: Utilities.cxx:244
Double_t sum
Definition: plot.C:31
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
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
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: Utilities.cxx:303
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: Utilities.cxx:256
static constexpr Double_t sr
Definition: Munits.h:164