RecoParticleFlow/bin/offsetStack.cc

0001 #include "TString.h"
0002 #include "TFile.h"
0003 #include "TProfile.h"
0004 #include "TH1.h"
0005 #include "THStack.h"
0006 #include "TCanvas.h"
0007 #include "TStyle.h"
0008 #include "TLegend.h"
0009 #include "TLatex.h"
0010
0011 #include <iostream>
0012 #include <map>
0013
0014 bool getHists(std::map<TString, TH1D*>& hists, TFile*& file, TString var, int var_val, float r);
0015 void setStyle();
0016 void setStack(THStack*& stack, std::map<TString, TH1D*>& hists);
0017 void split(const TString& str, TString& sub1, TString& sub2);
0018
0019 int main(int argc, char* argv[]) {
0020   if (argc < 5) {
0021     std::cout << "Please provide a label:filename, variable, deltaR value, and an out directory." << std::endl;
0022     std::cout << "You may also include a second file at the end." << std::endl;
0023     std::cout << "Example: offsetStack label1:file1.root npv 0.4 plots" << std::endl;
0024     std::cout << "Example: offsetStack label1:file1.root npv 0.4 plots label2:file2.root" << std::endl;
0025     return -1;
0026   }
0027
0028   TString label1, fname1;
0029   split(argv[1], label1, fname1);
0030   TString label = label1;
0031
0032   TString var = argv[2];
0033   float r = std::stof(argv[3]);
0034   TString outdir = argv[4];
0035
0036   TFile* file1 = TFile::Open(fname1);
0037   if (!file1) {
0038     std::cout << "Invalid file1: " << fname1 << std::endl;
0039     return -1;
0040   }
0041
0042   TH1F* h = (TH1F*)file1->FindObjectAny(var);
0043   if (!h) {
0044     std::cout << "Could not find " + var + " hist" << std::endl;
0045     return -1;
0046   }
0047   int avg = h->GetMean() + 0.5;
0048
0049   std::map<TString, TH1D*> hists1 = {{"chm", 0}, {"chu", 0}, {"nh", 0}, {"ne", 0}, {"hfh", 0}, {"hfe", 0}, {"lep", 0}};
0050   if (!getHists(hists1, file1, var, avg, r))
0051     return -1;
0052
0053   setStyle();
0054   TCanvas c("c", "c", 600, 600);
0055
0056   THStack* stack1 = new THStack(
0057       "stack1", Form(";#eta;<Offset Energy_{T}(#pi%.1f^{2})>/<%s> [GeV]", r, var == "npv" ? "N_{PV}" : "#mu"));
0058   setStack(stack1, hists1);
0059   stack1->Draw("hist");
0060   TString legPF = "F";
0061
0062   TLegend* leg = new TLegend(.4, .67, .65, .92);
0063
0064   //file2 included//
0065   if (argc > 5) {
0066     TString label2, fname2;
0067     split(argv[5], label2, fname2);
0068
0069     TFile* file2 = TFile::Open(fname2);
0070     if (!file2) {
0071       std::cout << "Invalid file2: " << fname2 << std::endl;
0072       return -1;
0073     }
0074
0075     TH1F* h2 = (TH1F*)file2->FindObjectAny(var);
0076     if (!h2) {
0077       std::cout << "Could not find " + var + " hist in file2" << std::endl;
0078       return -1;
0079     }
0080     int avg2 = h2->GetMean() + 0.5;
0081
0082     std::map<TString, TH1D*> hists2 = hists1;
0083     if (!getHists(hists2, file2, var, avg2, r))
0084       return -1;
0085
0086     THStack* stack2 = new THStack("stack2", "stack2");
0087     setStack(stack2, hists2);
0088     stack2->Draw("samepe");
0089
0090     legPF = "PF";
0091     leg->SetHeader(Form("#bf{Markers: %s, Histograms: %s}", label2.Data(), label1.Data()));
0092     label = label1 + "_vs_" + label2;
0093
0094     //Draw Markers for EM Deposits and Hadronic Deposits in two separate regions//
0095     TH1D* hfe_clone = (TH1D*)hists2["hfe"]->Clone("hfe_clone");
0096     TH1D* hfh_clone = (TH1D*)hists2["hfh"]->Clone("hfh_clone");
0097
0098     THStack* cloneStack = new THStack("cloneStack", "cloneStack");
0099     cloneStack->Add(hists2["ne"]);
0100     cloneStack->Add(hfe_clone);
0101     cloneStack->Add(hists2["nh"]);
0102     cloneStack->Add(hfh_clone);
0103     cloneStack->Draw("samepe");
0104
0105     hists2["ne"]->SetAxisRange(-2.9, 2.9);
0106     hists2["hfe"]->SetAxisRange(-5, -2.6);
0107     hists2["nh"]->SetAxisRange(-2.9, 2.9);
0108     hists2["hfh"]->SetAxisRange(-5, -2.6);
0109     hists2["chu"]->SetAxisRange(-2.9, 2.9);
0110     hists2["chm"]->SetAxisRange(-2.9, 2.9);
0111
0112     hfe_clone->SetAxisRange(2.6, 5);
0113     hfh_clone->SetAxisRange(2.6, 5);
0114   }
0115
0116   leg->SetBorderSize(0);
0117   leg->SetFillColor(0);
0118   leg->SetFillStyle(0);
0119   leg->SetTextSize(0.04);
0120   leg->SetTextFont(42);
0121
0122   leg->AddEntry(hists1["ne"], "Photons", legPF);
0123   leg->AddEntry(hists1["hfe"], "EM Deposits", legPF);
0124   leg->AddEntry(hists1["nh"], "Neutral Hadrons", legPF);
0125   leg->AddEntry(hists1["hfh"], "Hadronic Deposits", legPF);
0126   leg->AddEntry(hists1["chu"], "Unassoc. Charged Hadrons", legPF);
0127   leg->AddEntry(hists1["chm"], "Assoc. Charged Hadrons", legPF);
0128
0129   leg->Draw();
0130
0131   TLatex text;
0132   text.SetNDC();
0133
0134   text.SetTextSize(0.065);
0135   text.SetTextFont(61);
0136   text.DrawLatex(0.2, 0.87, "CMS");
0137
0138   text.SetTextSize(0.045);
0139   text.SetTextFont(42);
0140   text.DrawLatex(1 - label.Length() / 41., 0.96, label);
0141
0142   c.Print(outdir + "/stack_" + label + ".pdf");
0143 }
0144
0145 bool getHists(std::map<TString, TH1D*>& hists, TFile*& file, TString var, int var_val, float r) {
0146   for (auto& pair : hists) {
0147     TString name = Form("p_offset_eta_%s%i_%s", var.Data(), var_val, pair.first.Data());
0148     TProfile* p = (TProfile*)file->FindObjectAny(name);
0149     if (!p) {
0150       std::cout << "Could not find " + name << std::endl;
0151       return false;
0152     }
0153     pair.second = p->ProjectionX(pair.first);
0154     pair.second->Scale(r * r / 2 / var_val);
0155
0156     const double* xbins = p->GetXaxis()->GetXbins()->GetArray();
0157     for (int i = 1, n = p->GetNbinsX(); i <= n; i++) {
0158       pair.second->SetBinContent(i, pair.second->GetBinContent(i) / (xbins[i] - xbins[i - 1]));
0159       pair.second->SetBinError(i, pair.second->GetBinError(i) / (xbins[i] - xbins[i - 1]));
0160     }
0161   }
0162   return true;
0163 }
0164
0165 void setStyle() {
0166   gStyle->SetPadTopMargin(0.05);
0167   gStyle->SetPadBottomMargin(0.1);
0168   gStyle->SetPadLeftMargin(0.16);
0169   gStyle->SetPadRightMargin(0.02);
0170
0171   gStyle->SetOptStat(0);
0172   gStyle->SetOptTitle(0);
0173
0174   gStyle->SetTitleFont(42, "XYZ");
0175   gStyle->SetTitleSize(0.05, "XYZ");
0176   gStyle->SetTitleXOffset(0.9);
0177   gStyle->SetTitleYOffset(1.4);
0178
0179   gStyle->SetLabelFont(42, "XYZ");
0180   gStyle->SetLabelOffset(0.007, "XYZ");
0181   gStyle->SetLabelSize(0.04, "XYZ");
0182
0183   gStyle->SetPadTickX(1);
0184   gStyle->SetPadTickY(1);
0185 }
0186
0187 void setStack(THStack*& stack, std::map<TString, TH1D*>& hists) {
0188   stack->Add(hists["ne"]);
0189   stack->Add(hists["hfe"]);
0190   stack->Add(hists["nh"]);
0191   stack->Add(hists["hfh"]);
0192   stack->Add(hists["chu"]);
0193   stack->Add(hists["chm"]);
0194
0195   hists["ne"]->SetMarkerStyle(kMultiply);
0196   hists["hfe"]->SetMarkerStyle(kOpenStar);
0197   hists["nh"]->SetMarkerStyle(kOpenDiamond);
0198   hists["hfh"]->SetMarkerStyle(kOpenTriangleUp);
0199   hists["chu"]->SetMarkerStyle(kOpenCircle);
0200   hists["chm"]->SetMarkerStyle(kOpenCircle);
0201
0202   hists["ne"]->SetFillColor(kBlue);
0203   hists["hfe"]->SetFillColor(kViolet + 2);
0204   hists["nh"]->SetFillColor(kGreen);
0205   hists["hfh"]->SetFillColor(kPink + 6);
0206   hists["chu"]->SetFillColor(kRed - 9);
0207   hists["chm"]->SetFillColor(kRed);
0208
0209   hists["ne"]->SetLineColor(kBlack);
0210   hists["hfe"]->SetLineColor(kBlack);
0211   hists["nh"]->SetLineColor(kBlack);
0212   hists["hfh"]->SetLineColor(kBlack);
0213   hists["chu"]->SetLineColor(kBlack);
0214   hists["chm"]->SetLineColor(kBlack);
0215 }
0216
0217 void split(const TString& str, TString& sub1, TString& sub2) {
0218   int pos = str.First(':');
0219   sub1 = str(0, pos);
0220   sub2 = str(pos + 1, str.Length());
0221 }