diff --git a/NanoAnalysis/python/LHEAngProbFiller.py b/NanoAnalysis/python/LHEAngProbFiller.py index 8d5c3a7aff8..c3dad8c4b0a 100644 --- a/NanoAnalysis/python/LHEAngProbFiller.py +++ b/NanoAnalysis/python/LHEAngProbFiller.py @@ -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 = { @@ -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.") @@ -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]) diff --git a/NanoAnalysis/python/MELAProbHelper.py b/NanoAnalysis/python/MELAProbHelper.py index 1d48d7735f4..9a749b8b3b5 100644 --- a/NanoAnalysis/python/MELAProbHelper.py +++ b/NanoAnalysis/python/MELAProbHelper.py @@ -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 = { @@ -27,7 +31,6 @@ def __init__(self, MELA, MELASettings, ModuleContext): "Prod": None, "Dec": None, "Couplings": None, - # "isgen": None, "context": None, "computeprop": None, "propscheme": "FixedWidth", @@ -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) @@ -69,14 +75,15 @@ 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. @@ -84,35 +91,26 @@ def bookProbs(self, 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: @@ -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) @@ -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 diff --git a/NanoAnalysis/python/RecoProbFiller.py b/NanoAnalysis/python/RecoProbFiller.py index 3db8aebae58..6a6e5a54ecd 100644 --- a/NanoAnalysis/python/RecoProbFiller.py +++ b/NanoAnalysis/python/RecoProbFiller.py @@ -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') diff --git a/NanoAnalysis/python/nanoZZ4lAnalysis.py b/NanoAnalysis/python/nanoZZ4lAnalysis.py index 6b8d9d4194f..014f4ad2f12 100644 --- a/NanoAnalysis/python/nanoZZ4lAnalysis.py +++ b/NanoAnalysis/python/nanoZZ4lAnalysis.py @@ -254,15 +254,20 @@ post_sequence.append(mcTruthAnalyzer()) # Gen final state etc. # from ZZAnalysis.NanoAnalysis.genExtraFiller import * # post_sequence.append(genExtraFiller(mela)) Gen-level angles (not to be confused with LHE-level angles, filled by LHEAngProbFiller - + + if runMELA: + from ZZAnalysis.NanoAnalysis.LHEAngProbFiller import * + if NANOVERSION >= 15: + from ZZAnalysis.NanoAnalysis.LHEFiller import * + post_sequence.append(LHEFiller()) + post_sequence.append(LHEAngProbFiller(mela, NANOVERSION, melaSettings)) + if ADD_ALLEVENTS: # Add modules that produce the variables to be stored for all events at the beginning from ZZAnalysis.NanoAnalysis.genFiller import * from ZZAnalysis.NanoAnalysis.cloneBranches import * - from ZZAnalysis.NanoAnalysis.LHEAngProbFiller import * pre_sequence = [puWeight(LEPTON_SETUP, DATA_TAG), weights, genFiller(mela, dump=False), - LHEAngProbFiller(mela, NANOVERSION, melaSettings), cloneBranches(treeName='AllEvents', varlist=['run', 'luminosityBlock', 'event', 'GenDressedLepton_*', @@ -274,7 +279,6 @@ 'puWeight*', 'overallEventWeight', 'Pileup_nTrueInt', - 'LHEMela*', 'GenJet*', *(['ggH_NNLOPS_Weight'] if APPLY_QCD_GGF_UNCERT else []), ], @@ -282,11 +286,6 @@ continueFor = postPresel ), ] + pre_sequence - if NANOVERSION >= 15: - from ZZAnalysis.NanoAnalysis.LHEFiller import * - insertBefore(pre_sequence, 'LHEAngProbFiller', LHEFiller()) - - else : # Add them at the end, so that they are run only for selected events post_sequence.extend([puWeight(LEPTON_SETUP,DATA_TAG), @@ -337,6 +336,7 @@ #'keep PV*', #'keep Flag*', *(['keep MET_pt'] if NANOVERSION <=12 else ['keep PFMET_pt']), + *(['keep LHEMela*'] if runMELA else []), ] if IsMC: diff --git a/NanoAnalysis/test/prod/pyFragments/exampleProbabilities.py b/NanoAnalysis/test/prod/pyFragments/exampleProbabilities.py index 590c88218bf..154f69fc6ec 100644 --- a/NanoAnalysis/test/prod/pyFragments/exampleProbabilities.py +++ b/NanoAnalysis/test/prod/pyFragments/exampleProbabilities.py @@ -30,11 +30,14 @@ # - "useconstant": Boolean. This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true. -# - "match_mX": Boolean. If true, will set the Higgs mass to match the invariant mass of the daughter particles in each event. +# - "match_mX": Boolean. If true, will set the Higgs mass to match the invariant mass of the daughter particles in each event. # - "lepton_interference":"DefaultLeptonInterf", # - "ispm4l": Boolean. When True, causes a call to computePM4l(), which is a probability useful in calculating a signal-background discrimimant. -# - "dividep": NOTE: WIP. When finished, will take the "Name" variable of one probability and divide all probabilities with divdep=True by this probability. This is typically most useful for normalizing to a native probability. - +# - "dividep": Take the "Name" variable of one probability and divide all probabilities with divdep=True by this probability. This is typically most useful for normalizing to a native probability. Supported only for context="LHE". +# - "addPAux": also store the result of MELA.getPAux() - cf. slide 5 of https://indico.cern.ch/event/1619963/#11-hzz4l-legacy-signal-strengh +# - "addPmavjj": also store the result of computeDijetConvBW(False) - cf. slide 8 of https://indico.cern.ch/event/1619963/contributions/6826626/attachments/3190935/5678923/STXSupdate.pdf +# - "addPmavjj_true": also store the result of computeDijetConvBW(False) + #As a couple of general notes: # - append=True should always be set for each instance of setConf. diff --git a/NanoAnalysis/test/runLocal.py b/NanoAnalysis/test/runLocal.py index 9f87a3ac924..76cc4a768fb 100755 --- a/NanoAnalysis/test/runLocal.py +++ b/NanoAnalysis/test/runLocal.py @@ -18,7 +18,7 @@ #SampleToRun = "MC2023postBPix" #SampleToRun = "Data2024" #SampleToRun = "MC2024" -# SampleToRun = "MELA_Test" +#SampleToRun = "MELA_Test" #SampleToRun = "forNanoDoc" # To prepare variable lists with inspectNanoFile.py ### Obsolete Run2 samples @@ -221,14 +221,10 @@ setConf("NANOVERSION", 15) setConf("store","root://cms-xrd-global.cern.ch/") # Add probabilities - # import prod.pyFragments.LHE_Probs_ggH0PM_M125 import prod.pyFragments.EFT_RecoProbs # import prod.pyFragments.contextTest - # import prod.pyFragments.DefaultProbs - # import prod.pyFragments.LHE_Probs_ggH0PM_M125 - # import prod.pyFragments.EFT_RecoProbs - # import prod.pyFragments.contextTest - # import prod.pyFragments.DefaultProbs + import prod.pyFragments.DefaultProbs + import prod.pyFragments.LHE_Probs_ggH0PM_M125 setConf("fileNames", [ "/store/mc/RunIII2024Summer24NanoAODv15/GluGluH-Hto2Zto4L_Par-M-125_TuneCP5_13p6TeV_powheg-jhugen-pythia8/NANOAODSIM/150X_mcRun3_2024_realistic_v2-v2/110000/a9e03ff9-2146-4aff-bd26-69abcb98359f.root" # 10000 evts ])