CovMxSysts.cxx
Go to the documentation of this file.
2 #include "NuXAna/Vars/NusVars.h"
3 
5 
6 #include "TKey.h"
7 #include "TVectorD.h"
8 
9 #include <iostream>
10 
11 using std::cout;
12 using std::endl;
13 using std::make_pair;
14 using std::make_unique;
15 using std::map;
16 using std::pair;
17 using std::string;
18 using std::unique_ptr;
19 
20 namespace ana {
21 
22  // -------------------------------------------------------------------------------------
23  /// NCSyst
24  NCSyst::NCSyst(string name, string label, string sampleName,
25  map<int, pair<TH1*, TH1*> > shifts)
26  : NuISyst(name, label, sampleName, shifts)
27  {}
28 
29  // -------------------------------------------------------------------------------------
30  /// Implementation of standard ISyst::Shift function for NCSyst
31  void NCSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
32 
33  // Find the energy estimate for this slice
34  double calE = kNus20Energy(sr);
35 
36  // Find the signal or background oscillation channel
38  if (sr->mc.nnu == 0)
39  chan = SystChannel::kBkg;
40  else if (!sr->mc.nu[0].iscc)
41  chan = SystChannel::kSig;
42  else if (sr->mc.nu[0].iscc)
43  chan = SystChannel::kBkg;
44  if (chan == SystChannel::kUnknown) {
45  cout << "Error: Unknown oscillation channel." << endl;
46  abort();
47  }
48 
49  weight *= WeightFor(chan, sigma, calE);
50 
51  } // function NCSyst::Shift
52 
53  // // -------------------------------------------------------------------------------------
54  // /// SaveTo implementation for NCSyst
55  // void NCSyst::SaveTo(TDirectory* dir) const {
56 
57  // TDirectory* tmp = gDirectory;
58  // dir->cd();
59  // TObjString("NCSyst").Write("type");
60  // tmp->cd();
61  // NuISyst::SaveTo(dir);
62 
63  // } // function NCSyst::SaveTo
64 
65  // -------------------------------------------------------------------------------------
66  /// LoadFrom implementation for NCSyst
67  unique_ptr<NCSyst> NCSyst::LoadFrom(TDirectory* dir) {
68 
69  DontAddDirectory guard;
70 
71  // Check this is a NueSyst
72  TObjString* tag = (TObjString*)dir->Get("type");
73  assert(tag);
74  string sTag = tag->GetString().Data();
75  assert(sTag == "NCSyst");
76 
77  // Load names
78  TObjString* name = (TObjString*)dir->Get("name");
79  assert(name);
80  string sName = name->GetString().Data();
81  TObjString* label = (TObjString*)dir->Get("label");
82  assert(label);
83  string sLabel = label->GetString().Data();
84  TObjString* sampleName = (TObjString*)dir->Get("samplename");
85  assert(sampleName);
86  string sSampleName = sampleName->GetString().Data();
87 
88  map<int, pair<TH1*, TH1*>> shifts;
89 
90  // Load all shifts
91  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
92  TKey* key;
93  while ((key = (TKey*)next())) {
94  TDirectory* sigmaDir = (TDirectory*)key->ReadObj();
95  const TVectorD vSigma = *(TVectorD*)sigmaDir->Get("sigma");
96  TH1* hSignal = (TH1*)sigmaDir->Get("sig");
97  TH1* hBackground = (TH1*)sigmaDir->Get("bkg");
98  shifts[vSigma[0]] = make_pair(hSignal, hBackground);
99  }
100 
101  return make_unique<NCSyst>(sName, sLabel, sSampleName, shifts);
102 
103  } // function NumuSyst::LoadFrom
104 
105 } // namespace ana
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
NCSyst(std::string name, std::string label, std::string sample_name="", std::map< int, std::pair< TH1 *, TH1 * > > shifts={})
NCSyst.
Definition: CovMxSysts.cxx:24
const Var weight
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Var kNus20Energy([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kNEARDET &&sr->hdr.det!=caf::kFARDET) return(double) sr->slc.calE;double pars[4][6]={{1.049, 0.795, 0.8409, 0.17, 0.82,-1.00},{1.025, 0.797, 0.9162, 0.53,-0.26,-1.48},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00}};int detCur=(sr->hdr.det==caf::kFARDET)+((sr->spill.isRHC)<< 1);double e_EM=ana::kEME_2020(sr);double e_Had=ana::kHADE_2020(sr);double e_Cal=sr->slc.calE;return(e_EM/pars[detCur][0]+e_Had/pars[detCur][1])/(pars[detCur][2]+pars[detCur][3]*std::pow(e_Cal+pars[detCur][4], 2)*std::exp(pars[detCur][5]*e_Cal));})
Definition: NusVars.h:64
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
virtual double WeightFor(SystChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
Definition: CovMxSysts.cxx:83
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const char * label
SystChannel
Definition: CovMxSysts.h:12
caf::StandardRecord * sr
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
static std::unique_ptr< NCSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NCSyst.
Definition: CovMxSysts.cxx:67
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Implementation of standard ISyst::Shift function for NCSyst.
Definition: CovMxSysts.cxx:31
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
void next()
Definition: show_event.C:84
enum BeamMode string