JetMET Corrections

[lcorcodilos]

This thread is meant to start organizing the JME correction implementation as discussed at the meeting today (Jan 7th, 2021). We agreed to prioritize implementing whatever we can from CMSSW so that we can limit our work to "building the payload" in an RDataFrame friendly way.

@dmajumder please feel free to drop links to your old code here. Even just some include lines and a few lines that call methods from those headers would be useful in pointing us in the right direction.

[dmajumder]

Hi @lcorcodilos @mroguljic here's my implementation of JEC, JEC uncertainties, and JER.
https://github.com/dmajumder/VLQAna/blob/CMSSW_8_0_X_NewB2GAnaFW/src/JetMaker.cc

The main CMSSW classes are:
CondFormats/JetMETObjects/interface/JetCorrectorParameters.h
CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h
CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h

At the moment, I suggested that Matej implements the JEC uncertainties and the JER smearing and uncertainties. With nanoAODv7, we do not need to update the JEC itself.

The key parts of the JEC uncertainty implementation are:

• Creating the object to retrieve JECs from a payload: https://github.com/dmajumder/VLQAna/blob/CMSSW_8_0_X_NewB2GAnaFW/src/JetMaker.cc#L137
• Evaluating JEC uncertainty for a particular jet: https://github.com/dmajumder/VLQAna/blob/CMSSW_8_0_X_NewB2GAnaFW/src/JetMaker.cc#L265-L268

As for the JER, the values are here, taken from the JER TWiki:
https://github.com/dmajumder/VLQAna/blob/CMSSW_8_0_X_NewB2GAnaFW/interface/ApplyJER.h
numbers were hard-coded but this can be a nice lookup table somewhere.
The JER are applied to the jets here using the hybrid method:
https://github.com/dmajumder/VLQAna/blob/CMSSW_8_0_X_NewB2GAnaFW/src/JetMaker.cc#L236-L249

There's a lot going on in the above codes besides this and I can walk you through the main parts if needed.
You will see that the above is very similar tot the NanoAODTools implementation.

[mroguljic]

Hi everyone,

I started writing a class to handle JES shifts, following Devdatta's examples. The first goal was to get an object which can read uncertainties, given a jet, and to use this object to see if it returns sensible numbers. A very simple class is defined here:

• https://github.com/mroguljic/TIMBER/blob/master/TIMBER/Framework/src/JESUncShifter.cc [1]
• https://github.com/mroguljic/TIMBER/blob/master/TIMBER/Framework/include/JESUncShifter.h [2]

Constructor takes the payload name and creates an object to read it and there is a method which uses this object and returns the uncertainty for a given pt and eta.

I tried using this class with a simple TIMBER script

• https://github.com/mroguljic/X_YH_4b/blob/master/JECtesting.py#L16 [3]
  and calling a function which would read the uncertainty for a given jet
• https://github.com/mroguljic/TIMBER/blob/master/TIMBER/Framework/JESShiftTesting.cc [4]
  For the moment, dummy values of eta and pT are used instead of values for given jets.

When running [3], warning messages appear, indicating that jet eta wasn't set (but it was on L13 of [1]), thus giving an invalid uncertainty value of -999

%MSG-e JetCorrectionUncertainty::: 
 jet eta is not set
%MSG
%MSG-e SimpleJetCorrectionUncertainty: 
 bin variables out of range

Furthermore, there is also this error, which I guess is causing the problem:
Error in <TCling__LazyFunctionCreatorAutoloadForModule>: Failed to load library /afs/cern.ch/work/m/mrogulji/X_YH_4b/CMSSW_11_1_5/lib/slc7_amd64_gcc820/libCondFormatsJetMETObjects.so

Since I'm working on lxplus, I'm using CMSSW_11_1_5 per these instructions https://github.com/lcorcodilos/TIMBER/blob/master/doxysetup/Python3.md with the virtualenv. In order to use CMSSW software with TIMBER when doing JEC magic, from the CMSSW area, I did this:

git cms-init
git cms-addpkg CondFormats/JetMETObjects
scram b -j4

Perhaps there is something wrong with my virtualenv and/or the CMSSW setup, maybe someone else could try running this?

[mroguljic]

From what I gather, this could also be solved by using a standalone JEC reader class instead of trying to compile CMSSW and integrate it with TIMBER. Maybe we could ask JetMET or NanoAODTools team (or someone else?) if standalone libraries exist and are ready to be used. I see two benefits:

1. TIMBER doesn't need to rely on CMSSW
2. No python version restrictions and workarounds

Maybe this is a better long-term solution?

[lcorcodilos]

I agree but I think the key is to get things working in CMSSW (and thus, python 2.7) first so that we understand what is needed.

Then we could write a small utility that knows how to check for the scripts we've grabbed as standalone for a given CMSSW. I'm thinking something along the lines of python updateFromCMSSW.py --version=11_0_1. Then a user can update to the CMSSW version they want and we can keep the master branch of TIMBER updated with some widely used CMSSW version.

Of course, the issue is that we don't want 100s of files just to cover dependencies from CMSSW which I think gives us more reason to figure out what we need first.

[lcorcodilos]

@mroguljic I'm not sure if you eventually realized this but I was getting a seg fault because the .txt file had empty lines at the end and JetCorrectorParameters.cc didn't know how to handle it. Sort of silly but valuable info nonetheless!

Though, I also re-organized a bit your testing code. Primarily, I got rid of the testJESShifter step which you use for the Define because then the txt file will be opened for every single event and we just want to open it once. I'll paste below what I did instead.

If you've made any other progress, let me know. I'll work on this a bit more today and share my progress afterwards here.

# JECtesting.py
import ROOT
...
CompileCpp('JESUncShifter shifter = JESUncShifter("Summer19UL18_V5_MC_Uncertainty_AK8PFPuppi.txt");')
a.Define("JESUnc",'shifter.GetUncVal(0.12, 347.81, 1)') 
...

[mroguljic]

I haven't realized it, but this explains why the seg fault went away with changing the .txt file!

The problem with opening the file every time was bugging me, and I couldn't think of the solution. I didn't even realize that one could instantiate the shifter object and use it with RDataFrame, so that's solved.

[lcorcodilos]

This is actually exactly why the different modules have been built as classes 😄. The Correction() class does the instantiating external to the dataframe (basically runs the CompileCpp line) and then analyzer.AddCorrection() does the part of building the per-event eval via a Define() call. In other words, we want Correction() and AddCorrection() to eventually replace these two lines.

[lcorcodilos]

I'd like to just summarize the to-do list of what is needed and where to find the info on it since I've been digging around for a bit now. I can edit these if something looks incorrect (or more needs to be added).

• JES uncertainty. The CMSSW interfacing is essentially implemented now by JESUncShifter. We just need a way to package it in TIMBER and I think the most obvious would be through the Correction class. The issue is that the Correction class is implicitely for event weights and not "calibrations". However, it does everything that a Calibration() class would do at the moment. I have to think a bit more about it but what we want to avoid is tricking someone into thinking they can do a.AddCorrection(myJEScorrection) and get the right thing out. Maybe it makes sense to rename AddCorrection() to AddEventWeight(). Down the line, a dedicated Calibration() class could be used to fork the processing so you have three different sets of nodes, all with different, but that are acted on in parallel when you do a Define or Cut.
• JES re-correct the jet. This seems simple enough that I think we should just do it now while we have a sense of what's going on 😄
• JER uncertainty. Algorithm described here. Requires a dedicated smearer (they don't provide it for some reason...) The good news is that I already converted this from NanoAOD-tools - though it wasn't exactly planned... here. This just needs to be tested then.
• JMS uncertainty. Probably the easiest one. Just needs these flat values hardcoded.
• JMR uncertainty. Uses the same smearing procedure as JER (again, I already coded it up in C++ but it's untested).
• PUPPI SD mass correction. Low priority but there should be a module that deals with this. The logic in this is non-negligible. IIRC, one needs the uncorrected soft drop subjets to reconstruct the AK8 soft drop mass on which to apply this correction. More details here, though pending Run2 recommendations.

[lcorcodilos]

There infrastructural changes that I think should be made but only after we have things tested and working.

• ExternalTools -> ext
• Proper header files in include/ to match classes in src/ (it's otherwise a mess right now)
• Documentation page for modules intended for users to use (since there will be a mix of forward facing and internal .cc and .h files)

[lcorcodilos]

Actually decided to do the first of these two ahead of time but they will be lumped in with the JetMET modules.

[lcorcodilos]

There has been quite a bit of unforeseen work to get some of the infrastructure in place for these calibrations/uncertainties to work as we'd like. I've been doing the work on the JMEdev branch and will chronicle some of the work here along with a to-do of other items that need to be considered but aren't directly the modules themselves.

Done Already (directly JME related)

• Write JetRecalibrator class which handles the interfacing to the CMSSW based tools.
• Write JMEpaths class which handles the interfacing to the JME txt files
• Write JES_weight class which is the user-facing module to access the corrections (including recalibrations) and uncertainties by global tag, jet flavor, and uncertainty type ("" for the total and "" for a specific correction type if one would like to factorize the corrections as separate parameters).

Still needs to be done (directly JME related)

• Edit/rewrite JetSmearer using knowledge learned from JetRecalibrator.
• Write JER_weight.
• Write JMR_weight.
• Write JMS_weight.

Done Already (indirectly JME related)

• Moved TIMBER/Framework/ExternalTools/to TIMBER/Framework/ext/
• Organized TIMBER/Framework/src and TIMBER/Framework/include so that declarations and implementations are split (roughly - there are still some outstanding where they make sense).
• Added ability to create collections inside the RDataFrame. Can be done in one swoop with analyzer() option createAllCollections or manually with CreateCollection() or CreateAllCollections() but will be done on the fly as the user needs them otherwise. This is done by reading the Define or Cut action and parsing for +"s". Taking the example of "FatJet" (branches named FatJet_pt, FatJet_eta, etc), the collection FatJets will be built. TIMBER is doing two things - (1) writing and compiling the code for a struct that describes all FatJet attributes and (2) creating the actual object (vector of FatJet structs) that can be "stored" in the RDataFrame per-event.
  • The reason for this addition is so the new JME modules can take in the jet structs directly, thus simplifying the C++ code (and the actual use of the methods in the code).
  • This has penalties! Because there is more being compiled and built, memory usage increases (though not beyond values that are reasonable). It also increases processing time substantially. The difference between calculating the mean of FatJet_pt[0] and FatJets[0].pt is nearly a factor of 10 6 (with some help from the internet). If this is not the dominant process in the per-event loop, the differences are much smaller (factor of ~3 in one test where I used the collections everywhere I could). The collections should be reserved for using as input to C++ modules that benefit from having fewer positional arguments. There is a more efficient way to handle the construction of the collections (see issue Optimize SofA collections with DOD #36 ) but this is where I leave it for now.
  • Added ability to silence Define and Cut calls (so CreateCollection() is in the background in the terminal).
• Kept RunChain as attribute of analyzer so that contents can be easily accessed.
• Added dedicated method ReorderCollection(). This is meant to be used when JECs affect the pt of jets and a re-ordering is needed.
• Added CalibrateVar() method which will actually do the multiplication of a variable by the calibration weight (looping through uncertainty variations as well).
  • NOTE: This means a module feeds back {nominal, nominal+up, nominal-down}!
• Added Node.GetBaseNode() so that the top most parent can be accessed from a Node (ie. outside the analyzer).
• Added ModuleWorker class to handle all of the functionality shared by Correction and the new Calibration (both of which now inherit from ModuleWorker).
• Added Calibration class which doesn't do much different from ModuleWorker at the moment (it's just not Correction).
• To common.h, added TempDir (handles temporary directory storage) and ReadTarFile (opens and streams tarball contents - took a day to get working!) This was necessary because the JECs come in tarballs and untarring 800+ files and holding them in TIMBER is undesirable.
  • This required adding libarchive as a dependency. Added to setup.sh the recipe (if it doesn't exist) to download and build it inside of TIMBER/bin
• Added TIMBER/Utilities/CollectionGen.py which is used to build the C++ strings to create the collection structs.

Still needs to be done but lower priority (indirectly JME related)

• Need to add a function that returns a new ordering of pt vectors where the pt has been reweighted (ie. return indices of vector items sorted in descending order).
• HTML scraper to identify, download, and maintain (ie. update if "old") the JEC files from the Twiki