Skip to content
Snippets Groups Projects
charmonium_to_dimuon.py 11.58 KiB
###############################################################################
# (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration           #
#                                                                             #
# This software is distributed under the terms of the GNU General Public      #
# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
#                                                                             #
# In applying this licence, CERN does not waive the privileges and immunities #
# granted to it by virtue of its status as an Intergovernmental Organization  #
# or submit itself to any jurisdiction.                                       #
###############################################################################
"""
Make dimuon lines for charmonium decays. Shared by B&Q and B2CC.
"""

import Functors as F
from GaudiKernel.SystemOfUnits import GeV, MeV
from Moore.config import Hlt2Line, register_line_builder
from PyConf import configurable
from RecoConf.algorithms_thor import ParticleCombiner, ParticleFilter
from RecoConf.standard_particles import make_ismuon_long_muon

from Hlt2Conf.lines.bandq.builders.prefilters import make_prefilters
from Hlt2Conf.lines.config_pid import nopid_muons

# old: mass window [MeV] around the Jpsi and Psi2s,
# old: applied to Jpsi and Psi2s lines [2018: 120 MeV]
# new: exactly define the mass upper, lower limit of Jpsi.
# new: Easier for users to define asymmetric mass window considering the FSR.
_JPSI_PDG_MASS_ = 3096.9 * MeV
_PSI2S_PDG_MASS_ = 3686.1 * MeV
_MASSWINDOW_LOW_JPSI_ = 150 * MeV
_MASSWINDOW_HIGH_JPSI_ = 150 * MeV
_MASSWINDOW_LOW_PSI2S_ = 150 * MeV
_MASSWINDOW_HIGH_PSI2S_ = 150 * MeV
_MASSMIN_JPSI = _JPSI_PDG_MASS_ - _MASSWINDOW_LOW_JPSI_
_MASSMAX_JPSI = _JPSI_PDG_MASS_ + _MASSWINDOW_HIGH_JPSI_
_MASSMIN_PSI2S = _PSI2S_PDG_MASS_ - _MASSWINDOW_LOW_PSI2S_
_MASSMAX_PSI2S = _PSI2S_PDG_MASS_ + _MASSWINDOW_HIGH_PSI2S_

_PIDMU_JPSI = -5.0
_PIDMU_PSI2S = -5.0
_PIDMU_UPSILON = -5.0

full_lines = {}
turbo_lines = {}


# define a new muon combiner, to make the cut on muons more configurable.
@configurable
def make_charmonium_muons(
    make_particles=make_ismuon_long_muon,
    name="charmonium_muons_{hash}",
    minPt_muon=None,
    minP_muon=None,
    minIPChi2_muon=None,
    minIP_muon=None,
    minPIDmu=None,
    maxIPChi2_muon=None,
):
    code = F.require_all(F.ISMUON)

    if minPt_muon is not None:
        code &= F.PT > minPt_muon
    if minP_muon is not None:
        code &= F.P > minP_muon
    if minIPChi2_muon is not None:
        code &= F.OWNPVIPCHI2 > minIPChi2_muon
    if minIP_muon is not None:
        code &= F.OWNPVIP > minIP_muon

    # ignoring minPIDmu cut here, thus ignoring it in ALL lines in this file
    if not nopid_muons() and minPIDmu is not None:
        code &= F.require_all(F.PID_MU > minPIDmu)

    if maxIPChi2_muon is not None:
        code &= F.OWNPVIPCHI2 < maxIPChi2_muon

    return ParticleFilter(make_particles(), name=name, Cut=F.FILTER(code))


# the new dimuon base for charmonium reconstruction
@configurable
def make_charmonium_dimuon_base(
    name="charmonium_dimuon_base_{hash}",
    DecayDescriptor="J/psi(1S) -> mu+ mu-",
    maxVertexChi2=25,
    maxDOCAChi2=None,
    minPt_muon=None,
    minP_muon=None,
    minIPChi2_muon=None,
    minIP_muon=None,
    minPIDmu=None,
    minPt_dimuon=None,
    minMass_dimuon=None,
    maxMass_dimuon=None,
):
    # get the long muons
    muons = make_charmonium_muons(
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minIPChi2_muon=minIPChi2_muon,
        minIP_muon=minIP_muon,
        minPIDmu=minPIDmu,
    )

    combination_code = F.ALL
    if minMass_dimuon is not None:
        combination_code &= F.MASS > minMass_dimuon
    if maxMass_dimuon is not None:
        combination_code &= F.MASS < maxMass_dimuon
    if maxDOCAChi2 is not None:
        combination_code &= F.SDOCACHI2(1, 2) < maxDOCAChi2

    # require that the muons come from the same vertex
    vertex_code = F.require_all(F.CHI2DOF < maxVertexChi2)
    if minPt_dimuon is not None:
        vertex_code &= F.PT > minPt_dimuon

    return ParticleCombiner(
        name=name,
        Inputs=[muons, muons],
        DecayDescriptor=DecayDescriptor,
        CombinationCut=combination_code,
        CompositeCut=vertex_code,
    )


# @configurable
# def make_charmonium_samesign_dimuon_base(name='charmonium_samesign_dimuon_base_{hash}'):
#    return make_charmonium_dimuon_base(name=name, DecayDescriptor='J/psi(1S) -> mu+ mu+')


@configurable
def make_charmonium_samesign_dimuon_base(
    name="charmonium_samesign_dimuon_base_{hash}",
    DecayDescriptor="[J/psi(1S) -> mu+ mu+]cc",
    maxVertexChi2=25,
    maxDOCAChi2=None,
    minPt_muon=None,
    minP_muon=None,
    minIPChi2_muon=None,
    minIP_muon=None,
    minPIDmu=None,
    minPt_dimuon=None,
    minMass_dimuon=None,
    maxMass_dimuon=None,
):
    # get the long muons
    muons = make_charmonium_muons(
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minIPChi2_muon=minIPChi2_muon,
        minIP_muon=minIP_muon,
        minPIDmu=minPIDmu,
    )

    combination_code = F.ALL
    if minMass_dimuon is not None:
        combination_code &= F.MASS > minMass_dimuon
    if maxMass_dimuon is not None:
        combination_code &= F.MASS < maxMass_dimuon
    if maxDOCAChi2 is not None:
        combination_code &= F.SDOCACHI2(1, 2) < maxDOCAChi2

    # require that the muons come from the same vertex
    vertex_code = F.require_all(F.CHI2DOF < maxVertexChi2)
    if minPt_dimuon is not None:
        vertex_code &= F.PT > minPt_dimuon

    return ParticleCombiner(
        name=name,
        Inputs=[muons, muons],
        DecayDescriptor=DecayDescriptor,
        CombinationCut=combination_code,
        CompositeCut=vertex_code,
    )


@configurable
def make_charmonium_dimuon(
    name="charmonium_dimuon_{hash}",
    DecayDescriptor="J/psi(1S) -> mu+ mu-",
    minPt_dimuon=None,
    minPt_muon=300 * MeV,
    minP_muon=None,
    minIPChi2_muon=None,
    maxVertexChi2=25,
    minPIDmu=-5,
    ownpvdls_min=None,
    maxDOCAChi2=None,
    minMass_dimuon=None,
    maxMass_dimuon=None,
):
    make_particles = make_charmonium_dimuon_base(
        DecayDescriptor=DecayDescriptor,
        maxVertexChi2=maxVertexChi2,
        minIPChi2_muon=minIPChi2_muon,
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minPIDmu=minPIDmu,
        maxDOCAChi2=maxDOCAChi2,
        minMass_dimuon=minMass_dimuon,
        maxMass_dimuon=maxMass_dimuon,
    )

    code = F.require_all(F.CHI2DOF < maxVertexChi2)
    if minPt_dimuon is not None:
        code &= F.PT > minPt_dimuon
    if ownpvdls_min is not None:
        code &= F.OWNPVDLS > ownpvdls_min

    return ParticleFilter(make_particles, name=name, Cut=F.FILTER(code))


@configurable
def make_charmonium_samesign_dimuon(
    name="charmonium_samesign_dimuon_{hash}",
    DecayDescriptor="[J/psi(1S) -> mu+ mu+]cc",
    minPt_dimuon=None,
    minPt_muon=300 * MeV,
    minP_muon=None,
    maxVertexChi2=25,
    minIPChi2_muon=None,
    minPIDmu=-1,
    ownpvdls_min=None,
    maxDOCAChi2=None,
    minMass_dimuon=None,
    maxMass_dimuon=None,
):
    make_particles = make_charmonium_samesign_dimuon_base(
        DecayDescriptor=DecayDescriptor,
        maxVertexChi2=maxVertexChi2,
        minIPChi2_muon=minIPChi2_muon,
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minPIDmu=minPIDmu,
        maxDOCAChi2=maxDOCAChi2,
        minMass_dimuon=minMass_dimuon,
        maxMass_dimuon=maxMass_dimuon,
    )

    code = F.require_all(F.CHI2DOF < maxVertexChi2)
    if minPt_dimuon is not None:
        code &= F.PT > minPt_dimuon
    if ownpvdls_min is not None:
        code &= F.OWNPVDLS > ownpvdls_min

    return ParticleFilter(make_particles, name=name, Cut=F.FILTER(code))


@configurable
def make_jpsi(
    name="jpsi_{hash}",
    minMass_dimuon=_MASSMIN_JPSI,
    maxMass_dimuon=_MASSMAX_JPSI,
    minPt_muon=300 * MeV,
    minP_muon=None,
    minPt_Jpsi=None,
    maxVertexChi2=25,
    minPIDmu=_PIDMU_JPSI,
):
    code = F.ALL
    if minPt_Jpsi is not None:
        code &= F.PT > minPt_Jpsi

    dimuon = make_charmonium_dimuon(
        DecayDescriptor="J/psi(1S) -> mu+ mu-",
        minPt_dimuon=minPt_Jpsi,
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minPIDmu=minPIDmu,
        minMass_dimuon=minMass_dimuon,
        maxMass_dimuon=maxMass_dimuon,
        maxVertexChi2=maxVertexChi2,
    )

    return ParticleFilter(dimuon, name=name, Cut=F.FILTER(code))


@configurable
def make_psi2s(
    name="psi2s_{hash}",
    minMass_dimuon=_MASSMIN_PSI2S,
    maxMass_dimuon=_MASSMAX_PSI2S,
    minPt_muon=300 * MeV,
    minP_muon=None,
    minPt_Psi2S=None,
    maxPt_Psi2S=None,
    maxVertexChi2=25,
    minPIDmu=_PIDMU_PSI2S,
):
    code = F.ALL
    if minPt_Psi2S is not None:
        code &= F.PT > minPt_Psi2S
    if maxPt_Psi2S is not None:
        code &= F.PT < maxPt_Psi2S

    dimuon = make_charmonium_dimuon(
        DecayDescriptor="psi(2S) -> mu+ mu-",
        minPt_dimuon=minPt_Psi2S,
        minPt_muon=minPt_muon,
        minP_muon=minP_muon,
        minPIDmu=minPIDmu,
        minMass_dimuon=minMass_dimuon,
        maxMass_dimuon=maxMass_dimuon,
        maxVertexChi2=maxVertexChi2,
    )

    return ParticleFilter(dimuon, name=name, Cut=F.FILTER(code))


@configurable
def make_jpsi_tight(name="jpsi_tight_{hash}"):
    return make_jpsi(minPt_Jpsi=3000 * MeV, minPt_muon=650 * MeV, minP_muon=10 * GeV)


@configurable
def make_psi2s_tight(name="psi2s_tight_{hash}"):
    return make_psi2s(minPt_Psi2S=3000 * MeV, minPt_muon=650 * MeV, minP_muon=10 * GeV)


#############
# define lines
##############
@register_line_builder(turbo_lines)
@configurable
def JpsiToMuMu_line(
    name="Hlt2_JpsiToMuMu", prescale=0.3
):  # 230k Jpsi per pb-1 in 2024. Still 70k Jpsi per pb-1 after prescaling.
    line_alg = make_jpsi(
        minPt_muon=650 * MeV,
        minPIDmu=-2.0,
        minMass_dimuon=_MASSMIN_JPSI - 60,
        maxMass_dimuon=_MASSMAX_JPSI + 60,
    )
    return Hlt2Line(name=name, algs=make_prefilters() + [line_alg], prescale=prescale)


@register_line_builder(turbo_lines)
@configurable
def JpsiToMuMu_highpT_line(
    name="Hlt2_JpsiToMuMu_highpT",
    stream="bandq",
):  # for keeping the statistic for high-pT Jpsi samples
    line_alg = make_jpsi(
        minPt_muon=700 * MeV,
        minPIDmu=-2.0,
        minPt_Jpsi=5000 * MeV,  # A typical threshold above which NRQCD works.
        minP_muon=3000 * MeV,
        minMass_dimuon=_MASSMIN_JPSI,
        maxMass_dimuon=_MASSMAX_JPSI,
    )
    return Hlt2Line(name=name, algs=make_prefilters() + [line_alg], stream=stream)


@register_line_builder(turbo_lines)
@configurable
def Psi2SToMuMu_line(name="Hlt2_Psi2SToMuMu", prescale=1):
    line_alg = make_psi2s(
        minPt_muon=650 * MeV,
        minPIDmu=0.0,
        minMass_dimuon=_MASSMIN_PSI2S - 60,
        maxMass_dimuon=_MASSMAX_PSI2S + 60,
    )
    return Hlt2Line(name=name, algs=make_prefilters() + [line_alg], prescale=prescale)


@register_line_builder(full_lines)
@configurable
def JpsiToMuMuTight_line(name="Hlt2_DiMuonJPsiTightFull", prescale=1, persistreco=True):
    line_alg = make_jpsi_tight()
    return Hlt2Line(
        name=name,
        algs=make_prefilters() + [line_alg],
        prescale=prescale,
        persistreco=persistreco,
    )


@register_line_builder(full_lines)
@configurable
def Psi2SToMuMuTight_line(
    name="Hlt2_DiMuonPsi2STightFull", prescale=1, persistreco=True
):
    line_alg = make_psi2s_tight()
    return Hlt2Line(
        name=name,
        algs=make_prefilters() + [line_alg],
        prescale=prescale,
        persistreco=persistreco,
    )


all_lines = {}
all_lines.update(turbo_lines)
all_lines.update(full_lines)