Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
figmaker.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file figmaker.C
1 
2 #include "sPhenixStyle.C"
3 
4 void myText(Double_t x,Double_t y,Color_t color, const char *text, Double_t tsize=0.04) {
5 
6  TLatex l; l.SetTextAlign(12); l.SetTextSize(tsize);
7  l.SetNDC();
8  l.SetTextColor(color);
9  l.DrawLatex(x,y,text);
10 }
11 
12 TH1F* proj(TH2F* h2);
13 
14 TH1F* FBratio(TH1F* h){
15 
16  const int Nbin = h->GetXaxis()->GetNbins();
17  TH1F* hfb = new TH1F("temp32","",Nbin/2,0,Nbin/2);
18 
19  for(int i=0; i<Nbin/2; i++){
20  int b1 = i+1;
21  int b2 = Nbin-i;
22  float ratio = h->GetBinContent(b1)/h->GetBinContent(b2);
23  float err = sqrt( pow(h->GetBinError(b1)/h->GetBinContent(b1),2)+ pow(h->GetBinError(b2)/h->GetBinContent(b2),2))*pow(ratio,2);
24  hfb->SetBinContent(i,ratio);
25  hfb->SetBinError(i,err);
26  }
27  return hfb;
28 }
29 
30 
31 
32 void figmaker(){
33 
34 
36 
37  int totalEvents = 0;
38 
39  std::ifstream infile("../condor/runList.txt");
40  if (!infile) { std::cerr << "Error: Unable to open the file\n"; return 1;}
41  std::vector<int> runList;
42  int num;
43  while (infile >> num) runList.push_back(num);
44 
45  int Nruns= runList.size();
46 
47  for (int ir=0; ir<2; ir++){
48 
49  int run = runList[ir];
50 
51  TFile* fin = new TFile(Form("../condor/combine_out/out_%d.root",run));
52 
53  TH2F* h_emcal_mbd_correlation = (TH2F*) fin->Get("h_emcal_mbd_correlation" );
54  TH2F* h_ihcal_mbd_correlation = (TH2F*) fin->Get("h_ihcal_mbd_correlation" );
55  TH2F* h_ohcal_mbd_correlation = (TH2F*) fin->Get("h_ohcal_mbd_correlation" );
56  TH2F* h_emcal_hcal_correlation = (TH2F*) fin->Get("h_emcal_hcal_correlation");
57  TH1F* h_InvMass = (TH1F*) fin->Get("h_InvMass" );
58  TH2F* h_cemc_etaphi = (TH2F*) fin->Get("h_cemc_etaphi_wQA" );
59  TH2F* h_ihcal_etaphi = (TH2F*) fin->Get("h_ihcal_etaphi_wQA" );
60  TH2F* h_ohcal_etaphi = (TH2F*) fin->Get("h_ohcal_etaphi_wQA" );
61  TH2F* h_zdc_emcal_correlation = (TH2F*) fin->Get("h_zdc_emcal_correlation");
62  TH1F* hzdcNorthcalib = (TH1F*) fin->Get("hzdcNorthcalib");
63  TH1F* hzdcSouthcalib = (TH1F*) fin->Get("hzdcSouthcalib");
64  TH2F* h_etaphi_clus = (TH2F*) fin->Get("h_etaphi_clus" );
65  TH1F* hvtx_z_raw = (TH1F*) fin->Get("hvtx_z_raw");
66  TProfile2D* h_cemc_etaphi_time = (TProfile2D*) fin->Get("h_cemc_etaphi_time");
67  TProfile2D* h_ihcal_etaphi_time = (TProfile2D*) fin->Get("h_ihcal_etaphi_time");
68  TProfile2D* h_ohcal_etaphi_time = (TProfile2D*) fin->Get("h_ohcal_etaphi_time");
69 
70  TH2F* h_cemc_e_chi2 = (TH2F*) fin->Get("h_cemc_e_chi2");
71  TH2F* h_ohcal_e_chi2 = (TH2F*) fin->Get("h_ohcal_e_chi2");
72  TH2F* h_ihcal_e_chi2 = (TH2F*) fin->Get("h_ihcal_e_chi2");
73 
74  TCanvas* c1 = new TCanvas("c1", "c1", 400, 400);
75  h_emcal_mbd_correlation ->Draw("COLZ");
76  h_emcal_mbd_correlation ->SetXTitle("#Sigma #it{E}^{EMCal} [Arb]");
77  h_emcal_mbd_correlation ->SetYTitle("#Sigma #it{E}^{MBD} [Arb]");
78  h_emcal_mbd_correlation->GetXaxis()->SetNdivisions(505);
79  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
80  myText(0.22, 0.85, 1, Form("run %d", run));
81  gPad->SetRightMargin(0.15);
82  gPad->SetLogz();
83 
84  c1->SaveAs(Form("../plots/emcal_mbd_correlation_%d.pdf",run));
85 
86 
87  TCanvas* c2 = new TCanvas("c2", "c2", 400, 400);
88  h_ihcal_mbd_correlation ->Draw("COLZ");
89  h_ihcal_mbd_correlation ->SetXTitle("#Sigma #it{E}^{ihcal} [Arb]");
90  h_ihcal_mbd_correlation ->SetYTitle("#Sigma #it{E}^{MBD} [Arb]");
91  h_ihcal_mbd_correlation->GetXaxis()->SetNdivisions(505);
92  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
93  myText(0.22, 0.85, 1, Form("run %d", run));
94  gPad->SetRightMargin(0.15);
95  gPad->SetLogz();
96 
97  c2->SaveAs(Form("../plots/ihcal_mbd_correlation_%d.pdf",run));
98 
99  TCanvas* c3 = new TCanvas("c3", "c3", 400, 400);
100  h_ohcal_mbd_correlation ->Draw("COLZ");
101  h_ohcal_mbd_correlation ->SetXTitle("#Sigma #it{E}^{ohcal} [Arb]");
102  h_ohcal_mbd_correlation ->SetYTitle("#Sigma #it{E}^{MBD} [Arb]");
103  h_ohcal_mbd_correlation->GetXaxis()->SetNdivisions(505);
104  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
105  myText(0.22, 0.85, 1, Form("run %d", run));
106  gPad->SetRightMargin(0.15);
107  gPad->SetLogz();
108 
109  c3->SaveAs(Form("../plots/ohcal_mbd_correlation_%d.pdf",run));
110 
111 
112  TCanvas* c4 = new TCanvas("c4", "c4", 400, 400);
113  h_emcal_hcal_correlation->Draw("COLZ");
114  h_emcal_hcal_correlation->SetXTitle("#Sigma #it{E}^{EMCal} [Arb]");
115  h_emcal_hcal_correlation->SetYTitle("#Sigma #it{E}^{HCal} [Arb]");
116  h_emcal_hcal_correlation->GetXaxis()->SetNdivisions(505);
117  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
118  myText(0.22, 0.85, 1, Form("run %d", run));
119  gPad->SetRightMargin(0.15);
120  gPad->SetLogz();
121 
122  c4->SaveAs(Form("../plots/emcal_hcal_correlation_%d.pdf",run));
123 
124 
125  TCanvas* c5 = new TCanvas("c5", "c5", 400, 400);
126  h_InvMass ->Draw("");
127  h_InvMass ->SetXTitle("#it{M}_{#gamma#gamma}");
128  h_InvMass ->SetYTitle("counts");
129  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
130  myText(0.22, 0.85, 1, Form("run %d", run));
131 
132  c5->SaveAs(Form("../plots/InvMass_%d.pdf",run));
133 
134 
135  TCanvas* c6 = new TCanvas("c6", "c6", 400, 400);
136  h_cemc_etaphi ->Draw("COLZ");
137  h_cemc_etaphi ->SetXTitle("#it{#eta}_{i} EMCal");
138  h_cemc_etaphi ->SetYTitle("#it{#phi}_{i} EMCal");
139  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
140  myText(0.22, 0.85, 1, Form("run %d", run));
141  myText(0.22, 0.80,1,"EMCal");
142  gPad->SetRightMargin(0.15);
143 
144  c6->SaveAs(Form("../plots/cemc_etaphi_%d.pdf",run));
145 
146 
147 
148 
149  TCanvas* c7 = new TCanvas("c7", "c7", 400, 400);
150  h_ihcal_etaphi ->Draw("COLZ");
151  h_ihcal_etaphi ->SetXTitle("#it{#eta}_{i} iHcal");
152  h_ihcal_etaphi ->SetYTitle("#it{#phi}_{i} iHcal");
153  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
154  myText(0.22, 0.85, 1, Form("run %d", run));
155  myText(0.22, 0.80,1,"iHCal");
156  gPad->SetRightMargin(0.15);
157 
158  c7->SaveAs(Form("../plots/ihcal_etaphi_%d.pdf",run));
159 
160 
161  TCanvas* c8 = new TCanvas("c8", "c8", 400, 400);
162  h_ohcal_etaphi ->Draw("COLZ");
163  h_ohcal_etaphi ->SetXTitle("#it{#eta}_{i} oHcal");
164  h_ohcal_etaphi ->SetYTitle("#it{#phi}_{i} oHcal");
165  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
166  myText(0.22, 0.85, 1, Form("run %d", run));
167  myText(0.22, 0.80,1,"oHCal");
168  gPad->SetRightMargin(0.15);
169 
170  c8->SaveAs(Form("../plots/ohcal_etaphi_%d.pdf",run));
171 
172 
173  TH1F* h_emcal_proj = (TH1F*) proj(h_cemc_etaphi)->Clone("h_emcal_proj");
174 
175  TCanvas* c9 = new TCanvas("c9", "c9", 400, 400);
176  h_emcal_proj->Draw("hist");
177  h_emcal_proj->SetYTitle("N^{twr}(E_{T} > 1 GeV)");
178 
179  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
180  myText(0.22, 0.85, 1, Form("run %d", run));
181  myText(0.22, 0.80,1,"EMCal");
182 
183  c9->SaveAs(Form("../plots/cemc_proj_%d.pdf",run));
184 
185 
186  TH1F* h_ohcal_proj = (TH1F*) proj(h_ohcal_etaphi)->Clone("h_ohcal_proj");
187 
188  TCanvas* c10 = new TCanvas("c10", "c10", 400, 400);
189  h_ohcal_proj->Draw("hist");
190  h_ohcal_proj->SetYTitle("N^{twr}(E_{T} > 1 GeV)");
191  h_ohcal_proj->GetYaxis()->SetRangeUser(0, h_ohcal_proj->GetMaximum()*1.05);
192 
193  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
194  myText(0.22, 0.85, 1, Form("run %d", run));
195  myText(0.22, 0.80,1,"oHCal");
196 
197  c10->SaveAs(Form("../plots/ohcal_proj_%d.pdf",run));
198 
199 
200  TH1F* h_ihcal_proj = (TH1F*) proj(h_ihcal_etaphi)->Clone("h_ihcal_proj");
201 
202  TCanvas* c11 = new TCanvas("c11", "c10", 400, 400);
203  h_ihcal_proj->Draw("hist");
204  h_ihcal_proj->SetYTitle("N^{twr}(E_{T} > 1 GeV)");
205  h_ihcal_proj->SetXTitle("#eta_{i}");
206 
207  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
208  myText(0.22, 0.85, 1, Form("run %d", run));
209  myText(0.22, 0.80,1,"iHCal");
210 
211  c11->SaveAs(Form("../plots/ihcal_proj_%d.pdf",run));
212 
213 
214  TH1F* h_fb_ratio_emcal = FBratio(h_emcal_proj);
215  TH1F* h_fb_ratio_ihcal = FBratio(h_ihcal_proj);
216  TH1F* h_fb_ratio_ohcal = FBratio(h_ohcal_proj);
217 
218  TCanvas* c12 = new TCanvas("c12", "c12", 400, 400);
219  h_fb_ratio_emcal->Draw("ex0");
220  h_fb_ratio_emcal->SetYTitle("N^{twr}(#eta_{i})/N^{twr}(#eta_{N-i})");
221  h_fb_ratio_emcal->SetXTitle("#eta_{i}");
222  h_fb_ratio_emcal->GetYaxis()->SetRangeUser(0.1,2);
223 
224  h_fb_ratio_ohcal->Draw("ex0 same");
225  h_fb_ratio_ohcal->SetLineColor(kBlue);
226  h_fb_ratio_ohcal->SetMarkerColor(kBlue);
227  h_fb_ratio_ohcal->SetMarkerStyle(22);
228 
229  h_fb_ratio_ihcal->Draw("ex0 same");
230  h_fb_ratio_ihcal->SetLineColor(kRed);
231  h_fb_ratio_ihcal->SetMarkerColor(kRed);
232  h_fb_ratio_ihcal->SetMarkerStyle(33);
233 
234  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
235  myText(0.22, 0.85, 1, Form("run %d", run));
236  myText(0.22, 0.80,1,"EMCal");
237  myText(0.32,0.80,kBlue,"oHCal");
238  myText(0.42,0.80,kRed,"iHCal");
239 
240  c12->SaveAs(Form("../plots/h_fb_ratio_emcal_%d.pdf",run));
241 
242 
243  TCanvas* c13 = new TCanvas("c13", "c13", 400, 400);
244  h_zdc_emcal_correlation->Draw("COLZ");
245  h_zdc_emcal_correlation->SetXTitle("#Sigma #it{E}^{EMCal} [Arb]");
246  h_zdc_emcal_correlation->SetYTitle("#Sigma #it{E}^{ZDC} [Arb]");
247  h_zdc_emcal_correlation->GetXaxis()->SetNdivisions(505);
248  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
249  myText(0.22, 0.85, 1, Form("run %d", run));
250  gPad->SetRightMargin(0.15);
251  gPad->SetLogz();
252 
253  c13->SaveAs(Form("../plots/zdc_emcal_correlation_%d.pdf",run));
254 
255 
256  TCanvas* c14 = new TCanvas("c14", "c14", 400, 400);
257  hzdcNorthcalib->Draw();
258  hzdcNorthcalib->SetLineColor(kBlue);
259  hzdcNorthcalib->GetXaxis()->SetRangeUser(0.0,12000);
260  hzdcNorthcalib->SetXTitle("#Sigma #it{E}^{ZDC Side}");
261  hzdcNorthcalib->SetYTitle("Events");
262  hzdcNorthcalib->GetXaxis()->SetNdivisions(505);
263 
264  hzdcSouthcalib->Draw("same");
265  hzdcSouthcalib->SetLineColor(kRed);
266  gPad->SetLogy();
267 
268  myText(0.7, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
269  myText(0.7, 0.85, 1, Form("run %d", run));
270  myText(0.75, 0.80, kBlue,"North");
271  myText(0.65, 0.80, kRed,"South");
272 
273  c14->SaveAs(Form("../plots/zdc_e_northSouth_%d.pdf",run));
274 
275 
276 
277 
278  TCanvas* c15 = new TCanvas("c15", "c15", 400, 400);
279  h_etaphi_clus ->Draw("COLZ");
280  h_etaphi_clus ->SetXTitle("#it{#eta}_{i} EMCal Clusters");
281  h_etaphi_clus ->SetYTitle("#it{#phi}_{i} EMCal Clusters");
282  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
283  myText(0.22, 0.85, 1, Form("run %d", run));
284  gPad->SetRightMargin(0.15);
285 
286  c15->SaveAs(Form("../plots/etaphi_clus_%d.pdf",run));
287 
288  totalEvents += hvtx_z_raw->GetEntries();
289  int events = hvtx_z_raw->GetEntries();
290 
291  TCanvas* c16 = new TCanvas("c16", "c16", 400, 400);
292  hvtx_z_raw ->Draw("COLZ");
293  hvtx_z_raw ->SetXTitle("MBD Vertex #it{z} [cm]");
294  hvtx_z_raw ->SetYTitle("Events");
295  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
296  myText(0.22, 0.85, 1, Form("run %d", run));
297  myText(0.22, 0.80, 1, Form("events %d", events));
298 
299  c16->SaveAs(Form("../plots/vtx_z_%d.pdf",run));
300 
301 
302  TCanvas* c17 = new TCanvas("c17", "c17", 400, 400);
303  hzdcNorthcalib->Draw();
304  hzdcNorthcalib->SetLineColor(kBlue);
305  hzdcNorthcalib->GetXaxis()->SetRangeUser(10,300);
306  hzdcNorthcalib->SetXTitle("#Sigma #it{E}^{ZDC Side}");
307  hzdcNorthcalib->SetYTitle("Events");
308 
309  TGraph* gr_1n = new TGraph();
310  gr_1n->SetPoint(0,100,0);
311  gr_1n->SetPoint(1,100,1e7);
312  gr_1n->SetLineStyle(7);
313  gr_1n->Draw("l");
314 
315  hzdcSouthcalib->Draw("same");
316  hzdcSouthcalib->SetLineColor(kRed);
317  gPad->SetLogy();
318 
319  myText(0.7, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
320  myText(0.7, 0.85, 1, Form("run %d", run));
321  myText(0.75, 0.80, kBlue,"North");
322  myText(0.65, 0.80, kRed,"South");
323 
324  c17->SaveAs(Form("../plots/zdc_e_northSouth_1n_%d.pdf",run));
325 
326 
327  TCanvas* c18 = new TCanvas("c18", "c18", 400, 400);
328  h_cemc_etaphi_time->Draw("COLZ");
329  h_cemc_etaphi_time->SetXTitle("#it{#eta}_{i} EMCal");
330  h_cemc_etaphi_time->SetYTitle("#it{#phi}_{i} EMCal");
331  h_cemc_etaphi_time->GetXaxis()->SetNdivisions(505);
332 
333  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
334  myText(0.22, 0.85, 1, Form("run %d", run));
335  myText(0.22, 0.80, 1, "mean hit peak time EMCal");
336  gPad->SetRightMargin(0.15);
337 
338  c18->SaveAs(Form("../plots/cemc_etaphi_time%d.pdf",run));
339 
340  TCanvas* c19 = new TCanvas("c19", "c19", 400, 400);
341  h_ohcal_etaphi_time->Draw("COLZ");
342  h_ohcal_etaphi_time->SetXTitle("#it{#eta}_{i} oHcal");
343  h_ohcal_etaphi_time->SetYTitle("#it{#phi}_{i} oHcal");
344  h_ohcal_etaphi_time->GetXaxis()->SetNdivisions(505);
345 
346  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
347  myText(0.22, 0.85, 1, Form("run %d", run));
348  myText(0.22, 0.80, 1, "mean hit peak time Outter HCal");
349  gPad->SetRightMargin(0.15);
350 
351  c19->SaveAs(Form("../plots/ohcal_etaphi_time%d.pdf",run));
352 
353 
354  TCanvas* c20 = new TCanvas("c20", "c20", 400, 400);
355  h_ihcal_etaphi_time->Draw("COLZ");
356  h_ihcal_etaphi_time->SetXTitle("#it{#eta}_{i} iHcal");
357  h_ihcal_etaphi_time->SetYTitle("#it{#phi}_{i} iHcal");
358  h_ihcal_etaphi_time->GetXaxis()->SetNdivisions(505);
359 
360  myText(0.22, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
361  myText(0.22, 0.85, 1, Form("run %d", run));
362  myText(0.22, 0.80, 1, "mean hit peak time Inner HCal");
363  gPad->SetRightMargin(0.15);
364 
365  c20->SaveAs(Form("../plots/ihcal_etaphi_time%d.pdf",run));
366 
367  TCanvas* c21 = new TCanvas("c21", "c21", 400, 400);
368  h_ihcal_e_chi2->Draw("COLZ");
369  h_ihcal_e_chi2->SetXTitle("#i{E} [GeV] iHcal");
370  h_ihcal_e_chi2->SetYTitle("chi2 iHcal");
371  h_ihcal_e_chi2->GetXaxis()->SetNdivisions(505);
372  h_ihcal_e_chi2->GetXaxis()->SetRangeUser(-1,2);
373  gPad->SetLogy();
374  gPad->SetLogz();
375 
376  myText(0.52, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
377  myText(0.52, 0.85, 1, Form("run %d", run));
378  myText(0.52, 0.80, 1, "Inner HCal");
379  gPad->SetRightMargin(0.15);
380 
381  c21->SaveAs(Form("../plots/ihcal_e_chi2%d.pdf",run));
382 
383  TCanvas* c22 = new TCanvas("c22", "c22", 400, 400);
384  h_ohcal_e_chi2->Draw("COLZ");
385  h_ohcal_e_chi2->SetXTitle("#i{E} [GeV] oHcal");
386  h_ohcal_e_chi2->SetYTitle("chi2 oHcal");
387  h_ohcal_e_chi2->GetXaxis()->SetNdivisions(505);
388  h_ohcal_e_chi2->GetXaxis()->SetRangeUser(-1,7);
389  gPad->SetLogy();
390  gPad->SetLogz();
391 
392  myText(0.52, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
393  myText(0.52, 0.85, 1, Form("run %d", run));
394  myText(0.52, 0.80, 1, "Outer HCal");
395  gPad->SetRightMargin(0.15);
396 
397  c22->SaveAs(Form("../plots/ohcal_e_chi2%d.pdf",run));
398 
399 
400  TCanvas* c23 = new TCanvas("c23", "c23", 400, 400);
401  h_cemc_e_chi2->Draw("COLZ");
402  h_cemc_e_chi2->SetXTitle("#i{E} [GeV] EMCal ");
403  h_cemc_e_chi2->SetYTitle("chi2 oHcal");
404  h_cemc_e_chi2->GetXaxis()->SetNdivisions(505);
405  h_cemc_e_chi2->GetXaxis()->SetRangeUser(-1,15);
406  gPad->SetLogy();
407  gPad->SetLogz();
408 
409  myText(0.52, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
410  myText(0.52, 0.85, 1, Form("run %d", run));
411  myText(0.52, 0.80, 1, "EMCal");
412  gPad->SetRightMargin(0.15);
413 
414  c22->SaveAs(Form("../plots/cemc_e_chi2%d.pdf",run));
415 
416  }
417 
418  for (int ir=0; ir<Nruns; ir++){
419  cout << runList[ir] << ",";
420  }
421  cout << endl;
422  cout << "total events=" << totalEvents <<endl;
423  cout << endl;
424 
425 }
426 
427 
428 
429 
430 
431 
432 TH1F* proj(TH2F* h2){
433 
434  int x = h2->GetXaxis()->GetNbins();
435  int y = h2->GetYaxis()->GetNbins();
436  TH1F* h = (TH1F*) h2->ProjectionX("temp");
437  h->Reset();
438  for (int ix=1; ix<x+1; ix++){
439  float sum = 0;
440  float cc = 0;
441  for(int iy=1; iy<y+1; iy++){
442  float bc = h2->GetBinContent(ix,iy);
443  if (bc==0) cc+=1;
444  sum += bc;
445  }
446  if (cc==y) continue;
447  float sum_cor = sum*y/(y-cc);
448  h->SetBinContent(ix,sum_cor);
449  h->SetBinError(ix,sqrt(sum_cor));
450  }
451 
452 
453  return h;
454 }
455