Fittingsidebandfittest.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <vector>
4 #include <cmath>
5 #include <stdint.h>
6 #include "TFile.h"
7 #include "TTree.h"
8 #include "TCanvas.h"
9 #include "TF1.h"
10 #include "TH1F.h"
11 #include "TH2F.h"
12 #include <sstream>
13 #include <math.h>
14 #include "TF2.h"
15 #include "TH2.h"
16 #include "TCutG.h"
17 #include "TMath.h"
18 #include "TCanvas.h"
19 #include "TStyle.h"
20 #include "TRandom.h"
21 #include "TColor.h"
22 #include "TF3.h"
23 #include "TH3.h"
24 #include "TH3F.h"
25 #include "TMinuit.h"
26 #include "TLegend.h"
27 #include "THStack.h"
28 
29 //Global Constants used in the Fitting procedure
30 
31 //These bins, try to fit over the peak region
32 //and not the low statistics sidebands
33 const int xfitbins_low = 1;
34 const int xfitbins_high = 9;
35 const int yfitbins_low = 2;
36 const int yfitbins_high = 12;
37 const int xbins = 9;
38 const int ybins = 25;
39 
40 //Arrays to hold the histogram bin content
41 //Need the fitter to have access to this information without
42 //needing to open the histograms multiple times
43 float h_data[xbins][ybins];
44 float h_sig[xbins][ybins];
45 float h_numu[xbins][ybins];
46 float h_nc[xbins][ybins];
47 float h_anti[xbins][ybins];
49 
50 //These hold errors
51 float e_data[xbins][ybins];
52 float e_sig[xbins][ybins];
53 float e_numu[xbins][ybins];
54 float e_nc[xbins][ybins];
55 float e_anti[xbins][ybins];
57 
58 //Constant that holds POT ratio for data and MC
59 double pot_ratio = 1.4;
60 
61 
62 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
63 {
64  double chisq = 0;
65  double Data = 0;
66  double Signal = 0;
67  double Numu = 0;
68  double NC = 0;
69  double Anti = 0;
70  double Other = 0;
71  double Error = 0;
72  int dof = 0;
73  for( int i = 0; i < xfitbins_high; ++i){
74  for( int j = yfitbins_low; j < yfitbins_high; ++j){
75  Data = h_data[i][j];
76  Signal = h_sig[i][j];
77  Numu = h_numu[i][j];
78  NC = h_nc[i][j];
79  Anti = h_anti[i][j];
80  Other = h_other[i][j];
81  Error = e_data[i][j] * e_data[i][j] +
82  e_sig[i][j] * e_sig[i][j] +
83  e_numu[i][j] * e_numu[i][j] +
84  e_nc[i][j] * e_nc[i][j] +
85  e_anti[i][j] * e_anti[i][j] +
86  e_other[i][j] * e_other[i][j];
87  double Delta = 0;
88 
89  if( Data < 85 || pot_ratio * (Signal + Numu + NC + Anti + Other) < 85)
90  continue;
91  if( Error == 0) continue;
92  Delta = (Data - (par[0]*Numu + par[1]*NC + Signal + Anti + Other));
93  dof++;
94  chisq += Delta*Delta/Error;
95  } //end loop over y bins
96  }//end loop over x bins
97  std::cout << dof << std::endl;
98  f = chisq;
99 
100 }
101 
102 void MakePlots( TH2F data, TH2F signal, TH2F numu, TH2F nc,
103  TH2F anti, TH2F other, double scaling[3], std::string title)
104 {
105  auto data1 = (TH2F*)data.Clone("d1");
106  auto sig1 = (TH2F*)signal.Clone("d2");
107  auto numu1 = (TH2F*)numu.Clone("d3");
108  auto nc1 = (TH2F*)nc.Clone("d4");
109  auto anti1 = (TH2F*)anti.Clone("d5");
110  auto other1 = (TH2F*)other.Clone("d6");
111 
112 
113  auto data2 = (TH2F*)data.Clone("s1");
114  auto sig2 = (TH2F*)signal.Clone("s2");
115  auto numu2 = (TH2F*)numu.Clone("s3");
116  auto nc2 = (TH2F*)nc.Clone("s4");
117  auto anti2 = (TH2F*)anti.Clone("s5");
118  auto other2 = (TH2F*)other.Clone("s6");
119 
120  numu2->Scale(scaling[0]);
121  nc2->Scale(scaling[1]);
122  sig2->Scale(scaling[2]);
123 
124  other1->Add(sig1);
125  other1->Add(numu1);
126  other1->Add(nc1);
127  other1->Add(anti1);
128 
129  other2->Add(sig2);
130  other2->Add(numu2);
131  other2->Add(nc2);
132  other2->Add(anti2);
133 
134  data1->Divide(other1);
135  data2->Divide(other2);
136 
137 
138  std::string savename = title + "_before_fitting.png";
139 
140  TCanvas *c1 = new TCanvas("c1","c1",695,700);
141  c1->SetRightMargin(0.15);
142  data1->Draw("COLZ");
143  c1->SaveAs(savename.c_str());
144  c1->Close();
145  savename = title + "_after_fitting.png";
146  TCanvas *c2 = new TCanvas("c2","c2",695,700);
147  c2->SetRightMargin(0.15);
148  data2->Draw("COLZ");
149  c2->SaveAs(savename.c_str());
150  c2->Close();
151 
152  TFile *outf = new TFile("ScalingHistogram.root", "recreate");
153  other2->Divide(other1);
154  other2->Write("weights");
155  outf->Close();
156 
157 
158 }
159 
160 
161 void RunFitter( TH2F data, TH2F signal,
162  TH2F nc, TH2F anti,
163  TH2F other, TH2F numu, std::string title)
164 {
165  //Fill the arrays with bin content
166  for( int i = 0; i < xbins; ++i){
167  for( int j = 0; j < ybins; ++j){
168  h_data[i][j] = data.GetBinContent(i+1,j+1);
169  h_sig[i][j] = signal.GetBinContent(i+1,j+1);
170  h_numu[i][j] = numu.GetBinContent(i+1,j+1);
171  h_nc[i][j] = nc.GetBinContent(i+1,j+1);
172  h_anti[i][j] = anti.GetBinContent(i+1,j+1);
173  h_other[i][j] = other.GetBinContent(i+1,j+1);
174  e_data[i][j] = data.GetBinError(i+1,j+1);
175  e_sig[i][j] = signal.GetBinError(i+1,j+1);
176  e_numu[i][j] = numu.GetBinError(i+1,j+1);
177  e_nc[i][j] = nc.GetBinError(i+1,j+1);
178  e_anti[i][j] = anti.GetBinError(i+1,j+1);
179  e_other[i][j] = other.GetBinError(i+1,j+1);
180  //std::cout << h_data[i][j] << std::endl;
181 
182  } // end y bin loop
183  } //end x bin loop
184 
185  //Setup the Fit
186  TMinuit *gMinuit = new TMinuit(2);
187  gMinuit->SetFCN(fcn);
188  Double_t arglist[6];
189  Int_t ierflg = 0;
190  arglist[0] = 1;
191  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
192  double vstart[3] = { 1, 1, 1};
193  double step[3] = {0.1, 0.1, 0.1};
194  gMinuit->mnparm(0, "NumuCC Scaling", vstart[0], step[0],0,3,ierflg);
195  gMinuit->mnparm(1, "NC Scaling" , vstart[1], step[1],0,3,ierflg);
196  gMinuit->mnparm(2, "Signal Scaling", vstart[2], step[2],0,3,ierflg);
197 
198  arglist[0] = 0;
199  arglist[1] = 0;
200 
201  gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
202  gMinuit->mnexcm("MIGRAD" , arglist, 3, ierflg);
203 
204  //Grab the scaling parameters
205  Double_t scaling_factor[3];
206  Double_t err[3];
207  for( int i = 0; i < 3; ++i){
208  gMinuit->GetParameter(i, scaling_factor[i], err[i]);
209  }
210 
211  MakePlots(data,signal,numu,nc,anti,other,scaling_factor, title);
212 }
213 
214 
216 {
217  TFile *file_sideband = new TFile("SidebandFitTestHistograms.root", "READ");
218 
219  if( !(file_sideband) ){
220  std::cout << "SidebandFitTestHistograms.root does not exist." <<
221  " Run Plotsidebandfittest.C first" << std::endl;
222  return;
223  }
224 
225  //Load Fake Data From the ROOT File
226  TH2F* data_reweighted = (TH2F*)file_sideband->Get("data_reweighted");
227  TH2F* data_tufts = (TH2F*)file_sideband->Get("data_tufts");
228 
229  //Now Load MC Components
230  TH2F* mc_sig = (TH2F*)file_sideband->Get("mc_sig");
231  TH2F* mc_numu = (TH2F*)file_sideband->Get("mc_numu");
232  TH2F* mc_nc = (TH2F*)file_sideband->Get("mc_nc");
233  TH2F* mc_anti = (TH2F*)file_sideband->Get("mc_anti");
234  TH2F* mc_other = (TH2F*)file_sideband->Get("mc_other");
235 
236 
237  RunFitter( *(data_tufts), *(mc_sig),
238  *(mc_nc), *(mc_anti),
239  *(mc_other), *(mc_numu), "data_tufts_comparison");
240 
241 
242  RunFitter( *(data_reweighted), *(mc_sig),
243  *(mc_nc), *(mc_anti),
244  *(mc_other), *(mc_numu), "data_reweighted_comparison");
245 
246 }
247 
float e_numu[xbins][ybins]
float e_other[xbins][ybins]
const int xbins
float h_numu[xbins][ybins]
double pot_ratio
float e_anti[xbins][ybins]
Int_t par
Definition: SimpleIterate.C:24
void RunFitter(TH2F data, TH2F signal, TH2F nc, TH2F anti, TH2F other, TH2F numu, std::string title)
const int xfitbins_high
const int yfitbins_high
float h_other[xbins][ybins]
float e_nc[xbins][ybins]
void scaling(TH1D *hIn, const double shape_scale)
const XML_Char const XML_Char * data
Definition: expat.h:268
c2
Definition: demo5.py:33
float e_sig[xbins][ybins]
void MakePlots(TH2F data, TH2F signal, TH2F numu, TH2F nc, TH2F anti, TH2F other, double scaling[3], std::string title)
const int yfitbins_low
void Fittingsidebandfittest()
TFile * outf
Definition: testXsec.C:51
const double j
Definition: BetheBloch.cxx:29
float h_nc[xbins][ybins]
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
OStream cout
Definition: OStream.cxx:6
const int ybins
float h_sig[xbins][ybins]
const int NC
enum BeamMode nc
c1
Definition: demo5.py:24
const int xfitbins_low
float h_anti[xbins][ybins]
float e_data[xbins][ybins]
float h_data[xbins][ybins]
enum BeamMode string