Public Member Functions | Public Attributes | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
ana::nueccinc::NueCCIncTemplateFitter Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-12-03/NDAna/nuecc_inc/NueCCIncTemplateFitter.h"

Public Member Functions

 NueCCIncTemplateFitter (TH3F *da, std::vector< TH3F * > templates, std::vector< TH3F * > analysis, TH3F *cv, std::vector< TH3F * > systs_up, std::vector< TH3F * > systs_down, std::vector< TH3F * > systs_oneside, std::vector< TH3F * > multiverse, int numbins)
 
virtual ~NueCCIncTemplateFitter ()
 
void PrintValueHolders (int binx, int biny)
 
void FCNMinuit ()
 
std::vector< TH3F * > doFullFit ()
 
std::vector< TH1F * > getPIDFitTemplates (int binx, int biny)
 
std::vector< std::pair< TH2F *, TH2F * > > getFitNormalizationAndErrors ()
 
std::vector< TH1F * > getPIDTemplates (int binx, int biny)
 
TH2F * getCovarianceMatrixTemplateBins (int binx, int biny)
 
TH2F * getCorrelationMatrixTemplateBins (int binx, int biny)
 
TH2F * getCovarianceBetweenSysts (int binx, int biny)
 
TH2F * getCorrelationBetweenSysts (int binx, int biny)
 

Public Attributes

std::pair< double, double > nue_wei
 
std::pair< double, double > numu_wei
 
std::pair< double, double > nc_wei
 
double fit_covariance_matrix [3][3]
 
double chisquared
 

Private Member Functions

void FillValueHolders (int binx, int biny)
 
bool CheckSystsSize ()
 
std::vector< TH3F * > GetAnalysisTemplates ()
 
std::string GetAnalysisBinCaption (int binx, int biny)
 
bool CheckIfBinShouldBeUsed (int binx, int biny)
 
std::vector< TH1F * > CalculateOneSigmaShift (int binx, int biny)
 
void FillCovarianceMatrix (int binx, int biny, bool makePlots)
 
void FillStatisticCovarianceMatrix (int binx, int biny)
 
double GetCovarianceMatrixValue (int binx, int biny)
 
double GetStatisticCovarianceMatrixValue (int binx, int biny)
 
double myFunction3Vars (double numu_par, double nue_par, double nc_par)
 
double myFunction2Vars (double bkgd_par, double nue_par)
 
void SetBinWeightsAndErrors2D (int binx, int biny)
 
void PropagateFitUncertainty3D (int binx, int biny)
 
std::vector< std::pair< double, double > > GetTemplateWeightsAndErrors (int binx, int biny)
 
std::vector< TH1F * > GetTemplates (int binx, int biny, std::string name)
 
TH2F * FillCovarianceBetweenSysts (int binx, int biny)
 

Static Private Member Functions

static void fcn3Var (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 
static void fcn2Var (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 

Private Attributes

TH3F * hData3D
 
std::vector< TH3F * > hTemplates3D
 
std::vector< TH3F * > hAnalysis3D
 
TH3F * hCV3D
 
std::vector< TH3F * > hSystsUp3D
 
std::vector< TH3F * > hSystsDown3D
 
std::vector< TH3F * > hSystsOneSided
 
std::vector< TH3F * > hSystsMultiverse
 
const int NumberTemplateBins
 
TH2F * hSignalWeights2D
 
TH2F * hNumuCCWeights2D
 
TH2F * hNCWeights2D
 
TH2F * hSignalError2D
 
TH2F * hNumuCCError2D
 
TH2F * hNCError2D
 
TH2F * hChiSquaredResults2D
 
double DataValues [numPIDBins]
 
double SignalLike_Values [numPIDBins]
 
double NC_Values [numPIDBins]
 
double NumuLike_Values [numPIDBins]
 
double Other_Values [numPIDBins]
 
double CovarianceMatrix [numPIDBins][numPIDBins]
 
double StatCovarianceMatrix [numPIDBins][numPIDBins]
 

Detailed Description

Definition at line 64 of file NueCCIncTemplateFitter.h.

Constructor & Destructor Documentation

ana::nueccinc::NueCCIncTemplateFitter::NueCCIncTemplateFitter ( TH3F *  da,
std::vector< TH3F * >  templates,
std::vector< TH3F * >  analysis,
TH3F *  cv,
std::vector< TH3F * >  systs_up,
std::vector< TH3F * >  systs_down,
std::vector< TH3F * >  systs_oneside,
std::vector< TH3F * >  multiverse,
int  numbins 
)

Definition at line 28 of file NueCCIncTemplateFitter.cxx.

References make_pair().

32  :hData3D(da),
33  hTemplates3D(templates),
34  hAnalysis3D(analysis),
35  hCV3D(cv),
36  hSystsUp3D(systs_up),
37  hSystsDown3D(systs_down),
38  hSystsOneSided(systs_oneside),
39  hSystsMultiverse(multiverse),
40  NumberTemplateBins(numbins)
41  {
42  hSignalWeights2D = (TH2F*)hCV3D->Project3D("yx");
43  hSignalWeights2D->SetName("hSignalWeights2D");
44  hNumuCCWeights2D = (TH2F*)hCV3D->Project3D("yx");
45  hNumuCCWeights2D->SetName("hNumuCCWeights2D");
46  hNCWeights2D = (TH2F*)hCV3D->Project3D("yx");
47  hNCWeights2D->SetName("hNCWeights2D");
48  hSignalError2D = (TH2F*)hCV3D->Project3D("yx");
49  hSignalError2D->SetName("hSignalError2D");
50  hNumuCCError2D = (TH2F*)hCV3D->Project3D("yx");
51  hNumuCCError2D->SetName("hNumuCCError2D");
52  hNCError2D = (TH2F*)hCV3D->Project3D("yx");
53  hNCError2D->SetName("hNCError2D");
54  hChiSquaredResults2D = (TH2F*)hCV3D->Project3D("yx");
55  hChiSquaredResults2D->SetName("hChiSquaredResults");
56 
57  nue_wei = std::make_pair(0,0);
58  numu_wei = std::make_pair(0,0);
59  nc_wei = std::make_pair(0,0);
60  chisquared = 0;
61  }
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const std::string cv[Ncv]
virtual ana::nueccinc::NueCCIncTemplateFitter::~NueCCIncTemplateFitter ( )
inlinevirtual

Definition at line 75 of file NueCCIncTemplateFitter.h.

References nue_wei.

75 {};

Member Function Documentation

std::vector< TH1F * > ana::nueccinc::NueCCIncTemplateFitter::CalculateOneSigmaShift ( int  binx,
int  biny 
)
private

Definition at line 165 of file NueCCIncTemplateFitter.cxx.

References bin, stan::math::fabs(), hCV3D, hSystsDown3D, hSystsOneSided, hSystsUp3D, MECModelEnuComparisons::i, makeTrainCVSamples::int, and ana::UniqueName().

Referenced by FillCovarianceBetweenSysts(), and FillCovarianceMatrix().

167  {
168  //Function to calculate shift from one-sided and two sided systematics
169  std::vector<TH1F*> hResults;
170 
171 
172  TH1F* hCV =
173  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
174  binx,binx,biny,biny);
175 
176  //One Sided Systematics
177  for(int i = 0; i < (int)hSystsOneSided.size();i++){
178  TH1F* hHolder =
179  (TH1F*)hSystsOneSided[i]->ProjectionZ(ana::UniqueName().c_str(),
180  binx,binx,biny,biny);
181  hHolder->Add(hCV,-1);
182  hResults.push_back(hHolder);
183  }
184 
185  //Two Sided Systematics
186  for(int i = 0; i < (int)hSystsUp3D.size(); i++){
187  TH1F* hHolderUp =
188  (TH1F*)hSystsUp3D[i]->ProjectionZ(ana::UniqueName().c_str(),
189  binx,binx,biny,biny);
190  TH1F* hHolderDwn =
191  (TH1F*)hSystsDown3D[i]->ProjectionZ(ana::UniqueName().c_str(),
192  binx,binx,biny,biny);
193 
194  hHolderUp->Add(hCV,-1);
195  hHolderDwn->Add(hCV,-1);
196 
197  for(int bin = 0; bin <= hHolderUp->GetXaxis()->GetNbins();bin++){
198  if(fabs(hHolderDwn->GetBinContent(bin)) >
199  fabs(hHolderUp->GetBinContent(bin)))
200  hHolderUp->SetBinContent(bin,hHolderDwn->GetBinContent(bin));
201  }
202  hResults.push_back((TH1F*)hHolderUp->Clone(ana::UniqueName().c_str()));
203  }
204 
205  return hResults;
206  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
float bin[41]
Definition: plottest35.C:14
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
bool ana::nueccinc::NueCCIncTemplateFitter::CheckIfBinShouldBeUsed ( int  binx,
int  biny 
)
private

Definition at line 141 of file NueCCIncTemplateFitter.cxx.

References bin, om::cout, allTimeWatchdog::endl, hCV3D, NC_Values, NumberTemplateBins, NumuLike_Values, Other_Values, and SignalLike_Values.

Referenced by doFullFit().

142  {
143  //Subtract one to compare with vector entries
144  int signal_region_bin = 0;
146  signal_region_bin =
147  hCV3D->GetZaxis()->FindBin(sliding_signal_region[biny])-1;
148  else
149  signal_region_bin =
150  hCV3D->GetZaxis()->FindBin(signal_region) - 1;
151  std::cout << signal_region_bin << std::endl;
152 
153  float signal_amount = 0;
154  float bkgd_amount = 0;
155  for(int bin = signal_region_bin; bin < NumberTemplateBins; bin++)
156  {
157  signal_amount += SignalLike_Values[bin];
158  bkgd_amount += (NC_Values[bin]+NumuLike_Values[bin]+Other_Values[bin]);
159  }
160 
161  return ( ((signal_amount > 100) && (signal_amount/bkgd_amount > 0.25)) ||
162  ((signal_amount/bkgd_amount > 0.2) && (signal_amount > 300)));
163  }
const bool useSlidingSignalRegion
const float sliding_signal_region[13]
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
bool ana::nueccinc::NueCCIncTemplateFitter::CheckSystsSize ( )
private

Definition at line 96 of file NueCCIncTemplateFitter.cxx.

References hCV3D, hData3D, hSystsDown3D, and hSystsUp3D.

Referenced by doFullFit().

97  {
98  return ((hSystsUp3D.size() == hSystsDown3D.size()) &&
99  (hData3D->GetXaxis()->GetNbins() ==
100  hCV3D->GetXaxis()->GetNbins())
101  &&
102  (hData3D->GetYaxis()->GetNbins() ==
103  hCV3D->GetYaxis()->GetNbins())
104  &&
105  (hData3D->GetZaxis()->GetNbins() ==
106  hCV3D->GetZaxis()->GetNbins()));
107  }
std::vector< TH3F * > ana::nueccinc::NueCCIncTemplateFitter::doFullFit ( )

Definition at line 1122 of file NueCCIncTemplateFitter.cxx.

References make_true_q0q3_plots::caption, om::cerr, CheckIfBinShouldBeUsed(), CheckSystsSize(), hadd_many_files::counter, om::cout, allTimeWatchdog::endl, stan::math::fabs(), fcn2Var(), fcn3Var(), FCNMinuit(), FillCovarianceMatrix(), FillStatisticCovarianceMatrix(), FillValueHolders(), fit_covariance_matrix, stan::math::floor(), stan::math::fmin(), GetAnalysisBinCaption(), GetAnalysisTemplates(), GetXaxis(), GetYaxis(), MECModelEnuComparisons::i, make_pair(), nc_wei, nue_wei, numu_wei, cet::pow(), PropagateFitUncertainty3D(), moon_position_table_new3::second, SetBinWeightsAndErrors2D(), and string.

1123  {
1124  std::vector<TH3F*> badentry;
1125  bool CheckSysts = NueCCIncTemplateFitter::CheckSystsSize();
1126  if(!CheckSysts){
1127  std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
1128  return badentry;
1129  }
1130 
1131  std::vector<TH3F*> hResults =
1133 
1134  //Loop over all analysis bins (x and y axis bins)
1135  for(int xbin = 1; xbin <= hResults[0]->GetXaxis()->GetNbins(); xbin++){
1136  for(int ybin = 1; ybin <= hResults[0]->GetYaxis()->GetNbins(); ybin++){
1137 
1138  std::string caption =
1140 
1142 
1143  bool b_fit_this_bin =
1145 
1146  bool check_high_lo_energy = true;
1147  if(hResults[0]->GetYaxis()->GetBinLowEdge(ybin) == 6.0 ||
1148  hResults[0]->GetXaxis()->GetBinLowEdge(xbin) == 0.75)
1149  check_high_lo_energy = false;
1150 
1151  if(b_fit_this_bin && check_high_lo_energy){
1152  std::cout << "Doing Fit in this bin" << std::endl;
1156 
1157  //Now Perform the Fit
1158  //Fit 1: Get Good Starting Fit Values, Quickly
1159  //Three Parameter fit
1160 
1161  //Prepare Starting Seeds
1162  std::vector<float> nue_starting_points;
1163  std::vector<float> numu_starting_points;
1164  std::vector<float> nc_starting_points;
1165 
1166  for(float fillVal = 0.6; fillVal < 1.8; fillVal+=0.15){
1167  nue_starting_points.push_back(fillVal);
1168  numu_starting_points.push_back(fillVal);
1169  nc_starting_points.push_back(fillVal);
1170  }
1171 
1172 
1173  std::random_shuffle(nue_starting_points.begin(),
1174  nue_starting_points.end());
1175  std::random_shuffle(numu_starting_points.begin(),
1176  numu_starting_points.end());
1177  std::random_shuffle(nc_starting_points.begin(),
1178  nc_starting_points.end());
1179 
1180  //Try 3 Parameter Fit First
1181  std::vector<std::pair<float,std::vector<float>>>
1182  pairFitResultStartPts3Var;
1183 
1184  std::vector<std::vector<double>> vUncertainties;
1185  std::vector<std::vector<double>> vNormalizations;
1186 
1187  std::cout << nue_starting_points.size() <<
1188  "\t" << nc_starting_points.size() << "\t"
1189  << numu_starting_points.size() << std::endl;
1190 
1191  int counter = 0;
1192 
1193  for(auto nuestart : nue_starting_points){
1194  for(auto numustart : numu_starting_points){
1195  for(auto ncstart : nc_starting_points){
1196  counter++;
1197  TMinuit *gMinuit = new TMinuit(3);
1198  gMinuit->SetPrintLevel(-1);
1199  gMinuit->SetFCN(fcn3Var);
1200 
1201  Double_t arglist[4];
1202  Int_t ierflg = 0;
1203  arglist[0] = 1;
1204  gMinuit->mnexcm("SET ERR", arglist,1,ierflg);
1205 
1206  //Define Fit Parameters
1207  //0.15 is the approximate error on the parameter
1208  //0,0 refers to the lower and upper bounds on the parameters
1209  gMinuit->mnparm(0, "Numu-Like Scaling",
1210  numustart, 0.15, 0, 0, ierflg);
1211  gMinuit->mnparm(1, "Nue-Like Scaling",
1212  nuestart, 0.15, 0, 0, ierflg);
1213  gMinuit->mnparm(2, "NC Scaling",
1214  ncstart, 0.15, 0, 0, ierflg);
1215 
1216  arglist[0] = 0; //number of iterations 0 == no limit
1217  arglist[1] = 1; //1 == standard minimization strategy
1218  //SIMPLEX should provide a good starting point the next fits
1219  gMinuit->Command("SIMPLEX");
1220 
1221  //Fit 2: Improved Fit + Better Eror Estimation
1222  arglist[1] = 1; //2 == try to improve minimum (slower)
1223  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
1224 
1225  gMinuit->Command("SIMPLEX");
1226  gMinuit->Command("HESSE");
1227  gMinuit->Command("MIGRAD");
1228  gMinuit->Command("MINOS");
1229 
1230 
1231  //Grab Final Results
1232  double vstart[3];
1233  double verror[3];
1234  double fParNewStart;
1235  double fParNewErr;
1236  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1237  vstart[0] = fParNewStart;
1238  verror[0] = fParNewErr;
1239  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1240  vstart[1] = fParNewStart;
1241  verror[1] = fParNewErr;
1242  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1243  vstart[2] = fParNewStart;
1244  verror[2] = fParNewErr;
1245 
1246  //Before Ending Fits,
1247  //Make Sure There Was No Issue With the Fit
1248  bool troubleFitting = false;
1249  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1250  Int_t npari =0, nparx = 0, istat = 0;
1251  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1252 
1253  Double_t corr_mat_test[3][3];
1254  gMinuit->mnemat(*corr_mat_test,3);
1255 
1256  //Turn this test error matix into a correlation matrix
1257  double corr_01 =
1258  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1259  pow(corr_mat_test[1][1],0.5));
1260  double corr_02=
1261  corr_mat_test[0][2]/((pow(corr_mat_test[0][0],0.5) *
1262  pow(corr_mat_test[2][2],0.5)));
1263  double corr_12=
1264  corr_mat_test[1][2]/((pow(corr_mat_test[1][1],0.5) *
1265  pow(corr_mat_test[2][2],0.5)));
1266 
1267  if(istat != 3 ||
1268  vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1269  vstart[0] > 2 || vstart[1] > 2 || vstart[2] > 2 ||
1270  fabs(corr_01) > 0.97 || fabs(corr_02) > 0.97 ||
1271  fabs(corr_12) > 0.97)
1272  troubleFitting = true;
1273 
1274  std::vector<double> vUncert = {verror[0]/vstart[0],
1275  verror[1]/vstart[1],
1276  verror[2]/vstart[2]};
1277  std::vector<double> vNormal = {vstart[0],vstart[1],vstart[2]};
1278 
1279  std::vector<float> vStarts = {nuestart,numustart,ncstart};
1280 
1281  if(!troubleFitting){
1282  pairFitResultStartPts3Var.push_back(std::make_pair
1283  (fmin,vStarts));
1284  vUncertainties.push_back(vUncert);
1285  vNormalizations.push_back(vNormal);
1286  }
1287 
1288  } //loop over nc starting values
1289  } //loop over numu starting values
1290  }//loop over nue starting values
1291 
1292  if(printVerbose){
1293  for(uint i = 0; i < pairFitResultStartPts3Var.size(); i++){
1294  std::cout << pairFitResultStartPts3Var[i].first << "\t"
1295  << vNormalizations[i][0] << "\t"
1296  << vUncertainties[i][0] << "\t"
1297  << vNormalizations[i][1] << "\t"
1298  << vUncertainties[i][1] << "\t"
1299  << vNormalizations[i][0] << "\t"
1300  << vUncertainties[i][2] << "\t" << std::endl;
1301  }
1302  }
1303 
1304  bool shouldSwitchTo2VarFit = false;
1305  //Find the best fits
1306  int saveIndex = 0;
1307  if(pairFitResultStartPts3Var.size() > 0){
1308  //now we can loop through all of the items in the vector
1309  for(uint i = 0; i < pairFitResultStartPts3Var.size(); i++)
1310  {
1311  if(floor(pairFitResultStartPts3Var[saveIndex].first*1e3) >
1312  floor(pairFitResultStartPts3Var[i].first*1e3))
1313  {
1314  saveIndex = i;
1315  chisquared = pairFitResultStartPts3Var[saveIndex].first;
1316  }
1317  if(floor(pairFitResultStartPts3Var[saveIndex].first*1e3) ==
1318  floor(pairFitResultStartPts3Var[i].first*1e3))
1319  {
1320  if(vUncertainties[i][0] <
1321  vUncertainties[saveIndex][0] &&
1322  vUncertainties[i][1] <
1323  vUncertainties[saveIndex][1] &&
1324  vUncertainties[i][2] <
1325  vUncertainties[saveIndex][2])
1326  {
1327  saveIndex = i;
1328  chisquared = pairFitResultStartPts3Var[saveIndex].first;
1329  }
1330  }
1331  }
1332  }
1333  else shouldSwitchTo2VarFit = true;
1334 
1335 
1336  //Put 3 Var Fit Here
1337  if(!shouldSwitchTo2VarFit){
1338  //Fit 1: Get Good Starting Fit Values, Quickly
1339  //Three Parameter fit
1340  TMinuit *gMinuit = new TMinuit(3);
1341  gMinuit->SetPrintLevel(-1);
1342  if(printVerbose)gMinuit->SetPrintLevel(1);
1343  gMinuit->SetFCN(fcn3Var);
1344  Double_t arglist[4];
1345  Int_t ierflg = 0;
1346  arglist[0] = 1;
1347  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1348 
1349 
1350  double vstart[3] = {0.8,0.9,1.15};
1351  double verror[3] = {0,0,0};
1352  double vstep[3] = {0.15,0.15,0.15};
1353 
1354  //Define Fit Parameters
1355  gMinuit->mnparm(0, "Numu-Like Scaling",
1356  pairFitResultStartPts3Var[saveIndex].second[1],
1357  vstep[0], 0, 0, ierflg);
1358  gMinuit->mnparm(1, "Nue-Like Scaling",
1359  pairFitResultStartPts3Var[saveIndex].second[0],
1360  vstep[1], 0, 0, ierflg);
1361  gMinuit->mnparm(2, "NC Scaling",
1362  pairFitResultStartPts3Var[saveIndex].second[2],
1363  vstep[2], 0, 0, ierflg);
1364 
1365  arglist[0] = 0; //number of iterations 0 == no limit
1366  arglist[1] = 1; //1 == standard minimization strategy
1367 
1368  //SIMPLEX should provide a good starting point for further fits
1369  gMinuit->Command("SIMPLEX");
1370 
1371  //Fit 2: Improved Fit + Better Eror Estimation
1372  arglist[1] = 2; //2 == try to improve minimum (slower)
1373  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
1374  gMinuit->Command("SIMPLEX");
1375 
1376  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1377  Int_t npari =0, nparx = 0, istat = 0;
1378  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1379 
1380  gMinuit->Command("SIMPLEX");
1381  gMinuit->Command("HESSE");
1382  gMinuit->Command("MIGRAD");
1383  gMinuit->Command("MINOS");
1384 
1385  //Grab Final Results
1386  double fParNewStart;
1387  double fParNewErr;
1388  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1389  vstart[0] = fParNewStart;
1390  verror[0] = fParNewErr;
1391  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1392  vstart[1] = fParNewStart;
1393  verror[1] = fParNewErr;
1394  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1395  vstart[2] = fParNewStart;
1396  verror[2] = fParNewErr;
1397 
1398  //Before Ending Fits, Make Sure There Was No Issue With the Fit
1399  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1400 
1401  Double_t corr_mat_test[3][3];
1402  gMinuit->mnemat(*corr_mat_test,3);
1403 
1404  //Turn this test error matix into a correlation matrix
1405  double corr_01 =
1406  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1407  pow(corr_mat_test[1][1],0.5));
1408  double corr_02=
1409  corr_mat_test[0][2]/((pow(corr_mat_test[0][0],0.5) *
1410  pow(corr_mat_test[2][2],0.5)));
1411  double corr_12=
1412  corr_mat_test[1][2]/((pow(corr_mat_test[1][1],0.5) *
1413  pow(corr_mat_test[2][2],0.5)));
1414 
1415 
1416  if(istat != 3 ||
1417  vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1418  vstart[0] > 2. || vstart[1] > 2 || vstart[2] > 2.0 ||
1419  fabs(corr_01) > 0.98 || fabs(corr_02) > 0.98 ||
1420  fabs(corr_12) > 0.98)
1421  shouldSwitchTo2VarFit = true;
1422 
1423  //Get the covariance matrix
1424  Double_t emat[3][3];
1425  gMinuit->mnemat(*emat,3);
1426 
1427  nue_wei = std::make_pair (vstart[1],verror[1]);
1428  numu_wei = std::make_pair (vstart[0],verror[0]);
1429  nc_wei = std::make_pair (vstart[2],verror[2]);
1430  fit_covariance_matrix[0][0] = emat[0][0];
1431  fit_covariance_matrix[0][1] = emat[0][1];
1432  fit_covariance_matrix[0][2] = emat[0][2];
1433  fit_covariance_matrix[1][0] = emat[1][0];
1434  fit_covariance_matrix[1][1] = emat[1][1];
1435  fit_covariance_matrix[1][2] = emat[1][2];
1436  fit_covariance_matrix[2][0] = emat[2][0];
1437  fit_covariance_matrix[2][1] = emat[2][1];
1438  fit_covariance_matrix[2][2] = emat[2][2];
1439  std::cout << shouldSwitchTo2VarFit <<
1440  "\tMade it Here" << std::endl;
1441  }
1442 
1443  //Now If there were issues with the 3 Parameter Fit,
1444  //switch to a 2 ParameterFit
1445  std::vector<std::pair<float,std::vector<float>>>
1446  pairFitResultStartPts2Var;
1447 
1448  std::vector<std::vector<double>> vUncertainties2Var;
1449  std::vector<std::vector<double>> vNormalizations2Var;
1450 
1451  std::vector<float> signal_starting_points;
1452  std::vector<float> background_starting_points;
1453 
1454  for(float fillVal = 0.6; fillVal < 1.8; fillVal+=0.15){
1455  signal_starting_points.push_back(fillVal);
1456  background_starting_points.push_back(fillVal);
1457  }
1458 
1459  if(shouldSwitchTo2VarFit){
1460  //There Was trouble with the original fits,
1461  //most likely this was due to extremely
1462  //high correlation (approx 1)
1463  //between two of the fitting parameters
1464  //Switch to a fit that adjusts all backgrounds
1465  for(auto nuestart : signal_starting_points){
1466  for(auto bkgstart : background_starting_points){
1467  bool isGoodResult = true;
1468  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1469  Int_t npari =0, nparx = 0, istat = 0;
1470  Double_t arglist[4];
1471  Int_t ierflg = 0;
1472  TMinuit *jMinuit = new TMinuit(2);
1473  jMinuit->SetPrintLevel(-1);
1474  jMinuit->SetFCN(fcn2Var);
1475  arglist[0] = 1;
1476  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1477 
1478  jMinuit->mnparm(0, "Background Scaling",bkgstart,
1479  0.15, 0,0, ierflg);
1480  jMinuit->mnparm(1, "Nue CC Scaling",
1481  nuestart, 0.15, 0,0, ierflg);
1482 
1483  arglist[0] = 0; //number of iterations 0 == no limit
1484  arglist[1] = 1;
1485  jMinuit->Command("SIMPLEX");
1486  jMinuit->Command("HESSE");
1487 
1488  //Fit 2: Improved Fit + Better Eror Estimation
1489  arglist[1] = 2; //2 == try to improve minimum (slower)
1490  jMinuit->mnexcm("SET STR", arglist,1,ierflg);
1491 
1492  jMinuit->Command("SIMPLEX");
1493  jMinuit->Command("HESSE");
1494  jMinuit->Command("MIGRAD");
1495  //jMinuit->Command("MINOS");
1496 
1497 
1498  //Before Ending Fits,Make Sure There Was No Issue With the Fit
1499  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1500 
1501 
1502  //errors are
1503  // ISTAT:a status integerindicating how good is the covariance
1504  //*-* matrix:0= not calculated at all
1505  //*-* 1= approximation only, not accurate
1506  //*-* 2= full matrix, but forced positive-def
1507  //*-* 3= full accurate covariance matrix
1508  if(istat != 3)
1509  isGoodResult = false;
1510 
1511 
1512  //Check for high correlations between fit parameters,
1513  //can show an issue with the fit in this bin
1514  Double_t corr_mat_test[2][2];
1515  jMinuit->mnemat(*corr_mat_test,2);
1516 
1517 
1518  //Turn this test error matix into a correlation matrix
1519  double corr_01 =
1520  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1521  pow(corr_mat_test[1][1],0.5));
1522 
1523  if(fabs(corr_01) > 0.97)
1524  isGoodResult = false;
1525 
1526 
1527 
1528  std::vector<double> vUncert;
1529  std::vector<double> vNormal;
1530  for(int i = 0; i < 2; i++){
1531  double fParNormalization,fParError;
1532  jMinuit->GetParameter(i,fParNormalization,fParError);
1533  vUncert.push_back(fParError);
1534  vNormal.push_back(fParNormalization);
1535  if(fParNormalization < 0 || fParNormalization > 2 ||
1536  fParError/fParNormalization > 1)
1537  isGoodResult = false;
1538  if(!isGoodResult)break;
1539  }
1540 
1541  std::vector<float> vStart = {nuestart,bkgstart};
1542 
1543  if(isGoodResult)
1544  {pairFitResultStartPts2Var.push_back(std::make_pair
1545  (fmin,vStart));
1546  vUncertainties2Var.push_back(vUncert);
1547  vNormalizations2Var.push_back(vNormal);
1548 
1549  }
1550 
1551 
1552  }//background starting values loop
1553  } //nue starting values loop
1554  }//do 2Var Fit
1555 
1556  if(printVerbose){
1557  for(uint i = 0; i < pairFitResultStartPts2Var.size(); i++){
1558  std::cout << pairFitResultStartPts2Var[i].first << "\t"
1559  << vNormalizations2Var[i][0] << "\t"
1560  << vUncertainties2Var[i][0] << "\t"
1561  << vNormalizations2Var[i][1] << "\t"
1562  << vUncertainties2Var[i][1] << std::endl;
1563  }
1564  }
1565 
1566  //Grab best fit by ChiSquared Value
1567  //If those are equal sort by predicted uncertainties
1568  //All Uncertainties need to be better to be considered the best fit
1569 
1570  TH1F* hUncertaintyAverage = new TH1F("hUncertaintyAverage",
1571  ";Signal Fit Uncertainty;Count",
1572  100,0,2.0);
1573 
1574  saveIndex = 0;
1575 
1576 
1577  if(shouldSwitchTo2VarFit){
1578  for(uint i = 0; i < pairFitResultStartPts2Var.size(); i++)
1579  {
1580  hUncertaintyAverage->Fill(vUncertainties2Var[i][1]);
1581  if(floor(pairFitResultStartPts2Var[saveIndex].first*1e5) >
1582  floor(pairFitResultStartPts2Var[i].first*1e5))
1583  {
1584  saveIndex = i;
1585  chisquared = pairFitResultStartPts2Var[saveIndex].first;
1586  shouldSwitchTo2VarFit = true;
1587  }
1588  }
1589  }
1590 
1591  std::cout << "SaveIndex " << saveIndex
1592  << "************" << std::endl;
1593  std::cout << pairFitResultStartPts2Var.size() << std::endl;
1594 
1595 
1596  if(shouldSwitchTo2VarFit){
1597  //There Was trouble with the original fits,
1598  //most likely this was due to extremely high correlation (approx 1)
1599  //between two of the fitting parameters
1600  //Switch to a fit that adjusts all backgrounds
1601  double vstart[3] = {0.8,0.9,1.15};
1602  double verror[3] = {0,0,0};
1603  double vstep[3] = {0.15,0.15,0.15};//{0.01,0.01,0.01
1604 
1605  Double_t arglist[4];
1606  Int_t ierflg = 0;
1607  TMinuit *jMinuit = new TMinuit(2);
1608  jMinuit->SetPrintLevel(-1);
1609  if(printVerbose)jMinuit->SetPrintLevel(1);
1610  jMinuit->SetFCN(fcn2Var);
1611  arglist[0] = 1;
1612  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1613 
1614  double nue_start = 1.0;
1615  double bkgd_start = 1.0;
1616 
1617  if(pairFitResultStartPts2Var.size() > 0){
1618  nue_start = pairFitResultStartPts2Var[saveIndex].second[1];
1619  bkgd_start = pairFitResultStartPts2Var[saveIndex].second[0];
1620  std::cout << "Made it here" << std::endl;
1621  std::cout << nue_start << "\t" << bkgd_start << std::endl;
1622  }
1623 
1624 
1625  jMinuit->mnparm(0, "Background Scaling",
1626  bkgd_start,
1627  vstep[0], 0, 0, ierflg);
1628  jMinuit->mnparm(1, "Nue CC Scaling",
1629  nue_start,
1630  vstep[1], 0, 0, ierflg);
1631 
1632  arglist[0] = 0; //number of iterations 0 == no limit
1633  arglist[1] = 1;
1634 
1635  jMinuit->Command("SIMPLEX");
1636  jMinuit->Command("HESSE");
1637 
1638  //Fit 2: Improved Fit + Better Eror Estimation
1639  arglist[1] = 2; //2 == try to improve minimum (slower)
1640  jMinuit->mnexcm("SET STR", arglist,1,ierflg);
1641 
1642  jMinuit->Command("SIMPLEX");
1643  for(int i = 0; i < 2; i++){
1644  double fParNewStart = 0;
1645  double fParNewErr = 0;
1646  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1647  std::cout << "SIMPLEX" << "\t" << fParNewStart << "\t"
1648  << fParNewErr << std::endl;
1649  }
1650  jMinuit->Command("MIGRAD");
1651  for(int i = 0; i < 2; i++){
1652  double fParNewStart = 0;
1653  double fParNewErr = 0;
1654  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1655  std::cout << "MIGRAD" << "\t" << fParNewStart << "\t"
1656  << fParNewErr << std::endl;
1657  }
1658  //jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
1659  jMinuit->Command("HESSE");
1660  for(int i = 0; i < 2; i++){
1661  double fParNewStart = 0;
1662  double fParNewErr = 0;
1663  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1664  std::cout << "SEEK" << "\t" << fParNewStart << "\t"
1665  << fParNewErr << std::endl;
1666  }
1667  //jMinuit->Command("MINOS");
1668  for(int i = 0; i < 2; i++){
1669  double fParNewStart = 0;
1670  double fParNewErr = 0;
1671  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1672  std::cout << "MINOS" << "\t" << fParNewStart << "\t"
1673  << fParNewErr << std::endl;
1674  }
1675 
1676  double fParNewStart;
1677  double fParNewErr;
1678  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1679  vstart[0] = fParNewStart;
1680  verror[0] = fParNewErr;
1681  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1682  vstart[2] = fParNewStart;
1683  verror[2] = fParNewErr;
1684  jMinuit->GetParameter(1,fParNewStart,fParNewErr);
1685  vstart[1] = fParNewStart;
1686  verror[1] = fParNewErr;
1687  //Get the covariance matrix
1688  Double_t emat[2][2];
1689  jMinuit->mnemat(*emat,2);
1690 
1691  nue_wei = std::make_pair (vstart[1],verror[1]);
1692  numu_wei = std::make_pair (vstart[0],verror[0]);
1693  nc_wei = std::make_pair (vstart[2],verror[2]);
1694  fit_covariance_matrix[0][0] = emat[0][0];
1695  fit_covariance_matrix[0][1] = emat[0][1];
1696  fit_covariance_matrix[0][2] = 0;
1697  fit_covariance_matrix[1][0] = emat[1][0];
1698  fit_covariance_matrix[1][1] = emat[1][1];
1699  fit_covariance_matrix[1][2] = 0;
1700  fit_covariance_matrix[2][0] = emat[0][0];
1701  fit_covariance_matrix[2][1] = emat[0][1];
1702  fit_covariance_matrix[2][2] = emat[0][0];
1703  std::cout << shouldSwitchTo2VarFit <<
1704  "\tMade it Here, 2Var Fit true" << std::endl;
1705  }//2var results
1706 
1707 
1708  }//fit this bin
1709  else{
1710  std::cout << "Skipping Fit in " << caption << std::endl;
1711 
1712  nue_wei = std::make_pair (1,0);
1713  numu_wei = std::make_pair (1,0);
1714  nc_wei = std::make_pair (1,0);
1715  chisquared = 0;
1716  fit_covariance_matrix[0][0] = 0;
1717  fit_covariance_matrix[0][1] = 0;
1718  fit_covariance_matrix[0][2] = 0;
1719  fit_covariance_matrix[1][0] = 0;
1720  fit_covariance_matrix[1][1] = 0;
1721  fit_covariance_matrix[1][2] = 0;
1722  fit_covariance_matrix[2][0] = 0;
1723  fit_covariance_matrix[2][1] = 0;
1724  fit_covariance_matrix[2][2] = 0;
1725  }
1726 
1727  std::cout << "*********Fit for bin " << caption << " : " << std::endl;
1728  std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
1729  numu_wei.second << std::endl;
1730  std::cout << "NueCC Scaling Factor= " << nue_wei.first << " +- " <<
1731  nue_wei.second << std::endl;
1732  std::cout << "NC Scaling Factor= " << nc_wei.first << " +- " <<
1733  nc_wei.second << std::endl;
1734 
1737  }//loop over ybins
1738  }//loop over xbins
1739 
1740  //Now That we have gone over every Bin and propagated uncertainty
1741  //Return the final results
1742  std::vector<TH3F*> hResultsPostFit =
1744 
1745  for(uint i = 0; i < hResults.size(); i++){
1746  std::cout << i << "\t" << hResults[i]->Integral() << "\t"
1747  << hResultsPostFit[i]->Integral() << std::endl;
1748  }
1749 
1750 
1751  return hResultsPostFit;
1752  }
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
correl_xv GetYaxis() -> SetDecimals()
constexpr T pow(T x)
Definition: pow.h:75
static void fcn3Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
OStream cerr
Definition: OStream.cxx:7
std::string GetAnalysisBinCaption(int binx, int biny)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
correl_xv GetXaxis() -> SetDecimals()
void FillCovarianceMatrix(int binx, int biny, bool makePlots)
OStream cout
Definition: OStream.cxx:6
void FillStatisticCovarianceMatrix(int binx, int biny)
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
static void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
enum BeamMode string
unsigned int uint
void ana::nueccinc::NueCCIncTemplateFitter::fcn2Var ( Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag 
)
staticprivate

Definition at line 774 of file NueCCIncTemplateFitter.cxx.

References myFunction2Vars().

Referenced by doFullFit().

777  {
778  f = Fitter_obj->myFunction2Vars(par[0],par[1]);
779  }
Int_t par
Definition: SimpleIterate.C:24
NueCCIncTemplateFitter * Fitter_obj
double myFunction2Vars(double bkgd_par, double nue_par)
void ana::nueccinc::NueCCIncTemplateFitter::fcn3Var ( Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag 
)
staticprivate

Definition at line 711 of file NueCCIncTemplateFitter.cxx.

References myFunction3Vars().

Referenced by doFullFit().

714  {
715  f = Fitter_obj->myFunction3Vars(par[0],par[1],par[2]);
716  }
Int_t par
Definition: SimpleIterate.C:24
double myFunction3Vars(double numu_par, double nue_par, double nc_par)
NueCCIncTemplateFitter * Fitter_obj
void ana::nueccinc::NueCCIncTemplateFitter::FCNMinuit ( )

Definition at line 1103 of file NueCCIncTemplateFitter.cxx.

Referenced by doFullFit().

1104  {
1105  Fitter_obj = this;
1106  }
NueCCIncTemplateFitter * Fitter_obj
TH2F * ana::nueccinc::NueCCIncTemplateFitter::FillCovarianceBetweenSysts ( int  binx,
int  biny 
)
private

Definition at line 1004 of file NueCCIncTemplateFitter.cxx.

References CalculateOneSigmaShift(), hadd_many_files::counter, covariance(), hCV3D, hSystsMultiverse, MECModelEnuComparisons::i, makeTrainCVSamples::int, febshutoff_auto::integral, calib::j, make_pair(), gen_hdf5record::names, and ana::UniqueName().

Referenced by getCorrelationBetweenSysts(), and getCovarianceBetweenSysts().

1005  {
1006  //Make Sure this bin isn't excluded by the phase space cuts
1007  //If it is, don't worrky about calculating anything
1008  //Just checking one of the inner bins
1009  TH2F* hResultBad = new TH2F(ana::UniqueName().c_str(),";",1,0,1,1,0,1);
1010  if(hCV3D->Integral(binx,binx,biny,biny,1,hCV3D->GetZaxis()->GetNbins()) ==
1011  0) return hResultBad;
1012  //Systematic Samples
1013  //Turn Systamtic samples that are not generated from a multiverse
1014  //into a "one sigma" shift that we can randomly draw from in a moment
1015  std::vector<TH1F*> hSystsShifts =
1017 
1018  std::vector<std::pair<std::string,std::vector<double>>> hSysts_Pairs;
1019 
1020  std::vector<std::string> names = {"Cherenkov","Cal. Shape",
1021  "Light","Cal. Norm.",
1022  "Horn Current", "Beam Size",
1023  "Beam Shift X", "Beam Shift Y",
1024  "Horn 1 Shift X", "Horn 1 Shift Y",
1025  "Horn 2 Shift X", "Horn 2 Shift Y",
1026  "Target Z Shift", "Magnetic Field",
1027  "Horn Wather Thickness",
1028  "PPFX","GENIE"};
1029 
1030  TRandom *r0 = new TRandom();
1031 
1032  //Randomly draw from these one sigma shift samples
1033  //Draw 1000 from each sample
1034  for(int j = 0; j < (int)hSystsShifts.size(); j++){
1035  std::vector<double> vPseudoExp;
1036  for(int i = 0; i < 1000; i++){
1037  TH1F* hClone = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
1038  binx,binx,biny,biny);
1039  double wei = r0->Gaus(0,1);
1040  hClone->Add(hSystsShifts[j],wei);
1041  vPseudoExp.push_back(hClone->Integral(1,-1));
1042  delete hClone;
1043  }//loop over systematic samples
1044  hSysts_Pairs.push_back(std::make_pair(names[j],vPseudoExp));
1045  }//loop over pseudo experiments
1046 
1047  std::vector<double> vPPFX_Values;
1048  std::vector<double> vGenie_Values;
1049 
1050  for(int systNum = 0; systNum < (int)hSystsMultiverse.size();systNum++){
1051  double integral =
1052  hSystsMultiverse[systNum]->Integral(binx,binx,biny,biny,1,-1);
1053  if(systNum < 300)vGenie_Values.push_back(integral);
1054  else vPPFX_Values.push_back(integral);
1055  }//loop over multiverse universes
1056 
1057  hSysts_Pairs.push_back(std::make_pair("ppfx", vPPFX_Values));
1058  hSysts_Pairs.push_back(std::make_pair("genie", vGenie_Values));
1059 
1060  //Calculate Covariance
1061  //Multiverses are done first
1062  //Followed by the "multiverse" we created from the other shifts
1063  //above
1064  const int systNum = (int)hSysts_Pairs.size();
1065  double syst_covariance[systNum][systNum];
1066  for(int i = 0; i < systNum; i++){
1067  for(int j = 0; j < systNum; j++){
1068  double numerator_term = 0;
1069  double counter = hSysts_Pairs[i].second.size();
1070  for(int nUniv = 0; nUniv < counter; nUniv++){
1071  double term_i = hSysts_Pairs[i].second[nUniv] -
1072  hCV3D->Integral(binx,binx,biny,biny,1,-1);
1073  double term_j = hSysts_Pairs[j].second[nUniv] -
1074  hCV3D->Integral(binx,binx,biny,biny,1,-1);
1075  numerator_term += term_i*term_j * (1.0/(counter - 1));
1076  }
1077  double covariance = numerator_term;
1078  syst_covariance[i][j] = covariance;
1079  }
1080  }
1081 
1082 
1083  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),";",
1084  systNum,0,systNum,systNum,0,systNum);
1085 
1086  //Move the covariance matrix to the final result
1087  for(int i = 0; i < systNum; i++){
1088  for(int j = 0; j < systNum; j++){
1089  hResult->SetBinContent(i+1,j+1,syst_covariance[i][j]);
1090  hResult->GetXaxis()->SetBinLabel(i+1,names[i].c_str());
1091  hResult->GetYaxis()->SetBinLabel(j+1,names[j].c_str());
1092  }
1093  }
1094 
1095  hResult->GetZaxis()->SetTitle("Covariance (Events^{2})");
1096 
1097  return hResult;
1098  }
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< TH1F * > CalculateOneSigmaShift(int binx, int biny)
const double j
Definition: BetheBloch.cxx:29
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void ana::nueccinc::NueCCIncTemplateFitter::FillCovarianceMatrix ( int  binx,
int  biny,
bool  makePlots = false 
)
private

Definition at line 208 of file NueCCIncTemplateFitter.cxx.

References ana::AutoPlaceLegend(), demo5::c1, chisquared::c3, CalculateOneSigmaShift(), plot_validation_datamc::Clone(), covariance(), CovarianceMatrix, GetAnalysisBinCaption(), ana::MultiverseCorrelation::GetCorrelationMatrix(), ana::MultiverseCorrelation::GetCovarianceMatrix(), hCV3D, hSystsDown3D, hSystsMultiverse, hSystsOneSided, hSystsUp3D, MECModelEnuComparisons::i, std::isnan(), calib::j, kBlue, kRed, MECModelEnuComparisons::leg, make_pair(), makePlots(), NumberTemplateBins, ana::nueccinc::numPIDBins, output, std::sqrt(), string, ana::UniqueName(), and ymax.

Referenced by doFullFit(), getCorrelationMatrixTemplateBins(), and getCovarianceMatrixTemplateBins().

210  {
211  //Make Sure this bin isn't excluded by the phase space cuts
212  //If it is, don't worrky about calculating anything
213  //Just checking one of the inner bins
214  if(hCV3D->Integral(binx,binx,biny,biny,1,hCV3D->GetZaxis()->GetNbins()) ==
215  0) return;
216 
217  std::string binTitle =
219 
220 
221  //Systematic Samples
222  //Turn Systamtic samples that are not generated from a multiverse
223  //into a "one sigma" shift that we can randomly draw from in a moment
224  std::vector<TH1F*> hSystsShifts =
226 
227  std::map<std::string,TH2F*> mapCovariances;
228  std::map<std::string,TH2F*> mapCorrelations;
229  //Order of the systematics is Ckv, Calibration Shape, Light Levels, Cal
230  std::vector<std::string> vSystNames = {"ckv", "calib_shape", "light",
231  "calib", //"horn_current",
232  //"spot_size", "shiftx",
233  //"shifty", "horn1x", "horn1y",
234  //"horn2x", "horn2y", "targetz",
235  //"mag_field", "hornwater",
236  "multiverse"};
237 
238  //TRandom *r0 = new TRandom();
239 
240  int singlesided_counter = 0;
241  int doublesided_counter = 0;
242  for(uint systNum = 0; systNum < hSystsShifts.size(); systNum++){
243  //First draw systematic around the cv for the template
244  if((vSystNames[systNum] == "ckv" ||
245  vSystNames[systNum] == "calib_shape") && makePlots){
246  hSystsOneSided[singlesided_counter]->SetLineColorAlpha(kRed-9,0.5);
247  TH1F* hSyst =
248  (TH1F*)hSystsOneSided[singlesided_counter]->
249  ProjectionZ(ana::UniqueName().c_str(),
250  binx,binx,biny,biny);
251  TH1F* hCV1D =
252  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
253  binx,binx,biny,biny);
254  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
255  hCV1D->GetXaxis()->SetTitle("Template Bins");
256  hCV1D->SetTitle(binTitle.c_str());
257  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
258  c1->SetLeftMargin(0.15);
259  c1->SetBottomMargin(0.15);
260  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
261  hCV1D->Draw("hist");
262  hSyst->Draw("hist same");
263  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
264  leg->AddEntry(hSyst, "Sigma Shift", "fl");
265  leg->AddEntry(hCV1D, "Central Value", "fl");
266  leg->Draw();
267  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
268  "standard_template", vSystNames[systNum].c_str(),
269  binx,biny));
270  c1->Close();
271  delete leg;
272  delete c1;
273  singlesided_counter++;
274  }
275  else if(makePlots){
276  hSystsUp3D[doublesided_counter]->SetLineColorAlpha(kRed-9,0.5);
277  hSystsDown3D[doublesided_counter]->SetLineColorAlpha(kBlue-9,0.5);
278  TH1F* hSystsUp =
279  (TH1F*)hSystsUp3D[doublesided_counter]->
280  ProjectionZ(ana::UniqueName().c_str(),
281  binx,binx,biny,biny);
282  TH1F* hSystsDown =
283  (TH1F*)hSystsDown3D[doublesided_counter]->
284  ProjectionZ(ana::UniqueName().c_str(),
285  binx,binx,biny,biny);
286  TH1F* hCV1D =
287  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
288  binx,binx,biny,biny);
289  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
290  hCV1D->GetXaxis()->SetTitle("Template Bins");
291  hCV1D->SetTitle(binTitle.c_str());
292  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
293  c1->SetLeftMargin(0.15);
294  c1->SetBottomMargin(0.15);
295  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
296  hCV1D->Draw("hist");
297  hSystsUp->Draw("hist same");
298  hSystsDown->Draw("hist same");
299  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
300  leg->AddEntry(hSystsUp, "Sigma Up", "fl");
301  leg->AddEntry(hSystsDown, "Sigma Down", "fl");
302  leg->AddEntry(hCV1D, "Central Value", "fl");
303  leg->Draw();
304  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
305  "standard_template", vSystNames[systNum].c_str(),
306  binx,biny));
307  c1->Close();
308  delete leg;
309  delete c1;
310  doublesided_counter++;
311  }
312 
313  //Randomly draw from the one sigma shift samples
314  //Draw 1000 from each sample
315  std::vector<TH1F*> vPseudoExp;
316  for(uint i = 0; i < 1; i++){
317  TH1F* hClone = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
318  binx,binx,biny,biny);
319  double wei = 1;
320  hClone->Add(hSystsShifts[systNum],wei);
321  vPseudoExp.push_back(hClone);
322  }//loop over pseudo experiments
323 
324 
325  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
326  binx,binx,biny,biny);
327  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
328  hCV1D->GetXaxis()->SetTitle("Template Bins");
329  hCV1D->GetYaxis()->CenterTitle();
330  hCV1D->GetXaxis()->CenterTitle();
331  hCV1D->SetTitle(binTitle.c_str());
332 
333 
334  //Draw the Pseudo experiments
335  if(makePlots){
336  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
337  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
338  hCV1D->Draw("hist");
339  for(uint i = 0; i < vPseudoExp.size(); i++){
340  vPseudoExp[i]->SetLineColorAlpha(kBlue-9,0.50);
341  vPseudoExp[i]->Draw("hist same");
342  }
343  hCV1D->Draw("hist same");
344  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
345  "standard_template_pseudo_exp",
346  vSystNames[systNum].c_str(),binx,biny));
347  c1->Close();
348  }
349 
350  TH2F* hCovHolder =
351  new TH2F(ana::UniqueName().c_str(),
352  ";Template Bins;Template Bins; Covariance (Events^{2})",
355  TH2F* hCorHolder =
356  new TH2F(ana::UniqueName().c_str(),
357  ";Template Bins;Template Bins; Correlation",
360  hCovHolder->SetTitle(binTitle.c_str());
361  hCorHolder->SetTitle(binTitle.c_str());
362 
363  //Calculate Covariance Matrix For the Single Systematic
364  for(uint iu = 0; iu < vPseudoExp.size(); iu++){
365  for(int i = 0; i < numPIDBins; i++){
366  double xi=vPseudoExp[iu]->GetBinContent(i+1);
367  double ximean=hCV1D->GetBinContent(i+1);
368  for(int j = i; j < numPIDBins; j++){
369  double xj=vPseudoExp[iu]->GetBinContent(j+1);
370  double xjmean=hCV1D->GetBinContent(j+1);
371  double current_value = hCovHolder->GetBinContent(i+1,j+1);
372  double this_value = (xi-ximean)*(xj-xjmean);
373  hCovHolder->SetBinContent(i+1,j+1,current_value+this_value);
374  }
375  }
376 
377  }
378 
379  //Only calculate with the 1 sigma error band
380  const double N=(double)2;
381  hCovHolder->Scale(1.0/(N-1.0));
382  for(int i=0; i<numPIDBins; i++){
383  for(int j=i; j<numPIDBins; j++){
384  hCovHolder->SetBinContent(j+1,i+1,
385  hCovHolder->GetBinContent(i+1,j+1));
386  }
387  }
388 
389  for(int i=0; i<numPIDBins; i++){
390  for(int j=0; j<numPIDBins; j++){
391  double covariance = hCovHolder->GetBinContent(i+1,j+1);
392  double variance_i = sqrt(hCovHolder->GetBinContent(i+1,i+1));
393  double variance_j = sqrt(hCovHolder->GetBinContent(j+1,j+1));
394  double correlation = covariance/(variance_i*variance_j);
395  if(isnan(correlation)) correlation = 0;
396  hCorHolder->SetBinContent(i+1,j+1,
397  correlation);
398  }
399  }
400 
401 
402  mapCovariances.insert(std::make_pair(vSystNames[systNum],
403  hCovHolder));
404  mapCorrelations.insert(std::make_pair(vSystNames[systNum],
405  hCorHolder));
406  }
407 
408 
409  //Now For the Multiverse
410  //GENIE and Shape Only PPFX
411  {
412 
413  TH2F* hCovHolder =
414  new TH2F(ana::UniqueName().c_str(),
415  ";Template Bins;Template Bins; Covariance (Events^{2})",
416  numPIDBins,0,numPIDBins,
417  numPIDBins,0,numPIDBins);
418  TH2F* hCorHolder =
419  new TH2F(ana::UniqueName().c_str(),
420  ";Template Bins;Template Bins;Correlation",
421  numPIDBins,0,numPIDBins,
422  numPIDBins,0,numPIDBins);
423 
424  hCovHolder->SetTitle(binTitle.c_str());
425  hCorHolder->SetTitle(binTitle.c_str());
426 
427  TH2F* hCovHolderPPFX =
428  new TH2F(ana::UniqueName().c_str(),
429  ";Template Bins;Template Bins; Covariance (Events^{2})",
430  numPIDBins,0,numPIDBins,
431  numPIDBins,0,numPIDBins);
432  TH2F* hCorHolderPPFX =
433  new TH2F(ana::UniqueName().c_str(),
434  ";Template Bins;Template Bins;Correlation",
435  numPIDBins,0,numPIDBins,
436  numPIDBins,0,numPIDBins);
437 
438  hCovHolderPPFX->SetTitle(binTitle.c_str());
439  hCorHolderPPFX->SetTitle(binTitle.c_str());
440 
441  TH1D* hCV_Tester = (TH1D*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
442  binx,binx,biny,biny);
443 
444  std::vector<TH1D*> hMultiverseHistsGENIE;
445  std::vector<TH1D*> hMultiverseHistsPPFX;
446 
447  //Let's try to area normalize these
448  for(uint i = 0; i < hSystsMultiverse.size()- 100; i++){
449  TH1D* hHolder=
450  (TH1D*)hSystsMultiverse[i]->ProjectionZ
451  (ana::UniqueName().c_str(),binx,binx,biny,biny);
452  float ymax = hHolder->Integral();
453  float cvmax = hCV_Tester->Integral();
454  hHolder->Scale(cvmax/ymax);
455  hMultiverseHistsGENIE.push_back(hHolder);
456  }
457 
458  //PPFX Universes
459  for(uint i = hSystsMultiverse.size()-100;
460  i < hSystsMultiverse.size(); i++){
461  TH1D* hHolder=
462  (TH1D*)hSystsMultiverse[i]->ProjectionZ
463  (ana::UniqueName().c_str(),binx,binx,biny,biny);
464  float ymax = hHolder->Integral();
465  float cvmax = hCV_Tester->Integral();
466  hHolder->Scale(cvmax/ymax);
467  hMultiverseHistsPPFX.push_back(hHolder);
468  }
469 
470  ana::MultiverseCorrelation* multi_correlation_calc_genie =
471  new ana::MultiverseCorrelation(hMultiverseHistsGENIE);
472 
473  TH2F* hCov =
474  (TH2F*)multi_correlation_calc_genie->GetCovarianceMatrix();
475  TH2F* hCor =
476  (TH2F*)multi_correlation_calc_genie->GetCorrelationMatrix();
477 
478  ana::MultiverseCorrelation* multi_correlation_calc_ppfx =
479  new ana::MultiverseCorrelation(hMultiverseHistsPPFX);
480 
481  TH2F* hCovPPFX =
482  (TH2F*)multi_correlation_calc_ppfx->GetCovarianceMatrix();
483  TH2F* hCorPPFX =
484  (TH2F*)multi_correlation_calc_ppfx->GetCorrelationMatrix();
485 
486  for(int i = 1; i <= hCov->GetXaxis()->GetNbins(); i++){
487  for(int j = 1; j <= hCov->GetXaxis()->GetNbins(); j++){
488  hCovHolder->SetBinContent(i,j,hCov->GetBinContent(i,j));
489  hCorHolder->SetBinContent(i,j,hCor->GetBinContent(i,j));
490  hCovHolderPPFX->SetBinContent(i,j,hCovPPFX->GetBinContent(i,j));
491  hCorHolderPPFX->SetBinContent(i,j,hCorPPFX->GetBinContent(i,j));
492  }
493  }
494 
495 
496  if(makePlots){
497  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
498  binx,binx,biny,biny);
499  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
500  hCV1D->GetYaxis()->SetTitle("Template Bins");
501  hCV1D->GetYaxis()->CenterTitle();
502  hCV1D->GetXaxis()->CenterTitle();
503  hCV1D->SetTitle(binTitle.c_str());
504 
505  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
506  hCV1D->Draw("hist");
507  for(uint i = 0; i < hMultiverseHistsPPFX.size(); i++){
508  hMultiverseHistsPPFX[i]->SetLineColorAlpha(kBlue-9,0.50);
509  hMultiverseHistsPPFX[i]->Draw("hist same");
510  }
511  hCV1D->Draw("hist same");
512  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
513  "standard_template_pseudo_exp",
514  "ppfx",binx,biny));
515  c1->Close();
516  }
517 
518  if(makePlots){
519  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
520  binx,binx,biny,biny);
521  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
522  hCV1D->GetXaxis()->SetTitle("Template Bins");
523  hCV1D->SetTitle(binTitle.c_str());
524  hCV1D->GetYaxis()->CenterTitle();
525  hCV1D->GetXaxis()->CenterTitle();
526  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
527  hCV1D->Draw("hist");
528  for(uint i = 0; i < hMultiverseHistsGENIE.size()-1;i++){
529 
530  hMultiverseHistsGENIE[i]->SetLineColorAlpha(kBlue-9,0.50);
531  hMultiverseHistsGENIE[i]->Draw("hist same");
532  }
533  hCV1D->Draw("hist same");
534  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
535  "standard_template_pseudo_exp",
536  "genie",binx,biny));
537  c1->Close();
538  }
539 
540 
541  mapCovariances.insert(std::make_pair("genie",
542  hCovHolder));
543  mapCorrelations.insert(std::make_pair("genie",
544  hCorHolder));
545  mapCovariances.insert(std::make_pair("ppfx",
546  hCovHolderPPFX));
547  mapCorrelations.insert(std::make_pair("ppfx",
548  hCorHolderPPFX));
549  delete hCorPPFX;
550  delete hCovPPFX;
551  delete multi_correlation_calc_ppfx;
552  for(auto ele :hMultiverseHistsPPFX)
553  delete ele;
554  delete hCor;
555  delete hCov;
556  for(auto ele :hMultiverseHistsGENIE)
557  delete ele;
558  delete multi_correlation_calc_genie;
559 
560  }
561 
562  //Now Plot All Covariance Matrices Seperately
563  if(makePlots){
564  for(auto pair: mapCovariances){
565  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
566  pair.second->GetXaxis()->CenterTitle();
567  pair.second->GetYaxis()->CenterTitle();
568  pair.second->GetZaxis()->CenterTitle();
569  pair.second->GetZaxis()->SetTitle("Covariance (Events^{2})");
570  pair.second->SetTitle(binTitle.c_str());
571  c3->SetLeftMargin(0.1);
572  c3->SetRightMargin(0.15);
573  c3->SetBottomMargin(0.1);
574  pair.second->GetZaxis()->SetMaxDigits(3);
575  pair.second->Draw("COLZ");
576  c3->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
577  "covariance_single",
578  pair.first.c_str(),binx,biny));
579  c3->Close();
580  delete c3;
581  }
582  }//if makePlots
583 
584  if(makePlots){
585  for(auto pair: mapCorrelations){
586  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
587  pair.second->GetXaxis()->CenterTitle();
588  pair.second->GetYaxis()->CenterTitle();
589  pair.second->GetZaxis()->CenterTitle();
590  pair.second->SetTitle(binTitle.c_str());
591  pair.second->GetZaxis()->SetTitle("Correlation");
592  pair.second->GetZaxis()->SetRangeUser(-1,1);
593  c3->SetLeftMargin(0.1);
594  c3->SetRightMargin(0.15);
595  c3->SetBottomMargin(0.1);
596  pair.second->GetZaxis()->SetMaxDigits(3);
597  pair.second->Draw("COLZ");
598  c3->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
599  "correlation_single",
600  pair.first.c_str(),binx,biny));
601  c3->Close();
602  delete c3;
603  }
604  }//if makePlots
605 
606  //Now start adding up each of the individual covariance matrices
607  TH2F* hCovariance =
608  (TH2F*)mapCovariances["genie"]->Clone(ana::UniqueName().c_str());
609 
610 
611 
612 
613  std::string output = covarianceDir + "covariance_addition";
614  for(auto pair: mapCovariances){
615  if(pair.first == "genie")continue;
616  output+="_";
617  output+=pair.first;
618  hCovariance->Add(pair.second);
619  if(makePlots){
620  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
621  c3->SetLeftMargin(0.1);
622  c3->SetRightMargin(0.15);
623  c3->SetBottomMargin(0.1);
624  hCovariance->SetTitle(binTitle.c_str());
625  hCovariance->GetZaxis()->SetMaxDigits(3);
626  hCovariance->GetZaxis()->SetTitle("Covariance (Events^{2})");
627  hCovariance->Draw("COLZ");
628  c3->SaveAs(Form("%s_%i_%i.png",output.c_str(),binx,biny));
629  c3->Close();
630  delete c3;
631  }
632  }
633 
634  TH2F* hCorrelation =
635  (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
636 
637  for(int i=0; i<numPIDBins; i++){
638  for(int j=0; j<numPIDBins; j++){
639  double covariance = hCovariance->GetBinContent(i+1,j+1);
640  double variance_i = sqrt(hCovariance->GetBinContent(i+1,i+1));
641  double variance_j = sqrt(hCovariance->GetBinContent(j+1,j+1));
642  double correlation = covariance/(variance_i*variance_j);
643  if(isnan(correlation)) correlation = 0;
644  hCorrelation->SetBinContent(i+1,j+1,
645  correlation);
646  }
647  }
648 
649  if(makePlots){
650  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
651  hCorrelation->GetXaxis()->CenterTitle();
652  hCorrelation->GetYaxis()->CenterTitle();
653  hCorrelation->GetZaxis()->CenterTitle();
654  hCorrelation->SetTitle(binTitle.c_str());
655  hCorrelation->GetZaxis()->SetTitle("Correlation");
656  hCorrelation->GetZaxis()->SetRangeUser(-1,1);
657  c3->SetLeftMargin(0.1);
658  c3->SetRightMargin(0.1);
659  c3->SetBottomMargin(0.1);
660  hCorrelation->GetYaxis()->SetMaxDigits(3);
661  hCorrelation->Draw("COLZ");
662  c3->SaveAs(Form("%s%s_%i_%i.png",covarianceDir.c_str(),
663  "correlation_all",binx,biny));
664  c3->Close();
665  delete c3;
666  }
667 
668 
669  //Move the covariance matrix to the holder so that it can be used
670  //in the fit
671  for(int i = 0; i < NumberTemplateBins; i++){
672  for(int j = 0; j < NumberTemplateBins; j++){
673  CovarianceMatrix[i][j] = hCovariance->GetBinContent(i+1,j+1);
674  if((i == j) && hCovariance->GetBinContent(i+1,j+1) < 1)
675  CovarianceMatrix[i][j] = 10;
676  }
677  }
678 
679  for(uint i = 0; i < hSystsShifts.size(); i++){
680  delete hSystsShifts[i];
681  }
682  hSystsShifts.clear();
683  }
ofstream output
enum BeamMode kRed
Calculate bin to bin correlation matrix from the universe technique.
TH2D * GetCovarianceMatrix()
Get covariance matrix.
void makePlots()
Definition: makePlots.C:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
std::string GetAnalysisBinCaption(int binx, int biny)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< TH1F * > CalculateOneSigmaShift(int binx, int biny)
double CovarianceMatrix[numPIDBins][numPIDBins]
const std::string covarianceDir
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Double_t ymax
Definition: plot.C:25
TH2D * GetCorrelationMatrix()
Get correlation matrix.
const double j
Definition: BetheBloch.cxx:29
c1
Definition: demo5.py:24
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
enum BeamMode kBlue
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string
unsigned int uint
void ana::nueccinc::NueCCIncTemplateFitter::FillStatisticCovarianceMatrix ( int  binx,
int  biny 
)
private

Definition at line 685 of file NueCCIncTemplateFitter.cxx.

References hCV3D, hData3D, MECModelEnuComparisons::i, calib::j, NumberTemplateBins, cet::pow(), and StatCovarianceMatrix.

Referenced by doFullFit().

687  {
688  for(int i = 1; i <= NumberTemplateBins; i++){
689  for(int j = 1; j <= NumberTemplateBins; j++){
690  if(i == j){
691  StatCovarianceMatrix[i-1][j-1] =
692  pow(hData3D->GetBinError(binx,biny,i),2) +
693  pow(hCV3D->GetBinError(binx,biny,i),2);
694  }
695  else
696  StatCovarianceMatrix[i-1][j-1] = 0;
697  }
698  }
699  }
constexpr T pow(T x)
Definition: pow.h:75
const double j
Definition: BetheBloch.cxx:29
double StatCovarianceMatrix[numPIDBins][numPIDBins]
void ana::nueccinc::NueCCIncTemplateFitter::FillValueHolders ( int  binx,
int  biny 
)
private

Definition at line 69 of file NueCCIncTemplateFitter.cxx.

References DataValues, hData3D, hTemplates3D, MECModelEnuComparisons::i, NC_Values, NumberTemplateBins, NumuLike_Values, Other_Values, and SignalLike_Values.

Referenced by doFullFit(), and PrintValueHolders().

70  {
71 
72  for(int i = 0; i < NumberTemplateBins; i++){
73  DataValues[i] = ((double)hData3D->GetBinContent(binx,biny,i+1));
75  (double)hTemplates3D[1]->GetBinContent(binx,
76  biny,i+1)
77  + (double)hTemplates3D[2]->GetBinContent(binx,
78  biny,
79  i+1)
80  + (double)hTemplates3D[7]->GetBinContent(binx,
81  biny,
82  i+1);
83  NumuLike_Values[i] =
84  (double)hTemplates3D[3]->GetBinContent(binx,
85  biny,i+1)
86  + (double)hTemplates3D[4]->GetBinContent(binx,
87  biny,
88  i+1);
89  NC_Values[i] =
90  (double)hTemplates3D[5]->GetBinContent(binx,biny,i+1);
91  Other_Values[i] =
92  (double)hTemplates3D[6]->GetBinContent(binx,biny,i+1);
93  }//loop over template bins
94  }
std::string ana::nueccinc::NueCCIncTemplateFitter::GetAnalysisBinCaption ( int  binx,
int  biny 
)
private

Definition at line 119 of file NueCCIncTemplateFitter.cxx.

References make_true_q0q3_plots::caption, hCV3D, and string.

Referenced by doFullFit(), FillCovarianceMatrix(), getCorrelationMatrixTemplateBins(), getCovarianceMatrixTemplateBins(), getPIDFitTemplates(), and getPIDTemplates().

121  {
122  std::stringstream low_X;
123  low_X << std::fixed << std::setprecision(2) <<
124  hCV3D->GetXaxis()->GetBinLowEdge(binx);
125  std::stringstream hi_X;
126  hi_X << std::fixed << std::setprecision(2) <<
127  hCV3D->GetXaxis()->GetBinUpEdge(binx);
128 
129  std::stringstream low_Y;
130  low_Y << std::fixed << std::setprecision(2) <<
131  hCV3D->GetYaxis()->GetBinLowEdge(biny);
132  std::stringstream hi_Y;
133  hi_Y << std::fixed << std::setprecision(2) <<
134  hCV3D->GetYaxis()->GetBinUpEdge(biny);
135 
136  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
137  hi_X.str() + " " + low_Y.str() + " #leq E_{e} < " + hi_Y.str();
138  return caption;
139  }
enum BeamMode string
std::vector< TH3F * > ana::nueccinc::NueCCIncTemplateFitter::GetAnalysisTemplates ( )
private

Definition at line 109 of file NueCCIncTemplateFitter.cxx.

References hAnalysis3D, MECModelEnuComparisons::i, and ana::UniqueName().

Referenced by doFullFit(), and getFitNormalizationAndErrors().

110  {
111  std::vector<TH3F*> hResults;
112  for(unsigned int i = 0; i < hAnalysis3D.size(); i++)
113  hResults.push_back((TH3F*)
114  hAnalysis3D[i]->Clone
115  (ana::UniqueName().c_str()));
116  return hResults;
117  }
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH2F * ana::nueccinc::NueCCIncTemplateFitter::getCorrelationBetweenSysts ( int  binx,
int  biny 
)

Definition at line 1929 of file NueCCIncTemplateFitter.cxx.

References FillCovarianceBetweenSysts(), MECModelEnuComparisons::i, calib::j, std::sqrt(), ana::UniqueName(), and febshutoff_auto::val.

1930  {
1931  TH2F* hCovariance =
1933 
1934  TH2F* hResult = (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
1935 
1936  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1937  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1938  double val = hCovariance->GetBinContent(i,j);
1939  double dia_i = sqrt(hCovariance->GetBinContent(i,i));
1940  double dia_j = sqrt(hCovariance->GetBinContent(j,j));
1941 
1942  hResult->SetBinContent(i,j,val/(dia_i*dia_j));
1943  }//loop over y bins
1944  }//loop over x bins
1945 
1946  hResult->GetZaxis()->SetTitle("Correlation");
1947 
1948  return hResult;
1949  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double j
Definition: BetheBloch.cxx:29
TH2F * FillCovarianceBetweenSysts(int binx, int biny)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TH2F * ana::nueccinc::NueCCIncTemplateFitter::getCorrelationMatrixTemplateBins ( int  binx,
int  biny 
)

Definition at line 1870 of file NueCCIncTemplateFitter.cxx.

References make_true_q0q3_plots::caption, FillCovarianceMatrix(), GetAnalysisBinCaption(), GetCovarianceMatrixValue(), GetTemplates(), GetXaxis(), MECModelEnuComparisons::i, calib::j, std::sqrt(), string, ana::UniqueName(), and febshutoff_auto::val.

1872  {
1874  std::vector<TH1F*> hPIDAxis =
1875  NueCCIncTemplateFitter::GetTemplates(binx,biny,"nominal");
1876 
1877  TH2F* hCovariance =
1878  new TH2F(ana::UniqueName().c_str(),
1879  ";Template Bins;Template Bins;Covariance (Events^{2})",
1880  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1881  hPIDAxis[0]->GetXaxis()->GetNbins(),
1882  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1883  hPIDAxis[0]->GetXaxis()->GetNbins());
1884 
1885  TH2F* hCorrelation =
1886  new TH2F(ana::UniqueName().c_str(),
1887  ";Template Bins;Template Bins;Correlation",
1888  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1889  hPIDAxis[0]->GetXaxis()->GetNbins(),
1890  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1891  hPIDAxis[0]->GetXaxis()->GetNbins());
1892 
1893  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1894  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1895  double val =
1897  hCovariance->SetBinContent(i,j,val);
1898  }//loop over y bins
1899  }//loop over x bins
1900 
1901  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1902  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1903  double val = hCovariance->GetBinContent(i,j);
1904  double dia_i = sqrt(hCovariance->GetBinContent(i,i));
1905  double dia_j = sqrt(hCovariance->GetBinContent(j,j));
1906 
1907  hCorrelation->SetBinContent(i,j,val/(dia_i*dia_j));
1908  }//loop over y bins
1909  }//loop over x bins
1910 
1911  std::string caption =
1913  hCorrelation->SetTitle(caption.c_str());
1914 
1915  return ((TH2F*)hCorrelation->Clone("hCorrelation"));
1916  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string GetAnalysisBinCaption(int binx, int biny)
correl_xv GetXaxis() -> SetDecimals()
const double j
Definition: BetheBloch.cxx:29
void FillCovarianceMatrix(int binx, int biny, bool makePlots)
std::vector< TH1F * > GetTemplates(int binx, int biny, std::string name)
double GetCovarianceMatrixValue(int binx, int biny)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string
TH2F * ana::nueccinc::NueCCIncTemplateFitter::getCovarianceBetweenSysts ( int  binx,
int  biny 
)

Definition at line 1919 of file NueCCIncTemplateFitter.cxx.

References FillCovarianceBetweenSysts().

1920  {
1921  TH2F* hResult =
1923 
1924 
1925  return hResult;
1926 
1927  }
TH2F * FillCovarianceBetweenSysts(int binx, int biny)
TH2F * ana::nueccinc::NueCCIncTemplateFitter::getCovarianceMatrixTemplateBins ( int  binx,
int  biny 
)

Definition at line 1838 of file NueCCIncTemplateFitter.cxx.

References make_true_q0q3_plots::caption, FillCovarianceMatrix(), GetAnalysisBinCaption(), GetCovarianceMatrixValue(), GetTemplates(), GetXaxis(), MECModelEnuComparisons::i, calib::j, string, ana::UniqueName(), and febshutoff_auto::val.

1840  {
1842 
1843  std::vector<TH1F*> hPIDAxis =
1844  NueCCIncTemplateFitter::GetTemplates(binx,biny,"nominal");
1845 
1846  TH2F* hCovariance =
1847  new TH2F(ana::UniqueName().c_str(),//"hCovariance_holder",
1848  ";Template Bins;Template Bins;Covariance (Events^{2})",
1849  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1850  hPIDAxis[0]->GetXaxis()->GetNbins(),
1851  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1852  hPIDAxis[0]->GetXaxis()->GetNbins());
1853 
1854  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1855  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1856  double val =
1858  hCovariance->SetBinContent(i,j,val);
1859  }//loop over y bins
1860  }//loop over x bins
1861 
1862  std::string caption =
1864  hCovariance->SetTitle(caption.c_str());
1865 
1866  return ((TH2F*)hCovariance->Clone("hCovariance"));
1867  //ana::UniqueName().c_str()));
1868  }
std::string GetAnalysisBinCaption(int binx, int biny)
correl_xv GetXaxis() -> SetDecimals()
const double j
Definition: BetheBloch.cxx:29
void FillCovarianceMatrix(int binx, int biny, bool makePlots)
std::vector< TH1F * > GetTemplates(int binx, int biny, std::string name)
double GetCovarianceMatrixValue(int binx, int biny)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string
double ana::nueccinc::NueCCIncTemplateFitter::GetCovarianceMatrixValue ( int  binx,
int  biny 
)
private

Definition at line 701 of file NueCCIncTemplateFitter.cxx.

References CovarianceMatrix.

Referenced by getCorrelationMatrixTemplateBins(), getCovarianceMatrixTemplateBins(), myFunction2Vars(), and myFunction3Vars().

702  {
703  return CovarianceMatrix[binx][biny];
704  }
double CovarianceMatrix[numPIDBins][numPIDBins]
std::vector< std::pair< TH2F *, TH2F * > > ana::nueccinc::NueCCIncTemplateFitter::getFitNormalizationAndErrors ( )

Definition at line 1787 of file NueCCIncTemplateFitter.cxx.

References GetAnalysisTemplates(), GetTemplateWeightsAndErrors(), make_pair(), moon_position_table_new3::second, and ana::UniqueName().

1788  {
1789  std::vector<std::pair<TH2F*,TH2F*>> hResults;
1790 
1791  std::vector<TH3F*> vAnalysisBinning =
1793 
1794  TH2F* hTemplate = (TH2F*)vAnalysisBinning.front()->Project3D("yx");
1795  hTemplate->SetName("hFitWeightsErrorTemplate");
1796 
1797  TH2F* hSignalNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1798  TH2F* hSignalError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1799 
1800  TH2F* hNumuNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1801  TH2F* hNumuError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1802 
1803  TH2F* hNCNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1804  TH2F* hNCError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1805 
1806  for(int xbin = 1; xbin <= hTemplate->GetXaxis()->GetNbins(); xbin++){
1807  for(int ybin = 1; ybin <= hTemplate->GetYaxis()->GetNbins(); ybin++){
1808  std::vector<std::pair<double,double>> fit_results =
1810  hSignalNorm->SetBinContent(xbin,ybin,fit_results[0].first);
1811  hSignalError->SetBinContent(xbin,ybin,fit_results[0].second);
1812  hNumuNorm->SetBinContent(xbin,ybin,fit_results[1].first);
1813  hNumuError->SetBinContent(xbin,ybin,fit_results[1].second);
1814  hNCNorm->SetBinContent(xbin,ybin,fit_results[2].first);
1815  hNCError->SetBinContent(xbin,ybin,fit_results[2].second);
1816  }//loop over ybins
1817  }//loop over xbins
1818 
1819  hResults.push_back(std::make_pair(hSignalNorm,hSignalError));
1820  hResults.push_back(std::make_pair(hNumuNorm,hNumuError));
1821  hResults.push_back(std::make_pair(hNCNorm,hNCError));
1822 
1823  return hResults;
1824  }
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< std::pair< double, double > > GetTemplateWeightsAndErrors(int binx, int biny)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::vector< TH1F * > ana::nueccinc::NueCCIncTemplateFitter::getPIDFitTemplates ( int  binx,
int  biny 
)

Definition at line 1755 of file NueCCIncTemplateFitter.cxx.

References bin, make_true_q0q3_plots::caption, GetAnalysisBinCaption(), GetTemplates(), GetTemplateWeightsAndErrors(), nc, ana::nue, ana::nuebar, ana::numubar, fhicl::other, and string.

1757  {
1758  std::vector<TH1F*> hResults =
1759  NueCCIncTemplateFitter::GetTemplates(binx, biny,"fit");
1760  std::string caption =
1762  std::vector<std::pair<double,double>> vWeights =
1764  for(int bin = 1; bin <= hResults[0]->GetXaxis()->GetNbins(); bin++){
1765  double nue = hResults[1]->GetBinContent(bin) *vWeights[0].first;
1766  double nuebar = hResults[2]->GetBinContent(bin) *vWeights[0].first;
1767  double numu = hResults[3]->GetBinContent(bin) *vWeights[1].first;
1768  double numubar= hResults[4]->GetBinContent(bin) *vWeights[1].first;
1769  double nc = hResults[5]->GetBinContent(bin) *vWeights[2].first;
1770  double other = hResults[6]->GetBinContent(bin);
1771  double nuenf = hResults[7]->GetBinContent(bin) *vWeights[0].first;
1772 
1773  hResults[0]->SetBinContent(bin,nue+nuebar+numu+numubar+nc+other+nuenf);
1774  hResults[1]->SetBinContent(bin,nue);
1775  hResults[2]->SetBinContent(bin,nuebar);
1776  hResults[3]->SetBinContent(bin,numu);
1777  hResults[4]->SetBinContent(bin,numubar);
1778  hResults[5]->SetBinContent(bin,nc);
1779  hResults[6]->SetBinContent(bin,other);
1780  hResults[7]->SetBinContent(bin,nuenf);
1781  }//loop over template bins
1782  hResults[0]->SetTitle(caption.c_str());
1783  return hResults;
1784  }
std::string GetAnalysisBinCaption(int binx, int biny)
std::vector< std::pair< double, double > > GetTemplateWeightsAndErrors(int binx, int biny)
float bin[41]
Definition: plottest35.C:14
std::vector< TH1F * > GetTemplates(int binx, int biny, std::string name)
enum BeamMode nc
enum BeamMode string
std::vector< TH1F * > ana::nueccinc::NueCCIncTemplateFitter::getPIDTemplates ( int  binx,
int  biny 
)

Definition at line 1827 of file NueCCIncTemplateFitter.cxx.

References make_true_q0q3_plots::caption, GetAnalysisBinCaption(), GetTemplates(), and string.

1829  {
1830  std::vector<TH1F*> hResults =
1831  NueCCIncTemplateFitter::GetTemplates(binx, biny,"nominal");
1832  std::string caption =
1834  hResults[0]->SetTitle(caption.c_str());
1835  return hResults;
1836  }
std::string GetAnalysisBinCaption(int binx, int biny)
std::vector< TH1F * > GetTemplates(int binx, int biny, std::string name)
enum BeamMode string
double ana::nueccinc::NueCCIncTemplateFitter::GetStatisticCovarianceMatrixValue ( int  binx,
int  biny 
)
private

Definition at line 705 of file NueCCIncTemplateFitter.cxx.

References StatCovarianceMatrix.

Referenced by myFunction2Vars(), and myFunction3Vars().

707  {
708  return StatCovarianceMatrix[binx][biny];
709  }
double StatCovarianceMatrix[numPIDBins][numPIDBins]
std::vector< TH1F * > ana::nueccinc::NueCCIncTemplateFitter::GetTemplates ( int  binx,
int  biny,
std::string  name 
)
private

Definition at line 993 of file NueCCIncTemplateFitter.cxx.

References hTemplates3D, MECModelEnuComparisons::i, and ana::UniqueName().

Referenced by getCorrelationMatrixTemplateBins(), getCovarianceMatrixTemplateBins(), getPIDFitTemplates(), and getPIDTemplates().

995  {
996  std::vector<TH1F*> hResults;
997  for(uint i = 0; i < hTemplates3D.size(); i++){
998  hResults.push_back((TH1F*)hTemplates3D[i]->ProjectionZ
999  (ana::UniqueName().c_str(),binx,binx,biny,biny));
1000  }
1001  return hResults;
1002  }
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
unsigned int uint
std::vector< std::pair< double, double > > ana::nueccinc::NueCCIncTemplateFitter::GetTemplateWeightsAndErrors ( int  binx,
int  biny 
)
private

Definition at line 975 of file NueCCIncTemplateFitter.cxx.

References hNCError2D, hNCWeights2D, hNumuCCError2D, hNumuCCWeights2D, hSignalError2D, hSignalWeights2D, and make_pair().

Referenced by getFitNormalizationAndErrors(), and getPIDFitTemplates().

976  {
977  std::vector<std::pair<double,double>> vResult;
978 
979  vResult.push_back(std::make_pair(hSignalWeights2D->
980  GetBinContent(binx,biny),
982  GetBinContent(binx,biny)));
983  vResult.push_back(std::make_pair(hNumuCCWeights2D->
984  GetBinContent(binx,biny),
986  GetBinContent(binx,biny)));
987  vResult.push_back(std::make_pair(hNCWeights2D->GetBinContent(binx,biny),
988  hNCError2D->GetBinContent(binx,biny)));
989 
990  return vResult;
991  }
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
double ana::nueccinc::NueCCIncTemplateFitter::myFunction2Vars ( double  bkgd_par,
double  nue_par 
)
private

Definition at line 781 of file NueCCIncTemplateFitter.cxx.

References DataValues, GetCovarianceMatrixValue(), GetStatisticCovarianceMatrixValue(), MECModelEnuComparisons::i, calib::j, nc, NC_Values, NumberTemplateBins, ana::nueccinc::numPIDBins, NumuLike_Values, fhicl::other, Other_Values, and SignalLike_Values.

Referenced by fcn2Var().

783  {
784  float data = 0;
785  float signal = 0;
786  float numu = 0;
787  float nc = 0;
788  float other = 0;
789  // Covariance Matrix Inversion First
790  TMatrixF Vij_stat = TMatrix(numPIDBins,numPIDBins);
791  TMatrixF Vij = TMatrix(numPIDBins,numPIDBins);
792  for(int i =0; i < numPIDBins; i++){
793  for(int j = 0; j < NumberTemplateBins; j++){
794  Vij_stat[i][j] =
796  Vij[i][j] =
798  }
799  }
800  //Inversion
801  TMatrixD H3 = Vij_stat + Vij;
802  TDecompSVD svd(H3);
803  TMatrixD VijInv = svd.Invert();
804 
805  TMatrixD FitVector = TMatrixD(1,NumberTemplateBins);
806  TMatrixD FitVector2 = TMatrixD(NumberTemplateBins,1);
807 
808  //Calculate Chi Sq
809  for(int i =0; i < NumberTemplateBins; i++){
810  data = DataValues[i];
811  signal = SignalLike_Values[i];
812  numu = NumuLike_Values[i];
813  nc = NC_Values[i];
814  other = Other_Values[i];
815 
816  if( (data < 30) || (signal + numu + nc + other < 30)){
817  FitVector[0][i] = 0;
818  FitVector2[i][0] = 0;
819  continue;
820  }
821 
822  FitVector[0][i] = ( data - (nue_par*signal + bkgd_par*(numu + nc)
823  + other));
824 
825  FitVector2[i][0] = ( data -(nue_par*signal + bkgd_par*(numu + nc)
826  + other));
827  }
828 
829  TMatrixD Chi = FitVector * VijInv;
830  TMatrixD Chi2 = Chi * FitVector2;
831  float chisquare = Chi2(0,0);
832  return chisquare;
833  }
double GetStatisticCovarianceMatrixValue(int binx, int biny)
const XML_Char const XML_Char * data
Definition: expat.h:268
const double j
Definition: BetheBloch.cxx:29
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
enum BeamMode nc
double GetCovarianceMatrixValue(int binx, int biny)
double ana::nueccinc::NueCCIncTemplateFitter::myFunction3Vars ( double  numu_par,
double  nue_par,
double  nc_par 
)
private

Definition at line 719 of file NueCCIncTemplateFitter.cxx.

References DataValues, GetCovarianceMatrixValue(), GetStatisticCovarianceMatrixValue(), MECModelEnuComparisons::i, calib::j, nc, NC_Values, NumberTemplateBins, NumuLike_Values, fhicl::other, Other_Values, and SignalLike_Values.

Referenced by fcn3Var().

722  {
723  float data = 0;
724  float signal = 0;
725  float numu = 0;
726  float nc = 0;
727  float other = 0;
728  // Covariance Matrix Inversion First
729  TMatrixF Vij_stat = TMatrix(NumberTemplateBins,NumberTemplateBins);
730  TMatrixF Vij = TMatrix(NumberTemplateBins,NumberTemplateBins);
731  for(int i =0; i < NumberTemplateBins; i++){
732  for(int j = 0; j < NumberTemplateBins; j++){
733  Vij_stat[i][j] =
735  Vij[i][j] =
737  }
738  }
739  //Inversion
740  TMatrixD H3 = Vij_stat + Vij;
741  TDecompSVD svd(H3);
742  TMatrixD VijInv = svd.Invert();
743 
744  TMatrixD FitVector = TMatrixD(1,NumberTemplateBins);
745  TMatrixD FitVector2 = TMatrixD(NumberTemplateBins,1);
746 
747  //Calculate Chi Sq
748  for(int i =0; i < NumberTemplateBins; i++){
749  data = Fitter_obj->DataValues[i];
750  signal = Fitter_obj->SignalLike_Values[i];
751  numu = Fitter_obj->NumuLike_Values[i];
752  nc = Fitter_obj->NC_Values[i];
753  other = Fitter_obj->Other_Values[i];
754 
755  if( (data < 30) || (signal + numu + nc + other < 30)){
756  FitVector[0][i] = 0;
757  FitVector2[i][0] = 0;
758  continue;
759  }
760 
761  FitVector[0][i] = ( data - (nue_par*signal + numu_par*numu +
762  nc_par*nc + other));
763 
764  FitVector2[i][0] = ( data -(nue_par*signal + numu_par*numu +
765  nc_par*nc + other));
766  }
767 
768  TMatrixD Chi = FitVector * VijInv;
769  TMatrixD Chi2 = Chi * FitVector2;
770  float chisquare = Chi2(0,0);
771  return chisquare;
772  }
double GetStatisticCovarianceMatrixValue(int binx, int biny)
const XML_Char const XML_Char * data
Definition: expat.h:268
const double j
Definition: BetheBloch.cxx:29
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
enum BeamMode nc
double GetCovarianceMatrixValue(int binx, int biny)
NueCCIncTemplateFitter * Fitter_obj
void ana::nueccinc::NueCCIncTemplateFitter::PrintValueHolders ( int  binx,
int  biny 
)

Definition at line 1108 of file NueCCIncTemplateFitter.cxx.

References om::cout, DataValues, allTimeWatchdog::endl, FillValueHolders(), MECModelEnuComparisons::i, NC_Values, ana::nueccinc::numPIDBins, NumuLike_Values, Other_Values, and SignalLike_Values.

1109  {
1110  std::cout << "Values For Bins: " << binx << "\t" << biny << std::endl;
1112 
1113  for(unsigned int i = 0; i < numPIDBins; i++){
1114  std::cout << "Data: " << DataValues[i] << "\tSignal: " <<
1115  SignalLike_Values[i] << "\tNumu: " << NumuLike_Values[i] <<
1116  "\tNC: " << NC_Values[i] << "\tOther: " <<
1117  Other_Values[i] << std::endl;
1118  }
1119  }
OStream cout
Definition: OStream.cxx:6
void ana::nueccinc::NueCCIncTemplateFitter::PropagateFitUncertainty3D ( int  binx,
int  biny 
)
private

Definition at line 856 of file NueCCIncTemplateFitter.cxx.

References om::cout, allTimeWatchdog::endl, fit_covariance_matrix, hAnalysis3D, makeTrainCVSamples::int, nc_wei, nue_wei, numu_wei, cet::pow(), and std::sqrt().

Referenced by doFullFit().

857  {
858  for(int numChn = 0; numChn < (int)hAnalysis3D.size(); numChn++){
859  for(int zbin = 1;
860  zbin <= hAnalysis3D[0]->GetZaxis()->GetNbins();zbin++){
861  float nsig = hAnalysis3D[1]->GetBinContent(binx,biny,zbin);
862  float nnuebar = hAnalysis3D[2]->GetBinContent(binx,biny,zbin);
863  float nnumu = hAnalysis3D[3]->GetBinContent(binx,biny,zbin);
864  float nnumubar = hAnalysis3D[4]->GetBinContent(binx,biny,zbin);
865  float nnc = hAnalysis3D[5]->GetBinContent(binx,biny,zbin);
866  float nother = hAnalysis3D[6]->GetBinContent(binx,biny,zbin);
867  float nnf = hAnalysis3D[7]->GetBinContent(binx,biny,zbin);
868  if(numChn == 0){
869  //Add all Channels to total
870  float binCont = nsig * nue_wei.first;
871  binCont += nnuebar * nue_wei.first;
872  binCont += nnumu * numu_wei.first;
873  binCont += nnumubar * numu_wei.first;
874  binCont += nnc * nc_wei.first;
875  binCont += nother;
876  binCont += nnf * nue_wei.first;
877 
878  std::cout << hAnalysis3D[0]->GetBinContent(binx,biny,zbin)
879  << "\t" << binCont << std::endl;
880 
881  hAnalysis3D[0]->SetBinContent(binx,biny,zbin,binCont);
882 
883 
884 
885  //MINOS Error calculations already take the correlations between
886  //parameters into account
887  float binErr2 =
888  pow(fit_covariance_matrix[0][0],2)*pow((nnumu+nnumubar),2) +
889  pow(fit_covariance_matrix[1][1],2)* pow((nsig+nnuebar+nnf),2) +
890  pow(fit_covariance_matrix[2][2],2)* pow((nnc),2) +
891  nsig * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
892  nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
893  nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
894  nother;
895 
896  if(printVerbose)
897  std::cout << binx << "\t" << biny << "\t" << zbin << "\t"
898  << sqrt(binErr2) <<std::endl;
899 
900  if(fit_covariance_matrix[0][1] == 0 &&
901  fit_covariance_matrix[0][2] == 0 &&
902  fit_covariance_matrix[1][2] == 0 && nue_wei.second == 0 &&
903  numu_wei.second == 0 && nc_wei.second == 0){
904  binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
905  nnc + nother;
906  }
907  hAnalysis3D[0]->SetBinError(binx,biny,zbin,sqrt(binErr2));
908  }
909  if(numChn == 1){
910  //Signal Estimate
911  float binCont = nsig * nue_wei.first;
912 
913  //Just uncertainty on Signal estimate
914  //Signal + Background (taken into account in TMinuit internally
915  //when calculating the internal error matrix)
916 
917  //Need to add wrong sign estimate uncertainty
918  //Lets just add this in quadrature to the
919  //signal estimate + background estimate uncertainties
920  //We add uncertainties on the signal template estimate and
921  //statistical uncertainties on the WS component
922  float binErr = sqrt( pow(nue_wei.second,2)*pow(nsig,2) +
923  nsig * pow(nue_wei.first,2) +
924  pow(nue_wei.second,2)*pow(nnuebar,2) +
925  nnuebar * pow(nue_wei.first,2));
926 
927  if(printVerbose)
928  std::cout << binx << "\t" << biny << "\t" << zbin << "\t"
929  << nsig << "\t" << "\t" << sqrt(nsig) << "\t"
930  << binErr << std::endl;
931 
932  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
933  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
934  }
935  if(numChn == 2){
936  float binCont = nnuebar * nue_wei.first;
937  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
938  nnuebar * pow(nue_wei.first,2));
939  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
940  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
941  }
942  if(numChn == 3){
943  float binCont = nnumu * numu_wei.first;
944  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
945  nnumu * pow(numu_wei.first,2));
946  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
947  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
948  }
949  if(numChn == 4){
950  float binCont = nnumubar * numu_wei.first;
951  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
952  nnumubar * pow(numu_wei.first,2));
953  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
954  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
955  }
956  if(numChn == 5){
957  float binCont = nnc * nc_wei.first;
958  float binErr = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
959  nnc * pow(nc_wei.first,2));
960  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
961  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
962  }
963  if(numChn == 7){
964  float binCont = nnf * nue_wei.first;
965  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
966  nnf * pow(nue_wei.first,2));
967  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
968  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
969  }
970  }//loop over all Z axis bins (not electron kinematic bin)
971  }//loop over all interaction channels
972  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
OStream cout
Definition: OStream.cxx:6
void ana::nueccinc::NueCCIncTemplateFitter::SetBinWeightsAndErrors2D ( int  binx,
int  biny 
)
private

Definition at line 835 of file NueCCIncTemplateFitter.cxx.

References hChiSquaredResults2D, hNCError2D, hNCWeights2D, hNumuCCError2D, hNumuCCWeights2D, hSignalError2D, hSignalWeights2D, nc_wei, nue_wei, and numu_wei.

Referenced by doFullFit().

836  {
837  //Weights are good here
838  double signal_weis = nue_wei.first;
839  double nummu_weis = numu_wei.first;
840  double nc_weis = nc_wei.first;
841 
842  double signal_err = nue_wei.second;
843  double nummu_err = numu_wei.second;
844  double nc_err = nc_wei.second;
845 
846  hSignalWeights2D->SetBinContent(binx,biny,signal_weis);
847  hNumuCCWeights2D->SetBinContent(binx,biny,nummu_weis);
848  hNCWeights2D->SetBinContent(binx,biny,nc_weis);
849 
850  hSignalError2D->SetBinContent(binx,biny,signal_err);
851  hNumuCCError2D->SetBinContent(binx,biny,nummu_err);
852  hNCError2D->SetBinContent(binx,biny,nc_err);
853  hChiSquaredResults2D->SetBinContent(binx,biny,chisquared);
854  }

Member Data Documentation

double ana::nueccinc::NueCCIncTemplateFitter::chisquared

Definition at line 81 of file NueCCIncTemplateFitter.h.

double ana::nueccinc::NueCCIncTemplateFitter::CovarianceMatrix[numPIDBins][numPIDBins]
private

Definition at line 145 of file NueCCIncTemplateFitter.h.

Referenced by FillCovarianceMatrix(), and GetCovarianceMatrixValue().

double ana::nueccinc::NueCCIncTemplateFitter::DataValues[numPIDBins]
private
double ana::nueccinc::NueCCIncTemplateFitter::fit_covariance_matrix[3][3]

Definition at line 80 of file NueCCIncTemplateFitter.h.

Referenced by doFullFit(), and PropagateFitUncertainty3D().

std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hAnalysis3D
private

Definition at line 125 of file NueCCIncTemplateFitter.h.

Referenced by GetAnalysisTemplates(), and PropagateFitUncertainty3D().

TH2F* ana::nueccinc::NueCCIncTemplateFitter::hChiSquaredResults2D
private

Definition at line 138 of file NueCCIncTemplateFitter.h.

Referenced by SetBinWeightsAndErrors2D().

TH3F* ana::nueccinc::NueCCIncTemplateFitter::hCV3D
private
TH3F* ana::nueccinc::NueCCIncTemplateFitter::hData3D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hNCError2D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hNCWeights2D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hNumuCCError2D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hNumuCCWeights2D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hSignalError2D
private
TH2F* ana::nueccinc::NueCCIncTemplateFitter::hSignalWeights2D
private
std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hSystsDown3D
private
std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hSystsMultiverse
private

Definition at line 130 of file NueCCIncTemplateFitter.h.

Referenced by FillCovarianceBetweenSysts(), and FillCovarianceMatrix().

std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hSystsOneSided
private

Definition at line 129 of file NueCCIncTemplateFitter.h.

Referenced by CalculateOneSigmaShift(), and FillCovarianceMatrix().

std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hSystsUp3D
private
std::vector<TH3F*> ana::nueccinc::NueCCIncTemplateFitter::hTemplates3D
private

Definition at line 124 of file NueCCIncTemplateFitter.h.

Referenced by FillValueHolders(), and GetTemplates().

double ana::nueccinc::NueCCIncTemplateFitter::NC_Values[numPIDBins]
private
std::pair<double,double> ana::nueccinc::NueCCIncTemplateFitter::nc_wei
std::pair<double,double> ana::nueccinc::NueCCIncTemplateFitter::nue_wei
const int ana::nueccinc::NueCCIncTemplateFitter::NumberTemplateBins
private
std::pair<double,double> ana::nueccinc::NueCCIncTemplateFitter::numu_wei
double ana::nueccinc::NueCCIncTemplateFitter::NumuLike_Values[numPIDBins]
private
double ana::nueccinc::NueCCIncTemplateFitter::Other_Values[numPIDBins]
private
double ana::nueccinc::NueCCIncTemplateFitter::SignalLike_Values[numPIDBins]
private
double ana::nueccinc::NueCCIncTemplateFitter::StatCovarianceMatrix[numPIDBins][numPIDBins]
private

The documentation for this class was generated from the following files: