Functions | Variables
Calculate_covariances.C File Reference
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include <TRandom.h>
#include <TMatrixD.h>
#include <TRandom3.h>
#include <vector>
#include <iostream>
#include <string>

Go to the source code of this file.

Functions

TH1D * GetWanted1DFromENu (TH1D *h1In, double &IntegralVal, std::string name)
 
TH1D * GetWanted1DFromQ2 (TH1D *h1In, double &IntegralVal, std::string name)
 
TH1D * GetWanted1DFromMuKin (TH2D *h2In, double &IntegralVal, std::string name)
 
TMatrixD calc_cov_matrix (std::vector< TH1D * > &vhist, TH1D *hIncv)
 
TMatrixD calc_corr_matrix (TMatrixD &mcov)
 
TMatrixD calc_fer_matrix (TMatrixD &mcov, TH1D *hIncv)
 
TH1D * GetSymSyst (TH1D *hInshiftp, TH1D *hInshiftn, TH1D *hIncv, std::string SystName, std::string VarName)
 
TH2D * GetSymSyst (TH2D *hInshiftp, TH2D *hInshiftn, TH2D *hIncv, std::string SystName, std::string VarName)
 
void AddShift (TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
 
void AddOneSideShift (TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
 
void scaling (TH1D *hIn, const double shape_scale)
 
void MakeCov (std::string histName, TFile *fIn, TH1D *hIncv, std::vector< TH1D * > &hIngenie, std::vector< TH1D * > &hInppfx, std::vector< TH1D * > &hInsym, std::vector< TH1D * > &hInasym, int thismode)
 
void MakeCovUniverses (std::string VarName, TFile *fIn, TH1D *hIncv, std::vector< TH1D * > &hInuniv, int thismode)
 
void MakeCovSymSyst (std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInsym, int thismode)
 
void MakeCovASymSyst (std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInasym, int thismode)
 
void Calculate_covariances ()
 
std::vector< intstarting (const char *fnameIn, const char *fnameonly1DIn, const char *fnameOut, const char *cmode, const char *ccomp)
 
void Calculate_covariances (const char *fnameIn, const char *fnameonly1DIn, const char *fnameOut, const char *cmode, const char *ccomp)
 
int main (int argc, const char *argv[])
 

Variables

const int Nppfx = 100
 
const int Ngenie = 1000
 
const int Nuniv = Ngenie*Nppfx
 
const int Nbfoc = 22
 
const std::string bfoc [Nbfoc]
 
const int Ndsyst = 12
 
const std::string dsyst [Ndsyst]
 
const int Ncv = 3
 
const std::string cv [Ncv] = {"cv","nom_good_seventh","neutron_nom"}
 
const int Nwanted_MuKinBins = 172
 
const int NMuCosBins = 13
 
const int MuTperCosBins [NMuCosBins] = {6,7,8,9,9,9,13,14,17,20,20,20,20}
 
const int Nwanted_ENuBins = 18
 
const int Nwanted_Q2Bins = 15
 
const int Nsymsyst = 16
 
const std::string symsyst [Nsymsyst]
 
const int basesymsyst_rdm [Nsymsyst]
 
const int Nasymsyst = 2
 
const std::string asymsyst [Nasymsyst] = {"CalibShape","Cherenkov"}
 
const int baseasymsyst_rdm [Nasymsyst] = {26000000,27000000}
 
TRandom3 * r3 = new TRandom3(0)
 

Function Documentation

void AddOneSideShift ( TH1D *  hUniv,
double &  shift,
TH1D *  hIncv,
TH1D *  hSystIn 
)

Definition at line 419 of file Calculate_covariances.C.

References plot_xsec_1d::Nbins.

Referenced by MakeCov(), and MakeCovASymSyst().

419  {
420  int Nbins = hUniv->GetXaxis()->GetNbins();
421  for(int ii=1;ii<=Nbins;ii++){
422  double univ_val = hUniv->GetBinContent(ii);
423  double cv_val = hIncv->GetBinContent(ii);
424  double syst_val = hSystIn->GetBinContent(ii);
425  double shift_val = syst_val - cv_val;
426  shift_val *= shift;
427  univ_val += shift_val;
428  hUniv->SetBinContent(ii,univ_val);
429  }
430 }
void AddShift ( TH1D *  hUniv,
double &  shift,
TH1D *  hIncv,
TH1D *  hSystIn 
)

Definition at line 407 of file Calculate_covariances.C.

References plot_xsec_1d::Nbins.

Referenced by ana::SystMaker::GetNominal(), MakeCov(), and MakeCovSymSyst().

407  {
408  int Nbins = hUniv->GetXaxis()->GetNbins();
409  for(int ii=1;ii<=Nbins;ii++){
410  double univ_val = hUniv->GetBinContent(ii);
411  double cv_val = hIncv->GetBinContent(ii);
412  double syst_val = hSystIn->GetBinContent(ii);
413  double shift_val = syst_val - cv_val;
414  shift_val *= shift;
415  univ_val += shift_val;
416  hUniv->SetBinContent(ii,univ_val);
417  }
418 }
TMatrixD calc_corr_matrix ( TMatrixD mcov)

Definition at line 352 of file Calculate_covariances.C.

References std::sqrt().

Referenced by MakeCov(), MakeCovASymSyst(), MakeCovSymSyst(), and MakeCovUniverses().

352  {
353  TMatrixD cormx(mcov.GetNrows(),mcov.GetNcols());
354 
355  for(int ii=0; ii<mcov.GetNrows(); ii++){
356  for(int jj=0; jj<mcov.GetNcols(); jj++){
357  cormx[ii][jj] = mcov[ii][jj]/(sqrt(mcov[ii][ii])*sqrt(mcov[jj][jj]));
358  }
359  }
360  return cormx;
361 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
TMatrixD calc_cov_matrix ( std::vector< TH1D * > &  vhist,
TH1D *  hIncv 
)

Definition at line 327 of file Calculate_covariances.C.

Referenced by MakeCov(), MakeCovASymSyst(), MakeCovSymSyst(), and MakeCovUniverses().

327  {
328 
329  unsigned int N = vhist.size();
330  unsigned int Nbinsx = vhist[0]->GetNbinsX();
331  TMatrixD covmx(Nbinsx,Nbinsx);
332  //universes
333  for(int iu=0; iu<N; iu++){
334  for(int ii=0; ii< Nbinsx; ii++){
335  double xi=vhist[iu]->GetBinContent(ii+1);
336  double ximean=hIncv->GetBinContent(ii+1);
337  for(int jj=0; jj< Nbinsx; jj++){
338  double xj=vhist[iu]->GetBinContent(jj+1);
339  double xjmean=hIncv->GetBinContent(jj+1);
340  covmx[ii][jj]+=(xi-ximean)*(xj-xjmean) ;
341  }
342  }
343  }//end of univ
344 
345  for(int ii=0; ii< Nbinsx; ii++){
346  for(int jj=0; jj< Nbinsx; jj++){
347  covmx[ii][jj]/=(N-1.0);
348  }
349  }
350  return covmx;
351 }
TMatrixD calc_fer_matrix ( TMatrixD mcov,
TH1D *  hIncv 
)

Definition at line 362 of file Calculate_covariances.C.

Referenced by MakeCovSymSyst().

362  {
363  TMatrixD fermx(mcov.GetNrows(),mcov.GetNcols());
364  for(int ii=0; ii<mcov.GetNrows(); ii++){
365  for(int jj=0; jj<mcov.GetNcols(); jj++){
366  if(hIncv->GetBinContent(ii+1)!=0 && hIncv->GetBinContent(jj+1)!=0){
367  fermx[ii][jj] = mcov[ii][jj] / ( (hIncv->GetBinContent(ii+1)) * (hIncv->GetBinContent(jj+1)) );
368  }
369  }
370  }
371  return fermx;
372 }
void Calculate_covariances ( )

Definition at line 616 of file Calculate_covariances.C.

References asymsyst, om::cout, allTimeWatchdog::endl, exit(), Nasymsyst, Nsymsyst, and symsyst.

Referenced by main(), and starting().

616  {
617  std::cout<<"Usage: Calculate_covariances [file name In (2D)] [file name In (1D)] [file name Out] [mode] [comp]"<<std::endl;
618  std::cout<<"mode = total, shape or notppfx"<<std::endl;
619  std::cout<<"comp = all"<<std::endl;
620  std::cout<<" GENIE"<<std::endl;
621  std::cout<<" PPFX"<<std::endl;
622  for(int ii=0;ii<Nsymsyst;ii++){
623  std::cout<<" "<<symsyst[ii]<<std::endl;
624  }
625  for(int ii=0;ii<Nasymsyst;ii++){
626  std::cout<<" "<<asymsyst[ii]<<std::endl;
627  }
628  exit (1);
629 }
const std::string asymsyst[Nasymsyst]
const int Nasymsyst
OStream cout
Definition: OStream.cxx:6
const std::string symsyst[Nsymsyst]
exit(0)
const int Nsymsyst
void Calculate_covariances ( const char *  fnameIn,
const char *  fnameonly1DIn,
const char *  fnameOut,
const char *  cmode,
const char *  ccomp 
)

Definition at line 91 of file Calculate_covariances.C.

References asymsyst, bfoc, om::cout, cv, dsyst, allTimeWatchdog::endl, makeBrightnessMap::fOut, GetSymSyst(), GetWanted1DFromENu(), GetWanted1DFromMuKin(), GetWanted1DFromQ2(), MECModelEnuComparisons::i, MakeCov(), MakeCovASymSyst(), MakeCovSymSyst(), MakeCovUniverses(), Nasymsyst, Nbfoc, Ncv, Ndsyst, Ngenie, Nppfx, Nsymsyst, scaling(), starting(), and symsyst.

91  {
92 
93  //mode:
94  std::vector<int> ioption = starting(fnameIn, fnameonly1DIn,fnameOut, cmode,ccomp);
95  int imode = ioption[0];
96  int icomp = ioption[1];
97  std::cout<<"=> check: (imode,icomp) = ("<<imode<<", "<<icomp<<")"<<std::endl;
98  bool is_shape = false;
99  if(imode==1)is_shape = true;
100 
101  TFile* fOut = new TFile(fnameOut,"recreate");
102 
103  //Declarations:
104  TFile* file1D = new TFile(fnameonly1DIn,"read");
105  TFile* file2D = new TFile(fnameIn,"read");
106  //ENu
107  std::vector<TH1D*> h1cv_ENu, h1cv_ENu_input;
108  std::vector<TH1D*> h1ppfx_ENu, h1genie_ENu, h1ppfx_ENu_input, h1genie_ENu_input;
109  std::vector<TH1D*> h1bfoc_ENu, h1dsyst_ENu, h1bfoc_ENu_input, h1dsyst_ENu_input;
110  //Q2
111  std::vector<TH1D*> h1cv_Q2, h1cv_Q2_input;
112  std::vector<TH1D*> h1ppfx_Q2, h1genie_Q2, h1ppfx_Q2_input, h1genie_Q2_input;
113  std::vector<TH1D*> h1bfoc_Q2, h1dsyst_Q2, h1bfoc_Q2_input, h1dsyst_Q2_input;
114  //MuKin
115  std::vector<TH2D*> h2cv_MuKin_input;
116  std::vector<TH1D*> h1cv_MuKin;
117  std::vector<TH2D*> h2ppfx_MuKin_input, h2genie_MuKin_input;
118  std::vector<TH1D*> h1ppfx_MuKin, h1genie_MuKin;
119  std::vector<TH2D*> h2bfoc_MuKin_input, h2dsyst_MuKin_input;
120  std::vector<TH1D*> h1bfoc_MuKin, h1dsyst_MuKin;
121 
122  //Getting
123  //CV
124  double integral_cv_ENu, integral_cv_Q2, integral_cv_MuKin;
125  double integral_syst_ENu, integral_syst_Q2, integral_syst_MuKin;
126  double integral_aux;
127  for(unsigned i=0;i<Ncv;i++){
128  h1cv_ENu_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_ENu" ,cv[i].c_str(),cv[i].c_str())) );
129  h1cv_Q2_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_Q2" ,cv[i].c_str(),cv[i].c_str())) );
130  h2cv_MuKin_input.push_back( (TH2D*)file2D->Get(Form("%s/xsec_%s_2DProj",cv[i].c_str(),cv[i].c_str())) );
131  if(i==0){
132  h1cv_ENu.push_back( GetWanted1DFromENu( h1cv_ENu_input[i], integral_cv_ENu, Form("henu_%s", cv[i].c_str())) );
133  h1cv_Q2.push_back( GetWanted1DFromQ2( h1cv_Q2_input[i] , integral_cv_Q2, Form("hq2_%s", cv[i].c_str())) );
134  h1cv_MuKin.push_back( GetWanted1DFromMuKin( h2cv_MuKin_input[i], integral_cv_MuKin,Form("hmukin_%s",cv[i].c_str())) );
135  }
136  else{ //we currently do not use the cv for "nom_good_seventh" and "neutron_nom"
137  h1cv_ENu.push_back( GetWanted1DFromENu( h1cv_ENu_input[i], integral_aux, Form("henu_%s", cv[i].c_str())) );
138  h1cv_Q2.push_back( GetWanted1DFromQ2( h1cv_Q2_input[i] , integral_aux, Form("hq2_%s", cv[i].c_str())) );
139  h1cv_MuKin.push_back( GetWanted1DFromMuKin( h2cv_MuKin_input[i], integral_aux, Form("hmukin_%s",cv[i].c_str())) );
140  }
141  }
142 
143  //PPFX
144  for(unsigned i=0;i<Nppfx;i++){
145  h1ppfx_ENu_input.push_back( (TH1D*)file1D->Get(Form("ppfx%d/xsec_ppfx%d_ENu" ,i,i)));
146  h1ppfx_Q2_input.push_back( (TH1D*)file1D->Get(Form("ppfx%d/xsec_ppfx%d_Q2" ,i,i)));
147  h2ppfx_MuKin_input.push_back( (TH2D*)file2D->Get(Form("ppfx%d/xsec_ppfx%d_2DProj" ,i,i)));
148  h1ppfx_ENu.push_back( GetWanted1DFromENu( h1ppfx_ENu_input[i], integral_syst_ENu, Form("henu_ppfx%d", i)) );
149  h1ppfx_Q2.push_back( GetWanted1DFromQ2( h1ppfx_Q2_input[i], integral_syst_Q2, Form("hq2_ppfx%d" , i)) );
150  h1ppfx_MuKin.push_back( GetWanted1DFromMuKin(h2ppfx_MuKin_input[i],integral_syst_MuKin, Form("hmukin_ppfx%d",i)) );
151  if(is_shape){
152  scaling(h1ppfx_ENu[i] ,integral_cv_ENu / integral_syst_ENu);
153  scaling(h1ppfx_Q2[i] ,integral_cv_Q2 / integral_syst_Q2);
154  scaling(h1ppfx_MuKin[i],integral_cv_MuKin/ integral_syst_MuKin);
155  }
156  }
157  h1ppfx_ENu_input.clear();
158  h1ppfx_Q2_input.clear();
159  h2ppfx_MuKin_input.clear();
160 
161  //GENIE
162  for(unsigned i=0;i<Ngenie;i++){
163  h1genie_ENu_input.push_back( (TH1D*)file1D->Get(Form("genie%d/xsec_genie%d_ENu" ,i,i)));
164  h1genie_Q2_input.push_back( (TH1D*)file1D->Get(Form("genie%d/xsec_genie%d_Q2" ,i,i)));
165  h2genie_MuKin_input.push_back((TH2D*)file2D->Get(Form("genie%d/xsec_genie%d_2DProj",i,i)));
166  h1genie_ENu.push_back( GetWanted1DFromENu( h1genie_ENu_input[i], integral_syst_ENu, Form("henu_genie%d", i)) );
167  h1genie_Q2.push_back( GetWanted1DFromQ2( h1genie_Q2_input[i], integral_syst_Q2, Form("hq2_genie%d", i)) );
168  h1genie_MuKin.push_back( GetWanted1DFromMuKin(h2genie_MuKin_input[i],integral_syst_MuKin, Form("hmukin_genie%d",i)) );
169  if(is_shape){
170  scaling(h1genie_ENu[i] ,integral_cv_ENu / integral_syst_ENu);
171  scaling(h1genie_Q2[i] ,integral_cv_Q2 / integral_syst_Q2);
172  scaling(h1genie_MuKin[i],integral_cv_MuKin/ integral_syst_MuKin);
173  }
174  }
175  h1genie_ENu_input.clear();
176  h1genie_Q2_input.clear();
177  h2genie_MuKin_input.clear();
178 
179  //BEAM FOC and DET SYST
180  for(unsigned i=0;i<Nbfoc;i++){
181  h1bfoc_ENu_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_ENu" ,bfoc[i].c_str(),bfoc[i].c_str())) );
182  h1bfoc_Q2_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_Q2" ,bfoc[i].c_str(),bfoc[i].c_str())) );
183  h2bfoc_MuKin_input.push_back( (TH2D*)file2D->Get(Form("%s/xsec_%s_2DProj",bfoc[i].c_str(),bfoc[i].c_str())) );
184  }
185  for(unsigned i=0;i<Ndsyst;i++){
186  h1dsyst_ENu_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_ENu" ,dsyst[i].c_str(),dsyst[i].c_str())) );
187  h1dsyst_Q2_input.push_back( (TH1D*)file1D->Get(Form("%s/xsec_%s_Q2" ,dsyst[i].c_str(),dsyst[i].c_str())) );
188  h2dsyst_MuKin_input.push_back( (TH2D*)file2D->Get(Form("%s/xsec_%s_2DProj",dsyst[i].c_str(),dsyst[i].c_str())) );
189  }
190 
191  //Getting systematics histograms:
192  //Symmetricals:
193  std::vector<TH1D*> h1symsyst_ENu_input, h1symsyst_Q2_input;
194  std::vector<TH2D*> h2symsyst_MuKin_input;
195  std::vector<TH1D*> h1symsyst_ENu, h1symsyst_Q2, h1symsyst_MuKin;
196  //Asymmetricals:
197  //one side syst cal shape
198  std::vector<TH1D*> h1asymsyst_ENu_input, h1asymsyst_Q2_input;
199  std::vector<TH2D*> h2asymsyst_MuKin_input;
200  std::vector<TH1D*> h1asymsyst_ENu, h1asymsyst_Q2, h1asymsyst_MuKin;
201  for(unsigned i=0;i<Nsymsyst;i++){
202  //Foc: 11 systematics
203  if(i<11){
204  int ifoc = 2*i;
205  h1symsyst_ENu_input.push_back( GetSymSyst(h1bfoc_ENu_input[ifoc] ,h1bfoc_ENu_input[ifoc+1] ,h1cv_ENu_input[0] ,symsyst[i], "enu") );
206  h1symsyst_Q2_input.push_back( GetSymSyst(h1bfoc_Q2_input[ifoc] ,h1bfoc_Q2_input[ifoc+1] ,h1cv_Q2_input[0] ,symsyst[i], "q2") );
207  h2symsyst_MuKin_input.push_back( GetSymSyst(h2bfoc_MuKin_input[ifoc],h2bfoc_MuKin_input[ifoc+1],h2cv_MuKin_input[0],symsyst[i], "mukin") );
208  }
209  //Det sym syst: 5 systematics
210  else{
211  int idsyst = 2*(i-11);
212  h1symsyst_ENu_input.push_back( GetSymSyst(h1dsyst_ENu_input[idsyst] ,h1dsyst_ENu_input[idsyst+1] ,h1cv_ENu_input[0] ,symsyst[i], "enu") );
213  h1symsyst_Q2_input.push_back( GetSymSyst(h1dsyst_Q2_input[idsyst] ,h1dsyst_Q2_input[idsyst+1] ,h1cv_Q2_input[0] ,symsyst[i], "q2") );
214  h2symsyst_MuKin_input.push_back( GetSymSyst(h2dsyst_MuKin_input[idsyst],h2dsyst_MuKin_input[idsyst+1],h2cv_MuKin_input[0],symsyst[i], "mukin") );
215  }
216  h1symsyst_ENu.push_back( GetWanted1DFromENu( h1symsyst_ENu_input[i], integral_syst_ENu, Form("henu_%s" ,symsyst[i].c_str())) );
217  h1symsyst_Q2.push_back( GetWanted1DFromQ2( h1symsyst_Q2_input[i], integral_syst_Q2, Form("hq2_%s" ,symsyst[i].c_str())) );
218  h1symsyst_MuKin.push_back( GetWanted1DFromMuKin( h2symsyst_MuKin_input[i], integral_syst_MuKin, Form("hmukin_%s" ,symsyst[i].c_str())) );
219  if(is_shape){
220  scaling(h1symsyst_ENu[i] ,integral_cv_ENu / integral_syst_ENu);
221  scaling(h1symsyst_Q2[i] ,integral_cv_Q2 / integral_syst_Q2);
222  scaling(h1symsyst_MuKin[i],integral_cv_MuKin/ integral_syst_MuKin);
223  }
224  }
225  //Det sym syst: 2 systematics
226  for(unsigned i=0;i<Nasymsyst;i++){
227  h1asymsyst_ENu_input.push_back( h1dsyst_ENu_input[10+i]);
228  h1asymsyst_Q2_input.push_back( h1dsyst_Q2_input[10+i]);
229  h2asymsyst_MuKin_input.push_back(h2dsyst_MuKin_input[10+i]);
230  h1asymsyst_ENu.push_back( GetWanted1DFromENu( h1asymsyst_ENu_input[i], integral_syst_ENu, Form("henu_%s" ,asymsyst[i].c_str())) );
231  h1asymsyst_Q2.push_back( GetWanted1DFromQ2( h1asymsyst_Q2_input[i], integral_syst_Q2, Form("hq2_%s" ,asymsyst[i].c_str())) );
232  h1asymsyst_MuKin.push_back( GetWanted1DFromMuKin( h2asymsyst_MuKin_input[i],integral_syst_MuKin, Form("hmukin_%s" ,asymsyst[i].c_str())) );
233  if(is_shape){
234  scaling(h1asymsyst_ENu[i] ,integral_cv_ENu / integral_syst_ENu);
235  scaling(h1asymsyst_Q2[i] ,integral_cv_Q2 / integral_syst_Q2);
236  scaling(h1asymsyst_MuKin[i],integral_cv_MuKin/ integral_syst_MuKin);
237  }
238  }
239 
240  //Making universes:
241  fOut->mkdir("check");
242  //Total
243  if(icomp==0){
244  MakeCov("enu" , fOut, h1cv_ENu[0] , h1genie_ENu , h1ppfx_ENu , h1symsyst_ENu , h1asymsyst_ENu , imode);
245  MakeCov("q2" , fOut, h1cv_Q2[0] , h1genie_Q2 , h1ppfx_Q2 , h1symsyst_Q2 , h1asymsyst_Q2 , imode);
246  MakeCov("mukin", fOut, h1cv_MuKin[0], h1genie_MuKin, h1ppfx_MuKin, h1symsyst_MuKin, h1asymsyst_MuKin, imode);
247  }
248  //Only GENIE and PPFX
249  if(icomp==1){
250  MakeCovUniverses("enu_genie" , fOut, h1cv_ENu[0] , h1genie_ENu , imode);
251  MakeCovUniverses("q2_gene" , fOut, h1cv_Q2[0] , h1genie_Q2 , imode);
252  MakeCovUniverses("mukin_genie", fOut, h1cv_MuKin[0], h1genie_MuKin, imode);
253  }
254  if(icomp==2){
255  MakeCovUniverses("enu_ppfx" , fOut, h1cv_ENu[0] , h1ppfx_ENu , imode);
256  MakeCovUniverses("q2_ppfx" , fOut, h1cv_Q2[0] , h1ppfx_Q2 , imode);
257  MakeCovUniverses("mukin_ppfx" , fOut, h1cv_MuKin[0], h1ppfx_MuKin , imode);
258  }
259  //Only symmetric systematics:
260  int local_idx = -1;
261  if(icomp>=3 && icomp<=18){
262  local_idx = icomp -3;
263  MakeCovSymSyst("enu",local_idx,fOut,h1cv_ENu[0] ,h1symsyst_ENu[local_idx] ,imode);
264  MakeCovSymSyst("q2",local_idx,fOut,h1cv_Q2[0] ,h1symsyst_Q2[local_idx] ,imode);
265  MakeCovSymSyst("mukin",local_idx,fOut,h1cv_MuKin[0],h1symsyst_MuKin[local_idx],imode);
266  }
267  //Only asymmetric systematics:
268  if(icomp>=19 && icomp<=20){
269  local_idx = icomp -19;
270  MakeCovASymSyst("enu",local_idx,fOut,h1cv_ENu[0] ,h1asymsyst_ENu[local_idx] ,imode);
271  MakeCovASymSyst("q2",local_idx,fOut,h1cv_Q2[0] ,h1asymsyst_Q2[local_idx] ,imode);
272  MakeCovASymSyst("mukin",local_idx,fOut,h1cv_MuKin[0],h1asymsyst_MuKin[local_idx],imode);
273  }
274 
275 
276  fOut->Close();
277  std::cout<<"Ending "<<std::endl;
278 }
const int Nbfoc
const int Ngenie
TH1D * GetWanted1DFromENu(TH1D *h1In, double &IntegralVal, std::string name)
void MakeCovASymSyst(std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInasym, int thismode)
const std::string dsyst[Ndsyst]
TH1D * GetSymSyst(TH1D *hInshiftp, TH1D *hInshiftn, TH1D *hIncv, std::string SystName, std::string VarName)
void MakeCov(std::string histName, TFile *fIn, TH1D *hIncv, std::vector< TH1D * > &hIngenie, std::vector< TH1D * > &hInppfx, std::vector< TH1D * > &hInsym, std::vector< TH1D * > &hInasym, int thismode)
void scaling(TH1D *hIn, const double shape_scale)
const std::string asymsyst[Nasymsyst]
const std::string cv[Ncv]
const int Nasymsyst
TH1D * GetWanted1DFromMuKin(TH2D *h2In, double &IntegralVal, std::string name)
TH1D * GetWanted1DFromQ2(TH1D *h1In, double &IntegralVal, std::string name)
OStream cout
Definition: OStream.cxx:6
const std::string symsyst[Nsymsyst]
std::vector< int > starting(const char *fnameIn, const char *fnameonly1DIn, const char *fnameOut, const char *cmode, const char *ccomp)
const int Nppfx
const int Ndsyst
void MakeCovSymSyst(std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInsym, int thismode)
const int Ncv
const int Nsymsyst
const std::string bfoc[Nbfoc]
void MakeCovUniverses(std::string VarName, TFile *fIn, TH1D *hIncv, std::vector< TH1D * > &hInuniv, int thismode)
TH1D * GetSymSyst ( TH1D *  hInshiftp,
TH1D *  hInshiftn,
TH1D *  hIncv,
std::string  SystName,
std::string  VarName 
)

Definition at line 374 of file Calculate_covariances.C.

References std::abs(), cv, and plot_xsec_1d::Nbins.

Referenced by Calculate_covariances().

374  {
375  TH1D* houtput = (TH1D*)hIncv->Clone(Form("hSymSist%s_%s",SystName.c_str(),VarName.c_str() ));
376  int Nbins = hIncv->GetXaxis()->GetNbins();
377  for(int ii=1;ii<=Nbins;ii++){
378  double pp = hInshiftp->GetBinContent(ii);
379  double nn = hInshiftn->GetBinContent(ii);
380  double cv = hIncv->GetBinContent(ii);
381  double diff_pp = std::abs(pp-cv);
382  double diff_nn = std::abs(nn-cv);
383  if(diff_pp>diff_nn)houtput->SetBinContent(ii,pp);
384  else houtput->SetBinContent(ii,2.*cv-nn);
385  }
386  return houtput;
387 }
float abs(float number)
Definition: d0nt_math.hpp:39
const std::string cv[Ncv]
TH2D * GetSymSyst ( TH2D *  hInshiftp,
TH2D *  hInshiftn,
TH2D *  hIncv,
std::string  SystName,
std::string  VarName 
)

Definition at line 389 of file Calculate_covariances.C.

References std::abs(), cv, plot_xsec_2d::NbinsX, and plot_xsec_2d::NbinsY.

389  {
390  TH2D* houtput = (TH2D*)hIncv->Clone(Form("hSymSist%s_%s",SystName.c_str(),VarName.c_str() ));
391  int NbinsX = hIncv->GetXaxis()->GetNbins();
392  int NbinsY = hIncv->GetYaxis()->GetNbins();
393  for(int ii=1;ii<=NbinsX;ii++){
394  for(int jj=1;jj<=NbinsY;jj++){
395  double pp = hInshiftp->GetBinContent(ii,jj);
396  double nn = hInshiftn->GetBinContent(ii,jj);
397  double cv = hIncv->GetBinContent(ii,jj);
398  double diff_pp = std::abs(pp-cv);
399  double diff_nn = std::abs(nn-cv);
400  if(diff_pp>diff_nn)houtput->SetBinContent(ii,jj,pp);
401  else houtput->SetBinContent(ii,jj,2.*cv-nn);
402  }
403  }
404  return houtput;
405 }
float abs(float number)
Definition: d0nt_math.hpp:39
const std::string cv[Ncv]
TH1D * GetWanted1DFromENu ( TH1D *  h1In,
double &  IntegralVal,
std::string  name 
)

Definition at line 299 of file Calculate_covariances.C.

References hadd_many_files::counter, and Nwanted_ENuBins.

Referenced by Calculate_covariances().

299  {
300  TH1D* h1Wanted_ENu = new TH1D(name.c_str(),"",Nwanted_ENuBins,0,Nwanted_ENuBins);
301  int counter = 1;
302  int Nbinsx = h1In->GetXaxis()->GetNbins();
303  IntegralVal = 0;
304  for(int ii=1;ii<=Nbinsx;ii++){
305  if(ii<=3 || ii==Nbinsx)continue;
306  IntegralVal += ( h1In->GetBinContent(ii) ) * ( h1In->GetXaxis()->GetBinWidth(ii) );
307  h1Wanted_ENu->SetBinContent(counter,h1In->GetBinContent(ii));
308  h1Wanted_ENu->SetBinError( counter,0);
309  counter++;
310  }
311  return h1Wanted_ENu;
312 }
const XML_Char * name
Definition: expat.h:151
const int Nwanted_ENuBins
TH1D * GetWanted1DFromMuKin ( TH2D *  h2In,
double &  IntegralVal,
std::string  name 
)

Definition at line 279 of file Calculate_covariances.C.

References hadd_many_files::counter, MuTperCosBins, and Nwanted_MuKinBins.

Referenced by Calculate_covariances().

279  {
280  TH1D* h1Wanted_MuKin = new TH1D(name.c_str(),"",Nwanted_MuKinBins,0,Nwanted_MuKinBins);
281  int counter = 1;
282  int NbinX = h2In->GetXaxis()->GetNbins();
283  int NbinY = h2In->GetYaxis()->GetNbins();
284  IntegralVal = 0;
285  for(int ii=1;ii<=NbinX;ii++){
286  if(ii==1)continue;
287  for(int jj=1;jj<=NbinY;jj++){
288  if(jj==1 || jj==NbinY)continue;
289  if( (jj-1) > MuTperCosBins[ii-2] )continue;
290  IntegralVal += ( h2In->GetBinContent(ii,jj) ) * ( h2In->GetXaxis()->GetBinWidth(ii) ) * ( h2In->GetYaxis()->GetBinWidth(jj) );
291  h1Wanted_MuKin->SetBinContent(counter,h2In->GetBinContent(ii,jj));
292  h1Wanted_MuKin->SetBinError( counter,0);
293  counter++;
294  }
295  }
296  return h1Wanted_MuKin;
297 }
const XML_Char * name
Definition: expat.h:151
const int MuTperCosBins[NMuCosBins]
const int Nwanted_MuKinBins
TH1D * GetWanted1DFromQ2 ( TH1D *  h1In,
double &  IntegralVal,
std::string  name 
)

Definition at line 313 of file Calculate_covariances.C.

References hadd_many_files::counter, and Nwanted_Q2Bins.

Referenced by Calculate_covariances().

313  {
314  TH1D* h1Wanted_Q2 = new TH1D(name.c_str(),"",Nwanted_Q2Bins,0,Nwanted_Q2Bins);
315  int counter = 1;
316  int Nbinsx = h1In->GetXaxis()->GetNbins();
317  IntegralVal = 0;
318  for(int ii=1;ii<=Nbinsx;ii++){
319  if(ii==Nbinsx)continue;
320  IntegralVal += ( h1In->GetBinContent(ii) ) * ( h1In->GetXaxis()->GetBinWidth(ii) );
321  h1Wanted_Q2->SetBinContent(counter,h1In->GetBinContent(ii));
322  h1Wanted_Q2->SetBinError( counter,0);
323  counter++;
324  }
325  return h1Wanted_Q2;
326 }
const XML_Char * name
Definition: expat.h:151
const int Nwanted_Q2Bins
int main ( int  argc,
const char *  argv[] 
)

Definition at line 689 of file Calculate_covariances.C.

References Calculate_covariances().

689  {
690  if(argc==1)Calculate_covariances();
691  Calculate_covariances(argv[1],argv[2],argv[3],argv[4],argv[5]);
692  return 0;
693 }
void Calculate_covariances()
void MakeCov ( std::string  histName,
TFile *  fIn,
TH1D *  hIncv,
std::vector< TH1D * > &  hIngenie,
std::vector< TH1D * > &  hInppfx,
std::vector< TH1D * > &  hInsym,
std::vector< TH1D * > &  hInasym,
int  thismode 
)

Definition at line 545 of file Calculate_covariances.C.

References std::abs(), AddOneSideShift(), AddShift(), baseasymsyst_rdm, basesymsyst_rdm, calc_corr_matrix(), calc_cov_matrix(), om::cout, allTimeWatchdog::endl, Nasymsyst, Ngenie, Nppfx, Nsymsyst, r3, plot_timing_data::swap, and Write().

Referenced by Calculate_covariances().

545  {
546  std::cout<<"=> Making covariances"<<std::endl;
547  //General:
548  std::vector<TH1D*> vout;
549  int univ_id = -1;
550  double this_shift = 0;
551 
552  for(int ii=0;ii<Ngenie;ii++){
553  for(int jj=0;jj<Nppfx;jj++){
554  univ_id++;
555  if(univ_id%5000==0)std::cout<<"=> Shifting var: "<<univ_id<<" from "<<(Nppfx*Ngenie)<<std::endl;
556  vout.push_back( (TH1D*)hIncv->Clone(Form("h%s_u%d",histName.c_str(),univ_id)) );
557  //genie
558  vout[univ_id]->Add(hIngenie[ii]);
559  vout[univ_id]->Add(hIncv,-1);
560  //ppfx
561  if(thismode == 0 || thismode == 1){
562  if(ii==0 && jj==0)std::cout<<"=> Doing ppfx here"<<std::endl;
563  vout[univ_id]->Add(hInppfx[jj]);
564  vout[univ_id]->Add(hIncv,-1);
565  }
566  //sym syst
567  for(unsigned kk=0;kk<Nsymsyst;kk++){
568  r3->SetSeed(basesymsyst_rdm[kk]+univ_id);
569  this_shift = r3->Gaus(0.0,1.0);
570  AddShift(vout[univ_id],this_shift,hIncv, hInsym[kk]);
571  }
572  //one side syst cal shape:
573  for(unsigned kk=0;kk<Nasymsyst;kk++){
574  r3->SetSeed(baseasymsyst_rdm[kk]+univ_id);
575  this_shift = std::abs(r3->Gaus(0.0,1.0));
576  AddOneSideShift(vout[univ_id],this_shift,hIncv, hInasym[kk]);
577  }
578  }
579  } //end of new universes
580 
581  std::cout<<"We generated "<<vout.size()<< " of new universes for "<<histName<<std::endl;
582  TMatrixD cov = calc_cov_matrix(vout,hIncv);
583  TMatrixD cor = calc_corr_matrix(cov);
584  int nofbins = hIncv->GetXaxis()->GetNbins();
585 
586  TH2D* hcov = new TH2D(Form("hcov_%s",histName.c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
587  TH2D* hcor = new TH2D(Form("hcor_%s",histName.c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
588  for(int ii=1;ii<=nofbins;ii++){
589  for(int jj=1;jj<=nofbins;jj++){
590  hcov->SetBinContent(ii,jj,cov[ii-1][jj-1]);
591  hcor->SetBinContent(ii,jj,cor[ii-1][jj-1]);
592  }
593  }
594  std::cout<<"Storing results for "<<histName<<std::endl;
595  fIn->cd();
596  hIncv->Write(Form("hcv_%s",histName.c_str()));
597  hcov->Write( Form("cov_%s",histName.c_str()));
598  hcor->Write( Form("cor_%s",histName.c_str()));
599 
600  //storing sym systematics
601  fIn->mkdir(("check_"+histName).c_str());
602  fIn->cd(("check_"+histName).c_str());
603  for(int ii=0;ii<hInsym.size();ii++ )hInsym[ii]->Write();
604  for(int ii=0;ii<hInasym.size();ii++)hInasym[ii]->Write();
605 
606  //Clearing vectors:
607  hIngenie.clear();
608  hInppfx.clear();
609  hInsym.clear();
610  hInasym.clear();
611  vout.clear();
612  std::vector<TH1D*>().swap(vout);
613  vout.shrink_to_fit();
614 }
const int basesymsyst_rdm[Nsymsyst]
const int Ngenie
TMatrixD calc_cov_matrix(std::vector< TH1D * > &vhist, TH1D *hIncv)
const int baseasymsyst_rdm[Nasymsyst]
void AddOneSideShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
TRandom3 * r3
float abs(float number)
Definition: d0nt_math.hpp:39
const int Nasymsyst
OStream cout
Definition: OStream.cxx:6
const int Nppfx
TMatrixD calc_corr_matrix(TMatrixD &mcov)
void AddShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
const int Nsymsyst
gm Write()
void MakeCovASymSyst ( std::string  varName,
int  iVar,
TFile *  fIn,
TH1D *  hIncv,
TH1D *  hInasym,
int  thismode 
)

Definition at line 510 of file Calculate_covariances.C.

References std::abs(), AddOneSideShift(), asymsyst, baseasymsyst_rdm, calc_corr_matrix(), calc_cov_matrix(), om::cout, allTimeWatchdog::endl, Nuniv, and r3.

Referenced by Calculate_covariances().

510  {
511  if(thismode==2)return; //this calculation is not for "not ppfx mode"
512  std::cout<<"=> Making covariances asymmetric systematics: "<< varName<<", "<<asymsyst[iSyst] <<", mode: "<< thismode <<std::endl;
513 
514  //General:
515  std::vector<TH1D*> vout;
516  int univ_id = -1;
517  double this_shift = 0;
518  for(int ii=0;ii<Nuniv;ii++){
519  univ_id++;
520  vout.push_back( (TH1D*)hIncv->Clone(Form("h%s_%s_u%d",varName.c_str(),asymsyst[iSyst].c_str(),univ_id)) );
521  r3->SetSeed(baseasymsyst_rdm[iSyst]+univ_id);
522  this_shift = std::abs(r3->Gaus(0.0,1.0));
523  AddOneSideShift(vout[ii],this_shift,hIncv, hInasym);
524  } //end of new universes
525 
526  std::cout<<"We generated "<<vout.size()<< " of new universes for "<<varName<<", "<<asymsyst[iSyst]<<std::endl;
527  TMatrixD cov = calc_cov_matrix(vout,hIncv);
528  TMatrixD cor = calc_corr_matrix(cov);
529  int nofbins = hIncv->GetXaxis()->GetNbins();
530 
531  TH2D* hcov = new TH2D(Form("hcov_%s_%s",varName.c_str(),asymsyst[iSyst].c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
532  TH2D* hcor = new TH2D(Form("hcor_%s_%s",varName.c_str(),asymsyst[iSyst].c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
533  for(int ii=1;ii<=nofbins;ii++){
534  for(int jj=1;jj<=nofbins;jj++){
535  hcov->SetBinContent(ii,jj,cov[ii-1][jj-1]);
536  hcor->SetBinContent(ii,jj,cor[ii-1][jj-1]);
537  }
538  }
539  std::cout<<"Storing results for "<<asymsyst[iSyst]<<std::endl;
540  fIn->cd();
541  hcov->Write( Form("cov_%s_%s",varName.c_str(),asymsyst[iSyst].c_str()));
542  hcor->Write( Form("cor_%s_%s",varName.c_str(),asymsyst[iSyst].c_str()));
543 }
TMatrixD calc_cov_matrix(std::vector< TH1D * > &vhist, TH1D *hIncv)
const int baseasymsyst_rdm[Nasymsyst]
void AddOneSideShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
TRandom3 * r3
float abs(float number)
Definition: d0nt_math.hpp:39
const std::string asymsyst[Nasymsyst]
OStream cout
Definition: OStream.cxx:6
const int Nuniv
TMatrixD calc_corr_matrix(TMatrixD &mcov)
void MakeCovSymSyst ( std::string  varName,
int  iVar,
TFile *  fIn,
TH1D *  hIncv,
TH1D *  hInsym,
int  thismode 
)

Definition at line 467 of file Calculate_covariances.C.

References AddShift(), basesymsyst_rdm, calc_corr_matrix(), calc_cov_matrix(), calc_fer_matrix(), om::cout, allTimeWatchdog::endl, Nuniv, r3, and symsyst.

Referenced by Calculate_covariances().

467  {
468  if(thismode==2)return; //this calculation is not for the "not ppfx mode"
469  std::cout<<"=> Making covariances symmetric systematics: "<< varName<<", "<<symsyst[iSyst]<<", mode: "<< thismode <<std::endl;
470 
471  //General:
472  std::vector<TH1D*> vout;
473  int univ_id = -1;
474  double this_shift = 0;
475 
476  for(int ii=0;ii<Nuniv;ii++){
477  univ_id++;
478  vout.push_back( (TH1D*)hIncv->Clone(Form("h%s_%s_u%d",varName.c_str(),symsyst[iSyst].c_str(),univ_id)) );
479  r3->SetSeed(basesymsyst_rdm[iSyst]+univ_id);
480  this_shift = r3->Gaus(0.0,1.0);
481  AddShift(vout[ii],this_shift,hIncv, hInsym);
482  } //end of new universes
483 
484  std::cout<<"We generated "<<vout.size()<< " of new universes for "<<varName<<", "<<symsyst[iSyst]<<std::endl;
485  TMatrixD cov = calc_cov_matrix(vout,hIncv);
486  TMatrixD cor = calc_corr_matrix(cov);
487  TMatrixD fer = calc_fer_matrix(cov,hIncv);
488  int nofbins = hIncv->GetXaxis()->GetNbins();
489 
490  TH2D* hcov = new TH2D(Form("hcov_%s_%s",varName.c_str(),symsyst[iSyst].c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
491  TH2D* hcor = new TH2D(Form("hcor_%s_%s",varName.c_str(),symsyst[iSyst].c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
492  TH2D* hfer = new TH2D(Form("hfer_%s_%s",varName.c_str(),symsyst[iSyst].c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
493  for(int ii=1;ii<=nofbins;ii++){
494  for(int jj=1;jj<=nofbins;jj++){
495  hcov->SetBinContent(ii,jj,cov[ii-1][jj-1]);
496  hcor->SetBinContent(ii,jj,cor[ii-1][jj-1]);
497  hfer->SetBinContent(ii,jj,fer[ii-1][jj-1]);
498  }
499  }
500  std::cout<<"Storing results"<<std::endl;
501  fIn->cd();
502  hcov->Write( Form("cov_%s_%s",varName.c_str(),symsyst[iSyst].c_str()));
503  hcor->Write( Form("cor_%s_%s",varName.c_str(),symsyst[iSyst].c_str()));
504  fIn->cd("check");
505  hInsym->Write();
506  hIncv->Write();
507  hfer->Write();
508 }
const int basesymsyst_rdm[Nsymsyst]
TMatrixD calc_cov_matrix(std::vector< TH1D * > &vhist, TH1D *hIncv)
TRandom3 * r3
TMatrixD calc_fer_matrix(TMatrixD &mcov, TH1D *hIncv)
OStream cout
Definition: OStream.cxx:6
const std::string symsyst[Nsymsyst]
const int Nuniv
TMatrixD calc_corr_matrix(TMatrixD &mcov)
void AddShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
void MakeCovUniverses ( std::string  VarName,
TFile *  fIn,
TH1D *  hIncv,
std::vector< TH1D * > &  hInuniv,
int  thismode 
)

Definition at line 435 of file Calculate_covariances.C.

References calc_corr_matrix(), calc_cov_matrix(), om::cout, allTimeWatchdog::endl, and makeTrainCVSamples::int.

Referenced by Calculate_covariances().

435  {
436  if(thismode==2)return; //this calculation is not for the "not ppfx mode"
437  std::cout<<"=> Making covariances universes: "<< histName<<", mode: "<< thismode <<std::endl;
438 
439  std::vector<TH1D*> vout;
440  int univ_id = -1;
441  int NNuniv = int(hInuniv.size());
442  for(int ii=0;ii<NNuniv;ii++){
443  univ_id++;
444  vout.push_back( (TH1D*)hIncv->Clone(Form("h%s_u%d",histName.c_str(),univ_id)) );
445  vout[ii]->Add(hInuniv[ii]);
446  vout[ii]->Add(hIncv,-1);
447  } //end of new universes
448 
449  std::cout<<"We generated "<<vout.size()<< " of new universes for "<<histName<<std::endl;
450  TMatrixD cov = calc_cov_matrix(vout,hIncv);
451  TMatrixD cor = calc_corr_matrix(cov);
452  int nofbins = hIncv->GetXaxis()->GetNbins();
453 
454  TH2D* hcov = new TH2D(Form("hcov_%s",histName.c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
455  TH2D* hcor = new TH2D(Form("hcor_%s",histName.c_str()),"",nofbins,0,nofbins,nofbins,0,nofbins);
456  for(int ii=1;ii<=nofbins;ii++){
457  for(int jj=1;jj<=nofbins;jj++){
458  hcov->SetBinContent(ii,jj,cov[ii-1][jj-1]);
459  hcor->SetBinContent(ii,jj,cor[ii-1][jj-1]);
460  }
461  }
462  std::cout<<"Storing results for "<<histName<<std::endl;
463  fIn->cd();
464  hcov->Write( Form("cov_%s",histName.c_str()));
465  hcor->Write( Form("cor_%s",histName.c_str()));
466 }
TMatrixD calc_cov_matrix(std::vector< TH1D * > &vhist, TH1D *hIncv)
OStream cout
Definition: OStream.cxx:6
TMatrixD calc_corr_matrix(TMatrixD &mcov)
void scaling ( TH1D *  hIn,
const double  shape_scale 
)
std::vector< int > starting ( const char *  fnameIn,
const char *  fnameonly1DIn,
const char *  fnameOut,
const char *  cmode,
const char *  ccomp 
)

Definition at line 631 of file Calculate_covariances.C.

References asymsyst, Calculate_covariances(), om::cout, allTimeWatchdog::endl, Nasymsyst, Nsymsyst, string, and symsyst.

Referenced by Calculate_covariances().

631  {
632 
633  std::vector<int> voption(2,-1);
634  std::cout<<"Starting Calculate_covariances()"<<std::endl;
635  std::cout<<"Input:"<<std::endl;
636  std::cout<<"file name In (2D): "<< fnameIn <<std::endl;
637  std::cout<<"file name In (1D): "<< fnameonly1DIn<<std::endl;
638  std::cout<<"file name Out : "<< fnameOut <<std::endl;
639  std::cout<<"mode : "<< cmode <<std::endl;
640  std::cout<<"component : "<< ccomp <<std::endl;
641 
642  //mode
643  if(std::string(cmode)=="total"){
644  std::cout<<"Calculating total covariance"<<std::endl;
645  voption[0] = 0;
646  }
647  else if(std::string(cmode)=="shape"){
648  std::cout<<"Calculating shape covariance"<<std::endl;
649  voption[0] = 1;
650  }
651  else if(std::string(cmode)=="notppfx"){
652  std::cout<<"Calculating covariance without ppfx"<<std::endl;
653  voption[0] = 2;
654  }
655  else{
657  }
658  //component
659  if(std::string(ccomp)=="all"){
660  std::cout<<"Calculating all components"<<std::endl;
661  voption[1] = 0;
662  }
663  if(std::string(ccomp)=="GENIE"){
664  std::cout<<"Calculating "<<ccomp<<" systematic"<<std::endl;
665  voption[1] = 1;
666  }
667  if(std::string(ccomp)=="PPFX"){
668  std::cout<<"Calculating "<<ccomp<<" systematic"<<std::endl;
669  voption[1] = 2;
670  }
671  for(int ii=0;ii<Nsymsyst;ii++){
672  if(std::string(ccomp) == symsyst[ii]){
673  std::cout<<"Calculating "<<ccomp<<" systematic"<<std::endl;
674  voption[1] = 3 + ii;
675  }
676  }
677  for(int ii=0;ii<Nasymsyst;ii++){
678  if(std::string(ccomp) == asymsyst[ii]){
679  std::cout<<"Calculating "<<ccomp<<" systematic"<<std::endl;
680  voption[1] = 19 + ii;
681  }
682  }
683  if(voption[1] == -1)Calculate_covariances();
684 
685  return voption;
686 }
void Calculate_covariances()
const std::string asymsyst[Nasymsyst]
const int Nasymsyst
OStream cout
Definition: OStream.cxx:6
const std::string symsyst[Nsymsyst]
const int Nsymsyst
enum BeamMode string

Variable Documentation

const std::string asymsyst[Nasymsyst] = {"CalibShape","Cherenkov"}

Definition at line 63 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), MakeCovASymSyst(), and starting().

const int baseasymsyst_rdm[Nasymsyst] = {26000000,27000000}

Definition at line 64 of file Calculate_covariances.C.

Referenced by MakeCov(), and MakeCovASymSyst().

const int basesymsyst_rdm[Nsymsyst]
Initial value:
= {10000000,11000000,12000000,13000000,14000000,
15000000,16000000,17000000,18000000,19000000,
20000000,21000000,22000000,23000000,24000000,
25000000}

Definition at line 56 of file Calculate_covariances.C.

Referenced by MakeCov(), and MakeCovSymSyst().

const std::string bfoc[Nbfoc]
Initial value:
= {"2kA_Up","2kA_Dw",
"02mmBeamSpotSize_Up","02mmBeamSpotSize_Dw",
"1mmBeamShiftX_Up","1mmBeamShiftX_Dw",
"1mmBeamShiftY_Up","1mmBeamShiftY_Dw",
"3mmHorn1X_Up","3mmHorn1X_Dw",
"3mmHorn1Y_Up","3mmHorn1Y_Dw",
"3mmHorn2X_Up","3mmHorn2X_Dw",
"3mmHorn2Y_Up","3mmHorn2Y_Dw",
"7mmTargetZ_Up","7mmTargetZ_Dw",
"MagneticFieldinDecayPipe_Up","MagneticFieldinDecayPipe_Dw",
"1mmHornWater_Up","1mmHornWater_Dw"}

Definition at line 19 of file Calculate_covariances.C.

Referenced by Calculate_covariances().

const std::string cv[Ncv] = {"cv","nom_good_seventh","neutron_nom"}
const std::string dsyst[Ndsyst]
Initial value:
= {"angleup","angledw",
"calibpos","calibneg",
"lightup","lightdown",
"muoneup","muonedw",
"neutron_Up","neutron_Dw",
"calibshift",
"ckv"}

Definition at line 31 of file Calculate_covariances.C.

Referenced by Calculate_covariances().

const int MuTperCosBins[NMuCosBins] = {6,7,8,9,9,9,13,14,17,20,20,20,20}

Definition at line 44 of file Calculate_covariances.C.

Referenced by GetWanted1DFromMuKin().

const int Nasymsyst = 2

Definition at line 62 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), MakeCov(), and starting().

const int Nbfoc = 22

Definition at line 18 of file Calculate_covariances.C.

Referenced by Calculate_covariances().

const int Ncv = 3

Definition at line 38 of file Calculate_covariances.C.

Referenced by Calculate_covariances().

const int Ndsyst = 12

Definition at line 30 of file Calculate_covariances.C.

Referenced by Calculate_covariances().

const int Ngenie = 1000

Definition at line 16 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), and MakeCov().

const int NMuCosBins = 13

Definition at line 43 of file Calculate_covariances.C.

const int Nppfx = 100

Definition at line 15 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), and MakeCov().

const int Nsymsyst = 16

Definition at line 50 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), MakeCov(), and starting().

const int Nuniv = Ngenie*Nppfx
const int Nwanted_ENuBins = 18

Definition at line 46 of file Calculate_covariances.C.

Referenced by GetWanted1DFromENu().

const int Nwanted_MuKinBins = 172

Definition at line 42 of file Calculate_covariances.C.

Referenced by GetWanted1DFromMuKin().

const int Nwanted_Q2Bins = 15

Definition at line 47 of file Calculate_covariances.C.

Referenced by GetWanted1DFromQ2().

TRandom3* r3 = new TRandom3(0)
const std::string symsyst[Nsymsyst]
Initial value:
= {"FocCurrMag","FocBeamSpotSize","FocBeamShiftX","FocBeamShiftY",
"FocHorn1X", "FocHorn1Y", "FocHorn2X", "FocHorn2Y",
"FocTargetZ", "FocBFieldDPIP", "FocHornWaterIC",
"MuAngle", "CalibShift", "light", "MuE","Neu"}

Definition at line 51 of file Calculate_covariances.C.

Referenced by Calculate_covariances(), MakeCovSymSyst(), and starting().