Skip to content

Commit dff7f40

Browse files
mverwemverwe
authored
Merge pull request #10 from mverwe/master
CS variations of alpha
2 parents 351ab34 + 3339d6f commit dff7f40

File tree

2 files changed

+340
-0
lines changed

2 files changed

+340
-0
lines changed

‎plot/plotJetEnergyScaleCS.C

+143
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
//Quick dirty plotting macro for jet response from trees of JetToyHI framework
2+
3+
#include <iostream>
4+
#include <vector>
5+
#include <string>
6+
#include <algorithm>
7+
#include <fstream>
8+
9+
#include "TCanvas.h"
10+
#include "TFile.h"
11+
#include "TH3F.h"
12+
#include "TH1F.h"
13+
#include "TH1D.h"
14+
#include "TGraphAsymmErrors.h"
15+
#include "TLegend.h"
16+
17+
using namespace std;
18+
19+
int markerColor[11] = {1,kRed+1,4,kGreen+3,kOrange+7,kGray+2,kRed+5,kAzure+10,kGreen+4,kYellow+2,kGreen};
20+
int markerStyle[8] = {20,21,33,34,24,25,27,28};
21+
22+
TH1F *DrawFrame(double xmin = 0., double xmax = 1., double ymin = 0., double ymax = 1., TString xTitle = "x", TString yTitle = "y", bool setMargins = true);
23+
24+
void plotJetEnergyScaleCS(TString str = "JetToyHIResult.root") {
25+
26+
TFile *f = new TFile(str.Data());
27+
TTree *tr = dynamic_cast<TTree*>(f->Get("jetTree"));
28+
29+
tr->SetLineWidth(2);
30+
31+
int nEvt = tr->GetEntriesFast();
32+
33+
const int ncs = 5;
34+
35+
TString strLeg[ncs] = {"CS #alpha=-2",
36+
"CS #alpha=-1",
37+
"CS #alpha=0",
38+
"CS #alpha=1",
39+
"CS #alpha=2",
40+
};
41+
42+
43+
TH1 *hResponsePt[ncs] = {0};
44+
TH1 *hResponseM[ncs] = {0};
45+
46+
//----------------------------------------------------------------
47+
// plot the jet pT response
48+
//----------------------------------------------------------------
49+
TCanvas *c1 = new TCanvas("c1","c1",450,400);
50+
gPad->SetLogy();
51+
tr->SetLineColor(1);
52+
tr->Draw("csJetPt/sigJetPt>>hCSJetResponsePt(100,0.,2.)","sigJetPt>120. && abs(sigJetEta)<2.3");
53+
hResponsePt[0] = dynamic_cast<TH1*>(gDirectory->Get("hCSJetResponsePt"));
54+
55+
for(int i = 0; i<ncs; ++i) {
56+
tr->Draw(Form("csJet_%dPt/sigJetPt>>hCS%dJetResponsePt(100,0.,2.)",i,i),"sigJetPt>120. && abs(sigJetEta)<2.3","same");
57+
hResponsePt[i] = dynamic_cast<TH1*>(gDirectory->Get(Form("hCS%dJetResponsePt",i)));
58+
}
59+
60+
for(int i = 0; i<ncs; ++i) {
61+
if(!hResponsePt[i]) continue;
62+
hResponsePt[i]->Sumw2();
63+
if(nEvt>0) hResponsePt[i]->Scale(1./(double)nEvt);
64+
}
65+
66+
TH1F *fr1 = DrawFrame(0.,2.,5e-5,5e-1,"p_{T,sub}/p_{T,gen}","N/N_{evt}");
67+
TLegend *leg1 = new TLegend(0.23,0.6,0.45,0.93);
68+
leg1->SetFillColor(10);
69+
leg1->SetBorderSize(0);
70+
leg1->SetFillStyle(0);
71+
leg1->SetTextSize(0.05);
72+
for(int i = 0; i<ncs; ++i) {
73+
if(!hResponsePt[i]) continue;
74+
hResponsePt[i]->SetLineColor(markerColor[i]);
75+
hResponsePt[i]->SetMarkerColor(markerColor[i]);
76+
hResponsePt[i]->SetMarkerStyle(markerStyle[i]);
77+
hResponsePt[i]->Draw("same");
78+
leg1->AddEntry(hResponsePt[i],strLeg[i].Data(),"p");
79+
}
80+
leg1->Draw();
81+
82+
//----------------------------------------------------------------
83+
// plot the jet mass response
84+
//----------------------------------------------------------------
85+
TCanvas *c4 = new TCanvas("c4","c4",450,400);
86+
gPad->SetLogy();
87+
tr->Draw("csJetM/sigJetM>>hCSJetResponseM(100,0.,5.)","sigJetPt>120. && abs(sigJetEta)<2.3");
88+
hResponseM[0] = dynamic_cast<TH1*>(gDirectory->Get("hCSJetResponseM"));
89+
90+
for(int i =0; i<nsk; ++i) {
91+
tr->Draw(Form("csJet_%dM/sigJetM>>hCS%dJetResponseM(100,0.,2.)",i,i),"sigJetPt>120. && abs(sigJetEta)<2.3","same");
92+
hResponseM[i] = dynamic_cast<TH1*>(gDirectory->Get(Form("hCS%dJetResponseM",i)));
93+
}
94+
95+
for(int i = 0; i<ncs; ++i) {
96+
if(!hResponseM[i]) continue;
97+
hResponseM[i]->Sumw2();
98+
if(nEvt>0) hResponseM[i]->Scale(1./(double)nEvt);
99+
}
100+
101+
TH1F *fr2 = DrawFrame(0.,5.,5e-5,5e-1,"M_{sub}/M_{gen}","N/N_{evt}");
102+
TLegend *leg2 = new TLegend(0.5,0.6,0.93,0.93);
103+
leg2->SetFillColor(10);
104+
leg2->SetBorderSize(0);
105+
leg2->SetFillStyle(0);
106+
leg2->SetTextSize(0.05);
107+
for(int i = 0; i<ncs; ++i) {
108+
if(!hResponseM[i]) continue;
109+
hResponseM[i]->SetLineColor(markerColor[i]);
110+
hResponseM[i]->SetMarkerColor(markerColor[i]);
111+
hResponseM[i]->SetMarkerStyle(markerStyle[i]);
112+
hResponseM[i]->Draw("same");
113+
leg2->AddEntry(hResponseM[i],strLeg[i].Data(),"p");
114+
}
115+
leg2->Draw();
116+
117+
}
118+
119+
TH1F *DrawFrame(double xmin, double xmax, double ymin, double ymax, TString xTitle, TString yTitle, bool setMargins) {
120+
121+
if(setMargins) {
122+
gPad->SetLeftMargin(0.22);
123+
gPad->SetBottomMargin(0.15);
124+
gPad->SetRightMargin(0.1);//0.05);
125+
gPad->SetTopMargin(0.05);
126+
}
127+
128+
TH1F *frame = gPad->DrawFrame(xmin,ymin,xmax,ymax);
129+
frame->SetXTitle(xTitle.Data());
130+
frame->SetYTitle(yTitle.Data());
131+
frame->GetXaxis()->SetLabelSize(0.05);
132+
frame->GetYaxis()->SetLabelSize(0.05);
133+
frame->GetXaxis()->SetTitleSize(0.06);
134+
frame->GetYaxis()->SetTitleSize(0.06);
135+
frame->GetXaxis()->SetTitleOffset(1.0);
136+
frame->GetYaxis()->SetTitleOffset(1.3);
137+
frame->GetXaxis()->CenterTitle(true);
138+
frame->GetYaxis()->CenterTitle(true);
139+
140+
gPad->SetTicks(1,1);
141+
142+
return frame;
143+
}

‎runCSVariations.cc

+197
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
#include <iostream>
2+
#include <chrono>
3+
4+
#include "TFile.h"
5+
#include "TTree.h"
6+
7+
#include "fastjet/PseudoJet.hh"
8+
#include "fastjet/ClusterSequenceArea.hh"
9+
10+
#include "include/ProgressBar.h"
11+
12+
#include "PU14/EventMixer.hh"
13+
#include "PU14/CmdLine.hh"
14+
#include "PU14/PU14.hh"
15+
16+
#include "include/jetCollection.hh"
17+
#include "include/csSubtractor.hh"
18+
#include "include/csSubtractorFullEvent.hh"
19+
#include "include/skSubtractor.hh"
20+
#include "include/softDropGroomer.hh"
21+
#include "include/treeWriter.hh"
22+
#include "include/jetMatcher.hh"
23+
#include "include/randomCones.hh"
24+
#include "include/Angularity.hh"
25+
26+
using namespace std;
27+
using namespace fastjet;
28+
29+
// ./runCSVariations -hard /eos/project/j/jetquenching/JetWorkshop2017/samples/pythia8/dijet120/PythiaEventsTune14PtHat120_0.pu14 -pileup /eos/project/j/jetquenching/JetWorkshop2017/samples/thermal/Mult7000/ThermalEventsMult7000PtAv1.20_0.pu14 -nev 10
30+
31+
int main (int argc, char ** argv) {
32+
33+
auto start_time = std::chrono::steady_clock::now();
34+
35+
CmdLine cmdline(argc,argv);
36+
// inputs read from command line
37+
int nEvent = cmdline.value<int>("-nev",1); // first argument: command line option; second argument: default value
38+
//bool verbose = cmdline.present("-verbose");
39+
40+
std::cout << "will run on " << nEvent << " events" << std::endl;
41+
42+
// Uncomment to silence fastjet banner
43+
ClusterSequence::set_fastjet_banner_stream(NULL);
44+
45+
//to write info to root tree
46+
treeWriter trw("jetTree");
47+
48+
//Jet definition
49+
double R = 0.4;
50+
double ghostRapMax = 6.0;
51+
double ghost_area = 0.005;
52+
int active_area_repeats = 1;
53+
fastjet::GhostedAreaSpec ghost_spec(ghostRapMax, active_area_repeats, ghost_area);
54+
fastjet::AreaDefinition area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec);
55+
fastjet::JetDefinition jet_def(antikt_algorithm, R);
56+
57+
double jetRapMax = 3.0;
58+
fastjet::Selector jet_selector = SelectorAbsRapMax(jetRapMax);
59+
60+
Angularity width(1.,1.,R);
61+
Angularity pTD(0.,2.,R);
62+
63+
ProgressBar Bar(cout, nEvent);
64+
Bar.SetStyle(-1);
65+
66+
EventMixer mixer(&cmdline); //the mixing machinery from PU14 workshop
67+
68+
// loop over events
69+
int iev = 0;
70+
unsigned int entryDiv = (nEvent > 200) ? nEvent / 200 : 1;
71+
while ( mixer.next_event() && iev < nEvent )
72+
{
73+
// increment event number
74+
iev++;
75+
76+
Bar.Update(iev);
77+
Bar.PrintWithMod(entryDiv);
78+
79+
std::vector<fastjet::PseudoJet> particlesMerged = mixer.particles();
80+
81+
std::vector<double> eventWeight;
82+
eventWeight.push_back(mixer.hard_weight());
83+
eventWeight.push_back(mixer.pu_weight());
84+
85+
// cluster hard event only
86+
std::vector<fastjet::PseudoJet> particlesBkg, particlesSig;
87+
SelectorIsHard().sift(particlesMerged, particlesSig, particlesBkg); // this sifts the full event into two vectors of PseudoJet, one for the hard event, one for the underlying event
88+
89+
//---------------------------------------------------------------------------
90+
// jet clustering
91+
//---------------------------------------------------------------------------
92+
93+
// run the clustering, extract the signal jets
94+
fastjet::ClusterSequenceArea csSig(particlesSig, jet_def, area_def);
95+
jetCollection jetCollectionSig(sorted_by_pt(jet_selector(csSig.inclusive_jets(10.))));
96+
97+
//---------------------------------------------------------------------------
98+
// background subtraction
99+
//---------------------------------------------------------------------------
100+
101+
//run jet-by-jet constituent subtraction on mixed (hard+UE) event
102+
const int ncs = 5;
103+
double alpha[ncs] = {-2.,-1.,0.,1.,2.};
104+
std::vector<jetCollection> jetCollectionCSs;
105+
std::vector<double> rho;
106+
std::vector<double> rhom;
107+
for(int ics = 0; ics<ncs; ++ics) {
108+
csSubtractor csSub(R, alpha[ics], -1, 0.005,ghostRapMax,jetRapMax);
109+
csSub.setInputParticles(particlesMerged);
110+
jetCollection jetCollectionCS(csSub.doSubtraction());
111+
112+
if(ics==3) {
113+
//Background densities used by constituent subtraction
114+
rho.push_back(csSub.getRho());
115+
rhom.push_back(csSub.getRhoM());
116+
}
117+
jetCollectionCSs.push_back(jetCollectionCS);
118+
}
119+
120+
121+
//---------------------------------------------------------------------------
122+
// Groom the jets
123+
//---------------------------------------------------------------------------
124+
125+
126+
//SoftDrop grooming classic for signal jets (zcut=0.1, beta=0)
127+
softDropGroomer sdgSig(0.1, 0.0, R);
128+
jetCollection jetCollectionSigSD(sdgSig.doGrooming(jetCollectionSig));
129+
jetCollectionSigSD.addVector("zgSigSD", sdgSig.getZgs());
130+
jetCollectionSigSD.addVector("ndropSigSD", sdgSig.getNDroppedBranches());
131+
jetCollectionSigSD.addVector("dr12SigSD", sdgSig.getDR12());
132+
133+
//SoftDrop grooming classic for CS jets (zcut=0.1, beta=0)
134+
std::vector<jetCollection> jetCollectionCSSDs;
135+
136+
for(int ics = 0; ics<ncs; ++ics) {
137+
softDropGroomer sdgCS(0.1, 0.0, R);
138+
jetCollection jetCollectionCSSD(sdgCS.doGrooming(jetCollectionCSs[ics]));
139+
jetCollectionCSSD.addVector(Form("csJetSD_%dzg",ics), sdgCS.getZgs());
140+
jetCollectionCSSD.addVector(Form("csJetSD_%dndrop",ics), sdgCS.getNDroppedBranches());
141+
jetCollectionCSSD.addVector(Form("csJetSD_%ddr12",ics), sdgCS.getDR12());
142+
143+
jetCollectionCSSDs.push_back(jetCollectionCSSD);
144+
//std::cout << "n CS jets: " << jetCollectionCSs[ics].getJet().size() <<
145+
// " n CSSD jets: " << jetCollectionCSSDs[ics].getJet().size() << std::endl;
146+
}
147+
148+
149+
//match the CS jets to signal jets
150+
for(int ics = 0; ics<ncs; ++ics) {
151+
jetMatcher jmCS(R);
152+
jmCS.setBaseJets(jetCollectionCSs[ics]);
153+
jmCS.setTagJets(jetCollectionSig);
154+
jmCS.matchJets();
155+
156+
jmCS.reorderedToTag(jetCollectionCSs[ics]);
157+
jmCS.reorderedToTag(jetCollectionCSSDs[ics]);
158+
}
159+
160+
//---------------------------------------------------------------------------
161+
// write tree
162+
//---------------------------------------------------------------------------
163+
164+
//Give variable we want to write out to treeWriter.
165+
//Only vectors of the types 'jetCollection', and 'double', 'int', 'fastjet::PseudoJet' are supported
166+
167+
trw.addCollection("sigJet", jetCollectionSig);
168+
trw.addCollection("sigJetSD", jetCollectionSigSD);
169+
170+
for(int ics = 0; ics<ncs; ++ics) {
171+
trw.addCollection(Form("csJet_%d",ics), jetCollectionCSs[ics]);
172+
trw.addCollection(Form("csJetSD_%d",ics), jetCollectionCSSDs[ics]);
173+
}
174+
trw.addCollection("csRho", rho);
175+
trw.addCollection("csRhom", rhom);
176+
177+
trw.addCollection("eventWeight", eventWeight);
178+
179+
trw.fillTree();
180+
181+
}//event loop
182+
183+
Bar.Update(nEvent);
184+
Bar.Print();
185+
Bar.PrintLine();
186+
187+
TTree *trOut = trw.getTree();
188+
189+
TFile *fout = new TFile("JetToyHIResultCSVariations.root","RECREATE");
190+
trOut->Write();
191+
fout->Write();
192+
fout->Close();
193+
194+
double time_in_seconds = std::chrono::duration_cast<std::chrono::milliseconds>
195+
(std::chrono::steady_clock::now() - start_time).count() / 1000.0;
196+
std::cout << "runFromFile: " << time_in_seconds << std::endl;
197+
}

0 commit comments

Comments
 (0)