DecompWeights.cxx
Go to the documentation of this file.
2 
4 
5 #include <cassert>
6 #include <iostream>
7 
8 #include "TDirectory.h"
9 #include "TObjString.h"
11 
12 #include "TFile.h"
13 #include "TH1.h"
14 #include "TKey.h"
15 
16 namespace ana
17 {
19  const std::string& shortname,
20  const Var& decompVar)
21  : fHistoFname(fname), fShortName(shortname), fDecompVar(decompVar), fHistos()
22  {
23  }
24 
25  //----------------------------------------------------------------------
27  {
28  for(int l = 0; l < kNumComponents; ++l){
29  for(int k = 0; k < kIsOscillated; ++k){
30  delete fHistos[l][k];
31  fHistos[l][k] = 0;
32  }
33  }
34  }
35 
36  //----------------------------------------------------------------------
38  {
39  // Already been called
40  if(fHistos[0][0]) return;
41 
42  //open file
43  TFile fin (fHistoFname.c_str(),"read");
44  if(fin.IsZombie()){
45  std::cerr << "Warning: couldn't open " << fHistoFname << std::endl;
46  return;
47  }
48 
49  //Check that desired systematic is in file
50  // Loop through the detector folder, over all histograms
51 
52  TIter iterHist(gDirectory->GetListOfKeys());
53  TKey* keyHist;
54  int foundHisto=0;
55  bool isSeparatedByFlavor = false;
56 
57  while((keyHist = (TKey*)iterHist())) {
58  TString histName = keyHist->GetName(); // Get a histogram name
59  if(histName.Contains(fShortName)){ //is it the desired syst
60  foundHisto++;
61 
62  int flav = 0;
63  if(histName.Contains("numu")) flav = kNumu;
64  if(histName.Contains("nue")) flav = kNue;
65  if(histName.Contains("antinumu")) flav = kAntiNumu;
66  if(histName.Contains("antinue")) flav = kAntiNue;
67  if(histName.Contains("nc")) flav = kNC;
68 
69  int osc = 0;
70  if(histName.Contains("beam")) osc = kBeam;
71  if(histName.Contains("app")) osc = kAppearance;
72 
73  //store relevant histograms
74  fHistos[flav][osc] = (TH1D*) fin.Get(histName)->Clone();
75  // disassociate it from the file it came from
76  fHistos[flav][osc]->SetDirectory(0);
77  }
78  if(histName.Contains("numu")){
79  isSeparatedByFlavor=true;
80  }
81  }
82 
83  if (foundHisto==0) {
84  std::cerr << "Decomp Ratios for "<< fShortName
85  << " not in this file; aborting" << std::endl;
86  abort();
87  }
88 
89  if (!isSeparatedByFlavor) {
90  std::cerr << "Decomp Ratios for "<< fShortName
91  << " doesn't have required flavor information; aborting"
92  << std::endl;
93  abort();
94  }
95  fin.Close();
96  }
97 
98  //----------------------------------------------------------------------
100  {
102 
103  // Leave weight unaltered
104  if(sr->mc.nnu == 0) return 1;
105 
106  //Set the different weights for p/m sigma and neutrino flavor/energy
107 
108  double decompbin = fDecompVar(sr); // True neutrino energy
109 
110  //Find the component
111  int flav = 0;
112  switch(sr->mc.nu[0].pdg){
113  case 14: flav = kNumu; break;
114  case 12: flav = kNue; break;
115  case -14: flav = kAntiNumu; break;
116  case -12: flav = kAntiNue; break;
117  default: std::cerr << "Wrong nu flavor; ignore" << std::endl; return 1;
118  }
119  if(!sr->mc.nu[0].iscc)
120  flav = kNC;
121 
122  int osc = 0;
123  if(sr->mc.nu[0].pdgorig != sr->mc.nu[0].pdg)
124  osc = kAppearance;
125  else
126  osc = kBeam;
127 
128  //Find the right histogram
129  TH1D* h = fHistos[flav][osc];
130  if (h == 0){
131  std::cerr << fShortName+": Can't find desired histogram; ignore"<< std:: endl;
132  return 1;
133  }
134  //Only use weights if the energy is in the range of the histogram
135  if (decompbin >= h->GetXaxis()->GetXmin() &&
136  decompbin < h->GetXaxis()->GetXmax() ){
137  return h->GetBinContent(h->FindBin(decompbin));
138  }
139 
140  return 1;
141  }
142 }
TString fin
Definition: Style.C:24
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * fHistos[kNumComponents][kIsOscillated]
Definition: DecompWeights.h:34
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
std::string fHistoFname
Definition: DecompWeights.h:30
OStream cerr
Definition: OStream.cxx:7
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
DecompWeightFunc(const std::string &fname, const std::string &shortname, const Var &decompVar)
double GetWeight(const caf::SRProxy *sr) const
caf::StandardRecord * sr
Oscillation probability calculators.
Definition: Calcs.h:5
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
std::string fShortName
Definition: DecompWeights.h:31
void InitializeHistograms() const
enum BeamMode string