Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 14 additions & 9 deletions NanoAnalysis/python/LHEAngProbFiller.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ def __init__(self, MELA, NANOVERSION, MELASettings = None):
self.MELASettings = MELASettings
self.sortedSettings = []
self.ProbHelper = MELAProbHelper(self.MELA, self.MELASettings, "LHE")
print("***LHEAngProbFiller: set for: ", self.ProbHelper.names, flush=True)

# if MELASettings != None:
# defaults = {
Expand Down Expand Up @@ -72,9 +73,9 @@ def __init__(self, MELA, NANOVERSION, MELASettings = None):

def beginFile(self, inputFile, outputFile, inputTree, wrappedOutputTree):
self.out = wrappedOutputTree
self.out.branch("LHEMela_qH", "F", title="The mass of the Higgs candidate as reconstructed by the 4-leptons at LHE-level.")
self.out.branch("LHEMela_mZ1", "F", title="The mass of the first decay particle as reconstructed by 2 of the LHE-level leptons.")
self.out.branch("LHEMela_mZ2", "F", title="The mass of the second decay particle as reconstructed by 2 of the LHE-level leptons.")
# self.out.branch("LHEMela_qH", "F", title="The mass of the Higgs candidate as reconstructed by the 4-leptons at LHE-level.")
# self.out.branch("LHEMela_mZ1", "F", title="The mass of the first decay particle as reconstructed by 2 of the LHE-level leptons.")
# self.out.branch("LHEMela_mZ2", "F", title="The mass of the second decay particle as reconstructed by 2 of the LHE-level leptons.")
self.out.branch("LHEMela_costheta1", "F", limitedPrecision=16, title="In the Higgs' rest frame, theta_1 is the angle between the momentum of Z1 and the momentum of one of its decay products.")
self.out.branch("LHEMela_costheta2", "F", limitedPrecision=16, title="In the Higgs' rest frame, theta_2 is the angle between the momentum of Z2 and the momentum of one of its decay products.")
self.out.branch("LHEMela_Phi", "F", limitedPrecision=16, title="In the Higgs' rest frame, phi is the angle between the planes formed by the decay products of the two Z bosons.")
Expand Down Expand Up @@ -156,17 +157,21 @@ def analyze(self, event):
# self.MELA.setInputEvent(daughters, associated, mothers, 1)
self.MELA.setInputEvent(daughters, None, None, 0)
qH, mZ1, mZ2, costheta1, costheta2, Phi, costhetastar, Phi1 = self.MELA.computeDecayAngles()
self.out.fillBranch("LHEMela_costheta1", costheta1)
self.out.fillBranch("LHEMela_costheta2", costheta2)
self.out.fillBranch("LHEMela_Phi", Phi)
self.out.fillBranch("LHEMela_Phi1", Phi1)
self.out.fillBranch("LHEMela_costhetastar", costhetastar)
else:
qH, mZ1, mZ2, costheta1, costheta2, Phi, costhetastar, Phi1 = 0.,0.,0.,-999.,-999.,-999.,-999.,-999.
if len(daughters.toList()) != 4:
print(f"WARNING: LHEAngProbFiller: {len(daughters.toList())} LHE-leptons were selected for this event (4 expected)!")
elif abs(hMass - daughters.MTotal()) < 0.01:
else :
print(f"WARNING: LHEAngProbFiller: The invariant mass of the four LHE-leptons, {daughters.MTotal()}, is too different from the mass of the LHE-Higgs {hMass}! Expected a difference of less than 0.01, obtained a difference of ", hMass - daughters.MTotal())

# self.out.fillBranch("LHEMela_qH", qH)
# self.out.fillBranch("LHEMela_mZ1", mZ1)
# self.out.fillBranch("LHEMela_mZ2", mZ2)
self.out.fillBranch("LHEMela_costheta2", costheta2)
self.out.fillBranch("LHEMela_Phi", Phi)
self.out.fillBranch("LHEMela_Phi1", Phi1)
self.out.fillBranch("LHEMela_costhetastar", costhetastar)

self.MELA.resetInputEvent()

self.ProbHelper.fillProbs([daughters], [associated], [mothers])
Expand Down
139 changes: 69 additions & 70 deletions NanoAnalysis/python/MELAProbHelper.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,10 @@ def __init__(self, MELA, MELASettings, ModuleContext):
self.MELA = MELA
self.sortedSettings = []
self.ModuleContext = ModuleContext
self.names = []

if ModuleContext not in ["LHE", "Reco"] :
raise valueError("MELAProbHelper: invalid ModuleContext", ModuleContext)

if MELASettings != None:
defaults = {
Expand All @@ -27,7 +31,6 @@ def __init__(self, MELA, MELASettings, ModuleContext):
"Prod": None,
"Dec": None,
"Couplings": None,
# "isgen": None,
"context": None,
"computeprop": None,
"propscheme": "FixedWidth",
Expand All @@ -44,9 +47,12 @@ def __init__(self, MELA, MELASettings, ModuleContext):
}


### Split into reco and LHE probs
### Split into reco and LHE probs
for prob in MELASettings:
if (prob["context"] != self.ModuleContext) and (prob["context"] != "Any"): continue
if prob["context"] not in ["LHE", "Reco", "Any"] :
raise valueError(f"MELAProbHelper: invalid context {prob['context']} for {prob['Name']}")

if (prob["context"] != ModuleContext) and (prob["context"] != "Any"): continue

### Merge specific settings with defaults and check for unsupported values
fprob=copy.deepcopy(defaults)
Expand All @@ -69,50 +75,42 @@ def __init__(self, MELA, MELASettings, ModuleContext):
self.sortedSettings.append(fprob)

### Add index of probability to be used for dividep.
names = [d["Name"] for d in self.sortedSettings]
self.names = [d["Name"] for d in self.sortedSettings]
for prob in self.sortedSettings:
dp = prob["dividep"]
prob["dividep_eval"] = (None if dp == None else f"aCand.P_{dp}") # string to be evaluated to extract denominator
# print("***MELAPROBHELPER: ", prob["dividep_eval"])
prob["dividep_idx"] = (-1 if dp == None else names.index(dp)) # prob index, more efficient but would require keeping probs for all cands

print(f"***MELAProbHelper: probs: {names}", flush=True)
if dp == None :
prob["dividep_idx"] = -1
elif ModuleContext != "LHE" :
raise(ValueError(f"MELAProbHelper: for {prob['Name']}: 'dividep' is supported only for 'context=LHE'"))
else:
prob["dividep_idx"] = self.names.index(dp) # Index of probability to be used as denominator.

def bookProbs(self, wrappedOutputTree):
#this needs a per-module lenVar.
self.out = wrappedOutputTree
# print("***MELAHELPER: ", len(self.sortedSettings))
if len(self.sortedSettings) !=0 :
for p, prob in enumerate(self.sortedSettings):
print(prob["branchname"])
if self.ModuleContext == "LHE":
self.out.branch(prob["branchname"], "F", lenVar = None, limitedPrecision=16, title="User-defined LHE-level probability")
# print("***MELAProbHelper: wrote out ", prob["branchname"])
if prob.get("ispm4l", False):
self.out.branch(prob["branchname"]+"_ScaleUp", "F", lenVar = None, limitedPrecision=16, title="User-defined LHE-level m4l probability with Scale uncertainties up")
self.out.branch(prob["branchname"]+"_ScaleDown", "F", lenVar = None, limitedPrecision=16, title="User-defined LHE-level m4l probability with Scale uncertainties down")
self.out.branch(prob["branchname"]+"_SystUp", "F", lenVar = None, limitedPrecision=16, title="User-defined LHE-level m4l probability with Systematic uncertainties up")
self.out.branch(prob["branchname"]+"_SystDown", "F", lenVar = None, limitedPrecision=16, title="User-defined LHE-level m4l probability with Systematic uncertainties down")
if prob["computeprop"]:
self.out.branch(prob["branchname"]+"_prop", "F", limitedPrecision=16, title="User-defined weight to translate from POWHEG complex propagator scheme to JHUGen Breit-Wigner scheme")

else:
# print("***MELAProbHelper", prob["branchname"])
self.out.branch(prob["branchname"], "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined Reco-level probability")
# print("***MELAProbHelper: wrote out ", prob["branchname"])
if prob.get("ispm4l", False):
self.out.branch(prob["branchname"]+"_ScaleUp", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined Reco-level m4l probability with Scale uncertainties up")
self.out.branch(prob["branchname"]+"_ScaleDown", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined Reco-level m4l probability with Scale uncertainties down")
self.out.branch(prob["branchname"]+"_SystUp", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined Reco-level m4l probability with Systematic uncertainties up")
self.out.branch(prob["branchname"]+"_SystDown", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined Reco-level m4l probability with Systematic uncertainties down")
if prob["computeprop"]:
self.out.branch(prob["branchname"]+"_prop", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined weight to translate from POWHEG complex propagator scheme to JHUGen Breit-Wigner scheme")
if prob["addPAux"]:
self.out.branch(prob["branchname"]+"_aux", "F", lenVar = "nZZCand", limitedPrecision=16, title="User-defined auxiliary probability")
if prob["addPmavjj"]:
self.out.branch(prob["branchname"]+"_mavjj", "F", lenVar = "nZZCand", title="User-defined mavjj probability")
if prob["addPmavjj_true"]:
self.out.branch(prob["branchname"]+"_mvajj_true", "F", lenVar = "nZZCand", title="User-defined mvajj_true probability")
# print(prob["branchname"])
if self.ModuleContext == "LHE":
lenVar_ = None
else:
lenVar_ = "nZZCand"

self.out.branch(prob["branchname"], "F", lenVar=lenVar_, limitedPrecision=16, title=f"User-defined {self.ModuleContext}-level probability")
if prob.get("ispm4l", False):
self.out.branch(prob["branchname"]+"_ScaleUp", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined {self.ModuleContext}-level m4l probability with Scale uncertainties up")
self.out.branch(prob["branchname"]+"_ScaleDown", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined {self.ModuleContext}-level m4l probability with Scale uncertainties down")
self.out.branch(prob["branchname"]+"_SystUp", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined {self.ModuleContext}-level m4l probability with Systematic uncertainties up")
self.out.branch(prob["branchname"]+"_SystDown", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined {self.ModuleContext}-level m4l probability with Systematic uncertainties down")
if prob["computeprop"]:
self.out.branch(prob["branchname"]+"_prop", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined weight to translate from POWHEG complex propagator scheme to JHUGen Breit-Wigner scheme")
if prob["addPAux"]:
self.out.branch(prob["branchname"]+"_aux", "F", lenVar=lenVar_, limitedPrecision=16, title="User-defined auxiliary probability")
if prob["addPmavjj"]:
self.out.branch(prob["branchname"]+"_mavjj", "F", lenVar=lenVar_, title="User-defined mavjj probability")
if prob["addPmavjj_true"]:
self.out.branch(prob["branchname"]+"_mvajj_true", "F", lenVar=lenVar_, title="User-defined mvajj_true probability")
def fillProbs(self, candDaughters, candAssociated, candMothers):
if len(self.sortedSettings) == 0: return True
else:
Expand Down Expand Up @@ -141,8 +139,7 @@ def fillProbs(self, candDaughters, candAssociated, candMothers):
### Configure MELA for the event.
self.MELA.setProcess(MELA_Process, MELA_MatrixElement, MELA_Production)
self.MELA.differentiate_HWW_HZZ = MELA_separatewwzz
self.MELA.setMelaLeptonInterference(MELA_leptoninterference)

self.MELA.setMelaLeptonInterference(MELA_leptoninterference)

# Define arrays to fill with the probabilites for each candidate
probVec = [-999.]*len(candDaughters)
Expand Down Expand Up @@ -208,47 +205,49 @@ def fillProbs(self, candDaughters, candAssociated, candMothers):

# Handle divideP
if self.ModuleContext == "LHE":
vprob[iprob] = probVec[0]#probVec[iCand]

vprob[iprob] = probVec[0] # Note: iCand is always 0 for LHE
if MELA_divideP_idx != -1:
den = vprob[MELA_divideP_idx]
if den != 0:
probVec[iCand] /= den
probVec[0] /= den
else:
print("***MELAProbHelper: Protecting against division by 0.")
# else:
# print("***MELAProbHelper: Cannot run divideP for reco-level probability!")

# Write out all branches for this probability
if self.ModuleContext == "LHE":
self.out.fillBranch(MELA_branchname, probVec[0])
if MELA_ispm4l:
self.out.fillBranch(MELA_branchname+"_ScaleUp", probVec_ScaleUp[0])
self.out.fillBranch(MELA_branchname+"_ScaleDown", probVec_ScaleDown[0])
self.out.fillBranch(MELA_branchname+"_SystUp", probVec_SystUp[0])
self.out.fillBranch(MELA_branchname+"_SystDown", probVec_SystDown[0])
if MELA_computeprop:
self.out.fillBranch(MELA_branchname+"_prop", probPropVec[0])
else:
self.out.fillBranch(MELA_branchname, probVec)

if self.ModuleContext == "LHE": # Branches are not collections in this case
probVec=probVec[0]
if MELA_ispm4l:
self.out.fillBranch(MELA_branchname+"_ScaleUp", probVec_ScaleUp)
self.out.fillBranch(MELA_branchname+"_ScaleDown", probVec_ScaleDown)
self.out.fillBranch(MELA_branchname+"_SystUp", probVec_SystUp)
self.out.fillBranch(MELA_branchname+"_SystDown", probVec_SystDown)

if MELA_computeprop:
self.out.fillBranch(MELA_branchname+"_prop", probPropVec)

probVec_ScaleUp = probVec_ScaleUp[0]
probVec_ScaleDown = probVec_ScaleDown[0]
probVec_SystUp = probVec_SystUp[0]
probVec_SystDown = probVec_SystDown[0]
if MELA_computeprop and (MELA_prod or MELA_dec):
probPropVec = probPropVec[0]
if MELA_addPAux:
self.out.fillBranch(MELA_branchname+"_aux", probVec_PAux)

probVec_PAux = probVec_PAux[0]
if MELA_addPmavjj:
self.out.fillBranch(MELA_branchname+"_mavjj", probVec_mavjj)
probVec_mavjj = probVec_mavjj[0]
if MELA_addPmavjj_true:
self.out.fillBranch(MELA_branchname+"_mvajj_true", probVec_mvajj_true)

probVec_mvajj_true = probVec_mvajj_true[0]

self.out.fillBranch(MELA_branchname, probVec)

if MELA_ispm4l:
self.out.fillBranch(MELA_branchname+"_ScaleUp", probVec_ScaleUp)
self.out.fillBranch(MELA_branchname+"_ScaleDown", probVec_ScaleDown)
self.out.fillBranch(MELA_branchname+"_SystUp", probVec_SystUp)
self.out.fillBranch(MELA_branchname+"_SystDown", probVec_SystDown)

if MELA_computeprop and (MELA_prod or MELA_dec):
self.out.fillBranch(MELA_branchname+"_prop", probPropVec)

if MELA_addPAux:
self.out.fillBranch(MELA_branchname+"_aux", probVec_PAux)
if MELA_addPmavjj:
self.out.fillBranch(MELA_branchname+"_mavjj", probVec_mavjj)
if MELA_addPmavjj_true:
self.out.fillBranch(MELA_branchname+"_mvajj_true", probVec_mvajj_true)

return True

Expand Down
6 changes: 2 additions & 4 deletions NanoAnalysis/python/RecoProbFiller.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,17 +20,15 @@ def __init__(self, MELA, NANOVERSION, MELASettings = None):
self.sortedSettings = []
self.MELASettings = MELASettings
self.ProbHelper = MELAProbHelper(self.MELA, self.MELASettings, "Reco")
print("***RecoProbFiller: set for: ", self.ProbHelper.names, flush=True)


def beginFile(self, inputFile, outputFile, inputTree, wrappedOutputTree):
self.out = wrappedOutputTree
if self.MELASettings != None:
# print("***RecoProbFiller: Booking Probs")
self.ProbHelper.bookProbs(wrappedOutputTree)





def analyze(self, event):
if event.nZZCand==0: return True
cands = Collection(event, 'ZZCand')
Expand Down
Loading