muonid_opt.C
Go to the documentation of this file.
3 #include "CAFAna/Core/HistAxis.h"
4 #include "CAFAna/Vars/Vars.h"
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Cut.h"
7 #include "CAFAna/Core/MultiVar.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 
16 
19 
20 #include "TFile.h"
21 #include "CAFAna/Core/Spectrum.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 
26 
28 
33 
34 #include "CAFAna/Analysis/Plots.h"
35 #include "TLegend.h"
36 #include "TCanvas.h"
37 #include "TColor.h"
38 
39 #include <stdio.h>
40 #include <string.h>
41 
42 using namespace ana;
43 
44 void muonid_opt(std::string fname = "fMuonIDOpt.root")
45 {
46  TFile * fout = new TFile(fname.c_str(),"RECREATE");
47  if (!fout || fout == NULL || fout->IsZombie())
48  throw runtime_error("Could not open output file.");
49 
50  std::vector<std::string> dataset_labels = {"nominal", "lightdown", "lightup", "ckv", "calibneg", "calibpos", "calibshape"};
51 
52  std::vector<ana::SpectrumLoaderBase*> loaders;
53 
54  // Nominal removes muonid_training. Others are systematic shifts.
55  const vector<string> datasets =
56  { "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",
57  "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",
58  "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",
59  "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",
60  "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",
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-pos-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-func_v1"
63  };
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  const std::vector<Var> vars = {
72  kMuonID
73  };
74 
75  const std::vector<std::string> vars_labels = {"muonid"};
76 
77  const std::vector<string> vars_titles = {
78  "MuonID"
79  };
80 
81  const std::vector<Binning> vars_bins = {
82  Binning::Simple(100, -1, 1)
83  };
84 
85  // Preselection is quality, true detector fiducial, and tight containment cut.
86  const Cut kPreselection = kNumuMyQuality && kIsFiducialLoose && kNumuTightContainND;
87  const Cut kSignalCut = kIsNumuCCCut && kTrueFiducialLoose;
88 
89  std::vector<string> cuts_labels = {
90  "selection",
91  "sel_signal",
92  "sel_bkgd",
93  "all_signal",
94  "all_bkgd",
95  "all_numucc",
96  "all_notnumucc"
97  };
98 
99  std::vector<Cut> cuts = {
100  kPreselection, // Preselection
101  kPreselection && kSignalCut, // True && Preselection
102  kPreselection && !kSignalCut, // Bkgd && Preselection
103  kSignalCut, // All true signal
104  !kSignalCut, // All true bkgd
105  kIsNumuCCCut && kTrueDetector,// All numucc (inc. outside fiducial)
106  !kIsNumuCCCut && kTrueDetector// All !numucc (inc. outside fiducial)
107  };
108 
109  // Split truth by interaction mode as well
110  bool splitByInteractionMode = true;
111  if (splitByInteractionMode){
112  std::map<std::string, Cut> interaction_mode_cuts = {
113  {"_qe", kIsQE},
114  {"res", kIsRes},
115  {"dis", kIsDIS},
116  {"mec", kIsDytmanMEC},
117  {"oth", !kIsQE && !kIsRes && !kIsDIS && !kIsDytmanMEC}
118  };
119  for (auto intcut : interaction_mode_cuts){
120  cuts_labels.push_back("sel_sig_" + intcut.first);
121  cuts.push_back(kPreselection && kSignalCut && intcut.second);
122 
123  cuts_labels.push_back("all_sig_" + intcut.first);
124  cuts.push_back(kSignalCut && intcut.second);
125 
126  cuts_labels.push_back("all_numucc_" + intcut.first);
127  cuts.push_back(kIsNumuCCCut && kTrueDetector && intcut.second);
128  }
129  }
130 
131  assert(cuts.size() == cuts_labels.size());
132 
133  const Var basicWeight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
134 
135  SpectrumHandler * basicSpecs = new SpectrumHandler();
136  SpectrumHandler * xsecSysts = new SpectrumHandler();
137 
138  basicSpecs->SetLoaders(datasets, dataset_labels);
139  basicSpecs->SetSpillCuts(kStandardDQCuts && kTightBeamQualityCuts); // Since all loaders initially defined in this handler, this only needs to run once.
140 
141  // Set the multi-universe systematic spectrums to use the same nominal loader object as in basicSpecs
142  xsecSysts->SetLoader(basicSpecs->GetLoader("nominal"), "xsec");
143 
144  // Basic spectra with systematic datasets
145  basicSpecs->SetVars(vars, vars_labels, vars_titles, vars_bins);
146  basicSpecs->SetCuts(cuts, cuts_labels);
147  basicSpecs->SetWeights(basicWeight);
148  basicSpecs->SetShifts(kNoShift);
149 
150  // XSec Universe Spectrums
152  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
153  std::vector<std::string> genie_labels;
154  for (unsigned int i = 0; i < 1000; i++)
155  genie_labels.push_back(std::to_string(i));
156  xsecSysts->SetVars(vars, vars_labels, vars_titles, vars_bins);
157  xsecSysts->SetCuts(cuts, cuts_labels);
158  xsecSysts->SetWeights(basicWeight);
159  xsecSysts->SetShifts(genie_shifts, genie_labels);
160 
161  basicSpecs->CreateSpectrums();
162  xsecSysts->CreateSpectrums();
163 
164  basicSpecs->Go();
165  xsecSysts->Go();
166 
167  basicSpecs->SaveSpectrums(fout);
168  xsecSysts->SaveSpectrums(fout);
169 
170  fout->Close();
171 
172  return;
173 
174 }//end of void
const Cut kIsQE
Definition: TruthCuts.cxx:104
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
SpectrumLoaderBase * GetLoader(const std::string label)
const Cut kIsRes
Definition: TruthCuts.cxx:111
bool SaveSpectrums(TFile *f)
bool SetVars(const std::vector< Var > v, const std::vector< std::string > vl, const std::vector< std::string > vt, const std::vector< Binning > vb)
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)
const Cut kIsDIS
Definition: TruthCuts.cxx:118
bool SetSpillCuts(const std::vector< SpillCut > spillcuts)
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
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 kMuonID
Definition: NDXSecMuonPID.h:32
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
const Cut kPreselection
const std::map< std::pair< std::string, std::string >, Variable > vars
const Cut kIsNumuCCCut
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
bool SetLoaders(const std::vector< std::string > &ds, const std::vector< std::string > &dsnames)
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)
void muonid_opt(string outFilename="muonid_opt_hists.root")
Definition: muonid_opt.C:31
bool SetLoader(const std::string ds, const std::string dslab="")
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
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