bin_composition_pie_chart.C
Go to the documentation of this file.
1 #include "OscLib/IOscCalc.h"
3 
5 #include "CAFAna/Vars/FitVars.h"
9 
11 
12 #include "Utilities/func/MathUtil.h"
13 
14 #include "TCanvas.h"
15 #include "TH1.h"
16 #include "TLatex.h"
17 #include "TLegend.h"
18 #include "TMath.h"
19 #include "TPie.h"
20 #include "TPolyLine.h"
21 
22 using namespace ana;
23 
24 const double explode = .05;
25 const double R = .35;
26 
27 const int cols[3] = {kRed-7, kBlue-7, kGray};
28 const int text_cols[3] = {kRed, kBlue, kGray+2};
29 const int internal_line_col = kGray+2;
30 
31 const std::string component_labels[3] = {"Signal", "Beam bkg", "Cosmic bkg"};
32 const std::string sample_labels[4] = {"Low CVN", "Mid CVN", "High CVN", "Peripheral"};
33 
34 void leg(double x0, double y0, double x1, double y1)
35 {
36  TLegend* leg = new TLegend(x0, y0, x1, y1);
37  leg->SetFillStyle(0);
38  TH1* dummy = new TH1F("", "", 1, 0, 1);
39  dummy->SetLineWidth(3);
40  for(int i = 0; i < 3; ++i){
41  dummy->SetFillColor(cols[i]);
42  leg->AddEntry(dummy->Clone(), component_labels[i].c_str(), "bf");
43  }
44  leg->Draw();
45 }
46 
47 void pie_arc(double cx, double cy, double a0, double a1, double r,
48  int width, int lineCol, int fillCol, double fudgeR=0)
49 {
50  TPolyLine* l = new TPolyLine();
51  l->SetLineWidth(width);
52  l->SetLineColor(lineCol);
53 
54  l->SetPoint(0, cx+fudgeR*r*cos(a0), cy+fudgeR*r*sin(a0));
55  std::cout << cx << " " << cy << std::endl;
56  for(int i = 0; i <= 100; ++i){
57  const double a = a0 + (a1-a0)*i/100.;
58  l->SetPoint(i+1, cx+r*cos(a), cy+r*sin(a));
59  }
60  l->SetPoint(102, cx+fudgeR*r*cos(a1), cy+fudgeR*r*sin(a1));
61  l->SetPoint(103, cx+fudgeR*r*cos(a0), cy+fudgeR*r*sin(a0));
62 
63  l->SetFillColor(fillCol);
64  l->SetFillStyle(1001);
65 
66  if(fillCol >= 0) l->Draw("f");
67  l->Draw("l");
68 }
69 
70 double IntInBin(TH1* h, int bin)
71 {
72  int EBinsPerSelBin = 9;
73  int idx1 = 1 + EBinsPerSelBin*bin;
74  int idx2 = EBinsPerSelBin*(bin+1);
75  if (bin <= 2)
76  return h->Integral(idx1,idx2);
77  return h->GetBinContent(30); // Peripheral is merged to bin 30
78 }
79 
81 {
82  Spectrum *cos = LoadFromFile<Spectrum>(fname,name).release();
83  TH1 *ret = cos->ToTH1(kAna2017Livetime,kLivetime);
84  std::cout << ret->Integral() << " " << kAna2017Livetime << std::endl;
85  return ret;
86 }
87 
88 
89 void bin_composition_pie_chart(std::string fname="/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nom.root",
90  std::string predName="pred_nom2017",
91  std::string cname="/nova/ana/nu_e_ana/Ana2017/Predictions/cosmic/cosmic_spectra_prediction.root",
92  std::string cosName="All/cosm_merged_4bin")
93 {
94 
95  PredictionExtrap *core = LoadFromFile<PredictionExtrap>(fname,predName+"/predFid/predCore").release();
96  PredictionNoExtrap *fdmc = LoadFromFile<PredictionNoExtrap>(fname,predName+"/predFid/predNoExtrap").release();
98  PredictionNoExtrap *rock = LoadFromFile<PredictionNoExtrap>(fname,predName+"/predRock").release();
99 
100  double nonsScale = 1./10.45;
101  double swapScale = 1./12.91;
102  PredictionAddRock *pred = new PredictionAddRock(fid,rock,nonsScale,swapScale);
103 
104  // Close to one of the SA best fit points
105  double dCP = 1.5;
106  double ssth23 = 0.404;
107  double dmsq32 = 0.0027;
108 
110 
111  kFitDeltaInPiUnits.SetValue(calc,dCP);
112  kFitDmSq32.SetValue(calc,dmsq32);
113  kFitSinSqTheta23.SetValue(calc,ssth23);
114 
115  double pot = 9.48e20; // From Nitish, update to official once committed
116 
117  TH1 *psig = pred->PredictComponent(calc, Flavors::kNuMuToNuE,
119  TH1 *pbb = pred->Predict(calc).ToTH1(pot);
120  pbb->Add(psig,-1); // beam bkgd = total - sig
121  TH1 *pcos = PredictCosmic(cname, cosName);
122 
123 
124  new TCanvas("a", "b", 600, 600);
125 
126  double vals[] = { IntInBin(psig,0), IntInBin(pbb,0), IntInBin(pcos,0),
127  IntInBin(psig,1), IntInBin(pbb,1), IntInBin(pcos,1),
128  IntInBin(psig,2), IntInBin(pbb,2), IntInBin(pcos,2),
129  IntInBin(psig,3), IntInBin(pbb,3), IntInBin(pcos,3) };
130 
131  double tot = 0;
132  for(int i = 0; i < 12; ++i) tot += vals[i];
133 
134  std::cout << tot << " "<< psig->Integral() +pbb->Integral() +pcos->Integral() << std::endl;
135 
136  double tots[4];
137  tots[0] = vals[0]+vals[1]+vals[2];
138  tots[1] = vals[3]+vals[4]+vals[5];
139  tots[2] = vals[6]+vals[7]+vals[8];
140  tots[3] = vals[9]+vals[10]+vals[11];
141 
142  double accum = 0;
143 
144  double tot_accum = 0;
145  for(int i = 0; i < 4; ++i){
146  const double a0 = tot_accum/tot * 2*TMath::Pi();
147  tot_accum += tots[i];
148  const double a1 = tot_accum/tot * 2*TMath::Pi();
149  const double cx = .5+explode*cos((a0+a1)/2);
150  const double cy = .5+explode*sin((a0+a1)/2);
151 
152  for(int j = i*3; j < (i+1)*3; ++j){
153  const int col = cols[j%3];
154 
155  const double a0j = accum/tot * 2*TMath::Pi();
156  accum += vals[j];
157  const double a1j = accum/tot * 2*TMath::Pi();
158  pie_arc(cx, cy, a0j, a1j, R, 3, internal_line_col, col, 0.02);
159  }
160 
161  const std::string label = sample_labels[i];
162 
163  TLatex* ltx = new TLatex(cx+(R+.03)*cos((a0+a1)/2),
164  cy+(R+.03)*sin((a0+a1)/2),
165  label.c_str());
166  ltx->SetTextAlign(22);
167  double textAng = (a0+a1)/2*180/TMath::Pi()-90;
168  if(textAng < -90 || textAng > +90) textAng += 180;
169  ltx->SetTextAngle(textAng);
170 
171  ltx->Draw();
172 
173  // Draw the thick lines
174  pie_arc(cx, cy, a0, a1, R, 5, kBlack, -1);
175  }
176 
177  leg(0, .8, .25, 1);
178 
179  gPad->Print("pie.png");
180  gPad->Print("pie.pdf");
181 
182  // Begin bars
183 
184  new TCanvas;
185  TH1* htot = new TH1F("", ";;Expected Events", 17, 0, 17);
186  TH1* hcos = new TH1F("", "", 17, 0, 17);
187  TH1* hbkg = new TH1F("", "", 17, 0, 17);
188  for(int i = 0; i < 4; ++i){
189  for(int j = 0; j < 3; ++j){
190  htot->SetBinContent(4*i+j+2, tots[i]);
191  hcos->SetBinContent(4*i+j+2, vals[i*3+2]);
192  hbkg->SetBinContent(4*i+j+2, vals[i*3+1] + vals[i*3+2]);
193  }
194  }
195  htot->SetLineWidth(5);
196 
197  for(int i = 0; i < 4; ++i){
198  htot->GetXaxis()->SetBinLabel(3+4*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 
208  htot->SetFillColor(cols[0]);
209  htot->DrawClone();
210  htot->SetFillStyle(0);
211 
212  hbkg->SetLineWidth(3);
213  hbkg->SetLineColor(internal_line_col);
214  hbkg->SetFillColor(cols[1]);
215  hbkg->Draw("same");
216 
217  hcos->SetLineWidth(3);
218  hcos->SetLineColor(internal_line_col);
219  hcos->SetFillColor(cols[2]);
220  hcos->Draw("same");
221 
222  htot->Draw("same");
223 
224  leg(.125, .65, .375, .85);
225 
226  gPad->Print("bars.png");
227  gPad->Print("bars.pdf");
228 
229  // Begin mosaic
230 
231  new TCanvas("c", "d", 1000, 600);
232 
233  const double padding = .1;
234  const double internal_padding = .025;
235 
236  const double active_width = 1-2*padding-3*internal_padding;
237  const double active_height = 1-2*padding;
238 
239  // Fraction of total number of events in each sample
240  double frac[4] = {0,};
241  for(int i = 0; i < 4; ++i){
242  for(int j = 0; j < 3; ++j){
243  frac[i] += vals[i*3+j]/tot_accum;
244  }
245  }
246 
247  accum = 0;
248  for(int i = 0; i < 4; ++i){
249  const double x0 = padding + active_width*accum;
250  accum += frac[i];
251  const double x1 = padding + active_width*accum;
252  accum += internal_padding;
253  TBox* b = new TBox(x0, padding, x1, 1-padding);
254  b->SetLineWidth(5);
255  b->SetFillStyle(0);
256 
257  // Total in this column
258  const double tot = vals[i*3] + vals[i*3+1] + vals[i*3+2];
259  double accum2 = 0;
260  for(int j = 2; j >= 0; --j){
261  const double y0 = padding + active_height*accum2;
262  accum2 += vals[i*3+j]/tot;
263  const double y1 = padding + active_height*accum2;
264 
265  TBox* b = new TBox(x0, y0, x1, y1);
266  b->SetFillColor(cols[j]);
267  b->SetLineColor(internal_line_col);
268  b->SetLineWidth(3);
269  b->Draw("l");
270 
271  if(i == 0){
272  // This line is usefull to put the labels vertically
273  // on the left of the figure
274  //TLatex* ltx = new TLatex(padding*.75, (y0+y1)/2, component_labels[j].c_str());
275 
276  double labelsplit = 0.03;
277  if (j == 2){
278  TLatex* ltx1 = new TLatex((x0+x1)/2, (y0+y1)/2+labelsplit, "Cosmic");
279  ltx1->SetTextAlign(22);
280  ltx1->Draw();
281  TLatex* ltx2 = new TLatex((x0+x1)/2, (y0+y1)/2-labelsplit, "bkg");
282  ltx2->SetTextAlign(22);
283  ltx2->Draw();
284  }
285  else if (j == 1){
286  TLatex* ltx1 = new TLatex((x0+x1)/2, (y0+y1)/2+labelsplit, "Beam");
287  ltx1->SetTextAlign(22);
288  ltx1->Draw();
289  TLatex* ltx2 = new TLatex((x0+x1)/2, (y0+y1)/2-labelsplit, "bkg");
290  ltx2->SetTextAlign(22);
291  ltx2->Draw();
292  }
293  else{
294  TLatex* ltx = new TLatex((x0+x1)/2, (y0+y1)/2, component_labels[j].c_str());
295  ltx->SetTextAlign(22);
296  //ltx->SetTextAngle(90);
297  //ltx->SetTextColor(text_cols[j]);
298  ltx->Draw();
299  }
300  }
301  } // end for j
302 
303  TLatex* ltx = new TLatex((x0+x1)/2, padding/2, sample_labels[i].c_str());
304  ltx->SetTextAlign(22);
305  ltx->Draw();
306 
307  // This is the outline we computed above. Draw it after the rest
308  b->Draw("l");
309  }
310 
311  gPad->Print("mosaic.png");
312  gPad->Print("mosaic.pdf");
313 }
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double ssth23
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:148
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
const int text_cols[3]
const int internal_line_col
TH1 * PredictCosmic(std::string fname, std::string name)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const double explode
const char * label
double IntInBin(TH1 *h, int bin)
Charged-current interactions.
Definition: IPrediction.h:39
const double kAna2017Livetime
Definition: Exposures.h:200
const int cols[3]
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
TH1F * a1
Definition: f2_nu.C:476
const double R
double dCP
Int_t col[ntarg]
Definition: Style.C:29
const double a
#define pot
const std::string sample_labels[4]
const double j
Definition: BetheBloch.cxx:29
double frac(double x)
Fractional part.
float bin[41]
Definition: plottest35.C:14
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
const std::string component_labels[3]
T sin(T number)
Definition: d0nt_math.hpp:132
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double dmsq32
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const hit & b
Definition: hits.cxx:21
void pie_arc(double cx, double cy, double a0, double a1, double r, int width, int lineCol, int fillCol, double fudgeR=0)
void leg(double x0, double y0, double x1, double y1)
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
Prediction that just uses FD MC, with no extrapolation.
Take the output of an extrapolation and oscillate it as required.
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
enum BeamMode kBlue
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")
virtual Spectrum Predict(osc::IOscCalc *calc) const override
void rock(std::string suffix="full")
Definition: rock.C:28
TH1F * hbkg
Definition: Xdiff_gwt.C:58
enum BeamMode string