HCal/macros/calLinearity_ShowerCalib.C

0001 #include "TSystem.h"
0002 #include "TFile.h"
0003 #include "TH1F.h"
0004 #include "TF1.h"
0005 #include "TMath.h"
0006 #include "TCanvas.h"
0007 #include "TGraphAsymmErrors.h"
0008 #include "TLine.h"
0009 #include "TLegend.h"
0010 #include "TString.h"
0011
0012 #include <iostream>
0013 #include <string>
0014 #include <cstdio>
0015 #include <cstring>
0016
0017 using namespace std;
0018
0019 float ErrorAdd(float x, float y)
0020 {
0021     return sqrt(x*x+y*y);
0022 }
0023
0024 float ErrTimes(float x, float y, float dx, float dy)
0025 {
0026     return x*y*ErrorAdd(dx/x,dy/y);
0027 }
0028
0029 float ErrDiv(float x, float y, float dx, float dy)
0030 {
0031     return x/y*ErrorAdd(dx/x,dy/y);
0032 }
0033
0034 void calLinearity_ShowerCalib()
0035 {
0036   TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
0037   TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
0038   float mBeamEnergy[6] = {4.0,8.0,12.0,16.0,24.0,30.0};
0039   string runID[6] = {"1241","0422","0571","1683","1666","1684"};
0040
0041   TFile *File_InPut[6];
0042   TH2F *h_mAsymmEnergy_pion_ShowerCalib[6];
0043   TH1F *h_mEnergy_pion_ShowerCalib[6];
0044
0045   for(int i_energy = 0; i_energy < 6; ++i_energy)
0046   {
0047     string inputfile = Form("/sphenix/user/xusun/software/data/beam/ShowerCalibAna/Proto4ShowerInfoRAW_%s.root",runID[i_energy].c_str());
0048     File_InPut[i_energy] = TFile::Open(inputfile.c_str());
0049     string title = Form("%1.0f GeV & #pi^{-}",mBeamEnergy[i_energy]);
0050
0051     h_mAsymmEnergy_pion_ShowerCalib[i_energy] = (TH2F*)File_InPut[i_energy]->Get("h_mAsymmEnergy_pion_ShowerCalib")->Clone();
0052     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->SetTitle(title.c_str());
0053     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetTitle("Energy (GeV)");
0054     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->CenterTitle();
0055     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetNdivisions(505);
0056     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetTitle("E_{Asymm}");
0057     h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->CenterTitle();
0058
0059     h_mEnergy_pion_ShowerCalib[i_energy] = (TH1F*)h_mAsymmEnergy_pion_ShowerCalib[i_energy]->ProjectionY()->Clone();
0060     h_mEnergy_pion_ShowerCalib[i_energy]->SetStats(0);
0061     h_mEnergy_pion_ShowerCalib[i_energy]->SetTitle(title.c_str());
0062     h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetTitle("Energy (GeV)");
0063     h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->CenterTitle();
0064     h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetNdivisions(505);
0065     // h_mEnergy_pion_ShowerCalib[i_energy]->SetLineColor(2);
0066     h_mEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetRangeUser(0,1.2*h_mEnergy_pion_ShowerCalib[i_energy]->GetMaximum());
0067   }
0068
0069   TCanvas *c_energy_asymm = new TCanvas("c_energy_asymm","c_energy_asymm",800,1200);
0070   c_energy_asymm->Divide(2,3);
0071   for(int i_pad = 0; i_pad < 6; ++i_pad)
0072   {
0073     c_energy_asymm->cd(i_pad+1);
0074     c_energy_asymm->cd(i_pad+1)->SetLeftMargin(0.15);
0075     c_energy_asymm->cd(i_pad+1)->SetBottomMargin(0.15);
0076     c_energy_asymm->cd(i_pad+1)->SetTicks(1,1);
0077     c_energy_asymm->cd(i_pad+1)->SetGrid(0,0);
0078     h_mAsymmEnergy_pion_ShowerCalib[i_pad]->Draw("colz");
0079   }
0080   c_energy_asymm->SaveAs("./figures/sPHENIX_CollMeeting/c_energy_asymm_ShowerCalib.eps");
0081
0082   TF1 *f_landau[6];
0083   float mpv[6], err_mpv[6];
0084   float sigma[6], err_sigma[6];
0085   TCanvas *c_showercalib = new TCanvas("c_showercalib","c_showercalib",800,1200);
0086   c_showercalib->Divide(2,3);
0087   for(int i_pad = 0; i_pad < 6; ++i_pad)
0088   {
0089     c_showercalib->cd(i_pad+1);
0090     c_showercalib->cd(i_pad+1)->SetLeftMargin(0.15);
0091     c_showercalib->cd(i_pad+1)->SetBottomMargin(0.15);
0092     c_showercalib->cd(i_pad+1)->SetTicks(1,1);
0093     c_showercalib->cd(i_pad+1)->SetGrid(0,0);
0094     h_mEnergy_pion_ShowerCalib[i_pad]->Draw();
0095     string FuncName = Form("f_landau_%1.0fGeV",mBeamEnergy[i_pad]);
0096     f_landau[i_pad] = new TF1(FuncName.c_str(),"landau",0,100);
0097     f_landau[i_pad]->SetParameter(0,1.0);
0098     f_landau[i_pad]->SetParameter(1,mBeamEnergy[i_pad]);
0099     f_landau[i_pad]->SetParameter(2,1.0);
0100     f_landau[i_pad]->SetRange(0.5*mBeamEnergy[i_pad],mBeamEnergy[i_pad]*2.0);
0101     if(i_pad > 3) f_landau[i_pad]->SetRange(0.75*mBeamEnergy[i_pad],mBeamEnergy[i_pad]*1.9);
0102     h_mEnergy_pion_ShowerCalib[i_pad]->Fit(f_landau[i_pad],"R");
0103
0104     mpv[i_pad] = f_landau[i_pad]->GetParameter(1);
0105     err_mpv[i_pad] = f_landau[i_pad]->GetParError(1);
0106     g_linearity->SetPoint(i_pad,mBeamEnergy[i_pad],mpv[i_pad]);
0107     g_linearity->SetPointError(i_pad,0.0,0.0,err_mpv[i_pad],err_mpv[i_pad]);
0108
0109     sigma[i_pad] = f_landau[i_pad]->GetParameter(2);
0110     err_sigma[i_pad] = f_landau[i_pad]->GetParError(2);
0111     float val_resolution = sigma[i_pad]/mpv[i_pad];
0112     float err_resolution = ErrDiv(sigma[i_pad],mpv[i_pad],err_sigma[i_pad],err_mpv[i_pad]);
0113     g_resolution->SetPoint(i_pad,mBeamEnergy[i_pad],val_resolution);
0114     g_resolution->SetPointError(i_pad,0.0,0.0,err_resolution,err_resolution);
0115
0116     cout << "input energy is " << mBeamEnergy[i_pad] << ", calibrated energy is " << mpv[i_pad] << " +/- " << err_mpv[i_pad] << ", sigma is " << sigma[i_pad] << endl;
0117   }
0118   c_showercalib->SaveAs("./figures/sPHENIX_CollMeeting/c_showercalib.eps");
0119
0120   TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
0121   c_Linearity->cd();
0122   c_Linearity->cd()->SetLeftMargin(0.15);
0123   c_Linearity->cd()->SetBottomMargin(0.15);
0124   c_Linearity->cd()->SetTicks(1,1);
0125   c_Linearity->cd()->SetGrid(0,0);
0126
0127   TH1F *h_play = new TH1F("h_play","h_play",100,-0.5,99.5);
0128   for(int i_bin = 0; i_bin < 100; ++i_bin)
0129   {
0130     h_play->SetBinContent(i_bin+1,-100.0);
0131     h_play->SetBinError(i_bin+1,1.0);
0132   }
0133   h_play->SetTitle("");
0134   h_play->SetStats(0);
0135
0136   h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
0137   h_play->GetXaxis()->CenterTitle();
0138   h_play->GetXaxis()->SetRangeUser(0,40);
0139
0140   h_play->GetYaxis()->SetTitle("Calibrated Energy (GeV)");
0141   h_play->GetYaxis()->CenterTitle();
0142   h_play->GetYaxis()->SetRangeUser(0,40);
0143   h_play->DrawCopy("pE");
0144
0145   TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
0146   string inputlinearity = "linearity.txt";
0147   FILE *flinear = fopen(inputlinearity.c_str(), "r");
0148   if(flinear == NULL)
0149   {
0150     perror("Error opening file!");
0151   }
0152   else
0153   {
0154     float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
0155     char line[1024];
0156     int line_counter = 0;
0157     while(fgets(line,1024,flinear))
0158     {
0159       sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
0160       // cout << "x = " << x << ", y = " << y << ", err_x_low = " << err_x_low << ", err_x_high = " << err_x_high << ", err_y_low = " << err_y_low << ", err_y_high = " << err_y_high << endl;
0161       g_linearity_old->SetPoint(line_counter,x,y);
0162       g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
0163       line_counter++;
0164     }
0165   }
0166
0167   g_linearity_old->SetMarkerStyle(21);
0168   g_linearity_old->SetMarkerColor(2);
0169   g_linearity_old->SetMarkerSize(1.8);
0170   g_linearity_old->Draw("pE same");
0171
0172   g_linearity->SetMarkerStyle(20);
0173   g_linearity->SetMarkerColor(kGray+2);
0174   g_linearity->SetMarkerSize(2.0);
0175   g_linearity->Draw("pE same");
0176
0177   TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
0178   l_unity->SetLineColor(4);
0179   l_unity->SetLineStyle(2);
0180   l_unity->SetLineWidth(2);
0181   l_unity->Draw("l same");
0182
0183   TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
0184   leg_linear->SetBorderSize(0);
0185   leg_linear->SetFillColor(0);
0186   leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
0187   leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
0188   leg_linear->AddEntry(l_unity,"unity","l");
0189   leg_linear->Draw("same");
0190
0191   c_Linearity->SaveAs("./figures/sPHENIX_CollMeeting/c_Linearity.eps");
0192
0193   TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
0194   c_Resolution->cd();
0195   c_Resolution->cd()->SetLeftMargin(0.15);
0196   c_Resolution->cd()->SetBottomMargin(0.15);
0197   c_Resolution->cd()->SetTicks(1,1);
0198   c_Resolution->cd()->SetGrid(0,0);
0199
0200   h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
0201   h_play->GetXaxis()->CenterTitle();
0202   h_play->GetXaxis()->SetRangeUser(0,40);
0203
0204   h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
0205   h_play->GetYaxis()->CenterTitle();
0206   h_play->GetYaxis()->SetRangeUser(0,0.7);
0207   h_play->DrawCopy("pE");
0208
0209   TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
0210   string inputresolution = "resolution.txt";
0211   FILE *fres = fopen(inputresolution.c_str(), "r");
0212   if(fres == NULL)
0213   {
0214     perror("Error opening file!");
0215   }
0216   else
0217   {
0218     float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
0219     char line[1024];
0220     int line_counter = 0;
0221     while(fgets(line,1024,fres))
0222     {
0223       sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
0224       // cout << "x = " << x << ", y = " << y << ", err_x_low = " << err_x_low << ", err_x_high = " << err_x_high << ", err_y_low = " << err_y_low << ", err_y_high = " << err_y_high << endl;
0225       g_resolution_old->SetPoint(line_counter,x,y);
0226       g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
0227       line_counter++;
0228     }
0229   }
0230
0231   g_resolution_old->SetMarkerStyle(21);
0232   g_resolution_old->SetMarkerColor(2);
0233   g_resolution_old->SetMarkerSize(1.8);
0234   g_resolution_old->Draw("pE same");
0235
0236   g_resolution->SetMarkerStyle(20);
0237   g_resolution->SetMarkerColor(kGray+2);
0238   g_resolution->SetMarkerSize(2.0);
0239   g_resolution->Draw("pE same");
0240
0241   TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
0242   leg_res->SetBorderSize(0);
0243   leg_res->SetFillColor(0);
0244   leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
0245   leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
0246   leg_res->Draw("same");
0247
0248   c_Resolution->SaveAs("./figures/sPHENIX_CollMeeting/c_Resolution.eps");
0249
0250 }