JetBkgdSub/offline/JetReso.C

0001 #include <iostream>
0002 #include <vector>
0003
0004 using namespace std;
0005
0006
0007 Int_t JetReso(TString infile="matched.root", TString outfile="resohists.root")
0008 {
0009
0010     //=======================================================
0011     //================ Input File ===========================
0012     //=======================================================
0013     std::cout << "Input file: " << infile.Data() << std::endl;
0014
0015     // open input file
0016     TFile fin(infile.Data(), "READ");
0017     if(!fin.IsOpen())
0018     {
0019         std::cout << "Error: could not open input file" << std::endl;
0020         exit(-1);
0021     }
0022
0023     // get tree
0024     TTree *tree = (TTree*)fin.Get("T");
0025     if(!tree)
0026     {
0027         std::cout << "Error: could not find tree" << std::endl;
0028         exit(-1);
0029     }
0030
0031     // get branches
0032     int centrality = 0;
0033     double rho_area = 0;
0034     double rho_mult = 0;
0035     float weight = 0;
0036     std::vector<float> *matched_pt_iter = 0;
0037     std::vector<float> *matched_pt_area = 0;
0038     std::vector<float> *matched_pt_mult = 0;
0039     std::vector<float> *matched_pt_iter_unsub = 0;
0040     std::vector<float> *matched_pt_area_unsub = 0;
0041     std::vector<float> *matched_pt_mult_unsub = 0;
0042     std::vector<float> *matched_truth_pt_iter = 0;
0043     std::vector<float> *matched_truth_pt_area = 0;
0044     std::vector<float> *matched_truth_pt_mult = 0;
0045     std::vector<float> *unmatched_pt_iter = 0;
0046     std::vector<float> *unmatched_pt_area = 0;
0047     std::vector<float> *unmatched_pt_mult = 0;
0048     std::vector<float> *unmatched_pt_iter_unsub = 0;
0049     std::vector<float> *unmatched_pt_area_unsub = 0;
0050     std::vector<float> *unmatched_pt_mult_unsub = 0;
0051     std::vector<float> *missed_truth_pt_iter = 0;
0052     std::vector<float> *missed_truth_pt_area = 0;
0053     std::vector<float> *missed_truth_pt_mult = 0;
0054
0055     tree->SetBranchAddress("weight", &weight);
0056     tree->SetBranchAddress("centrality", &centrality);
0057     tree->SetBranchAddress("rho_area", &rho_area);
0058     tree->SetBranchAddress("rho_mult", &rho_mult);
0059     tree->SetBranchAddress("matched_pt_iter", &matched_pt_iter);
0060     tree->SetBranchAddress("matched_pt_area", &matched_pt_area);
0061     tree->SetBranchAddress("matched_pt_mult", &matched_pt_mult);
0062     tree->SetBranchAddress("matched_pt_iter_unsub", &matched_pt_iter_unsub);
0063     tree->SetBranchAddress("matched_pt_area_unsub", &matched_pt_area_unsub);
0064     tree->SetBranchAddress("matched_pt_mult_unsub", &matched_pt_mult_unsub);
0065     tree->SetBranchAddress("matched_truth_pt_iter", &matched_truth_pt_iter);
0066     tree->SetBranchAddress("matched_truth_pt_area", &matched_truth_pt_area);
0067     tree->SetBranchAddress("matched_truth_pt_mult", &matched_truth_pt_mult);
0068     tree->SetBranchAddress("unmatched_pt_iter", &unmatched_pt_iter);
0069     tree->SetBranchAddress("unmatched_pt_area", &unmatched_pt_area);
0070     tree->SetBranchAddress("unmatched_pt_mult", &unmatched_pt_mult);
0071     tree->SetBranchAddress("unmatched_pt_iter_unsub", &unmatched_pt_iter_unsub);
0072     tree->SetBranchAddress("unmatched_pt_area_unsub", &unmatched_pt_area_unsub);
0073     tree->SetBranchAddress("unmatched_pt_mult_unsub", &unmatched_pt_mult_unsub);
0074     tree->SetBranchAddress("missed_truth_pt_iter", &missed_truth_pt_iter);
0075     tree->SetBranchAddress("missed_truth_pt_area", &missed_truth_pt_area);
0076     tree->SetBranchAddress("missed_truth_pt_mult", &missed_truth_pt_mult);
0077
0078     //=======================================================
0079     // declare histograms
0080     //=======================================================
0081     const double pt_range[] = {10,15,20,25,30,35,40,45,50,60,80};
0082     const int pt_N = sizeof(pt_range)/sizeof(double) - 1;
0083
0084     const int resp_N = 100;
0085     double resp_bins[resp_N+1];
0086     for(int i = 0; i < resp_N+1; i++){
0087         resp_bins[i] = 2.0/resp_N * i;
0088     }
0089
0090     const double cent_bins[] = {-1, 0, 10, 30, 50, 80}; //the first bin is for inclusive in centrality/pp events
0091     const int cent_N = sizeof(cent_bins)/sizeof(double) - 1;
0092
0093     TH3D * h3_reso_area = new TH3D("h3_reso_area", "h3_reso_area", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
0094     TH3D * h3_reso_mult = new TH3D("h3_reso_mult", "h3_reso_mult", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
0095     TH3D * h3_reso_iter = new TH3D("h3_reso_iter", "h3_reso_iter", pt_N, pt_range, resp_N, resp_bins, cent_N, cent_bins);
0096
0097     // fill histograms
0098     int nentries = tree->GetEntries();
0099     for (int iEvent = 0; iEvent < nentries; iEvent++){
0100         // get entry
0101         tree->GetEntry(iEvent);
0102
0103         for (unsigned int ijet =0; ijet < matched_pt_iter->size(); ijet++){
0104             h3_reso_iter->Fill(matched_truth_pt_iter->at(ijet), matched_pt_iter->at(ijet)/matched_truth_pt_iter->at(ijet), centrality, weight);
0105             h3_reso_iter->Fill(matched_truth_pt_iter->at(ijet), matched_pt_iter->at(ijet)/matched_truth_pt_iter->at(ijet), -1, weight);
0106         }
0107         for (unsigned int ijet =0; ijet < matched_pt_area->size(); ijet++){
0108             h3_reso_area->Fill(matched_truth_pt_area->at(ijet), matched_pt_area->at(ijet)/matched_truth_pt_area->at(ijet), centrality, weight);
0109             h3_reso_area->Fill(matched_truth_pt_area->at(ijet), matched_pt_area->at(ijet)/matched_truth_pt_area->at(ijet), -1, weight);
0110         }
0111         for (unsigned int ijet =0; ijet < matched_pt_mult->size(); ijet++){
0112             h3_reso_mult->Fill(matched_truth_pt_mult->at(ijet), matched_pt_mult->at(ijet)/matched_truth_pt_mult->at(ijet), centrality, weight);
0113             h3_reso_mult->Fill(matched_truth_pt_mult->at(ijet), matched_pt_mult->at(ijet)/matched_truth_pt_mult->at(ijet), -1, weight);
0114         }
0115     }
0116
0117     //=======================================================
0118     // JES and JER
0119     //=======================================================
0120     TH1D * h_jes_area[cent_N];
0121     TH1D * h_jer_area[cent_N];
0122     TH1D * h_jes_mult[cent_N];
0123     TH1D * h_jer_mult[cent_N];
0124     TH1D * h_jes_iter[cent_N];
0125     TH1D * h_jer_iter[cent_N];
0126
0127     for (int icent = 1; icent <= cent_N; icent++){
0128         h_jes_area[icent-1] = new TH1D(Form("h_jes_area_%d", icent), Form("h_jes_area_%d", icent), pt_N, pt_range);
0129         h_jer_area[icent-1] = new TH1D(Form("h_jer_area_%d", icent), Form("h_jer_area_%d", icent), pt_N, pt_range);
0130         h_jes_mult[icent-1] = new TH1D(Form("h_jes_mult_%d", icent), Form("h_jes_mult_%d", icent), pt_N, pt_range);
0131         h_jer_mult[icent-1] = new TH1D(Form("h_jer_mult_%d", icent), Form("h_jer_mult_%d", icent), pt_N, pt_range);
0132         h_jes_iter[icent-1] = new TH1D(Form("h_jes_iter_%d", icent), Form("h_jes_iter_%d", icent), pt_N, pt_range);
0133         h_jer_iter[icent-1] = new TH1D(Form("h_jer_iter_%d", icent), Form("h_jer_iter_%d", icent), pt_N, pt_range);
0134     }
0135
0136     for(int icent = 0; icent < cent_N; ++icent)
0137     {
0138         // area
0139         for (int i = 0; i < pt_N; ++i)
0140         {
0141             // fit function
0142             TF1 *func = new TF1("func","gaus",0, 1.2);
0143             h3_reso_area->GetXaxis()->SetRange(i+1,i+1);
0144             h3_reso_area->GetZaxis()->SetRange(icent+1,icent+1);
0145             TH1D * h_temp = (TH1D*)h3_reso_area->Project3D("y");
0146             h_temp->SetName(Form("h_jes_area_%d_%d", icent, i));
0147             h_temp->Fit(func,"","",0,1.2);
0148             h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
0149             float dsigma = func -> GetParError(2);
0150             float denergy = func -> GetParError(1);
0151             float sigma = func -> GetParameter(2);
0152             float energy = func -> GetParameter(1);
0153             float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
0154             h_jes_area[icent]->SetBinContent(i+1, energy);
0155             h_jes_area[icent]->SetBinError(i+1, denergy);
0156             h_jer_area[icent]->SetBinContent(i+1, sigma/energy);
0157             h_jer_area[icent]->SetBinError(i+1, djer);
0158         }
0159
0160         // mult
0161         for (int i = 0; i < pt_N; ++i)
0162         {
0163             // fit function
0164             TF1 *func = new TF1("func","gaus",0, 1.2);
0165             h3_reso_mult->GetXaxis()->SetRange(i+1,i+1);
0166             h3_reso_mult->GetZaxis()->SetRange(icent+1,icent+1);
0167             TH1D * h_temp = (TH1D*)h3_reso_mult->Project3D("y");
0168             h_temp->SetName(Form("h_jes_mult_%d_%d", icent, i));
0169             h_temp->Fit(func,"","",0,1.2);
0170             h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
0171             float dsigma = func -> GetParError(2);
0172             float denergy = func -> GetParError(1);
0173             float sigma = func -> GetParameter(2);
0174             float energy = func -> GetParameter(1);
0175             float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
0176             h_jes_mult[icent]->SetBinContent(i+1, energy);
0177             h_jes_mult[icent]->SetBinError(i+1, denergy);
0178             h_jer_mult[icent]->SetBinContent(i+1, sigma/energy);
0179             h_jer_mult[icent]->SetBinError(i+1, djer);
0180         }
0181
0182         // iter
0183         for (int i = 0; i < pt_N; ++i)
0184         {
0185             // fit function
0186             TF1 *func = new TF1("func","gaus",0, 1.2);
0187             h3_reso_iter->GetXaxis()->SetRange(i+1,i+1);
0188             h3_reso_iter->GetZaxis()->SetRange(icent+1,icent+1);
0189             TH1D * h_temp = (TH1D*)h3_reso_iter->Project3D("y");
0190             h_temp->SetName(Form("h_jes_iter_%d_%d", icent, i));
0191             h_temp->Fit(func,"","",0,1.2);
0192             h_temp->Fit(func,"","",func->GetParameter(1)-1.5*func->GetParameter(2),func->GetParameter(1)+1.5*func->GetParameter(2));
0193             float dsigma = func -> GetParError(2);
0194             float denergy = func -> GetParError(1);
0195             float sigma = func -> GetParameter(2);
0196             float energy = func -> GetParameter(1);
0197             float djer = dsigma/energy + (sigma*denergy)/(energy*energy);
0198             h_jes_iter[icent]->SetBinContent(i+1, energy);
0199             h_jes_iter[icent]->SetBinError(i+1, denergy);
0200             h_jer_iter[icent]->SetBinContent(i+1, sigma/energy);
0201             h_jer_iter[icent]->SetBinError(i+1, djer);
0202         }
0203
0204
0205     }
0206
0207     //=======================================================
0208     //================ Output File ==========================
0209     //=======================================================
0210     std::cout << "Output file: " << outfile.Data() << std::endl;
0211     TFile * fout = new TFile(outfile.Data(), "RECREATE");
0212     h3_reso_area->Write();
0213     h3_reso_mult->Write();
0214     h3_reso_iter->Write();
0215     for (int icent = 1; icent <= cent_N; icent++){
0216         h_jes_area[icent-1]->Write();
0217         h_jer_area[icent-1]->Write();
0218         h_jes_mult[icent-1]->Write();
0219         h_jer_mult[icent-1]->Write();
0220         h_jes_iter[icent-1]->Write();
0221         h_jer_iter[icent-1]->Write();
0222     }
0223     fout->Close();
0224
0225     // close input file
0226     fin.Close();
0227     return 0;
0228
0229
0230
0231 }