Calculate_covariances.C
Go to the documentation of this file.
1 
2 
3 #include "TFile.h"
4 #include "TH1D.h"
5 #include "TH2D.h"
6 #include <TRandom.h>
7 #include <TMatrixD.h>
8 #include <TRandom3.h>
9 
10 #include <vector>
11 #include <iostream>
12 #include <string>
13 
14 //Constants:
15 const int Nppfx = 100;
16 const int Ngenie = 1000;
17 const int Nuniv = Ngenie*Nppfx;
18 const int Nbfoc = 22;
19 const std::string bfoc[Nbfoc] = {"2kA_Up","2kA_Dw",
20  "02mmBeamSpotSize_Up","02mmBeamSpotSize_Dw",
21  "1mmBeamShiftX_Up","1mmBeamShiftX_Dw",
22  "1mmBeamShiftY_Up","1mmBeamShiftY_Dw",
23  "3mmHorn1X_Up","3mmHorn1X_Dw",
24  "3mmHorn1Y_Up","3mmHorn1Y_Dw",
25  "3mmHorn2X_Up","3mmHorn2X_Dw",
26  "3mmHorn2Y_Up","3mmHorn2Y_Dw",
27  "7mmTargetZ_Up","7mmTargetZ_Dw",
28  "MagneticFieldinDecayPipe_Up","MagneticFieldinDecayPipe_Dw",
29  "1mmHornWater_Up","1mmHornWater_Dw"};
30 const int Ndsyst = 12;
31 const std::string dsyst[Ndsyst] = {"angleup","angledw",
32  "calibpos","calibneg",
33  "lightup","lightdown",
34  "muoneup","muonedw",
35  "neutron_Up","neutron_Dw",
36  "calibshift",
37  "ckv"};
38 const int Ncv = 3;
39 const std::string cv[Ncv] = {"cv","nom_good_seventh","neutron_nom"};
40 
41 
42 const int Nwanted_MuKinBins = 172;
43 const int NMuCosBins = 13;
44 const int MuTperCosBins[NMuCosBins] = {6,7,8,9,9,9,13,14,17,20,20,20,20};
45 
46 const int Nwanted_ENuBins = 18;
47 const int Nwanted_Q2Bins = 15;
48 
49 //sysm syst
50 const int Nsymsyst = 16;
51 const std::string symsyst[Nsymsyst] = {"FocCurrMag","FocBeamSpotSize","FocBeamShiftX","FocBeamShiftY",
52  "FocHorn1X", "FocHorn1Y", "FocHorn2X", "FocHorn2Y",
53  "FocTargetZ", "FocBFieldDPIP", "FocHornWaterIC",
54  "MuAngle", "CalibShift", "light", "MuE","Neu"};
55 
56 const int basesymsyst_rdm[Nsymsyst] = {10000000,11000000,12000000,13000000,14000000,
57  15000000,16000000,17000000,18000000,19000000,
58  20000000,21000000,22000000,23000000,24000000,
59  25000000};
60 
61 //One side syst
62 const int Nasymsyst = 2;
63 const std::string asymsyst[Nasymsyst] = {"CalibShape","Cherenkov"};
64 const int baseasymsyst_rdm[Nasymsyst] = {26000000,27000000};
65 
66 TH1D* GetWanted1DFromENu(TH1D* h1In, double& IntegralVal, std::string name);
67 TH1D* GetWanted1DFromQ2(TH1D* h1In , double& IntegralVal, std::string name);
68 TH1D* GetWanted1DFromMuKin(TH2D* h2In, double& IntegralVal, std::string name);
69 TMatrixD calc_cov_matrix(std::vector<TH1D*>& vhist, TH1D* hIncv);
71 TMatrixD calc_fer_matrix(TMatrixD& mcov, TH1D* hIncv);
72 TH1D* GetSymSyst(TH1D* hInshiftp,TH1D* hInshiftn,TH1D* hIncv,std::string SystName,std::string VarName);
73 TH2D* GetSymSyst(TH2D* hInshiftp,TH2D* hInshiftn,TH2D* hIncv,std::string SystName,std::string VarName);
74 void AddShift( TH1D* hUniv,double& shift,TH1D* hIncv, TH1D* hSystIn);
75 void AddOneSideShift(TH1D* hUniv,double& shift,TH1D* hIncv, TH1D* hSystIn);
76 void scaling(TH1D* hIn, const double shape_scale);
77 
78 void MakeCov(std::string histName, TFile* fIn, TH1D* hIncv, std::vector<TH1D*>& hIngenie, std::vector<TH1D*>& hInppfx,
79  std::vector<TH1D*>& hInsym, std::vector<TH1D*>& hInasym, int thismode);
80 void MakeCovUniverses(std::string VarName, TFile* fIn, TH1D* hIncv, std::vector<TH1D*>& hInuniv, int thismode);
81 void MakeCovSymSyst( std::string varName, int iVar, TFile* fIn, TH1D* hIncv, TH1D* hInsym , int thismode);
82 void MakeCovASymSyst( std::string varName, int iVar, TFile* fIn, TH1D* hIncv, TH1D* hInasym , int thismode);
83 
84 //General variables
85 TRandom3 *r3 = new TRandom3(0);
86 //
88 std::vector<int> starting(const char* fnameIn, const char* fnameonly1DIn,const char* fnameOut, const char* cmode, const char* ccomp);
89 
90 //Compiling: g++ -o Calculate_covariances Calculate_covariances.C `root-config --cflags --glibs`
91 void Calculate_covariances(const char* fnameIn, const char* fnameonly1DIn,const char* fnameOut, const char* cmode, const char* ccomp){
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 }
279 TH1D* GetWanted1DFromMuKin(TH2D* h2In, double& IntegralVal, std::string name){
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 }
298 
299 TH1D* GetWanted1DFromENu(TH1D* h1In, double& IntegralVal, std::string name){
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 }
313 TH1D* GetWanted1DFromQ2(TH1D* h1In, double& IntegralVal, std::string name){
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 }
327 TMatrixD calc_cov_matrix(std::vector<TH1D*>& vhist, TH1D* hIncv){
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 }
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 }
362 TMatrixD calc_fer_matrix(TMatrixD& mcov, TH1D* hIncv){
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 }
373 
374 TH1D* GetSymSyst(TH1D* hInshiftp,TH1D* hInshiftn,TH1D* hIncv,std::string SystName,std::string VarName){
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 }
388 
389 TH2D* GetSymSyst(TH2D* hInshiftp,TH2D* hInshiftn,TH2D* hIncv,std::string SystName,std::string VarName){
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 }
406 
407 void AddShift(TH1D* hUniv,double& shift,TH1D* hIncv, TH1D* hSystIn){
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 }
419 void AddOneSideShift(TH1D* hUniv,double& shift,TH1D* hIncv, TH1D* hSystIn){
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 }
431 void scaling(TH1D* hIn, const double shape_scale){
432  hIn->Scale(shape_scale);
433 }
434 
435 void MakeCovUniverses(std::string histName, TFile* fIn, TH1D* hIncv, std::vector<TH1D*>& hInuniv, int thismode){
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 }
467 void MakeCovSymSyst( std::string varName, int iSyst, TFile* fIn, TH1D* hIncv, TH1D* hInsym, int thismode){
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 }
509 
510 void MakeCovASymSyst(std::string varName,int iSyst, TFile* fIn, TH1D* hIncv, TH1D* hInasym, int thismode){
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 }
544 
545 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){
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 }
615 
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 }
630 
631 std::vector<int> starting(const char* fnameIn, const char* fnameonly1DIn,const char* fnameOut, const char* cmode, const char* ccomp){
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 }
687 
688 # ifndef __CINT__
689 int main(int argc, const char* argv[]){
690  if(argc==1)Calculate_covariances();
691  Calculate_covariances(argv[1],argv[2],argv[3],argv[4],argv[5]);
692  return 0;
693 }
694 # endif
const int Nbfoc
const int basesymsyst_rdm[Nsymsyst]
const XML_Char * name
Definition: expat.h:151
const int MuTperCosBins[NMuCosBins]
const int Ngenie
TMatrixD calc_cov_matrix(std::vector< TH1D * > &vhist, TH1D *hIncv)
TH1D * GetWanted1DFromENu(TH1D *h1In, double &IntegralVal, std::string name)
const int baseasymsyst_rdm[Nasymsyst]
T sqrt(T number)
Definition: d0nt_math.hpp:156
int main(int argc, const char *argv[])
void AddOneSideShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
void MakeCovASymSyst(std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInasym, int thismode)
const std::string dsyst[Ndsyst]
void Calculate_covariances()
TRandom3 * r3
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)
float abs(float number)
Definition: d0nt_math.hpp:39
void scaling(TH1D *hIn, const double shape_scale)
const std::string asymsyst[Nasymsyst]
TMatrixD calc_fer_matrix(TMatrixD &mcov, TH1D *hIncv)
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]
const int Nwanted_ENuBins
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const int NMuCosBins
const int Nuniv
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
TMatrixD calc_corr_matrix(TMatrixD &mcov)
void MakeCovSymSyst(std::string varName, int iVar, TFile *fIn, TH1D *hIncv, TH1D *hInsym, int thismode)
const int Nwanted_Q2Bins
exit(0)
const int Ncv
void AddShift(TH1D *hUniv, double &shift, TH1D *hIncv, TH1D *hSystIn)
const int Nwanted_MuKinBins
const int Nsymsyst
const std::string bfoc[Nbfoc]
gm Write()
void MakeCovUniverses(std::string VarName, TFile *fIn, TH1D *hIncv, std::vector< TH1D * > &hInuniv, int thismode)