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  //----------------------------------------------------------------------
154  double LogLikelihoodDerivative(double e, double o, double dedx)
155  {
156  double ret = 2*dedx;
157  if(o) ret -= 2*o*dedx / e;
158  return ret;
159  }
160 
161  //----------------------------------------------------------------------
162  double LogLikelihoodDerivative(const TH1D* eh, const TH1D* oh,
163  const std::vector<double>& dedx)
164  {
165  const double* ea = eh->GetArray();
166  const double* oa = oh->GetArray();
167 
168  double ret = 0;
169  for(unsigned int i = 0; i < dedx.size(); ++i){
170  ret += LogLikelihoodDerivative(ea[i], oa[i], dedx[i]);
171  }
172  return ret;
173  }
174 
175  //----------------------------------------------------------------------
176  double Chi2CovMx(const TVectorD* e, const TVectorD* o, const TMatrixD* covmxinv)
177  {
178  assert (e->GetNrows() == o->GetNrows());
179 
180  TVectorD diff = *o - *e;
181  return diff * ((*covmxinv) * diff); // operator* for two TVectorDs is the "dot product" (i.e., v1 * v2 = v1^{trans}v1)
182  }
183 
184  //----------------------------------------------------------------------
185  double Chi2CovMx(const TH1* e, const TH1* o, const TMatrixD* covmxinv)
186  {
187  TVectorD eVec(e->GetNbinsX());
188  TVectorD oVec(o->GetNbinsX());
189  for (int bin = 1; bin <= e->GetNbinsX(); bin++)
190  eVec[bin-1] = e->GetBinContent(bin);
191  for (int bin = 1; bin <= o->GetNbinsX(); bin++)
192  oVec[bin-1] = o->GetBinContent(bin);
193 
194  return Chi2CovMx(&eVec, &oVec, covmxinv);
195  }
196 
197  //----------------------------------------------------------------------
199  int nbinsx, double xmin, double xmax,
200  int nbinsy, double ymin, double ymax)
201  {
202  DontAddDirectory guard;
203 
204  // How wide the bins will be once we're done
205  const double xwidth = (xmax-xmin)/(nbinsx-1);
206  const double ywidth = (ymax-ymin)/(nbinsy-1);
207 
208  // Move the bin edges so that the limits occur at the centres
209  xmin -= xwidth/2; ymin -= ywidth/2;
210  xmax += xwidth/2; ymax += ywidth/2;
211 
212  return new TH2F(UniqueName().c_str(), title.c_str(),
213  nbinsx, xmin, xmax,
214  nbinsy, ymin, ymax);
215  }
216 
217 
218  //----------------------------------------------------------------------
219  std::unique_ptr<TMatrixD> SymmMxInverse(const TMatrixD& mx)
220  {
221  // check if there are any null rows/columns.
222  // if there are, they make the matrix singular.
223  // we will remove them temporarily,
224  // invert the matrix, then put them back afterwards.
225  std::set<int> nullRows;
226  for (auto row = mx.GetRowLwb(); row <= mx.GetRowUpb(); row++)
227  {
228  bool rowIsNull = true;
229  for (auto col = mx.GetColLwb(); col <= mx.GetColUpb(); col++)
230  {
231  if (mx[row][col] != 0.)
232  {
233  rowIsNull = false;
234  break;
235  }
236  }
237 
238  if (rowIsNull)
239  nullRows.insert(row);
240  }
241 
242  std::cerr << " Notice: covariance matrix '" << mx.GetName() << "' has " << nullRows.size() << " null rows.\n"
243  << "They will be removed before inverting and added back afterwards." << std::endl;
244 
245  // create a new matrix for inverting, skipping the null rows
246  auto invMx = std::make_unique<TMatrixD>(mx.GetRowLwb(), mx.GetRowUpb() - nullRows.size(),
247  mx.GetColLwb(), mx.GetColUpb() - nullRows.size());
248  unsigned int skippedRows = 0;
249  for (auto row = mx.GetRowLwb(); row <= mx.GetRowUpb(); row++)
250  {
251  if (nullRows.find(row) != nullRows.end())
252  {
253  skippedRows++;
254  continue;
255  }
256  unsigned int skippedCols = 0;
257  for (auto col = mx.GetColLwb(); col <= mx.GetColUpb(); col++)
258  {
259  // since we assumed the matrix is symmetric,
260  // we can just use the null rows list here
261  if (nullRows.find(col) != nullRows.end())
262  {
263  skippedCols++;
264  continue;
265  }
266 
267  (*invMx)[col-skippedCols][row-skippedRows] = (*invMx)[row-skippedRows][col-skippedCols] = mx[row][col];
268  }
269  }
270 
271  invMx->Invert();
272 
273  // put back the empty rows if there were any
274  if (nullRows.size())
275  {
276  skippedRows = 0;
277  auto retMx = std::make_unique<TMatrixD>(mx.GetRowLwb(), mx.GetRowUpb(),
278  mx.GetColLwb(), mx.GetColUpb());
279  for (auto row = mx.GetRowLwb(); row <= mx.GetRowUpb(); row++)
280  {
281  if (nullRows.find(row) != nullRows.end())
282  {
283  skippedRows++;
284  continue;
285  }
286 
287  unsigned int skippedCols = skippedRows;
288  for (auto col = row; col <= mx.GetColUpb(); col++)
289  {
290  if (nullRows.find(col) != nullRows.end())
291  {
292  skippedCols++;
293  continue;
294  }
295 
296  (*retMx)[col][row] = (*retMx)[row][col] = (*invMx)[row-skippedRows][col-skippedCols];
297  }
298  }
299 
300  return retMx;
301  }
302 
303  return invMx;
304  }
305 
306  // Helper functions for MakeTHND().
307  namespace{
308  // Eventually the bin parameters will all be unpacked and we just pass them
309  // on to the regular constructor.
310  template<class T, class... A> T* MakeHist(A... a)
311  {
312  DontAddDirectory guard;
313  return new T(a...);
314  }
315 
316  // This function consumes bins from the start of the argument list and
317  // pushes their translations onto the list of arguments at the end.
318  template<class T, class... A> T* MakeHist(const Binning& firstBin,
319  A... args)
320  {
321  if(firstBin.IsSimple())
322  return MakeHist<T>(args...,
323  firstBin.NBins(), firstBin.Min(), firstBin.Max());
324  else
325  return MakeHist<T>(args...,
326  firstBin.NBins(), &firstBin.Edges().front());
327  }
328  }
329 
330  // Concrete instantiations. MakeHist() requires us to put the bin arguments
331  // first...
332  //----------------------------------------------------------------------
333  TH1D* MakeTH1D(const char* name, const char* title, const Binning& bins)
334  {
335  return MakeHist<TH1D>(bins, name, title);
336  }
337 
338  //----------------------------------------------------------------------
339  TH1F* MakeTH1F(const char* name, const char* title, const Binning& bins)
340  {
341  return MakeHist<TH1F>(bins, name, title);
342  }
343 
344  //----------------------------------------------------------------------
345  TH2D* MakeTH2D(const char* name, const char* title,
346  const Binning& binsx,
347  const Binning& binsy)
348  {
349  return MakeHist<TH2D>(binsx, binsy, name, title);
350  }
351 
352  //----------------------------------------------------------------------
353  TH2F* MakeTH2F(const char* name, const char* title,
354  const Binning& binsx,
355  const Binning& binsy)
356  {
357  return MakeHist<TH2F>(binsx, binsy, name, title);
358  }
359 
360  //----------------------------------------------------------------------
361 
362  TH3D* MakeTH3D(const char* name, const char* title,
363  const Binning& binsx,
364  const Binning& binsy,
365  const Binning& binsz)
366  {
367  return MakeHist<TH3D>(name, title,
368  binsx.NBins(), &binsx.Edges().front(),
369  binsy.NBins(), &binsy.Edges().front(),
370  binsz.NBins(), &binsz.Edges().front());
371  }
372 
373  //----------------------------------------------------------------------
374  TH2* ToTH2(const Spectrum& s, double exposure, ana::EExposureType expotype,
375  const Binning& binsx, const Binning& binsy, ana::EBinType bintype)
376  {
377  DontAddDirectory guard;
378 
379  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
380  return ToTH2Helper(std::move(h1), binsx, binsy, bintype);
381 
382  }
383 
384  //----------------------------------------------------------------------
385  TH2* ToTH2(const Ratio& r,
386  const Binning& binsx, const Binning& binsy)
387  {
388  DontAddDirectory guard;
389 
390  std::unique_ptr<TH1> h1(r.ToTH1());
391  return ToTH2Helper(std::move(h1), binsx, binsy);
392 
393  }
394 
395  //----------------------------------------------------------------------
396 
397  TH2* ToTH2Helper(std::unique_ptr<TH1> h1,
398  const Binning& binsx, const Binning& binsy,
399  ana::EBinType bintype)
400  {
401  // Make sure it's compatible with having been made with this binning
402  assert(h1->GetNbinsX() == binsx.NBins()*binsy.NBins());
403 
404  TH2* ret = MakeTH2D("", UniqueName().c_str(), binsx, binsy);
405 
406  for(int i = 0; i < h1->GetNbinsX(); ++i){
407  const double val = h1->GetBinContent(i+1);
408  const double err = h1->GetBinError(i+1);
409 
410  const int ix = i / binsy.NBins();
411  const int iy = i % binsy.NBins();
412 
413  ret->SetBinContent(ix+1, iy+1, val);
414  ret->SetBinError (ix+1, iy+1, err);
415  }
416 
417  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
418 
419  return ret;
420  }
421 
422  //----------------------------------------------------------------------
423 
424  TH3* ToTH3(const Spectrum& s, double exposure, ana::EExposureType expotype,
425  const Binning& binsx, const Binning& binsy, const Binning& binsz,
426  ana::EBinType bintype)
427  {
428  DontAddDirectory guard;
429 
430  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
431 
432  return ToTH3Helper(std::move(h1), binsx, binsy, binsz, bintype);
433  }
434 
435  //----------------------------------------------------------------------
436 
437  TH3* ToTH3(const Ratio& r,
438  const Binning& binsx, const Binning& binsy, const Binning& binsz)
439  {
440  DontAddDirectory guard;
441 
442  std::unique_ptr<TH1> h1(r.ToTH1());
443 
444  return ToTH3Helper(std::move(h1), binsx, binsy, binsz);
445  }
446 
447  //----------------------------------------------------------------------
448  TH3* ToTH3Helper(std::unique_ptr<TH1> h1,
449  const Binning& binsx,
450  const Binning& binsy,
451  const Binning& binsz,
452  ana::EBinType bintype)
453  {
454 
455  const int nx = binsx.NBins();
456  const int ny = binsy.NBins();
457  const int nz = binsz.NBins();
458 
459  // Make sure it's compatible with having been made with this binning
460  assert(h1->GetNbinsX() == nx*ny*nz);
461 
462  TH3* ret;
463 
464  // If all three axes are simple, we can call a simpler constructor
465  if(binsx.IsSimple() && binsy.IsSimple() && binsz.IsSimple()){
466  ret = new TH3F(UniqueName().c_str(), "",
467  nx, binsx.Min(), binsx.Max(),
468  ny, binsy.Min(), binsy.Max(),
469  nz, binsz.Min(), binsz.Max());
470 
471  }
472  else{
473  // If >= 1 axis has a not simple binning, you have to pass the
474  // edges vector to the TH3 constructor
475  ret = new TH3F(UniqueName().c_str(), "",
476  nx, &binsx.Edges().front(),
477  ny, &binsy.Edges().front(),
478  nz, &binsz.Edges().front());
479  }
480 
481  for(int i = 0; i < h1->GetNbinsX(); ++i){
482  const double val = h1->GetBinContent(i+1);
483  const double err = h1->GetBinError(i+1);
484 
485  const int nynz = ny*nz;
486  const int nmodnynz = i%nynz;
487  const int ix = i/nynz;
488  const int iy = nmodnynz/nz;
489  const int iz = i%nz;
490 
491  ret->SetBinContent(ix+1, iy+1, iz+1, val);
492  ret->SetBinError (ix+1, iy+1, iz+1, err);
493  }
494 
495  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
496 
497  return ret;
498 
499  }
500 
501  //----------------------------------------------------------------------
502  std::vector<std::string> Wildcard(const std::string& wildcardString)
503  {
504  // Expand environment variables and wildcards like the shell would
505  wordexp_t p;
506  const int status = wordexp(wildcardString.c_str(), &p, WRDE_SHOWERR);
507 
508  if(status != 0){
509  std::cerr << "Wildcard string '" << wildcardString
510  << "' returned error " << status << " from wordexp()."
511  << std::endl;
512  return {};
513  }
514 
515  std::vector<std::string> fileList;
516 
517  for(unsigned int i = 0; i < p.we_wordc; ++i){
518  // Check the file exists before adding it
519  struct stat sb;
520  if(stat(p.we_wordv[i], &sb) == 0)
521  fileList.push_back(p.we_wordv[i]);
522  }
523 
524  wordfree(&p);
525 
526  return fileList;
527  }
528 
529  //----------------------------------------------------------------------
531  {
532 #ifdef NOVACMAKE
533  const char* pub = getenv("NOVASOFT_DIR");
534  const char* priv = getenv("NOVASOFT_DIR");
535 #else
536  const char* pub = getenv("FW_BASE");
537  const char* priv = getenv("FW_RELEASE_BASE");
538 #endif
539 
540  if(priv && priv != std::string(".")){
541  const std::string ret = std::string(priv)+"/"+Dir;
542  struct stat junk;
543  if(stat(ret.c_str(), &junk) == 0) return ret;
544  }
545  assert(pub);
546  return std::string(pub)+"/"+Dir;
547  }
548 
549  //----------------------------------------------------------------------
551  {
552  return FindPackageDir("CAFAna");
553  }
554 
555  //----------------------------------------------------------------------
556  std::vector<std::string> LoadFileList(const std::string& listfile)
557  {
558  std::vector<std::string> ret;
559 
560  std::ifstream is(listfile);
561  if(!is.good()){
562  std::cerr << "Can't open file list '" << listfile << "'. Aborting." << std::endl;
563  abort();
564  }
565 
566  while(!is.eof()){
568  is >> fname;
569  if(!fname.empty()) ret.push_back(fname);
570  }
571  return ret;
572  }
573 
574  //----------------------------------------------------------------------
575  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir)
576  {
577  std::map<std::string, std::string> ret;
578 
579  TIter next(dir->GetListOfKeys());
580  while(TObject* key = next()){
581  TObject* obj = dir->Get(key->GetName());
582  assert(obj);
583 
584  const TString className = obj->ClassName();
585  if(className == "TObjString"){
586  ret[key->GetName()] = ((TObjString*)obj)->GetString();
587  }
588  else{
589  std::cerr << "Unexpected object " << key->GetName() << " of type " << className << " while looking for metadata. Ignoring" << std::endl;
590  }
591  }
592 
593  return ret;
594  }
595 
596  //----------------------------------------------------------------------
597  void CombineMetadata(std::map<std::string, std::string>& base,
598  const std::map<std::string, std::string>& add,
599  std::set<std::string>& mask)
600  {
601  for(auto it: add){
602  const std::string& key = it.first;
603 
604  // Needs special handling anyway, leave it blank.
605  if(key == "parents") continue;
606 
607  // Accumulate the runs list
608  if(key == "runs"){
609  const std::string& r1 = base[key];
610  const std::string& r2 = it.second;
611 
612  assert(!r2.empty());
613 
614  // when concatenating MC files,
615  // it's possible to get duplicates
616  // (run+subrun the same, only differing
617  // in cycle number, which isn't part
618  // of this metadata field)
619  // don't repeat, because SAM gets angry.
620  if (r1.find(r2.substr(1, r2.size() - 2)) != std::string::npos)
621  continue;
622 
623  if(r1.empty()){
624  base[key] = r2;
625  continue;
626  }
627 
628  // "[foo]" + "[bar]"
629  std::string sum = r1+&r2[1]; // "[foo]bar]"
630  sum[r1.size()-1] = ','; // "[foo,bar]"
631  base[key] = sum;
632  continue;
633  }
634 
635  if(base.find(key) == base.end()){
636  // If it's new, add it
637  base[key] = it.second;
638  }
639  else{
640  if(key == "simulated.number_of_spills" ||
641  key == "event_count" ||
642  key == "online.totalevents"){
643  // These two fields should be accumulated
644  base[key] = TString::Format("%d",
645  atoi(base[key].c_str()) +
646  atoi(it.second.c_str())).Data();
647  }
648  else{
649  // If it's a clash, record it
650  if(base[key] != it.second) mask.insert(key);
651  }
652  }
653  }
654  }
655 
656 
657  //----------------------------------------------------------------------
658  void WriteCAFMetadata(TDirectory* dir,
659  const std::map<std::string, std::string>& meta)
660  {
661  TDirectory* tmp = gDirectory;
662  dir->cd();
663 
664  for(auto it: meta){
665  TObjString str(it.second.c_str());
666  str.Write(it.first.c_str());
667  }
668 
669  dir->Save();
670 
671  tmp->cd();
672  }
673 
674  //----------------------------------------------------------------------
675 
677  {
678  // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
679  // CAF files and fitting the resulting distribution to a gaussian.
680  // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
681  // and containment cuts.
682  // CAF Files used:
683  // NearDet:
684  // prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
685  // FarDet:
686  // prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08
687  // Unknown for NDOS
688 
689  switch(det){
690  case caf::kFARDET: return TVector3(-6.833e-05,
691  +6.388e-02,
692  +9.980e-1).Unit();
693  case caf::kNDOS: return TVector3(-3.845e-4,
694  +6.517e-2,
695  +9.967e-1).Unit();
696  case caf::kNEARDET: return TVector3(-8.424e-04,
697  -6.174e-02,
698  +9.981e-1).Unit();
699  default:
700  std::cerr<< "Aborting. Bad geometry configuration: " << det << "\n";
701  abort();
702  }
703 
704  }
705 
706 
707  //----------------------------------------------------------------------
709  {
710  return (getenv("_CONDOR_SCRATCH_DIR") != 0);
711  }
712 
713 
714  //----------------------------------------------------------------------
715  size_t JobNumber()
716  {
717  if (!RunningOnGrid())
718  throw std::runtime_error("Not running on grid!");
719 
720  return std::stoi(std::string(getenv("PROCESS")));
721  }
722 
723  //----------------------------------------------------------------------
725  {
726  // I would be much much happier to do this in proper code, but I'm not sure
727  // how, there's no samweb C++ API?
728  return system(TString::Format("samweb list-definitions --defname %s | grep %s", def.c_str(), def.c_str()).Data()) == 0;
729  }
730 
731  //----------------------------------------------------------------------
732  bool AlmostEqual(double a, double b)
733  {
734  if(a == 0 && b == 0) return true;
735 
736  return fabs(a-b)/std::max(a, b) < .0001; // allow 0.01% error
737  }
738 
739 
740  //----------------------------------------------------------------------
742  {
743  static bool first = true;
744  static bool onsite = false;
745 
746  if (first && unauth) {
747  first = false;
748  char chostname[255];
749  gethostname(chostname, 255);
750  std::string hostname = chostname;
751 
752  if ( hostname.find("fnal.gov") != std::string::npos ){
753  onsite = true;
754  std::cout << "Using unauthenticated xrootd access (port 1095) while on-site, hostname: " << hostname << std::endl;
755  }
756  else {
757  onsite = false;
758  std::cout << "Using authenticated xrootd access (port 1094) access while off-site, hostname: " << hostname << std::endl;
759  }
760  }
761 
762  if(loc.rfind("/pnfs/", 0) == 0){ // ie begins with
763  if ( onsite && unauth )
764  loc = std::string("root://fndca1.fnal.gov:1095//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
765  else
766  loc = std::string("root://fndca1.fnal.gov:1094//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
767  }
768  return loc;
769  }
770 
771  //----------------------------------------------------------------------
772  std::size_t Hash(const std::vector<std::size_t> &vals)
773  {
774  std::size_t seed = 0;
775 
776  for (auto v : vals)
777  seed ^= std::hash<std::size_t>{}(v) + 0x9e3779b9 + (seed<<6) + (seed>>2);
778  return seed;
779  }
780 
781  //----------------------------------------------------------------------
782  std::size_t Hash(const caf::SRSpillProxy *sr)
783  {
784  return Hash({sr->run, sr->subrun, sr->evt});
785  }
786 
787  //----------------------------------------------------------------------
788  std::size_t Hash(const caf::SRProxy *sr)
789  {
790  return Hash({Hash(&(sr->spill)), std::size_t(sr->hdr.cycle), sr->hdr.subevt});
791  }
792 
793  //----------------------------------------------------------------------
794  FitToFourier::FitToFourier(TH1* h, double xlo, double xhi, int NOsc)
795  : fHist(h), fxlo(xlo), fxhi(xhi), fNOsc(NOsc)
796  {
797  }
798 
799  //----------------------------------------------------------------------
800  FitToFourier::~FitToFourier()
801  {
802  }
803 
804  //----------------------------------------------------------------------
805  double FitToFourier::operator()(double *x, double *par) const
806  {
807  double x0 = x[0];
808  double val = par[0];
809  for (int i = 1; i <= fNOsc; i++)
810  val += par[2*i-1]*sin(i*M_PI*x0) + par[2*i]*cos(i*M_PI*x0);
811  return val;
812  }
813 
814  //----------------------------------------------------------------------
815  TF1* FitToFourier::Fit() const
816  {
817  std::vector<double> s(fNOsc, 0.);
818  std::vector<double> c(fNOsc, 0.);
819  int nBins = 0;
820  for(int i = 1; i <= fHist->GetNbinsX(); ++i){
821  const double x = M_PI * fHist->GetXaxis()->GetBinCenter(i);
822  const double y = fHist->GetBinContent(i);
823 
824  if(y == 0) continue;
825  ++nBins;
826 
827  for(int n = 0; n <= fNOsc; ++n){
828  s[n] += y * sin(n*x);
829  c[n] += y * cos(n*x);
830  }
831  }
832 
833  for(int n = 0; n <= fNOsc; ++n){
834  s[n] *= 2./nBins;
835  c[n] *= 2./nBins;
836  }
837 
838  TF1* f = new TF1(UniqueName().c_str(), this, fxlo, fxhi, 2*fNOsc+1);
839 
840  f->SetParameter(0, c[0]/2);
841  for(int n = 1; n <= fNOsc; ++n){
842  f->SetParameter(n*2-1, s[n]);
843  f->SetParameter(n*2, c[n]);
844  }
845 
846  // Because ROOT is having problems drawing f if I don't
847  double min = fHist->GetMinimum();
848  double max = fHist->GetMaximum();
849  f->GetYaxis()->SetRangeUser(0.8*min, 1.2*max);
850  return f;
851  }
852 }
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: Utilities.cxx:345
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:1986
Near Detector underground.
Definition: SREnums.h:10
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()
Definition: Utilities.cxx:715
TH3D * MakeTH3D(const char *name, const char *title, const Binning &binsx, const Binning &binsy, const Binning &binsz)
Definition: Utilities.cxx:362
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: CutFlow_header.h:5
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:122
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:564
bool RunningOnGrid()
Definition: Utilities.cxx:708
TH3 * ToTH3Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Definition: Utilities.cxx:448
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:1980
fVtxDx Fit("f")
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1969
TH2 * ToTH2Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:397
void CombineMetadata(std::map< std::string, std::string > &base, const std::map< std::string, std::string > &add, std::set< std::string > &mask)
Definition: Utilities.cxx:597
Divide bin contents by bin widths.
Definition: Utilities.h:45
OStream cerr
Definition: OStream.cxx:7
std::vector< std::string > Wildcard(const std::string &wildcardString)
Definition: Utilities.cxx:502
std::unique_ptr< TMatrixD > SymmMxInverse(const TMatrixD &mx)
Definition: Utilities.cxx:219
Float_t tmp
Definition: plot.C:36
stan::math::var LogLikelihood(const std::vector< stan::math::var > &exp, const TH1 *obs)
Variant that handles the prediction in the form of Stan vars.
Definition: StanUtils.cxx:11
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:741
std::string FindCAFAnaDir()
Definition: Utilities.cxx:550
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
FloatingExceptionOnNaN(bool enable=true)
Definition: Utilities.cxx:65
list fileList
Definition: cafExposure.py:30
bool IsSimple() const
Definition: Binning.h:29
Double_t ymax
Definition: plot.C:25
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Definition: Utilities.cxx:575
EBinType
Definition: Utilities.h:42
double LogLikelihoodDerivative(double e, double o, double dedx)
Definition: Utilities.cxx:154
const XML_Char * s
Definition: expat.h:262
double Min() const
Definition: Binning.h:26
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:244
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:90
Int_t col[ntarg]
Definition: Style.C:29
const double a
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: Utilities.h:50
std::string getenv(std::string const &name)
Class for parsing *.fcl files for their metadata information.
void WriteCAFMetadata(TDirectory *dir, const std::map< std::string, std::string > &meta)
Definition: Utilities.cxx:658
bool AlmostEqual(double a, double b)
Definition: Utilities.cxx:732
caf::Proxy< int > cycle
Definition: SRProxy.h:224
float bin[41]
Definition: plottest35.C:14
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:530
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:374
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
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Definition: Utilities.cxx:333
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double Max() const
Definition: Binning.h:27
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH2F * MakeTH2F(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: Utilities.cxx:353
TH1F * h1
UInt_t par
Definition: AnaPlotMaker.h:47
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::size_t Hash(const std::vector< std::size_t > &vals)
Definition: Utilities.cxx:772
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:724
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
int NBins() const
Definition: Binning.h:25
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
TVector3 NuMIBeamDirection(caf::Det_t det)
Definition: Utilities.cxx:676
Double_t ymin
Definition: plot.C:24
std::vector< std::string > LoadFileList(const std::string &listfile)
Definition: Utilities.cxx:556
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
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
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Definition: Utilities.cxx:424
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:30
TH1F * MakeTH1F(const char *name, const char *title, const Binning &bins)
Definition: Utilities.cxx:339
double Chi2CovMx(const TVectorD *e, const TVectorD *o, const TMatrixD *covmxinv)
Definition: Utilities.cxx:176
static constexpr Double_t sr
Definition: Munits.h:164
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax)
Definition: Utilities.cxx:198