Dear ROOTers,
I am using ROOT 5.34.24 on Mac OSX 10.11.3 to analyse different *.root event files crated by MadGraph 5 2.2.3 and found that my code is running very slow if I change “arbitrary” pieces of code. I am able to remove “faulty” or unused parts of the code and then my machine processes roughly 1k events per second. However, if I add certain code (even in untriggered if-loops), the processing speed reduces to roughly 5 events per second.
Just to give you one of many examples where this happens (see also code of files at the bottom of this post): If I add the line
anywhere after the particle loop in ang_distr.C ended and with all terms of the sum being well defined integers, everything works fine. However, just adding any other integer, e.g.
massively slows the code down.
Another example: If I add cout << "hello" << endl; just anywhere after the particle loop, the code also slows down. Uncommenting some of the pieces in ang_distr.C yields the same effect.
I was not able to track down the cause of this issue, but I think (and hope) I might just be overseeing something obvious.
Because I cannot give a more specific description of this problem, I copied the full (but minimal) working code below in a state, where it is fast.
My code consists of 4 files: the main *.C file ang_distr.C, the corresponding *.h header file ang_distr.h and two header files holding some functions I defined, angles.h and my_func.h.
Thank you very much for any help and apologies for this long post!
Cheers,
Patrick
Edit:
I thought it might be helpful to summarise the structure of the main routine Bool_t ang_distr::Process(Long64_t entry) {} of the file ang_distr.C to give you a better overview. The routine is organised as follows (pseudo code):
Reset histograms and define most variables
if(Particle_size > 0)
{
for (int iii=0; iii < Particle_size; iii++) { // for particles in event
{
find incoming particles
find outgoing particles
}
recreate final_state // slowing down appears here...
recreate intermediate particles // ... and/or here...
choose scenario to study (no quarks, 2 quarks, ...) // ...and/or here
} ang_distr.C:
#define ang_distr_cxx
#include "angles.h"
#include "my_func.h"
#include "ang_distr.h"
#include <iostream>
#include <cmath>
#include <TH2.h>
#include <TStyle.h>
#include <TLorentzVector.h>
#include <TCanvas.h>
#include <TClonesArray.h>
void ang_distr::Begin(TTree * /*tree*/) {
CosPrimeCut = 0.5;
N_Events = 1;
hCphistar->Reset();
hCphiprime->Reset();
hCallinvmass->Reset();
hCWinvmass->Reset();
TString option = GetOption();
}
void ang_distr::SlaveBegin(TTree * /*tree*/) {
TString option = GetOption();
}
Bool_t ang_distr::Process(Long64_t entry) {
GetEntry(entry);
final_state = 0;
N_incoming = 0; // counts incoming particles
N_outgoing = 0; // counts all outgoing particles
N_out_q = 0; // counts outgoing quarks only
N_out_Z = 0; // counts outgoing Z bosons only
N_out_Wp = 0; // counts outgoing W+ bosons only
em_out = 0; // dummy index to count if there is an outgoing electron
ep_out = 0; // dummy index to count if there is an outgoing positron
mum_out = 0; // dummy index to count if there is an outgoing muon
mup_out = 0; // dummy index to count if there is an outgoing anti-muon
ve_out = 0; // dummy index to count if there is an outgoing neutrino_e
a_ve_out = 0; // dummy index to count if there is an outgoing anti_neutrino_e
vm_out = 0; // dummy index to count if there is an outgoing neutrino_mu
a_vm_out = 0; // dummy index to count if there is an outgoing anti_neutrino_mu
Ecm = 0.; // initialise the center of mass energy
if(Particle_size > 0) { // Particle size start
for (int iii=0; iii < Particle_size; iii++) { // Particle loop start
if (Particle_Status[iii] == -1) { // find incoming particle loop start
switch (N_incoming) { // set incoming 4-momenta start
case 0: incoming1.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_incoming = 1;
break;
case 1: incoming2.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_incoming = 2;
break;
default: std::cout << "Too many incoming particles." << std::endl;
return 0;
} // set incoming 4-momenta end
} // find incoming particle loop end
if (Particle_Status[iii] == 1) { // find outgoing particle loop start
switch (Particle_PID[iii]) { // set outgoing 4-momenta start
case 11: electron.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
em_out = 1;
break;
case -11: positron.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
ep_out = 2;
break;
case 12: neutrino_e.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
ve_out = 3;
break;
case -12: anti_neutrino_e.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
a_ve_out = 4;
break;
case 13: muon.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
mum_out = 5;
break;
case -13: anti_muon.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
mup_out = 6;
break;
case 14: neutrino_mu.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
vm_out = 10;
break;
case -14: anti_neutrino_mu.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
a_vm_out = 8;
break;
case 23:
switch (N_out_Z) {
case 0: outgoing_Z1.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_Z++;
N_outgoing++;
break;
case 1: outgoing_Z2.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_Z++;
N_outgoing++;
break;
default: cout << "Error in outgoing Z assignments." << endl;
return 0;
}
break;
case 24:
switch (N_out_Wp) {
case 0: outgoing_Wp1.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_Wp++;
N_outgoing++;
break;
case 1: outgoing_Wp2.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_Wp++;
N_outgoing++;
break;
default: cout << "Error in outgoing W+ assignments." << endl;
return 0;
}
break;
case -24: outgoing_Wm.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_outgoing++;
break;
case 1: case 2: case 3: case 4: case 5: case 6: case -1: case -2: case -3: case -4: case -5: case -6:
switch (N_out_q) {
case 0: quark1.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_q++;
N_outgoing++;
break;
case 1: quark2.SetPxPyPzE(Particle_Px[iii],Particle_Py[iii],Particle_Pz[iii],Particle_E[iii]);
N_out_q++;
N_outgoing++;
break;
default: cout << "Error in outgoing quark assignments." << endl;
return 0;
}
break;
default: cout << "Outgoing particle neither electron, muon, neutrino, Z, W+- nor quark." << endl;
return 0;
} // set outgoing 4-momenta end
} // find outgoing particle loop end
} // Particle loop end
// All TLorentzVectors set. Now on to the calculations
// Recreate the final state (not yet optimised) and assign the Parent TLorentzVectors
//final_state = em_out + ep_out + ve_out + a_ve_out + mum_out + mup_out + vm_out + a_vm_out;
final_state = 14;
int abc = em_out + ep_out + ve_out + a_ve_out + em_out; // adding any other integer variable slows it down
//cout << "hello" << endl; // uncommenting this also slows down the code
switch (final_state) {
case 14: // w+w- > zz > e+ e- mu+ mu-
Parent1 = positron + electron;
Parent2 = muon + anti_muon;
cms = Parent1 + Parent2;
Ecm = cms.M();
break;
case 16: // w+z > w+z > e+ ve mu+ mu-
Parent1 = positron + neutrino_e;
Parent2 = muon + anti_muon;
cms = Parent1 + Parent2;
Ecm = cms.M();
break;
case 18: // w+w- > w+w- > e+ ve mu- vm~
Parent1 = positron + neutrino_e;
Parent2 = muon + anti_neutrino_mu;
cms = Parent1 + Parent2;
Ecm = cms.M();
break;
case 21: // w+w+ > w+w+ > e+ ve mu+ vm
Parent1 = positron + neutrino_e;
Parent2 = anti_muon + neutrino_mu;
cms = Parent1 + Parent2;
Ecm = cms.M();
break;
default: cout << "Error in final state reconstruction." << endl;
return 0;
}
// Find the initial intermediate particle from "incoming - outgoing" (the first pair of W's from pp > (W+W- > ZZ) > jj ZZ)
if (N_out_q == 0) {
intermediate1 = incoming1;
intermediate2 = incoming2;
}
else if ( (incoming1.Angle(quark1.Vect()) < incoming1.Angle(quark2.Vect())) && (N_out_q > 0) ) {
intermediate1 = incoming1 - quark1;
intermediate2 = incoming2 - quark2;
}
else {
intermediate1 = incoming1 - quark2;
intermediate2 = incoming2 - quark1;
}
// Incoming and outgoing particles set, now calculate cos(theta_v) and cos(theta*) depending on the final states
if (N_out_q == 2) {
cout << "q2-Event No. " << N_Events << endl;
//if ( (abs(calc_cos_v(Parent1, intermediate1, cms)) < CosPrimeCut) && (abs(quark1.PseudoRapidity() - quark2.PseudoRapidity()) > 4.) && (quark1.E() > 300.) && (quark2.E() > 300.) && (quark1.Pt() > 30.) && (quark2.Pt() > 30.) && (quark1.PseudoRapidity() < 4.5) && (quark2.PseudoRapidity() < 4.5) && is_near_Z_res(electron,positron) && is_near_Z_res(muon,anti_muon) && (Ecm > 0) ) {
if ( N_out_q == 2//(abs(calc_cos_v(Parent1, intermediate1, cms)) < CosPrimeCut) &&
//(abs(quark1.PseudoRapidity() - quark2.PseudoRapidity()) > 4.) &&
//(quark1.E() > 300.) && (quark2.E() > 300.) && (quark1.Pt() > 30.) && (quark2.Pt() > 30.) &&
//(quark1.PseudoRapidity() < 4.5) && (quark2.PseudoRapidity() < 4.5) &&
//is_near_Z_res(electron,positron) && is_near_Z_res(muon,anti_muon) &&
//(Ecm > 0)
) {
//hCphiprime->Fill(calc_cos_v(Parent1, intermediate1, cms));
//hCphistar->Fill(calc_cos_star(electron, positron));
hCallinvmass->Fill(cms.M());
hCWinvmass->Fill(intermediate1.M());
}
}
else if (N_out_q == 0) {
if (N_outgoing == 2) {
cout << "2-Event No. " << N_Events << endl;
if ( (abs(calc_cos_v(outgoing_Z1, intermediate1, cms)) < CosPrimeCut) && (Ecm > 0) ) {
hCphiprime->Fill(calc_cos_v(outgoing_Z1, intermediate1, cms));
hCallinvmass->Fill(cms.M());
}
}
else if (N_outgoing == 4) {
cout << "4-Event No. " << N_Events << endl;
if ( (abs(calc_cos_v(Parent1, intermediate1, cms)) < CosPrimeCut) && (Ecm > 0) ) {
hCphiprime->Fill(calc_cos_v(Parent1, intermediate1, cms));
hCphistar->Fill(calc_cos_star(muon, anti_muon));
hCallinvmass->Fill((cms).M());
}
}
else {cout << "Error: Wrong amount of outgoing particles in zero quark case." << endl;}
}
else {cout << "Error: Wrong number of quarks (should be 0 or 2, is " << N_out_q << ")." << endl;}
} // Particle size end
N_Events++;
return kTRUE;
}
void ang_distr::SlaveTerminate() {
}
void ang_distr::Terminate() {
ofstream myfile;
//myfile.open ("surviving_events.txt", ios_base::app);
//myfile << N_surv_events << " out of " << N_Events << " events survived the cut." << std::endl;
//myfile.close();
TCanvas *can1 = new TCanvas("can1");
can1->cd();
hCphiprime->Draw();
can1->SaveAs("data/rawdata/DELtheta_prime-wpwm_to_zz_to_epemmupmum-a0-0.5-0.C");
TCanvas *can2 = new TCanvas("can2");
can2->cd();
hCphistar->Draw();
can2->SaveAs("data/rawdata/DELtheta_star-wpwm_to_zz_to_epemmupmum-a0-0.5-0.C");
TCanvas *can3 = new TCanvas("can3");
can3->cd();
hCallinvmass->Draw();
can3->SaveAs("data/rawdata/DELall_inv-wpwm_to_zz_to_epemmupmum-a0-0.5-0.C");
TCanvas *can4 = new TCanvas("can4");
can4->cd();
hCWinvmass->Draw();
can4->SaveAs("data/rawdata/DELW_inv-wpwm_to_zz_to_epemmupmum-a0-0.5-0.C");
}ang_distr.h
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Wed Apr 13 18:40:13 2016 by ROOT version 5.34/34
// from TTree LHEF/Analysis tree
// found on file: w+w-_to_zz-a=1.root
//////////////////////////////////////////////////////////
#ifndef ang_distr_h
#define ang_distr_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TLegend.h>
#include <TTree.h>
#include <TBranch.h>
#include <TLorentzVector.h>
// Header file for the classes stored in the TTree if any.
#include <TClonesArray.h>
#include <TObject.h>
// Fixed size dimensions of array or collections stored in the TTree if any.
TH1F *hCphistar = new TH1F("hstar", "hCphistar", 20, -1, 1);
TH1F *hCphiprime = new TH1F("hprime", "hCphiprime", 20, -1, 1);
TH1F *hCallinvmass = new TH1F("hallinvmass", "hCallinvmass", 20, 0, 1200);
TH1F *hCWinvmass = new TH1F("hWinvmass", "hCWinvmass", 20, -600, 0);
//TH1F *hCphistar = new TH1F("hCphistar", "hCphistar", 20, -1, 1);
//TH1F *hCphiprime = new TH1F("hCphiprime", "hCphiprime", 20, -1, 1);
//TH1F *hCallinvmass = new TH1F("hCallinvmass", "hCallinvmass", 20, 0, 1200);
//TH1F *hCWinvmass = new TH1F("hCWinvmass", "hCWinvmass", 20, -600, 0);
Double_t CosPrimeCut; // Cut on cos(theta_V)
Double_t Ecm; // Center of Mass energy
Int_t final_state; // final state indicator
Int_t N_Events; // total number of events of the *.root file
Int_t N_surv_events; // number of events surviving the applied cuts
Int_t N_incoming; // counts incoming particles
Int_t N_outgoing; // counts all outgoing particles
Int_t N_out_q; // counts outgoing quarks only
Int_t N_out_Z; // counts outgoing Z bosons only
Int_t N_out_Wp; // counts outgoing W+ bosons only
Int_t em_out; // dummy index to count if there is an outgoing electron
Int_t ep_out; // dummy index to count if there is an outgoing positron
Int_t mum_out; // dummy index to count if there is an outgoing muon
Int_t mup_out; // dummy index to count if there is an outgoing anti-muon
Int_t ve_out; // dummy index to count if there is an outgoing neutrino_e
Int_t a_ve_out; // dummy index to count if there is an outgoing anti_neutrino_e
Int_t vm_out; // dummy index to count if there is an outgoing neutrino_mu
Int_t a_vm_out; // dummy index to count if there is an outgoing anti_neutrino_mu
TLorentzVector outgoing_Z1;
TLorentzVector outgoing_Z2;
TLorentzVector outgoing_Wp1;
TLorentzVector outgoing_Wp2;
TLorentzVector outgoing_Wm;
TLorentzVector outgoing_gen1;
TLorentzVector outgoing_gen2;
TLorentzVector outgoing_gen3;
TLorentzVector quark1;
TLorentzVector quark2;
TLorentzVector electron;
TLorentzVector positron;
TLorentzVector muon;
TLorentzVector anti_muon;
TLorentzVector neutrino_e;
TLorentzVector anti_neutrino_e;
TLorentzVector neutrino_mu;
TLorentzVector anti_neutrino_mu;
TLorentzVector incoming1;
TLorentzVector incoming2;
TLorentzVector intermediate1;
TLorentzVector intermediate2;
TLorentzVector Parent1;
TLorentzVector Parent2;
TLorentzVector cms;
// Fixed size dimensions of array or collections stored in the TTree if any.
const Int_t kMaxEvent = 1;
const Int_t kMaxRwgt = 1;
const Int_t kMaxParticle = 20;
class ang_distr : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Declaration of leaf types
Int_t Event_;
UInt_t Event_fUniqueID[kMaxEvent]; //[Event_]
UInt_t Event_fBits[kMaxEvent]; //[Event_]
Long64_t Event_Number[kMaxEvent]; //[Event_]
Int_t Event_Nparticles[kMaxEvent]; //[Event_]
Int_t Event_ProcessID[kMaxEvent]; //[Event_]
Double_t Event_Weight[kMaxEvent]; //[Event_]
Double_t Event_ScalePDF[kMaxEvent]; //[Event_]
Double_t Event_CouplingQED[kMaxEvent]; //[Event_]
Double_t Event_CouplingQCD[kMaxEvent]; //[Event_]
Int_t Event_size;
Int_t Rwgt_;
UInt_t Rwgt_fUniqueID[kMaxRwgt]; //[Rwgt_]
UInt_t Rwgt_fBits[kMaxRwgt]; //[Rwgt_]
Double_t Rwgt_Weight[kMaxRwgt]; //[Rwgt_]
Int_t Rwgt_size;
Int_t Particle_;
UInt_t Particle_fUniqueID[kMaxParticle]; //[Particle_]
UInt_t Particle_fBits[kMaxParticle]; //[Particle_]
Int_t Particle_PID[kMaxParticle]; //[Particle_]
Int_t Particle_Status[kMaxParticle]; //[Particle_]
Int_t Particle_Mother1[kMaxParticle]; //[Particle_]
Int_t Particle_Mother2[kMaxParticle]; //[Particle_]
Int_t Particle_ColorLine1[kMaxParticle]; //[Particle_]
Int_t Particle_ColorLine2[kMaxParticle]; //[Particle_]
Double_t Particle_Px[kMaxParticle]; //[Particle_]
Double_t Particle_Py[kMaxParticle]; //[Particle_]
Double_t Particle_Pz[kMaxParticle]; //[Particle_]
Double_t Particle_E[kMaxParticle]; //[Particle_]
Double_t Particle_M[kMaxParticle]; //[Particle_]
Double_t Particle_PT[kMaxParticle]; //[Particle_]
Double_t Particle_Eta[kMaxParticle]; //[Particle_]
Double_t Particle_Phi[kMaxParticle]; //[Particle_]
Double_t Particle_Rapidity[kMaxParticle]; //[Particle_]
Double_t Particle_LifeTime[kMaxParticle]; //[Particle_]
Double_t Particle_Spin[kMaxParticle]; //[Particle_]
Int_t Particle_size;
// List of branches
TBranch *b_Event_; //!
TBranch *b_Event_fUniqueID; //!
TBranch *b_Event_fBits; //!
TBranch *b_Event_Number; //!
TBranch *b_Event_Nparticles; //!
TBranch *b_Event_ProcessID; //!
TBranch *b_Event_Weight; //!
TBranch *b_Event_ScalePDF; //!
TBranch *b_Event_CouplingQED; //!
TBranch *b_Event_CouplingQCD; //!
TBranch *b_Event_size; //!
TBranch *b_Rwgt_; //!
TBranch *b_Rwgt_fUniqueID; //!
TBranch *b_Rwgt_fBits; //!
TBranch *b_Rwgt_Weight; //!
TBranch *b_Rwgt_size; //!
TBranch *b_Particle_; //!
TBranch *b_Particle_fUniqueID; //!
TBranch *b_Particle_fBits; //!
TBranch *b_Particle_PID; //!
TBranch *b_Particle_Status; //!
TBranch *b_Particle_Mother1; //!
TBranch *b_Particle_Mother2; //!
TBranch *b_Particle_ColorLine1; //!
TBranch *b_Particle_ColorLine2; //!
TBranch *b_Particle_Px; //!
TBranch *b_Particle_Py; //!
TBranch *b_Particle_Pz; //!
TBranch *b_Particle_E; //!
TBranch *b_Particle_M; //!
TBranch *b_Particle_PT; //!
TBranch *b_Particle_Eta; //!
TBranch *b_Particle_Phi; //!
TBranch *b_Particle_Rapidity; //!
TBranch *b_Particle_LifeTime; //!
TBranch *b_Particle_Spin; //!
TBranch *b_Particle_size; //!
ang_distr(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~ang_distr() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(ang_distr,0);
};
#endif
#ifdef ang_distr_cxx
void ang_distr::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("Event", &Event_, &b_Event_);
fChain->SetBranchAddress("Event.fUniqueID", Event_fUniqueID, &b_Event_fUniqueID);
fChain->SetBranchAddress("Event.fBits", Event_fBits, &b_Event_fBits);
fChain->SetBranchAddress("Event.Number", Event_Number, &b_Event_Number);
fChain->SetBranchAddress("Event.Nparticles", Event_Nparticles, &b_Event_Nparticles);
fChain->SetBranchAddress("Event.ProcessID", Event_ProcessID, &b_Event_ProcessID);
fChain->SetBranchAddress("Event.Weight", Event_Weight, &b_Event_Weight);
fChain->SetBranchAddress("Event.ScalePDF", Event_ScalePDF, &b_Event_ScalePDF);
fChain->SetBranchAddress("Event.CouplingQED", Event_CouplingQED, &b_Event_CouplingQED);
fChain->SetBranchAddress("Event.CouplingQCD", Event_CouplingQCD, &b_Event_CouplingQCD);
fChain->SetBranchAddress("Event_size", &Event_size, &b_Event_size);
fChain->SetBranchAddress("Rwgt", &Rwgt_, &b_Rwgt_);
fChain->SetBranchAddress("Rwgt.fUniqueID", &Rwgt_fUniqueID, &b_Rwgt_fUniqueID);
fChain->SetBranchAddress("Rwgt.fBits", &Rwgt_fBits, &b_Rwgt_fBits);
fChain->SetBranchAddress("Rwgt.Weight", &Rwgt_Weight, &b_Rwgt_Weight);
fChain->SetBranchAddress("Rwgt_size", &Rwgt_size, &b_Rwgt_size);
fChain->SetBranchAddress("Particle", &Particle_, &b_Particle_);
fChain->SetBranchAddress("Particle.fUniqueID", Particle_fUniqueID, &b_Particle_fUniqueID);
fChain->SetBranchAddress("Particle.fBits", Particle_fBits, &b_Particle_fBits);
fChain->SetBranchAddress("Particle.PID", Particle_PID, &b_Particle_PID);
fChain->SetBranchAddress("Particle.Status", Particle_Status, &b_Particle_Status);
fChain->SetBranchAddress("Particle.Mother1", Particle_Mother1, &b_Particle_Mother1);
fChain->SetBranchAddress("Particle.Mother2", Particle_Mother2, &b_Particle_Mother2);
fChain->SetBranchAddress("Particle.ColorLine1", Particle_ColorLine1, &b_Particle_ColorLine1);
fChain->SetBranchAddress("Particle.ColorLine2", Particle_ColorLine2, &b_Particle_ColorLine2);
fChain->SetBranchAddress("Particle.Px", Particle_Px, &b_Particle_Px);
fChain->SetBranchAddress("Particle.Py", Particle_Py, &b_Particle_Py);
fChain->SetBranchAddress("Particle.Pz", Particle_Pz, &b_Particle_Pz);
fChain->SetBranchAddress("Particle.E", Particle_E, &b_Particle_E);
fChain->SetBranchAddress("Particle.M", Particle_M, &b_Particle_M);
fChain->SetBranchAddress("Particle.PT", Particle_PT, &b_Particle_PT);
fChain->SetBranchAddress("Particle.Eta", Particle_Eta, &b_Particle_Eta);
fChain->SetBranchAddress("Particle.Phi", Particle_Phi, &b_Particle_Phi);
fChain->SetBranchAddress("Particle.Rapidity", Particle_Rapidity, &b_Particle_Rapidity);
fChain->SetBranchAddress("Particle.LifeTime", Particle_LifeTime, &b_Particle_LifeTime);
fChain->SetBranchAddress("Particle.Spin", Particle_Spin, &b_Particle_Spin);
fChain->SetBranchAddress("Particle_size", &Particle_size, &b_Particle_size);
}
Bool_t ang_distr::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
#endif // #ifdef ang_distr_cxxmy_func.h
#ifndef my_func_H
#define my_func_H
#include <cmath>
#include <iostream>
#include <algorithm>
#include <vector>
bool is_near_Z_res(TLorentzVector LV1, TLorentzVector LV2)
{
TLorentzVector LV_tot = LV1 + LV2;
if (LV_tot.M() > 81.0 && LV_tot.M() < 101.0) {return true;}
//if (LV_tot.M() > 84.406 && LV_tot.M() < 99.37) {return true;}
else {return false;}
}
#endifangles.h
#ifndef angles_H
#define angles_H
#include "TLorentzVector.h"
#include <iostream>
// Calculate cos(theta*) (angle between fermion and gauge boson boost)
inline Double_t calc_cos_star(TLorentzVector LV1,TLorentzVector LV2)
{
Double_t cos_star;
TLorentzVector cmstmp = LV1 + LV2;
LV1.Boost(-cmstmp.BoostVector());
cos_star=cos(LV1.Angle(cmstmp.Vect()));
return cos_star;
}
// Calculate cos(theta_v) (angle between Out1 and Intermediate1 in restframe)
inline Double_t calc_cos_v(TLorentzVector Out1, TLorentzVector Intermediate1, TLorentzVector cms)
{
Double_t cos_v;
Out1.Boost(-cms.BoostVector());
Intermediate1.Boost(-cms.BoostVector());
//cms.Boost(-cms.BoostVector());
cos_v=cos(Out1.Angle(Intermediate1.Vect()));
return cos_v;
}
#endif