#!/usr/bin/env python import sys from array import array if len(sys.argv) < 4: print(" Usage: DELPHES.py input_file type output_file") sys.exit(1) #import os # add Delphes lib path to system LD_LIBRARY_PATH evironment variable #pathDelphesLib = "" #os.environ['LD_LIBRARY_PATH'] = pathDelphesLib+':'+os.environ['LD_LIBRARY_PATH'] import ROOT # 必须在Delphes的根目录下运行 ROOT.gSystem.Load("libDelphes.so") #pathDelphes = "" #ROOT.gInterpreter.AddIncludePath(pathDelphes) # declare classes #ROOT.gInterpreter.Declare('#include "classes/DelphesClasses.h"') #ROOT.gInterpreter.Declare('#include "external/ExRootAnalysis/ExRootTreeReader.h"') print("input") inputFile = sys.argv[1] #sys.exit(0) # Create chain of root trees chain = ROOT.TChain("Delphes") chain.Add(inputFile) sampleType = sys.argv[2] # Create object of class ExRootTreeReader treeReader = ROOT.ExRootTreeReader(chain) numberOfEntries = treeReader.GetEntries() # Get pointers to branches used in this analysis #branchParticle = treeReader.UseBranch('Particle') branchJet = treeReader.UseBranch("Jet") branchElectron = treeReader.UseBranch("Electron") branchWeight = treeReader.UseBranch("Weight") branchEvent = treeReader.UseBranch("Event") output_filename = sys.argv[3] outputFile = ROOT.TFile(output_filename,"UPDATE") # Book histograms histInvMass = ROOT.TH1F("mass_eplus_eminus", "Invariant mass of e+ e-", 50, 0.0, 150.0) histInvMass.Sumw2(True) # Book TTree treeSlimmed = ROOT.TTree("tree","slimmed tree") mass = array('d', [0]) treeSlimmed.Branch('massReco', mass, 'massReco/D') xsec = array('d',[0]) treeSlimmed.Branch('xsec', xsec, 'xsec/D') N = array('d',[0]) treeSlimmed.Branch('N', N, 'N/I') eventWeight = array('d',[0]) treeSlimmed.Branch('eventWeight', eventWeight, 'eventWeight/D') # Loop over all events for entry in range(0, numberOfEntries): # Load selected branches with data from specified event treeReader.ReadEntry(entry) ## main MC event weight weight = branchEvent[0].Weight electrons = [] positrons = [] #for electron in branchElectron for i in range(branchElectron.GetEntries()) : electron = branchElectron.At(i) charge = electron.Charge pt = electron.PT eta = electron.Eta phi = electron.Phi m = 0.5/1000 #0.5MeV ## looseCut = 0.3, mediumCut = 0.2, tightCut = 0.1 IsoCut = 0.3 passIso = electron.IsolationVar < IsoCut if(passIso and electron.Charge < 0 ) : electrons.append(electron) if(passIso and electron.Charge > 0 ) : positrons.append(electron) if(len(electrons)==1 and len(positrons)==1): massReco = (electrons[0].P4() + positrons[0].P4()).M() histInvMass.Fill(massReco, weight/numberOfEntries) mass[0] = massReco eventWeight[0] = weight/numberOfEntries xsec[0] = weight N[0] = numberOfEntries treeSlimmed.Fill() outputFile.mkdir(sampleType) outputFile.cd(sampleType) histInvMass.Write() treeSlimmed.Write() outputFile.Close()
Baidu
map