Randomly selecting 3 events from TTree, which has 1 million events

ROOT
1

ROOT Version: 6.08.06:
Compiler: 5.2.1

Hi, I wonder if I can get help from this website about PyRoot

I currently have a root file called qcd.root (I wish I could have uploaded in this forum, but the size of file is 3.1 GB)

It is a root file, which as been generated through following steps

Use PythiaRun to generate 100000 pp->jj+X QCD events, output file being qcd.hepmc

Use DelphesHepmc executable with CMS card file does jet clustering and generates qcd.root

Therefore,qcd.root is a 100000 QCD Jet clustered events file, which has Delphes TTree, and within there, branches of Event, Weight, Track, Jet, and so on.

My primary goal is to randomly sample N events from the sample of N events and histogram all jets with |y| < 0.5 and pT > 1 TeV.

I currently have my code, which works, but lack completeness.

`from ROOT import gSystem, TFile, TH1F, TCanvas
from histutil import *
#-------------------------------------------------------------------------------
rootfile = 'Jet.root'
def main():
   
    # read all entries and fill the histograms
    myfile = TFile("qcd.root")
    mytree = myfile.Delphes
    hfile = TFile(rootfile, "recreate")
    
    if not hfile.IsOpen():
        hutil.error("randomhist.py",
                    "can't open rootfile Event.root")
    gSystem.Load('libDelphes')   # for i in range(1000):
    hfile.cd()
    c1 = TCanvas('c1', 'Histogram', 200, 10, 700, 900)
   # s = str(i)
    hpx = TH1F("hpx","px distribution 1st", 100, 100, 4200)
    for event in mytree:
        for jet in event.Jet:
            hpx.Fill(jet.PT)
        hpx.Write()
    myfile.Close()
    hfile.Close()
    
#-------------------------------------------------------------------------------
try:
    main()
except KeyboardInterrupt:
    print 'ciao!'
`

For this code, I have 1 main question.

In my current code, the for loop for event goes over ENTIRE 1 MILLION events and creates that number of hpx. However, I want to Randomly select only 100 event from my tree and and then write out 1 px.

How may I do so?

Thank you for help!
In addition, if you notice something wrong or inefficient about my code, please let me know!

2

Hi @trenta_coollime,

The way to do it is to generate some random numbers, e.g. with ROOT.gRandom (see here) to select the 100 events you want to fill your histogram with.
The way you tell the TTree that you want only those 100 random events is by looping a bit differently, so instead of doing:

for event in mytree:  # this loops over everything

you do

event = mytree.event
for i in get_random_entry_numbers():`
  mytree.GetEntry(i)
  for jet in event.Jet:
  ...

The function get_random_entry_numbers() could be a generator of random numbers with ROOT.gRandom.

The main point here is to is call GetEntry() on the tree to specifically get the entries you want.

Cheers,
Enric

1 Like
3

Thank you for the help.

However, the code written above did not work for me.
It says that TTree does not have objectification for event

I looked at my qcd.root to find source of problem and its structure is like

qcd.root ->
Delphes Tree ->
Branches of “Event”, “Weight”, “Particle”, “Jet”, and so on

Would there be a way to use the method for my Root file’s structure?

Thank you!
(In addition, how may I specify the random events number, not necessarily to be 100, but also 3 or 5 or any other arbitrary number)

4

HI @trenta_coollime,

Can you try with the following code:

for i in get_random_entry_numbers():`
  mytree.GetEntry(i)
  jets = getattr(myTree, 'event/Jet')
  for jet in jets:
  ...

If this does not work, can you paste here the result of myTree.Print(), please?

Cheers,
Enric

5

someFile.txt (116.0 KB)
I attached the .txt output for mytree.Print() output.

For your second suggestion, it claims

Traceback (most recent call last):
  File "randomhist.py", line 36, in <module>
    main()
  File "randomhist.py", line 24, in main
    for i in get_random_entry_numbers():
NameError: global name 'get_random_entry_numbers' is not defined

My code at the moment looks like

from ROOT import gSystem, TFile, TH1F, TCanvas, gRandom
from histutil import *
from array import array
#-------------------------------------------------------------------------------
rootfile = 'Jet.root'
xbins = array('f',[56.0,64.0,74.0,84.0,97.0,114.,133.,153.,174.,196.,220.,245.,272.,300.,330.,362.,395.,430.,468.,507.,548.,592.,638.,686.,737.,790.,846.,905.,967.,1032.,1101.,1172.,1248.,1327.,1410.,1497.,1588.,1684.,1784.,1890.,2000.,2116.,2238.,2366.,2500.,2640.,2787.,2941.,3103.,3273.,3450.,3637.,3832.,4037.])

def main():
    gSystem.Load('libDelphes')   # for i in range(1000):
    # read all entries and fill the histograms
    myfile = TFile("qcd.root")
    mytree = myfile.Delphes
    hfile = TFile(rootfile, "recreate")
    
    if not hfile.IsOpen():
        hutil.error("randomhist.py",
                    "can't open rootfile Event.root")
    
    hfile.cd()
   # c1 = TCanvas('c1', 'Histogram', 200, 10, 700, 900)
   # s = str(i)
    hpx = TH1F("hpx","px distribution 1st",53,xbins)
    
    for i in get_random_entry_numbers():
        mytree.GetEntry(i)
        jets = getattr(myTree, 'event/Jet')
        for jet in jets:
            if jet.PT > 1000 and abs(jet.Eta)<.5:
                hpx.Fill(jet.PT)
        hpx.Draw()
    myfile.Close()
    hfile.Close()
    
#-------------------------------------------------------------------------------
try:
    main()
except KeyboardInterrupt:
    print 'ciao!'

I appreciate your help
Thank you.

6

Hi @trenta_coollime ,

I see in the result of Print that you have a branch called Jet that is an array of jets, I guess that is the one you want to use.

The function get_random_entry_numbers is something you need to implement. As I suggested above, it could be a generator of random numbers with ROOT.gRandom. For the sake of making it work, you can start with something like:

for i in [0,1,2]:
  mytree.GetEntry(i)
  jets = mytree.Jet
  for jet in jets:
  ...

or alternatively, this:

for i in [0,1,2]:
  mytree.GetEntry(i)
  jets = getattr(mytree, 'Jet')
  for jet in jets:
  ...

Cheers,
Enric

1 Like
7

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.