File indexing completed on 2025-08-06 08:15:48
0001 #include <string>
0002 #include <TString.h>
0003 #include <TFile.h>
0004 #include <TH2F.h>
0005 #include <TH1F.h>
0006 #include <TF1.h>
0007 #include <TCanvas.h>
0008 #include <TGraphAsymmErrors.h>
0009 #include <iostream>
0010 #include <fstream>
0011
0012 float ErrorAdd(float x, float y)
0013 {
0014 return sqrt(x*x+y*y);
0015 }
0016
0017 float ErrTimes(float x, float y, float dx, float dy)
0018 {
0019 return x*y*ErrorAdd(dx/x,dy/y);
0020 }
0021
0022 float ErrDiv(float x, float y, float dx, float dy)
0023 {
0024 return x/y*ErrorAdd(dx/x,dy/y);
0025 }
0026
0027 void extractMIPEnergy()
0028 {
0029 string inputenergy[4] = {"5GeV","8GeV","12GeV","60GeV"};
0030 string runId[4] = {"1087","0422","0571","0821"};
0031 float mean_energy[4] = {5.0,8.0,12.0,60.0};
0032
0033 TH2F *h_mAsymmEnergy_pion[4];
0034 TH1F *h_mMIPEnergy_pion[4];
0035 TFile *File_InPut[4];
0036
0037 for(int i_energy = 0; i_energy < 4; ++i_energy)
0038 {
0039 string inputfile = Form("/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_%s.root",runId[i_energy].c_str());
0040 File_InPut[i_energy] = TFile::Open(inputfile.c_str());
0041 h_mAsymmEnergy_pion[i_energy] = (TH2F*)File_InPut[i_energy]->Get("h_mAsymmEnergy_pion");
0042 h_mMIPEnergy_pion[i_energy] = (TH1F*)h_mAsymmEnergy_pion[i_energy]->ProjectionY()->Clone();
0043 }
0044
0045 TF1 *f_gaus_MIP[4];
0046 float MIP_mean[4];
0047 float MIP_width[4];
0048 TCanvas *c_MIPEnergy = new TCanvas("c_MIPEnergy","c_MIPEnergy",1500,500);
0049 c_MIPEnergy->Divide(3,1);
0050 for(int i_pad = 0; i_pad < 3; ++i_pad)
0051 {
0052 c_MIPEnergy->cd(i_pad+1);
0053 c_MIPEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
0054 c_MIPEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
0055 c_MIPEnergy->cd(i_pad+1)->SetTicks(1,1);
0056 c_MIPEnergy->cd(i_pad+1)->SetGrid(0,0);
0057 h_mMIPEnergy_pion[i_pad]->SetTitle(inputenergy[i_pad].c_str());
0058 h_mMIPEnergy_pion[i_pad]->GetXaxis()->SetTitle("Tower Calibrated Energy (GeV)");
0059 h_mMIPEnergy_pion[i_pad]->Draw("hE");
0060 string FuncName = Form("f_gaus_MIP_%d",i_pad);
0061 f_gaus_MIP[i_pad] = new TF1(FuncName.c_str(),"gaus",0,5);
0062 f_gaus_MIP[i_pad]->SetParameter(0,1.0);
0063 f_gaus_MIP[i_pad]->SetParameter(1,1.0);
0064 f_gaus_MIP[i_pad]->SetParameter(2,0.1);
0065 f_gaus_MIP[i_pad]->SetRange(0.4,0.9);
0066 h_mMIPEnergy_pion[i_pad]->Fit(f_gaus_MIP[i_pad],"NQR");
0067 f_gaus_MIP[i_pad]->SetLineColor(2);
0068 f_gaus_MIP[i_pad]->SetLineStyle(2);
0069 f_gaus_MIP[i_pad]->SetLineWidth(2);
0070 f_gaus_MIP[i_pad]->Draw("l same");
0071 MIP_mean[i_pad] = f_gaus_MIP[i_pad]->GetParameter(1);
0072 MIP_width[i_pad] = f_gaus_MIP[i_pad]->GetParameter(2);
0073 cout << "MIP for " << inputenergy[i_pad] << " is: " << MIP_mean[i_pad] << " +/- " << MIP_width[i_pad] << endl;
0074 }
0075
0076 int i_MIP = 1;
0077 cout << "default MIP for all energy chosen to be " << inputenergy[i_MIP] << " is: " << MIP_mean[i_MIP] << " +/- " << MIP_width[i_MIP] << endl;
0078
0079 ofstream File_OutPut("MIPEnergy.txt");
0080 File_OutPut << "default MIP at " << inputenergy[i_MIP] << " is: " << MIP_mean[i_MIP] << " +/- " << MIP_width[i_MIP] << endl;
0081 File_OutPut.close();
0082
0083 c_MIPEnergy->SaveAs("../figures/HCAL_ShowerCalib/c_MIPEnergy.eps");
0084
0085 #if 0
0086 TF1 *f_gaus_Reco[4];
0087 float Reco_mean[4];
0088 float err_mean[4];
0089 float Reco_width[4];
0090 float err_width[4];
0091 float Reco_resolution[4];
0092 float err_resolution[4];
0093 TCanvas *c_RecoEnergy = new TCanvas("c_RecoEnergy","c_RecoEnergy",1500,500);
0094 c_RecoEnergy->Divide(3,1);
0095 for(int i_pad = 0; i_pad < 3; ++i_pad)
0096 {
0097 c_RecoEnergy->cd(i_pad+1);
0098 c_RecoEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
0099 c_RecoEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
0100 c_RecoEnergy->cd(i_pad+1)->SetTicks(1,1);
0101 c_RecoEnergy->cd(i_pad+1)->SetGrid(0,0);
0102 h_mMIPEnergy_pion[i_pad]->Draw("hE");
0103
0104 string FuncName = Form("f_gaus_Reco_%d",i_pad);
0105 f_gaus_Reco[i_pad] = new TF1(FuncName.c_str(),"gaus",0,100);
0106 f_gaus_Reco[i_pad]->SetParameter(0,1.0);
0107 f_gaus_Reco[i_pad]->SetParameter(1,mean_energy[i_pad]/3.0);
0108 f_gaus_Reco[i_pad]->SetParameter(2,1.0);
0109 f_gaus_Reco[i_pad]->SetRange(MIP_mean[i_MIP]+2.0*MIP_width[i_MIP],0.5*mean_energy[i_pad]);
0110 h_mMIPEnergy_pion[i_pad]->Fit(f_gaus_Reco[i_pad],"NQR");
0111 f_gaus_Reco[i_pad]->SetLineColor(2);
0112 f_gaus_Reco[i_pad]->SetLineStyle(2);
0113 f_gaus_Reco[i_pad]->SetLineWidth(2);
0114 f_gaus_Reco[i_pad]->Draw("l same");
0115
0116 Reco_mean[i_pad] = f_gaus_Reco[i_pad]->GetParameter(1);
0117 err_mean[i_pad] = f_gaus_Reco[i_pad]->GetParError(1);
0118 Reco_width[i_pad] = f_gaus_Reco[i_pad]->GetParameter(2);
0119 err_width[i_pad] = f_gaus_Reco[i_pad]->GetParError(2);
0120 cout << "Reco for " << inputenergy[i_pad] << " is: " << Reco_mean[i_pad] << " +/- " << Reco_width[i_pad] << endl;
0121
0122 Reco_resolution[i_pad] = Reco_width[i_pad]/Reco_mean[i_pad];
0123 err_resolution[i_pad] = ErrDiv(Reco_width[i_pad],Reco_mean[i_pad],err_width[i_pad],err_mean[i_pad]);
0124 }
0125
0126 TGraphAsymmErrors *g_lieanrity = new TGraphAsymmErrors();
0127 TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
0128 for(int i_point = 0; i_point < 3; ++i_point)
0129 {
0130 g_lieanrity->SetPoint(i_point,mean_energy[i_point],Reco_mean[i_point]);
0131 g_lieanrity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
0132
0133 g_resolution->SetPoint(i_point,mean_energy[i_point],Reco_resolution[i_point]);
0134 g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
0135 }
0136
0137 TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
0138 c_Linearity->cd();
0139 c_Linearity->cd()->SetLeftMargin(0.15);
0140 c_Linearity->cd()->SetBottomMargin(0.15);
0141 c_Linearity->cd()->SetTicks(1,1);
0142 c_Linearity->cd()->SetGrid(0,0);
0143
0144 TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
0145 for(int i_bin = 0; i_bin < 100; ++i_bin)
0146 {
0147 h_play->SetBinContent(i_bin+1,-10.0);
0148 h_play->SetBinError(i_bin+1,1.0);
0149 }
0150 h_play->SetTitle("");
0151 h_play->SetStats(0);
0152 h_play->GetXaxis()->SetTitle("input Energy (GeV)");
0153 h_play->GetXaxis()->CenterTitle();
0154 h_play->GetXaxis()->SetNdivisions(505);
0155 h_play->GetXaxis()->SetRangeUser(0.0,20.0);
0156
0157 h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
0158 h_play->GetYaxis()->CenterTitle();
0159 h_play->GetYaxis()->SetRangeUser(0.0,10.0);
0160 h_play->DrawCopy("pE");
0161
0162 g_lieanrity->SetMarkerStyle(24);
0163 g_lieanrity->SetMarkerColor(2);
0164 g_lieanrity->SetMarkerSize(1.2);
0165 g_lieanrity->Draw("PE same");
0166
0167 TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
0168 f_pol->SetParameter(0,0.0);
0169 f_pol->SetParameter(1,1.0);
0170 f_pol->SetRange(0.0,20.0);
0171 g_lieanrity->Fit(f_pol,"NR");
0172 f_pol->SetLineColor(2);
0173 f_pol->SetLineStyle(2);
0174 f_pol->SetLineWidth(2);
0175 f_pol->Draw("l same");
0176
0177 TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
0178 string inputresolution = "resolution.txt";
0179 FILE *fres = fopen(inputresolution.c_str(), "r");
0180 if(fres == NULL)
0181 {
0182 perror("Error opening file!");
0183 }
0184 else
0185 {
0186 float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
0187 char line[1024];
0188 int line_counter = 0;
0189 while(fgets(line,1024,fres))
0190 {
0191 sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
0192
0193 g_resolution_old->SetPoint(line_counter,x,y);
0194 g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
0195 line_counter++;
0196 }
0197 }
0198 TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
0199 c_Resolution->cd();
0200 c_Resolution->cd()->SetLeftMargin(0.15);
0201 c_Resolution->cd()->SetBottomMargin(0.15);
0202 c_Resolution->cd()->SetTicks(1,1);
0203 c_Resolution->cd()->SetGrid(0,0);
0204
0205 h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
0206 h_play->GetYaxis()->CenterTitle();
0207 h_play->GetYaxis()->SetRangeUser(0.0,0.8);
0208 h_play->DrawCopy("pE");
0209
0210 g_resolution->SetMarkerStyle(20);
0211 g_resolution->SetMarkerColor(kGray+2);
0212 g_resolution->SetMarkerSize(2.0);
0213 g_resolution->Draw("pE same");
0214
0215 g_resolution_old->SetMarkerStyle(21);
0216 g_resolution_old->SetMarkerColor(2);
0217 g_resolution_old->SetMarkerSize(1.8);
0218 g_resolution_old->Draw("pE same");
0219
0220 TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
0221 leg_res->SetBorderSize(0);
0222 leg_res->SetFillColor(0);
0223 leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
0224 leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
0225 leg_res->Draw("same");
0226 #endif
0227 }
