BeamSysts.cxx
Go to the documentation of this file.
2 
4 
6 
7 #include "TFile.h"
8 #include "TH1.h"
9 #include "TKey.h"
10 
11 namespace ana
12 {
13 
15  {
16  // if (sr->mc.nnu != 1)
17  // return 1.;
18 
19  return GetWeight(sr, kPlus);
20  }
21 
22  // Smoothened versions of PPFX Flux CV and Multiverse weights
24  {
25  if(UnivIdx >= 100){
26  std::cerr << "Error : Universe Index too large!" << std::endl;
27  abort();
28  }
29 
30  static NuTruthVar* Cache[100] = {0};
31  std::string UnivIdxStr = TString::Format("%02d", UnivIdx).Data();
32  if(!Cache[UnivIdx])
33  Cache[UnivIdx] = new NuTruthVar(BeamWeightFunc(FindCAFAnaDir()+"/data/beam/ppfx_smooth_multiverse_weights_2020.root",
34  "ppfx_univ"+UnivIdxStr));
35 
36  return *Cache[UnivIdx];
37  }
38 
39  BeamSyst* GetPPFXUnivSyst(int UnivIdx)
40  {
41  if(UnivIdx >= 100){
42  std::cerr << "Error : Universe Index too large!" << std::endl;
43  abort();
44  }
45 
46  static BeamSyst* Cache[100] = {0};
47  std::string UnivIdxStr = TString::Format("%02d", UnivIdx).Data();
48  if(!Cache[UnivIdx])
49  Cache[UnivIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_smooth_multiverse_weights_2020.root"),
50  "ppfx_univ"+UnivIdxStr, "PPFX Universe "+UnivIdxStr);
51 
52  return Cache[UnivIdx];
53 
54  }
55 
56 
57  const BeamSyst kPPFXCVSyst((FindCAFAnaDir()+"/data/beam/ppfx_smooth_multiverse_weights_2020.root"),
58  "ppfx_cv", "PPFX CV");
59 
60  // preliminary PC systematics for Hadron production
61  // docdb: 22532
62 
64  {
65  if(PCIdx >= 50){
66  std::cerr << "Error : Component Index too large!" << std::endl;
67  abort();
68  }
69 
70  static BeamSyst* Cache[50] = {0};
71  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
72  if(!Cache[PCIdx])
73  Cache[PCIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_pc_shifts_smooth_2017.root"),
74  "ppfx_pc"+PCIdxStr, "PPFX Universe "+PCIdxStr);
75 
76  return Cache[PCIdx];
77  }
78 
79  // PCA systs on full set of flux systematics including PPFX Hadron production and Beam Focussing
81  {
82  if(PCIdx >= 50){
83  std::cerr << "Error : Component Index too large!" << std::endl;
84  abort();
85  }
86 
87  static BeamSyst* Cache[50] = {0};
88  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
89  if(!Cache[PCIdx])
90  Cache[PCIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_hadp_beam_pc_shifts_fn_2018.root"),
91  "ppfx_hadp_beam_pc"+PCIdxStr, "Flux Component "+PCIdxStr);
92 
93  return Cache[PCIdx];
94  }
95 
96  // PCA systs on full set of flux systematics including PPFX Hadron production and Beam Focussing for 2020 using the nuTree, not recTree
98  {
99  if(PCIdx >= 50){
100  std::cerr << "Error : Component Index too large!" << std::endl;
101  abort();
102  }
103 
104  static BeamSyst* Cache[50] = {0};
105  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
106  if(!Cache[PCIdx])
107  Cache[PCIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_hadp_beam_pc_shifts_fn_2020_nutree.root"),
108  "ppfx_hadp_beam_pc"+PCIdxStr, "Flux Component "+PCIdxStr);
109 
110  return Cache[PCIdx];
111  }
112 
113  // PPFX systs for ND SBAna, uses diagonalisation in F and N basis instead of the F/N basis for the standard analysis
115  {
116  if(PCIdx >= 10){
117  std::cerr << "Error : Component Index too large!" << std::endl;
118  abort();
119  }
120 
121  static BeamSyst* Cache[10] = {0};
122  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
123  if(!Cache[PCIdx])
124  Cache[PCIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_pc_shifts_nd_sbana_2017.root"),
125  "ppfx_ndsb_ana_pc"+PCIdxStr, "Hadron Production ND Component "+PCIdxStr);
126 
127  return Cache[PCIdx];
128  }
129 
130  // PCA systs on full set of flux systematics including PPFX Hadron production and Beam Focussing using an ND-only covariance matrix
131  // May be useful for ND analyses
133  {
134  if(PCIdx >= 50){
135  std::cerr << "Error : Component Index too large!" << std::endl;
136  abort();
137  }
138 
139  static BeamSyst* Cache[50] = {0};
140  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
141  if(!Cache[PCIdx])
142  Cache[PCIdx] = new BeamSyst((FindCAFAnaDir()+"/data/beam/ppfx_hadp_beam_pc_shifts_nd_2018.root"),
143  "ppfx_hadp_beam_nd_pc"+PCIdxStr, "Flux Component "+PCIdxStr);
144 
145  return Cache[PCIdx];
146  }
147 
148  /// Horn Current +/-2kA
149  const BeamSyst
150  kBeamHornCurrent((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
151  "2kA",
152  "2kAHornCurrent","+/- 2kA Horn Current");
153 
154  ///Beam Spot Size 1.3 +/- 0.2 mm both XY
155  const BeamSyst
156  kBeamSpotSize((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
157  "0.2mmBeamSpotSize",
158  "0p2mmBeamSpotSize"," 1.3 +/- 0.2 mm Spot Size");
159 
160  ///Beam position on target +-1 mm, X/Y separately
161  const BeamSyst
162  kBeamPosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
163  "1mmBeamShiftX","Beam Position X");
164 
165  const BeamSyst
166  kBeamPosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
167  "1mmBeamShiftY","Beam Position Y");
168 
169  ///Horn 1 and 2 position +-3mm in X and Y separately
170  const BeamSyst
171  kBeamH1PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
172  "3mmHorn1X","+/- 3mm Horn 1 X Position");
173  const BeamSyst
174  kBeamH1PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
175  "3mmHorn1Y","+/- 3mm Horn 1 Y Position");
176  const BeamSyst
177  kBeamH2PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
178  "3mmHorn2X","+/- 3mm Horn 2 X Position");
179  const BeamSyst
180  kBeamH2PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
181  "3mmHorn2Y","+/- 3mm Horn 2 Y Position");
182 
183  ///Target z position shift +/-7mm
184  const BeamSyst
185  kBeamTarget((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
186  "7mmTargetZ","+/- 7mm Target z Position");
187 
188  ///Constant magnetic field in decay pipe
189  const BeamSyst
190  kBeamMagField((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
191  "Magnetic Field in Decay Pipe",
192  "MagnFieldDecayPipe","Magnetic Field in Decay Pipe");
193 
194  ///New vs old NuMI target
195  const BeamSyst
196  kBeamNewTarget((FindCAFAnaDir()+"/data/beam/ratioNewToOldTgt.root"),
197  "NewTarget", "New vs Old NuMI Target");
198 
199  ///New target at 160 kA to new target at 200 kA
200  const BeamSyst
201  kBeamNewTgt160to200kA((FindCAFAnaDir()+"/data/beam/ratioNewTgt160to200kA.root"),
202  "NewTgt160to200kA", "New target 160 kA/200 kA");
203 
204  ///New target at 180 kA to new target at 200 kA
205  const BeamSyst
206  kBeamNewTgt180to200kA((FindCAFAnaDir()+"/data/beam/ratioNewTgt180to200kA.root"),
207  "NewTgt180to200kA", "New target 180 kA/200 kA");
208 
209  ///Water layer on horn 1: +/- 1mm
210  const BeamSyst
211  kBeamGeomWater((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
212  "1mmHornWater","+/- 1mm water on Horn 1");
213 
214  ///All Beam Transport systematics combined in quadratures
215  const BeamSyst
216  kBeamAllTransport((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),
217  "totErr",
218  "beamTransportComb","Combined Beam Transport Systematics");
219 
220  //----------------------------------------------------------------------
222  const std::string& histName)
223  : fFileName(fname), fHistName(histName), fHistos()
224  {
225  }
226 
227  //----------------------------------------------------------------------
229  {
230  for(int l = 0; l < kNumFluxType; ++l){
231  for(int i = 0; i < kNumDets; ++i){
232  for(int j = 0; j < kNumSigns; ++j){
233  for(int k = 0; k < kNumFlavors; ++k){
234  delete fHistos[l][i][j][k];
235  fHistos[l][i][j][k] = 0;
236  }
237  }
238  }
239  }
240  }
241 
242  //----------------------------------------------------------------------
244  {
245  // Already been called
246  if(fHistos[0][0][0][0]) return;
247 
248  //open file
249  TFile fin (fFileName.c_str(),"read");
250  if(fin.IsZombie()){
251  std::cerr << "Warning: couldn't open " << fFileName << std::endl;
252  return;
253  }
254 
255  //Check that desired systematic is in file
256  // Loop through the detector folder, over all histograms
257 
258  TIter iterHist(gDirectory->GetListOfKeys());
259  TKey* keyHist;
260  int foundHisto=0;
261  bool isSeparatedByFlavor = false;
262 
263  while((keyHist = (TKey*)iterHist())) {
264  TString histName = keyHist->GetName(); // Get a histogram name
265  if(histName.Contains(fHistName)){ //is it the desired syst
266  foundHisto++;
267 
268  int flux=0;
269  if(histName.Contains("FHC")) flux = kFHC;
270  if(histName.Contains("RHC")) flux = kRHC;
271 
272  int det = 0;
273  if(histName.Contains("ND")) det = kND;
274  if(histName.Contains("FD")) det = kFD;
275 
276  int sign = 0;
277  if(histName.Contains("minus")) sign = kMinus;
278  if(histName.Contains("plus" )) sign = kPlus;
279 
280  int flav = 0;
281  if(histName.Contains("numu")) flav = kNumu;
282  if(histName.Contains("nue")) flav = kNue;
283  if(histName.Contains("anumu")) flav = kAntiNumu;
284  if(histName.Contains("anue")) flav = kAntiNue;
285 
286  //store relevant histograms
287  fHistos[flux][det][sign][flav] = (TH1D*) fin.Get(histName)->Clone();
288  // disassociate it from the file it came from
289  fHistos[flux][det][sign][flav]->SetDirectory(0);
290  }
291  if(histName.Contains("numu")){
292  isSeparatedByFlavor=true;
293  }
294  }
295 
296  if (foundHisto==0) {
297  std::cerr << "Beam systematic histogram " << fHistName
298  << " is not in file " << fFileName <<"; aborting" << std::endl;
299  abort();
300  }
301 
302  if (!isSeparatedByFlavor) {
303  std::cerr << "Beam systematic "<< fHistName
304  << " doesn't have required flavor information; aborting"
305  << std::endl;
306  abort();
307  }
308  fin.Close();
309  }
310 
311  //----------------------------------------------------------------------
313  ESign sign) const
314  {
316 
317  int flux = kFHC;
318  if(nu->isRHC) flux = kRHC;
319 
320  // Need different histogram per detector
321  int det = 0;
322  switch (nu->det)
323  {
324  case caf::kNEARDET : det = kND; break;
325  case caf::kFARDET : det = kFD; break;
326  default: std::cerr <<"Wrong detector; ignore" << std::endl; return 1;
327  }
328 
329  //Set the different weights for p/m sigma and neutrino flavor/energy
330 
331  double energy = nu->E; // True neutrino energy
332 
333  //Find the neutrino flavor
334  int flav = 0;
335  switch(nu->pdgorig)
336  {
337  case 14: flav = kNumu; break;
338  case 12: flav = kNue; break;
339  case -14: flav = kAntiNumu; break;
340  case -12: flav = kAntiNue; break;
341  default: std::cerr << "Wrong nu flavor; ignore" << std::endl; return 1;
342  }
343 
344  //Find the right histogram for the detector and sign
345  TH1D* h = fHistos[flux][det][sign][flav];
346  if (h == 0){
347  std::cerr << fHistName+": Can't find desired histogram; ignore"<< std:: endl;
348  return 1;
349  }
350  //Only use weights if the energy is in the range of the histogram
351  if (energy > h->GetXaxis()->GetXmin() &&
352  energy < h->GetXaxis()->GetXmax() ){
353  return h->Interpolate(energy);
354  }
355 
356  //Main use case is for PPFX smoothened weights.
357  //Return the weight in the last bin if the energy is larger than the range of the histogram
358  //Generalized for all BeamWeightFuncs and BeamSysts
359  if (energy > h->GetXaxis()->GetXmax())
360  return h->GetBinContent(h->GetNbinsX());
361 
362  return 1;
363  }
364 
365  //----------------------------------------------------------------------
367  const std::string& shortname,
368  const std::string& latexname)
369  : ISyst(shortname, latexname),
370  BeamSystOrWeightBase(fname, shortname)
371  {
372  }
373  //----------------------------------------------------------------------
375  const std::string& histname,
376  const std::string& shortname,
377  const std::string& latexname)
378  : ISyst(shortname, latexname),
379  BeamSystOrWeightBase(fname, histname)
380  {
381  }
382 
383  //----------------------------------------------------------------------
385  caf::SRNeutrinoProxy* sr, double& weight) const
386  {
387  if(sigma == 0) return;
388 
389  // Need different histogram for plus or minus sigma
390  ESign sign = kPlus;
391  if(sigma<0) sign = kMinus;
392 
393  weight *=1+(GetWeight(sr, sign)-1.)*TMath::Abs(sigma);
394  weight = std::max(0., weight);
395  }
396 
397 }
BeamSyst * GetPPFXPrincipals(int PCIdx)
Definition: BeamSysts.cxx:63
Near Detector underground.
Definition: SREnums.h:10
T max(const caf::Proxy< T > &a, T b)
TString fin
Definition: Style.C:24
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const BeamSyst kBeamNewTarget((FindCAFAnaDir()+"/data/beam/ratioNewToOldTgt.root"),"NewTarget","New vs Old NuMI Target")
New vs old NuMI target.
Definition: BeamSysts.h:119
const BeamSyst kBeamTarget((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"7mmTargetZ","+/- 7mm Target z Position")
Target z position shift +/-7mm.
Definition: BeamSysts.h:113
Beam systematics:
Definition: BeamSysts.h:37
const Var weight
const BeamSyst kBeamHornCurrent((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"2kA","2kAHornCurrent","+/- 2kA Horn Current")
Horn Current +/-2kA.
Definition: BeamSysts.h:97
NuTruthVar & GetPPFXFluxUnivWgtSmooth_NT(int UnivIdx)
Definition: BeamSysts.cxx:23
OStream cerr
Definition: OStream.cxx:7
BeamSyst * GetPPFXPrincipalsSBAna(int PCIdx)
Definition: BeamSysts.cxx:114
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Loaders::FluxType flux
BeamSyst * GetFluxPrincipals2018(int PCIdx)
Definition: BeamSysts.cxx:80
const BeamSyst kBeamGeomWater((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmHornWater","+/- 1mm water on Horn 1")
Water layer on horn 1: +/- 1mm.
Definition: BeamSysts.h:128
const BeamSyst kBeamNewTgt160to200kA((FindCAFAnaDir()+"/data/beam/ratioNewTgt160to200kA.root"),"NewTgt160to200kA","New target 160 kA/200 kA")
New target at 160 kA to new target at 200 kA.
Definition: BeamSysts.h:122
const BeamSyst kBeamH1PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1X","+/- 3mm Horn 1 X Position")
Horn 1 and 2 position +-3mm in X and Y separately.
Definition: BeamSysts.h:107
BeamSystOrWeightBase(const std::string &fname, const std::string &shortName)
Definition: BeamSysts.cxx:221
BeamSyst * GetPPFXUnivSyst(int UnivIdx)
Definition: BeamSysts.cxx:39
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
double energy
Definition: plottest35.C:25
const BeamSyst kBeamAllTransport((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"totErr","beamTransportComb","Combined Beam Transport Systematics")
All Beam Transport systematics combined in quadratures.
Definition: BeamSysts.h:131
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
double operator()(const caf::SRNeutrinoProxy *sr) const
Definition: BeamSysts.cxx:14
BeamSyst * GetFluxPrincipals2020(int PCIdx)
Definition: BeamSysts.cxx:97
double sigma(TH1F *hist, double percentile)
const BeamSyst kBeamPosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftY","Beam Position Y")
Definition: BeamSysts.h:104
const BeamSyst kBeamNewTgt180to200kA((FindCAFAnaDir()+"/data/beam/ratioNewTgt180to200kA.root"),"NewTgt180to200kA","New target 180 kA/200 kA")
New target at 180 kA to new target at 200 kA.
Definition: BeamSysts.h:125
const BeamSyst kPPFXCVSyst((FindCAFAnaDir()+"/data/beam/ppfx_smooth_multiverse_weights_2020.root"),"ppfx_cv","PPFX CV")
Definition: BeamSysts.h:69
BeamSyst * GetFluxPrincipalsND2018(int PCIdx)
Definition: BeamSysts.cxx:132
void TruthShift(double sigma, caf::SRNeutrinoProxy *sr, double &weight) const override
Definition: BeamSysts.cxx:384
const BeamSyst kBeamH1PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1Y","+/- 3mm Horn 1 Y Position")
Definition: BeamSysts.h:108
BeamWeightFunc(const std::string &fname, const std::string &histname)
Definition: BeamSysts.h:57
double GetWeight(const caf::SRNeutrinoProxy *nu, ESign sign) const
Definition: BeamSysts.cxx:312
const BeamSyst kBeamSpotSize((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"0.2mmBeamSpotSize","0p2mmBeamSpotSize"," 1.3 +/- 0.2 mm Spot Size")
Beam Spot Size 1.3 +/- 0.2 mm both XY.
Definition: BeamSysts.h:100
const BeamSyst kBeamH2PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2X","+/- 3mm Horn 2 X Position")
Definition: BeamSysts.h:109
const BeamSyst kBeamPosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftX","Beam Position X")
Beam position on target +-1 mm, X/Y separately.
Definition: BeamSysts.h:103
void InitializeHistograms() const
Definition: BeamSysts.cxx:243
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const BeamSyst kBeamH2PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2Y","+/- 3mm Horn 2 Y Position")
Definition: BeamSysts.h:110
const BeamSyst kBeamMagField((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"Magnetic Field in Decay Pipe","MagnFieldDecayPipe","Magnetic Field in Decay Pipe")
Constant magnetic field in decay pipe.
Definition: BeamSysts.h:116
Template for Var and SpillVar.
BeamSyst(const std::string &fname, const std::string &shortname, const std::string &latexname)
Definition: BeamSysts.cxx:366
TH1D * fHistos[kNumFluxType][kNumDets][kNumSigns][kNumFlavors]
Definition: BeamSysts.h:33
def sign(x)
Definition: canMan.py:197
enum BeamMode string