fit_thresh_corrs.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TDecompChol.h"
3 #include "TDecompLU.h"
4 #include "TDecompSVD.h"
5 #include "TF1.h"
6 #include "TFile.h"
7 #include "TGraph.h"
8 #include "TH2.h"
9 #include "TMatrixD.h"
10 #include "TVectorD.h"
11 #include "TProfile.h"
12 #include "TTree.h"
13 
14 // Order of fit in each dimension
15 const int Na = 12; // W
16 const int Nb = 6; // Cells
17 const int Nc = 4; // Position in module
18 const int Nflag = 3;
19 
20 // NearDet settings
21 //const int cellMax = 32*3;
22 //const double cutW = 190;
23 //const double cutWMuonX = 90;
24 //const double cellCorrMin = 0.5;
25 //const double cellCorrMax = 2.0;
26 //const double wCorrMin = 0.5;
27 //const double wCorrMax = 2.0;
28 
29 // TestBeam settings
30 const int cellMax = 32*2;
31 const double cutW = 114;
32 const double cutWMuonX = 90;//Needs to be set but is unused. Value is irrelevant.
33 const double cellCorrMin = 0.5;
34 const double cellCorrMax = 2.0;
35 const double wCorrMin = 0.5;
36 const double wCorrMax = 2.0;
37 
38 // FarDet settings, uncomment as appropriate
39 
40 //int cellMax = 32*12;
41 //double cutW = 750;
42 //const double cellCorrMin = 0.5;//2.0;
43 //const double cellCorrMax = 2.0;//2.6;
44 //const double wCorrMin = 0.5;//1.8;
45 //const double wCorrMax = 2.0;//3.0;
46 
47 TH1F* singlecellcorr(TH2* h2, int bin, TString name)
48 {
49 
50  int nbinsx = h2->GetXaxis()->GetNbins();
51  int minx = h2->GetXaxis()->GetXmin();
52  int maxx = h2->GetXaxis()->GetXmax();
53  TH1F* h1 = new TH1F (name, name+";W (cm);", nbinsx, minx, maxx);//h2->GetXaxis()->GetXmin(), h2->GetXaxis()->GetXmax);
54 
55  for(int i=0;i<nbinsx;i++)
56  {
57  h1->SetBinContent(i, h2->GetBinContent(i, bin));
58  }
59  return h1;
60 
61 }
62 
64 {
65  // Linearize three dimensions into one
67  // The A matrix has a symmetry property. Build up this vector instead
68  TVectorD Avec(Na*Nb*Nc*8);
69  TVectorD Nvec(Na*Nb*Nc);
70  TMatrixD Aflag(Na*Nb*Nc+Nflag, Nflag);
71  TVectorD vec(Na*Nb*Nc+Nflag);
72 
73  // Precompute values of pow(w/1000, a+d). 1000 is to scale to approximately
74  // 0-1 so that matrix entries will all be comparably sized.
75  double wads[h2->GetNbinsX()+2][Na*2];
76  for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw)
77  for(int ad = 0; ad < Na*2; ++ad)
78  wads[iw][ad] = pow(h2->GetXaxis()->GetBinCenter(iw)/1000, ad);
79 
80  // Precompute values of pow(cell/1000., b+e)
81  double cbes[32*12][Nb*2];
82  for(int cell = 0; cell < 32*12; ++cell)
83  for(int be = 0; be < Nb*2; ++be)
84  cbes[cell][be] = pow(cell/1000., be);
85 
86  // Precompute values of pow(mode, c+f)
87  double mcfs[32][Nc*2];
88  for(int mod = 0; mod < 32; ++mod)
89  for(int cf = 0; cf < Nc*2; ++cf)
90  mcfs[mod][cf] = pow(mod, cf);
91 
92  for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
93  // std::cout << iw << std::endl;
94  // const double w = h2->GetXaxis()->GetBinCenter(iw);
95  for(int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
96  // One for underflow, and one because we added one earlier to put a white
97  // space around
98  const int cell = ic-2;
99  const int mod = cell%32;
100 
101  // I don't understand why start and end behave differently
102  // physically. But empirically they do.
103  int isStart = int(mod == 0);
104  int isEnd = int(mod == 31);
105  int isMid = int(mod == 15 || mod == 16);
106 
107  const double vi = h2->GetBinContent(iw, ic);
108  // Don't fit to zeros
109  if(!vi) continue;
110 
111  for(int a = 0; a < Na*2; ++a){
112  const double wa = wads[iw][a]; //pow(w/1000, a);
113  for(int b = 0; b < Nb*2; ++b){
114  const double cb = cbes[cell][b]; //pow(cell/1000., b);
115  for(int c = 0; c < Nc*2; ++c){
116  const double mc = mcfs[mod][c]; //pow(mod/10., c);
117 
118  Avec[Nb*Nc*4*a+Nc*2*b+c] += wa*cb*mc;
119 
120  if(a < Na && b < Nb && c < Nc){
121  vec[Nb*Nc*a+Nc*b+c] += vi*wa*cb*mc;
122 
123  if(isStart) Aflag[Nb*Nc*a+Nc*b+c][0] += wa*cb*mc;
124  if(isEnd) Aflag[Nb*Nc*a+Nc*b+c][1] += wa*cb*mc;
125  if(isMid) Aflag[Nb*Nc*a+Nc*b+c][2] += wa*cb*mc;
126  }
127  } // end for c
128  } // end for b
129  } // end for a
130 
131  // Only need to do diagonals. Can't have both flags at once
132  if(isStart) Aflag[Na*Nb*Nc+0][0] += 1;
133  if(isEnd) Aflag[Na*Nb*Nc+1][1] += 1;
134  if(isMid) Aflag[Na*Nb*Nc+2][2] += 1;
135 
136  vec[Na*Nb*Nc+0] += vi*isStart;
137  vec[Na*Nb*Nc+1] += vi*isEnd;
138  vec[Na*Nb*Nc+2] += vi*isMid;
139  } // end for ic
140  } // end for iw
141 
142  // Expand Avec out into the cells of A
143  for(int a = 0; a < Na; ++a){
144  for(int b = 0; b < Nb; ++b){
145  for(int c = 0; c < Nc; ++c){
146  for(int d = 0; d < Na; ++d){
147  for(int e = 0; e < Nb; ++e){
148  for(int f = 0; f < Nc; ++f){
149  A[Nb*Nc*d+Nc*e+f][Nb*Nc*a+Nc*b+c] = Avec[Nb*Nc*4*(a+d)+Nc*2*(b+e)+c+f];
150  }
151  }
152  }
153  }
154  }
155  }
156 
157  for(int flag = 0; flag < Nflag; ++flag){
158  for(int abc = 0; abc < Na*Nb*Nc+Nflag; ++abc){
159  const double val = Aflag[abc][flag];
160  A[Na*Nb*Nc+flag][abc] = val;
161  A[abc][Na*Nb*Nc+flag] = val;
162  }
163  }
164 
165  // A.Print();
166  // vec.Print();
167 
168  TDecompSVD dc(A);
169  bool junk;
170  TVectorD p = dc.Solve(vec, junk);
171 
172  return p;
173 }
174 
176 {
177  //std::cout << "SHOULD HAVE CHANGES" << std::endl;
178  // Output from ThresholdAna_module
179 
180  // TFile* fin = new TFile("fdThresh_1.root");
181  // TFile* fin = new TFile("/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/threshold_2ndAna_fulldetector_ndmc_idealconditions.root"); // threshold and shadowing file. Output file of ThresholdAna_module.cc
182  // TFile* fin = new TFile("/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/consolidated/threshana_hadded.root");
183  // TFile* fin = new TFile("/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/threshold_2ndAna_fulldetector_ndmc_idealconditions.root");
184  TFile* fin = new TFile("/pnfs/nova/persistent/users/kmulder/testbeam_rel_calib_28062019/thresholdCorr_mc/thresholdAna_mc_hadded.root");
185  TTree* tr = new TTree("tree", "tree");
186  int nViews = 2;//CHANGE TO 2 FOR FD, 4 FOR ND
187  int view;
188  tr->Branch("view", &view);
189  int wpow, cellpow, modpow, startpow, midpow, endpow;
190  double coeff;
191  tr->Branch("wpow", &wpow);
192  tr->Branch("cellpow", &cellpow);
193  tr->Branch("modpow", &modpow);
194  tr->Branch("startpow", &startpow);
195  tr->Branch("midpow", &midpow);
196  tr->Branch("endpow", &endpow);
197  tr->Branch("coeff", &coeff);
198 
199  TFile* f_correction = new TFile("threshold_comparison.root", "RECREATE"); // this file will be used to compare threshold and shielding correction factors.
200  //TFile* fout = new TFile("thresh_fits.root", "RECREATE");
201 
202  double maxW = -1;
203  const TString viewStrArr[4]={"xview","yview","muxview","muyview"};
204  const TString viewStrUserArr[4]={"X View","Y View","Muon Catcher X View","Muon Catcher Y View"};
205  for(view = 0; view < nViews; ++view){
206  const TString viewStr = viewStrArr[view];
207  const TString viewStrUser = viewStrUserArr[view];
208 
209  TH2* h2 = 0; // Wait until maxW is set
210 
211  for(int cell = 0; cell < cellMax; ++cell){
212  TH2* h = (TH2*)fin->Get(TString::Format("thresh/ratioVsW_%s_cell%d", viewStr.Data(), cell));
213  if(!h) continue;
214 
215  maxW = h->GetXaxis()->GetXmax();
216 
217  if(!h2) h2 = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
218 
219  TProfile* pr = h->ProfileX();
220 
221  for(int n = 0; n < pr->GetNbinsX(); ++n){
222  // Weird stuff happens right at the edges. Don't include it in the fit,
223  // because it compromises the main body.
224  if(fabs(h2->GetXaxis()->GetBinCenter(n)) > cutW && view!=2) continue;
225  if(fabs(h2->GetXaxis()->GetBinCenter(n)) > cutWMuonX && view==2) continue;
226 
227  // +2, one for overflow, one to put a white cell padding around
228  h2->SetBinContent(n, cell+2, pr->GetBinContent(n));
229  }
230  }
231 
232 
233  TH2* l2 = 0;
234  for(int cell = 0; cell < cellMax; ++cell){
235  TH2* l = (TH2*)fin->Get(TString::Format("thresh/fPEpercmVsW_%s_cell%d", viewStr.Data(), cell));
236  if(!l) continue;
237 
238  maxW = l->GetXaxis()->GetXmax();
239 
240  if(!l2) l2 = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
241 
242  TProfile* prl = l->ProfileX();
243  for(int n = 0; n < prl->GetNbinsX(); ++n){
244  // Weird stuff happens right at the edges. Don't include it in the fit,
245  // because it compromises the main body.
246  if(fabs(l2->GetXaxis()->GetBinCenter(n)) > cutW && view!=2) continue;
247  if(fabs(l2->GetXaxis()->GetBinCenter(n)) > cutWMuonX && view==2) continue;
248 
249  // +2, one for overflow, one to put a white cell padding around
250  l2->SetBinContent(n, cell+2, prl->GetBinContent(n));
251  }
252  }
253  new TCanvas;
254  l2->GetXaxis()->CenterTitle();
255  l2->GetYaxis()->CenterTitle();
256  l2->Draw("colz");
257  // Don't need to see zeros
258  //l2->SetMinimum(.5);
259  //l2->SetMaximum(3.0);
260  gPad->Print("plots/thresh_PEpercm_"+viewStr+".eps");
261 
262 
263  new TCanvas;
264  h2->GetXaxis()->CenterTitle();
265  h2->GetYaxis()->CenterTitle();
266  h2->Draw("colz");
267  // Don't need to see zeros
268  h2->SetMinimum(.5);
269  h2->SetMaximum(3.0);
270  gPad->Print("plots/thresh_corr_"+viewStr+".eps");
271 
272  new TCanvas;
273  TString newplotname = "Cell31";
274  TH1F* cell31plot = singlecellcorr(h2, 31, newplotname);
275  cell31plot->Draw();
276  gPad->Print("plots/thresh_corr_cell31"+viewStr+".eps");
277 
278  const TVectorD p = solve(h2);
279 
280 
281  // Save coefficients to tree
282  startpow = 0;
283  midpow = 0;
284  endpow = 0;
285  for(int a = 0; a < Na; ++a){
286  for(int b = 0; b < Nb; ++b){
287  for(int c = 0; c < Nc; ++c){
288  wpow = a;
289  cellpow = b;
290  modpow = c;
291  coeff = p[Nb*Nc*a+Nc*b+c];
292  tr->Fill();
293  }
294  }
295  }
296  wpow = 0;
297  cellpow = 0;
298  modpow = 0;
299  startpow = 1;
300  coeff = p[Na*Nb*Nc+0];
301  tr->Fill();
302  startpow = 0;
303  endpow = 1;
304  coeff = p[Na*Nb*Nc+1];
305  tr->Fill();
306  endpow = 0;
307  midpow = 1;
308  coeff = p[Na*Nb*Nc+2];
309  tr->Fill();
310 
311  // Make a map of what the fit looks like
312  TH2* h2fit = new TH2F("", viewStrUser+";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
313  for(int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
314  const double w = h2->GetXaxis()->GetBinCenter(iw);
315  for(int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
316  const int cell = ic-2;
317  const int mod = cell%32;
318  int isStart = int(mod == 0);
319  int isEnd = int(mod == 31);
320  int isMid = int(mod == 15 || mod == 16);
321 
322  double v = 0;
323  for(int a = 0; a < Na; ++a){
324  for(int b = 0; b < Nb; ++b){
325  for(int c = 0; c < Nc; ++c){
326  v += p[Nb*Nc*a+Nc*b+c]*pow(w/1000, a)*pow(cell/1000., b)*pow(mod, c);
327  }
328  }
329  }
330  if(isStart) v += p[Na*Nb*Nc+0];
331  if(isEnd) v += p[Na*Nb*Nc+1];
332  if(isMid) v += p[Na*Nb*Nc+2];
333  if(h2->GetBinContent(iw, ic)) h2fit->SetBinContent(iw, ic, v);
334  }
335  }
336 
337  new TCanvas;
338  h2fit->GetXaxis()->CenterTitle();
339  h2fit->GetYaxis()->CenterTitle();
340  h2fit->Draw("colz");
341  h2fit->SetMinimum(h2->GetMinimum());
342  h2fit->SetMaximum(h2->GetMaximum());
343  gPad->Print("plots/thresh_corr_fit_"+viewStr+".eps");
344 
345  new TCanvas;
346  TH1* fx = h2fit->ProjectionX("a"+viewStr);
347  int nbinsy = 0;
348  for(int i = 0; i < h2->GetNbinsY()+2; ++i){
349  if(h2fit->GetBinContent(h2->GetNbinsX()/2+1, i)) ++nbinsy;
350  }
351  fx->Scale(1./nbinsy);
352  fx->SetLineColor(kRed);
353  fx->GetYaxis()->SetTitle("Correction factor");
354  fx->GetYaxis()->CenterTitle();
355  fx->SetTitle(viewStrUser);
356  fx->GetYaxis()->SetRangeUser(wCorrMin, wCorrMax);//1.1, 1.3
357  fx->Draw("l");
358  TH1* dx = h2->ProjectionX("b"+viewStr);
359  dx->Scale(1./nbinsy);
360  dx->SetLineColor(kBlack);
361  dx->Draw("same");
362 
363  gPad->Print("plots/thresh_corr_proj_W_"+viewStr+".eps");
364  fx->SetName("a"+viewStr);
365  fx->Write();
366 
367  new TCanvas;
368  TH1* fy = h2fit->ProjectionY("c"+viewStr);
369  int nbinsx = 0;
370  for(int i = 0; i < h2->GetNbinsX()+2; ++i){
371  std::cout << "xbin : " << i << ", ybins/2 +1 : " << h2->GetNbinsY()/2+1 << ", content: " << h2fit->GetBinContent(i, h2->GetNbinsY()/2+1) << std::endl;
372  if(h2fit->GetBinContent(i, h2->GetNbinsY()/2+1)) ++nbinsx;
373  }
374  std::cout << nbinsx << " nbinsx" << std::endl;
375  fy->Scale(1./nbinsx);
376  fy->SetLineColor(kRed);
377  fy->GetYaxis()->SetTitle("Correction factor");
378  fy->GetYaxis()->CenterTitle();
379  fy->SetTitle(viewStrUser);
380  fy->GetYaxis()->SetRangeUser(cellCorrMin, cellCorrMax);//1.1, 1.25
381  fy->Draw("l");
382  TH1* dy = h2->ProjectionY("d"+viewStr);
383  dy->Scale(1./nbinsx);
384  dy->SetLineColor(kBlack);
385  dy->Draw("same");
386 
387  gPad->Print("plots/thresh_corr_proj_cell_"+viewStr+".eps");
388  fy->SetName("c"+viewStr);
389  fy->Write();
390 
391  new TCanvas;
392  TH2* ratio = (TH2*)h2->Clone();
393  ratio->Divide(h2fit);
394  ratio->Draw("colz");
395  ratio->SetMinimum(.97);
396  ratio->SetMaximum(1.03);
397  gPad->Print("plots/thresh_corr_ratio_"+viewStr+".eps");
398 
399  new TCanvas;
400  TH1* px = ratio->ProjectionX(viewStr+"_px");
401  px->Scale(1./nbinsy);
402  px->GetYaxis()->SetRangeUser(.98, 1.02);
403  px->Draw();
404  gPad->Print("plots/thresh_corr_ratio_W_"+viewStr+".eps");
405 
406  new TCanvas;
407  TH1* py = ratio->ProjectionY(viewStr+"_py");
408  py->Scale(1./nbinsx);
409  py->GetYaxis()->SetRangeUser(.98, 1.02);
410  py->Draw();
411  gPad->Print("plots/thresh_corr_ratio_cell_"+viewStr+".eps");
412  } // end for view
413 // f_correction->Close();
414 
415  TFile* fout = new TFile("thresh_fits.root", "RECREATE");
416  tr->SetDirectory(fout);
417  tr->Write();
418  fout->Close();
419 
420  std::cout << "Wrote results to thresh_fits.root" << std::endl;
421 }
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const int Na
TH1F * singlecellcorr(TH2 *h2, int bin, TString name)
const char * p
Definition: xmltok.h:285
TH1 * ratio(TH1 *h1, TH1 *h2)
constexpr T pow(T x)
Definition: pow.h:75
const double cutW
const int Nc
const double wCorrMin
double dy[NP][NC]
double dx[NP][NC]
const double a
Float_t d
Definition: plot.C:236
Eigen::VectorXd vec
TVectorD solve(TH2 *h2)
float bin[41]
Definition: plottest35.C:14
double coeff(double W, int m, int l, bool real)
const double cellCorrMin
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
static const double A
Definition: Units.h:82
const double cutWMuonX
Int_t nbinsx
Definition: plot.C:23
const int cellMax
const double wCorrMax
const hit & b
Definition: hits.cxx:21
void fit_thresh_corrs()
const double cellCorrMax
const int Nb
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
const int Nflag