NueCCIncMRECorrection.cxx
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
6 #include "TObjString.h"
7 #include "TVectorD.h"
8 #include "TVector3.h"
9 #include "TH1F.h"
10 #include "TH3F.h"
11 #include "TH3.h"
12 #include "TCanvas.h"
13 #include <cassert>
14 #include <vector>
15 #include <string>
16 #include <iostream>
17 #include "TSVDUnfold.h"
18 #include "TDecompSVD.h"
19 
20 #include "/cvmfs/nova.opensciencegrid.org/externals/roounfold/v1_1_1b/NULL-e17--prof/NULL/src/RooUnfoldResponse.h"
21 #include "/cvmfs/nova.opensciencegrid.org/externals/roounfold/v1_1_1b/NULL-e17--prof/NULL/src/RooUnfoldBayes.h"
22 #include "/cvmfs/nova.opensciencegrid.org/externals/roounfold/v1_1_1b/NULL-e17--prof/NULL/src/RooUnfoldErrors.h"
23 #include "/cvmfs/nova.opensciencegrid.org/externals/roounfold/v1_1_1b/NULL-e17--prof/NULL/src/RooUnfoldSvd.h"
24 
25 
26 namespace ana{
27  namespace nueccinc{
30  HistAxis recoAxis1D,HistAxis recoAxis2D,
31  HistAxis trueAxis1D,HistAxis trueAxis2D,
32  Cut preselection, Cut full_selection,
33  SystShifts shift, SystShifts shiftMC, Var wei, Var weiMC)
34  :fDataPre1D(lDA, recoAxis1D, preselection, shift, wei),
35  fDataPre2D(lDA, recoAxis2D, preselection, shift, wei),
36  fDataFull1D(lDA, recoAxis1D, full_selection, shift, wei),
37  fDataFull2D(lDA, recoAxis2D, full_selection, shift, wei),
38  fMCPre1D(lMC, recoAxis1D, preselection, shiftMC, weiMC),
39  fMCPre2D(lMC, recoAxis2D, preselection, shiftMC, weiMC),
40  fMCFull1D(lMC, recoAxis1D, full_selection, shiftMC, weiMC),
41  fMCFull2D(lMC, recoAxis2D, full_selection, shiftMC, weiMC),
42  fUnfoldingPre1D(lMC,recoAxis1D,trueAxis1D,
43  preselection,shiftMC,weiMC),
44  fUnfoldingPre2D(lMC,recoAxis2D,trueAxis2D,
45  preselection,shiftMC,weiMC),
46  fUnfoldingFull1D(lMC,recoAxis1D,trueAxis1D,
47  full_selection,shiftMC,weiMC),
48  fUnfoldingFull2D(lMC,recoAxis2D,trueAxis2D,
49  full_selection,shiftMC,weiMC)
50  {}
51 
53  {
54  if(dimension == "1d"){
55  TH1D* hDataPre = (TH1D*)fDataPre1D.ToTH1(6.72829e20);
56  TH1D* hDataFull = (TH1D*)fDataFull1D.ToTH1(6.72829e20);
57  TH1D* hMCPre = (TH1D*)fMCPre1D.ToTH1(6.72829e20);
58  TH1D* hMCFull = (TH1D*)fMCFull1D.ToTH1(6.72829e20);
59  TH2D* hResponsePre = (TH2D*)fUnfoldingPre1D.ToTH2(6.72829e20);
60  TH2D* hResponseFull = (TH2D*)fUnfoldingFull1D.ToTH2(6.72829e20);
61 
62  TH1D* hRecoMCPre =
63  (TH1D*)hResponsePre->ProjectionX(ana::UniqueName().c_str());
64  TH1D* hRecoMCFull =
65  (TH1D*)hResponseFull->ProjectionX(ana::UniqueName().c_str());
66  TH1D* hTrueMCPre =
67  (TH1D*)hResponsePre->ProjectionY(ana::UniqueName().c_str());
68  TH1D* hTrueMCFull =
69  (TH1D*)hResponseFull->ProjectionY(ana::UniqueName().c_str());
70 
71 
72 
73  RooUnfoldResponse unf_pre(hRecoMCPre,hTrueMCPre,
74  hResponsePre, "Pre_Response");
75 
76  RooUnfoldResponse unf_full(hRecoMCFull,hTrueMCFull,
77  hResponseFull,"Full_Response");
78 
79  RooUnfoldSvd unf_da_pre(&unf_pre, hDataPre, 3,
80  ana::UniqueName().c_str());
81 
82  RooUnfoldSvd unf_mc_pre(&unf_pre, hMCPre, 3,
83  ana::UniqueName().c_str());
84 
85  RooUnfoldSvd unf_da_full(&unf_full, hDataFull, 3,
86  ana::UniqueName().c_str());
87 
88  RooUnfoldSvd unf_mc_full(&unf_full, hMCFull, 3,
89  ana::UniqueName().c_str());
90 
91  TH1F* hData_Correction_denom =
92  (TH1F*)unf_da_pre.Hreco(RooUnfold::ErrorTreatment::kCovariance);
93  hData_Correction_denom->SetName(ana::UniqueName().c_str());
94 
95  TH1F* hMC_Correction_denom =
96  (TH1F*)unf_mc_pre.Hreco(RooUnfold::ErrorTreatment::kCovariance);
97  hMC_Correction_denom->SetName(ana::UniqueName().c_str());
98 
99  TH1F* hData_Correction_numerator =
100  (TH1F*)unf_da_full.Hreco(RooUnfold::ErrorTreatment::kCovariance);
101  hData_Correction_numerator->SetName(ana::UniqueName().c_str());
102 
103  TH1F* hMC_Correction_numerator =
104  (TH1F*)unf_mc_full.Hreco(RooUnfold::ErrorTreatment::kCovariance);
105  hMC_Correction_numerator->SetName(ana::UniqueName().c_str());
106 
107  hData_Correction_numerator->Divide(hData_Correction_denom);
108  hMC_Correction_numerator->Divide(hMC_Correction_denom);
109 
110  hData_Correction_numerator->Divide(hMC_Correction_numerator);
111 
112  return hData_Correction_numerator;
113 
114  }
115  else{
116  TH1D* hDataPre = (TH1D*)fDataPre2D.ToTH1(6.72829e20);
117  TH1D* hDataFull = (TH1D*)fDataFull2D.ToTH1(6.72829e20);
118  TH1D* hMCPre = (TH1D*)fMCPre2D.ToTH1(6.72829e20);
119  TH1D* hMCFull = (TH1D*)fMCFull2D.ToTH1(6.72829e20);
120  TH2D* hResponsePre = (TH2D*)fUnfoldingPre2D.ToTH2(6.72829e20);
121  TH2D* hResponseFull = (TH2D*)fUnfoldingFull2D.ToTH2(6.72829e20);
122 
123  TH1D* hRecoMCPre =
124  (TH1D*)hResponsePre->ProjectionX(ana::UniqueName().c_str());
125  TH1D* hRecoMCFull =
126  (TH1D*)hResponseFull->ProjectionX(ana::UniqueName().c_str());
127  TH1D* hTrueMCPre =
128  (TH1D*)hResponsePre->ProjectionY(ana::UniqueName().c_str());
129  TH1D* hTrueMCFull =
130  (TH1D*)hResponseFull->ProjectionY(ana::UniqueName().c_str());
131 
132 
133 
134  RooUnfoldResponse unf_pre(hRecoMCPre,hTrueMCPre,
135  hResponsePre, "Pre_Response");
136 
137  RooUnfoldResponse unf_full(hRecoMCFull,hTrueMCFull,
138  hResponseFull,"Full_Response");
139 
140  RooUnfoldBayes unf_da_pre(&unf_pre, hDataPre, 2,
141  ana::UniqueName().c_str());
142 
143  RooUnfoldBayes unf_mc_pre(&unf_pre, hMCPre, 2,
144  ana::UniqueName().c_str());
145 
146  RooUnfoldBayes unf_da_full(&unf_full, hDataFull, 2,
147  ana::UniqueName().c_str());
148 
149  RooUnfoldBayes unf_mc_full(&unf_full, hMCFull, 2,
150  ana::UniqueName().c_str());
151 
152  TH1F* hData_Correction_denom =
153  (TH1F*)unf_da_pre.Hreco(RooUnfold::ErrorTreatment::kCovariance);
154  hData_Correction_denom->SetName(ana::UniqueName().c_str());
155 
156  TH1F* hMC_Correction_denom =
157  (TH1F*)unf_mc_pre.Hreco(RooUnfold::ErrorTreatment::kCovariance);
158  hMC_Correction_denom->SetName(ana::UniqueName().c_str());
159 
160  TH1F* hData_Correction_numerator =
161  (TH1F*)unf_da_full.Hreco(RooUnfold::ErrorTreatment::kCovariance);
162  hData_Correction_numerator->SetName(ana::UniqueName().c_str());
163 
164  TH1F* hMC_Correction_numerator =
165  (TH1F*)unf_mc_full.Hreco(RooUnfold::ErrorTreatment::kCovariance);
166  hMC_Correction_numerator->SetName(ana::UniqueName().c_str());
167 
168  hData_Correction_numerator->Divide(hData_Correction_denom);
169  hMC_Correction_numerator->Divide(hMC_Correction_denom);
170 
171  hData_Correction_numerator->Divide(hMC_Correction_numerator);
172 
173  return hData_Correction_numerator;
174  }
175  }
176 
177 
179  {
180  std::vector<Spectrum> hResult;
181  if(dimName == "1d"){
182  hResult.push_back(fDataPre1D);
183  }
184  if(dimName == "2d"){
185  hResult.push_back(fDataPre2D);
186  }
187  return hResult[0];
188  }
189 
190 
192  {
193  TH1F* hResult = NueCCIncMRECorrection::CalcEfficiencyCorrection(dimension);
194  return hResult;
195  }
196 
198  {
199  Spectrum preselected = NueCCIncMRECorrection::Preselected(dimName);
200 
201  TH1F* hResult = (TH1F*)preselected.ToTH1(6.72829e20);
202  return hResult;
203  }
204 
206  {
207  std::vector<Spectrum> hResult;
208  if(dimName == "1d"){
209  hResult.push_back(fDataFull1D);
210  }
211  if(dimName == "2d"){
212  hResult.push_back(fDataFull2D);
213  }
214  return hResult[0];
215  }
216 
217 
219  {
220  Spectrum selected = NueCCIncMRECorrection::Selected(dimName);
221 
222  TH1F* hResult = (TH1F*)selected.ToTH1(6.72829e20);
223  return hResult;
224  }
225 
226 
228  {
229  Spectrum preselected = NueCCIncMRECorrection::Preselected(dimName);
230  Spectrum selected = NueCCIncMRECorrection::Selected(dimName);
231 
232  TH1F* hResult = (TH1F*)selected.ToTH1(6.72829e20);
233  hResult->Divide((TH1F*)preselected.ToTH1(6.72829e20));
234  return hResult;
235  }
236 
237 
238 
239  void NueCCIncMRECorrection::SaveTo(TDirectory* dir) const{
240  TDirectory *tmp = gDirectory;
241  dir->cd();
242 
243  TObjString("NueCCIncMRECorrection").Write("type");
244 
245  fDataPre1D.SaveTo(dir->mkdir("fDataPre1D"));
246  fDataPre2D.SaveTo(dir->mkdir("fDataPre2D"));
247  fDataFull1D.SaveTo(dir->mkdir("fDataFull1D"));
248  fDataFull2D.SaveTo(dir->mkdir("fDataFull2D"));
249  fMCPre1D.SaveTo(dir->mkdir("fMCPre1D"));
250  fMCPre2D.SaveTo(dir->mkdir("fMCPre2D"));
251  fMCFull1D.SaveTo(dir->mkdir("fMCFull1D"));
252  fMCFull2D.SaveTo(dir->mkdir("fMCFull2D"));
253  fUnfoldingPre1D.SaveTo(dir->mkdir("fUnfoldingPre1D"));
254  fUnfoldingPre2D.SaveTo(dir->mkdir("fUnfoldingPre2D"));
255  fUnfoldingFull1D.SaveTo(dir->mkdir("fUnfoldingFull1D"));
256  fUnfoldingFull2D.SaveTo(dir->mkdir("fUnfoldingFull2D"));
257  tmp->cd();
258  }
259 
260  std::unique_ptr<NueCCIncMRECorrection>
262 
263  Spectrum* da_pre_1D =
264  ana::LoadFrom<Spectrum>(dir->GetDirectory("fDataPre1D")).release();
265  Spectrum* da_pre_2D =
266  ana::LoadFrom<Spectrum>(dir->GetDirectory("fDataPre2D")).release();
267  Spectrum* da_full_1D =
268  ana::LoadFrom<Spectrum>(dir->GetDirectory("fDataFull1D")).release();
269  Spectrum* da_full_2D =
270  ana::LoadFrom<Spectrum>(dir->GetDirectory("fDataFull2D")).release();
271 
272 
273  Spectrum* mc_pre_1D =
274  ana::LoadFrom<Spectrum>(dir->GetDirectory("fMCPre1D")).release();
275  Spectrum* mc_pre_2D =
276  ana::LoadFrom<Spectrum>(dir->GetDirectory("fMCPre2D")).release();
277  Spectrum* mc_full_1D =
278  ana::LoadFrom<Spectrum>(dir->GetDirectory("fMCFull1D")).release();
279  Spectrum* mc_full_2D =
280  ana::LoadFrom<Spectrum>(dir->GetDirectory("fMCFull2D")).release();
281 
282  Spectrum* unfolding_pre_1D =
283  ana::LoadFrom<Spectrum>(dir->GetDirectory("fUnfoldingPre1D")).release();
284  Spectrum* unfolding_pre_2D =
285  ana::LoadFrom<Spectrum>(dir->GetDirectory("fUnfoldingPre2D")).release();
286  Spectrum* unfolding_full_1D =
287  ana::LoadFrom<Spectrum>(dir->GetDirectory("fUnfoldingFull1D")).release();
288  Spectrum* unfolding_full_2D =
289  ana::LoadFrom<Spectrum>(dir->GetDirectory("fUnfoldingFull2D")).release();
290 
291 
292 
293 
294  return std::unique_ptr<NueCCIncMRECorrection>
295  (
296  new NueCCIncMRECorrection(*da_pre_1D,*da_pre_2D,*da_full_1D,*da_full_2D,
297  *mc_pre_1D,*mc_pre_2D,*mc_full_1D,*mc_full_2D,
298  *unfolding_pre_1D, *unfolding_pre_2D,
299  *unfolding_full_1D, *unfolding_full_2D)
300  );
301 
302  }
303 
304  }//namespace nueccinc
305 }//namespace ana
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1F * getPreselected1D(std::string dimName)
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:148
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
TH1F * CalcEfficiencyCorrection(std::string dimension)
Float_t tmp
Definition: plot.C:36
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
NueCCIncMRECorrection(SpectrumLoaderBase &lDA, SpectrumLoaderBase &lMC, HistAxis recoAxis1D, HistAxis recoAxis2D, HistAxis trueAxis1D, HistAxis trueAxis2D, Cut preselection, Cut full_selection, SystShifts shift, SystShifts shiftMC, Var wei, Var weiMC)
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
Spectrum Preselected(std::string dimName)
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
TH1F * getEfficiencyCorrection(std::string dimension)
TH1F * getEfficiency1D(std::string dimName)
static std::unique_ptr< NueCCIncMRECorrection > LoadFrom(TDirectory *dir)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Spectrum Selected(std::string dimName)