ForwardCalo/JetEnergyResolution/fjetResol.C

0001 /*
0002    Auther: Chong Kim (ckim.phenix@gmail.com)
0003
0004    - Target: study resolution of 1st/2nd leading jets reconstructed in fsPHENIX calorimeter
0005    - Environments:
0006      a. Requires root file contains reconstructed forward jets, which is results of
0007         /macros/macros/g4simulations/Fun4All_G4_fsPHENIX.C
0008         * modified Fun4All_G4_fsPHENIX.C used for this code can be found at:
0009         /gpfs/mnt/gpfs04/sphenix/user/ckim/fjets/Fun4All_G4_fsPHENIX.C
0010      b. Uses root branch "ntp_truthjet"
0011
0012    - How code works:
0013      a. List up 1st/2nd leading jets using true pT (gpt) which satisfying min pT cut (cutPT)
0014         (function searchLeadingJets)
0015      b. Study target variable (ex. energy) using collected 1st/2nd leading jets in each event
0016         For now only energy is checked but any variable in "ntp_truthjet" can be studied
0017      c. Loop through multiple events and root files
0018
0019    - Parameters (ex. Pythia version, beam energy, Anti_kT radius, etc) can be arranged via arrays in top
0020    - Tested by using output files in "/gpfs/mnt/gpfs04/sphenix/user/ckim/fjets/Ana/results_1121" at Mar.2, 2017
0021 */
0022
0023 #include <TCanvas.h>
0024 #include <TFile.h>
0025 #include <TF1.h>
0026 #include <TGraph.h>
0027 #include <TGraphErrors.h>
0028 #include <TH1.h>
0029 #include <TH2.h>
0030 #include <TLatex.h>
0031 #include <TLegend.h>
0032 #include <TNtuple.h>
0033 #include <TString.h>
0034 #include <TStyle.h>
0035 #include <iostream>
0036 using namespace std;
0037
0038 std::vector<int> searchLeadingJets(TNtuple *T, float cutPT, const char *targetV = "gpt");
0039
0040 //=========================================================
0041 void fjetResol(
0042         const char* path   = "/gpfs/mnt/gpfs04/sphenix/user/ckim/fjets/Ana/results_1121",
0043         const char* var    = "e",
0044         const char* suffix = "",
0045         bool print = true
0046         )
0047 {
0048
0049     //Parameters
0050     const int   setPy[]    = {8/*, 6*/}; //Pythia8 or 6
0051     const int   setBeamE[] = {200/*, 510*/};
0052     const int   setR[]     = {4, 6}; //Anti-kT radius
0053     const int   setMinPT[] = {/*3, 4,*/ 5, 10, 15}; //Minimum parton pT in pythia cfg
0054     const float setEta[]   = {1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6};
0055     const float setE[]     = {20,25,30,35,40,45,50, 60,70,80,90,100,110,120,130,140,150,160, 180,200}; //True jet E
0056
0057     const int nsetPy    = sizeof(setPy)   /sizeof(setPy[0]);
0058     const int nsetBeamE = sizeof(setBeamE)/sizeof(setBeamE[0]);
0059     const int nsetR     = sizeof(setR)    /sizeof(setR[0]);
0060     const int nsetMinPT = sizeof(setMinPT)/sizeof(setMinPT[0]);
0061     const int nsetEta   = sizeof(setEta)  /sizeof(setEta[0]) - 1;
0062     const int nsetE     = sizeof(setE)    /sizeof(setE[0]) - 1;
0063
0064     //----------------------------------------------
0065
0066     //Histograms to be filled
0067     TH2F *h2[nsetPy][nsetBeamE][nsetR][nsetEta];
0068     TH1F *h1[nsetPy][nsetBeamE][nsetR][nsetEta][nsetE];
0069     for (int a=0; a<nsetPy;    a++)
0070     for (int b=0; b<nsetBeamE; b++)
0071     for (int c=0; c<nsetR;     c++)
0072     for (int x=0; x<nsetEta;   x++)
0073     {
0074         const char* h2Name  = Form("sanity_py%i_%i_r0%i_eta%i", setPy[a],setBeamE[b],setR[c],x);
0075         h2[a][b][c][x] = new TH2F(h2Name, "", 100,0,b==0?100:200, 60,0,2);
0076         for (int y=0; y<nsetE; y++)
0077         {
0078             const char* h1Name = Form("h1_%s_py%i_%i_r0%i_eta%i_e%i", var,setPy[a],setBeamE[b],setR[c],x,y);
0079             h1[a][b][c][x][y] = new TH1F(h1Name, "", 150,-1.5,1.5);
0080         }
0081     }
0082
0083     //----------------------------------------------
0084
0085     //Iterate through each file and tree
0086     for (int a=0; a<nsetPy;    a++)
0087     for (int b=0; b<nsetBeamE; b++)
0088     for (int c=0; c<nsetR;     c++)
0089     for (int d=0; d<nsetMinPT; d++)
0090     {
0091         //Link input
0092         const char *finName = Form("%s/pythia%i_%i_r0%i_pT%i.root", path,setPy[a],setBeamE[b],setR[c],setMinPT[d]);
0093         TFile   *F = TFile::Open(finName);
0094         TNtuple *T = (TNtuple*)F->Get("ntp_truthjet");
0095         cout <<Form("Open %s from %s...", T->GetName(), F->GetName()) <<endl;
0096
0097         //Search 1st/2nd leading jets in this file
0098         std::vector<int> jetList = searchLeadingJets(T, (float)setMinPT[d]);
0099         //for (unsigned int i=0; i<jetList.size(); i++) cout <<i <<" " <<jetList[i] <<endl;
0100
0101         //Link branches
0102         float event, geta, gphi, ge, eta, phi, e;
0103         T->SetBranchAddress("event", &event);
0104         T->SetBranchAddress("geta",  &geta);
0105         T->SetBranchAddress("gphi",  &gphi);
0106         T->SetBranchAddress("ge",    &ge);
0107         T->SetBranchAddress("eta",   &eta);
0108         T->SetBranchAddress("phi",   &phi);
0109         T->SetBranchAddress("e",     &e);
0110
0111         //Iterate through for leading jets found
0112         for (unsigned int i=0; i<jetList.size(); i++)
0113         {
0114             T->GetEntry(jetList[i]);
0115
0116             int id_eta = -1;
0117             for (int x=0; x<nsetEta; x++) { if (geta > setEta[x] && geta < setEta[x+1]) id_eta = x; }
0118             if (id_eta == -1) continue;
0119
0120             h2[a][b][c][id_eta]->Fill(ge, e/ge); //Sanity plots
0121
0122             int id_e = -1;
0123             for (int y=0; y<nsetE; y++) { if (ge > setE[y] && ge < setE[y+1]) id_e = y; }
0124             if (id_e == -1) continue;
0125
0126             if (!strcmp(var, "e")) h1[a][b][c][id_eta][id_e]->Fill(e/ge);
0127         }//i, leading jets
0128
0129         //Cleanup
0130         T->ResetBranchAddresses();
0131         T->Delete();
0132         F->Delete();
0133     }//d, c, b, a
0134
0135     //----------------------------------------------
0136
0137     gStyle->SetOptStat("e");
0138
0139     //Sanity check plots
0140     TCanvas *c1[nsetPy][nsetBeamE][nsetR];
0141     for (int a=0; a<nsetPy;    a++)
0142     for (int b=0; b<nsetBeamE; b++)
0143     for (int c=0; c<nsetR;     c++)
0144     {
0145         const char* c1Name = Form("sanity_py%i_%i_r0%i", setPy[a],setBeamE[b],setR[c]);
0146         c1[a][b][c] = new TCanvas(c1Name, c1Name, 1280, 800);
0147         c1[a][b][c]->Divide(nsetEta/2, 2, 0.015, 0.015, 10);
0148         for (int x=0; x<nsetEta; x++)
0149         {
0150             c1[a][b][c]->cd(x+1)->SetGrid();
0151             h2[a][b][c][x]->SetTitle(Form("%2.1f < #eta < %2.1f;True e;Reco e/True e", setEta[x], setEta[x+1]));
0152             h2[a][b][c][x]->GetYaxis()->SetTitleOffset(1.45);
0153             h2[a][b][c][x]->Draw("colz");
0154         }
0155         if (print) c1[a][b][c]->Print(Form("%s%s.png", c1[a][b][c]->GetName(),suffix));
0156     }//c, b, a
0157
0158     //Get resolutions
0159     TCanvas *c2 = new TCanvas(Form("resol_%s", var), var, 1280, 800);
0160     c2->Divide(nsetEta/2, 2, 0.015, 0.015, 10);
0161     TF1 *f1Gaus = new TF1("f1Gaus", "[0]*TMath::Gaus(x, [1], [2])", -0.2, 0.2);
0162     TLegend *leg1 = new TLegend(0.4, 0.6/*0.9 - 0.15*nsetPy*/, 0.9, 0.9);
0163     for (int a=0; a<nsetPy;    a++)
0164     for (int b=0; b<nsetBeamE; b++)
0165     for (int c=0; c<nsetR;     c++)
0166     for (int x=0; x<nsetEta;   x++)
0167     {
0168         TGraphErrors *g1 = new TGraphErrors();
0169         g1->SetName(Form("g1_%s_py%i_%i_r0%i_eta%i_e%i", var,setPy[a],setBeamE[b],setR[c],x));
0170         g1->SetLineColor(2 * (b+1));
0171         g1->SetMarkerColor(2 * (b+1));
0172         g1->SetMarkerSize(1.);
0173         g1->SetMarkerStyle(20 + 4*c);
0174
0175         for (int y=0; y<nsetE; y++)
0176         {
0177             if (h1[a][b][c][x][y]->GetEntries() < 50)
0178             {
0179                 h1[a][b][c][x][y]->Delete();
0180                 continue;
0181             }
0182             f1Gaus->SetParameters(
0183                     h1[a][b][c][x][y]->GetMaximum(),
0184                     h1[a][b][c][x][y]->GetMean(),
0185                     h1[a][b][c][x][y]->GetRMS()
0186                     );
0187             h1[a][b][c][x][y]->Fit(
0188                     f1Gaus->GetName(), "QR0", "",
0189                     h1[a][b][c][x][y]->GetMean() - h1[a][b][c][x][y]->GetRMS() * 3,
0190                     h1[a][b][c][x][y]->GetMean() + h1[a][b][c][x][y]->GetRMS() * 3
0191                     );
0192             g1->SetPoint(g1->GetN(), (setE[y]+setE[y+1])/2, f1Gaus->GetParameter(2));
0193             g1->SetPointError(g1->GetN()-1, 0, f1Gaus->GetParError(2));
0194         }//y (energy)
0195
0196         c2->cd(x+1)->SetGrid();
0197         if (a==0 && b==0 && c==0)
0198         {
0199             const char *yTitle = "#sigma (Reco e / True e)";
0200             const char *gTitle = Form("%2.1f < #eta < %2.1f;True e;%s", setEta[x],setEta[x+1], yTitle);
0201             gPad->DrawFrame(20-5, 0.05, 200+5, 0.3, gTitle);
0202             gPad->Modified();
0203         }
0204         if (a==0 && x==0)
0205         {
0206             leg1->AddEntry(g1, Form("PYTHIA%i, #sqrt{s} = %i, R = 0.%i", setPy[a], setBeamE[b], setR[c]), "p");
0207         }
0208         g1->Draw("pe same");
0209         if (a==(nsetPy-1) && b==(nsetBeamE-1) && c==(nsetR-1) && x==0) leg1->Draw("same");
0210     }//x (eta), c, b, a
0211     if (print) c2->Print(Form("%s%s.png", c2->GetName(),suffix));
0212
0213     return;
0214 }//Main
0215
0216 //==============================================================================
0217 std::vector<int> searchLeadingJets(TNtuple *T, float cutPT, const char *targetV)
0218 {
0219     std::vector<int> jetList; //Return (event # of target jets)
0220
0221     float event, gpt, var;
0222     T->SetBranchAddress("event", &event);
0223     if (strcmp("gpt", targetV))
0224     {
0225         T->SetBranchAddress("gpt", &gpt);
0226         cout <<"Search leading jets by using " <<targetV <<endl;
0227     }
0228     T->SetBranchAddress(targetV, &var); //gpt or ge (truth info)
0229
0230     //Note: high1_v should always higher than high2_v
0231     int evtCheck = -1;
0232     int high1_n  = -1; //1st leading jet's entry # in the tree
0233     int high2_n  = -1;
0234     float high1_v = -1; //1st leading jet's value
0235     float high2_v = -2;
0236
0237     const int allEntries = T->GetEntries();
0238     for (int a=0; a<allEntries; a++)
0239     {
0240         /*
0241            Iterate through all entries with same event number,
0242            find 1st and 2nd leading jets of the event,
0243            then save them when next event (different event #) found
0244         */
0245
0246         T->GetEntry(a);
0247
0248         //pT cut
0249         float tempPT = (strcmp("gpt", targetV))?gpt:var;
0250         if (tempPT < cutPT) continue;
0251
0252         if (evtCheck != event || a == allEntries) //Save jet entries
0253         {
0254             if (high1_n != -1 && high2_n != -1)
0255             {
0256                 if (high1_n < high2_n)
0257                 {
0258                     jetList.push_back(high1_n);
0259                     jetList.push_back(high2_n);
0260                 }
0261                 else
0262                 {
0263                     jetList.push_back(high2_n);
0264                     jetList.push_back(high1_n);
0265                 }
0266             }
0267
0268             //Update dummy indices for next iteration
0269             evtCheck = event;
0270             high1_n = a;
0271             high2_n = -1;
0272             high1_v = var;
0273             high2_v = -1;
0274         }
0275         else //Iterate & Update
0276         {
0277             if (var > high1_v)
0278             {
0279                 high2_n = high1_n;
0280                 high2_v = high1_v;
0281                 high1_n = a;
0282                 high1_v = var;
0283             }
0284             else if (var > high2_v)
0285             {
0286                 high2_n = a;
0287                 high2_v = var;
0288             }
0289         }
0290     }//a, nEvents
0291
0292     T->ResetBranchAddresses();
0293     return jetList;
0294 }//searchLeadingJets