Public Member Functions | Private Member Functions | Private Attributes | List of all members
protoana::PDSPThinSliceFitter Class Reference

#include <PDSPThinSliceFitter.h>

Public Member Functions

 PDSPThinSliceFitter (std::string fcl_file, std::string output_file)
 
void FillMCEvents ()
 
void BuildMCSamples ()
 
void SaveMCSamples ()
 
void GetNominalFluxes ()
 
void BuildDataHists ()
 
void InitializeMCSamples ()
 
void CompareDataMC (std::string extra_name, TDirectory *xsec_dir, TDirectory *plot_dir, bool post_fit=false)
 
void ScaleMCToData ()
 
void ScaleDataToNorm ()
 
void RunFitAndSave ()
 
 ~PDSPThinSliceFitter ()
 

Private Member Functions

void NormalFit ()
 
void Pulls ()
 
void Configure (std::string fcl_file)
 
void DefineFitFunction ()
 
void MakeMinimizer ()
 
void ParameterScans ()
 
void DoThrows (const TH1D &pars, const TMatrixD *cov)
 
void Do1DShifts (const TH1D &pars, bool prefit=false)
 
void SetBestFit ()
 
void GetCurrentTruthHists (std::map< int, std::vector< TH1 * >> &throw_hists, std::map< int, std::vector< TH1 * >> &throw_inc_hists, std::map< int, std::vector< TH1 * >> &throw_xsec_hists)
 
void PlotThrows (std::map< int, std::vector< TH1 * >> &throw_hists, std::map< int, std::vector< TH1 * >> &truth_throw_hists, std::map< int, std::vector< TH1 * >> &truth_inc_hists, std::map< int, std::vector< TH1 * >> &truth_xsec_hists)
 
void BuildFakeDataXSecs (bool use_scales=true)
 
void BuildDataFromToy ()
 
double CalcChi2SystTerm ()
 
double CalcRegTerm ()
 
void MakeThrowsTree (TTree &tree, std::vector< double > &branches)
 
std::vector< double > GetBestFitParsVec ()
 
double BetheBloch (double energy, double mass)
 
double densityEffect (double beta, double gamma)
 

Private Attributes

ThinSliceDriverfThinSliceDriver
 
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
 
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
 
ThinSliceDataSet fDataSet
 
std::map< int, boolfIsSignalSample
 
TFile fMCFile
 
TTree * fMCTree
 
TFile fDataFile
 
TTree * fDataTree
 
TFile fOutputFile
 
ROOT::Math::Functor fFitFunction
 
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
 
std::vector< double > fMinimizerInitVals
 
TH1D fPreFitParsNormal
 
std::map< int, double > fNominalFluxes
 
std::map< int, double > fFakeFluxes
 
std::map< int, std::vector< std::vector< double > > > fFluxesBySample
 
std::map< int, std::vector< std::vector< double > > > fFakeFluxesBySample
 
std::map< int, std::vector< int > > fFluxParsToSamples
 
double fDataFlux
 
double fMCDataScale = 1.
 
std::map< int, std::vector< double > > fSignalParameters
 
std::map< int, std::vector< std::string > > fSignalParameterNames
 
size_t fTotalSignalParameters = 0
 
std::map< int, double > fFluxParameters
 
std::map< int, std::stringfFluxParameterNames
 
size_t fTotalFluxParameters = 0
 
std::map< std::string, ThinSliceSystematicfSystParameters
 
std::vector< std::stringfSystParameterNames
 
std::vector< double > fParLimits
 
std::vector< double > fParLimitsUp
 
size_t fTotalSystParameters = 0
 
std::map< std::string, size_t > fCovarianceBins
 
bool fAddSystTerm
 
bool fAddRegTerm
 
double fRegFactor = 0.
 
TMatrixD * fCovMatrix
 
TMatrixD * fCovMatrixDisplay
 
TDecompChol * fInputChol
 
std::map< std::string, double > fToyValues
 
TRandom3 fRNG
 
std::map< int, std::vector< double > > fFakeDataScales
 
std::map< int, std::vector< double > > fBestFitSignalPars
 
std::map< std::string, ThinSliceSystematicfBestFitSystPars
 
std::map< int, double > fBestFitFluxPars
 
std::map< int, TH1 * > fNominalXSecs
 
std::map< int, TH1 * > fNominalIncs
 
std::map< int, TH1 * > fBestFitXSecs
 
std::map< int, TH1 * > fBestFitIncs
 
std::map< int, TH1 * > fFakeDataXSecs
 
std::map< int, TH1 * > fFakeDataIncs
 
std::map< int, TH1D * > fBestFitSelectionHists
 
std::map< int, std::vector< double > > fBestFitTruthVals
 
std::vector< ThinSliceEventfEvents
 
std::vector< ThinSliceEventfFakeDataEvents
 
std::string fMCFileName
 
std::string fDataFileName
 
std::string fTreeName
 
std::vector< fhicl::ParameterSetfSelectionSets
 
std::vector< fhicl::ParameterSetfSampleSets
 
std::map< int, std::stringfFluxTypes
 
int fMaxCalls
 
int fMaxIterations
 
size_t fNFitSteps = 0
 
unsigned int fNScanSteps
 
double fTolerance
 
double fLowerLimit
 
double fUpperLimit
 
std::vector< std::pair< int, int > > fPlotStyle
 
bool fPlotRebinned
 
bool fRandomStart
 
std::string fDriverName
 
std::string fAnalysis
 
fhicl::ParameterSet fAnalysisOptions
 
double fPitch
 
std::string fSliceMethod
 
bool fMultinomial
 
bool fDoFakeData
 
bool fDoThrows
 
bool fDoScans
 
bool fDo1DShifts
 
bool fDoSysts
 
bool fRunHesse
 
bool fFixVariables
 
std::map< std::string, double > fSystsToFix
 
std::string fFakeDataRoutine
 
bool fDoFluctuateStats
 
bool fSplitMC
 
int fSplitVal = 0
 
bool fFillIncidentInFunction = false
 
bool fFitUnderOverflow = false
 
bool fUseFakeSamples = false
 
bool fFitFlux
 
size_t fNThrows
 
size_t fMaxRethrows
 
std::string fFitType = "Normal"
 
size_t fNPulls
 
bool fDoScaleDataToNorm
 
double fDataNorm
 
std::vector< double > fIncidentRecoBins
 
std::vector< double > fTrueIncidentBins
 
std::vector< double > fBeamEnergyBins
 
std::vector< int > fIncidentSamples
 
std::vector< int > fMeasurementSamples
 
bool fDrawXSecUnderflow
 
std::map< int, std::vector< double > > fSignalBins
 

Detailed Description

Definition at line 30 of file PDSPThinSliceFitter.h.

Constructor & Destructor Documentation

protoana::PDSPThinSliceFitter::PDSPThinSliceFitter ( std::string  fcl_file,
std::string  output_file 
)

Definition at line 41 of file PDSPThinSliceFitter.cxx.

43  : fOutputFile(output_file.c_str(), "RECREATE") {
44  Configure(fcl_file);
45  gStyle->SetOptStat(0);
46  gROOT->SetBatch(1);
47 
48  try {
51  }
52  catch (const std::runtime_error & e) {
54  throw e;
55  }
56 
57 }
void Configure(std::string fcl_file)
const double e
static ThinSliceDriverRegistry * Instance()
ThinSliceDriver * GetDriver(const std::string &name, const fhicl::ParameterSet &extra_options)
fhicl::ParameterSet fAnalysisOptions
protoana::PDSPThinSliceFitter::~PDSPThinSliceFitter ( )

Definition at line 59 of file PDSPThinSliceFitter.cxx.

59  {
60  fOutputFile.Close();
61 }

Member Function Documentation

double protoana::PDSPThinSliceFitter::BetheBloch ( double  energy,
double  mass 
)
inlineprivate

Definition at line 190 of file PDSPThinSliceFitter.h.

190  {
191  double K,rho,Z,A, charge, me, I, gamma, /*momentum ,*/wmax, pitch;
192  long double beta;
193  K = 0.307;
194  rho = 1.4;
195  charge = 1;
196  Z = 18;
197  A = 39.948;
198  I = pow(10,-6)*10.5*18; //MeV
199  me = 0.51; //MeV me*c^2
200  pitch = 1;
201 
202  //momentum = sqrt( pow(energy,2) - pow(massicle,2));
203  //beta = momentum/sqrt(pow(massicle,2) + pow(momentum,2));
204  //gamma = 1/sqrt(1 - pow(beta,2));
205 
206  gamma = (energy + mass) / mass;
207  beta = sqrt( 1 - 1/pow(gamma,2));
208 
209  wmax = 2*me*pow(beta,2)*pow(gamma,2)/(1+2*gamma*me/mass + pow(me,2)/pow(mass,2));
210 
211 
212  double dEdX;
213  //multiply by rho to have dEdX MeV/cm in LAr
214 
215  dEdX = pitch*(rho*K*Z*pow(charge,2))/(A*pow(beta,2))*(0.5*log(2*me*pow(gamma,2)*pow(beta,2)*wmax/pow(I,2)) - pow(beta,2) - densityEffect( beta, gamma )/2 );
216 
217  return dEdX;
218  };
double beta(double KE, const simb::MCParticle *part)
constexpr T pow(T x)
Definition: pow.h:72
double dEdX(double KE, const simb::MCParticle *part)
double densityEffect(double beta, double gamma)
double gamma(double KE, const simb::MCParticle *part)
#define A
Definition: memgrp.cpp:38
void protoana::PDSPThinSliceFitter::BuildDataFromToy ( )
private

Definition at line 1361 of file PDSPThinSliceFitter.cxx.

1361  {
1362  if (!fAddSystTerm) {
1364  = "Error: Need to include systematic term and covariance matrix";
1365  throw std::runtime_error(message);
1366  }
1367 
1368  std::vector<double> vals, init_vals;
1369  for (auto it = fSignalParameters.begin();
1370  it != fSignalParameters.end(); ++it) {
1371  for (size_t i = 0; i < it->second.size(); ++i) {
1372  vals.push_back(it->second.at(i));
1373  init_vals.push_back(it->second.at(i));
1374  }
1375  }
1376  for (auto it = fFluxParameters.begin();
1377  it != fFluxParameters.end(); ++it) {
1378  vals.push_back(it->second);
1379  init_vals.push_back(it->second);
1380  }
1381 
1382  std::vector<size_t> bins;
1383  std::vector<std::string> names;
1384  for (auto it = fSystParameters.begin();
1385  it != fSystParameters.end(); ++it) {
1386  bins.push_back(fCovarianceBins[it->first]);
1387  names.push_back(it->first);
1388  init_vals.push_back(it->second.GetValue());
1389  }
1390 
1391  size_t base_par = fTotalSignalParameters + fTotalFluxParameters;
1392 
1393  if (fAnalysisOptions.get<bool>("SetToy", false)) {
1394  std::vector<double> toy_vals
1395  = fAnalysisOptions.get<std::vector<double>>("ToyVals");
1396  for (double & v : toy_vals) {
1397  vals.push_back(v);
1398  }
1399 
1400  for (size_t i = 0; i < vals.size(); ++i) {
1401  std::cout << vals[i] << " ";
1402  if (i >= base_par) {
1403  std::cout << names[i - base_par];
1404  fToyValues[names[i - base_par]] = vals[i];
1405  }
1406  std::cout << std::endl;
1407  }
1408  }
1409  else {
1410 
1411  bool rethrow = true;
1412  while (rethrow) {
1413  bool all_pos = true;
1414  TVectorD rand(fTotalSystParameters);
1415  for (size_t i = 0; i < fTotalSystParameters; ++i) {
1416  rand[i] = fRNG.Gaus();
1417  }
1418  TVectorD rand_times_chol = fInputChol->GetU()*rand;
1419 
1420  for (size_t i = 0; i < fTotalSystParameters; ++i) {
1421  double val = rand_times_chol[bins[i]] + init_vals[base_par + i];
1422  //if (val < fParLimits[base_par + i]) {
1423  if (val < fSystParameters[names[i]].GetGenThrowLimit()) {
1424  all_pos = false;
1425  std::cout << "Rethrowing " << i << " " << val <<
1426  fSystParameters[names[i]].GetGenThrowLimit() <<
1427  std::endl;
1428  }
1429  //if (val > fParLimitsUp[base_par + i]) {
1430  if (val > fSystParameters[names[i]].GetGenThrowLimitUp()) {
1431  all_pos = false;
1432  std::cout << "Rethrowing " << i << " " << val <<
1433  fSystParameters[names[i]].GetGenThrowLimitUp() <<
1434  std::endl;
1435  }
1436  }
1437 
1438  if (all_pos) {
1439  for (size_t i = 0; i < fTotalSystParameters; ++i) {
1440  if (fSystsToFix.find(names[i]) != fSystsToFix.end()) {
1441 
1442  }
1443  double val = (fSystsToFix.find(names[i]) == fSystsToFix.end() ?
1444  rand_times_chol[bins[i]] + init_vals[base_par + i] :
1445  fSystsToFix[names[i]]);
1446 
1447  vals.push_back(val);
1448  fToyValues[names[i]] = val;
1449  }
1450  }
1451  rethrow = !all_pos;
1452  }
1453 
1454  std::cout << "Vals: " << vals.size() << " total " <<
1455  base_par +
1456  fTotalSystParameters << std::endl;
1457  for (size_t i = 0; i < vals.size(); ++i) {
1458  std::cout << vals[i] << " ";
1459  if (i >= base_par)
1460  std::cout << names[i - base_par];
1461  std::cout << std::endl;
1462  }
1463  }
1464  fFillIncidentInFunction = true;
1465  fUseFakeSamples = true;
1466  fFitFunction(&vals[0]);
1468  BuildFakeDataXSecs(false);
1469 
1470  fUseFakeSamples = false;
1471  //Refill the hists for comparisons
1472  fFitFunction(&init_vals[0]);
1473  fFillIncidentInFunction = false;
1474 
1475 
1476 }
std::string string
Definition: nybbler.cc:12
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< int, double > fFluxParameters
std::map< std::string, double > fToyValues
std::map< std::string, size_t > fCovarianceBins
std::map< int, std::vector< double > > fSignalParameters
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
void FillHistsFromSamples(const std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, double &flux)
std::map< std::string, double > fSystsToFix
void BuildFakeDataXSecs(bool use_scales=true)
static std::vector< std::string > const names
Definition: FragmentType.hh:8
fhicl::ParameterSet fAnalysisOptions
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::BuildDataHists ( )

Definition at line 668 of file PDSPThinSliceFitter.cxx.

668  {
669  //Create the data hists
670  fDataSet = ThinSliceDataSet(fIncidentRecoBins, fSelectionSets);
671 
672  TFile inputFile((!fDoFakeData ? fDataFileName.c_str() : fMCFileName.c_str()),
673  "OPEN");
674  TTree * tree = (TTree*)inputFile.Get(fTreeName.c_str());
675  if (!fDoFakeData) {
677  }
678  else if (fFakeDataRoutine == "Toy") {
680  }
681  else if (fFakeDataRoutine == "Asimov") {
682  std::cout << "Doing Asimov" << std::endl;
683  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
684  std::vector<ThinSliceSample> & samples_vec = it->second[0];
685  fFakeDataScales[it->first] = std::vector<double>(samples_vec.size(), 1.);
686  }
687 
688  std::vector<double> vals;
689  for (auto it = fSignalParameters.begin();
690  it != fSignalParameters.end(); ++it) {
691  for (size_t i = 0; i < it->second.size(); ++i) {
692  vals.push_back(it->second[i]);
693  }
694  }
695  for (auto it = fFluxParameters.begin();
696  it != fFluxParameters.end(); ++it) {
697  vals.push_back(it->second);
698  }
699 
700  for (auto it = fSystParameters.begin();
701  it != fSystParameters.end(); ++it) {
702  vals.push_back(it->second.GetValue());
703  }
705  fUseFakeSamples = true;
706  fFitFunction(&vals[0]);
708  BuildFakeDataXSecs(false);
709 
710  fUseFakeSamples = false;
711  //Refill the hists for comparisons
712  fFitFunction(&vals[0]);
713  fFillIncidentInFunction = false;
714  }
715  else {
716  for (auto it = fFakeSamples.begin(); it != fFakeSamples.end(); ++it) {
717  for (size_t i = 0; i < it->second.size(); ++i) {
718  for (size_t j = 0; j < it->second[i].size(); ++j) {
719  it->second[i][j].Reset();
720  }
721  }
722  }
723  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
724  std::vector<ThinSliceSample> & samples_vec = it->second[0];
725  fFakeDataScales[it->first] = std::vector<double>(samples_vec.size(), 0.);
726  }
729  fSplitVal);
730  BuildFakeDataXSecs(false);
731  }
732 
733  inputFile.Close();
734 
735  fOutputFile.cd();
736  if (fDoScaleDataToNorm) {
737  ScaleDataToNorm();
738  }
739  TDirectory * out = (TDirectory *)fOutputFile.mkdir("Data");
740  out->cd();
741  fDataSet.GetIncidentHist().Write();
742  std::map<int, TH1 *> & selected_hists = fDataSet.GetSelectionHists();
743  for (auto it = selected_hists.begin(); it != selected_hists.end(); ++it) {
744  it->second->Write();
745  }
746 
748 }
virtual void BuildDataHists(TTree *tree, ThinSliceDataSet &data_set, double &flux, int split_val=0)=0
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::map< std::string, ThinSliceSystematic > fSystParameters
std::vector< double > fBeamEnergyBins
std::map< int, double > fFluxParameters
std::vector< double > fIncidentRecoBins
std::map< int, TH1 * > & GetSelectionHists()
std::map< int, std::vector< double > > fSignalParameters
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
void FillHistsFromSamples(const std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, double &flux)
static QFile inputFile
std::vector< fhicl::ParameterSet > fSelectionSets
std::map< int, std::vector< double > > fFakeDataScales
std::map< int, bool > fIsSignalSample
void BuildFakeDataXSecs(bool use_scales=true)
QTextStream & endl(QTextStream &s)
virtual void BuildFakeData(TTree *tree, std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, const std::map< int, bool > &signal_sample_checks, ThinSliceDataSet &data_set, double &flux, std::map< int, std::vector< double >> &sample_scales, std::vector< double > &beam_energy_bins, int split_val=0)=0
void protoana::PDSPThinSliceFitter::BuildFakeDataXSecs ( bool  use_scales = true)
private

Definition at line 2654 of file PDSPThinSliceFitter.cxx.

2654  {
2655  //First, set all samples to fake data scales
2656  if (use_scales) {
2657  for (auto it = fFakeSamples.begin(); it != fFakeSamples.end(); ++it) {
2658  for (size_t i = 0; i < it->second.size(); ++i) {
2659  for (size_t j = 0; j < it->second[i].size(); ++j) {
2660  std::cout << "Fake scale " << it->first << " " << i << " " << j << " " << fFakeDataScales[it->first][j] << std::endl;
2661  it->second[i][j].SetFactorAndScale(fFakeDataScales[it->first][j]);
2662  }
2663  }
2664  }
2665  }
2666 
2667  std::map<int, TH1D*> fake_data_totals;
2668 
2669  for (auto s : fMeasurementSamples) {
2670  //auto & samples_vec_2D = fSamples[s];
2671  auto & samples_vec_2D = fFakeSamples[s];
2672  std::vector<double> & bins = fSignalBins[s];
2673  std::string xsec_name = "FakeData" +
2674  samples_vec_2D[0][1].GetName() + "XSec";
2675  TH1D * temp_xsec = new TH1D(xsec_name.c_str(), "",
2676  fSignalBins[s].size() - 1, &bins[0]);
2677  for (size_t i = 0; i < samples_vec_2D.size(); ++i) {
2678  for (size_t j = 1; j < samples_vec_2D[i].size() - 1; ++j) {
2679  temp_xsec->AddBinContent(j, samples_vec_2D[i][j].GetVariedFlux());
2680  //temp_xsec->AddBinContent(j, samples_vec_2D[i][j].GetVariedFlux()*
2681  // fFakeDataScales[s][j]);
2682  }
2683  }
2684 
2685  std::string inc_name = "FakeData" +
2686  samples_vec_2D[0][1].GetName() + "Inc";
2687  TH1D * temp_inc = new TH1D(inc_name.c_str(), "", fSignalBins[s].size() - 1,
2688  &fSignalBins[s][0]);
2689  for (size_t i = 0; i < fIncidentSamples.size(); ++i) {
2690  //auto & vec_2D = fSamples[fIncidentSamples[i]];
2691  auto & vec_2D = fFakeSamples[fIncidentSamples[i]];
2692  for (size_t j = 0; j < vec_2D.size(); ++j) {
2693  auto & samples_vec = vec_2D[j];
2694  for (size_t k = 0; k < samples_vec.size(); ++k) {
2695  if (fSliceMethod == "E") {
2696  samples_vec[k].FillESliceHist(*temp_inc);
2697  }
2698  else {
2699  samples_vec[k].FillHistFromIncidentEnergies(*temp_inc);
2700  }
2701  }
2702  }
2703  }
2704 
2705  std::string tot_name = "FakeData" +
2706  samples_vec_2D[0][1].GetName() + "Tot";
2707  fake_data_totals[s] = (TH1D*)temp_xsec->Clone(tot_name.c_str());
2708  temp_xsec->Divide(temp_inc);
2709 
2710  if (fSliceMethod == "E") {
2711  for (int i = 1; i <= temp_xsec->GetNbinsX(); ++i) {
2712  temp_xsec->SetBinContent(i, -1.*log(1. - temp_xsec->GetBinContent(i)));
2713  }
2714  temp_xsec->Scale(1.E27*39.948/(1.4 * 6.022E23));
2715  for (int i = 1; i <= temp_xsec->GetNbinsX(); ++i) {
2716 
2717  double bethe_val = BetheBloch(temp_xsec->GetBinCenter(i), 139.57);
2718 
2719  temp_xsec->SetBinContent(i, (bethe_val*
2720  temp_xsec->GetBinContent(i)/
2721  temp_xsec->GetBinWidth(i)));
2722  }
2723  }
2724  else if (fSliceMethod == "Alt") {
2725  for (int i = 1; i <= temp_xsec->GetNbinsX(); ++i) {
2726  temp_xsec->SetBinContent(i, -1.*log(1. - temp_xsec->GetBinContent(i)));
2727  }
2728  temp_xsec->Scale(1.E27*39.948/(fPitch * 1.4 * 6.022E23));
2729  }
2730  else {
2731  for (int i = 1; i <= temp_xsec->GetNbinsX(); ++i) {
2732  temp_xsec->SetBinContent(i, -1.*log(1. - temp_xsec->GetBinContent(i)));
2733  }
2734  temp_xsec->Scale(1.E27/ (fPitch * 1.4 * 6.022E23 / 39.948 ));
2735  }
2736  fFakeDataXSecs[s] = temp_xsec;
2737  fFakeDataIncs[s] = temp_inc;
2738 
2739  fFakeDataXSecs[s]->SetDirectory(0);
2740  fFakeDataIncs[s]->SetDirectory(0);
2741  fake_data_totals[s]->SetDirectory(0);
2742 
2743  //std::cout << "Xsec vals " << fFakeDataXSecs[s] << std::endl;
2744  //for (int i = 1; i <= temp_xsec->GetNbinsX(); ++i) {
2745  // std::cout << temp_xsec->GetBinContent(i) << " ";
2746  //}
2747  //std::cout << std::endl;
2748  }
2749  fOutputFile.cd();
2750  TDirectory * out = (TDirectory *)fOutputFile.mkdir("FakeDataXSecs");
2751  out->cd();
2752  for (auto it = fFakeDataXSecs.begin(); it != fFakeDataXSecs.end(); ++it) {
2753  it->second->Write();
2754  fFakeDataIncs[it->first]->Write();
2755  fake_data_totals[it->first]->Write();
2756  }
2757 
2758 }
std::vector< int > fMeasurementSamples
std::map< int, TH1 * > fFakeDataXSecs
std::string string
Definition: nybbler.cc:12
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
std::map< int, TH1 * > fFakeDataIncs
double BetheBloch(double energy, double mass)
std::map< int, std::vector< double > > fSignalBins
std::map< int, std::vector< double > > fFakeDataScales
static QCString * s
Definition: config.cpp:1042
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::BuildMCSamples ( )

Definition at line 318 of file PDSPThinSliceFitter.cxx.

318  {
319  //Open the MC file and set branches
320  TFile fMCFile(fMCFileName.c_str(), "OPEN");
321  fMCTree = (TTree*)fMCFile.Get(fTreeName.c_str());
322 
323  //FillMCEvents();
325  fSplitMC);
332  fMCFile.Close();
333 
334  if (fDoSysts) {
337  fOutputFile);
338  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
339  for (size_t i = 0; i < it->second.size(); ++i) {
340  for (size_t j = 0; j < it->second[i].size(); ++j) {
341  fFakeSamples[it->first][i][j].SetSystematicSplines(
342  it->second[i][j].GetAllSplines());
343  }
344  }
345  }
346  }
347 
348 }
std::vector< ThinSliceEvent > fEvents
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
virtual void BuildMCSamples(const std::vector< ThinSliceEvent > &events, std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, const std::map< int, bool > &signal_sample_checks, std::map< int, double > &nominal_fluxes, std::map< int, std::vector< std::vector< double >>> &fluxes_by_sample, std::vector< double > &beam_energy_bins)=0
virtual void FillMCEvents(TTree *tree, std::vector< ThinSliceEvent > &events, std::vector< ThinSliceEvent > &fake_data_events, int &split_val, const bool &do_split)=0
std::map< std::string, ThinSliceSystematic > fSystParameters
std::vector< double > fBeamEnergyBins
std::vector< ThinSliceEvent > fFakeDataEvents
std::map< int, std::vector< std::vector< double > > > fFakeFluxesBySample
std::map< int, double > fFakeFluxes
virtual void SetupSysts(const std::vector< ThinSliceEvent > &events, std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, const std::map< int, bool > &signal_sample_checks, std::vector< double > &beam_energy_bins, const std::map< std::string, ThinSliceSystematic > &pars, TFile &output_file)=0
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
std::map< int, std::vector< std::vector< double > > > fFluxesBySample
std::map< int, bool > fIsSignalSample
std::map< int, double > fNominalFluxes
double protoana::PDSPThinSliceFitter::CalcChi2SystTerm ( )
private

Definition at line 2760 of file PDSPThinSliceFitter.cxx.

2760  {
2761  double result = 0.;
2762  for (auto it = fSystParameters.begin(); it != fSystParameters.end(); ++it) {
2763  for (auto it2 = fSystParameters.begin();
2764  it2 != fSystParameters.end(); ++it2) {
2765  double val_1 = it->second.GetValue();
2766  double val_2 = it2->second.GetValue();
2767  double central_1 = it->second.GetCentral();
2768  double central_2 = it2->second.GetCentral();
2769  int bin_1 = fCovarianceBins[it->first];
2770  int bin_2 = fCovarianceBins[it2->first];
2771 
2772  result += (val_1 - central_1)*(val_2 - central_2)*
2773  (*fCovMatrix)[bin_1][bin_2];
2774  }
2775  }
2776  return result;
2777 }
static QCString result
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< std::string, size_t > fCovarianceBins
double protoana::PDSPThinSliceFitter::CalcRegTerm ( )
private

Definition at line 2801 of file PDSPThinSliceFitter.cxx.

2801  {
2802  double result = 0.;
2803  for (auto it = fSignalParameters.begin(); it != fSignalParameters.end(); ++it) {
2804  for (size_t i = 0; i < it->second.size()-1; ++i) {
2805  result += std::pow((it->second[i] - it->second[i+1]), 2);
2806  }
2807  }
2808 
2809  return fRegFactor*result;
2810 }
static QCString result
constexpr T pow(T x)
Definition: pow.h:72
std::map< int, std::vector< double > > fSignalParameters
void protoana::PDSPThinSliceFitter::CompareDataMC ( std::string  extra_name,
TDirectory *  xsec_dir,
TDirectory *  plot_dir,
bool  post_fit = false 
)

Definition at line 769 of file PDSPThinSliceFitter.cxx.

771  {
776  plot_dir, fPlotRebinned, post_fit);
777 
778  xsec_dir->cd();
779 
780  //Get the incident histogram from all of the relevant samples
781  for (size_t i = 0; i < fMeasurementSamples.size(); ++i) {
782  int sample_ID = fMeasurementSamples[i];
783  auto & samples_vec_2D
784  = fSamples[sample_ID];
785 
786  std::vector<double> signal_bins = fSignalBins[sample_ID];
787  if (fDrawXSecUnderflow) signal_bins.insert(signal_bins.begin(), 0.);
788 
789  std::string signal_name = extra_name;
790  signal_name += samples_vec_2D[0][1].GetName() + "Hist";
791  TH1D signal_hist(signal_name.c_str(), "", signal_bins.size() - 1,
792  &signal_bins[0]);
793 
794  for (size_t j = 0; j < samples_vec_2D.size(); ++j) {
795  auto & samples_vec = samples_vec_2D[j];
796  for (size_t k = (fDrawXSecUnderflow? 0 : 1);
797  k < samples_vec.size()-1; ++k) {
798  ThinSliceSample & sample = samples_vec[k];
799  if (fDrawXSecUnderflow) {
800  signal_hist.AddBinContent(k+1, sample.GetVariedFlux());
801  }
802  else {
803  signal_hist.AddBinContent(k, sample.GetVariedFlux());
804  }
805  }
806  }
807  signal_hist.Write();
808 
809  //Get the incident histogram from all of the relevant samples
810  std::string inc_name = extra_name;
811  inc_name += "TotalIncident" + samples_vec_2D[0][0].GetName();
812 
813  TH1D * total_incident_hist = new TH1D(inc_name.c_str(), "",
814  signal_bins.size() - 1,
815  &signal_bins[0]);
816  std::map<int, std::vector<TH1D*>> temp_hists;
817  for (size_t i = 0; i < fIncidentSamples.size(); ++i) {
818  auto & vec_2D = fSamples[fIncidentSamples[i]];
819  temp_hists[fIncidentSamples[i]] = std::vector<TH1D*>();
820  for (size_t j = 0; j < vec_2D.size(); ++j) {
821  auto & samples_vec = vec_2D[j];
822  for (size_t k = 0; k < samples_vec.size(); ++k) {
823  if (j == 0) {
824  std::string name = extra_name;
825  name += "Incident" + samples_vec_2D[0][1].GetName() +
826  samples_vec[k].GetName() + std::to_string(k);
827  std::string title = samples_vec[k].GetSelectionHists().begin()->second->GetTitle();
828  temp_hists[fIncidentSamples[i]].push_back(
829  new TH1D(name.c_str(),
830  title.c_str(),
831  signal_bins.size() - 1,
832  &signal_bins[0]));
833  }
834  if (fSliceMethod == "E") {
835  samples_vec[k].FillESliceHist(*total_incident_hist);
836  samples_vec[k].FillESliceHist(
837  *(temp_hists[fIncidentSamples[i]][k]));
838  }
839  else {
840  samples_vec[k].FillHistFromIncidentEnergies(*total_incident_hist);
841  samples_vec[k].FillHistFromIncidentEnergies(
842  *(temp_hists[fIncidentSamples[i]][k]));
843  }
844  }
845  }
846  }
847  total_incident_hist->Write();
848  if (post_fit) {
849  fBestFitIncs[sample_ID] = total_incident_hist;
850  }
851  else {
852  fNominalIncs[sample_ID] = total_incident_hist;
853  }
854 
855  std::string xsec_name = extra_name + samples_vec_2D[0][1].GetName() +
856  "XSec";
857  TH1D * xsec_hist = (TH1D*)signal_hist.Clone(xsec_name.c_str());
858  xsec_hist->Sumw2();
859  xsec_hist->Divide(total_incident_hist);
860  if (fSliceMethod == "E") {
861  for (int i = 1; i <= xsec_hist->GetNbinsX(); ++i) {
862  xsec_hist->SetBinContent(i, -1.*log(1. - xsec_hist->GetBinContent(i)));
863  }
864  xsec_hist->Scale(1.E27*39.948/(1.4 * 6.022E23));
865  for (int i = 1; i <= xsec_hist->GetNbinsX(); ++i) {
866 
867  double bethe_val = BetheBloch(xsec_hist->GetBinCenter(i), 139.57);
868 
869  xsec_hist->SetBinContent(i, (bethe_val*
870  xsec_hist->GetBinContent(i)/
871  xsec_hist->GetBinWidth(i)));
872  }
873  }
874  else if (fSliceMethod == "Alt") {
875  for (int i = 1; i <= xsec_hist->GetNbinsX(); ++i) {
876  xsec_hist->SetBinContent(i, -1.*log(1. - xsec_hist->GetBinContent(i)));
877  }
878  xsec_hist->Scale(1.E27*39.948/(fPitch * 1.4 * 6.022E23));
879  }
880  else {
881  for (int i = 1; i <= xsec_hist->GetNbinsX(); ++i) {
882  xsec_hist->SetBinContent(i, -1.*log(1. - xsec_hist->GetBinContent(i)));
883  }
884  xsec_hist->Scale(1.E27/ (fPitch * 1.4 * 6.022E23 / 39.948 ));
885  }
886  xsec_hist->Write();
887  if (post_fit) {
888  fBestFitXSecs[sample_ID] = xsec_hist;
889  }
890  else {
891  fNominalXSecs[sample_ID] = xsec_hist;
892  }
893 
894  std::string stack_name = extra_name;
895  stack_name += "IncidentStack" + samples_vec_2D[0][1].GetName();
896  THStack inc_stack(stack_name.c_str(), ";Reconstructed KE (MeV)");
897  int iColor = 0;
898  for (auto it = temp_hists.begin(); it != temp_hists.end(); ++it) {
899  for (size_t i = 0; i < it->second.size(); ++i) {
900  std::pair<int, int> color_fill =
902  it->second[i]->SetFillColor(color_fill.first);
903  it->second[i]->SetFillStyle(color_fill.second);
904  it->second[i]->SetLineColor(1);
905  inc_stack.Add(it->second[i]);
906  ++iColor;
907  }
908  }
909  inc_stack.Write();
910  }
911 
912  //Get the Best fit selection hists
913  std::map<int, std::vector<TH1*>> temp_map;
915  fPlotRebinned);
916  for (auto it = temp_map.begin(); it != temp_map.end(); ++it) {
917  fBestFitSelectionHists[it->first] = (TH1D*)it->second[0];
918  }
919 
920  //Get the best fit truth vals
921  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
922  fBestFitTruthVals[it->first]
923  = std::vector<double>(it->second[0].size(), 0.);
924 
925  for (size_t i = 0; i < it->second[0].size(); ++i) {
926  double best_fit_val_i = 0.;
927  for (size_t j = 0; j < it->second.size(); ++j) {
928  best_fit_val_i += it->second[j][i].GetVariedFlux();
929  }
930 
931  fBestFitTruthVals[it->first][i] = best_fit_val_i;
932  }
933  }
934 }
static QCString name
Definition: declinfo.cpp:673
std::vector< int > fMeasurementSamples
std::pair< int, int > GetColorAndStyle(size_t i, const std::vector< std::pair< int, int >> &plot_style)
std::map< int, TH1 * > fBestFitXSecs
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
void CompareDataMC(std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, ThinSliceDataSet &data_set, TFile &output_file, std::vector< std::pair< int, int >> plot_style, int nPars, TDirectory *plot_dir, bool plot_rebinned=false, bool post_fit=false)
std::string string
Definition: nybbler.cc:12
std::vector< std::pair< int, int > > fPlotStyle
std::map< int, TH1D * > fBestFitSelectionHists
std::map< int, TH1 * > fNominalXSecs
std::map< int, std::vector< double > > fBestFitTruthVals
std::map< int, TH1 * > fBestFitIncs
std::map< int, TH1 * > fNominalIncs
virtual void GetCurrentHists(std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, ThinSliceDataSet &data_set, std::map< int, std::vector< TH1 * >> &hists, bool plot_rebinned)=0
double BetheBloch(double energy, double mass)
std::map< int, std::vector< double > > fSignalBins
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
void protoana::PDSPThinSliceFitter::Configure ( std::string  fcl_file)
private

Definition at line 1902 of file PDSPThinSliceFitter.cxx.

1902  {
1903  fhicl::ParameterSet pset;
1904 
1905  char const* fhicl_env = getenv("FHICL_FILE_PATH");
1906  std::string search_path;
1907 
1908  if (!fhicl_env) {
1909  std::cerr << "Expected environment variable FHICL_FILE_PATH is missing " <<
1910  "or empty: using \".\"\n";
1911  search_path = ".";
1912  }
1913  else {
1914  search_path = std::string{fhicl_env};
1915  }
1916 
1917  cet::filepath_first_absolute_or_lookup_with_dot lookupPolicy{search_path};
1918  pset = fhicl::ParameterSet::make(fcl_file, lookupPolicy);
1919 
1920  //Getting the configurable parameters
1921  fMCFileName = pset.get<std::string>("MCFileName");
1922  fDataFileName = pset.get<std::string>("DataFileName");
1923  fTreeName = pset.get<std::string>("TreeName");
1924 
1925  fSelectionSets = pset.get<std::vector<fhicl::ParameterSet>>("Selections");
1926 
1927  fIncidentRecoBins = pset.get<std::vector<double>>("IncidentRecoBins");
1928  fTrueIncidentBins = pset.get<std::vector<double>>("TrueIncidentBins");
1929  fBeamEnergyBins = pset.get<std::vector<double>>("BeamEnergyBins");
1930  fIncidentSamples = pset.get<std::vector<int>>("IncidentSamples");
1931  fMeasurementSamples = pset.get<std::vector<int>>("MeasurementSamples");
1932 
1933  fDrawXSecUnderflow = pset.get<bool>("DrawXSecUnderflow");
1934 
1935  fSampleSets = pset.get<std::vector<fhicl::ParameterSet>>("Samples");
1936  std::vector<std::pair<int, std::string>> temp_vec
1937  = pset.get<std::vector<std::pair<int, std::string>>>("FluxTypes");
1938  fFluxTypes = std::map<int, std::string>(temp_vec.begin(), temp_vec.end());
1939 
1940  fFitFlux = pset.get<bool>("FitFlux");
1941  if (fFitFlux) {
1942  for (size_t i = 0; i < temp_vec.size() - 1; ++i) {
1943  fFluxParameterNames[temp_vec[i].first] = "par_" + temp_vec[i].second +
1944  "_flux";
1945  fFluxParameters[temp_vec[i].first] = 1.;
1946  }
1947  fTotalFluxParameters = temp_vec.size() - 1;
1948  }
1949 
1950  fMaxCalls = pset.get<int>("MaxCalls", 1e9);
1951  fMaxIterations = pset.get<int>("MaxIterations", 1e6);
1952  fNScanSteps = pset.get<unsigned int>("NScanSteps") + 1;
1953  fTolerance = pset.get<double>("Tolerance");
1954  fLowerLimit = pset.get<double>("LowerLimit");
1955  fUpperLimit = pset.get<double>("UpperLimit");
1956  fPlotStyle = pset.get<std::vector<std::pair<int,int>>>("PlotStyle");
1957  fPlotRebinned = pset.get<bool>("PlotRebinned");
1958  fRandomStart = pset.get<bool>("RandomStart");
1959 
1960  fDriverName = pset.get<std::string>("DriverName");
1961  fAnalysisOptions = pset.get<fhicl::ParameterSet>("AnalysisOptions");
1962  fPitch = fAnalysisOptions.get<double>("WirePitch");
1963  fSliceMethod = fAnalysisOptions.get<std::string>("SliceMethod");
1964  fMultinomial = fAnalysisOptions.get<bool>("Multinomial", true);
1965  fDoFakeData = pset.get<bool>("DoFakeData");
1966  if (fDoFakeData) {
1967  fFakeDataRoutine = fAnalysisOptions.get<std::string>("FakeDataRoutine");
1968  }
1969  fSplitMC = pset.get<bool>("SplitMC");
1970  fDoThrows = pset.get<bool>("DoThrows");
1971  fDoScans = pset.get<bool>("DoScans");
1972  fRunHesse = pset.get<bool>("RunHesse");
1973  fDo1DShifts = pset.get<bool>("Do1DShifts");
1974  fDoSysts = pset.get<bool>("DoSysts");
1975  fFixVariables = pset.get<bool>("FixVariables", false);
1976  fFitUnderOverflow = pset.get<bool>("FitUnderOverflow", false);
1977  if (fFixVariables) {
1978  std::vector<std::pair<std::string, double>> temp_vec
1979  = pset.get<std::vector<std::pair<std::string, double>>>("SystsToFix");
1980  fSystsToFix = std::map<std::string, double>(temp_vec.begin(), temp_vec.end());
1981  }
1982  fDoFluctuateStats = pset.get<bool>("FluctuateStats");
1983 
1984  fNThrows = pset.get<size_t>("NThrows");
1985  fMaxRethrows = pset.get<size_t>("MaxRethrows");
1986 
1987  fFitType = pset.get<std::string>("FitType");
1988  fNPulls = pset.get<size_t>("NPulls");
1989 
1990  fDoScaleDataToNorm = pset.get<bool>("ScaleDataToNorm", false);
1991  fDataNorm = pset.get<double>("DataNorm", 1.);
1992 
1993  fRNG = TRandom3(pset.get<int>("RNGSeed", 0));
1994 
1995  fAddRegTerm = pset.get<bool>("AddRegTerm", false);
1996  fRegFactor = pset.get<double>("RegFactor", 0.);
1997 
1998  //Initialize systematics
1999  if (fDoSysts) {
2000  fAddSystTerm = pset.get<bool>("AddSystTerm");
2001  std::vector<fhicl::ParameterSet> par_vec
2002  = pset.get<std::vector<fhicl::ParameterSet>>("Systematics");
2003  for (size_t i = 0; i < par_vec.size(); ++i) {
2004  ThinSliceSystematic syst(par_vec[i]);
2005  fSystParameters[par_vec[i].get<std::string>("Name")] = syst;
2007  }
2008 
2009  if (fAddSystTerm) {
2010  std::vector<std::pair<std::string, int>> temp_vec
2011  = pset.get<std::vector<std::pair<std::string, int>>>("CovarianceBins");
2013  = std::map<std::string, size_t>(temp_vec.begin(), temp_vec.end());
2014 
2015  TFile cov_file(pset.get<std::string>("CovarianceFile").c_str());
2016  fCovMatrix = (TMatrixD*)cov_file.Get(
2017  pset.get<std::string>("CovarianceMatrix").c_str());
2018  fCovMatrixDisplay = (TMatrixD*)fCovMatrix->Clone();
2019  if (!fCovMatrix->IsSymmetric()) {
2020  std::string message = "PDSPThinSliceFitter::Configure: ";
2021  message += "Error. Input covariance matrix ";
2022  message += "is not symmetric";
2023  throw std::runtime_error(message);
2024  }
2025 
2026  fInputChol = new TDecompChol(*fCovMatrix);
2027  fInputChol->Decompose();
2028 
2029  fCovMatrix->Invert();
2030  }
2031  }
2032 }
std::vector< int > fMeasurementSamples
std::vector< double > fTrueIncidentBins
std::string string
Definition: nybbler.cc:12
static ParameterSet make(intermediate_table const &tbl)
Definition: ParameterSet.cc:68
std::map< std::string, ThinSliceSystematic > fSystParameters
std::vector< double > fBeamEnergyBins
std::map< int, double > fFluxParameters
std::vector< std::pair< int, int > > fPlotStyle
std::vector< double > fIncidentRecoBins
std::vector< fhicl::ParameterSet > fSampleSets
std::map< std::string, size_t > fCovarianceBins
std::string getenv(std::string const &name)
Definition: getenv.cc:15
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::map< int, std::string > fFluxParameterNames
std::vector< fhicl::ParameterSet > fSelectionSets
std::map< std::string, double > fSystsToFix
std::map< int, std::string > fFluxTypes
fhicl::ParameterSet fAnalysisOptions
void protoana::PDSPThinSliceFitter::DefineFitFunction ( )
private

use samples

Definition at line 1742 of file PDSPThinSliceFitter.cxx.

1742  {
1743  ///use samples
1744  fFitFunction = ROOT::Math::Functor(
1745  [&](double const * coeffs) {
1746 
1747  //Set all the parameters
1748  //coeffs are ordered according to
1749  //keys of signal par map
1750  //----vector for given signal
1751  //then keys of flux par map
1752  //then values of syst parameters
1753  size_t par_position = 0;
1754  for (auto it = fSignalParameters.begin();
1755  it != fSignalParameters.end(); ++it) {
1756  for (size_t i = 0; i < it->second.size(); ++i) {
1757  it->second.at(i) = coeffs[par_position];
1758  ++par_position;
1759  }
1760  }
1761  for (auto it = fFluxParameters.begin();
1762  it != fFluxParameters.end(); ++it) {
1763  it->second = coeffs[par_position];
1764  ++par_position;
1765  }
1766 
1767  for (auto it = fSystParameters.begin();
1768  it != fSystParameters.end(); ++it) {
1769  it->second.SetValue(coeffs[par_position]);
1770  //std::cout << it->first << " " << it->second.GetValue() << std::endl;
1771  ++par_position;
1772  }
1773 
1774  //Refill the hists
1775  if (!fUseFakeSamples) {
1781  }
1782  else {
1788  }
1789 
1790  //Scale to Data
1791  if (!fUseFakeSamples) {
1792  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
1793  for (size_t i = 0; i < it->second.size(); ++i) {
1794  for (size_t j = 0; j < it->second[i].size(); ++j) {
1795  it->second[i][j].ScaleToDataMC();
1796  }
1797  }
1798  }
1799  }
1800 
1801 
1802  double total_varied_flux = 0.;
1803  double total_nominal_flux = 0.;
1804  std::vector<double> varied_flux(fBeamEnergyBins.size() - 1, 0.);
1805  std::vector<double> nominal_flux(fBeamEnergyBins.size() - 1, 0.);
1806 
1807  //Go through and determine the scaled flux
1808  for (auto it = fFluxesBySample.begin();
1809  it != fFluxesBySample.end(); ++it) {
1810  int sample_ID = it->first;
1811  std::vector<std::vector<ThinSliceSample>> & samples_2D
1812  = fSamples[sample_ID];
1813  for (size_t i = 0; i < samples_2D.size(); ++i) {
1814  std::vector<ThinSliceSample> & samples = samples_2D[i];
1815  for (size_t j = 0; j < samples.size(); ++j) {
1816  ThinSliceSample & sample = samples.at(j);
1817  int flux_type = sample.GetFluxType();
1818  if (flux_type == 1) {
1819  nominal_flux[i] += (!fUseFakeSamples ?
1820  fFluxesBySample[sample_ID][i][j] :
1821  fFakeFluxesBySample[sample_ID][i][j]);
1822  varied_flux[i] += sample.GetVariedFlux();
1823  }
1824  total_nominal_flux += (!fUseFakeSamples ?
1825  fFluxesBySample[sample_ID][i][j] :
1826  fFakeFluxesBySample[sample_ID][i][j]);
1827  total_varied_flux += sample.GetVariedFlux();
1828 
1829  }
1830  }
1831  }
1832 
1833  //double total_flux_factor
1834  // = (fMultinomial ? (total_nominal_flux/total_varied_flux) : 1.);
1835  double total_flux_factor = (total_nominal_flux/total_varied_flux);
1836 
1837  double nominal_primary = 0., varied_primary = 0.;
1838  for (size_t i = 0; i < nominal_flux.size(); ++i) {
1839  nominal_primary += nominal_flux[i];
1840  varied_primary += varied_flux[i];
1841  }
1842 
1843  std::vector<double> flux_factor = nominal_flux;
1844  for (size_t i = 0; i < flux_factor.size(); ++i) {
1845  flux_factor[i] = (nominal_flux[i] > 1.e-7 ?
1846  flux_factor[i] / varied_flux[i] :
1847  0.)*(varied_primary/nominal_primary);
1848  }
1849 
1850  if (!fUseFakeSamples) {
1851  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
1852  for (size_t i = 0; i < it->second.size(); ++i) {
1853  for (size_t j = 0; j < it->second[i].size(); ++j) {
1854  int flux_type = it->second[i][j].GetFluxType();
1855  //double total = 0.;
1856  //std::cout << "Varied flux1 : " << std::setprecision(8) <<it->second[i][j].GetVariedFlux() << std::endl;
1857  //std::cout << "scaling by " << total_flux_factor << " " << (flux_type == 1 ? flux_factor[i] : 1.) <<
1858  // " " << total_flux_factor*(flux_type == 1 ? flux_factor[i] : 1.) << std::endl;
1859 
1860  it->second[i][j].SetFactorAndScale(
1861  total_flux_factor*(flux_type == 1 ? flux_factor[i] : 1.));
1862  //std::cout << "Varied flux2 : " << std::setprecision(8) <<it->second[i][j].GetVariedFlux() << std::endl;
1863  }
1864  }
1865  }
1866  }
1867  else {
1868  for (auto it = fFakeSamples.begin(); it != fFakeSamples.end(); ++it) {
1869  for (size_t i = 0; i < it->second.size(); ++i) {
1870  for (size_t j = 0; j < it->second[i].size(); ++j) {
1871  int flux_type = it->second[i][j].GetFluxType();
1872  //double total = 0.;
1873  //std::cout << "Varied flux1 : " << std::setprecision(8) <<it->second[i][j].GetVariedFlux() << std::endl;
1874  //std::cout << "scaling by " << total_flux_factor << " " << (flux_type == 1 ? flux_factor[i] : 1.) <<
1875  // " " << total_flux_factor*(flux_type == 1 ? flux_factor[i] : 1.) << std::endl;
1876  it->second[i][j].SetFactorAndScale(
1877  total_flux_factor*(flux_type == 1 ? flux_factor[i] : 1.));
1878  //std::cout << "Varied flux2 : " << std::setprecision(8) <<it->second[i][j].GetVariedFlux() << std::endl;
1879  }
1880  }
1881  }
1882  }
1883 
1884  std::pair<double, size_t> chi2_points
1886  ++fNFitSteps;
1887 
1888  double syst_chi2 = (fAddSystTerm ? CalcChi2SystTerm() : 0.);
1889 
1890  //Do I need to do this in ThinSliceDriver?
1891  //double reg_term = (fAddRegTerm ? CalcRegTerm() : 0.);
1892  //
1893 
1894  //std::cout << (chi2_points.first + syst_chi2) << std::endl;
1895  return (chi2_points.first + syst_chi2 /* + reg_term */);
1896  },
1898  std::cout << "Done F2" << std::endl;
1899 
1900 }
std::vector< ThinSliceEvent > fEvents
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
virtual void RefillMCSamples(const std::vector< ThinSliceEvent > &events, std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, const std::map< int, bool > &signal_sample_checks, std::vector< double > &beam_energy_bins, const std::map< int, std::vector< double >> &signal_pars, const std::map< int, double > &flux_pars, const std::map< std::string, ThinSliceSystematic > &syst_pars, bool fit_under_over, bool fill_incident=false)=0
std::map< std::string, ThinSliceSystematic > fSystParameters
std::vector< double > fBeamEnergyBins
std::map< int, double > fFluxParameters
std::vector< ThinSliceEvent > fFakeDataEvents
std::map< int, std::vector< std::vector< double > > > fFakeFluxesBySample
virtual std::pair< double, size_t > CalculateChi2(std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, ThinSliceDataSet &data_set)=0
std::map< int, std::vector< double > > fSignalParameters
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
std::map< int, std::vector< std::vector< double > > > fFluxesBySample
std::map< int, bool > fIsSignalSample
QTextStream & endl(QTextStream &s)
double protoana::PDSPThinSliceFitter::densityEffect ( double  beta,
double  gamma 
)
inlineprivate

Definition at line 220 of file PDSPThinSliceFitter.h.

220  {
221  double lar_C = 5.215, lar_x0 = 0.201, lar_x1 = 3, lar_a = 0.196, lar_k = 3;
222  long double x = log10(beta * gamma);
223 
224  if( x >= lar_x1 ) return 2*log(10)*x - lar_C;
225 
226  else if ( lar_x0 <= x && x < lar_x1) return 2*log(10)*x - lar_C + lar_a * pow(( lar_x1 - x ) , lar_k );
227 
228  else return 0.; //if x < lar_x0
229  };
double beta(double KE, const simb::MCParticle *part)
constexpr T pow(T x)
Definition: pow.h:72
double gamma(double KE, const simb::MCParticle *part)
list x
Definition: train.py:276
void protoana::PDSPThinSliceFitter::Do1DShifts ( const TH1D &  pars,
bool  prefit = false 
)
private

Definition at line 1692 of file PDSPThinSliceFitter.cxx.

1692  {
1693  //std::string dir_name = (prefit ? "PreFit1DShifts" : "1DShifts");
1694  TDirectory * shift_dir = fOutputFile.mkdir((prefit ? "PreFit1DShifts" : "1DShifts"));
1695  fFillIncidentInFunction = true;
1696 
1697  //Iterate over the bins in the parameter hist
1698  for (int i = 1; i <= pars.GetNbinsX(); ++i) {
1699 
1700  std::string par_name = pars.GetXaxis()->GetBinLabel(i);
1701 
1702  std::vector<double> vals;
1703  if (pars.GetBinError(i) < 1.e-9) continue;
1704  //Set to +1 sigma for parameter i
1705  for (int j = 1; j <= pars.GetNbinsX(); ++j) {
1706  if (j == i) {
1707  vals.push_back(pars.GetBinContent(j) + pars.GetBinError(i));
1708  }
1709  else {
1710  vals.push_back(pars.GetBinContent(j));
1711  }
1712  }
1713  fFitFunction(&vals[0]);
1714 
1715  shift_dir->cd();
1716  std::string dir_name = par_name + "_PlusSigma";
1717  TDirectory * plus_dir = shift_dir->mkdir(dir_name.c_str());
1718  TDirectory * xsec_dir = plus_dir->mkdir("XSec");
1719  CompareDataMC(dir_name, xsec_dir, plus_dir);
1720 
1721  vals.clear();
1722 
1723  //Set to -1 sigma for parameter i
1724  for (int j = 1; j <= pars.GetNbinsX(); ++j) {
1725  if (j == i) {
1726  vals.push_back(pars.GetBinContent(j) - pars.GetBinError(i));
1727  }
1728  else {
1729  vals.push_back(pars.GetBinContent(j));
1730  }
1731  }
1732  fFitFunction(&vals[0]);
1733 
1734  shift_dir->cd();
1735  dir_name = par_name + "_MinusSigma";
1736  TDirectory * minus_dir = shift_dir->mkdir(dir_name.c_str());
1737  xsec_dir = minus_dir->mkdir("XSec");
1738  CompareDataMC(dir_name, xsec_dir, minus_dir);
1739  }
1740 }
std::string string
Definition: nybbler.cc:12
void CompareDataMC(std::string extra_name, TDirectory *xsec_dir, TDirectory *plot_dir, bool post_fit=false)
void protoana::PDSPThinSliceFitter::DoThrows ( const TH1D &  pars,
const TMatrixD *  cov 
)
private

Definition at line 1523 of file PDSPThinSliceFitter.cxx.

1523  {
1524  std::vector<std::vector<double>> vals;
1525  TVectorD rand(pars.GetNbinsX());
1526 
1527  TDecompChol chol(*cov);
1528  bool success = chol.Decompose();
1529  if (!success) {
1530  std::cout << "Error" << std::endl;
1531  return;
1532  }
1533 
1534  fOutputFile.mkdir("Throws");
1535  fOutputFile.cd("Throws");
1536 
1537  TTree throws_tree("tree", "");
1538  std::vector<double> throws_branches;
1539  std::vector<TVectorD *> throws_arrays;
1540  std::vector<TVectorD *> throws_xsec_arrays;
1541 
1542  MakeThrowsTree(throws_tree, throws_branches);
1543  //MakeThrowsArrays(throws_arrays);
1544  for (auto it = fSignalParameters.begin(); it != fSignalParameters.end(); ++it) {
1545  for (size_t i = 0; i < it->second.size(); ++i) {
1546  throws_arrays.push_back(new TVectorD(fNThrows));
1547  throws_xsec_arrays.push_back(new TVectorD(fNThrows));
1548  }
1549  }
1550 
1551  for (auto it = fFluxParameters.begin(); it != fFluxParameters.end(); ++it) {
1552  throws_arrays.push_back(new TVectorD(fNThrows));
1553  }
1554 
1555  for (auto it = fSystParameters.begin(); it != fSystParameters.end(); ++it) {
1556  throws_arrays.push_back(new TVectorD(fNThrows));
1557  }
1558 
1559 
1560 
1561  TCanvas cPars("cParametersThrown", "");
1562  cPars.SetTicks();
1563 
1564  //Build up the map of selection hists
1565  std::map<int, std::vector<TH1*>> throw_hists;
1566  for (auto it = fDataSet.GetSelectionHists().begin();
1567  it != fDataSet.GetSelectionHists().end(); ++it) {
1568  throw_hists[it->first] = std::vector<TH1*>();
1569  }
1570 
1571  std::map<int, std::vector<TH1*>> truth_throw_hists;
1572  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
1573  truth_throw_hists[it->first] = std::vector<TH1*>();
1574  }
1575 
1576  std::map<int, std::vector<TH1*>> truth_inc_throw_hists,
1577  truth_xsec_throw_hists;
1578  for (auto it = fMeasurementSamples.begin();
1579  it != fMeasurementSamples.end(); ++it) {
1580  truth_inc_throw_hists[*it] = std::vector<TH1*>();
1581  truth_xsec_throw_hists[*it] = std::vector<TH1*>();
1582  }
1583 
1584  /* for (auto it = fCombinedMeasurements.begin();
1585  it != fCombinedMeasurements.end(); ++it) {
1586  truth_combined_inc_throw_hists[*it] = std::vector<TH1*>();
1587  truth_combined_xsec_throw_hists[*it] = std::vector<TH1*>();
1588  } */
1589 
1590  //int n_signal_flux_pars = fTotalSignalParameters + fTotalFluxParameters;
1591  for (size_t i = 0; i < fNThrows; ++i) {
1592  if (! (i%100) ) {
1593  std::cout << "Throw " << i << "/" << fNThrows << std::endl;
1594  //auto start = std::chrono::high_resolution_clock::now();
1595  }
1596  vals.push_back(std::vector<double>(pars.GetNbinsX(), 1.));
1597 
1598  bool rethrow = true;
1599  size_t nRethrows = 0;
1600  while (rethrow && nRethrows < fMaxRethrows) {
1601  bool all_pos = true;
1602  for (int j = 1; j <= pars.GetNbinsX(); ++j) {
1603  rand[j-1] = fRNG.Gaus();
1604  }
1605 
1606  TVectorD rand_times_chol = chol.GetU()*rand;
1607 
1608  for (int j = 1; j <= pars.GetNbinsX(); ++j) {
1609  //std::cout << rand_times_chol[j-1] << " ";
1610  vals.back()[j-1] = pars.GetBinContent(j) + rand_times_chol[j-1];
1611 
1612  //Configure this to truncate at 0 or rethrow
1613  if (vals.back()[j-1] < fParLimits[j-1]) {
1614  all_pos = false;
1615  std::cout << "Rethrowing " << j-1 << " " << vals.back()[j-1] << " " << fParLimits[j-1] << std::endl;
1616  }
1617  if (vals.back()[j-1] > fParLimitsUp[j-1]) {
1618  all_pos = false;
1619  std::cout << "Rethrowing " << j-1 << " " << vals.back()[j-1] << " " << fParLimitsUp[j-1] << std::endl;
1620  }
1621  rethrow = !all_pos;
1622 
1623  }
1624  //std::cout << std::endl;
1625  ++nRethrows;
1626  }
1627 
1628  for (size_t j = 0; j < vals.back().size(); ++j) {
1629  throws_branches[j] = vals.back()[j];
1630  (*throws_arrays[j])[i] = vals.back()[j];
1631  }
1632  throws_tree.Fill();
1633 
1634  //Applies the variations according to the thrown parameters
1635  fFillIncidentInFunction = true;
1636  fFitFunction(&(vals.back()[0]));
1638  GetCurrentTruthHists(truth_throw_hists,
1639  truth_inc_throw_hists,
1640  truth_xsec_throw_hists);
1641  int xsec_bin = 0;
1642  for(auto it = truth_xsec_throw_hists.begin();
1643  it != truth_xsec_throw_hists.end(); ++it) {
1644  TH1 * xsec_hist = it->second.back();
1645  if (fSliceMethod == "E") {
1646  for (int j = 1; j <= xsec_hist->GetNbinsX(); ++j) {
1647  xsec_hist->SetBinContent(j, -1.*log(1. - xsec_hist->GetBinContent(j)));
1648  }
1649  xsec_hist->Scale(1.E27*39.948/(1.4 * 6.022E23));
1650  for (int j = 1; j <= xsec_hist->GetNbinsX(); ++j) {
1651 
1652  double bethe_val = BetheBloch(xsec_hist->GetBinCenter(j), 139.57);
1653 
1654  xsec_hist->SetBinContent(i, (bethe_val*
1655  xsec_hist->GetBinContent(j)/
1656  xsec_hist->GetBinWidth(j)));
1657  }
1658  }
1659  else if (fSliceMethod == "Alt") {
1660  for (int j = 1; j <= xsec_hist->GetNbinsX(); ++j) {
1661  xsec_hist->SetBinContent(j, -1.*log(1. - xsec_hist->GetBinContent(j)));
1662  }
1663  xsec_hist->Scale(1.E27*39.948/(fPitch * 1.4 * 6.022E23));
1664  }
1665  else {
1666  for (int j = 1; j <= xsec_hist->GetNbinsX(); ++j) {
1667  xsec_hist->SetBinContent(j, -1.*log(1. - xsec_hist->GetBinContent(j)));
1668  }
1669  xsec_hist->Scale(1.E27/ (fPitch * 1.4 * 6.022E23 / 39.948 ));
1670  }
1671  for (int j = 1; j <= xsec_hist->GetNbinsX(); ++j) {
1672  (*throws_xsec_arrays[xsec_bin])[i] = xsec_hist->GetBinContent(j);
1673  ++xsec_bin;
1674  }
1675  }
1676  }
1677 
1678  throws_tree.Write();
1679  for (size_t i = 0; i < throws_arrays.size(); ++i) {
1680  std::string name = "array" + std::to_string(i);
1681  throws_arrays[i]->Write(name.c_str());
1682  }
1683  for (size_t i = 0; i < throws_xsec_arrays.size(); ++i) {
1684  std::string name = "xsec_array" + std::to_string(i);
1685  throws_xsec_arrays[i]->Write(name.c_str());
1686  }
1687 
1688  PlotThrows(throw_hists, truth_throw_hists,
1689  truth_inc_throw_hists, truth_xsec_throw_hists);
1690 }
static QCString name
Definition: declinfo.cpp:673
std::vector< int > fMeasurementSamples
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::string string
Definition: nybbler.cc:12
std::map< std::string, ThinSliceSystematic > fSystParameters
void MakeThrowsTree(TTree &tree, std::vector< double > &branches)
std::map< int, double > fFluxParameters
void GetCurrentTruthHists(std::map< int, std::vector< TH1 * >> &throw_hists, std::map< int, std::vector< TH1 * >> &throw_inc_hists, std::map< int, std::vector< TH1 * >> &throw_xsec_hists)
std::vector< double > fParLimitsUp
std::map< int, TH1 * > & GetSelectionHists()
std::map< int, std::vector< double > > fSignalParameters
virtual void GetCurrentHists(std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, ThinSliceDataSet &data_set, std::map< int, std::vector< TH1 * >> &hists, bool plot_rebinned)=0
double BetheBloch(double energy, double mass)
void PlotThrows(std::map< int, std::vector< TH1 * >> &throw_hists, std::map< int, std::vector< TH1 * >> &truth_throw_hists, std::map< int, std::vector< TH1 * >> &truth_inc_hists, std::map< int, std::vector< TH1 * >> &truth_xsec_hists)
std::vector< double > fParLimits
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::FillMCEvents ( )

Definition at line 350 of file PDSPThinSliceFitter.cxx.

350  {
351 
352  std::cout << "Filling MC Events" << std::endl;
353 
354  int sample_ID, selection_ID, event, run, subrun;
355  int true_beam_PDG;
356  double true_beam_interactingEnergy, reco_beam_interactingEnergy;
357  double true_beam_endP, true_beam_mass, true_beam_endZ;
358  double reco_beam_endZ, true_beam_startP;
359  double beam_inst_P;
360  std::vector<double> * reco_beam_incidentEnergies = 0x0,
361  * true_beam_incidentEnergies = 0x0,
362  * true_beam_traj_Z = 0x0,
363  * true_beam_traj_KE = 0x0,
364  * reco_daughter_track_thetas = 0x0,
365  * reco_daughter_track_scores = 0x0;
366  std::vector<int> * true_beam_slices = 0x0;
367  fMCTree->SetBranchAddress("event", &event);
368  fMCTree->SetBranchAddress("subrun", &subrun);
369  fMCTree->SetBranchAddress("run", &run);
370 
371  fMCTree->SetBranchAddress("true_beam_PDG", &true_beam_PDG);
372 
373  fMCTree->SetBranchAddress("new_interaction_topology", &sample_ID);
374  fMCTree->SetBranchAddress("selection_ID", &selection_ID);
375  fMCTree->SetBranchAddress("true_beam_interactingEnergy",
376  &true_beam_interactingEnergy);
377  fMCTree->SetBranchAddress("true_beam_endP", &true_beam_endP);
378  fMCTree->SetBranchAddress("true_beam_endZ", &true_beam_endZ);
379  fMCTree->SetBranchAddress("true_beam_mass", &true_beam_mass);
380  fMCTree->SetBranchAddress("reco_beam_interactingEnergy",
381  &reco_beam_interactingEnergy);
382  fMCTree->SetBranchAddress("reco_beam_endZ", &reco_beam_endZ);
383  fMCTree->SetBranchAddress("reco_beam_incidentEnergies",
384  &reco_beam_incidentEnergies);
385  fMCTree->SetBranchAddress("true_beam_incidentEnergies",
386  &true_beam_incidentEnergies);
387  fMCTree->SetBranchAddress("true_beam_slices",
388  &true_beam_slices);
389  fMCTree->SetBranchAddress("true_beam_startP", &true_beam_startP);
390  fMCTree->SetBranchAddress("true_beam_traj_Z", &true_beam_traj_Z);
391  fMCTree->SetBranchAddress("true_beam_traj_KE", &true_beam_traj_KE);
392  std::vector<double> * calibrated_dQdX = 0x0, * beam_EField = 0x0,
393  * track_pitch = 0x0;
394  fMCTree->SetBranchAddress("reco_beam_calibrated_dQdX_SCE", &calibrated_dQdX);
395  fMCTree->SetBranchAddress("reco_beam_EField_SCE", &beam_EField);
396  fMCTree->SetBranchAddress("reco_beam_TrkPitch_SCE", &track_pitch);
397  fMCTree->SetBranchAddress("beam_inst_P", &beam_inst_P);
398  fMCTree->SetBranchAddress("reco_daughter_allTrack_Theta", &reco_daughter_track_thetas);
399  fMCTree->SetBranchAddress("reco_daughter_PFP_trackScore_collection",
400  &reco_daughter_track_scores);
401 
402  std::vector<double> * g4rw_alt_primary_plus_sigma_weight = 0x0,
403  * g4rw_alt_primary_minus_sigma_weight = 0x0,
404  * g4rw_full_primary_plus_sigma_weight = 0x0,
405  * g4rw_full_primary_minus_sigma_weight = 0x0;
406  std::vector<std::vector<double>> * g4rw_primary_grid_weights = 0x0,
407  * g4rw_full_grid_weights = 0x0,
408  * g4rw_full_grid_proton_weights = 0x0;
409  fMCTree->SetBranchAddress("g4rw_alt_primary_plus_sigma_weight",
410  &g4rw_alt_primary_plus_sigma_weight);
411  fMCTree->SetBranchAddress("g4rw_alt_primary_minus_sigma_weight",
412  &g4rw_alt_primary_minus_sigma_weight);
413  fMCTree->SetBranchAddress("g4rw_full_primary_plus_sigma_weight",
414  &g4rw_full_primary_plus_sigma_weight);
415  fMCTree->SetBranchAddress("g4rw_full_primary_minus_sigma_weight",
416  &g4rw_full_primary_minus_sigma_weight);
417  fMCTree->SetBranchAddress("g4rw_full_grid_weights", &g4rw_full_grid_weights);
418  fMCTree->SetBranchAddress("g4rw_full_grid_proton_weights", &g4rw_full_grid_proton_weights);
419 
420  fMCTree->SetBranchAddress("g4rw_primary_grid_weights",
421  &g4rw_primary_grid_weights);
422 
423  std::vector<std::vector<double>> * daughter_dQdXs = 0x0,
424  * daughter_resRanges = 0x0,
425  * daughter_EFields = 0x0;
426  fMCTree->SetBranchAddress(
427  "reco_daughter_allTrack_calibrated_dQdX_SCE", &daughter_dQdXs);
428  fMCTree->SetBranchAddress(
429  "reco_daughter_allTrack_resRange_SCE", &daughter_resRanges);
430  fMCTree->SetBranchAddress(
431  "reco_daughter_allTrack_EField_SCE", &daughter_EFields);
432  bool has_pi0_shower;
433  fMCTree->SetBranchAddress("has_shower_dist_energy", &has_pi0_shower);
434 
435  int nentries = fMCTree->GetEntries();
436  if (fSplitMC) {
437  fSplitVal = fMCTree->GetEntries()/2;
438  std::cout << "Note: Splitting MC in half. " <<
439  fSplitVal << "/" << fMCTree->GetEntries() <<std::endl;
440  nentries = fSplitVal;
441  }
442  for (int i = 0; i < nentries; ++i) {
443  fMCTree->GetEntry(i);
444 
445  fEvents.push_back(ThinSliceEvent(event, subrun, run));
446  fEvents.back().SetSampleID(sample_ID);
447  fEvents.back().SetSelectionID(selection_ID);
448  fEvents.back().SetTrueInteractingEnergy(true_beam_interactingEnergy);
449  fEvents.back().SetRecoInteractingEnergy(reco_beam_interactingEnergy);
450  fEvents.back().SetTrueEndP(true_beam_endP);
451  fEvents.back().SetTrueEndZ(true_beam_endZ);
452  fEvents.back().SetTrueStartP(true_beam_startP);
453  fEvents.back().SetTrueMass(true_beam_mass);
454  fEvents.back().SetRecoEndZ(reco_beam_endZ);
455 
456  fEvents.back().SetRecoIncidentEnergies(*reco_beam_incidentEnergies);
457  fEvents.back().SetTrueIncidentEnergies(*true_beam_incidentEnergies);
458  fEvents.back().SetTrueTrajZ(*true_beam_traj_Z);
459  fEvents.back().SetTrueTrajKE(*true_beam_traj_KE);
460  fEvents.back().SetTrueSlices(*true_beam_slices);
461  fEvents.back().SetdQdXCalibrated(*calibrated_dQdX);
462  fEvents.back().SetEField(*beam_EField);
463  fEvents.back().SetTrackPitch(*track_pitch);
464  fEvents.back().SetBeamInstP(beam_inst_P);
465  fEvents.back().SetPDG(true_beam_PDG);
466  fEvents.back().SetRecoDaughterTrackThetas(*reco_daughter_track_thetas);
467  fEvents.back().SetRecoDaughterTrackScores(*reco_daughter_track_scores);
468  fEvents.back().SetHasPi0Shower(has_pi0_shower);
469  fEvents.back().MakeG4RWBranch("g4rw_alt_primary_plus_sigma_weight",
470  *g4rw_alt_primary_plus_sigma_weight);
471  fEvents.back().MakeG4RWBranch("g4rw_alt_primary_minus_sigma_weight",
472  *g4rw_alt_primary_minus_sigma_weight);
473  fEvents.back().MakeG4RWBranch("g4rw_full_primary_plus_sigma_weight",
474  *g4rw_full_primary_plus_sigma_weight);
475  fEvents.back().MakeG4RWBranch("g4rw_full_primary_minus_sigma_weight",
476  *g4rw_full_primary_minus_sigma_weight);
477  for (size_t j = 0; j < daughter_dQdXs->size(); ++j) {
478  fEvents.back().AddRecoDaughterTrackdQdX((*daughter_dQdXs)[j]);
479  fEvents.back().AddRecoDaughterTrackResRange((*daughter_resRanges)[j]);
480  fEvents.back().AddRecoDaughterEField((*daughter_EFields)[j]);
481  }
482 
483  for (size_t j = 0; j < g4rw_primary_grid_weights->size(); ++j) {
484  std::string name_full = "g4rw_full_grid_weights_" + std::to_string(j);
485  fEvents.back().MakeG4RWBranch(name_full, (*g4rw_full_grid_weights)[j]);
486 
487  std::string name_primary = "g4rw_primary_grid_weights_" +
488  std::to_string(j);
489  fEvents.back().MakeG4RWBranch(name_primary,
490  (*g4rw_primary_grid_weights)[j]);
491  }
492  fEvents.back().MakeG4RWBranch("g4rw_full_grid_proton_weights",
493  (*g4rw_full_grid_proton_weights)[0]);
494  }
495 
496  for (int i = fSplitVal; i < fMCTree->GetEntries(); ++i) {
497  fMCTree->GetEntry(i);
498 
499  fFakeDataEvents.push_back(ThinSliceEvent(event, subrun, run));
500  fFakeDataEvents.back().SetSampleID(sample_ID);
501  fFakeDataEvents.back().SetSelectionID(selection_ID);
502  fFakeDataEvents.back().SetTrueInteractingEnergy(true_beam_interactingEnergy);
503  fFakeDataEvents.back().SetRecoInteractingEnergy(reco_beam_interactingEnergy);
504  fFakeDataEvents.back().SetTrueEndP(true_beam_endP);
505  fFakeDataEvents.back().SetTrueEndZ(true_beam_endZ);
506  fFakeDataEvents.back().SetTrueStartP(true_beam_startP);
507  fFakeDataEvents.back().SetTrueMass(true_beam_mass);
508  fFakeDataEvents.back().SetRecoEndZ(reco_beam_endZ);
509 
510  fFakeDataEvents.back().SetRecoIncidentEnergies(*reco_beam_incidentEnergies);
511  fFakeDataEvents.back().SetTrueIncidentEnergies(*true_beam_incidentEnergies);
512  fFakeDataEvents.back().SetTrueTrajZ(*true_beam_traj_Z);
513  fFakeDataEvents.back().SetTrueTrajKE(*true_beam_traj_KE);
514  fFakeDataEvents.back().SetTrueSlices(*true_beam_slices);
515  fFakeDataEvents.back().SetdQdXCalibrated(*calibrated_dQdX);
516  fFakeDataEvents.back().SetEField(*beam_EField);
517  fFakeDataEvents.back().SetTrackPitch(*track_pitch);
518  fFakeDataEvents.back().SetBeamInstP(beam_inst_P);
519  fFakeDataEvents.back().SetPDG(true_beam_PDG);
520  fFakeDataEvents.back().SetRecoDaughterTrackThetas(*reco_daughter_track_thetas);
521  fFakeDataEvents.back().SetRecoDaughterTrackScores(*reco_daughter_track_scores);
522  fFakeDataEvents.back().SetHasPi0Shower(has_pi0_shower);
523  fFakeDataEvents.back().MakeG4RWBranch("g4rw_alt_primary_plus_sigma_weight",
524  *g4rw_alt_primary_plus_sigma_weight);
525  fFakeDataEvents.back().MakeG4RWBranch("g4rw_alt_primary_minus_sigma_weight",
526  *g4rw_alt_primary_minus_sigma_weight);
527  fFakeDataEvents.back().MakeG4RWBranch("g4rw_full_primary_plus_sigma_weight",
528  *g4rw_full_primary_plus_sigma_weight);
529  fFakeDataEvents.back().MakeG4RWBranch("g4rw_full_primary_minus_sigma_weight",
530  *g4rw_full_primary_minus_sigma_weight);
531  for (size_t j = 0; j < daughter_dQdXs->size(); ++j) {
532  fFakeDataEvents.back().AddRecoDaughterTrackdQdX((*daughter_dQdXs)[j]);
533  fFakeDataEvents.back().AddRecoDaughterTrackResRange((*daughter_resRanges)[j]);
534  fFakeDataEvents.back().AddRecoDaughterEField((*daughter_EFields)[j]);
535  }
536 
537  for (size_t j = 0; j < g4rw_primary_grid_weights->size(); ++j) {
538  std::string name_full = "g4rw_full_grid_weights_" + std::to_string(j);
539  fFakeDataEvents.back().MakeG4RWBranch(name_full, (*g4rw_full_grid_weights)[j]);
540 
541  std::string name_primary = "g4rw_primary_grid_weights_" +
542  std::to_string(j);
543  fFakeDataEvents.back().MakeG4RWBranch(name_primary,
544  (*g4rw_primary_grid_weights)[j]);
545  }
546  fFakeDataEvents.back().MakeG4RWBranch("g4rw_full_grid_proton_weights",
547  (*g4rw_full_grid_proton_weights)[0]);
548  }
549 
550  std::cout << "Filled MC Events" << std::endl;
551 }
std::vector< ThinSliceEvent > fEvents
std::string string
Definition: nybbler.cc:12
std::vector< ThinSliceEvent > fFakeDataEvents
auto selection_ID
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
QTextStream & endl(QTextStream &s)
std::vector< double > protoana::PDSPThinSliceFitter::GetBestFitParsVec ( )
private

Definition at line 1478 of file PDSPThinSliceFitter.cxx.

1478  {
1479  std::vector<double> results;
1480 
1481  for (auto it = fSignalParameters.begin();
1482  it != fSignalParameters.end(); ++it) {
1483  for (size_t i = 0; i < it->second.size(); ++i) {
1484  results.push_back(it->second.at(i));
1485  }
1486  }
1487  for (auto it = fFluxParameters.begin();
1488  it != fFluxParameters.end(); ++it) {
1489  results.push_back(it->second);
1490  }
1491 
1492  for (auto it = fSystParameters.begin();
1493  it != fSystParameters.end(); ++it) {
1494  results.push_back(it->second.GetValue());
1495  }
1496 
1497  return results;
1498 }
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< int, double > fFluxParameters
std::map< int, std::vector< double > > fSignalParameters
void protoana::PDSPThinSliceFitter::GetCurrentTruthHists ( std::map< int, std::vector< TH1 * >> &  throw_hists,
std::map< int, std::vector< TH1 * >> &  throw_inc_hists,
std::map< int, std::vector< TH1 * >> &  throw_xsec_hists 
)
private

Definition at line 2591 of file PDSPThinSliceFitter.cxx.

2594  {
2595  //Loop over the samples
2596  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
2597  //Get the number of bins from the first entry of the beam energy bins
2598  std::vector<std::vector<ThinSliceSample>> & samples_vec_2D = it->second;
2599  size_t nBins = samples_vec_2D[0].size();
2600  std::string name = it->second[0][0].GetName() + "Throw" +
2601  std::to_string(throw_hists[it->first].size());
2602  TH1D * temp_hist = new TH1D(name.c_str(), "", nBins, 0, nBins);
2603  for (size_t i = 0; i < samples_vec_2D.size(); ++i) {
2604  std::vector<ThinSliceSample> & samples_vec = samples_vec_2D[i];
2605  for (size_t j = 0; j < it->second[i].size(); ++j) {
2606  temp_hist->AddBinContent(j+1, samples_vec[j].GetVariedFlux());
2607  }
2608  }
2609  throw_hists[it->first].push_back(temp_hist);
2610  }
2611 
2612  for (auto it = throw_inc_hists.begin(); it != throw_inc_hists.end(); ++it) {
2613  int s = it->first;
2614  auto & samples_vec_2D = fSamples[s];
2615  const std::vector<double> & bins = fSignalBins[s];
2616  std::string name = samples_vec_2D[0][0].GetName();
2617  name += "IncidentThrow" +
2618  std::to_string(throw_inc_hists[it->first].size());
2619  TH1D * temp_inc_hist = new TH1D(name.c_str(), "", bins.size() - 1, &bins[0]);
2620 
2621 
2622  name = samples_vec_2D[0][0].GetName();
2623  name += "XSecThrow" +
2624  std::to_string(throw_inc_hists[it->first].size());
2625  TH1D * temp_xsec_hist = new TH1D(name.c_str(), "", bins.size() - 1,
2626  &bins[0]);
2627  for (auto i_s : fIncidentSamples) {
2628  auto & incident_vec_2D = fSamples[i_s];
2629  for (size_t i = 0; i < incident_vec_2D.size(); ++i) {
2630  for (size_t j = 0; j < incident_vec_2D[i].size(); ++j) {
2631  if (fSliceMethod == "E") {
2632  incident_vec_2D[i][j].FillESliceHist(*temp_inc_hist);
2633  }
2634  else {
2635  incident_vec_2D[i][j].FillHistFromIncidentEnergies(*temp_inc_hist);
2636  }
2637  }
2638  }
2639  }
2640  throw_inc_hists[s].push_back(temp_inc_hist);
2641 
2642  for (int i = 1; i <= temp_xsec_hist->GetNbinsX(); ++i) {
2643  temp_xsec_hist->SetBinContent(
2644  i, throw_hists[s].back()->GetBinContent(i+1));
2645  }
2646  temp_xsec_hist->Divide(temp_inc_hist);
2647  throw_xsec_hists[s].push_back(temp_xsec_hist);
2648  }
2649 }
static QCString name
Definition: declinfo.cpp:673
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::string string
Definition: nybbler.cc:12
std::map< int, std::vector< double > > fSignalBins
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
static QCString * s
Definition: config.cpp:1042
void protoana::PDSPThinSliceFitter::GetNominalFluxes ( )
void protoana::PDSPThinSliceFitter::InitializeMCSamples ( )

Definition at line 194 of file PDSPThinSliceFitter.cxx.

194  {
195  for (size_t i = 0; i < fSampleSets.size(); ++i) {
196  fhicl::ParameterSet sample_set = fSampleSets[i];
197  std::string sample_name = sample_set.get<std::string>("Name");
198  int sample_ID = sample_set.get<int>("ID");
199 
200  if (fSamples.find(sample_ID) == fSamples.end()) {
201  fSamples[sample_ID] = std::vector<std::vector<ThinSliceSample>>();
202  fFakeSamples[sample_ID] = std::vector<std::vector<ThinSliceSample>>();
203  fFluxesBySample[sample_ID] = std::vector<std::vector<double>>();
204  fFakeFluxesBySample[sample_ID] = std::vector<std::vector<double>>();
205  }
206  fIsSignalSample[sample_ID] = sample_set.get<bool>("IsSignal");
207 
208  int flux_type = sample_set.get<int>("FluxType");
209  std::vector<double> bins
210  = sample_set.get<std::vector<double>>("SignalBins");
211  fSignalBins[sample_ID] = bins;
212 
213  for (size_t j = 1; j < fBeamEnergyBins.size(); ++j) {
214  fSamples[sample_ID].push_back(std::vector<ThinSliceSample>());
215  fFakeSamples[sample_ID].push_back(std::vector<ThinSliceSample>());
216  fFluxesBySample[sample_ID].push_back(std::vector<double>());
217  fFakeFluxesBySample[sample_ID].push_back(std::vector<double>());
218  if (sample_set.get<bool>("IsSignal")) {
219 
220  if (j == 1) {
221  fSignalParameters[sample_ID] = std::vector<double>();
222  fSignalParameterNames[sample_ID] = std::vector<std::string>();
223  }
224 
225  ThinSliceSample underflow_sample(
226  sample_name + "Underflow",
227  flux_type, fSelectionSets,
229  fSamples[sample_ID].back().push_back(underflow_sample);
230 
231  ThinSliceSample fake_underflow_sample(
232  sample_name + "FakeUnderflow",
233  flux_type, fSelectionSets,
235  fFakeSamples[sample_ID].back().push_back(fake_underflow_sample);
236 
237  if (j == 1 && fFitUnderOverflow) {
238  fSignalParameters[sample_ID].push_back(1.);
239 
240  std::string par_name = "par_" + sample_name + "_underflow";
241  fSignalParameterNames[sample_ID].push_back(par_name);
243  }
244 
245  fFluxesBySample[sample_ID].back().push_back(0.);
246  fFakeFluxesBySample[sample_ID].back().push_back(0.);
247 
248  for (size_t k = 1; k < bins.size(); ++k) {
249  ThinSliceSample sample(sample_name, flux_type, fSelectionSets,
251  j, true,
252  {bins[k-1], bins[k]});
253 
254  std::string fake_name = sample_name + "Fake";
255  ThinSliceSample fake_sample(fake_name, flux_type, fSelectionSets,
257  j, true,
258  {bins[k-1], bins[k]});
259  fSamples[sample_ID].back().push_back(sample);
260  fFakeSamples[sample_ID].back().push_back(fake_sample);
261 
262  if (j == 1) {
263  fSignalParameters[sample_ID].push_back(1.);
264 
265  std::string par_name = "par_" + sample_name + "_" +
266  PreciseToString(bins[k-1]) + "_" +
267  PreciseToString(bins[k]);
268  fSignalParameterNames[sample_ID].push_back(par_name);
270  }
271  fFluxesBySample[sample_ID].back().push_back(0.);
272  fFakeFluxesBySample[sample_ID].back().push_back(0.);
273  }
274 
275  ThinSliceSample overflow_sample(sample_name + "Overflow",
276  flux_type, fSelectionSets,
278  fSamples[sample_ID].back().push_back(overflow_sample);
279 
280  ThinSliceSample fake_overflow_sample(sample_name + "FakeOverflow",
281  flux_type, fSelectionSets,
283  if (j == 1 && fFitUnderOverflow) {
284  fSignalParameters[sample_ID].push_back(1.);
285 
286  std::string par_name = "par_" + sample_name + "_overflow";
287  fSignalParameterNames[sample_ID].push_back(par_name);
289  }
290 
291  fFakeSamples[sample_ID].back().push_back(fake_overflow_sample);
292  fFluxesBySample[sample_ID].back().push_back(0.);
293  fFakeFluxesBySample[sample_ID].back().push_back(0.);
294  }
295  else {
296  ThinSliceSample sample(sample_name, flux_type, fSelectionSets,
298  fSamples[sample_ID].back().push_back(sample);
299 
300  std::string fake_name = sample_name + "Fake";
301  ThinSliceSample fake_sample(fake_name, flux_type, fSelectionSets,
303  fFakeSamples[sample_ID].back().push_back(fake_sample);
304  fFluxesBySample[sample_ID].back().push_back(0.);
305  fFakeFluxesBySample[sample_ID].back().push_back(0.);
306  if (fFluxParameters.find(flux_type) != fFluxParameters.end() &&
307  j == 1) {
308  fFluxParsToSamples[flux_type].push_back(sample_ID);
309  }
310  }
311  }
312  }
313 
314  MakeMinimizer();
315 }
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::vector< double > fTrueIncidentBins
std::map< int, std::vector< int > > fFluxParsToSamples
std::string string
Definition: nybbler.cc:12
std::vector< double > fBeamEnergyBins
std::map< int, double > fFluxParameters
std::vector< double > fIncidentRecoBins
std::map< int, std::vector< std::vector< double > > > fFakeFluxesBySample
std::vector< fhicl::ParameterSet > fSampleSets
std::map< int, std::vector< std::string > > fSignalParameterNames
std::map< int, std::vector< double > > fSignalParameters
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::map< int, std::vector< std::vector< ThinSliceSample > > > fFakeSamples
std::map< int, std::vector< std::vector< double > > > fFluxesBySample
std::map< int, std::vector< double > > fSignalBins
std::vector< fhicl::ParameterSet > fSelectionSets
std::string PreciseToString(const double val, const int n=2)
std::map< int, bool > fIsSignalSample
void protoana::PDSPThinSliceFitter::MakeMinimizer ( )
private

Definition at line 64 of file PDSPThinSliceFitter.cxx.

64  {
65  fMinimizer = std::unique_ptr<ROOT::Math::Minimizer>
66  (ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"));
67 
68  fMinimizer->SetMaxFunctionCalls(fMaxCalls);
69  fMinimizer->SetMaxIterations(fMaxIterations);
70  fMinimizer->SetTolerance(fTolerance);
71 
73  TH1D parsHist("preFitPars", "", total_parameters, 0, total_parameters);
74  fPreFitParsNormal = TH1D("preFitParsNormal", "",
75  total_parameters, 0, total_parameters);
76 
77  size_t n_par = 0;
78  for (auto it = fSignalParameters.begin();
79  it != fSignalParameters.end(); ++it) {
80  for (size_t i = 0; i < it->second.size(); ++i) {
81  if (fRandomStart) {
82  it->second[i] = fRNG.Gaus(1., .1);
83  }
84  fMinimizer->SetVariable(n_par, fSignalParameterNames[it->first][i].c_str(),
85  it->second[i], 0.01);
86  fMinimizer->SetVariableLimits(n_par, fLowerLimit, fUpperLimit);
87  fMinimizerInitVals.push_back(it->second[i]);
88  parsHist.SetBinContent(n_par+1, it->second[i]);
89  parsHist.SetBinError(n_par+1, 0.);
90  //parsHist.GetXaxis()->SetBinLabel(
91  // n_par+1, fSignalParameterNames[it->first][i].c_str());
92 
93  fPreFitParsNormal.SetBinContent(n_par+1, it->second[i]);
94  fPreFitParsNormal.SetBinError(n_par+1, 0.);
95  fPreFitParsNormal.GetXaxis()->SetBinLabel(
96  n_par+1, fSignalParameterNames[it->first][i].c_str());
97  fParLimits.push_back(0.);
98  fParLimitsUp.push_back(100.);
99  ++n_par;
100 
101  }
102  }
103 
104  for (auto it = fFluxParameters.begin();
105  it != fFluxParameters.end(); ++it) {
106  if (fRandomStart) {
107  it->second = fRNG.Gaus(1., .1);
108  }
109  fMinimizer->SetVariable(n_par, fFluxParameterNames[it->first].c_str(),
110  it->second, 0.01);
111  fMinimizerInitVals.push_back(it->second);
112  fMinimizer->SetVariableLimits(n_par, fLowerLimit, fUpperLimit);
113  parsHist.SetBinContent(n_par+1, it->second);
114  parsHist.SetBinError(n_par+1, 0.);
115  //parsHist.GetXaxis()->SetBinLabel(
116  // n_par+1, fFluxParameterNames[it->first].c_str());
117  fPreFitParsNormal.SetBinContent(n_par+1, it->second);
118  fPreFitParsNormal.SetBinError(n_par+1, 0.);
119  fPreFitParsNormal.GetXaxis()->SetBinLabel(
120  n_par+1, fFluxParameterNames[it->first].c_str());
121  fParLimits.push_back(0.);
122  fParLimitsUp.push_back(100.);
123  ++n_par;
124  }
125 
126  for (auto it = fSystParameters.begin(); it != fSystParameters.end(); ++it) {
127  std::cout << "Adding parameter " << it->second.GetName().c_str() <<
128  " " << it->second.GetCentral() << std::endl;
129  double val = it->second.GetCentral();
130  if (fRandomStart) {
131  it->second.SetValue(fRNG.Gaus(1., .1)*it->second.GetCentral());
132  val = it->second.GetValue();
133  }
134  fMinimizer->SetVariable(n_par, it->second.GetName().c_str(),
135  val, 0.01);
136  fMinimizerInitVals.push_back(val);
137  fMinimizer->SetVariableLimits(n_par, it->second.GetLowerLimit(),
138  it->second.GetUpperLimit());
139  if (fFixVariables &&
140  fSystsToFix.find(it->first) != fSystsToFix.end()) {
141  std::cout << "Fixing variable " << it->second.GetName() <<
142  " to value " << fSystsToFix.at(it->first) <<
143  std::endl;
144  fMinimizer->SetVariableValue(n_par, fSystsToFix.at(it->first));
145  fMinimizer->FixVariable(n_par);
146  }
147 
148 
149  fParLimits.push_back(it->second.GetThrowLimit());
150  fParLimitsUp.push_back(it->second.GetThrowLimitUp());
151  if (abs(it->second.GetCentral()) < 1.e-5) {
152  parsHist.SetBinContent(n_par+1, val + 1.);
153  if (fAddSystTerm) {
154  size_t cov_bin = fCovarianceBins[it->first];
155  parsHist.SetBinError(n_par+1, sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin]));
156  std::cout << "1 Set Error: " << it->first << " " <<
157  sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin]) << std::endl;
158  }
159  }
160  else {
161  parsHist.SetBinContent(n_par+1,
162  val/it->second.GetCentral());
163  if (fAddSystTerm) {
164  size_t cov_bin = fCovarianceBins[it->first];
165  parsHist.SetBinError(n_par+1,
166  sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin])/it->second.GetCentral());
167  std::cout << "2 Set Error: " << it->first << " " <<
168  sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin]) << " " << it->second.GetCentral() << " " <<
169  sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin])/it->second.GetCentral() << std::endl;
170  }
171  }
172 
173  //parsHist.GetXaxis()->SetBinLabel(
174  // n_par+1, it->second.GetName().c_str());
175 
176  fPreFitParsNormal.SetBinContent(n_par+1, it->second.GetValue());
177  fPreFitParsNormal.GetXaxis()->SetBinLabel(
178  n_par+1, it->second.GetName().c_str());
179  if (fAddSystTerm) {
180  size_t cov_bin = fCovarianceBins[it->first];
181  fPreFitParsNormal.SetBinError(n_par+1, sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin]));
182  }
183  ++n_par;
184  }
185 
186  parsHist.SetMarkerColor(kBlue);
187  parsHist.SetMarkerStyle(20);
188  fOutputFile.cd();
189  parsHist.Write();
190  fPreFitParsNormal.Write();
191 
192 }
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< int, double > fFluxParameters
std::vector< double > fParLimitsUp
std::map< int, std::vector< std::string > > fSignalParameterNames
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
std::map< std::string, size_t > fCovarianceBins
T abs(T value)
const double e
std::map< int, std::vector< double > > fSignalParameters
std::vector< double > fMinimizerInitVals
std::map< int, std::string > fFluxParameterNames
std::map< std::string, double > fSystsToFix
std::vector< double > fParLimits
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::MakeThrowsTree ( TTree &  tree,
std::vector< double > &  branches 
)
private

Definition at line 2779 of file PDSPThinSliceFitter.cxx.

2779  {
2780  for (auto it = fSignalParameters.begin(); it != fSignalParameters.end(); ++it) {
2781  for (size_t i = 0; i < it->second.size(); ++i) {
2782  branches.push_back(0.);
2783  tree.Branch(fSignalParameterNames[it->first][i].c_str(),
2784  &branches.back());
2785  }
2786  }
2787 
2788  for (auto it = fFluxParameters.begin(); it != fFluxParameters.end(); ++it) {
2789  branches.push_back(0.);
2790  tree.Branch(fFluxParameterNames[it->first].c_str(),
2791  &branches.back());
2792  }
2793 
2794  for (auto it = fSystParameters.begin(); it != fSystParameters.end(); ++it) {
2795  branches.push_back(0.);
2796  tree.Branch(it->second.GetName().c_str(), &branches.back());
2797  }
2798 }
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< int, double > fFluxParameters
std::map< int, std::vector< std::string > > fSignalParameterNames
std::map< int, std::vector< double > > fSignalParameters
std::map< int, std::string > fFluxParameterNames
void protoana::PDSPThinSliceFitter::NormalFit ( )
private

Definition at line 936 of file PDSPThinSliceFitter.cxx.

936  {
937  std::cout << "Running Fit" << std::endl;
938  bool found_minimum = false;
939  try {
940  found_minimum = fMinimizer->Minimize();
941  }
942  catch (const std::exception & e) {
943  std::cerr << e.what() << std::endl;
944  std::cerr << "exiting safely" << std::endl;
945  found_minimum = 0;
946  }
947 
948  TVector output_fit_status(1);
949  output_fit_status[0] = fMinimizer->Status();
950 
951  if (!found_minimum) {
952  std::cout << "Failed to find minimum" << std::endl;
953  fOutputFile.cd();
954  output_fit_status.Write("fit_status");
955 
956  if (fMinimizer->Status() == 3) {
957  std::cout << "Edm too large: " << fMinimizer->Edm() << std::endl;
958  }
959 
960  if (fDoScans) ParameterScans();
961  }
962  else {
963  std::vector<double> vals;
964  std::cout << "Found minimimum. Fit status: " << fMinimizer->Status() << std::endl;
965  double chi2_syst = (fAddSystTerm ? CalcChi2SystTerm() : 0.);
966  double chi2_stat = fThinSliceDriver->CalculateChi2(fSamples, fDataSet).first;
967  //double chi2_reg = (fAddRegTerm ? CalcRegTerm() : 0.);
968  std::cout << "chi2 Syst: " << chi2_syst << std::endl;
969  std::cout << "chi2 Stat: " << chi2_stat << std::endl;
970 
971 
972  if (fRunHesse) {
973  std::cout << "Running Hesse" << std::endl;
974  bool hesse_good = fMinimizer->Hesse();
975  std::cout << "hesse good? " << hesse_good << std::endl;
976  TVector output_hesse_status(1);
977 
978  fOutputFile.cd();
979  output_hesse_status[0] = (hesse_good ? 1 : 0);
980  output_hesse_status.Write("hesse_status");
981  }
982  /*
983  std::cout <<
984  fTotalSignalParameters << " " <<
985  fTotalFluxParameters << " " <<
986  fTotalSystParameters << std::endl;*/
988  std::cout << total_parameters << std::endl;
989  for (size_t i = 0; i < total_parameters; ++i) {
990  std::cout << fMinimizer->VariableName(i) << " " << fMinimizer->X()[i] <<
991  std::endl;
992  vals.push_back(fMinimizer->X()[i]);
993  }
994 
995  TH2D covHist("covHist", "", total_parameters, 0,
996  total_parameters, total_parameters, 0,
997  total_parameters);
998  TH2D corrHist("corrHist", "", total_parameters, 0,
999  total_parameters, total_parameters, 0,
1000  total_parameters);
1001 
1002  TH1D parsHist("postFitPars", "", total_parameters, 0,
1003  total_parameters);
1004  TH1D parsHist_normal_val("postFitParsNormal", "",
1005  total_parameters, 0, total_parameters);
1006  TH1D * toyParsHist = 0x0,
1007  * toyParsHist_normal = 0x0;
1008  if (fFakeDataRoutine == "Toy") {
1009  toyParsHist = (TH1D*)parsHist.Clone("toyPars");
1010  toyParsHist->Reset();
1011  toyParsHist_normal = (TH1D*)parsHist.Clone("toyPars_normal");
1012  toyParsHist_normal->Reset();
1013  }
1014 
1015 
1016 
1017  int n_par = 0;
1018  for (auto it = fSignalParameters.begin();
1019  it != fSignalParameters.end(); ++it) {
1020  for (size_t i = 0; i < it->second.size(); ++i) {
1021  parsHist.SetBinContent(n_par+1, fMinimizer->X()[n_par]);
1022  parsHist.SetBinError(n_par+1,
1023  sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1024  parsHist.GetXaxis()->SetBinLabel(
1025  n_par+1, fSignalParameterNames[it->first][i].c_str());
1026 
1027  parsHist_normal_val.SetBinContent(n_par+1, fMinimizer->X()[n_par]);
1028  parsHist_normal_val.SetBinError(
1029  n_par+1, sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1030  parsHist_normal_val.GetXaxis()->SetBinLabel(
1031  n_par+1, fSignalParameterNames[it->first][i].c_str());
1032 
1033  if (fFakeDataRoutine == "Toy")
1034  toyParsHist->SetBinContent(n_par+1, -999.);
1035 
1036 
1037  //covHist.GetXaxis()->SetBinLabel(
1038  // n_par+1, fSignalParameterNames[it->first][i].c_str());
1039  //corrHist.GetXaxis()->SetBinLabel(
1040  // n_par+1, fSignalParameterNames[it->first][i].c_str());
1041  //covHist.GetYaxis()->SetBinLabel(
1042  // n_par+1, fSignalParameterNames[it->first][i].c_str());
1043  //corrHist.GetYaxis()->SetBinLabel(
1044  // n_par+1, fSignalParameterNames[it->first][i].c_str());
1045 
1046  ++n_par;
1047  }
1048  }
1049  for (auto it = fFluxParameters.begin();
1050  it != fFluxParameters.end(); ++it) {
1051  parsHist.SetBinContent(n_par+1, fMinimizer->X()[n_par]);
1052  parsHist.SetBinError(n_par+1,
1053  sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1054  parsHist.GetXaxis()->SetBinLabel(
1055  n_par+1, fFluxParameterNames[it->first].c_str());
1056 
1057  parsHist_normal_val.SetBinContent(n_par+1, fMinimizer->X()[n_par]);
1058  parsHist_normal_val.SetBinError(
1059  n_par+1, sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1060  parsHist_normal_val.GetXaxis()->SetBinLabel(
1061  n_par+1, fFluxParameterNames[it->first].c_str());
1062 
1063 
1064  if (fFakeDataRoutine == "Toy")
1065  toyParsHist->SetBinContent(n_par+1, -999.);
1066 
1067  //covHist.GetXaxis()->SetBinLabel(
1068  // n_par+1, fFluxParameterNames[it->first].c_str());
1069  //corrHist.GetXaxis()->SetBinLabel(
1070  // n_par+1, fFluxParameterNames[it->first].c_str());
1071  //covHist.GetYaxis()->SetBinLabel(
1072  // n_par+1, fFluxParameterNames[it->first].c_str());
1073  //corrHist.GetYaxis()->SetBinLabel(
1074  // n_par+1, fFluxParameterNames[it->first].c_str());
1075 
1076  ++n_par;
1077  }
1078 
1079  for (auto it = fSystParameters.begin();
1080  it != fSystParameters.end(); ++it) {
1081 
1082  if (abs(it->second.GetCentral()) < 1.e-5) {
1083  parsHist.SetBinContent(n_par+1, fMinimizer->X()[n_par] + 1.);
1084  parsHist.SetBinError(n_par+1,
1085  sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1086  if (fFakeDataRoutine == "Toy")
1087  toyParsHist->SetBinContent(n_par+1, fToyValues[it->first] + 1.);
1088 
1089  }
1090  else {
1091  parsHist.SetBinContent(n_par+1, fMinimizer->X()[n_par]/it->second.GetCentral());
1092  parsHist.SetBinError(n_par+1,
1093  sqrt(fMinimizer->CovMatrix(n_par, n_par))/it->second.GetCentral());
1094  if (fFakeDataRoutine == "Toy")
1095  toyParsHist->SetBinContent(n_par+1, fToyValues[it->first]/it->second.GetCentral());
1096  }
1097 
1098  if (fFakeDataRoutine == "Toy") {
1099  toyParsHist_normal->SetBinContent(n_par+1, fToyValues[it->first]);
1100  size_t cov_bin = fCovarianceBins[it->first];
1101  toyParsHist_normal->SetBinError(
1102  n_par+1, sqrt((*fCovMatrixDisplay)[cov_bin][cov_bin]));
1103  }
1104 
1105  parsHist.GetXaxis()->SetBinLabel(
1106  n_par+1, it->second.GetName().c_str());
1107 
1108  parsHist_normal_val.SetBinContent(n_par+1, fMinimizer->X()[n_par]);
1109  parsHist_normal_val.SetBinError(
1110  n_par+1, sqrt(fMinimizer->CovMatrix(n_par, n_par)));
1111  parsHist_normal_val.GetXaxis()->SetBinLabel(
1112  n_par+1, it->second.GetName().c_str());
1113 
1114 
1115  //covHist.GetXaxis()->SetBinLabel(
1116  // n_par+1, it->second.GetName().c_str());
1117  //corrHist.GetXaxis()->SetBinLabel(
1118  // n_par+1, it->second.GetName().c_str());
1119  //covHist.GetYaxis()->SetBinLabel(
1120  // n_par+1, it->second.GetName().c_str());
1121  //corrHist.GetYaxis()->SetBinLabel(
1122  // n_par+1, it->second.GetName().c_str());
1123 
1124  ++n_par;
1125  }
1126 
1127  TMatrixD * cov = new TMatrixD(total_parameters, total_parameters);
1128 
1129  //std::cout << "First cov: " << std::endl;
1130  for (size_t i = 0; i < total_parameters; ++i) {
1131  for (size_t j = 0; j < total_parameters; ++j) {
1132  double cov_val = fMinimizer->CovMatrix(i, j);
1133  covHist.SetBinContent(i+1, j+1, cov_val);
1134  corrHist.SetBinContent(i+1, j+1, fMinimizer->Correlation(i, j));
1135  //std::cout << fMinimizer->Correlation(i, j) << " ";
1136  (*cov)(i, j) = fMinimizer->CovMatrix(i, j);
1137  }
1138  //std::cout << std::endl;
1139  }
1140 
1141  if (fFixVariables && fAddSystTerm) {
1143  for (auto it = fSystParameters.begin(); it != fSystParameters.end();
1144  ++it) {
1145  if (fSystsToFix.find(it->first) != fSystsToFix.end()) {
1146  size_t cov_bin = fCovarianceBins[it->first];
1147  covHist.SetBinContent(n_par+1, n_par+1,
1148  (*fCovMatrixDisplay)[cov_bin][cov_bin]);
1149  (*cov)(n_par, n_par) = (*fCovMatrixDisplay)[cov_bin][cov_bin];
1150  }
1151  ++n_par;
1152  }
1153 
1154  /*
1155  std::cout << "New cov: " << std::endl;
1156  for (size_t i = 0; i < total_parameters; ++i) {
1157  for (size_t j = 0; j < total_parameters; ++j) {
1158 
1159  std::cout << covHist.GetBinContent(i+1, j+1) << " ";
1160  }
1161  std::cout << std::endl;
1162  }*/
1163  }
1164 
1165  fOutputFile.cd();
1166  output_fit_status.Write("fit_status");
1167  parsHist.SetFillColor(kRed);
1168  parsHist.SetFillStyle(3144);
1169  parsHist.SetMarkerColor(kRed);
1170  parsHist.SetMarkerStyle(20);
1171  parsHist.Write();
1172  parsHist_normal_val.Write();
1173 
1174  TVectorD chi2_syst_out(1);
1175  chi2_syst_out[0] = chi2_syst;
1176  chi2_syst_out.Write("chi2_syst");
1177  TVectorD chi2_stat_out(1);
1178  chi2_stat_out[0] = chi2_stat;
1179  chi2_stat_out.Write("chi2_stat");
1180 
1184 
1185  //Drawing pre + post fit pars
1186  TCanvas cPars("cParameters", "");
1187  cPars.SetTicks();
1188  parsHist.GetYaxis()->SetTitle("Parameter Values");
1189  parsHist.SetTitleSize(.05, "Y");
1190  //parsHist.Draw("e2");
1191  TH1D * preFitParsHist = (TH1D*)fOutputFile.Get("preFitPars");
1192  preFitParsHist->SetTitle(";Parameter;Parameter Value");
1193  preFitParsHist->GetXaxis()->SetTitleOffset(.8);
1194  preFitParsHist->GetYaxis()->SetTitleOffset(.8);
1195  preFitParsHist->SetTitleSize(.05, "XY");
1196  preFitParsHist->SetFillColor(kBlue);
1197  preFitParsHist->SetFillStyle(3001);
1198  //preFitParsHist->Draw("pe2 same");
1199  preFitParsHist->SetMaximum(3.);
1200  preFitParsHist->SetMinimum(-1.);
1201  preFitParsHist->Draw("pe2");
1202  parsHist.Draw("e2 same");
1203  TLine l(0., 1., parsHist.GetXaxis()->GetBinUpEdge(parsHist.GetNbinsX()), 1.);
1204  l.SetLineColor(kBlack);
1205  l.Draw("same");
1206  TLegend leg(.15, .65, .45, .85);
1207  leg.AddEntry(preFitParsHist, "Pre-Fit #pm 1 #sigma", "pf");
1208  leg.AddEntry(&parsHist, "Post-Fit #pm 1 #sigma", "pf");
1209 
1210  if (fFakeDataRoutine == "Toy") {
1211  toyParsHist->SetMarkerStyle(20);
1212  toyParsHist->SetMarkerColor(kBlack);
1213  toyParsHist->Draw("same p");
1214  leg.AddEntry(toyParsHist, "Toy Value", "p");
1215  toyParsHist->Write();
1216  toyParsHist_normal->Write();
1217  }
1218  preFitParsHist->SetLineWidth(0);
1219  preFitParsHist->Draw("p same");
1220  parsHist.SetLineWidth(0);
1221  parsHist.Draw("p same");
1222 
1223  /*
1224  if (fAddSystTerm) {
1225  //std::string chi2_str = "#chi^{2}_{Syst}";
1226  //chi2_str += std::to_string(chi2_syst);
1227  //leg.AddEntry((TObject*)0x0, chi2_str.c_str(), "");
1228  TString chi2_str;
1229  //chi2_str.Form("#chi^{2}_{Syst} = %.2f", chi2_syst);
1230  chi2_str.Form("#chi^{2} = %.2f", chi2_syst);
1231  leg.AddEntry((TObject*)0x0, chi2_str, "");
1232  }*/
1233  leg.Draw();
1234  cPars.Write();
1235 
1236  covHist.SetTitle(";Parameter;Parameter");
1237  covHist.SetTitleSize(.05, "XY");
1238  covHist.GetXaxis()->SetTitleOffset(.8);
1239  covHist.GetYaxis()->SetTitleOffset(.8);
1240  covHist.Write();
1241  cov->Write("CovMatrix");
1242  corrHist.SetTitle(";Parameter;Parameter");
1243  corrHist.SetTitleSize(.05, "XY");
1244  corrHist.GetXaxis()->SetTitleOffset(.8);
1245  corrHist.GetYaxis()->SetTitleOffset(.8);
1246  corrHist.SetMaximum(1.);
1247  corrHist.SetMinimum(-1.);
1248  corrHist.Write();
1249 
1250  TCanvas cCorr("cCorr", "");
1251  corrHist.Draw("colz");
1252  cCorr.Write();
1253 
1254  if (fDoScans)
1255  ParameterScans();
1256 
1257  fFillIncidentInFunction = true;
1258  fFitFunction(&vals[0]);
1259 
1260  //SetBestFit(); //Deprecated
1261 
1262  //save post fit stacks
1263  TDirectory * top_dir = fOutputFile.GetDirectory("");
1264  TDirectory * xsec_dir = fOutputFile.mkdir("PostFitXSec");
1265  std::string extra_name = "PostFit";
1266  CompareDataMC(extra_name, xsec_dir, top_dir, true);
1267  if (fDoThrows)
1268  DoThrows(parsHist_normal_val, cov);
1269 
1270  if (fDo1DShifts) {
1271  if (fAddSystTerm)
1273  Do1DShifts(parsHist_normal_val);
1274  }
1275  }
1276 }
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::string string
Definition: nybbler.cc:12
std::map< std::string, ThinSliceSystematic > fSystParameters
std::map< int, double > fFluxParameters
std::map< std::string, double > fToyValues
std::map< int, double > fBestFitFluxPars
static QStrList * l
Definition: config.cpp:1044
virtual std::pair< double, size_t > CalculateChi2(std::map< int, std::vector< std::vector< ThinSliceSample >>> &samples, ThinSliceDataSet &data_set)=0
std::map< int, std::vector< double > > fBestFitSignalPars
std::map< std::string, ThinSliceSystematic > fBestFitSystPars
std::map< int, std::vector< std::string > > fSignalParameterNames
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
std::map< std::string, size_t > fCovarianceBins
T abs(T value)
const double e
std::map< int, std::vector< double > > fSignalParameters
std::map< int, std::string > fFluxParameterNames
void CompareDataMC(std::string extra_name, TDirectory *xsec_dir, TDirectory *plot_dir, bool post_fit=false)
std::map< std::string, double > fSystsToFix
void Do1DShifts(const TH1D &pars, bool prefit=false)
void DoThrows(const TH1D &pars, const TMatrixD *cov)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::ParameterScans ( )
private

Definition at line 1500 of file PDSPThinSliceFitter.cxx.

1500  {
1501  fOutputFile.cd();
1502  TDirectory * out = (TDirectory *)fOutputFile.mkdir("Scans");
1503  out->cd();
1504 
1506  std::cout << "Total parameters " << total_parameters << std::endl;
1507 
1508  double * x = new double[fNScanSteps] {};
1509  double * y = new double[fNScanSteps] {};
1510  for (size_t i = 0; i < total_parameters; ++i) {
1511  std::cout << "\tParameter " << fMinimizer->VariableName(i) << std::endl;
1512  bool scanned = fMinimizer->Scan(i, fNScanSteps, x, y);
1513  if (scanned) {
1514  TGraph gr(fNScanSteps - 1, x, y);
1515  gr.Write(fMinimizer->VariableName(i).c_str());
1516  }
1517  }
1518 
1519  delete[] x;
1520  delete[] y;
1521 }
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
list x
Definition: train.py:276
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::PlotThrows ( std::map< int, std::vector< TH1 * >> &  throw_hists,
std::map< int, std::vector< TH1 * >> &  truth_throw_hists,
std::map< int, std::vector< TH1 * >> &  truth_inc_hists,
std::map< int, std::vector< TH1 * >> &  truth_xsec_hists 
)
private

Definition at line 2045 of file PDSPThinSliceFitter.cxx.

2049  {
2050  std::map<int, TH1*> data_hists
2051  = (fPlotRebinned ?
2054 
2055  //Build best fit hists and get bins for covariance
2056  //std::map<int, TH1D*> best_fit_selection_hists;
2057  int nBins = 0;
2058  for (auto it = data_hists.begin(); it != data_hists.end(); ++it ) {
2059  //TH1D * best_fit_hist = (TH1D*)it->second->Clone();
2060  //best_fit_hist->Reset();
2061  //for (auto it2 = fSamples.begin(); it2 != fSamples.end(); ++it2) {
2062  // for (size_t i = 0; i < it2->second.size(); ++i) {
2063  // for (size_t j = 0; j < it2->second[i].size(); ++j) {
2064  // it2->second[i][j].SetFactorToBestFit();
2065  // best_fit_hist->Add(
2066  // (TH1D*)(fPlotRebinned ?
2067  // it2->second[i][j].GetRebinnedSelectionHist(it->first) :
2068  // it2->second[i][j].GetSelectionHist(it->first)));
2069  // }
2070  // }
2071  //}
2072  //best_fit_selection_hists[it->first] = best_fit_hist;
2073  //nBins += best_fit_hist->GetNbinsX();
2074  nBins += it->second->GetNbinsX();
2075  }
2076 
2077  TH2D selection_cov("SelectionCov", "", nBins, 0, nBins, nBins, 0, nBins);
2078 
2079  nBins = 0;
2080  std::map<int, size_t> sample_bins;
2081  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
2082  nBins += it->second[0].size();
2083  sample_bins[it->first] = it->second[0].size();
2084  }
2085 
2086  //std::map<int, std::vector<double>> best_fit_truth;
2087  std::map<int, std::vector<double>> best_fit_errs;
2088 
2089  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
2090  //best_fit_truth[it->first]
2091  // = std::vector<double>(sample_bins[it->first], 0.);
2092  best_fit_errs[it->first]
2093  = std::vector<double>(sample_bins[it->first], 0.);
2094 
2095  //for (size_t i = 0; i < sample_bins[it->first]; ++i) {
2096  // double best_fit_val_i = 0.;
2097  // for (size_t j = 0; j < it->second.size(); ++j) {
2098  // best_fit_val_i += it->second[j][i].GetVariedFlux();
2099  // }
2100 
2101  // best_fit_truth[it->first][i] = best_fit_val_i;
2102  //}
2103  }
2104 
2105  TH2D interaction_cov("interaction_cov", "", nBins, 0, nBins, nBins, 0, nBins);
2106  std::map<int, std::vector<double>> best_fit_inc_truth;
2107  std::map<int, std::vector<double>> best_fit_xsec_truth;
2108  std::map<int, std::vector<double>> best_fit_inc_errs;
2109  std::map<int, std::vector<double>> best_fit_xsec_errs;
2110 
2111  nBins = 0;
2112  std::map<int, size_t> xsec_bins;
2113  for (auto it = fBestFitIncs.begin(); it != fBestFitIncs.end(); ++it) {
2114  int s = it->first;
2115  nBins += it->second->GetNbinsX();
2116  xsec_bins[s] = it->second->GetNbinsX();
2117 
2118  best_fit_inc_truth[s] = std::vector<double>(xsec_bins[s], 0.);
2119  best_fit_xsec_truth[s] = std::vector<double>(xsec_bins[s], 0.);
2120  best_fit_inc_errs[s] = std::vector<double>(xsec_bins[s], 0.);
2121  best_fit_xsec_errs[s] = std::vector<double>(xsec_bins[s], 0.);
2122 
2123  for (size_t i = 0; i < xsec_bins[s]; ++i) {
2124  best_fit_inc_truth[s][i] = it->second->GetBinContent(i+1);
2125  best_fit_xsec_truth[s][i] = fBestFitXSecs[s]->GetBinContent(i+1);
2126  }
2127  }
2128 
2129  //TH2D incident_cov("incident_cov", "", nBins, 0, nBins, nBins, 0, nBins);
2130  TH2D xsec_cov("xsec_cov", "", nBins, 0, nBins, nBins, 0, nBins);
2131  TH2D xsec_corr("xsec_corr", "", nBins, 0, nBins, nBins, 0, nBins);
2132  xsec_cov.SetTitle(";Cross Section Bin;Cross section Bin");
2133  xsec_cov.SetTitleSize(.05, "XY");
2134  xsec_cov.GetXaxis()->SetTitleOffset(.8);
2135  xsec_corr.SetTitle(";Cross Section Bin;Cross section Bin");
2136  xsec_corr.SetTitleSize(.05, "XY");
2137  xsec_corr.GetXaxis()->SetTitleOffset(.8);
2138  xsec_corr.SetMaximum(1.);
2139  xsec_corr.SetMinimum(-1.);
2140  TMatrixD xsec_cov_matrix(nBins, nBins);
2141 
2142  std::map<int, std::vector<double>> mean_xsecs;
2143  for (size_t z = 0; z < fNThrows; ++z) {
2144  for (auto it = truth_xsec_hists.begin(); it != truth_xsec_hists.end(); ++it) {
2145  std::vector<TH1 *> xsec_hists_i = it->second;
2146  if (z == 0)
2147  mean_xsecs[it->first] = std::vector<double>(xsec_bins[it->first], 0.);
2148  for (size_t i = 0; i < xsec_bins[it->first]; ++i) {
2149  mean_xsecs[it->first][i] += xsec_hists_i[z]->GetBinContent(i+1)/fNThrows;
2150  }
2151  }
2152  }
2153 
2154  for (size_t z = 0; z < fNThrows; ++z) {
2155  int bin_i = 1;
2156  for (auto it = fBestFitSelectionHists.begin();
2157  it != fBestFitSelectionHists.end(); ++it) {
2158  TH1D * best_fit = it->second;
2159  int selection_ID = it->first;
2160  std::vector<TH1*> & temp_throws = throw_hists[selection_ID];
2161  for (int i = 1; i <= best_fit->GetNbinsX(); ++i) {
2162  double best_fit_val_i = best_fit->GetBinContent(i);
2163  int bin_j = 1;
2164  for (auto it2 = fBestFitSelectionHists.begin();
2165  it2 != fBestFitSelectionHists.end(); ++it2) {
2166 
2167  TH1D * best_fit_2 = it2->second;
2168  int selection_ID_2 = it2->first;
2169  std::vector<TH1*> & temp_throws_2 = throw_hists[selection_ID_2];
2170  for (int j = 1; j <= best_fit_2->GetNbinsX(); ++j) {
2171  double best_fit_val_j = best_fit_2->GetBinContent(j);
2172  double val = (best_fit_val_i - temp_throws[z]->GetBinContent(i))*
2173  (best_fit_val_j - temp_throws_2[z]->GetBinContent(j));
2174  selection_cov.SetBinContent(
2175  bin_i, bin_j, (val/temp_throws.size() +
2176  selection_cov.GetBinContent(bin_i, bin_j)));
2177  ++bin_j;
2178  }
2179  }
2180  ++bin_i;
2181  }
2182  }
2183 
2184  bin_i = 1;
2185  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
2186  std::vector<TH1 *> throw_hists_i = truth_throw_hists[it->first];
2187 
2188  for (size_t i = 0; i < sample_bins[it->first]; ++i) {
2189  double best_fit_val_i = fBestFitTruthVals[it->first][i];
2190 
2191  int bin_j = 1;
2192  for (auto it2 = fSamples.begin(); it2 != fSamples.end(); ++it2) {
2193  std::vector<TH1 *> throw_hists_j = truth_throw_hists[it2->first];
2194  for (size_t j = 0; j < sample_bins[it2->first]; ++j) {
2195  double best_fit_val_j = fBestFitTruthVals[it2->first][j];
2196 
2197  double val
2198  = (throw_hists_i[z]->GetBinContent(i+1) - best_fit_val_i)*
2199  (throw_hists_j[z]->GetBinContent(j+1) - best_fit_val_j);
2200  interaction_cov.SetBinContent(
2201  bin_i, bin_j,
2202  (interaction_cov.GetBinContent(bin_i, bin_j) +
2203  val/throw_hists_i.size()));
2204  if (bin_i == bin_j && (z == fNThrows - 1)) {
2205  best_fit_errs[it->first][i]
2206  = sqrt(interaction_cov.GetBinContent(bin_i, bin_j));
2207  }
2208  ++bin_j;
2209  }
2210  }
2211 
2212  ++bin_i;
2213  }
2214  }
2215 
2216 
2217 
2218 
2219  bin_i = 1;
2220  for (auto it = truth_inc_hists.begin(); it != truth_inc_hists.end(); ++it) {
2221  std::vector<TH1 *> xsec_hists_i = truth_xsec_hists[it->first];
2222  for (size_t i = 0; i < xsec_bins[it->first]; ++i) {
2223  //double best_fit_xsec_i = mean_xsecs[it->first][i];
2224  double best_fit_xsec_i = best_fit_xsec_truth[it->first][i];
2225 
2226  int bin_j = 1;
2227  for (auto it2 = truth_inc_hists.begin(); it2 != truth_inc_hists.end();
2228  ++it2) {
2229  std::vector<TH1 *> xsec_hists_j = truth_xsec_hists[it2->first];
2230  for (size_t j = 0; j < xsec_bins[it2->first]; ++j) {
2231  //double best_fit_xsec_j = mean_xsecs[it2->first][j];
2232  double best_fit_xsec_j = best_fit_xsec_truth[it2->first][j];
2233 
2234  double val
2235  = (xsec_hists_i[z]->GetBinContent(i+1) - best_fit_xsec_i)*
2236  (xsec_hists_j[z]->GetBinContent(j+1) - best_fit_xsec_j);
2237  xsec_cov.SetBinContent(
2238  bin_i, bin_j,
2239  (xsec_cov.GetBinContent(bin_i, bin_j) +
2240  val/fNThrows));
2241  if (bin_i == bin_j && (z == fNThrows - 1)) {
2242  best_fit_xsec_errs[it->first][i]
2243  = sqrt(xsec_cov.GetBinContent(bin_i, bin_j));
2244  }
2245  ++bin_j;
2246  }
2247  }
2248  ++bin_i;
2249  }
2250  }
2251  }
2252 
2253  for (int i = 0; i < nBins; ++i) {
2254  for (int j = 0; j < nBins; ++j) {
2255  xsec_cov_matrix[i][j] = xsec_cov.GetBinContent(i+1, j+1);
2256  double corr_val =
2257  xsec_cov.GetBinContent(i+1, j+1)/
2258  sqrt(xsec_cov.GetBinContent(i+1, i+1)*
2259  xsec_cov.GetBinContent(j+1, j+1));
2260  xsec_corr.SetBinContent(i+1, j+1, corr_val);
2261  }
2262  }
2263  xsec_cov_matrix.Invert();
2264 
2265 
2266  fOutputFile.cd("Throws");
2267  selection_cov.Write();
2268  interaction_cov.Write();
2269  xsec_cov.Write();
2270  xsec_corr.Write();
2271 
2272  TCanvas cXSecCorr("cXSecCorr", "");
2273  xsec_corr.Draw("colz");
2274  cXSecCorr.Write();
2275 
2276 
2277 
2278  double nominal_xsec_chi2 = 0.;
2279  double fake_xsec_chi2 = 0.;
2280  int bin_i = 0;
2281  for (auto it = best_fit_xsec_truth/*mean_xsecs*/.begin(); it != best_fit_xsec_truth/*mean_xsecs*/.end();
2282  ++it) {
2283  for (size_t i = 0; i < it->second.size(); ++i) {
2284  double measured_val_i = it->second[i];
2285  double mc_val_i = fNominalXSecs[it->first]->GetBinContent(i+1);
2286  double fake_val_i = ((fDoFakeData /*&& fFakeDataRoutine != "Toy"*/) ?
2287  fFakeDataXSecs[it->first]->GetBinContent(i+1) :
2288  0.);
2289  int bin_j = 0;
2290  for (auto it2 = best_fit_xsec_truth/*mean_xsecs*/.begin();
2291  it2 != best_fit_xsec_truth/*mean_xsecs*/.end();
2292  ++it2) {
2293  for (size_t j = 0; j < it2->second.size(); ++j) {
2294  double measured_val_j = it2->second[j];
2295  double mc_val_j = fNominalXSecs[it2->first]->GetBinContent(j+1);
2296  double fake_val_j = ((fDoFakeData /*&& fFakeDataRoutine != "Toy"*/) ?
2297  fFakeDataXSecs[it2->first]->GetBinContent(j+1) :
2298  0.);
2299  nominal_xsec_chi2 += ((measured_val_i - mc_val_i)*
2300  xsec_cov_matrix[bin_i][bin_j]*
2301  (measured_val_j - mc_val_j));
2302  if (fDoFakeData /*&& fFakeDataRoutine != "Toy"*/) {
2303  fake_xsec_chi2 += ((measured_val_i - fake_val_i)*
2304  xsec_cov_matrix[bin_i][bin_j]*
2305  (measured_val_j - fake_val_j));
2306  }
2307  ++bin_j;
2308  }
2309  //++bin_j;
2310  }
2311  ++bin_i;
2312  }
2313  //++bin_i;
2314  }
2315 
2316  int bin_count = 0;
2317  for (auto it = data_hists.begin(); it != data_hists.end(); ++it) {
2318  int selection_ID = it->first;
2319  std::vector<TH1*> hists = throw_hists.at(selection_ID);
2320 
2321  std::string canvas_name = "cThrow" +
2322  fDataSet.GetSelectionName(selection_ID);
2323  TCanvas cThrow(canvas_name.c_str(), "");
2324  cThrow.SetTicks();
2325 
2326  std::string name = "Throw" + fDataSet.GetSelectionName(selection_ID);
2327  auto data_hist = it->second;
2328  std::vector<double> xs, xs_width;
2329  std::vector<double> ys, errs;
2330  for (int i = 1;
2331  i <= fBestFitSelectionHists[it->first]->GetNbinsX(); ++i) {
2332  ys.push_back(
2333  fBestFitSelectionHists[it->first]->GetBinContent(i));
2334  errs.push_back(
2335  sqrt(selection_cov.GetBinContent(bin_count+i, bin_count+i)));
2336  xs.push_back(data_hist->GetBinCenter(i));
2337  xs_width.push_back(data_hist->GetBinWidth(i)/2.);
2338  }
2339 
2340  TGraphAsymmErrors throw_gr(data_hist->GetNbinsX(),
2341  &xs[0], &ys[0],
2342  &xs_width[0], &xs_width[0], &errs[0], &errs[0]);
2343 
2344  throw_gr.SetTitle(fDataSet.GetSelectionName(selection_ID).c_str());
2345  throw_gr.GetXaxis()->SetTitle(data_hist->GetXaxis()->GetTitle());
2346  throw_gr.SetFillStyle(3144);
2347  throw_gr.SetFillColor(kRed);
2348  data_hist->Draw();
2349  throw_gr.Draw("same a2");
2350  data_hist->Draw("same e1");
2351  fOutputFile.cd("Throws");
2352  cThrow.Write();
2353 
2354  bin_count += data_hist->GetNbinsX();
2355  }
2356 
2357  bin_count = 0;
2358  for (auto it = truth_throw_hists.begin(); it != truth_throw_hists.end(); ++it) {
2359  int sample_ID = it->first;
2360 
2361  std::vector<double> xs, xs_width;
2362  for (size_t i = 0; i < sample_bins[it->first]; ++i) {
2363  xs.push_back(i + 0.5);
2364  xs_width.push_back(.5);
2365  }
2366 
2367  std::string name = "hNominal" + fSamples[sample_ID][0][0].GetName();
2368  TH1D temp_nominal(name.c_str(), "", xs.size(), 0, xs.size());
2369  std::string dummy_name = "hDummy" + fSamples[sample_ID][0][0].GetName();
2370  TH1D dummy(dummy_name.c_str(), "", xs.size(), 0, xs.size());
2371  if (fIsSignalSample[sample_ID]) {
2372  dummy.GetXaxis()->SetBinLabel(1, "Underflow");
2373  dummy.GetXaxis()->SetBinLabel(dummy.GetNbinsX(), "Overflow");
2374  for (int i = 2; i < dummy.GetNbinsX(); ++i) {
2375  std::ostringstream low_stream, high_stream;
2376  low_stream.precision(2);
2377  high_stream.precision(2);
2378  low_stream << std::fixed <<
2379  fSamples[sample_ID][0][i-1].RangeLowEnd();
2380  high_stream
2381  << std::fixed <<
2382  fSamples[sample_ID][0][i-1].RangeHighEnd();
2383  std::string label = low_stream.str();
2384  label += " - ";
2385  label += high_stream.str();
2386  //std::string label
2387  // = std::to_string(fSamples[sample_ID][0][i-1].RangeLowEnd());
2388  //label += " - ";
2389  //label += std::to_string(fSamples[sample_ID][0][i-1].RangeHighEnd());
2390  dummy.GetXaxis()->SetBinLabel(i, label.c_str());
2391  }
2392  }
2393  else {
2394  dummy.GetXaxis()->SetBinLabel(1, "");
2395  }
2396 
2397 
2398  std::vector<std::vector<ThinSliceSample>> & samples_vec_2D
2399  = fSamples[sample_ID];
2400  for (size_t i = 0; i < samples_vec_2D.size(); ++i) {
2401  for (size_t j = 0; j < samples_vec_2D[i].size(); ++j) {
2402  temp_nominal.AddBinContent(j+1, samples_vec_2D[i][j].GetNominalFlux());
2403  }
2404  }
2405 
2406  double max = -999.;
2407  for (size_t i = 0; i < sample_bins[it->first]; ++i) {
2408  if ((fBestFitTruthVals[sample_ID][i] + best_fit_errs[sample_ID][i]) > max)
2409  max = (fBestFitTruthVals[sample_ID][i] + best_fit_errs[sample_ID][i]);
2410 
2411  if (temp_nominal.GetBinContent(i+1) > max)
2412  max = temp_nominal.GetBinContent(i+1);
2413  }
2414 
2415  fOutputFile.cd("Throws");
2416  std::string canvas_name = "cTruthThrow" + fSamples[sample_ID][0][0].GetName();
2417  TCanvas cThrow(canvas_name.c_str(), "");
2418  cThrow.SetTicks();
2419  TGraphAsymmErrors throw_gr(xs.size(),
2420  &xs[0], &fBestFitTruthVals[it->first][0],
2421  &xs_width[0], &xs_width[0],
2422  &best_fit_errs[it->first][0],
2423  &best_fit_errs[it->first][0]);
2424  throw_gr.SetFillStyle(3144);
2425  throw_gr.SetFillColor(kRed);
2426  throw_gr.SetMinimum(0.);
2427  throw_gr.SetMaximum(1.5*max);
2428  dummy.SetMaximum(1.5*max);
2429  dummy.Draw();
2430  //throw_gr.Draw("a2");
2431  throw_gr.Draw("same 2");
2432  throw_gr.Draw("p");
2433 
2434  temp_nominal.SetMarkerColor(kBlue);
2435  temp_nominal.SetMarkerStyle(20);
2436  temp_nominal.Draw("same p");
2437 
2438  TLegend leg;
2439  leg.AddEntry(&throw_gr, "Throws", "lpf");
2440  leg.AddEntry(&temp_nominal, "Nominal", "p");
2441 
2442  if (fDoFakeData && fFakeDataRoutine != "Toy") {
2443  name = "hVaried" + fSamples[sample_ID][0][0].GetName();
2444  TH1D * temp_varied = (TH1D*)temp_nominal.Clone(name.c_str());
2445  for (size_t i = 0; i < xs.size(); ++i) {
2446  temp_varied->SetBinContent(
2447  i+1, temp_varied->GetBinContent(i+1)*fFakeDataScales[sample_ID][i]);
2448  }
2449  temp_varied->SetMarkerColor(kBlack);
2450  temp_varied->SetMarkerStyle(20);
2451  temp_varied->Draw("same p");
2452  leg.AddEntry(temp_varied, "Fake Data", "p");
2453  }
2454 
2455  leg.Draw();
2456  cThrow.Write();
2457 
2458  bin_count += xs.size();
2459  }
2460 
2461  for (auto it = best_fit_xsec_truth.begin(); it != best_fit_xsec_truth.end();
2462  ++it) {
2463  int sample_ID = it->first;
2464 
2465  std::vector<double> xs, xs_width;
2466  for (size_t i = 0; i < xsec_bins[sample_ID]; ++i) {
2467  xs.push_back(i + 0.5);
2468  xs_width.push_back(.5);
2469  }
2470 
2471  fOutputFile.cd("Throws");
2472  std::string canvas_name = "cXSecThrow" + fSamples[sample_ID][0][0].GetName();
2473  std::string gr_name = "grXSecThrow" + fSamples[sample_ID][0][0].GetName();
2474  TCanvas cThrow(canvas_name.c_str(), "");
2475  cThrow.SetTicks();
2476 
2477  std::string dummy_name = "hDummyXSec" + fSamples[sample_ID][0][0].GetName();
2478  TH1D dummy(dummy_name.c_str(), "", xs.size(), 0, xs.size());
2479  for (int i = 1; i <= dummy.GetNbinsX(); ++i) {
2480  std::ostringstream low_stream, high_stream;
2481  low_stream.precision(2);
2482  high_stream.precision(2);
2483  low_stream << std::fixed <<
2484  fSamples[sample_ID][0][i].RangeLowEnd();
2485  high_stream
2486  << std::fixed <<
2487  fSamples[sample_ID][0][i].RangeHighEnd();
2488  std::string label = low_stream.str();
2489  label += " - ";
2490  label += high_stream.str();
2491  //std::string label
2492  // = std::to_string(fSamples[sample_ID][0][i-1].RangeLowEnd());
2493  //label += " - ";
2494  //label += std::to_string(fSamples[sample_ID][0][i-1].RangeHighEnd());
2495  dummy.GetXaxis()->SetBinLabel(i, label.c_str());
2496  }
2497 
2498  TGraphAsymmErrors throw_gr(xs.size(),
2499  &xs[0], &best_fit_xsec_truth/*mean_xsecs*/[it->first][0],
2500  &xs_width[0], &xs_width[0],
2501  &best_fit_xsec_errs[it->first][0],
2502  &best_fit_xsec_errs[it->first][0]);
2503  //TGraph mean_gr(xs.size(), &xs[0], &mean_xsecs[it->first][0]);
2504  double max = -999.;
2505  for (size_t i = 0; i < best_fit_xsec_truth[it->first].size(); ++i) {
2506  if ((best_fit_xsec_truth[it->first][i] +
2507  best_fit_xsec_errs[it->first][0]) > max) {
2508  max = (best_fit_xsec_truth[it->first][i] +
2509  best_fit_xsec_errs[it->first][0]);
2510  }
2511  }
2512  //throw_gr.SetFillStyle(3144);
2513  //throw_gr.SetFillColor(kRed);
2514  throw_gr.SetMarkerStyle(20);
2515  throw_gr.SetMinimum(0.);
2516  dummy.SetMaximum(1.5*max);
2517  dummy.SetMinimum(0.);
2518  dummy.Draw();
2519  dummy.GetXaxis()->SetTitle("Kinetic Energy (MeV)");
2520  dummy.GetYaxis()->SetTitle("#sigma (mb)");
2521  //throw_gr.Draw("same 2");
2522  throw_gr.Draw("same 2");
2523  //mean_gr.Draw("same 2");
2524  throw_gr.Draw("p");
2525 
2526  std::vector<double> nominal_xsec_vals;
2527  for (size_t i = 0; i < xs.size(); ++i) {
2528  nominal_xsec_vals.push_back(
2529  fNominalXSecs[it->first]->GetBinContent(i+1));
2530  }
2531  TGraph nominal_gr(xs.size(), &xs[0], &nominal_xsec_vals[0]);
2532  nominal_gr.SetMarkerColor(kBlue);
2533  nominal_gr.SetMarkerStyle(20);
2534  nominal_gr.Draw("same p");
2535 
2536  TLegend leg;
2537  leg.AddEntry(&throw_gr, "Measured", "lpf");
2538  //leg.AddEntry(&mean_gr, "Mean", "lpf");
2539  leg.AddEntry(&nominal_gr, "Nominal", "p");
2540 
2541  if (fDoFakeData /*&& fFakeDataRoutine != "Toy"*/) {
2542  //std::cout << "Plotting fake data" << std::endl;
2543  std::vector<double> fake_xsec_vals;
2544  //std::cout << xs.size() << std::endl;
2545  for (size_t i = 0; i < xs.size(); ++i) {
2546 
2547  //std::cout << fFakeDataXSecs[it->first] << std::endl;
2548  //std::cout << "Testing" << std::endl;
2549  //std::cout << fFakeDataXSecs[it->first]->GetNbinsX() << std::endl;
2550  fake_xsec_vals.push_back(
2551  fFakeDataXSecs[it->first]->GetBinContent(i+1));
2552  //std::cout << "Adding " << fake_xsec_vals.back() << std::endl;
2553  }
2554  TGraph * fake_gr = new TGraph(xs.size(), &xs[0], &fake_xsec_vals[0]);
2555  fake_gr->Write();
2556  fake_gr->SetMarkerColor(kRed);
2557  fake_gr->SetMarkerStyle(20);
2558  fake_gr->Draw("same p");
2559  leg.AddEntry(fake_gr, "Fake Data", "p");
2560  }
2561 
2562  TString chi2_str;
2563  chi2_str.Form("Nominal #chi^{2} = %.2f", nominal_xsec_chi2);
2564  leg.AddEntry((TObject*)0x0, chi2_str, "");
2565  //std::string chi2_str = "Nominal #chi^{2} = " +
2566  // std::to_string(nominal_xsec_chi2);
2567  //leg.AddEntry((TObject*)0x0, chi2_str.c_str(), "");
2568  if (fDoFakeData /*&& fFakeDataRoutine != "Toy"*/) {
2569  //std::string fake_chi2_str = "Fake Data #chi^{2} = " +
2570  // std::to_string(fake_xsec_chi2);
2571  //leg.AddEntry((TObject*)0x0, fake_chi2_str.c_str(), "");
2572  TString fake_chi2_str;
2573  fake_chi2_str.Form("Fake Data #chi^{2} = %.2f", fake_xsec_chi2);
2574  leg.AddEntry((TObject*)0x0, fake_chi2_str, "");
2575  }
2576  leg.Draw("same");
2577  gPad->RedrawAxis();
2578  cThrow.Write();
2579  throw_gr.Write(gr_name.c_str());
2580  }
2581 
2582  TVectorD fake_chi2_out(1);
2583  fake_chi2_out[0] = fake_xsec_chi2;
2584  fake_chi2_out.Write("FakeDataChi2");
2585 
2586  TVectorD nominal_chi2_out(1);
2587  nominal_chi2_out[0] = nominal_xsec_chi2;
2588  nominal_chi2_out.Write("NominalChi2");
2589 }
static QCString name
Definition: declinfo.cpp:673
std::map< int, TH1 * > fBestFitXSecs
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::map< int, TH1 * > fFakeDataXSecs
std::string string
Definition: nybbler.cc:12
std::map< int, TH1D * > fBestFitSelectionHists
std::map< int, TH1 * > fNominalXSecs
std::map< int, TH1 * > & GetSelectionHists()
std::map< int, std::vector< double > > fBestFitTruthVals
std::map< int, TH1 * > fBestFitIncs
std::map< int, TH1 * > & GetRebinnedSelectionHists()
static int max(int a, int b)
auto selection_ID
std::string & GetSelectionName(int id)
cet::LibraryManager dummy("noplugin")
std::map< int, std::vector< double > > fFakeDataScales
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
std::map< int, bool > fIsSignalSample
static constexpr double ys
Definition: Units.h:103
static QCString * s
Definition: config.cpp:1042
void protoana::PDSPThinSliceFitter::Pulls ( )
private

Definition at line 1278 of file PDSPThinSliceFitter.cxx.

1278  {
1280 
1281  //TH2D hPulls("Pulls", "", total_parameters, 0, total_parameters, 100, -5, 5 );
1282  TTree pulls_tree("pulls_tree", "");
1283  std::vector<double> pulls;
1284  pulls_tree.Branch("pulls", &pulls);
1285  std::vector<double> vals;
1286  pulls_tree.Branch("vals", &vals);
1287 
1289  for (size_t i = 0; i < fNPulls; ++i) {
1290  std::cout << "Fitting for pulls: " << i << "/" << fNPulls << std::endl;
1291  fMinimizer->SetVariableValues(&fMinimizerInitVals[0]);
1292  std::cout << "\tGenerating fluctuation" << std::endl;
1294  std::cout << "\tDone" << std::endl;
1295  int fit_status = fMinimizer->Minimize();
1296  if (fit_status) {
1297  std::cout << "Fit successful" << std::endl;
1298  pulls.clear();
1299  vals.clear();
1300  for (size_t j = 0; j < total_parameters; ++j) {
1301  double pull_val = (fMinimizer->X()[j] - fMinimizerInitVals[j])/
1302  sqrt(fMinimizer->CovMatrix(j, j));
1303  pulls.push_back(pull_val);
1304  vals.push_back(fMinimizer->X()[j]);
1305  //hPulls.Fill(j+.5, pull_val);
1306  }
1307  pulls_tree.Fill();
1308  }
1309  }
1310  fOutputFile.cd();
1311  pulls_tree.Write();
1312  //hPulls.Write();
1313 
1314 }
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
std::vector< double > fMinimizerInitVals
QTextStream & endl(QTextStream &s)
void protoana::PDSPThinSliceFitter::RunFitAndSave ( )

Definition at line 1316 of file PDSPThinSliceFitter.cxx.

1316  {
1317  TFile fMCFile(fMCFileName.c_str(), "OPEN");
1318  fMCTree = (TTree*)fMCFile.Get(fTreeName.c_str());
1319 
1320 
1322  fMinimizer->SetFunction(fFitFunction);
1323 
1324  BuildDataHists();
1325  if (fDoFluctuateStats) {
1328  }
1329 
1330  ScaleMCToData();
1331 
1332  SaveMCSamples();
1333 
1334  TDirectory * top_dir = fOutputFile.GetDirectory("");
1335  TDirectory * xsec_dir = fOutputFile.mkdir("PreFitXSec");
1336  std::string extra_name = "PreFit";
1337  CompareDataMC(extra_name, xsec_dir, top_dir);
1338 
1339  if (fFitType == "Normal" /*|| fFitType == "Toy"*/) {
1340  NormalFit();
1341  }
1342  else if (fFitType == "Pulls") {
1343  Pulls();
1344  }
1345  else if (fFitType == "None") {
1346 
1347  }
1348  else {
1349  std::string message = "Error: invalid fit type given -- ";
1350  message += fFitType;
1351  throw std::runtime_error(message);
1352  }
1353 
1354 
1355  if (fDoSysts)
1357 
1358  fMCFile.Close();
1359 }
std::string string
Definition: nybbler.cc:12
std::unique_ptr< ROOT::Math::Minimizer > fMinimizer
virtual void WrapUpSysts(TFile &output_file)=0
void CompareDataMC(std::string extra_name, TDirectory *xsec_dir, TDirectory *plot_dir, bool post_fit=false)
void protoana::PDSPThinSliceFitter::SaveMCSamples ( )

Definition at line 751 of file PDSPThinSliceFitter.cxx.

751  {
752  fOutputFile.cd();
753  fOutputFile.mkdir("MC_Samples");
754  fOutputFile.cd("MC_Samples");
755  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
756  for (size_t i = 0; i < it->second.size(); ++i) {
757  auto vec = it->second.at(i);
758  for (size_t j = 0; j < vec.size(); ++j) {
759  const std::map<int, TH1 *> & hists = vec[j].GetSelectionHists();
760  for (auto it2 = hists.begin(); it2 != hists.end(); ++it2) {
761  it2->second->Write();
762  }
763  }
764  }
765  }
766 
767 }
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
void protoana::PDSPThinSliceFitter::ScaleDataToNorm ( )

Definition at line 553 of file PDSPThinSliceFitter.cxx.

553  {
554  std::map<int, TH1 *> & selected_hists = fDataSet.GetSelectionHists();
555  for (auto it = selected_hists.begin(); it != selected_hists.end(); ++it) {
556  it->second->Scale(fDataNorm/fDataFlux);
557  }
558 
559  for (auto it = fFakeFluxesBySample.begin();
560  it != fFakeFluxesBySample.end(); ++it) {
561  for (size_t i = 0; i < it->second.size(); ++i) {
562  for (size_t j = 0; j < it->second[i].size(); ++j) {
563  it->second[i][j] *= fDataNorm/fDataFlux;
564  }
565  }
566  }
567 
568  fDataFlux = fDataNorm;
569 
570 
571 }
std::map< int, TH1 * > & GetSelectionHists()
std::map< int, std::vector< std::vector< double > > > fFakeFluxesBySample
void protoana::PDSPThinSliceFitter::ScaleMCToData ( )

Definition at line 573 of file PDSPThinSliceFitter.cxx.

573  {
574  double total_nominal = 0.;
575 
576  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
577  for (size_t i = 0; i < it->second.size(); ++i) {
578  for (size_t j = 0; j < it->second[i].size(); ++j) {
579  total_nominal += it->second[i][j].GetNominalFlux();
580  }
581  }
582  }
583 
584  /*
585  for (auto it = fNominalFluxes.begin(); it != fNominalFluxes.end(); ++it) {
586  total_nominal += it->second;
587  }*/
588 
589  fMCDataScale = fDataFlux/total_nominal;
590  for (auto it = fNominalFluxes.begin(); it != fNominalFluxes.end(); ++it) {
591  it->second *= fMCDataScale;
592  }
593 
594  for (auto it = fFluxesBySample.begin();
595  it != fFluxesBySample.end(); ++it) {
596  for (size_t i = 0; i < it->second.size(); ++i) {
597  for (size_t j = 0; j < it->second[i].size(); ++j) {
598  it->second[i][j] *= fMCDataScale;
599  }
600  }
601  }
602 
603  std::cout << "MC Data Scale: " << fMCDataScale << std::endl;
604 
605  double new_total_mc = 0., new_total_data = 0., mc_flux = 0.;
606  double data_bins = 0., mc_bins = 0.;
607  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
608  for (size_t i = 0; i < it->second.size(); ++i) {
609  for (size_t j = 0; j < it->second[i].size(); ++j) {
610  it->second[i][j].SetDataMCScale(fMCDataScale);
611 
612  mc_flux += it->second[i][j].GetNominalFlux();
613  const std::map<int, TH1 *> & hists
614  = it->second[i][j].GetSelectionHists();
615  for (auto it2 = hists.begin(); it2 != hists.end(); ++it2) {
616  new_total_mc += it2->second->Integral();
617  for (int k = 1; k <= it2->second->GetNbinsX(); ++k) {
618  mc_bins += it2->second->GetBinContent(k);
619  }
620  }
621  }
622  }
623  }
624  std::map<int, TH1 *> & selected_hists = fDataSet.GetSelectionHists();
625  for (auto it = selected_hists.begin(); it != selected_hists.end(); ++it) {
626  new_total_data += it->second->Integral();
627  for (int k = 1; k <= it->second->GetNbinsX(); ++k) {
628  data_bins += it->second->GetBinContent(k);
629  }
630  }
631 
632  /*
633  std::cout << "Data: " << std::setprecision(20) << fDataFlux << std::endl;
634  std::cout << "MC: " << mc_flux << std::endl;
635  std::cout << "MC hists: " << new_total_mc << std::endl;
636  std::cout << "Data hists: " << new_total_data << std::endl;
637  std::cout << "MC bins: " << mc_bins << std::endl;
638  std::cout << "Data bins: " << data_bins << std::endl;
639 
640  std::cout << "Data/MC: " << data_bins/mc_bins << std::endl;*/
641  double new_factor = data_bins/mc_bins;
642 
643  //mc_bins = 0.;
644  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
645  for (size_t i = 0; i < it->second.size(); ++i) {
646  for (size_t j = 0; j < it->second[i].size(); ++j) {
647 
648  it->second[i][j].ExtraFactor(new_factor);
649 
650  /*
651  const std::map<int, TH1 *> & hists
652  = it->second[i][j].GetSelectionHists();
653  for (auto it2 = hists.begin(); it2 != hists.end(); ++it2) {
654  //new_total_mc += it2->second->Integral();
655  for (int k = 1; k <= it2->second->GetNbinsX(); ++k) {
656  mc_bins += it2->second->GetBinContent(k);
657  }
658  }*/
659  }
660  }
661  }
662 
663 
664  //std::cout << "New MC bins: " << mc_bins << std::endl;
665 
666 }
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples
std::map< int, TH1 * > & GetSelectionHists()
std::map< int, std::vector< std::vector< double > > > fFluxesBySample
QTextStream & endl(QTextStream &s)
std::map< int, double > fNominalFluxes
void protoana::PDSPThinSliceFitter::SetBestFit ( )
private

Definition at line 2034 of file PDSPThinSliceFitter.cxx.

2034  {
2035  for (auto it = fSamples.begin(); it != fSamples.end(); ++it) {
2036  for (size_t i = 0; i < it->second.size(); ++i) {
2037  for (size_t j = 0; j < it->second[i].size(); ++j) {
2038  it->second[i][j].SetBestFit();
2039  }
2040  }
2041  }
2042 }
std::map< int, std::vector< std::vector< ThinSliceSample > > > fSamples

Member Data Documentation

bool protoana::PDSPThinSliceFitter::fAddRegTerm
private

Definition at line 126 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fAddSystTerm
private

Definition at line 126 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fAnalysis
private

Definition at line 162 of file PDSPThinSliceFitter.h.

fhicl::ParameterSet protoana::PDSPThinSliceFitter::fAnalysisOptions
private

Definition at line 163 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fBeamEnergyBins
private

Definition at line 184 of file PDSPThinSliceFitter.h.

std::map<int, double> protoana::PDSPThinSliceFitter::fBestFitFluxPars
private

Definition at line 137 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fBestFitIncs
private

Definition at line 139 of file PDSPThinSliceFitter.h.

std::map<int, TH1D*> protoana::PDSPThinSliceFitter::fBestFitSelectionHists
private

Definition at line 142 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<double> > protoana::PDSPThinSliceFitter::fBestFitSignalPars
private

Definition at line 135 of file PDSPThinSliceFitter.h.

std::map<std::string, ThinSliceSystematic> protoana::PDSPThinSliceFitter::fBestFitSystPars
private

Definition at line 136 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<double> > protoana::PDSPThinSliceFitter::fBestFitTruthVals
private

Definition at line 143 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fBestFitXSecs
private

Definition at line 139 of file PDSPThinSliceFitter.h.

std::map<std::string, size_t> protoana::PDSPThinSliceFitter::fCovarianceBins
private

Definition at line 125 of file PDSPThinSliceFitter.h.

TMatrixD* protoana::PDSPThinSliceFitter::fCovMatrix
private

Definition at line 128 of file PDSPThinSliceFitter.h.

TMatrixD * protoana::PDSPThinSliceFitter::fCovMatrixDisplay
private

Definition at line 128 of file PDSPThinSliceFitter.h.

TFile protoana::PDSPThinSliceFitter::fDataFile
private

Definition at line 81 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fDataFileName
private

Definition at line 149 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fDataFlux
private

Definition at line 109 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fDataNorm
private

Definition at line 182 of file PDSPThinSliceFitter.h.

ThinSliceDataSet protoana::PDSPThinSliceFitter::fDataSet
private

Definition at line 77 of file PDSPThinSliceFitter.h.

TTree* protoana::PDSPThinSliceFitter::fDataTree
private

Definition at line 82 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDo1DShifts
private

Definition at line 167 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoFakeData
private

Definition at line 167 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoFluctuateStats
private

Definition at line 171 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoScaleDataToNorm
private

Definition at line 181 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoScans
private

Definition at line 167 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoSysts
private

Definition at line 167 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDoThrows
private

Definition at line 167 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fDrawXSecUnderflow
private

Definition at line 186 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fDriverName
private

Definition at line 161 of file PDSPThinSliceFitter.h.

std::vector<ThinSliceEvent> protoana::PDSPThinSliceFitter::fEvents
private

Definition at line 145 of file PDSPThinSliceFitter.h.

std::vector<ThinSliceEvent> protoana::PDSPThinSliceFitter::fFakeDataEvents
private

Definition at line 145 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fFakeDataIncs
private

Definition at line 140 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fFakeDataRoutine
private

Definition at line 170 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<double> > protoana::PDSPThinSliceFitter::fFakeDataScales
private

Definition at line 134 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fFakeDataXSecs
private

Definition at line 140 of file PDSPThinSliceFitter.h.

std::map<int, double> protoana::PDSPThinSliceFitter::fFakeFluxes
private

Definition at line 105 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<std::vector<double> > > protoana::PDSPThinSliceFitter::fFakeFluxesBySample
private

Definition at line 106 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<std::vector<ThinSliceSample> > > protoana::PDSPThinSliceFitter::fFakeSamples
private

Definition at line 75 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fFillIncidentInFunction = false
private

Definition at line 174 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fFitFlux
private

Definition at line 177 of file PDSPThinSliceFitter.h.

ROOT::Math::Functor protoana::PDSPThinSliceFitter::fFitFunction
private

Definition at line 84 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fFitType = "Normal"
private

Definition at line 179 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fFitUnderOverflow = false
private

Definition at line 175 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fFixVariables
private

Definition at line 168 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<std::vector<double> > > protoana::PDSPThinSliceFitter::fFluxesBySample
private

Definition at line 106 of file PDSPThinSliceFitter.h.

std::map<int, std::string> protoana::PDSPThinSliceFitter::fFluxParameterNames
private

Definition at line 117 of file PDSPThinSliceFitter.h.

std::map<int, double> protoana::PDSPThinSliceFitter::fFluxParameters
private

Definition at line 116 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<int> > protoana::PDSPThinSliceFitter::fFluxParsToSamples
private

Definition at line 108 of file PDSPThinSliceFitter.h.

std::map<int, std::string> protoana::PDSPThinSliceFitter::fFluxTypes
private

Definition at line 153 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fIncidentRecoBins
private

Definition at line 184 of file PDSPThinSliceFitter.h.

std::vector<int> protoana::PDSPThinSliceFitter::fIncidentSamples
private

Definition at line 185 of file PDSPThinSliceFitter.h.

TDecompChol* protoana::PDSPThinSliceFitter::fInputChol
private

Definition at line 129 of file PDSPThinSliceFitter.h.

std::map<int, bool> protoana::PDSPThinSliceFitter::fIsSignalSample
private

Definition at line 78 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fLowerLimit
private

Definition at line 157 of file PDSPThinSliceFitter.h.

int protoana::PDSPThinSliceFitter::fMaxCalls
private

Definition at line 154 of file PDSPThinSliceFitter.h.

int protoana::PDSPThinSliceFitter::fMaxIterations
private

Definition at line 154 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fMaxRethrows
private

Definition at line 178 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fMCDataScale = 1.
private

Definition at line 110 of file PDSPThinSliceFitter.h.

TFile protoana::PDSPThinSliceFitter::fMCFile
private

Definition at line 79 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fMCFileName
private

Definition at line 148 of file PDSPThinSliceFitter.h.

TTree* protoana::PDSPThinSliceFitter::fMCTree
private

Definition at line 80 of file PDSPThinSliceFitter.h.

std::vector<int> protoana::PDSPThinSliceFitter::fMeasurementSamples
private

Definition at line 185 of file PDSPThinSliceFitter.h.

std::unique_ptr<ROOT::Math::Minimizer> protoana::PDSPThinSliceFitter::fMinimizer
private

Definition at line 85 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fMinimizerInitVals
private

Definition at line 86 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fMultinomial
private

Definition at line 166 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fNFitSteps = 0
private

Definition at line 155 of file PDSPThinSliceFitter.h.

std::map<int, double> protoana::PDSPThinSliceFitter::fNominalFluxes
private

Definition at line 105 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fNominalIncs
private

Definition at line 138 of file PDSPThinSliceFitter.h.

std::map<int, TH1*> protoana::PDSPThinSliceFitter::fNominalXSecs
private

Definition at line 138 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fNPulls
private

Definition at line 180 of file PDSPThinSliceFitter.h.

unsigned int protoana::PDSPThinSliceFitter::fNScanSteps
private

Definition at line 156 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fNThrows
private

Definition at line 178 of file PDSPThinSliceFitter.h.

TFile protoana::PDSPThinSliceFitter::fOutputFile
private

Definition at line 83 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fParLimits
private

Definition at line 123 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fParLimitsUp
private

Definition at line 123 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fPitch
private

Definition at line 164 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fPlotRebinned
private

Definition at line 159 of file PDSPThinSliceFitter.h.

std::vector<std::pair<int, int> > protoana::PDSPThinSliceFitter::fPlotStyle
private

Definition at line 158 of file PDSPThinSliceFitter.h.

TH1D protoana::PDSPThinSliceFitter::fPreFitParsNormal
private

Definition at line 87 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fRandomStart
private

Definition at line 160 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fRegFactor = 0.
private

Definition at line 127 of file PDSPThinSliceFitter.h.

TRandom3 protoana::PDSPThinSliceFitter::fRNG
private

Definition at line 133 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fRunHesse
private

Definition at line 167 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<std::vector<ThinSliceSample> > > protoana::PDSPThinSliceFitter::fSamples
private

Definition at line 75 of file PDSPThinSliceFitter.h.

std::vector<fhicl::ParameterSet> protoana::PDSPThinSliceFitter::fSampleSets
private

Definition at line 152 of file PDSPThinSliceFitter.h.

std::vector<fhicl::ParameterSet> protoana::PDSPThinSliceFitter::fSelectionSets
private

Definition at line 151 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<double> > protoana::PDSPThinSliceFitter::fSignalBins
private

Definition at line 187 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<std::string> > protoana::PDSPThinSliceFitter::fSignalParameterNames
private

Definition at line 113 of file PDSPThinSliceFitter.h.

std::map<int, std::vector<double> > protoana::PDSPThinSliceFitter::fSignalParameters
private

Definition at line 112 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fSliceMethod
private

Definition at line 165 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fSplitMC
private

Definition at line 172 of file PDSPThinSliceFitter.h.

int protoana::PDSPThinSliceFitter::fSplitVal = 0
private

Definition at line 173 of file PDSPThinSliceFitter.h.

std::vector<std::string> protoana::PDSPThinSliceFitter::fSystParameterNames
private

Definition at line 122 of file PDSPThinSliceFitter.h.

std::map<std::string, ThinSliceSystematic> protoana::PDSPThinSliceFitter::fSystParameters
private

Definition at line 121 of file PDSPThinSliceFitter.h.

std::map<std::string, double> protoana::PDSPThinSliceFitter::fSystsToFix
private

Definition at line 169 of file PDSPThinSliceFitter.h.

ThinSliceDriver* protoana::PDSPThinSliceFitter::fThinSliceDriver
private

Definition at line 74 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fTolerance
private

Definition at line 157 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fTotalFluxParameters = 0
private

Definition at line 118 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fTotalSignalParameters = 0
private

Definition at line 114 of file PDSPThinSliceFitter.h.

size_t protoana::PDSPThinSliceFitter::fTotalSystParameters = 0
private

Definition at line 124 of file PDSPThinSliceFitter.h.

std::map<std::string, double> protoana::PDSPThinSliceFitter::fToyValues
private

Definition at line 131 of file PDSPThinSliceFitter.h.

std::string protoana::PDSPThinSliceFitter::fTreeName
private

Definition at line 150 of file PDSPThinSliceFitter.h.

std::vector<double> protoana::PDSPThinSliceFitter::fTrueIncidentBins
private

Definition at line 184 of file PDSPThinSliceFitter.h.

double protoana::PDSPThinSliceFitter::fUpperLimit
private

Definition at line 157 of file PDSPThinSliceFitter.h.

bool protoana::PDSPThinSliceFitter::fUseFakeSamples = false
private

Definition at line 176 of file PDSPThinSliceFitter.h.


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