vertex_optimize.C
Go to the documentation of this file.
2 
4 #include "CAFAna/Core/HistAxis.h"
5 #include "CAFAna/Vars/Vars.h"
6 #include "CAFAna/Core/Binning.h"
7 #include "CAFAna/Core/Cut.h"
8 #include "CAFAna/Core/Ratio.h"
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Cuts/TruthCuts.h"
13 #include "3FlavorAna/Cuts/NueCutsFirstAna.h"
14 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
21 
22 
23 
24 #include "TFile.h"
25 #include "CAFAna/Core/Spectrum.h"
26 #include "TH1.h"
27 #include "TH2.h"
28 
29 #include "CAFAna/Analysis/Plots.h"
30 #include "TLegend.h"
31 #include "TCanvas.h"
32 #include "TColor.h"
33 
38 
39 #include <stdio.h>
40 #include <string.h>
41 
42 using namespace ana;
43 
44 //------------------------------------------------------------
45 
46 void vertex_optimize(std::string filename = "fVertexOpt.root")
47 {
48  TFile * fout = new TFile(filename.c_str(),"RECREATE");
49  if (!fout || fout == NULL || fout->IsZombie())
50  throw runtime_error("Could not open output file.");
51 
52  std::vector<std::string> dataset_labels = {"nominal", "lightdown", "lightup", "ckv", "calibneg", "calibpos", "calibshape"};
53 
54  std::vector<ana::SpectrumLoaderBase*> loaders;
55 
56  const std::vector<std::string> datasets =
57  { "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1_minus_muonid_training",
58  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1",
59  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1",
60  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1",
61  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1",
62  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1",
63  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1"
64  };
65 
66  // Print all the inputs for easy checking
67  std::cout<<"\nDataset names:\n";
68  for(auto &id: datasets)
69  std::cout<<id<<"\n";
70 
71  // Preselection here are just quality and containment. Fiducial become the optimization function
72  const Cut preselMinusFiducial =
73  kNumuMyQuality && // ensures reconstruction
74  kNumuTightContainND; // containment
75 
76  // Selection is simply presel && MuonID.
77  const Cut selectionCut =
78  preselMinusFiducial && kMuonIDCut; // Muon selection > 0.24 (reoptimized)
79 
80  // Full signal cut goes from
81  const Cut kIsTruthSig = kIsNumuCC && kTrueDetector; // true signal events
82 
83  const std::vector<Cut> cuts = {
84  preselMinusFiducial,
85  preselMinusFiducial && kIsTruthSig,
86  preselMinusFiducial && !kIsTruthSig,
87  selectionCut,
88  selectionCut && kIsTruthSig,
89  selectionCut && !kIsTruthSig,
90  kIsTruthSig
91  };
92 
93  const std::vector<std::string> cut_labels = {
94  "presel",
95  "presel_sig",
96  "presel_bkg",
97  "selection",
98  "sel_sig",
99  "sel_bkg",
100  "signal"
101  };
102 
103  assert(cuts.size() == cut_labels.size());
104 
105  // all the vars we want
106  const int mvar = 3;
107 
108  std::vector<HistAxis> taxes = {
109  HistAxis("Vertex in X Direction (cm)", Binning::Simple( 20, -200.0, 200.0), kMuStartX),
110  HistAxis("Vertex in Y Direction (cm)", Binning::Simple( 20, -200.0, 200.0), kMuStartY),
111  HistAxis("Vertex in Z Direction (cm)", Binning::Simple( 64, 0.0, 1280.0), kMuStartZ)
112  };
113 
114 
115  std::vector<std::string> taxes_labels = {"vtxx", "vtxy", "vtxz"};
116 
117  const Var basicWeight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
118 
119  assert(taxes.size() == taxes_labels.size());
120 
121  SpectrumHandler * basicSpecs = new SpectrumHandler();
122  SpectrumHandler * xsecSysts = new SpectrumHandler();
123 
124  basicSpecs->SetLoaders(datasets, dataset_labels);
125  basicSpecs->SetSpillCuts(kStandardDQCuts && kTightBeamQualityCuts); // Since all loaders initially defined in this handler, this only needs to run once.
126 
127  // Set the multi-universe systematic spectrums to use the same nominal loader object as in basicSpecs
128  xsecSysts->SetLoader(basicSpecs->GetLoader("nominal"), "xsec");
129 
130  // Basic spectra with systematic datasets
131  basicSpecs->SetHistAxes(taxes, taxes_labels);
132  basicSpecs->SetCuts(cuts, cut_labels);
133  basicSpecs->SetWeights(basicWeight);
134  basicSpecs->SetShifts(kNoShift);
135 
136  // XSec Universe Spectrums
137  const unsigned int nUniverses = 100;
139  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
140  std::vector<std::string> genie_labels;
141  for (unsigned int i = 0; i < nUniverses; i++)
142  genie_labels.push_back(std::to_string(i));
143  xsecSysts->SetHistAxes(taxes, taxes_labels);
144  xsecSysts->SetCuts(cuts, cut_labels);
145  xsecSysts->SetWeights(basicWeight);
146  xsecSysts->SetShifts(genie_shifts, genie_labels);
147 
148  basicSpecs->CreateSpectrums();
149  xsecSysts->CreateSpectrums();
150 
151  basicSpecs->Go();
152  xsecSysts->Go();
153 
154  basicSpecs->SaveSpectrums(fout);
155  xsecSysts->SaveSpectrums(fout);
156 
157  fout->Close();
158 
159  return;
160 }//end of plot_data_mc
161 
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
const Cut kNumuTightContainND([](const caf::SRProxy *sr){if(sr->vtx.nelastic< 1) return false;int ibesttrk=kBestTrack(sr);for(unsigned int i=0;i< sr->vtx.elastic[0].fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic[0].fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(int i=0;i< int(sr->trk.kalman.ntracks);++i){if(i==ibesttrk) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275|| sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;return((sr->trk.kalman.tracks[ibesttrk].stop.Z()< 1275 ||sr->trk.kalman.tracks[ibesttrk].trkyposattrans< 55) &&sr->trk.kalman.tracks[ibesttrk].trkfwdcellnd > 5 &&sr->trk.kalman.tracks[ibesttrk].trkbakcellnd > 10);})
Definition: NumuCCIncCuts.h:29
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
bool SaveSpectrums(TFile *f)
const Cut kNumuMyQuality([](const caf::SRProxy *sr){return(sr->trk.kalman.ntracks > 0 && sr->slc.nhit > 20 && sr->slc.ncontplanes > 4);})
Definition: NumuCCIncCuts.h:25
bool SetShifts(const std::vector< SystShifts > s, const std::vector< std::string > sl)
string filename
Definition: shutoffs.py:106
std::vector< std::string > cut_labels
const SpillCut kTightBeamQualityCuts([](const caf::SRSpillProxy *s){if(s->ismc) return true; if(s->trigger==2) return true;if(s->spilltimesec==0 &&s->deltaspilltimensec==0 &&s->widthx==0) return bool(s->isgoodspill);if(std::abs(s->deltaspilltimensec) > 0.5e9) return false;if(s->spillpot< 2e12) return false;if(s->hornI< -202|| s->hornI >-198) return false;if(s->posx< -2.00|| s->posx >+2.00) return false;if(s->posy< -2.00|| s->posy >+2.00) return false;return kBeamWidthCut(s);})
Definition: SpillCuts.h:10
const Var kMuStartX([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;return sr->trk.kalman.tracks[ibesttrk].start.X();})
Definition: NumuCCIncVars.h:47
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
void vertex_optimize(std::string filename="fVertexOpt.root")
const Cut kMuonIDCut
bool SetHistAxes(const std::vector< HistAxis > h, const std::vector< std::string > hl)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
assert(nhit_max >=nhit_nbins)
bool SetLoader(const std::string ds, const std::string dslab="")
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Var kMuStartZ([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;return sr->trk.kalman.tracks[ibesttrk].start.Z();})
Definition: NumuCCIncVars.h:49
const Var kMuStartY([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;return sr->trk.kalman.tracks[ibesttrk].start.Y();})
Definition: NumuCCIncVars.h:48
const Cut kTrueDetector
bool SetCuts(const std::vector< Cut > c, const std::vector< std::string > cl)
std::vector< Dataset > datasets
Definition: Periods.h:3
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
bool SetWeights(const std::vector< Var > w, const std::vector< std::string > wl)
enum BeamMode string