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