FCSurface.cxx
Go to the documentation of this file.
1 #include "CAFAna/FC/FCSurface.h"
2 
3 #include "CAFAna/FC/FCPoint.h"
4 #include "CAFAna/Core/Progress.h"
6 
7 #include "TCanvas.h"
8 #include "TFile.h"
9 #include "TH2.h"
10 #include "TTree.h"
11 #include "TVirtualPad.h"
12 
13 #include <cassert>
14 #include <cmath>
15 #include <iostream>
16 
17 namespace ana
18 {
19  //----------------------------------------------------------------------
20  FCSurface::FCSurface(int nbinsx, double xmin, double xmax,
21  int nbinsy, double ymin, double ymax)
22  {
23  DontAddDirectory guard;
24 
25  fBinning = new TH2F(UniqueName().c_str(), "",
26  nbinsx, xmin, xmax, nbinsy, ymin, ymax);
27 
28  fBins.resize((nbinsx+2)*(nbinsy+2));
29  fPrintedWarning = false;
30  }
31 
32  //----------------------------------------------------------------------
34  {
35  delete fBinning;
36  }
37 
38  //----------------------------------------------------------------------
40  {
41  if(pt.Up()<0)
42  return;
43 
44  int bin;
45  if(plot == "delta")
46  bin = fBinning->FindBin(pt.TrueDelta(), 0);
47  else if(plot == "ssth23")
48  bin = fBinning->FindBin(pt.TrueTh23(), 0);
49  else if(plot == "th13")
50  bin = fBinning->FindBin(pt.TrueTh13(), 0);
51  else if(plot == "dmsq")
52  bin = fBinning->FindBin(pt.TrueDmsq(), 0);
53  else if(plot == "deltassth23")
54  bin = fBinning->FindBin(pt.TrueDelta(), pt.TrueTh23());
55  else if(plot == "ssth23dmsq32") {
56  bin = fBinning->FindBin(pt.TrueTh23(), pt.TrueDmsq());
57  }
58  else if(plot == "deltath13")
59  bin = fBinning->FindBin(pt.TrueTh13(), pt.TrueDelta());
60  else if(plot == "mass")
61  bin = fBinning->FindBin(pt.ProfDmsq(), 0);
62  else {
63  std::cout << "I don't know what plot you're trying to make." << std::endl;
64  std::cout << "You might have to add the option to FCSurface::Add(). " << std::endl;
65  exit(1);
66  }
67 
68  // For deltassth23, since the first and last x bins cover the same range
69  // in sin(delta_CP), we want to merge these two columns of bins
70  if(plot == "deltassth23"){
71  // Get binwidth so we can decide whether trueDelta falls in either
72  // of these bins
73  double binWidth = fBinning->GetXaxis()->GetBinWidth(1);
74  // If trueDelta falls into the first bin, also add it to the last
75  if(pt.TrueDelta() > -1*binWidth/2 && pt.TrueDelta() < binWidth/2) {
76  int mergeBin = fBinning->FindBin(pt.TrueDelta()+2, pt.TrueTh23());
77  fBins[mergeBin].Add(pt);
78  }
79  // If trueDelta falls into the last bin, also add it to the first
80  if(pt.TrueDelta() > 2-binWidth/2 && pt.TrueDelta() < 2+binWidth/2) {
81  int mergeBin = fBinning->FindBin(pt.TrueDelta()-2, pt.TrueTh23());
82  fBins[mergeBin].Add(pt);
83  }
84  }
85 
86  fBins[bin].Add(pt);
87  }
88 
89  //----------------------------------------------------------------------
91  {
92  Progress prog("Filling FCSurface from collection");
93  const unsigned int N = coln.NPoints();
94 
95  for(unsigned int n = 0; n < N; ++n){
96  if(n%100000 == 0) prog.SetProgress(double(n)/N);
97  Add(coln.Point(n), plot);
98  }
99  prog.Done();
100  }
101 
102  //----------------------------------------------------------------------
104  {
105  DontAddDirectory guard;
106 
107  TH2* hRet = new TH2F(*fBinning);
108  hRet->SetName(UniqueName().c_str());
109 
110  // Don't do overflow and underflow
111  for(int x = 1; x < hRet->GetNbinsX()+1; ++x){
112  for(int y = 1; y < hRet->GetNbinsY()+1; ++y){
113  const int bin = hRet->GetBin(x, y);
114  if(fBins[bin].Empty()){
115  double crit2D = -2*log(1-sig);
116  hRet->SetBinContent(bin,crit2D);
117  if (!fPrintedWarning) std::cout << "WARNING, NO FCPOINTS IN THIS BIN, fall back on 2D Gaussian up value." << std::endl;
118  fPrintedWarning = true;
119  continue;
120  }
121  const double up = fBins[bin].UpValueForSignificance(sig);
122  assert(!std::isinf(up) && !std::isnan(up));
123  hRet->SetBinContent(bin, up);
124  }
125  }
126 
127  hRet->SetMinimum(0);
128 
129  return hRet;
130  }
131 
132  //----------------------------------------------------------------------
134  {
135  TVirtualPad* c = gPad;
136  if(!c) c = new TCanvas;
137 
138  const int X = fBinning->GetNbinsX();
139  const int Y = fBinning->GetNbinsY();
140  c->Divide(X, Y);
141 
142  for(int x = 1; x < X+1; ++x){
143  for(int y = 1; y < Y+1; ++y){
144  c->cd(x+(y-1)*X);
145  const int bin = fBinning->GetBin(x, y);
146  fBins[bin].Draw();
147  }
148  }
149  }
150 
151  //----------------------------------------------------------------------
152  const FCBin* FCSurface::GetBin(int x, int y) const
153  {
154  return &fBins[fBinning->GetBin(x, y)];
155  }
156 
157  //----------------------------------------------------------------------
159  {
160  TFile fout(fname.c_str(), "RECREATE");
161  fBinning->Write("bins");
162 
163  TTree tr("fcsurf", "fcsurf");
164  FCBin::InitToTree(&tr);
165  for(const FCBin& bin : fBins) bin.ToTree(&tr);
166  tr.Write();
167  }
168 
169  //----------------------------------------------------------------------
170  std::unique_ptr<FCSurface> FCSurface::LoadFromFile(const std::string& fname)
171  {
172  DontAddDirectory guard;
173 
174  std::unique_ptr<FCSurface> ret(new FCSurface);
175 
176  TFile* fin = TFile::Open(pnfs2xrootd(fname).c_str());
177  assert(!fin->IsZombie());
178 
179  ret->fBinning = (TH2F*)(fin->Get("bins")->Clone());
180 
181  TTree* tr = (TTree*)fin->Get("fcsurf");
182  assert(tr);
183 
185 
186  for(int i = 0; i < tr->GetEntries(); ++i)
187  ret->fBins.push_back(FCBin::FromTree(tr, i));
188 
189  delete fin;
190 
191  return ret;
192  }
193 
194 } // namespace
TString fin
Definition: Style.C:24
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
double TrueTh23() const
Definition: FCPoint.h:33
double TrueTh13() const
Definition: FCPoint.h:34
std::vector< FCBin > fBins
Definition: FCSurface.h:38
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const FCBin * GetBin(int x, int y) const
Definition: FCSurface.cxx:152
unsigned int NPoints() const
Definition: FCCollection.h:26
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Float_t Y
Definition: plot.C:38
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
Definition: FCSurface.cxx:170
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
Double_t ymax
Definition: plot.C:25
bool fPrintedWarning
Definition: FCSurface.h:39
A collection of Feldman-Cousins experiments at the same oscillation parameters.
Definition: FCBin.h:13
TH2F * fBinning
Definition: FCSurface.h:37
void SaveToFile(const std::string &fname) const
Definition: FCSurface.cxx:158
double TrueDelta() const
Definition: FCPoint.h:32
static FCBin FromTree(TTree *tr, int idx)
Definition: FCBin.cxx:125
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double TrueDmsq() const
Definition: FCPoint.h:35
static void InitFromTree(TTree *tr)
Definition: FCBin.cxx:107
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Int_t nbinsx
Definition: plot.C:23
static void InitToTree(TTree *tr)
Definition: FCBin.cxx:89
exit(0)
double ProfDmsq() const
Definition: FCPoint.h:36
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
Double_t ymin
Definition: plot.C:24
double Up() const
Definition: FCPoint.h:38
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t X
Definition: plot.C:38
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const FCPoint & Point(int n) const
Definition: FCCollection.h:27
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12