bin_composition_pie_chart.C
Go to the documentation of this file.
3 
5 #include "CAFAna/Vars/FitVars.h"
7 
9 
10 #include "Utilities/func/MathUtil.h"
11 
12 #include "TCanvas.h"
13 #include "TH1.h"
14 #include "TLatex.h"
15 #include "TLegend.h"
16 #include "TMath.h"
17 #include "TPie.h"
18 #include "TPolyLine.h"
19 #include "TRandom3.h"
20 
21 using namespace ana;
22 
23 const double explode = .05;
24 const double R = .35;
25 
26 const int cols[3] = {kRed-7, kBlue-7, kGray};
27 const int text_cols[3] = {kRed, kBlue, kGray+2};
28 const int internal_line_col = kGray+2;
29 
30 const std::string component_labels[3] = {"Signal", "Beam bkg", "Cosmic bkg"};
31 const std::string sample_labels[3] = {"Low CVN", "High CVN", "Peripheral"};
32 
33 void leg(double x0, double y0, double x1, double y1)
34 {
35  TLegend* leg = new TLegend(x0, y0, x1, y1);
36  leg->SetFillStyle(0);
37  TH1* dummy = new TH1F("", "", 1, 0, 1);
38  dummy->SetLineWidth(3);
39  for(int i = 0; i < 3; ++i){
40  dummy->SetFillColor(cols[i]);
41  leg->AddEntry(dummy->Clone(), component_labels[i].c_str(), "bf");
42  }
43  leg->Draw();
44 }
45 
46 void pie_arc(double cx, double cy, double a0, double a1, double r,
47  int width, int lineCol, int fillCol, double fudgeR=0)
48 {
49  TPolyLine* l = new TPolyLine();
50  l->SetLineWidth(width);
51  l->SetLineColor(lineCol);
52 
53  l->SetPoint(0, cx+fudgeR*r*cos(a0), cy+fudgeR*r*sin(a0));
54  for(int i = 0; i <= 100; ++i){
55  const double a = a0 + (a1-a0)*i/100.;
56  l->SetPoint(i+1, cx+r*cos(a), cy+r*sin(a));
57  }
58  l->SetPoint(102, cx+fudgeR*r*cos(a1), cy+fudgeR*r*sin(a1));
59  l->SetPoint(103, cx+fudgeR*r*cos(a0), cy+fudgeR*r*sin(a0));
60 
61  l->SetFillColor(fillCol);
62  l->SetFillStyle(1001);
63 
64  if(fillCol >= 0) l->Draw("f");
65  l->Draw("l");
66 }
67 
68 double IntInBin(TH1* h, int bin)
69 {
70  int EBinsPerSelBin = 9;
71  int idx1 = 1 + EBinsPerSelBin*bin;
72  int idx2 = EBinsPerSelBin*(bin+1);
73  if (bin <= 2)
74  return h->Integral(idx1,idx2);
75  return h->GetBinContent(21); // Peripheral is merged to bin 30
76 }
77 
79 {
80  Spectrum *cos = LoadFromFile<Spectrum>(fname,name).release();
81  TH1 *ret = cos->ToTH1(kAna2017Livetime,kLivetime);
82  return ret;
83 }
84 
85 
87 {
88  // Need a calc before we can do anything, seed at 2017 best fit
89  double dCP = 1.21;
90  double ssth23 = 0.446;
91  double dmsq32 = 0.00244;
92 
94 
95  kFitDeltaInPiUnits.SetValue(calc,dCP);
96  kFitDmSq32.SetValue(calc,dmsq32);
97  kFitSinSqTheta23.SetValue(calc,ssth23);
98 
99  // Pull the prediction
100  std::string decompStr = isFHC ? "combo" : "prop";
101  std::string beamStr = isFHC ? "fhc" : "rhc";
102  const IPrediction *pred =
103  GetNuePrediction2018(decompStr, calc, false, beamStr, true);
104 
105  double pot = 9.48e20; // From Nitish, update to official once committed
106 
107  TH1 *psig = pred->PredictComponent(calc, Flavors::kNuMuToNuE,
109  if (!isFHC)
110  psig = pred->PredictComponent(calc, Flavors::kNuMuToNuE,
112  TH1 *pbb = pred->Predict(calc).ToTH1(pot);
113  pbb->Add(psig,-1); // beam bkgd = total - sig
114 
115  // And the cosmic prediction
116  TH1D *pcos = GetNueCosmics2018(beamStr).first;
117 
118 
119  TCanvas *canPie = new TCanvas("a", "b", 600, 600);
120 
121  double vals[] = { IntInBin(psig,0), IntInBin(pbb,0), IntInBin(pcos,0),
122  IntInBin(psig,1), IntInBin(pbb,1), IntInBin(pcos,1),
123  IntInBin(psig,2), IntInBin(pbb,2), IntInBin(pcos,2)};
124 
125  double tot = 0;
126  for(int i = 0; i < 9; ++i) tot += vals[i];
127 
128  double tots[3];
129  tots[0] = vals[0]+vals[1]+vals[2];
130  tots[1] = vals[3]+vals[4]+vals[5];
131  tots[2] = vals[6]+vals[7]+vals[8];
132 
133 
134  double accum = 0;
135 
136  double tot_accum = 0;
137  for(int i = 0; i < 3; ++i){
138  const double a0 = tot_accum/tot * 2*TMath::Pi();
139  tot_accum += tots[i];
140  const double a1 = tot_accum/tot * 2*TMath::Pi();
141  const double cx = .5+explode*cos((a0+a1)/2);
142  const double cy = .5+explode*sin((a0+a1)/2);
143 
144  for(int j = i*3; j < (i+1)*3; ++j){
145  const int col = cols[j%3];
146 
147  const double a0j = accum/tot * 2*TMath::Pi();
148  accum += vals[j];
149  const double a1j = accum/tot * 2*TMath::Pi();
150  pie_arc(cx, cy, a0j, a1j, R, 3, internal_line_col, col, 0.02);
151  }
152 
153  const std::string label = sample_labels[i];
154 
155  TLatex* ltx = new TLatex(cx+(R+.03)*cos((a0+a1)/2),
156  cy+(R+.03)*sin((a0+a1)/2),
157  label.c_str());
158  ltx->SetTextAlign(22);
159  double textAng = (a0+a1)/2*180/TMath::Pi()-90;
160  if(textAng < -90 || textAng > +90) textAng += 180;
161  ltx->SetTextAngle(textAng);
162 
163  ltx->Draw();
164 
165  // Draw the thick lines
166  pie_arc(cx, cy, a0, a1, R, 5, kBlack, -1);
167  }
168 
169  leg(0, .8, .25, 1);
170 
171  std::string appenStr = isFHC ? "_FHC" : "_RHC";
172  gPad->Print(("pie"+appenStr+".png").c_str());
173  gPad->Print(("pie"+appenStr+".pdf").c_str());
174 
175  // Begin bars
176 
177  TRandom3 *rand = new TRandom3(); // Poisson fake data, until we have data
178  TCanvas *canBars = new TCanvas();
179  Double_t xbins[] = {0,1,4,5,8,9,12,13};
180  TH1* htot = new TH1F("", ";;Expected Events", 7, xbins);
181  TH1* htot_nofill = new TH1F("", ";;Expected Events", 7, xbins);
182  TH1* hdata = new TH1F("", ";;Expected Events", 7, xbins);
183  TH1* hcos = new TH1F("", "", 7, xbins);
184  TH1* hbkg = new TH1F("", "", 7, xbins);
185  for(int i = 0; i < 3; ++i){
186  htot->SetBinContent(2*i+2, tots[i]);
187  htot_nofill->SetBinContent(2*i+2, tots[i]);
188  hcos->SetBinContent(2*i+2, vals[i*3+2]);
189  hbkg->SetBinContent(2*i+2, vals[i*3+1] + vals[i*3+2]);
190  hdata->SetBinContent(2*i+2, rand->Poisson(tots[i])); // Poisson for now
191  }
192  htot->SetLineWidth(5);
193  htot_nofill->SetLineWidth(5);
194  hdata->SetLineWidth(5);
195  hdata->SetMarkerStyle(8);
196 
197  for(int i = 0; i < 3; ++i){
198  htot->GetXaxis()->SetBinLabel(2+2*i, sample_labels[i].c_str());
199  }
200 
201  htot->GetXaxis()->SetLabelSize(1.5*htot->GetYaxis()->GetTitleSize());
202  htot->GetXaxis()->LabelsOption("h");
203 
204  htot->GetXaxis()->SetTickLength(0);
205 
206  htot->GetYaxis()->CenterTitle();
207  double max = 1.40*htot->GetMaximum();
208  htot->SetMaximum(max);
209 
210  htot->SetFillColor(cols[0]);
211  htot->Draw("hist");
212  //htot->DrawClone();
213  //htot->SetFillStyle(0);
214 
215  hbkg->SetLineWidth(3);
216  hbkg->SetLineColor(internal_line_col);
217  hbkg->SetFillColor(cols[1]);
218  hbkg->Draw("same");
219 
220  hcos->SetLineWidth(3);
221  hcos->SetLineColor(internal_line_col);
222  hcos->SetFillColor(cols[2]);
223  hcos->Draw("same");
224  htot_nofill->Draw("same,hist");
225  hdata->Draw("E,P,same");
226 
227  double xlo = isFHC ? 0.125 : 0.60;
228  TLegend* legBars = new TLegend(xlo, 0.60, xlo+0.25, 0.85);
229  legBars->SetFillStyle(0);
230  TH1* dummy = new TH1F("", "", 1, 0, 1);
231  dummy->SetLineWidth(3);
232  for(int i = 0; i < 3; ++i){
233  dummy->SetFillColor(cols[i]);
234  legBars->AddEntry(dummy->Clone(), component_labels[i].c_str(), "bf");
235  }
236  legBars->AddEntry(hdata, "Data", "lep");
237  legBars->Draw();
238 
239  gPad->Print(("bars"+appenStr+".png").c_str());
240  gPad->Print(("bars"+appenStr+".pdf").c_str());
241 
242  // Begin mosaic
243 
244  TCanvas *canMosaic = new TCanvas("c", "d", 1000, 600);
245 
246  const double padding = .1;
247  const double internal_padding = .025;
248 
249  const double active_width = 1-2*padding-2*internal_padding;
250  const double active_height = 1-2*padding;
251 
252  // Fraction of total number of events in each sample
253  double frac[3] = {0,};
254  for(int i = 0; i < 3; ++i){
255  for(int j = 0; j < 3; ++j){
256  frac[i] += vals[i*3+j]/tot_accum;
257  }
258  }
259 
260  accum = 0;
261  for(int i = 0; i < 3; ++i){
262  const double x0 = padding + active_width*accum;
263  accum += frac[i];
264  const double x1 = padding + active_width*accum;
265  accum += internal_padding;
266  TBox* b = new TBox(x0, padding, x1, 1-padding);
267  b->SetLineWidth(5);
268  b->SetFillStyle(0);
269 
270  // Total in this column
271  const double tot = vals[i*3] + vals[i*3+1] + vals[i*3+2];
272  double accum2 = 0;
273  for(int j = 2; j >= 0; --j){
274  const double y0 = padding + active_height*accum2;
275  accum2 += vals[i*3+j]/tot;
276  const double y1 = padding + active_height*accum2;
277 
278  TBox* b = new TBox(x0, y0, x1, y1);
279  b->SetFillColor(cols[j]);
280  b->SetLineColor(internal_line_col);
281  b->SetLineWidth(3);
282  b->Draw("l");
283 
284  if(i == 0){
285  // This line is usefull to put the labels vertically
286  // on the left of the figure
287  //TLatex* ltx = new TLatex(padding*.75, (y0+y1)/2, component_labels[j].c_str());
288 
289  double labelsplit = 0.03;
290  if (j == 2){
291  TLatex* ltx1 =
292  new TLatex((x0+x1)/2, (y0+y1)/2, "Cosmic bkg");
293  ltx1->SetTextAlign(22);
294  ltx1->Draw();
295  }
296  else if (j == 1){
297  TLatex* ltx1 =
298  new TLatex((x0+x1)/2, (y0+y1)/2, "Beam bkg");
299  ltx1->SetTextAlign(22);
300  ltx1->Draw();
301  }
302  else{
303  TLatex* ltx =
304  new TLatex((x0+x1)/2, (y0+y1)/2, component_labels[j].c_str());
305  ltx->SetTextAlign(22);
306  //ltx->SetTextAngle(90);
307  //ltx->SetTextColor(text_cols[j]);
308  ltx->Draw();
309  }
310  }
311  } // end for j
312 
313  TLatex* ltx = new TLatex((x0+x1)/2, padding/2, sample_labels[i].c_str());
314  ltx->SetTextAlign(22);
315  ltx->Draw();
316 
317  // This is the outline we computed above. Draw it after the rest
318  b->Draw("l");
319  }
320 
321  gPad->Print(("mosaic"+appenStr+".png").c_str());
322  gPad->Print(("mosaic"+appenStr+".pdf").c_str());
323 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:85
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
double ssth23
Antineutrinos-only.
Definition: IPrediction.h:50
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:16
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:67
static bool isFHC
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:115
const int text_cols[3]
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
void pie_arc(double cx, double cy, double a0, double a1, double r, int width, int lineCol, int fillCol, double fudgeR=0)
const char * label
Charged-current interactions.
Definition: IPrediction.h:39
void bin_composition_pie_chart(std::string fname="/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nom.root", std::string predName="pred_nom2017", std::string cname="/nova/ana/nu_e_ana/Ana2017/Predictions/cosmic/cosmic_spectra_prediction.root", std::string cosName="All/cosm_merged_4bin")
General interface to any calculator that lets you set the parameters.
const double kAna2017Livetime
Definition: Exposures.h:200
const double explode
TH1F * a1
Definition: f2_nu.C:476
double dCP
Int_t col[ntarg]
Definition: Style.C:29
const double a
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalculator *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
#define pot
const int cols[3]
const double j
Definition: BetheBloch.cxx:29
double frac(double x)
Fractional part.
float bin[41]
Definition: plottest35.C:14
void leg(double x0, double y0, double x1, double y1)
Neutrinos-only.
Definition: IPrediction.h:49
TH1 * PredictCosmic(std::string fname, std::string name)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
double dmsq32
const int xbins
Definition: MakePlots.C:82
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
const hit & b
Definition: hits.cxx:21
const double R
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const int internal_line_col
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
const std::string sample_labels[4]
const std::string component_labels[3]
double IntInBin(TH1 *h, int bin)
h
Definition: demo3.py:41
TH1F * hbkg
Definition: Xdiff_gwt.C:58