PlotPreSelectionPlots.C File Reference

Go to the source code of this file.

Enumerations
enum	CutSide {   kLow, kHigh, kBoth, kNoCut,   kLow, kHigh, kBoth, kNoCut,   kLow, kHigh, kBoth, kNeither }
enum	FomType {   kSOverSB, kSOverSSB, kBinnedSOverSSB, kBinnedSOverSSBSq,   kEfficiency, kPurity, kEffTimesPurity, kSOverSB,   kSOverSSB, kEfficiency, kPurity, kEffTimesPurity,   kSOverSB, kSOverSSB, kEfficiency, kPurity,   kEffTimesPurity }

Functions
void	printValues (std::vector< double > theseValues, std::vector< double > theseFOMValues)
void	getMaximumFOM (TH1D *h)
double	calcSOverSB (double sig, double bg)
double	calcSOverSSB (double sig, double bg, double fracUn)
double	calcEfficiency (double selSig, double totSig)
double	calcPurity (double selSig, double selTot)
double	getFOM (FomType thisFomType, double selSig, double selBg, double totSig, double selTot, double fracUn)
double	getFracUncertainty (TH1D *h, std::vector< float > fracUncertainties)
TH1D *	getFOMHist (std::string fomName, CutSide thisCutSide, FomType thisFomType, TH1D *all, TH1D *sig, TH1D *trueSig, TH2D *energyCorr=nullptr, std::vector< float > fracUncertainties={})
void	PlotPreSelectionPlots ()

Enumeration Type Documentation

enum CutSide

Enumerator
kLow
kHigh
kBoth
kNoCut
kLow
kHigh
kBoth
kNoCut
kLow
kHigh
kBoth
kNeither

Definition at line 5 of file PlotPreSelectionPlots.C.

             {
  kLow,
  kHigh,
  kBoth,
  kNoCut
};

enum FomType

Enumerator
kSOverSB
kSOverSSB
kBinnedSOverSSB
kBinnedSOverSSBSq
kEfficiency
kPurity
kEffTimesPurity
kSOverSB
kSOverSSB
kEfficiency
kPurity
kEffTimesPurity
kSOverSB
kSOverSSB
kEfficiency
kPurity
kEffTimesPurity

Definition at line 12 of file PlotPreSelectionPlots.C.

             {
  kSOverSB,       // signal over sqrt background
  kSOverSSB,      // signal over sqrt (signal plus background)
  kEfficiency,    // efficiency (total true selected / total true)
                  //  -- n.b. "total true" is after the Concat cuts, so not quite a true purity
  kPurity,        // purity (total true selected / total selected)
  kEffTimesPurity // efficiency times purity
};

Function Documentation

double calcEfficiency	(	double	selSig,
		double	totSig
	)

Definition at line 70 of file PlotPreSelectionPlots.C.

Referenced by getFOM().

                                                   {
  if (totSig == 0) return 0;
  else return selSig/totSig;
}

double calcPurity	(	double	selSig,
		double	selTot
	)

Definition at line 76 of file PlotPreSelectionPlots.C.

Referenced by getFOM().

                                               {
  if (selTot == 0) return 0;
  else return selSig/selTot;
}

double calcSOverSB	(	double	sig,
		double	bg
	)

Definition at line 58 of file PlotPreSelectionPlots.C.

References std::sqrt().

Referenced by getFOM().

                                         {
  if (sig == 0) return 0;
  else return sig/std::sqrt(bg);
}

double calcSOverSSB	(	double	sig,
		double	bg,
		double	fracUn
	)

Definition at line 64 of file PlotPreSelectionPlots.C.

References cet::pow(), and std::sqrt().

Referenced by getFOM().

                                                         {
  if (sig == 0) return 0;
  else return sig/std::sqrt(sig+bg+std::pow(fracUn,2));
}

double getFOM	(	FomType	thisFomType,
		double	selSig,
		double	selBg,
		double	totSig,
		double	selTot,
		double	fracUn
	)

Definition at line 81 of file PlotPreSelectionPlots.C.

References calcEfficiency(), calcPurity(), calcSOverSB(), calcSOverSSB(), kEfficiency, kEffTimesPurity, kPurity, kSOverSB, and kSOverSSB.

Referenced by getFOMHist().

                                                                                                            {

  switch (thisFomType){
    case kSOverSB:
      return calcSOverSB(selSig, selBg);
    case kSOverSSB:
      return calcSOverSSB(selSig, selBg, fracUn);
    case kEfficiency:
      return calcEfficiency(selSig, totSig);
    case kPurity:
      return calcPurity(selSig, selTot);
    case kEffTimesPurity:
      return calcEfficiency(selSig, totSig) * calcPurity(selSig, selTot);
    default:
      throw std::logic_error("thisFomType not accepted");
  }

}

TH1D* getFOMHist	(	std::string	fomName,
		CutSide	thisCutSide,
		FomType	thisFomType,
		TH1D *	all,
		TH1D *	sig,
		TH1D *	trueSig,
		TH2D *	energyCorr = nullptr,
		std::vector< float >	fracUncertainties = {}
	)

Definition at line 121 of file PlotPreSelectionPlots.C.

References getFOM(), getFracUncertainty(), MECModelEnuComparisons::i, kBoth, kHigh, kLow, and proj.

Referenced by PlotPreSelectionPlots().

                                                      {}){

  TH1D* fomReturnHist = new TH1D(fomName.c_str(),
                                 ";;FOM",
                                 all->GetNbinsX(),
                                 all->GetBinLowEdge(1),
                                 all->GetBinLowEdge(all->GetNbinsX()+1));

  double selSig, totSig, selTot, selBg;
  if (thisCutSide == CutSide::kBoth) 
    return fomReturnHist;               
  else if (thisCutSide == CutSide::kLow){
    for (int i = 1; i < fomReturnHist->GetNbinsX()+1; i++){
      selTot = all->Integral(i, fomReturnHist->GetNbinsX());
      selSig = sig->Integral(i, fomReturnHist->GetNbinsX());
      selBg  = selTot - selSig;
      totSig = trueSig->Integral();

      TH1D* proj = nullptr;
      if (energyCorr != nullptr){
        proj = energyCorr->ProjectionX("proj", i, fomReturnHist->GetNbinsX());
      }
      double fracUncertainty = (proj == nullptr) ? 0 : getFracUncertainty(proj, fracUncertainties);

      fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, fracUncertainty));
    }
  }
  else if (thisCutSide == CutSide::kHigh){
    for (int i = fomReturnHist->GetNbinsX()+1; i > 0; i--){
      selTot = all->Integral(1,i);
      selSig = sig->Integral(1,i);
      selBg  = selTot - selSig;
      totSig = trueSig->Integral();

      TH1D* proj = nullptr;
      if (energyCorr != nullptr){
        proj = energyCorr->ProjectionX("proj", 1, i);
      }
      double fracUncertainty = (proj == nullptr) ? 0 : getFracUncertainty(proj, fracUncertainties);

      fomReturnHist->SetBinContent(i, getFOM(thisFomType, selSig, selBg, totSig, selTot, fracUncertainty));

    }
  }

  return fomReturnHist;
}

double getFracUncertainty	(	TH1D *	h,
		std::vector< float >	fracUncertainties
	)

Definition at line 100 of file PlotPreSelectionPlots.C.

References om::cout, allTimeWatchdog::endl, and MECModelEnuComparisons::i.

Referenced by getFOMHist().

                                                                      {

  if (h->GetNbinsX() != fracUncertainties.size()+1){
    std::cout << "hBins: " << h->GetNbinsX()
              << " fracUncertainties: " << fracUncertainties.size() << std::endl;
  }

  TH1D* hShifted = (TH1D*)h->Clone("hShifted");

  for (int i = 2; i < h->GetNbinsX()+1; ++i){
    hShifted->SetBinContent(i,h->GetBinContent(i)*(1+fracUncertainties.at(i-2))); 
  }


  double totUnc = (hShifted->Integral()/h->Integral())-1;

  return totUnc;

}

void getMaximumFOM

(

TH1D *

h

)

Definition at line 38 of file PlotPreSelectionPlots.C.

References MECModelEnuComparisons::i, and printValues().

Referenced by PlotPreSelectionPlots().

                           {
  
  std::vector<double> cutValues;
  std::vector<double> FOMValues;
  
  double maxFOM = h->GetMaximum();
  for (int i = 1; i < h->GetNbinsX()+1; ++i)
  {
    if (h->GetBinContent(i) == maxFOM)
    {
        cutValues.push_back(h->GetBinLowEdge(i));
        FOMValues.push_back(h->GetBinContent(i));
    }
  }

  printValues(cutValues, FOMValues);

}

void PlotPreSelectionPlots

(

)

Definition at line 177 of file PlotPreSelectionPlots.C.

References allInOneTrainingPlots::axis, demo5::c1, om::cout, Draw(), EffTimesPurHist, allTimeWatchdog::endl, getFOMHist(), getMaximumFOM(), MECModelEnuComparisons::i, kBoth, kEffTimesPurity, kGreen, kHigh, kLow, kNoCut, kOrange, kRed, kSOverSB, kSOverSSB, kViolet, MECModelEnuComparisons::leg, APDHVSetting::mod, Scale(), SetFillColor(), SetLineColor(), SetLineWidth(), SetTitle(), SOverSBHist, SOverSSBHist, std::sqrt(), string, thisPad, and topPad.

                            {

  TFile* covmxFile = new TFile("/pnfs/nova/persistent/analysis/nux/nus20/cmf_covmx/2020-02-14/covariance_matrix_hist_all.root", 
                               "read");

  TH2D* covmx = (TH2D*)covmxFile->Get("covar/covarianceMatrix");
  std::vector<float> fracUncertainty;

  for (int i = 93; i < 93+38; ++i){
    double mod = 1;
    if (i > 93+0 && i < 93+20)
      mod = 1;
    fracUncertainty.push_back(std::sqrt(covmx->GetBinContent(i,i))*mod);
  }

  std::vector< std::string > cutLevels = {
    "NoCut" ,
    "Quality",
    "QualityVtx",
    "QualityVtxFiducial",
    "FullSelection"
  };

  // first up define some vectors for stuff we're going to loop over
  std::vector< std::string > histNames1D = {
    "kHitsPerPlane",
    "kNSliceHits",
    "kPartPtp",
    "kIsElastic",
    "kNFuzzyProng",
    "kNContPlanes",
    "kVtxX",
    "kVtxY",
    "kVtxZ",
    "kDistAllTop",
    "kDistAllBottom",
    "kDistAllEast",
    "kDistAllWest",
    "kDistAllFront",
    "kDistAllBack",
    "kCVNnc",
    "kCVNnc_looseptp",
    "kCVNnc_oldpresel",
    "kCaloE"
  };

  std::vector< std::string > histTitles1D = {
    ";Hits Per Plane;Events;",
    ";Number Hits in Slice;Events",
    ";Beam-Transverse Momentum Fraction; Events",
    ";Number of Elastic Vertices;Events",
    ";Number of Fuzzy Prongs;Events",
    ";Number of Continuous Planes;Events",
    ";Vertex x (cm);Events",
    ";Vertex y (cm);Events",
    ";Vertex z (cm);Events",
    ";Minimum Distance To Top;Events",
    ";Minimum Distance To Bottom;Events",
    ";Minimum Distance To East;Events",
    ";Minimum Distance To West;Events",
    ";Minimum Distance To Front;Events",
    ";Minimum Distance To Back;Events",
    ";CVNnc Score;Events",
    ";CVNnc Score (looseptp);Events",
    ";CVNnc Score (oldpresel);Events",
    ";Calorimetric Energy (GeV);Events"
  };

  std::vector< CutSide > histCutSide = {
    kHigh  , // kHitsPerPlane
    kLow   , // kNSliceHits
    kHigh  , // kPartPtp
    kLow   , // kIsElastic
    kLow   , // kNFuzzyProng
    kLow   , // kNContPlanes
    kBoth  , // kVtxX
    kBoth  , // kVtxY
    kBoth  , // kVtxZ
    kLow   , // kDistAllTopHist
    kLow   , // kDistAllBottomHist
    kLow   , // kDistAllEastHist
    kLow   , // kDistAllWestHist
    kLow   , // kDistAllFrontHist
    kLow   , // kDistAllBackHist
    kLow   , // kCVNnc
    kLow   , // kCVNnc_looseptpHist
    kLow   , // kCVNnc_oldpreselHist
    kNoCut // kCaloE 
  };

  std::vector< bool > plotFOM = {
    true,
    true,
    true,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    false,
    true,
    true,
    true,
    false
  };

  std::vector< std::vector< double > > currentCuts = {
    {0.0}       ,  // kHitsPerPlane
    {0.0}       ,  // kNSliceHits
    {1.0}       ,  // kPartPtp
    {1}         ,  // kIsElastic
    {1}         ,  // kNFuzzyProng
    {2}         ,  // kNContPlanes
    {-100, 100} ,  // kVtxX
    {-100, 100} ,  // kVtxY
    {150 , 1000},  // kVtxZ
    {20}        ,  // kDistAllTopHist
    {20}        ,  // kDistAllBottomHist
    {20}        ,  // kDistAllEastHist
    {20}        ,  // kDistAllWestHist
    {150}       ,  // kDistAllFrontHist
    {50}        ,  // kDistAllBackHist
    {0}         ,  // kCVNnc
    {0.42}      ,  // kCVNnc_looseptpHist
    {0}         ,  // kCVNnc_oldpreselHist
    {0,0}          // kCaloE
  };

  std::vector< std::string > histNames2D = {
    "kHitsPerPlaneNSliceHitsHist",
    "kHitsPerPlanePartPtpHist",
    "kNSliceHitsPartPtpHist"
  };

  std::vector< std::string > histTitles2D = {
    ";Hits Per Plane; Number of Hits in Slice",
    ";Hits Per Plane; kPartPtp",
    ";Number of Hits in Slice; kPartPtp"
  };

  double maxVal   = 1.5;
  int    textSize = 28;
  // fist deal with the 1d histograms
  for (int iH = 0; iH < histNames1D.size(); iH++){

    TCanvas* c1Ratio = new TCanvas("c1Ratio", "c1Ratio", 600, 600);
    TPad *topPad     = new TPad("topPad", "", 0.005, 0.3, 0.995, 0.995);
    TPad *bottomPad  = new TPad("bottomPad", "", 0.005, 0.005, 0.995, 0.3);
    topPad   ->SetTopMargin(0.1);
    topPad   ->SetBottomMargin(0.005);
    bottomPad->SetTopMargin(0.005);
    bottomPad->SetBottomMargin(0.35);
    bottomPad->SetGridy();
    topPad   ->Draw();
    bottomPad->Draw();

    TCanvas* c1 = new TCanvas("c1", "c1", 600, 500);
    TPad *thisPad = new TPad("topPad", "", 0.005, 0.005, 0.995, 0.995);
    thisPad ->SetTopMargin(0.1);
    thisPad ->SetBottomMargin(0.15);
    thisPad->Draw();
    thisPad->cd();

    for (int icl = 0; icl < cutLevels.size(); ++icl){
      std::cout << std::string("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)) << std::endl;
      TH1D* histSim                = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)).c_str());
      TH1D* histSimTrueNC          = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)+"TrueNC").c_str());
      TH1D* histSimTrueNCNoCut     = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+"NoCutTrueNC").c_str());
      TH1D* histDatNoCut           = (TH1D*)_file0->Get(("dat"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)).c_str());

      // get 2D histograms for variables correlated with energy
      std::string name = histNames1D.at(iH);
      std::string simEnergyString = "sim"+name+"CaloEHist"+cutLevels.at(icl);
      TH2D* histSimEnergy = (TH2D*)_file0->Get(simEnergyString.c_str());

      double scalef = 1;
      if (histDatNoCut != nullptr)
        scalef = (double)histDatNoCut->Integral()/histSim->Integral();

      histSim                -> Sumw2();
      histSim                -> Scale(scalef);
      histSimTrueNC          -> Sumw2();
      histSimTrueNC          -> Scale(scalef);
      histSimTrueNCNoCut     -> Sumw2();
      histSimTrueNCNoCut     -> Scale(scalef);

      histSim ->SetLineColor(kGray+1);
      histSimTrueNC->SetFillColor(kAzure+1);
      histSimTrueNC->SetLineColor(kAzure+1);
      histSimTrueNC->SetMarkerColor(kAzure+1);

      histSim ->SetLineWidth(2);
      histSim ->GetYaxis()->SetRangeUser(0.0001, histSim->GetMaximum()*maxVal);
      histSim ->SetTitle(histTitles1D.at(iH).c_str());
      histSim ->GetYaxis()->SetTitleOffset(1.25);
      histSim ->GetYaxis()->SetMaxDigits(3);
      histSim ->GetYaxis()->SetNdivisions(515);
      histSim ->GetXaxis()->SetTitleFont(43);
      histSim ->GetXaxis()->SetTitleSize(textSize);
      histSim ->GetXaxis()->CenterTitle();
      histSim ->GetYaxis()->SetTitleFont(43);
      histSim ->GetYaxis()->SetTitleSize(textSize);
      histSim ->GetYaxis()->CenterTitle();
      histSim ->GetYaxis()->SetLabelFont(43);
      histSim ->GetYaxis()->SetLabelSize(textSize);
      histSim ->GetXaxis()->SetLabelFont(43);
      histSim ->GetXaxis()->SetLabelSize(textSize);
      bool isDrawData = false;
      if ((histNames1D.at(iH).find("CVN") == std::string::npos && 
            histNames1D.at(iH).find("CaloE") == std::string::npos) || 
          (histNames1D.at(iH).find("CVN") != std::string::npos &&
            cutLevels.at(icl).find("QualityVtxFiducial") != std::string::npos) ||
          (histNames1D.at(iH).find("CVN") != std::string::npos &&
            cutLevels.at(icl).find("FullSelection") != std::string::npos) ||
          (histNames1D.at(iH).find("CaloE") != std::string::npos &&
            cutLevels.at(icl).find("FullSelection") != std::string::npos))
        isDrawData = true;

      if(isDrawData){
        c1Ratio->cd();
        topPad->cd();
        histSim ->GetYaxis() ->SetTitleOffset(1.3);
        histSim ->GetYaxis()->SetMaxDigits(3);
        histSim ->GetYaxis()->SetNdivisions(515);

      }

      histSim       ->Draw("hist");
      histSimTrueNC ->Draw("same hist");
      histSim       ->Draw("same hist");
      if (isDrawData){
        histDatNoCut ->SetMarkerStyle(20);
        histDatNoCut ->SetMarkerSize(0.6);

        histDatNoCut ->GetXaxis()->SetTitleFont(43);
        histDatNoCut ->GetXaxis()->SetTitleSize(textSize);
        histDatNoCut ->GetXaxis()->CenterTitle();
        histDatNoCut ->GetYaxis()->SetTitleFont(43);
        histDatNoCut ->GetYaxis()->SetTitleSize(textSize);
        histDatNoCut ->GetYaxis()->CenterTitle();
        histDatNoCut ->GetYaxis()->SetLabelFont(43);
        histDatNoCut ->GetYaxis()->SetLabelSize(textSize);
        histDatNoCut ->GetXaxis()->SetLabelFont(43);
        histDatNoCut ->GetXaxis()->SetLabelSize(textSize);
        histDatNoCut ->SetTitle(histTitles1D.at(iH).c_str());


        histDatNoCut ->Draw("same p");
        bottomPad->cd();
        TH1D* ratioHist = (TH1D*)histDatNoCut->Clone("ratioHist");
        ratioHist->Divide(histSim);
        ratioHist->GetYaxis()->SetRangeUser(0.7, 1.29);
        ratioHist->GetYaxis()->SetNdivisions(505);
        ratioHist->GetYaxis()->SetTitle("Data/Sim");
        ratioHist->GetYaxis()->SetTitleOffset(1.3);
        ratioHist->GetXaxis()->SetTitleOffset(3.0);
        ratioHist->Draw("p");
        topPad->cd();
      }
      TH1D* EffTimesPurHist = (TH1D*)getFOMHist("EffTimesPur"+histNames1D.at(iH)+cutLevels.at(icl), 
                                                histCutSide.at(iH), 
                                                FomType::kEffTimesPurity, 
                                                histSim, 
                                                histSimTrueNC, 
                                                histSimTrueNCNoCut);

      TH1D* SOverSSBHist    = (TH1D*)getFOMHist("SOverSSB"+histNames1D.at(iH)+cutLevels.at(icl), 
                                                histCutSide.at(iH), 
                                                FomType::kSOverSSB,
                                                histSim, 
                                                histSimTrueNC, 
                                                histSimTrueNCNoCut,
                                                histSimEnergy,
                                                fracUncertainty);

      TH1D* SOverSBHist     = (TH1D*)getFOMHist("SOverSB"+histNames1D.at(iH)+cutLevels.at(icl), 
                                                histCutSide.at(iH), 
                                                FomType::kSOverSB,
                                                histSim, 
                                                histSimTrueNC, 
                                                histSimTrueNCNoCut);

      // define new axis on the right hand side for the EffTimesPurHist
      TLegend* legFOM = new TLegend(0.5, 0.70, 0.8, 0.89);
      if (plotFOM.at(iH) == true){
 
        std::cout 
          << "Printing cut values: " 
          << "\n-- Efficiency x Purity: "
          << std::endl;
        getMaximumFOM(EffTimesPurHist);
        std::cout 
          << "\n-- SOverSSB:"
          << std::endl;
        getMaximumFOM(SOverSSBHist);
         std::cout 
          << "\n-- SOverSB:"
          << std::endl;
        getMaximumFOM(SOverSBHist);
     
        double scaleval;
        if (isDrawData)
          scaleval = histDatNoCut->GetMaximum();
        else
          scaleval = histSim->GetMaximum()/maxVal;

        EffTimesPurHist->GetYaxis()->SetRangeUser(0,EffTimesPurHist->GetMaximum()*maxVal);
        EffTimesPurHist->Scale(scaleval/EffTimesPurHist->GetMaximum());
        EffTimesPurHist->SetLineColor(kGreen-3);
        EffTimesPurHist->SetLineWidth(2);
        EffTimesPurHist->Draw("same hist");

        SOverSSBHist->Scale(1./SOverSSBHist->GetMaximum() * scaleval);
        SOverSSBHist->SetLineColor(kViolet+1);
        SOverSSBHist->SetLineWidth(2);
        SOverSSBHist->Draw("same hist");

        SOverSBHist->Scale(0.8/SOverSBHist->GetMaximum() * scaleval);
        SOverSBHist->SetLineColor(kOrange+1);
        SOverSBHist->SetLineWidth(2);
        SOverSBHist->Draw("same hist");

        legFOM->SetFillStyle(0);
        legFOM->SetLineWidth(0);
        legFOM->SetTextFont(43);
        legFOM->SetTextSize(textSize*0.75);
        legFOM->AddEntry(EffTimesPurHist, "Eff #times Purity", "l");
        legFOM->AddEntry(SOverSBHist    , "S/#sqrt{B}", "l");
        legFOM->AddEntry(SOverSSBHist   , "S/#sqrt{S+B+#sigma^{2}_{syst}}", "l");

        gPad->RedrawAxis();

        //draw an axis on the right side
        TGaxis *axis = new TGaxis(EffTimesPurHist->GetBinLowEdge(EffTimesPurHist->GetNbinsX()+1), 0.,
                                  EffTimesPurHist->GetBinLowEdge(EffTimesPurHist->GetNbinsX()+1), histSim->GetMaximum(),
                                  0.0001,maxVal,510,"+L");

        axis->SetTitleFont(43);
        axis->SetTitleSize(textSize);
        axis->SetLabelFont(43);
        axis->SetLabelSize(textSize);
        axis->SetTitle("Figure Of Merit");
        axis->SetTitleOffset(1.4);
        axis->CenterTitle();
        axis->Draw();
      }
      // draw boxes
      TBox* box  = new TBox();
      TBox* box2 = new TBox(-5000, -5000, -5000, -5000); 
      box->SetY1(0);
      box->SetY2(histSim->GetMaximum());
      if (histCutSide.at(iH) == kLow){
        box->SetX1(histSim->GetBinLowEdge(1));
        box->SetX2(currentCuts.at(iH).at(0));
      }
      if (histCutSide.at(iH) == kHigh){
        box->SetX1(currentCuts.at(iH).at(0));
        box->SetX2(histSim->GetBinLowEdge(histSim->GetNbinsX()+1));
      }
      if (histCutSide.at(iH) == kBoth){
        box->SetX1(histSim->GetBinLowEdge(1));
        box->SetX2(currentCuts.at(iH).at(0));
        box2->SetX1(currentCuts.at(iH).at(1));
        box2->SetX2(histSim->GetBinLowEdge(histSim->GetNbinsX()+1));
        box2->SetLineWidth(0);
        box2->SetFillColorAlpha(kGray, 0.4);
        box2->SetY1(0);
        box2->SetY2(histSim->GetMaximum());
        box2->Draw("same");
      }
      box->SetLineWidth(0);
      box->SetFillColorAlpha(kGray, 0.4);
      if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
        box->Draw("same");

      //gPad->Update();
      gPad->RedrawAxis();
      gPad->Modified();

      if (plotFOM.at(iH) == true)
        legFOM->Draw("same");

      TLegend* leg = new TLegend(0.2, 0.70, 0.5, 0.89);
      leg->AddEntry(histSim , "All Events"    , "l" );
      leg->AddEntry(histSimTrueNC, "True NC Events", "f" );
      if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
        leg->AddEntry(box       , "Proposed Cuts"  , "lf");
      if (isDrawData)
        leg->AddEntry(histDatNoCut, "Data", "p");
      leg->SetLineWidth(0);
      leg->SetFillStyle(0);
      leg->SetTextFont(43);
      leg->SetTextSize(textSize*0.75);
      leg->Draw("same");

      if (!isDrawData){
        c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".png").c_str());
        c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".pdf").c_str());
      }
      else{ 
        c1Ratio->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".png").c_str());
        c1Ratio->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+".pdf").c_str());
      }
      thisPad->cd();

      // get out of detector plots
      TH1D* histSimTrueVtxOutOfDet = (TH1D*)_file0->Get(("sim"+histNames1D.at(iH)+"Hist"+cutLevels.at(icl)+"TrueVtxOutOfDet" ).c_str());
      if (histSimTrueVtxOutOfDet == nullptr) continue;
      histSimTrueVtxOutOfDet -> Sumw2();
      histSimTrueVtxOutOfDet -> Scale(scalef);
      histSimTrueVtxOutOfDet -> SetLineColor(kRed+1);
      histSimTrueVtxOutOfDet -> SetFillColor(kRed+1);
      histSimTrueVtxOutOfDet -> SetLineWidth(2);
      histSim                -> Draw("hist");
      histSimTrueVtxOutOfDet -> Draw("same hist");
      histSim                -> Draw("same hist");

      if (currentCuts.at(iH).at(0) != 0 || currentCuts.at(iH).size() > 1)
        box->Draw("same");
      if (histCutSide.at(iH) == kBoth)
        box2->Draw("same");

      //gPad->Update();
      gPad->RedrawAxis();
      gPad->Modified();

      TLegend *leg2 = new TLegend(0.2, 0.7, 0.5, 0.89);
      leg2->AddEntry(histSim, "All Events", "l");
      leg2->AddEntry(histSimTrueVtxOutOfDet, "Rock Events", "f");
      leg2->AddEntry(box       , "Proposed Cuts"  , "lf");
      leg2->SetLineWidth(0);
      leg2->SetFillStyle(0);
      leg2->SetTextFont(43);
      leg2->SetTextSize(textSize*0.75);

      leg2->Draw("same");

      c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+"containment.png").c_str());
      c1->SaveAs((histNames1D.at(iH)+cutLevels.at(icl)+"containment.pdf").c_str());
    }
  }

  for (int iH = 0; iH < histNames2D.size(); iH++){

    TCanvas* c1   = new TCanvas("c1", "c1", 600, 400);
    TPad *thisPad = new TPad("topPad", "", 0.005, 0.005, 0.995, 0.995);
    thisPad ->SetTopMargin(0.1);
    thisPad ->SetBottomMargin(0.15);
    thisPad->Draw();

 
    c1->SetRightMargin(0.18);
    c1->cd();
    thisPad->cd();

    for (int icl = 0; icl < cutLevels.size(); ++icl){
      TH2D* histSim  = (TH2D*)_file0->Get(("sim"+histNames2D.at(iH)+cutLevels.at(icl)).c_str());
      TH2D* histSimTrueNC = (TH2D*)_file0->Get(("sim"+histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC").c_str());

      histSim       -> SetTitle(("Background"+histTitles2D.at(iH)).c_str());
      histSimTrueNC -> SetTitle(("True NC in Detector"+histTitles2D.at(iH)).c_str());

      gStyle      -> SetTitleFont(43);
      gStyle      -> SetTitleFontSize(0.06);

      histSim ->GetXaxis()->SetTitleFont(43);
      histSim ->GetXaxis()->SetTitleSize(textSize);
      histSim ->GetXaxis()->CenterTitle();
      histSim ->GetYaxis()->SetTitleFont(43);
      histSim ->GetYaxis()->SetTitleSize(textSize);
      histSim ->GetYaxis()->CenterTitle();
      histSimTrueNC ->GetXaxis()->SetTitleFont(43);
      histSimTrueNC ->GetXaxis()->SetTitleSize(textSize);
      histSimTrueNC ->GetXaxis()->CenterTitle();
      histSimTrueNC ->GetYaxis()->SetTitleFont(43);
      histSimTrueNC ->GetYaxis()->SetTitleSize(textSize);
      histSimTrueNC ->GetYaxis()->CenterTitle();

      gPad->Modified();
      gPad->Update();

      histSim->Add(histSimTrueNC, -1);
      histSim->Draw("colz");
      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+".png").c_str());
      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+".pdf").c_str());

      histSimTrueNC->Draw("colz");
      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC.png").c_str());
      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"TrueNC.pdf").c_str());
      
      histSim->Draw("box");
      histSimTrueNC->SetLineColor(kAzure+1);
      histSimTrueNC->Draw("box same");

      TLegend* leg = new TLegend(0.6, 0.75, 0.75, 0.85);
      leg->AddEntry(histSim, "All Events");
      leg->AddEntry(histSimTrueNC, "True NC Events");
      leg->SetLineWidth(0);
      leg->SetFillStyle(0);
      leg->SetTextFont(43);
      leg->SetTextSize(textSize*0.75);
      leg->Draw("same");

      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"together.png").c_str());
      c1->SaveAs((histNames2D.at(iH)+cutLevels.at(icl)+"together.pdf").c_str());
    }
  }
}

void printValues	(	std::vector< double >	theseValues,
		std::vector< double >	theseFOMValues
	)

Definition at line 21 of file PlotPreSelectionPlots.C.

References om::cout, allTimeWatchdog::endl, and MECModelEnuComparisons::i.

Referenced by getMaximumFOM().

                                                                                 {

  std::cout << "cut value: ";
  for (int i = 0; i < theseValues.size(); ++i){
    std::cout << " " << theseValues.at(i);
  }
  std::cout << std::endl;

  std::cout << "FOM value: ";
  for (int i = 0; i < theseFOMValues.size(); ++i){
    std::cout << " " << theseFOMValues.at(i);
  }
  std::cout << "\n" << std::endl;

}