How to write a code to read entry in a Delphes output file event by event?

Hi,

I want to write a code using python to : In a Delphes output root file , using a TChain to add the root file , and read it event by event, and for each event , read the Jet , electron , muon …etc. ( In Delphes , it suggest us to use DelphesAnalysis but I find it is hard to use, so I want to write something I can control)

for example, like the example in Delphes( but use ExRootAnalysis and C++) ,

    TClonesArray *branchElectron = treeReader->UseBranch("Electron");
    
    for(entry = 0; entry < allEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);

    // Loop over all electrons in event
    for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
    {
      electron = (Electron*) branchElectron->At(i);
      particle = (GenParticle*) electron->Particle.GetObject();
    }
  }

I tried it , for example , get the root file by:

chain = TChain('Delphes') chain.Add('zhh.root')
and get the Branch by:

but I don’t know how to read the Branch event by event.

Also I have read the example in root tutorials,

[code] t1->SetBranchAddress(“px”,&px);
t1->SetBranchAddress(“py”,&py);
t1->SetBranchAddress(“pz”,&pz);
t1->SetBranchAddress(“random”,&random);
t1->SetBranchAddress(“ev”,&ev);

//create two histograms
TH1F *hpx = new TH1F(“hpx”,“px distribution”,100,-3,3);
TH2F *hpxpy = new TH2F(“hpxpy”,“py vs px”,30,-3,3,30,-3,3);

//read all entries and fill the histograms
Long64_t nentries = t1->GetEntries();
for (Long64_t i=0;i<nentries;i++) {
t1->GetEntry(i);
hpx->Fill(px);
hpxpy->Fill(px,py);
}[/code]

It seems first in Python there is no reference or address so I don’t know how to use “SetBranchAddress” and second it seems not event by event.

Is there anyone can help me ?

Best,
Li

Hi Li,

you can proceed as follows:

[....]
#create two histograms
hpx   = ROOT.TH1F("hpx","px distribution",100,-3,3);
hpxpy = ROOT.TH2F("hpxpy","py vs px",30,-3,3,30,-3,3);
   
# read all entries and fill the histograms
myfile = ROOT.TFile("zhh.root")
mytree = myfile.Delphes
for event in mytree:
    jet = event.Jet
    hpx.Fill(jet.px)
    hpxpy.Fill(jet.px, jet.py)

It is not clear if you have only one jet or many. In case you have many

for event in mytree:
    for jet in event.Jet:
        hpx.Fill(jet.px)
        hpxpy.Fill(jet.px, jet.py)

Cheers,
Danilo

[quote=“dpiparo”]Hi Li,

you can proceed as follows:

[....]
#create two histograms
hpx   = ROOT.TH1F("hpx","px distribution",100,-3,3);
hpxpy = ROOT.TH2F("hpxpy","py vs px",30,-3,3,30,-3,3);
   
# read all entries and fill the histograms
myfile = ROOT.TFile("zhh.root")
mytree = myfile.Delphes
for event in mytree:
    jet = event.Jet
    hpx.Fill(jet.px)
    hpxpy.Fill(jet.px, jet.py)

It is not clear if you have only one jet or many. In case you have many

for event in mytree:
    for jet in event.Jet:
        hpx.Fill(jet.px)
        hpxpy.Fill(jet.px, jet.py)

Cheers,
Danilo[/quote]

Hi, Danilo
It is amazing for the code looks so nice. But I have some more questions:

  1. sometimes I have more files, for example , zhh0.root zhh1,root zhh2.root …etc .I usually deal with multi-files by TChain , how to deal with this issue with TFile ? Or how to write a TChain version ?

  2. I add :

for event in mytree: print event.Jet.GetEntries() for jet in event.Jet: hpx.Fill(jet.PT) print jet.Track.GetEntries()
I think the second line , I add, can get the number of jet information, am I right?

  1. Physically a jet contains tracks,etc… so

[code]
for(j = 0; j < jet->Constituents.GetEntriesFast(); ++j)
{
object = jet->Constituents.At(j);

    // Check if the constituent is accessible
    if(object == 0) continue;

    if(object->IsA() == GenParticle::Class())
    {
      particle = (GenParticle*) object;
      // cout << "    GenPart pt: " << particle->PT << ", eta: " << particle->Eta << ", phi: " << particle->Phi << endl;
      momentum += particle->P4();
    }
    else if(object->IsA() == Track::Class())
    {
      track = (Track*) object;
      // cout << "    Track pt: " << track->PT << ", eta: " << track->Eta << ", phi: " << track->Phi << endl;
      momentum += track->P4();
    }
    else if(object->IsA() == Tower::Class())
    {
      tower = (Tower*) object;
      // cout << "    Tower pt: " << tower->ET << ", eta: " << tower->Eta << ", phi: " << tower->Phi << endl;
      momentum += tower->P4();
    }
  }
  [/code]

The IsA() method provide a way to identify the institutions of Jet, how can I deal with it in the python ?

Thanks again for the brief code. And this is why I want to learn python coding, they are so readable.

Best,
Li

Hi Li,

glad to hear that. Thank Python :slight_smile:

Now, your points:

  1. Of course you can use the TChain. My bad for using a TFile only.
  2. It depends on the type of the object in the branch. If it is a TClonesArray or if it implements the GetEntries method, yes.
  3. The same as in C++? Did you try and encounter problems?

Cheers,
Danilo

All work: IsA(), isinstance(), type(). But! No need your case: all branch call P4(), so all object has P4(), yes? Alors, always good.

-Dom

Good point Dom!

Danilo

[quote=“dpiparo”]Hi Li,

glad to hear that. Thank Python :slight_smile:

Now, your points:

  1. Of course you can use the TChain. My bad for using a TFile only.
  2. It depends on the type of the object in the branch. If it is a TClonesArray or if it implements the GetEntries method, yes.
  3. The same as in C++? Did you try and encounter problems?

Cheers,
Danilo[/quote]

Hi Danilo,
Yes , my problem is that,

import ROOT
ROOT.gSystem.Load('libDelphes')

chain = ROOT.TChain('Delphes')
chain.Add('zhh.root')

for event in chain:
    #print event.Jet.GetEntries()
    for jet in event.Jet:
        #hpx.Fill(jet.PT)
        for constituent in jet.Constituents:
            if( constituent.IsA() == Track ):
            #if(isinstance(constituent,Track)):
                print ' it is a track '

There is a bug:

ui@ui:~/tools/package/Delphes-3.3.2$ python test.py TClass::TClass:0: RuntimeWarning: no dictionary for class CompBase is available Traceback (most recent call last): File "test.py", line 42, in <module> if( constituent.IsA() == Track ): NameError: name 'Track' is not defined

The problem is that I do not know how to point the class Track which in C++ , it is :

the Track::Class() give something to compare . How can I do such thing in python?

Best,
Li

All work: IsA(), isinstance(), type(). But! No need your case: all branch call P4(), so all object has P4(), yes? Alors, always good.

-Dom[/quote]

Hi Dom,
Thanks for your advise , the point I want to know the constituents in jet is Track , Tower or other things is to reconstruct vertex using Track, for example.

I tried this:

It seems that it can find the right Track object , but it is a little ugly and low efficiency and may be wrong ( if there is no Track in the event ? ) , so I want to write a right way like:

But I don’t know how to write “the Track class” in the python.

Best,
Li

[code]type(constituent) == ROOT.Track[/code]
[code]constituent.IsA() == ROOT.Track.Class()[/code]
first expect faster (2 & 3 need extra call comparison).

-Dom

first expect faster (2 & 3 need extra call comparison).

-Dom

[quote=“Dominique”]isinstance(constituent, ROOT.Track)

first expect faster (2 & 3 need extra call comparison).

-Dom[/quote]

Hi Dom,
Thanks a lot!
I am so surprise that I forget that…
It is surprise that I thought ExRootAnalysis ( used in Delphes) should have re-wrote the Track class ( maybe a subclass ), so the Track class in the Delphes Class should not be equal to the ROOT.Track . Now it seems ExRootAnalysis do not re-write it.

Again, thanks a lot! Now I can write a python version to analysis the Delphes output file.

Best,
Li

Hi Li,

is the missing dictionary issue also fixed?

Danilo

[quote=“dpiparo”]Hi Li,

is the missing dictionary issue also fixed?

Danilo[/quote]

Hi Danilo,
no. I didn’t. Using

import ROOT
ROOT.gSystem.Load('libDelphes')
chain = ROOT.TChain('Delphes')
chain.Add('zhh.root')

for event in chain:
    #print event.Jet.GetEntries()
    for jet in event.Jet:
        pass

with " for event in chain " there is a warning , and using

import ROOT
ROOT.gSystem.Load('libDelphes')
myfile = ROOT.TFile("zhh.root")
mytree = myfile.Delphes

for event in mytree:
    #print event.Jet.GetEntries()
    for jet in event.Jet:
        pass

there is a same warning, too.

Best,
Li

Hi Li,

that is a consistent behaviour.
What one should find out is where the dictionary for “CompBase” is located (it is not a ROOT class*). Do you have a hint?

Cheers,
Danilo

cp3.irmp.ucl.ac.be/projects/del … sClasses.h

[quote=“dpiparo”]Hi Li,

that is a consistent behaviour.
What one should find out is where the dictionary for “CompBase” is located (it is not a ROOT class*). Do you have a hint?

Cheers,
Danilo

cp3.irmp.ucl.ac.be/projects/del … sClasses.h[/quote]

Hi,
I guess it’s from Delphes ,

ui@ui:~/tools/package/Delphes$ grep -nr "CompBase" ./
Binary file ./pileup2root matches
Binary file ./DelphesSTDHEP matches
Binary file ./DelphesHepMC matches
Binary file ./root2pileup matches
Binary file ./libDelphes.so matches
Binary file ./lhco2root matches
./tmp/classes/ClassesDict.h:65:extern G__linked_taginfo G__ClassesDictLN_CompBase;
Binary file ./tmp/classes/DelphesClasses.o matches
./tmp/classes/ClassesDict.cc:6559:   G__memvar_setup((void*)(&Candidate::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6677:   G__memvar_setup((void*)(&GenParticle::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6762:   G__memvar_setup((void*)(&Photon::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6786:   G__memvar_setup((void*)(&Electron::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6809:   G__memvar_setup((void*)(&Muon::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6852:   G__memvar_setup((void*)(&Jet::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6885:   G__memvar_setup((void*)(&Track::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6906:   G__memvar_setup((void*)(&Tower::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:6925:   G__memvar_setup((void*)(&HectorHit::fgCompare),85,0,0,G__get_linked_tagnum(&G__ClassesDictLN_CompBase),-1,-2,1,"fgCompare=",0,"!");
./tmp/classes/ClassesDict.cc:7013:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7047:   G__memfunc_setup("GetCompare",999,G__ClassesDict_522_0_3, 85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 3);
./tmp/classes/ClassesDict.cc:7200:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7359:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7381:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7409:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7437:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7459:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7487:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7509:   G__memfunc_setup("GetCompare",999,(G__InterfaceMethod) NULL,85, G__get_linked_tagnum(&G__ClassesDictLN_CompBase), -1, 0, 0, 1, 1, 9, "", (char*)NULL, (void*) NULL, 1);
./tmp/classes/ClassesDict.cc:7717:G__linked_taginfo G__ClassesDictLN_CompBase = { "CompBase" , 99 , -1 };
./tmp/classes/ClassesDict.cc:7789:  G__ClassesDictLN_CompBase.tagnum = -1 ;
./tmp/classes/ClassesDict.cc:7864:   G__get_linked_tagnum_fwd(&G__ClassesDictLN_CompBase);
Binary file ./tmp/classes/ClassesDict.o matches
Binary file ./root2lhco matches
Binary file ./stdhep2pileup matches
Binary file ./libDelphesNoFastJet.so matches
Binary file ./DelphesLHEF matches
Binary file ./Example1 matches
./classes/DelphesClasses.cc:35:CompBase *GenParticle::fgCompare = 0;
./classes/DelphesClasses.cc:36:CompBase *Photon::fgCompare = CompPT<Photon>::Instance();
./classes/DelphesClasses.cc:37:CompBase *Electron::fgCompare = CompPT<Electron>::Instance();
./classes/DelphesClasses.cc:38:CompBase *Muon::fgCompare = CompPT<Muon>::Instance();
./classes/DelphesClasses.cc:39:CompBase *Jet::fgCompare = CompPT<Jet>::Instance();
./classes/DelphesClasses.cc:40:CompBase *Track::fgCompare = CompPT<Track>::Instance();
./classes/DelphesClasses.cc:41:CompBase *Tower::fgCompare = CompE<Tower>::Instance();
./classes/DelphesClasses.cc:42:CompBase *HectorHit::fgCompare = CompE<HectorHit>::Instance();
./classes/DelphesClasses.cc:43:CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
./classes/SortableObject.h:36:class CompBase
./classes/SortableObject.h:39:  virtual ~CompBase() { }
./classes/SortableObject.h:53:  virtual const CompBase *GetCompare() const = 0;
./classes/SortableObject.h:63:class CompE: public CompBase
./classes/SortableObject.h:89:class CompPT: public CompBase
./classes/SortableObject.h:115:class CompMomentumPt: public CompBase
./classes/SortableObject.h:141:class CompET: public CompBase
./classes/SortableObject.h:167:class CompDeltaR: public CompBase
./classes/DelphesClasses.h:160:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:161:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:253:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:254:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:288:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:289:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:321:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:322:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:381:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:382:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:425:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:426:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:454:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:455:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:480:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:481:  const CompBase *GetCompare() const { return fgCompare; }
./classes/DelphesClasses.h:579:  static CompBase *fgCompare; //!
./classes/DelphesClasses.h:580:  const CompBase *GetCompare() const { return fgCompare; }
Binary file ./hepmc2pileup matches

Best,
Li

Hi Li,

after loading libDelphes at interpreter commandline, what does this command return?

gSystem->Load("libDelphes.so");
TClass* c = TClass::GetClass("CompBase");
cout << c << endl;

I think you are using ROOT5. I believe this problem would be solved by ROOT6: did you try it?

Cheers,
Danilo

[quote=“dpiparo”]Hi Li,

after loading libDelphes at interpreter commandline, what does this command return?

gSystem->Load("libDelphes.so");
TClass* c = TClass::GetClass("CompBase");
cout << c << endl;

I think you are using ROOT5. I believe this problem would be solved by ROOT6: did you try it?

Cheers,
Danilo[/quote]

Hi Danilo,

ROOT 5.34/34 (v5-34-34@v5-34-34, Oct 02 2015, 16:30:37 on linuxx8664gcc)

CINT/ROOT C/C++ Interpreter version 5.18.00, July 2, 2010
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
root [0] gSystem->Load("libDelphes.so");
root [1] TClass* c = TClass::GetClass("CompBase");
Warning in <TClass::TClass>: no dictionary for class CompBase is available
root [2] cout << c << endl;
0x14f28e0
root [3] 

This is the result.
And I am using CheckMATE so I have to use root 5. I do not try root 6 yet.

Best,
Li

Hi Li,

this could mean that the class is not properly selected in the selection files.
I see on github github.com/delphes/delphes/blob … sLinkDef.h that the class does not appear in the linkdef. Could you try to add it (#pragma link C++ class CompBase+;), rebuild delphes and try again?

Cheers,
D

[quote=“dpiparo”]Hi Li,

this could mean that the class is not properly selected in the selection files.
I see on github github.com/delphes/delphes/blob … sLinkDef.h that the class does not appear in the linkdef. Could you try to add it (#pragma link C++ class CompBase+;), rebuild delphes and try again?

Cheers,
D[/quote]

Hi,
After adding the code and rebuilt,

ui@ui:~/tools/package/delphes$ root -l
root [0] gSystem->Load("libDelphes.so");
root [1] TClass* c = TClass::GetClass("CompBase");
root [2] cout << c << endl;
0x11c5a40
root [3] 

amazing, No warning. You are so good at coding. Can you tell me the reason of this warning ?

Best,
Li