A problem with LHEF tree and "if" operator

ROOT
1

Good day!
I have no ideas anymore what to do with that…
I’m using the following scrpit.

[code]#define test_JHU_4e_cxx
#include “test_JHU_4e.h”
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void test_JHU_4e::Loop()
{
// In a ROOT session, you can do:
// Root > .L test_JHU_4e.C
// Root > test_JHU_4e t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries
//

// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, Insert statements like:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus(“branchname”,1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
//by b_branchname->GetEntry(ientry); //read only this branch
gROOT->ProcessLine(".x style.C");
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
//gStyle->SetOptStat(00111111);

hfile = new TFile(“hist_test_JHU_4e.root”, “RECREATE”, “ROOT histograms”);

//TH1::SetDefaultSumw2(kTRUE);

TH1F *h_Higgs_Minv = new TH1F(“h_Higgs_Minv”, “pp->x0->ZZ->4l 8 TeV”, 100, 0.0, 500.0);
TH1F *h_4l_PT = new TH1F(“h_4l_PT”, “pp->x0->ZZ->4l 8 TeV”, 100, 0.0, 500.0);
TH1F *h_4l_PT_log = new TH1F(“h_4l_PT_log”, “pp->x0->ZZ->4l 8 TeV”, 100, 0.0, 500.0);
TH1F *h_Z1_Minv = new TH1F(“h_Z1_Minv”, “pp->x0->ZZ->4l 8 TeV”, 100, 0.0, 150.0);
TH1F *h_Z2_Minv = new TH1F(“h_Z2_Minv”, “pp->x0->ZZ->4l 8 TeV”, 100, 0.0, 150.0);
TH1F *h_cost_star = new TH1F(“h_cost_star”, “pp->x0->ZZ->4l 8 TeV”, 15, -1.0, 1.0);
TH1F *h_cost1 = new TH1F(“h_cost1”, “pp->x0->ZZ->4l 8 TeV”, 15, -1.0, 1.0);
TH1F *h_cost2 = new TH1F(“h_cost2”, “pp->x0->ZZ->4l 8 TeV”, 15, -1.0, 1.0);
TH1F *h_fi = new TH1F(“h_fi”, “pp->x0->ZZ->4l 8 TeV”, 15, -3.14, 3.14);
TH1F *h_fi1 = new TH1F(“h_fi1”, “pp->x0->ZZ->4l 8 TeV”, 15, -3.14, 3.14);

h_Z1_Minv->SetLineColor(kRed);
h_Z1_Minv->GetXaxis()->SetTitle(“Z1 Minv [GeV]”);
h_Z1_Minv->GetYaxis()->SetTitle(“Normalized entries”);

h_Z2_Minv->SetLineColor(kRed);
h_Z2_Minv->GetXaxis()->SetTitle(“Z2 Minv [GeV]”);
h_Z2_Minv->GetYaxis()->SetTitle(“Normalized entries”);

h_Higgs_Minv->SetLineColor(kRed);
h_Higgs_Minv->GetXaxis()->SetTitle(“Higgs Minv [GeV]”);
h_Higgs_Minv->GetYaxis()->SetTitle(“Normalized entries”);

h_4l_PT->SetLineColor(kRed);
h_4l_PT_log->SetLineColor(kRed);
h_4l_PT->GetXaxis()->SetTitle(“Higgs PT [GeV]”);
h_4l_PT_log->GetXaxis()->SetTitle(“Higgs PT [GeV]”);
h_4l_PT->GetYaxis()->SetTitle(“Normalized entries”);
h_4l_PT_log->GetYaxis()->SetTitle(“Events / 5 GeV”);

h_cost_star->SetLineColor(kBlue);
h_cost_star->GetXaxis()->SetTitle(“cos(#theta^{*})”);
h_cost_star->GetYaxis()->SetTitle(“Normalized entries”);

h_cost1->SetLineColor(kBlue);
h_cost1->GetXaxis()->SetTitle(“cos(#theta_{1})”);
h_cost1->GetYaxis()->SetTitle(“Normalized entries”);

h_cost2->SetLineColor(kBlue);
h_cost2->GetXaxis()->SetTitle(“cos(#theta_{2})”);
h_cost2->GetYaxis()->SetTitle(“Normalized entries”);

h_fi->SetLineColor(kBlue);
h_fi->GetXaxis()->SetTitle("#Phi");
h_fi->GetYaxis()->SetTitle(“Normalized entries”);

h_fi1->SetLineColor(kBlue);
h_fi1->GetXaxis()->SetTitle("#Phi_{1}");
h_fi1->GetYaxis()->SetTitle(“Normalized entries”);

double cost_star;
double cost1, cost2;
double fi, fi1, cosfi, cosfi1;

int i;

double MZ = 91.1876;
double Minv, Z1_Minv, Z2_Minv, Z3_Minv, Z4_Minv;
double scale, d1, d2, d3, d4;
double pt1=0, pt2=0, pt3=0;
int n_pt1, n_pt2, n_pt3;

TVector3 nz;
nz[0] = 0.0; nz[1] = 0.0; nz[2] = 1.0;
TVector3 z1_rfr_Z2, parton, z1, z2, v1p, v2p, v3p, v4p, n1p, n2p, nscp, z2_rfr_Z1;

TLorentzVector Z1, Z2, Z3, Z4, v1, v2, v3, v4, vH, q, vZ_tmp, v_tmp, Z2_rfr_Z1, Z1_rfr_Z2;

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;
if(!(jentry%1000)) cout << jentry << endl;

 v1.SetXYZT(Particle_Px[5], Particle_Py[5], Particle_Pz[5], Particle_E[5]);
 v2.SetXYZT(Particle_Px[6], Particle_Py[6], Particle_Pz[6], Particle_E[6]);
 v3.SetXYZT(Particle_Px[7], Particle_Py[7], Particle_Pz[7], Particle_E[7]);
 v4.SetXYZT(Particle_Px[8], Particle_Py[8], Particle_Pz[8], Particle_E[8]);
//cout << Particle_PID[5] << Particle_PID[6] << Particle_PID[7] << Particle_PID[8] << endl; 
//cin.get();

	 Z1 = v1 + v2;
     Z2 = v3 + v4;
     Z3 = v1 + v4;
     Z4 = v2 + v3;
     Z1_Minv = Z1.Mag();
     Z2_Minv = Z2.Mag();
     Z3_Minv = Z3.Mag();
     Z4_Minv = Z4.Mag();

     d1 = TMath::Abs(Z1_Minv - MZ);
     d2 = TMath::Abs(Z2_Minv - MZ);
     d3 = TMath::Abs(Z3_Minv - MZ);
     d4 = TMath::Abs(Z4_Minv - MZ);

     if(TMath::Min(d3, d4) < TMath::Min(d1, d2)) {

        Z1 = Z3;
        Z2 = Z4;
        d1 = d3;
        d2 = d4;

		v_tmp = v2;
	    v2 = v4;
	    v4 = v_tmp;
     }


     if(d2 < d1) {

        vZ_tmp = Z2;
        Z2 = Z1;
        Z1 = vZ_tmp;

	    v_tmp = v1;
	    v1 = v3;
	    v3 = v_tmp;
	    v_tmp = v2;
	    v2 = v4;
	    v4 = v_tmp;
}
  Z1_Minv = Z1.Mag();
  Z2_Minv = Z2.Mag();
  h_Z1_Minv->Fill(Z1_Minv);
  h_Z2_Minv->Fill(Z2_Minv);

  vH = Z1 + Z2;

  q.SetPxPyPzE( 0, 0, ( vH.E() + vH.Pz() ) / 2., ( vH.E() + vH.Pz() ) / 2. );
  q.Boost( -( vH.BoostVector() ) );

  parton = q.Vect().Unit();

  Z1.Boost( -( vH.BoostVector() ) ); 
  Z2.Boost( -( vH.BoostVector() ) );

  z1 = Z1.Vect().Unit();
  z2 = Z2.Vect().Unit();

  cost_star = TMath::Cos( parton.Angle( z1 ) );

  TLorentzVector vv1, vv2, vv3, vv4;
  vv1 = v1;
  vv2 = v2;
  vv3 = v3;
  vv4 = v4;
  vv1.Boost( -( vH.BoostVector() ) );
  vv2.Boost( -( vH.BoostVector() ) );
  vv3.Boost( -( vH.BoostVector() ) );
  vv4.Boost( -( vH.BoostVector() ) );

  v1p = vv1.Vect();
  v2p = vv2.Vect();
  v3p = vv3.Vect();
  v4p = vv4.Vect();

  n1p = v1p.Cross( v2p ).Unit();
  n2p = v3p.Cross( v4p ).Unit();
  nscp = nz.Cross( z1 ).Unit();
  fi = ( z1.Dot( n1p.Cross( n2p ) ) /
		  TMath::Abs( z1.Dot( n1p.Cross( n2p ) ) ) *
		  TMath::ACos( -n1p.Dot( n2p ) ) );
  fi1 = ( z1.Dot( n1p.Cross( nscp ) ) /
		   TMath::Abs( z1.Dot( n1p.Cross( nscp ) ) ) *
		   TMath::ACos( n1p.Dot( nscp ) ) );

  Z2_rfr_Z1 = Z2;
  Z2_rfr_Z1.Boost( -( Z1.BoostVector() ) ); 
  z2_rfr_Z1 = Z2_rfr_Z1.Vect();

  Z1_rfr_Z2 = Z1;
  Z1_rfr_Z2.Boost( -( Z2.BoostVector() ) );
  z1_rfr_Z2 = Z1_rfr_Z2.Vect();
  vv1.Boost( -( Z1.BoostVector() ) );
  vv3.Boost( -( Z2.BoostVector() ) );


  cost1 = - ( z2_rfr_Z1.Dot( vv1.Vect() ) /
		     TMath::Abs( z2_rfr_Z1.Mag() * vv1.Vect().Mag() ) );
  cost2 = - ( z1_rfr_Z2.Dot( vv3.Vect() ) /
		     TMath::Abs( z1_rfr_Z2.Mag() * vv3.Vect().Mag() ) );

	 h_Higgs_Minv->Fill(vH.M());
	 h_4l_PT->Fill(vH.Pt());
     h_4l_PT_log->Fill(vH.Pt());
	 h_cost_star->Fill(cost_star);
	 h_cost1->Fill(cost1);
	 h_cost2->Fill(cost2);
	 h_fi1->Fill(fi1);
     h_fi->Fill(fi);
 };
 
 	h_Z1_Minv->Sumw2();
    h_Z2_Minv->Sumw2();
    h_cost_star->Sumw2();
	h_cost1->Sumw2();
	h_cost2->Sumw2();
	h_fi->Sumw2();
    h_fi1->Sumw2();
	
 scale = 1.0/h_Higgs_Minv->GetEntries();
 h_Higgs_Minv->Scale(scale);
 h_Higgs_Minv->SetMinimum(0);
 scale = 1.0/h_4l_PT->GetEntries();
 h_4l_PT->Scale(scale);
 h_4l_PT->SetMinimum(0);
 scale = 1.0/h_Z1_Minv->GetEntries();
 h_Z1_Minv->Scale(scale);
 h_Z1_Minv->SetMinimum(0);
 scale = 1.0/h_Z2_Minv->GetEntries();
 h_Z2_Minv->Scale(scale);
 h_Z2_Minv->SetMinimum(0);
 scale = 1.0/h_cost_star->GetEntries();
 h_cost_star->Scale(scale);
 h_cost_star->SetMinimum(0);
 h_cost1->Scale(scale);
 h_cost1->SetMinimum(0);
 scale = 1.0/h_cost2->GetEntries();
 h_cost2->Scale(scale);
 h_cost2->SetMinimum(0);
 scale = 1.0/h_fi1->GetEntries();
 h_fi1->Scale(scale);
 h_fi1->SetMinimum(0);
 scale = 1.0/h_fi->GetEntries();
 h_fi->Scale(scale);
 h_fi->SetMinimum(0);

 hfile->Write();

} //end main loop
[/code]
This script plots some cinematic distribusions. Input root file contains LHEF standart tree. Runtime is about ~ 60 seconds.
Now I want to make some particle selection and add these lines at the begining of the Loop:

i=Particle_PID[5]+Particle_PID[7]; if(i/4==6) continue;
Now this scrpit have runtime about 60 hours!
This simple selection just ruin it! Also I have increased calculation time per event with increasing number of running event. Maybe it’s some king of overflow, I don’t know. I’ve tried a thousand variation of “if” operators and noone works. The problem cause, I suppose, when I want to decline more that one possible conditions.
Also if I try

i=Particle_PID[5]+Particle_PID[6]; if(i/4==6) continue;
The code works fine. Particle_PID[5]=11 or 13 (lepton or muon), Particle_PID[6]=-11 or -13 (antilepton or antimuon from the same Z boson!). Particle_PID[7] is the lepton of muon from the second Z-boson.

Help me please, this is a very annoying long-story bug, I can provide any information you need!

Best Regards,
Nikita

2

Hi,

I am afraid that the question is not related to ROOT.
You could try adding some output or reformulate this condition to kick off debugging.

Best,
Danilo

3

Thanks for the reply, Danilo
Can you, please, tell me some more about debug? How can I use it at ROOT?
I’m using ROOT v5-34-00 at Windows 7x64.

Best regards,
Nikita

4

Hi,

there are plenty of debuggers out there, pretty good ones by Microsoft actually. For this case a printout on screen (cout) would do.

Cheers,
Danilo

5

Hi Nikita,

could you please provide some kind of minimum code snippet reproducing your problem?
Or at the very least please tell us how to run your script from your original post and where one could take needed input files / headers (if any).

6

[quote=“yus”]Hi Nikita,

could you please provide some kind of minimum code snippet reproducing your problem?
Or at the very least please tell us how to run your script from your original post and where one could take needed input files / headers (if any).[/quote]
I’ll try to provide some “brief” version of my code:

void program::Loop()
{
#defining input file, creating histograms, creating variables etc#

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;

//reading some parameters
     v1.SetXYZT(Particle_Px[5], Particle_Py[5], Particle_Pz[5], Particle_E[5]);
     v2.SetXYZT(Particle_Px[6], Particle_Py[6], Particle_Pz[6], Particle_E[6]);
     v3.SetXYZT(Particle_Px[7], Particle_Py[7], Particle_Pz[7], Particle_E[7]);
     v4.SetXYZT(Particle_Px[8], Particle_Py[8], Particle_Pz[8], Particle_E[8]);

#calculating some values#

#filling histograms#
};
#normalizing histograms and adding errorbars#
}

This code works fine and it’s very fast. But if I want to make some parameters selection it works awful (VERY slow).
I mean if I add some logical operators in the beggining of the loop (if, else, switch etc) which needs to decline 2 or more variants the calculating time increases at ~ 1000 times.
But, if I add the selection at the fitting stage at the end of the loop, like that:

void program::Loop()
{
#defining input file, creating histograms, creating variables etc#

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;

//reading some parameters
     v1.SetXYZT(Particle_Px[5], Particle_Py[5], Particle_Pz[5], Particle_E[5]);
     v2.SetXYZT(Particle_Px[6], Particle_Py[6], Particle_Pz[6], Particle_E[6]);
     v3.SetXYZT(Particle_Px[7], Particle_Py[7], Particle_Pz[7], Particle_E[7]);
     v4.SetXYZT(Particle_Px[8], Particle_Py[8], Particle_Pz[8], Particle_E[8]);

#calculating some values#

if(selection id passed){
#filling histograms#
};

};
#normalizing histograms and adding errorbars#
}

all works fine.
So, maybe the problem is caused by huge amount of code which logical operators must handle during the loop?

7

Hi Nikita,

it’s still very difficult to tell without any code snippet that we can actually run. If you could slightly modify your snippet from last post adding input files (which should be in some public folder on e.g. lxplus) and telling us exactly how to run your code, that would be great.