fit_thresh_corrs_fb.C
Go to the documentation of this file.
1 
2 #include "TCanvas.h"
3 #include "TDecompChol.h"
4 #include "TDecompLU.h"
5 #include "TDecompSVD.h"
6 #include "TF1.h"
7 #include "TFile.h"
8 #include "TGraph.h"
9 #include "TH2.h"
10 #include "TMatrixD.h"
11 #include "TVectorD.h"
12 #include "TProfile.h"
13 #include "TTree.h"
14 #include "TStyle.h"
15 
16 ///Commentary
17 //==============
18 // Ryan Nichol changed this code on 20th March 2019 to remove the loop over position in the moudle and replace the start, mid and end flags
19 // with a flag for each of the 32 positions in the module.
20 // Also at this time we changes the arbitrary 1000 scale factor to dividing by cutW and cellMax (respectively)
21 // Finally we changed the order of the polynomial in W from 12->15
22 // We also chnaged the output format to store some more useful numbers
23 
24 // Order of fit in each dimension
25 const int Na = 15; // W
26 const int Nb = 6; // Cells
27 const int Nflag = 32;
28 
29 //CHANGE
30 // NearDet settings, uncomment as appropriate
31 // const int cellMax = 32*3;
32 // const double cutW = 190;
33 const double cutWMuonX = 80;//DON'T uncomment this line
34 // const double cellCorrMin = 0.5;
35 // const double cellCorrMax = 2.0;
36 // const double wCorrMin = 0.5;
37 // const double wCorrMax = 2.0;
38 
39 // FarDet settings, uncomment as appropriate
40 int cellMax = 32*12;
41 double cutW = 750;
42 const double cellCorrMin = 2.0;
43 const double cellCorrMax = 2.6;
44 const double wCorrMin = 1.8;
45 const double wCorrMax = 3.0;
46 
47 //CHANGE true if using fibre brightness
48 const bool fFiberBrightnessMode = true;
49 const bool fMuCFiberBrightnessMode = false;
50 
52 {
53  // Linearize three dimensions into one
55 
56  // The A matrix has a symmetry property. Build up this vector instead
57  TVectorD Avec(Na*Nb*8); //Why is this 8 here?
58  TVectorD Nvec(Na*Nb);
59  TMatrixD Aflag(Na*Nb+Nflag, Nflag);
60  TVectorD vec(Na*Nb+Nflag);
61 
62  // Precompute values of pow(w/cutW, a+d). 1000 is to scale to approximately
63  // 0-1 so that matrix entries will all be comparably sized.
64  double wads[h2->GetNbinsX()+2][Na*2];
65  for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw)
66  for(int ad = 0; ad < Na*2; ++ad)
67  wads[iw][ad] = pow(h2->GetXaxis()->GetBinCenter(iw)/cutW, ad);
68 
69  // Precompute values of pow(cell/cellMax, b+e)
70  double cbes[cellMax][Nb*2];
71  for(int cell = 0; cell < cellMax; ++cell)
72  for(int be = 0; be < Nb*2; ++be)
73  cbes[cell][be] = pow(cell/double(cellMax), be);
74 
75  for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
76  // std::cout << iw << std::endl;
77  // const double w = h2->GetXaxis()->GetBinCenter(iw);
78  for(int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
79  // One for underflow, and one because we added one earlier to put a white
80  // space around
81  // Argh Chris's choice of silly binning is very annoying
82  const int cell = ic-2;
83  const int mod = cell%32;
84 
85  const double vi = h2->GetBinContent(iw, ic);
86  // Don't fit to zeros
87  if(!vi) continue;
88 
89  for(int a = 0; a < Na*2; ++a){
90  const double wa = wads[iw][a]; //pow(w/cutW, a);
91  for(int b = 0; b < Nb*2; ++b){
92  const double cb = cbes[cell][b]; //pow(cell/cellMax, b);
93 
94  Avec[Nb*2*a+b] += wa*cb;
95 
96  if(a < Na && b < Nb ){
97  vec[Nb*a+b] += vi*wa*cb;
98  Aflag[Nb*a+b][mod] += wa*cb;
99  }
100  } // end for b
101  } // end for a
102 
103  // Now need to deal with the diagonals for the flags for cells within a module dependence
104  Aflag[Na*Nb+mod][mod] += 1;
105  vec[Na*Nb+mod] += vi;
106  } // end for ic
107  } // end for iw
108 
109  // Expand Avec out into the cells of A
110  for(int a = 0; a < Na; ++a){
111  for(int b = 0; b < Nb; ++b){
112  for(int d = 0; d < Na; ++d){
113  for(int e = 0; e < Nb; ++e){
114  A[Nb*d+e][Nb*a+b] = Avec[Nb*2*(a+d)+(b+e)];
115  }
116  }
117  }
118  }
119 
120  for(int flag = 0; flag < Nflag; ++flag){
121  for(int abc = 0; abc < Na*Nb+Nflag; ++abc){
122  const double val = Aflag[abc][flag];
123  A[Na*Nb+flag][abc] = val;
124  A[abc][Na*Nb+flag] = val;
125  }
126  }
127 
128  // A.Print();
129  // vec.Print();
130 
131  TDecompSVD dc(A);
132  bool junk;
133  TVectorD p = dc.Solve(vec, junk);
134 
135  return p;
136 }
137 
138 void fit_thresh_corrs_fb(std::string thresh_input_file, std::string thresh_fit_file, int nBrightnesses=12, int minfb=0)
139 {
140  gStyle->SetOptStat(0);
141 
142  // Output from ThresholdAna_module
143  // CHANGE
144  TFile* fin = new TFile(thresh_input_file.c_str());
145  TFile* fout = new TFile(thresh_fit_file.c_str(), "recreate");
146  TTree* tr = new TTree("tree", "tree");
147 
148  int nViews = 2;//4 //CHANGE FD: 2, ND with FB: 2, ND without FB (i.e. for muon catcher): 4
149  // int nBrightnesses = 4;//CHANGE
150  // int minfb = 0; //CHANGE :I recommend doing this in two batches: fb0-4, fb5-8, i.e. nBrightnesses = 5 or 4, minfb = 0 or 5
151  //if(fFiberBrightnessMode) nBrightnesses = 9; //I recommend doing this in two batches: fb0-4, fb5-8
152 
153  int view;
154  int fiberbrightness;
155  int wpow, cellpow, modpow;//, startpow, midpow, endpow;
156  double coeff;
157  double dcellMax,dcutW, dcutWMax;
158  int NaTree, NbTree, NflagTree;
159 
160  tr->Branch("fiberbrightness", &fiberbrightness);
161  tr->Branch("view", &view);
162  // tr->Branch("Na",const_cast<int&>(Na));
163  // tr->Branch("Nb",const_cast<int&>(Nb));
164  // tr->Branch("Nflag",const_cast<int&>(Nflag));
165  tr->Branch("Na",&NaTree);
166  tr->Branch("Nb",&NbTree);
167  tr->Branch("Nflag",&NflagTree);
168  tr->Branch("wpow", &wpow);
169  tr->Branch("cellpow", &cellpow);
170  tr->Branch("modpow", &modpow);
171  tr->Branch("dcutW",&dcutW);
172  tr->Branch("dcutWMax",&dcutWMax);
173  tr->Branch("dcellMax",&dcellMax);
174  tr->Branch("coeff", &coeff);
175 
176  NaTree = Na;
177  NbTree = Nb;
178  NflagTree = Nflag;
179  dcellMax=cellMax;
180  dcutWMax = cutW;
181 
182  // TFile* f_correction = new TFile("threshold_correction.root", "RECREATE"); // this file will be used to compare threshold and shielding correction factors.
183  // TFile* fout = new TFile("thresh_fits_thirdanalysis_realandupdatedlightlevels_nd_fb5-8.root", "RECREATE"); //CHANGE
184 
185  double maxW = -1;
186  const TString viewStrArr[4]={"xview","yview","muxview","muyview"};
187  const TString viewStrUserArr[4]={"X View","Y View","Muon Catcher X View","Muon Catcher Y View"};
188  for(fiberbrightness = minfb; fiberbrightness < nBrightnesses+minfb; ++fiberbrightness){
189 
190  for(view = 0; view < nViews; ++view){
191 
192  if(view > 1 && !fMuCFiberBrightnessMode && fiberbrightness > 0) continue;
193  if(view == 2) dcutW = cutWMuonX;
194  else dcutW = cutW;
195 
196  const TString viewStr = viewStrArr[view];
197  const TString viewStrUser = viewStrUserArr[view];
198 
199  TString fbViewStr = viewStr;
200  TString fbViewStrUser = viewStrUser;
201 
202  TH2* h2 = 0; // Wait until maxW is set
203 
204  if((fFiberBrightnessMode && view <= 1) || (fMuCFiberBrightnessMode && view > 1)){
205  fbViewStr = TString::Format("%s_fb%i",viewStr.Data(), fiberbrightness);
206  fbViewStrUser = TString::Format("%s, Fiber Brightness:%i",viewStrUser.Data(), fiberbrightness);
207  }
208 
209  for(int cell = 0; cell < cellMax; ++cell){
210  //std::cout << "Looking for " << TString::Format("thresh/ratioVsW_%s_cell%d", fbViewStr.Data(), cell) << std::endl;
211 
212  TH2* h = (TH2*)fin->Get(TString::Format("thresh/ratioVsW_%s_cell%d", fbViewStr.Data(), cell));
213 
214  if(!h) {
215  std::cout << "Can't find thresh/ratioVsW_" << fbViewStr.Data() << "_cell" << cell << std::endl;
216  continue;
217  }
218 
219  maxW = h->GetXaxis()->GetXmax();
220 
221  if(!h2){
222  h2 = new TH2F("", fbViewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
223  std::cout << "Created h2 with title " << fbViewStrUser.Data() << std::endl;
224  }
225 
226  TProfile* pr = h->ProfileX();
227 
228  for(int n = 0; n < pr->GetNbinsX(); ++n){
229  // Weird stuff happens right at the edges. Don't include it in the fit,
230  // because it compromises the main body.
231  if(fabs(h2->GetXaxis()->GetBinCenter(n)) > cutW) continue;
232  if(h2->GetXaxis()->GetBinCenter(n) > cutWMuonX && view==2) continue;
233 
234  // +2, one for overflow, one to put a white cell padding around
235  h2->SetBinContent(n, cell+2, pr->GetBinContent(n));
236  }
237  }
238 
239  if(!h2) std::cout << "Can't find h2!" << std::endl;
240  else std::cout << "Drawing h2..." << std::endl;
241  new TCanvas;
242  h2->GetXaxis()->CenterTitle();
243  h2->GetYaxis()->CenterTitle();
244  h2->Draw("colz");
245  // Don't need to see zeros
246  h2->SetMinimum(.5);
247  //gPad->Print("plots/thresh_corr_"+fbViewStr+".png");
248  gPad->Print("plots/thresh_corr_"+fbViewStr+".png");
249 
250  std::cout << "Solving h2..." << std::endl;
251  const TVectorD p = solve(h2);
252 
253 
254  // Save coefficients to tree
255  std::cout << "Saving coefficients..." << std::endl;
256  for(int a = 0; a < Na; ++a){
257  for(int b = 0; b < Nb; ++b){
258  wpow = a;
259  cellpow = b;
260  coeff = p[Nb*a+b];
261  tr->Fill();
262  }
263  }
264  wpow = 0;
265  cellpow = 0;
266  for(modpow=0;modpow<32;modpow++) {
267  coeff = p[Na*Nb+modpow];
268  tr->Fill();
269  }
270 
271  // Make a map of what the fit looks like
272  std::cout << "Saving a map..." << std::endl;
273  TH2* h2fit = new TH2F("", fbViewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
274  for(int iw = 1; iw < h2->GetNbinsX()+1; ++iw){
275  const double w = h2->GetXaxis()->GetBinCenter(iw);
276  for(int ic = 1; ic < h2->GetNbinsY()+1; ++ic){
277  const int cell = ic-2;
278  const int mod = cell%32;
279 
280  double v = 0;
281  for(int a = 0; a < Na; ++a){
282  for(int b = 0; b < Nb; ++b){
283  v += p[Nb*a+b]*pow(w/dcutWMax, a)*pow(cell/dcellMax, b);
284  }
285  }
286  //Now add in the poisiton within module addition
287  v += p[Na*Nb+mod];
288  if(h2->GetBinContent(iw, ic)) h2fit->SetBinContent(iw, ic, v);
289  }
290  }
291 
292  std::cout << "Drawing the map..." << std::endl;
293  new TCanvas;
294  h2fit->GetXaxis()->CenterTitle();
295  h2fit->GetYaxis()->CenterTitle();
296  h2fit->Draw("colz");
297  h2fit->SetMinimum(h2->GetMinimum());
298  h2fit->SetMaximum(h2->GetMaximum());
299  //h2fit->SetMinimum(0.5);
300  //h2fit->SetMaximum(1.4);
301  //gPad->Print("plots/thresh_corr_fit_"+fbViewStr+".png");
302  gPad->Print("plots/thresh_corr_fit_"+fbViewStr+".png");
303 
304  std::cout << "And the map's x projection..." << std::endl;
305  new TCanvas;
306  TH1* fx = h2fit->ProjectionX("a"+fbViewStr);
307  int nbinsy = 0;
308  for(int i = 0; i < h2->GetNbinsY()+2; ++i){
309  if(h2fit->GetBinContent(h2->GetNbinsX()/2+1, i)) ++nbinsy;
310  }
311  fx->Scale(1./nbinsy);
312  fx->SetLineColor(kRed);
313  fx->GetYaxis()->SetTitle("Correction factor");
314  fx->GetYaxis()->CenterTitle();
315  fx->SetMinimum(1.0);//1.1);
316  fx->SetMaximum(2.0);//was 1.3
317  fx->SetTitle(fbViewStrUser);
318  //fx->GetYaxis()->SetRangeUser(wCorrMin, wCorrMax);
319 
320  fx->Draw("l");
321  TH1* dx = h2->ProjectionX("b"+fbViewStr);
322  dx->Scale(1./nbinsy);
323  dx->SetLineColor(kBlack);
324  dx->Draw("same");
325 
326  //gPad->Print("plots/thresh_corr_proj_W_"+fbViewStr+".png");
327  gPad->Print("plots/thresh_corr_proj_W_"+fbViewStr+".png");
328  fx->SetName("a"+fbViewStr);
329  fx->Write();
330 
331  std::cout << "And the map's y projection..." << std::endl;
332  new TCanvas;
333  TH1* fy = h2fit->ProjectionY("c"+fbViewStr);
334  int nbinsx = 0;
335  for(int i = 0; i < h2->GetNbinsX()+2; ++i){
336  if(h2fit->GetBinContent(i, h2->GetNbinsY()/2+1)) ++nbinsx;
337  }
338  // fy->Sumw2();
339  fy->Scale(1./nbinsx);
340  fy->SetLineColor(kRed);
341  fy->GetYaxis()->SetTitle("Correction factor");
342  fy->GetYaxis()->CenterTitle();
343  fy->SetMinimum(1.0);//1.1
344  fy->SetMaximum(2.0);//was 1.6
345  fy->SetTitle(fbViewStrUser);
346  //fy->GetYaxis()->SetRangeUser(cellCorrMin, cellCorrMax);
347  fy->Draw("l");
348  TH1* dy = h2->ProjectionY("d"+fbViewStr);
349  dy->Scale(1./nbinsx);
350  dy->SetLineColor(kBlack);
351  dy->Draw("same");
352 
353  //gPad->Print("plots/thresh_corr_proj_cell_"+fbViewStr+".png");
354  gPad->Print("plots/thresh_corr_proj_cell_"+fbViewStr+".png");
355  fy->SetName("c"+fbViewStr);
356  fy->Write();
357 
358  new TCanvas;
359  TH2* ratio = (TH2*)h2->Clone();
360  ratio->Divide(h2fit);
361  ratio->Draw("colz");
362  ratio->SetMinimum(.97);
363  ratio->SetMaximum(1.03);
364  //gPad->Print("plots/thresh_corr_ratio_"+fbViewStr+".png");
365  gPad->Print("plots/thresh_corr_ratio_"+fbViewStr+".png");
366 
367  new TCanvas;
368  TH1* px = ratio->ProjectionX(fbViewStr+"_px");
369  px->Scale(1./nbinsy);
370  px->GetYaxis()->SetRangeUser(.98, 1.02);
371  px->Draw();
372  //gPad->Print("plots/thresh_corr_ratio_W_"+fbViewStr+".png");
373  gPad->Print("plots/thresh_corr_ratio_W_"+fbViewStr+".png");
374 
375  new TCanvas;
376  TH1* py = ratio->ProjectionY(fbViewStr+"_py");
377  py->Scale(1./nbinsx);
378  py->GetYaxis()->SetRangeUser(.98, 1.02);
379  py->Draw();
380  //gPad->Print("plots/thresh_corr_ratio_cell_"+fbViewStr+".png");
381  gPad->Print("plots/thresh_corr_ratio_cell_"+fbViewStr+".png");
382  }// end for view
383 
384  } // end for fibre brightness
385  // f_correction->Close();
386 
387  //TFile* fout = new TFile("thresh_fits_fb5-8.root", "RECREATE");
388  tr->Write();
389  fout->Close();
390  //std::cout << "Wrote results to thresh_fits_fd_gain100_fb0-4.root" << std::endl;
391 
392 }
TString fin
Definition: Style.C:24
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int cellMax
const char * p
Definition: xmltok.h:285
TH1 * ratio(TH1 *h1, TH1 *h2)
constexpr T pow(T x)
Definition: pow.h:75
double dy[NP][NC]
double dx[NP][NC]
const double wCorrMax
void fit_thresh_corrs_fb(std::string thresh_input_file, std::string thresh_fit_file, int nBrightnesses=12, int minfb=0)
const double a
const bool fMuCFiberBrightnessMode
Float_t d
Definition: plot.C:236
Eigen::VectorXd vec
double coeff(double W, int m, int l, bool real)
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
const bool fFiberBrightnessMode
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const double cutWMuonX
static const double A
Definition: Units.h:82
double cutW
const double wCorrMin
Int_t nbinsx
Definition: plot.C:23
const hit & b
Definition: hits.cxx:21
const double cellCorrMax
const int Nb
const double cellCorrMin
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
TVectorD solve(TH2 *h2)
const int Nflag
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
const int Na
Commentary.