ATLAS Offline Software
Loading...
Searching...
No Matches
Classes | Namespaces | Functions
JetCalibTools_PlotJMSFactors.cxx File Reference
#include "JetCalibTools/JetCalibrationTool.h"
#include "xAODJet/Jet.h"
#include "xAODJet/JetContainer.h"
#include "xAODJet/JetAuxContainer.h"
#include "xAODRootAccess/TEvent.h"
#include "xAODRootAccess/TStore.h"
#include "AsgTools/StandaloneToolHandle.h"
#include "CxxUtils/checker_macros.h"
#include <vector>
#include "TString.h"
#include "TH2D.h"
#include "TCanvas.h"

Go to the source code of this file.

Classes

class  jet::JetFourMomAccessor
 JetFourMomAccessor is an extension of JetAttributeAccessor::AccessorWrapper<xAOD::JetFourMom_t> AccessorWrapper<xAOD::JetFourMom_t> purpose is to provide a direct and simple access to JetFourMom_t attributes (which are internally saved as 4 floats inside jets). More...

Namespaces

namespace  jet

Functions

double mTA (const double trackMass, const double trackPt, const double caloPt)
int main ATLAS_NOT_THREAD_SAFE (int argc, char *argv[])

Function Documentation

◆ ATLAS_NOT_THREAD_SAFE()

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

Definition at line 47 of file JetCalibTools_PlotJMSFactors.cxx.

48{
49 // Check argument usage
50 if (argc != 5 && argc != 6)
51 {
52 std::cout << "USAGE: " << argv[0] << " <JetCollection> <ConfigFile> <OutputFile> <mass type> (dev mode switch)" << std::endl;
53 std::cout << "\tMass types: \"calo\", \"TA\", \"comb\" (capitalization is ignored)" << std::endl;
54 return 1;
55 }
56
57 // Store the arguments
58 const TString jetAlgo = argv[1];
59 const TString config = argv[2];
60 const TString outFile = argv[3];
61 const TString massType = argv[4];
62 const bool isDevMode = argc > 5 && (TString(argv[5]) == "true" || TString(argv[5]) == "dev");
63
64 // Derived information
65 const bool outFileIsExtensible = outFile.EndsWith(".pdf") || outFile.EndsWith(".ps");
66 const bool outFileIsRoot = outFile.EndsWith(".root");
67 const bool doCombMass = !massType.CompareTo("comb",TString::kIgnoreCase);
68 const bool doCaloMass = doCombMass || !massType.CompareTo("calo",TString::kIgnoreCase);
69 const bool doTAMass = doCombMass || !massType.CompareTo("ta",TString::kIgnoreCase);
70
71 // Assumed constants
72 const TString calibSeq = "JMS"; // only want to apply the JMS here
73 const bool isData = false; // doesn't actually matter for JMS, which is always active
74 const float massForScan = 80.385e3; // W-boson
75 const float etaForScan = 0; // central jets
76 const float etaRange = (doCombMass || doTAMass) ? 2 : 3; // how far to go in |eta| (-etaRange,etaRange)
77 const float trackMassFactor = 2./3.; // Recall: mTA = ptCalo * (mTrack/ptTrack), this multiplies mTrack
78 const float trackPtFactor = 2./3.; // Recall: mTA = ptCalo * (mTrack/ptTrack), this multiplies ptTrack
79
80 // Accessor strings
81 const TString startingScaleString = "JetEtaJESScaleMomentum";
82 const TString caloMassScaleString = "JetJMSScaleMomentumCalo";
83 const TString taMassScaleString = "JetJMSScaleMomentumTA";
84 const TString taMassFloatString = "JetTrackAssistedMassCalibrated";
85 const TString detectorEtaString = "DetectorEta";
86 const TString trackMassString = "TrackSumMass";
87 const TString trackPtString = "TrackSumPt";
88
89 // Accessors
90 jet::JetFourMomAccessor startingScale(startingScaleString.Data());
91 jet::JetFourMomAccessor caloMassScale(caloMassScaleString.Data());
92 jet::JetFourMomAccessor taMassScale(taMassScaleString.Data());
93 SG::AuxElement::Accessor<float> mTAfloat(taMassFloatString.Data());
94 SG::AuxElement::Accessor<float> detectorEta(detectorEtaString.Data());
95 SG::AuxElement::Accessor<float> trackMass(trackMassString.Data());
96 SG::AuxElement::Accessor<float> trackPt(trackPtString.Data());
97
98
99 // Create the calib tool
100 asg::StandaloneToolHandle<IJetCalibrationTool> calibTool;
101 calibTool.setTypeAndName("JetCalibrationTool/MyJetCalibTool");
102 {
103 if (calibTool.initialize().isFailure())
104 {
105 std::cout << "Failed to make ana tool" << std::endl;
106 return 2;
107 }
108 if (calibTool.setProperty("JetCollection",jetAlgo.Data()).isFailure())
109 {
110 std::cout << "Failed to set JetCollection: " << jetAlgo.Data() << std::endl;
111 return 3;
112 }
113 if (calibTool.setProperty("ConfigFile",config.Data()).isFailure())
114 {
115 std::cout << "Failed to set ConfigFile: " << config.Data() << std::endl;
116 return 3;
117 }
118 if (calibTool.setProperty("CalibSequence",calibSeq.Data()).isFailure())
119 {
120 std::cout << "Failed to set CalibSequence: " << calibSeq.Data() << std::endl;
121 return 3;
122 }
123 if (calibTool.setProperty("IsData",isData).isFailure())
124 {
125 std::cout << "Failed to set IsData: " << (isData ? std::string("true") : std::string("false")) << std::endl;
126 return 3;
127 }
128 if (isDevMode && calibTool.setProperty("DEVmode",isDevMode).isFailure())
129 {
130 std::cout << "Failed to set DEVmode" << std::endl;
131 return 4;
132 }
133 if (calibTool.retrieve().isFailure())
134 {
135 std::cout << "Failed to initialize the JetCalibTool" << std::endl;
136 return 5;
137 }
138 }
139
140
141 // Build a jet container and a jet for us to manipulate later
142 xAOD::TEvent event;
143 xAOD::TStore store;
144 xAOD::JetContainer* jets = new xAOD::JetContainer();
145 jets->setStore(new xAOD::JetAuxContainer());
146 jets->push_back(new xAOD::Jet());
147 xAOD::Jet* jet = jets->at(0);
148
149 xAOD::JetContainer* calibJets = new xAOD::JetContainer();
150 calibJets->setStore(new xAOD::JetAuxContainer());
151 calibJets->push_back(new xAOD::Jet());
152 xAOD::Jet* calibJet = calibJets->at(0);
153
154 // Make the histograms to fill
155 std::vector<TH2D*> hists_pt_eta;
156 std::vector<TH2D*> hists_pt_mpt;
157 if (doCaloMass)
158 {
159 hists_pt_eta.push_back(new TH2D("JMS_calo_pt_eta",Form("JMS (calo) for jets with mass=%.1f GeV",massForScan/1.e3),1150,200,2500,20*etaRange,-etaRange,etaRange));
160 hists_pt_mpt.push_back(new TH2D("JMS_calo_pt_mpt",Form("JMS (calo): for jets with #eta=%.1f",etaForScan),1150,200,2500,100,0,1));
161 }
162 if (doTAMass)
163 {
164 hists_pt_eta.push_back(new TH2D("JMS_TA_pt_eta",Form("JMS (TA) for jets with mass=%.1f GeV (TA factor = %.1f)",massForScan/1.e3,trackMassFactor/trackPtFactor),1150,200,2500,20*etaRange,-etaRange,etaRange));
165 hists_pt_mpt.push_back(new TH2D("JMS_TA_pt_mpt",Form("JMS (TA) for jets with #eta=%.1f",etaForScan),1150,200,2500,100,0,1));
166 }
167 if (doCombMass)
168 {
169 hists_pt_eta.push_back(new TH2D("JMS_comb_pt_eta",Form("JMS (comb) for jets with mass=%.1f GeV",massForScan/1.e3),1150,200,2500,20*etaRange,-etaRange,etaRange));
170 hists_pt_mpt.push_back(new TH2D("JMS_comb_pt_mpt",Form("JMS (comb) for jets with eta=%.1f",etaForScan),1150,200,2500,100,0,1));
171 }
172
173 // Fill the pt vs eta histogram
174 for (int xBin = 1; xBin <= hists_pt_eta.at(0)->GetNbinsX(); ++xBin)
175 {
176 const double pt = hists_pt_eta.at(0)->GetXaxis()->GetBinCenter(xBin)*1.e3;
177 for (int yBin = 1; yBin <= hists_pt_eta.at(0)->GetNbinsY(); ++yBin)
178 {
179 const double eta = hists_pt_eta.at(0)->GetYaxis()->GetBinCenter(yBin);
180
181 // Set the main 4-vector and scale 4-vector
182 jet->setJetP4(xAOD::JetFourMom_t(pt,eta,0,massForScan));
183 detectorEta(*jet) = eta;
184 trackMass(*jet) = trackMassFactor*massForScan;
185 trackPt(*jet) = trackPtFactor*pt;
186 startingScale.setAttribute(*jet,xAOD::JetFourMom_t(pt,eta,0,massForScan));
187 // Repeat for jet to calibrate
188 calibJet->setJetP4(xAOD::JetFourMom_t(pt,eta,0,massForScan));
189 detectorEta(*calibJet) = eta;
190 trackMass(*calibJet) = trackMassFactor*massForScan;
191 trackPt(*calibJet) = trackPtFactor*pt;
192 startingScale.setAttribute(*calibJet,xAOD::JetFourMom_t(pt,eta,0,massForScan));
193
194 // Jet kinematics set, now apply calibration
195 if (calibTool->modify(*calibJets).isFailure()) {
196 std::cout << "Failed to apply jet calibration" << std::endl;
197 return 6;
198 }
199
200 // Calculate the scale factors
201 const double JMS = calibJet->m()/startingScale(*jet).mass();
202 const double JMScalo = doCombMass ? caloMassScale(*calibJet).mass()/startingScale(*jet).mass() : JMS;
203 const double JMSTA = (doCombMass ? taMassScale(*calibJet).mass() : mTAfloat(*calibJet))/mTA(trackMass(*jet),trackPt(*jet),startingScale(*jet).pt());
204
205 // JMS retrieved, fill the plot(s)
206 size_t plotIndex = 0;
207 if (doCaloMass)
208 hists_pt_eta.at(plotIndex++)->SetBinContent(xBin,yBin,JMScalo);
209 if (doTAMass)
210 hists_pt_eta.at(plotIndex++)->SetBinContent(xBin,yBin,JMSTA);
211 if (doCombMass)
212 hists_pt_eta.at(plotIndex++)->SetBinContent(xBin,yBin,JMS);
213 }
214 }
215
216 // Fill the pt vs m/pt histogram
217 for (int xBin = 1; xBin <= hists_pt_mpt.at(0)->GetNbinsX(); ++xBin)
218 {
219 const double pt = hists_pt_mpt.at(0)->GetXaxis()->GetBinCenter(xBin)*1.e3;
220 for (int yBin = 1; yBin <= hists_pt_mpt.at(0)->GetNbinsY(); ++yBin)
221 {
222 const double mpt = hists_pt_mpt.at(0)->GetYaxis()->GetBinCenter(yBin);
223 const double mass = pt*mpt;
224
225 // Set the main 4-vector and scale 4-vector
226 jet->setJetP4(xAOD::JetFourMom_t(pt,etaForScan,0,mass));
227 detectorEta(*jet) = etaForScan;
228 trackMass(*jet) = trackMassFactor*mass;
229 trackPt(*jet) = trackPtFactor*pt;
230 startingScale.setAttribute(*jet,xAOD::JetFourMom_t(pt,etaForScan,0,mass));
231 // Repeat for jet to calibrate
232 calibJet->setJetP4(xAOD::JetFourMom_t(pt,etaForScan,0,mass));
233 detectorEta(*calibJet) = etaForScan;
234 trackMass(*calibJet) = trackMassFactor*mass;
235 trackPt(*calibJet) = trackPtFactor*pt;
236 startingScale.setAttribute(*calibJet,xAOD::JetFourMom_t(pt,etaForScan,0,mass));
237
238 // Jet kinematics set, now apply calibration
239 if(calibTool->modify(*calibJets).isFailure()){
240 std::cout << "Failed to apply jet calibration" << std::endl;
241 return 6;
242 }
243
244 // Calculate the scale factors
245 const double JMS = calibJet->m()/startingScale(*jet).mass();
246 const double JMScalo = doCombMass ? caloMassScale(*calibJet).mass()/startingScale(*jet).mass() : JMS;
247 const double JMSTA = (doCombMass ? taMassScale(*calibJet).mass() : mTAfloat(*calibJet))/mTA(trackMass(*jet),trackPt(*jet),startingScale(*jet).pt());
248
249 // JMS retrieved, fill the plot(s)
250 size_t plotIndex = 0;
251 if (doCaloMass)
252 hists_pt_mpt.at(plotIndex++)->SetBinContent(xBin,yBin,JMScalo);
253 if (doTAMass)
254 hists_pt_mpt.at(plotIndex++)->SetBinContent(xBin,yBin,JMSTA);
255 if (doCombMass)
256 hists_pt_mpt.at(plotIndex++)->SetBinContent(xBin,yBin,JMS);
257 }
258 }
259
260
261 // Make the plots
262
263 // First the canvas
264 TCanvas* canvas = new TCanvas("canvas");
265 canvas->SetMargin(0.07,0.13,0.1,0.10);
266 canvas->SetFillStyle(4000);
267 canvas->SetFillColor(0);
268 canvas->SetFrameBorderMode(0);
269 canvas->cd();
270 canvas->SetLogx(true);
271
272 // Now labels/etc
273 for (TH2* hist : hists_pt_eta)
274 {
275 hist->SetStats(false);
276 hist->GetXaxis()->SetTitle("Jet #it{p}_{T} [GeV]");
277 hist->GetXaxis()->SetTitleOffset(1.35);
278 hist->GetXaxis()->SetMoreLogLabels();
279 hist->GetYaxis()->SetTitle("#eta");
280 hist->GetYaxis()->SetTitleOffset(0.9);
281 hist->GetZaxis()->SetTitle("m_{JES+JMS} / m_{JES}");
282 //hist->GetZaxis()->SetRangeUser(0.4,1.1);
283 }
284 for (TH2* hist : hists_pt_mpt)
285 {
286 hist->SetStats(false);
287 hist->GetXaxis()->SetTitle("Jet #it{p}_{T} [GeV]");
288 hist->GetXaxis()->SetTitleOffset(1.35);
289 hist->GetXaxis()->SetMoreLogLabels();
290 hist->GetYaxis()->SetTitle("m / #it{p}_{T}");
291 hist->GetYaxis()->SetTitleOffset(0.9);
292 hist->GetZaxis()->SetTitle("m_{JES+JMS} / m_{JES}");
293 //hist->GetZaxis()->SetRangeUser(0.4,1.1);
294 }
295
296 // Now write them out
297 if (outFileIsExtensible)
298 {
299 canvas->Print(outFile+"[");
300 for (size_t iHist = 0; iHist < hists_pt_eta.size(); ++iHist)
301 {
302 hists_pt_eta.at(iHist)->Draw("colz");
303 canvas->Print(outFile);
304 hists_pt_mpt.at(iHist)->Draw("colz");
305 canvas->Print(outFile);
306 }
307 canvas->Print(outFile+"]");
308 }
309 else if (outFileIsRoot)
310 {
311 TFile* outRootFile = new TFile(outFile,"RECREATE");
312 std::cout << "Writing to output ROOT file: " << outFile << std::endl;
313 outRootFile->cd();
314 for (size_t iHist = 0; iHist < hists_pt_eta.size(); ++iHist)
315 {
316 hists_pt_eta.at(iHist)->Write();
317 hists_pt_mpt.at(iHist)->Write();
318 }
319 outRootFile->Close();
320 }
321 else
322 {
323 unsigned int counter = 1;
324 for (size_t iHist = 0; iHist < hists_pt_eta.size(); ++iHist)
325 {
326 hists_pt_eta.at(iHist)->Draw("colz");
327 canvas->Print(Form("%u-fixMass-%s",counter,outFile.Data()));
328 hists_pt_mpt.at(iHist)->Draw("colz");
329 canvas->Print(Form("%u-fixMass-%s",counter++,outFile.Data()));
330 }
331 }
332
333 // Clean up a bit (although the program is about to end...)
334 for (TH2D* hist : hists_pt_eta)
335 delete hist;
336 hists_pt_eta.clear();
337
338 for (TH2D* hist : hists_pt_mpt)
339 delete hist;
340 hists_pt_mpt.clear();
341
342 return 0;
343}
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
mTA
double mTA(const double trackMass, const double trackPt, const double caloPt)
Definition JetCalibTools_PlotJMSFactors.cxx:40
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.
SG::AuxElement::Accessor
SG::Accessor< T, ALLOC > Accessor
Definition AuxElement.h:572
asg::StandaloneToolHandle
an "initializing" ToolHandle for stand-alone applications
Definition StandaloneToolHandle.h:44
asg::StandaloneToolHandle::setProperty
StatusCode setProperty(const std::string &name, T2 &&value)
Definition StandaloneToolHandle.h:105
asg::StandaloneToolHandle::retrieve
StatusCode retrieve()
initialize the tool, will succeed if the tool was already initialized
Definition StandaloneToolHandle.h:147
asg::StandaloneToolHandle::setTypeAndName
void setTypeAndName(const std::string &typeAndName)
Definition StandaloneToolHandle.h:101
asg::StandaloneToolHandle::initialize
StatusCode initialize()
initialize the tool, will fail if the tool was already initialized
Definition StandaloneToolHandle.h:158
jet::JetFourMomAccessor
JetFourMomAccessor is an extension of JetAttributeAccessor::AccessorWrapper<xAOD::JetFourMom_t> Acces...
Definition JetCalibTools_PlotJESFactors.cxx:32
string
STL class.
xAOD::Jet_v1::setJetP4
void setJetP4(const JetFourMom_t &p4)
Definition Jet_v1.cxx:171
xAOD::Jet_v1::m
virtual double m() const
The invariant mass of the particle.
Definition Jet_v1.cxx:59
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
DQPostProcessTest.outFile
outFile
Comment Out Those You do not wish to run.
Definition DQPostProcessTest.py:36
LArCellBinning.etaRange
etaRange
Filling Eta range.
Definition LArCellBinning.py:128
SGTest::store
TestStore store
Definition TestStore.cxx:23
TauTrackVars::trackPt
bool trackPt(const xAOD::TauJet &, const xAOD::TauTrack &track, float &out)
Definition ConstituentLoaderTauTrack.cxx:113
std
STL namespace.
test_pyathena.counter
counter
Definition test_pyathena.py:15
test_pyathena.pt
pt
Definition test_pyathena.py:11
xAOD::Jet
Jet_v1 Jet
Definition of the current "jet version".
Definition Event/xAOD/xAODJet/xAODJet/Jet.h:17
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
xAOD::JetFourMom_t
ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< double > > JetFourMom_t
Base 4 Momentum type for Jet.
Definition JetTypes.h:17

◆ mTA()

double mTA ( const double trackMass,
const double trackPt,
const double caloPt )

Definition at line 40 of file JetCalibTools_PlotJMSFactors.cxx.

41{
42 return trackMass * caloPt / trackPt;
43}
Generated on for ATLAS Offline Software by doxygen 1.14.0