From 5df6f5b1b2f6410daa8954fbfecff7726116d03c Mon Sep 17 00:00:00 2001
From: Ganrong Wang <ganrong.wang@cern.ch>
Date: Mon, 10 Jul 2023 23:12:09 +0200
Subject: [PATCH] Add B2JpsiX line

---
 .../StrippingB2CC/StrippingB2JpsiX.py         | 287 ++++++++++++++++++
 1 file changed, 287 insertions(+)
 create mode 100755 Phys/StrippingSelections/python/StrippingSelections/StrippingB2CC/StrippingB2JpsiX.py

diff --git a/Phys/StrippingSelections/python/StrippingSelections/StrippingB2CC/StrippingB2JpsiX.py b/Phys/StrippingSelections/python/StrippingSelections/StrippingB2CC/StrippingB2JpsiX.py
new file mode 100755
index 000000000..b89ef64ea
--- /dev/null
+++ b/Phys/StrippingSelections/python/StrippingSelections/StrippingB2CC/StrippingB2JpsiX.py
@@ -0,0 +1,287 @@
+###############################################################################
+# (c) Copyright 2000-2019 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.                                       #
+###############################################################################
+__author__  = 'Ganrong Wang'
+__date__    = '06/07/2023'
+
+__all__ = ( 'B2JpsiXConf', 'default_config' )
+
+"""
+(1) Bd,s->Jpsi X, X=eta, omega, eta_p or phi, with X->pi+ pi- pi0, pi0 including resolved pi0 and merged pi0
+"""
+from Gaudi.Configuration import *
+from GaudiConfUtils.ConfigurableGenerators import FilterDesktop, CombineParticles
+from PhysSelPython.Wrappers import Selection, DataOnDemand, MergedSelection, AutomaticData
+from StrippingConf.StrippingLine import StrippingLine
+from StrippingUtils.Utils import LineBuilder
+from GaudiKernel.SystemOfUnits import MeV
+
+default_config = {
+    'NAME'         : 'B2JpsiX',
+    'BUILDERTYPE'  : 'B2JpsiXConf',
+    'CONFIG'       :
+        {
+        'Pion_Trghp'                     : 0.5 
+        , 'Pion_minipchi2'                 : 6.0 
+        , 'Pion_pidk'                      : 5 
+        , 'Pi0_res_pt'                     : 300 
+        , 'Pi0_res_cl'                     : -1000
+        , 'Pi0_mer_pt'                     : 300 
+        , 'jpsi_mintree_pt'                : 500.0
+        , 'jpsi_mm_min'                    : 2996.916
+        , 'jpsi_mm_max'                    : 3196.916 
+        , 'jpsi_bpvdls_max'                : 3 
+        , 'jpsi_bpvdls_min'                : -3
+        , 'jpsi_mintree_pidmu'             : 0.0 
+        , 'X_apt'                          : 1000
+        , 'X_bpvvdz'                       : 0 
+        , 'X_vfaspf'                       : 9 
+        , 'X_bpvdira'                      : 0.95
+        , 'X_bpvvdchi2'                    : 25
+        , 'X_adamass'                      : 200 
+        , 'B_mass_min'                     : 4500
+        , 'B_mass_max'                     : 6500
+        , 'B_bpvipchi2'                    : 25
+        , 'B_bpvdira'                      : 0.9995
+        , 'B_bpvvdchi2'                    : 64
+        , 'B_mipchi2'                      : 9 
+        , 'B_vchi2vdof'                    : 10
+        , 'BMassWindow'                    : 1500
+        , 'Prescale'                       : 1.0 
+        , 'Postscale'                      : 1.0 
+        },
+    'WGs'     : [ 'B2CC' ],
+    'STREAMS' : [ 'Dimuon' ]
+    }
+
+
+
+class B2JpsiXConf(LineBuilder) :
+    """
+    Builder for R_X measurements
+    """
+
+    # now just define keys. Default values are fixed later
+    __configuration_keys__ = default_config['CONFIG'].keys()
+
+    def __init__(self, name, config):
+        LineBuilder.__init__(self, name, config)
+
+        self._name = name
+
+        name_B2Jpsieta   = name+"_eta2pipipi0"
+        name_B2Jpsiomega = name+"_omega2pipipi0"
+        name_B2Jpsietap  = name+"_etaprime2pipipi0"
+        name_B2Jpsiphi   = name+"_phi2pipipi0"
+
+        from StandardParticles import StdLoosePions as Pions
+        from StandardParticles import StdLoosePhi2KK as Phis
+        from StandardParticles import StdLooseJpsi2MuMu as JPsis
+
+        # 1 : Make the particles we will be actually using
+        SelPions  = self._filterHadron( name = "PionsFor" + self._name, sel = Pions, params = config )
+        SelPi0s   = self._myPi0s("PiZerosFor"+ self._name, params = config)
+
+        # 2 : Dileptons
+        SelJPsis  = self._filterJpsi( name = "JPsisFor" + self._name, sel = JPsis, params = config )
+        SelEta    = self._X2PiPiPi0( name="EtaFor"+self._name, Pions=SelPions, Pi0=SelPi0s, params = config )
+        SelOmega  = self._X2PiPiPi0( name="OmegasFor"+self._name, Pions=SelPions, Pi0=SelPi0s, params = config )
+        SelEtap   = self._X2PiPiPi0( name="EtapFor"+self._name, Pions=SelPions, Pi0=SelPi0s, params = config )
+        SelPhi    = self._X2PiPiPi0( name="PhiFor"+self._name, Pions=SelPions, Pi0=SelPi0s, params = config )
+
+
+        # 3 : Combine
+        SelB2Jpsieta   = self._makeB2LLX(name_B2Jpsieta  , XPart=SelEta,   lhcbJPsi=SelJPsis, params=config)
+        SelB2Jpsiomega = self._makeB2LLX(name_B2Jpsiomega, XPart=SelOmega, lhcbJPsi=SelJPsis, params=config)
+        SelB2Jpsietap  = self._makeB2LLX(name_B2Jpsietap,  XPart=SelEtap,  lhcbJPsi=SelJPsis, params=config)
+        SelB2Jpsiphi   = self._makeB2LLX(name_B2Jpsiphi,   XPart=SelPhi,   lhcbJPsi=SelJPsis, params=config)
+
+        # 4 : Declare Lines
+        SPDFilter = {
+            'Code'      : " ( recSummary(LHCb.RecSummary.nSPDhits,'Raw/Spd/Digits') < 600 )" ,
+            'Preambulo' : [ "from LoKiNumbers.decorators import *", "from LoKiCore.basic import LHCb" ]
+            }
+        self.B2JpsietaLine = StrippingLine(name_B2Jpsieta+"Line", 
+                                     prescale = config['Prescale'], 
+                                     postscale = 1, 
+                                     selection = SelB2Jpsieta, 
+                                     RelatedInfoTools = self._RelInfoTools(SelB2Jpsieta),
+                                     FILTER = SPDFilter, 
+                                     RequiredRawEvents = [], 
+                                     MDSTFlag = False )
+
+        self.B2JpsiomegaLine = StrippingLine(name_B2Jpsiomega+"Line", 
+                                     prescale = config['Prescale'], 
+                                     postscale = 1,  
+                                     selection = SelB2Jpsiomega, 
+                                     RelatedInfoTools = self._RelInfoTools(SelB2Jpsiomega),
+                                     FILTER = SPDFilter, 
+                                     RequiredRawEvents = [], 
+                                     MDSTFlag = False )
+
+        self.B2JpsietapLine = StrippingLine(name_B2Jpsietap+"Line", 
+                                     prescale = config['Prescale'], 
+                                     postscale = 1,  
+                                     selection = SelB2Jpsietap, 
+                                     RelatedInfoTools = self._RelInfoTools(SelB2Jpsietap),
+                                     FILTER = SPDFilter, 
+                                     RequiredRawEvents = [], 
+                                     MDSTFlag = False )
+
+        self.B2JpsiphiLine = StrippingLine(name_B2Jpsiphi+"Line", 
+                                     prescale = config['Prescale'], 
+                                     postscale = 1,  
+                                     selection = SelB2Jpsiphi, 
+                                     RelatedInfoTools = self._RelInfoTools(SelB2Jpsiphi),
+                                     FILTER = SPDFilter, 
+                                     RequiredRawEvents = [], 
+                                     MDSTFlag = False )
+
+
+        # 5 : register Line
+        self.registerLine( self.B2JpsietaLine   )
+        self.registerLine( self.B2JpsiomegaLine )
+        self.registerLine( self.B2JpsietapLine  )
+        self.registerLine( self.B2JpsiphiLine   )
+
+
+#####################################################
+    def _RelInfoTools(self, selection) :
+        """
+        Return related information for the given selection
+        """
+        # Use defaults where ever possible
+        _decay1 = "Bottom -> (J/psi(1S) -> ^[l+]CC [l-]CC) X"
+        _decay2 = "Bottom -> (J/psi(1S) -> [l+]CC ^[l-]CC) X"
+        RelInfo = [{ "Type":"RelInfoConeVariables", "Location":"ConeIsoInfo" },
+                   { "Type":"RelInfoVertexIsolation", "Location":"VtxIsoInfo" },
+                   { "Type":"RelInfoVertexIsolationBDT", "Location":"VtxIsoInfoBDT" },
+                   { "Type":"RelInfoBs2MuMuBIsolations", "Location":"BSMUMUVARIABLES"},
+                   # { "Type":"RelInfoBs2MuMuTrackIsolations", "Location":"BSMUMUTrackVARIABLES"},
+                   { "Type":"RelInfoBs2MuMuTrackIsolations", "DaughterLocations":{_decay1:"Muon1MuMuTrkIso", _decay2:"Muon2MuMuTrackIso"} },
+                   { "Type":"RelInfoMuonIsolation", "DaughterLocations":{_decay1:"Muon1Iso", _decay2:"Muon2Iso"} },
+                   ]
+        return RelInfo
+
+#####################################################
+    def _filterHadron( self, name, sel, params ):
+        """
+        Filter for all hadronic final states
+        """
+        # requires all basic particles to have IPCHI2 > KaonIPCHI2
+        # and hadron PT > KaonPT
+        # Added a ghost probability cut to reduce the runI rate
+        _Code = "(TRGHP < %(Pion_Trghp)s) &" \
+                "(MIPCHI2DV(PRIMARY) > %(Pion_minipchi2)s) & " \
+                "((HASRICH) & (PIDK <  %(Pion_pidk)s))" % params
+        # Actually implement the stuff
+        _Filter = FilterDesktop(Code = _Code)
+        return Selection( name, Algorithm = _Filter, RequiredSelections = [ sel ] )
+
+#####################################################
+
+    def _filterJpsi(self, name, sel, params):
+        """
+		Filter Jpsi from Jpsi2mumu
+		"""
+        _jpsi_cuts = " (MINTREE('mu+'==ABSID,PT) > %(jpsi_mintree_pt)s *MeV ) & (MM >  %(jpsi_mm_min)s) & (MM <  %(jpsi_mm_max)s) & ((BPVDLS> %(jpsi_bpvdls_max)s) | (BPVDLS< %(jpsi_bpvdls_min)s)) & (MINTREE('mu+'==ABSID,PIDmu) >  %(jpsi_mintree_pidmu)s)" % params
+        _filter_jpsi = FilterDesktop(Code= _jpsi_cuts)
+
+        return Selection( name, Algorithm = _filter_jpsi, RequiredSelections = [ sel ] )
+
+
+#####################################################
+    def _myPi0s(self, name, params):
+        """
+        Filter Pi0 from Std Pi0
+        """
+        from StandardParticles import StdLooseResolvedPi0 as _pi0resolved
+        from StandardParticles import StdLooseMergedPi0 as _pi0merged
+
+        _pi0_res_cuts = "(PT > %(Pi0_res_pt)s *MeV) & (CL > %(Pi0_res_cl)s)" % params
+        _pi0_mer_cuts = "(PT > %(Pi0_mer_pt)s *MeV)" % params
+
+        _filter_pi0resolved = FilterDesktop(Code=_pi0_res_cuts)
+        _filter_pi0merged = FilterDesktop(Code=_pi0_mer_cuts)
+        _selpi0resolved = Selection("Selection_"+name+"_pi0resolved", RequiredSelections=[_pi0resolved], Algorithm=_filter_pi0resolved)
+        _selpi0merged = Selection("Selection_"+name+"_pi0merged", RequiredSelections=[_pi0merged], Algorithm=_filter_pi0merged)
+        _sel = MergedSelection("Selection_"+name+"_pi0", RequiredSelections=[_selpi0resolved,_selpi0merged])
+        return _sel
+
+#####################################################
+    def _X2PiPiPi0( self, name, Pions, Pi0, params):
+        """
+        Make X -> pi+ pi- pi0 , X= eta,omega(782),eta_prime,phi(1020)
+        """
+        _X_Com_cuts = "(APT> %(X_apt)s*MeV) " % params #PreVertexCuts
+        _X_Mon_cuts = "(BPVVDZ>%(X_bpvvdz)s) & (VFASPF(VCHI2)<%(X_vfaspf)s) & (BPVDIRA>%(X_bpvdira)s) & (BPVVDCHI2>%(X_bpvvdchi2)s)" % params #postVertexCuts
+
+        _X2pipipizero = CombineParticles()
+
+        if (name.startswith("EtaFor")):
+            _X2pipipizero.DecayDescriptor = "eta -> pi+ pi- pi0"
+            _X_Com_cuts += "& (ADAMASS('eta') < %(X_adamass)s *MeV)" % params
+
+        if (name.startswith("OmegasFor")):
+            _X2pipipizero.DecayDescriptor = "omega(782) -> pi+ pi- pi0"
+            _X_Com_cuts += "& (ADAMASS('omega(782)') < %(X_adamass)s *MeV)" % params
+
+        if (name.startswith("EtapFor")):
+            _X2pipipizero.DecayDescriptor = "eta_prime -> pi+ pi- pi0"
+            _X_Com_cuts += "& (ADAMASS('eta_prime') < %(X_adamass)s *MeV)" % params
+
+        if (name.startswith("PhiFor")):
+            _X2pipipizero.DecayDescriptor = "phi(1020) -> pi+ pi- pi0"
+            _X_Com_cuts += "& (ADAMASS('phi(1020)') < %(X_adamass)s *MeV)" % params
+
+        _X2pipipizero.CombinationCut = _X_Com_cuts
+        _X2pipipizero.MotherCut = _X_Mon_cuts
+        _XConf = _X2pipipizero.configurable("Combine_"+name+"_PiPiPi0")
+        _selX2PIPIPIZERO = Selection( "Selection_"+name+"_X2pipipizero",
+                                          Algorithm = _XConf,
+                                          RequiredSelections = [ Pions, Pi0 ] )
+        return _selX2PIPIPIZERO
+
+
+#####################################################
+    def _makeB2LLX( self, name, XPart, lhcbJPsi, params):
+
+        """
+        CombineParticles / Selection for the B
+        """
+
+        #PreVertexCuts
+        B_com_cuts = "in_range( %(B_mass_min)s,AM, %(B_mass_max)s)" % params
+        #PostVertexCuts
+        B_mon_cuts = "(BPVIPCHI2() < %(B_bpvipchi2)s ) "\
+                 "& (BPVDIRA > %(B_bpvdira)s) "\
+                 "& (BPVVDCHI2 > %(B_bpvvdchi2)s) "\
+                 "& (MAXTREE(ISBASIC,MIPCHI2DV(PRIMARY)) > %(B_mipchi2)s )" \
+				 "& (VFASPF(VCHI2/VDOF) < %(B_vchi2vdof)s)"% params
+
+        if 'eta2pipipi0' in name:
+            _Decays =  ["B_s0 -> J/psi(1S) eta"]
+        if 'omega2pipipi0' in name:
+            _Decays =  ["B_s0 -> J/psi(1S) omega(782)"]
+        if 'etaprime2pipipi0' in name:
+            _Decays =  ["B_s0 -> J/psi(1S) eta_prime"]
+        if 'phi2pipipi0' in name:
+            _Decays =  ["B_s0 -> J/psi(1S) phi(1020)"]
+
+
+        _Combine = CombineParticles(DecayDescriptors = _Decays,
+                                    CombinationCut = B_com_cuts,
+                                    MotherCut = B_mon_cuts )
+
+        # Do it this way to prevent multiple candidates coming from the same J/psi candidate being in XPart and dilepton
+        selB2jpsiX = Selection("AllExceptJpsi"+name, Algorithm = _Combine, RequiredSelections = [ lhcbJPsi, XPart ])
+        return selB2jpsiX
+
+#####################################################
-- 
GitLab