TUtil.hh

Go to the documentation of this file.

// March 28 2011
// S. Jindariani ([email protected])
// Y. Gao ([email protected])
// K. Burkett ([email protected])
 
 
#ifndef ZZ_COMMON
#define ZZ_COMMON
#include <string>
#include <vector>
#include <complex>
// MelaIO class
#include "MelaIO.h"
// Couplings classes
#include "TCouplings.hh"
// MCFM utilities
#include "TMCFMUtils.hh"
// Mod_Parameters
#include "TModParameters.hh"
// NNPDF Driver for JHUGen
#include "TNNPDFDriver.hh"
// Mod Kinematics
#include "TModKinematics.hh"
// JHUGenMELA
#include "TModJHUGen.hh"
#include "TModJHUGenMELA.hh"
// Higgs + 0 jet
#include "TModHiggsMatEl.hh"
#include "TModGravitonMatEl.hh"
#include "TModZprimeMatEl.hh"
// Higgs + 1/2 jets
#include "TModHiggsJJMatEl.hh"
#include "TModHiggsJMatEl.hh"
// VH
#include "TModVHiggsMatEl.hh"
// ttH
#include "TModTTBHMatEl.hh"
// ROOT includes
#include "TLorentzVector.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TFile.h"
#include "TF1.h"
 
 
namespace TUtil{
  extern bool forbidMassiveLeptons;
  extern bool forbidMassiveJets;
  extern TVar::FermionMassRemoval LeptonMassScheme;
  extern TVar::FermionMassRemoval JetMassScheme;
 
  void applyLeptonMassCorrection(bool flag=true);
  void applyJetMassCorrection(bool flag=true);
  void setLeptonMassScheme(TVar::FermionMassRemoval scheme=TVar::ConserveDifermionMass);
  void setJetMassScheme(TVar::FermionMassRemoval scheme=TVar::ConserveDifermionMass);
  // This version makes the masses of p1 and p2 to be m1 and m2, leaving p1+p2 unchanged.
  void constrainedRemovePairMass(TLorentzVector& p1, TLorentzVector& p2, double m1=0, double m2=0);
  // This version simply scales momentum to match energy for the desired mass.
  void scaleMomentumToEnergy(const TLorentzVector& massiveJet, TLorentzVector& masslessJet, double mass=0);
  // Function that has generic removal features
  std::pair<TLorentzVector, TLorentzVector> removeMassFromPair(
    TLorentzVector const& jet1, int const& jet1Id,
    TLorentzVector const& jet2, int const& jet2Id,
    double m1=0, double m2=0
    );
  // Function that adjusts top daughter kinematics
  void adjustTopDaughters(SimpleParticleCollection_t& daughters); // Daughters are arranged as b, Wf, Wfb
  // Compute a fake jet from the massless jets
  void computeFakeJet(TLorentzVector const& realJet, TLorentzVector const& others, TLorentzVector& fakeJet); // Input massive + higgs -> output massless fake jet
 
  // TLorentzVector::Boost in complex plane
  std::pair<TLorentzVector, TLorentzVector> ComplexBoost(TVector3 const& beta, TLorentzVector const& p4);
 
  void computeAngles(
    float& costhetastar,
    float& costheta1,
    float& costheta2,
    float& Phi,
    float& Phi1,
    TLorentzVector Z1_lept1, int Z1_lept1Id,
    TLorentzVector Z1_lept2, int Z1_lept2Id,
    TLorentzVector Z2_lept1, int Z2_lept1Id,
    TLorentzVector Z2_lept2, int Z2_lept2Id
    );
  void computeAnglesCS(
    float const& pbeam,
    float& costhetastar,
    float& costheta1,
    float& costheta2,
    float& Phi,
    float& Phi1,
    TLorentzVector Z1_lept1, int Z1_lept1Id,
    TLorentzVector Z1_lept2, int Z1_lept2Id,
    TLorentzVector Z2_lept1, int Z2_lept1Id,
    TLorentzVector Z2_lept2, int Z2_lept2Id
    );
  // Angles of associated production
  void computeVBFAngles(
    float& costhetastar,
    float& costheta1,
    float& costheta2,
    float& Phi,
    float& Phi1,
    float& Q2V1,
    float& Q2V2,
    TLorentzVector p4M11, int Z1_lept1Id,
    TLorentzVector p4M12, int Z1_lept2Id,
    TLorentzVector p4M21, int Z2_lept1Id,
    TLorentzVector p4M22, int Z2_lept2Id,
    TLorentzVector jet1, int jet1Id,
    TLorentzVector jet2, int jet2Id,
    TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
    TLorentzVector* injet2=0, int injet2Id=0
    );
  void computeVBFAngles_ComplexBoost(
    float& costhetastar,
    float& costheta1_real, float& costheta1_imag,
    float& costheta2_real, float& costheta2_imag,
    float& Phi,
    float& Phi1,
    float& Q2V1,
    float& Q2V2,
    TLorentzVector p4M11, int Z1_lept1Id,
    TLorentzVector p4M12, int Z1_lept2Id,
    TLorentzVector p4M21, int Z2_lept1Id,
    TLorentzVector p4M22, int Z2_lept2Id,
    TLorentzVector jet1, int jet1Id,
    TLorentzVector jet2, int jet2Id,
    TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
    TLorentzVector* injet2=0, int injet2Id=0
    );
  void computeVHAngles(
    float& costhetastar,
    float& costheta1,
    float& costheta2,
    float& Phi,
    float& Phi1,
    float& m1,
    float& m2,
    TLorentzVector p4M11, int Z1_lept1Id,
    TLorentzVector p4M12, int Z1_lept2Id,
    TLorentzVector p4M21, int Z2_lept1Id,
    TLorentzVector p4M22, int Z2_lept2Id,
    TLorentzVector jet1, int jet1Id,
    TLorentzVector jet2, int jet2Id,
    TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
    TLorentzVector* injet2=0, int injet2Id=0
    );
  // hs: cos theta*
  // h[x]: cos theta_X
  void computeTTHAngles(
    // ttH system
    float& hs,
    float& hincoming,
    float& hTT,
    float& PhiTT,
    float& Phi1,
 
    // tt system
    float& hbb,
    float& hWW,
    float& Phibb,
    float& Phi1bb,
 
    // Wplus system
    float& hWplusf,
    float& PhiWplusf,
 
    // Wminus system
    float& hWminusf,
    float& PhiWminusf,
 
    // H->ZZ/ZA/AA vectors. For A decays, leave second vector (0,0,0,0) with id=-9000
    // Z/A1
    TLorentzVector p4M11, int Z1_lept1Id,
    TLorentzVector p4M12, int Z1_lept2Id,
    // Z/A2
    TLorentzVector p4M21, int Z2_lept1Id,
    TLorentzVector p4M22, int Z2_lept2Id,
    // TT system
    TLorentzVector b, int bId,
    TLorentzVector Wplusf, int WplusfId,
    TLorentzVector Wplusfb, int WplusfbId,
    TLorentzVector bbar, int bbarId,
    TLorentzVector Wminusf, int WminusfId,
    TLorentzVector Wminusfb, int WminusfbId,
    TLorentzVector* injet1=0, int injet1Id=0, // Gen. partons in lab frame
    TLorentzVector* injet2=0, int injet2Id=0
  );
 
  // Parameter settings
  void SetEwkCouplingParameters(double ext_Gf, double ext_aemmz, double ext_mW, double ext_mZ, double ext_xW, int ext_ewscheme);
  void SetMass(double inmass, int ipart);
  void SetDecayWidth(double inwidth, int ipart);
  double GetMass(int ipart);
  double GetDecayWidth(int ipart);
  double GetMass(const MELAParticle* part);
  double GetDecayWidth(const MELAParticle* part);
  void GetMassWidth(int ipart, double& m, double& ga);
  void GetMassWidth(const MELAParticle* part, double& m, double& ga);
  void SetCKMElements(double* invckm_ud, double* invckm_us, double* invckm_cd, double* invckm_cs, double* invckm_ts, double* invckm_tb, double* invckm_ub=0, double* invckm_cb=0, double* invckm_td=0);
  void SetMadgraphCKMElements(double ckmlambda=0.2265, double ckma=0.79, double ckmrho=0.141, double ckmeta=0);
  double GetCKMElement(int iquark, int jquark);
  std::complex<double> GetMadgraphCKMElement(int iquark, int jquark);
  double InterpretScaleScheme(const TVar::Production& production, const TVar::MatrixElement& matrixElement, const TVar::EventScaleScheme& scheme, TLorentzVector p[mxpart]);
  void SetAlphaS(double& Q_ren, double& Q_fac, double multiplier_ren, double multiplier_fac, int mynloop, int mynflav, std::string mypartons); // Q_ren/fac -> Q_ren/fac * multiplier_ren/fac
  void GetAlphaS(double* alphas_, double* alphasmz_); // Get last alpha_s value set
 
  // chooser.f split into 3 different functions
  bool MCFM_chooser(
    const TVar::Process& process, const TVar::Production& production, const TVar::LeptonInterference& leptonInterf,
    const TVar::VerbosityLevel& verbosity,
    const TVar::simple_event_record& mela_event
    );
  bool MCFM_SetupParticleCouplings(
    const TVar::Process& process, const TVar::Production& production,
    const TVar::VerbosityLevel& verbosity,
    const TVar::simple_event_record& mela_event,
    std::vector<int>* partOrder, std::vector<int>* apartOrder
    );
  TString GetMCFMParticleLabel(const int& pid, bool useQJ, bool useExtendedConventions);
 
  // JHUGen-specific wrappers
  void InitJHUGenMELA(const char* pathtoPDFSet, int PDFMember, double collider_sqrts);
  void SetJHUGenHiggsMassWidth(double MReso, double GaReso);
  void SetJHUGenDistinguishWWCouplings(bool doAllow);
  void ResetAmplitudeIncludes();
 
  // Spin-0 couplings
  void SetMadgraphSpinZeroCouplings(SpinZeroCouplings const* Hcouplings);
  void SetMCFMSpinZeroCouplings(bool useBSM, SpinZeroCouplings const* Hcouplings, bool forceZZ);
  void SetMCFMaTQGCCouplings(bool useBSM, aTQGCCouplings const* couplings);
  void SetJHUGenSpinZeroVVCouplings(double Hvvcoupl[SIZE_HVV][2], double Hvvpcoupl[SIZE_HVV][2], double Hvpvpcoupl[SIZE_HVV][2], int Hvvcoupl_cqsq[SIZE_HVV_CQSQ], double HvvLambda_qsq[SIZE_HVV_LAMBDAQSQ][SIZE_HVV_CQSQ], bool useWWcoupl);
  void SetJHUGenSpinZeroGGCouplings(double Hggcoupl[SIZE_HGG][2]);
  void SetJHUGenSpinZeroQQCouplings(double Hqqcoupl[SIZE_HQQ][2]);
  // Spin-1 couplings
  void SetJHUGenSpinOneCouplings(double Zqqcoupl[SIZE_ZQQ][2], double Zvvcoupl[SIZE_ZVV][2]);
  // Spin-2 couplings
  void SetJHUGenSpinTwoCouplings(double Gacoupl[SIZE_GGG][2], double Gvvcoupl[SIZE_GVV][2], double Gvvpcoupl[SIZE_GVV][2], double Gvpvpcoupl[SIZE_GVV][2], double qLeftRightcoupl[SIZE_GQQ][2]);
  // Vprime / contact couplings
  void SetJHUGenVprimeContactCouplings(double Zpffcoupl[SIZE_Vpff][2], double Wpffcoupl[SIZE_Vpff][2]);
  // AZff
  void SetMCFMAZffCouplings(bool useBSM, AZffCouplings const* Zcouplings);
  void SetJHUGenAZffCouplings(bool needAZff, double AZffcoupl[SIZE_AZff][2]);
 
  // PS cuts, unused
  bool MCFM_masscuts(double s[][mxpart], const TVar::Process& process);
  bool MCFM_smalls(double s[][mxpart], int npart);
 
  // ME computations
  double SumMatrixElementPDF(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement, const TVar::LeptonInterference& leptonInterf,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    TVar::VerbosityLevel verbosity
    );
  double JHUGenMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    TVar::VerbosityLevel verbosity
    );
  double MadgraphMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    TVar::VerbosityLevel verbosity
    );
  double HJJMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    TVar::VerbosityLevel verbosity
    );
  double VHiggsMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    bool includeHiggsDecay,
    TVar::VerbosityLevel verbosity
    );
  double TTHiggsMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    int topDecay, int topProcess,
    TVar::VerbosityLevel verbosity
    );
  double BBHiggsMatEl(
    const TVar::Process& process, const TVar::Production& production, const TVar::MatrixElement& matrixElement,
    TVar::event_scales_type* event_scales, MelaIO* RcdME,
    const double& EBEAM,
    int botProcess,
    TVar::VerbosityLevel verbosity
    );
 
  int WipeMEArray(const TVar::Process& process, const TVar::Production& production, const int id[mxpart], double msq[nmsq][nmsq], const TVar::VerbosityLevel& verbosity);
  bool CheckPartonMomFraction(const TLorentzVector& p0, const TLorentzVector& p1, double xx[2], const double& EBEAM, const TVar::VerbosityLevel& verbosity);
  void ComputePDF(const TLorentzVector& p0, const TLorentzVector& p1, double fx1[nmsq], double fx2[nmsq], const double& EBEAM, const TVar::VerbosityLevel& verbosity);
  double SumMEPDF(const TLorentzVector& p0, const TLorentzVector& p1, double msq[nmsq][nmsq], MelaIO* RcdME, const double& EBEAM, const TVar::VerbosityLevel& verbosity);
 
  // Propagator reweighting
  double ResonancePropagator(double const& sqrts, TVar::ResonancePropagatorScheme scheme);
 
  // Boost the particles with or without associated ones to pT=0 frame and return std::vectors filled with (id, momentum) pairs
  void GetBoostedParticleVectors(
    MELACandidate* melaCand,
    TVar::simple_event_record& mela_event,
    TVar::VerbosityLevel verbosity=TVar::DEBUG
    );
 
  // Convert vectors of simple particles to MELAParticles and create a MELACandidate
  // The output lists could be members of TEvtProb directly.
  MELACandidate* ConvertVectorFormat(
    // Inputs
    SimpleParticleCollection_t* pDaughters, // Cannot be 0
    SimpleParticleCollection_t* pAssociated, // Allowed to be 0
    SimpleParticleCollection_t* pMothers, // Allowed to be 0
    bool isGen,
    // Outputs
    std::vector<MELAParticle*>* particleList,
    std::vector<MELACandidate*>* candList
    );
  // Convert the vector of three body decay daughters (as simple particles) to MELAParticles and create a MELAThreeBodyDecayCandidate
  // The output lists could be members of TEvtProb directly.
  MELAThreeBodyDecayCandidate* ConvertThreeBodyDecayCandidate(
    // Input
    SimpleParticleCollection_t* tbdDaughters,
    // Outputs
    std::vector<MELAParticle*>* particleList,
    std::vector<MELAThreeBodyDecayCandidate*>* tbdCandList
    );
  void PrintCandidateSummary(MELACandidate* cand);
  void PrintCandidateSummary(TVar::simple_event_record* cand);
 
}
 
#endif