EnergySysts2018.cxx
Go to the documentation of this file.
1 #include <map>
2 
3 #include "TRandom3.h"
4 
7 
8 namespace ana
9 {
10  const NeutronVisEScaleSyst2018 kNeutronVisEScaleProngsSyst2018(true);
11  const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false);
12 
13  //----------------------------------------------------------------------
15  Shift(double sigma, caf::SRProxy* sr, double& weight) const
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  }
155 }
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2018
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2003
const Var weight
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:213
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2095
const NeutronVisEScaleSyst2018 kNeutronVisEScaleProngsSyst2018(true)
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2077
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2002
Double_t scale
Definition: plot.C:25
if(dump)
Float_t E
Definition: plot.C:20
caf::StandardRecord * sr
caf::Proxy< float > hadcalE
Definition: SRProxy.h:168
double sigma(TH1F *hist, double percentile)
TRandom * RNG(const caf::SRProxy *sr) const
Get the random number generator.
Definition: EvtRNGSyst.cxx:17
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2017
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2101
caf::Proxy< float > calE
Definition: SRProxy.h:1248
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2105