Looping over tree branches to find index of leading jet

asidley · July 31, 2019, 12:29pm

Hi, extreme newbie here. I am not sure if this is a ROOT problem per se or just a c++ problem, since I’m very new to the language. I am writing a function that goes through a Chain of trees from 5 files, and looks for the index of the entry that passes the tight selection criteria for being the leading jet. This function is written inside the Macro generated by the Makeclass() method (which I have called testMC.C). The code compiles fine, but always returns -1, even though it should return the index of the leading jet (and I am certain that at least one jet passes the selection criteria).

Int_t testMC::JetSelection() {
  //finds the index of leading jet based on tight selection criteria                                                                                                                         
  Int_t jc = -1;
  for (int i=0; i<nJet; i++) {
  fChain->GetEntry(i);

    if ((*jetEta)[i] <= 2.6) {
      if ((*jetNHF)[i] < 0.9 && (*jetNEF)[i] < 0.9 && (*jetNConstituents)[i] > 1 && (*jetMuonEnF)[i] < 0.8 && (*jetCHF)[i] > 0 && (*jetNCharged)[i] > 0 && (*jetCEF)[i] < 0.8) {
        jc = i;
      }
        if (jc != -1) {
          break;
      }
    }
    else if ((*jetEta)[i] > 2.6 && (*jetEta)[i] <= 2.7) {
      if ((*jetNHF)[i] < 0.9 && (*jetNEF)[i] < 0.99 && (*jetMuonEnF)[i] < 0.8 && (*jetNCharged)[i]>0 && (*jetCEF)[i] < 0.8) {
          jc = i;
      }
          if (jc != -1) {
          break;
          }
      }
    else if ((*jetEta)[i] > 2.7 && (*jetEta)[i] <= 3.0) {
      if (((*jetNEF)[i] < 0.99 && (*jetNEF)[i] > 0.2) && (*jetNNeutral)[i] >2) {
              jc = i;
      }
      if (jc != -1) {
              break;
            }
          }
    else if ((*jetEta)[i] > 3.0 && (*jetEta)[i] <= 5.0) {
      if ((*jetNHF)[i] > 0.2 && (*jetNEF)[i] < 0.9 && (*jetNNeutral)[i] > 10) {
                jc = i;
      }
      if (jc!= -1) {
                break;
              }
    }
  }
  return jc;
}

Wile_E_Coyote · July 31, 2019, 12:35pm

Well,

  for (int i=0; i<nJet; i++) {
  fChain->GetEntry(i);

means that you read tree entries 0, 1, ..., (nJet - 1) and I suspect that it is not really what you want.

asidley · July 31, 2019, 12:38pm

Ah you are right. Copy-paste error. But I removed that line and still have the same problem

Wile_E_Coyote · July 31, 2019, 12:49pm

I think you need to attach your full source code.

asidley · July 31, 2019, 12:58pm

This is the full Macro. The header class is just the default generated by MakeClass(), with the directories of my files added.

#define testMC_cxx
#include "testMC.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
Int_t testMC::JetSelection() {
  //finds the index of leading jet based on tight selection criteria                                                                                                                         
  Int_t jc = -1;
  for (int i=0; i<nJet; i++) {
    // fChain->GetEntry(i);                                                                                                                                                                  

    if ((*jetEta)[i] <= 2.6) {
      if ((*jetNHF)[i] < 0.9 && (*jetNEF)[i] < 0.9 && (*jetNConstituents)[i] > 1 && (*jetMuonEnF)[i] < 0.8 && (*jetCHF)[i] > 0 && (*jetNCharged)[i] > 0 && (*jetCEF)[i] < 0.8) {
        jc = i;
      }
        if (jc != -1) {
          break;
      }
    }
    else if ((*jetEta)[i] > 2.6 && (*jetEta)[i] <= 2.7) {
      if ((*jetNHF)[i] < 0.9 && (*jetNEF)[i] < 0.99 && (*jetMuonEnF)[i] < 0.8 && (*jetNCharged)[i]>0 && (*jetCEF)[i] < 0.8) {
          jc = i;
      }
          if (jc != -1) {
          break;
          }
      }
    else if ((*jetEta)[i] > 2.7 && (*jetEta)[i] <= 3.0) {
      if (((*jetNEF)[i] < 0.99 && (*jetNEF)[i] > 0.2) && (*jetNNeutral)[i] >2) {
              jc = i;
      }
      if (jc != -1) {
              break;
            }
          }
    else if ((*jetEta)[i] > 3.0 && (*jetEta)[i] <= 5.0) {
      if ((*jetNHF)[i] > 0.2 && (*jetNEF)[i] < 0.9 && (*jetNNeutral)[i] > 10) {
                jc = i;
      }
      if (jc!= -1) {
                break;
              }
    }
  }
return jc 
}
void testMC::LeadingLoop(Int_t leadindex) {
  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntriesFast();

  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
  }

  TFile f("leadingfile", "CREATE");
  TH1F* lh_jetpt = new TH1F("lh_jetpt", "leading jet pt histogram;jetpt;number of entries", 100, 0, 430);
  TH1F* lh_jeteta = new TH1F("lh_jeteta", "leading jet eta histogram;eta;number of entries", 100, -6, 6);
  TH1F* lh_jetphi = new TH1F("lh_jetphi", "leading jet phi histogram;phi;number of entries", 100, -3.5, 3.5);
  TH1F* lh_phoet = new TH1F("lh_phoet", "leading photon et histogram;Et;number of entries", 100, 0, 430);
  TH1F* lh_phoeta = new TH1F("lh_phoeta", "leading photon eta histogram;eta;number of entries", 100, -3.05, 3.05);
  TH1F* lh_phophi = new TH1F("lh_phophi", "leading photon phi histogram;phi;number of entries", 100, -4, 4);
  TH1F* lh_elept = new TH1F("lh_elept", "leading electron pt histogram;pt;number of entries", 100, 0, 400);
  TH1F* lh_eleeta = new TH1F("lh_eleeta", "leading electron eta histogram;eta;number of entries", 100, -3.5, 3.5);
  TH1F* lh_elephi = new TH1F("lh_elephi", "leading electron phi histogram;phi;number of entries", 100, -3.5, 3.5);                                                                                                                                     
  TH1F* lh_mupt = new TH1F("lh_mupt", "leading muon pt histogram;pt;number of entries", 100, 0, 2000);
  TH1F* lh_mueta = new TH1F("lh_mueta", "leading muon eta histogram;eta;number of entries", 100, 0, 4);
  TH1F* lh_muphi = new TH1F("lh_muphi", "leading muon phi histogram;phi;number of entries", 100, -4, 4);


for (int i=0; i<nentries; i++) {
    fChain->GetEntry(i);
    lh_jetpt->Fill((*jetPt)[leadindex]);
    lh_jeteta->Fill((*jetEta)[leadindex]);
    lh_jetphi->Fill((*jetPhi)[leadindex]);
    lh_phoet->Fill((*phoEt)[leadindex]);
    lh_phoeta->Fill((*phoEta)[leadindex]);
    lh_phophi->Fill((*phoPhi)[leadindex]);
    lh_elept->Fill((*elePt)[leadindex]);
    lh_eleeta->Fill((*eleEta)[leadindex]);
    lh_elephi->Fill((*elePhi)[leadindex]);
    lh_mupt->Fill((*muPt)[leadindex]);
    lh_mueta->Fill((*muEta)[leadindex]);
    lh_muphi->Fill((*muPhi)[leadindex]);

 }
    lh_jetpt->Write("lh_jetpt");
    lh_jeteta->Write("lh_jeteta");
    lh_jetphi->Write("lh_jetphi");
    lh_phoet->Write("lh_phoet");
    lh_phoeta->Write("lh_phoeta");
    lh_phophi->Write("lh_phophi");
    lh_elept->Write("lh_elept");
    lh_eleeta->Write("lh_eleeta");
    lh_elephi->Write("lh_elephi");
    lh_mupt->Write("lh_mupt");
    lh_mueta->Write("lh_mueta");
    lh_muphi->Write("lh_muphi");
  f.Close();
}

Wile_E_Coyote · July 31, 2019, 1:08pm

Try:

void testMC::LeadingLoop() {
  if (fChain == 0) return;
  
  TFile f("leadingfile.root", "RECREATE");
  TH1F* lh_jetpt = new TH1F("lh_jetpt", "leading jet pt histogram;jetpt;number of entries", 100, 0, 430);
  TH1F* lh_jeteta = new TH1F("lh_jeteta", "leading jet eta histogram;eta;number of entries", 100, -6, 6);
  TH1F* lh_jetphi = new TH1F("lh_jetphi", "leading jet phi histogram;phi;number of entries", 100, -3.5, 3.5);
  TH1F* lh_phoet = new TH1F("lh_phoet", "leading photon et histogram;Et;number of entries", 100, 0, 430);
  TH1F* lh_phoeta = new TH1F("lh_phoeta", "leading photon eta histogram;eta;number of entries", 100, -3.05, 3.05);
  TH1F* lh_phophi = new TH1F("lh_phophi", "leading photon phi histogram;phi;number of entries", 100, -4, 4);
  TH1F* lh_elept = new TH1F("lh_elept", "leading electron pt histogram;pt;number of entries", 100, 0, 400);
  TH1F* lh_eleeta = new TH1F("lh_eleeta", "leading electron eta histogram;eta;number of entries", 100, -3.5, 3.5);
  TH1F* lh_elephi = new TH1F("lh_elephi", "leading electron phi histogram;phi;number of entries", 100, -3.5, 3.5);
  TH1F* lh_mupt = new TH1F("lh_mupt", "leading muon pt histogram;pt;number of entries", 100, 0, 2000);
  TH1F* lh_mueta = new TH1F("lh_mueta", "leading muon eta histogram;eta;number of entries", 100, 0, 4);
  TH1F* lh_muphi = new TH1F("lh_muphi", "leading muon phi histogram;phi;number of entries", 100, -4, 4);
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    
    Int_t leadindex = JetSelection();
    if (leadindex != -1) {
      lh_jetpt->Fill((*jetPt)[leadindex]);
      lh_jeteta->Fill((*jetEta)[leadindex]);
      lh_jetphi->Fill((*jetPhi)[leadindex]);
      lh_phoet->Fill((*phoEt)[leadindex]);
      lh_phoeta->Fill((*phoEta)[leadindex]);
      lh_phophi->Fill((*phoPhi)[leadindex]);
      lh_elept->Fill((*elePt)[leadindex]);
      lh_eleeta->Fill((*eleEta)[leadindex]);
      lh_elephi->Fill((*elePhi)[leadindex]);
      lh_mupt->Fill((*muPt)[leadindex]);
      lh_mueta->Fill((*muEta)[leadindex]);
      lh_muphi->Fill((*muPhi)[leadindex]);
    }
  }
  
#if 0 /* 0 or 1 */
  f.cd(); // change the current directory to "f"
  lh_jetpt->Write("lh_jetpt");
  lh_jeteta->Write("lh_jeteta");
  lh_jetphi->Write("lh_jetphi");
  lh_phoet->Write("lh_phoet");
  lh_phoeta->Write("lh_phoeta");
  lh_phophi->Write("lh_phophi");
  lh_elept->Write("lh_elept");
  lh_eleeta->Write("lh_eleeta");
  lh_elephi->Write("lh_elephi");
  lh_mupt->Write("lh_mupt");
  lh_mueta->Write("lh_mueta");
  lh_muphi->Write("lh_muphi");
#else /* 0 or 1 */
  f.Write(); // automatically saves all "owned" histograms
#endif /* 0 or 1 */
  f.Close(); // automatically deletes all "owned" histograms
}

asidley · July 31, 2019, 1:14pm

Calling this function gives this error which I’m not sure what it means…

 *** Break *** segmentation violation



===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
    gdb.printing.register_pretty_printer(gdb.current_objfile(),
    gdb.printing.register_pretty_printer(gdb.current_objfile(),
#0  0x00007f7cfaffda3c in waitpid () from /lib64/libc.so.6
#1  0x00007f7cfaf7bde2 in do_system () from /lib64/libc.so.6
#2  0x00007f7cfbb6ccf4 in TUnixSystem::StackTrace() () from /usr/lib64/root/libCore.so.6.18
#3  0x00007f7cfbb6eb3c in TUnixSystem::DispatchSignals(ESignals) () from /usr/lib64/root/libCore.so.6.18
#4  <signal handler called>
#5  0x00007f7cfc2cca3c in ?? ()
#6  0x0000000000000000 in ?? ()
===========================================================


The lines below might hint at the cause of the crash.
You may get help by asking at the ROOT forum http://root.cern.ch/forum
Only if you are really convinced it is a bug in ROOT then please submit a
report at http://root.cern.ch/bugs Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5  0x00007f7cfc2cca3c in ?? ()
#6  0x0000000000000000 in ?? ()
===========================================================

asidley · July 31, 2019, 1:17pm

Also, that function depends on the JetSelection() function working properly, and it still doesn’t

Wile_E_Coyote · July 31, 2019, 2:03pm

You can try to add some protection, e.g. something like:

    // ...
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    
    // just a precaution ... assuming that all these pointers reference some std::vector
    if ( (!jetPt) || (!jetEta) || (!jetPhi) || (!phoEt) || (!phoEta) || (!phoPhi) || (!elePt) || (!eleEta) || (!elePhi) || (!muPt) || (!muEta) || (!muPhi) ||
         (jetPt->size() < nJet) || (jetEta->size() < nJet) || (jetPhi->size() < nJet) || (phoEt->size() < nJet) || (phoEta->size() < nJet) || (phoPhi->size() < nJet) || (elePt->size() < nJet) || (eleEta->size() < nJet) || (elePhi->size() < nJet) || (muPt->size() < nJet) || (muEta->size() < nJet) || (muPhi->size() < nJet) ||
         (!jetNHF) || (!jetNEF) || (!jetNConstituents) || (!jetMuonEnF) || (!jetCHF) || (!jetNCharged) ||  (!jetCEF) ||  (!jetNNeutral) ||
         (jetNHF->size() < nJet) || (jetNEF->size() < nJet) || (jetNConstituents->size() < nJet) || (jetMuonEnF->size() < nJet) || (jetCHF->size() < nJet) || (jetNCharged->size() < nJet) ||  (jetCEF->size() < nJet) ||  (jetNNeutral->size() < nJet) ) continue;
    
    Int_t leadindex = JetSelection();
    // ...

asidley · August 1, 2019, 12:04pm

Thank you, this actually worked! I’m still confused though, because JetSelection() always returns -1, yet LeadingLoop(), which uses the output of JetSelection() as a vector index, still manages to find the leading jet and create histograms with its variables…