EICMCAnalysis/jets/plot_jets.C

0001 /* ROOT includes */
0002 #include <TROOT.h>
0003 #include <TFile.h>
0004 #include <TTree.h>
0005 #include <TCut.h>
0006 #include <TF1.h>
0007 #include <TH1F.h>
0008 #include <TH2F.h>
0009 #include <TGraphErrors.h>
0010 #include <TProfile.h>
0011 #include <TCanvas.h>
0012 #include <TString.h>
0013
0014 using namespace std;
0015
0016 int
0017 plot_jets(TString jetfilename)
0018 {
0019   TFile *fin = new TFile(jetfilename, "OPEN");
0020
0021   TTree* jets = (TTree*)fin->Get("jets");
0022
0023   /* basic cut that is applied to every plot */
0024   TCut cut_base("1");
0025
0026   /* select different ranges in pseudorapidity */
0027   vector< TCut > v_cuts_eta;
0028   v_cuts_eta.push_back( TCut("(jet_truth_eta>-4&&jet_truth_eta<-2)") );
0029   v_cuts_eta.push_back( TCut("(jet_truth_eta>-1&&jet_truth_eta<1)") );
0030   v_cuts_eta.push_back( TCut("(jet_truth_eta>2&&jet_truth_eta<4)") );
0031
0032   /* select different ranges in energy */
0033   vector< TCut > v_cuts_energy;
0034   v_cuts_energy.push_back( TCut("(jet_truth_e>0&&jet_truth_e<10)") );
0035   v_cuts_energy.push_back( TCut("(jet_truth_e>10&&jet_truth_e<20)") );
0036   v_cuts_energy.push_back( TCut("(jet_truth_e>20&&jet_truth_e<50)") );
0037
0038   /* temporary cancas where objects that need to be drawn end up */
0039   TCanvas *ctemp = new TCanvas("**Scratch**");
0040
0041
0042   /* PLOT: True jet energies vs pseudorapidity */
0043   ctemp->cd();
0044
0045   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);
0046   jets->Draw("jet_truth_e : jet_truth_eta >> h_jets_energy_vs_eta","","colz");
0047
0048   TCanvas *c4 = new TCanvas();
0049   h_jets_energy_vs_eta->Draw("colz");
0050   c4->SetLogz();
0051   c4->Print("plots/plot_truth_jet_energy_vs_Eta.eps");
0052
0053   /* PLOT: Histogram of jet energies in different pseudorapidity ranges */
0054   ctemp->cd();
0055
0056   TH1F* h_jets_energy_eta1 = new TH1F("h_jets_energy_eta1",";E_{jet}^{truth} (GeV);counts",100,0,50);
0057   TH1F* h_jets_energy_eta2 = (TH1F*)h_jets_energy_eta1->Clone("h_jets_energy_eta2");
0058   TH1F* h_jets_energy_eta3 = (TH1F*)h_jets_energy_eta1->Clone("h_jets_energy_eta3");
0059
0060   h_jets_energy_eta1->SetLineColor(kBlue);
0061   h_jets_energy_eta2->SetLineColor(kRed);
0062   h_jets_energy_eta3->SetLineColor(kGreen+1);
0063
0064   jets->Draw("jet_truth_e >> h_jets_energy_eta1", v_cuts_eta.at(0), "");
0065   jets->Draw("jet_truth_e >> h_jets_energy_eta2", v_cuts_eta.at(1), "");
0066   jets->Draw("jet_truth_e >> h_jets_energy_eta3", v_cuts_eta.at(2), "");
0067
0068   TCanvas *c2 = new TCanvas();
0069   h_jets_energy_eta2->Draw("");
0070   h_jets_energy_eta1->Draw("same");
0071   h_jets_energy_eta3->Draw("same");
0072   c2->Print("plots/plot_hist_etruth_etaranges.eps");
0073
0074   /* PLOT: Jet energy resolution vs energy- all pseudorapidities */
0075   ctemp->cd();
0076
0077   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);
0078   jets->Draw("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_e >> h_eres_vs_e", cut_base);
0079
0080   TCanvas *c5 = new TCanvas();
0081   h_eres_vs_e->Draw("colz");
0082   c5->Print("plots/plot_eres_vs_e_all.eps");
0083
0084   /* PLOT: Jet energy resolution vs eta- all energies */
0085   ctemp->cd();
0086
0087   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);
0088   jets->Draw("(jet_smear_e-jet_truth_e)/jet_truth_e:jet_truth_eta >> h_eres_vs_eta", cut_base);
0089
0090   TCanvas *c6 = new TCanvas();
0091   h_eres_vs_eta->Draw("colz");
0092   c6->Print("plots/plot_eres_vs_eta_all.eps");
0093
0094   //TProfile* hprof_eres_vs_eta = h_eres_vs_eta->ProfileX();
0095   //hprof_eres_vs_eta->GetYaxis()->SetTitle( h_eres_vs_eta->GetYaxis()->GetTitle() );
0096   //
0097   //TCanvas *c7 = new TCanvas();
0098   //hprof_eres_vs_eta->Draw();
0099
0100
0101   /* PLOT: Jet energy resolution vs energy in different regions of pseudorapidity */
0102   ctemp->cd();
0103
0104   TH1F* h_frame_emean = new TH1F("h_frame_emean", "", 10,2.5,52.5);
0105   h_frame_emean->GetYaxis()->SetRangeUser(-0.4,0.4);
0106   h_frame_emean->GetXaxis()->SetTitle("E_{jet}^{truth} [GeV]");
0107   h_frame_emean->GetYaxis()->SetTitle("(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}");
0108
0109   TH1F* h_frame_esigma = new TH1F("h_frame_esigma", "", 10,2.5,52.5);
0110   h_frame_esigma->GetYaxis()->SetRangeUser(0.1,0.4);
0111   h_frame_esigma->GetXaxis()->SetTitle("E_{jet}^{truth} [GeV]");
0112   h_frame_esigma->GetYaxis()->SetTitle("#sigma ( (E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth} )");
0113
0114   for ( unsigned etaidx = 0; etaidx < v_cuts_eta.size(); etaidx++ )
0115     {
0116       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);
0117       Int_t nbins_x_etacut = h_eres_etacut->GetXaxis()->GetNbins();
0118
0119       TGraphErrors* g_emean_etacut = new TGraphErrors( nbins_x_etacut );
0120       TGraphErrors* g_esigma_etacut = new TGraphErrors( nbins_x_etacut );
0121
0122       TCanvas *c_evseta_1 = new TCanvas(TString::Format( "etarange_%d", etaidx ), (TString::Format( "Pseudorapidity range: %s", v_cuts_eta.at(etaidx).GetTitle() )));
0123       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) );
0124       h_eres_etacut->Draw("COLZ");
0125       c_evseta_1->Print(TString::Format( "plots/plot_%s.eps", c_evseta_1->GetName() ));
0126
0127       /* Loop over all x-bins in 2D histogram, create Y projections, and get standard deviation */
0128       for ( Int_t i = 0; i < nbins_x_etacut; i++ )
0129         {
0130           ctemp->cd();
0131
0132           TH1D* h_proj = h_eres_etacut->ProjectionY("py",i+1,i+1);
0133
0134           /* skip projections with too few entries */
0135           if ( h_proj->GetEntries() < 100 )
0136             continue;
0137
0138       /* Fit gaussian to projected histogram */
0139           h_proj->Fit("gaus","Q");
0140           TF1* fit = h_proj->GetFunction("gaus");
0141
0142           TCanvas *cstore = new TCanvas();
0143
0144           TString proj_name = TString::Format("Projection_etacut%d_bincenter_%.01fGeV", etaidx, h_eres_etacut->GetXaxis()->GetBinCenter(i+1) );
0145
0146           cstore->SetName(proj_name);
0147           cstore->SetTitle(proj_name);
0148           h_proj->DrawClone("clone");
0149       cstore->Print(TString::Format( "plots/plot_%s.eps", cstore->GetName() ));
0150
0151           cout << "RMS: " << h_eres_etacut->GetXaxis()->GetBinCenter(i+1) << " --> " <<  h_proj->GetRMS() << endl;
0152
0153           if (fit)
0154             {
0155               cout << "FIT: " << h_eres_etacut->GetXaxis()->GetBinCenter(i+1) << " --> " <<  fit->GetParameter(2) << endl;
0156               g_emean_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(1));
0157               g_emean_etacut->SetPointError(i, 0, fit->GetParError(1));
0158
0159               g_esigma_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(2));
0160               g_esigma_etacut->SetPointError(i, 0, fit->GetParError(2));
0161             }
0162
0163           //g_esigma_etacut->SetPoint(i, h_eres_etacut->GetXaxis()->GetBinCenter(i+1), h_proj->GetRMS());
0164         }
0165
0166       /* Draw resolution graph */
0167       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() )));
0168       h_frame_emean->Draw();
0169       g_emean_etacut->Draw("Psame");
0170       c_evseta_2->Print(TString::Format( "plots/plot_%s.eps", c_evseta_2->GetName() ));
0171
0172       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() )));
0173       h_frame_esigma->Draw();
0174       g_esigma_etacut->Draw("Psame");
0175       c_evseta_3->Print(TString::Format( "plots/plot_%s.eps", c_evseta_3->GetName() ));
0176     }
0177
0178   /* PLOT: Jet energy resolution vs pseudorapidity per energy range */
0179     ctemp->cd();
0180
0181   TH1F* h_frame_energycut_emean = new TH1F("h_frame_energycut_emean", "", 32,-4,4);
0182   h_frame_energycut_emean->GetYaxis()->SetRangeUser(-0.4,0.4);
0183   h_frame_energycut_emean->GetXaxis()->SetTitle("#eta_{jet}^{truth}");
0184   h_frame_energycut_emean->GetYaxis()->SetTitle("(E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth}");
0185
0186   TH1F* h_frame_energycut_esigma = new TH1F("h_frame_energycut_esigma", "", 32,-4,4);
0187   h_frame_energycut_esigma->GetYaxis()->SetRangeUser(0.1,0.4);
0188   h_frame_energycut_esigma->GetXaxis()->SetTitle("#eta_{jet}^{truth}");
0189   h_frame_energycut_esigma->GetYaxis()->SetTitle("#sigma ( (E_{jet}^{smear}-E_{jet}^{truth}) / E_{jet}^{truth} )");
0190
0191   for ( unsigned energyidx = 0; energyidx < v_cuts_energy.size(); energyidx++ )
0192     {
0193       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);
0194       Int_t nbins_x_energycut = h_eres_energycut->GetXaxis()->GetNbins();
0195
0196       TGraphErrors* g_emean_energycut = new TGraphErrors( nbins_x_energycut );
0197       TGraphErrors* g_esigma_energycut = new TGraphErrors( nbins_x_energycut );
0198
0199       TCanvas *c_evseta_1 = new TCanvas(TString::Format( "energyrange_%d", energyidx ), (TString::Format( "Energy range: %s", v_cuts_energy.at(energyidx).GetTitle() )));
0200       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) );
0201       h_eres_energycut->Draw("COLZ");
0202       c_evseta_1->Print(TString::Format( "plots/plot_%s.eps", c_evseta_1->GetName() ));
0203
0204       /* Loop over all x-bins in 2D histogram, create Y projections, and get standard deviation */
0205       for ( Int_t i = 0; i < nbins_x_energycut; i++ )
0206         {
0207           ctemp->cd();
0208
0209           TH1D* h_proj = h_eres_energycut->ProjectionY("py",i+1,i+1);
0210
0211           /* skip projections with too few entries */
0212           if ( h_proj->GetEntries() < 100 )
0213             continue;
0214
0215       /* Fit gaussian to projected histogram */
0216           h_proj->Fit("gaus","Q");
0217           TF1* fit = h_proj->GetFunction("gaus");
0218
0219           TCanvas *cstore = new TCanvas();
0220
0221           TString proj_name = TString::Format("Projection_energycut%d_bincenter_%.01fEta", energyidx, h_eres_energycut->GetXaxis()->GetBinCenter(i+1) );
0222
0223           cstore->SetName(proj_name);
0224           cstore->SetTitle(proj_name);
0225           h_proj->DrawClone("clone");
0226       cstore->Print(TString::Format( "plots/plot_%s.eps", cstore->GetName() ));
0227
0228           cout << "RMS: " << h_eres_energycut->GetXaxis()->GetBinCenter(i+1) << " --> " <<  h_proj->GetRMS() << endl;
0229
0230           if (fit)
0231             {
0232               cout << "FIT: " << h_eres_energycut->GetXaxis()->GetBinCenter(i+1) << " --> " <<  fit->GetParameter(2) << endl;
0233               g_emean_energycut->SetPoint(i, h_eres_energycut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(1));
0234               g_emean_energycut->SetPointError(i, 0, fit->GetParError(1));
0235
0236               g_esigma_energycut->SetPoint(i, h_eres_energycut->GetXaxis()->GetBinCenter(i+1), fit->GetParameter(2));
0237               g_esigma_energycut->SetPointError(i, 0, fit->GetParError(2));
0238             }
0239         }
0240
0241       /* Draw resolution graph */
0242       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() )));
0243       h_frame_energycut_emean->Draw();
0244       g_emean_energycut->Draw("Psame");
0245       c_evseta_2->Print(TString::Format( "plots/plot_%s.eps", c_evseta_2->GetName() ));
0246
0247       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() )));
0248       h_frame_energycut_esigma->Draw();
0249       g_esigma_energycut->Draw("Psame");
0250       c_evseta_3->Print(TString::Format( "plots/plot_%s.eps", c_evseta_3->GetName() ));
0251     }
0252
0253   return 0;
0254 }