Read lhe and make histos

an.py

import sys
#sys.path.insert(0, '/gwpool/users/tecedor/.local/lib/python3.9/site-packages/')
sys.path.insert(0, '/gwpool/users/gpizzati/.local/lib/python3.9/site-packages/')
import ROOT
import pylhe
import awkward as ak
import numpy as np
#arr = pylhe.to_awkward(pylhe.read_lhe_with_attributes('/gwpool/users/tecedor/genproductions/GridpackConfig/bin/MadGraph5_aMCatNLO/lhe/HHjj_ewk_smhloop0_dim6_3ops_15904236_10.lhe'))

files = []
files = list(map(lambda k: '/gwpool/users/tecedor/genproductions/GridpackConfig/bin/MadGraph5_aMCatNLO/lhe/' + k, files))

index = int(sys.argv[0])
file = files[index]

print(file, index)

arr = pylhe.to_awkward(pylhe.read_lhe_with_attributes(file))

#arr = pylhe.to_awkward(pylhe.read_lhe_with_attributes('HHjj_ewk_smhloop0_dim6_3ops_15904236_10.lhe'))

#norm = 1 / 200
#ak.sum(200 * arr.weights.values[:, 0] * norm)

#-------------------Higgs variables-----------------------------

mHH=(arr.particles.vector[:, -4]+arr.particles.vector[:, -3]).M
deltaPhiHH = (arr.particles.vector[:, -4]).deltaphi(arr.particles.vector[:, -3])
deltaEtaHH = (arr.particles.vector[:, -4]).deltaeta(arr.particles.vector[:, -3])

firstIsFirst_H = np.array(arr.particles.vector[:, -4].pt > arr.particles.vector[:, -3].pt)

ptH1  = np.array(arr.particles.vector[:, -4].pt)
phiH1 = np.array(arr.particles.vector[:, -4].phi)
etaH1 = np.array(arr.particles.vector[:, -4].eta)
ptH1 [~firstIsFirst_H] = np.array(arr.particles.vector[:, -3].pt) [~firstIsFirst_H]
phiH1[~firstIsFirst_H] = np.array(arr.particles.vector[:, -3].phi)[~firstIsFirst_H]
etaH1[~firstIsFirst_H] = np.array(arr.particles.vector[:, -3].eta)[~firstIsFirst_H]

ptH2  = np.array(arr.particles.vector[:, -3].pt)
phiH2 = np.array(arr.particles.vector[:, -3].phi)
etaH2 = np.array(arr.particles.vector[:, -3].eta)
ptH2 [~firstIsFirst_H] = np.array(arr.particles.vector[:, -4].pt) [~firstIsFirst_H]
phiH2[~firstIsFirst_H] = np.array(arr.particles.vector[:, -4].phi)[~firstIsFirst_H]
etaH2[~firstIsFirst_H] = np.array(arr.particles.vector[:, -4].eta)[~firstIsFirst_H]

#-------------------Jet variables-------------------------------

mJJ=(arr.particles.vector[:, -2]+arr.particles.vector[:, -1]).M

firstIsFirst_J = np.array(arr.particles.vector[:, -2].pt > arr.particles.vector[:, -1].pt)

deltaPhiJJ = (arr.particles.vector[:, -2]).deltaphi(arr.particles.vector[:, -1])
deltaEtaJJ = (arr.particles.vector[:, -2]).deltaeta(arr.particles.vector[:, -1])

ptJ1 =np.array( arr.particles.vector[:, -2].pt)
phiJ1  =np.array( arr.particles.vector[:, -2].phi)
etaJ1 =np.array( arr.particles.vector[:, -2].eta)
ptJ1 [~firstIsFirst_J] =np.array( arr.particles.vector[:, -1].pt) [~firstIsFirst_J]
phiJ1 [~firstIsFirst_J] = np.array(arr.particles.vector[:, -1].phi) [~firstIsFirst_J]
etaJ1 [~firstIsFirst_J] = np.array(arr.particles.vector[:, -1].eta) [~firstIsFirst_J]

ptJ2 =np.array( arr.particles.vector[:, -1].pt)
phiJ2  =np.array( arr.particles.vector[:, -1].phi)
etaJ2 = np.array(arr.particles.vector[:, -1].eta)
ptJ2 [~firstIsFirst_J] =np.array( arr.particles.vector[:, -2].pt) [~firstIsFirst_J]
phiJ2 [~firstIsFirst_J] = np.array(arr.particles.vector[:, -2].phi) [~firstIsFirst_J]
etaJ2 [~firstIsFirst_J] = np.array(arr.particles.vector[:, -2].eta) [~firstIsFirst_J]

makeDatacard.py

import ROOT



rootFile = '.root'

f = ROOT.TFile(rootFile)
operator = 'cH'
binName = 'mHH_ ' + operator

header = f'''
## Shape input card
imax 1 number of channels
jmax * number of background
kmax * number of nuisance parameters
----------------------------------------------------------------------------------------------------
bin         {binName}
observation 0
shapes  *           * output_hist.root     histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
shapes  data_obs           * output_hist.root     histo_Data
bin     {binName} {binName} {binName}
'''




rows = [['process'], ['process'], ['rate']]

shapes = [f'quad_{operator}', f'sm_lin_quad_{operator}', 'sm']


for i, shape in enumerate(shapes):
    h = f.Get(f'histo_{shape}')
    tot = h.Integral()
    rows[0].append(shape)
    rows[1].append(str(i+1))
    rows[2].append(str(round(tot, 5)))

txt = header 
for row in rows:
    txt += '\t'.join(row) + '\n'

with open('datacard.txt', 'w') as file:
    file.write(txt)

reader.py

import sys
sys.path.insert(0, '/gwpool/users/tecedor/.local/lib/python3.9/site-packages/')
import ROOT
ROOT.EnableImplicitMT()
ROOT.gROOT.SetBatch(True)
ROOT.TH1.SetDefaultSumw2(True)
import argparse 
files = ['output_' + str(i) + '.root' for i in range(0,200)]

path = '/gwpool/users/tecedor/condor_analyzer/root/'
files = list(map(lambda k: path + k, files))
n_files = len(files)


df = ROOT.RDataFrame("Events", files)

weights = {
        "w_sm": f"weights[0]/{n_files}",
        "w_lin_cH": f"(0.5*(weights[2]/{n_files} - weights[1]/{n_files}))",
        "w_quad_cH": f"(-w_sm + 0.5*(weights[2]/{n_files} + weights[1]/{n_files}))",
        "w_sm_lin_quad_cH": f"(weights[2]/{n_files})",
        "w_lin_cHW": f"(0.5*(weights[4]/{n_files} - weights[3]/{n_files}))",
        "w_quad_cHW": f"(-w_sm + 0.5*(weights[4]/{n_files} + weights[3]/{n_files}))",
        "w_sm_lin_quad_cHW": f"(weights[4]/{n_files})",
        "w_lin_cHq3": f"(0.5*(weights[6]/{n_files} - weights[5]/{n_files}))",
        "w_quad_cHq3": f"(-w_sm + 0.5*(weights[6]/{n_files} + weights[5]/{n_files}))",
        "w_sm_lin_quad_cHq3": f"(weights[6]/{n_files})",
        }

def defaultParser():

    parser = argparse.ArgumentParser()

    parser.add_argument("-v",
                        type = str,
                        )

    parser.add_argument("-l",
                        type = int,
                        )

    parser.add_argument("-xmin",
                        type = int,
                        )

    parser.add_argument("-xmax",
                        type = int,
                        )
    return parser

parser = defaultParser()
args = parser.parse_args()
print(args.v, args.l, args.xmin, args.xmax)

variable = args.v
lum = args.l
xmax = args.xmax
xmin = args.xmin


for key in weights.keys():
    df = df.Define(key, weights[key])


histos = []
for w in weights.keys():
    weight = f'{lum} * 1000.0 * ' + w
    df = df.Redefine(w, weight)
    histos.append(df.Histo1D(('histo_' + '_'.join(w.split('_')[1:]), "", 40, int(args.xmin),int(args.xmax)), variable, w))

f = ROOT.TFile('output_hist.root', 'RECREATE')
f.cd()

savedData = False
for histo in histos:
    histo.Write()
    if not savedData:
        h = histo.Clone()
        for i in range(0, h.GetNbinsX() +1):
            h.SetBinContent(i, 0)
        h.SetName('histo_Data')
        h.Write()
        savedData = True


f.Close()

"""
#h=df.Histo1D("mHH","w")
c = ROOT.TCanvas("c","c",1000,1000)
colors = [ROOT.kBlack, ROOT.kBlue, ROOT.kRed]
for h in histos:
    h.SetLineColor(colors[histos.index(h)])
    h.Draw('same')
    print(h.Integral())
c.SetLogy()
c.Draw()
c.Print("test.png")
"""


"""

for f in files:
    file = ROOT.TFile(f)
    f.Get(


df = ROOT.RDataFrame("Events", files)
"""