ISurface.cxx
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TGraph.h"
3 #include "TKey.h"
4 #include "TMarker.h"
5 #include "TH2.h"
6 #include "TPad.h"
7 #include "TROOT.h"
8 #include "TStyle.h"
9 #include "TVectorD.h"
10 
11 #include "OscLib/IOscCalc.h"
12 
13 #include "CAFAna/Core/Utilities.h"
14 #include "CAFAna/Fit/ISurface.h"
15 
16 namespace ana
17 {
18  //---------------------------------------------------------------------
19  void ISurface::Draw() const
20  {
21  EnsureAxes();
22 
23  fHist->Draw("colz same");
24 
25  // colz obliterated them
26  gPad->RedrawAxis();
27 
28  gPad->Update();
29  }
30 
31  //---------------------------------------------------------------------
32  void ISurface::DrawBestFit(Color_t color, Int_t marker) const
33  {
34  EnsureAxes();
35 
36  TMarker *mark = new TMarker(fBestFitX, fBestFitY, marker);
37  mark->SetMarkerSize(1.5);
38  mark->SetMarkerColor(color);
39  mark->Draw();
40  gPad->Update();
41  }
42 
43  //----------------------------------------------------------------------
44  void ISurface::DrawContour(TH2 *fc, Style_t style, Color_t color,
45  double minchi)
46  {
47  EnsureAxes();
48 
49  std::vector<TGraph *> gs = GetGraphs(fc, minchi);
50 
51  for (TGraph *g: gs)
52  {
53  g->SetLineWidth(3);//2);
54  g->SetLineStyle(style);
55  g->SetLineColor(color);
56  g->Draw("l");
57  }
58 
59  gPad->Update();
60  }
61 
62  //---------------------------------------------------------------------
63  void ISurface::EnsureAxes() const
64  {
65  // Could have a file temporarily open
66  DontAddDirectory guard;
67 
68  // If this pad has already been drawn in, already has axes
69  if (gPad && !gPad->GetListOfPrimitives()->IsEmpty()) return;
70 
71  const TAxis *ax = fHist->GetXaxis();
72  const TAxis *ay = fHist->GetYaxis();
73  const double Nx = ax->GetNbins();
74  const double Ny = ay->GetNbins();
75 
76  // Axes with limits where the user originally requested them, which we
77  // adjusted to be the centres of the first and last bins.
78  TH2 *axes = new TH2C(UniqueName().c_str(),
79  TString::Format(";%s;%s",
80  ax->GetTitle(), ay->GetTitle()),
81  Nx - 1, BinCenterX(1), BinCenterX(Nx),
82  Ny - 1, BinCenterY(1), BinCenterY(Ny));
83  axes->Draw();
84 
85  if(fHist){
86  // "colz same" will reuse axis's min and max, so set them helpfully here
87  axes->SetMinimum(fHist->GetMinimum());
88  axes->SetMaximum(fHist->GetMaximum());
89  }
90 
91  axes->SetTitle(fHist->GetTitle());
92  axes->GetXaxis()->SetLabelSize(ax->GetLabelSize());
93  axes->GetYaxis()->SetLabelSize(ay->GetLabelSize());
94  axes->GetXaxis()->SetLabelOffset(ax->GetLabelOffset());
95  axes->GetYaxis()->SetLabelOffset(ay->GetLabelOffset());
96  axes->GetXaxis()->CenterTitle();
97  axes->GetYaxis()->CenterTitle();
98 
99  if(fLogX) gPad->SetLogx();
100  if(fLogY) gPad->SetLogy();
101 
102  gPad->Update();
103  }
104 
105  //----------------------------------------------------------------------
106  std::vector<TGraph*> ISurface::GetGraphs(TH2* fc, double minchi)
107  {
108  std::vector<TGraph*> ret;
109 
110  if(minchi < 0) minchi = fBestLikelihood;
111  std::unique_ptr<TH2> surf((TH2*)fHist->Clone(UniqueName().c_str()));
112  surf->Add(fc, -1);
113 
114  TVirtualPad* bak = gPad;
115 
116  const bool wasbatch = gROOT->IsBatch();
117  gROOT->SetBatch(); // User doesn't want to see our temporary canvas
118  TCanvas tmp;
119 
120  gStyle->SetOptStat(0);
121 
122  const double level = minchi-fBestLikelihood;
123  surf->SetContour(1, &level);
124  surf->Draw("cont list");
125 
126  tmp.Update();
127  tmp.Paint();
128 
129  gROOT->SetBatch(wasbatch);
130  gPad = bak;
131 
132  // The graphs we need (contained inside TLists, contained inside
133  // TObjArrays) are in the list of specials. But we need to be careful about
134  // types, because other stuff can get in here too (TDatabasePDG for
135  // example).
136  TCollection* specs = gROOT->GetListOfSpecials();
137 
138  TIter nextSpec(specs);
139  while(TObject* spec = nextSpec()){
140  if(!spec->InheritsFrom(TObjArray::Class())) continue;
141  TObjArray* conts = (TObjArray*)spec;
142 
143  if(conts->IsEmpty()) continue;
144 
145  if(!conts->At(0)->InheritsFrom(TList::Class())) continue;
146  TList* cont = (TList*)conts->At(0);
147 
148  TIter nextObj(cont);
149  // Contour could be split into multiple pieces
150  while(TObject* obj = nextObj()){
151  if(!obj->InheritsFrom(TGraph::Class())) continue;
152 
153  ret.push_back((TGraph*)obj->Clone(UniqueName().c_str()));
154  } // end for obj
155  } // end for spec
156 
157  return ret;
158  }
159 
160  //----------------------------------------------------------------------
161  TH2 * ISurface::ToTH2(double minchi) const
162  {
163  // Could have a file temporarily open
164  DontAddDirectory guard;
165 
166  TH2 *ret = new TH2F(*fHist);
167 
168  if (minchi >= 0)
169  {
170  for (int x = 0; x < ret->GetNbinsX() + 2; ++x)
171  {
172  for (int y = 0; y < ret->GetNbinsY() + 2; ++y)
173  {
174  ret->SetBinContent(x, y, ret->GetBinContent(x, y) + fBestLikelihood - minchi);
175  }
176  }
177  }
178 
179  return ret;
180  }
181 
182  //----------------------------------------------------------------------
183  void ISurface::SetTitle(const char *str)
184  {
185  fHist->SetTitle(str);
186  }
187 
188  //----------------------------------------------------------------------
189  void ISurface::SaveToHelper(TDirectory *dir) const
190  {
191  TDirectory *oldDir = gDirectory;
192 
193  dir->cd();
194 
195  TVectorD v(3);
196  v[0] = fBestLikelihood;
197  v[1] = fBestFitX;
198  v[2] = fBestFitY;
199  v.Write("minValues");
200 
201  fHist->Write("hist");
202 
203  TVectorD s(fSeedValues.size(), &fSeedValues[0]);
204  s.Write("seeds");
205 
206  TObjString(fLogX ? "yes" : "no").Write("logx");
207  TObjString(fLogY ? "yes" : "no").Write("logy");
208 
209  if (oldDir)
210  oldDir->cd();
211  }
212 
213  //----------------------------------------------------------------------
215  {
216 
217  const TVectorD v = *(TVectorD *) dir->Get("minValues");
218  const TVectorD s = *(TVectorD *) dir->Get("seeds");
219 
220  surf.fHist = (TH2F *) dir->Get("hist");
221  surf.fBestLikelihood = v[0];
222  surf.fBestFitX = v[1];
223  surf.fBestFitY = v[2];
224 
225  for (int idx = 0; idx < s.GetNrows(); ++idx)
226  surf.fSeedValues.push_back(s[idx]);
227 
228  }
229 
230  //----------------------------------------------------------------------
231  /// Helper function for the gaussian approximation surfaces
232  TH2* Flat(double level, const ISurface& s)
233  {
234  TH2* h = s.ToTH2();
235 
236  for(int x = 0; x < h->GetNbinsX()+2; ++x)
237  for(int y = 0; y < h->GetNbinsY()+2; ++y)
238  h->SetBinContent(x, y, level);
239 
240  return h;
241  }
242 
243  //----------------------------------------------------------------------
244  double ISurface::BinCenterX(int bin) const
245  {
246  const TAxis* ax = fHist->GetXaxis();
247  return fLogX ? ax->GetBinCenterLog(bin) : ax->GetBinCenter(bin);
248  }
249 
250  //----------------------------------------------------------------------
251  double ISurface::BinCenterY(int bin) const
252  {
253  const TAxis* ax = fHist->GetYaxis();
254  return fLogY ? ax->GetBinCenterLog(bin) : ax->GetBinCenter(bin);
255  }
256 
257 } // namespace ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
void SetTitle(const char *str)
Definition: ISurface.cxx:183
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
Float_t tmp
Definition: plot.C:36
std::vector< TGraph * > GetGraphs(TH2 *fc, double minchi=-1)
For expert use, custom painting of contours.
Definition: ISurface.cxx:106
const XML_Char * s
Definition: expat.h:262
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
TH2F * fHist
Definition: ISurface.h:66
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
bool fLogX
Definition: ISurface.h:67
double fBestFitY
Definition: ISurface.h:65
void SaveToHelper(TDirectory *dir) const
dir should already be the appropriate sub-directory
Definition: ISurface.cxx:189
float bin[41]
Definition: plottest35.C:14
bool fLogY
Definition: ISurface.h:67
std::vector< double > fSeedValues
Definition: ISurface.h:68
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
void EnsureAxes() const
Definition: ISurface.cxx:63
TDirectory * dir
Definition: macro.C:5
double fBestFitX
Definition: ISurface.h:64
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
static void FillSurfObj(ISurface &surf, TDirectory *dir)
Definition: ISurface.cxx:214
TH2 * Flat(double level, const ISurface &s)
Helper function for the gaussian approximation surfaces.
Definition: ISurface.cxx:232
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
double fBestLikelihood
Definition: ISurface.h:63
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
surf
Definition: demo5.py:35