ATLAS Offline Software
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

 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
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
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 }

◆ 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 }
SGTest::store
TestStore store
Definition: TestStore.cxx:23
mTA
double mTA(const double trackMass, const double trackPt, const double caloPt)
Definition: JetCalibTools_PlotJMSFactors.cxx:40
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:68
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
LArCellBinning.etaRange
etaRange
Filling Eta range.
Definition: LArCellBinning.py:128
muonBucketDump.outRootFile
outRootFile
Definition: muonBucketDump.py:7
plotmaker.hist
hist
Definition: plotmaker.py:148
test_pyathena.pt
pt
Definition: test_pyathena.py:11
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
xAOD::Jet_v1::setJetP4
void setJetP4(const JetFourMom_t &p4)
Definition: Jet_v1.cxx:171
xAOD::JetAuxContainer_v1
Temporary container used until we have I/O for AuxStoreInternal.
Definition: JetAuxContainer_v1.h:37
jet
Definition: JetCalibTools_PlotJESFactors.cxx:23
asg::StandaloneToolHandle::initialize
StatusCode initialize()
initialize the tool, will fail if the tool was already initialized
Definition: StandaloneToolHandle.h:158
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
asg::StandaloneToolHandle::setProperty
StatusCode setProperty(const std::string &name, T2 &&value)
Definition: StandaloneToolHandle.h:105
asg::StandaloneToolHandle
an "initializing" ToolHandle for stand-alone applications
Definition: StandaloneToolHandle.h:44
LArCellNtuple.argv
argv
Definition: LArCellNtuple.py:152
DeMoAtlasDataLoss.canvas
dictionary canvas
Definition: DeMoAtlasDataLoss.py:187
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
IJetCalibrationTool::modify
virtual StatusCode modify(xAOD::JetContainer &jets) const override final
Apply calibration to a jet container (for IJetModifier interface).
Definition: IJetCalibrationTool.h:33
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
asg::StandaloneToolHandle::retrieve
StatusCode retrieve()
initialize the tool, will succeed if the tool was already initialized
Definition: StandaloneToolHandle.h:147
TauGNNUtils::Variables::Track::trackPt
bool trackPt(const xAOD::TauJet &, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:472
DQPostProcessTest.outFile
outFile
Comment Out Those You do not wish to run.
Definition: DQPostProcessTest.py:37
xAOD::JetFourMom_t
ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< double > > JetFourMom_t
Base 4 Momentum type for Jet.
Definition: JetTypes.h:17
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
asg::StandaloneToolHandle::setTypeAndName
void setTypeAndName(const std::string &typeAndName)
Definition: StandaloneToolHandle.h:101
xAOD::Jet_v1
Class describing a jet.
Definition: Jet_v1.h:57
xAOD::TStore
A relatively simple transient store for objects created in analysis.
Definition: TStore.h:44
xAOD::Jet_v1::m
virtual double m() const
The invariant mass of the particle.
Definition: Jet_v1.cxx:59
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.
test_pyathena.counter
counter
Definition: test_pyathena.py:15
xAOD::TEvent
Tool for accessing xAOD files outside of Athena.
Definition: Control/xAODRootAccess/xAODRootAccess/TEvent.h:84