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<Nruns; 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 
384  TCanvas* c22 = new TCanvas("c22", "c22", 400, 400);
385  h_ohcal_e_chi2->Draw("COLZ");
386  h_ohcal_e_chi2->SetXTitle("#i{E} [GeV] oHcal");
387  h_ohcal_e_chi2->SetYTitle("chi2 oHcal");
388  h_ohcal_e_chi2->GetXaxis()->SetNdivisions(505);
389  h_ohcal_e_chi2->GetXaxis()->SetRangeUser(-1,7);
390  gPad->SetLogy();
391  gPad->SetLogz();
392 
393  myText(0.52, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
394  myText(0.52, 0.85, 1, Form("run %d", run));
395  myText(0.52, 0.80, 1, "Outer HCal");
396  gPad->SetRightMargin(0.15);
397 
398  c22->SaveAs(Form("../plots/ohcal_e_chi2%d.pdf",run));
399 
400 
401  TCanvas* c23 = new TCanvas("c23", "c23", 400, 400);
402  h_cemc_e_chi2->Draw("COLZ");
403  h_cemc_e_chi2->SetXTitle("#i{E} [GeV] EMCal ");
404  h_cemc_e_chi2->SetYTitle("chi2 EMCal");
405  h_cemc_e_chi2->GetXaxis()->SetNdivisions(505);
406  h_cemc_e_chi2->GetXaxis()->SetRangeUser(-1,15);
407  gPad->SetLogy();
408  gPad->SetLogz();
409 
410  myText(0.52, 0.9, 1, "#it{#bf{sPHENIX}} Internal");
411  myText(0.52, 0.85, 1, Form("run %d", run));
412  myText(0.52, 0.80, 1, "EMCal");
413  gPad->SetRightMargin(0.15);
414 
415  c23->SaveAs(Form("../plots/cemc_e_chi2%d.pdf",run));
416 
417  }
418 
419  for (int ir=0; ir<Nruns; ir++){
420  cout << runList[ir] << ",";
421  }
422  cout << endl;
423  cout << "total events=" << totalEvents <<endl;
424  cout << endl;
425 
426 }
427 
428 
429 
430 
431 
432 
433 TH1F* proj(TH2F* h2){
434 
435  int x = h2->GetXaxis()->GetNbins();
436  int y = h2->GetYaxis()->GetNbins();
437  TH1F* h = (TH1F*) h2->ProjectionX("temp");
438  h->Reset();
439  for (int ix=1; ix<x+1; ix++){
440  float sum = 0;
441  float cc = 0;
442  for(int iy=1; iy<y+1; iy++){
443  float bc = h2->GetBinContent(ix,iy);
444  if (bc==0) cc+=1;
445  sum += bc;
446  }
447  if (cc==y) continue;
448  float sum_cor = sum*y/(y-cc);
449  h->SetBinContent(ix,sum_cor);
450  h->SetBinError(ix,sqrt(sum_cor));
451  }
452 
453 
454  return h;
455 }
456