Public Member Functions | Protected Member Functions | Private Attributes | List of all members
ana::NeutronVisEScaleSyst2018 Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-02-24/3FlavorAna/Systs/EnergySysts2018.h"

Inheritance diagram for ana::NeutronVisEScaleSyst2018:
ana::EvtRNGSyst ana::ISyst

Public Member Functions

 NeutronVisEScaleSyst2018 (bool useProngs, double threshE=0.040, double sf=3.6, double prob=0.33)
 
void Shift (double sigma, caf::SRProxy *sr, double &weight) const override
 Perform the systematic shift. More...
 
virtual const std::stringShortName () const final
 The name printed out to the screen. More...
 
virtual const std::stringLatexName () const final
 The name used on plots (ROOT's TLatex syntax) More...
 
virtual void TruthShift (double sigma, caf::SRNeutrinoProxy *nu, double &weight) const
 
virtual bool IsGenieReweight () const
 GENIE reweights can only provide +/-1,2sigma. More...
 

Protected Member Functions

TRandom * RNG (const caf::SRProxy *sr) const
 Get the random number generator. More...
 

Private Attributes

bool fUseProngs
 
double fThreshE
 
double fScale
 
double fProb
 

Detailed Description

Definition at line 9 of file EnergySysts2018.h.

Constructor & Destructor Documentation

ana::NeutronVisEScaleSyst2018::NeutronVisEScaleSyst2018 ( bool  useProngs,
double  threshE = 0.040,
double  sf = 3.6,
double  prob = 0.33 
)
inline

Definition at line 12 of file EnergySysts2018.h.

References Shift(), sigma(), sr, and ana::weight.

13  : EvtRNGSyst(std::string("NeutronEvis") + (useProngs ? "Prongs" : "Primaries") + "Syst2018", "Neutron visible energy systematic 2018"),
14  fUseProngs(useProngs), fThreshE(threshE), fScale(sf), fProb(prob)
15  {}
EvtRNGSyst(const std::string &shortName, const std::string &latexName)
Definition: EvtRNGSyst.cxx:12
enum BeamMode string

Member Function Documentation

virtual bool ana::ISyst::IsGenieReweight ( ) const
inlinevirtualinherited

GENIE reweights can only provide +/-1,2sigma.

Reimplemented in ana::SummedSyst.

Definition at line 56 of file ISyst.h.

56 {return false;}
virtual const std::string& ana::ISyst::LatexName ( ) const
inlinefinalvirtualinherited

The name used on plots (ROOT's TLatex syntax)

Definition at line 30 of file ISyst.h.

References ana::ISyst::fLatexName, ana::ISyst::Shift(), sigma(), sr, and ana::weight.

Referenced by ana::PredictionInterp::DebugPlotColz(), GetGENIEShiftLabels(), ana::NuISyst::SaveTo(), SystsGENIEAna(), and WriteSystName().

30 {return fLatexName;}
std::string fLatexName
Definition: ISyst.h:60
TRandom * ana::EvtRNGSyst::RNG ( const caf::SRProxy sr) const
protectedinherited

Get the random number generator.

Definition at line 17 of file EvtRNGSyst.cxx.

References ana::EvtRNGSyst::fLastSeed, ana::EvtRNGSyst::fRng, samweb_client.utility::hash, and ana::Hash().

Referenced by Shift().

18  {
19  auto hash = Hash(sr);
20  if (fLastSeed != hash)
21  {
22  fRng->SetSeed(hash);
23  fLastSeed = hash;
24  }
25 
26  return fRng.get();
27  }
std::unique_ptr< TRandom > fRng
Definition: EvtRNGSyst.h:27
std::size_t fLastSeed
Definition: EvtRNGSyst.h:26
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:460
void ana::NeutronVisEScaleSyst2018::Shift ( double  sigma,
caf::SRProxy sr,
double &  weight 
) const
overridevirtual

Perform the systematic shift.

Override this function if your systematic depends on non-SRNeutrino quantities. If it is SRNeutrino-only, implement the other function, and let this default forward to you when necessary.

Parameters
sigmaNumber of sigma to shift record by
srThe record to inspect and alter
weightScale this weight for reweighting systematics

Reimplemented from ana::ISyst.

Definition at line 15 of file EnergySysts2018.cxx.

References caf::Proxy< caf::SRSlice >::calE, E, caf::Proxy< caf::SRVertexBranch >::elastic, caf::Proxy< caf::StandardRecord >::energy, fProb, fScale, fThreshE, fUseProngs, caf::Proxy< caf::SRElastic >::fuzzyk, caf::Proxy< caf::SRNumuEnergy >::hadcalE, if(), std::isnan(), caf::Proxy< caf::SRElastic >::IsValid, caf::Proxy< caf::StandardRecord >::mc, caf::Proxy< caf::SRTruthBranch >::nnu, caf::Proxy< caf::SRTruthBranch >::nu, caf::Proxy< caf::SREnergyBranch >::numu, caf::Proxy< caf::SRFuzzyK >::png, caf::Proxy< caf::SRFuzzyK >::png2d, ana::EvtRNGSyst::RNG(), scale, caf::Proxy< caf::StandardRecord >::slc, and caf::Proxy< caf::StandardRecord >::vtx.

Referenced by NeutronVisEScaleSyst2018().

16  {
17  if (sr->mc.nnu != 1)
18  return;
19 
20  // only do anything if at least 1 MeV of neutron-caused visE
21  if (sr->mc.nu[0].visEinslcNeutronBirks < 0.001)
22  return;
23 
24  // todo: what should sigma apply to?
25  // here I assume it's the scale factor, but it could reasonably be
26  // any of the three parameters, or some combination of all of them
27  double scale = 1 + sigma * (fScale - 1);
28 
29  std::vector<caf::SRProngProxy*> prongs;
30  if (sr->vtx.elastic.IsValid){
31  for (caf::SRFuzzyKProngProxy & png : sr->vtx.elastic.fuzzyk.png)
32  prongs.push_back(static_cast<caf::SRProngProxy*>(&png)); // this is a slicing operation, but we don't need the 3D part of it below anyway
33  for (caf::SRProngProxy & png : sr->vtx.elastic.fuzzyk.png2d)
34  prongs.push_back(&png);
35  }
36  double addedE = 0;
37  if (fUseProngs){
38  for (auto & png : prongs){
39  // is this right? do we want to search the whole mother list, or just the immediate mother? hmm.
40  // for now, no choice -- motherlist is filled with 0s until ART 2.11 comes in with its fix to ROOT
41  // bool neutronFound = false;
42  //// std::cout << " prong mother list size = " << png->truth.motherlist.size() << std::endl;
43  // for (const auto & m : png->truth.motherlist)
44  // {
45  //// std::cout << " ancestor pdg = " << m << std::endl;
46  // if (m == 2112)
47  // {
48  // neutronFound = true;
49  // break;
50  // }
51  // }
52  // if (!neutronFound)
53  // continue;
54  if (png->truth.motherpdg != 2112)
55  continue;
56 
57 
58  // yes: adjust all neutrons below some threshold.
59  // skip the ones above it.
60  double visE = png->truth.visEinslc; //visEinslcBirks;
61  if (visE > fThreshE)
62  continue;
63 
64  double moreE = (RNG(sr)->Uniform() < fProb) ? visE * (scale - 1) : 0;
65  addedE += moreE;
66 
67  png->calE += moreE;
68  if (png->calE < 0)
69  png->calE = 0;
70  png->weightedCalE += moreE;
71  if (png->weightedCalE < 0)
72  png->weightedCalE = 0;
73  } // for (png)
74  } // if (fUseProngs)
75  else{
76  std::map<int, double> addedEbyTrkID;
77  for (auto & prim : sr->mc.nu[0].prim){
78  if (prim.pdg != 2112)
79  continue;
80 
81  // neutrons won't leave any energy directly (probably) but best be safe
82  double visE = prim.visEinslcBirks + prim.daughterVisEinslcBirks;
83  // concession to rapidly changing CAF structure
84  if (std::isnan(visE))
85  visE = prim.visEinslc + prim.daughterVisEinslc;
86  if (visE > fThreshE)
87  continue;
88 
89  double moreE = (RNG(sr)->Uniform() < fProb) ? visE * (scale - 1) : 0;
90  addedE += moreE;
91  addedEbyTrkID[prim.trkID] += moreE;
92  }
93 
94  // try to redistribute energy to the prongs.
95  // first: find any prongs that are directly associated with the primary neutrons
96  // that we tagged above.
97  std::vector<caf::SRProngProxy*> otherNeutronProngs; // for the next loop
98  for (auto & png : prongs){
99  int motherTrk = 0;
100  for (const auto & mother : png->truth.motherlist){
101  if (addedEbyTrkID.find(mother) != addedEbyTrkID.end()){
102  motherTrk = mother;
103  break;
104  }
105  }
106  if (!motherTrk){
107  // save other neutron prongs for a second pass.
108  // note: threshold originally CAME from calE, so better respect it here!
109  if (png->truth.motherpdg == 2112 && png->calE <= fThreshE)
110  otherNeutronProngs.push_back(png);
111  continue;
112  }
113  // todo: what if there are multiple prongs that descend from this neutron??
114  double E = addedEbyTrkID.at(motherTrk);
115 
116  png->calE += E;
117  if (png->calE < 0)
118  png->calE = 0;
119  png->weightedCalE += E;
120  if (png->weightedCalE < 0)
121  png->weightedCalE = 0;
122 
123  addedEbyTrkID.erase(motherTrk);
124  }
125 
126  // now distribute any energy left to the other neutron-descendants...
127  // since we don't know which ones they "should" be,
128  // do something dumb and assign them "randomly"
129  for (auto & png : otherNeutronProngs){
130  if (addedEbyTrkID.empty())
131  break;
132 
133  png->calE += addedEbyTrkID.begin()->second;
134  if (png->calE < 0)
135  png->calE = 0;
136  png->weightedCalE += addedEbyTrkID.begin()->second;
137  if (png->weightedCalE < 0)
138  png->weightedCalE = 0;
139 
140  addedEbyTrkID.erase(addedEbyTrkID.begin());
141 // std::cout << " added energy vector size is now: " << addedEbyTrkID.size() << std::endl;
142  }
143 
144  }
145 
146  // add the neutron energy back to the energy estimator inputs.
147  // * for numu, we posit that the reco muon is never affected. so we only adjust Ehad.
148  // * for nue, we've already modified neutron-descended prongs above,
149  // but most of those won't be judged "electromagnetic" by CVN,
150  // so we care about the "hadronic" part, which is just slice calE - (sum of "EM" prong calEs).
151  // Thus modifying slice calE is good enough.
152  sr->energy.numu.hadcalE += addedE; // this goes into the numu energy estimator function
153  sr->slc.calE += addedE; // this guy is used by the nue energy estimator function
154  }
caf::Proxy< std::vector< int > > motherlist
Definition: SRProxy.h:639
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:214
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
Proxy for caf::SRFuzzyKProng.
Definition: SRProxy.h:2003
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
caf::Proxy< int > motherpdg
Definition: SRProxy.h:641
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
caf::Proxy< caf::SRParticleTruth > truth
Definition: SRProxy.h:1932
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
Double_t scale
Definition: plot.C:25
if(dump)
caf::Proxy< float > weightedCalE
Definition: SRProxy.h:1936
Float_t E
Definition: plot.C:20
caf::Proxy< float > hadcalE
Definition: SRProxy.h:169
double sigma(TH1F *hist, double percentile)
TRandom * RNG(const caf::SRProxy *sr) const
Get the random number generator.
Definition: EvtRNGSyst.cxx:17
Proxy for caf::SRProng.
Definition: SRProxy.h:1897
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
caf::Proxy< float > calE
Definition: SRProxy.h:1292
caf::Proxy< float > calE
Definition: SRProxy.h:1907
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
caf::Proxy< float > visEinslc
Definition: SRProxy.h:653
virtual const std::string& ana::ISyst::ShortName ( ) const
inlinefinalvirtualinherited
virtual void ana::ISyst::TruthShift ( double  sigma,
caf::SRNeutrinoProxy nu,
double &  weight 
) const
inlinevirtualinherited

For systematics that deal only with the neutrino truth and not any reconstruction/PID details. Systematics defined this way will work on nuTree-derived spectra too (e.g. denominators of efficiencies).

Reimplemented in ana::BeamSyst, demo::DemoSyst1, ana::GenericSystComponentScale< T >, ana::GenericSystComponentScale< T >, ana::GenericSystComponentScale< T >, ana::ReinteractionSyst, and ana::NOvARwgtSyst.

Definition at line 46 of file ISyst.h.

Referenced by ana::ISyst::Shift().

49  {
50  // Implement this function if your systematic depends only
51  // SRNeutrino. Left blank by default, since systematics using other
52  // information can do nothing sensible to the nuTree.
53  }

Member Data Documentation

double ana::NeutronVisEScaleSyst2018::fProb
private

Definition at line 23 of file EnergySysts2018.h.

Referenced by Shift().

double ana::NeutronVisEScaleSyst2018::fScale
private

Definition at line 22 of file EnergySysts2018.h.

Referenced by Shift().

double ana::NeutronVisEScaleSyst2018::fThreshE
private

Definition at line 21 of file EnergySysts2018.h.

Referenced by Shift().

bool ana::NeutronVisEScaleSyst2018::fUseProngs
private

Definition at line 20 of file EnergySysts2018.h.

Referenced by Shift().


The documentation for this class was generated from the following files: