10 #include "TDirectory.h" 15 #include "TObjString.h" 38 fBackup = TH1::AddDirectoryStatus();
39 TH1::AddDirectory(
false);
51 const char*
s =
getenv(
"IFDH_SILENT");
53 if(!fSet) setenv(
"IFDH_SILENT",
"1", 1);
60 if(!fSet) unsetenv(
"IFDH_SILENT");
68 feclearexcept(FE_INVALID);
70 fegetexceptflag(&
fBackup, FE_INVALID);
73 feenableexcept(FE_INVALID);
75 fedisableexcept(FE_INVALID);
81 fesetexceptflag(&
fBackup, FE_INVALID);
85 std::unique_ptr<TMatrixD>
CalcCovMx(
const std::vector<TArrayD*> & binSets,
int firstBin,
int lastBin)
87 if (binSets.size() < 2)
88 return std::unique_ptr<TMatrixD>(
nullptr);
91 lastBin = binSets[0]->GetSize() - 1;
93 int nBins = lastBin - firstBin + 1;
95 std::vector<double> binMeans(nBins);
96 for(
const auto & binSet : binSets )
98 for ( decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++ )
99 binMeans[binIdx] += (*binSet)[binIdx];
101 for (decltype(lastBin) binIdx = firstBin; binIdx <= lastBin; binIdx++)
102 binMeans[binIdx] /= binSets.size();
105 auto covmx = std::make_unique<TMatrixD>(
nBins,
nBins);
107 for(
unsigned int hist_idx = 0; hist_idx < binSets.size(); ++hist_idx )
110 for( decltype(nBins)
i = 0;
i <
nBins;
i++ )
112 double xi = (*(binSets[hist_idx]))[
i];
113 for( decltype(nBins) k =
i; k <
nBins; k++ )
115 double xk = (*(binSets[hist_idx]))[k];
116 (*covmx)[
i][k] += (xi - binMeans[
i]) * (xk - binMeans[k]);
118 (*covmx)[k][
i] = (*covmx)[
i][k];
124 (*covmx) *= 1./(binSets.size()-1);
132 assert(eh->GetNbinsX() == oh->GetNbinsX());
136 int bufferBins = useOverflow? 2 : 1;
138 for(
int i = 0;
i < eh->GetNbinsX()+bufferBins; ++
i){
139 const double e = eh->GetBinContent(
i);
140 const double o = oh->GetBinContent(
i);
151 int nbinsy,
double ymin,
double ymax,
bool ylog)
155 if(xlog){xmin =
log(xmin); xmax =
log(xmax);}
156 if(ylog){ymin =
log(ymin); ymax =
log(ymax);}
159 const double xwidth = (xmax-
xmin)/(nbinsx-1);
160 const double ywidth = (ymax-
ymin)/(nbinsy-1);
163 xmin -= xwidth/2; ymin -= ywidth/2;
164 xmax += xwidth/2; ymax += ywidth/2;
166 std::vector<double> xedges(nbinsx+1);
167 std::vector<double> yedges(nbinsy+1);
170 xedges[
i] = xmin + (xmax-
xmin)*
i/
double(nbinsx);
171 if(xlog) xedges[
i] =
exp(xedges[
i]);
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]);
178 return new TH2F(
UniqueName().c_str(), title.c_str(),
180 nbinsy, &yedges.front());
187 const char*
pub =
getenv(
"NOVASOFT_DIR");
190 const char* pub =
getenv(
"FW_BASE");
191 const char* priv =
getenv(
"FW_RELEASE_BASE");
197 if(
stat(ret.c_str(), &junk) == 0)
return ret;
212 std::vector<std::string>
ret;
214 std::ifstream is(listfile);
223 if(!fname.empty()) ret.push_back(fname);
231 std::map<std::string, std::string>
ret;
233 TIter
next(dir->GetListOfKeys());
235 TObject* obj = dir->Get(
key->GetName());
238 const TString className = obj->ClassName();
239 if(className ==
"TObjString"){
240 ret[
key->GetName()] = ((TObjString*)obj)->GetString();
243 std::cerr <<
"Unexpected object " <<
key->GetName() <<
" of type " << className <<
" while looking for metadata. Ignoring" <<
std::endl;
252 const std::map<std::string, std::string>&
add,
253 std::set<std::string>&
mask)
259 if(key ==
"parents")
continue;
274 if (r1.find(r2.substr(1, r2.size() - 2)) != std::string::npos)
284 sum[r1.size()-1] =
',';
289 if(base.find(key) == base.end()){
291 base[
key] =
it.second;
294 if(key ==
"simulated.number_of_spills" ||
295 key ==
"event_count" ||
296 key ==
"online.totalevents"){
299 atoi(base[key].c_str()) +
300 atoi(
it.second.c_str())).Data();
304 if(base[key] !=
it.second) mask.insert(key);
313 const std::map<std::string, std::string>&
meta)
315 TDirectory*
tmp = gDirectory;
319 TObjString
str(
it.second.c_str());
320 str.Write(
it.first.c_str());
353 std::cerr<<
"Aborting. Bad geometry configuration: " << det <<
"\n";
362 static bool cache_set =
false;
364 cache = (
getenv(
"_CONDOR_SCRATCH_DIR") != 0);
374 static int cache = -1;
379 cache = std::atoi(env);
386 std::cout <<
"Stride() called, but CAFANA_STRIDE is not set (--stride not passed?)" <<
std::endl;
398 static int cache = -1;
403 cache = std::atoi(env);
410 std::cout <<
"Offset() called, but CAFANA_OFFSET is not set (--offset not passed?)" <<
std::endl;
422 static int cache = 0;
427 cache = std::atoi(env);
460 std::size_t
Hash(
const std::vector<std::size_t> &vals)
462 std::size_t
seed = 0;
465 seed ^= std::hash<std::size_t>{}(
v) + 0x9e3779b9 + (seed<<6) + (seed>>2);
483 : fHist(h), fxlo(xlo), fxhi(xhi), fNOsc(NOsc)
505 std::vector<double>
s(
fNOsc, 0.);
506 std::vector<double>
c(
fNOsc, 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);
516 s[
n] += y *
sin(
n*x);
517 c[
n] += y *
cos(
n*x);
528 f->SetParameter(0, c[0]/2);
530 f->SetParameter(
n*2-1, s[
n]);
531 f->SetParameter(n*2, c[n]);
537 f->GetYaxis()->SetRangeUser(0.8*min, 1.2*max);
caf::Proxy< caf::SRSpill > spill
Near Detector underground.
Det_t
Which NOvA detector?
size_t JobNumber()
What's the process number for a grid job?
Far Detector at Ash River.
Cuts and Vars for the 2020 FD DiF Study.
std::map< std::string, double > xmax
bool RunningOnGrid()
Is this a grid (condor) job?
caf::Proxy< caf::SRHeader > hdr
Proxy for caf::StandardRecord.
void CombineMetadata(std::map< std::string, std::string > &base, const std::map< std::string, std::string > &add, std::set< std::string > &mask)
base += add
FitToFourier(TH1 *h, double xlo, double xhi, int NOsc)
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
std::string FindCAFAnaDir()
const XML_Char int const XML_Char int const XML_Char * base
double operator()(double *x, double *par) const
FloatingExceptionOnNaN(bool enable=true)
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.
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
int Limit()
Value passed to –limit, or -1 if not specified.
Prototype Near Detector on the Surface.
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.
std::string getenv(std::string const &name)
caf::Proxy< unsigned int > evt
fvar< T > exp(const fvar< T > &x)
void WriteCAFMetadata(TDirectory *dir, const std::map< std::string, std::string > &meta)
Write map of metadata parameters into a CAF file.
size_t Stride(bool allow_default)
Value passed to –stride, or 1 if not specified.
std::string FindPackageDir(std::string Dir)
caf::Proxy< unsigned int > subrun
static float min(const float a, const float b, const float c)
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
size_t Offset(bool allow_default)
Value passed to –offset, or 0 if not specified.
caf::Proxy< unsigned int > run
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Extract map of metadata parameters from a CAF file.
assert(nhit_max >=nhit_nbins)
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
~FloatingExceptionOnNaN()
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Prevent histograms being added to the current directory.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
add("abs", expr_type(int_type()), expr_type(int_type()))
std::string UniqueName()
Return a different string each time, for creating histograms.