PredictionAddRock.cxx
Go to the documentation of this file.
3 
7 
8 #include "TDirectory.h"
9 #include "TObjString.h"
10 
11 #include "TH1.h"
12 #include "TVectorD.h"
13 
14 namespace ana
15 {
17 
18  //----------------------------------------------------------------------
20  PredictionNoExtrap *predRock,
21  double scaleNonSwap,
22  double scaleFluxSwap)
23  : fPredFid(predFid),
24  fPredRock(predRock),
25  fScaleNonSwap(scaleNonSwap),
26  fScaleFluxSwap(scaleFluxSwap)
27  {
28  }
29 
30  //----------------------------------------------------------------------
32  {
33  delete fPredFid;
34  delete fPredRock;
35  }
36 
37  //----------------------------------------------------------------------
39  {
40  return PredictComponent(calc,
43  Sign::kBoth);
44  }
45 
46  //----------------------------------------------------------------------
48  const SystShifts& syst) const
49  {
50  return PredictComponentSyst(calc, syst,
53  Sign::kBoth);
54  }
55 
56 
57  //----------------------------------------------------------------------
59  Flavors::Flavors_t flav,
61  Sign::Sign_t sign) const
62  {
63  return PredictComponentSyst(calc,kNoShift,flav,curr,sign);
64  }
65 
66  //----------------------------------------------------------------------
68  const SystShifts& syst,
69  Flavors::Flavors_t flav,
71  Sign::Sign_t sign) const
72  {
73  // Clone the systs
74  SystShifts systClone = syst;
75 
76  // First Check for RockScaleSyst
77  double systScale = 1 + systClone.GetShift(&kRockScaleSyst);
78  if(systScale < 0) systScale = 0;
79 
80  // IPredictions don't support predict component syst
81  systClone.SetShift(&kRockScaleSyst,0);
82 
83  // Strategy: assume PredictComponentSyst is properly implemented for the
84  // fiducial prediction, so start with that. Treat the rock as a series of
85  // delta's to that distribution, all predicted without syst shifts.
86  Spectrum ret = fPredFid->PredictComponentSyst(calc,systClone,
87  flav,curr,sign);
88 
89  Spectrum rock = RockComponent(calc, flav, curr, sign);
90  rock.Scale(systScale);
91  ret += rock;
92 
93  return ret;
94  }
95 
96  //----------------------------------------------------------------------
98  Flavors::Flavors_t flav,
100  Sign::Sign_t sign) const
101  {
102  Spectrum ret(0, {}, {}, 0, 0);
103 
104  if(curr & Current::kNC){
105  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
106 
107  // Predict the NC part
108  const double ncScale = (fScaleFluxSwap * fScaleNonSwap) /
110  ret = fPredRock->PredictComponent(calc, flav, Current::kNC, sign);
111  ret.Scale(ncScale);
112  }
113  else{
114  // Otherwise, still need to calculate something to get a suitably-binned
115  // Spectrum, but shouldn't keep in the total.
116  ret = fPredRock->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth);
117  ret.Clear();
118  }
119 
120 
121  if(curr & Current::kCC){
122  if (flav & Flavors::kNuEToNuE){
123  Spectrum delta = fPredRock->PredictComponent(calc, Flavors::kNuEToNuE,
124  Current::kCC, sign);
125  delta.Scale(fScaleNonSwap);
126  ret += delta;
127  }
128  if (flav & Flavors::kNuEToNuMu){
129  Spectrum delta = fPredRock->PredictComponent(calc, Flavors::kNuEToNuMu,
130  Current::kCC, sign);
131  delta.Scale(fScaleFluxSwap);
132  ret += delta;
133  }
134  // If you're astute, you'll realize we don't do anything for tauswap.
135  // We never had tauswap rock files made, so it doesn't make sense to sum
136  // them in.
137 
138  if (flav & Flavors::kNuMuToNuE){
139  Spectrum delta = fPredRock->PredictComponent(calc, Flavors::kNuMuToNuE,
140  Current::kCC, sign);
141  delta.Scale(fScaleFluxSwap);
142  ret += delta;
143  }
144  if (flav & Flavors::kNuMuToNuMu){
145  Spectrum delta = fPredRock->PredictComponent(calc,Flavors::kNuMuToNuMu,
146  Current::kCC, sign);
147  delta.Scale(fScaleNonSwap);
148  ret += delta;
149  }
150  }
151 
152  return ret;
153  }
154 
155  //----------------------------------------------------------------------
158  const SystShifts& shift,
159  double pot,
160  std::unordered_map<const ISyst*, std::vector<double>>& dchi) const
161  {
162  if(dchi.count(&kRockScaleSyst) &&
163  shift.GetShift(&kRockScaleSyst) > -1){
164  // Fiducial part doesn't know about the rock systematic
165  dchi.erase(dchi.find(&kRockScaleSyst));
166  SystShifts systClone = shift;
167  systClone.SetShift(&kRockScaleSyst, 0);
168  fPredFid->Derivative(calc, systClone, pot, dchi);
169 
170  // Add that back in manually
171  TH1D* hrock = RockComponent(calc, Flavors::kAll, Current::kBoth, Sign::kBoth).ToTH1(pot);
172  std::vector<double>& dchiRock = dchi[&kRockScaleSyst];
173  const int N = hrock->GetNbinsX()+2;
174  dchiRock.resize(N);
175  for(int i = 0; i < N; ++i)
176  dchiRock[i] = hrock->GetBinContent(i);
177  HistCache::Delete(hrock);
178  }
179  else{
180  // Otherwise we can just shortcut
181  fPredFid->Derivative(calc, shift, pot, dchi);
182  }
183  }
184 
185  //----------------------------------------------------------------------
187  {
188  // Start with the fiducial prediction, then add rock later
190 
191  // Now, we need to add in the rock estimate.
192  // Here, we need to scale down rock components because at generation time,
193  // we keep ~10% of spills - those that leave an FLSHit in the detector
194  // There are different rock repition rates for swap / nonswap files.
195  // If the neutrino changes flavor, scale by the fluxswap scale, otherwise
196  // use the nonswap scale
197  if (from == to){
199  delta.Scale(fScaleFluxSwap);
200  ret += delta;
201  }
202  else{
204  delta.Scale(fScaleNonSwap);
205  ret += delta;
206  }
207 
208  // If from / to are an invalid pair, trust that fPredFid and fPredRock
209  // will properly assert.
210  return ret;
211  }
212 
213  //----------------------------------------------------------------------
214  //nc
216  {
218 
220  double scale = (fScaleFluxSwap*fScaleNonSwap) /
222  Rock.Scale(scale); // NC's come from both swap and fluxswap files
223 
224  return Fid + Rock;
225  }
226 
228  {
229  Spectrum Fid = fPredFid->ComponentNC();
230 
231  Spectrum Rock = fPredRock->ComponentNC();
232  double scale = (fScaleFluxSwap*fScaleNonSwap) /
234  Rock.Scale(scale); // NC's come from both swap and fluxswap files
235 
236  return Fid + Rock;
237  }
238 
240  {
242 
244  double scale = (fScaleFluxSwap*fScaleNonSwap) /
246  Rock.Scale(scale); // NC's come from both swap and fluxswap files
247 
248  return Fid + Rock;
249  }
250  //end nc
251  //----------------------------------------------------------------------
252  void PredictionAddRock::SaveTo(TDirectory* dir) const
253  {
254  TDirectory* tmp = gDirectory;
255 
256  dir->cd();
257 
258  TObjString("PredictionAddRock").Write("type");
259 
260  fPredFid->SaveTo(dir->mkdir("predFid"));
261  fPredRock->SaveTo(dir->mkdir("predRock"));
262 
263  TVectorD scales(2);
264  scales[0] = fScaleNonSwap;
265  scales[1] = fScaleFluxSwap;
266  scales.Write("fRockScales");
267 
268  tmp->cd();
269  }
270 
271  //----------------------------------------------------------------------
272  std::unique_ptr<PredictionAddRock> PredictionAddRock::LoadFrom(TDirectory* dir)
273  {
274  assert(dir->GetDirectory("predFid"));
275  assert(dir->GetDirectory("predRock"));
276  auto predFid = ana::LoadFrom<IPrediction>(dir->GetDirectory("predFid")).release();
277  auto predRock = ana::LoadFrom<PredictionNoExtrap>(dir->GetDirectory("predRock")).release();
278 
279  TVectorD* scales = (TVectorD*)dir->Get("fRockScales");
280  double scaleNonSwap = (*scales)[0];
281  double scaleFluxSwap = (*scales)[1];
282 
283 
284  return std::unique_ptr<PredictionAddRock>(new PredictionAddRock(predFid,predRock,scaleNonSwap,scaleFluxSwap));
285  }
286 }
OscillatableSpectrum ComponentCC(int from, int to) const override
virtual void Derivative(osc::IOscCalculator *calc, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dchi) const override
virtual Spectrum ComponentNC() const
Definition: IPrediction.h:111
Spectrum ComponentNC() const override
virtual void SaveTo(TDirectory *dir) const override
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
Spectrum ComponentNCTotal() const override
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
Spectrum ComponentNCTotal() const override
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
double delta
Definition: runWimpSim.h:98
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
virtual void SaveTo(TDirectory *dir) const
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
virtual void Derivative(osc::IOscCalculator *calc, const SystShifts &shift, double pot, std::unordered_map< const ISyst *, std::vector< double >> &dchi) const
Definition: IPrediction.h:92
Float_t tmp
Definition: plot.C:36
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:734
virtual Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
static std::unique_ptr< PredictionAddRock > LoadFrom(TDirectory *dir)
virtual Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Spectrum ComponentNCAnti() const override
Double_t scale
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
virtual Spectrum ComponentNCAnti() const
Definition: IPrediction.h:113
OscillatableSpectrum ComponentCC(int from, int to) const override
#define pot
T GetShift(const ISyst *syst) const
virtual void SaveTo(TDirectory *dir) const override
def sign(x)
Definition: canMan.py:204
const SystShifts kNoShift
Definition: SystShifts.h:112
(&#39; survival&#39;)
Definition: IPrediction.h:19
virtual Spectrum ComponentNCTotal() const
Definition: IPrediction.h:109
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
PredictionNoExtrap * fPredRock
TDirectory * dir
Definition: macro.C:5
virtual OscillatableSpectrum ComponentCC(int from, int to) const
Definition: IPrediction.h:103
string syst
Definition: plotSysts.py:176
Neutral-current interactions.
Definition: IPrediction.h:40
Spectrum ComponentNCAnti() const override
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Prediction that just uses FD MC, with no extrapolation.
Spectrum RockComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
virtual Spectrum Predict(osc::IOscCalculator *calc) const override
Spectrum with true energy information, allowing it to be oscillated
void Scale(double scale) const
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const override
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
void rock(std::string suffix="full")
Definition: rock.C:28
Spectrum ComponentNC() const override
const DummyRockScaleSyst kRockScaleSyst