Mela Class Reference

#include <Mela.h>

Collaboration diagram for Mela:

Public Member Functions
	Mela (double LHCsqrts_=13., double mh_=125., TVar::VerbosityLevel verbosity_=TVar::ERROR)
	the MELA constructor More...
	Mela (const Mela &other)
	Copy constructor for MELA. More...
	~Mela ()
	MELA destructor. More...
void	build (double mh_)
	This is the actual building of the tool that occurs in each instance of the Mela::Mela constructor. More...
void	setProcess (TVar::Process myModel, TVar::MatrixElement myME, TVar::Production myProduction)
	Sets the process, matrix element, and production that MELA is to use for this event. Calls ZZMatrixElement::set_Process, which calls TEvtProb::SetProcess. More...
void	setVerbosity (TVar::VerbosityLevel verbosity_=TVar::ERROR)
	Sets the verbosity for MELA outside of the initial constructor. More...
TVar::VerbosityLevel	getVerbosity ()
	Gets the current verbosity level for MELA. More...
void	setMelaLeptonInterference (TVar::LeptonInterference myLepInterf=TVar::DefaultLeptonInterf)
	Sets the MELA Lepton Interference. More...
void	setRemoveLeptonMasses (bool MasslessLeptonSwitch=true)
	either permits or forbids massive leptons. More...
void	setRemoveJetMasses (bool MasslessLeptonSwitch=true)
	either permits or forbids massive jets. More...
void	setMelaPrimaryHiggsMass (double myHiggsMass)
	Sets the mass of the "primary" higgs. More...
void	setMelaHiggsMass (double myHiggsMass, int index=0)
	Sets the mass of your chosen Higgs. More...
void	setMelaHiggsWidth (double myHiggsWidth=-1, int index=0)
	Sets the width of your chosen Higgs. More...
void	setMelaHiggsMassWidth (double myHiggsMass, double myHiggsWidth, int index)
	a combination of setMelaHiggsMass and setMelaHiggsWidth. More...
void	setRenFacScaleMode (TVar::EventScaleScheme renormalizationSch, TVar::EventScaleScheme factorizationSch, double ren_sf, double fac_sf)
	Sets the renormalization and the factorization schemes. More...
void	setCandidateDecayMode (TVar::CandidateDecayMode mode)
	Sets the decay mode for your event. More...
void	setCurrentCandidateFromIndex (unsigned int icand)
	Switches the candidate that you are working on to another candidate based off of an index. More...
void	setCurrentCandidate (MELACandidate *cand)
	Switches the candidate that you are working on to another candidate object specified. More...
void	setInputEvent (SimpleParticleCollection_t *pDaughters, SimpleParticleCollection_t *pAssociated=0, SimpleParticleCollection_t *pMothers=0, bool isGen=false)
	Sets the input event for MELA. MELA cannot run without this. More...
void	resetInputEvent ()
	Resets the event in preparation for the next iteration of the event loop. More...
void	setTempCandidate (SimpleParticleCollection_t *pDaughters, SimpleParticleCollection_t *pAssociated=0, SimpleParticleCollection_t *pMothers=0, bool isGen=false)
	Sets a temporary MELA candidate, by setting melaCand in Xcal2 to a temporary candidate without pushing this candidate to the candList of Xcal2. More...
void	appendTopCandidate (SimpleParticleCollection_t *TopDaughters)
	Adds a top quark as a MELA candidate. More...
void	resetMass (double inmass, int ipart)
	Resets the mass for a particle that is an electroweak parameter according to its id. More...
void	resetWidth (double inwidth, int ipart)
	Resets the width for a particle that is an electroweak parameter according to its id. More...
void	resetQuarkMasses ()
	Resets the masses of each quark to their original values. More...
void	resetMCFM_EWKParameters (double ext_Gf, double ext_aemmz, double ext_mW, double ext_mZ, double ext_xW, int ext_ewscheme=3)
	Resets the electroweak parameters back to their defaults. More...
double	getPrimaryMass (int ipart)
	A function to get the current primary EW/QCD parameters from MELA. More...
double	getPrimaryWidth (int ipart)
	A function to get the current primary EW/QCD parameters from MELA. More...
double	getHiggsWidthAtPoleMass (double mass)
	Returns the width of the Higgs at a given pole mass as a calculation. More...
MelaIO *	getIORecord ()
	Returns the MELAIO object, and by consequence, the entire parton-by-parton matrix element record. More...
MELACandidate *	getCurrentCandidate ()
	Gets the current MELA top-level (input) candList object. More...
int	getCurrentCandidateIndex ()
	Returns the index of the current MELA candidate - returns -1 if there is no candidate to be found. More...
int	getNCandidates ()
	Returns the size of the candidate list TEvtProb::candList. More...
std::vector< MELATopCandidate_t * > *	getTopCandidateCollection ()
	Same as getNCandidates but specifically for Top Quark Candidates. More...
void	getConstant (float &prob)
	This function returns a multiplicative constant to the matrix element calculation. More...
void	getPAux (float &prob)
RooSpin::modelMeasurables	getMeasurablesRRV ()
	Returns a RooSpin::modelMeasureables object containing many observable quantities. More...
void	computeDecayAngles (float &qH, float &m1, float &m2, float &costheta1, float &costheta2, float &Phi, float &costhetastar, float &Phi1)
	computes the decay angles for gg -> H -> ZZ as defined in Figure 1 of arXiv:1001.3396 More...
void	computeVBFAngles (float &Q2V1, float &Q2V2, float &costheta1, float &costheta2, float &Phi, float &costhetastar, float &Phi1)
	computes the decay angles for Vector Boson Fusion (VBF) -> H -> ZZ as defined in Figure 1 of arXiv:1208.4018 More...
void	computeVBFAngles_ComplexBoost (float &Q2V1, float &Q2V2, float &costheta1_real, float &costheta1_imag, float &costheta2_real, float &costheta2_imag, float &Phi, float &costhetastar, float &Phi1)
void	computeVHAngles (float &mVstar, float &mV, float &costheta1, float &costheta2, float &Phi, float &costhetastar, float &Phi1)
	computes the decay angles for Vector Boson Fusion (VBF) -> H -> ZZ as defined in Figure 1 of arXiv:1208.4018 More...
void	computeTTHAngles (int topDecay, float &mT1, float &mW1, float &mT2, float &mW2, float &costheta1, float &costheta2, float &Phi, float &costhetastar, float &Phi1, float &hbb, float &hWW, float &Phibb, float &Phi1bb, float &hWplusf, float &PhiWplusf, float &hWminusf, float &PhiWminusf)
void	computeP_selfDspin0 (double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2], float &prob, bool useConstant=true)
	This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 0. More...
void	computeP_selfDspin1 (double selfDZqqcoupl_input[SIZE_ZQQ][2], double selfDZvvcoupl_input[SIZE_ZVV][2], float &prob, bool useConstant=true)
	This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1. More...
void	computeP_selfDspin1 (double selfDZvvcoupl_input[SIZE_ZVV][2], float &prob, bool useConstant=true)
	This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1, but sets default values for Zqq couplings. More...
void	computeP_selfDspin2 (double selfDGggcoupl_input[SIZE_GGG][2], double selfDGqqcoupl_input[SIZE_GQQ][2], double selfDGvvcoupl_input[SIZE_GVV][2], float &prob, bool useConstant=true)
	This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 2. More...
void	computeP_selfDspin2 (double selfDGggcoupl_input[SIZE_GGG][2], double selfDGvvcoupl_input[SIZE_GVV][2], float &prob, bool useConstant=true)
	This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1, but sets default values for Zqq couplings. More...
void	computeP (float &prob, bool useConstant=true)
	Computes the probability for the probabilities on the decay side of things using the constituent daughter 4 vectors along with any jets that could exist. More...
void	computeD_CP (TVar::MatrixElement myME, TVar::Process myType, float &prob)
	computes the value of D_CP More...
void	computeProdDecP (double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2], double selfDHwwcoupl_input[nSupportedHiggses][SIZE_HVV][2], double selfDaTQGCcoupl_input[SIZE_ATQGC][2], double selfDAZffcoupl_input[SIZE_AZff][2], float &prob, bool useConstant=true)
	computes the combined production and decay probability while taking in coupling arrays More...
void	computeProdDecP (float &prob, bool useConstant=true)
	computes the combined production and decay probability, with couplings set beforehand More...
void	computeProdP (double selfDHggcoupl_input[SIZE_HGG][2], double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2], double selfDHwwcoupl_input[nSupportedHiggses][SIZE_HVV][2], float &prob, bool useConstant=true)
	computes Production side probabilities while taking in coupling arrays More...
void	computeProdP (float &prob, bool useConstant=true)
	computes Production side probabilities using couplings set beforehand More...
void	computeProdP_VH (double selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2], float &prob, bool includeHiggsDecay=false, bool useConstant=true)
void	computeProdP_VH (float &prob, bool includeHiggsDecay=false, bool useConstant=true)
void	computeProdP_ttH (float &prob, int topProcess=2, int topDecay=0, bool useConstant=true)
void	compute4FermionWeight (float &w)
void	getXPropagator (TVar::ResonancePropagatorScheme scheme, float &prop)
void	computePM4l (TVar::SuperMelaSyst syst, float &prob)
void	computeDijetConvBW (float &prob, bool useTrueBW=false)
void	computeD_gg (TVar::MatrixElement myME, TVar::Process myType, float &prob)
std::vector< TLorentzVector >	calculate4Momentum (double Mx, double M1, double M2, double theta, double theta1, double theta2, double Phi1, double Phi)

Static Public Member Functions
static void	cleanLinkedFiles ()

Public Attributes
TRandom3	melaRandomNumber
RooRealVar *	mzz_rrv
RooRealVar *	z1mass_rrv
RooRealVar *	z2mass_rrv
RooRealVar *	costhetastar_rrv
RooRealVar *	costheta1_rrv
RooRealVar *	costheta2_rrv
RooRealVar *	phi_rrv
RooRealVar *	phi1_rrv
RooRealVar *	Y_rrv
RooRealVar *	upFrac_rrv
RooAbsPdf *	pdf
ScalarPdfFactory_HVV *	ggSpin0Model
VectorPdfFactory *	spin1Model
TensorPdfFactory_ppHVV *	spin2Model
RooqqZZ_JHU_ZgammaZZ_fast *	qqZZmodel
SuperMELA *	super
double	selfDHggcoupl [nSupportedHiggses][SIZE_HGG][2]
double	selfDHg4g4coupl [nSupportedHiggses][SIZE_HGG][2]
double	selfDHqqcoupl [nSupportedHiggses][SIZE_HQQ][2]
double	selfDHbbcoupl [nSupportedHiggses][SIZE_HQQ][2]
double	selfDHttcoupl [nSupportedHiggses][SIZE_HQQ][2]
double	selfDHb4b4coupl [nSupportedHiggses][SIZE_HQQ][2]
double	selfDHt4t4coupl [nSupportedHiggses][SIZE_HQQ][2]
double	selfDHzzcoupl [nSupportedHiggses][SIZE_HVV][2]
double	selfDHwwcoupl [nSupportedHiggses][SIZE_HVV][2]
double	selfDHzzLambda_qsq [nSupportedHiggses][SIZE_HVV_LAMBDAQSQ][SIZE_HVV_CQSQ]
double	selfDHwwLambda_qsq [nSupportedHiggses][SIZE_HVV_LAMBDAQSQ][SIZE_HVV_CQSQ]
int	selfDHzzCLambda_qsq [nSupportedHiggses][SIZE_HVV_CQSQ]
int	selfDHwwCLambda_qsq [nSupportedHiggses][SIZE_HVV_CQSQ]
bool	differentiate_HWW_HZZ
double	selfDHzzpcoupl [SIZE_HVV][2]
double	selfDHzpzpcoupl [SIZE_HVV][2]
double	selfDZpffcoupl [SIZE_Vpff][2]
double	selfDHwwpcoupl [SIZE_HVV][2]
double	selfDHwpwpcoupl [SIZE_HVV][2]
double	selfDWpffcoupl [SIZE_Vpff][2]
double	selfDM_Zprime
double	selfDGa_Zprime
double	selfDM_Wprime
double	selfDGa_Wprime
double	selfDZqqcoupl [SIZE_ZQQ][2]
double	selfDZvvcoupl [SIZE_ZVV][2]
double	selfDGqqcoupl [SIZE_GQQ][2]
double	selfDGggcoupl [SIZE_GGG][2]
double	selfDGvvcoupl [SIZE_GVV][2]
double	selfDGvvpcoupl [SIZE_GVV][2]
double	selfDGvpvpcoupl [SIZE_GVV][2]
double	selfDaTQGCcoupl [SIZE_ATQGC][2]
double	selfDAZffcoupl [SIZE_AZff][2]

Protected Member Functions
void	printLogo () const
void	setSpinZeroCouplings ()
void	setSpinOneCouplings ()
void	setSpinTwoCouplings ()
void	setATQGCCouplings ()
void	setAZffCouplings ()
bool	configureAnalyticalPDFs ()
void	reset_SelfDCouplings ()
void	reset_PAux ()
void	reset_CandRef ()
void	constructDggr (float bkg_VAMCFM_noscale, float ggzz_VAMCFM_noscale, float ggHZZ_prob_pure_noscale, float ggHZZ_prob_int_noscale, float widthScale, float &myDggr)
void	getPConstantHandles ()
void	deletePConstantHandles ()
MelaPConstant *	getPConstantHandle (TVar::MatrixElement me_, TVar::Production prod_, TVar::Process proc_, TString relpath, TString spname, const bool useSqrts=false)
void	deletePConstantHandle (MelaPConstant *&handle)
void	computeConstant (float &prob)
void	setConstant ()
float	getConstant_JHUGenUndecayed ()
float	getConstant_4l ()
float	getConstant_2l2q ()
float	getConstant_4q ()
float	getConstant_FourFermionDecay (const int &decid)

Protected Attributes
double	LHCsqrts
TVar::Process	myModel_
TVar::MatrixElement	myME_
TVar::Production	myProduction_
TVar::LeptonInterference	myLepInterf_
TVar::VerbosityLevel	myVerbosity_
ZZMatrixElement *	ZZME
SuperDijetMela *	superDijet
float	auxiliaryProb
MELACandidate *	melaCand
MelaPConstant *	pAvgSmooth_JHUGen_JQCD_HSMHiggs [TVar::nFermionMassRemovalSchemes-1]
MelaPConstant *	pAvgSmooth_JHUGen_JJQCD_HSMHiggs [TVar::nFermionMassRemovalSchemes-1]
MelaPConstant *	pAvgSmooth_JHUGen_JJVBF_HSMHiggs [TVar::nFermionMassRemovalSchemes-1]
MelaPConstant *	pAvgSmooth_JHUGen_Had_ZH_HSMHiggs [TVar::nFermionMassRemovalSchemes-1]
MelaPConstant *	pAvgSmooth_JHUGen_Had_WH_HSMHiggs [TVar::nFermionMassRemovalSchemes-1]
MelaPConstant *	pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4mu
MelaPConstant *	pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4e
MelaPConstant *	pAvgSmooth_JHUGen_ZZGG_HSMHiggs_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_HSMHiggs_4mu
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_HSMHiggs_4e
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_S_HSMHiggs_4mu
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_S_HSMHiggs_4e
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_S_HSMHiggs_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_S_HSMHiggs_4mu
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_S_HSMHiggs_4e
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_S_HSMHiggs_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_S_HSMHiggs_4mu
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_S_HSMHiggs_4e
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_S_HSMHiggs_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_ZZQQB_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_ZZQQB_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_JJVBF_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_Had_ZH_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_Had_WH_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_JJQCD_bkgZZ_4mu
MelaPConstant *	pAvgSmooth_MCFM_JJQCD_bkgZZ_4e
MelaPConstant *	pAvgSmooth_MCFM_JJQCD_bkgZZ_2mu2e
MelaPConstant *	pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q

Detailed Description

Definition at line 48 of file Mela.h.

Constructor & Destructor Documentation

Mela() [1/2]

Mela::Mela	(	double	LHCsqrts_ = 13.,
		double	mh_ = 125.,
		TVar::VerbosityLevel	verbosity_ = TVar::ERROR
	)

the MELA constructor

Parameters

[in]	LHCsqrts_	The luminosity of your collider in TeV. Default is LHC at 13 TeV.
[in]	mh_	The mass of the Higgs in GeV. Default is 125 GeV
[in]	verbosity_	The verbosity of MELA that you desire, as defined in TVar::ERROR

Definition at line 94 of file Mela.cc.

    :
  melaRandomNumber(35797),
  LHCsqrts(LHCsqrts_),
  myVerbosity_(verbosity_),
  ZZME(0),
  auxiliaryProb(0.),
  melaCand(0)
{
  this->printLogo();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Start Mela constructor" << endl;
  build(mh_);
  if (myVerbosity_>=TVar::DEBUG) MELAout << "End Mela constructor" << endl;
}

Mela() [2/2]

Mela::Mela

(

const Mela &

other

)

Copy constructor for MELA.

Parameters

[in]

other

another MELA instance

Definition at line 111 of file Mela.cc.

                            :
melaRandomNumber(35797),
LHCsqrts(other.LHCsqrts),
myVerbosity_(other.myVerbosity_),
ZZME(0),
auxiliaryProb(0.),
melaCand(0)
{
  double mh_ = other.ZZME->get_PrimaryHiggsMass();
  build(mh_);
}

~Mela()

Mela::~Mela

(

)

MELA destructor.

Definition at line 122 of file Mela.cc.

           {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Begin Mela destructor" << endl;
 
  //setRemoveLeptonMasses(false); // Use Run 1 scheme for not removing lepton masses. Notice the switch itself is defined as an extern, so it has to be set to default value at the destructor!
  setRemoveLeptonMasses(true); // Use Run 2 scheme for removing lepton masses. Notice the switch itself is defined as an extern, so it has to be set to default value at the destructor!
 
  // Delete the derived RooFit objects first...
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying analytical PDFs" << endl;
  delete ggSpin0Model;
  delete spin1Model;
  delete spin2Model;
  delete qqZZmodel;
  // ...then delete the observables.
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying analytical PDFs observables" << endl;
  delete mzz_rrv;
  delete z1mass_rrv;
  delete z2mass_rrv;
  delete costhetastar_rrv;
  delete costheta1_rrv;
  delete costheta2_rrv;
  delete phi_rrv;
  delete phi1_rrv;
  delete Y_rrv;
  delete upFrac_rrv;
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying SuperDijetMELA" << endl;
  delete superDijet;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying SuperMELA" << endl;
  delete super;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying ZZME" << endl;
  delete ZZME;
 
  // Delete ME constant handles
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela destructor: Destroying PConstant handles" << endl;
  deletePConstantHandles();
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "End Mela destructor" << endl;
}

Member Function Documentation

appendTopCandidate()

void Mela::appendTopCandidate

(

SimpleParticleCollection_t *

TopDaughters

)

Adds a top quark as a MELA candidate.

See also

Wrapper for ZZMatrixElement::append_TopCandidate, which is a wrapper for TEvtProb::AppendTopCandidate

Parameters

[in]

TopDaughters

A SimpleParticleCollection_t of the top's daughters, or decay products

Definition at line 363 of file Mela.cc.

{ ZZME->append_TopCandidate(TopDaughters); }

build()

void Mela::build

(

double

mh_

)

This is the actual building of the tool that occurs in each instance of the Mela::Mela constructor.

Parameters

mh_	This is the mass of the Higgs in GeV

Definition at line 175 of file Mela.cc.

                          {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Start Mela::build" << endl;
  //setRemoveLeptonMasses(false); // Use Run 1 scheme for not removing fermion masses
  setRemoveLeptonMasses(true); // Use Run 2 scheme for removing fermion masses to compute MEs that expect massless fermions properly
 
  const double maxSqrts = 8.;
 
  // Create symlinks to the required files, if these are not already present (do nothing otherwise)
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Create symlinks to the required files if these are not already present:" << endl;
 
  const string MELAPKGPATH = TVar::GetMELAPath();
  if (myVerbosity_>=TVar::DEBUG)  MELAout << "\t- MELA package path: " << MELAPKGPATH << endl;
 
  const string mcfmWarning = MELAPKGPATH + "data/ffwarn.dat"; symlink(mcfmWarning.c_str(), "ffwarn.dat");
  const string mcfm_brsm_o = MELAPKGPATH + "data/br.sm1"; symlink(mcfm_brsm_o.c_str(), "br.sm1");
  const string mcfm_brsm_t = MELAPKGPATH + "data/br.sm2"; symlink(mcfm_brsm_t.c_str(), "br.sm2");
  const string mcfmInput1 = MELAPKGPATH + "data/input.DAT"; symlink(mcfmInput1.c_str(), "input.DAT");
  const string mcfmInput2 = MELAPKGPATH + "data/process.DAT"; symlink(mcfmInput2.c_str(), "process.DAT");
  if (myVerbosity_>=TVar::DEBUG) MELAout << "\t- MCFM symlinks are done" << endl;
  mkdir("Pdfdata", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  const string mcfmInput3 = MELAPKGPATH + "data/Pdfdata/cteq6l1.tbl"; symlink(mcfmInput3.c_str(), "Pdfdata/cteq6l1.tbl");
  const string mcfmInput4 = MELAPKGPATH + "data/Pdfdata/cteq6l.tbl"; symlink(mcfmInput4.c_str(), "Pdfdata/cteq6l.tbl");
  if (myVerbosity_>=TVar::DEBUG) MELAout << "\t- PDF symlinks are done" << endl;
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Create variables used in anaMELA" << endl;
  mzz_rrv = new RooRealVar("mzz", "m_{ZZ}", mh_, 0., 1000.);
  z1mass_rrv = new RooRealVar("z1mass", "m_{Z1}", 0., 160.);
  z2mass_rrv = new RooRealVar("z2mass", "m_{Z2}", 0., 200.);
  costhetastar_rrv = new RooRealVar("costhetastar", "cos#theta^{*}", -1., 1.);
  costheta1_rrv = new RooRealVar("costheta1", "cos#theta_{1}", -1., 1.);
  costheta2_rrv = new RooRealVar("costheta2", "cos#theta_{2}", -1., 1.);
  phi_rrv= new RooRealVar("phi", "#Phi", -TMath::Pi(), TMath::Pi());
  phi1_rrv= new RooRealVar("phi1", "#Phi_{1}", -TMath::Pi(), TMath::Pi());
  Y_rrv = new RooRealVar("Yzz", "#Y_{ZZ}", 0, -4, 4);
  upFrac_rrv = new RooRealVar("upFrac", "fraction up-quarks", .5, 0., 1.);
 
  RooSpin::modelMeasurables measurables_;
  measurables_.h1 = costheta1_rrv;
  measurables_.h2 = costheta2_rrv;
  measurables_.Phi = phi_rrv;
  measurables_.m1 = z1mass_rrv;
  measurables_.m2 = z2mass_rrv;
  measurables_.m12 = mzz_rrv;
  measurables_.hs = costhetastar_rrv;
  measurables_.Phi1 = phi1_rrv;
  measurables_.Y = Y_rrv;
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Create anaMELA PDF factories" << endl;
  ggSpin0Model = new ScalarPdfFactory_HVV(measurables_, false, RooSpin::kVdecayType_Zll, RooSpin::kVdecayType_Zll); // RooSpin::kVdecayType_Zll,RooSpin::kVdecayType_Zll==ZZ4l
  spin1Model = new VectorPdfFactory(z1mass_rrv, z2mass_rrv, costhetastar_rrv, costheta1_rrv, costheta2_rrv, phi_rrv, phi1_rrv, mzz_rrv);
  spin2Model = new TensorPdfFactory_ppHVV(measurables_, RooSpin::kVdecayType_Zll, RooSpin::kVdecayType_Zll);
  qqZZmodel = new RooqqZZ_JHU_ZgammaZZ_fast("qqZZmodel", "qqZZmodel", *z1mass_rrv, *z2mass_rrv, *costheta1_rrv, *costheta2_rrv, *phi_rrv, *costhetastar_rrv, *phi1_rrv, *mzz_rrv, *upFrac_rrv);
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Paths for ZZMatrixElement" << endl;
  const string path_HiggsWidthFile = MELAPKGPATH + "data/HiggsTotalWidth_YR3.txt";
  if (myVerbosity_>=TVar::DEBUG) MELAout << "\t- Cross section/width file: " << path_HiggsWidthFile << endl;
  const string path_nnpdf = MELAPKGPATH + "data/Pdfdata/NNPDF30_lo_as_0130.LHgrid";
  char path_nnpdf_c[] = "Pdfdata/NNPDF30_lo_as_0130.LHgrid";
  int pdfmember = 0;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "\t- Linking NNPDF path " << path_nnpdf.c_str() << " -> " << path_nnpdf_c << endl;
  symlink(path_nnpdf.c_str(), path_nnpdf_c);
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Start ZZMatrixElement" << endl;
  ZZME = new ZZMatrixElement(path_nnpdf_c, pdfmember, path_HiggsWidthFile.substr(0, path_HiggsWidthFile.length()-23).c_str(), 1000.*LHCsqrts/2., myVerbosity_);
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Set ZZMatrixElement masses" << endl;
  setMelaPrimaryHiggsMass(mh_);
  setMelaHiggsMass(mh_, 0); setMelaHiggsMass(-1., 1);
  setMelaHiggsWidth(-1., 0); setMelaHiggsWidth(-1., 1);
  setMelaLeptonInterference(TVar::DefaultLeptonInterf);
  setCandidateDecayMode(TVar::CandidateDecay_ZZ); // Default decay mode is ZZ at the start
 
 
  /***** CONSTANTS FOR MATRIX ELEMENTS *****/
  getPConstantHandles();
 
 
  /***** SuperMELA *****/
  // Deactivate generation messages
  RooMsgService::instance().getStream(1).removeTopic(NumIntegration);
  RooMsgService::instance().setStreamStatus(1, kFALSE);
  RooMsgService::instance().setStreamStatus(0, kFALSE);// silence also the error messages, but should really be looked at.
 
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Start superMELA" << endl;
  int superMELA_LHCsqrts = LHCsqrts;
  if (superMELA_LHCsqrts > maxSqrts) superMELA_LHCsqrts = maxSqrts;
  super = new SuperMELA(mh_, "4mu", superMELA_LHCsqrts); // preliminary intialization, we adjust the flavor later
  char cardpath[500];
  sprintf(cardpath, "data/CombinationInputs/SM_inputs_%dTeV/inputs_4mu.txt", superMELA_LHCsqrts);
  string cardfile = MELAPKGPATH + cardpath;
  super->SetPathToCards(cardfile.substr(0, cardfile.length()-14).c_str());
  super->SetVerbosity((myVerbosity_>=TVar::DEBUG));
  super->init();
 
 
  /***** SuperDijetMela *****/
  float superDijetSqrts = LHCsqrts;
  if (superDijetSqrts<13.) superDijetSqrts=13.; // Dijet resolution does not exist for 7 or 8 TeV
  superDijet = new SuperDijetMela(superDijetSqrts, myVerbosity_);
 
  // Initialize the couplings to 0 and end Mela constructor
  reset_SelfDCouplings();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "End Mela::build" << endl;
}

calculate4Momentum()

std::vector< TLorentzVector > Mela::calculate4Momentum	(	double	Mx,
		double	M1,
		double	M2,
		double	theta,
		double	theta1,
		double	theta2,
		double	Phi1,
		double	Phi
	)

Definition at line 524 of file Mela.cc.

                                                                                                                                                      {
    return ZZME->Calculate4Momentum(Mx, M1, M2, theta, theta1, theta2, Phi1, Phi);
}

compute4FermionWeight()

void Mela::compute4FermionWeight

(

float &

w

)

Definition at line 1807 of file Mela.cc.

                                        { // Lepton interference using JHUGen
  reset_PAux();
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    bool hasFailed=false;
    int id_original[2][2];
    for (int iv=0; iv<2; iv++){
      MELAParticle* Vi = melaCand->getSortedV(iv);
      int ndau=Vi->getNDaughters();
      if (ndau!=2 || !(PDGHelpers::isAZBoson(Vi->id) || PDGHelpers::isAPhoton(Vi->id))){ w=1; hasFailed=true; break; } // Veto WW, ZG, GG
      for (int ivd=0; ivd<2; ivd++) id_original[iv][ivd]=Vi->getDaughter(ivd)->id;
    }
    if (
      !PDGHelpers::isALepton(id_original[0][0])
      ||
      !PDGHelpers::isALepton(id_original[0][1])
      ||
      !PDGHelpers::isALepton(id_original[1][0])
      ||
      !PDGHelpers::isALepton(id_original[1][1])
      ){
      if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeWeight: Function is not implemented for decay states other than 4l/2l2l." << endl;
      w=0;
      hasFailed=true;
    }
    /*
    if (
    (id_original[0][0]==0 && id_original[0][1]==0)
    ||
    (id_original[1][0]==0 && id_original[1][1]==0)
    ){ // Return 1 if any pairs of quarks are unknown
    w=1;
    hasFailed=true;
    }
    */
 
    if (!hasFailed){
      float dXsec_HZZ_JHU, dXsec_HZZ_JHU_interf; // temporary prob
 
      // Calculate dXsec ratio by directly modifying the candidate
      computeP(dXsec_HZZ_JHU, false);
      for (int ivd=0; ivd<2; ivd++) melaCand->getSortedV(1)->getDaughter(ivd)->id=id_original[0][0]*(1-2*ivd);
      computeP(dXsec_HZZ_JHU_interf, false);
      for (int ivd=0; ivd<2; ivd++) melaCand->getSortedV(1)->getDaughter(ivd)->id=id_original[1][ivd];
 
      w = dXsec_HZZ_JHU_interf / dXsec_HZZ_JHU;
      // protect against anomalously large weights
      if (w>5.) w=25./w;
    }
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
}

computeD_CP()

void Mela::computeD_CP	(	TVar::MatrixElement	myME,
		TVar::Process	myType,
		float &	prob
	)

computes the value of D_CP

See also

calls Mela::computeP_selfDspin0

Parameters

[in]	myME	the matrix element you want to use, set as a TVar::MatrixElement
[in]	myType	the process you want to use, set as a TVar::Process
[out]	prob	This is the value of the output probability - edited by reference

Definition at line 1339 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeD_CP" << endl;
  double coupl_mix[nSupportedHiggses][SIZE_HVV][2] ={ { { 0 } } };
  double coupl_1[nSupportedHiggses][SIZE_HVV][2] ={ { { 0 } } };
  double coupl_2[nSupportedHiggses][SIZE_HVV][2] ={ { { 0 } } };
 
  switch (myType){
  case TVar::D_g1g4:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][3][0] =2.521;
    coupl_1[0][0][0] =1.;
    coupl_2[0][3][0] =2.521;
    break;
  case TVar::D_g1g4_pi_2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][3][1] =2.521;
    coupl_1[0][0][0] =1.;
    coupl_2[0][3][1] =2.521;
    break;
  case TVar::D_g1g2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][1][0] = 1.638;
    coupl_1[0][0][0] =1.;
    coupl_2[0][1][0] = 1.638;
    break;
  case TVar::D_g1g2_pi_2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][1][1] = 1.638;
    coupl_1[0][0][0] =1.;
    coupl_2[0][1][1] = 1.638;
    break;
  case TVar::D_g1g1prime2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][11][0] = 12046.01;
    coupl_1[0][0][0] =1.;
    coupl_2[0][11][0] = 12046.01;
    break;
  case TVar::D_zzzg:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][4][0] = 0.0688;
    coupl_1[0][0][0] =1.;
    coupl_2[0][4][0] = 0.0688;
    break;
  case TVar::D_zzgg:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][7][0] = -0.0898;
    coupl_1[0][0][0] =1.;
    coupl_2[0][7][0] = -0.0898;
    break;
  case TVar::D_zzzg_PS:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][6][0] = 0.0855;
    coupl_1[0][0][0] =1.;
    coupl_2[0][6][0] = 0.0855;
    break;
  case TVar::D_zzgg_PS:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][9][0] = -0.0907;
    coupl_1[0][0][0] =1.;
    coupl_2[0][9][0] = -0.0907;
    break;
  case TVar::D_zzzg_g1prime2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][30][0] = -7591.914;
    coupl_1[0][0][0] =1.;
    coupl_2[0][30][0] = -7591.914;
    break;
  case TVar::D_zzzg_g1prime2_pi_2:
    coupl_mix[0][0][0] =1.;
    coupl_mix[0][30][1] = -7591.914;
    coupl_1[0][0][0] =1.;
    coupl_2[0][30][1] = -7591.914;
    break;
  default:
    MELAout <<"Error: Not supported!"<<endl;
  }
 
  float pMix, p1, p2;
  setProcess(TVar::SelfDefine_spin0, myME, TVar::ZZGG);
  computeP_selfDspin0(coupl_mix, pMix, true);
  computeP_selfDspin0(coupl_1, p1, true);
  computeP_selfDspin0(coupl_2, p2, true);
  prob = pMix- p1- p2;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeD_CP" << endl;
}

computeD_gg()

void Mela::computeD_gg	(	TVar::MatrixElement	myME,
		TVar::Process	myType,
		float &	prob
	)

Definition at line 1934 of file Mela.cc.

   {
  prob=-99;
  if (myME != TVar::MCFM || myType != TVar::D_gg10){
    MELAout << "Only support MCFM and D_gg10"<<endl;
    return;
  }
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    float bkg_VAMCFM, ggzz_VAMCFM_noscale, ggHZZ_prob_pure_noscale, ggHZZ_prob_int_noscale, bkgHZZ_prob_noscale;
    float ggScale=0;
    setProcess(TVar::bkgZZ, myME, TVar::ZZGG); computeP(ggzz_VAMCFM_noscale, false);
    setProcess(TVar::HSMHiggs, myME, TVar::ZZGG); computeP(ggHZZ_prob_pure_noscale, false);
    setProcess(TVar::bkgZZ_SMHiggs, myME, TVar::ZZGG); computeP(bkgHZZ_prob_noscale, false); setConstant(); getConstant(ggScale);
    if (ggScale>0.){
      bkgHZZ_prob_noscale /= ggScale;
      ggHZZ_prob_pure_noscale /= ggScale;
      ggzz_VAMCFM_noscale /= ggScale;
    }
    ggHZZ_prob_int_noscale = bkgHZZ_prob_noscale - ggHZZ_prob_pure_noscale -  ggzz_VAMCFM_noscale;
 
    setProcess(TVar::bkgZZ, myME, TVar::ZZQQB); computeP(bkg_VAMCFM, true);
 
    constructDggr(bkg_VAMCFM, ggzz_VAMCFM_noscale, ggHZZ_prob_pure_noscale, ggHZZ_prob_int_noscale, 10., prob); // Use 10 for Dgg10
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
}

computeDecayAngles()

void Mela::computeDecayAngles	(	float &	qH,
		float &	m1,
		float &	m2,
		float &	costheta1,
		float &	costheta2,
		float &	Phi,
		float &	costhetastar,
		float &	Phi1
	)

computes the decay angles for gg -> H -> ZZ as defined in Figure 1 of arXiv:1001.3396

See also

Calls TUtil::computeAngles to fully calculate angles

Use TUtil::computeAngles if you would like to calculate the angles through inputting your own 4-vectors and quantities rather than through the MELA event loop

Attention

This function requires there to be an input event already defined

Warning

This function will edit all values inplace. Every value is technically an output.

The selection for m1 and m2 are based upon which of the leptons is closer to the z mass. $m_1 \leq m_2$.

Parameters

[out]	qH	The mass of the "Higgs" as reconstructed by the constituent 4 leptons
[out]	m1	The mass of the first decay particle as reconstructed by 2 of the leptons
[out]	m2	The mass of the first decay particle as reconstructed by 2 of the leptons
[out]	costheta1	$\cos(\theta_1)$
[out]	costheta2	$\cos(\theta_2)$
[out]	Phi	$\phi$
[out]	costhetastar	$\cos(\theta*)$
[out]	Phi1	$\phi_1$

Definition at line 558 of file Mela.cc.

 {
  using TVar::simple_event_record;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeDecayAngles" << endl;
 
  qH=0; m1=0; m2=0; costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
 
    int partIncCode=TVar::kNoAssociated; // Only use associated partons in the pT=0 frame boost
    simple_event_record mela_event;
    mela_event.AssociationCode=partIncCode;
    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
    if (daughters.size()<2 || daughters.size()>4 || mela_event.intermediateVid.size()!=2){
      if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeDecayAngles: Number of daughters " << daughters.size() << " or number of intermediate Vs " << mela_event.intermediateVid << " not supported!" << endl;
      return;
    }
 
    // Make sure there are exactly 4 daughters, null or not
    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
    else if (daughters.size()==2){
      daughters.push_back(SimpleParticle_t(-9000, nullVector));
      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
    }
    qH = (daughters.at(0).second+daughters.at(1).second+daughters.at(2).second+daughters.at(3).second).M();
    m1 = (daughters.at(0).second+daughters.at(1).second).M();
    m2 = (daughters.at(2).second+daughters.at(3).second).M();
 
    TUtil::computeAngles(
      costhetastar, costheta1, costheta2, Phi, Phi1,
      daughters.at(0).second, daughters.at(0).first,
      daughters.at(1).second, daughters.at(1).first,
      daughters.at(2).second, daughters.at(2).first,
      daughters.at(3).second, daughters.at(3).first
    );
 
    // Protect against NaN
    if (!std::isfinite(costhetastar)) costhetastar=0;
    if (!std::isfinite(costheta1)) costheta1=0;
    if (!std::isfinite(costheta2)) costheta2=0;
    if (!std::isfinite(Phi)) Phi=0;
    if (!std::isfinite(Phi1)) Phi1=0;
  }
  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeDecayAngles: No possible melaCand in TEvtProb to compute angles." << endl;
 
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeDecayAngles" << endl;
}

computeDijetConvBW()

void Mela::computeDijetConvBW	(	float &	prob,
		bool	useTrueBW = false
	)

Definition at line 2792 of file Mela.cc.

                                                        {
  melaCand = getCurrentCandidate();
  prob = superDijet->GetConvBW(myProduction_, melaCand, useTrueBW);
  reset_CandRef();
}

computeP()

void Mela::computeP	(	float &	prob,
		bool	useConstant = true
	)

Computes the probability for the probabilities on the decay side of things using the constituent daughter 4 vectors along with any jets that could exist.

See also

Mela::computeDecayAngles

Attention

All matrix elements (JHUGen, MCFM, or analytic), productions, and processes can be used here

Parameters

[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1183 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeP" << endl;
  reset_PAux();
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
    float mZZ=0, mZ1=0, mZ2=0, costheta1=0, costheta2=0, Phi=0, costhetastar=0, Phi1=0;
 
    if (myME_ == TVar::ANALYTICAL){
      computeDecayAngles(
        mZZ, mZ1, mZ2,
        costheta1, costheta2, Phi,
        costhetastar, Phi1
        );
      costhetastar_rrv->setVal(costhetastar);
      costheta1_rrv->setVal(costheta1);
      costheta2_rrv->setVal(costheta2);
      phi_rrv->setVal(Phi);
      phi1_rrv->setVal(Phi1);
      z1mass_rrv->setVal(mZ1);
      z2mass_rrv->setVal(mZ2);
      mzz_rrv->setVal(mZZ);
      Y_rrv->setConstant(true); // Just to avoid integrating over this variable unnecessarily
 
      bool computeAnaMELA = configureAnalyticalPDFs();
      if (computeAnaMELA){
        if (myProduction_==TVar::ZZINDEPENDENT){
          RooAbsPdf* integral = (RooAbsPdf*)pdf->createIntegral(RooArgSet(*costhetastar_rrv, *phi1_rrv));
          prob = integral->getVal();
          delete integral;
        }
        else prob = pdf->getVal();
      }
      else if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeP: The specified anaMELA configuration is not valid!" << endl;
 
      Y_rrv->setConstant(false);
    }
    else if (myME_ == TVar::JHUGen || myME_ == TVar::MCFM){
      setAZffCouplings();
      if (!(myME_ == TVar::MCFM  && myProduction_ == TVar::ZZINDEPENDENT &&  (myModel_ == TVar::bkgZZ || myModel_ == TVar::bkgWW || myModel_ == TVar::bkgZGamma || myModel_ == TVar::bkgGammaGamma))){
        if (myME_ == TVar::MCFM || myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
        else if (myModel_ == TVar::SelfDefine_spin1) setSpinOneCouplings();
        else if (myModel_ == TVar::SelfDefine_spin2) setSpinTwoCouplings();
        ZZME->computeXS(prob);
      }
      else{
        computeDecayAngles(
          mZZ, mZ1, mZ2,
          costheta1, costheta2, Phi,
          costhetastar, Phi1
          );
 
        if (myVerbosity_>=TVar::DEBUG){ // Notify first
          MELAout << "Mela::computeP: Condition (myME_ == TVar::MCFM  && myProduction_ == TVar::ZZINDEPENDENT &&  myModel_ == TVar::bkgZZ/WW/ZGamma/ZJJ/GammaGamma)." << endl;
          vector<TLorentzVector> pDauVec = calculate4Momentum(mZZ, mZ1, mZ1, acos(costhetastar), acos(costheta1), acos(costheta2), Phi1, Phi);
          MELAout
            << "\tOriginal mZZ=" << mZZ << " "
            << "m1=" << mZ1 << " "
            << "m2=" << mZ2 << " "
            << "h1=" << costheta1 << " "
            << "h2=" << costheta2 << " "
            << "Phi=" << Phi << " "
            << "hs=" << costhetastar << " "
            << "Phi1=" << Phi1 << endl;
          MELAout << "\tfor daughters:" << endl;
          for (int iv=0; iv<2; iv++){
            for (int idau=0; idau<min(2, melaCand->getSortedV(iv)->getNDaughters()); idau++){
              MELAout
                << "id=" << melaCand->getSortedV(iv)->getDaughter(idau)->id << " "
                << "x=" << pDauVec.at(2*iv+idau).X() << " "
                << "y=" << pDauVec.at(2*iv+idau).Y() << " "
                << "z=" << pDauVec.at(2*iv+idau).Z() << " "
                << "t=" << pDauVec.at(2*iv+idau).T() << endl;
            }
          }
 
        }
 
        prob = 0;
        int gridsize_hs = 5;
        double hs_min = 0; //-1.;
        double hs_max = 1;
        double hs_step = (hs_max - hs_min) / double(gridsize_hs);
 
        int gridsize_phi1 = 5;
        double phi1_min = 0; //-TMath::Pi();
        double phi1_max = TMath::Pi();
        double phi1_step = (phi1_max - phi1_min) / double(gridsize_phi1);
 
        for (int i_hs = 0; i_hs < gridsize_hs + 1; i_hs++){
          double hs_val = hs_min + i_hs * hs_step;
          for (int i_phi1 = 0; i_phi1 < gridsize_phi1 + 1; i_phi1++){
            double phi1_val = phi1_min + i_phi1 * phi1_step;
            float temp_prob=0;
 
            // Get identical 4-vectors
            SimpleParticleCollection_t daughters;
            vector<TLorentzVector> pDauVec = calculate4Momentum(mZZ, mZ1, mZ2, acos(hs_val), acos(costheta1), acos(costheta2), phi1_val, Phi);
            for (int iv=0; iv<2; iv++){
              for (int idau=0; idau<min(2, melaCand->getSortedV(iv)->getNDaughters()); idau++){
                SimpleParticle_t tmpPair(melaCand->getSortedV(iv)->getDaughter(idau)->id, pDauVec.at(2*iv+idau));
                daughters.push_back(tmpPair);
              }
            }
            if (myVerbosity_>=TVar::DEBUG){ // Summarize the integrated particles
              MELAout << "Mela::computeP: hs, Phi1 are now " << hs_val << " " << phi1_val << endl;
              unsigned int idau=1;
              for (SimpleParticle_t const& tmpPart:daughters){
                MELAout << "Dau " << idau << " "
                  << "id=" << tmpPart.first << " "
                  << "x=" << tmpPart.second.X() << " "
                  << "y=" << tmpPart.second.Y() << " "
                  << "z=" << tmpPart.second.Z() << " "
                  << "t=" << tmpPart.second.T() << endl;
                idau++;
              }
            }
            vector<MELAParticle*> partList_tmp;
            vector<MELACandidate*> candList_tmp;
            MELACandidate* cand_tmp = TUtil::ConvertVectorFormat(
              &daughters,
              0,
              0,
              false,
              &partList_tmp,
              &candList_tmp
              );
            if (myVerbosity_>=TVar::ERROR && !cand_tmp) MELAerr << "Mela::computeP: Failed to construct temporary candidate!" << endl;
            setCurrentCandidate(cand_tmp);
            if (myVerbosity_>=TVar::DEBUG && cand_tmp){ MELAout << "Mela::computeP: ZZINDEPENDENT calculation produces candidate:" << endl; TUtil::PrintCandidateSummary(cand_tmp); }
            // calculate the ME
            ZZME->computeXS(temp_prob);
            // Delete the temporary particles
            for (MELACandidate* tmpPart:candList_tmp) delete tmpPart; // Only one candidate should really be here
            for (MELAParticle* tmpPart:partList_tmp) delete tmpPart;
            setCurrentCandidate(melaCand);
            prob += temp_prob;
          }
        }
        prob =  prob / float((gridsize_hs + 1) * (gridsize_phi1 +1));
      }
    }
 
    if (useConstant) computeConstant(prob);
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeP" << endl;
}

computeP_selfDspin0()

void Mela::computeP_selfDspin0	(	double	selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		float &	prob,
		bool	useConstant = true
	)

This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 0.

See also

Mela::computeP, Mela::selfDHqqcoupl, Mela::selfDHggcoupl, Mela::selfDHttcoupl, Mela::selfDHbbcoupl, Mela::selfDHzzcoupl, Mela::selfDHwwcoupl

Warning

I (Mohit Srivastav) would personally avoid using this function, as it obfuscates the couplings you are really using .

It is always a good idea to set couplings yourself so that you know what you're really doing.

Parameters

[in]	selfDHvvcoupl_input	an array corresponding to the couplings of Mela::selfDHzzcoupl (between the Higgs and the Z Boson), Mela::selfDHwwcoupl using the indices in TCouplingsBase.hh
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1077 of file Mela.cc.

   {
  //remove this line, it's here for counting purposes.  Zpffcoupl
  // Don't set these, and you will get 0.
  if (myME_==TVar::JHUGen){
    for (int jh=0; jh<(int)nSupportedHiggses; jh++){
      // JHUGen ME production side is Hqq=1, Hgg=1
      selfDHqqcoupl[jh][0][0] = 1.0;
      selfDHggcoupl[jh][0][0] = 1.0;
    }
  }
  else if (myME_==TVar::MCFM){
    for (int jh=0; jh<(int)nSupportedHiggses; jh++){
      // MCFM ME production side is Htt=Hbb=1
      selfDHttcoupl[jh][0][0] = 1.0;
      selfDHbbcoupl[jh][0][0] = 1.0;
    }
  }
  for (int jh=0; jh<(int)nSupportedHiggses; jh++){
    for (int im=0; im<2; im++){
      for (int ic=0; ic<SIZE_HVV; ic++){
        selfDHzzcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im];
        selfDHwwcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im]; // Just for extra protection since differentiate_HWW_HZZ is set to false.
      }
    }
  }
 
  computeP(
    prob,
    useConstant
    );
}

computeP_selfDspin1() [1/2]

void Mela::computeP_selfDspin1	(	double	selfDZqqcoupl_input[SIZE_ZQQ][2],
		double	selfDZvvcoupl_input[SIZE_ZVV][2],
		float &	prob,
		bool	useConstant = true
	)

This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1.

See also

Mela::computeP, Mela::selfDHqqcoupl, Mela::selfDHggcoupl, Mela::selfDHttcoupl, Mela::selfDHbbcoupl, Mela::selfDHzzcoupl, Mela::selfDHwwcoupl

Warning

I (Mohit Srivastav) would personally avoid using this function, as it obfuscates the couplings you are really using .

It is always a good idea to set couplings yourself so that you know what you're really doing.

Parameters

[in]	selfDZqqcoupl_input	an array corresponding to the couplings of Mela::selfDZqqcoupl (between the Higgs and quarks) using the indices in TCouplingsBase.hh
[in]	selfDZvvcoupl_input	an array corresponding to the couplings of Mela::selfDZvvcoupl (between the Higgs and the Z/W Bosons) using the indices in TCouplingsBase.hh
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1112 of file Mela.cc.

   {
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_ZQQ; ic++) selfDZqqcoupl[ic][im] = selfDZqqcoupl_input[ic][im];
    for (int ic=0; ic<SIZE_ZVV; ic++) selfDZvvcoupl[ic][im] = selfDZvvcoupl_input[ic][im];
  }
 
  computeP(
    prob,
    useConstant
    );
}

computeP_selfDspin1() [2/2]

void Mela::computeP_selfDspin1	(	double	selfDZvvcoupl_input[SIZE_ZVV][2],
		float &	prob,
		bool	useConstant = true
	)

This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1, but sets default values for Zqq couplings.

See also

Mela::computeP, Mela::selfDHqqcoupl, Mela::selfDHggcoupl, Mela::selfDHttcoupl, Mela::selfDHbbcoupl, Mela::selfDHzzcoupl, Mela::selfDHwwcoupl

Warning

I (Mohit Srivastav) would personally avoid using this function, as it obfuscates the couplings you are really using .

It is always a good idea to set couplings yourself so that you know what you're really doing.

Parameters

[in]	selfDZvvcoupl_input	an array corresponding to the couplings of Mela::selfDZvvcoupl (between the Higgs and the Z/W Bosons) using the indices in TCouplingsBase.hh
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1128 of file Mela.cc.

   {
  // Initialize the quark_left_right couplings
  selfDZqqcoupl[0][0]=1.0;
  selfDZqqcoupl[1][0]=1.0;
 
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_ZVV; ic++) selfDZvvcoupl[ic][im] = selfDZvvcoupl_input[ic][im];
  }
 
  computeP(
    prob,
    useConstant
    );
}

computeP_selfDspin2() [1/2]

void Mela::computeP_selfDspin2	(	double	selfDGggcoupl_input[SIZE_GGG][2],
		double	selfDGqqcoupl_input[SIZE_GQQ][2],
		double	selfDGvvcoupl_input[SIZE_GVV][2],
		float &	prob,
		bool	useConstant = true
	)

This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 2.

See also

Mela::computeP, Mela::selfDGggcoupl, Mela::selfDGqqcoupl, Mela::selfDGvvcoupl

Warning

I (Mohit Srivastav) would personally avoid using this function, as it obfuscates the couplings you are really using .

It is always a good idea to set couplings yourself so that you know what you're really doing.

Parameters

[in]	selfDGggcoupl_input	an array corresponding to the couplings of Mela::selfDGggcoupl (between the "graviton" and gluons) using the indices in TCouplingsBase.hh
[in]	selfDGqqcoupl_input	an array corresponding to the couplings of Mela::selfDGqqcoupl (between the "graviton" and quarks) using the indices in TCouplingsBase.hh
[in]	selfDGvvcoupl_input	an array corresponding to the couplings of Mela::selfDGvvcoupl (between the "graviton" and the W/Z bosons) using the indices in TCouplingsBase.hh
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1146 of file Mela.cc.

   {
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_GGG; ic++) selfDGggcoupl[ic][im] = selfDGggcoupl_input[ic][im];
    for (int ic=0; ic<SIZE_GQQ; ic++) selfDGqqcoupl[ic][im] = selfDGqqcoupl_input[ic][im];
    for (int ic=0; ic<SIZE_GVV; ic++) selfDGvvcoupl[ic][im] = selfDGvvcoupl_input[ic][im];
  }
 
  computeP(
    prob,
    useConstant
    );
}

computeP_selfDspin2() [2/2]

void Mela::computeP_selfDspin2	(	double	selfDGggcoupl_input[SIZE_GGG][2],
		double	selfDGvvcoupl_input[SIZE_GVV][2],
		float &	prob,
		bool	useConstant = true
	)

This is a function that calls Mela::computeP with preset quark, gluon, and vector boson couplings for spin 1, but sets default values for Zqq couplings.

See also

Mela::computeP, Mela::selfDGggcoupl, Mela::selfDGqqcoupl, Mela::selfDGvvcoupl

Warning

I (Mohit Srivastav) would personally avoid using this function, as it obfuscates the couplings you are really using .

It is always a good idea to set couplings yourself so that you know what you're really doing.

Parameters

[in]	selfDGggcoupl_input	an array corresponding to the couplings of Mela::selfDGggcoupl (between the "graviton" and gluons) using the indices in TCouplingsBase.hh
[in]	selfDGvvcoupl_input	an array corresponding to the couplings of Mela::selfDGvvcoupl (between the "graviton" and the W/Z bosons) using the indices in TCouplingsBase.hh
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1164 of file Mela.cc.

   {
  selfDGqqcoupl[0][0]=1.0; // Set these incorrectly, and you see left-right asymmetries in qqG (or nothing)
  selfDGqqcoupl[1][0]=1.0;
 
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_GGG; ic++) selfDGggcoupl[ic][im] = selfDGggcoupl_input[ic][im];
    for (int ic=0; ic<SIZE_GVV; ic++) selfDGvvcoupl[ic][im] = selfDGvvcoupl_input[ic][im];
  }
 
  computeP(
    prob,
    useConstant
    );
}

computePM4l()

void Mela::computePM4l	(	TVar::SuperMelaSyst	syst,
		float &	prob
	)

Definition at line 1864 of file Mela.cc.

                                                         {
  reset_PAux();
  prob=-99;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    bool hasFailed=false;
    int id_original[2][2];
    for (int iv=0; iv<2; iv++){
      MELAParticle* Vi = melaCand->getSortedV(iv);
      int ndau=Vi->getNDaughters();
      if (ndau!=2 || !(PDGHelpers::isAZBoson(Vi->id) || PDGHelpers::isAPhoton(Vi->id))){ hasFailed=true; break; } // Veto WW, ZG, GG
      for (int ivd=0; ivd<2; ivd++) id_original[iv][ivd]=Vi->getDaughter(ivd)->id;
    }
 
    if (!hasFailed){
      if (abs(id_original[0][0])==11 && abs(id_original[1][0])==11 && abs(id_original[0][1])==11 && abs(id_original[1][1])==11) super->SetDecayChannel("4e");
      else if (abs(id_original[0][0])==13 && abs(id_original[1][0])==13 && abs(id_original[0][1])==13 && abs(id_original[1][1])==13) super->SetDecayChannel("4mu");
      else if (
        (abs(id_original[0][0])==11 && abs(id_original[0][1])==11 && abs(id_original[1][0])==13 && abs(id_original[1][1])==13)
        ||
        (abs(id_original[0][0])==13 && abs(id_original[0][1])==13 && abs(id_original[1][0])==11 && abs(id_original[1][1])==11)
        ) super->SetDecayChannel("2e2mu");
      else{ if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computePM4l: SuperMELA is currently not implemented for decay states other than 4e. 4mu, 2e2mu." << endl; hasFailed=true; }
    }
 
    if (!hasFailed){
      double mZZ = melaCand->m();
      // currently only supported signal is ggH(0+), only supported background is summed paramterization
      if (syst == TVar::SMSyst_None){
        std::pair<double, double> m4lP = super->M4lProb(mZZ);
        if (myModel_ == TVar::HSMHiggs) prob = m4lP.first;
        else if (myModel_ == TVar::bkgZZ) prob = m4lP.second;
      }
      else{
        //systematics for p(m4l)
        float mZZtmp=mZZ;
        float meanErr=float(super->GetSigShapeSystematic("meanCB"));
        float sigmaErr=float(super->GetSigShapeSystematic("sigmaCB"));
        float sigmaCB=float(super->GetSigShapeParameter("sigmaCB"));
        if (syst == TVar::SMSyst_ScaleUp) mZZtmp = mZZ*(1.0+meanErr);
        else if (syst == TVar::SMSyst_ScaleDown) mZZtmp = mZZ*(1.0-meanErr);
        else if (syst == TVar::SMSyst_ResUp || syst ==  TVar::SMSyst_ResDown) mZZtmp=melaRandomNumber.Gaus(mZZ, sigmaErr*sigmaCB);
 
        std::pair<double, double> m4lP = super->M4lProb(mZZtmp);
        if (myModel_ == TVar::HSMHiggs) prob = m4lP.first;
        else if (myModel_ == TVar::bkgZZ) prob = m4lP.second;
      }
    }
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
}

computeProdDecP() [1/2]

void Mela::computeProdDecP	(	double	selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		double	selfDHwwcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		double	selfDaTQGCcoupl_input[SIZE_ATQGC][2],
		double	selfDAZffcoupl_input[SIZE_AZff][2],
		float &	prob,
		bool	useConstant = true
	)

computes the combined production and decay probability while taking in coupling arrays

Attention

The JHUGen Matrix element is not supported

See also

Calls Mela::computeProdDecP(prob, useConstant)

Mela::selfDHvvcoupl, Mela::selfDHwwcoupl, Mela::selfDaTQGCcoupl, Mela::selfDAZffcoupl

Parameters

[in]	selfDHvvcoupl_input	input for the couplings set in Mela::selfDHvvcoupl
[in]	selfDHwwcoupl_input	input for the couplings set in Mela::selfDHwwcoupl
[in]	selfDaTQGCcoupl_input	input for the couplings set in Mela::selfDaTQGCcoupl
[in]	selfDAZffcoupl_input	input for the couplings set in Mela::selfDAZffcoupl
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1430 of file Mela.cc.

   {
  for (int jh=0; jh<(int)nSupportedHiggses; jh++){
    for (int im=0; im<2; im++){
      for (int ic=0; ic<SIZE_HVV; ic++){
        selfDHzzcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im];
        selfDHwwcoupl[jh][ic][im] = selfDHwwcoupl_input[jh][ic][im]; // Just for extra protection since differentiate_HWW_HZZ is set to false.
      }
    }
  }
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_ATQGC; ic++) selfDaTQGCcoupl[ic][im] = selfDaTQGCcoupl_input[ic][im];
  }
  for (int im=0; im<2; im++){
    for (int ic=0; ic<SIZE_AZff; ic++) selfDAZffcoupl[ic][im] = selfDAZffcoupl_input[ic][im];
  }
  computeProdDecP(
    prob,
    useConstant
    );
}

computeProdDecP() [2/2]

void Mela::computeProdDecP	(	float &	prob,
		bool	useConstant = true
	)

computes the combined production and decay probability, with couplings set beforehand

Attention

The JHUGen Matrix element is not supported

The following production modes are supported using the MCFM Matrix Element (as named in TVar::Production) TVar::Had_WH TVar::Had_ZH TVar::Had_WH_S TVar::Had_ZH_S TVar::Had_WH_TU TVar::Had_ZH_TU TVar::Lep_ZH TVar::Lep_WH TVar::Lep_ZH_S TVar::Lep_WH_S TVar::Lep_ZH_TU TVar::Lep_WH_TU TVar::JJVBF TVar::JJEW TVar::JJEWQCD TVar::JJQCD TVar::JJVBF_S TVar::JJEW_S TVar::JJEWQCD_S TVar::JJQCD_S TVar::JJVBF_TU TVar::JJEW_TU TVar::JJEWQCD_TU TVar::JJQCD_TU

Parameters

[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1457 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeProdDecP" << endl;
  reset_PAux();
  melaCand = getCurrentCandidate();
 
  bool hasFailed = false;
  if (myME_ != TVar::MCFM){
    MELAout << "Mela::computeProdDecP ME is not supported for ME " << myME_ << endl;
    hasFailed = true;
  }
  if (
    !(
    myProduction_==TVar::Had_WH || myProduction_==TVar::Had_ZH
    || myProduction_==TVar::Had_WH_S || myProduction_==TVar::Had_ZH_S
    || myProduction_==TVar::Had_WH_TU || myProduction_==TVar::Had_ZH_TU
    || myProduction_==TVar::Lep_ZH || myProduction_==TVar::Lep_WH
    || myProduction_==TVar::Lep_ZH_S || myProduction_==TVar::Lep_WH_S
    || myProduction_==TVar::Lep_ZH_TU || myProduction_==TVar::Lep_WH_TU
    || myProduction_==TVar::JJVBF || myProduction_==TVar::JJEW || myProduction_==TVar::JJEWQCD || myProduction_==TVar::JJQCD
    || myProduction_==TVar::JJVBF_S || myProduction_==TVar::JJEW_S || myProduction_==TVar::JJEWQCD_S || myProduction_==TVar::JJQCD_S
    || myProduction_==TVar::JJVBF_TU || myProduction_==TVar::JJEW_TU || myProduction_==TVar::JJEWQCD_TU || myProduction_==TVar::JJQCD_TU
    )
    ){
    MELAout << "Mela::computeProdDecP production mode is not supported for production " << myProduction_ << endl;
    hasFailed = true;
  }
  if (!melaCand) hasFailed=true;
  if (hasFailed) prob=0;
  else{
    setSpinZeroCouplings();
    setATQGCCouplings();
    ZZME->computeProdXS_VVHVV(prob);
    if (useConstant) computeConstant(prob);
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeProdDecP" << endl;
}

computeProdP() [1/2]

void Mela::computeProdP	(	double	selfDHggcoupl_input[SIZE_HGG][2],
		double	selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		double	selfDHwwcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		float &	prob,
		bool	useConstant = true
	)

computes Production side probabilities while taking in coupling arrays

Attention

The MCFM Matrix element is not supported

See also

Calls Mela::computeProdP(prob, useConstant)

Mela:selfDHggcoupl, Mela::selfDHvvcoupl, Mela::selfDHwwcoupl

Parameters

[in]	selfDHggcoupl_input	input for the couplings set in Mela::setlfDHggcoupl
[in]	selfDHvvcoupl_input	input for the couplings set in Mela::selfDHvvcoupl
[in]	selfDHwwcoupl_input	input for the couplings set in Mela::selfDHwwcoupl
[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1501 of file Mela.cc.

   {
  for (int im=0; im<2; im++){ for (int ic=0; ic<SIZE_HGG; ic++) selfDHggcoupl[0][ic][im] = selfDHggcoupl_input[ic][im]; }
  for (int jh=0; jh<(int)nSupportedHiggses; jh++){
    for (int im=0; im<2; im++){
      for (int ic=0; ic<SIZE_HVV; ic++){
        selfDHzzcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im];
        selfDHwwcoupl[jh][ic][im] = selfDHwwcoupl_input[jh][ic][im]; // Just for extra protection since differentiate_HWW_HZZ is set to false.
      }
    }
  }
  computeProdP(
    prob,
    useConstant
    );
}

computeProdP() [2/2]

void Mela::computeProdP	(	float &	prob,
		bool	useConstant = true
	)

computes Production side probabilities using couplings set beforehand

Attention

The MCFM Matrix element is not supported

Parameters

[out]	prob	This is the value of the output probability - edited by reference
[in]	useConstant	This turns on the calculation of a corrective constant to different probabilities through Mela::getConstant. If you would like the "pure" MELA calculation to be run, set useConstant to false. By default true.

Definition at line 1522 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeProdP" << endl;
  if (myProduction_ == TVar::ttH || myProduction_ == TVar::bbH) computeProdP_ttH(prob, 2, 0, useConstant);
  else if (myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH || myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH || myProduction_ == TVar::GammaH) computeProdP_VH(prob, false, useConstant);
  else{
    reset_PAux();
 
    melaCand = getCurrentCandidate();
    if (melaCand){
      TLorentzVector nullFourVector(0, 0, 0, 0);
      bool isJet2Fake = false;
      MELACandidate* candOriginal = melaCand;
 
      unsigned int firstJetIndex=0;
      TLorentzVector jet1, higgs;
      TLorentzVector jet1massless(0, 0, 0, 0);
      TLorentzVector jet2massless(0, 0, 0, 0);
      higgs=melaCand->p4;
      if (myProduction_ == TVar::JJQCD || myProduction_ == TVar::JJVBF){
        int njets=0;
        {
          int ip=0;
          for (MELAParticle* tmpPart:melaCand->getAssociatedJets()){
            if (tmpPart->passSelection){
              njets++;
              if (njets==1){
                firstJetIndex = ip;
                jet1 = tmpPart->p4;
              }
            }
            ip++;
          }
        }
        if (njets==1){
          TUtil::scaleMomentumToEnergy(jet1, jet1massless);
          TUtil::computeFakeJet(jet1massless, higgs, jet2massless);
          isJet2Fake=true;
 
          // Scale jet2massless pz if necessary
          const double threshold = 1000.*LHCsqrts/2.;
          TLorentzVector pTotal = higgs+jet1massless+jet2massless;
          double sysZ = pTotal.Z();
          if (fabs(sysZ)>threshold){
            double maxpz2 = threshold - higgs.Z() - jet1massless.Z();
            if (fabs(maxpz2)>0.){
              double ratio = jet2massless.Z()/maxpz2;
              double absp=sqrt(pow(jet2massless.Pt(), 2)+pow(jet2massless.Z()*ratio, 2));
              if (myVerbosity_>=TVar::INFO) MELAout << "Mela::computeProdP, isJet2Fake=true case: Rescaling pz of fake jet by " << ratio << " and energy = " << absp << "." << endl;
              jet2massless.SetXYZT(jet2massless.X(), jet2massless.Y(), jet2massless.Z()*ratio, absp);
            }
            else{
              if (myVerbosity_>=TVar::INFO) MELAout << "Mela::computeProdP, isJet2Fake=true case: Unable to rescaling pz of fake jet since max(|pz|)<0. Setting to 0 with appropriate energy = pT = " << jet2massless.Pt() << "." << endl;
              jet2massless.SetXYZT(jet2massless.X(), jet2massless.Y(), 0., jet2massless.Pt());
            }
          }
        }
      }
 
      if (isJet2Fake){ // Do the integration
        MELACandidate* candCopy = melaCand->shallowCopy();
        MELAParticle* firstJet = candCopy->getAssociatedJet(firstJetIndex);
        firstJet->p4.SetXYZT(jet1massless.X(), jet1massless.Y(), jet1massless.Z(), jet1massless.T()); // Re-assign momenta of the first jet. Be careful, it changes candOriginal as well!
        MELAParticle fakeJet(0, jet2massless); // jet2massless is the unknown jet
        candCopy->addAssociatedJet(&fakeJet);
        setCurrentCandidate(candCopy);
 
        if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
        ZZME->computeProdXS_JJH(prob); // Higgs + 2 jets: SBF or WBF main probability
 
        int nGrid=11;
        std::vector<double> etaArray;
        std::vector<double> pArray;
        double eta_max = 10;
        if (jet2massless.Pt()>0.) eta_max = max(eta_max, 1.2*fabs(jet2massless.Eta()));
        double eta_min = -eta_max;
 
        for (int iter=0; iter<nGrid; iter++){
          float prob_temp=0;
 
          double jet2temp_eta = ((double)iter)*(eta_max-eta_min) / (((double)nGrid) - 1.) + eta_min;
          etaArray.push_back(jet2temp_eta);
          double jet2temp_sinh_eta = TMath::SinH(jet2temp_eta);
          double jet2temp_pz = jet2massless.Pt()*jet2temp_sinh_eta;
          fakeJet.p4.SetZ(jet2temp_pz);
          fakeJet.p4.SetX(jet2massless.X()); fakeJet.p4.SetY(jet2massless.Y()); fakeJet.p4.SetT(fakeJet.p4.P());
 
          // Skip case with invalid pz
          const double threshold = 1000.*LHCsqrts/2.;
          TLorentzVector pTotal = higgs+jet1massless+fakeJet.p4;
          double sys = (pTotal.T()+fabs(pTotal.Z()))/2.;
          if (fabs(sys)<threshold){
            if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
            ZZME->computeProdXS_JJH(prob_temp);
          }
          pArray.push_back((double)prob_temp);
        }
 
        const double grid_precision = 0.15;
        int ctr_iter=0;
        for (int iG=0; iG<nGrid-1; iG++){ // For each spacing, first compare the average of end points to spline value
          if (pArray[iG]==pArray[iG+1]) continue;
          if (etaArray[iG]==etaArray[iG+1]) continue;
 
          ctr_iter++;
 
          TGraph interpolator(nGrid, etaArray.data(), pArray.data());
          double derivative_first = (pArray[1]-pArray[0])/(etaArray[1]-etaArray[0]);
          double derivative_last = (pArray[nGrid-1]-pArray[nGrid-2])/(etaArray[nGrid-1]-etaArray[nGrid-2]);
          TSpline3 spline("spline", &interpolator, "b1e1", derivative_first, derivative_last);
          double x_middle = (etaArray[iG]+etaArray[iG+1])*0.5;
          double y_middle = (pArray[iG]+pArray[iG+1])*0.5;
          double y_sp = spline.Eval(x_middle);
          if (y_sp<0) y_sp = 0;
 
          std::vector<double>::iterator gridIt;
 
          if (fabs(y_sp-y_middle)<grid_precision*fabs(y_middle) || fabs(etaArray[iG+1]-etaArray[iG])<1e-3){
            gridIt = pArray.begin()+iG+1;
            pArray.insert(gridIt, y_sp);
            gridIt = etaArray.begin()+iG+1;
            etaArray.insert(gridIt, x_middle);
            iG++; // Pass to next bin
          }
          else{
            float prob_temp=0;
 
            double jet2temp_eta = x_middle;
            gridIt = etaArray.begin()+iG+1;
            etaArray.insert(gridIt, x_middle);
            double jet2temp_sinh_eta = TMath::SinH(jet2temp_eta);
            double jet2temp_pz = jet2massless.Pt()*jet2temp_sinh_eta;
            fakeJet.p4.SetZ(jet2temp_pz);
            fakeJet.p4.SetX(jet2massless.X()); fakeJet.p4.SetY(jet2massless.Y()); fakeJet.p4.SetT(fakeJet.p4.P());
 
            // Skip case with invalid pz
            const double threshold = 1000.*LHCsqrts/2.;
            TLorentzVector pTotal = higgs+jet1massless+fakeJet.p4;
            double sys = (pTotal.T()+fabs(pTotal.Z()))/2.;
            if (fabs(sys)<threshold){
              if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
              ZZME->computeProdXS_JJH(prob_temp);
            }
            gridIt = pArray.begin()+iG+1;
            pArray.insert(gridIt, (double)prob_temp);
            iG--; // Do not pass to next bin, repeat until precision is achieved.
          }
          nGrid++;
        }
 
        if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela::computeProdP: Number of iterations for JVBF eta integration: " << ctr_iter << endl;
 
        auxiliaryProb = 0;
        int iGFirst=0, iGLast=nGrid-1;
        for (int iG=1; iG<nGrid; iG++){
          if (pArray[iG]>0 && pArray[iG-1]==0){
            iGFirst = iG-1;
            break;
          }
        }
        for (int iG=nGrid-2; iG>=0; iG--){
          if (pArray[iG]>0 && pArray[iG+1]==0){
            iGLast = iG+1;
            break;
          }
        }
        double dEtaGrid = etaArray[iGLast] - etaArray[iGFirst];
        for (int iG=iGFirst; iG<iGLast-1; iG++){
          double dEta = etaArray[iG+1] - etaArray[iG];
          double sumProb = pArray[iG]+pArray[iG+1];
          sumProb *= 0.5;
          dEta = dEta/dEtaGrid;
          double addProb = sumProb*dEta;
          auxiliaryProb += (float)addProb;
        }
 
        firstJet->p4.SetXYZT(jet1.X(), jet1.Y(), jet1.Z(), jet1.T()); // Re-assign momenta of the first jet to original. Be careful, it changes candOriginal as well!
        delete candCopy; // Delete the shallow copy
        setCurrentCandidate(candOriginal);
        melaCand = getCurrentCandidate();
        if (myVerbosity_>=TVar::DEBUG){
          if (melaCand!=candOriginal) MELAerr << "Mela::computeProdP: melaCand!=candOriginal at the end of the fake jet scenario!" << endl;
        }
 
        if (fabs(prob)>0) auxiliaryProb /= prob;
      }
      else{
        if (myProduction_ == TVar::JJQCD || myProduction_ == TVar::JJVBF){
          if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
          ZZME->computeProdXS_JJH(prob); // Higgs + 2 jets: SBF or WBF
        }
        else if (myProduction_ == TVar::JQCD){
          // No anomalous couplings are implemented in HJ
          ZZME->computeProdXS_JH(prob); // Higgs + 1 jet; only SM is supported for now.
        }
      }
      if (useConstant) computeConstant(prob);
    }
 
    reset_SelfDCouplings();
    reset_CandRef();
  }
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeProdP" << endl;
}

computeProdP_ttH()

void Mela::computeProdP_ttH	(	float &	prob,
		int	topProcess = 2,
		int	topDecay = 0,
		bool	useConstant = true
	)

Definition at line 1778 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeProdP_ttH" << endl;
  reset_PAux();
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
    ZZME->computeProdXS_ttH(prob,topProcess, topDecay);
    if (useConstant) computeConstant(prob);
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeProdP_ttH" << endl;
}

computeProdP_VH() [1/2]

void Mela::computeProdP_VH	(	double	selfDHvvcoupl_input[nSupportedHiggses][SIZE_HVV][2],
		float &	prob,
		bool	includeHiggsDecay = false,
		bool	useConstant = true
	)

Definition at line 1730 of file Mela.cc.

   {
  // Dedicated function for VH ME
  selfDHggcoupl[0][0][0] = 1.0; // Don't set this, and you might get 0 in the future for ggVH.
  for (int jh=0; jh<(int)nSupportedHiggses; jh++){
    for (int im=0; im<2; im++){
      for (int ic=0; ic<SIZE_HVV; ic++){
  //remove this line, it's here for counting purposes.  Zpffcoupl
        selfDHzzcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im];
        selfDHwwcoupl[jh][ic][im] = selfDHvvcoupl_input[jh][ic][im]; // Just for extra protection since differentiate_HWW_HZZ is set to false.
      }
    }
  }
 
  computeProdP_VH(
    prob,
    includeHiggsDecay,
    useConstant
    );
}

computeProdP_VH() [2/2]

void Mela::computeProdP_VH	(	float &	prob,
		bool	includeHiggsDecay = false,
		bool	useConstant = true
	)

Definition at line 1754 of file Mela.cc.

   {
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeProdP_VH" << endl;
  reset_PAux();
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    if (myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH || myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH || myProduction_ == TVar::GammaH){
      if (myModel_ == TVar::SelfDefine_spin0) setSpinZeroCouplings();
      ZZME->computeProdXS_VH(prob, includeHiggsDecay); // VH
 
      if (useConstant) computeConstant(prob);
    }
  }
 
  reset_SelfDCouplings();
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeProdP_VH" << endl;
}

computeTTHAngles()

void Mela::computeTTHAngles	(	int	topDecay,
		float &	mT1,
		float &	mW1,
		float &	mT2,
		float &	mW2,
		float &	costheta1,
		float &	costheta2,
		float &	Phi,
		float &	costhetastar,
		float &	Phi1,
		float &	hbb,
		float &	hWW,
		float &	Phibb,
		float &	Phi1bb,
		float &	hWplusf,
		float &	PhiWplusf,
		float &	hWminusf,
		float &	PhiWminusf
	)

Definition at line 895 of file Mela.cc.

 {
  using TVar::simple_event_record;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeVBFAngles" << endl;
 
  mT1=mT2=mW1=mW2
  =costheta1=costheta2=Phi=costhetastar=Phi1
  =hbb=hWW=Phibb=Phi1bb
  =hWplusf=PhiWplusf
  =hWminusf=PhiWminusf=0;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
 
    int partIncCode;
    int nRequested_Tops=1;
    int nRequested_Antitops=1;
    if (topDecay>0) partIncCode=TVar::kUseAssociated_UnstableTops; // Look for unstable tops
    else partIncCode=TVar::kUseAssociated_StableTops; // Look for stable tops
 
    simple_event_record mela_event;
    mela_event.AssociationCode=partIncCode;
    mela_event.nRequested_Tops=nRequested_Tops;
    mela_event.nRequested_Antitops=nRequested_Antitops;
    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
    SimpleParticleCollection_t& mothers = mela_event.pMothers;
    //SimpleParticleCollection_t& aparts = mela_event.pAssociated;
    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
    SimpleParticleCollection_t topDaughters; topDaughters.reserve(3);
    SimpleParticleCollection_t antitopDaughters; antitopDaughters.reserve(3);
 
    if (topDecay>0){
      // Daughters are assumed to have been ordered as b, Wf, Wfb already.
      for (unsigned int itd=0; itd<mela_event.pTopDaughters.at(0).size(); itd++) topDaughters.push_back(mela_event.pTopDaughters.at(0).at(itd));
      for (unsigned int itd=0; itd<mela_event.pAntitopDaughters.at(0).size(); itd++) antitopDaughters.push_back(mela_event.pAntitopDaughters.at(0).at(itd));
    }
    else{
      for (unsigned int itop=0; itop<mela_event.pStableTops.size(); itop++) topDaughters.push_back(mela_event.pStableTops.at(itop));
      for (unsigned int itop=0; itop<mela_event.pStableAntitops.size(); itop++) antitopDaughters.push_back(mela_event.pStableAntitops.at(itop));
    }
 
    if (topDecay==0 && (mela_event.pStableTops.size()<1 || mela_event.pStableAntitops.size()<1)){
      if (myVerbosity_>=TVar::ERROR) MELAerr
        << "TUtil::TTHiggsMatEl: Number of stable tops (" << mela_event.pStableTops.size() << ")"
        << " and number of stable antitops (" << mela_event.pStableAntitops.size() << ")"
        << " in ttH process are not 1!" << endl;
      return;
    }
    else if (topDecay>0 && (mela_event.pTopDaughters.size()<1 || mela_event.pAntitopDaughters.size()<1)){
      if (myVerbosity_>=TVar::ERROR) MELAerr
        << "Mela::computeTTHAngles: Number of set of top daughters (" << mela_event.pTopDaughters.size() << ")"
        << " and number of set of antitop daughters (" << mela_event.pAntitopDaughters.size() << ")"
        << " in ttH process are not 1!" << endl;
      return;
    }
    else if (topDecay>0 && (mela_event.pTopDaughters.at(0).size()!=3 || mela_event.pAntitopDaughters.at(0).size()!=3)){
      if (myVerbosity_>=TVar::ERROR) MELAerr
        << "Mela::computeTTHAngles: Number of top daughters (" << mela_event.pTopDaughters.at(0).size() << ")"
        << " and number of antitop daughters (" << mela_event.pAntitopDaughters.at(0).size() << ")"
        << " in ttH process are not 3!" << endl;
      return;
    }
    if (topDaughters.size()<3){ for (size_t ip=topDaughters.size(); ip<3; ip++) topDaughters.emplace_back(-9000, nullVector); }
    if (antitopDaughters.size()<3){ for (size_t ip=antitopDaughters.size(); ip<3; ip++) antitopDaughters.emplace_back(-9000, nullVector); }
 
    // Make sure there are exactly 4 daughters, null or not
    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
      SimpleParticle_t& firstPart = daughters.at(0);
      firstPart.first=25;
      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
      daughters.erase(daughters.begin()+4, daughters.end());
    }
    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
    else if (daughters.size()==2){
      daughters.push_back(SimpleParticle_t(-9000, nullVector));
      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
    }
 
    // Compute masses
    {
      TLorentzVector pT;
      TLorentzVector pW;
      for (size_t ip=0; ip<topDaughters.size(); ip++){
        if (ip>0){ pT += topDaughters.at(ip).second; pW += topDaughters.at(ip).second; }
        else pT += topDaughters.at(ip).second;
      }
      mT1=pT.M();
      mW1=pW.M();
    }
    {
      TLorentzVector pT;
      TLorentzVector pW;
      for (size_t ip=0; ip<antitopDaughters.size(); ip++){
        if (ip>0){ pT += antitopDaughters.at(ip).second; pW += antitopDaughters.at(ip).second; }
        else pT += antitopDaughters.at(ip).second;
      }
      mT2=pT.M();
      mW2=pW.M();
    }
 
    TUtil::computeTTHAngles(
      costhetastar, costheta1, costheta2, Phi, Phi1,
      hbb, hWW, Phibb, Phi1bb,
      hWplusf, PhiWplusf,
      hWminusf, PhiWminusf,
 
      daughters.at(0).second, daughters.at(0).first,
      daughters.at(1).second, daughters.at(1).first,
      daughters.at(2).second, daughters.at(2).first,
      daughters.at(3).second, daughters.at(3).first,
 
      topDaughters.at(0).second, topDaughters.at(0).first,
      topDaughters.at(1).second, topDaughters.at(1).first,
      topDaughters.at(2).second, topDaughters.at(2).first,
 
      antitopDaughters.at(0).second, antitopDaughters.at(0).first,
      antitopDaughters.at(1).second, antitopDaughters.at(1).first,
      antitopDaughters.at(2).second, antitopDaughters.at(2).first,
 
      &(mothers.at(0).second), mothers.at(0).first,
      &(mothers.at(1).second), mothers.at(1).first
    );
 
    // Protect against NaN
    if (!std::isfinite(costhetastar)) costhetastar=0;
    if (!std::isfinite(costheta1)) costheta1=0;
    if (!std::isfinite(costheta2)) costheta2=0;
    if (!std::isfinite(Phi)) Phi=0;
    if (!std::isfinite(Phi1)) Phi1=0;
 
    if (!std::isfinite(hbb)) hbb=0;
    if (!std::isfinite(hWW)) hWW=0;
    if (!std::isfinite(Phibb)) Phibb=0;
    if (!std::isfinite(Phi1bb)) Phi1bb=0;
 
    if (!std::isfinite(hWplusf)) hWplusf=0;
    if (!std::isfinite(PhiWplusf)) PhiWplusf=0;
 
    if (!std::isfinite(hWminusf)) hWminusf=0;
    if (!std::isfinite(PhiWminusf)) PhiWminusf=0;
 
    if (myVerbosity_>=TVar::DEBUG) MELAout
      << "Mela::computeTTHAngles: (h1, h2, Phi, hs, Phi1) = "
      << costheta1 << ", " << costheta2 << ", " << Phi << ", "
      << costhetastar << ", " << Phi1 << endl;
  }
  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeTTHAngles: No possible melaCand in TEvtProb to compute angles." << endl;
 
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeTTHAngles" << endl;
}

computeVBFAngles()

void Mela::computeVBFAngles	(	float &	Q2V1,
		float &	Q2V2,
		float &	costheta1,
		float &	costheta2,
		float &	Phi,
		float &	costhetastar,
		float &	Phi1
	)

computes the decay angles for Vector Boson Fusion (VBF) -> H -> ZZ as defined in Figure 1 of arXiv:1208.4018

See also

Calls TUtil::computeVBFAngles to fully calculate angles

Use TUtil::computeVBFAngles if you would like to calculate the angles through inputting your own 4-vectors and quantities rather than through the MELA event loop

Attention

This function requires there to be an input event already defined

Warning

This function will edit all values inplace. Every value is technically an output.

The selection for m1 and m2 are based upon which of the leptons is closer to the z mass. $m_1 \leq m_2$.

Parameters

[out]	Q2V1	The mass of the first decay particle as reconstructed by 2 of the leptons
[out]	m2V2	The mass of the first decay particle as reconstructed by 2 of the leptons
[out]	costheta1	$\cos(\theta_1)$
[out]	costheta2	$\cos(\theta_2)$
[out]	Phi	$\phi$
[out]	costhetastar	$\cos(\theta*)$
[out]	Phi1	$\phi_1$

Definition at line 620 of file Mela.cc.

 {
  using TVar::simple_event_record;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeVBFAngles" << endl;
 
  Q2V1=0; Q2V2=0; costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
 
    int nRequested_AssociatedJets=2;
    int partIncCode=TVar::kUseAssociated_Jets; // Only use associated partons in the pT=0 frame boost
    simple_event_record mela_event;
    mela_event.AssociationCode=partIncCode;
    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
    SimpleParticleCollection_t& mothers = mela_event.pMothers;
    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
    if ((int)aparts.size()!=nRequested_AssociatedJets){ if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVBFAngles: Number of associated particles is not 2!" << endl; return; }
 
    // Make sure there are exactly 4 daughters, null or not
    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
      SimpleParticle_t& firstPart = daughters.at(0);
      firstPart.first=25;
      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
      daughters.erase(daughters.begin()+4, daughters.end());
    }
    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
    else if (daughters.size()==2){
      daughters.push_back(SimpleParticle_t(-9000, nullVector));
      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
    }
 
    if (myVerbosity_>=TVar::DEBUG){
      MELAout << "Mela::computeVBFAngles: Giving the following particles to TUtil::computeVBFAngles:" << endl;
      for (unsigned int i=0; i<std::min(daughters.size(), (SimpleParticleCollection_t::size_type) 4); i++) MELAout << daughters.at(i) << endl;
      for (unsigned int i=0; i<std::min(aparts.size(), (SimpleParticleCollection_t::size_type) 2); i++) MELAout << aparts.at(i) << endl;
      for (unsigned int i=0; i<std::min(mothers.size(), (SimpleParticleCollection_t::size_type) 2); i++) MELAout << mothers.at(i) << endl;
    }
 
    TUtil::computeVBFAngles(
      costhetastar, costheta1, costheta2, Phi, Phi1, Q2V1, Q2V2,
      daughters.at(0).second, daughters.at(0).first,
      daughters.at(1).second, daughters.at(1).first,
      daughters.at(2).second, daughters.at(2).first,
      daughters.at(3).second, daughters.at(3).first,
      aparts.at(0).second, aparts.at(0).first,
      aparts.at(1).second, aparts.at(1).first,
      &(mothers.at(0).second), mothers.at(0).first,
      &(mothers.at(1).second), mothers.at(1).first
    );
 
    // Protect against NaN
    if (!std::isfinite(costhetastar)) costhetastar=0;
    if (!std::isfinite(costheta1)) costheta1=0;
    if (!std::isfinite(costheta2)) costheta2=0;
    if (!std::isfinite(Phi)) Phi=0;
    if (!std::isfinite(Phi1)) Phi1=0;
    if (!std::isfinite(Q2V1)) Q2V1=0;
    if (!std::isfinite(Q2V2)) Q2V2=0;
 
    if (myVerbosity_>=TVar::DEBUG) MELAout
      << "Mela::computeVBFAngles: (Q2_1, Q2_2, h1, h2, Phi, hs, Phi1) = "
      << Q2V1 << ", " << Q2V2 << ", "
      << costheta1 << ", " << costheta2 << ", " << Phi << ", "
      << costhetastar << ", " << Phi1 << endl;
  }
  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVBFAngles: No possible melaCand in TEvtProb to compute angles." << endl;
 
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeVBFAngles" << endl;
}

computeVBFAngles_ComplexBoost()

void Mela::computeVBFAngles_ComplexBoost	(	float &	Q2V1,
		float &	Q2V2,
		float &	costheta1_real,
		float &	costheta1_imag,
		float &	costheta2_real,
		float &	costheta2_imag,
		float &	Phi,
		float &	costhetastar,
		float &	Phi1
	)

Definition at line 702 of file Mela.cc.

 {
  using TVar::simple_event_record;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeVBFAngles_ComplexBoost" << endl;
 
  Q2V1=0; Q2V2=0; costheta1_real=0; costheta2_real=0; costheta1_imag=0; costheta2_imag=0; Phi=0; costhetastar=0; Phi1=0;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
 
    int nRequested_AssociatedJets=2;
    int partIncCode=TVar::kUseAssociated_Jets; // Only use associated partons in the pT=0 frame boost
    simple_event_record mela_event;
    mela_event.AssociationCode=partIncCode;
    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
    SimpleParticleCollection_t& mothers = mela_event.pMothers;
    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
    if ((int) aparts.size()!=nRequested_AssociatedJets){ if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVBFAngles_ComplexBoost: Number of associated particles is not 2!" << endl; return; }
 
    // Make sure there are exactly 4 daughters, null or not
    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
      SimpleParticle_t& firstPart = daughters.at(0);
      firstPart.first=25;
      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
      daughters.erase(daughters.begin()+4, daughters.end());
    }
    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
    else if (daughters.size()==2){
      daughters.push_back(SimpleParticle_t(-9000, nullVector));
      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
    }
 
    TUtil::computeVBFAngles_ComplexBoost(
      costhetastar, costheta1_real, costheta1_imag, costheta2_real, costheta2_imag, Phi, Phi1, Q2V1, Q2V2,
      daughters.at(0).second, daughters.at(0).first,
      daughters.at(1).second, daughters.at(1).first,
      daughters.at(2).second, daughters.at(2).first,
      daughters.at(3).second, daughters.at(3).first,
      aparts.at(0).second, aparts.at(0).first,
      aparts.at(1).second, aparts.at(1).first,
      &(mothers.at(0).second), mothers.at(0).first,
      &(mothers.at(1).second), mothers.at(1).first
    );
 
    // Protect against NaN
    if (!std::isfinite(costhetastar)) costhetastar=0;
    if (!std::isfinite(costheta1_real)) costheta1_real=0;
    if (!std::isfinite(costheta2_real)) costheta2_real=0;
    if (!std::isfinite(costheta1_imag)) costheta1_imag=0;
    if (!std::isfinite(costheta2_imag)) costheta2_imag=0;
    if (!std::isfinite(Phi)) Phi=0;
    if (!std::isfinite(Phi1)) Phi1=0;
    if (!std::isfinite(Q2V1)) Q2V1=0;
    if (!std::isfinite(Q2V2)) Q2V2=0;
 
    if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela::computeVBFAngles_ComplexBoost: result = " << Q2V1 << ", " << Q2V2
                                           << ", " << costheta1_real << " + " << costheta1_imag << "i, "
                                           << costheta2_real << " + " << costheta2_imag << "i, " << Phi << ", "
                                           << costhetastar << ", " << Phi1 << endl;
    if (myVerbosity_>=TVar::DEBUG) MELAout
      << "Mela::computeVBFAngles_ComplexBoost: (Q2_1, Q2_2, h1, h2, Phi, hs, Phi1) = "
      << Q2V1 << ", " << Q2V2 << ", "
      << costheta1_real << " + " << costheta1_imag << "i, "
      << costheta2_real << " + " << costheta2_imag << "i, "
      << Phi << ", " << costhetastar << ", " << Phi1 << endl;
  }
  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVBFAngles_ComplexBoost: No possible melaCand in TEvtProb to compute angles." << endl;
 
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeVBFAngles_ComplexBoost" << endl;
}

computeVHAngles()

void Mela::computeVHAngles	(	float &	mVstar,
		float &	mV,
		float &	costheta1,
		float &	costheta2,
		float &	Phi,
		float &	costhetastar,
		float &	Phi1
	)

computes the decay angles for Vector Boson Fusion (VBF) -> H -> ZZ as defined in Figure 1 of arXiv:1208.4018

See also

Calls TUtil::computeVHAngles to fully calculate angles

Use TUtil::computeVHAngles if you would like to calculate the angles through inputting your own 4-vectors and quantities rather than through the MELA event loop

Attention

This function requires there to be an input event already defined

Warning

This function will edit all values inplace. Every value is technically an output.

The selection for m1 and m2 are based upon which of the leptons is closer to the z mass. $m_1 \leq m_2$.

Parameters

[out]	mVstar	The mass of the virtual Z boson in ZH production
[out]	mV	The mass of the Z in decay
[out]	costheta1	$\cos(\theta_1)$
[out]	costheta2	$\cos(\theta_2)$
[out]	Phi	$\phi$
[out]	costhetastar	$\cos(\theta*)$
[out]	Phi1	$\phi_1$

Definition at line 786 of file Mela.cc.

 {
  using TVar::simple_event_record;
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: Begin computeVHAngles" << endl;
 
  mVstar = 0; mV = 0; costheta1=0; costheta2=0; Phi=0; costhetastar=0; Phi1=0;
 
  melaCand = getCurrentCandidate();
  if (melaCand){
    TLorentzVector const nullVector(0, 0, 0, 0);
 
    if (!(myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH || myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH || myProduction_ == TVar::GammaH)){
      if (myVerbosity_>=TVar::ERROR) MELAerr << "Mela::computeVHAngles: Production is not supported! " << ProductionName(myProduction_) << endl;
      return;
    }
 
    int nRequested_AssociatedJets=0;
    int nRequested_AssociatedLeptons=0;
    int nRequested_AssociatedPhotons=0;
    int AssociationVCompatibility=0;
    int partIncCode=TVar::kNoAssociated; // Just to avoid warnings
    if (myProduction_ == TVar::Had_ZH || myProduction_ == TVar::Had_WH){ // Only use associated partons
      partIncCode=TVar::kUseAssociated_Jets;
      nRequested_AssociatedJets=2;
    }
    else if (myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Lep_WH){ // Only use associated leptons(+)neutrinos
      partIncCode=TVar::kUseAssociated_Leptons;
      nRequested_AssociatedLeptons=2;
    }
    else if (myProduction_ == TVar::GammaH){ // Only use associated photon
      partIncCode=TVar::kUseAssociated_Photons;
      nRequested_AssociatedPhotons=1;
    }
    if (myProduction_==TVar::Lep_WH || myProduction_==TVar::Had_WH) AssociationVCompatibility=24;
    else if (myProduction_==TVar::Lep_ZH || myProduction_==TVar::Had_ZH) AssociationVCompatibility=23;
    else if (myProduction_==TVar::GammaH) AssociationVCompatibility=22;
    simple_event_record mela_event;
    mela_event.AssociationCode=partIncCode;
    mela_event.AssociationVCompatibility=AssociationVCompatibility;
    mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
    mela_event.nRequested_AssociatedLeptons=nRequested_AssociatedLeptons;
    mela_event.nRequested_AssociatedPhotons=nRequested_AssociatedPhotons;
    TUtil::GetBoostedParticleVectors(melaCand, mela_event, myVerbosity_);
    SimpleParticleCollection_t& mothers = mela_event.pMothers;
    SimpleParticleCollection_t& aparts = mela_event.pAssociated;
    SimpleParticleCollection_t& daughters = mela_event.pDaughters;
 
    if ((aparts.size()<(unsigned int) (nRequested_AssociatedJets+nRequested_AssociatedLeptons) && myProduction_!=TVar::GammaH) || (aparts.size()<(unsigned int) nRequested_AssociatedPhotons && myProduction_==TVar::GammaH)){
      if (myVerbosity_>=TVar::ERROR){
        MELAerr << "Mela::computeVHAngles: Number of associated particles (" << aparts.size() << ") is less than ";
        if (myProduction_!=TVar::GammaH) MELAerr << (nRequested_AssociatedJets+nRequested_AssociatedLeptons);
        else MELAerr << nRequested_AssociatedPhotons;
        MELAerr << endl;
      }
      return;
    }
 
    // Make sure there are exactly 4 daughters, null or not
    if (daughters.size()>4){ // Unsupported size, default to undecayed Higgs
      SimpleParticle_t& firstPart = daughters.at(0);
      firstPart.first=25;
      for (auto it=daughters.cbegin()+1; it!=daughters.cend(); it++){ firstPart.second = firstPart.second + it->second; }
      daughters.erase(daughters.begin()+4, daughters.end());
    }
    if (daughters.size()%2==1){ for (unsigned int ipar=daughters.size(); ipar<4; ipar++) daughters.push_back(SimpleParticle_t(-9000, nullVector)); }
    else if (daughters.size()==2){
      daughters.push_back(SimpleParticle_t(-9000, nullVector));
      daughters.insert(daughters.begin()+1, SimpleParticle_t(-9000, nullVector));
    }
 
    TUtil::computeVHAngles(
      costhetastar, costheta1, costheta2, Phi, Phi1, mVstar, mV,
      daughters.at(0).second, daughters.at(0).first,
      daughters.at(1).second, daughters.at(1).first,
      daughters.at(2).second, daughters.at(2).first,
      daughters.at(3).second, daughters.at(3).first,
      aparts.at(0).second, aparts.at(0).first,
      aparts.at(1).second, aparts.at(1).first,
      &(mothers.at(0).second), mothers.at(0).first,
      &(mothers.at(1).second), mothers.at(1).first
    );
 
    // Protect against NaN
    if (!std::isfinite(costhetastar)) costhetastar=0;
    if (!std::isfinite(costheta1)) costheta1=0;
    if (!std::isfinite(costheta2)) costheta2=0;
    if (!std::isfinite(Phi)) Phi=0;
    if (!std::isfinite(Phi1)) Phi1=0;
    if (!std::isfinite(mVstar)) mVstar=0;
    if (!std::isfinite(mV)) mV=0;
 
    if (myVerbosity_>=TVar::DEBUG) MELAout
      << "Mela::computeVHAngles: (mVstar, mV, h1, h2, Phi, hs, Phi1) = "
      << mVstar << ", " << mV << ", " << costheta1 << ", " << costheta2 << ", " << Phi << ", "
      << costhetastar << ", " << Phi1 << endl;
  }
  else if (myVerbosity_>=TVar::DEBUG) MELAerr << "Mela::computeVHAngles: No possible melaCand in TEvtProb to compute angles." << endl;
 
  reset_CandRef();
  if (myVerbosity_>=TVar::DEBUG) MELAout << "Mela: End computeVHAngles" << endl;
}

getConstant()

void Mela::getConstant

(

float &

prob

)

This function returns a multiplicative constant to the matrix element calculation.

Attention

if useConstant=True in the following functions the constant will be applied to the matrix element:

• computeP_selfDspin0
• computeP_selfDspin1
• computeP_selfDspin1
• computeP_selfDspin2
• computeP_selfDspin2
• computeP
• computeProdDecP
• computeProdDecP
• computeProdP
• computeProdP
• computeProdP_VH
• computeProdP_VH
• computeProdP_ttH
• computePM4l

Parameters

[out]

prob

This is the number that will be multiplied inplace by the constant.

Definition at line 2168 of file Mela.cc.

{ prob = getIORecord()->getMEConst(); }

getCurrentCandidate()

MELACandidate * Mela::getCurrentCandidate

(

)

Gets the current MELA top-level (input) candList object.

See also

Wrapper for ZZMatrixElement::get_CurrentCandidate, which is a wrapper for TEvtProb::GetCurrentCandidate, which returns the protected element TEvtProb::melaCand

Definition at line 546 of file Mela.cc.

{ return ZZME->get_CurrentCandidate(); }

getCurrentCandidateIndex()

int Mela::getCurrentCandidateIndex

(

)

Returns the index of the current MELA candidate - returns -1 if there is no candidate to be found.

See also

Wrapper for ZZMatrixElement::get_CurrentCandidateIndex, which is a wrapper for TEvtProb::GetCurrentCandidateIndex

Definition at line 547 of file Mela.cc.

{ return ZZME->get_CurrentCandidateIndex(); }

getHiggsWidthAtPoleMass()

double Mela::getHiggsWidthAtPoleMass

(

double

mass

)

Returns the width of the Higgs at a given pole mass as a calculation.

See also

Wrapper for ZZMatrixElement::get_HiggsWidthAtPoleMass, which is a wrapper for TEvtProb::GetHiggsWidthAtPoleMass, which is a wrapper for MELAHXSWidth::HiggsWidth

Attention

This function does an independent calculation for the mass of the Higgs at a certain mass, and does not rely on any other set MELA values

Parameters

[in]

mass

The value of the Higgs pole mass you want in GeV

Definition at line 517 of file Mela.cc.

{ return ZZME->get_HiggsWidthAtPoleMass(mass); }

getIORecord()

MelaIO * Mela::getIORecord

(

)

Returns the MELAIO object, and by consequence, the entire parton-by-parton matrix element record.

See also

wrapper for ZZMatrixElement::get_IORecord, which is a wrapper for TEvtProb::GetIORecord, which is a wrapper for MelaIO::getRef

Definition at line 544 of file Mela.cc.

{ return ZZME->get_IORecord(); }

getMeasurablesRRV()

RooSpin::modelMeasurables Mela::getMeasurablesRRV

(

)

Returns a RooSpin::modelMeasureables object containing many observable quantities.

Attention

The parameters in question are:

• $\cos(\theta_1)$
• $\cos(\theta_2)$
• $\phi$
• mass of the first Z boson ($M_{Z1}$)
• mass of the second Z boson ($M_{Z2}$)
• mass of the full particle ($M_{ZZ}$)
• $\cos(\theta*)$
• $\phi_1$
• Y

Remarks

these angles are defined in Figure 1 of arXiv:1001.3396 as well as (described more generally) in Figure 1 of arXiv:1208.4018

Definition at line 528 of file Mela.cc.

                                               {
  RooSpin::modelMeasurables measurables;
  measurables.h1 = costheta1_rrv;
  measurables.h2 = costheta2_rrv;
  measurables.Phi = phi_rrv;
  measurables.m1 = z1mass_rrv;
  measurables.m2 = z2mass_rrv;
  measurables.m12 = mzz_rrv;
  measurables.hs = costhetastar_rrv;
  measurables.Phi1 = phi1_rrv;
  measurables.Y = Y_rrv;
  return measurables;
}

getNCandidates()

int Mela::getNCandidates

(

)

Returns the size of the candidate list TEvtProb::candList.

See also

Wrapper for ZZMatrixElement::get_NCandidates, which is a wrapper for TEvtProb::GetNCandidates, which returns the size of TEvtProb::candList

Definition at line 548 of file Mela.cc.

{ return ZZME->get_NCandidates(); }

getPAux()

void Mela::getPAux

(

float &

prob

)

Definition at line 554 of file Mela.cc.

{ prob = auxiliaryProb; }

getPrimaryMass()

double Mela::getPrimaryMass

(

int

ipart

)

A function to get the current primary EW/QCD parameters from MELA.

See also

Wrapper for ZZMatrixElement::get_PrimaryMass, which is a Wrapper for TEvtProb::GetPrimaryMass, which calls TUtil::GetMass.

Warning

(Higgs mass/width used in the Matrix Element could be different). Refer to setMelaHiggsMassWidth or setMELAHiggsMass for explicit references to those.

Attention

You can input negative ids to edit antiparticles (i.e. 24 -> W+, -24 -> W-)

In MCFM both the particle and its antiparticle's mass are changed at the same time. In JHUGen this is not the case - they are separate entries.

The particles that you can change when using the MCFM matrix element are as as follow:

absolute value of ipart	Particle Name
8	t Prime quark (4th generation)
7	b Prime quark (4th generation)
6	Top Quark
5	Bottom Quark
4	Charm Quark
3	Strange Quark
2	Up Quark
1	Down Quark
11	Electron
13	Muon
15	Tau
23	Z Boson
24	W Boson
25	Higgs Boson

The particles that you can change when using the JHUGen matrix element are as follow:

ipart	Particle Name	Sign Sensitive
0 OR 21	Gluon	No
1	Down Quark	Yes
2	Up Quark	Yes
3	Strange Quark	Yes
4	Charm Quark	Yes
5	Bottom Quark	Yes
6	Top Quark	Yes
11	Electron	Yes
22	Photon	No
23	Z Boson	No
24	W Boson	Yes
13	Muon	Yes
15	Tau	Yes
12	Electron Neutrino	Yes
14	Muon Neutrino	Yes
16	Tau Neutrino	Yes
25	Higgs Boson	No
32	Z Prime	No
33	Z Prime 2	No
34	W Prime	Yes

Parameters

[in]

ipart

the particle you would like to get the mass of.

Definition at line 515 of file Mela.cc.

{ return ZZME->get_PrimaryMass(ipart); }

getPrimaryWidth()

double Mela::getPrimaryWidth

(

int

ipart

)

A function to get the current primary EW/QCD parameters from MELA.

See also

Wrapper for ZZMatrixElement::get_PrimaryMass, which is a Wrapper for TEvtProb::GetPrimaryMass, which calls TUtil::GetMass.

Warning

(Higgs mass/width used in the Matrix Element could be different). Refer to setMelaHiggsMassWidth or setMELAHiggsMass for explicit references to those.

Attention

You can input negative ids to edit antiparticles (i.e. 24 -> W+, -24 -> W-)

In MCFM both the particle and its antiparticle's mass are changed at the same time. In JHUGen this is not the case - they are separate entries.

The particles that you can change when using the MCFM matrix element are as as follow:

absolute value of ipart	Particle Name
15	Tau
23	Z Boson
24	W Boson
25	Higgs Boson

The particles that you can change when using the JHUGen matrix element are as follow:

ipart	Particle Name	Sign Sensitive
0 OR 21	Gluon	No
1	Down Quark	Yes
2	Up Quark	Yes
3	Strange Quark	Yes
4	Charm Quark	Yes
5	Bottom Quark	Yes
6	Top Quark	Yes
11	Electron	Yes
22	Photon	No
23	Z Boson	No
24	W Boson	Yes
13	Muon	Yes
15	Tau	Yes
12	Electron Neutrino	Yes
14	Muon Neutrino	Yes
16	Tau Neutrino	Yes
25	Higgs Boson	No
32	Z Prime	No
33	Z Prime 2	No
34	W Prime	Yes

Parameters

[in]

ipart

the particle you would like to get the width of.

Definition at line 516 of file Mela.cc.

{ return ZZME->get_PrimaryWidth(ipart); }

getTopCandidateCollection()

std::vector< MELATopCandidate_t * > * Mela::getTopCandidateCollection

(

)

Same as getNCandidates but specifically for Top Quark Candidates.

See also

Wrapper for ZZMatrixElement::get_TopCandidateCollection, which is a wrapper for TEvtProb::GetTopCandidates, which returns TEvtProb::topCandList

Definition at line 549 of file Mela.cc.

{ return ZZME->get_TopCandidateCollection(); }

getVerbosity()

TVar::VerbosityLevel Mela::getVerbosity

(

)

Gets the current verbosity level for MELA.

Returns

a TVar::VerbosityLevel describing the verbosity level for MELA. This is a number from 0 to 5, and corresponds to the values in TVar::VerbosityLevel.

Definition at line 333 of file Mela.cc.

{ return myVerbosity_; }

getXPropagator()

void Mela::getXPropagator	(	TVar::ResonancePropagatorScheme	scheme,
		float &	prop
	)

Definition at line 1799 of file Mela.cc.

                                                                          {
  prop=0.;
  melaCand = getCurrentCandidate();
  if (melaCand) ZZME->get_XPropagator(scheme, prop);
  reset_CandRef();
}

resetInputEvent()

void Mela::resetInputEvent

(

)

Resets the event in preparation for the next iteration of the event loop.

See also

Wrapper for ZZMatrixElement::reset_InputEvent, which is a wrapper for TEvtProb::ResetInputEvent

Attention

Without resetting the input event at the end of each for loop, behavior could be unexpected!

It is important to call this at the end of every event loop iteration to clean up TEvtProb!

Definition at line 356 of file Mela.cc.

{ ZZME->reset_InputEvent(); }

resetMass()

void Mela::resetMass	(	double	inmass,
		int	ipart
	)

Resets the mass for a particle that is an electroweak parameter according to its id.

See also

Wrapper for ZZMatrixElement::reset_Mass, which is a wrapper for TEvtProb::ResetMass, which is a wrapper for TUtil::SetMass

This finally interfaces with either the function SetMass in mod_parameters.F90 for JHUGen, accessed through .TModParameters.hh

in mod_parameters.F90 the conversion from ipart to the values in JHUGen's internal code happens with the function convertLHEreverse(Part)

or the masses in the JHUGen-MCFM library at /src/Inc/masses.F or at TMCFM::spinzerohiggs_anomcoupl (for tPrime and bPrime)

Warning

It is not recommended to use this function to edit the mass of the Higgs in MCFM. Please use either setMelaHiggsMass or setMelaHiggsMassWidth for that.

Attention

You can input negative ids to edit antiparticles (i.e. 24 -> W+, -24 -> W-)

In MCFM both the particle and its antiparticle's mass are changed at the same time. In JHUGen this is not the case - they are separate entries.

The particles that you can change when using the MCFM matrix element are as as follow:

absolute value of ipart	Particle Name
8	t Prime quark (4th generation)
7	b Prime quark (4th generation)
6	Top Quark
5	Bottom Quark
4	Charm Quark
3	Strange Quark
2	Up Quark
1	Down Quark
11	Electron
13	Muon
15	Tau
23	Z Boson
24	W Boson
25	Higgs Boson

The particles that you can change when using the JHUGen matrix element are as follow:

ipart	Particle Name	Sign Sensitive
0 OR 21	Gluon	No
1	Down Quark	Yes
2	Up Quark	Yes
3	Strange Quark	Yes
4	Charm Quark	Yes
5	Bottom Quark	Yes
6	Top Quark	Yes
11	Electron	Yes
22	Photon	No
23	Z Boson	No
24	W Boson	Yes
13	Muon	Yes
15	Tau	Yes
12	Electron Neutrino	Yes
14	Muon Neutrino	Yes
16	Tau Neutrino	Yes
25	Higgs Boson	No
32	Z Prime	No
33	Z Prime 2	No
34	W Prime	Yes

Parameters

[in]	inmass	the mass that you want in GeV
[in]	ipart	the particle whose mass you wish to change

Definition at line 508 of file Mela.cc.

{ ZZME->reset_Mass(inmass, ipart); }

resetMCFM_EWKParameters()

void Mela::resetMCFM_EWKParameters	(	double	ext_Gf,
		double	ext_aemmz,
		double	ext_mW,
		double	ext_mZ,
		double	ext_xW,
		int	ext_ewscheme = 3
	)

Resets the electroweak parameters back to their defaults.

See also

Wrapper for ZZMatrixElement::reset_MCFM_EWKParameters, which is a wrapper for TEvtProb::ResetMCFM_EWKParameters, which calls TUtil::SetEwkCouplingParameters

Definition at line 511 of file Mela.cc.

                                                                                                                                {
  ZZME->reset_MCFM_EWKParameters(ext_Gf, ext_aemmz, ext_mW, ext_mZ, ext_xW, ext_ewscheme);
}

resetQuarkMasses()

void Mela::resetQuarkMasses

(

)

Resets the masses of each quark to their original values.

See also

Wrapper for ZZMatrixElement::reset_QuarkMasses, which is a wrapper for TEvtProb::ResetQuarkMasses, which calls TEvtProb::ResetMass

Warning

You should run this command each time you would like to reset the mass of a quark to mitigate unexpected behavior

Attention

The quarks that are reset are as follow, with their default masses in GeV listed as well:

Quark	Mass
Down	0.001
Up	0.005
Strange	0.1
Charm	1.275
Bottom	4.75
Top	173.2
B Prime	1e5
T Prime	1e5

Definition at line 510 of file Mela.cc.

{ ZZME->reset_QuarkMasses(); }

resetWidth()

void Mela::resetWidth	(	double	inwidth,
		int	ipart
	)

Resets the width for a particle that is an electroweak parameter according to its id.

See also

Wrapper for ZZMatrixElement::reset_Width, which is a wrapper for TEvtProb::ResetWidth, which is a wrapper for TUtil::SetDecayWidth

This finally interfaces with either the function SetDecayWidth in mod_parameters.F90 for JHUGen, accessed through .TModParameters.hh

in mod_parameters.F90 the conversion from ipart to the values in JHUGen's internal code happens with the function convertLHEreverse(Part)

or the masses in the JHUGen-MCFM library at /src/Inc/masses.F or at TMCFM::spinzerohiggs_anomcoupl (for tPrime and bPrime)

Warning

It is not recommended to use this function to edit the mass of the Higgs in MCFM. Please use either setMelaHiggsWidth or setMelaHiggsMassWidth for that.

Attention

You can input negative ids to edit antiparticles (i.e. 24 -> W+, -24 -> W-)

In MCFM both the particle and its antiparticle's mass are changed at the same time. In JHUGen this is not the case - they are separate entries.

The particles that you can change when using the MCFM matrix element are as as follow:

absolute value of ipart	Particle Name
15	Tau
23	Z Boson
24	W Boson
25	Higgs Boson

The particles that you can change when using the JHUGen matrix element are as follow:

ipart	Particle Name	Sign Sensitive
0 OR 21	Gluon	No
1	Down Quark	Yes
2	Up Quark	Yes
3	Strange Quark	Yes
4	Charm Quark	Yes
5	Bottom Quark	Yes
6	Top Quark	Yes
11	Electron	Yes
22	Photon	No
23	Z Boson	No
24	W Boson	Yes
13	Muon	Yes
15	Tau	Yes
12	Electron Neutrino	Yes
14	Muon Neutrino	Yes
16	Tau Neutrino	Yes
25	Higgs Boson	No
32	Z Prime	No
33	Z Prime 2	No
34	W Prime	Yes

Parameters

[in]	inmass	the mass that you want in GeV
[in]	ipart	the particle whose mass you wish to change

Definition at line 509 of file Mela.cc.

{ ZZME->reset_Width(inwidth, ipart); }

setCandidateDecayMode()

void Mela::setCandidateDecayMode

(

TVar::CandidateDecayMode

mode

)

Sets the decay mode for your event.

See also

Wrapper for ZZMatrixElement::set_CandidateDecayMode, which is a wrapper for TEvtProb::SetCandidateDecayMode

Parameters

[in]

mode

The decay mode you would like picked from TVar::CandidateDecayMode

Definition at line 340 of file Mela.cc.

{ ZZME->set_CandidateDecayMode(mode); }

setCurrentCandidate()

void Mela::setCurrentCandidate

(

MELACandidate *

cand

)

Switches the candidate that you are working on to another candidate object specified.

See also

Wrapper for ZZMatrixElement::set_CurrentCandidate, which is a wrapper for TEvtProb::SetCurrentCandidate.

Parameters

[in]

cand

The MELACandidate object you would like to switch to

Definition at line 342 of file Mela.cc.

{ ZZME->set_CurrentCandidate(cand); }

setCurrentCandidateFromIndex()

void Mela::setCurrentCandidateFromIndex

(

unsigned int

icand

)

Switches the candidate that you are working on to another candidate based off of an index.

See also

Wrapper for ZZMatrixElement::set_CurrentCandidateFromIndex, which is a wrapper for TEvtProb::SetCurrentCandidateFromIndex.

Parameters

[in]

icand

The index of the candidate that you would like to switch to

Definition at line 341 of file Mela.cc.

{ ZZME->set_CurrentCandidateFromIndex(icand); }

setInputEvent()

void Mela::setInputEvent	(	SimpleParticleCollection_t *	pDaughters,
		SimpleParticleCollection_t *	pAssociated = 0,
		SimpleParticleCollection_t *	pMothers = 0,
		bool	isGen = false
	)

Sets the input event for MELA. MELA cannot run without this.

See also

Wrapper for ZZMatrixElement::set_InputEvent, which is a wrapper for TEvtProb::SetInputEvent, which calls TUtil::ConvertVectorFormat

Attention

An input event must be set for each event in an event loop, otherwise MELA will throw a segmentation error

Parameters

[in]	pDaughters	A SimpleParticleCollection_t of particle daughters (decay products)
[in]	pAssociated	A SimpleParticleCollection_t of associated particles (i.e. jets), by default 0 (no jets)
[in]	pMothers	A SimpleParticleCollection_t of particle mothers (i.e. gluons), by default 0 (reco data contains no mother information)
[in]	isGen	A boolean signifying whether the event in question is a Gen event or a reco event, by default false (reco)

Definition at line 343 of file Mela.cc.

   {
  ZZME->set_InputEvent(
    pDaughters,
    pAssociated,
    pMothers,
    isGen
    );
}

setMelaHiggsMass()

void Mela::setMelaHiggsMass	(	double	myHiggsMass,
		int	index = 0
	)

Sets the mass of your chosen Higgs.

See also

Wrapper for ZZMatrixElement::set_mHiggs.

Parameters

[in]	myHiggsMass	This is the mass of the Higgs that you would like, in GeV
[in]	index	This is either 0 or 1, depending on if you want to change the first or second resonance. By default it is 0.

Definition at line 336 of file Mela.cc.

{ ZZME->set_mHiggs(myHiggsMass, index); }

setMelaHiggsMassWidth()

void Mela::setMelaHiggsMassWidth	(	double	myHiggsMass,
		double	myHiggsWidth,
		int	index
	)

a combination of setMelaHiggsMass and setMelaHiggsWidth.

See also

Wrapper for ZZMatrixElement::set_mHiggs_wHiggs.

Parameters

[in]	myHiggsMass	This is the mass of the Higgs that you would like, in GeV
[in]	myHiggsWidth	This is the mass of the Higgs that you would like, in GeV
[in]	index	This is either 0 or 1, depending on if you want to change the first or second resonance. By default it is 0.

Definition at line 338 of file Mela.cc.

{ ZZME->set_mHiggs_wHiggs(myHiggsMass, myHiggsWidth, index); }

setMelaHiggsWidth()

void Mela::setMelaHiggsWidth	(	double	myHiggsWidth = -1,
		int	index = 0
	)

Sets the width of your chosen Higgs.

See also

Wrapper for ZZMatrixElement::set_wHiggs

Parameters

[in]	myHiggsWidth	This is the mass of the Higgs that you would like, in GeV
[in]	index	This is either 0 or 1, depending on if you want to change the first or second resonance. By default it is 0.

Definition at line 337 of file Mela.cc.

{ ZZME->set_wHiggs(myHiggsWidth, index); }

setMelaLeptonInterference()

void Mela::setMelaLeptonInterference

(

TVar::LeptonInterference

myLepInterf = TVar::DefaultLeptonInterf

)

Sets the MELA Lepton Interference.

Parameters

[in]

myLepInterf

sets the myLepInterf_ variable to one found in TVar::LeptonInterface through ZZMatrixElement::set_LeptonInterface. TVar::DefaultLeptonInterf by default.

Definition at line 339 of file Mela.cc.

{ myLepInterf_=myLepInterf; ZZME->set_LeptonInterference(myLepInterf); }

setMelaPrimaryHiggsMass()

void Mela::setMelaPrimaryHiggsMass

(

double

myHiggsMass

)

Sets the mass of the "primary" higgs.

See also

Wrapper for the function ZZMatrixElement::set_PrimaryHiggsMass, which is a wrapper for TEvtProb::SetPrimaryHiggsMass

Parameters

[in]

myHiggsMass

This is the mass of the Higgs that you would like, in GeV

Attention

The primary Higgs is the first resonance - nominally MELA can have 2 resonances when working through JHUGen-MCFM

Remarks

This function is effectively the same as setMelaHiggsMass(mass, 0)

Definition at line 335 of file Mela.cc.

{ ZZME->set_PrimaryHiggsMass(myHiggsMass); }

setProcess()

void Mela::setProcess	(	TVar::Process	myModel,
		TVar::MatrixElement	myME,
		TVar::Production	myProduction
	)

Sets the process, matrix element, and production that MELA is to use for this event. Calls ZZMatrixElement::set_Process, which calls TEvtProb::SetProcess.

Attention

Remember to set the process for each event, otherwise the MELA event loop will throw a segmentation error.

Parameters

[in]	myModel	a TVar for the Process you would like, as defined in TVar::Process
[in]	myME	a TVar for the matrix element you would like, as defined in TVar::MatrixElement
[in]	myProduction	a TVar for the production mode you would like, as defined in TVar::Production

Definition at line 310 of file Mela.cc.

                                                                                             {
  myME_ = myME;
  myProduction_ = myProduction;
  // In case s-channel processes are passed for JHUGen ME, flip them back to JHUGen-specific productions.
  if (myME_==TVar::JHUGen){
    if (myProduction_==TVar::Had_ZH_S) myProduction_=TVar::Had_ZH;
    else if (myProduction_==TVar::Had_WH_S) myProduction_=TVar::Had_WH;
    else if (myProduction_==TVar::Lep_ZH_S) myProduction_=TVar::Lep_ZH;
    else if (myProduction_==TVar::Lep_WH_S) myProduction_=TVar::Lep_WH;
    else if (myProduction_==TVar::JJVBF_S) myProduction_=TVar::JJVBF;
    else if (myProduction_==TVar::JJQCD_S) myProduction_=TVar::JJQCD;
  }
  myModel_ = myModel;
  if (ZZME!=0) ZZME->set_Process(myModel_, myME_, myProduction_);
}

setRemoveJetMasses()

void Mela::setRemoveJetMasses

(

bool

MasslessLeptonSwitch = true

)

either permits or forbids massive jets.

See also

Wrapper for the function TUtil::applyJetMassCorrection

Parameters

[in]

MasslessLeptonSwitch

Whether you would like your jets to be massive, by default true.

Definition at line 520 of file Mela.cc.

{ TUtil::applyJetMassCorrection(MasslessLeptonSwitch); }

setRemoveLeptonMasses()

void Mela::setRemoveLeptonMasses

(

bool

MasslessLeptonSwitch = true

)

either permits or forbids massive leptons.

See also

Wrapper for the function TUtil::applyLeptonMassCorrection

Parameters

[in]

MasslessLeptonSwtich

Whether you would like your leptons to be massive, by default true.

Definition at line 519 of file Mela.cc.

{ TUtil::applyLeptonMassCorrection(MasslessLeptonSwitch); }

setRenFacScaleMode()

void Mela::setRenFacScaleMode	(	TVar::EventScaleScheme	renormalizationSch,
		TVar::EventScaleScheme	factorizationSch,
		double	ren_sf,
		double	fac_sf
	)

Sets the renormalization and the factorization schemes.

See also

Wrapper for ZZMatrixElement::set_RenFacScaleMode, which is a wrapper for TEvtProb::SetRenFacScaleMode, which edits the TVar::event_scales_type struct in TEvtProb.event_scales

Parameters

[in]	renormalizationSch	This is the renormalization scheme that you are picking from TVar::EventScaleScheme
[in]	factorizationSch	This is the factorization scheme that you are picking from TVar::EventScaleScheme
[in]	ren_sf	This is the renormalization scale factor that you would like
[in]	fac_sf	This is the scale factor for the factorization scale that you would like

Definition at line 521 of file Mela.cc.

                                                                                                                                         {
  ZZME->set_RenFacScaleMode(renormalizationSch, factorizationSch, ren_sf, fac_sf);
}

setTempCandidate()

void Mela::setTempCandidate	(	SimpleParticleCollection_t *	pDaughters,
		SimpleParticleCollection_t *	pAssociated = 0,
		SimpleParticleCollection_t *	pMothers = 0,
		bool	isGen = false
	)

Sets a temporary MELA candidate, by setting melaCand in Xcal2 to a temporary candidate without pushing this candidate to the candList of Xcal2.

See also

Wrapper for ZZMatrixElement::set_TempCandidate

Parameters

[in]	pDaughters	A SimpleParticleCollection_t of particle daughters (decay products)
[in]	pAssociated	A SimpleParticleCollection_t of associated particles (i.e. jets), by default 0 (no jets)
[in]	pMothers	A SimpleParticleCollection_t of particle mothers (i.e. gluons), by default 0 (reco data contains no mother information)
[in]	isGen	A boolean signifying whether the event in question is a Gen event or a reco event, by default false (reco)

Definition at line 357 of file Mela.cc.

   { ZZME->set_TempCandidate(pDaughters, pAssociated, pMothers); }

setVerbosity()

void Mela::setVerbosity

(

TVar::VerbosityLevel

verbosity_ = TVar::ERROR

)

Sets the verbosity for MELA outside of the initial constructor.

Parameters

[in]

verbosity_

The verbosity of MELA that you desire, as defined in TVar::ERROR.

Definition at line 325 of file Mela.cc.

                                                    {
  myVerbosity_=verbosity_;
  if (ZZME) ZZME->set_Verbosity(myVerbosity_);
  if (super) super->SetVerbosity((myVerbosity_>=TVar::DEBUG));
  if (superDijet) superDijet->SetVerbosity(myVerbosity_);
  if (ggSpin0Model) ggSpin0Model->setVerbosity(myVerbosity_);
  if (spin2Model) spin2Model->setVerbosity(myVerbosity_);
}

Member Data Documentation

costheta1_rrv

RooRealVar* Mela::costheta1_rrv

Definition at line 939 of file Mela.h.

costheta2_rrv

RooRealVar* Mela::costheta2_rrv

Definition at line 940 of file Mela.h.

costhetastar_rrv

RooRealVar* Mela::costhetastar_rrv

Definition at line 938 of file Mela.h.

ggSpin0Model

ScalarPdfFactory_HVV* Mela::ggSpin0Model

Definition at line 947 of file Mela.h.

melaRandomNumber

TRandom3 Mela::melaRandomNumber

Definition at line 934 of file Mela.h.

mzz_rrv

RooRealVar* Mela::mzz_rrv

Definition at line 935 of file Mela.h.

pdf

RooAbsPdf* Mela::pdf

Definition at line 946 of file Mela.h.

phi1_rrv

RooRealVar* Mela::phi1_rrv

Definition at line 942 of file Mela.h.

phi_rrv

RooRealVar* Mela::phi_rrv

Definition at line 941 of file Mela.h.

qqZZmodel

RooqqZZ_JHU_ZgammaZZ_fast* Mela::qqZZmodel

Definition at line 950 of file Mela.h.

spin1Model

VectorPdfFactory* Mela::spin1Model

Definition at line 948 of file Mela.h.

spin2Model

TensorPdfFactory_ppHVV* Mela::spin2Model

Definition at line 949 of file Mela.h.

super

SuperMELA* Mela::super

Definition at line 952 of file Mela.h.

upFrac_rrv

RooRealVar* Mela::upFrac_rrv

Definition at line 944 of file Mela.h.

Y_rrv

RooRealVar* Mela::Y_rrv

Definition at line 943 of file Mela.h.

z1mass_rrv

RooRealVar* Mela::z1mass_rrv

Definition at line 936 of file Mela.h.

z2mass_rrv

RooRealVar* Mela::z2mass_rrv

Definition at line 937 of file Mela.h.

---

The documentation for this class was generated from the following files:

• MELA/interface/Mela.h
• MELA/src/Mela.cc