Merge branch 'master' into v1.0

Conflicts:
	CMakeLists.txt
	README.md
	TTbar_FullyLeptonic/TTbar_FullyLeptonic.cc

CMakeLists.txt

@@ -15,10 +15,20 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-cmake_minimum_required(VERSION 2.8)
+cmake_minimum_required(VERSION 3.4)
 project(MoMEMta-Tutorials)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -pedantic -Wextra")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -pedantic -Wextra")
+
+# Require c++11 *at least*, use default compiler standard if possible
+if (CMAKE_CXX_STANDARD_COMPUTED_DEFAULT STRLESS "11" OR
+        CMAKE_CXX_STANDARD_COMPUTED_DEFAULT STREQUAL "98")
+    set(CMAKE_CXX_STANDARD 11)
+    set(CMAKE_CXX_STANDARD_REQUIRED ON)
+endif()
+
+# Stick to the standard 
+set(CMAKE_CXX_EXTENSIONS OFF) 
 
 # Find dependencices
 

Paper_configs/MELA_Higgs_JP_sig.lua

@@ -0,0 +1,68 @@
+load_modules("MatrixElements/pp_x0_ZZ_4mu/build/libme_pp_x0_ZZ_4mu.so")
+
+lepton_plus_1 = declare_input("lepton_plus_1")
+lepton_plus_2 = declare_input("lepton_plus_2")
+lepton_minus_1 = declare_input("lepton_minus_1")
+lepton_minus_2 = declare_input("lepton_minus_2")
+
+parameters = {
+    energy = 13000.,
+    higgs_width = 6.3823e-3, -- default width, can be changed from analysis script
+    -- Path to param card to be changed from the analysis script;
+    -- There is only card prepared for each set of signal couplings.
+    param_card = "", 
+}
+
+inputs = {
+    lepton_plus_1.reco_p4,
+    lepton_minus_1.reco_p4,
+    lepton_plus_2.reco_p4,
+    lepton_minus_2.reco_p4,
+    }
+
+-- Build the partonic initial state
+BuildInitialState.boost = {
+    do_transverse_boost = true,
+    particles = inputs,
+}
+
+-- Call matrix element on fully defined partonic event
+MatrixElement.h_4f = {
+    pdf = "CT10nlo",
+    pdf_scale = 125.,
+    matrix_element = "pp_x0_ZZ_4mu_HC_UFO_P1_Sigma_HC_UFO_gg_mupmupmummum",
+    matrix_element_parameters = {
+        card = parameter("param_card"),
+    },
+    override_parameters = {
+        mdl_WX0 = parameter("higgs_width"),
+    },
+    initialState = "boost::partons",
+    -- Matrix element expects mu(+)mu(+)mu(-)mu(-),
+    -- whereas inputs are given as mu(+)mu(-)mu(+)mu(-),
+    -- hence indexing has to be reordered:
+    particles = {
+        inputs = inputs,
+        ids = {
+            {
+                pdg_id = -13,
+                me_index = 1,
+            },
+            {
+                pdg_id = 13,
+                me_index = 3,
+            },
+            {
+                pdg_id = -13,
+                me_index = 2,
+            },
+            {
+                pdg_id = 13,
+                me_index = 4,
+            },
+        }
+    },
+}
+
+-- Define quantity to be returned to MoMEMta
+integrand("h_4f::output")

Paper_configs/MELA_ZZ_bkg.lua

@@ -0,0 +1,58 @@
+load_modules("MatrixElements/pp_4mu/build/libpp_4mu.so")
+
+lepton_plus_1 = declare_input("lepton_plus_1")
+lepton_plus_2 = declare_input("lepton_plus_2")
+lepton_minus_1 = declare_input("lepton_minus_1")
+lepton_minus_2 = declare_input("lepton_minus_2")
+
+parameters = {
+    energy = 13000.,
+}
+
+inputs = {
+    lepton_plus_1.reco_p4,
+    lepton_minus_1.reco_p4,
+    lepton_plus_2.reco_p4,
+    lepton_minus_2.reco_p4,
+    }
+
+-- Build the partonic initial state
+BuildInitialState.boost = {
+    do_transverse_boost = true,
+    particles = inputs,
+}
+
+-- Call matrix element on fully defined partonic event
+MatrixElement.pp_4f = {
+    pdf = "CT10nlo",
+    pdf_scale = 91.,
+    matrix_element = "pp_4mu_sm_P1_Sigma_sm_uux_mupmummupmum",
+    matrix_element_parameters = {
+        card = "MatrixElements/pp_4mu/Cards/param_card.dat"
+    },
+    initialState = "boost::partons",
+    particles = {
+        inputs = inputs,
+        ids = {
+            {
+                pdg_id = -13,
+                me_index = 1,
+            },
+            {
+                pdg_id = 13,
+                me_index = 2,
+            },
+            {
+                pdg_id = -13,
+                me_index = 3,
+            },
+            {
+                pdg_id = 13,
+                me_index = 4,
+            },
+        }
+    },
+}
+
+-- Define quantity to be returned to MoMEMta
+integrand("pp_4f::output")

Paper_configs/SingleTop_TWminus.lua

@@ -0,0 +1,189 @@
+load_modules('MatrixElements/ME_gb_tw_mup_mum/libme_ME_gb_tw_mup_mum.so')
+
+cuba = {
+    -- Require 2% precision on weight
+    relative_accuracy = 0.02,
+    -- But limit number of evaluations if precision is not reached
+    max_eval = 500000,
+    n_start = 50000,
+}
+
+parameters = {
+    energy = 13000.,
+    top_mass = 173.,
+    top_width = 1.491500,
+    W_mass = 80.419002,
+    W_width = 2.047600,
+
+    ME_card = 'MatrixElements/ME_gb_tw_mup_mum/Cards/param_card.dat',
+}
+
+local lep_p = declare_input("lepton+")
+local lep_m = declare_input("lepton-")
+local bjet = declare_input("jet")
+
+-- Dummy massless 4-vector we'll use to integrate over its degrees of freedom
+local dummy_neutrino = declare_input("dummy_neutrino")
+
+-- Generate W and top invariant masses according to their Breit-Wigner
+BreitWignerGenerator.flatter_s12 = {
+    ps_point = add_dimension(),
+    mass = parameter('W_mass'),
+    width = parameter('W_width')
+}
+BreitWignerGenerator.flatter_s123 = {
+    ps_point = add_dimension(),
+    mass = parameter('top_mass'),
+    width = parameter('top_width')
+}
+BreitWignerGenerator.flatter_s45 = {
+    ps_point = add_dimension(),
+    mass = parameter('W_mass'),
+    width = parameter('W_width')
+}
+
+-- Extremely simplified transfer functions on lepton and b-jet energies
+GaussianTransferFunctionOnEnergy.tf_lep_m = {
+    ps_point = add_dimension(),
+    reco_particle = lep_m.reco_p4,
+    sigma = 0.02,
+    sigma_range = 5.,
+}
+lep_m.set_gen_p4("tf_lep_m::output")
+
+GaussianTransferFunctionOnEnergy.tf_lep_p = {
+    ps_point = add_dimension(),
+    reco_particle = lep_p.reco_p4,
+    sigma = 0.02,
+    sigma_range = 5.,
+}
+lep_p.set_gen_p4("tf_lep_p::output")
+
+GaussianTransferFunctionOnEnergy.tf_bjet = {
+    ps_point = add_dimension(),
+    reco_particle = bjet.reco_p4,
+    sigma = 0.15,
+    sigma_range = 5.,
+}
+bjet.set_gen_p4("tf_bjet::output")
+
+-- We integrate over the neutrino's Phi angle: use a "flat" transfer function (1 everywhere)
+-- The initial "reco" 4-vector is random, it only needs to have the right mass
+FlatTransferFunctionOnPhi.neutrino_phi_fixed = {
+    ps_point = add_dimension(),
+    reco_particle = dummy_neutrino.reco_p4,
+}
+
+-- Change of variables: align Top and W (from Top) propagators with grid!, i.e.:
+--      from the neutrinos's Theta angle and momentum,
+--      to the W and Top masses.
+-- The remaining d.o.f. (neutrino Phi) was taken care of by the "flat" T.F.
+-- The module "spits out" a completely specified neutrino 4-vector, using as input the generated Top and W masses
+SecondaryBlockB.sb_b = {
+    p1 = 'neutrino_phi_fixed::output',
+    p2 = lep_p.gen_p4,
+    p3 = bjet.gen_p4,
+    s12 = 'flatter_s12::s',
+    s123 = 'flatter_s123::s',
+}
+
+StandardPhaseSpace.phaseSpaceOut = {
+    particles = {lep_p.gen_p4, bjet.gen_p4, lep_m.gen_p4}
+}
+
+-- Iterate over solutions from secondary block
+Looper.looper_sb_b = {
+    solutions = "sb_b::solutions",
+    path = Path("met_substracted", "blockb", "looper_blockb")
+}
+
+    -- Define X = MET - neutrino (from W from Top)
+    -- We'll use X to fix the other neutrino's transverse momentum
+    VectorLinearCombinator.met_substracted = {
+        inputs = {'met::p4', 'looper_sb_b::particles/1'},
+        coefficients = {1., -1.}
+    }
+
+    -- "Main" change of variables:
+    -- Change 
+    --      from the other neutrino's d.o.f.
+    --      to the W's mass
+    -- Take care of momentum conservation in the process (i.e., also removes the 2 Bjorken-x variables).
+    -- The module "spits out" the other neutrino's completely specified 4-vector
+    BlockB.blockb = {
+        p2 = lep_m.gen_p4,
+        pT_is_met = true,
+        met = 'met_substracted::output',
+        s12 = 'flatter_s45::s',
+    }
+
+    -- Iterate over solutions from main block
+    Looper.looper_blockb = {
+        solutions = "blockb::solutions",
+        path = Path("boost", "twMinus", "integrand")
+    }
+
+        -- Final state parton event is now fully specified
+        full_inputs = {lep_p.gen_p4, bjet.gen_p4, lep_m.gen_p4, 'looper_sb_b::particles/1', 'looper_blockb::particles/1'}
+
+        -- Build initial state using momentum conservation
+        BuildInitialState.boost = {
+            do_transverse_boost = true,
+            particles = full_inputs
+        }
+
+        -- All terms that have to be multiplied together with the matrix element to define the integrand
+        jacobians = {'flatter_s12::jacobian', 'flatter_s123::jacobian', 'flatter_s45::jacobian', 'neutrino_phi_fixed::TF_times_jacobian', 'looper_blockb::jacobian', 'looper_sb_b::jacobian',
+                    'tf_lep_m::TF_times_jacobian', 'tf_lep_p::TF_times_jacobian', 'tf_bjet::TF_times_jacobian', 'phaseSpaceOut::phase_space'}
+
+        MatrixElement.twMinus = {
+          pdf = 'CT10nlo',
+          pdf_scale = parameter('top_mass'),
+
+          matrix_element = 'ME_gb_tw_mup_mum_sm_no_b_mass',
+          matrix_element_parameters = {
+              card =  parameter('ME_card'),
+          },
+
+          initialState = 'boost::partons',
+
+          particles = {
+            inputs = full_inputs,
+            ids = { 
+              {
+                pdg_id = -13,
+                me_index = 1,
+              },
+
+              {
+                pdg_id = 5,
+                me_index = 3,
+              },
+
+              {
+                pdg_id = 13,
+                me_index = 4,
+              },
+
+              {
+                pdg_id = 14,
+                me_index = 2,
+              },
+
+              {
+                pdg_id = -14,
+                me_index = 5,
+              }
+            }
+          },
+
+          jacobians = jacobians
+        }
+
+        -- Sum answer over the two nested loops over solution
+        DoubleLooperSummer.integrand = {
+            input = "twMinus::output"
+        }
+
+-- End of loops, define integrand to return to integrator
+integrand("integrand::sum")