ATLAS Offline Software
testResolution.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "JetUncertainties/JetUncertaintiesTool.h"
6 #include "JetUncertainties/Helpers.h"
7 
8 #include "xAODJet/Jet.h"
9 #include "xAODJet/JetContainer.h"
10 #include "xAODJet/JetAuxContainer.h"
11 #include "xAODEventInfo/EventInfo.h"
12 #include "xAODEventInfo/EventInfoContainer.h"
13 #include "xAODEventInfo/EventInfoAuxContainer.h"
14 
15 #include "PATInterfaces/SystematicRegistry.h"
16 #include "PATInterfaces/SystematicVariation.h"
17 #include "PATInterfaces/SystematicsUtil.h"
18 
19 #ifdef ROOTCORE
20 # include "xAODRootAccess/TEvent.h"
21 # include "xAODRootAccess/TStore.h"
22 #endif // ROOTCORE
23 
24 #include "OptionHelper.h"
25 
26 #include "TH1D.h"
27 #include "TCanvas.h"
28 #include "TF1.h"
29 #include "TLatex.h"
30 
31 
32 void setJetKinematics(xAOD::Jet& jet, double pt, double eta, double phi, double mass)
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 }
49 
50 int main (int argc, char* argv[])
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  if (argc == 8) optHelper.Initialize(jet::utils::vectorize<TString>(argv[7],";"));
78  else optHelper.Initialize(std::vector<TString>());
79 
80  if (!outFile.EndsWith(".pdf"))
81  {
82  printf("Only pdf output files are currently supported\n");
83  exit(1);
84  }
85 
86  if (jet::utils::vectorize<TString>(component,",").size() != 2)
87  {
88  printf("Bad component formatting, got \"%s\"\n",component.Data());
89  exit(1);
90  }
91  TString compName = jet::utils::vectorize<TString>(component,",").at(0);
92  float shift = atof(jet::utils::vectorize<TString>(component,",").at(1));
93 
94  bool isData = false;
95  if (!isDataStr.CompareTo("true",TString::kIgnoreCase))
96  isData = true;
97  else if (!isDataStr.CompareTo("false",TString::kIgnoreCase))
98  isData = false;
99  else
100  {
101  printf("Unable to determine whether to configure for data or MC: %s\n",isDataStr.Data());
102  exit(1);
103  }
104 
105  JetUncertaintiesTool* tool = new JetUncertaintiesTool("JESprov");
106  if (tool->setProperty("JetDefinition",jetDef.Data()).isFailure())
107  exit(1);
108  if (tool->setProperty("MCType",mcType.Data()).isFailure())
109  exit(1);
110  if (tool->setProperty("ConfigFile",config.Data()).isFailure())
111  exit(1);
112  if (tool->setProperty("IsData",isData).isFailure())
113  exit(1);
114  if (optHelper.GetPath() != "")
115  if (tool->setProperty("Path",optHelper.GetPath().Data()).isFailure())
116  exit(1);
117  if (tool->setScaleToGeV().isFailure())
118  exit(1);
119  if (tool->initialize().isFailure())
120  exit(1);
121 
122 
123  // Build a jet container and a jet for us to manipulate later
124  xAOD::TEvent event;
125  xAOD::TStore store;
126  xAOD::JetContainer* jets = new xAOD::JetContainer();
127  jets->setStore(new xAOD::JetAuxContainer());
128  jets->push_back(new xAOD::Jet());
129  xAOD::Jet* jet = jets->at(0);
130 
131  // Build an EventInfo object for us to manipulate later
132  xAOD::EventInfoContainer* eInfos = new xAOD::EventInfoContainer();
133  eInfos->setStore(new xAOD::EventInfoAuxContainer());
134  eInfos->push_back(new xAOD::EventInfo());
135  xAOD::EventInfo* eInfo = eInfos->at(0);
136 
137  // Ensure that the specified component is a valid systematic
138  // This is a +1sigma variation
139  CP::SystematicVariation variation(Form("JET_%s",compName.Data()),shift);
140  if (!tool->isAffectedBySystematic(variation))
141  {
142  printf("The specified variation was not recognized: JET_%s\n",compName.Data());
143  exit(1);
144  }
145  CP::SystematicSet syst;
146  syst.insert(variation);
147  if (tool->applySystematicVariation(syst) != StatusCode::SUCCESS)
148  {
149  printf("Failed to apply systematic variation\n");
150  exit(1);
151  }
152 
153  // Get info on the variation
154  const size_t compIndex = tool->getComponentIndex("JET_"+compName);
155  const std::set<jet::CompScaleVar::TypeEnum> scaleVars = tool->getComponentScaleVars(compIndex);
156  const jet::CompScaleVar::TypeEnum scaleVar = scaleVars.size() == 1 ? *(scaleVars.begin()) : jet::CompScaleVar::UNKNOWN;
157  const jet::JetTopology::TypeEnum topology = tool->getComponentTopology(compIndex);
158 
159  printf("Component is index %zu, with ScaleVar %s and topology %s\n",compIndex,jet::CompScaleVar::enumToString(scaleVar).Data(),jet::JetTopology::enumToString(topology).Data());
160 
161  // Prepare the canvas
162  TCanvas canvas("canvas");
163  canvas.SetMargin(0.12,0.04,0.15,0.04);
164  canvas.SetFillStyle(4000);
165  canvas.SetFillColor(0);
166  canvas.SetFrameBorderMode(0);
167  canvas.cd();
168  canvas.Print(outFile+"[");
169 
170 
171  // Prepare to add labels
172  TLatex tex;
173  tex.SetNDC();
174  tex.SetTextFont(42);
175  tex.SetTextSize(0.04);
176 
177 
178  // Define the jet that we want to smear repeatedly
179  // Center the jet around pT of 1 TeV and mass of 100 GeV (easy numbers)
180  // eta_bins = [0.0,0.2,0.7,1.3,1.8,2.5,3.2,3.5,4.5]
181  const double pT = 1000;
182  const double mass = 100;
183  const double phi = 0;
184  const std::vector<double> etaVals = optHelper.GetFixedEtaVals();
185  const int numSmear = 100000;
186 
187  for (const double eta : etaVals)
188  {
189  // Prepare a histogram to fill repeatedly as we smear the value
190  TH1D smearPt(Form("smearPt_eta%.2f",eta),"",100,500,1500);
191  TH1D smearMass(Form("smearMass_eta%.2f",eta),"",100,50,150);
192  smearPt.Sumw2();
193  smearMass.Sumw2();
194  smearPt.SetStats(0);
195  smearMass.SetStats(0);
196  smearPt.GetXaxis()->SetTitle("#it{p}_{T} [GeV]");
197  smearMass.GetXaxis()->SetTitle("#it{m} [GeV]");
198  smearPt.GetYaxis()->SetTitle("Number of events");
199  smearMass.GetYaxis()->SetTitle("Number of events");
200 
201  // Now smear repeatedly from the same starting jet
202  for (int iSmear = 0; iSmear < numSmear; ++iSmear)
203  {
204  setJetKinematics(*jet,pT,eta,phi,mass);
205  if (tool->applyCorrection(*jet,*eInfo) != CP::CorrectionCode::Ok)
206  {
207  printf("Error while smearing, iteration %d\n",iSmear);
208  exit(1);
209  }
210  // Fill the histograms with the smeared jet
211  smearPt.Fill(jet->pt());
212  smearMass.Fill(jet->m());
213  }
214 
215  // Get info on the expected values
216  setJetKinematics(*jet,pT,eta,phi,mass);
217  const double nomData = tool->getNominalResolutionData(*jet,scaleVar,topology);
218  const double nomMC = tool->getNominalResolutionMC(*jet,scaleVar,topology);
219  const double uncert = tool->getUncertainty(compIndex,*jet,*eInfo,scaleVar);
220 
221  const double fullUncMC = jet::CompScaleVar::isRelResolutionType(scaleVar) ? nomMC * uncert : uncert;
222  const double fullUncData = jet::CompScaleVar::isRelResolutionType(scaleVar) ? nomData * uncert : uncert;
223  const double smearMC = sqrt(pow(nomMC + fabs(shift*fullUncMC),2) - pow(nomMC,2));
224  const double smearData = nomData != JESUNC_ERROR_CODE ? sqrt(pow(nomData + fabs(shift*fullUncData),2) - pow(nomData,2)) : 0;
225 
226  // Fit a Gaussian
227  TF1 fitPt("fitPt","gaus");
228  smearPt.Fit(&fitPt,"E");
229  TF1 fitMass("fitMass","gaus");
230  smearMass.Fit(&fitMass,"E");
231 
232  // Draw and print the plot
233  smearPt.Draw();
234  tex.DrawLatex(0.70,0.9,"Gaussian fit results");
235  tex.DrawLatex(0.70,0.85,Form("#mu = %.0f GeV",fitPt.GetParameter(1)));
236  tex.DrawLatex(0.70,0.80,Form("#sigma = %.1f GeV",fitPt.GetParameter(2)));
237  tex.DrawLatex(0.70,0.75,Form("#mu/#sigma = %.1f%%",fitPt.GetParameter(2)/fitPt.GetParameter(1)*100));
238  tex.DrawLatex(0.15,0.9,"Expectation from tool");
239  tex.DrawLatex(0.15,0.85,"#sigma_{smear}^{2} = (#sigma_{nom} + |N#delta#sigma|)^{2} - (#sigma_{nom})^{2}");
240  if (nomData != JESUNC_ERROR_CODE)
241  tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = %.1f%%",100*nomData));
242  else
243  tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = N/A"));
244  tex.DrawLatex(0.15,0.75,Form("#sigma_{nom}^{MC}/#it{p}_{T} = %.1f%%",100*nomMC));
245  tex.DrawLatex(0.15,0.70,Form("#delta#sigma/#it{p}_{T} = %.1f%%",100*uncert));
246  tex.DrawLatex(0.15,0.65,Form("N_{sigma} = %+.1f",shift));
247  if (nomData != JESUNC_ERROR_CODE)
248  tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = %.1f%%",100*smearData));
249  else
250  tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = N/A"));
251  tex.DrawLatex(0.15,0.50,Form("#sigma_{smear}^{MC}/#it{p}_{T} = %.1f%%",100*smearMC));
252  canvas.Print(outFile);
253 
254 
255 
256  smearMass.Draw();
257  tex.DrawLatex(0.70,0.9,"Gaussian fit results");
258  tex.DrawLatex(0.70,0.85,Form("#mu = %.0f GeV",fitMass.GetParameter(1)));
259  tex.DrawLatex(0.70,0.80,Form("#sigma = %.1f GeV",fitMass.GetParameter(2)));
260  tex.DrawLatex(0.70,0.75,Form("#mu/#sigma = %.1f%%",fitMass.GetParameter(2)/fitMass.GetParameter(1)*100));
261  tex.DrawLatex(0.15,0.9,"Expectation from tool");
262  tex.DrawLatex(0.15,0.85,"#sigma_{smear}^{2} = (#sigma_{nom} + |N#delta#sigma|)^{2} - (#sigma_{nom})^{2}");
263  if (nomData != JESUNC_ERROR_CODE)
264  tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = %.1f%%",100*nomData));
265  else
266  tex.DrawLatex(0.15,0.80,Form("#sigma_{nom}^{data}/#it{p}_{T} = N/A"));
267  tex.DrawLatex(0.15,0.75,Form("#sigma_{nom}^{MC}/#it{p}_{T} = %.1f%%",100*nomMC));
268  tex.DrawLatex(0.15,0.70,Form("#delta#sigma/#it{p}_{T} = %.1f%%",100*uncert));
269  tex.DrawLatex(0.15,0.65,Form("N_{sigma} = %+.1f",shift));
270  if (nomData != JESUNC_ERROR_CODE)
271  tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = %.1f%%",100*smearData));
272  else
273  tex.DrawLatex(0.15,0.55,Form("#sigma_{smear}^{data}/#it{p}_{T} = N/A"));
274  tex.DrawLatex(0.15,0.50,Form("#sigma_{smear}^{MC}/#it{p}_{T} = %.1f%%",100*smearMC));
275  canvas.Print(outFile);
276  }
277 
278  // Done printing
279  canvas.Print(outFile+"]");
280 
281  return 0;
282 }
jet::CompMassDef::CombMassTop
@ CombMassTop
Definition: UncertaintyEnum.h:79
SGTest::store
TestStore store
Definition: TestStore.cxx:23
CalculateHighPtTerm.pT
pT
Definition: ICHEP2016/CalculateHighPtTerm.py:57
Jet.h
AddEmptyComponent.compName
compName
Definition: AddEmptyComponent.py:32
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
xAOD::JetAttributeAccessor::AccessorWrapper::setAttribute
void setAttribute(SG::AuxElement &p, const TYPE &v) const
Definition: JetAccessors.h:54
Data
@ Data
Definition: BaseObject.h:11
EventInfoAuxContainer.h
jet::CompMassDef::CaloMass
@ CaloMass
Definition: UncertaintyEnum.h:74
test_pyathena.pt
pt
Definition: test_pyathena.py:11
CP::SystematicSet
Class to wrap a set of SystematicVariations.
Definition: SystematicSet.h:31
xAOD::EventInfoContainer
EventInfoContainer_v1 EventInfoContainer
Define the latest version of the container.
Definition: EventInfoContainer.h:17
OptionHelper.h
jet::OptionHelper::GetFixedEtaVals
std::vector< double > GetFixedEtaVals() const
Definition: OptionHelper.h:550
CP::SystematicVariation
Definition: SystematicVariation.h:47
jet::CompScaleVar::isRelResolutionType
bool isRelResolutionType(const TypeEnum type)
Definition: UncertaintyEnum.cxx:379
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
EventInfoContainer.h
Helpers.h
jet::CompMassDef::CombMassHbb
@ CombMassHbb
Definition: UncertaintyEnum.h:78
JESUNC_ERROR_CODE
#define JESUNC_ERROR_CODE
Definition: Reconstruction/Jet/JetUncertainties/JetUncertainties/Helpers.h:23
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
jet::CompScaleVar::enumToString
TString enumToString(const TypeEnum type)
Definition: UncertaintyEnum.cxx:250
CalculateHighPtTerm.jetDef
dictionary jetDef
Definition: ICHEP2016/CalculateHighPtTerm.py:28
xAOD::JetAuxContainer_v1
Temporary container used until we have I/O for AuxStoreInternal.
Definition: JetAuxContainer_v1.h:37
jet::CompMassDef::CombMassWZ
@ CombMassWZ
Definition: UncertaintyEnum.h:77
jet
Definition: JetCalibTools_PlotJESFactors.cxx:23
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
jet::CompMassDef::getJetScaleString
TString getJetScaleString(const TypeEnum type)
Definition: UncertaintyEnum.cxx:225
jet::CompMassDef::TAMass
@ TAMass
Definition: UncertaintyEnum.h:75
SystematicRegistry.h
JetUncertaintiesTool
Definition: JetUncertaintiesTool.h:44
LArCellNtuple.argv
argv
Definition: LArCellNtuple.py:152
DeMoAtlasDataLoss.canvas
dictionary canvas
Definition: DeMoAtlasDataLoss.py:187
TEvent.h
jet::JetTopology::enumToString
TString enumToString(const TypeEnum type)
Definition: UncertaintyEnum.cxx:609
jet::OptionHelper::Initialize
bool Initialize(const std::vector< TString > &options)
Definition: OptionHelper.h:242
DQHistogramMergeRegExp.argc
argc
Definition: DQHistogramMergeRegExp.py:20
jet::JetFourMomAccessor
JetFourMomAccessor is an extension of JetAttributeAccessor::AccessorWrapper<xAOD::JetFourMom_t> Acces...
Definition: JetCalibTools_PlotJESFactors.cxx:32
jet::OptionHelper::GetPath
TString GetPath() const
Definition: OptionHelper.h:40
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
calibdata.exit
exit
Definition: calibdata.py:236
jet::CompScaleVar::UNKNOWN
@ UNKNOWN
Definition: UncertaintyEnum.h:93
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:37
xAOD::EventInfoAuxContainer_v1
Auxiliary information about the pileup events.
Definition: EventInfoAuxContainer_v1.h:31
xAOD::JetFourMom_t
ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< double > > JetFourMom_t
Base 4 Momentum type for Jet.
Definition: JetTypes.h:17
CP::SystematicSet::insert
void insert(const SystematicVariation &systematic)
description: insert a systematic into the set
Definition: SystematicSet.cxx:88
AtlCoolConsole.tool
tool
Definition: AtlCoolConsole.py:453
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
jet::OptionHelper
Definition: OptionHelper.h:23
EventInfo.h
xAOD::EventInfo_v1
Class describing the basic event information.
Definition: EventInfo_v1.h:43
setJetKinematics
void setJetKinematics(xAOD::Jet &jet, double pt, double eta, double phi, double mass)
Definition: testResolution.cxx:32
jet::CompScaleVar::TypeEnum
TypeEnum
Definition: UncertaintyEnum.h:91
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
CP::CorrectionCode::Ok
@ Ok
The correction was done successfully.
Definition: CorrectionCode.h:38
xAOD::TStore
A relatively simple transient store for objects created in analysis.
Definition: TStore.h:44
JetContainer.h
optHelper
jet::OptionHelper optHelper
Definition: MakeUncertaintyPlots.cxx:52
main
int main(int argc, char *argv[])
Definition: testResolution.cxx:50
jet::CompMassDef::CombMassQCD
@ CombMassQCD
Definition: UncertaintyEnum.h:76
JetAuxContainer.h
defineDB.jets
list jets
Definition: JetTagCalibration/share/defineDB.py:24
xAOD::JetContainer
JetContainer_v1 JetContainer
Definition of the current "jet container version".
Definition: JetContainer.h:17
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
jet::JetTopology::TypeEnum
TypeEnum
Definition: UncertaintyEnum.h:208
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
SystematicVariation.h
TStore.h
xAOD::TEvent
Tool for accessing xAOD files outside of Athena.
Definition: Control/xAODRootAccess/xAODRootAccess/TEvent.h:84
JetUncertaintiesTool.h
SystematicsUtil.h
Generated on Sun Dec 22 2024 21:19:14 for ATLAS Offline Software by   doxygen 1.8.18