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