Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetReso.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetReso.C
1 #include <iostream>
2 #include <vector>
3 
4 using namespace std;
5 
6 
7 Int_t JetReso(TString infile="matched.root", TString outfile="resohists.root")
8 {
9 
10  //=======================================================
11  //================ Input File ===========================
12  //=======================================================
13  std::cout << "Input file: " << infile.Data() << std::endl;
14 
15  // open input file
16  TFile fin(infile.Data(), "READ");
17  if(!fin.IsOpen())
18  {
19  std::cout << "Error: could not open input file" << std::endl;
20  exit(-1);
21  }
22 
23  // get tree
24  TTree *tree = (TTree*)fin.Get("T");
25  if(!tree)
26  {
27  std::cout << "Error: could not find tree" << std::endl;
28  exit(-1);
29  }
30 
31  // get branches
32  int centrality = 0;
33  double rho_area = 0;
34  double rho_mult = 0;
35  float weight = 0;
36  std::vector<float> *matched_pt_iter = 0;
37  std::vector<float> *matched_pt_area = 0;
38  std::vector<float> *matched_pt_mult = 0;
39  std::vector<float> *matched_pt_iter_unsub = 0;
40  std::vector<float> *matched_pt_area_unsub = 0;
41  std::vector<float> *matched_pt_mult_unsub = 0;
42  std::vector<float> *matched_truth_pt_iter = 0;
43  std::vector<float> *matched_truth_pt_area = 0;
44  std::vector<float> *matched_truth_pt_mult = 0;
45  std::vector<float> *unmatched_pt_iter = 0;
46  std::vector<float> *unmatched_pt_area = 0;
47  std::vector<float> *unmatched_pt_mult = 0;
48  std::vector<float> *unmatched_pt_iter_unsub = 0;
49  std::vector<float> *unmatched_pt_area_unsub = 0;
50  std::vector<float> *unmatched_pt_mult_unsub = 0;
51  std::vector<float> *missed_truth_pt_iter = 0;
52  std::vector<float> *missed_truth_pt_area = 0;
53  std::vector<float> *missed_truth_pt_mult = 0;
54 
55  tree->SetBranchAddress("weight", &weight);
56  tree->SetBranchAddress("centrality", &centrality);
57  tree->SetBranchAddress("rho_area", &rho_area);
58  tree->SetBranchAddress("rho_mult", &rho_mult);
59  tree->SetBranchAddress("matched_pt_iter", &matched_pt_iter);
60  tree->SetBranchAddress("matched_pt_area", &matched_pt_area);
61  tree->SetBranchAddress("matched_pt_mult", &matched_pt_mult);
62  tree->SetBranchAddress("matched_pt_iter_unsub", &matched_pt_iter_unsub);
63  tree->SetBranchAddress("matched_pt_area_unsub", &matched_pt_area_unsub);
64  tree->SetBranchAddress("matched_pt_mult_unsub", &matched_pt_mult_unsub);
65  tree->SetBranchAddress("matched_truth_pt_iter", &matched_truth_pt_iter);
66  tree->SetBranchAddress("matched_truth_pt_area", &matched_truth_pt_area);
67  tree->SetBranchAddress("matched_truth_pt_mult", &matched_truth_pt_mult);
68  tree->SetBranchAddress("unmatched_pt_iter", &unmatched_pt_iter);
69  tree->SetBranchAddress("unmatched_pt_area", &unmatched_pt_area);
70  tree->SetBranchAddress("unmatched_pt_mult", &unmatched_pt_mult);
71  tree->SetBranchAddress("unmatched_pt_iter_unsub", &unmatched_pt_iter_unsub);
72  tree->SetBranchAddress("unmatched_pt_area_unsub", &unmatched_pt_area_unsub);
73  tree->SetBranchAddress("unmatched_pt_mult_unsub", &unmatched_pt_mult_unsub);
74  tree->SetBranchAddress("missed_truth_pt_iter", &missed_truth_pt_iter);
75  tree->SetBranchAddress("missed_truth_pt_area", &missed_truth_pt_area);
76  tree->SetBranchAddress("missed_truth_pt_mult", &missed_truth_pt_mult);
77 
78  //=======================================================
79  // declare histograms
80  //=======================================================
81  const double pt_range[] = {10,15,20,25,30,35,40,45,50,60,80};
82  const int pt_N = sizeof(pt_range)/sizeof(double) - 1;
83 
84  const int resp_N = 100;
85  double resp_bins[resp_N+1];
86  for(int i = 0; i < resp_N+1; i++){
87  resp_bins[i] = 2.0/resp_N * i;
88  }
89 
90  const double cent_bins[] = {-1, 0, 10, 30, 50, 80}; //the first bin is for inclusive in centrality/pp events
91  const int cent_N = sizeof(cent_bins)/sizeof(double) - 1;
92 
93  TH3D * h3_reso_area = new TH3D("h3_reso_area", "h3_reso_area", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
94  TH3D * h3_reso_mult = new TH3D("h3_reso_mult", "h3_reso_mult", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
95  TH3D * h3_reso_iter = new TH3D("h3_reso_iter", "h3_reso_iter", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
96 
97  // fill histograms
98  int nentries = tree->GetEntries();
99  for (int iEvent = 0; iEvent < nentries; iEvent++){
100  // get entry
101  tree->GetEntry(iEvent);
102 
103  for (unsigned int ijet =0; ijet < matched_pt_iter->size(); ijet++){
104  h3_reso_iter->Fill(matched_truth_pt_iter->at(ijet), matched_pt_iter->at(ijet)/matched_truth_pt_iter->at(ijet), centrality, weight);
105  h3_reso_iter->Fill(matched_truth_pt_iter->at(ijet), matched_pt_iter->at(ijet)/matched_truth_pt_iter->at(ijet), -1, weight);
106  }
107  for (unsigned int ijet =0; ijet < matched_pt_area->size(); ijet++){
108  h3_reso_area->Fill(matched_truth_pt_area->at(ijet), matched_pt_area->at(ijet)/matched_truth_pt_area->at(ijet), centrality, weight);
109  h3_reso_area->Fill(matched_truth_pt_area->at(ijet), matched_pt_area->at(ijet)/matched_truth_pt_area->at(ijet), -1, weight);
110  }
111  for (unsigned int ijet =0; ijet < matched_pt_mult->size(); ijet++){
112  h3_reso_mult->Fill(matched_truth_pt_mult->at(ijet), matched_pt_mult->at(ijet)/matched_truth_pt_mult->at(ijet), centrality, weight);
113  h3_reso_mult->Fill(matched_truth_pt_mult->at(ijet), matched_pt_mult->at(ijet)/matched_truth_pt_mult->at(ijet), -1, weight);
114  }
115  }
116 
117  //=======================================================
118  // JES and JER
119  //=======================================================
120  TH1D * h_jes_area[cent_N];
121  TH1D * h_jer_area[cent_N];
122  TH1D * h_jes_mult[cent_N];
123  TH1D * h_jer_mult[cent_N];
124  TH1D * h_jes_iter[cent_N];
125  TH1D * h_jer_iter[cent_N];
126 
127  for (int icent = 1; icent <= cent_N; icent++){
128  h_jes_area[icent-1] = new TH1D(Form("h_jes_area_%d", icent), Form("h_jes_area_%d", icent), pt_N, pt_range);
129  h_jer_area[icent-1] = new TH1D(Form("h_jer_area_%d", icent), Form("h_jer_area_%d", icent), pt_N, pt_range);
130  h_jes_mult[icent-1] = new TH1D(Form("h_jes_mult_%d", icent), Form("h_jes_mult_%d", icent), pt_N, pt_range);
131  h_jer_mult[icent-1] = new TH1D(Form("h_jer_mult_%d", icent), Form("h_jer_mult_%d", icent), pt_N, pt_range);
132  h_jes_iter[icent-1] = new TH1D(Form("h_jes_iter_%d", icent), Form("h_jes_iter_%d", icent), pt_N, pt_range);
133  h_jer_iter[icent-1] = new TH1D(Form("h_jer_iter_%d", icent), Form("h_jer_iter_%d", icent), pt_N, pt_range);
134  }
135 
136  for(int icent = 0; icent < cent_N; ++icent)
137  {
138  // area
139  for (int i = 0; i < pt_N; ++i)
140  {
141  // fit function
142  TF1 *func = new TF1("func","gaus",0, 1.2);
143  h3_reso_area->GetXaxis()->SetRange(i+1,i+1);
144  h3_reso_area->GetZaxis()->SetRange(icent+1,icent+1);
145  TH1D * h_temp = (TH1D*)h3_reso_area->Project3D("y");
146  h_temp->SetName(Form("h_jes_area_%d_%d", icent, i));
147  h_temp->Fit(func,"","",0,1.2);
148  h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
149  float dsigma = func -> GetParError(2);
150  float denergy = func -> GetParError(1);
151  float sigma = func -> GetParameter(2);
152  float energy = func -> GetParameter(1);
153  float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
154  h_jes_area[icent]->SetBinContent(i+1, energy);
155  h_jes_area[icent]->SetBinError(i+1, denergy);
156  h_jer_area[icent]->SetBinContent(i+1, sigma/energy);
157  h_jer_area[icent]->SetBinError(i+1, djer);
158  }
159 
160  // mult
161  for (int i = 0; i < pt_N; ++i)
162  {
163  // fit function
164  TF1 *func = new TF1("func","gaus",0, 1.2);
165  h3_reso_mult->GetXaxis()->SetRange(i+1,i+1);
166  h3_reso_mult->GetZaxis()->SetRange(icent+1,icent+1);
167  TH1D * h_temp = (TH1D*)h3_reso_mult->Project3D("y");
168  h_temp->SetName(Form("h_jes_mult_%d_%d", icent, i));
169  h_temp->Fit(func,"","",0,1.2);
170  h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
171  float dsigma = func -> GetParError(2);
172  float denergy = func -> GetParError(1);
173  float sigma = func -> GetParameter(2);
174  float energy = func -> GetParameter(1);
175  float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
176  h_jes_mult[icent]->SetBinContent(i+1, energy);
177  h_jes_mult[icent]->SetBinError(i+1, denergy);
178  h_jer_mult[icent]->SetBinContent(i+1, sigma/energy);
179  h_jer_mult[icent]->SetBinError(i+1, djer);
180  }
181 
182  // iter
183  for (int i = 0; i < pt_N; ++i)
184  {
185  // fit function
186  TF1 *func = new TF1("func","gaus",0, 1.2);
187  h3_reso_iter->GetXaxis()->SetRange(i+1,i+1);
188  h3_reso_iter->GetZaxis()->SetRange(icent+1,icent+1);
189  TH1D * h_temp = (TH1D*)h3_reso_iter->Project3D("y");
190  h_temp->SetName(Form("h_jes_iter_%d_%d", icent, i));
191  h_temp->Fit(func,"","",0,1.2);
192  h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
193  float dsigma = func -> GetParError(2);
194  float denergy = func -> GetParError(1);
195  float sigma = func -> GetParameter(2);
196  float energy = func -> GetParameter(1);
197  float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
198  h_jes_iter[icent]->SetBinContent(i+1, energy);
199  h_jes_iter[icent]->SetBinError(i+1, denergy);
200  h_jer_iter[icent]->SetBinContent(i+1, sigma/energy);
201  h_jer_iter[icent]->SetBinError(i+1, djer);
202  }
203 
204 
205  }
206 
207  //=======================================================
208  //================ Output File ==========================
209  //=======================================================
210  std::cout << "Output file: " << outfile.Data() << std::endl;
211  TFile * fout = new TFile(outfile.Data(), "RECREATE");
212  h3_reso_area->Write();
213  h3_reso_mult->Write();
214  h3_reso_iter->Write();
215  for (int icent = 1; icent <= cent_N; icent++){
216  h_jes_area[icent-1]->Write();
217  h_jer_area[icent-1]->Write();
218  h_jes_mult[icent-1]->Write();
219  h_jer_mult[icent-1]->Write();
220  h_jes_iter[icent-1]->Write();
221  h_jer_iter[icent-1]->Write();
222  }
223  fout->Close();
224 
225  // close input file
226  fin.Close();
227  return 0;
228 
229 
230 
231 }