ATLAS Offline Software
Loading...
Searching...
No Matches
Functions
testResolution.cxx File Reference
#include "JetUncertainties/JetUncertaintiesTool.h"
#include "JetUncertainties/Helpers.h"
#include "xAODJet/Jet.h"
#include "xAODJet/JetContainer.h"
#include "xAODJet/JetAuxContainer.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODEventInfo/EventInfoContainer.h"
#include "xAODEventInfo/EventInfoAuxContainer.h"
#include "PATInterfaces/SystematicRegistry.h"
#include "PATInterfaces/SystematicVariation.h"
#include "PATInterfaces/SystematicsUtil.h"
#include "OptionHelper.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TLatex.h"

Go to the source code of this file.

Functions

void setJetKinematics (xAOD::Jet &jet, double pt, double eta, double phi, double mass)
int main (int argc, char *argv[])

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 50 of file testResolution.cxx.

51{
52 StatusCode::enableFailure();
53 jet::OptionHelper optHelper;
54
55
56 if (argc != 7 && argc != 8)
57 {
58 printf("Expected arguments:\n");
59 printf("\t1. Output file (pdf)\n");
60 printf("\t2. Jet definition\n");
61 printf("\t3. MC type\n");
62 printf("\t4. Config file\n");
63 printf("\t5. Component name and variation\n");
64 printf("\t\tExample: \"FourVecResUnc,+1\"\n");
65 printf("\t6. IsData (\"true\" or \"false\")\n");
66 printf("\t7. Options (optional argument), semi-colon delimited, examples:\n");
67 printf("\t\tisDijet=false\n");
68 printf("\t\tisLargeR=false\n");
69 exit(1);
70 }
71 TString outFile = argv[1];
72 TString jetDef = argv[2];
73 TString mcType = argv[3];
74 TString config = argv[4];
75 TString component = argv[5];
76 TString isDataStr = argv[6];
77 try{
78 if (argc == 8) optHelper.Initialize(jet::utils::vectorize<TString>(argv[7],";"));
79 else optHelper.Initialize(std::vector<TString>());
80 } catch (const std::exception& e) {
81 std::cerr<<"Exception in optHelper.Initialize: " << e.what()<<std::endl;
82 exit(1);
83 } catch (...) {
84 std::cerr<<"Unknown exception in optHelper.Initialize"<<std::endl;
85 exit(1);
86 }
87
88 if (!outFile.EndsWith(".pdf"))
89 {
90 printf("Only pdf output files are currently supported\n");
91 exit(1);
92 }
93
94 if (jet::utils::vectorize<TString>(component,",").size() != 2)
95 {
96 printf("Bad component formatting, got \"%s\"\n",component.Data());
97 exit(1);
98 }
99 TString compName = jet::utils::vectorize<TString>(component,",").at(0);
100 float shift = atof(jet::utils::vectorize<TString>(component,",").at(1));
101
102 bool isData = false;
103 if (!isDataStr.CompareTo("true",TString::kIgnoreCase))
104 isData = true;
105 else if (!isDataStr.CompareTo("false",TString::kIgnoreCase))
106 isData = false;
107 else
108 {
109 printf("Unable to determine whether to configure for data or MC: %s\n",isDataStr.Data());
110 exit(1);
111 }
112
113 JetUncertaintiesTool* tool = new JetUncertaintiesTool("JESprov");
114 if (tool->setProperty("JetDefinition",jetDef.Data()).isFailure())
115 exit(1);
116 if (tool->setProperty("MCType",mcType.Data()).isFailure())
117 exit(1);
118 if (tool->setProperty("ConfigFile",config.Data()).isFailure())
119 exit(1);
120 if (tool->setProperty("IsData",isData).isFailure())
121 exit(1);
122 if (optHelper.GetPath() != "")
123 if (tool->setProperty("Path",optHelper.GetPath().Data()).isFailure())
124 exit(1);
125 if (tool->setScaleToGeV().isFailure())
126 exit(1);
127 if (tool->initialize().isFailure())
128 exit(1);
129
130
131 // Build a jet container and a jet for us to manipulate later
132 xAOD::TEvent event;
133 xAOD::TStore store;
134 xAOD::JetContainer* jets = new xAOD::JetContainer();
135 jets->setStore(new xAOD::JetAuxContainer());
136 jets->push_back(new xAOD::Jet());
137 xAOD::Jet* jet = jets->at(0);
138
139 // Build an EventInfo object for us to manipulate later
140 xAOD::EventInfoContainer* eInfos = new xAOD::EventInfoContainer();
141 eInfos->setStore(new xAOD::EventInfoAuxContainer());
142 eInfos->push_back(new xAOD::EventInfo());
143 xAOD::EventInfo* eInfo = eInfos->at(0);
144
145 // Ensure that the specified component is a valid systematic
146 // This is a +1sigma variation
147 CP::SystematicVariation variation(Form("JET_%s",compName.Data()),shift);
148 if (!tool->isAffectedBySystematic(variation))
149 {
150 printf("The specified variation was not recognized: JET_%s\n",compName.Data());
151 exit(1);
152 }
153 CP::SystematicSet syst;
154 syst.insert(variation);
155 if (tool->applySystematicVariation(syst) != StatusCode::SUCCESS)
156 {
157 printf("Failed to apply systematic variation\n");
158 exit(1);
159 }
160
161 // Get info on the variation
162 const size_t compIndex = tool->getComponentIndex("JET_"+compName);
163 const std::set<jet::CompScaleVar::TypeEnum> scaleVars = tool->getComponentScaleVars(compIndex);
164 const jet::CompScaleVar::TypeEnum scaleVar = scaleVars.size() == 1 ? *(scaleVars.begin()) : jet::CompScaleVar::UNKNOWN;
165 const jet::JetTopology::TypeEnum topology = tool->getComponentTopology(compIndex);
166
167 printf("Component is index %zu, with ScaleVar %s and topology %s\n",compIndex,jet::CompScaleVar::enumToString(scaleVar).Data(),jet::JetTopology::enumToString(topology).Data());
168
169 // Prepare the canvas
170 TCanvas canvas("canvas");
171 canvas.SetMargin(0.12,0.04,0.15,0.04);
172 canvas.SetFillStyle(4000);
173 canvas.SetFillColor(0);
174 canvas.SetFrameBorderMode(0);
175 canvas.cd();
176 canvas.Print(outFile+"[");
177
178
179 // Prepare to add labels
180 TLatex tex;
181 tex.SetNDC();
182 tex.SetTextFont(42);
183 tex.SetTextSize(0.04);
184
185
186 // Define the jet that we want to smear repeatedly
187 // Center the jet around pT of 1 TeV and mass of 100 GeV (easy numbers)
188 // eta_bins = [0.0,0.2,0.7,1.3,1.8,2.5,3.2,3.5,4.5]
189 const double pT = 1000;
190 const double mass = 100;
191 const double phi = 0;
192 const std::vector<double> etaVals = optHelper.GetFixedEtaVals();
193 const int numSmear = 100000;
194
195 for (const double eta : etaVals)
196 {
197 // Prepare a histogram to fill repeatedly as we smear the value
198 TH1D smearPt(Form("smearPt_eta%.2f",eta),"",100,500,1500);
199 TH1D smearMass(Form("smearMass_eta%.2f",eta),"",100,50,150);
200 smearPt.Sumw2();
201 smearMass.Sumw2();
202 smearPt.SetStats(0);
203 smearMass.SetStats(0);
204 smearPt.GetXaxis()->SetTitle("#it{p}_{T} [GeV]");
205 smearMass.GetXaxis()->SetTitle("#it{m} [GeV]");
206 smearPt.GetYaxis()->SetTitle("Number of events");
207 smearMass.GetYaxis()->SetTitle("Number of events");
208
209 // Now smear repeatedly from the same starting jet
210 for (int iSmear = 0; iSmear < numSmear; ++iSmear)
211 {
212 setJetKinematics(*jet,pT,eta,phi,mass);
213 if (tool->applyCorrection(*jet,*eInfo) != CP::CorrectionCode::Ok)
214 {
215 printf("Error while smearing, iteration %d\n",iSmear);
216 exit(1);
217 }
218 // Fill the histograms with the smeared jet
219 smearPt.Fill(jet->pt());
220 smearMass.Fill(jet->m());
221 }
222
223 // Get info on the expected values
224 setJetKinematics(*jet,pT,eta,phi,mass);
225 const double nomData = tool->getNominalResolutionData(*jet,scaleVar,topology);
226 const double nomMC = tool->getNominalResolutionMC(*jet,scaleVar,topology);
227 const double uncert = tool->getUncertainty(compIndex,*jet,*eInfo,scaleVar);
228
229 const double fullUncMC = jet::CompScaleVar::isRelResolutionType(scaleVar) ? nomMC * uncert : uncert;
230 const double fullUncData = jet::CompScaleVar::isRelResolutionType(scaleVar) ? nomData * uncert : uncert;
231 const double smearMC = sqrt(pow(nomMC + fabs(shift*fullUncMC),2) - pow(nomMC,2));
232 const double smearData = nomData != JESUNC_ERROR_CODE ? sqrt(pow(nomData + fabs(shift*fullUncData),2) - pow(nomData,2)) : 0;
233
234 // Fit a Gaussian
235 TF1 fitPt("fitPt","gaus");
236 smearPt.Fit(&fitPt,"E");
237 TF1 fitMass("fitMass","gaus");
238 smearMass.Fit(&fitMass,"E");
239
240 // Draw and print the plot
241 smearPt.Draw();
242 tex.DrawLatex(0.70,0.9,"Gaussian fit results");
243 tex.DrawLatex(0.70,0.85,Form("#mu = %.0f GeV",fitPt.GetParameter(1)));
244 tex.DrawLatex(0.70,0.80,Form("#sigma = %.1f GeV",fitPt.GetParameter(2)));
245 tex.DrawLatex(0.70,0.75,Form("#mu/#sigma = %.1f%%",fitPt.GetParameter(2)/fitPt.GetParameter(1)*100));
246 tex.DrawLatex(0.15,0.9,"Expectation from tool");
247 tex.DrawLatex(0.15,0.85,"#sigma_{smear}^{2} = (#sigma_{nom} + |N#delta#sigma|)^{2} - (#sigma_{nom})^{2}");
248 if (nomData != JESUNC_ERROR_CODE)
249 tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = %.1f%%",100*nomData));
250 else
251 tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = N/A"));
252 tex.DrawLatex(0.15,0.75,Form("#sigma_{nom}^{MC}/#it{p}_{T} = %.1f%%",100*nomMC));
253 tex.DrawLatex(0.15,0.70,Form("#delta#sigma/#it{p}_{T} = %.1f%%",100*uncert));
254 tex.DrawLatex(0.15,0.65,Form("N_{sigma} = %+.1f",shift));
255 if (nomData != JESUNC_ERROR_CODE)
256 tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = %.1f%%",100*smearData));
257 else
258 tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = N/A"));
259 tex.DrawLatex(0.15,0.50,Form("#sigma_{smear}^{MC}/#it{p}_{T} = %.1f%%",100*smearMC));
260 canvas.Print(outFile);
261
262
263
264 smearMass.Draw();
265 tex.DrawLatex(0.70,0.9,"Gaussian fit results");
266 tex.DrawLatex(0.70,0.85,Form("#mu = %.0f GeV",fitMass.GetParameter(1)));
267 tex.DrawLatex(0.70,0.80,Form("#sigma = %.1f GeV",fitMass.GetParameter(2)));
268 tex.DrawLatex(0.70,0.75,Form("#mu/#sigma = %.1f%%",fitMass.GetParameter(2)/fitMass.GetParameter(1)*100));
269 tex.DrawLatex(0.15,0.9,"Expectation from tool");
270 tex.DrawLatex(0.15,0.85,"#sigma_{smear}^{2} = (#sigma_{nom} + |N#delta#sigma|)^{2} - (#sigma_{nom})^{2}");
271 if (nomData != JESUNC_ERROR_CODE)
272 tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = %.1f%%",100*nomData));
273 else
274 tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = N/A"));
275 tex.DrawLatex(0.15,0.75,Form("#sigma_{nom}^{MC}/#it{p}_{T} = %.1f%%",100*nomMC));
276 tex.DrawLatex(0.15,0.70,Form("#delta#sigma/#it{p}_{T} = %.1f%%",100*uncert));
277 tex.DrawLatex(0.15,0.65,Form("N_{sigma} = %+.1f",shift));
278 if (nomData != JESUNC_ERROR_CODE)
279 tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = %.1f%%",100*smearData));
280 else
281 tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = N/A"));
282 tex.DrawLatex(0.15,0.50,Form("#sigma_{smear}^{MC}/#it{p}_{T} = %.1f%%",100*smearMC));
283 canvas.Print(outFile);
284 }
285
286 // Done printing
287 canvas.Print(outFile+"]");
288
289 return 0;
290}
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
Data
@ Data
Definition BaseObject.h:11
optHelper
jet::OptionHelper optHelper
Definition MakeUncertaintyPlots.cxx:52
CP::CorrectionCode::Ok
@ Ok
The correction was done successfully.
Definition CorrectionCode.h:38
CP::SystematicSet
Class to wrap a set of SystematicVariations.
Definition SystematicSet.h:31
CP::SystematicSet::insert
void insert(const SystematicVariation &systematic)
description: insert a systematic into the set
Definition SystematicSet.cxx:87
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
xAOD::TEvent
Tool for accessing xAOD files outside of Athena.
Definition Control/xAODRootAccess/xAODRootAccess/TEvent.h:59
xAOD::TStore
A relatively simple transient store for objects created in analysis.
Definition TStore.h:45
Base_Fragment.mass
mass
Definition Sherpa_i/share/common/Base_Fragment.py:59
CxxUtils::atof
double atof(std::string_view str)
Converts a string into a double / float.
Definition Control/CxxUtils/Root/StringUtils.cxx:91
DQPostProcessTest.outFile
outFile
Comment Out Those You do not wish to run.
Definition DQPostProcessTest.py:36
SGTest::store
TestStore store
Definition TestStore.cxx:23
conifer::pow
constexpr int pow(int x)
Definition conifer.h:27
jet::CompScaleVar::isRelResolutionType
bool isRelResolutionType(const TypeEnum type)
Definition UncertaintyEnum.cxx:379
jet::CompScaleVar::enumToString
TString enumToString(const TypeEnum type)
Definition UncertaintyEnum.cxx:250
jet::JetTopology::enumToString
TString enumToString(const TypeEnum type)
Definition UncertaintyEnum.cxx:609
jet::utils::vectorize
bool vectorize(const TString &str, const TString &sep, std::vector< T > &result)
Definition Reconstruction/Jet/JetUncertainties/JetUncertainties/Helpers.h:147
python.base_data.config
config
Definition base_data.py:20
xAOD::Jet
Jet_v1 Jet
Definition of the current "jet version".
Definition Event/xAOD/xAODJet/xAODJet/Jet.h:17
xAOD::EventInfoContainer
EventInfoContainer_v1 EventInfoContainer
Define the latest version of the container.
Definition EventInfoContainer.h:17
xAOD::EventInfo
EventInfo_v1 EventInfo
Definition of the latest event info version.
Definition IEventInfoCnvTool.h:18
xAOD::EventInfoAuxContainer
EventInfoAuxContainer_v1 EventInfoAuxContainer
Define the latest version of the auxiliary container.
Definition EventInfoAuxContainer.h:16
xAOD::JetAuxContainer
JetAuxContainer_v1 JetAuxContainer
Definition of the current jet auxiliary container.
Definition JetAuxContainer.h:22
xAOD::JetContainer
JetContainer_v1 JetContainer
Definition of the current "jet container version".
Definition JetContainer.h:17
setJetKinematics
void setJetKinematics(xAOD::Jet &jet, double pt, double eta, double phi, double mass)
Definition testResolution.cxx:32

◆ setJetKinematics()

void setJetKinematics ( xAOD::Jet & jet,
double pt,
double eta,
double phi,
double mass )

Definition at line 32 of file testResolution.cxx.

33{
34 static jet::JetFourMomAccessor scaleCalo(jet::CompMassDef::getJetScaleString(jet::CompMassDef::CaloMass).Data());
35 static jet::JetFourMomAccessor scaleTA(jet::CompMassDef::getJetScaleString(jet::CompMassDef::TAMass).Data());
36 static jet::JetFourMomAccessor scaleCombQCD(jet::CompMassDef::getJetScaleString(jet::CompMassDef::CombMassQCD).Data());
37 static jet::JetFourMomAccessor scaleCombWZ(jet::CompMassDef::getJetScaleString(jet::CompMassDef::CombMassWZ).Data());
38 static jet::JetFourMomAccessor scaleCombHbb(jet::CompMassDef::getJetScaleString(jet::CompMassDef::CombMassHbb).Data());
39 static jet::JetFourMomAccessor scaleCombTop(jet::CompMassDef::getJetScaleString(jet::CompMassDef::CombMassTop).Data());
40
41 jet.setJetP4( xAOD::JetFourMom_t(pt,eta,phi,mass));
42 scaleCalo.setAttribute( jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
43 scaleTA.setAttribute( jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
44 scaleCombQCD.setAttribute(jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
45 scaleCombWZ.setAttribute( jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
46 scaleCombHbb.setAttribute(jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
47 scaleCombTop.setAttribute(jet, xAOD::JetFourMom_t(pt,eta,phi,mass));
48}
jet::JetFourMomAccessor
JetFourMomAccessor is an extension of JetAttributeAccessor::AccessorWrapper<xAOD::JetFourMom_t> Acces...
Definition JetCalibTools_PlotJESFactors.cxx:32
jet::CompMassDef::getJetScaleString
TString getJetScaleString(const TypeEnum type)
Definition UncertaintyEnum.cxx:225
xAOD::JetFourMom_t
ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< double > > JetFourMom_t
Base 4 Momentum type for Jet.
Definition JetTypes.h:17
Generated on for ATLAS Offline Software by doxygen 1.14.0