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

#include <ProtoDUNEFit.h>

Public Member Functions

 ProtoDUNEFit ()
 
 ProtoDUNEFit (std::string configPath)
 
virtual ~ProtoDUNEFit ()
 
void BuildWorkspace (TString Outputfile, int analysis=-1)
 
bool ApplySystematicToSample (RooStats::HistFactory::Sample &sample, TH1 *histo, std::vector< TH1 * > systvec, bool hasnormfactor, bool isnorm)
 
bool BuildSystThenApplyToSample (RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, protoana::HistType this_histType, size_t iChan, size_t iTopo, size_t iTruthBin=999)
 
bool BuildSignalSystThenApplyToSample (RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iChan, size_t iTopo, size_t iTruthBin)
 
bool BuildBackgroundSystThenApplyToSample (RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iChan, size_t iTopo)
 
bool BuildIncidentSystThenApplyToSample (RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iTopo)
 
std::vector< std::vector< std::pair< TH1 *, TH1 * > > > BuildIncidentSignalSyst (size_t iTopo)
 
bool ApplyBuiltSystToSample (TH1 *histo, TH1 *high_hist, TH1 *low_hist, RooStats::HistFactory::Sample &sample, std::string syst_name, std::string syst_type, bool hasnormfactor)
 
std::vector< TH1 * > DrawXSecs (RooFitResult *fitresult=0x0)
 
std::vector< TH2 * > DrawSmearingMatrix (RooFitResult *fitresult, bool doPostFit=false)
 

Private Member Functions

bool Configure (std::string configPath)
 
bool FillHistogramVectors_Pions ()
 
bool FillSidebandHistograms_Pions ()
 
void ScaleMCToData (bool data_is_mc=false)
 
void ScaleMuonContent ()
 
void AddSamplesAndChannelsToMeasurement (RooStats::HistFactory::Measurement &meas)
 
void AddIncidentSamplesAndChannelsToMeasurement (RooStats::HistFactory::Measurement &meas)
 
void AddSidebandSamplesAndChannelsToMeasurement (RooStats::HistFactory::Measurement &meas)
 
void DecorateEfficiency (TGraphAsymmErrors *eff, std::string x_title="E_{true} at vertex [MeV]")
 

Private Attributes

std::string _RecoTreeName
 
std::string _Minimizer
 
std::string _TruthTreeName
 
std::vector< double > _RecoBinning
 
std::vector< double > _TruthBinning
 
std::vector< double > _SidebandBinning
 
std::vector< int > _SignalTopology
 
std::vector< int > _BackgroundTopology
 
std::vector< int > _IncidentTopology
 
std::vector< std::string_DataFileNames
 
std::vector< std::string_MCFileNames
 
std::vector< std::string_DataControlSampleFiles
 
std::vector< std::string_MCControlSampleFiles
 
std::vector< std::string_SidebandTopologyName
 
std::vector< std::string_SystFileNames
 
std::vector< std::string_SystToConsider
 
std::vector< std::string_SystType
 
std::vector< std::string_BackgroundTopologyName
 
std::vector< std::string_SignalTopologyName
 
std::vector< std::string_ChannelNames
 
std::vector< std::string_IncidentMCFileNames
 
std::vector< std::string_IncidentDataFileNames
 
std::vector< std::string_IncidentTopologyName
 
std::vector< int > _SidebandTopology
 
int _FitStrategy
 
int _NToys
 
double _IgnoreStatisticalErrorBelow
 
double _IgnoreSystematicErrorBelow
 
bool _EnableMinosError
 
bool _DoAsimovFit
 
bool _EnableStatisticalError
 
bool _EnableSystematicError
 
bool _NormalisedSystematic
 
bool _FitInReco
 
bool _UseComputedSysts
 
std::vector< TH1 * > _bkghistos
 
std::vector< TH1 * > _sighistos
 
std::vector< TH1 * > _truthsighistos
 
std::vector< TH1 * > _datahistos
 
std::vector< TH1 * > _syst_hists
 
std::vector< size_t > _bkg_chan_index
 
std::vector< size_t > _sig_chan_index
 
std::vector< size_t > _bkg_topo_index
 
std::vector< size_t > _sig_topo_index
 
std::vector< size_t > _sig_truth_index
 
std::vector< size_t > _inc_sig_topo_index
 
std::vector< size_t > _inc_bkg_topo_index
 
std::vector< TH1 * > _incbkghistos
 
std::vector< TH1 * > _incsighistos
 
std::vector< TH1 * > _incdatahistos
 
std::vector< TH1 * > _sideband_hists_mc
 
std::vector< TH1 * > _sideband_hists_data
 
std::vector< TGraphAsymmErrors * > _efficiencyGraphs
 
TGraphAsymmErrors * _incidentEfficiency
 
TH1 * _incidentEfficiencyNum
 
TH1 * _incidentEfficiencyDenom
 
std::vector< TH1 * > _interactingEfficiencyDenoms
 
std::vector< TH1 * > _interactingEfficiencyNums
 
std::vector< TGraphAsymmErrors * > _interactingEfficiencies
 
bool _AddIncidentToMeasurement
 
bool _DoNegativeReco
 
bool _DoScaleMCToData
 
bool _AddSidebandsToMeasurement
 
bool _DoScaleMuonContent
 
bool _RandSigPriors
 
TRandom3 rand
 
std::vector< int > _AddBackgroundFactors
 
std::vector< int > _AddIncidentBackgroundFactors
 
std::vector< int > _enable_bkg_factor
 
std::vector< int > _enable_inc_bkg_factor
 
bool _OnlyDrawXSecs
 
bool _DoDrawXSecs
 
bool _DataIsMC
 
double _EndZCut
 
double _WirePitch
 
double _ScaleFactor = 1.
 
double _IncidentScaleFactor = 1.
 
double _MuonScaleFactor = 1.
 
double _PionScaleFactor = 1.
 

Detailed Description

Definition at line 36 of file ProtoDUNEFit.h.

Constructor & Destructor Documentation

protoana::ProtoDUNEFit::ProtoDUNEFit ( )

Definition at line 46 of file ProtoDUNEFit.cxx.

46  {
47  //********************************************************************
48 
49 }
protoana::ProtoDUNEFit::ProtoDUNEFit ( std::string  configPath)

Definition at line 52 of file ProtoDUNEFit.cxx.

52  {
53  //********************************************************************
54 
55  Configure(configPath);
56 
57 }
bool Configure(std::string configPath)
protoana::ProtoDUNEFit::~ProtoDUNEFit ( )
virtual

Definition at line 60 of file ProtoDUNEFit.cxx.

60  {
61  //********************************************************************
62 
63 }

Member Function Documentation

void protoana::ProtoDUNEFit::AddIncidentSamplesAndChannelsToMeasurement ( RooStats::HistFactory::Measurement &  meas)
private

Definition at line 789 of file ProtoDUNEFit.cxx.

789  {
790  //********************************************************************
791 
792  TString channelname("ChannelIncident");
793  RooStats::HistFactory::Channel channel(channelname.Data());
794  // Add data to channel
796  for(unsigned int j=0; j < _incdatahistos.size(); j++){
797  // Clone histogram
798  TH1D* htemp = (TH1D*)(_incdatahistos.at(j)->Clone());
799 
800  mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Adding Incident dataset " << htemp->GetName() << " to channel " << channelname.Data() << " with " << htemp->Integral() << " events";
801  data.SetHisto(htemp);
802  channel.SetData(data);
803  }
804 
805  // Get histograms with the systematics
806 
807  std::vector<TH1*> systvec;
809  //Just use the first for now
811  std::cout << "Got " << systvec.size() << " syst hists" << std::endl;
812  }
813 
814  // Add bkg samples to channel
815  for(unsigned int j=0; j < _incbkghistos.size(); j++){
816  // Clone histogram
817  TH1D* htemp = (TH1D*)(_incbkghistos.at(j)->Clone());
818 
819  TString hname(htemp->GetName());
820  //if(htemp->GetEntries() == 0 || htemp->Integral() == 0) continue;
821 
822  mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Adding Incident MC sample " << hname.Data() << " to channel " << channelname.Data() << " with " << htemp->Integral() << " events";
823  TString samplename = hname + TString("_sample");
824  RooStats::HistFactory::Sample sample(samplename.Data());
825  sample.SetNormalizeByTheory(true);
826 
827  // Check to enable statistical uncertainty
829  //mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Enable statistical uncertainty";
830  sample.ActivateStatError();
831  TH1* staterrorhisto = protoana::ProtoDUNEFitUtils::GetStatsSystHistogram(htemp);
832  RooStats::HistFactory::StatError& staterror = sample.GetStatError();
833  staterror.SetUseHisto();
834  staterror.SetErrorHist(staterrorhisto);
835  }
836 
837 
838 
840  if (_UseComputedSysts) {
841  ApplySystematicToSample(sample, htemp, systvec, false, false);
842  }
843  else{
844  BuildIncidentSystThenApplyToSample(sample, htemp, false, false,
846  }
847  }
848 
849 
850  // Set histogram for sample
851  sample.SetHisto(htemp);
852 
853  std::cout << "Check: " << j << " " << hname << " " <<
855  if (_enable_inc_bkg_factor[j]) {
856  if (htemp->Integral() > 0) {
857  TString poiname = hname;
858  poiname.ReplaceAll("MC","POI");
859  poiname.ReplaceAll("_Histo","");
860  poiname.ReplaceAll(".","");
861  poiname.ReplaceAll("-","_");
862  poiname.ReplaceAll(channelname.Data(),"");
863  poiname.ReplaceAll("__","_");
864  meas.SetPOI(poiname.Data());
865  //sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 2.0);
866  if (_RandSigPriors) {
867  sample.AddNormFactor(poiname.Data(), rand.Gaus(1.0, .1), 0., 2.);
868  }
869  else {
870  sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 2.0);
871  }
872  mf::LogInfo("BuildMeasurement") << "Sample " << sample.GetName() << " has normalisation parameter " << poiname.Data();
873  }
874  }
875  // Add sample to channel
876  channel.AddSample(sample);
877  }
878 
879  //Needed to make the systs for the bin-by-bin incident sample
880  std::vector<std::vector<std::pair<TH1*, TH1*>>> inc_sig_systs;
882  std::cout << "Incident systs for topo: " << _IncidentTopology[0] << " " <<
884  inc_sig_systs = BuildIncidentSignalSyst(0);
885  }
886 
887  // Add signal samples to channel
888  for(unsigned int j=0; j < _incsighistos.size(); j++){
889  // Clone histogram
890  TH1D* htemp = (TH1D*)(_incsighistos.at(j)->Clone());
891 
892  TString hname(htemp->GetName());
893  //if(htemp->GetEntries() == 0 || htemp->Integral() == 0) continue;
894 
895  mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Adding Incident MC sample" << hname.Data() << " to channel " << channelname.Data() << " with " << htemp->Integral() << " events";
896  TString samplename = hname + TString("_sample");
897  RooStats::HistFactory::Sample sample(samplename.Data());
898  sample.SetNormalizeByTheory(true);
899 
900  // Check to enable statistical uncertainty
902  //mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Enable statistical uncertainty";
903  sample.ActivateStatError();
904  TH1* staterrorhisto = protoana::ProtoDUNEFitUtils::GetStatsSystHistogram(htemp);
905  RooStats::HistFactory::StatError& staterror = sample.GetStatError();
906  staterror.SetUseHisto();
907  staterror.SetErrorHist(staterrorhisto);
908  }
909 
910 
912  if (_UseComputedSysts) {
913  ApplySystematicToSample(sample, htemp, systvec, false, false);
914  }
915  else{
916  for (size_t k = 0; k < _SystToConsider.size(); ++k) {
917  std::cout << "Applying syst " << _SystToConsider[k] << " to " <<
918  htemp->GetName() << " " << inc_sig_systs[k][j].first <<
919  inc_sig_systs[k][j].second << std::endl;
920  ApplyBuiltSystToSample(htemp, inc_sig_systs[k][j].first,
921  inc_sig_systs[k][j].second, sample,
922  _SystToConsider[k], _SystType[k], false);
923  }
924  }
925  }
926 
927 
928  // Set histogram for sample
929  sample.SetHisto(htemp);
930 
931  // Add POI to this sample
932  TString poiname = hname;
933  poiname.ReplaceAll("MC","POI");
934  poiname.ReplaceAll("_Histo","");
935  poiname.ReplaceAll(".0","");
936  poiname.ReplaceAll("-","_");
937  meas.SetPOI(poiname.Data()); // AddPOI would also work
938  //sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
939  if (_RandSigPriors) {
940  sample.AddNormFactor(poiname.Data(), rand.Gaus(1.0, .1), 0., 100.);
941  }
942  else {
943  sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
944  }
945  mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Sample " << sample.GetName() << " has normalisation parameter " << poiname.Data();
946 
947  // Add sample to channel
948  channel.AddSample(sample);
949  }
950 
951  // Statistical uncertainty less than 1% is ignored
952  channel.SetStatErrorConfig(_IgnoreStatisticalErrorBelow,"Poisson"); // Poisson or Gaussian
953 
954  // Add channel to measurement
955  mf::LogInfo("AddIncidentSamplesAndChannelsToMeasurement") << "Adding channel " << channel.GetName() << " to measurement " << meas.GetName();
956  meas.AddChannel(channel);
957 
958 }
std::vector< std::string > _SystFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::vector< std::pair< TH1 *, TH1 * > > > BuildIncidentSignalSyst(size_t iTopo)
bool ApplySystematicToSample(RooStats::HistFactory::Sample &sample, TH1 *histo, std::vector< TH1 * > systvec, bool hasnormfactor, bool isnorm)
TH1 * GetStatsSystHistogram(TH1 *nominal)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< TH1 * > _incsighistos
Definition: ProtoDUNEFit.h:129
uint8_t channel
Definition: CRTFragment.hh:201
std::vector< TH1 * > _incdatahistos
Definition: ProtoDUNEFit.h:129
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< int > _IncidentTopology
Definition: ProtoDUNEFit.h:107
Fw2dFFT::Data Data
double _IgnoreStatisticalErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
bool BuildIncidentSystThenApplyToSample(RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iTopo)
std::vector< int > _enable_inc_bkg_factor
Definition: ProtoDUNEFit.h:147
std::vector< TH1 * > _incbkghistos
Definition: ProtoDUNEFit.h:129
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
std::vector< size_t > _inc_bkg_topo_index
Definition: ProtoDUNEFit.h:127
QTextStream & endl(QTextStream &s)
bool ApplyBuiltSystToSample(TH1 *histo, TH1 *high_hist, TH1 *low_hist, RooStats::HistFactory::Sample &sample, std::string syst_name, std::string syst_type, bool hasnormfactor)
std::vector< TH1 * > GetSystHistograms(std::string name)
void protoana::ProtoDUNEFit::AddSamplesAndChannelsToMeasurement ( RooStats::HistFactory::Measurement &  meas)
private

Definition at line 616 of file ProtoDUNEFit.cxx.

616  {
617  //********************************************************************
618 
619  const int nmcchannels = _MCFileNames.size();
620  for(int i=0; i < nmcchannels; i++){
621  TString channelname = Form("Channel%s", _ChannelNames[i].c_str());
622  RooStats::HistFactory::Channel channel(channelname.Data());
623  // Add data to channel
625  for(unsigned int j=0; j < _datahistos.size(); j++){
626  // Clone histogram
627  TH1D* htemp = (TH1D*)(_datahistos.at(j)->Clone());
628 
629  TString hname(htemp->GetName());
630  if(hname.Contains(channelname.Data()) && hname.Contains("Data")){
631  mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Adding dataset " <<
632  hname.Data() << " to channel " << channelname.Data() << " with " <<
633  htemp->Integral() << " events";
634 
635  data.SetHisto(htemp);
636  channel.SetData(data);
637  }
638  }
639 
640  // Get histograms with the systematics
641 
642  std::vector<TH1*> systvec;
644  std::cout << "Attempting to get systs from " << _SystFileNames[i] << std::endl;
646  std::cout << "Got " << systvec.size() << " systematic histograms" << std::endl;
647  }
648 
649  // Add bkg samples to channel
650  for(unsigned int j=0; j < _bkghistos.size(); j++){
651  // Clone histogram
652  TH1D* htemp = (TH1D*)(_bkghistos.at(j)->Clone());
653 
654  TString hname(htemp->GetName());
655  //if(htemp->GetEntries() == 0 || htemp->Integral() == 0) continue;
656  if(hname.Contains(channelname.Data()) && hname.Contains("MC")){
657  mf::LogInfo("AddSamplesAndChannelsToMeasurement") <<
658  "Adding MC sample " << hname.Data() << " to channel " <<
659  channelname.Data() << " with " << htemp->Integral() << " events " <<
660  htemp->GetEntries();
661 
662  TString samplename = hname + TString("_sample");
663  RooStats::HistFactory::Sample sample(samplename.Data());
664  sample.SetNormalizeByTheory(true);
665 
666  // Check to enable statistical uncertainty
668  //mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Enable statistical uncertainty";
669  sample.ActivateStatError();
670  TH1* staterrorhisto = protoana::ProtoDUNEFitUtils::GetStatsSystHistogram(htemp);
671  RooStats::HistFactory::StatError& staterror = sample.GetStatError();
672  staterror.SetUseHisto();
673  staterror.SetErrorHist(staterrorhisto);
674  }
675 
677  if (_UseComputedSysts) {
678  ApplySystematicToSample(sample, htemp, systvec, false, false);
679  }
680  else {
681  BuildBackgroundSystThenApplyToSample(sample, htemp, false, false,
683  //BuildSystThenApplyToSample(sample, htemp, false, false,
684  // kBackground, _bkg_chan_index[j], _bkg_topo_index[j]);
685  }
686  }
687 
688  // Set histogram for sample
689  sample.SetHisto(htemp);
690 
691  if (_enable_bkg_factor[j]) {
692  if (htemp->Integral() > 0) {
693  TString poiname = hname;
694  poiname.ReplaceAll("MC","POI");
695  poiname.ReplaceAll("_Histo","");
696  poiname.ReplaceAll(".","");
697  poiname.ReplaceAll("-","_");
698  poiname.ReplaceAll(channelname.Data(),"");
699  poiname.ReplaceAll("__","_");
700  meas.SetPOI(poiname.Data());
701  //sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100./*2.0*/);
702  if (_RandSigPriors) {
703  sample.AddNormFactor(poiname.Data(), rand.Gaus(1.0, .1), 0., 100.);
704  }
705  else {
706  sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
707  }
708  mf::LogInfo("BuildMeasurement") << "Sample " << sample.GetName() << " has normalisation parameter " << poiname.Data();
709  }
710  }
711 
712  // Add sample to channel
713  channel.AddSample(sample);
714  }
715  }
716 
717  // Add signal samples to channel
718  for(unsigned int j=0; j < _sighistos.size(); j++){
719  // Clone histogram
720  TH1D* htemp = (TH1D*)(_sighistos.at(j)->Clone());
721 
722  TString hname(htemp->GetName());
723  //if(htemp->GetEntries() == 0 || htemp->Integral() == 0) continue;
724  if(hname.Contains(channelname.Data()) && hname.Contains("MC")){
725  mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Adding MC sample" << hname.Data() << " to channel " << channelname.Data() << " with " << htemp->Integral() << " events";
726  TString samplename = hname + TString("_sample");
727  RooStats::HistFactory::Sample sample(samplename.Data());
728  sample.SetNormalizeByTheory(true);
729 
730  // Check to enable statistical uncertainty
732  //mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Enable statistical uncertainty";
733  sample.ActivateStatError();
734  TH1* staterrorhisto = protoana::ProtoDUNEFitUtils::GetStatsSystHistogram(htemp);
735  RooStats::HistFactory::StatError& staterror = sample.GetStatError();
736  staterror.SetUseHisto();
737  staterror.SetErrorHist(staterrorhisto);
738  }
739 
741  if (_UseComputedSysts) { //Make this a vector?
742  ApplySystematicToSample(sample, htemp, systvec, true, _NormalisedSystematic);
743  }
744  else {
745  BuildSignalSystThenApplyToSample(sample, htemp, true,
747  _sig_truth_index[j]);
748  }
749  }
750 
751  // Set histogram for sample
752  sample.SetHisto(htemp);
753 
754  // Add POI to this sample
755  TString poiname = hname;
756  poiname.ReplaceAll("MC","POI");
757  poiname.ReplaceAll("_Histo","");
758  poiname.ReplaceAll(".0","");
759  poiname.ReplaceAll("-","_");
760  poiname.ReplaceAll(channelname.Data(),"");
761  poiname.ReplaceAll("__","_");
762  poiname.ReplaceAll("_0_1000000","");
763  meas.SetPOI(poiname.Data()); // AddPOI would also work
764 
765  if (_RandSigPriors) {
766  sample.AddNormFactor(poiname.Data(), rand.Gaus(1.0, .1), 0., 100.);
767  }
768  else {
769  sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
770  }
771  mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Sample " << sample.GetName() << " has normalisation parameter " << poiname.Data();
772 
773  // Add sample to channel
774  channel.AddSample(sample);
775  }
776  }
777 
778  // Statistical uncertainty less than 1% is ignored
779  channel.SetStatErrorConfig(_IgnoreStatisticalErrorBelow,"Poisson"); // Poisson or Gaussian
780 
781  // Add channel to measurement
782  mf::LogInfo("AddSamplesAndChannelsToMeasurement") << "Adding channel " << channel.GetName() << " to measurement " << meas.GetName();
783  meas.AddChannel(channel);
784  }
785 
786 }
std::vector< std::string > _SystFileNames
Definition: ProtoDUNEFit.h:109
bool BuildSignalSystThenApplyToSample(RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iChan, size_t iTopo, size_t iTruthBin)
bool ApplySystematicToSample(RooStats::HistFactory::Sample &sample, TH1 *histo, std::vector< TH1 * > systvec, bool hasnormfactor, bool isnorm)
TH1 * GetStatsSystHistogram(TH1 *nominal)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< size_t > _sig_truth_index
Definition: ProtoDUNEFit.h:127
uint8_t channel
Definition: CRTFragment.hh:201
std::vector< size_t > _sig_topo_index
Definition: ProtoDUNEFit.h:127
std::vector< int > _enable_bkg_factor
Definition: ProtoDUNEFit.h:147
std::vector< TH1 * > _bkghistos
Definition: ProtoDUNEFit.h:124
std::vector< size_t > _sig_chan_index
Definition: ProtoDUNEFit.h:126
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
Fw2dFFT::Data Data
double _IgnoreStatisticalErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< TH1 * > _datahistos
Definition: ProtoDUNEFit.h:124
std::vector< size_t > _bkg_chan_index
Definition: ProtoDUNEFit.h:126
bool BuildBackgroundSystThenApplyToSample(RooStats::HistFactory::Sample &sample, TH1 *histo, bool hasnormfactor, bool isnorm, size_t iChan, size_t iTopo)
std::vector< size_t > _bkg_topo_index
Definition: ProtoDUNEFit.h:127
std::vector< TH1 * > _sighistos
Definition: ProtoDUNEFit.h:124
QTextStream & endl(QTextStream &s)
std::vector< TH1 * > GetSystHistograms(std::string name)
void protoana::ProtoDUNEFit::AddSidebandSamplesAndChannelsToMeasurement ( RooStats::HistFactory::Measurement &  meas)
private

Definition at line 961 of file ProtoDUNEFit.cxx.

962  {
963 //********************************************************************
964 
965  std::string message_source = "AddSidebandSamplesAndChannelsToMeasurement";
966 
967  //Fix this
968  for (size_t i = 0; i < _MCControlSampleFiles.size(); ++i) { // remove this loop -- not needed
969  TString toponame = Form("Topo%s", _SidebandTopologyName[i].c_str());
970  //Replace Topology Name here with 'Sideband Name'? maybe APA2
971  TString channelname = Form("SidebandChannelAPA2");
972  //TString channelname = Form("SidebandChannel%s",
973  // _SidebandTopologyName[i].c_str());
974  RooStats::HistFactory::Channel channel(channelname.Data());
975 
976  // Add data to channel
978  for (size_t j = 0; j < _sideband_hists_data.size(); ++j) {
979  TH1D* htemp = (TH1D*)_sideband_hists_data[j]->Clone();
980  TString hname(htemp->GetName());
981 
982  if(hname.Contains(channelname.Data()) && hname.Contains("Data")){
983  mf::LogInfo(message_source) <<
984  "Adding sideband dataset " << hname.Data() << " to channel " <<
985  channelname.Data() << " with " << htemp->Integral() << " events";
986  data.SetHisto(htemp);
987  channel.SetData(data);
988  }
989  }
990 
991  // Add bkg samples to channel
992  for (size_t j = 0; j < _sideband_hists_mc.size(); ++j) {
993  TH1D* htemp = (TH1D*)_sideband_hists_mc[j]->Clone();
994  mf::LogInfo(message_source) << "sideband " << j << " " <<
995  htemp->GetName() << " " <<
996  htemp->GetEntries() << " " <<
997  htemp->Integral();
998  TString hname(htemp->GetName());
999  if (/*htemp->GetEntries() == 0 || */htemp->Integral() == 0) continue;
1000  if(hname.Contains(channelname.Data()) && hname.Contains("MC")){
1001  mf::LogInfo(message_source) << "Adding MC sideband sample " <<
1002  hname.Data() << " to channel " << channelname.Data() <<
1003  " with " << htemp->Integral() << " events";
1004 
1005  TString samplename = hname + TString("_sample");
1006  RooStats::HistFactory::Sample sample(samplename.Data());
1007  sample.SetNormalizeByTheory(true);
1008 
1009  // Check to enable statistical uncertainty
1011  //mf::LogInfo(message_source) << "Enable statistical uncertainty";
1012  sample.ActivateStatError();
1013  TH1* staterrorhisto =
1015  RooStats::HistFactory::StatError& staterror = sample.GetStatError();
1016  staterror.SetUseHisto();
1017  staterror.SetErrorHist(staterrorhisto);
1018  }
1019 
1020  //if(_EnableSystematicError){
1021  //ApplySystematicToSample(sample, htemp, systvec, false, false);
1022  //}
1023 
1024  // Set histogram for sample
1025  sample.SetHisto(htemp);
1026 
1027  // Add bkg normalisation parameter to this sample
1028  //if(hname.Contains(toponame.Data())){
1029  if(j == 0) {
1030  TString poiname = hname;
1031  poiname.ReplaceAll("MC","POI");
1032  poiname.ReplaceAll("_Histo","");
1033  poiname.ReplaceAll("_Hist","");
1034  poiname.ReplaceAll(channelname.Data(),"");
1035  poiname.ReplaceAll("SidebandChannel","");
1036  poiname.ReplaceAll(".","");
1037  poiname.ReplaceAll("-","_");
1038  poiname.ReplaceAll("__","_");
1039  meas.SetPOI(poiname.Data());
1040  //sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
1041  if (_RandSigPriors) {
1042  sample.AddNormFactor(poiname.Data(), rand.Gaus(1.0, .1), 0., 100.);
1043  }
1044  else {
1045  sample.AddNormFactor(poiname.Data(), 1.0, 0.0, 100.0);
1046  }
1047  mf::LogInfo(message_source) << "Sample " << sample.GetName() <<
1048  " has normalisation parameter " << poiname.Data();
1049  }
1050 
1051  // Add sample to channel
1052  channel.AddSample(sample);
1053  }
1054  }
1055 
1056  // Statistical uncertainty less than 1% is ignored
1057  // Poisson or Gaussian
1058  channel.SetStatErrorConfig(_IgnoreStatisticalErrorBelow,"Poisson");
1059 
1060  // Add channel to measurement
1061  mf::LogInfo(message_source) << "Adding sideband channel " <<
1062  channel.GetName() << " to measurement " << meas.GetName();
1063  meas.AddChannel(channel);
1064  }
1065 
1066 }
std::vector< std::string > _SidebandTopologyName
Definition: ProtoDUNEFit.h:109
TH1 * GetStatsSystHistogram(TH1 *nominal)
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
uint8_t channel
Definition: CRTFragment.hh:201
std::vector< TH1 * > _sideband_hists_mc
Definition: ProtoDUNEFit.h:130
Fw2dFFT::Data Data
double _IgnoreStatisticalErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< TH1 * > _sideband_hists_data
Definition: ProtoDUNEFit.h:130
std::vector< std::string > _MCControlSampleFiles
Definition: ProtoDUNEFit.h:109
bool protoana::ProtoDUNEFit::ApplyBuiltSystToSample ( TH1 *  histo,
TH1 *  high_hist,
TH1 *  low_hist,
RooStats::HistFactory::Sample &  sample,
std::string  syst_name,
std::string  syst_type,
bool  hasnormfactor 
)

Definition at line 2167 of file ProtoDUNEFit.cxx.

2171  {
2172 
2173  double nominal_integral = histo->Integral();
2174  double high_integral = high_hist->Integral();
2175  double high_val = high_integral/nominal_integral;
2176  double low_integral = low_hist->Integral();
2177  double low_val = low_integral/nominal_integral;
2178  std::string hname = histo->GetName();
2179 
2180  if ((fabs(1. - low_val) < _IgnoreSystematicErrorBelow) ||
2181  (fabs(high_val - 1.) < _IgnoreSystematicErrorBelow)) {
2182  std::cout << "WARNING::Ignoring systematic uncertainty " << syst_name <<
2183  " for sample " << hname << " as it is below the threshold " <<
2184  _IgnoreSystematicErrorBelow << " , syst fraction = " <<
2185  1. - low_val << " , " << high_val - 1. << std::endl;
2186  return false;
2187  }
2188 
2189  if (nominal_integral < 1.e-4 ||
2190  high_integral < 1.e-4 ||
2191  low_integral < 1.e-4 ) {
2192  return false;
2193  }
2194 
2195  //need this above
2196  if (syst_type == "NormOnly") {
2197  return false;
2198  }
2199 
2200  // Mixed (norm and/or shape) systematic case
2201  if(!hasnormfactor){ // Sample without normalisation parameter
2202 
2203  // Single bin samples and without norm factor are considered as norm syst
2205  sample.AddOverallSys(syst_name.c_str(), low_val, high_val);//systname from above
2206 
2207  std::cout << "INFO::Adding single bin and without norm parameter " <<
2208  "norm systematic " << syst_name << " , to sample " << hname <<
2209  " , with high value = " << high_val <<
2210  " , and low value = " << low_val <<
2211  " , with nominal integral = " << nominal_integral <<
2212  " , high histogram integral = " << high_integral <<
2213  " , low histogram integral = " << low_integral <<
2214  " , and syst fraction = " << 1. - low_val <<
2215  " , " << high_val - 1. << std::endl;
2216  return true;
2217  }
2218  else{ // multi-bin region and without norm parameter
2219  RooStats::HistFactory::HistoSys syst;
2220  syst.SetName(syst_name.c_str());
2221  syst.SetHistoHigh(high_hist);
2222  syst.SetHistoLow(low_hist);
2223  sample.AddHistoSys(syst);
2224 
2225  std::cout << "INFO::Adding multi-bin and without norm parameter " <<
2226  "shape+norm systematic " << syst_name <<
2227  " , to sample " << hname <<
2228  " , with nominal integral = " << nominal_integral <<
2229  " , high histogram integral = " << high_integral <<
2230  " , low histogram integral = " << low_integral <<
2231  ", and syst fraction = " << 1. - low_val <<
2232  " , " << high_val - 1. << std::endl;
2233  return true;
2234  }
2235  }
2236 
2237  return true;
2238 }
std::string string
Definition: nybbler.cc:12
const double e
double _IgnoreSystematicErrorBelow
Definition: ProtoDUNEFit.h:119
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::ApplySystematicToSample ( RooStats::HistFactory::Sample &  sample,
TH1 *  histo,
std::vector< TH1 * >  systvec,
bool  hasnormfactor,
bool  isnorm 
)

Definition at line 1526 of file ProtoDUNEFit.cxx.

1526  {
1527 //********************************************************************
1528 
1529  if(!histo){
1530  std::cerr << "ERROR::No input histogram found! Will not apply systematics!" << std::endl;
1531  return false;
1532  }
1533 
1534  if(systvec.empty()){
1535  std::cout << "INFO::No detector systematic vector found! Will not apply systematics!" << std::endl;
1536  return false;
1537  }
1538 
1539  TString hname(histo->GetName());
1540  hname.ReplaceAll("_Histo","");
1541 
1542  if(_SystToConsider.size() != _SystType.size()){
1543  std::cerr << "Vectors SystToConsider and SystType do not have the same size. Will not apply systematics. Check your fcl file." << std::endl;
1544  return false;
1545  }
1546 
1547  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
1548  TString systname(_SystToConsider[i].c_str());
1549  TString systtype(_SystType[i].c_str());
1550 
1551  TH1* highhist = NULL; TH1* lowhist = NULL;
1552  for (size_t j = 0; j < systvec.size(); ++j) {
1553  TH1* systhisto = (TH1*)systvec[j];
1554  TString systhistoname(systhisto->GetName());
1555 
1556  if (!systhistoname.Contains(hname.Data())) {
1557  continue;
1558  }
1559  if (!systhistoname.Contains(systname.Data())) {
1560  continue;
1561  }
1562 
1563  if(systhistoname.Contains("LOW") || systhistoname.Contains("Low") || systhistoname.Contains("low")){
1564  lowhist = systhisto;
1565  }
1566  if(systhistoname.Contains("HIGH") || systhistoname.Contains("High") || systhistoname.Contains("high")){
1567  highhist = systhisto;
1568  }
1569  }
1570 
1571  if(!highhist || !lowhist){
1572  std::cerr << "ERROR::Stage1: Systematic histograms not found! " <<
1573  "Will not apply systematic " << systname.Data() <<
1574  " to histogram " << histo->GetName() << ". Will skip!" <<
1575  std::endl;
1576  continue;
1577  }
1578 
1580  histo, highhist, "UP");
1582  histo, lowhist, "DOWN");
1583 
1584  if(!highsyst || !lowsyst){
1585  std::cerr << "ERROR::Stage2: Systematic histograms not found! " <<
1586  "Will not apply systematic " << systname.Data() <<
1587  " to histogram " << histo->GetName() << ". Will skip!" << std::endl;
1588  continue;
1589  }
1590 
1591  Double_t highval = highsyst->Integral()/histo->Integral();
1592  Double_t lowval = lowsyst->Integral()/histo->Integral();
1593 
1594  if (fabs(1. - lowval) < _IgnoreSystematicErrorBelow ||
1595  fabs(1. - highval) < _IgnoreSystematicErrorBelow) {
1596 
1597  std::cout << "WARNING::Ignoring systematic uncertainty " << systname
1598  << " for sample " << hname << " as it is below the threshold "
1599  << _IgnoreSystematicErrorBelow << " , syst fraction = "
1600  << 1. - lowval << " , " << highval - 1. << std::endl;
1601 
1602  continue;
1603  }
1604 
1605  // Case for normalisation systematic only
1606  if(systtype == "NormOnly"){
1607  if(isnorm && hname.Contains("Signal")){
1608  if(highval > 0. && lowval > 0.){
1609  std::cout << "INFO::Normalised systematic " << systname.Data() << " is considered." << std::endl;
1610  highsyst->Scale(1./highval);
1611  lowsyst->Scale(1./lowval);
1612  }
1613  }
1614  sample.AddOverallSys(systname.Data(), lowval, highval);
1615  std::cout << "INFO::Adding norm systematic by default with name " << systname << " , to sample " << hname << " , with high value = " << highval << " , and low value = " << lowval << " , with nominal integral = " << histo->Integral() << " , high histogram integral = " << highsyst->Integral() << " , low histogram integral = " << lowsyst->Integral() << " , and syst fraction = " << 1.-lowsyst->Integral()/histo->Integral() << " , " << highsyst->Integral()/histo->Integral()-1. << std::endl;
1616  continue;
1617  }
1618 
1619  // Mixed (norm and/or shape) systematic case
1620  if(!hasnormfactor){ // Sample without normalisation parameter
1621  if(protoana::ProtoDUNEFitUtils::IsSingleBinHisto(histo)){ // Single bin samples and without norm factor are considered as norm syst
1622  sample.AddOverallSys(systname.Data(), lowval, highval);
1623  std::cout << "INFO::Adding single bin and without norm parameter norm systematic " << systname << " , to sample " << hname << " , with high value = " << highval << " , and low value = " << lowval << " , with nominal integral = " << histo->Integral() << " , high histogram integral = " << highsyst->Integral() << " , low histogram integral = " << lowsyst->Integral() << " , and syst fraction = " << 1.-lowsyst->Integral()/histo->Integral() << " , " << highsyst->Integral()/histo->Integral()-1. << std::endl;
1624  }
1625  else{ // multi-bin region and without norm parameter
1626  RooStats::HistFactory::HistoSys syst;
1627  syst.SetName(systname.Data());
1628  syst.SetHistoHigh(highsyst);
1629  syst.SetHistoLow(lowsyst);
1630  sample.AddHistoSys(syst);
1631  std::cout << "INFO::Adding multi-bin and without norm parameter shape+norm systematic " << systname << " , to sample " << hname << " , with nominal integral = " << histo->Integral() << " , high histogram integral = " << highsyst->Integral() << " , low histogram integral = " << lowsyst->Integral() << ", and syst fraction = " << 1.-lowsyst->Integral()/histo->Integral() << " , " << highsyst->Integral()/histo->Integral()-1. << std::endl;
1632  }
1633  }
1634  else{ // Sample with normalisation parameter
1635  // Apply normalised systematic uncertainty
1636  if(isnorm && hname.Contains("Signal")){
1637  if(highval > 0. && lowval > 0.){
1638  std::cout << "INFO::Normalised systematic " << systname.Data() << " is considered." << std::endl;
1639  highsyst->Scale(1./highval);
1640  lowsyst->Scale(1./lowval);
1641  }
1642  }
1643 
1644  if(protoana::ProtoDUNEFitUtils::IsSingleBinHisto(histo)){ // Single bin samples and with norm factor are considered as norm syst
1645  sample.AddOverallSys(systname.Data(), lowval, highval);
1646  std::cout << "INFO::Adding norm systematic " << systname << " , to sample " << hname << " , with high value = " << highval << " , and low value = " << lowval << " , with nominal integral = " << histo->Integral() << " , high histogram integral = " << highsyst->Integral() << " , low histogram integral = " << lowsyst->Integral() << " , and syst fraction = " << 1.-lowsyst->Integral()/histo->Integral() << " , " << highsyst->Integral()/histo->Integral()-1. << std::endl;
1647  }
1648  else{ // multi-bin region and with norm parameter
1649  RooStats::HistFactory::HistoSys syst;
1650  syst.SetName(systname.Data());
1651  syst.SetHistoHigh(highsyst);
1652  syst.SetHistoLow(lowsyst);
1653  sample.AddHistoSys(syst);
1654  std::cout << "INFO::Adding shape+norm systematic " << systname << " , to sample " << hname << " , with nominal integral = " << histo->Integral() << " , high histogram integral = " << highsyst->Integral() << " , low histogram integral = " << lowsyst->Integral() << ", and syst fraction = " << 1.-lowsyst->Integral()/histo->Integral() << " , " << highsyst->Integral()/histo->Integral()-1. << std::endl;
1655  }
1656  }
1657  }
1658 
1659  return true;
1660 
1661 }
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
TH1 * GetSystematicHistoFromNominal(TH1 *nominal, TH1 *syst, TString name)
double _IgnoreSystematicErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::BuildBackgroundSystThenApplyToSample ( RooStats::HistFactory::Sample &  sample,
TH1 *  histo,
bool  hasnormfactor,
bool  isnorm,
size_t  iChan,
size_t  iTopo 
)

Definition at line 1811 of file ProtoDUNEFit.cxx.

1814  {
1815 
1816  TH1 * low_hist = 0x0;
1817  TH1 * high_hist = 0x0;
1818 
1819  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
1820  std::string syst_name = _SystToConsider[i];
1821  std::string syst_type = _SystType[i];
1822 
1823  double tmin = _TruthBinning[0];
1824  double tmax = _TruthBinning.back();
1825  if(_FitInReco){
1826  tmin = 0.0;
1827  tmax = 1000000.0;
1828  }
1829 
1830  low_hist =
1833  _ChannelNames[iChan], _BackgroundTopologyName[iTopo],
1834  _BackgroundTopology[iTopo], _EndZCut, tmin, tmax, _DoNegativeReco, -1,
1835  syst_name);
1836 
1837  high_hist =
1840  _ChannelNames[iChan], _BackgroundTopologyName[iTopo],
1841  _BackgroundTopology[iTopo], _EndZCut, tmin, tmax, _DoNegativeReco, 1,
1842  syst_name);
1843 
1844  if (!(low_hist && high_hist)) {
1845  continue;
1846  }
1847 
1848  std::cout << "Adding systs: " << low_hist->GetName() << " " <<
1849  high_hist->GetName() << std::endl;
1850  _syst_hists.push_back((TH1*)low_hist->Clone());
1851  _syst_hists.push_back((TH1*)high_hist->Clone());
1852 
1853  ApplyBuiltSystToSample(histo, high_hist, low_hist, sample, syst_name,
1854  syst_type, hasnormfactor);
1855  }
1856 
1857  return true;
1858 }
TH1 * FillMCBackgroundHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double endZ_cut, double minval=0.0, double maxval=100000.0, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::string string
Definition: nybbler.cc:12
std::vector< std::string > _BackgroundTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< int > _BackgroundTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
QTextStream & endl(QTextStream &s)
bool ApplyBuiltSystToSample(TH1 *histo, TH1 *high_hist, TH1 *low_hist, RooStats::HistFactory::Sample &sample, std::string syst_name, std::string syst_type, bool hasnormfactor)
std::vector< std::vector< std::pair< TH1 *, TH1 * > > > protoana::ProtoDUNEFit::BuildIncidentSignalSyst ( size_t  iTopo)

Definition at line 1918 of file ProtoDUNEFit.cxx.

1918  {
1919 
1920  //std::vector<TH1 *> low_hists, high_hists;
1921  std::vector<std::vector<std::pair<TH1*, TH1*>>> results;
1922 
1923  TH1 * low_hist = 0x0;
1924  TH1 * high_hist = 0x0;
1925 
1926  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
1927  std::string syst_name = _SystToConsider[i];
1928  std::string syst_type = _SystType[i];
1929 
1930  low_hist =
1933  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1935 
1936  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1937  low_hist->Add(
1940  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1942  }
1943 
1944 
1945  high_hist =
1948  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1950 
1951  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1952  high_hist->Add(
1955  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1957  }
1958 
1959  if (!(low_hist && high_hist)) {
1960  return results;
1961  }
1962 
1963  //Create a new vector for this syst
1964  results.push_back(std::vector<std::pair<TH1*, TH1*>>());
1965 
1966  //Split up the histos
1967  for (int j = 1; j <= low_hist->GetNbinsX(); j++) {
1968  TString low_name = Form("%s_IncBin%i", low_hist->GetName(),j);
1969  TString high_name = Form("%s_IncBin%i", high_hist->GetName(),j);
1970 
1971  TH1* temp_low_hist = (TH1*)low_hist->Clone();
1972  temp_low_hist->Reset();
1973  temp_low_hist->SetName(low_name.Data());
1974 
1975  TH1* temp_high_hist = (TH1*)high_hist->Clone();
1976  temp_high_hist->Reset();
1977  temp_high_hist->SetName(high_name.Data());
1978 
1979  for (int k = 1; k <= low_hist->GetNbinsX(); k++) {
1980  if (k == j) {
1981  temp_low_hist->SetBinContent(k, low_hist->GetBinContent(k));
1982  temp_high_hist->SetBinContent(k, high_hist->GetBinContent(k));
1983  }
1984  else {
1985  temp_low_hist->SetBinContent(k, 0.0);
1986  temp_high_hist->SetBinContent(k, 0.0);
1987  }
1988  }
1989  results.back().push_back({temp_high_hist, temp_low_hist});
1990  //low_hists.push_back(temp_low_hist);
1991  //high_hists.push_back(temp_high_hist);
1992  std::cout << "Adding systs: " << temp_low_hist->GetName() << " " <<
1993  temp_high_hist->GetName() << std::endl;
1994  _syst_hists.push_back((TH1*)temp_low_hist->Clone());
1995  _syst_hists.push_back((TH1*)temp_high_hist->Clone());
1996  }
1997  }
1998 
1999  std::cout << "Systs: " << _SystToConsider.size() << " " << results.size() << std::endl;
2000 
2001  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
2002  std::cout << i << " " << _SystToConsider[i] << " " << results[i].size() << std::endl;
2003  for (size_t j = 0; j < results[i].size(); ++j) {
2004  std::cout << results[i][j].first << " " << results[i][j].second << std::endl;
2005  }
2006  }
2007 
2008  return results;
2009 
2010 }
TH1 * FillMCIncidentHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string topo, int toponum, double reco_beam_endZ_cut, double minval=0., double maxval=100000., bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::string string
Definition: nybbler.cc:12
struct vector vector
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< int > _IncidentTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::BuildIncidentSystThenApplyToSample ( RooStats::HistFactory::Sample &  sample,
TH1 *  histo,
bool  hasnormfactor,
bool  isnorm,
size_t  iTopo 
)

Definition at line 1860 of file ProtoDUNEFit.cxx.

1862  {
1863 
1864  TH1 * low_hist = 0x0;
1865  TH1 * high_hist = 0x0;
1866 
1867  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
1868  std::string syst_name = _SystToConsider[i];
1869  std::string syst_type = _SystType[i];
1870 
1871  low_hist =
1874  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1876 
1877  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1878  low_hist->Add(
1881  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1883  }
1884 
1885 
1886  high_hist =
1889  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1891 
1892  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1893  high_hist->Add(
1896  _IncidentTopologyName[iTopo], _IncidentTopology[iTopo], 0., 10000.,
1898  }
1899 
1900  if (!(low_hist && high_hist)) {
1901  continue;
1902  }
1903 
1904  std::cout << "Adding systs: " << low_hist->GetName() << " " <<
1905  high_hist->GetName() << std::endl;
1906  _syst_hists.push_back((TH1*)low_hist->Clone());
1907  _syst_hists.push_back((TH1*)high_hist->Clone());
1908 
1909  ApplyBuiltSystToSample(histo, high_hist, low_hist, sample, syst_name,
1910  syst_type, hasnormfactor);
1911  }
1912 
1913  return true;
1914 }
TH1 * FillMCIncidentHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string topo, int toponum, double reco_beam_endZ_cut, double minval=0., double maxval=100000., bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::string string
Definition: nybbler.cc:12
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< int > _IncidentTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
QTextStream & endl(QTextStream &s)
bool ApplyBuiltSystToSample(TH1 *histo, TH1 *high_hist, TH1 *low_hist, RooStats::HistFactory::Sample &sample, std::string syst_name, std::string syst_type, bool hasnormfactor)
bool protoana::ProtoDUNEFit::BuildSignalSystThenApplyToSample ( RooStats::HistFactory::Sample &  sample,
TH1 *  histo,
bool  hasnormfactor,
bool  isnorm,
size_t  iChan,
size_t  iTopo,
size_t  iTruthBin 
)

Definition at line 1762 of file ProtoDUNEFit.cxx.

1765  {
1766 
1767  TH1 * low_hist = 0x0;
1768  TH1 * high_hist = 0x0;
1769 
1770  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
1771  std::string syst_name = _SystToConsider[i];
1772  std::string syst_type = _SystType[i];
1773 
1774  if (iTruthBin > _TruthBinning.size() - 1) {
1775  std::cout << "Error! Requesting truth bin " << iTruthBin <<
1776  " from binning vector of size " <<
1777  _TruthBinning.size() << std::endl;
1778  return false;
1779  }
1780 
1781  low_hist =
1784  _ChannelNames[iChan], _SignalTopologyName[iTopo],
1785  _SignalTopology[iTopo], _TruthBinning[iTruthBin-1],
1786  _TruthBinning[iTruthBin], _EndZCut, _DoNegativeReco, -1, syst_name);
1787 
1788  high_hist =
1791  _ChannelNames[iChan], _SignalTopologyName[iTopo],
1792  _SignalTopology[iTopo], _TruthBinning[iTruthBin-1],
1793  _TruthBinning[iTruthBin], _EndZCut, _DoNegativeReco, 1, syst_name);
1794 
1795  if (!(low_hist && high_hist)) {
1796  continue;
1797  }
1798 
1799  std::cout << "Adding systs: " << low_hist->GetName() << " " <<
1800  high_hist->GetName() << std::endl;
1801  _syst_hists.push_back((TH1*)low_hist->Clone());
1802  _syst_hists.push_back((TH1*)high_hist->Clone());
1803 
1804  ApplyBuiltSystToSample(histo, high_hist, low_hist, sample, syst_name,
1805  syst_type, hasnormfactor);
1806  }
1807 
1808  return true;
1809 }
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _SignalTopology
Definition: ProtoDUNEFit.h:107
TH1 * FillMCSignalHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double minval, double maxval, double endZ_cut, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
QTextStream & endl(QTextStream &s)
bool ApplyBuiltSystToSample(TH1 *histo, TH1 *high_hist, TH1 *low_hist, RooStats::HistFactory::Sample &sample, std::string syst_name, std::string syst_type, bool hasnormfactor)
bool protoana::ProtoDUNEFit::BuildSystThenApplyToSample ( RooStats::HistFactory::Sample &  sample,
TH1 *  histo,
bool  hasnormfactor,
bool  isnorm,
protoana::HistType  this_histType,
size_t  iChan,
size_t  iTopo,
size_t  iTruthBin = 999 
)

Definition at line 2012 of file ProtoDUNEFit.cxx.

2015  {
2016 
2017  TH1 * low_hist = 0x0;
2018  TH1 * high_hist = 0x0;
2019 
2020  for (size_t i = 0; i < _SystToConsider.size(); ++i) {
2021  std::string syst_name = _SystToConsider[i];
2022  std::string syst_type = _SystType[i];
2023 
2024  switch (this_histType) {
2025  case kSignal: {
2026 
2027  if (iTruthBin > _TruthBinning.size() - 1) {
2028  std::cout << "Error! Requesting truth bin " << iTruthBin <<
2029  " from binning vector of size " <<
2030  _TruthBinning.size() << std::endl;
2031  return false;
2032  }
2033 
2034  low_hist =
2037  _ChannelNames[iChan], _SignalTopologyName[iTopo],
2038  _SignalTopology[iTopo], _TruthBinning[iTruthBin-1],
2039  _TruthBinning[iTruthBin], _EndZCut, _DoNegativeReco, -1, syst_name);
2040 
2041  high_hist =
2044  _ChannelNames[iChan], _SignalTopologyName[iTopo],
2045  _SignalTopology[iTopo], _TruthBinning[iTruthBin-1],
2046  _TruthBinning[iTruthBin], _EndZCut, _DoNegativeReco, 1, syst_name);
2047 
2048  break;
2049  }
2050  case kBackground: {
2051 
2052  double tmin = _TruthBinning[0];
2053  double tmax = _TruthBinning.back();
2054  if(_FitInReco){
2055  tmin = 0.0;
2056  tmax = 1000000.0;
2057  }
2058 
2059  low_hist =
2062  _ChannelNames[iChan], _BackgroundTopologyName[iTopo],
2063  _BackgroundTopology[iTopo], _EndZCut, tmin, tmax, _DoNegativeReco, -1,
2064  syst_name);
2065 
2066  high_hist =
2069  _ChannelNames[iChan], _BackgroundTopologyName[iTopo],
2070  _BackgroundTopology[iTopo], _EndZCut, tmin, tmax, _DoNegativeReco, 1,
2071  syst_name);
2072  break;
2073  }
2074  case kIncident: {
2075  break;
2076  /*low_hist =
2077  protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions(
2078  _IncidentMCFileNames[0], _RecoTreeName, _RecoBinning,
2079  _ChannelNames[0], _IncidentTopologyName[iTopo], _IncidentTopology[iTopo],
2080  _EndZCut);*/
2081  }
2082  default: {
2083  return false;
2084  }
2085  }
2086 
2087  if (!(low_hist && high_hist)) {
2088  continue;
2089  }
2090 
2091  //ApplyBuiltSystToSample(histo, high_hist, low_hist, sample, syst_name,
2092  //syst_type, hasnormfactor);
2093 
2094  double nominal_integral = histo->Integral();
2095  double high_integral = high_hist->Integral();
2096  double high_val = high_integral/nominal_integral;
2097  double low_integral = low_hist->Integral();
2098  double low_val = low_integral/nominal_integral;
2099  std::string hname = histo->GetName();
2100 
2101  if ((fabs(1. - low_val) < _IgnoreSystematicErrorBelow) ||
2102  (fabs(high_val - 1.) < _IgnoreSystematicErrorBelow)) {
2103  std::cout << "WARNING::Ignoring systematic uncertainty " << syst_name <<
2104  " for sample " << hname << " as it is below the threshold " <<
2105  _IgnoreSystematicErrorBelow << " , syst fraction = " <<
2106  fabs(1. - low_val) << " , " << fabs(high_val - 1.) << std::endl;
2107  continue;
2108  }
2109 
2110  if (nominal_integral < 1.e-4 ||
2111  high_integral < 1.e-4 ||
2112  low_integral < 1.e-4 ) {
2113  continue;
2114  }
2115 
2116  _syst_hists.push_back((TH1*)low_hist->Clone());
2117  _syst_hists.push_back((TH1*)high_hist->Clone());
2118 
2119 
2120  //need this above
2121  if (syst_type == "NormOnly") {
2122  continue;
2123  }
2124 
2125  // Mixed (norm and/or shape) systematic case
2126  if(!hasnormfactor){ // Sample without normalisation parameter
2127 
2128  // Single bin samples and without norm factor are considered as norm syst
2130  sample.AddOverallSys(syst_name.c_str(), low_val, high_val);//systname from above
2131 
2132  std::cout << "INFO::Adding single bin and without norm parameter " <<
2133  "norm systematic " << syst_name << " , to sample " << hname <<
2134  " , with high value = " << high_val <<
2135  " , and low value = " << low_val <<
2136  " , with nominal integral = " << nominal_integral <<
2137  " , high histogram integral = " << high_integral <<
2138  " , low histogram integral = " << low_integral <<
2139  " , and syst fraction = " << 1. - low_val <<
2140  " , " << high_val - 1. << std::endl;
2141  continue;
2142  }
2143  else{ // multi-bin region and without norm parameter
2144  RooStats::HistFactory::HistoSys syst;
2145  syst.SetName(syst_name.c_str());
2146  syst.SetHistoHigh(high_hist);
2147  syst.SetHistoLow(low_hist);
2148  sample.AddHistoSys(syst);
2149 
2150  std::cout << "INFO::Adding multi-bin and without norm parameter " <<
2151  "shape+norm systematic " << syst_name <<
2152  " , to sample " << hname <<
2153  " , with nominal integral = " << nominal_integral <<
2154  " , high histogram integral = " << high_integral <<
2155  " , low histogram integral = " << low_integral <<
2156  ", and syst fraction = " << 1. - low_val <<
2157  " , " << high_val - 1. << std::endl;
2158  continue;
2159  }
2160  }
2161 
2162  }
2163 
2164  return true;
2165 }
TH1 * FillMCBackgroundHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double endZ_cut, double minval=0.0, double maxval=100000.0, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
std::vector< std::string > _BackgroundTopologyName
Definition: ProtoDUNEFit.h:109
const double e
std::vector< int > _BackgroundTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
double _IgnoreSystematicErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _SignalTopology
Definition: ProtoDUNEFit.h:107
TH1 * FillMCSignalHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double minval, double maxval, double endZ_cut, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
QTextStream & endl(QTextStream &s)
void protoana::ProtoDUNEFit::BuildWorkspace ( TString  Outputfile,
int  analysis = -1 
)

Definition at line 66 of file ProtoDUNEFit.cxx.

66  {
67  //********************************************************************
68 
69  bool hfilled = false;
70  bool hfilled_sidebands = false;
71  // Pion analysis
72  if(analysis == 1) {
73  if (_DoScaleMCToData) {
75  }
76  if (_DoScaleMuonContent) {
78  }
79 
80  hfilled = FillHistogramVectors_Pions();
82  hfilled_sidebands =
84  if (!hfilled_sidebands) return;
85  }
86  }
87 
88  if (!hfilled) return;
89 
90  if (_OnlyDrawXSecs) {
91  mf::LogInfo("BuildWorkspace") << "Only drawing XSecs. Exiting";
92 
93  TFile *f = new TFile(Outputfile.Data(), "RECREATE");
94  TDirectory *HistoDir = f->mkdir("OriginalHistograms");
95  HistoDir->cd();
96 
97  std::vector<TH1 *> pre_xsecs = DrawXSecs();
98  for (size_t k = 0; k < pre_xsecs.size(); ++k) {
99  pre_xsecs[k]->Write();
100  }
101 
102  std::vector<TH2 *> pre_smear = DrawSmearingMatrix(0x0);
103  for (size_t k = 0; k < pre_smear.size(); ++k) {
104  pre_smear[k]->Write();
105  }
106 
107  for(unsigned int k = 0; k < _sighistos.size(); k++){
108  _sighistos[k]->Write();
109  }
110  for(unsigned int k = 0; k < _bkghistos.size(); k++){
111  _bkghistos[k]->Write();
112  }
113  for(unsigned int k = 0; k < _datahistos.size(); k++){
114  _datahistos[k]->Write();
115  }
116  for(unsigned int k = 0; k < _truthsighistos.size(); k++){
117  _truthsighistos[k]->Write();
118  }
119 
120  // Decorate and save efficiency graphs
121  for(unsigned int i=0; i < _efficiencyGraphs.size(); i++){
122  TGraphAsymmErrors* effgraph = (TGraphAsymmErrors*)_efficiencyGraphs.at(i);
123 
124  TCanvas* ceffgraph = new TCanvas(Form("ceffgraph%i",i),Form("Efficiency for channel %i",i));
125 
126  DecorateEfficiency( effgraph );
127  effgraph->Draw("*a");
128 
129  ceffgraph->Write();
130  }
131 
132  for(unsigned int k = 0; k < _incsighistos.size(); k++){
133  _incsighistos[k]->Write();
134  }
135  for(unsigned int k = 0; k < _incbkghistos.size(); k++){
136  _incbkghistos[k]->Write();
137  }
138  for(unsigned int k = 0; k < _incdatahistos.size(); k++){
139  _incdatahistos[k]->Write();
140  }
141 
142  DecorateEfficiency(_incidentEfficiency, "E_{true} at slice [MeV]");
143  _incidentEfficiency->Write();
144  _incidentEfficiencyNum->Write();
145  _incidentEfficiencyDenom->Write();
146  ////////////////////////////
147 
148  //Interacting efficiency
149  for( size_t i = 0; i < _interactingEfficiencyDenoms.size(); ++i ){
150  _interactingEfficiencyDenoms[i]->Write();
151  _interactingEfficiencyNums[i]->Write();
152 
154  _interactingEfficiencies[i]->Write();
155 
156  }
157 
158 
159  f->Close();
160 
161  return;
162  }
163 
164  // Get pion flux
165  //TH1* pionflux_ereco_histo = protoana::ProtoDUNESelectionUtils::FillMCFlux_Pions(_MCFileNames[0], _TruthTreeName, _TruthBinning, 1);
166  //TH1* pionflux_etruth_histo = protoana::ProtoDUNESelectionUtils::FillMCFlux_Pions(_MCFileNames[0], _TruthTreeName, _TruthBinning, 2);
167  //TH1* pionflux_eff_histo = (TH1*)pionflux_ereco_histo->Clone();
168  //pionflux_eff_histo->Divide(pionflux_etruth_histo);
169  //pionflux_eff_histo->SetNameTitle("pionflux_eff_histo","Efficiency for incident particles");
170  //TH1* pionrecoflux_ereco_histo = protoana::ProtoDUNESelectionUtils::FillMCFlux_Pions(_MCFileNames[0], _TruthTreeName, _RecoBinning, 3);
171  //TH1* pionrecoflux_eff_histo = (TH1*)pionrecoflux_ereco_histo->Clone();
172  //pionrecoflux_eff_histo->Divide(pionflux_etruth_histo);
173  //pionrecoflux_eff_histo->SetNameTitle("pionrecoflux_eff_histo","Efficiency for incident particles");
174 
175  // Create measurement object
176  RooStats::HistFactory::Measurement meas("ProtoDUNEFitExample","ProtoDUNE fit example");
177 
178  // Since this is not LHC, don't care about luminosity. Set to constant
179  meas.SetLumi(1.0);
180  meas.SetLumiRelErr(0.001);
181  meas.AddConstantParam("Lumi");
182 
183  // SetExportOnly to false meaning that we will fit the model
184  meas.SetExportOnly(false);
185 
186  // Add samples and channels to measurement
190 
191  // Sidebands
194  }
195 
196 
197  // Print info about meas
198  meas.PrintTree();
199 
200  // Export meas to workspace
201  RooStats::HistFactory::HistoToWorkspaceFactoryFast h2w(meas);
202  RooWorkspace* ws = h2w.MakeCombinedModel(meas);
203 
204  // Define interpolation code: 6th order polynomial interpolation and linear extrapolation
206 
207  // Save initial snapshot
208  protoana::ProtoDUNEFitUtils::SaveSnapshot(ws, Form("%s_initial_snapshot",ws->GetName()));
209 
210  // Fit and/or generate data using the workspace
212 
213  // Fit options - must be defined first
214  fitandgen->SetFitStrategy(_FitStrategy);
215  fitandgen->SetMinimiser(_Minimizer);
217  fitandgen->EnableMinosError();
218 
219  RooFitResult *fitresult = NULL;
220  TTree *toys_tree = NULL;
221  TString treename("protodUNE_dataresults");
222 
223  // Plots before fit
224  std::vector<TString> truebinsnameVec;
225  std::vector<TString> sidebandNames;
226 
227  for(size_t l = 0; l < _BackgroundTopologyName.size(); l++){
228  TString str = Form("%s", _BackgroundTopologyName[l].c_str());
229  truebinsnameVec.push_back(str);
230  }
231 
232  if(!_FitInReco){
233  for(size_t l = 1; l < _TruthBinning.size(); l++){
234  TString str = Form("Signal %.1f-%.1f", _TruthBinning[l-1], _TruthBinning[l]);
235  truebinsnameVec.push_back(str);
236  }
237  }
238 
239  for (size_t l = 0; l < _SidebandTopologyName.size(); ++l) {
240  TString str = Form("%s", _SidebandTopologyName[l].c_str());
241  sidebandNames.push_back(str);
242  }
243 
244  std::vector< TString > incidentNameVec;
245  for (size_t l = 1; l < _IncidentTopologyName.size(); ++l) {
246  TString str = Form("%s", _IncidentTopologyName[l].c_str());
247  std::cout << "NameVec " << l << " " << str << std::endl;
248  incidentNameVec.push_back(str);
249  }
250  for (size_t l = 1; l < _TruthBinning.size(); ++l) {
251  incidentNameVec.push_back(Form("%s_%.1f-%.1f",
253  std::cout << incidentNameVec.back() << std::endl;
254  }
255 
256  std::vector<TCanvas*> bfplots =
258  ws, "beforefit", "Poisson", "ratio", truebinsnameVec, _RecoBinning,
259  incidentNameVec, sidebandNames, _SidebandBinning, "Before Fit", _DoNegativeReco);
260 
261  RooAbsData *asimovdata = ws->data("asimovData");
262  std::vector<TCanvas*> bfAsimovplots =
264  ws, "asimov", "Poisson", "ratio", truebinsnameVec, _RecoBinning,
265  incidentNameVec, sidebandNames, _SidebandBinning, "Asimov Dataset", _DoNegativeReco, asimovdata);
266 
267  // ----------------------------------------------------------------------------------------------------
268  // Check if this is MC toys case
269  // ----------------------------------------------------------------------------------------------------
270 
271  if (_NToys > 1/*0*/) {
272  toys_tree = fitandgen->GenerateAndFit(ws, _NToys);
273  std::cout << "Got tree " << toys_tree << std::endl;
274  toys_tree->SetNameTitle("protodUNE_mctoysresults", "protodUNE_mctoysresults");
275 
276  // For each set of plots need to clone the output tree to avoid ROOT crashes
277  TTree *mctoys_results1 = (TTree*)toys_tree->Clone();
278  TCanvas* nuisancecanvas = protoana::ProtoDUNEFitUtils::PlotNuisanceParameters(mctoys_results1, ws);
279 
280  TTree *mctoys_results2 = (TTree*)toys_tree->Clone();
281  TCanvas* pullscanvas = protoana::ProtoDUNEFitUtils::PlotParametersPull(mctoys_results2, ws);
282 
283  TTree *mctoys_results3 = (TTree*)toys_tree->Clone();
284  TCanvas* nuispullscanvas = protoana::ProtoDUNEFitUtils::PlotNuisanceParametersPull(mctoys_results3, ws);
285 
286  TTree *mctoys_results4 = (TTree*)toys_tree->Clone();
287  //TCanvas* avresultcanvas = protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys(mctoys_results4, ws, "Channel0", "SigTopo");
288  TCanvas* avresultcanvas = protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys(mctoys_results4, ws, "POI", "POI");
289  std::cout << "plotted ave" << std::endl;
290 
291  // Save workspace
293 
294  TFile *f = new TFile(Outputfile.Data(), "UPDATE");
295 
296  // Save all histograms used in the measurement
297  meas.writeToFile(f);
298 
299  toys_tree->Write();
300  nuisancecanvas->Write();
301  if (_NToys > 1) {
302  pullscanvas->Write();
303  nuispullscanvas->Write();
304  }
305  avresultcanvas->Write();
306 
307  for(unsigned int i=0; i < bfAsimovplots.size(); i++){
308  bfAsimovplots[i]->Write();
309  }
310  for(unsigned int i=0; i < bfplots.size(); i++){
311  bfplots[i]->Write();
312  }
313 
314  // Histograms below are also saved from the measurement
315  /*
316  // Save all input histograms
317  TDirectory *HistoDir = f->mkdir("OriginalHistograms");
318  HistoDir->cd();
319 
320  for(unsigned int k = 0; k < _sighistos.size(); k++){
321  _sighistos[k]->Write();
322  }
323  for(unsigned int k = 0; k < _bkghistos.size(); k++){
324  _bkghistos[k]->Write();
325  }
326  for(unsigned int k = 0; k < _datahistos.size(); k++){
327  _datahistos[k]->Write();
328  }
329  for(unsigned int k = 0; k < _truthsighistos.size(); k++){
330  _truthsighistos[k]->Write();
331  }
332 
333  // Decorate and save efficiency graphs
334  for(unsigned int i=0; i < _efficiencyGraphs.size(); i++){
335  TGraphAsymmErrors* effgraph = (TGraphAsymmErrors*)_efficiencyGraphs.at(i);
336 
337  TCanvas* ceffgraph = new TCanvas(Form("ceffgraph%i",i),Form("Efficiency for channel %i",i));
338 
339  TAxis *ax = effgraph->GetHistogram()->GetXaxis();
340  Double_t x1 = ax->GetBinLowEdge(1);
341  //Double_t x2 = ax->GetBinUpEdge(ax->GetNbins());
342  effgraph->GetHistogram()->GetXaxis()->Set(_TruthBinning.size(),x1,_TruthBinning.size());
343  for(unsigned int i = 1; i < _TruthBinning.size(); i++){
344  TString ibinstr = Form("%.1f-%.1f",_TruthBinning[i-1],_TruthBinning[i]);
345  effgraph->GetHistogram()->GetXaxis()->SetBinLabel(i, ibinstr.Data());
346  }
347 
348  effgraph->Draw("*a");
349  effgraph->SetMarkerStyle(20);
350  effgraph->SetMarkerColor(1);
351  effgraph->SetTitle("Efficiency");
352  effgraph->GetXaxis()->SetTitle("E_{true} at vertex [MeV]");
353  effgraph->GetYaxis()->SetTitle("Efficiency");
354  effgraph->GetYaxis()->SetTitleOffset(1.25);
355 
356  ceffgraph->Write();
357  }
358 
359  HistoDir->cd("..");
360  */
361 
362  f->Close();
363 
364  // Nothing else to do here
365  return;
366  }
367 
368  // ----------------------------------------------------------------------------------------------------
369  // Continue with either asimov or data fit
370  // ----------------------------------------------------------------------------------------------------
371 
372  else if (_NToys == 1) {
373  fitresult = fitandgen->GenerateAndFitOneToy(ws);
374  }
375  else if (_DoAsimovFit) {
376  fitresult = fitandgen->FitAsimovData(ws);
377  treename = "protodUNE_asimovresults";
378  }
379  else{
380  fitresult = fitandgen->FitData(ws);
381  }
382 
383  // Try different fit configurations if the fit failed
384  // Reset all values and errors
385  //protoana::ProtoDUNEFitUtils::ResetAllValuesAndErrors(ws);
386 
387  if(fitresult->status() != 0){
388  std::cout << "WARNING::Fit failed - trying with Scan" << std::endl;
389  fitandgen->SetAlgorithm("Scan");
390  fitresult = fitandgen->FitData(ws);
391  }
392 
393  if(fitresult->status() != 0){
394  if(_FitStrategy == 0){
395  std::cout << "WARNING::Fit failed with fit strategy 0 - trying with fit strategy 1." << std::endl;
396  fitandgen->SetAlgorithm(ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
397  fitandgen->SetFitStrategy(1);
398  fitresult = fitandgen->FitData(ws);
399  }
400  }
401 
402  if(fitresult->status() != 0){
403  std::cout << "WARNING::Fit failed - trying with improve" << std::endl;
404  fitandgen->SetAlgorithm("migradimproved");
405  fitandgen->SetMinimiser("Minuit");
406  fitresult = fitandgen->FitData(ws);
407  }
408 
409  // Save results from data or asimov fit
410  if(!fitresult){
411  std::cerr << "No fit resuls found. Will exit!" << std::endl;
412  return;
413  }
414 
415  // Import fit result in workspace
416  ws->import(*fitresult,kTRUE);
417 
418  // Print fit results on the screen
419  //fitresult->Print();
420 
421  std::vector<TCanvas*> afplots =
423  ws, "afterfit", "Poisson", "ratio", truebinsnameVec, _RecoBinning,
424  incidentNameVec, sidebandNames, _SidebandBinning, "After fit", _DoNegativeReco,
425  NULL, fitresult);
426 
427 
428  // Save post-fit workspace snapshot
429  protoana::ProtoDUNEFitUtils::SaveSnapshot(ws, Form("%s_postfit_snapshot",ws->GetName()));
430  protoana::ProtoDUNEFitUtils::SaveSnapshot(ws, Form("%s_postfitForPlots_snapshot",ws->GetName()));
431 
432  TH2* fitcov = protoana::ProtoDUNEFitUtils::GetFitCovariance(fitresult);
434 
435  toys_tree = fitandgen->RooFitResultToTTree(ws, fitresult);
436  toys_tree->SetNameTitle(treename.Data(),treename.Data());
437 
438  // Plot NLL
439  std::vector<TCanvas*> nllplots = protoana::ProtoDUNEFitUtils::PlotNLL(
440  ws, "PDFit", fitresult);
441 
442  TTree *mctoys_results1 = (TTree*)toys_tree->Clone();
443  TCanvas* nuisancecanvas = protoana::ProtoDUNEFitUtils::PlotNuisanceParameters(mctoys_results1, ws);
444 
445  TTree *mctoys_results2 = (TTree*)toys_tree->Clone();
446 
447  RooArgList prefit_POI = fitresult->floatParsInit();
448  TCanvas* avresultcanvas = protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys(mctoys_results2, ws, "POI", "POI", &prefit_POI);
449 
450  //Testing
451  /*
452  RooArgList prefit_POI = fitresult->floatParsInit();
453  TIterator* itr = prefit_POI.createIterator();
454  RooRealVar * var = 0x0;
455  std::cout << "Prefit!!" << std::endl;
456  while ( (var = (RooRealVar*)itr->Next()) ) {
457  std::string name = var->GetName();
458  if (name.find("POI") == std::string::npos) continue;
459  std::cout << var->GetName() << " " << var->getVal() << std::endl;
460 
461  }
462  */
463 
464 
465  //TTree *mctoys_results3 = (TTree*)toys_tree->Clone();
466  //TCanvas* avresultcanvas_inc = protoana::ProtoDUNEFitUtils::PlotAverageResultsFromToys(mctoys_results3, ws, "POI", "POI");
467 
468  // Fit fixing nuisance parameters
469  //protoana::ProtoDUNEFitUtils::LoadSnapshot(ws, Form("%s_postfitForPlots_snapshot",ws->GetName()));
470  //protoana::ProtoDUNEFitUtils::MakeNuisanceParamsConstant(ws,"");
471  //RooFitResult *fitresultnosysts = fitandgen->FitData(ws);
472  //if(fitresultnosysts){
473  //ws->import(*fitresultnosysts,kTRUE);
474  //}
475  //protoana::ProtoDUNEFitUtils::SaveSnapshot(ws, Form("%s_postfitNosysts_snapshot",ws->GetName()));
476 
477  // Save workspace
479 
480  TFile *f = new TFile(Outputfile.Data(), "UPDATE");
481 
482  // Save all histograms used in the measurement
483  meas.writeToFile(f);
484 
485  fitcov->Write();
486  fitcor->Write();
487  toys_tree->Write();
488  if(nuisancecanvas)
489  nuisancecanvas->Write();
490  if(avresultcanvas)
491  avresultcanvas->Write();
492  //if(avresultcanvas_inc)
493  //avresultcanvas_inc->Write();
494 
495  for(unsigned int i=0; i < bfAsimovplots.size(); i++){
496  bfAsimovplots[i]->Write();
497  }
498  for(unsigned int i=0; i < bfplots.size(); i++){
499  bfplots[i]->Write();
500  }
501  for(unsigned int i=0; i < afplots.size(); i++){
502  afplots[i]->Write();
503  }
504 
505  TDirectory *NLLDir = f->mkdir("NLLPlots");
506  NLLDir->cd();
507  for (size_t i = 0; i < nllplots.size(); ++i) {
508  nllplots[i]->Write();
509  }
510  NLLDir->cd("..");
511 
512  //pionflux_etruth_histo->Write();
513  //pionflux_ereco_histo->Write();
514  //pionflux_eff_histo->Write();
515  //pionrecoflux_ereco_histo->Write();
516  //pionrecoflux_eff_histo->Write();
517 
518 
519  TDirectory *HistoDir = f->mkdir("OriginalHistograms");
520  HistoDir->cd();
521 
522  for(unsigned int k = 0; k < _sighistos.size(); k++){
523  _sighistos[k]->Write();
524  }
525  for(unsigned int k = 0; k < _bkghistos.size(); k++){
526  _bkghistos[k]->Write();
527  }
528  for(unsigned int k = 0; k < _datahistos.size(); k++){
529  _datahistos[k]->Write();
530  }
531  for(unsigned int k = 0; k < _truthsighistos.size(); k++){
532  _truthsighistos[k]->Write();
533  }
534 
535  if (_syst_hists.size()) {
536  TDirectory * SystDir = f->mkdir("SystematicHists");
537  SystDir->cd();
538 
539  for (size_t i = 0; i < _syst_hists.size(); ++i) {
540  _syst_hists[i]->Write();
541  }
542  HistoDir->cd();
543  }
544  // Decorate and save efficiency graphs
545  for(unsigned int i=0; i < _efficiencyGraphs.size(); i++){
546  TGraphAsymmErrors* effgraph = (TGraphAsymmErrors*)_efficiencyGraphs.at(i);
547 
548  TCanvas* ceffgraph = new TCanvas(Form("ceffgraph%i",i),Form("Efficiency for channel %i",i));
549 
550  DecorateEfficiency( effgraph );
551  effgraph->Draw("*a");
552 
553  ceffgraph->Write();
554  }
555 
556  // Try drawing the xsecs
557  if (_DoDrawXSecs) {
558  std::vector<TH1 *> pre_xsecs = DrawXSecs();
559  for (size_t k = 0; k < pre_xsecs.size(); ++k) {
560  pre_xsecs[k]->Write();
561  }
562 
563  std::vector<TH1 *> post_xsecs = DrawXSecs(fitresult);
564  for (size_t k = 0; k < post_xsecs.size(); ++k) {
565  post_xsecs[k]->Write();
566  }
567  std::cout << "Finished XSec" << std::endl;
568  }
569 
570  for(unsigned int k = 0; k < _incsighistos.size(); k++){
571  _incsighistos[k]->Write();
572  }
573  for(unsigned int k = 0; k < _incbkghistos.size(); k++){
574  _incbkghistos[k]->Write();
575  }
576  for(unsigned int k = 0; k < _incdatahistos.size(); k++){
577  _incdatahistos[k]->Write();
578  }
579 
580  std::cout << "Wrote incident" << std::endl;
581 
582  //Incident efficiency
583  TCanvas* ceffgraph = new TCanvas( "ceffgraph", "Efficiency" );
584 
585  _incidentEfficiency->Draw("a");
586  DecorateEfficiency(_incidentEfficiency, "E_{true} at slice [MeV]");
587 
588  ceffgraph->Write();
589 
590  _incidentEfficiency->Write();
591  _incidentEfficiencyNum->Write();
592  _incidentEfficiencyDenom->Write();
593  ////////////////////////////
594 
595  //Interacting efficiency
596  for( size_t i = 0; i < _interactingEfficiencyDenoms.size(); ++i ){
597  _interactingEfficiencyDenoms[i]->Write();
598  _interactingEfficiencyNums[i]->Write();
599 
601  _interactingEfficiencies[i]->Write();
602 
603  }
604 
605  HistoDir->cd("..");
606 
607  std::cout << "Closing" << std::endl;
608  f->Close();
609  std::cout << "Closed" << std::endl;
610 
611  return;
612 
613 }
std::vector< TH1 * > _interactingEfficiencyNums
Definition: ProtoDUNEFit.h:138
std::vector< TGraphAsymmErrors * > _efficiencyGraphs
Definition: ProtoDUNEFit.h:132
std::vector< TCanvas * > PlotDatasetsAndPdfs(RooWorkspace *work, TString name, TString error, TString plottodraw, std::vector< TString > binnames, std::vector< double > recobins, std::vector< TString > incidentBinNames, std::vector< TString > sidebandBinNames, std::vector< double > sidebandBins, TString measurement="PDFit", bool doNegativeReco=false, RooAbsData *data=NULL, RooFitResult *result=NULL)
void SetInterpolationCode(RooWorkspace *ws, int code)
RooFitResult * FitData(RooWorkspace *ws, bool isWeighted=false)
void SaveWorkspace(RooWorkspace *ws, TString outFileName)
std::vector< std::string > _SidebandTopologyName
Definition: ProtoDUNEFit.h:109
TGraphAsymmErrors * _incidentEfficiency
Definition: ProtoDUNEFit.h:134
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< TH1 * > _truthsighistos
Definition: ProtoDUNEFit.h:124
std::vector< std::string > _BackgroundTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _incsighistos
Definition: ProtoDUNEFit.h:129
std::vector< TH1 * > _incdatahistos
Definition: ProtoDUNEFit.h:129
std::vector< TH1 * > _bkghistos
Definition: ProtoDUNEFit.h:124
static QStrList * l
Definition: config.cpp:1044
void AddIncidentSamplesAndChannelsToMeasurement(RooStats::HistFactory::Measurement &meas)
std::vector< TH1 * > _interactingEfficiencyDenoms
Definition: ProtoDUNEFit.h:137
std::vector< TH1 * > _syst_hists
Definition: ProtoDUNEFit.h:125
void AddSidebandSamplesAndChannelsToMeasurement(RooStats::HistFactory::Measurement &meas)
void AddSamplesAndChannelsToMeasurement(RooStats::HistFactory::Measurement &meas)
RooFitResult * FitAsimovData(RooWorkspace *w)
std::vector< double > _SidebandBinning
Definition: ProtoDUNEFit.h:105
TH2 * GetFitCovariance(RooFitResult *result)
TCanvas * PlotParametersPull(TTree *tree, RooWorkspace *ws)
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< TH2 * > DrawSmearingMatrix(RooFitResult *fitresult, bool doPostFit=false)
std::vector< TGraphAsymmErrors * > _interactingEfficiencies
Definition: ProtoDUNEFit.h:139
TCanvas * PlotNuisanceParameters(TTree *tree, RooWorkspace *ws)
std::vector< TH1 * > DrawXSecs(RooFitResult *fitresult=0x0)
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< TCanvas * > PlotNLL(RooWorkspace *work, TString name, RooFitResult *result, bool plotPLL=false)
TH2 * GetFitCorrelationMatrix(RooFitResult *result)
void SaveSnapshot(RooWorkspace *ws, TString snapshotname)
std::vector< TH1 * > _datahistos
Definition: ProtoDUNEFit.h:124
std::vector< TH1 * > _incbkghistos
Definition: ProtoDUNEFit.h:129
TTree * GenerateAndFit(RooWorkspace *ws, int nexp)
TCanvas * PlotAverageResultsFromToys(TTree *tree, RooWorkspace *ws, TString channelname, TString catname, RooArgList *PreFit_POI=0x0)
TTree * RooFitResultToTTree(RooWorkspace *ws, RooFitResult *res)
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
void ScaleMCToData(bool data_is_mc=false)
TCanvas * PlotNuisanceParametersPull(TTree *tree, RooWorkspace *ws)
RooFitResult * GenerateAndFitOneToy(RooWorkspace *ws)
std::vector< TH1 * > _sighistos
Definition: ProtoDUNEFit.h:124
void DecorateEfficiency(TGraphAsymmErrors *eff, std::string x_title="E_{true} at vertex [MeV]")
static QCString str
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::Configure ( std::string  configPath)
private

Definition at line 1664 of file ProtoDUNEFit.cxx.

1664  {
1665  //********************************************************************
1666 
1667  cet::filepath_first_absolute_or_lookup_with_dot cetfilelook("FHICL_FILE_PATH");
1668 
1669  // Parse configuration file
1670  auto const pset = fhicl::ParameterSet::make(configPath, cetfilelook);
1671 
1672  _RecoTreeName = pset.get<std::string>("RecoTreeName");
1673  _Minimizer = pset.get<std::string>("Minimizer");
1674  _TruthTreeName = pset.get<std::string>("TruthTreeName");
1675 
1676  _DataFileNames = pset.get<std::vector<std::string>>("DataFileNames");
1677  _MCFileNames = pset.get<std::vector<std::string>>("MCFileNames");
1678  _MCControlSampleFiles = pset.get<std::vector<std::string>>("MCControlSampleFiles");
1679  _SidebandTopologyName = pset.get<std::vector<std::string>>("SidebandTopologyName");
1680  _SidebandTopology = pset.get<std::vector<int>>("SidebandTopology");
1681  //_SidebandBinning = pset.get<std::pair<double,double>>("SidebandBinning");
1682  _SidebandBinning = pset.get<std::vector<double>>("SidebandBinning");
1683  _DataControlSampleFiles = pset.get<std::vector<std::string>>("DataControlSampleFiles");
1684  _IncidentMCFileNames = pset.get<std::vector<std::string>>("IncidentMCFileNames");
1685  _IncidentDataFileNames = pset.get<std::vector<std::string>>("IncidentDataFileNames");
1686  _SystFileNames = pset.get<std::vector<std::string>>("SystFileNames");
1687  _SystToConsider = pset.get<std::vector<std::string>>("SystToConsider");
1688  _SystType = pset.get<std::vector<std::string>>("SystType");
1689  _BackgroundTopologyName = pset.get<std::vector<std::string>>("BackgroundTopologyName");
1690  _SignalTopologyName = pset.get<std::vector<std::string>>("SignalTopologyName");
1691  _IncidentTopologyName = pset.get<std::vector<std::string>>("IncidentTopologyName");
1692  _ChannelNames = pset.get<std::vector<std::string>>("ChannelNames");
1693 
1694  _FitStrategy = pset.get<int>("FitStrategy");
1695  _NToys = pset.get<int>("NToys");
1696 
1697  _IgnoreStatisticalErrorBelow = pset.get<double>("IgnoreStatisticalErrorBelow");
1698  _IgnoreSystematicErrorBelow = pset.get<double>("IgnoreSystematicErrorBelow");
1699 
1700  _EnableMinosError = pset.get<bool>("EnableMinosError");
1701  _DoAsimovFit = pset.get<bool>("DoAsimovFit");
1702  _EnableStatisticalError = pset.get<bool>("EnableStatisticalError");
1703  _EnableSystematicError = pset.get<bool>("EnableSystematicError");
1704  _UseComputedSysts = pset.get<bool>("UseComputedSysts");
1705  _NormalisedSystematic = pset.get<bool>("NormalisedSystematic");
1706 
1707  _RecoBinning = pset.get< std::vector<double> >("RecoBinning");
1708  _TruthBinning = pset.get< std::vector<double> >("TruthBinning");
1709  _FitInReco = pset.get<bool>("FitInReco");
1710 
1711  _SignalTopology = pset.get< std::vector<int> >("SignalTopology");
1712  _BackgroundTopology = pset.get< std::vector<int> >("BackgroundTopology");
1713  _IncidentTopology = pset.get< std::vector<int> >("IncidentTopology");
1714 
1715  _AddIncidentToMeasurement = pset.get<bool>("AddIncidentToMeasurement");
1716  _AddSidebandsToMeasurement = pset.get<bool>("AddSidebandsToMeasurement");
1717  _DoNegativeReco = pset.get<bool>("DoNegativeReco");
1718  _EndZCut = pset.get<double>("EndZCut");
1719  _AddBackgroundFactors = pset.get<std::vector<int>>("AddBackgroundFactors");
1720  _AddIncidentBackgroundFactors = pset.get<std::vector<int>>("AddIncidentBackgroundFactors");
1721  //_AddIncidentBackgroundFactors = pset.get<bool>("AddIncidentBackgroundFactors");
1722 
1723  _DoScaleMCToData = pset.get<bool>("DoScaleMCToData");
1724  _DoScaleMuonContent = pset.get<bool>("DoScaleMuonContent");
1725  if (_DoScaleMuonContent)
1726  _MuonScaleFactor = pset.get<double>("MuonScaleFactor");
1727 
1728  _RandSigPriors = pset.get<bool>("RandSigPriors");
1729  //_StatFluctuation = pset.get<bool>("StatFluctuation");
1730  _DataIsMC = pset.get<bool>("DataIsMC");
1731  _OnlyDrawXSecs = pset.get<bool>("OnlyDrawXSecs");
1732  _DoDrawXSecs = pset.get<bool>("DoDrawXSecs");
1733  _WirePitch = pset.get<double>("WirePitch");
1734 
1735  return true;
1736 
1737 }
std::vector< std::string > _SystFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _IncidentDataFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _DataFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _SidebandTopologyName
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
static ParameterSet make(intermediate_table const &tbl)
Definition: ParameterSet.cc:68
std::vector< int > _AddBackgroundFactors
Definition: ProtoDUNEFit.h:147
std::vector< std::string > _BackgroundTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< int > _BackgroundTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _SystToConsider
Definition: ProtoDUNEFit.h:109
std::vector< double > _SidebandBinning
Definition: ProtoDUNEFit.h:105
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
std::string _TruthTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _AddIncidentBackgroundFactors
Definition: ProtoDUNEFit.h:147
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< int > _IncidentTopology
Definition: ProtoDUNEFit.h:107
double _IgnoreSystematicErrorBelow
Definition: ProtoDUNEFit.h:119
double _IgnoreStatisticalErrorBelow
Definition: ProtoDUNEFit.h:119
std::vector< std::string > _SystType
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< int > _SidebandTopology
Definition: ProtoDUNEFit.h:116
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _SignalTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _DataControlSampleFiles
Definition: ProtoDUNEFit.h:109
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCControlSampleFiles
Definition: ProtoDUNEFit.h:109
void protoana::ProtoDUNEFit::DecorateEfficiency ( TGraphAsymmErrors *  eff,
std::string  x_title = "E_{true} at vertex [MeV]" 
)
private

Definition at line 1739 of file ProtoDUNEFit.cxx.

1740  {
1741  TAxis * ax = eff->GetHistogram()->GetXaxis();
1742  double x1 = ax->GetBinLowEdge(1);
1743  /*eff->GetHistogram()->GetXaxis()*/ax->Set((_TruthBinning.size() - 1),
1744  x1, (_TruthBinning.size() - 1));
1745 
1746  for(size_t i = 1; i < _TruthBinning.size(); ++i) {
1747  TString label = Form("%.1f-%.1f", _TruthBinning[i-1], _TruthBinning[i]);
1748  /*eff->GetHistogram()->GetXaxis()*/ax->SetBinLabel(i, label.Data());
1749  }
1750 
1751  //eff->Draw("*a");
1752  eff->SetMarkerStyle(20);
1753  eff->SetMarkerColor(1);
1754  eff->SetTitle("Efficiency");
1755  eff->GetXaxis()->SetTitle(x_title.c_str());
1756  eff->GetYaxis()->SetTitle("Efficiency");
1757  eff->GetYaxis()->SetTitleOffset(1.25);
1758 
1759 }
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< TH2 * > protoana::ProtoDUNEFit::DrawSmearingMatrix ( RooFitResult *  fitresult,
bool  doPostFit = false 
)

Definition at line 2384 of file ProtoDUNEFit.cxx.

2384  {
2385 
2386  std::vector<TH2 *> smearing_hists;
2387 
2388  std::map<std::string, std::vector<double>> POI_vals;
2389 
2390  if (doPostFit) {
2391  RooArgList floatParsList = fitresult->floatParsFinal();
2392  TIterator* itr = floatParsList.createIterator();
2393  RooRealVar * var = 0x0;
2394  while ( (var = (RooRealVar*)itr->Next()) ) {
2395  std::string name = var->GetName();
2396  if (name.find("POI") == std::string::npos) continue;
2397  std::cout << var->GetName() << " " << var->getVal() << std::endl;
2398 
2399  bool found = false;
2400  for (size_t i = 0; i < _SignalTopologyName.size(); ++i) {
2401  if (name.find(_SignalTopologyName[i]) != std::string::npos) {
2402  POI_vals[_SignalTopologyName[i]].push_back(var->getVal());
2403  found = true;
2404  break;
2405  }
2406  }
2407  if (!found)
2408  POI_vals["INC"].push_back(var->getVal());
2409  }
2410  }
2411 
2412  std::string incident_name = (doPostFit ? "IncidentSignalPostFit2D" :
2413  "IncidentSignalPreFit2D");
2414  TH2 * inc_signal_hist = new TH2D(incident_name.c_str(), "",
2415  _TruthBinning.size() - 1, 0,
2416  _TruthBinning.size() - 1,
2417  (_DoNegativeReco ? _RecoBinning.size() :
2418  _RecoBinning.size()-1),
2419  (_DoNegativeReco ? -1 : 0),
2420  _RecoBinning.size() - 1);
2421 
2422  for (size_t i = 0; i < _TruthBinning.size()-1; ++i) {
2423  for (int j = 1; j <= inc_signal_hist->GetNbinsY(); ++j) {
2424  //Scale by the post fit parameter value
2425  inc_signal_hist->SetBinContent(i+1, j,
2426  _incsighistos[i]->GetBinContent(j)*(doPostFit ? POI_vals["INC"][i] : 1.));
2427 
2428  if (i == 0) {
2429  if (_DoNegativeReco && j == 1) {
2430  inc_signal_hist->GetYaxis()->SetBinLabel(j, "< 0");
2431  }
2432  else if (_DoNegativeReco) {
2433  TString jbinstr = Form("%.1f-%.1f", _RecoBinning[j-2], _RecoBinning[j-1]);
2434  inc_signal_hist->GetYaxis()->SetBinLabel(j, jbinstr);
2435  }
2436  else {
2437  TString jbinstr = Form("%.1f-%.1f", _RecoBinning[j-1], _RecoBinning[j]);
2438  inc_signal_hist->GetYaxis()->SetBinLabel(j, jbinstr);
2439  }
2440 
2441  }
2442  }
2443  TString ibinstr = Form("%.1f-%.1f", _TruthBinning[i], _TruthBinning[i+1]);
2444  inc_signal_hist->GetXaxis()->SetBinLabel(i+1, ibinstr);
2445  inc_signal_hist->SetTitle(";E_{True} (MeV);E_{Reco} (MeV)");
2446  }
2447  smearing_hists.push_back(inc_signal_hist);
2448 
2449  std::map<std::string, std::vector<TH1 *>> signal_hists_by_topo;
2450  std::map<std::string, std::vector<TH1 *>> mixed_hists_by_topo;
2451  for (std::string topo : _SignalTopologyName) {
2452  for (size_t i = 0; i < _sighistos.size(); ++i) {
2453  std::string name = _sighistos[i]->GetName();
2454  if (name.find("_" + topo + "_") == std::string::npos) continue;
2455  if (name.find("MC_Channel" + topo) == std::string::npos) {
2456  mixed_hists_by_topo[topo].push_back(_sighistos[i]);
2457  }
2458  else {
2459  signal_hists_by_topo[topo].push_back(_sighistos[i]);
2460  }
2461  }
2462  }
2463 
2464 
2465  for (std::string topo : _SignalTopologyName) {
2466  std::string signal_name = "Signal" + topo +
2467  (doPostFit ? "PostFit2D" : "PreFit2D");
2468  int nBinsY = inc_signal_hist->GetNbinsY();
2469  TH2 * signal_hist =
2470  new TH2D(signal_name.c_str(), "",
2471  inc_signal_hist->GetNbinsX(), 0,
2472  inc_signal_hist->GetNbinsX(),
2473  nBinsY,
2474  inc_signal_hist->GetYaxis()->GetBinLowEdge(1),
2475  (inc_signal_hist->GetYaxis()->GetBinLowEdge(nBinsY) +
2476  inc_signal_hist->GetYaxis()->GetBinWidth(nBinsY)));
2477 
2478 
2479 
2480  for (size_t i = 0; i < signal_hists_by_topo[topo].size(); ++i) {
2481  TH1 * hist = signal_hists_by_topo[topo][i];
2482  double scale = (doPostFit ? POI_vals[topo][i] : 1.);
2483  for (int j = 1; j <= inc_signal_hist->GetNbinsY(); ++j) {
2484  signal_hist->SetBinContent(
2485  i+1, j, scale*hist->GetBinContent(j));
2486  if (i == 0) {
2487  if (_DoNegativeReco && j == 1) {
2488  signal_hist->GetYaxis()->SetBinLabel(j, "< 0");
2489  }
2490  else if (_DoNegativeReco) {
2491  TString jbinstr = Form("%.1f-%.1f", _RecoBinning[j-2], _RecoBinning[j-1]);
2492  signal_hist->GetYaxis()->SetBinLabel(j, jbinstr);
2493  }
2494  else {
2495  TString jbinstr = Form("%.1f-%.1f", _RecoBinning[j-1], _RecoBinning[j]);
2496  signal_hist->GetYaxis()->SetBinLabel(j, jbinstr);
2497  }
2498 
2499  }
2500  }
2501  TString ibinstr = Form("%.1f-%.1f", _TruthBinning[i], _TruthBinning[i+1]);
2502  signal_hist->GetXaxis()->SetBinLabel(i+1, ibinstr);
2503  signal_hist->SetTitle(";E_{True} (MeV);E_{Reco} (MeV)");
2504  }
2505  smearing_hists.push_back(signal_hist);
2506  }
2507 
2508  return smearing_hists;
2509 }
static QCString name
Definition: declinfo.cpp:673
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
std::vector< TH1 * > _incsighistos
Definition: ProtoDUNEFit.h:129
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
int var
Definition: 018_def.c:9
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
std::vector< TH1 * > _sighistos
Definition: ProtoDUNEFit.h:124
QTextStream & endl(QTextStream &s)
std::vector< TH1 * > protoana::ProtoDUNEFit::DrawXSecs ( RooFitResult *  fitresult = 0x0)

Definition at line 2240 of file ProtoDUNEFit.cxx.

2240  {
2241 
2242  std::cout << "Drawing XSec" << std::endl;
2243 
2244  std::vector<TH1 *> xsecs;
2245 
2246  std::map<std::string, std::vector<double>> POI_vals;
2247 
2248  if (fitresult) {
2249  RooArgList floatParsList = fitresult->floatParsFinal();
2250  TIterator* itr = floatParsList.createIterator();
2251  RooRealVar * var = 0x0;
2252  while ( (var = (RooRealVar*)itr->Next()) ) {
2253  std::string name = var->GetName();
2254  if (name.find("POI") == std::string::npos) continue;
2255  std::cout << var->GetName() << " " << var->getVal() << std::endl;
2256 
2257  bool found = false;
2258  for (size_t i = 0; i < _SignalTopologyName.size(); ++i) {
2259  if (name.find(_SignalTopologyName[i]) != std::string::npos) {
2260  POI_vals[_SignalTopologyName[i]].push_back(var->getVal());
2261  found = true;
2262  break;
2263  }
2264  }
2265  if (!found)
2266  POI_vals["INC"].push_back(var->getVal());
2267  }
2268  }
2269 
2270  //Make the incident hists
2271 
2272  std::string incident_name = (fitresult ? "IncidentSignalPostFit" : "IncidentSignalPreFit");
2273  TH1 * inc_signal_hist = new TH1D(incident_name.c_str(), "",
2274  _TruthBinning.size()-1, 0,
2275  _TruthBinning.size()-1);
2276  //TH1 * inc_signal_hist_prefit = new TH1D("IncidentSignalPreFit", "",
2277  // _TruthBinning.size()-1, 0,
2278  // _TruthBinning.size()-1);
2279  for (size_t i = 0; i < _TruthBinning.size()-1; ++i) {
2280  std::cout << "Making bin for Truth Bin: " << _TruthBinning[i] << " " <<
2281  _TruthBinning[i+1] << std::endl;
2282  std::cout << "Inc hist: " << _incsighistos[i]->GetName() << std::endl;
2283 
2284  //Scale by the post fit parameter value
2285  inc_signal_hist->SetBinContent(i+1,
2286  _incsighistos[i]->Integral()*(fitresult ? POI_vals["INC"][i] : 1.));
2287  //inc_signal_hist_prefit->SetBinContent(i+1, _incsighistos[i]->Integral());
2288  }
2289 
2290  //xsecs.push_back((TH1*)inc_signal_hist_prefit->Clone("IncidentSignalPreFitNoEff"));
2291 
2292  //Efficiency correction
2293  for (int i = 0; i < _incidentEfficiency->GetN(); ++i) {
2294  if (_incidentEfficiency->GetY()[i] < 1.e-5) continue;
2295  inc_signal_hist->SetBinContent(i+1,
2296  inc_signal_hist->GetBinContent(i+1) / _incidentEfficiency->GetY()[i]);
2297 
2298  //inc_signal_hist_prefit->SetBinContent(i+1,
2299  // inc_signal_hist_prefit->GetBinContent(i+1) / _incidentEfficiency->GetY()[i]);
2300  }
2301 
2302  xsecs.push_back(inc_signal_hist);
2303 
2304  std::map<std::string, std::vector<TH1 *>> signal_hists_by_topo;
2305  std::map<std::string, std::vector<TH1 *>> mixed_hists_by_topo;
2306  for (std::string topo : _SignalTopologyName) {
2307  for (size_t i = 0; i < _sighistos.size(); ++i) {
2308  std::string name = _sighistos[i]->GetName();
2309  if (name.find("_" + topo + "_") == std::string::npos) continue;
2310  if (name.find("MC_Channel" + topo) == std::string::npos) {
2311  mixed_hists_by_topo[topo].push_back(_sighistos[i]);
2312  }
2313  else {
2314  signal_hists_by_topo[topo].push_back(_sighistos[i]);
2315  }
2316  }
2317  }
2318 
2319 
2320  for (std::string topo : _SignalTopologyName) {
2321  std::string signal_name = "Signal" + topo +
2322  (fitresult ? "PostFit" : "PreFit");
2323  TH1 * signal_hist = new TH1D(signal_name.c_str(), "",
2324  inc_signal_hist->GetNbinsX(), 0,
2325  inc_signal_hist->GetNbinsX());
2326  //TH1 * signal_hist_prefit = new TH1D(("PreFitSignal" + topo).c_str(), "",
2327  // inc_signal_hist->GetNbinsX(), 0,
2328  // inc_signal_hist->GetNbinsX());
2329 
2330  TGraphAsymmErrors * eff = 0x0;
2331  for (size_t i = 0; i < _interactingEfficiencies.size(); ++i) {
2332  std::string name = _interactingEfficiencies[i]->GetName();
2333  if (name.find(topo) != std::string::npos) {
2334  eff = _interactingEfficiencies[i];
2335  break;
2336  }
2337  }
2338 
2339  for (size_t i = 0; i < signal_hists_by_topo[topo].size(); ++i) {
2340  //correct by eff
2341  double from_signal = signal_hists_by_topo[topo][i]->Integral();
2342  from_signal = from_signal / (eff->GetY()[i] > 1.e-5 ? eff->GetY()[i] : 1.);
2343 
2344  //double from_mixed = mixed_hists_by_topo[topo][i]->Integral();
2345 
2346  double combined = (from_signal /*+ from_mixed*/) /*/ eff->GetY()[i]*/;
2347 
2348  signal_hist->SetBinContent(i+1, (fitresult ? POI_vals[topo][i] : 1.)*combined);
2349  //signal_hist_prefit->SetBinContent(i+1, combined);
2350  }
2351 
2352  xsecs.push_back(signal_hist);
2353  //xsecs.push_back(signal_hist_prefit);
2354 
2355 
2356  //Now divide the interacting hist by the incident hist
2357  std::string xsec_name = "XSEC_" + topo +
2358  (fitresult ? "_PostFit" : "_PreFit");
2359  TH1 * xsec = (TH1*)signal_hist->Clone(xsec_name.c_str());
2360  xsec->Sumw2();
2361  inc_signal_hist->Sumw2();
2362  xsec->Divide(inc_signal_hist);
2363  xsec->Scale(1.E27/ (_WirePitch * 1.4 * 6.022E23 / 39.948 ));
2364 
2365  for (size_t i = 1; i < _TruthBinning.size(); ++i) {
2366  TString ibinstr = Form("%.1f-%.1f",_TruthBinning[i-1],_TruthBinning[i]);
2367  xsec->GetXaxis()->SetBinLabel(i, ibinstr);
2368  }
2369 
2370  //TCanvas * xsec_canvas = new TCanvas(("Canvas_XSec_"+chan).c_str(), "xsec", 500, 400);
2371  xsec->SetMinimum(0.);
2372  xsec->SetMarkerColor(1);
2373  xsec->SetLineColor(1);
2374  xsec->SetMarkerStyle(20);
2375  xsec->SetMarkerSize(.75);
2376  xsec->SetLineWidth(2);
2377  xsecs.push_back(xsec);
2378  }
2379 
2380  return xsecs;
2381 
2382 }
static QCString name
Definition: declinfo.cpp:673
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
TGraphAsymmErrors * _incidentEfficiency
Definition: ProtoDUNEFit.h:134
std::string string
Definition: nybbler.cc:12
std::vector< TH1 * > _incsighistos
Definition: ProtoDUNEFit.h:129
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< TGraphAsymmErrors * > _interactingEfficiencies
Definition: ProtoDUNEFit.h:139
int var
Definition: 018_def.c:9
std::vector< TH1 * > _sighistos
Definition: ProtoDUNEFit.h:124
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::FillHistogramVectors_Pions ( )
private

Definition at line 1069 of file ProtoDUNEFit.cxx.

1069  {
1070 //********************************************************************
1071 
1072  const int nmcchannels = _MCFileNames.size();
1073  const int ndatachannels = _DataFileNames.size();
1074  const int nchannelnames = _ChannelNames.size();
1075  const int ninctoponames = _IncidentTopologyName.size();
1076 
1077  const int nbkgtopo = _BackgroundTopology.size();
1078  const int nbkgtoponames = _BackgroundTopologyName.size();
1079  const int nsigtopo = _SignalTopology.size();
1080  const int nsigtoponames = _SignalTopologyName.size();
1081  const int ninctopo = _IncidentTopology.size();
1082 
1083  if(nmcchannels != ndatachannels){
1084  mf::LogError("FillHistogramVectors_Pions") << "The number of data and MC channels is not the same. Check MCFileNames and DataFileNames in the fcl file. Time to die!";
1085  return false;
1086  }
1087 
1088  if(nmcchannels != nchannelnames || ndatachannels != nchannelnames){
1089  mf::LogError("FillHistogramVectors_Pions") << "The channel names do not correspond to the input files. Check MCFileNames, DataFileNames and ChannelNames in the fcl file. Time to die!";
1090  return false;
1091  }
1092 
1093  if(nbkgtopo != nbkgtoponames){
1094  mf::LogError("FillHistogramVectors_Pions") << "Background topologies and background names do not have the same length. Check BackgroundTopology and BackgroundTopologyName in the fcl file. Time to die!";
1095  return false;
1096  }
1097 
1098  if(nsigtopo != nsigtoponames){
1099  mf::LogError("FillHistogramVectors_Pions") << "Signal topologies and background name vectors do not have the same size. Check SignalTopology and SignalTopologyName in the fcl file. Time to die!";
1100  return false;
1101  }
1102 
1103  if(ninctopo != ninctoponames){
1104  mf::LogError("FillHistogramVectors_Pions") << "Incident topologies and name vectors do not have the same size. Check IncidentTopology and IncidentTopologyName in the fcl file. Time to die!";
1105  return false;
1106  }
1107 
1108  // Get total number of data and MC triggers
1109  //int nmctriggers = protoana::ProtoDUNESelectionUtils::GetNTriggers_Pions(_MCFileNames[0], _RecoTreeName);
1110  //int ndatatriggers = protoana::ProtoDUNESelectionUtils::GetNTriggers_Pions(_DataFileNames[0], _RecoTreeName, false);
1111  //double mcnorm = (double)ndatatriggers/nmctriggers;
1112  //mf::LogInfo("FillHistogramVectors_Pions") << "Total number of MC triggers = " << nmctriggers << ", total number of data triggers = " << ndatatriggers << " , data/MC = " << mcnorm;
1113 
1114  Double_t tmin = _TruthBinning[0];
1115  Double_t tmax = _TruthBinning[_TruthBinning.size()-1];
1116  if(_FitInReco){
1117  tmin = 0.0;
1118  tmax = 1000000.0;
1119  }
1120 
1121  for(int i=0; i < nmcchannels; i++){
1122  for(int j=0; j < nbkgtopo; j++){
1123  int topo = _BackgroundTopology[j];
1124  _bkghistos.push_back(
1127  _BackgroundTopologyName[j], topo, _EndZCut, tmin, tmax,
1128  _DoNegativeReco, 0, "", _ScaleFactor,
1131 
1132  //For systs later
1133  _bkg_chan_index.push_back(i);
1134  _bkg_topo_index.push_back(j);
1135 
1136  //_incbkghistos.push_back( protoana::ProtoDUNESelectionUtils::FillMCBackgroundHistogram_Pions(_MCFileNames[0], _RecoTreeName, _RecoBinning, i, topo, true) );
1137  }
1138 
1139  //TH1D* sigevenshisto = new TH1D(Form("sigevenshisto_channel%i",i),Form("sigevenshisto_channel%i",i),_TruthBinning.size()-1,0,_TruthBinning.size()-1);
1140  //TH1D* truevenshisto = new TH1D(Form("truevenshisto_channel%i",i),Form("truevenshisto_channel%i",i),_TruthBinning.size()-1,0,_TruthBinning.size()-1);
1141 
1142  for(int j=0; j < nsigtopo; j++){
1143  int topo = _SignalTopology[j];
1144  for(unsigned int k=1; k < _TruthBinning.size(); k++){
1145  if(_FitInReco){
1146  TH1* hsignal =
1149  _ChannelNames[i], _SignalTopologyName[j], topo, tmin, tmax,
1150  _EndZCut);
1151 
1152  // Make one histogram per bin
1153  for(int k=1; k <= hsignal->GetNbinsX(); k++){
1154  TString hname = Form("%s_RecoBin%i",hsignal->GetName(), k);
1155  TH1 *hsignal_clone = (TH1*)hsignal->Clone();
1156  hsignal_clone->SetName(hname.Data());
1157  hsignal_clone->SetBinContent(k, hsignal->GetBinContent(k));
1158  for(int kk=1; kk <= hsignal_clone->GetNbinsX(); kk++){
1159  if(kk == k) continue;
1160  hsignal_clone->SetBinContent(kk, 0.0);
1161  }
1162  _sighistos.push_back(hsignal_clone);
1163  }
1164  }
1165  else{
1166  _sighistos.push_back(
1169  _ChannelNames[i], _SignalTopologyName[j], topo,
1171  _DoNegativeReco, 0, "", _ScaleFactor,
1173 
1174  //For systs later
1175  _sig_chan_index.push_back(i);
1176  _sig_topo_index.push_back(j);
1177  _sig_truth_index.push_back(k);
1178 
1179  //sigevenshisto->SetBinContent(k, (_sighistos.back())->Integral() + sigevenshisto->GetBinContent(k));
1180  }
1181  }
1182  //_incsighistos.push_back( protoana::ProtoDUNESelectionUtils::FillMCSignalHistogram_Pions(_MCFileNames[0], _RecoTreeName, _RecoBinning, i, topo, 0, 0, true) );
1183  //TH1* inchistoall = protoana::ProtoDUNESelectionUtils::FillMCSignalHistogram_Pions(_MCFileNames[0], _RecoTreeName, _RecoBinning, i, topo, 0, 0, true);
1184  //for(Int_t k=1; k <= inchistoall->GetNbinsX(); k++){
1185  //TH1 *newhist = (TH1*)inchistoall->Clone();
1186  //newhist->SetNameTitle(Form("%s_RecoBin%i",inchistoall->GetName(),k), Form("%s in Reco Bin %i",inchistoall->GetName(),k));
1187  //for(Int_t kk=1; kk <= newhist->GetNbinsX(); kk++){
1188  // if(kk != k){
1189  // newhist->SetBinContent(kk,0);
1190  // newhist->SetBinError(kk,0);
1191  // }
1192  //}
1193  //_incsighistos.push_back(newhist);
1194  //}
1195  }
1196 
1197  //_truthsighistos.push_back( protoana::ProtoDUNESelectionUtils::FillMCTruthSignalHistogram_Pions( _MCFileNames[0], _TruthTreeName, _TruthBinning, i) );
1198  //truevenshisto->Add(_truthsighistos.back());
1199 
1200  //TGraphAsymmErrors* effgraph = new TGraphAsymmErrors(sigevenshisto,truevenshisto);
1201  //effgraph->SetNameTitle(Form("Efficiency_channel%i",i), Form("Efficiency for channel %i",i));
1202  //_efficiencyGraphs.push_back(effgraph);
1203 
1204  }
1205 
1206  for(int i=0; i < ndatachannels; i++){
1207  _datahistos.push_back(
1211  /*
1212  if (_StatFluctuation) {
1213  for (int j = 0; j <= _datahistos.back()->GetNbinsX(); ++j) {
1214  double new_val = rand.Poisson(_datahistos.back()->GetBinContent(j));
1215  _datahistos.back()->SetBinContent(j, new_val);
1216  }
1217  }
1218  */
1219  }
1220 
1221  for(int i=0; i < ninctopo; i++){
1222  if (i == 0) {
1223  for (size_t k = 1; k < _TruthBinning.size(); ++k) {
1224  TH1 * inchisto =
1231  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1232  inchisto->Add(
1239  }
1240  TString hname = Form("%s_TrueBin_%.1f-%.1f", inchisto->GetName(),
1241  _TruthBinning[k-1], _TruthBinning[k]);
1242  inchisto->SetName(hname.Data());
1243  _incsighistos.push_back(inchisto);
1244  _inc_sig_topo_index.push_back(i);
1245  }
1246  }
1247  else {
1248  TH1* inchisto =
1252  0., 10000., _DoNegativeReco, 0, "", _IncidentScaleFactor,
1254 
1255  for (size_t j = 1; j < _IncidentMCFileNames.size(); ++j) {
1256  inchisto->Add(
1260  0., 10000., _DoNegativeReco, 0, "", _IncidentScaleFactor,
1262  }
1263 
1264  inchisto->SetNameTitle(Form("MC_ChannelIncident_%s_Histo",
1265  _IncidentTopologyName[i].c_str()),
1266  Form("Incident MC for topology %s",
1267  _IncidentTopologyName[i].c_str()));
1268 
1269  _incbkghistos.push_back(inchisto);
1271  _inc_bkg_topo_index.push_back(i);
1272  }
1273  }
1274 
1275  TH1* incdatahisto =
1277  //_DataFileNames[0], _RecoTreeName, _RecoBinning, _ChannelNames[0],
1278  _IncidentDataFileNames[0], _RecoTreeName, _RecoBinning, ""/*_ChannelNames[0]*/,
1279  _EndZCut, _DoNegativeReco, true);
1280 
1281  for (size_t i = 1; i < _IncidentDataFileNames.size(); ++i) {
1282  incdatahisto->Add(
1284  //_DataFileNames[i], _RecoTreeName, _RecoBinning, _ChannelNames[i],
1285  _IncidentDataFileNames[i], _RecoTreeName, _RecoBinning, ""/*_ChannelNames[i]*/,
1286  _EndZCut, _DoNegativeReco, true));
1287  }
1288  _incdatahistos.push_back(incdatahisto);
1289  /*
1290  if (_StatFluctuation) {
1291  for (int j = 0; j <= _incdatahistos.back()->GetNbinsX(); ++j) {
1292  double new_val = rand.Poisson(_incdatahistos.back()->GetBinContent(j));
1293  _incdatahistos.back()->SetBinContent(j, new_val);
1294  }
1295  }
1296  */
1297 
1298  std::pair<TH1 *, TH1 *> inc_eff_num_denom =
1302 
1303  _incidentEfficiencyNum = inc_eff_num_denom.first;
1304  _incidentEfficiencyDenom = inc_eff_num_denom.second;
1305  for (int i = 1; i < _incidentEfficiencyDenom->GetNbinsX(); ++i) {
1306  if (_incidentEfficiencyDenom->GetBinContent(i) < 1.e-5)
1307  _incidentEfficiencyDenom->SetBinContent(i, 1.);
1308  }
1309 
1310  //_incidentEfficiency = new TGraphAsymmErrors(_incidentEfficiencyNum,
1311  // _incidentEfficiencyDenom);
1312 
1313  //Need to build up the hists for events the pass the selection
1314  TH1D * incident_numerator = new TH1D("incident_numerator", "",
1315  _TruthBinning.size()-1, 0,
1316  _TruthBinning.size()-1);
1317 
1318  for (size_t i = 0; i < _TruthBinning.size() - 1; ++i) {
1319  incident_numerator->SetBinContent(i+1, _incsighistos[i]->Integral());
1320  }
1321 
1322  //_incidentEfficiencyDenom->Scale(_ScaleFactor);
1323  _incidentEfficiency = new TGraphAsymmErrors(incident_numerator,
1325 
1326  _incidentEfficiency->SetNameTitle("MC_Incident_Efficiency", "Efficiency");
1327 
1328  std::map<std::string, std::vector<TH1 *>> signal_hists_by_topo;
1329  for (std::string topo : _SignalTopologyName) {
1330  for (size_t i = 0; i < _sighistos.size(); ++i) {
1331  std::string name = _sighistos[i]->GetName();
1332  //std::cout << name << std::endl;
1333  if (name.find("_" + topo + "_") == std::string::npos) continue;
1334  if (name.find("MC_Channel" + topo) != std::string::npos) {
1335  signal_hists_by_topo[topo].push_back(_sighistos[i]);
1336  }
1337  }
1338  }
1339 
1340  //Interacting Efficiencies
1341  for( int i = 0; i < nsigtopo; ++i ){
1342  int topo = _SignalTopology[i];
1343 
1344  std::pair< TH1 *, TH1 * > eff_num_denom =
1347  _ChannelNames[i], _SignalTopologyName[i], topo, _EndZCut,
1348  0, _ScaleFactor);
1349 
1350  _interactingEfficiencyNums.push_back(eff_num_denom.first);
1351  for (int i = 1; i < eff_num_denom.second->GetNbinsX(); ++i) {
1352  if (eff_num_denom.second->GetBinContent(i) < 1.e-5)
1353  eff_num_denom.second->SetBinContent(i, 1.);
1354  }
1355  _interactingEfficiencyDenoms.push_back(eff_num_denom.second);
1356 
1357  //TGraphAsymmErrors * eff = new TGraphAsymmErrors(eff_num_denom.first,
1358  // eff_num_denom.second);
1359 
1360  TH1D * interacting_numerator = new TH1D("incident_numerator", "",
1361  _TruthBinning.size()-1, 0,
1362  _TruthBinning.size()-1);
1363 
1364  std::cout << "n bins: " << eff_num_denom.first->GetNbinsX() << " " <<
1365  eff_num_denom.second->GetNbinsX() << " " <<
1366  interacting_numerator->GetNbinsX() << std::endl;
1367 
1368  for (size_t j = 0; j < signal_hists_by_topo[_SignalTopologyName[i]].size(); ++j) {
1369  interacting_numerator->SetBinContent(
1370  j+1, signal_hists_by_topo[_SignalTopologyName[i]][j]->Integral());
1371  }
1372  //eff_num_denom.second->Scale(_ScaleFactor);
1373  TGraphAsymmErrors * eff = new TGraphAsymmErrors(interacting_numerator, eff_num_denom.second);
1374 
1375  std::string name = "MC_Channel_" + _ChannelNames[i] + "_" +
1376  _SignalTopologyName[i] + "_Interacting_Efficiency";
1377 
1378  eff->SetNameTitle(name.c_str(), "Efficiency");
1379 
1380  _interactingEfficiencies.push_back(eff);
1381 
1382  }
1383 
1384 
1385 
1386  return true;
1387 
1388 }
static QCString name
Definition: declinfo.cpp:673
TH1 * FillMCIncidentHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string topo, int toponum, double reco_beam_endZ_cut, double minval=0., double maxval=100000., bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< TH1 * > _interactingEfficiencyNums
Definition: ProtoDUNEFit.h:138
std::vector< std::string > _IncidentDataFileNames
Definition: ProtoDUNEFit.h:109
std::vector< std::string > _DataFileNames
Definition: ProtoDUNEFit.h:109
TH1 * FillMCBackgroundHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double endZ_cut, double minval=0.0, double maxval=100000.0, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< std::string > _SignalTopologyName
Definition: ProtoDUNEFit.h:109
TGraphAsymmErrors * _incidentEfficiency
Definition: ProtoDUNEFit.h:134
TH1 * FillDataHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, double reco_beam_endZ_cut, bool doNegativeReco=false, bool IsIncidentHisto=false)
std::string string
Definition: nybbler.cc:12
std::vector< int > _AddBackgroundFactors
Definition: ProtoDUNEFit.h:147
std::vector< size_t > _sig_truth_index
Definition: ProtoDUNEFit.h:127
std::vector< std::string > _BackgroundTopologyName
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _incsighistos
Definition: ProtoDUNEFit.h:129
std::vector< TH1 * > _incdatahistos
Definition: ProtoDUNEFit.h:129
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< size_t > _sig_topo_index
Definition: ProtoDUNEFit.h:127
std::vector< int > _enable_bkg_factor
Definition: ProtoDUNEFit.h:147
std::vector< TH1 * > _bkghistos
Definition: ProtoDUNEFit.h:124
std::vector< size_t > _sig_chan_index
Definition: ProtoDUNEFit.h:126
std::vector< TH1 * > _interactingEfficiencyDenoms
Definition: ProtoDUNEFit.h:137
std::vector< int > _BackgroundTopology
Definition: ProtoDUNEFit.h:107
std::vector< std::string > _ChannelNames
Definition: ProtoDUNEFit.h:109
std::vector< double > _TruthBinning
Definition: ProtoDUNEFit.h:104
std::vector< std::string > _MCFileNames
Definition: ProtoDUNEFit.h:109
std::string _TruthTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _AddIncidentBackgroundFactors
Definition: ProtoDUNEFit.h:147
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::vector< TGraphAsymmErrors * > _interactingEfficiencies
Definition: ProtoDUNEFit.h:139
std::pair< TH1 *, TH1 * > GetMCIncidentEfficiency(std::string fileName, std::string treeName, std::vector< double > bins, double reco_beam_endZ_cut, int doSyst=0, double weight=1.)
std::vector< int > _IncidentTopology
Definition: ProtoDUNEFit.h:107
std::vector< size_t > _inc_sig_topo_index
Definition: ProtoDUNEFit.h:127
std::vector< std::string > _IncidentTopologyName
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
std::vector< int > _enable_inc_bkg_factor
Definition: ProtoDUNEFit.h:147
std::vector< TH1 * > _datahistos
Definition: ProtoDUNEFit.h:124
std::vector< size_t > _bkg_chan_index
Definition: ProtoDUNEFit.h:126
std::pair< TH1 *, TH1 * > GetMCInteractingEfficiency(std::string fileName, std::string treeName, std::vector< double > bins, std::string channel, std::string topo, int toponum, double endZ_cut, int doSyst=0, double weight=1.)
std::vector< int > _SignalTopology
Definition: ProtoDUNEFit.h:107
std::vector< TH1 * > _incbkghistos
Definition: ProtoDUNEFit.h:129
std::vector< size_t > _bkg_topo_index
Definition: ProtoDUNEFit.h:127
TH1 * FillMCSignalHistogram_Pions(std::string filename, std::string treename, std::vector< double > recoBins, std::string channel, std::string topo, int toponum, double minval, double maxval, double endZ_cut, bool doNegativeReco=false, int doSyst=0, std::string systName="", double weight=1.0, std::pair< double, double > PiMuScale={1., 1.})
std::vector< double > _RecoBinning
Definition: ProtoDUNEFit.h:104
std::vector< TH1 * > _sighistos
Definition: ProtoDUNEFit.h:124
std::vector< size_t > _inc_bkg_topo_index
Definition: ProtoDUNEFit.h:127
QTextStream & endl(QTextStream &s)
bool protoana::ProtoDUNEFit::FillSidebandHistograms_Pions ( )
private

Definition at line 1391 of file ProtoDUNEFit.cxx.

1391  {
1392 //********************************************************************
1393 
1394  std::string message_source = "FillSidebandHistograms_Pions";
1395 
1396  if (_MCControlSampleFiles.size() != _DataControlSampleFiles.size()) {
1397  mf::LogError(message_source) << "The number of data and MC control " <<
1398  "samples is not the same. Check MCControlSampleFiles and " <<
1399  "DataControlSampleFiles in the fcl file.";
1400  return false;
1401  }
1402 
1403  if (_SidebandTopologyName.size() != _SidebandTopology.size()) {
1404  mf::LogError(message_source) << "The number of Sideband topologies " <<
1405  "does not match the number of sideband topology names. Check " <<
1406  "SidebandTopologyName and SidebandTopology in the fcl file.";
1407  return false;
1408  }
1409 
1410  for (size_t i = 0; i < _SidebandTopologyName.size(); ++i) {
1411  _sideband_hists_mc.push_back(
1416  }
1417 
1418  std::cout << "GOT " << _sideband_hists_mc.size() << " SIDEBAND HISTS" <<
1419  std::endl;
1420 
1421  _sideband_hists_data.push_back(
1425 
1426  /*
1427  if (_StatFluctuation) {
1428  for (int j = 0; j <= _sideband_hists_data.back()->GetNbinsX(); ++j) {
1429  double new_val = rand.Poisson(_sideband_hists_data.back()->GetBinContent(j));
1430  _sideband_hists_data.back()->SetBinContent(j, new_val);
1431  }
1432  }
1433  */
1434 
1435  return true;
1436 }
std::vector< std::string > _SidebandTopologyName
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< TH1 * > _sideband_hists_mc
Definition: ProtoDUNEFit.h:130
std::vector< double > _SidebandBinning
Definition: ProtoDUNEFit.h:105
std::vector< int > _SidebandTopology
Definition: ProtoDUNEFit.h:116
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
TH1 * FillDataSidebandHistogram_Pions(std::string filename, std::string treename, std::string channel, double endZ_cut, std::vector< double > binning, double weight=1.)
std::vector< std::string > _DataControlSampleFiles
Definition: ProtoDUNEFit.h:109
std::vector< TH1 * > _sideband_hists_data
Definition: ProtoDUNEFit.h:130
TH1 * FillMCSidebandHistogram_Pions(std::string filename, std::string treename, std::string channel, std::string topo, int toponum, double endZ_cut, std::vector< double > binning, double weight=1.)
QTextStream & endl(QTextStream &s)
std::vector< std::string > _MCControlSampleFiles
Definition: ProtoDUNEFit.h:109
void protoana::ProtoDUNEFit::ScaleMCToData ( bool  data_is_mc = false)
private

Definition at line 1439 of file ProtoDUNEFit.cxx.

1439  {
1440 //********************************************************************
1441  //Get the number of incident pions from MC
1442  double nIncidentPionsMC = 0.;
1443  double nPrimaryPionsMC = 0.;
1444  for (size_t i = 0; i < _IncidentMCFileNames.size(); ++i) {
1445  //Get the tree
1446  TFile incidentFile(_IncidentMCFileNames[i].c_str(), "OPEN");
1447  TTree * incident_tree = (TTree*)incidentFile.Get(_RecoTreeName.c_str());
1448 
1449  //std::string cut = "(true_beam_PDG == 211 || true_beam_PDG == -13)";
1450  nIncidentPionsMC += incident_tree->GetEntries();
1451 
1452  std::string cut = "passBeamCut && primary_ends_inAPA3"; // && primary_passes_chi2";
1453 
1454  nPrimaryPionsMC += incident_tree->GetEntries(cut.c_str());
1455  incidentFile.Close();
1456  }
1457  std::cout << "Got " << nIncidentPionsMC << " incident MC Pions" << std::endl;
1458 
1459  //Get the number of incident pions from Data
1460  double nIncidentPionsData = 0.;
1461  double nPrimaryPionsData = 0.;
1462  for (size_t i = 0; i < _IncidentDataFileNames.size(); ++i) {
1463  //Get the tree
1464  TFile incidentFile(_IncidentDataFileNames[i].c_str(), "OPEN");
1465  TTree * incident_tree = (TTree*)incidentFile.Get(_RecoTreeName.c_str());
1466 
1467  //std::string cut = (data_is_mc ?
1468  // "(true_beam_PDG == 211 || true_beam_PDG == -13)" :
1469  // "data_BI_PDG_candidates[0] == 13");
1470  nIncidentPionsData += incident_tree->GetEntries();
1471 
1472  std::string cut = "passBeamCut && primary_ends_inAPA3";//&& primary_passes_chi2";
1473 
1474  nPrimaryPionsData += incident_tree->GetEntries(cut.c_str());
1475 
1476  incidentFile.Close();
1477  }
1478 
1479 
1480  //_ScaleFactor = nPrimaryPionsData/nPrimaryPionsMC;
1481  _IncidentScaleFactor = nIncidentPionsData/nIncidentPionsMC;
1483  std::cout << "Scale factor: " << _ScaleFactor << std::endl;
1484  std::cout << "Incident scale factor: " << _IncidentScaleFactor << std::endl;
1485 
1486  /*
1487  for(size_t i = 0; i < _sighistos.size(); i++){
1488  _sighistos[i]->Scale(scale_factor);
1489  }
1490  for(size_t i = 0; i < _bkghistos.size(); i++){
1491  _bkghistos[i]->Scale(scale_factor);
1492  }
1493  for(size_t i = 0; i < _incsighistos.size(); i++){
1494  _incsighistos[i]->Scale(scale_factor);
1495  }
1496  for(size_t i = 0; i < _incbkghistos.size(); i++){
1497  _incbkghistos[i]->Scale(scale_factor);
1498  }
1499  */
1500 
1501 }
std::vector< std::string > _IncidentDataFileNames
Definition: ProtoDUNEFit.h:109
std::string string
Definition: nybbler.cc:12
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
QTextStream & endl(QTextStream &s)
void protoana::ProtoDUNEFit::ScaleMuonContent ( )
private

Definition at line 1504 of file ProtoDUNEFit.cxx.

1504  {
1505  //Get the number of incident pions from MC
1506  double nPiMC = 0.;
1507  double nMuMC = 0.;
1508  for (size_t i = 0; i < _IncidentMCFileNames.size(); ++i) {
1509  //Get the tree
1510  TFile incidentFile(_IncidentMCFileNames[i].c_str(), "OPEN");
1511  TTree * incident_tree = (TTree*)incidentFile.Get(_RecoTreeName.c_str());
1512 
1513  nPiMC += incident_tree->GetEntries("true_beam_PDG == 211");
1514  nMuMC += incident_tree->GetEntries("true_beam_PDG == -13");
1515  incidentFile.Close();
1516  }
1517 
1518  _PionScaleFactor = ((nPiMC + nMuMC) - (_MuonScaleFactor*nMuMC)) / nPiMC;
1519  std::cout << "nMu & nPi: " << nMuMC << " " << nPiMC << std::endl;
1520  std::cout << "Muon & Pion Scale: " << _MuonScaleFactor << " " <<
1522 }
std::vector< std::string > _IncidentMCFileNames
Definition: ProtoDUNEFit.h:109
std::string _RecoTreeName
Definition: ProtoDUNEFit.h:102
QTextStream & endl(QTextStream &s)

Member Data Documentation

std::vector<int> protoana::ProtoDUNEFit::_AddBackgroundFactors
private

Definition at line 147 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_AddIncidentBackgroundFactors
private

Definition at line 147 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_AddIncidentToMeasurement
private

Definition at line 141 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_AddSidebandsToMeasurement
private

Definition at line 142 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_BackgroundTopology
private

Definition at line 107 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_BackgroundTopologyName
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_bkg_chan_index
private

Definition at line 126 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_bkg_topo_index
private

Definition at line 127 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_bkghistos
private

Definition at line 124 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_ChannelNames
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_DataControlSampleFiles
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_DataFileNames
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_datahistos
private

Definition at line 124 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DataIsMC
private

Definition at line 150 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DoAsimovFit
private

Definition at line 120 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DoDrawXSecs
private

Definition at line 149 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DoNegativeReco
private

Definition at line 141 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DoScaleMCToData
private

Definition at line 141 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_DoScaleMuonContent
private

Definition at line 143 of file ProtoDUNEFit.h.

std::vector<TGraphAsymmErrors*> protoana::ProtoDUNEFit::_efficiencyGraphs
private

Definition at line 132 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_enable_bkg_factor
private

Definition at line 147 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_enable_inc_bkg_factor
private

Definition at line 147 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_EnableMinosError
private

Definition at line 120 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_EnableStatisticalError
private

Definition at line 120 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_EnableSystematicError
private

Definition at line 121 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_EndZCut
private

Definition at line 151 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_FitInReco
private

Definition at line 121 of file ProtoDUNEFit.h.

int protoana::ProtoDUNEFit::_FitStrategy
private

Definition at line 118 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_IgnoreStatisticalErrorBelow
private

Definition at line 119 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_IgnoreSystematicErrorBelow
private

Definition at line 119 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_inc_bkg_topo_index
private

Definition at line 127 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_inc_sig_topo_index
private

Definition at line 127 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_incbkghistos
private

Definition at line 129 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_incdatahistos
private

Definition at line 129 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_IncidentDataFileNames
private

Definition at line 109 of file ProtoDUNEFit.h.

TGraphAsymmErrors* protoana::ProtoDUNEFit::_incidentEfficiency
private

Definition at line 134 of file ProtoDUNEFit.h.

TH1* protoana::ProtoDUNEFit::_incidentEfficiencyDenom
private

Definition at line 136 of file ProtoDUNEFit.h.

TH1* protoana::ProtoDUNEFit::_incidentEfficiencyNum
private

Definition at line 135 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_IncidentMCFileNames
private

Definition at line 109 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_IncidentScaleFactor = 1.
private

Definition at line 154 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_IncidentTopology
private

Definition at line 107 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_IncidentTopologyName
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_incsighistos
private

Definition at line 129 of file ProtoDUNEFit.h.

std::vector< TGraphAsymmErrors * > protoana::ProtoDUNEFit::_interactingEfficiencies
private

Definition at line 139 of file ProtoDUNEFit.h.

std::vector< TH1 * > protoana::ProtoDUNEFit::_interactingEfficiencyDenoms
private

Definition at line 137 of file ProtoDUNEFit.h.

std::vector< TH1 * > protoana::ProtoDUNEFit::_interactingEfficiencyNums
private

Definition at line 138 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_MCControlSampleFiles
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_MCFileNames
private

Definition at line 109 of file ProtoDUNEFit.h.

std::string protoana::ProtoDUNEFit::_Minimizer
private

Definition at line 102 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_MuonScaleFactor = 1.
private

Definition at line 156 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_NormalisedSystematic
private

Definition at line 121 of file ProtoDUNEFit.h.

int protoana::ProtoDUNEFit::_NToys
private

Definition at line 118 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_OnlyDrawXSecs
private

Definition at line 149 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_PionScaleFactor = 1.
private

Definition at line 156 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_RandSigPriors
private

Definition at line 144 of file ProtoDUNEFit.h.

std::vector<double> protoana::ProtoDUNEFit::_RecoBinning
private

Definition at line 104 of file ProtoDUNEFit.h.

std::string protoana::ProtoDUNEFit::_RecoTreeName
private

Definition at line 102 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_ScaleFactor = 1.
private

Definition at line 153 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_sideband_hists_data
private

Definition at line 130 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_sideband_hists_mc
private

Definition at line 130 of file ProtoDUNEFit.h.

std::vector<double> protoana::ProtoDUNEFit::_SidebandBinning
private

Definition at line 105 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_SidebandTopology
private

Definition at line 116 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_SidebandTopologyName
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_sig_chan_index
private

Definition at line 126 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_sig_topo_index
private

Definition at line 127 of file ProtoDUNEFit.h.

std::vector<size_t> protoana::ProtoDUNEFit::_sig_truth_index
private

Definition at line 127 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_sighistos
private

Definition at line 124 of file ProtoDUNEFit.h.

std::vector<int> protoana::ProtoDUNEFit::_SignalTopology
private

Definition at line 107 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_SignalTopologyName
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_syst_hists
private

Definition at line 125 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_SystFileNames
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_SystToConsider
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<std::string> protoana::ProtoDUNEFit::_SystType
private

Definition at line 109 of file ProtoDUNEFit.h.

std::vector<double> protoana::ProtoDUNEFit::_TruthBinning
private

Definition at line 104 of file ProtoDUNEFit.h.

std::vector<TH1*> protoana::ProtoDUNEFit::_truthsighistos
private

Definition at line 124 of file ProtoDUNEFit.h.

std::string protoana::ProtoDUNEFit::_TruthTreeName
private

Definition at line 102 of file ProtoDUNEFit.h.

bool protoana::ProtoDUNEFit::_UseComputedSysts
private

Definition at line 122 of file ProtoDUNEFit.h.

double protoana::ProtoDUNEFit::_WirePitch
private

Definition at line 151 of file ProtoDUNEFit.h.

TRandom3 protoana::ProtoDUNEFit::rand
private

Definition at line 145 of file ProtoDUNEFit.h.


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