Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-18 09:16:47 

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 <TProfile.h>
0010 #include <TLine.h>
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 extractLevelingEnergyLinearity()
0028 {
0029   string input_5GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_1087.root";
0030   TFile *File_5GeV = TFile::Open(input_5GeV.c_str());
0031   // TH2F *h_mAsymmEnergy_pion_5GeV = (TH2F*)File_5GeV->Get("h_mAsymmEnergy_electron");
0032   TH2F *h_mAsymmEnergy_pion_5GeV = (TH2F*)File_5GeV->Get("h_mAsymmEnergy_pion_leveling");
0033   TH1F *h_mEnergy_pion_5GeV = (TH1F*)h_mAsymmEnergy_pion_5GeV->ProjectionY()->Clone();
0034 
0035   string input_8GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0422.root";
0036   TFile *File_8GeV = TFile::Open(input_8GeV.c_str());
0037   // TH2F *h_mAsymmEnergy_pion_8GeV = (TH2F*)File_8GeV->Get("h_mAsymmEnergy_electron");
0038   TH2F *h_mAsymmEnergy_pion_8GeV = (TH2F*)File_8GeV->Get("h_mAsymmEnergy_pion_leveling");
0039   TH1F *h_mEnergy_pion_8GeV = (TH1F*)h_mAsymmEnergy_pion_8GeV->ProjectionY()->Clone();
0040 
0041   string input_12GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0571.root";
0042   TFile *File_12GeV = TFile::Open(input_12GeV.c_str());
0043   // TH2F *h_mAsymmEnergy_pion_12GeV = (TH2F*)File_12GeV->Get("h_mAsymmEnergy_electron");
0044   TH2F *h_mAsymmEnergy_pion_12GeV = (TH2F*)File_12GeV->Get("h_mAsymmEnergy_pion_leveling");
0045   TH1F *h_mEnergy_pion_12GeV = (TH1F*)h_mAsymmEnergy_pion_12GeV->ProjectionY()->Clone();
0046 
0047   string input_60GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0821.root";
0048   TFile *File_60GeV = TFile::Open(input_60GeV.c_str());
0049   // TH2F *h_mAsymmEnergy_pion_60GeV = (TH2F*)File_60GeV->Get("h_mAsymmEnergy_electron");
0050   TH2F *h_mAsymmEnergy_pion_60GeV = (TH2F*)File_60GeV->Get("h_mAsymmEnergy_pion");
0051   TH1F *h_mEnergy_pion_60GeV = (TH1F*)h_mAsymmEnergy_pion_60GeV->ProjectionY()->Clone();
0052 
0053   TCanvas *c_AsymmEnergy = new TCanvas("c_AsymmEnergy","c_AsymmEnergy",1600,400);
0054   c_AsymmEnergy->Divide(4,1);
0055   for(int i_pad = 0; i_pad < 4; ++i_pad)
0056   {
0057     c_AsymmEnergy->cd(i_pad+1);
0058     c_AsymmEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
0059     c_AsymmEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
0060     c_AsymmEnergy->cd(i_pad+1)->SetTicks(1,1);
0061     c_AsymmEnergy->cd(i_pad+1)->SetGrid(0,0);
0062   }
0063   c_AsymmEnergy->cd(1);
0064   h_mAsymmEnergy_pion_5GeV->Draw("colz");
0065   h_mAsymmEnergy_pion_5GeV->ProfileX("p_mAsymmEnergy_pion_5GeV",1,-1,"i");
0066   p_mAsymmEnergy_pion_5GeV->SetMarkerStyle(20);
0067   p_mAsymmEnergy_pion_5GeV->SetMarkerColor(1);
0068   p_mAsymmEnergy_pion_5GeV->SetMarkerSize(1.0);
0069   p_mAsymmEnergy_pion_5GeV->Draw("pE same");
0070   
0071   c_AsymmEnergy->cd(2);
0072   h_mAsymmEnergy_pion_8GeV->Draw("colz");
0073   h_mAsymmEnergy_pion_8GeV->ProfileX("p_mAsymmEnergy_pion_8GeV",1,-1,"i");
0074   p_mAsymmEnergy_pion_8GeV->SetMarkerStyle(20);
0075   p_mAsymmEnergy_pion_8GeV->SetMarkerColor(1);
0076   p_mAsymmEnergy_pion_8GeV->SetMarkerSize(1.0);
0077   p_mAsymmEnergy_pion_8GeV->Draw("pE same");
0078   
0079   c_AsymmEnergy->cd(3);
0080   h_mAsymmEnergy_pion_12GeV->Draw("colz");
0081   h_mAsymmEnergy_pion_12GeV->ProfileX("p_mAsymmEnergy_pion_12GeV",1,-1,"i");
0082   p_mAsymmEnergy_pion_12GeV->SetMarkerStyle(20);
0083   p_mAsymmEnergy_pion_12GeV->SetMarkerColor(1);
0084   p_mAsymmEnergy_pion_12GeV->SetMarkerSize(1.0);
0085   p_mAsymmEnergy_pion_12GeV->Draw("pE same");
0086   
0087   c_AsymmEnergy->cd(4);
0088   h_mAsymmEnergy_pion_60GeV->Draw("colz");
0089   h_mAsymmEnergy_pion_60GeV->ProfileX("p_mAsymmEnergy_pion_60GeV",1,-1,"i");
0090   p_mAsymmEnergy_pion_60GeV->SetMarkerStyle(20);
0091   p_mAsymmEnergy_pion_60GeV->SetMarkerColor(1);
0092   p_mAsymmEnergy_pion_60GeV->SetMarkerSize(1.0);
0093   p_mAsymmEnergy_pion_60GeV->Draw("pE same");
0094   
0095   float mean_energy[4] = {5.0,8.0,12.0,60.0};
0096   float val_mean[4];
0097   float err_mean[4];
0098   float val_sigma[3];
0099   float err_sigma[3];
0100   float val_resolution[3];
0101   float err_resolution[3];
0102 
0103   TCanvas *c_Energy = new TCanvas("c_Energy","c_Energy",1600,400);
0104   c_Energy->Divide(4,1);
0105   for(int i_pad = 0; i_pad < 4; ++i_pad)
0106   {
0107     c_Energy->cd(i_pad+1);
0108     c_Energy->cd(i_pad+1)->SetLeftMargin(0.15);
0109     c_Energy->cd(i_pad+1)->SetBottomMargin(0.15);
0110     c_Energy->cd(i_pad+1)->SetTicks(1,1);
0111     c_Energy->cd(i_pad+1)->SetGrid(0,0);
0112   }
0113   c_Energy->cd(1);
0114   h_mEnergy_pion_5GeV->Draw("hE");
0115   TF1 *f_gaus_5GeV = new TF1("f_gaus_5GeV","gaus",0,100);
0116   f_gaus_5GeV->SetParameter(0,1.0);
0117   f_gaus_5GeV->SetParameter(1,h_mEnergy_pion_5GeV->GetMean());
0118   f_gaus_5GeV->SetParameter(2,1.0);
0119   f_gaus_5GeV->SetRange(1.2,2.8);
0120   h_mEnergy_pion_5GeV->Fit(f_gaus_5GeV,"NR");
0121   f_gaus_5GeV->SetLineColor(2);
0122   f_gaus_5GeV->SetLineStyle(2);
0123   f_gaus_5GeV->SetLineWidth(2);
0124   f_gaus_5GeV->Draw("l same");
0125   val_mean[0]       = f_gaus_5GeV->GetParameter(1);
0126   err_mean[0]       = f_gaus_5GeV->GetParError(1);
0127   val_sigma[0]      = f_gaus_5GeV->GetParameter(2); // extract calibrated energy
0128   err_sigma[0]      = f_gaus_5GeV->GetParError(2);
0129   val_resolution[0] = val_sigma[0]/val_mean[0];
0130   err_resolution[0] = ErrDiv(val_sigma[0],val_mean[0],err_sigma[0],err_mean[0]);
0131   
0132   c_Energy->cd(2);
0133   h_mEnergy_pion_8GeV->Draw("hE");
0134   TF1 *f_gaus_8GeV = new TF1("f_gaus_8GeV","gaus",0,100);
0135   f_gaus_8GeV->SetParameter(0,1.0);
0136   f_gaus_8GeV->SetParameter(1,h_mEnergy_pion_8GeV->GetMean());
0137   f_gaus_8GeV->SetParameter(2,1.0);
0138   f_gaus_8GeV->SetRange(1.0,4.2);
0139   h_mEnergy_pion_8GeV->Fit(f_gaus_8GeV,"NR");
0140   f_gaus_8GeV->SetLineColor(2);
0141   f_gaus_8GeV->SetLineStyle(2);
0142   f_gaus_8GeV->SetLineWidth(2);
0143   f_gaus_8GeV->Draw("l same");
0144   val_mean[1]       = f_gaus_8GeV->GetParameter(1);
0145   err_mean[1]       = f_gaus_8GeV->GetParError(1);
0146   val_sigma[1]      = f_gaus_8GeV->GetParameter(2);
0147   err_sigma[1]      = f_gaus_8GeV->GetParError(2);
0148   val_resolution[1] = val_sigma[1]/val_mean[1];
0149   err_resolution[1] = ErrDiv(val_sigma[1],val_mean[1],err_sigma[1],err_mean[1]);
0150   
0151   c_Energy->cd(3);
0152   h_mEnergy_pion_12GeV->Draw("hE");
0153   TF1 *f_gaus_12GeV = new TF1("f_gaus_12GeV","gaus",0,100);
0154   f_gaus_12GeV->SetParameter(0,1.0);
0155   f_gaus_12GeV->SetParameter(1,h_mEnergy_pion_12GeV->GetMean());
0156   f_gaus_12GeV->SetParameter(2,1.0);
0157   f_gaus_12GeV->SetRange(1.5,7.0);
0158   h_mEnergy_pion_12GeV->Fit(f_gaus_12GeV,"NR");
0159   f_gaus_12GeV->SetLineColor(2);
0160   f_gaus_12GeV->SetLineStyle(2);
0161   f_gaus_12GeV->SetLineWidth(2);
0162   f_gaus_12GeV->Draw("l same");
0163   val_mean[2]       = f_gaus_12GeV->GetParameter(1);
0164   err_mean[2]       = f_gaus_12GeV->GetParError(1);
0165   val_sigma[2]      = f_gaus_12GeV->GetParameter(2);
0166   err_sigma[2]      = f_gaus_12GeV->GetParError(2);
0167   val_resolution[2] = val_sigma[2]/val_mean[2];
0168   err_resolution[2] = ErrDiv(val_sigma[2],val_mean[2],err_sigma[2],err_mean[2]);
0169   
0170   c_Energy->cd(4);
0171   h_mEnergy_pion_60GeV->Draw("hE");
0172   TF1 *f_gaus_60GeV = new TF1("f_gaus_60GeV","gaus",0,100);
0173   f_gaus_60GeV->SetParameter(0,1.0);
0174   f_gaus_60GeV->SetParameter(1,h_mEnergy_pion_60GeV->GetMean());
0175   f_gaus_60GeV->SetParameter(2,1.0);
0176   f_gaus_60GeV->SetRange(10.2,30.2);
0177   h_mEnergy_pion_60GeV->Fit(f_gaus_60GeV,"NR");
0178   f_gaus_60GeV->SetLineColor(2);
0179   f_gaus_60GeV->SetLineStyle(2);
0180   f_gaus_60GeV->SetLineWidth(2);
0181   f_gaus_60GeV->Draw("l same");
0182   val_mean[3]       = f_gaus_60GeV->GetParameter(1);
0183   err_mean[3]       = f_gaus_60GeV->GetParError(1);
0184 
0185   c_Energy->SaveAs("../figures/HCAL_ShowerCalib/c_Energy_LevelingCorr.eps");
0186   
0187   TGraphAsymmErrors *g_lieanrity = new TGraphAsymmErrors();
0188   TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
0189   for(int i_point = 0; i_point < 4; ++i_point)
0190   {
0191     g_lieanrity->SetPoint(i_point,mean_energy[i_point],val_mean[i_point]);
0192     g_lieanrity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
0193 
0194     g_resolution->SetPoint(i_point,mean_energy[i_point],val_resolution[i_point]);
0195     g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
0196   }
0197 
0198   TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
0199   c_Linearity->cd();
0200   c_Linearity->cd()->SetLeftMargin(0.15);
0201   c_Linearity->cd()->SetBottomMargin(0.15);
0202   c_Linearity->cd()->SetTicks(1,1);
0203   c_Linearity->cd()->SetGrid(0,0);
0204 
0205   TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
0206   for(int i_bin = 0; i_bin < 100; ++i_bin)
0207   {
0208     h_play->SetBinContent(i_bin+1,-10.0);
0209     h_play->SetBinError(i_bin+1,1.0);
0210   }
0211   h_play->SetTitle("");
0212   h_play->SetStats(0);
0213   h_play->GetXaxis()->SetTitle("input Energy (GeV)");
0214   h_play->GetXaxis()->CenterTitle();
0215   h_play->GetXaxis()->SetNdivisions(505);
0216   h_play->GetXaxis()->SetRangeUser(0.0,80.0);
0217 
0218   h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
0219   h_play->GetYaxis()->CenterTitle();
0220   h_play->GetYaxis()->SetRangeUser(0.0,40.0);
0221   h_play->DrawCopy("pE");
0222 
0223   g_lieanrity->SetMarkerStyle(24);
0224   g_lieanrity->SetMarkerColor(2);
0225   g_lieanrity->SetMarkerSize(1.2);
0226   g_lieanrity->Draw("PE same");
0227 
0228   TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
0229   f_pol->SetParameter(0,0.0);
0230   f_pol->SetParameter(1,1.0);
0231   f_pol->SetRange(0.0,80.0);
0232   g_lieanrity->Fit(f_pol,"NR");
0233   f_pol->SetLineColor(2);
0234   f_pol->SetLineStyle(2);
0235   f_pol->SetLineWidth(2);
0236   f_pol->Draw("l same");
0237 
0238   c_Linearity->SaveAs("../figures/HCAL_ShowerCalib/c_Linearity_LevelingCorr.eps");
0239   
0240   /*
0241   TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
0242   l_unity->SetLineColor(4);
0243   l_unity->SetLineStyle(2);
0244   l_unity->SetLineWidth(2);
0245   l_unity->Draw("l same");
0246   */
0247 
0248   TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
0249   c_Resolution->cd();
0250   c_Resolution->cd()->SetLeftMargin(0.15);
0251   c_Resolution->cd()->SetBottomMargin(0.15);
0252   c_Resolution->cd()->SetTicks(1,1);
0253   c_Resolution->cd()->SetGrid(0,0);
0254 
0255   h_play->GetYaxis()->SetTitle("$DeltaE/<E>");
0256   h_play->GetYaxis()->CenterTitle();
0257   h_play->GetYaxis()->SetRangeUser(0.0,0.8);
0258   h_play->DrawCopy("pE");
0259 
0260   g_resolution->SetMarkerStyle(24);
0261   g_resolution->SetMarkerColor(2);
0262   g_resolution->SetMarkerSize(1.2);
0263   g_resolution->Draw("PE same");
0264 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!