2 #include "TDecompChol.h" 4 #include "TDecompSVD.h" 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);
57 h1->SetBinContent(
i, h2->GetBinContent(
i, bin));
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);
81 double cbes[32*12][
Nb*2];
83 for(
int be = 0; be <
Nb*2; ++be)
87 double mcfs[32][Nc*2];
89 for(
int cf = 0; cf < Nc*2; ++cf)
92 for(
int iw = 0; iw < h2->GetNbinsX()+2; ++iw){
95 for(
int ic = 0; ic < h2->GetNbinsY()+2; ++ic){
98 const int cell = ic-2;
99 const int mod = cell%32;
103 int isStart =
int(mod == 0);
104 int isEnd =
int(mod == 31);
105 int isMid =
int(mod == 15 || mod == 16);
107 const double vi = h2->GetBinContent(iw, ic);
111 for(
int a = 0;
a < Na*2; ++
a){
112 const double wa = wads[iw][
a];
113 for(
int b = 0;
b < Nb*2; ++
b){
114 const double cb = cbes[
cell][
b];
115 for(
int c = 0;
c < Nc*2; ++
c){
116 const double mc = mcfs[
mod][
c];
118 Avec[Nb*Nc*4*
a+Nc*2*
b+
c] += wa*cb*
mc;
120 if(
a < Na &&
b < Nb &&
c < Nc){
121 vec[Nb*Nc*
a+Nc*
b+
c] += vi*wa*cb*
mc;
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;
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;
136 vec[Na*Nb*Nc+0] += vi*isStart;
137 vec[Na*Nb*Nc+1] += vi*isEnd;
138 vec[Na*Nb*Nc+2] += vi*isMid;
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];
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;
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");
188 tr->Branch(
"view", &view);
189 int wpow, cellpow, modpow, startpow, midpow, endpow;
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);
199 TFile* f_correction =
new TFile(
"threshold_comparison.root",
"RECREATE");
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];
215 maxW = h->GetXaxis()->GetXmax();
217 if(!h2) h2 =
new TH2F(
"", viewStrUser+
";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
219 TProfile*
pr = h->ProfileX();
221 for(
int n = 0;
n < pr->GetNbinsX(); ++
n){
224 if(
fabs(h2->GetXaxis()->GetBinCenter(
n)) >
cutW && view!=2)
continue;
225 if(
fabs(h2->GetXaxis()->GetBinCenter(
n)) >
cutWMuonX && view==2)
continue;
228 h2->SetBinContent(
n,
cell+2, pr->GetBinContent(
n));
238 maxW = l->GetXaxis()->GetXmax();
240 if(!l2) l2 =
new TH2F(
"", viewStrUser+
";W;Cell", 100, -maxW, +maxW, cellMax+2, -1, cellMax+1);
242 TProfile* prl = l->ProfileX();
243 for(
int n = 0;
n < prl->GetNbinsX(); ++
n){
246 if(
fabs(l2->GetXaxis()->GetBinCenter(
n)) >
cutW && view!=2)
continue;
247 if(
fabs(l2->GetXaxis()->GetBinCenter(
n)) >
cutWMuonX && view==2)
continue;
250 l2->SetBinContent(
n,
cell+2, prl->GetBinContent(
n));
254 l2->GetXaxis()->CenterTitle();
255 l2->GetYaxis()->CenterTitle();
260 gPad->Print(
"plots/thresh_PEpercm_"+viewStr+
".eps");
264 h2->GetXaxis()->CenterTitle();
265 h2->GetYaxis()->CenterTitle();
270 gPad->Print(
"plots/thresh_corr_"+viewStr+
".eps");
273 TString newplotname =
"Cell31";
276 gPad->Print(
"plots/thresh_corr_cell31"+viewStr+
".eps");
285 for(
int a = 0;
a <
Na; ++
a){
286 for(
int b = 0;
b <
Nb; ++
b){
287 for(
int c = 0;
c <
Nc; ++
c){
291 coeff = p[Nb*Nc*
a+Nc*
b+
c];
300 coeff = p[Na*
Nb*
Nc+0];
304 coeff = p[Na*
Nb*
Nc+1];
308 coeff = p[Na*
Nb*
Nc+2];
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);
323 for(
int a = 0;
a <
Na; ++
a){
324 for(
int b = 0;
b <
Nb; ++
b){
325 for(
int c = 0;
c <
Nc; ++
c){
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);
338 h2fit->GetXaxis()->CenterTitle();
339 h2fit->GetYaxis()->CenterTitle();
341 h2fit->SetMinimum(h2->GetMinimum());
342 h2fit->SetMaximum(h2->GetMaximum());
343 gPad->Print(
"plots/thresh_corr_fit_"+viewStr+
".eps");
346 TH1* fx = h2fit->ProjectionX(
"a"+viewStr);
348 for(
int i = 0;
i < h2->GetNbinsY()+2; ++
i){
349 if(h2fit->GetBinContent(h2->GetNbinsX()/2+1,
i)) ++nbinsy;
351 fx->Scale(1./nbinsy);
352 fx->SetLineColor(
kRed);
353 fx->GetYaxis()->SetTitle(
"Correction factor");
354 fx->GetYaxis()->CenterTitle();
355 fx->SetTitle(viewStrUser);
358 TH1*
dx = h2->ProjectionX(
"b"+viewStr);
359 dx->Scale(1./nbinsy);
360 dx->SetLineColor(kBlack);
363 gPad->Print(
"plots/thresh_corr_proj_W_"+viewStr+
".eps");
364 fx->SetName(
"a"+viewStr);
368 TH1* fy = h2fit->ProjectionY(
"c"+viewStr);
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;
375 fy->Scale(1./nbinsx);
376 fy->SetLineColor(
kRed);
377 fy->GetYaxis()->SetTitle(
"Correction factor");
378 fy->GetYaxis()->CenterTitle();
379 fy->SetTitle(viewStrUser);
382 TH1*
dy = h2->ProjectionY(
"d"+viewStr);
383 dy->Scale(1./nbinsx);
384 dy->SetLineColor(kBlack);
387 gPad->Print(
"plots/thresh_corr_proj_cell_"+viewStr+
".eps");
388 fy->SetName(
"c"+viewStr);
392 TH2*
ratio = (TH2*)h2->Clone();
393 ratio->Divide(h2fit);
395 ratio->SetMinimum(.97);
396 ratio->SetMaximum(1.03);
397 gPad->Print(
"plots/thresh_corr_ratio_"+viewStr+
".eps");
400 TH1* px = ratio->ProjectionX(viewStr+
"_px");
401 px->Scale(1./nbinsy);
402 px->GetYaxis()->SetRangeUser(.98, 1.02);
404 gPad->Print(
"plots/thresh_corr_ratio_W_"+viewStr+
".eps");
407 TH1* py = ratio->ProjectionY(viewStr+
"_py");
408 py->Scale(1./nbinsx);
409 py->GetYaxis()->SetRangeUser(.98, 1.02);
411 gPad->Print(
"plots/thresh_corr_ratio_cell_"+viewStr+
".eps");
415 TFile*
fout =
new TFile(
"thresh_fits.root",
"RECREATE");
416 tr->SetDirectory(fout);
fvar< T > fabs(const fvar< T > &x)
TH1F * singlecellcorr(TH2 *h2, int bin, TString name)
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
double coeff(double W, int m, int l, bool real)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)