UtilsExt.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/UtilsExt.h"
2 
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/Ratio.h"
5 #include "CAFAna/Core/Spectrum.h"
6 
7 #include "TDirectory.h"
8 #include "TH2.h"
9 #include "TH3.h"
10 
11 #include <iostream>
12 
13 #include "wordexp.h"
14 
15 #include "sys/stat.h"
16 
17 namespace ana
18 {
19  //----------------------------------------------------------------------
21  {
22  static int N = 0;
23  return TString::Format("cafanauniq%d", N++).Data();
24  }
25 
26  //----------------------------------------------------------------------
28  {
29  fBackup = TH1::AddDirectoryStatus();
30  TH1::AddDirectory(false);
31  }
32 
33  //----------------------------------------------------------------------
35  {
36  TH1::AddDirectory(fBackup);
37  }
38 
39  //----------------------------------------------------------------------
40  bool AlmostEqual(double a, double b, double eps)
41  {
42  if(a == 0 && b == 0) return true;
43 
44  return fabs(a-b) <= eps * std::max(fabs(a), fabs(b));
45  }
46 
47  // Helper functions for MakeTHND().
48  namespace{
49  // Eventually the bin parameters will all be unpacked and we just pass them
50  // on to the regular constructor.
51  template<class T, class... A> T* MakeHist(A... a)
52  {
53  DontAddDirectory guard;
54  return new T(a...);
55  }
56 
57  // This function consumes bins from the start of the argument list and
58  // pushes their translations onto the list of arguments at the end.
59  template<class T, class... A> T* MakeHist(const Binning& firstBin,
60  A... args)
61  {
62  if(firstBin.IsSimple())
63  return MakeHist<T>(args...,
64  firstBin.NBins(), firstBin.Min(), firstBin.Max());
65  else
66  return MakeHist<T>(args...,
67  firstBin.NBins(), &firstBin.Edges().front());
68  }
69  }
70 
71  // Concrete instantiations. MakeHist() requires us to put the bin arguments
72  // first...
73  //----------------------------------------------------------------------
74  TH1D* MakeTH1D(const char* name, const char* title, const Binning& bins)
75  {
76  return MakeHist<TH1D>(bins, name, title);
77  }
78 
79  //----------------------------------------------------------------------
80  TH1F* MakeTH1F(const char* name, const char* title, const Binning& bins)
81  {
82  return MakeHist<TH1F>(bins, name, title);
83  }
84 
85  //----------------------------------------------------------------------
86  TH2D* MakeTH2D(const char* name, const char* title,
87  const Binning& binsx,
88  const Binning& binsy)
89  {
90  return MakeHist<TH2D>(binsx, binsy, name, title);
91  }
92 
93  //----------------------------------------------------------------------
94  TH2F* MakeTH2F(const char* name, const char* title,
95  const Binning& binsx,
96  const Binning& binsy)
97  {
98  return MakeHist<TH2F>(binsx, binsy, name, title);
99  }
100 
101  //----------------------------------------------------------------------
102  TH3D* MakeTH3D(const char* name, const char* title,
103  const Binning& binsx,
104  const Binning& binsy,
105  const Binning& binsz)
106  {
107  return MakeHist<TH3D>(name, title,
108  binsx.NBins(), &binsx.Edges().front(),
109  binsy.NBins(), &binsy.Edges().front(),
110  binsz.NBins(), &binsz.Edges().front());
111  }
112 
113 
114  //----------------------------------------------------------------------
115  TH2* ToTH2(const Spectrum& s, double exposure, ana::EExposureType expotype,
116  const Binning& binsx, const Binning& binsy, ana::EBinType bintype)
117  {
118  DontAddDirectory guard;
119 
120  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
121  return ToTH2Helper(std::move(h1), binsx, binsy, bintype);
122  }
123 
124  //----------------------------------------------------------------------
125  TH2* ToTH2(const Ratio& r,
126  const Binning& binsx, const Binning& binsy)
127  {
128  DontAddDirectory guard;
129 
130  std::unique_ptr<TH1> h1(r.ToTH1());
131  return ToTH2Helper(std::move(h1), binsx, binsy);
132  }
133 
134  //----------------------------------------------------------------------
135  TH2* ToTH2Helper(std::unique_ptr<TH1> h1,
136  const Binning& binsx, const Binning& binsy,
137  ana::EBinType bintype)
138  {
139  // Make sure it's compatible with having been made with this binning
140  assert(h1->GetNbinsX() == binsx.NBins()*binsy.NBins());
141 
142  TH2* ret = MakeTH2D("", UniqueName().c_str(), binsx, binsy);
143 
144  for(int i = 0; i < h1->GetNbinsX(); ++i){
145  const double val = h1->GetBinContent(i+1);
146  const double err = h1->GetBinError(i+1);
147 
148  const int ix = i / binsy.NBins();
149  const int iy = i % binsy.NBins();
150 
151  ret->SetBinContent(ix+1, iy+1, val);
152  ret->SetBinError (ix+1, iy+1, err);
153  }
154 
155  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
156 
157  return ret;
158  }
159 
160  //----------------------------------------------------------------------
161 
162  TH3* ToTH3(const Spectrum& s, double exposure, ana::EExposureType expotype,
163  const Binning& binsx, const Binning& binsy, const Binning& binsz,
164  ana::EBinType bintype)
165  {
166  DontAddDirectory guard;
167 
168  std::unique_ptr<TH1> h1(s.ToTH1(exposure, expotype));
169 
170  return ToTH3Helper(std::move(h1), binsx, binsy, binsz, bintype);
171  }
172 
173  //----------------------------------------------------------------------
174 
175  TH3* ToTH3(const Ratio& r,
176  const Binning& binsx, const Binning& binsy, const Binning& binsz)
177  {
178  DontAddDirectory guard;
179 
180  std::unique_ptr<TH1> h1(r.ToTH1());
181 
182  return ToTH3Helper(std::move(h1), binsx, binsy, binsz);
183  }
184 
185  //----------------------------------------------------------------------
186  TH3* ToTH3Helper(std::unique_ptr<TH1> h1,
187  const Binning& binsx,
188  const Binning& binsy,
189  const Binning& binsz,
190  ana::EBinType bintype)
191  {
192  const int nx = binsx.NBins();
193  const int ny = binsy.NBins();
194  const int nz = binsz.NBins();
195 
196  // Make sure it's compatible with having been made with this binning
197  assert(h1->GetNbinsX() == nx*ny*nz);
198 
199  TH3* ret;
200 
201  // If all three axes are simple, we can call a simpler constructor
202  if(binsx.IsSimple() && binsy.IsSimple() && binsz.IsSimple()){
203  ret = new TH3F(UniqueName().c_str(), "",
204  nx, binsx.Min(), binsx.Max(),
205  ny, binsy.Min(), binsy.Max(),
206  nz, binsz.Min(), binsz.Max());
207  }
208  else{
209  // If >= 1 axis has a not simple binning, you have to pass the
210  // edges vector to the TH3 constructor
211  ret = new TH3F(UniqueName().c_str(), "",
212  nx, &binsx.Edges().front(),
213  ny, &binsy.Edges().front(),
214  nz, &binsz.Edges().front());
215  }
216 
217  for(int i = 0; i < h1->GetNbinsX(); ++i){
218  const double val = h1->GetBinContent(i+1);
219  const double err = h1->GetBinError(i+1);
220 
221  const int nynz = ny*nz;
222  const int nmodnynz = i%nynz;
223  const int ix = i/nynz;
224  const int iy = nmodnynz/nz;
225  const int iz = i%nz;
226 
227  ret->SetBinContent(ix+1, iy+1, iz+1, val);
228  ret->SetBinError (ix+1, iy+1, iz+1, err);
229  }
230 
231  if(bintype == ana::EBinType::kBinDensity) ret->Scale(1, "width");
232 
233  return ret;
234  }
235 
236  //----------------------------------------------------------------------
238  {
239  static bool first = true;
240  static bool onsite = false;
241 
242  if (first && unauth) {
243  first = false;
244  char chostname[255];
245  gethostname(chostname, 255);
246  std::string hostname = chostname;
247 
248  if ( hostname.find("fnal.gov") != std::string::npos ){
249  onsite = true;
250  std::cout << "Using unauthenticated xrootd access (port 1095) while on-site, hostname: " << hostname << std::endl;
251  }
252  else {
253  onsite = false;
254  std::cout << "Using authenticated xrootd access (port 1094) access while off-site, hostname: " << hostname << std::endl;
255  }
256  }
257 
258  if(loc.rfind("/pnfs/", 0) == 0){ // ie begins with
259  if ( onsite && unauth )
260  loc = std::string("root://fndca1.fnal.gov:1095//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
261  else
262  loc = std::string("root://fndca1.fnal.gov:1094//pnfs/fnal.gov/usr/")+&loc.c_str()[6];
263  }
264  return loc;
265  }
266 
267  //----------------------------------------------------------------------
268  std::vector<std::string> Wildcard(const std::string& wildcardString)
269  {
270  // Expand environment variables and wildcards like the shell would
271  wordexp_t p;
272  const int status = wordexp(wildcardString.c_str(), &p, WRDE_SHOWERR);
273 
274  if(status != 0){
275  std::cerr << "Wildcard string '" << wildcardString
276  << "' returned error " << status << " from wordexp()."
277  << std::endl;
278  return {};
279  }
280 
281  std::vector<std::string> fileList;
282 
283  for(unsigned int i = 0; i < p.we_wordc; ++i){
284  // Check the file exists before adding it
285  struct stat sb;
286  if(stat(p.we_wordv[i], &sb) == 0)
287  fileList.push_back(p.we_wordv[i]);
288  }
289 
290  wordfree(&p);
291 
292  return fileList;
293  }
294 
295  //----------------------------------------------------------------------
297  {
298  // I would be much much happier to do this in proper code, but I'm not sure
299  // how, there's no samweb C++ API?
300  return system(TString::Format("samweb list-definitions --defname %s | grep %s", def.c_str(), def.c_str()).Data()) == 0;
301  }
302 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
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
system("rm -rf microbeam.root")
int status
Definition: fabricate.py:1613
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
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:148
TH1F * MakeTH1F(const char *name, const char *title, const Binning &bins)
Definition: UtilsExt.cxx:80
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: UtilsExt.cxx:186
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
Definition: UtilsExt.cxx:74
TH2F * MakeTH2F(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: UtilsExt.cxx:94
Divide bin contents by bin widths.
Definition: UtilsExt.h:31
OStream cerr
Definition: OStream.cxx:7
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
void MakeHist(int n)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
list fileList
Definition: cafExposure.py:30
bool IsSimple() const
Definition: Binning.h:33
EBinType
Definition: UtilsExt.h:28
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: UtilsExt.cxx:268
double Min() const
Definition: Binning.h:30
bool AlmostEqual(double a, double b, double eps)
Definition: UtilsExt.cxx:40
const double a
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: UtilsExt.h:35
Represent the ratio between two spectra.
Definition: Ratio.h:8
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
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: UtilsExt.cxx:86
double Max() const
Definition: Binning.h:31
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1F * h1
static const double A
Definition: Units.h:82
bool SAMDefinitionExists(const std::string &def)
Definition: UtilsExt.cxx:296
int NBins() const
Definition: Binning.h:29
const hit & b
Definition: hits.cxx:21
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: UtilsExt.cxx:135
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
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: UtilsExt.cxx:102
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
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: UtilsExt.cxx:162
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
enum BeamMode string