Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Jet_reso.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Jet_reso.C
1 #include "sPhenixStyle.h"
2 #include "sPhenixStyle.C"
3 
4 
5 void Jet_reso()
6 {
7  SetsPhenixStyle();
8  TH1::SetDefaultSumw2();
9  TH2::SetDefaultSumw2();
10 
11  TChain * ct = new TChain("T");
12  //if you want to run one file use this
13  ct->Add("../macro/output.root");
14 
15  //if you want to combine multiple files use this
16  /*for(int i = 0; i < 100; i++){
17  ct->Add(Form("../macro/condor/output_%i.root",i));
18  }*/
19 
20  vector<float> *eta = 0;
21  vector<float> *phi = 0;
22  vector<float> *pt = 0;
23  vector<float> *e = 0;
24  vector<float> *subtracted_et = 0;
25  vector<float> *truthEta = 0;
26  vector<float> *truthPhi = 0;
27  vector<float> *truthPt = 0;
28  vector<float> *truthE = 0;
29  vector<float> *subseedEta = 0;
30  vector<float> *subseedPhi = 0;
31  vector<float> *subseedPt = 0;
32  vector<float> *subseedE = 0;
33  vector<int> *subseedCut = 0;
34  vector<float> *rawseedEta = 0;
35  vector<float> *rawseedPhi = 0;
36  vector<float> *rawseedPt = 0;
37  vector<float> *rawseedE = 0;
38  vector<int> *rawseedCut = 0;
39  float *totalCalo = 0;
40  int cent;
41 
42  ct->SetBranchAddress("eta",&eta);
43  ct->SetBranchAddress("phi",&phi);
44  ct->SetBranchAddress("pt",&pt);
45  ct->SetBranchAddress("e",&e);
46  ct->SetBranchAddress("subtracted_et",&subtracted_et);
47  ct->SetBranchAddress("truthEta",&truthEta);
48  ct->SetBranchAddress("truthPhi",&truthPhi);
49  ct->SetBranchAddress("truthPt",&truthPt);
50  ct->SetBranchAddress("truthE",&truthE);
51  ct->SetBranchAddress("cent", &cent);
52 
53  ct->SetBranchAddress("rawseedEta", &rawseedEta);
54  ct->SetBranchAddress("rawseedPhi", &rawseedPhi);
55  ct->SetBranchAddress("rawseedPt", &rawseedPt);
56  ct->SetBranchAddress("rawseedE", &rawseedE);
57  ct->SetBranchAddress("rawseedCut", &rawseedCut);
58  ct->SetBranchAddress("subseedEta", &subseedEta);
59  ct->SetBranchAddress("subseedPhi", &subseedPhi);
60  ct->SetBranchAddress("subseedPt", &subseedPt);
61  ct->SetBranchAddress("subseedE", &subseedE);
62  ct->SetBranchAddress("subseedCut", &subseedCut);
63 
64  const Float_t pt_range[] = {10,15,20,25,30,35,40,45,50,60,80};
65  const int pt_N = sizeof(pt_range)/sizeof(float) - 1;
66  const int resp_N = 60;
67  Float_t resp_bins[resp_N+1];
68  for(int i = 0; i < resp_N+1; i++){
69  resp_bins[i] = 1.2/resp_N * i;
70  }
71  const int eta_N = 40;
72  Float_t eta_bins[eta_N+1];
73  for(int i = 0; i < eta_N+1; i++){
74  eta_bins[i] = -1.1 + 2.2/eta_N * i;
75  }
76  const int phi_N = 40;
77  Float_t phi_bins[phi_N+1];
78  for(int i = 0; i < phi_N+1; i++){
79  phi_bins[i] = -TMath::Pi() + 2*TMath::Pi()/phi_N * i;
80  }
81  const int subet_N = 400;
82  Float_t subet_bins[subet_N+1];
83  for(int i = 0; i < subet_N+1; i++){
84  subet_bins[i] = i*0.5;
85  }
86  const float cent_bins[] = {-1, 0, 10, 30, 50, 80}; //the first bin is for inclusive in centrality/pp events
87  const int cent_N = sizeof(cent_bins)/sizeof(float) - 1;
88 
89  TH3F *h_MC_Reso_pt = new TH3F("h_MC_Reso", "" , pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);;
90  h_MC_Reso_pt->GetXaxis()->SetTitle("p_{T}^{truth} [GeV]");
91  h_MC_Reso_pt->GetYaxis()->SetTitle("p_{T}^{reco}/p_{T}^{truth}");
92 
93  TH2F *h_pt_reco = new TH2F("h_pt_reco","",subet_N,subet_bins, cent_N, cent_bins);
94  h_pt_reco->GetXaxis()->SetTitle("p_{T} [GeV]");
95  TH2F *h_pt_true = new TH2F("h_pt_true","",pt_N,pt_range, cent_N, cent_bins);
96  h_pt_true->GetXaxis()->SetTitle("p_{T} [GeV]");
97  TH2F *h_pt_true_matched = new TH2F("h_pt_true_matched","",subet_N,subet_bins, cent_N, cent_bins);
98  h_pt_true_matched->GetXaxis()->SetTitle("p_{T} [GeV]");
99  TH3F *h_eta_phi_reco = new TH3F("h_eta_phi_reco","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
100  h_eta_phi_reco->GetXaxis()->SetTitle("#eta");
101  h_eta_phi_reco->GetYaxis()->SetTitle("#phi");
102  TH3F *h_eta_phi_true = new TH3F("h_eta_phi_true","",eta_N, eta_bins, phi_N, phi_bins, cent_N, cent_bins);
103  TH3F *h_subEt_pt = new TH3F("h_subEt_pt","",pt_N,pt_range, subet_N, subet_bins, cent_N, cent_bins);
104 
105 
106  TGraphErrors *g_jes[cent_N];
107  TGraphErrors *g_jer[cent_N];
108  for(int icent = 0; icent < cent_N; icent++){
109  g_jes[icent] = new TGraphErrors(pt_N);
110  g_jer[icent] = new TGraphErrors(pt_N);
111  }
112 
113  int nentries = ct->GetEntries();
114  for(int i = 0; i < nentries; i++){
115  ct->GetEntry(i);
116 
117  int njets = truthPt->size();
118  int nrecojets = pt->size();
119 
120  for(int tj = 0; tj < njets; tj++){
121  //fill truth hists
122  h_pt_true->Fill(truthPt->at(tj), cent);
123  h_pt_true->Fill(truthPt->at(tj), -1);
124  h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), cent);
125  h_eta_phi_true->Fill(truthEta->at(tj),truthPhi->at(tj), -1);
126  //do reco to truth jet matching
127  float matchEta, matchPhi, matchPt, matchE, matchsubtracted_et, dR;
128  float dRMax = 100;
129  for(int rj = 0; rj < nrecojets; rj++){
130  float dEta = truthEta->at(tj) - eta->at(rj);
131  float dPhi = truthPhi->at(tj) - phi->at(rj);
132  while(dPhi > TMath::Pi()) dPhi -= 2*TMath::Pi();
133  while(dPhi < -TMath::Pi()) dPhi += 2*TMath::Pi();
134  dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
135  if(dR < dRMax){
136  matchEta = eta->at(rj);
137  matchPhi = phi->at(rj);
138  matchPt = pt->at(rj);
139  matchE = e->at(rj);
140  matchsubtracted_et = subtracted_et->at(rj);
141  dRMax = dR;
142  }
143  }
144 
145  if(dRMax > 0.3) continue;
146  h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), cent);
147  h_pt_true_matched->Fill(truthPt->at(tj), cent);
148  h_eta_phi_reco->Fill(matchEta,matchPhi, cent);
149  h_pt_reco->Fill(matchPt, cent);
150  h_MC_Reso_pt->Fill(truthPt->at(tj),matchPt/truthPt->at(tj), -1);
151  h_pt_true_matched->Fill(truthPt->at(tj), -1);
152  h_eta_phi_reco->Fill(matchEta,matchPhi, -1);
153  h_pt_reco->Fill(matchPt, -1);
154  h_subEt_pt->Fill(matchPt, matchsubtracted_et, cent);
155  h_subEt_pt->Fill(matchPt, matchsubtracted_et, -1);
156  }
157  }
158 
159  TCanvas *c = new TCanvas("c","c");
160  //c->Print("fits04.pdf("); //uncomment these to get a pdf with all the fits
161  TLegend *leg = new TLegend(.25,.2,.6,.5);
162  leg->SetFillStyle(0);
163  gStyle->SetOptFit(1);
164  for(int icent = 0; icent < cent_N; ++icent){
165  for (int i = 0; i < pt_N; ++i){
166  TF1 *func = new TF1("func","gaus",0,1.2);
167  h_MC_Reso_pt->GetXaxis()->SetRange(i+1,i+1);
168  h_MC_Reso_pt->GetZaxis()->SetRange(icent+1,icent+1);
169  TH1F *h_temp = (TH1F*) h_MC_Reso_pt->Project3D("y");
170  h_temp->Fit(func,"","",0,1.2);
171  h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
172  func->SetLineColor(kRed);
173  h_temp->Draw();
174  leg->AddEntry("",Form("%2.0f%%-%2.0f%%",h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+1),h_MC_Reso_pt->GetZaxis()->GetBinLowEdge(icent+2)),"");
175  leg->AddEntry("",Form("%2.0f < p_T < %2.0f GeV",h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+1),h_MC_Reso_pt->GetXaxis()->GetBinLowEdge(i+2)),"");
176  leg->Draw("SAME");
177  /*-------for calculating the JER uncertainty----*/
178  float dsigma = func -> GetParError(2);
179  float denergy = func -> GetParError(1);
180  float sigma = func -> GetParameter(2);
181  float energy = func -> GetParameter(1);
182 
183  float djer = dsigma/energy + sigma*denergy/pow(energy,2);//correct way to calculate jer
184  //accounting for the fact that the
185  //mean response and width are correlated
186  //c->Print("fits04.pdf");
187  leg->Clear();
188  g_jes[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(1));
189  g_jes[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(1));
190  g_jer[icent]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(2)/func->GetParameter(1));
191 
192  //g_jer[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(2));
193  g_jer[icent]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),djer);
194  }
195  }
196  //c->Print("fits04.pdf)");
197 
198 
199 
200  TFile *f_out = new TFile("hists.root","RECREATE");
201  for(int icent = 0; icent < cent_N; ++icent){
202  g_jes[icent]->Write(Form("g_jes_cent%i",icent));
203  g_jer[icent]->Write(Form("g_jer_cent%i",icent));
204  }
205 
206  h_pt_true->Write();
207  h_eta_phi_true->Write();
208  h_pt_true_matched->Write();
209  h_eta_phi_reco->Write();
210  h_pt_reco->Write();
211  h_subEt_pt->Write();
212 }
213