genie_syst.C
Go to the documentation of this file.
1 // Biswa B. uses tag release S16-08-24. Going to be in line with him.
2 // 2017/03/01 Jon requested Biswa to update the NDAna Vars and Cuts.
3 // Biswa commited changes the same day. Using development should be fime.
4 
5 #ifdef __CINT__
6 void genie_syst(
7  int ndatafiles = INT_MAX,
8  int nmcfiles = INT_MAX,
9  string dataflistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.b_nd_numi_fhc_period3_v1_goodruns.txt",
10  string mcflistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt"
11 )
12 {
13  std::cout << "Sorry, you must run in compiled mode" << std::endl;
14 }
15 #else
16 
17 #include "CAFAna/Core/HistAxis.h"
18 #include "CAFAna/Cuts/TruthCuts.h"
22 
23 using namespace ana;
24 
26  kRemidCut &&
28  kIsFiducial &&
30  kECut ;
31 
32 const Cut kMyIsNuMu(
33  [](const caf::SRProxy* sr)
34  {
35  return sr->mc.nu[0].pdg == 14;
36  });
37 
39  int ndatafiles = INT_MAX,
40  int nmcfiles = INT_MAX,
41  string dataflistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.b_nd_numi_fhc_period3_v1_goodruns.txt",
42  string mcflistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt"
43 )
44 {
45  /// load file names into vectors
46  vector<string> dataflist = LoadFileList(dataflistpn);
47  vector<string> mcflist = LoadFileList(mcflistpn);
48  if(!dataflist.size() || !mcflist.size())
49  {
50  cout << "Not able to load the file lists. Quit." << endl;
51  return;
52  }
53  /// only run over the specified number of files
54  dataflist.resize(min(min(ndatafiles, INT_MAX), (int)dataflist.size()));
55  mcflist.resize(min(min(nmcfiles, INT_MAX), (int)mcflist.size()));
56  cout << mcflist[0] << endl;
57 
58  /// loaders
59  SpectrumLoader dataloader(dataflist);
60  SpectrumLoader mcloader(mcflist);
61 
62  /// first background estimator try: TrivialBkgdEstimator
63  std::vector<Cut> trueBkgds;
64  //~ trueBkgds.push_back(!kIsNumuCC);
65  trueBkgds.push_back(!kMyIsNuMu);
66  TrivialBkgdEstimator bkgdest(
67  mcloader,
68  kRecoE_Custom_Axis,
69  kMyAllNumuCCCuts,
70  trueBkgds
71  );
72 
73  /// Dan's cross section analysis class
74  CrossSectionAnalysis xsecana(
75  mcloader,
76  dataloader,
77  kRecoE_Custom_Axis,
78  kTrueE_Custom_Axis,
79  kMyAllNumuCCCuts,
80  kMyIsNuMu,
81  &bkgdest,
82  TVector3(-176,-172,25),
83  TVector3(177,179,1150),
84  2
85  );
86 
87  /// save results
88  //~ TFile ouf("xsec_output.root","recreate");
89  //~ xsecana.SaveTo(gDirectory);
90  //~ ouf.Close();
91  dataloader.Go();
92  cout << ">>> data loader done" << endl;
93  mcloader.Go();
94  cout << ">>> MC loader done" << endl;
95 
96  TH1* h = xsecana.Result();
97  h->Draw();
98 }
99 
100 #endif
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kMyIsNuMu([](const caf::SRProxy *sr){return sr->mc.nu[0].pdg==14;})
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
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
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
void genie_syst(int ndatafiles=INT_MAX, int nmcfiles=INT_MAX, string dataflistpn="/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.b_nd_numi_fhc_period3_v1_goodruns.txt", string mcflistpn="/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt")
Definition: genie_syst.C:38
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
Definition: Utilities.cxx:210
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const Cut kIsFiducial([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return false;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;bool isfid=(sr->trk.kalman.tracks[ibesttrk].start.X()< vtxmax.X()&&sr->trk.kalman.tracks[ibesttrk].start.X() > vtxmin.X()&&sr->trk.kalman.tracks[ibesttrk].start.Y() > vtxmin.Y()&&sr->trk.kalman.tracks[ibesttrk].start.Y()< vtxmax.Y()&&sr->trk.kalman.tracks[ibesttrk].start.Z() > vtxmin.Z()&&sr->trk.kalman.tracks[ibesttrk].start.Z()< vtxmax.Z());return isfid;})
Definition: NumuCCIncCuts.h:27
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
const Cut kRemidCut
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Generic organizational class for a cross section analysis.
Just return the MC prediction for the background.
const Cut kMyAllNumuCCCuts
Definition: genie_syst.C:25