HCal/macros/calLinearity.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
0011 #include <iostream>
0012 #include <string>
0013 #include <cstdio>
0014 #include <cstring>
0015
0016 using namespace std;
0017
0018 float ErrorAdd(float x, float y)
0019 {
0020     return sqrt(x*x+y*y);
0021 }
0022
0023 float ErrTimes(float x, float y, float dx, float dy)
0024 {
0025     return x*y*ErrorAdd(dx/x,dy/y);
0026 }
0027
0028 float ErrDiv(float x, float y, float dx, float dy)
0029 {
0030     return x/y*ErrorAdd(dx/x,dy/y);
0031 }
0032
0033 void calLinearity()
0034 {
0035   gSystem->Load("libPrototype4.so");
0036   gSystem->Load("libProto4_HCalShowerCalib.so");
0037
0038   TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
0039   TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
0040   float mBeamEnergy[5] = {8.0,12.0,16.0,24.0,30.0};
0041
0042   string inputlist = "EnergyScan.list";
0043   FILE *fp = fopen(inputlist.c_str(), "r");
0044   if(fp == NULL)
0045   {
0046     perror("Error opening file!");
0047   }
0048   else
0049   {
0050     char line[1024];
0051     int line_counter = 0;
0052     while(fgets(line,1024,fp))
0053     {
0054       line[strcspn(line, "\n")] = '\0';
0055       string inputfile = line;
0056       cout << "read in: " << inputfile.c_str() << endl;
0057       TFile *File_HCAL = TFile::Open(inputfile.c_str());
0058       TH1F *h_BeamMom = (TH1F*)File_HCAL->Get("hBeam_Mom")->Clone();
0059       // float inputEnergy = TMath::Abs(h_BeamMom->GetMean()); // get input energy
0060       float inputEnergy = mBeamEnergy[line_counter]; // get input energy
0061
0062       TH1F *h_calib = (TH1F*)File_HCAL->Get("h_hcal_total_calib")->Clone();
0063       TF1 *f_gaus = new TF1("f_gaus","gaus",0,100);
0064       f_gaus->SetParameter(0,1.0);
0065       f_gaus->SetParameter(1,inputEnergy);
0066       f_gaus->SetParameter(2,1.0);
0067       h_calib->Fit(f_gaus,"NQ");
0068       float val_calib = f_gaus->GetParameter(1); // extract calibrated energy
0069       float err_calib = f_gaus->GetParError(1);
0070       g_linearity->SetPoint(line_counter,inputEnergy,val_calib);
0071       g_linearity->SetPointError(line_counter,0.0,0.0,err_calib,err_calib);
0072
0073       float val_sigma = f_gaus->GetParameter(2); // extract calibrated energy
0074       float err_sigma = f_gaus->GetParError(2);
0075       float val_resolution = val_sigma/val_calib;
0076       float err_resolution = ErrDiv(val_sigma,val_calib,err_sigma,err_calib);
0077       g_resolution->SetPoint(line_counter,inputEnergy,val_resolution);
0078       g_resolution->SetPointError(line_counter,0.0,0.0,err_resolution,err_resolution);
0079
0080       line_counter++;
0081       File_HCAL->Close();
0082       cout << "input energy is " << inputEnergy << ", calibrated energy is " << val_calib << " +/- " << err_calib << ", sigma is " << val_sigma << endl;
0083     }
0084   }
0085
0086   TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
0087   c_Linearity->cd();
0088   c_Linearity->cd()->SetLeftMargin(0.15);
0089   c_Linearity->cd()->SetBottomMargin(0.15);
0090   c_Linearity->cd()->SetTicks(1,1);
0091   c_Linearity->cd()->SetGrid(0,0);
0092
0093   TH1F *h_play = new TH1F("h_play","h_play",100,-0.5,99.5);
0094   for(int i_bin = 0; i_bin < 100; ++i_bin)
0095   {
0096     h_play->SetBinContent(i_bin+1,-100.0);
0097     h_play->SetBinError(i_bin+1,1.0);
0098   }
0099   h_play->SetTitle("");
0100   h_play->SetStats(0);
0101
0102   h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
0103   h_play->GetXaxis()->CenterTitle();
0104   h_play->GetXaxis()->SetRangeUser(0,40);
0105
0106   h_play->GetYaxis()->SetTitle("Calibrated Energy (GeV)");
0107   h_play->GetYaxis()->CenterTitle();
0108   h_play->GetYaxis()->SetRangeUser(0,40);
0109   h_play->DrawCopy("pE");
0110
0111   TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
0112   string inputlinearity = "linearity.txt";
0113   FILE *flinear = fopen(inputlinearity.c_str(), "r");
0114   if(flinear == NULL)
0115   {
0116     perror("Error opening file!");
0117   }
0118   else
0119   {
0120     float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
0121     char line[1024];
0122     int line_counter = 0;
0123     while(fgets(line,1024,flinear))
0124     {
0125       sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
0126       // 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;
0127       g_linearity_old->SetPoint(line_counter,x,y);
0128       g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
0129       line_counter++;
0130     }
0131   }
0132
0133   g_linearity_old->SetMarkerStyle(21);
0134   g_linearity_old->SetMarkerColor(2);
0135   g_linearity_old->SetMarkerSize(1.8);
0136   g_linearity_old->Draw("pE same");
0137
0138   g_linearity->SetMarkerStyle(20);
0139   g_linearity->SetMarkerColor(kGray+2);
0140   g_linearity->SetMarkerSize(2.0);
0141   g_linearity->Draw("pE same");
0142
0143   TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
0144   l_unity->SetLineColor(4);
0145   l_unity->SetLineStyle(2);
0146   l_unity->SetLineWidth(2);
0147   l_unity->Draw("l same");
0148
0149   TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
0150   leg_linear->SetBorderSize(0);
0151   leg_linear->SetFillColor(0);
0152   leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
0153   leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
0154   leg_linear->AddEntry(l_unity,"unity","l");
0155   leg_linear->Draw("same");
0156
0157   c_Linearity->SaveAs("/gpfs/mnt/gpfs04/sphenix/user/xusun/HCAL_BeamTest/figures/c_Linearity.eps");
0158
0159   TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
0160   c_Resolution->cd();
0161   c_Resolution->cd()->SetLeftMargin(0.15);
0162   c_Resolution->cd()->SetBottomMargin(0.15);
0163   c_Resolution->cd()->SetTicks(1,1);
0164   c_Resolution->cd()->SetGrid(0,0);
0165
0166   h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
0167   h_play->GetXaxis()->CenterTitle();
0168   h_play->GetXaxis()->SetRangeUser(0,40);
0169
0170   h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
0171   h_play->GetYaxis()->CenterTitle();
0172   h_play->GetYaxis()->SetRangeUser(0,0.7);
0173   h_play->DrawCopy("pE");
0174
0175   TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
0176   string inputresolution = "resolution.txt";
0177   FILE *fres = fopen(inputresolution.c_str(), "r");
0178   if(fres == NULL)
0179   {
0180     perror("Error opening file!");
0181   }
0182   else
0183   {
0184     float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
0185     char line[1024];
0186     int line_counter = 0;
0187     while(fgets(line,1024,fres))
0188     {
0189       sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
0190       // 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;
0191       g_resolution_old->SetPoint(line_counter,x,y);
0192       g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
0193       line_counter++;
0194     }
0195   }
0196
0197   g_resolution_old->SetMarkerStyle(21);
0198   g_resolution_old->SetMarkerColor(2);
0199   g_resolution_old->SetMarkerSize(1.8);
0200   g_resolution_old->Draw("pE same");
0201
0202   g_resolution->SetMarkerStyle(20);
0203   g_resolution->SetMarkerColor(kGray+2);
0204   g_resolution->SetMarkerSize(2.0);
0205   g_resolution->Draw("pE same");
0206
0207   TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
0208   leg_res->SetBorderSize(0);
0209   leg_res->SetFillColor(0);
0210   leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
0211   leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
0212   leg_res->Draw("same");
0213
0214   c_Resolution->SaveAs("/gpfs/mnt/gpfs04/sphenix/user/xusun/HCAL_BeamTest/figures/c_Resolution.eps");
0215 }