Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_jets.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_jets.C
1 /* ROOT includes */
2 #include <TROOT.h>
3 #include <TFile.h>
4 #include <TTree.h>
5 #include <TCut.h>
6 #include <TF1.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TGraphErrors.h>
10 #include <TProfile.h>
11 #include <TCanvas.h>
12 #include <TString.h>
13 
14 using namespace std;
15 
16 int
17 plot_jets(TString jetfilename)
18 {
19  TFile *fin = new TFile(jetfilename, "OPEN");
20 
21  TTree* jets = (TTree*)fin->Get("jets");
22 
23  /* basic cut that is applied to every plot */
24  TCut cut_base("1");
25 
26  /* select different ranges in pseudorapidity */
27  vector< TCut > v_cuts_eta;
28  v_cuts_eta.push_back( TCut("(jet_truth_eta>-4&&jet_truth_eta<-2)") );
29  v_cuts_eta.push_back( TCut("(jet_truth_eta>-1&&jet_truth_eta<1)") );
30  v_cuts_eta.push_back( TCut("(jet_truth_eta>2&&jet_truth_eta<4)") );
31 
32  /* select different ranges in energy */
33  vector< TCut > v_cuts_energy;
34  v_cuts_energy.push_back( TCut("(jet_truth_e>0&&jet_truth_e<10)") );
35  v_cuts_energy.push_back( TCut("(jet_truth_e>10&&jet_truth_e<20)") );
36  v_cuts_energy.push_back( TCut("(jet_truth_e>20&&jet_truth_e<50)") );
37 
38  /* temporary cancas where objects that need to be drawn end up */
39  TCanvas *ctemp = new TCanvas("**Scratch**");
40 
41 
42  /* PLOT: True jet energies vs pseudorapidity */
43  ctemp->cd();
44 
45  TH2F* h_jets_energy_vs_eta = new TH2F("h_jets_energy_vs_eta",";#eta_{jet}^{truth};E_{jet}^{truth}",36,-4.5,4.5,25,0,100);
46  jets->Draw("jet_truth_e : jet_truth_eta >> h_jets_energy_vs_eta","","colz");
47 
48  TCanvas *c4 = new TCanvas();
49  h_jets_energy_vs_eta->Draw("colz");
50  c4->SetLogz();
51  c4->Print("plots/plot_truth_jet_energy_vs_Eta.eps");
52 
53  /* PLOT: Histogram of jet energies in different pseudorapidity ranges */
54  ctemp->cd();
55 
56  TH1F* h_jets_energy_eta1 = new TH1F("h_jets_energy_eta1",";E_{jet}^{truth} (GeV);counts",100,0,50);
57  TH1F* h_jets_energy_eta2 = (TH1F*)h_jets_energy_eta1->Clone("h_jets_energy_eta2");
58  TH1F* h_jets_energy_eta3 = (TH1F*)h_jets_energy_eta1->Clone("h_jets_energy_eta3");
59 
60  h_jets_energy_eta1->SetLineColor(kBlue);
61  h_jets_energy_eta2->SetLineColor(kRed);
62  h_jets_energy_eta3->SetLineColor(kGreen+1);
63 
64  jets->Draw("jet_truth_e >> h_jets_energy_eta1", v_cuts_eta.at(0), "");
65  jets->Draw("jet_truth_e >> h_jets_energy_eta2", v_cuts_eta.at(1), "");
66  jets->Draw("jet_truth_e >> h_jets_energy_eta3", v_cuts_eta.at(2), "");
67 
68  TCanvas *c2 = new TCanvas();
69  h_jets_energy_eta2->Draw("");
70  h_jets_energy_eta1->Draw("same");
71  h_jets_energy_eta3->Draw("same");
72  c2->Print("plots/plot_hist_etruth_etaranges.eps");
73 
74  /* PLOT: Jet energy resolution vs energy- all pseudorapidities */
75  ctemp->cd();
76 
77  TH2F* h_eres_vs_e = new TH2F("h_eres_vs_e",";E_{jet}^{truth};(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}",2*12,2.5,62.5,2*20,-2,4);
78  jets->Draw("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_e >> h_eres_vs_e", cut_base);
79 
80  TCanvas *c5 = new TCanvas();
81  h_eres_vs_e->Draw("colz");
82  c5->Print("plots/plot_eres_vs_e_all.eps");
83 
84  /* PLOT: Jet energy resolution vs eta- all energies */
85  ctemp->cd();
86 
87  TH2F* h_eres_vs_eta = new TH2F("h_eres_vs_eta",";#eta_{jet}^{truth};(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}",2*18,-4.5,4.5,2*20,-2,4);
88  jets->Draw("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_eta >> h_eres_vs_eta", cut_base);
89 
90  TCanvas *c6 = new TCanvas();
91  h_eres_vs_eta->Draw("colz");
92  c6->Print("plots/plot_eres_vs_eta_all.eps");
93 
94  //TProfile* hprof_eres_vs_eta = h_eres_vs_eta->ProfileX();
95  //hprof_eres_vs_eta->GetYaxis()->SetTitle( h_eres_vs_eta->GetYaxis()->GetTitle() );
96  //
97  //TCanvas *c7 = new TCanvas();
98  //hprof_eres_vs_eta->Draw();
99 
100 
101  /* PLOT: Jet energy resolution vs energy in different regions of pseudorapidity */
102  ctemp->cd();
103 
104  TH1F* h_frame_emean = new TH1F("h_frame_emean", "", 10,2.5,52.5);
105  h_frame_emean->GetYaxis()->SetRangeUser(-0.4,0.4);
106  h_frame_emean->GetXaxis()->SetTitle("E_{jet}^{truth} [GeV]");
107  h_frame_emean->GetYaxis()->SetTitle("(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}");
108 
109  TH1F* h_frame_esigma = new TH1F("h_frame_esigma", "", 10,2.5,52.5);
110  h_frame_esigma->GetYaxis()->SetRangeUser(0.1,0.4);
111  h_frame_esigma->GetXaxis()->SetTitle("E_{jet}^{truth} [GeV]");
112  h_frame_esigma->GetYaxis()->SetTitle("#sigma ( (E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth} )");
113 
114  for ( unsigned etaidx = 0; etaidx < v_cuts_eta.size(); etaidx++ )
115  {
116  TH2F* h_eres_etacut = new TH2F(TString::Format("h_eres_etacut%d",etaidx), ";E_{jet}^{truth};(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}",10,2.5,52.5,20,-1,1);
117  Int_t nbins_x_etacut = h_eres_etacut->GetXaxis()->GetNbins();
118 
119  TGraphErrors* g_emean_etacut = new TGraphErrors( nbins_x_etacut );
120  TGraphErrors* g_esigma_etacut = new TGraphErrors( nbins_x_etacut );
121 
122  TCanvas *c_evseta_1 = new TCanvas(TString::Format( "etarange_%d", etaidx ), (TString::Format( "Pseudorapidity range: %s", v_cuts_eta.at(etaidx).GetTitle() )));
123  jets->Draw(TString::Format("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_e >> h_eres_etacut%d",etaidx), cut_base && v_cuts_eta.at(etaidx) );
124  h_eres_etacut->Draw("COLZ");
125  c_evseta_1->Print(TString::Format( "plots/plot_%s.eps", c_evseta_1->GetName() ));
126 
127  /* Loop over all x-bins in 2D histogram, create Y projections, and get standard deviation */
128  for ( Int_t i = 0; i < nbins_x_etacut; i++ )
129  {
130  ctemp->cd();
131 
132  TH1D* h_proj = h_eres_etacut->ProjectionY("py",i+1,i+1);
133 
134  /* skip projections with too few entries */
135  if ( h_proj->GetEntries() < 100 )
136  continue;
137 
138  /* Fit gaussian to projected histogram */
139  h_proj->Fit("gaus","Q");
140  TF1* fit = h_proj->GetFunction("gaus");
141 
142  TCanvas *cstore = new TCanvas();
143 
144  TString proj_name = TString::Format("Projection_etacut%d_bincenter_%.01fGeV", etaidx, h_eres_etacut->GetXaxis()->GetBinCenter(i+1) );
145 
146  cstore->SetName(proj_name);
147  cstore->SetTitle(proj_name);
148  h_proj->DrawClone("clone");
149  cstore->Print(TString::Format( "plots/plot_%s.eps", cstore->GetName() ));
150 
151  cout << "RMS: " << h_eres_etacut->GetXaxis()->GetBinCenter(i+1) << " --> " << h_proj->GetRMS() << endl;
152 
153  if (fit)
154  {
155  cout << "FIT: " << h_eres_etacut->GetXaxis()->GetBinCenter(i+1) << " --> " << fit->GetParameter(2) << endl;
156  g_emean_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(1));
157  g_emean_etacut->SetPointError(i, 0, fit->GetParError(1));
158 
159  g_esigma_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(2));
160  g_esigma_etacut->SetPointError(i, 0, fit->GetParError(2));
161  }
162 
163  //g_esigma_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), h_proj->GetRMS());
164  }
165 
166  /* Draw resolution graph */
167  TCanvas *c_evseta_2 = new TCanvas(TString::Format( "mean_vs_e_etarange%d", etaidx ), (TString::Format( "Pseudorapidity range: %s", v_cuts_eta.at(etaidx).GetTitle() )));
168  h_frame_emean->Draw();
169  g_emean_etacut->Draw("Psame");
170  c_evseta_2->Print(TString::Format( "plots/plot_%s.eps", c_evseta_2->GetName() ));
171 
172  TCanvas *c_evseta_3 = new TCanvas(TString::Format( "sigma_vs_e_etarange%d", etaidx), (TString::Format( "Pseudorapidity range: %s", v_cuts_eta.at(etaidx).GetTitle() )));
173  h_frame_esigma->Draw();
174  g_esigma_etacut->Draw("Psame");
175  c_evseta_3->Print(TString::Format( "plots/plot_%s.eps", c_evseta_3->GetName() ));
176  }
177 
178  /* PLOT: Jet energy resolution vs pseudorapidity per energy range */
179  ctemp->cd();
180 
181  TH1F* h_frame_energycut_emean = new TH1F("h_frame_energycut_emean", "", 32,-4,4);
182  h_frame_energycut_emean->GetYaxis()->SetRangeUser(-0.4,0.4);
183  h_frame_energycut_emean->GetXaxis()->SetTitle("#eta_{jet}^{truth}");
184  h_frame_energycut_emean->GetYaxis()->SetTitle("(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}");
185 
186  TH1F* h_frame_energycut_esigma = new TH1F("h_frame_energycut_esigma", "", 32,-4,4);
187  h_frame_energycut_esigma->GetYaxis()->SetRangeUser(0.1,0.4);
188  h_frame_energycut_esigma->GetXaxis()->SetTitle("#eta_{jet}^{truth}");
189  h_frame_energycut_esigma->GetYaxis()->SetTitle("#sigma ( (E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth} )");
190 
191  for ( unsigned energyidx = 0; energyidx < v_cuts_energy.size(); energyidx++ )
192  {
193  TH2F* h_eres_energycut = new TH2F(TString::Format("h_eres_energycut%d",energyidx), ";E_{jet}^{truth};(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}",32,-4,4,20,-1,1);
194  Int_t nbins_x_energycut = h_eres_energycut->GetXaxis()->GetNbins();
195 
196  TGraphErrors* g_emean_energycut = new TGraphErrors( nbins_x_energycut );
197  TGraphErrors* g_esigma_energycut = new TGraphErrors( nbins_x_energycut );
198 
199  TCanvas *c_evseta_1 = new TCanvas(TString::Format( "energyrange_%d", energyidx ), (TString::Format( "Energy range: %s", v_cuts_energy.at(energyidx).GetTitle() )));
200  jets->Draw(TString::Format("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_eta >> h_eres_energycut%d",energyidx), cut_base && v_cuts_energy.at(energyidx) );
201  h_eres_energycut->Draw("COLZ");
202  c_evseta_1->Print(TString::Format( "plots/plot_%s.eps", c_evseta_1->GetName() ));
203 
204  /* Loop over all x-bins in 2D histogram, create Y projections, and get standard deviation */
205  for ( Int_t i = 0; i < nbins_x_energycut; i++ )
206  {
207  ctemp->cd();
208 
209  TH1D* h_proj = h_eres_energycut->ProjectionY("py",i+1,i+1);
210 
211  /* skip projections with too few entries */
212  if ( h_proj->GetEntries() < 100 )
213  continue;
214 
215  /* Fit gaussian to projected histogram */
216  h_proj->Fit("gaus","Q");
217  TF1* fit = h_proj->GetFunction("gaus");
218 
219  TCanvas *cstore = new TCanvas();
220 
221  TString proj_name = TString::Format("Projection_energycut%d_bincenter_%.01fEta", energyidx, h_eres_energycut->GetXaxis()->GetBinCenter(i+1) );
222 
223  cstore->SetName(proj_name);
224  cstore->SetTitle(proj_name);
225  h_proj->DrawClone("clone");
226  cstore->Print(TString::Format( "plots/plot_%s.eps", cstore->GetName() ));
227 
228  cout << "RMS: " << h_eres_energycut->GetXaxis()->GetBinCenter(i+1) << " --> " << h_proj->GetRMS() << endl;
229 
230  if (fit)
231  {
232  cout << "FIT: " << h_eres_energycut->GetXaxis()->GetBinCenter(i+1) << " --> " << fit->GetParameter(2) << endl;
233  g_emean_energycut->SetPoint(i, h_eres_energycut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(1));
234  g_emean_energycut->SetPointError(i, 0, fit->GetParError(1));
235 
236  g_esigma_energycut->SetPoint(i, h_eres_energycut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(2));
237  g_esigma_energycut->SetPointError(i, 0, fit->GetParError(2));
238  }
239  }
240 
241  /* Draw resolution graph */
242  TCanvas *c_evseta_2 = new TCanvas(TString::Format( "mean_vs_e_energyrange%d", energyidx ), (TString::Format( "Energy range: %s", v_cuts_energy.at(energyidx).GetTitle() )));
243  h_frame_energycut_emean->Draw();
244  g_emean_energycut->Draw("Psame");
245  c_evseta_2->Print(TString::Format( "plots/plot_%s.eps", c_evseta_2->GetName() ));
246 
247  TCanvas *c_evseta_3 = new TCanvas(TString::Format( "sigma_vs_e_energyrange%d", energyidx), (TString::Format( "Energy range: %s", v_cuts_energy.at(energyidx).GetTitle() )));
248  h_frame_energycut_esigma->Draw();
249  g_esigma_energycut->Draw("Psame");
250  c_evseta_3->Print(TString::Format( "plots/plot_%s.eps", c_evseta_3->GetName() ));
251  }
252 
253  return 0;
254 }