NDFDCVNCutSelector.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Sample.h"
3 #include "CAFAna/Core/Registry.h"
5 #include "CAFAna/Cuts/QuantileCuts.h"
11 
14 
15 #include "NuXAna/Cuts/NusCuts18.h"
16 #include "NuXAna/Cuts/NusCuts20.h"
23 
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TMatrixT.h"
27 
28 // print bin contents of histogram
29 // ---------------------------------------------------------------------------
30 void PrintHist(TH1D* h){
31  std::cout << "hist name: " << h->GetName() << std::endl;
32  for (int i = 0; i < h->GetNbinsX(); ++i)
33  std::cout << h->GetBinContent(i+1) << std::endl;
34 }
35 
36 // apply projection to get selected events for a given cut value
37 //-----------------------------------------------------------------------------
38 TH1D* GetSelectedEvents(TH2D* h, std::string det, int cutBin){
39 
40  std::string histName = det+"proj"+std::to_string(cutBin);
41  TH1D* proj = h->ProjectionX(histName.c_str(), cutBin+1, h->GetNbinsY()+1);
42 
43  return proj;
44 }
45 
46 // concatenate the nd and fd histograms
47 // ----------------------------------------------------------------------------
48 TH1D* ConcatenateHists(TH1D* nd, TH1D* fd){
49 
50  int ndBins = nd->GetNbinsX();
51  int fdBins = fd->GetNbinsX();
52  int totBins = ndBins + fdBins;
53 
54  TH1D* h = new TH1D("allE", ";Energy Bin;", totBins, 0, totBins);
55 
56  for (int i = 0; i < totBins; ++i){
57  if (i < ndBins)
58  h->SetBinContent(i+1, nd->GetBinContent(i+1));
59  else
60  h->SetBinContent(i+1, fd->GetBinContent(i+1-ndBins));
61  }
62 
63  //PrintHist(nd);
64  //PrintHist(fd);
65  //PrintHist(h);
66 
67  return h;
68 }
69 
70 // function to find and remove empty rows in matrix
71 // ---------------------------------------------------------------------------
73 
74  // find empty rows and columns
75  std::vector<int> emptyRows;
76  for (int i = 0; i < mat->GetNrows();++i){
77  float rowTotal = 0;
78  for (int j = 0; j < mat->GetNcols(); ++j){
79  rowTotal += (*mat)(i,j);
80  }
81  if (rowTotal == 0)
82  emptyRows.push_back(i);
83  }
84  std::vector<int> emptyCols;
85  for (int j = 0; j < mat->GetNcols();++j){
86  float colTotal = 0;
87  for (int i = 0; i < mat->GetNrows(); ++i){
88  colTotal += (*mat)(i,j);
89  }
90  if (colTotal == 0)
91  emptyCols.push_back(j);
92  }
93 
94  std::vector<int>::iterator it;
95  TMatrixD out(mat->GetNrows()-emptyRows.size(), mat->GetNcols()-emptyCols.size());
96  int posi = 0;
97  for (int i = 0; i < mat->GetNrows(); ++i){
98  it = std::find(emptyRows.begin(), emptyRows.end(), i);
99  if (it == emptyRows.end()) continue;
100 
101  int posj = 0;
102  for (int j = 0; j < mat->GetNrows(); ++j){
103  it = std::find(emptyCols.begin(), emptyCols.end(), i);
104  if (it == emptyCols.end()) continue;
105 
106  out(posi, posj) = (*mat)(i,j);
107 
108  posj++;
109  }
110  posi++;
111  }
112 
113  return out;
114 }
115 
116 // covnert fractional covariance matrix to full covariance matrix
117 // including the statistical uncertainty
118 // ----------------------------------------------------------------------------
119 TMatrixD FracToFull(TMatrixD* fracCov, TH1D* spectra){
120 
121  TMatrixD inMat = *fracCov;
122 
123  int nBins = spectra->GetNbinsX();
124  TMatrixD fullCov(nBins, nBins);
125 
126  for (int i =0; i < nBins; ++i){
127  for (int j = 0; j < nBins; ++j){
128 
129  // multiply out systematic error
130  fullCov(i,j) = inMat(i,j) *
131  spectra->GetBinContent(i+1) *
132  spectra->GetBinContent(j+1);
133 
134  // add in variance from statistical uncertainties
135  if (i == j)
136  fullCov(i,j) = fullCov(i,j) + spectra->GetBinContent(i+1);
137 
138  }
139  }
140 
141  //fullCov.Print();
142 
143  return fullCov;
144 
145 }
146 
147 // main function
148 // ----------------------------------------------------------------------------
150 
151  // load covariance matrix
152  TFile* covmxFile = new TFile("/pnfs/nova/persistent/analysis/nux/nus20/selection/covmx.root");
153 
154  TMatrixD* fraccovmx = (TMatrixD*)covmxFile->Get("fraccovmx");
155 
156  std::string simEnergyStringFD = "FD/kCaloEkCVNnc_looseptpPreselQualityFiducialTotalMC";
157  std::string simEnergyStringND = "ND/kCaloEkCVNnc_looseptpPreselQualityFiducialTotalMC";
158 
159  TFile* spectraFile = new TFile(file.c_str(), "read");
160 
161  TH2D* hCVNEnergyND
162  = (TH2D*)spectraFile->Get(simEnergyStringND.c_str());
163 
164  TH2D* hCVNEnergyFD
165  = (TH2D*)spectraFile->Get(simEnergyStringFD.c_str());
166 
167  // now need to scale to some POT
168 
169 
170  TH2D* fomVals
171  = new TH2D("fomVals", ";ND CVN cut value;FD CVN cut value", 100, 0, 1, 100, 0, 1);
172 
173  TH1D* ndE; TH1D* fdE; TH1D* allE;
174 
175  for (int iND = 0; iND < 100; ++iND){
176  ndE = GetSelectedEvents(hCVNEnergyND, "nd", iND);
177 
178  //for (int j = 0; j < ndE->GetNbinsX(); ++j)
179  // ndE->SetBinContent(j+1,0);
180 
181  for (int iFD = 0; iFD < 100; ++iFD){
182  fdE = GetSelectedEvents(hCVNEnergyFD, "fd", iFD);
183 
184  // concatenate histograms and get the full matrix
185  allE = ConcatenateHists(ndE, fdE);
186 
187  TMatrixD fullcovmx = FracToFull(fraccovmx, allE);
188 
189  //TMatrixD zmfullcov = RemoveEmptyRowsAndColumns(&fullcovmx);
190 
191  //fullcovmx.Print();
192 
193  std::vector<double> fulldecorrelatedunc =
194  GetDecorrelatedUncertainty(&fullcovmx,
195  ndE->GetNbinsX(),
196  fdE->GetNbinsX());
197 
198  double fomVal = calcBinnedSOverSSB(allE, fulldecorrelatedunc);
199 
200  fomVals->SetBinContent(iND+1,iFD+1, fomVal);
201 
202  delete allE;
203 
204  }
205  }
206 
207  double val = 0;
208  float selBinI = 0;
209  float selBinJ = 0;
210  for (int i = 0; i < fomVals->GetNbinsX(); ++i){
211  for (int j = 0; j < fomVals->GetNbinsY(); ++j){
212  //if (i > 95){
213  // std::cout
214  // << fomVals->GetXaxis()->GetBinLowEdge(i+1) << ", "
215  // << fomVals->GetYaxis()->GetBinLowEdge(j+1) << ": "
216  // << fomVals->GetBinContent(i+1, j+1) << std::endl;
217  //}
218  if (fomVals->GetBinContent(i+1,j+1) > val){
219  val = fomVals->GetBinContent(i+1, j+1);
220  selBinI = fomVals->GetXaxis()->GetBinLowEdge(i+1);
221  selBinJ = fomVals->GetYaxis()->GetBinLowEdge(j+1);
222  }
223  }
224  }
225 
226  std::cout
227  << "FOM is maximised by" << val
228  << " at ND: " << selBinI
229  << ", FD: " << selBinJ << std::endl;
230 
231  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
232  c1->SetRightMargin(0.18);
233  fomVals->Draw("colz");
234  c1->SaveAs("fomVals.png");
235 
236 }
int nBins
Definition: plotROC.py:16
set< int >::iterator it
TH1D * ConcatenateHists(TH1D *nd, TH1D *fd)
double calcBinnedSOverSSB(TH1D *projSig, std::vector< double > totUnc)
Definition: FOMUtilities.h:161
TMatrixD RemoveEmptyRowsAndColumns(TMatrixD *mat)
const double j
Definition: BetheBloch.cxx:29
void NDFDCVNCutSelector(std::string file)
OStream cout
Definition: OStream.cxx:6
TMatrixD FracToFull(TMatrixD *fracCov, TH1D *spectra)
Float_t proj
Definition: plot.C:35
TFile * file
Definition: cellShifts.C:17
TH1D * GetSelectedEvents(TH2D *h, std::string det, int cutBin)
c1
Definition: demo5.py:24
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void PrintHist(TH1D *h)
Eigen::MatrixXd mat
std::vector< double > GetDecorrelatedUncertainty(TMatrixD *inmat, int ndbins, int fdbins)
enum BeamMode string