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 using std::cout;
10 using std::endl;
11 using std::make_pair;
12 using std::make_unique;
13 using std::map;
14 using std::pair;
15 using std::string;
16 using std::unique_ptr;
17 
18 namespace ana {
19 
20  // -------------------------------------------------------------------------------------
21  /// NCSyst
22  NCSyst::NCSyst(string name, string label, string sampleName,
23  map<int, pair<TH1*, TH1*> > shifts)
24  : NuISyst(name, label, sampleName, shifts)
25  {}
26 
27  // -------------------------------------------------------------------------------------
28  /// Implementation of standard ISyst::Shift function for NCSyst
29  void NCSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
30 
31  // Find the energy estimate for this slice
32  double calE = kNus20Energy(sr);
33 
34  // Find the signal or background oscillation channel
36  if (sr->mc.nnu == 0)
37  chan = SystChannel::kBkg;
38  else if (!sr->mc.nu[0].iscc)
39  chan = SystChannel::kSig;
40  else if (sr->mc.nu[0].iscc)
41  chan = SystChannel::kBkg;
42  if (chan == SystChannel::kUnknown) {
43  cout << "Error: Unknown oscillation channel." << endl;
44  abort();
45  }
46 
47  weight *= WeightFor(chan, sigma, calE);
48 
49  } // function NCSyst::Shift
50 
51  // -------------------------------------------------------------------------------------
52  /// SaveTo implementation for NCSyst
53  void NCSyst::SaveTo(TDirectory* dir) const {
54 
55  TDirectory* tmp = gDirectory;
56  dir->cd();
57  TObjString("NCSyst").Write("type");
58  tmp->cd();
59  NuISyst::SaveTo(dir);
60 
61  } // function NCSyst::SaveTo
62 
63  // -------------------------------------------------------------------------------------
64  /// LoadFrom implementation for NCSyst
65  unique_ptr<NCSyst> NCSyst::LoadFrom(TDirectory* dir) {
66 
67  DontAddDirectory guard;
68 
69  // Check this is a NueSyst
70  TObjString* tag = (TObjString*)dir->Get("type");
71  assert(tag);
72  string sTag = tag->GetString().Data();
73  assert(sTag == "NCSyst");
74 
75  // Load names
76  TObjString* name = (TObjString*)dir->Get("name");
77  assert(name);
78  string sName = name->GetString().Data();
79  TObjString* label = (TObjString*)dir->Get("label");
80  assert(label);
81  string sLabel = label->GetString().Data();
82  TObjString* sampleName = (TObjString*)dir->Get("samplename");
83  assert(sampleName);
84  string sSampleName = sampleName->GetString().Data();
85 
86  map<int, pair<TH1*, TH1*>> shifts;
87 
88  // Load all shifts
89  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
90  TKey* key;
91  while ((key = (TKey*)next())) {
92  TDirectory* sigmaDir = (TDirectory*)key->ReadObj();
93  const TVectorD vSigma = *(TVectorD*)sigmaDir->Get("sigma");
94  TH1* hSignal = (TH1*)sigmaDir->Get("sig");
95  TH1* hBackground = (TH1*)sigmaDir->Get("bkg");
96  shifts[vSigma[0]] = make_pair(hSignal, hBackground);
97  }
98 
99  return make_unique<NCSyst>(sName, sLabel, sSampleName, shifts);
100 
101  } // function NumuSyst::SaveTo
102 
103 } // 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:22
const Var weight
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
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:574
virtual double WeightFor(SystChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
Definition: CovMxSysts.cxx:80
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
Float_t tmp
Definition: plot.C:36
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:65
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Implementation of standard ISyst::Shift function for NCSyst.
Definition: CovMxSysts.cxx:29
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
virtual void SaveTo(TDirectory *dir) const override
SaveTo implementation for NCSyst.
Definition: CovMxSysts.cxx:53
assert(nhit_max >=nhit_nbins)
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
virtual void SaveTo(TDirectory *dir) const
SaveTo implementation for NuISyst.
Definition: CovMxSysts.cxx:51
void next()
Definition: show_event.C:84