diff --git a/CMakeLists.txt b/CMakeLists.txt index bba2d728e4539edaad441564cf4a22a176c8ff75..05b5d7acee56be31942497fd93f851a4ba63cbd6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ ############################################################################### -# (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration # +# (c) Copyright 2000-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". # @@ -62,7 +62,6 @@ lhcb_add_subdirectories( Phys/LoKiAlgo Phys/LoKiArrayFunctors Phys/LoKiFitters - Phys/LoKiJets Phys/LoKiPhys Phys/LoKiProtoParticles Phys/LoKiTracks diff --git a/Phys/LoKiJets/CMakeLists.txt b/Phys/LoKiJets/CMakeLists.txt deleted file mode 100644 index 3f9fb4b783847348502643e30fbd241992c89e39..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/CMakeLists.txt +++ /dev/null @@ -1,48 +0,0 @@ -############################################################################### -# (c) Copyright 2000-2021 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. # -############################################################################### -#[=======================================================================[.rst: -Phys/LoKiJets -------------- -#]=======================================================================] - -gaudi_add_module(LoKiJets - SOURCES - src/LoKiFastJetMaker.cpp - src/LoKiFastJetWithAreaMaker.cpp - src/LoKiJetMakerAlg.cpp - src/LoKiJetMakerWR2VtxAlg.cpp - src/LoKiJetParticleMaker.cpp - src/LoKiJets2JetsAlg.cpp - src/LoKiSeedConeJetMaker.cpp - src/LoKiSeedFinder.cpp - src/LoKiVVSeedFinder.cpp - LINK - AIDA::aida - Boost::headers - FastJet::FastJet - Gaudi::GaudiAlgLib - Gaudi::GaudiKernel - Gaudi::GaudiUtilsLib - LHCb::LHCbMathLib - LHCb::LoKiCoreLib - LHCb::PhysEvent - LHCb::RecEvent - Phys::DaVinciInterfacesLib - Phys::DaVinciKernelLib - Phys::DaVinciTypesLib - Phys::LoKiAlgo - Phys::LoKiPhysLib - ROOT::GenVector - ROOT::Hist - ROOT::MathCore -) - -gaudi_install(PYTHON) diff --git a/Phys/LoKiJets/doc/release.notes b/Phys/LoKiJets/doc/release.notes deleted file mode 100755 index 888ed48bdd2ff8c95c9bb32f8c3d3b7a48f3d8c1..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/doc/release.notes +++ /dev/null @@ -1,444 +0,0 @@ -! ----------------------------------------------------------------------------- -! Package : Phys/LoKiJets -! Responsible : Victor Coco victor.coco@cern.ch -! Purpose : Jet finding algorithm embedded into LHCb software -! ----------------------------------------------------------------------------- - -!========================= LoKiJets v6r15 2016-03-04 ========================== - -! 2016-02-14 - Phil Ilten - - Added option to turn off banner for HLT. - -!========================= LoKiJets v6r14 2015-11-25 ========================== - -! 2015-11-01 - Gerhard Raven - - replace endreq with endmsg - -!========================= LoKiJets v6r13 2015-10-07 ========================== -! 2015-08-12 - Gerhard Raven - - remove #include of obsolete Gaudi headers - -!========================= LoKiJets v6r12 2014-10-30 ========================= - -! 2014-10-20 - Chris Jones - - General downgrading of various Error and Warning message, to sanitise a bit - things for running in production, the Stripping. - -!========================= LoKiJets v6r11p1 2014-09-30 ========================= - -! 2014-09-18 - Marco Clemencic - - Fixed compilation problem with head version of ROOT. - -!========================= LoKiJets v6r11 2014-07-25 ========================= - -! 2014-06-18 - Wouter Hulsbergen - - Change default values of RecombinationScheme and Strategy in LokiFastJetmaker - to default values in fastjet. The change in Strategy saves a factor ~50 in - execution time. - -!========================= LoKiJets v6r10 2014-02-20 ========================= - -! 2014-02-06 - Marco Clemencic - - Reverted the workaround because it's now in Gaudi. - -! 2014-02-03 - Marco Clemencic - - Fixed problems with some incompatibilities between the new Plugin Service and - the old one. - Some templated tools are called using explicitly the old naming, and that was - not working with the new plugin service. - -! 2014-01-25 - Chris Jones - - Sanitise unneccessary info() printout for production. - - Replace cout (not allowed) with debug(). - - Protect debug() and verbose() messages behind if statements. - -!========================= LoKiJets v6r9 2013-10-29 ========================= - -! 2013-10-10 - Chris Jones - - Some minor clean ups and typo fixes. - -! 2013-10-10 - Cedric Potterat - - Fixing Memory Leak detected by Valgrind. - -!========================= LoKiJets v6r8p2 2013-06-13 ========================= - -! 2013-05-21 - Chris Jones - - Do not use cout (not allowed), replace with info() messages. - - Protect debug() with if ( msglevel(MSG::DEBUG) ) checks - - Fix clang 3.2 warnings. - -!========================= LoKiJets v6r8p1 2013-04-05 ========================= - -! 2013-03-22 - Cedric Potterat - - Fix Covertiy defect in LoKiSeedFinder.cpp - -!========================= LoKiJets v6r8 2012-11-30 ========================= - -! 2012-11-29 - Marco Clemencic - - Added CMake configuration file. - -! 2012-11-26 - Victor Coco - - LoKiFastJetWithAreaMaker now compatible with fastjet 3.0. - TODO: implement more easy config of the area type - -!========================= LoKiJets v6r7 2012-10-04 ========================= - -! 2012-09-27 - Murilo Rangel -- JetID cuts Updated also for no PV association - -! 2012-09-26 - Murilo Rangel -- Updated jec.root file and JetID cuts - -!========================= LoKiJets v6r6p1 2012-07-09 ========================= - -! 2012-06-29 - Chris Jones - - Use a public instance of LoKi::DistanceCalculator - -!========================= LoKiJets v6r6 2012-06-28 ========================= - -! 2012-06-27 - potterat - - correct the the jet 2 pv association. - -!========================= LoKiJets v6r5 2012-05-03 ========================= - -! 2012-04-27 - rangel -- Added more jet information - -!========================= LoKiJets v6r4 2012-03-01 ========================= - -! 2012-02-29 - rangel - - n90 small bug fix - -!========================= LoKiJets v6r3 2012-02-03 ========================= - -! 2012-01-28 - Chris Jones - - Fix compilation warnings introduced with the last commit - - Remove obsolete file LoKiJets_dll.cpp - -! 2012-01-27 - Victor Coco - - Add some jet ID variable in LoKiJetMakerAlg - - Clean LoKiFastJetMaker from outdated options - - Add JEC reference histogram in histo/jec.root - - Add JEC correction in LoKiJetMakerAlg - -! 2012-01-17 - Victor Coco - - Add some protection against jets made of daughter with no PV info + update - refernece point and endVertex of the jets that have associated PV - -!========================= LoKiJets v6r2 2011-12-15 ========================= - -! 2011-11-28 - Murilo Rangel - - Add first JetID variables in JetMakerAlg. - -! 2011-11-25 - Victor Coco - - Add structure for appending jet ID info in JetMakerAlg. To be continued... - -! 2011-11-23 - Chris Jones - - Fix various gcc 4.6 warnings - -!================== LoKiJets v6r1 2011-11-10 =============================== - -! 2011-11-07 - Cedric Potterat - - Setting the option to set the referencePoint of the Jet to the associated - vertex (if true) in LoKiJetMakerAlg - -! 2011-10-27 - Victor Coco - - Rethinking of the jet to vertex association. Modification in LoKiJetMakerAlg - -!================== LoKiJets v6r0 2011-04-29 =============================== -LoKiJets add new methods with multi PV, move codes in JetAccessories(MC) -and correction in the cmt/requirements of LoKiJets - -! 2011-04-19 - Cédric Potterat - - new LoKiJetMakerWR2VtxAlg.cpp - . new algorith to run the JetMaker algorithms on a RecVertx container - - new version of src/LoKiSeedConeJetMaker.{cpp,h}, src/LoKiSeedFinder.{cpp,h}, - src/LoKiVVSeedFinder.{cpp,h}, LoKiFastJetMaker.{cpp,h}, LoKiFastJetWithAreaMaker.cpp : - . update for the new IJetMaker.h fonction (with vertex) - . new algorithm for the SeedFinders - . add the possiblility to assign a PV to a FastJet (via option) - - move LoKiMCPartBHadronSeed.{cpp,h}, PartonicJets2HepMCJets.cpp, HepMCJets2Jets.cpp, - HepMCJets2HepMCJets.cpp, PartonicJets2HepMCJets.cpp to JetAccessoriesMC - - src/LoKiJetParticleMaker - . increase the reserver of the output vector from 20 to 100 - - cmt/requirements - . Increase to v6r0 - - -!================== LoKiJets v5r5 2011-04-21 =============================== - -! 2011-04-08 - Juan Palacios - - src/HepMCJets2HepMCJets.cpp - - src/Jets2Jets.cpp - - src/HepMCJets2Jets.cpp - - src/PartonicJets2HepMCJets.cpp - . Fix memory leaks. - - cmt/requirements - . Increase to v5r5. - -!================== LoKiJets v5r4 2011-02-20 =============================== - -! 2011-02-08 - Juan Palacios - - cmt/requirements - . Use Phys/LoKiAlgoMC. - -! 2011-02-02 - Cédric Potterat - - add src/LoKiSeedConeJetMaker - creat Jet with a cone algo around a seed - (seed created by a call to a LoKiJetMaker Algo) - - add src/LoKiSeedFinder - creat Seed for LoKiSeedConeJetMaker (Barcelona and Milano Algo) - can be call as a LoKiJetMaker alog, output Seed - - add src/LoKiVVSeedFinder - creat Seed for LoKiSeedConeJetMaker (Lausanne Algo) - can be call as a LoKiJetMaker alog, output Seed - - add src/LoKiMCPartBHadronSeed - creat Seed for LoKiSeedConeJetMaker (MCPart B Hadron is the Seed) - can be call as a LoKiJetMaker alog, output Seed - - Increase to v5r4. - -!================== LoKiJets v5r3 2011-02-04 =============================== - -! 2011-02-04 - Juan Palacios - - Add needed LoKi MC packages. - - Increase to v5r3. - -!================== LoKiJets v5r2p2 2010-07-12 =============================== - -! 2010-07-12 - Juan Palacios - - cmt/requirements - . Remove Phys/Particle2MCTruth. Not actually needed. - - src/HepMCJets2Jets.cpp - . Remove include of IP2MCP.h - -!================== LoKiJets v5r2p1 2010-07-12 =============================== - -! 2010-07-12 - Juan Palacios - - cmt/requirements - . Need Phys/Particle2MCTruth. Increase to v5r2p1. - -!================== LoKiJets v5r2 2010-04-30 ================================= - -! 2010-04-16 - Juan Palacios - - src/LoKiJetParticleMaker.cpp, src/LoKiJets2JetsAlg.cpp - . use Particle::Range instead of Particle::Container for inputs. - - cmt/requirements - . Increase to v5r2 - -!================== LoKiJets v5r1 2010-02-22 ================================= - -! 2010-02-04 - Victor Coco - - rename Configuration.py to JetMaker_Config.py in order to fix some warning - concerning configurables - -!================== LoKiJets v5r0 2010-02-02 ================================= - -! 2010-01-11 - Victor Coco - - remove stdint.h include from LoKiFastJetMaker (not recognised under windows) - -! 2010-01-08 - Victor Coco - - remove call to CommonParticles.Utils (not needed and unknown in windows - configuration --> package not build) - -! 2009-12-14 - Victor Coco - - add some python configuration for jet making + some simple option file. - Not definitive, just as starting point - - add all the needed classes to relate jets at the different level of simulation. - In the end they might go in a LoKiJetsMCTools package (except the Jets2Jets tool) - - Clean a bit, remove warning, remove completly Kt jet maker - - Update few things like the nominal way of calling StdJets on demand - - - -!================== LoKiJets v4r1 - -! 2008-11-09 - Vanya BELYAEV - - remove options directory with all opts-files - - add python directory with corresponding on-demand configuration of - standard Kt-Jets on-demand. - - to enable on-demand cratino of jets: - - import LoKiJets.StdKtJets_Configuration - - The location is defined to be "Phys/StdKtJets" - - The jets are made from: - - "Phys/Std/NoPIDPions" - and - "Phys/StdLooseAllPhotons" - -! 2008-06-26 - Vanya BELYAEV - - - LoKi::JetMaker - use "LoKi::FastJetMaker" as the default jet maker - - - KtJetMaker is declaread as obsolete - - - cmt/requirements - - version increment to v4r1 - -!========================== LoKiJets v4r0 2007-12-18 ===================== - -! 2007-12-18 - Patrick Koppenburg - - Small patch by Victor Coco - -! 2007-12-10 - Vanya BELYAEV - - Add fastjet package, but still KtJet as the default engine - new files: - - src/LoKiFastJetMaker.cpp - src/LoKiFastJetMaker.h - src/LoKiFastJetWithAreaMaker.cpp - - After successfull tests, the "Kt"-stuff will be removed - - - -!========================== LoKiJets v3r0 2007-12-10 ===================== - -! 2007-12-10 - Vanya BELYAEV - - - use ktjet from LCG_Interfaces - no local ktejet is used anymore - - postpone a bit the fastjet tools - -! 2007-10-16 - Vanya BELYAEV - - add new tool LoKi::KtJetyMaker - the simplest implementation - of IJetMaker interface using KtJet algorithm - - redesigne the algorithm to invoke th eJEt-Making tool - - add LoKi::JetParticleMaker - the implementation of "particle maker" - - cmt/requirements - increment the version to v3r0 - - -!========================== LoKiJets v2r4 2007-09-24 ===================== - -! 2007-08-22 - Patrick Koppenburg - - LoKiKtJetMaker.cpp - Replace all by particle -> setParticleID( m_Id ); by - particle -> setParticleID( LHCb::ParticleID (m_Id) ); - - Increment to v2r4 - -! 2007-07-17 - Patrick Koppenburg - - Removed _load file. Rlease when needed. - -!========================== LoKiJets v2r3p1 2007-04-19 =================== -! 2007-04-19 - P. Koppenburg - Remove controversial doxygen comment - -!========================== LoKiJets v2r3 2007-03-26 =================== - -! 2007-03-27 - Vanya BELYAEV - - src/LoKiKtJetMaker.cpp - use pointer instead of copy... - [somethng is definitely wrong inside KtJet and/or KtLorentzVector] - - cmt/requirements - version increment to v2r3 - -!========================== LoKiJets v2r2 2007-03-04 =================== - -! 2007-03-04 - Vanya BELYAEV - - - src/LoKiKtJetMaker.cpp - improve the selection of the associated particle - - cmt/requirements - version increment to v2r2 - -!========================== LoKiJets v2r1 2006-11-12 =================== -! 2006-11-12 - Vanya BELYAEV - - src/LoKiKtJetMaker.cpp - - many fixes, in particular, fix for the particle combiner, - better treatment of orting schemes, better error handling, - guard against empty input lists, etc... - - Finally it works more or less OK. - -- options/StdKtJets.opts, - options/StdMCKtJets.opts - options/StdHepKtJets.opts - - configuration files to run the corresponding jet-finding algorithms: - - - for regular Kt jets - $LOKIJETSOPTS/StdKtJets.opts - - - for Kt jets from MC-particles: - $LOKIJETSOPTS/StdMCKtJets.opts - - - for Kt jets from HepMC-particles: - $LOKIJETSOPTS/StdHepMCKtJets.opts - -- cmt/requirements - - version increment to v2r1 - - -!========================== LoKiJets v2r0 2006-10-18 =================== -! 2006-10-18 - Vanya BELYAEV - - - src/LoKiKtJetMaker.cpp - fix a bug with the creation of jets for particle ("isAtomic"-flag) - - (re)-tag as v2r0 ! - - -!========================== LoKiJets v2r0 2006-09-06 =================== - -!2006-09-06 - Vanya BELYAEV - - reimport the pckage under Phys has and adapt for DC06 software branch - - tool is temporary(?) removed - - cmt/requiremens - version change to v2r0 - - -! 2005-07-26 - Vanya BELYAEV - - src/KtJetTool.cpp - new tool which implements IJetMaker interface - - -!========================== LoKiJets v1r1 2005-06-13 =================== -! 2005-06-02 - P. Koppenburg - KtJets moved to External/ - -!========================== LoKiJets v1r0 2005-02-21 ========================== -! 2005-03-29 - Vanya BELYAEV - - - no new code, just SLOCCount report - -Categorizing files. -Computing results. - -SLOC Directory SLOC-by-Language (Sorted) -150 src_top_dir cpp=150 -36 cmt csh=20,sh=16 -0 CVS (none) -0 doc (none) -0 src_CVS (none) - -Totals grouped by language (dominant language first): -cpp: 150 (80.65%) -csh: 20 (10.75%) -sh: 16 (8.60%) - -Total Physical Source Lines of Code (SLOC) = 186 -Development Effort Estimate, Person-Years (Person-Months) = 0.03 (0.41) - (Basic COCOMO model, Person-Months = 2.4 * (KSLOC**1.05)) -Schedule Estimate, Years (Months) = 0.15 (1.78) - (Basic COCOMO model, Months = 2.5 * (person-months**0.38)) -Estimated Average Number of Developers (Effort/Schedule) = 0.23 -Total Estimated Cost to Develop = $ 4,620 - (average salary = $56,286/year, overhead = 2.40). -SLOCCount is Open Source Software/Free Software, licensed under the FSF GPL. -Please credit this data as "generated using 'SLOCCount' by David A. Wheeler." - - -! 2005-03-21 - Vanya BELYAEV - - New package - -! ----------------------------------------------------------------------------- -! The END -! ----------------------------------------------------------------------------- diff --git a/Phys/LoKiJets/histo/jec.root b/Phys/LoKiJets/histo/jec.root deleted file mode 100644 index 8d88e3fa9a125e663c71987be05e43aa773ca85c..0000000000000000000000000000000000000000 Binary files a/Phys/LoKiJets/histo/jec.root and /dev/null differ diff --git a/Phys/LoKiJets/python/LoKiJets/StdJets_Configuration.py b/Phys/LoKiJets/python/LoKiJets/StdJets_Configuration.py deleted file mode 100755 index ec0ac88bb648195f07919f771787009547b5e4f5..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/python/LoKiJets/StdJets_Configuration.py +++ /dev/null @@ -1,53 +0,0 @@ -#!/usr/bin/env python -############################################################################### -# (c) Copyright 2000-2018 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. # -############################################################################### -# ============================================================================= -# ============================================================================= -## @file LoKiJets/StdJets.py -# configuration file for 'Standard Jets made with FastJet' -# @author Victor COCO victor.coco@cern.ch -# @date 2009-09-09 -# ============================================================================= -""" -Configuration file for jets made with FastJet -""" -from __future__ import print_function -__author__ = "Victor COCO victor.coco@cern.ch" -# ============================================================================= -__all__ = ('StdJets', 'FastJet', 'locations') -# ============================================================================= -from Gaudi.Configuration import * -from Configurables import LoKi__JetMaker, LoKi__FastJetMaker -from CommonParticles.Utils import * - -## create the algorithm -StdJets = LoKi__JetMaker('StdJets') - -# select the jet maker -StdJets.JetMaker = 'LoKi__FastJetMaker' - -StdJets.addTool(LoKi__FastJetMaker) - -FastJet = getattr(StdJets, 'LoKi__FastJetMaker') - -## configure Data-On-Demand service -locations = updateDoD(StdJets) - -## ============================================================================ -if '__main__' == __name__: - - print(__doc__) - print(__author__) - print(locationsDoD(locations)) - -# ============================================================================= -# The END -# ============================================================================= diff --git a/Phys/LoKiJets/python/LoKiJets/__init__.py b/Phys/LoKiJets/python/LoKiJets/__init__.py deleted file mode 100755 index 0f1935074d33ef3cdecabacc28e80aa77b669a5f..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/python/LoKiJets/__init__.py +++ /dev/null @@ -1,44 +0,0 @@ -#!/usr/bin/env python -############################################################################### -# (c) Copyright 2000-2018 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. # -############################################################################### -# ============================================================================= -# @file -# Mandatory file to keep Python 'package' for Phys/LoKiJets -# -# This file is a part of LoKi project - -# "C++ ToolKit for Smart and Friendly Physics Analysis" -# -# The package has been designed with the kind help from -# Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, -# contributions and advices from G.Raven, J.van Tilburg, -# A.Golutvin, P.Koppenburg have been used in the design. -# -# @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl -# @date 2008-11-09 -# ============================================================================= -""" -Mandatory file to keep Python 'package' for Phys/LoKiJets - -This file is a part of LoKi project - -\"C++ ToolKit for Smart and Friendly Physics Analysis\" - -The package has been designed with the kind help from -Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, -contributions and advices from G.Raven, J.van Tilburg, -A.Golutvin, P.Koppenburg have been used in the design. -""" -# ============================================================================= -__author__ = " Vanya BELYAEV Ivan.Belyaev@nikhef.nl " -# ============================================================================= - -# ============================================================================= -# The END -# ============================================================================= diff --git a/Phys/LoKiJets/src/LoKiFastJetMaker.cpp b/Phys/LoKiJets/src/LoKiFastJetMaker.cpp deleted file mode 100644 index a0b752d7771651408af0524650aa13511c54f171..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiFastJetMaker.cpp +++ /dev/null @@ -1,198 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "LoKi/ILoKiSvc.h" - -#include "LoKiFastJetMaker.h" - -#include "fastjet/ClusterSequence.hh" -#include "fastjet/config.h" - -#include "DetDesc/IDetectorElement.h" - -/** - * Implementation file for class LoKi::FastJetMaker - * @author Vanya BELYAEV ibelyaev@physics.syr.edu - * @author Victor COCO coco@lapp.in2p3.fr - * @date 2007-10-22 - */ - -/* standard initialization of the tool - * @return status code - */ -StatusCode LoKi::FastJetMaker::initialize() { - StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) { return sc; } - // check the parameters - sc = check(); - if ( sc.isFailure() ) { return sc; } - // - svc<LoKi::ILoKiSvc>( "LoKiSvc", true ); - // - if ( !m_showBanner ) fastjet::ClusterSequence::set_fastjet_banner_stream( 0 ); - return sc; -} -// ============================================================================ -// prepare the input information -// ============================================================================ -fastjet::JetDefinition LoKi::FastJetMaker::prepare( const IJetMaker::Input& input, - std::vector<fastjet::PseudoJet>& jets ) const { - // input container of "particles" - jets.reserve( input.size() ); - for ( IJetMaker::Input::const_iterator ip = input.begin(); input.end() != ip; ++ip ) { - const LHCb::Particle* p = *ip; - if ( NULL == p ) { - Warning( "Invalid input particle" ).ignore(); - continue; - } - jets.push_back( makeJet( p, to_user_index( ip - input.begin() ) ) ); - } - - // if(m_type == fastjet::plugin_algorithm) - // { - // fastjet::JetDefinition::Plugin * plugin; - // double overlap_threshold = 0.5; - // plugin = new fastjet::SISConePlugin (m_r, overlap_threshold); - // fastjet::JetDefinition jet_def(plugin); - // delete plugin; - // return jet_def ; - - //} - // else{ - - fastjet::JetFinder finder = (fastjet::JetFinder)m_type; - fastjet::RecombinationScheme scheme = (fastjet::RecombinationScheme)m_recom; - fastjet::Strategy strategy = (fastjet::Strategy)m_strat; - fastjet::JetDefinition jet_def( finder, m_r, scheme, strategy ); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "fastjet::JetDefinition:" << endmsg; - debug() << jet_def.description() << endmsg; - } - - return jet_def; - // } -} -// =========================================================================== -// find the jets -// =========================================================================== -StatusCode LoKi::FastJetMaker::makeJets( const IJetMaker::Input& input_, // - const LHCb::RecVertex& /* vtx_ */, // - IJetMaker::Jets& jets_ ) const { - makeJets( input_, jets_ ).ignore(); - return StatusCode::SUCCESS; -} - -StatusCode LoKi::FastJetMaker::makeJets( const IJetMaker::Input& input_, IJetMaker::Jets& jets_ ) const { - - StatusCode sc = check(); - if ( sc.isFailure() ) { return Error( "Invalid configuration of fastjet" ); } - - // input container of "particles" - Jets_ inputs; - // prepare the input dat and define the jets - - fastjet::JetDefinition jetDef = prepare( input_, inputs ); - - if ( inputs.empty() ) { - IJetMaker::Jets output; - output.reserve( 0 ); - jets_ = output; - if ( msgLevel( MSG::DEBUG ) ) { m_nJets += output.size(); } - return StatusCode::SUCCESS; - } - // Jets found - Jets_ jets; - - fastjet::ClusterSequence* clusters = new fastjet::ClusterSequence( inputs, jetDef ); - - switch ( m_sort ) { - case 3: - jets = sorted_by_rapidity( clusters->inclusive_jets( m_ptmin ) ); - break; - case 2: - jets = sorted_by_pt( clusters->inclusive_jets( m_ptmin ) ); - break; - case 1: - jets = sorted_by_E( clusters->inclusive_jets( m_ptmin ) ); - break; - default: - jets = sorted_by_pt( clusters->inclusive_jets( m_ptmin ) ); - break; - } - - if ( jets.empty() ) { Warning( "No jets from fastjet::ClusterSequence", StatusCode::SUCCESS, 0 ).ignore(); } - - // - if ( NULL == m_combiner ) { m_combiner = tool<IParticleCombiner>( m_combinerName, this ); } - - // Get the default geometry FIXME, use functional framework - auto lhcb = getDet<IDetectorElement>( m_standardGeometry_address ); - if ( !lhcb ) { throw GaudiException( "Could not load geometry", name(), StatusCode::FAILURE ); } - auto& geometry = *lhcb->geometry(); - - IJetMaker::Jets output; - output.reserve( jets.size() ); - - LoKi::Point3D point = LoKi::Point3D( 0, 0, 0 ); - - for ( Jets_::iterator ijet = jets.begin(); jets.end() != ijet; ++ijet ) { - const Jet& jet = *ijet; - const Constituents& constituents = clusters->constituents( jet ); - if ( constituents.empty() ) { Warning( "Jet is 'empty'!" ).ignore(); } - - LHCb::Particle pJet; - LHCb::Vertex vJet; - LHCb::Particle::ConstVector daughters; - - pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); - - pJet.setReferencePoint( point ); - - for ( Constituents::const_iterator ic = constituents.begin(); constituents.end() != ic; ++ic ) { - const Jet& c = *ic; - // find the appropriate input particle - const int index = from_user_index( c.user_index() ); - if ( 0 > index || (int)inputs.size() <= index ) { - Warning( "Invalid index for a constituent!" ).ignore(); - continue; - } // CONTINUE - // get the appropriate particle: - const LHCb::Particle* p = input_[index]; - // add the particle into the vertex - daughters.push_back( p ); - } - if ( daughters.empty() ) { - Warning( "Empty list of of daughter particles, skip it" ).ignore(); - continue; - } - // use the tool - const StatusCode sc = m_combiner->combine( daughters, pJet, vJet, geometry ); - if ( sc.isFailure() ) { - Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - continue; - } - // redefine the momentum - pJet.setMomentum( Gaudi::LorentzVector( jet.px(), jet.py(), jet.pz(), jet.e() ) ); - // - output.push_back( pJet.clone() ); - } - - if ( msgLevel( MSG::DEBUG ) ) { m_nJets += output.size(); } - - jets_ = output; - - delete clusters; - - return StatusCode::SUCCESS; -} - -DECLARE_COMPONENT( LoKi::FastJetMaker ) diff --git a/Phys/LoKiJets/src/LoKiFastJetMaker.h b/Phys/LoKiJets/src/LoKiFastJetMaker.h deleted file mode 100644 index d635e8d3982f53bd9805284218c45e54d8902970..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiFastJetMaker.h +++ /dev/null @@ -1,290 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -#pragma once - -#include "GaudiAlg/GaudiTool.h" - -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" - -#include "Event/Particle.h" -#include "LoKi/Geometry.h" -#include "LoKi/Kinematics.h" - -#include "fastjet/JetDefinition.hh" -#include "fastjet/PseudoJet.hh" - -namespace LoKi { - - /** @class FastJetMaker - * - * The most trivial, FastJet based implementaion of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Vanya BELYAEV ibelyaev@physics.syr.edu - * @author Victor COCO coco@lapp.in2p3.fr - * @date 2007-10-22 - */ - class FastJetMaker : public virtual IJetMaker, public GaudiTool { - public: - /** The main method: jet-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* jetMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets jets ; - * - * // find the jets! - * StatusCode sc = jetMaker -> makeJets ( input , jets ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iJet = jets.begin() ; - * jets.end() != iJet ; ++iJet ) - * { - * // get the jet - * LHCb::Particle* jet = *iJet ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param jets container of output jets - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according Victor Coco's studies. - * - */ - FastJetMaker( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ) - // - , m_jetID( 98 ) - // - , m_type( 0 ) - , m_recom( fastjet::E_scheme ) - , m_r( 0.7 ) - , m_ptmin( 0 ) - // - , m_sort( 2 ) - // - , m_strat( fastjet::Best ) - , m_showBanner( false ) - // - , m_combinerName( "MomentumCombiner" ) - , m_combiner( 0 ) { - // - declareInterface<IJetMaker>( this ); - // - declareProperty( "JetID", m_jetID, "Particle ID for the Jet" ); - // - declareProperty( "Type", m_type, "The type of the algorithm [0:kt,1:Cambridge,2:anti-kt,default:kt]" ); - declareProperty( "Recombination", m_recom, "Recombination scheme: fastjet::RecombinationScheme" ); - declareProperty( "RParameter", m_r ); - declareProperty( "PtMin", m_ptmin ); - // - declareProperty( "Sort", m_sort, "Sorting Criteria for jets [3:eta,2:pt,1:E,default:pt]" ); - declareProperty( "Strategy", m_strat, "JetFinding strategy, see fastjet::Strategy " ); - declareProperty( "ShowBanner", m_showBanner = false, "Print the FastJet banner if true." ); - // define momentum combiner - declareProperty( "ParticleCombiner", m_combinerName ); - } - /// destructor - virtual ~FastJetMaker() {} - - public: - /** standard initialization of the tool - * @return status code - */ - StatusCode initialize() override; - - protected: - /// make the detailed check of all parameters - inline StatusCode check() const { - // verify the algorithm type - switch ( m_type ) { - case fastjet::kt_algorithm: - break; - case fastjet::cambridge_algorithm: - break; - case fastjet::antikt_algorithm: - break; - case fastjet::genkt_algorithm: - break; - case fastjet::cambridge_for_passive_algorithm: - break; - case fastjet::genkt_for_passive_algorithm: - break; - case fastjet::plugin_algorithm: - break; - default: - return Error( "Invalid JetFinder-algorithm is specified" ); - } - // verify the recombination sceme - switch ( m_recom ) { - case fastjet::E_scheme: - break; - case fastjet::pt_scheme: - break; - case fastjet::pt2_scheme: - break; - case fastjet::Et_scheme: - break; - case fastjet::Et2_scheme: - break; - case fastjet::BIpt_scheme: - break; - case fastjet::BIpt2_scheme: - break; - default: - return Error( "Invalid RecombinationScheme is specified" ); - } - // verify the strategy - switch ( m_strat ) { - /// experimental ... - case fastjet::N2MinHeapTiled: - break; - /// fastest from about 50..10^4 - case fastjet::N2Tiled: - break; - /// legacy - case fastjet::N2PoorTiled: - break; - /// fastest below 50 - case fastjet::N2Plain: - break; - /// worse even than the usual N^3 algorithms - case fastjet::N3Dumb: - break; - /// automatic selection of the best (based on N) - case fastjet::Best: - break; - /// best of the NlnN variants -- best overall for N>10^4 - case fastjet::NlnN: - break; - /// legacy N ln N using 3pi coverage of cylinder - case fastjet::NlnN3pi: - break; - /// legacy N ln N using 4pi coverage of cylinder - case fastjet::NlnN4pi: - break; - /// Chan's closest pair method (in a variant with 4pi coverage), - /// for use exclusively wit h the Cambridge algorithm - case fastjet::NlnNCam4pi: - break; - case fastjet::NlnNCam2pi2R: - break; - case fastjet::NlnNCam: - break; // 2piMultD - /// the plugin has been used... - case fastjet::plugin_strategy: - break; - default: - return Error( "Invalid Strategy is specified" ); - } - // check the minimum momentum - if ( 0 > m_ptmin ) { Warning( "PtMin is negative " ).ignore(); } - // - return StatusCode::SUCCESS; - } - - protected: - // prepare the input information - fastjet::JetDefinition prepare( const IJetMaker::Input& input, std::vector<fastjet::PseudoJet>& jets ) const; - // - int to_user_index( const int index ) const { return index + 10000; } - int from_user_index( const int index ) const { return index - 10000; } - - private: - // the default constructor is disabled - FastJetMaker(); - // the copy constructor is disabled - FastJetMaker( const FastJetMaker& ); - // the assignement operator is disabled - FastJetMaker& operator=( const FastJetMaker& ); - - protected: - Gaudi::Property<std::string> m_standardGeometry_address{this, "StandardGeometryTop", "/dd/Structure/LHCb"}; - - // proposed jet ID - int m_jetID; ///< proposed jet ID - // KtEvent flag - int m_type; ///< KtEvent flag/mode - // Angular distance scheme - int m_angle; ///< angular distance scheme - // Recombination scheme - int m_recom; ///< recombination scheme - // R-parameter - double m_r; ///< R-parameters - // ptMin-parameter - double m_ptmin; ///< pt-min parameter - // jet sorting criteria - int m_sort; ///< jet sorting criteria - // jet finding strategy - int m_strat; ///< jet finding strategy - // print the FastJet banner - bool m_showBanner; ///< print the FastJet banner - // combiner - std::string m_combinerName; - mutable IParticleCombiner* m_combiner; ///< combiner to be used - - private: - mutable Gaudi::Accumulators::StatCounter<> m_nJets{this, "#jets"}; - }; - - typedef fastjet::PseudoJet Jet; - typedef std::vector<Jet> Jets_; - typedef std::vector<Jet> Constituents; - /// trivial function which "converts" particle into the "jet" - inline Jet makeJet( const LHCb::Particle* p, const int index ) { - if ( 0 == p ) { return Jet(); } - const Gaudi::LorentzVector& v = p->momentum(); - Jet jet( v.Px(), v.Py(), v.Pz(), v.E() ); - jet.set_user_index( index ); - return jet; - } -} // End of namespace LoKi diff --git a/Phys/LoKiJets/src/LoKiFastJetWithAreaMaker.cpp b/Phys/LoKiJets/src/LoKiFastJetWithAreaMaker.cpp deleted file mode 100644 index 8955fe40e1f7abebcf3fd844ba15a250a00b296f..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiFastJetWithAreaMaker.cpp +++ /dev/null @@ -1,318 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "LoKiFastJetMaker.h" - -#include "fastjet/ClusterSequenceArea.hh" - -#include "DetDesc/IDetectorElement.h" - -namespace LoKi { - - /** @class FastJetWithAreaMaker - * - * The most trivial, FastJet based implementation of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Vanya BELYAEV ibelyaev@physics.syr.edu - * @author Victor COCO coco@lapp.in2p3.fr - * @date 2007-10-22 - */ - class FastJetWithAreaMaker : public LoKi::FastJetMaker { - public: - /** The main method: jet-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* jetMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets jets ; - * - * // find the jets! - * StatusCode sc = jetMaker -> makeJets ( input , jets ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iJet = jets.begin() ; - * jets.end() != iJet ; ++iJet ) - * { - * // get the jet - * LHCb::Particle* jet = *iJet ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param jets container of output jets - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according Victor Coco's studies. - * - */ - FastJetWithAreaMaker( const std::string& type, const std::string& name, const IInterface* parent ) - : LoKi::FastJetMaker( type, name, parent ) - // area type - , m_area_type( 20 ) ///< the type of area method - // the maximal pseudorapidity for ghosts - , m_ghost_etamax( 8 ) ///< the maximal pseudorapidity for ghosts - // number of repetitions for active area evaluation - , m_active_area_repeats( 2 ) ///< number of repetitions for active area evaluation - // the area of the ghost - , m_ghost_area( 0.01 ) ///< the area of the ghost - // the scatterer parameter for the grid - , m_grid_scatterer( 1.e-4 ) ///< the scatterer parameter for the grid - // the scatterer parameter for kt - , m_kt_scatterer( 0.1 ) ///< the scatterer parameter for kt - // the kt for ghosts - , m_mean_ghost_kt( 1.e-100 ) ///< the kt for ghosts - // pt/area strategy - , m_ptarea_strategy( 0 ) /// deprecated... to put back - // pt/area strategy - , m_ptarea_range( 2 ) ///< pt/area range - , m_effRfact( 1. ) { - // - declareProperty( "AreaType", m_area_type, "The method to get the area calculation" ); - // - declareProperty( "GhostEtaMax", m_ghost_etamax, "The maximal pseudorapidity for active area ghosts" ); - // - declareProperty( "AreaRepeat", m_active_area_repeats, "Repeate active area evaluation" ); - // - declareProperty( "GhostArea", m_ghost_area, "The area, associated with single ghost" ); - // - declareProperty( "GridScatterer", m_grid_scatterer, - "The fractional random fluctiatuon of ghosts on eta-phi grid" ); - // - declareProperty( "KtScatterer", m_kt_scatterer, "The fractional random fluctiatuon of ghosts kt" ); - // - declareProperty( "GhostKt", m_mean_ghost_kt, "The average transverse momentum of the ghosts" ); - // - declareProperty( "PtPerAreaStrategy", m_ptarea_strategy, "The strategy for pt/unit area evaluation" ); - // - declareProperty( "PtPerAreaRange", m_ptarea_range, "The range for pt/unit area evaluation" ); - // - declareProperty( "VoronoiEffectiveRFactor", m_effRfact, "Effective R Factor for Voronoi area determination" ); - // - }; - - private: - // the method to get the area - int m_area_type; ///< the method to get the area - // the maximal pseudorapidity for ghosts - double m_ghost_etamax; ///< the maximal pseudorapidity for ghosts - // number of repetitions for active area evaluation - int m_active_area_repeats; ///< number of repetitions for active area evaluation - // the area of the ghost - double m_ghost_area; ///< the area of the ghosts - // the scatterer parameter for the grid - double m_grid_scatterer; // the scatterer parameter for the grid - // the scatterer parameter for kt - double m_kt_scatterer; ///< the scatterer parameter for kt - // the kt for ghosts - double m_mean_ghost_kt; ///< the kt for ghosts - // pt/area strategy - unsigned int m_ptarea_strategy; ///< pt/area strategy - // pt/area range - double m_ptarea_range; ///< pt/area range - // effective R factor for voronoi area - double m_effRfact; ///< effective R factor for voronoi area - - mutable Gaudi::Accumulators::StatCounter<> m_nJets{this, "#jets"}; - }; -} // End of namespace LoKi - -/** @file - * Implementation file for class LoKi::FastJetWithAreaMaker - * @author Vanya BELYAEV ibelyaev@physics.syr.edu - * @author Victor COCO coco@lapp.in2p3.fr - * @date 2007-10-22 - */ - -// find the jets -StatusCode LoKi::FastJetWithAreaMaker::makeJets( const IJetMaker::Input& input_, const LHCb::RecVertex& /* vtx_ */, - IJetMaker::Jets& jets_ ) const { - makeJets( input_, jets_ ).ignore(); - return StatusCode::SUCCESS; -} - -StatusCode LoKi::FastJetWithAreaMaker::makeJets( const IJetMaker::Input& input_, IJetMaker::Jets& jets_ ) const { - - StatusCode sc = check(); - if ( sc.isFailure() ) { return Error( "Invalid configurtaion of fastjet" ); } - - // input data - Jets_ inputs; - - // prepare the input dat and define the jets - fastjet::JetDefinition jet_def = prepare( input_, inputs ); - - // Jets found - Jets_ jets; - - fastjet::AreaDefinition area_def; - - if ( m_area_type == -1 ) { - return Error( "Invalid configurtaion of area type in fastjet" ); - } else if ( m_area_type < 11 ) { - // specify the active area - fastjet::GhostedAreaSpec active_area( m_ghost_etamax, m_active_area_repeats, m_ghost_area, m_grid_scatterer, - m_kt_scatterer, m_mean_ghost_kt ); - fastjet::AreaDefinition myarea_def( fastjet::active_area, active_area ); - area_def = myarea_def; - - } else if ( m_area_type == 20 ) { - fastjet::VoronoiAreaSpec vor_area( m_effRfact ); - fastjet::AreaDefinition myarea_def( fastjet::voronoi_area, vor_area ); - area_def = myarea_def; - } else { - fastjet::AreaDefinition myarea_def( m_area_type ); - area_def = myarea_def; - } - - // clusterisation sequence - fastjet::ClusterSequenceArea clusters // fastjet::ClusterSequenceWithArea clusters - ( inputs, jet_def, area_def ); - - switch ( m_sort ) { - case 3: - jets = sorted_by_rapidity( clusters.inclusive_jets( m_ptmin ) ); - break; - case 2: - jets = sorted_by_pt( clusters.inclusive_jets( m_ptmin ) ); - break; - case 1: - jets = sorted_by_E( clusters.inclusive_jets( m_ptmin ) ); - break; - default: - jets = sorted_by_pt( clusters.inclusive_jets( m_ptmin ) ); - break; - } - - if ( jets.empty() ) { Warning( "No jets from fastjet::ClusterSequenceWithArea" ).ignore(); } - - // - if ( 0 == m_combiner ) { m_combiner = tool<IParticleCombiner>( m_combinerName, this ); } - - // Get the default geometry FIXME, use functional framework - auto lhcb = getDet<IDetectorElement>( m_standardGeometry_address ); - if ( !lhcb ) { throw GaudiException( "Could not load geometry", name(), StatusCode::FAILURE ); } - auto& geometry = *lhcb->geometry(); - - IJetMaker::Jets output; - output.reserve( jets.size() ); - - LoKi::Point3D point = LoKi::Point3D( 0, 0, 0 ); - - // get the "pt per unit area" estimate - // const double ptPerUnitArea = clusters.parabolic_pt_per_unit_area - // ( (fastjet::ClusterSequenceArea::mean_pt_strategies) m_ptarea_strategy , - // m_ptarea_range ) ; - - for ( Jets_::iterator ijet = jets.begin(); jets.end() != ijet; ++ijet ) { - const Jet& jet = *ijet; - const Constituents& constituents = clusters.constituents( jet ); - if ( constituents.empty() ) { Warning( "Jet is 'empty'!" ).ignore(); } - - LHCb::Particle pJet, pJetArea; - LHCb::Vertex vJet; - LHCb::Particle::ConstVector daughters; - - pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); - pJetArea.setParticleID( LHCb::ParticleID( -1 * m_jetID ) ); - - pJet.setReferencePoint( point ); - pJetArea.setReferencePoint( point ); - - // set the jet active area - pJet.addInfo( LHCb::Particle::JetActiveArea, clusters.area( jet ) ); - // set the jet active area uncertainty - pJet.addInfo( LHCb::Particle::JetActiveAreaError, clusters.area_error( jet ) ); - // get active area 4-vector - fastjet::PseudoJet area = clusters.area_4vector( jet ); - pJet.addInfo( LHCb::Particle::JetActiveAreaPx, area.px() ); - pJet.addInfo( LHCb::Particle::JetActiveAreaPy, area.py() ); - pJet.addInfo( LHCb::Particle::JetActiveAreaPz, area.pz() ); - pJet.addInfo( LHCb::Particle::JetActiveAreaE, area.e() ); - pJet.addInfo( LHCb::Particle::JetPtPerUnitArea, 0. ); - - for ( Constituents::const_iterator ic = constituents.begin(); constituents.end() != ic; ++ic ) { - const Jet& c = *ic; - // find the appropriate input particle - const int index = from_user_index( c.user_index() ); - if ( 0 > index || (int)inputs.size() <= index ) { - Warning( "Invalid index for a constituent!" ).ignore(); - continue; - } // CONTINUE - // get the appropriate particle: - const LHCb::Particle* p = input_[index]; - // add the particle into the vertex - daughters.push_back( p ); - } - if ( daughters.empty() ) { - Warning( "Empty list of of daughter particles, skip it" ).ignore(); - continue; - } - // use the tool - StatusCode sc = m_combiner->combine( daughters, pJet, vJet, geometry ); - if ( sc.isFailure() ) { - Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - continue; - } - // redefine the momentum - pJet.setMomentum( Gaudi::LorentzVector( jet.px(), jet.py(), jet.pz(), jet.e() ) ); - // - - output.push_back( pJet.clone() ); - } - - if ( msgLevel( MSG::DEBUG ) ) { m_nJets += output.size(); } - - jets_ = output; - - return StatusCode::SUCCESS; -} - -DECLARE_COMPONENT( LoKi::FastJetWithAreaMaker ) diff --git a/Phys/LoKiJets/src/LoKiJetMakerAlg.cpp b/Phys/LoKiJets/src/LoKiJetMakerAlg.cpp deleted file mode 100644 index e6e5f766dfa70a0b01b88d4f9baba824210efa35..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiJetMakerAlg.cpp +++ /dev/null @@ -1,358 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "LoKi/Algo.h" -#include "LoKi/ParticleContextCuts.h" -#include "LoKi/ParticleCuts.h" -#include "LoKi/VertexCuts.h" - -#include "Kernel/DaVinciAlgorithm.h" -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" -#include "Kernel/PVSentry.h" -#include "Kernel/Particle2Vertex.h" - -#include "Event/Particle.h" - -#include "AIDA/IHistogram2D.h" -#include "GaudiKernel/IHistogramSvc.h" -#include "boost/lexical_cast.hpp" -#include <GaudiUtils/Aida2ROOT.h> - -#include "Math/VectorUtil.h" -#include "TH2D.h" -#include "TMath.h" - -class TH2D; - -namespace LoKi { - - /** @class JetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Vanya BELYAEV belyaev@lapp.in2p3.fr - * @date 2005-03-21 - */ - class JetMaker : public LoKi::Algo { - public: - /** Standard constructor - * @param name instance name - * @param pSvc pointer to Service Locator - */ - JetMaker( const std::string& name, ISvcLocator* pSvc ) - : LoKi::Algo( name, pSvc ) - // - , m_makerName( "LoKi::FastJetMaker" ) - , m_maker( 0 ) - , m_associate2Vertex( false ) - , m_runJetID( true ) { - // - declareProperty( "JetMaker", m_makerName, "Type type/name of jet-maker tool (IJetMaker interface)" ); - declareProperty( "Associate2Vertex", m_associate2Vertex, "Jet reconstruction per vertex" ); - declareProperty( "JetID", m_runJetID, "Append Jet ID information to the jet" ); - declareProperty( "ApplyJetID", m_applyJetID = false, "Apply jet ID cuts" ); - declareProperty( "ApplyJEC", m_applyJEC = false ); - - declareProperty( "HistoPath", m_histo_path = "JEC" ); - // - } - - public: - /** standard execution of the algorithm - * @see LoKi::Algo - * @return status code - */ - StatusCode initialize() override; - StatusCode analyse() override; - - StatusCode appendJetIDInfo( LHCb::Particle* jet ); - StatusCode JEC( LHCb::Particle* jet ); - - private: - /// maker name - std::string m_makerName; // jet maker name - /// maker - const IJetMaker* m_maker; // jet maker to be used - /// associate two vertex? - bool m_associate2Vertex; // make jet per vertex - /// append jet ID info - bool m_runJetID; // append jet ID info - /// apply JEC - bool m_applyJEC; - /// apply JEC - bool m_applyJetID; - /// histograms for JEC - std::vector<TH2D*> m_histosJEC; - /// histo path - std::string m_histo_path; - - mutable Gaudi::Accumulators::StatCounter<> m_nJets{this, "#jets"}; - }; - -} // end of namespace LoKi - -/** @file - * Implementation file for class LoKi::JetMaker - * @date 2005-03-21 - * @author Vanya BELYAEV belyaev@lapp.in2p3.fr - */ - -/* standard execution of the algorithm - * @see LoKi::Algo - * @return status code - */ - -StatusCode LoKi::JetMaker::initialize() { - StatusCode sc = LoKi::Algo::initialize(); - if ( sc.isFailure() ) { return sc; } - // Read in the histograms: - if ( m_applyJEC ) { - for ( int i = 1; i < 4; ++i ) { - std::string histoname = "JEC_PV" + boost::lexical_cast<std::string>( i ); - AIDA::IHistogram2D* aida = get<AIDA::IHistogram2D>( histoSvc(), m_histo_path + histoname ); - if ( 0 == aida ) warning() << "Could not find AIDA::IHistogram2D* " << m_histo_path + histoname << "." << endmsg; - m_histosJEC.push_back( Gaudi::Utils::Aida2ROOT::aida2root( aida ) ); - } - } - // - return sc; -} - -StatusCode LoKi::JetMaker::analyse() { - using namespace LoKi; - using namespace LoKi::Types; - using namespace LoKi::Cuts; - - if ( m_associate2Vertex ) { - - // A cut to get the x position of the bestPV of input particles (would be better to code a VKEY functor) - LoKi::Types::Fun bestVertexVX = BPV( this, VX ); - LoKi::Types::Fun bestVertexVY = BPV( this, VY ); - LoKi::Types::Fun bestVertexVZ = BPV( this, VZ ); - - // A cut to check that a jet contains information able to link it to a PV - LoKi::Types::Cut withPVPointingInfo = - NINTREE( ( ABSID == 310 || ABSID == 3122 ) || - ( HASTRACK && static_cast<int>( LHCb::Track::Types::Downstream ) != TRTYPE ) ) > 0; - - Range part = select( "part", PALL ); - - // Loop over PV list and make jets out of appropriate inputs related to the PVs - const LHCb::RecVertex::Range pvs = this->primaryVertices(); - for ( LHCb::RecVertex::Range::const_iterator i_pv = pvs.begin(); pvs.end() != i_pv; i_pv++ ) { - IJetMaker::Input inputs; - for ( Range::const_iterator i_p = part.begin(); part.end() != i_p; i_p++ ) { - // Neutral inputs - if ( Q( *i_p ) == 0 ) { - // Take all neutrals exept V0s - if ( ABSID( *i_p ) != 310 && ABSID( *i_p ) != 3122 ) inputs.push_back( *i_p ); - // Keep only the V0s that match to the vertex - else if ( bestVertexVX( *i_p ) == VX( *i_pv ) && bestVertexVY( *i_p ) == VY( *i_pv ) && - bestVertexVZ( *i_p ) == VZ( *i_pv ) ) - inputs.push_back( *i_p ); - else - continue; - } - // Charged inputs - else { - // Take all downstream tracks - if ( static_cast<int>( LHCb::Track::Types::Downstream ) == TRTYPE( *i_p ) ) inputs.push_back( *i_p ); - // Keep only the tracks with velo segment that match to the vertex - else if ( bestVertexVX( *i_p ) == VX( *i_pv ) && bestVertexVY( *i_p ) == VY( *i_pv ) && - bestVertexVZ( *i_p ) == VZ( *i_pv ) ) - inputs.push_back( *i_p ); - else - continue; - } - } - - // input container of "particles" - IJetMaker::Jets jets; - - LoKi::Types::Fun mtf = LoKi::Cuts::INFO( 9003, -10. ); - LoKi::Types::Fun n90 = LoKi::Cuts::INFO( 9002, -10. ); - LoKi::Types::Fun cpf = LoKi::Cuts::INFO( 9006, -10. ); - LoKi::Types::Fun width = LoKi::Cuts::INFO( 9007, -10. ); - LoKi::Types::Fun nPVInfo = LoKi::Cuts::INFO( 9005, -10. ); - if ( 0 == m_maker ) { m_maker = tool<IJetMaker>( m_makerName, m_makerName, this ); } - - // make the jets - StatusCode sc = m_maker->makeJets( inputs.begin(), inputs.end(), jets ); - - if ( sc.isFailure() ) { return Error( "Error from jet maker", sc ); } - - // save all jets - - while ( !jets.empty() ) { - LHCb::Particle* jet = jets.back(); - if ( m_applyJEC ) this->JEC( jet ).ignore(); - if ( m_runJetID ) this->appendJetIDInfo( jet ).ignore(); - // If the jet contain info on PV, assign a PV and update the P2PV relation table - if ( withPVPointingInfo( jet ) ) { - jet->setReferencePoint( Gaudi::XYZPoint( ( *i_pv )->position() ) ); - LHCb::Vertex* vJet = new LHCb::Vertex(); - vJet->setPosition( ( *i_pv )->position() ); - vJet->setCovMatrix( ( *i_pv )->covMatrix() ); - vJet->setChi2( ( *i_pv )->chi2() ); - vJet->setNDoF( ( *i_pv )->nDoF() ); - vJet->setOutgoingParticles( jet->daughters() ); - jet->setEndVertex( vJet ); - this->relate( jet, *i_pv ); - } - if ( !( m_applyJetID && m_runJetID && ( mtf( jet ) > 0.75 || nPVInfo( jet ) < 2 ) ) ) { - sc = save( "jets", jet ); - if ( sc.isFailure() ) { return Error( "Error from save function in jet maker", sc ); } - } - jets.pop_back(); - unRelatePV( jet ); - delete jet->endVertex(); - delete jet; - } - } - } else { - - Range all = select( "all", LoKi::Cuts::PALL ); - // input container of "particles" - IJetMaker::Jets jets; - - if ( 0 == m_maker ) { m_maker = tool<IJetMaker>( m_makerName, m_makerName, this ); } - - // make the jets - StatusCode sc = m_maker->makeJets( all.begin(), all.end(), jets ); - - if ( sc.isFailure() ) { return Error( "Error from jet maker", sc ); } - - LoKi::Types::Fun mtf = LoKi::Cuts::INFO( 9003, -10. ); - LoKi::Types::Fun n90 = LoKi::Cuts::INFO( 9002, -10. ); - LoKi::Types::Fun cpf = LoKi::Cuts::INFO( 9006, -10. ); - LoKi::Types::Fun width = LoKi::Cuts::INFO( 9007, -10. ); - LoKi::Types::Fun nPVInfo = LoKi::Cuts::INFO( 9005, -10. ); - - // save all jets - while ( !jets.empty() ) { - LHCb::Particle* jet = jets.back(); - if ( m_runJetID ) this->appendJetIDInfo( jet ).ignore(); - if ( m_applyJetID && m_runJetID && ( mtf( jet ) > 0.75 || nPVInfo( jet ) < 2 ) ) { - jets.pop_back(); - delete jet; - continue; - } - - sc = save( "jets", jet ); - if ( sc.isFailure() ) { return Error( "Error from save function in jet maker", sc ); } - - jets.pop_back(); - delete jet; - } - } - - if ( msgLevel( MSG::DEBUG ) ) { m_nJets += selected( "jets" ).size(); } - setFilterPassed( true ); - - return StatusCode::SUCCESS; -} - -StatusCode LoKi::JetMaker::appendJetIDInfo( LHCb::Particle* jet ) { - // Get Jet Daughters - std::vector<const LHCb::Particle*> daughtersvector = jet->daughtersVector(); - std::vector<const LHCb::Particle*>::iterator idaughter = daughtersvector.begin(); - - double mtf; /// Highest pT track / Jet pT - double cpf; /// charged pT fraction - V0s are not included - double width; /// jet width - int n90; /// Number of items responsible for at least 90% of the jet momentum - int ntrk; /// Number of tracks - - float auxptmax = -1, sumpt = 0; - int iitems = 0; - double tpx = 0, tpy = 0; - std::vector<float> itemspt; - ntrk = n90 = width = 0; - - for ( ; idaughter != daughtersvector.end(); ++idaughter ) { - const LHCb::Particle* daughter = *idaughter; - if ( daughter->particleID().threeCharge() != 0 ) { - ntrk++; - auxptmax = daughter->momentum().Pt() > auxptmax ? daughter->momentum().Pt() : auxptmax; - tpx += daughter->momentum().Px(); - tpy += daughter->momentum().Py(); - } - iitems++; - float pt = daughter->momentum().Pt(); - sumpt += pt; - itemspt.push_back( pt ); - for ( int ii = 0; ii < iitems; ii++ ) { - if ( itemspt[ii] < pt ) { - float aux = itemspt[ii]; - itemspt[ii] = pt; - pt = aux; - } - } - width += ROOT::Math::VectorUtil::DeltaR( daughter->momentum(), jet->momentum() ) * daughter->momentum().Pt(); - } - - mtf = auxptmax / jet->momentum().Pt(); - mtf = 0 > mtf ? 0 : mtf; - mtf = 1 < mtf ? 1 : mtf; - cpf = std::sqrt( tpx * tpx + tpy * tpy ) / jet->momentum().Pt(); - width /= sumpt; - - sort( itemspt.begin(), itemspt.end() ); - float auxptsum = 0; - n90 = 0; - for ( int ii = iitems - 1; ii >= 0; ii-- ) { - auxptsum += itemspt[ii]; - n90++; - if ( auxptsum / sumpt > 0.9 ) break; - } - - LoKi::Types::Fun NsatCells = LoKi::Cuts::SUMTREE( LoKi::Cuts::INFO( 955, 0. ), LoKi::Cuts::Q == 0, 0. ); - LoKi::Types::Fun N_HasPVInfo = LoKi::Cuts::NINTREE( - ( LoKi::Cuts::ABSID == 310 || LoKi::Cuts::ABSID == 3122 ) || - ( LoKi::Cuts::HASTRACK && static_cast<int>( LHCb::Track::Types::Downstream ) != LoKi::Cuts::TRTYPE ) ); - - jet->addInfo( 9001, ntrk ); - jet->addInfo( 9002, n90 ); - jet->addInfo( 9003, mtf ); - jet->addInfo( 9006, cpf ); - jet->addInfo( 9007, width ); - jet->addInfo( 9004, NsatCells( jet ) ); - jet->addInfo( 9005, N_HasPVInfo( jet ) ); - - return StatusCode::SUCCESS; -} - -StatusCode LoKi::JetMaker::JEC( LHCb::Particle* jet ) { - const int PV = this->primaryVertices().size(); - const int usePV = ( PV > 3 ? 3 : PV ); - const TH2D* histo = m_histosJEC[usePV - 1]; - const int xbin = histo->GetXaxis()->FindFixBin( LoKi::Cuts::PT( jet ) / 1000. ); - const int ybin = histo->GetYaxis()->FindFixBin( LoKi::Cuts::ETA( jet ) ); - const double cor = histo->GetBinContent( xbin, ybin ); - // Store the uncorrected kinematics - - jet->addInfo( 9100, cor ); - jet->addInfo( 9101, PV ); - - jet->setMomentum( Gaudi::LorentzVector( cor * LoKi::Cuts::PX( jet ), cor * LoKi::Cuts::PY( jet ), - cor * LoKi::Cuts::PZ( jet ), cor * LoKi::Cuts::E( jet ) ) ); - - return StatusCode::SUCCESS; -} - -DECLARE_COMPONENT( LoKi::JetMaker ) diff --git a/Phys/LoKiJets/src/LoKiJetMakerWR2VtxAlg.cpp b/Phys/LoKiJets/src/LoKiJetMakerWR2VtxAlg.cpp deleted file mode 100644 index c878b599f16eeed619ab860f8250f2770d6cde92..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiJetMakerWR2VtxAlg.cpp +++ /dev/null @@ -1,184 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "LoKi/Algo.h" -#include "LoKi/ParticleCuts.h" - -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" - -#include "Event/Particle.h" - -namespace LoKi { - - /** @class JetMakerWR2Vtx - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Vanya BELYAEV belyaev@lapp.in2p3.fr - * @date 2005-03-21 - */ - class JetMakerWR2Vtx : public LoKi::Algo { - public: - /** Standard constructor - * @param name instance name - * @param pSvc pointer to Service Locator - */ - JetMakerWR2Vtx( const std::string& name, ISvcLocator* pSvc ) - : LoKi::Algo( name, pSvc ) - // - , m_makerName( "LoKi::FastJetMaker" ) - , m_maker( 0 ) - , m_VtxLoc( "" ) { - // - declareProperty( "JetMaker", m_makerName, "Type type/name of jet-maker tool (IJetMaker interface)" ); - declareProperty( "InputRecVtxLoc", m_VtxLoc, "RecVtx Input location" ); - - // - } - - /** standard execution of the algorithm - * @see LoKi::Algo - * @return status code - */ - StatusCode analyse() override; - - private: - /// maker name - std::string m_makerName; // jet maker name - /// maker - const IJetMaker* m_maker; // jet maker to be used - std::string m_VtxLoc; - - mutable Gaudi::Accumulators::StatCounter<> m_nJets{this, "#jets"}; - }; - -} // end of namespace LoKi - -/** @file - * Implementation file for class LoKi::JetMakerWR2Vtx - * @date 2005-03-21 - * @author Vanya BELYAEV belyaev@lapp.in2p3.fr - */ - -/* standard execution of the algorithm - * @see LoKi::Algo - * @return status code - */ -StatusCode LoKi::JetMakerWR2Vtx::analyse() { - using namespace LoKi; - using namespace LoKi::Types; - // select all input particles - Range all = select( "all", LoKi::Cuts::PALL ); - - // input container of "particles" - const std::vector<const LHCb::Particle*> input( all.begin(), all.end() ); - - // const LHCb::RecVertex vtx; - if ( 0 == m_maker ) { m_maker = tool<IJetMaker>( m_makerName, m_makerName, this ); } - - if ( m_VtxLoc == "" ) { - - IJetMaker::Jets jets; - - StatusCode sc = m_maker->makeJets( input, jets ); - if ( sc.isFailure() ) { return Error( "Error from seed maker", sc ); } - - int myj = 0; - while ( !jets.empty() ) { - - LHCb::Particle* jet = jets.back(); - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "SeedsP " << myj << "| nb of combin Seeds: " << jet->weight() - << " nb part in the jet: " << jet->daughtersVector().size() << " Pt: " << jet->pt() << endmsg; - debug() << " | m : " << jet->momentum() << endmsg; - debug() << " | eta: " << jet->momentum().eta() << " / phi: " << jet->momentum().phi() << endmsg; - debug() << " |Chi2: " << jet->endVertex()->chi2() << endmsg; - debug() << " |nDoF: " << jet->endVertex()->nDoF() << endmsg; - debug() << " |C/N : " << jet->endVertex()->chi2PerDoF() << endmsg; - } - - save( "jets", jet ).ignore(); - jets.pop_back(); - delete jet; - } - } - - if ( !exist<LHCb::RecVertex::Range>( m_VtxLoc ) ) return Error( "the VtxLoc do not exsit" ); - - const LHCb::RecVertex::Range PVs = get<LHCb::RecVertex::Range>( m_VtxLoc ); - - if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx size : " << PVs.size() << endmsg; - - int i = 0; - // std::vector<int> sjets; - for ( LHCb::RecVertex::Range::const_iterator ipv2 = PVs.begin(); ipv2 != PVs.end(); ++ipv2 ) { - // sjets.at(i) = 0; - - if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx: " << i << endmsg; - - IJetMaker::Jets jets; - - const LHCb::RecVertex vtx = **ipv2; - - // make the jets - StatusCode sc = m_maker->makeJets( input, vtx, jets ); - - if ( sc.isFailure() ) { return Error( "Error from seed maker", sc ); } - - // save all jets - int myj = 0; - - while ( !jets.empty() ) { - - LHCb::Particle* jet = jets.back(); - // sjets.at(i)++; - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "SeedsP " << myj << "| nb of combin Seeds: " << jet->weight() - << " nb part in the jet: " << jet->daughtersVector().size() << " Pt: " << jet->pt() << endmsg; - debug() << " | m : " << jet->momentum() << endmsg; - debug() << " | eta: " << jet->momentum().eta() << " / phi: " << jet->momentum().phi() << endmsg; - debug() << " |Chi2: " << jet->endVertex()->chi2() << endmsg; - debug() << " |nDoF: " << jet->endVertex()->nDoF() << endmsg; - debug() << " |C/N : " << jet->endVertex()->chi2PerDoF() << endmsg; - } - - save( "jets", jet ).ignore(); - jets.pop_back(); - delete jet; - } - // i++; - } - - // int kk=0; - // for(int jj = 0; jj != (int) sjets.size(); jj++) - // if( sjets.at(i) > 0) - // kk++; - - // if(kk > 1) std::cout<<"evt :"<<evt<<std::endl; - - if ( msgLevel( MSG::DEBUG ) ) { m_nJets += selected( "jets" ).size(); } - - setFilterPassed( true ); - - return StatusCode::SUCCESS; -} - -DECLARE_COMPONENT( LoKi::JetMakerWR2Vtx ) diff --git a/Phys/LoKiJets/src/LoKiJetParticleMaker.cpp b/Phys/LoKiJets/src/LoKiJetParticleMaker.cpp deleted file mode 100644 index 4aa1d22615caf7a58a67f0e7c447c4e1dd4a9916..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiJetParticleMaker.cpp +++ /dev/null @@ -1,120 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -// ============================================================================ -// include files -// ============================================================================ -// STD & STL -// ============================================================================ -#include <string> -// ============================================================================ -// GaudiAlg -// ============================================================================ -#include "GaudiAlg/GaudiTool.h" -// ============================================================================ -// DaVinci Kernel -// ============================================================================ -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleMaker.h" -// ============================================================================ -namespace LoKi { - // ========================================================================== - /** @class JetParticleMaker - * The simple implementation of interface IParticleMaker for building the - * jets using IJetMaker tool - * @see IParticleMaker - * @see IJetMaker - * @author Vanya BELYAEV ibelyaev@physics.syr.edu - * @date 2007-10-15 - */ - class JetParticleMaker : public virtual IParticleMaker, public GaudiTool { - public: - // ======================================================================== - /// make the particles - StatusCode makeParticles( LHCb::Particle::ConstVector& particles ) override; - /// initialize the tool - StatusCode initialize() override { - StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) { return sc; } - // locate the jet-finder - m_jetMaker = tool<IJetMaker>( m_jetMakerName, this ); - // check the inputs - if ( m_inputs.empty() ) { Warning( "No input locations are specified!" ).ignore(); } - return StatusCode::SUCCESS; - } - // ======================================================================== - /// The standard constructor - JetParticleMaker( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ), m_jetMakerName( "LoKi::FastJetMaker" ), m_jetMaker( 0 ) { - declareInterface<IParticleMaker>( this ); - - declareProperty( "JetMaker", m_jetMakerName, "The type/name of Jet-Maker tool" ); - declareProperty( "Inputs", m_inputs, "The list of input locations " ); - } - // ======================================================================== - private: - // ======================================================================== - // the jet maker name - std::string m_jetMakerName; ///< the jet-maker name - // the jet maker itself - const IJetMaker* m_jetMaker; ///< the jet-maker itself - typedef std::vector<std::string> Inputs; - // input locations - Inputs m_inputs; ///< the input locations - // counters - mutable Gaudi::Accumulators::SummingCounter<unsigned int> m_nbJetsCounter{this, "#jets"}; - // ======================================================================== - }; - // ========================================================================== -} // end of namespace LoKi -// ============================================================================ -// make the particles -// ============================================================================ -StatusCode LoKi::JetParticleMaker::makeParticles( LHCb::Particle::ConstVector& particles ) { - - IJetMaker::Input inputs; - // loop over all input locations - for ( Inputs::const_iterator i = m_inputs.begin(); m_inputs.end() != i; ++i ) { - if ( exist<LHCb::Particle::Range>( *i ) ) { - const LHCb::Particle::Range parts = get<LHCb::Particle::Range>( *i ); - inputs.insert( inputs.end(), parts.begin(), parts.end() ); - } else if ( exist<LHCb::Particle::Range>( ( *i ) + "/Particles" ) ) { - const LHCb::Particle::Range parts = get<LHCb::Particle::Range>( ( *i ) + "/Particles" ); - inputs.insert( inputs.end(), parts.begin(), parts.end() ); - } else { - return Error( "No valid location: " + ( *i ) ); - } - } - - // empty container of conmstituents? - if ( inputs.empty() ) { Warning( "Empty container of input particles " ).ignore(); } - - // prepare the container for output jets - IJetMaker::Jets jets; - jets.reserve( 100 ); - - // use Jet Maker tool: - StatusCode sc = m_jetMaker->makeJets( inputs, jets ); - if ( sc.isFailure() ) { return Error( "The error form JetMaker", sc ); } - - // add the jets to the output container: - particles.insert( particles.end(), jets.begin(), jets.end() ); - - // some statistics - m_nbJetsCounter += jets.size(); - - return StatusCode::SUCCESS; -} -// ============================================================================ -/// the factory -DECLARE_COMPONENT( LoKi::JetParticleMaker ) -// ============================================================================] -// The END -// ============================================================================ diff --git a/Phys/LoKiJets/src/LoKiJets2JetsAlg.cpp b/Phys/LoKiJets/src/LoKiJets2JetsAlg.cpp deleted file mode 100644 index 3361a3bf2d07e5c96ad97c5a51b3f3a7f1928c31..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiJets2JetsAlg.cpp +++ /dev/null @@ -1,144 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -// ============================================================================ -// Include files -// ============================================================================ -#include "GaudiAlg/GaudiAlgorithm.h" -// ============================================================================ -// DaVinci Kernel -// ============================================================================ -#include "Kernel/IJets2Jets.h" -// =========================================================================== -// Event -// ============================================================================ -#include "Event/Particle.h" -// ============================================================================ - -// ========================================================================== -/** @class Jets2JetsAlg - * @date 2009-10-29 - * @author Victor Coco victor.coco@cern.ch - * - */ -class Jets2JetsAlg : public GaudiAlgorithm { - // ======================================================================== - /// the friend factory for instantiation - // ======================================================================== -public: - // ======================================================================== - /** Standard constructor - * @param name instance name - * @param pSvc pointer to Service Locator - */ - Jets2JetsAlg( const std::string& name, ISvcLocator* pSvc ) - : GaudiAlgorithm( name, pSvc ) - // - , m_jetMatcherName() - , m_jetALocation() - , m_jetBLocation() - , m_output() - , m_jetMatcher( 0 ) { - // - declareProperty( "Jets2Jets", m_jetMatcherName, "Type type/name of jet-jetMatcher tool (IJets2Jets interface)" ); - declareProperty( "JetsALocation", m_jetALocation, "Location of the first set of jet" ); - declareProperty( "JetsBLocation", m_jetBLocation, "Location of the second set of jet" ); - declareProperty( "OutputTable", m_output, "Location of the output table" ); - // - } - // ======================================================================== -public: - // ======================================================================== - /** standard execution of the algorithm - * @see LoKi::Algo - * @return status code - */ - StatusCode initialize() override; - StatusCode execute() override; - // ======================================================================== - - // ======================================================================== -private: - // ======================================================================== - /// jetMatcher name - std::string m_jetMatcherName; ///< jet jetMatcher name - std::string m_jetALocation; ///< Location of the first set of jet - std::string m_jetBLocation; ///< Location of the second set of jet - std::string m_output; ///< Location of the second set of jet - /// jetMatcher - const IJets2Jets* m_jetMatcher; ///< jet jetMatcher to be used - // ======================================================================== -}; -// ========================================================================== -// end of namespace LoKi -// ============================================================================ -/** @file - * Implementation file for class Jets2JetsAlg - * Client for implementation of IJets2Jets, furnishing the sets of jets out of - * transiant event store location and recording the resulting output table - * - * @date 2009-10-29 - * @author Victor COCO victor.coco@cern.ch - */ - -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( Jets2JetsAlg ) - -StatusCode Jets2JetsAlg::initialize() { - StatusCode sc = GaudiAlgorithm::initialize(); - if ( sc.isFailure() ) return sc; - if ( m_jetALocation.empty() || m_jetBLocation.empty() ) return Error( "No input locations is specified" ); - if ( m_jetMatcherName.empty() || m_jetBLocation.empty() ) return Error( "No jet matching tool specified" ); - if ( m_output.empty() ) return Error( "No OutputTable is specified" ); - return sc; -} - -StatusCode Jets2JetsAlg::execute() { - // Get the jets - const LHCb::Particles* PartsA = get<LHCb::Particles>( m_jetALocation ); - const LHCb::Particles* PartsB = get<LHCb::Particles>( m_jetBLocation ); - if ( PartsA->size() == 0 || PartsB->size() == 0 ) { - Warning( "No jets to match in the event" ).ignore(); - setFilterPassed( false ); - return StatusCode::SUCCESS; - } - IJets2Jets::Jets* JetsA = new IJets2Jets::Jets(); - IJets2Jets::Jets* JetsB = new IJets2Jets::Jets(); - for ( LHCb::Particles::const_iterator ip = PartsA->begin(); PartsA->end() != ip; ip++ ) JetsA->push_back( *ip ); - for ( LHCb::Particles::const_iterator ip = PartsB->begin(); PartsB->end() != ip; ip++ ) JetsB->push_back( *ip ); - if ( 0 == m_jetMatcher ) m_jetMatcher = tool<IJets2Jets>( m_jetMatcherName, m_jetMatcherName, this ); - - // create the relation table - IJets2Jets::Table* table = new IJets2Jets::Table(); - put( table, m_output ); - - // call the tool - m_jetMatcher->makeRelation( *JetsA, *JetsB, *table ); - - // total number of established links - const size_t links = table->relations().size(); - - // make a statistics - if ( msgLevel( MSG::VERBOSE ) ) { counter( "# " + m_jetALocation + "->" + m_jetBLocation ) += links; } - // - if ( table->relations().empty() ) { - Warning( "The relation table is empty!" ).ignore(); - setFilterPassed( false ); - } - setFilterPassed( true ); - - return StatusCode::SUCCESS; -} - -// =========================================================================== -/// The factory -// ============================================================================ -// The END -// ============================================================================ diff --git a/Phys/LoKiJets/src/LoKiSeedConeJetMaker.cpp b/Phys/LoKiJets/src/LoKiSeedConeJetMaker.cpp deleted file mode 100644 index 4b7d0b3a465eaa4196e40d530ffd402112552592..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiSeedConeJetMaker.cpp +++ /dev/null @@ -1,508 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" - -#include "LoKi/Geometry.h" -#include "LoKi/ILoKiSvc.h" -#include "LoKi/Kinematics.h" - -#include "Event/Particle.h" -#include "Event/RecVertex.h" -#include "Kernel/IParticleTransporter.h" - -#include "DetDesc/IDetectorElement.h" - -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/ITupleTool.h" - -namespace LoKi { - - /** @class SeedConeJetMaker - * - * The most trivial, SeedConeJet based implementaion of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Cedric POTTERAT cedric.potterat@cern.ch - * @date 2011-01-31 - */ - class SeedConeJetMaker : public virtual IJetMaker, public GaudiTool { - public: - /** The main method: jet-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* jetMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets jets ; - * - * // find the jets! - * StatusCode sc = jetMaker -> makeJets ( input , jets ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iJet = jets.begin() ; - * jets.end() != iJet ; ++iJet ) - * { - * // get the jet - * LHCb::Particle* jet = *iJet ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param jets container of output jets - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according to EPFL/UB/CERN studies. - * - */ - SeedConeJetMaker( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ) - // - , m_jetID( 10098 ) - // - , m_r( 0.7 ) - , m_ptmin( 8000 ) - // - , m_sort( 4 ) - , m_sortSeed( 4 ) - - // - , m_combinerName( "MomentumCombiner" ) - , m_combiner( 0 ) - - , m_seedFinderName( "LoKi__SeedFinder" ) - , m_seedFinder( 0 ) - , m_uniquetrk( false ) - - { - // - declareInterface<IJetMaker>( this ); - // - declareProperty( "JetID", m_jetID, "Particle ID for the Jet" ); - declareProperty( "RParameter", m_r ); - declareProperty( "PtMin", m_ptmin ); - // - declareProperty( "SortJet", m_sort, - "Sorting Criteria for jets [0:none,1:pt,2:E, 3:eta,4:DauThenPt, default:DauThenPt]" ); - - declareProperty( "SortSeed", m_sort, - "Sorting Criteria for Seeds [0:none,1:pt,2:E, 3:eta,4:DauThenPt, default:DauThenPt]" ); - // define momentum combiner - declareProperty( "ParticleCombiner", m_combinerName ); - declareProperty( "SeedFinder", m_seedFinderName, "VVSeedFinder, UBSeedFinder, ..." ); - - declareProperty( "ConeUniqueTrk", m_uniquetrk, "Cone With unique trk ON/OFF" ); - } - - /** standard initialization of the tool - * @return status code - */ - StatusCode initialize() override; - - protected: - /// make the detailed check of all parameters - - inline StatusCode check() const { - if ( 0 > m_ptmin ) { Warning( "PtMin is negative " ).ignore(); } - return StatusCode::SUCCESS; - } - - protected: - int to_user_index( const int index ) const { return index + 10000; } - int from_user_index( const int index ) const { return index - 10000; } - - protected: - Gaudi::Property<std::string> m_standardGeometry_address{this, "StandardGeometryTop", "/dd/Structure/LHCb"}; - - // proposed jet ID - int m_jetID; ///< proposed jet ID - // R-parameter - double m_r; ///< R-parameters - // ptMin-parameter - double m_ptmin; ///< pt-min parameter - // jet sorting criteria - int m_sort; ///< jet sorting criteria - int m_sortSeed; ///< seed sorting criteria - - // combiner - std::string m_combinerName; - mutable IParticleCombiner* m_combiner; ///< combiner to be used - - std::string m_seedFinderName; - mutable IJetMaker* m_seedFinder; ///< combiner to be used - - bool m_uniquetrk; - - ITupleTool* m_tuple; - - IJetMaker::Jets JetCone( double const& Rmax, IJetMaker::Jets& seeds, IJetMaker::Input const& inputs, - IGeometryInfo const& geometry ) const; - - IJetMaker::Jets JetConePurged( double const& Rmax, IJetMaker::Jets& seeds, IJetMaker::Input const& inputs, - IGeometryInfo const& geometry ) const; - - private: - mutable Gaudi::Accumulators::StatCounter<> m_nJets{this, "#jets"}; - }; -} // namespace LoKi - -class sortDauPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - if ( obj1->weight() == obj2->weight() ) - if ( obj1->daughtersVector().size() == obj2->daughtersVector().size() ) - return obj1->pt() > obj2->pt(); - else - return obj1->daughtersVector().size() > obj2->daughtersVector().size(); - else - return obj1->weight() > obj2->weight(); - } -}; - -class sortE { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().E() > obj2->momentum().E(); - } -}; - -class sortEta { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().Eta() > obj2->momentum().Eta(); - } -}; - -class sortPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { return obj1->pt() > obj2->pt(); } -}; - -/* standard initialization of the tool - * @return status code - */ -StatusCode LoKi::SeedConeJetMaker::initialize() { - StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) { return sc; } - - svc<LoKi::ILoKiSvc>( "LoKiSvc", true ); - - // select your prefered SeedFinder - m_seedFinder = tool<IJetMaker>( m_seedFinderName, this ); - - m_combiner = tool<IParticleCombiner>( m_combinerName, this ); - - return sc; -} - -// =========================================================================== -// find the jets via a Seed Finder -// =========================================================================== - -StatusCode LoKi::SeedConeJetMaker::makeJets( const IJetMaker::Input& input_, IJetMaker::Jets& jets_ ) const { - - const LHCb::RecVertex::Container* verts = get<LHCb::RecVertex::Container>( LHCb::RecVertexLocation::Primary ); - if ( verts->size() != 1 ) return StatusCode::SUCCESS; - LHCb::RecVertex::Container::const_iterator iv = verts->begin(); - - const LHCb::RecVertex vtx = **iv; - - makeJets( input_, vtx, jets_ ).ignore(); - - return StatusCode::SUCCESS; -} - -StatusCode LoKi::SeedConeJetMaker::makeJets( const IJetMaker::Input& input_, const LHCb::RecVertex& vtx_, - IJetMaker::Jets& jets_ ) const { - - IJetMaker::Jets Seeds; - - // create seeds with the input particles - m_seedFinder->makeJets( input_, vtx_, Seeds ).ignore(); - - if ( Seeds.empty() ) { - Warning( "No Seeds from: " + m_seedFinderName ).ignore(); - return StatusCode::SUCCESS; - } - - // sort the seeds according to the selected creteria (SortSeed opt) - if ( m_sortSeed == 1 ) std::stable_sort( Seeds.begin(), Seeds.end(), sortPt() ); - if ( m_sortSeed == 2 ) std::stable_sort( Seeds.begin(), Seeds.end(), sortE() ); - if ( m_sortSeed == 3 ) std::stable_sort( Seeds.begin(), Seeds.end(), sortEta() ); - if ( m_sortSeed == 4 ) std::stable_sort( Seeds.begin(), Seeds.end(), sortDauPt() ); - - LHCb::Particle::Vector output; - LHCb::Particle::Vector outputP; - - output.reserve( Seeds.size() ); - outputP.reserve( Seeds.size() ); - - // Get the default geometry FIXME, use functional framework - auto lhcb = getDet<IDetectorElement>( m_standardGeometry_address ); - if ( !lhcb ) { throw GaudiException( "Could not load geometry", name(), StatusCode::FAILURE ); } - auto& geometry = *lhcb->geometry(); - - // create jets arrount the seeds directions, use the RParameter opt, PtJetMin and the jetID. - output = JetCone( m_r, Seeds, input_, geometry ); - // sort the jets according to the selected creteria (SortJet opt) - if ( m_sort == 1 ) std::stable_sort( output.begin(), output.end(), sortPt() ); - if ( m_sort == 2 ) std::stable_sort( output.begin(), output.end(), sortE() ); - if ( m_sort == 3 ) std::stable_sort( output.begin(), output.end(), sortEta() ); - if ( m_sort == 4 ) std::stable_sort( output.begin(), output.end(), sortDauPt() ); - int myj = 0; - for ( LHCb::Particle::Vector::iterator ijet = output.begin(); output.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "ConeJets " << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - } - - if ( m_uniquetrk ) { - outputP = JetConePurged( m_r, output, input_, geometry ); - if ( m_sort == 1 ) std::stable_sort( outputP.begin(), outputP.end(), sortPt() ); - if ( m_sort == 2 ) std::stable_sort( outputP.begin(), outputP.end(), sortE() ); - if ( m_sort == 3 ) std::stable_sort( outputP.begin(), outputP.end(), sortEta() ); - if ( m_sort == 4 ) std::stable_sort( outputP.begin(), outputP.end(), sortDauPt() ); - - } else - outputP = output; - - myj = 0; - for ( LHCb::Particle::Vector::iterator ijet = outputP.begin(); outputP.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "ConeJetsP" << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - } - - if ( msgLevel( MSG::INFO ) ) { m_nJets += outputP.size(); } - - jets_ = outputP; - - // Delete the pointer to the seed - for ( IJetMaker::Jets::iterator ijet = Seeds.begin(); Seeds.end() != ijet; ++ijet ) delete *ijet; - - return StatusCode::SUCCESS; -} - -//============================================================================= -// JETCONE: THE HEART OF EVERYTHING, BUILDS A JET AROUND A SEED -//============================================================================= -IJetMaker::Jets LoKi::SeedConeJetMaker::JetCone( double const& Rmax, IJetMaker::Jets& seeds, - IJetMaker::Input const& inputs, IGeometryInfo const& geometry ) const { - - LHCb::Particle::Vector pJetVec; - IJetMaker::Jets::const_iterator iseeds; - - // loop over the seeds - for ( iseeds = seeds.begin(); iseeds != seeds.end(); ++iseeds ) { - - const LHCb::Particle* pSeed = *iseeds; - - LHCb::Particle pJet; - LHCb::Vertex vJet; - LHCb::Particle::ConstVector daughters; - daughters.reserve( inputs.size() ); - - pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); - pJet.setReferencePoint( pSeed->referencePoint() ); - pJet.setEndVertex( pSeed->endVertex() ); - double phiSeed = 0; - double etaSeed = 0; - - if ( pSeed->endVertex() != NULL ) { - Gaudi::XYZVector seedd = Gaudi::XYZVector( Gaudi::XYZPoint( pSeed->endVertex()->position() ) - - Gaudi::XYZPoint( pSeed->referencePoint() ) ); - phiSeed = seedd.phi(); - etaSeed = seedd.eta(); - - } else { - phiSeed = pSeed->momentum().phi(); - etaSeed = pSeed->momentum().eta(); - } - - const double PI = 3.14159265; - - Gaudi::LorentzVector ptot = Gaudi::LorentzVector( 0., 0., 0., 0. ); - - LHCb::Particle::ConstVector::const_iterator ip; - // loop over the particles - for ( ip = inputs.begin(); ip != inputs.end(); ++ip ) { - - const LHCb::Particle* myPart = *ip; - Gaudi::LorentzVector p1 = myPart->momentum(); - - double phiPart = myPart->momentum().phi(); - double etaPart = myPart->momentum().eta(); - double Dphi = fabs( phiSeed - phiPart ); - if ( Dphi > PI ) { Dphi = 2 * PI - Dphi; }; - double R = sqrt( Dphi * Dphi + ( etaSeed - etaPart ) * ( etaSeed - etaPart ) ); - - // continue if not in DeltaR (Cone geom in eta-phi plan) - if ( R > Rmax ) continue; - ptot = ptot + p1; - - daughters.push_back( myPart ); - } - - if ( daughters.empty() ) Warning( "Empty list of of daughter particles, skip it" ).ignore(); - - // continue if the jet is below PtJetMin opt - if ( ptot.Pt() < m_ptmin ) continue; - // save the particles involved as the daugther of the jet - StatusCode sc = m_combiner->combine( daughters, pJet, vJet, geometry ); - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - - pJet.setMomentum( Gaudi::LorentzVector( ptot ) ); - pJetVec.push_back( pJet.clone() ); - } - // return your jets - return pJetVec; -} - -//============================================================================= -// JETCONE: THE HEART OF EVERYTHING, BUILDS A JET AROUND A SEED -//============================================================================= -IJetMaker::Jets LoKi::SeedConeJetMaker::JetConePurged( double const& Rmax, IJetMaker::Jets& seeds, - IJetMaker::Input const& inputs, - IGeometryInfo const& geometry ) const { - - LHCb::Particle::Vector pJetVec; - IJetMaker::Jets::const_iterator iseeds; - - LHCb::Particle::Vector inputs2; - - LHCb::Particle::ConstVector::const_iterator ip; - // loop over the particles - for ( ip = inputs.begin(); ip != inputs.end(); ++ip ) inputs2.push_back( ( *ip )->clone() ); - - // loop over the seeds - for ( iseeds = seeds.begin(); iseeds != seeds.end(); ++iseeds ) { - - const LHCb::Particle* pSeed = *iseeds; - - LHCb::Particle pJet; - LHCb::Vertex vJet; - LHCb::Particle::ConstVector daughters; - daughters.reserve( inputs.size() ); - - pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); - pJet.setReferencePoint( pSeed->referencePoint() ); - pJet.setEndVertex( pSeed->endVertex() ); - double phiSeed = 0; - double etaSeed = 0; - - if ( pSeed->endVertex() != NULL ) { - Gaudi::XYZVector seedd = Gaudi::XYZVector( Gaudi::XYZPoint( pSeed->endVertex()->position() ) - - Gaudi::XYZPoint( pSeed->referencePoint() ) ); - phiSeed = seedd.phi(); - etaSeed = seedd.eta(); - - } else { - phiSeed = pSeed->momentum().phi(); - etaSeed = pSeed->momentum().eta(); - } - - const double PI = 3.14159265; - - Gaudi::LorentzVector ptot = Gaudi::LorentzVector( 0., 0., 0., 0. ); - - LHCb::Particle::Vector inputs3; // part not used! - - LHCb::Particle::Vector::iterator ip; - // loop over the particles - for ( ip = inputs2.begin(); ip != inputs2.end(); ++ip ) { - - LHCb::Particle* myPart = *ip; - Gaudi::LorentzVector p1 = myPart->momentum(); - - double phiPart = myPart->momentum().phi(); - double etaPart = myPart->momentum().eta(); - double Dphi = fabs( phiSeed - phiPart ); - if ( Dphi > PI ) { Dphi = 2 * PI - Dphi; }; - double R = sqrt( Dphi * Dphi + ( etaSeed - etaPart ) * ( etaSeed - etaPart ) ); - - // continue if not in DeltaR (Cone geom in eta-phi plan) - if ( R > Rmax ) { - inputs3.push_back( myPart ); - } else { - ptot = ptot + p1; - - daughters.push_back( myPart ); - } - } - inputs2 = inputs3; // remove the particle already assingned to a jet - - if ( daughters.empty() ) Warning( "Empty list of of daughter particles, skip it" ).ignore(); - - // continue if the jet is below PtJetMin opt - if ( ptot.Pt() < m_ptmin ) continue; - // save the particles involved as the daugther of the jet - StatusCode sc = m_combiner->combine( daughters, pJet, vJet, geometry ); - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - - pJet.setMomentum( Gaudi::LorentzVector( ptot ) ); - pJetVec.push_back( pJet.clone() ); - } - // return your jets - return pJetVec; -} - -DECLARE_COMPONENT( LoKi::SeedConeJetMaker ) diff --git a/Phys/LoKiJets/src/LoKiSeedFinder.cpp b/Phys/LoKiJets/src/LoKiSeedFinder.cpp deleted file mode 100644 index 80ea8be2643e8444fa5e98d0fff5fd59f85ceaf0..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiSeedFinder.cpp +++ /dev/null @@ -1,969 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" - -#include "LoKi/Geometry.h" -#include "LoKi/ILoKiSvc.h" -#include "LoKi/Kinematics.h" - -#include "Event/Particle.h" -#include "Event/RecVertex.h" -#include "Kernel/IDistanceCalculator.h" -#include "Kernel/IParticleTransporter.h" -#include "Kernel/IVertexFit.h" - -#include "DetDesc/IDetectorElement.h" - -#include "GaudiAlg/GaudiTool.h" - -namespace LoKi { - - /** @class VVSeedFinder - * - * The SeedFinder implementaion of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Cedric POTTERAT cedric.potterat@cern.ch - * @date 2011-01-31 - */ - class SeedFinder : public virtual IJetMaker, public GaudiTool { - public: - /** The main method: seed-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* seedMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets seeds ; - * - * // find the jets! - * StatusCode sc = seedMaker -> makeJets ( input , seeds ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iSeed = seeds.begin() ; - * seeds.end() != iSeed ; ++iSeed ) - * { - * // get the jet - * LHCb::Particle* seed = *iSeed ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param seeds container of output seeds (of type Particle) - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according to EPFL/UB/CERN studies. - * - */ - SeedFinder( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ) - , m_seedID( 20097 ) - , m_r( 0.15 ) - , m_sort( 4 ) - , m_combinerName( "MomentumCombiner" ) - , m_combiner( 0 ) - , m_dist( 0 ) - , m_fitter( 0 ) - , m_PtTrackMin( 600.0 ) - , m_PTrackMin( 1000.0 ) - , m_IPmin( 0.1 ) - , m_Signif( 2.5 ) - , m_DMK0( 10.0 ) - , m_DtrakMaxK0( 0.5 ) - , m_TseedVtxMin( 1.0 ) - , m_TseedVtxMax( 200.0 ) - , m_TseedVtxMinAnyPV( 0.1 ) - , m_DtrakMax( 0.5 ) - , m_PtSeedsMin( 1000 ) - , m_SeedsMaxChi2DoF( 50. ) - , m_Triplets( true ) - , m_DRmin( 0.1 ) - , m_DRmax( 50. ) - , m_TrkChi2DoF( 2.5 ) - , m_PVveto( true ) - // , m_dirAng ( 5.) - - { - declareInterface<IJetMaker>( this ); - declareProperty( "SeedID", m_seedID, "Particle ID for the Seed" ); - declareProperty( "SeedRParameter", m_r ); - declareProperty( "Sort", m_sort, "Sorting Criteria for jets [0:none,1:pt,2:E,3:eta, default:Pt]" ); - declareProperty( "ParticleCombiner", m_combinerName ); - declareProperty( "SeedPtTrackMin", m_PtTrackMin, "pt of the track used for Vertexing" ); - declareProperty( "SeedPTrackMin", m_PTrackMin, "p of the track used for Vertexing" ); - declareProperty( "SeedIPmin", m_IPmin, "ip of the track used for Vertexing" ); - declareProperty( "SeedSignif", m_Signif, "signif oft he track used for Vertexing" ); - declareProperty( "SeedDMK0", m_DMK0, "mass window for Ks" ); - declareProperty( "SeedDtrakMaxK0", m_DtrakMaxK0, "dca window for Ks" ); - declareProperty( "SeedTseedVtxMin", m_TseedVtxMin, "min distance btw PV and the vtx" ); - declareProperty( "SeedTseedVtxMinAnyPV", m_TseedVtxMinAnyPV, "min distance btw any PV and the vtx" ); - declareProperty( "SeedTseedVtxMax", m_TseedVtxMax, "max distance btw PV and the vtx" ); - declareProperty( "SeedDtrakMax", m_DtrakMax, "dca window for vtx" ); - declareProperty( "SeedPtSeedsMin", m_PtSeedsMin, "min pt of the seeds" ); - declareProperty( "SeedVtxMaxChi2PerDoF", m_SeedsMaxChi2DoF, "max chi2 per dof for the vtx fit of the seed" ); - declareProperty( "SeedTriplets", m_Triplets, "built Vtx with 3 tracks" ); - declareProperty( "SeedDRmin", m_DRmin, "min positon in R of the vtx" ); - declareProperty( "SeedDRmax", m_DRmax, "max positon in R of the vtx" ); - declareProperty( "SeedTrkChi2PerDoF", m_TrkChi2DoF, "max chi2PerDoF for the track used for the vtx" ); - declareProperty( "vetoPV", m_PVveto, - "exclude vertex near to any PV with distance fixe by 'SeedTseedVtxMinAnyPV' and exclude the trk " - "associated to any PV to construct a secondary vtx" ); - // declareProperty ( "SeedMaxDirAngle" , m_dirAng , "Max direction angle btw the seed momentum - // and the direction vtx-> seed" ); - } - /** standard initialization of the tool - * @return status code - */ - StatusCode initialize() override; - - protected: - /// make the detailed check of all parameters - - inline StatusCode check() const { - if ( 0 > m_ptmin ) { Warning( "PtMin is negative " ).ignore(); } - return StatusCode::SUCCESS; - } - - protected: - Gaudi::Property<std::string> m_standardGeometry_address{this, "StandardGeometryTop", "/dd/Structure/LHCb"}; - - int to_user_index( const int index ) const { return index + 10000; } - int from_user_index( const int index ) const { return index - 10000; } - // proposed jet ID - int m_seedID; ///< proposed jet ID - // R-parameter - double m_r; ///< R-parameters - // ptMin-parameter - double m_ptmin; ///< pt-min parameter - // jet sorting criteria - int m_sort; ///< jet sorting criteria - // combiner - std::string m_combinerName; - mutable IParticleCombiner* m_combiner; ///< combiner to be used - - IDistanceCalculator* m_dist; - IVertexFit* m_fitter; - - double m_Rmax; - double m_PtTrackMin; - double m_PTrackMin; - double m_IPmin; - double m_Signif; - double m_DMK0; - double m_DtrakMaxK0; - double m_TseedVtxMin; - double m_TseedVtxMax; - double m_TseedVtxMinAnyPV; - double m_DtrakMax; - double m_PtSeedsMin; - double m_SeedsMaxChi2DoF; - bool m_Triplets; - double m_DRmin; - double m_DRmax; - double m_TrkChi2DoF; - bool m_PVveto; - double m_dirAng; - - void RemoveTracks( LHCb::Particle::ConstVector& particles, const LHCb::RecVertex PV ) const; - - double getDeltaR( LHCb::Particle* p1, LHCb::Particle* p2 ) const; - - private: - mutable Gaudi::Accumulators::StatCounter<> m_nSeeds{this, "#seeds"}; - }; -} // namespace LoKi - -class sortDauPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - if ( obj1->weight() == obj2->weight() ) - if ( obj1->daughtersVector().size() == obj2->daughtersVector().size() ) - return obj1->pt() > obj2->pt(); - else - return obj1->daughtersVector().size() > obj2->daughtersVector().size(); - else - return obj1->weight() > obj2->weight(); - } -}; - -class sortCHi2NDoF { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - - if ( obj1->daughtersVector().size() == obj2->daughtersVector().size() ) - return obj1->endVertex()->chi2PerDoF() < obj2->endVertex()->chi2PerDoF(); - else - return obj1->daughtersVector().size() > obj2->daughtersVector().size(); - } -}; - -class sortDirAng { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - - Gaudi::XYZVector seedd1 = Gaudi::XYZVector( Gaudi::XYZPoint( obj1->endVertex()->position() ) - - Gaudi::XYZPoint( obj1->referencePoint() ) ); - Gaudi::XYZVector seedd2 = Gaudi::XYZVector( Gaudi::XYZPoint( obj2->endVertex()->position() ) - - Gaudi::XYZPoint( obj2->referencePoint() ) ); - - double dotprod1 = std::fabs( std::acos( seedd1.Unit().Dot( obj1->momentum().Vect().Unit() ) ) ); - double dotprod2 = std::fabs( std::acos( seedd2.Unit().Dot( obj2->momentum().Vect().Unit() ) ) ); - - return dotprod1 < dotprod2; - } -}; - -class sortCOMBO { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - - Gaudi::XYZVector seedd1 = Gaudi::XYZVector( Gaudi::XYZPoint( obj1->endVertex()->position() ) - - Gaudi::XYZPoint( obj1->referencePoint() ) ); - Gaudi::XYZVector seedd2 = Gaudi::XYZVector( Gaudi::XYZPoint( obj2->endVertex()->position() ) - - Gaudi::XYZPoint( obj2->referencePoint() ) ); - - double dotprod1 = std::fabs( std::acos( seedd1.Unit().Dot( obj1->momentum().Vect().Unit() ) ) ); - double dotprod2 = std::fabs( std::acos( seedd2.Unit().Dot( obj2->momentum().Vect().Unit() ) ) ); - - double combo1 = 1000 * 1000 * ( dotprod1 * obj1->endVertex()->chi2PerDoF() ) / - ( obj1->daughtersVector().size() * obj1->pt() * obj1->pt() ); - double combo2 = 1000 * 1000 * ( dotprod2 * obj2->endVertex()->chi2PerDoF() ) / - ( obj2->daughtersVector().size() * obj2->pt() * obj2->pt() ); - - return combo1 < combo2; - } -}; - -class sortE { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().E() > obj2->momentum().E(); - } -}; - -class sortEta { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().Eta() > obj2->momentum().Eta(); - } -}; - -class sortPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { return obj1->pt() > obj2->pt(); } -}; - -/* standard initialization of the tool - * @return status code - */ -StatusCode LoKi::SeedFinder::initialize() { - const StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) { return sc; } - - svc<LoKi::ILoKiSvc>( "LoKiSvc", true ); - - m_dist = tool<IDistanceCalculator>( "LoKi::DistanceCalculator" ); - - m_combiner = tool<IParticleCombiner>( m_combinerName, this ); - - m_fitter = tool<IVertexFit>( "OfflineVertexFitter" ); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << " ------------------------------------------------------ " << endmsg; - debug() << " | | " << endmsg; - debug() << " | SeedsFinder Tool for Phys Particles | " << endmsg; - debug() << " | | " << endmsg; - debug() << " ------------------------------------------------------ " << endmsg; - debug() << " | | " << endmsg; - debug() << " | " << endmsg; - debug() << " | SeedsFinder parameters: " << endmsg; - debug() << " | SeedID " << m_seedID << endmsg; - debug() << " | SeedTriplets " << m_Triplets << endmsg; - debug() << " | PtTrackMin " << m_PtTrackMin << endmsg; - debug() << " | PTrackMin " << m_PTrackMin << endmsg; - debug() << " | TrkChi2PerDoF " << m_TrkChi2DoF << endmsg; - debug() << " | IPmin " << m_IPmin << endmsg; - debug() << " | Signif " << m_Signif << endmsg; - debug() << " | DMK0 " << m_DMK0 << endmsg; - debug() << " | TseedVtxMin " << m_TseedVtxMin << endmsg; - debug() << " | TseedVtxMax " << m_TseedVtxMax << endmsg; - debug() << " | DtrakMax " << m_DtrakMax << endmsg; - debug() << " | PtSeedsMin " << m_PtSeedsMin << endmsg; - debug() << " | SeedsVtxMaxChi2PerDoF " << m_SeedsMaxChi2DoF << endmsg; - debug() << " | " << endmsg; - debug() << " | | " << endmsg; - debug() << " ------------------------------------------------------ " << endmsg; - } - - return sc; -} - -// =========================================================================== -// find seeds -// =========================================================================== - -StatusCode LoKi::SeedFinder::makeJets( const IJetMaker::Input& input_, IJetMaker::Jets& jets_ ) const { - - const LHCb::RecVertex::Container* verts = get<LHCb::RecVertex::Container>( LHCb::RecVertexLocation::Primary ); - if ( verts->size() != 1 ) return StatusCode::SUCCESS; - LHCb::RecVertex::Container::const_iterator iv = verts->begin(); - - const LHCb::RecVertex vtx = **iv; - - makeJets( input_, vtx, jets_ ).ignore(); - - return StatusCode::SUCCESS; -} - -StatusCode LoKi::SeedFinder::makeJets( const IJetMaker::Input& input_, const LHCb::RecVertex& RecVert_, - IJetMaker::Jets& jets_ ) const { - - const LHCb::RecVertex::Container* verts = get<LHCb::RecVertex::Container>( LHCb::RecVertexLocation::Primary ); - LHCb::RecVertex::Container::const_iterator iv; - LHCb::VertexBase* RecVert = RecVert_.clone(); - - std::vector<LHCb::VertexBase*> PVs; - - for ( iv = verts->begin(); iv != verts->end(); iv++ ) - if ( ( *iv )->position() != RecVert->position() ) PVs.push_back( ( *iv )->clone() ); - - IJetMaker::Jets Seeds; - IJetMaker::Jets SeedsPurged; - - double ipl, chi2l; - StatusCode sc; - LHCb::Vertex vtx, vseedbest; - double dca12, dcae12; - // Kshort Mass - const double MK0 = 497.7; - - std::vector<Gaudi::XYZVector> Slopes_flag; - - LHCb::Particle::ConstVector PartIPK0Sub; - LHCb::Particle::ConstVector PartIP; - LHCb::Particle::ConstVector all; - LHCb::Particle::ConstVector KsDau; - - double tof = 0, tofe = 0; - - int testLiHi = 0; - int testClone = 0; - int testLong = 0; - int testCharged = 0; - - // Get the default geometry FIXME, use functional framework - auto lhcb = getDet<IDetectorElement>( m_standardGeometry_address ); - if ( !lhcb ) { throw GaudiException( "Could not load geometry", name(), StatusCode::FAILURE ); } - auto& geometry = *lhcb->geometry(); - - //----------------------------------------// - // test the particles // - //----------------------------------------// - - for ( LHCb::Particle::ConstVector::const_iterator ip = input_.begin(); input_.end() != ip; ++ip ) { - const LHCb::Particle* p = *ip; - if ( 0 == p ) { - Warning( "Invalid input particle" ).ignore(); - continue; - } - const LHCb::ProtoParticle* myProto = p->proto(); - if ( myProto != NULL ) { - const LHCb::Track* trk = myProto->track(); - if ( trk != NULL ) { - if ( trk->flags() == 1 || trk->flags() == 4 ) continue; - // == 1 : backward, == 4 clone - double CloneDist = trk->info( LHCb::Track::AdditionalInfo::CloneDist, 9999. ); - double TrkLi = trk->likelihood(); - double GhostProb = trk->ghostProbability(); - double LongTrk = static_cast<int>( trk->type() ); - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << "addinfo Likelihood = " << TrkLi << endmsg; - verbose() << "addinfo CloneDist = " << CloneDist << endmsg; - verbose() << "addinfo GhostProb = " << GhostProb << endmsg; - verbose() << "addinfo Longb = " << LongTrk << endmsg; - verbose() << "p Trk = " << trk->momentum() / Gaudi::Units::GeV << endmsg; - } - if ( CloneDist < 9999. ) testClone++; - if ( TrkLi < -60. ) testLiHi++; - if ( LongTrk == 3 ) testLong++; - //----------------------------------------// - // discard clones and Ghost // - //----------------------------------------// - if ( CloneDist < 9999. ) continue; - if ( TrkLi < -60. ) continue; - } - } - - bool m_FilterPart = true; - //----------------------------------------// - // check if 1 particle is not twice on // - // the TES with the slopes // - //----------------------------------------// - bool flag_same = false; - if ( m_FilterPart ) { - if ( Slopes_flag.size() != 0 ) { - for ( int k = 0; k < (int)Slopes_flag.size(); k++ ) { - if ( Slopes_flag.at( k ) == p->slopes() ) { - flag_same = true; - if ( msgLevel( MSG::VERBOSE ) ) - verbose() << "same slopes !! s1: " << Slopes_flag.at( k ) << " s2: " << p->slopes() << endmsg; - break; - } - } - } - if ( !flag_same ) { - Slopes_flag.push_back( p->slopes() ); - } else { - continue; - } - } - - //----------------------------------------// - // save all the part that arrived here // - // as the global container // - //----------------------------------------// - all.push_back( p ); - - if ( p->charge() == 0 ) continue; - testCharged++; - if ( p->proto() == NULL ) continue; - if ( p->proto()->track() == NULL ) continue; - if ( p->proto()->track()->type() != LHCb::Track::Types::Long ) continue; - if ( p->proto()->track()->chi2PerDoF() > m_TrkChi2DoF ) continue; - - sc = m_dist->distance( p, RecVert, ipl, chi2l, geometry ); - - if ( p->p() < m_PTrackMin ) continue; - if ( p->pt() < m_PtTrackMin ) continue; - if ( ipl < m_IPmin ) continue; - if ( sqrt( chi2l ) < m_Signif ) continue; - - //----------------------------------------// - // save all the charge, with a good Pt // - // and a good IP/IPE as the input // - // particles for the Cone Jet Algo // - //----------------------------------------// - PartIP.push_back( p ); - } - - if ( m_PVveto ) { - for ( iv = verts->begin(); iv != verts->end(); ++iv ) { RemoveTracks( PartIP, **iv ); } - } - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "nb of ghost (LiH) : " << testLiHi << endmsg; - debug() << "nb of clone : " << testClone << endmsg; - debug() << "nb of long : " << testLong << endmsg; - debug() << "nb of charged : " << testCharged << endmsg; - debug() << "nb of 'photon' : " << (int)all.size() - double( testCharged ) << endmsg; - debug() << "Particle ALL size : " << input_.size() << endmsg; - debug() << "Particle ALL size (JetCone) : " << all.size() << endmsg; - debug() << "Particle INPUT for JET : " << PartIP.size() << endmsg; - } - - int cntK0 = 0; - - //============================================================================= - // COMBINEPIPIMINK0: GET RIF OF K0S - //============================================================================= - for ( LHCb::Particle::ConstVector::const_iterator ip = PartIP.begin(); PartIP.end() != ip; ++ip ) { - const LHCb::Particle* p = *ip; - bool testK = false; - - for ( LHCb::Particle::ConstVector::const_iterator kp = KsDau.begin(); KsDau.end() != kp; ++kp ) - if ( *ip == *kp ) { - testK = true; - break; - } - - if ( !testK ) - for ( LHCb::Particle::ConstVector::const_iterator ip2 = ip + 1; PartIP.end() != ip2; ++ip2 ) { - const LHCb::Particle* p2 = *ip2; - - sc = m_dist->distance( ( p ), ( p2 ), dca12, dcae12, geometry ); - if ( !sc ) { warning() << "can't mesure the dist " << endmsg; } - if ( dca12 > m_DtrakMax ) continue; // dca too large bwt the tracks - - Gaudi::LorentzVector sum = ( p )->momentum() + ( p2 )->momentum(); - - //============================================================================= - // COMBINEPIPIMINK0: GET RIF OF K0S - //============================================================================= - if ( ( p )->particleID().abspid() == 211 && ( p2 )->particleID().abspid() == 211 && - ( ( p )->charge() ) * ( ( p2 )->charge() ) < 0 ) - if ( fabs( sum.M() - MK0 ) < m_DMK0 ) { - cntK0++; - testK = true; - - KsDau.push_back( p ); - KsDau.push_back( p2 ); - } - } - - if ( testK ) continue; - PartIPK0Sub.push_back( p ); - } - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "nb K0 : " << cntK0 << endmsg; - debug() << "Particle INPUT for JET no K0 : " << PartIPK0Sub.size() << endmsg; - } - - if ( PartIPK0Sub.size() < 2 ) { return Warning( "Not enough good part for seeding", StatusCode::SUCCESS ); } - - Gaudi::XYZPoint BL_P = Gaudi::XYZPoint( 0, 0, 0 ); - Gaudi::LorentzVector BL_M = Gaudi::LorentzVector( 0, 0, 1, 0 ); - - if ( exist<LHCb::Particle::Range>( "/Event/BeamLine" ) ) { - const LHCb::Particle::Range BL = get<LHCb::Particle::Range>( "/Event/BeamLine" ); - const LHCb::Particle* tmp = *( BL.begin() ); - BL_P = Gaudi::XYZPoint( tmp->referencePoint() ); - BL_M = Gaudi::LorentzVector( tmp->momentum() ); - // m_BeamLine->setMomentum( tmp->momentum() ); - if ( msgLevel( MSG::DEBUG ) ) debug() << "Beam line position " << BL_P << " direction " << BL_M << endmsg; - } else { - if ( msgLevel( MSG::DEBUG ) ) - debug() << "No Beam line found at " - << "/Event/BeamLine" << endmsg; - - BL_P = Gaudi::XYZPoint( RecVert->position() ); - - if ( msgLevel( MSG::DEBUG ) ) debug() << "Beam line position " << BL_P << " direction " << BL_M << endmsg; - } - - // COMBIN PART - - LHCb::Particle::ConstVector::const_iterator jp, kp, lp; - for ( jp = PartIPK0Sub.begin(); jp != PartIPK0Sub.end(); jp++ ) { - for ( kp = jp + 1; kp != PartIPK0Sub.end(); kp++ ) { - - sc = m_dist->distance( ( *jp ), ( *kp ), dca12, dcae12, geometry ); - if ( !sc ) { warning() << "can't mesure the dist " << endmsg; } - if ( dca12 > m_DtrakMax ) continue; // dca too large bwt the tracks - - Gaudi::LorentzVector sum = ( *jp )->momentum() + ( *kp )->momentum(); - - if ( m_Triplets ) { - - // add Try 3-tracks Seed otpion - for ( lp = ( kp + 1 ); lp != PartIPK0Sub.end(); lp++ ) { - - double dca13, dcae13; - m_dist->distance( ( *jp ), ( *lp ), dca13, dcae13, geometry ).ignore(); - if ( dca13 > m_DtrakMax ) continue; // dca too large btw the trks - - double dca23, dcae23; - m_dist->distance( ( *lp ), ( *kp ), dca23, dcae23, geometry ).ignore(); - if ( dca23 > m_DtrakMax ) continue; // dca too large btw the trks - - StatusCode scvtx = m_fitter->fit( vtx, **jp, **kp, **lp, geometry ); - if ( scvtx.isFailure() ) { - warning() << "VTX Fit failed" << endmsg; - continue; - } - - if ( vtx.chi2PerDoF() > m_SeedsMaxChi2DoF ) continue; - // cut on chi2 - double mydz = -1; - if ( scvtx ) mydz = vtx.position().z() - RecVert->position().z(); - if ( mydz < 0 ) continue; - // the vtx is before the PV - - // intersection of the beam line with the XY plane, - // find the lambda parameter of the line. - double lambda = ( vtx.position().z() - BL_P.z() ) / BL_M.z(); - - // find x and y of intersection point - double x = BL_P.x() + lambda * BL_M.x(); - double y = BL_P.y() + lambda * BL_M.y(); - - x -= vtx.position().x(); - y -= vtx.position().y(); - - double dr = sqrt( x * x + y * y ); - if ( dr < m_DRmin || dr > m_DRmax ) continue; - - // the vtx is too close or to far from the beam direction - - //--------------------------- - // Seed level cuts ----------- - //--------------------------- - LHCb::VertexBase* PPvtx2 = vtx.clone(); - m_dist->distance( PPvtx2, RecVert, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMin || tof > m_TseedVtxMax ) continue; - // the vtx is too close or to far from the PV (~time of flight) - bool PVveto = false; - if ( m_PVveto ) - for ( std::vector<LHCb::VertexBase*>::iterator pv = PVs.begin(); pv != PVs.end(); pv++ ) { - m_dist->distance( PPvtx2, *pv, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMinAnyPV ) { - PVveto = true; - break; - } - } - if ( PVveto ) continue; - - Gaudi::LorentzVector sum2 = sum + ( *lp )->momentum(); - if ( sum2.Pt() < m_PtSeedsMin ) continue; - // Pt too soft - - LHCb::Particle::ConstVector daughters2; - LHCb::Particle pSeed2; - LHCb::Vertex vSeed2; - - daughters2.push_back( *lp ); - daughters2.push_back( *kp ); - daughters2.push_back( *jp ); - pSeed2.setParticleID( LHCb::ParticleID( m_seedID ) ); - pSeed2.setReferencePoint( RecVert->position() ); - pSeed2.setMomentum( sum2 ); - - StatusCode sc = m_combiner->combine( daughters2, pSeed2, vSeed2, geometry ); - - pSeed2.setEndVertex( vtx.clone() ); - // save the trks as the daugthers - // remove - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - Seeds.push_back( pSeed2.clone() ); - // remove lp to test another part - - } // end loop lp - } - - // 2-seed - LHCb::Particle pSeed; - LHCb::Particle::ConstVector daughters; - LHCb::Vertex vSeed; - - StatusCode scvtx = m_fitter->fit( vtx, **jp, **kp, geometry ); - - if ( scvtx.isFailure() ) { - warning() << "VTX Fit failed" << endmsg; - continue; - } - - if ( vtx.chi2PerDoF() > m_SeedsMaxChi2DoF ) continue; - - double mydz = -1; - if ( scvtx ) mydz = vtx.position().z() - RecVert->position().z(); - if ( mydz < 0 ) continue; - // the vtx is before the PV - - // Gaudi::XYZPoint dv = Gaudi::XYZPoint(vtx.position() - m_BeamLine->referencePoint().position()); - - // intersection of the beam line with the XY plane, - // find the lambda parameter of the line. - double lambda = ( vtx.position().z() - BL_P.z() ) / BL_M.z(); - - // find x and y of intersection point - double x = BL_P.x() + lambda * BL_M.x(); - double y = BL_P.y() + lambda * BL_M.y(); - - x -= vtx.position().x(); - y -= vtx.position().y(); - - double dr = sqrt( x * x + y * y ); - if ( dr < m_DRmin || dr > m_DRmax ) continue; - // the vtx is too close or to far from the beam direction - - //--------------------------- - // Seed level cuts ----------- - //--------------------------- - LHCb::VertexBase* PPvtx2 = vtx.clone(); - m_dist->distance( PPvtx2, RecVert, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMin || tof > m_TseedVtxMax ) continue; - // the vtx is too close or to far from the PV (~time of flight) - - bool PVveto = false; - if ( m_PVveto ) - for ( std::vector<LHCb::VertexBase*>::iterator pv = PVs.begin(); pv != PVs.end(); pv++ ) { - m_dist->distance( PPvtx2, *pv, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMinAnyPV ) { - PVveto = true; - break; - } - } - if ( PVveto ) continue; - - if ( sum.Pt() < m_PtSeedsMin ) continue; - // Pt too soft - - daughters.push_back( *jp ); - daughters.push_back( *kp ); - - pSeed.setParticleID( LHCb::ParticleID( m_seedID ) ); - pSeed.setReferencePoint( RecVert->position() ); - pSeed.setMomentum( sum ); - StatusCode sc = m_combiner->combine( daughters, pSeed, vSeed, geometry ); - // save the trks as the daugthers - pSeed.setEndVertex( vtx.clone() ); - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - Seeds.push_back( pSeed.clone() ); - - } // end loop kp - } // end loop jp - - if ( Seeds.empty() ) { - Warning( "No Seeds from LoKi::SeedFinder" ).ignore(); - return StatusCode::SUCCESS; - } - - //============================================================================= - // SORT THE SEEDS - //============================================================================= - - if ( m_sort == 1 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortPt() ); - else if ( m_sort == 2 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortE() ); - else if ( m_sort == 3 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortEta() ); - else if ( m_sort == 4 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortDauPt() ); - else if ( m_sort == 5 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortCHi2NDoF() ); - else if ( m_sort == 6 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortDirAng() ); - else if ( m_sort == 7 ) - std::stable_sort( Seeds.begin(), Seeds.end(), sortCOMBO() ); - - // remove seed who share a track with a higher racked seed - int myj = 0; - - for ( IJetMaker::Jets::iterator ijet = Seeds.begin(); Seeds.end() != ijet; ++ijet ) { - myj++; - - bool samedau = false; - LHCb::Particle* Seed1 = *ijet; - if ( Seed1->daughters().size() < 2 ) continue; - - const LHCb::Particle* dau11( 0 ); - const LHCb::Particle* dau12( 0 ); - const LHCb::Particle* dau13( 0 ); - - dau11 = Seed1->daughters().at( 0 ).target(); - dau12 = Seed1->daughters().at( 1 ).target(); - - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << "Seeds " << myj << "| nb of combin Seeds: " << Seed1->weight() - << " nb part in the jet: " << Seed1->daughtersVector().size() << " Pt: " << Seed1->pt() << endmsg; - verbose() << " | m : " << Seed1->momentum() << endmsg; - verbose() << " | eta: " << Seed1->momentum().eta() << " / phi: " << Seed1->momentum().phi() << endmsg; - verbose() << " |dau1: " << dau11->pt() << endmsg; - verbose() << " |dau2: " << dau12->pt() << endmsg; - } - - if ( Seed1->daughters().size() > 2 ) { - dau13 = Seed1->daughters().at( 2 ).target(); - if ( msgLevel( MSG::VERBOSE ) ) verbose() << " |dau3: " << dau13->pt() << endmsg; - } - - Gaudi::XYZVector seedd = Gaudi::XYZVector( Gaudi::XYZPoint( Seed1->endVertex()->position() ) - - Gaudi::XYZPoint( Seed1->referencePoint() ) ); - - double dotprod = seedd.Unit().Dot( Seed1->momentum().Vect().Unit() ); - - double a = std::fabs( std::acos( dotprod ) ); - - double combo1 = 1000 * 1000 * ( a * Seed1->endVertex()->chi2PerDoF() ) / - ( Seed1->daughtersVector().size() * Seed1->pt() * Seed1->pt() ); - - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << " |PAng: " << a << endmsg; - verbose() << " |Chi2: " << Seed1->endVertex()->chi2() << endmsg; - verbose() << " |nDoF: " << Seed1->endVertex()->nDoF() << endmsg; - verbose() << " |C/N : " << Seed1->endVertex()->chi2PerDoF() << endmsg; - verbose() << " |COMB: " << combo1 << endmsg; - } - - for ( IJetMaker::Jets::iterator ijet2 = SeedsPurged.begin(); SeedsPurged.end() != ijet2; ++ijet2 ) { - - LHCb::Particle* Seed2 = *ijet2; - - const LHCb::Particle* dau21( 0 ); - const LHCb::Particle* dau22( 0 ); - const LHCb::Particle* dau23( 0 ); - - dau21 = Seed2->daughters().at( 0 ).target(); - dau22 = Seed2->daughters().at( 1 ).target(); - if ( Seed2->daughters().size() > 2 ) dau23 = Seed2->daughters().at( 2 ).target(); - - // double R = getDeltaR(Seed1, Seed2); - - // if(R < 0.5){//start the merging process - - if ( ( dau11 == dau21 ) || ( dau12 == dau22 ) || ( dau11 == dau22 ) || ( dau12 == dau21 ) ) { - samedau = true; - break; - } - if ( Seed1->daughters().size() > 2 ) - if ( ( dau13 == dau21 ) || ( dau13 == dau22 ) ) { - samedau = true; - break; - } - if ( Seed2->daughters().size() > 2 ) - if ( ( dau11 == dau23 ) || ( dau12 == dau23 ) ) { - samedau = true; - break; - } - if ( Seed2->daughters().size() > 2 && Seed1->daughters().size() > 2 ) - if ( dau13 == dau23 ) { - samedau = true; - break; - } - - // } - } - if ( samedau ) continue; - - SeedsPurged.push_back( Seed1 ); - } - - myj = 0; - for ( IJetMaker::Jets::iterator ijet = SeedsPurged.begin(); SeedsPurged.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "SeedsP " << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - - const LHCb::Particle* dau11( 0 ); - const LHCb::Particle* dau12( 0 ); - dau11 = pJet->daughters().at( 0 ).target(); - dau12 = pJet->daughters().at( 1 ).target(); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << " |dau1: " << dau11->pt() << endmsg; - debug() << " |dau2: " << dau12->pt() << endmsg; - } - - if ( pJet->daughters().size() > 2 ) { - const LHCb::Particle* dau13( 0 ); - dau13 = pJet->daughters().at( 2 ).target(); - if ( msgLevel( MSG::DEBUG ) ) debug() << " |dau3: " << dau13->pt() << endmsg; - } - - Gaudi::XYZVector seedd = Gaudi::XYZVector( Gaudi::XYZPoint( pJet->endVertex()->position() ) - - Gaudi::XYZPoint( pJet->referencePoint() ) ); - - double dotprod = seedd.Unit().Dot( pJet->momentum().Vect().Unit() ); - - double a = std::fabs( std::acos( dotprod ) ); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << " |dirA: " << a << endmsg; - debug() << " |Chi2: " << pJet->endVertex()->chi2() << endmsg; - debug() << " |nDoF: " << pJet->endVertex()->nDoF() << endmsg; - debug() << " |C/N : " << pJet->endVertex()->chi2PerDoF() << endmsg; - } - - double combo1 = 1000 * 1000 * ( a * pJet->endVertex()->chi2PerDoF() ) / - ( pJet->daughtersVector().size() * pJet->pt() * pJet->pt() ); - - if ( msgLevel( MSG::DEBUG ) ) debug() << " |COMB: " << combo1 << endmsg; - } - - if ( msgLevel( MSG::INFO ) ) { m_nSeeds += SeedsPurged.size(); } - - jets_ = SeedsPurged; - - return StatusCode::SUCCESS; -} - -//============================================================================= -// Remove trks form a PV -//============================================================================= - -void LoKi::SeedFinder::RemoveTracks( LHCb::Particle::ConstVector& particles, const LHCb::RecVertex PV ) const { - - // Remove all tracks from the PV - LHCb::Particle::ConstVector tmp; - SmartRefVector<LHCb::Track>::const_iterator iPV = PV.tracks().begin(); - int endkey = PV.tracks().back()->key(); - - for ( LHCb::Particle::ConstVector::const_iterator i = particles.begin(); i != particles.end(); ++i ) { - if ( ( *i )->proto()->track() == NULL ) continue; - const LHCb::Track* tk = ( *i )->proto()->track(); - - while ( ( ( *iPV )->key() < tk->key() ) && ( *iPV )->key() != endkey ) { iPV++; } - if ( ( *iPV )->key() == tk->key() ) { - if ( ( *iPV )->key() != endkey ) iPV++; - continue; - } - tmp.push_back( *i ); - } - // Copy back tmp to particles - particles = tmp; -} - -//============================================================================= -// GETDELTAR: CALCULATES DELTAR -//============================================================================= -double LoKi::SeedFinder::getDeltaR( LHCb::Particle* p1, LHCb::Particle* p2 ) const { - - const double PI = 3.14159265; - double R, Dphi, phi1, phi2, e1, e2; - phi1 = p1->momentum().phi(); - phi2 = p2->momentum().phi(); - e1 = p1->momentum().eta(); - e2 = p2->momentum().eta(); - Dphi = fabs( phi1 - phi2 ); - if ( Dphi > PI ) { Dphi = 2 * PI - Dphi; }; - R = sqrt( Dphi * Dphi + ( e1 - e2 ) * ( e1 - e2 ) ); - return R; -} - -DECLARE_COMPONENT( LoKi::SeedFinder ) diff --git a/Phys/LoKiJets/src/LoKiVVSeedFinder.cpp b/Phys/LoKiJets/src/LoKiVVSeedFinder.cpp deleted file mode 100644 index a92be7e96459aa6db27b9003e3dc59a65f503644..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiVVSeedFinder.cpp +++ /dev/null @@ -1,1118 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ - -#include "LoKi/Geometry.h" -#include "LoKi/ILoKiSvc.h" -#include "LoKi/Kinematics.h" - -#include "Kernel/IDistanceCalculator.h" -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" -#include "Kernel/IVertexFit.h" - -#include "Event/Particle.h" -#include "Event/RecVertex.h" -#include "Kernel/IParticleTransporter.h" -#include "LHCbMath/LHCbMath.h" - -#include "DetDesc/IDetectorElement.h" - -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/ITupleTool.h" - -namespace LoKi { - - /** @class VVSeedFinder - * - * The VVSeedFinder (EPFL) implementaion of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Cedric POTTERAT cedric.potterat@cern.ch - * @date 2011-01-31 - */ - class VVSeedFinder : public virtual IJetMaker, public GaudiTool { - public: - /** The main method: seed-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* seedMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets seeds ; - * - * // find the jets! - * StatusCode sc = seedMaker -> makeJets ( input , seeds ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iSeed = seeds.begin() ; - * seeds.end() != iSeed ; ++iSeed ) - * { - * // get the jet - * LHCb::Particle* seed = *iSeed ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param seeds container of output seeds (of type Particle) - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according to EPFL studies. - * @link http://cdsweb.cern.ch/record/1211038/files/LHCb-INT-2009-023.pdf LHCb-INT-2009-023 - Bay/Potterat - * @endlink - * @link http://cdsweb.cern.ch/record/1266883/files/CERN-THESIS-2010-074.pdf CERN_THESIS-2010-074 - Potterat - * @endlink - * - */ - VVSeedFinder( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ) - // - , m_seedID( 20098 ) - // - , m_r( 0.15 ) - // - , m_sort( 4 ) - // - , m_combinerName( "MomentumCombiner" ) - , m_combiner( 0 ) - , m_dist( 0 ) - , m_fitter( 0 ) - , m_PtTrackMin( 600.0 ) - , m_PTrackMin( 1000.0 ) - , m_IPmin( 0.1 ) - , m_Signif( 2.5 ) - , m_DMK0( 10.0 ) - , m_TseedVtxMin( 1.0 ) - , m_TseedVtxMax( 200.0 ) - , m_TseedVtxMinAnyPV( 0.1 ) - , m_DtrakMax( 0.5 ) - , m_PtSeedsMin( 1000. ) - , m_PtMergedSeedsMin( 1000. ) - , m_SeedsMaxChi2DoF( 50. ) - , m_Triplets( false ) - , m_DRmin( 0.1 ) - , m_DRmax( 10. ) - , m_TrkChi2DoF( 2.5 ) - , m_DeltaRSeeds( 0.5 ) - , m_preFilter( true ) - , m_jetFilter( true ) - , m_PVveto( true ) { - // - declareInterface<IJetMaker>( this ); - // - declareProperty( "SeedID", m_seedID, "Particle ID for the Seed" ); - declareProperty( "SeedRParameter", m_r ); - declareProperty( "Sort", m_sort, - "Sorting Criteria for jets [0:none,1:pt,2:E,3:eta, 4:ProtoSeed the Dau then Pt, default:4]" ); - // define momentum combiner - declareProperty( "ParticleCombiner", m_combinerName ); - declareProperty( "SeedPtTrackMin", m_PtTrackMin, "pt of the track used for Vertexing" ); - declareProperty( "SeedPTrackMin", m_PTrackMin, "p of the track used for Vertexing" ); - declareProperty( "SeedIPmin", m_IPmin, "ip of the track used for Vertexing" ); - declareProperty( "SeedSignif", m_Signif, "signif oft he track used for Vertexing" ); - declareProperty( "SeedDMK0", m_DMK0, "mass window for Ks" ); - declareProperty( "SeedTseedVtxMin", m_TseedVtxMin, "min distance btw PV and the vtx" ); - declareProperty( "SeedTseedVtxMax", m_TseedVtxMax, "max distance btw PV and the vtx" ); - declareProperty( "SeedTseedVtxMinAnyPV", m_TseedVtxMinAnyPV, "min distance btw any PV and the vtx" ); - - declareProperty( "SeedDtrakMax", m_DtrakMax, "dca window for vtx" ); - declareProperty( "SeedPtSeedsMin", m_PtSeedsMin, "min pt of the seeds" ); - declareProperty( "SeedVtxMaxChi2PerDoF", m_SeedsMaxChi2DoF, "max chi2 per dof for the vtx fit of the seed" ); - declareProperty( "PtMergedSeedsMin", m_PtMergedSeedsMin, "min pt of the merged seeds" ); - - declareProperty( "SeedTriplets", m_Triplets, "built Vtx with 3 tracks" ); - declareProperty( "SeedDRmin", m_DRmin, "min positon in R of the vtx" ); - declareProperty( "SeedDRmax", m_DRmax, "max positon in R of the vtx" ); - declareProperty( "SeedTrkChi2PerDoF", m_TrkChi2DoF, "max chi2PerDoF for the track used for the vtx" ); - declareProperty( "SeedsDeltaR", m_DeltaRSeeds, "min DR btw 2 seeds otherwise:merge" ); - declareProperty( "SeedpreFilter", m_preFilter, "pre filter the protoseeds" ); - declareProperty( "SeedFilter", m_jetFilter, "filter the protoseeds" ); - declareProperty( "vetoPV", m_PVveto, - "exclude vertex near to any PV with distance fixe by 'SeedTseedVtxMinAnyPV' and exclude the trk " - "associated to any PV to construct a secondary vtx" ); - } - /** standard initialization of the tool - * @return status code - */ - StatusCode initialize() override; - - protected: - /// make the detailed check of all parameters - - inline StatusCode check() const { - if ( 0 > m_ptmin ) { Warning( "PtMin is negative " ).ignore(); } - return StatusCode::SUCCESS; - } - - Gaudi::Property<std::string> m_standardGeometry_address{this, "StandardGeometryTop", "/dd/Structure/LHCb"}; - - int to_user_index( const int index ) const { return index + 10000; } - int from_user_index( const int index ) const { return index - 10000; } - // proposed jet ID - int m_seedID; ///< proposed jet ID - // R-parameter - double m_r; ///< R-parameters - // ptMin-parameter - double m_ptmin; ///< pt-min parameter - // jet sorting criteria - int m_sort; ///< jet sorting criteria - // combiner - std::string m_combinerName; - mutable IParticleCombiner* m_combiner; ///< combiner to be used - - IDistanceCalculator* m_dist; - IVertexFit* m_fitter; - - double m_Rmax; - double m_PtTrackMin; - double m_PTrackMin; - double m_IPmin; - double m_Signif; - double m_DMK0; - double m_TseedVtxMin; - double m_TseedVtxMax; - double m_TseedVtxMinAnyPV; - double m_DtrakMax; - double m_PtSeedsMin; - double m_PtMergedSeedsMin; - double m_SeedsMaxChi2DoF; - bool m_Triplets; - double m_DRmin; - double m_DRmax; - double m_TrkChi2DoF; - double m_DeltaRSeeds; - bool m_preFilter; - bool m_jetFilter; - bool m_PVveto; - - IJetMaker::Jets JetCone( double const& Rmax, IJetMaker::Jets& seeds, IJetMaker::Input const& inputs, - IGeometryInfo const& geometry ) const; - - IJetMaker::Jets FilterProtoJets( const double& DeltaRSeeds, IJetMaker::Jets& ProtoJets ) const; - - double getDeltaR( LHCb::Particle* p1, LHCb::Particle* p2 ) const; - - void RemoveTracks( LHCb::Particle::ConstVector& particles, const LHCb::RecVertex PV ) const; - - private: - mutable Gaudi::Accumulators::StatCounter<> m_nSeeds{this, "#seeds"}; - }; - -} // namespace LoKi - -class sortDauPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - if ( obj1->weight() == obj2->weight() ) - if ( obj1->daughtersVector().size() == obj2->daughtersVector().size() ) - return obj1->pt() > obj2->pt(); - else - return obj1->daughtersVector().size() > obj2->daughtersVector().size(); - else - return obj1->weight() > obj2->weight(); - } -}; - -class sortE { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().E() > obj2->momentum().E(); - } -}; - -class sortEta { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().Eta() > obj2->momentum().Eta(); - } -}; - -class sortPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { return obj1->pt() > obj2->pt(); } -}; - -/* standard initialization of the tool - * @return status code - */ -StatusCode LoKi::VVSeedFinder::initialize() { - StatusCode sc = GaudiTool::initialize(); - if ( sc.isFailure() ) { return sc; } - - svc<LoKi::ILoKiSvc>( "LoKiSvc", true ); - - m_dist = tool<IDistanceCalculator>( "LoKi::DistanceCalculator" ); - if ( !m_dist ) { - err() << "Unable to Retrieve LoKi::DistanceCalculator" << endmsg; - return StatusCode::FAILURE; - } - - m_combiner = tool<IParticleCombiner>( m_combinerName, this ); - if ( !m_combiner ) { - err() << "Unable to Retrieve Default ParticleCombiner" << endmsg; - return StatusCode::FAILURE; - } - - m_fitter = tool<IVertexFit>( "OfflineVertexFitter" ); - if ( !m_fitter ) { - err() << "Unable to Retrieve Default VertexFitter" << endmsg; - return StatusCode::FAILURE; - } - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << " ------------------------------------------------------ " << endmsg; - debug() << " | | " << endmsg; - debug() << " | \\ /\\ / | " << endmsg; - debug() << " | \\/ \\/SeedsFinder Tool for Phys Particles | " << endmsg; - debug() << " | | " << endmsg; - debug() << " | \\ / | " << endmsg; - debug() << " | \\/ EPFLausanne - 2009-11 | " << endmsg; - debug() << " | | A. Bay & C. Potterat | " << endmsg; - debug() << " | | | " << endmsg; - debug() << " | _ _ | " << endmsg; - debug() << " | | |_||_||_ | " << endmsg; - debug() << " | |__ | | ||__ | " << endmsg; - debug() << " ------------------------------------------------------ " << endmsg; - debug() << " | | " << endmsg; - debug() << " | " << endmsg; - debug() << " | VVSeedsFinder parameters: " << endmsg; - debug() << " | SeedID " << m_seedID << endmsg; - debug() << " | SeedTriplets " << m_Triplets << endmsg; - debug() << " | SeedRParameter " << m_r << endmsg; - debug() << " | PtTrackMin " << m_PtTrackMin << endmsg; - debug() << " | PTrackMin " << m_PTrackMin << endmsg; - debug() << " | SeedTrkChi2PerDoF " << m_TrkChi2DoF << endmsg; - debug() << " | IPmin " << m_IPmin << endmsg; - debug() << " | Signif " << m_Signif << endmsg; - debug() << " | DMK0 " << m_DMK0 << endmsg; - debug() << " | TseedVtxMin " << m_TseedVtxMin << endmsg; - debug() << " | TseedVtxMax " << m_TseedVtxMax << endmsg; - debug() << " | DtrakMax " << m_DtrakMax << endmsg; - debug() << " | DeltaRSeeds " << m_DeltaRSeeds << endmsg; - debug() << " | PtSeedsMin " << m_PtSeedsMin << endmsg; - debug() << " | PtMergedSeedsMin " << m_PtMergedSeedsMin << endmsg; - debug() << " | " << endmsg; - debug() << " | | " << endmsg; - debug() << " ------------------------------------------------------ " << endmsg; - } - - return StatusCode::SUCCESS; -} -// =========================================================================== -// find seeds -// =========================================================================== -StatusCode LoKi::VVSeedFinder::makeJets( const IJetMaker::Input& input_, IJetMaker::Jets& jets_ ) const { - - const LHCb::RecVertex::Container* verts = get<LHCb::RecVertex::Container>( LHCb::RecVertexLocation::Primary ); - if ( verts->size() != 1 ) return StatusCode::SUCCESS; - LHCb::RecVertex::Container::const_iterator iv = verts->begin(); - const LHCb::RecVertex vtx = **iv; - makeJets( input_, vtx, jets_ ).ignore(); - return StatusCode::SUCCESS; -} - -StatusCode LoKi::VVSeedFinder::makeJets( const IJetMaker::Input& input_, const LHCb::RecVertex& RecVert_, - IJetMaker::Jets& jets_ ) const { - - const LHCb::RecVertex::Container* verts = get<LHCb::RecVertex::Container>( LHCb::RecVertexLocation::Primary ); - LHCb::RecVertex::Container::const_iterator iv; - LHCb::VertexBase* RecVert = RecVert_.clone(); - - std::vector<LHCb::VertexBase*> PVs; - - for ( iv = verts->begin(); iv != verts->end(); iv++ ) - if ( ( *iv )->position() != RecVert->position() ) PVs.push_back( ( *iv )->clone() ); - - IJetMaker::Jets Seeds; - double ipl, chi2l; - StatusCode sc; - LHCb::Vertex vtx, vseedbest; - double dca12, dcae12; - // Kshort Mass - const double MK0 = 497.7; - - std::vector<Gaudi::XYZVector> Slopes_flag; - - LHCb::Particle::ConstVector PartIPK0Sub; - LHCb::Particle::ConstVector PartIP; - LHCb::Particle::ConstVector all; - LHCb::Particle::ConstVector KsDau; - - double tof = 0, tofe = 0; - - int testLiHi = 0; - int testClone = 0; - int testLong = 0; - int testCharged = 0; - - // Get the default geometry FIXME, use functional framework - auto lhcb = getDet<IDetectorElement>( m_standardGeometry_address ); - if ( !lhcb ) { throw GaudiException( "Could not load geometry", name(), StatusCode::FAILURE ); } - auto& geometry = *lhcb->geometry(); - - //----------------------------------------// - // test the particles // - //----------------------------------------// - - for ( LHCb::Particle::ConstVector::const_iterator ip = input_.begin(); input_.end() != ip; ++ip ) { - const LHCb::Particle* p = *ip; - if ( 0 == p ) { - Warning( "Invalid input particle" ).ignore(); - continue; - } - const LHCb::ProtoParticle* myProto = p->proto(); - if ( myProto != NULL ) { - const LHCb::Track* trk = myProto->track(); - if ( trk != NULL ) { - if ( trk->flags() == 1 || trk->flags() == 4 ) continue; - // == 1 : backward, == 4 clone - if ( trk->flags() == 1 ) continue; - double CloneDist = trk->info( LHCb::Track::AdditionalInfo::CloneDist, 9999. ); - double TrkLi = trk->likelihood(); - double GhostProb = trk->ghostProbability(); - double LongTrk = static_cast<int>( trk->type() ); - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << "addinfo Likelihood = " << TrkLi << endmsg; - verbose() << "addinfo CloneDist = " << CloneDist << endmsg; - verbose() << "addinfo GhostProb = " << GhostProb << endmsg; - verbose() << "addinfo Longb = " << LongTrk << endmsg; - verbose() << "p Trk = " << trk->momentum() / Gaudi::Units::GeV << endmsg; - } - if ( CloneDist < 9999. ) testClone++; - if ( TrkLi < -60. ) testLiHi++; - if ( LongTrk == 3 ) testLong++; - //----------------------------------------// - // discard clones and Ghost // - //----------------------------------------// - if ( CloneDist < 9999. ) continue; - if ( TrkLi < -60. ) continue; - } - } - - bool m_FilterPart = true; - //----------------------------------------// - // check if 1 particle is not twice on // - // the TES with the slopes // - //----------------------------------------// - bool flag_same = false; - if ( m_FilterPart ) { - if ( Slopes_flag.size() != 0 ) { - for ( int k = 0; k < (int)Slopes_flag.size(); k++ ) { - if ( Slopes_flag.at( k ) == p->slopes() ) { - flag_same = true; - if ( msgLevel( MSG::VERBOSE ) ) - verbose() << "same slopes !! s1: " << Slopes_flag.at( k ) << " s2: " << p->slopes() << endmsg; - break; - } - } - } - if ( !flag_same ) { - Slopes_flag.push_back( p->slopes() ); - } else { - continue; - } - } - - //----------------------------------------// - // save all the part that arrived here // - // as the global container // - //----------------------------------------// - all.push_back( p ); - - if ( p->charge() == 0 ) continue; - testCharged++; - if ( p->proto() == NULL ) continue; - if ( p->proto()->track() == NULL ) continue; - if ( p->proto()->track()->type() != LHCb::Track::Types::Long ) continue; - if ( p->proto()->track()->chi2PerDoF() > m_TrkChi2DoF ) continue; - - sc = m_dist->distance( p, RecVert, ipl, chi2l, geometry ); - - if ( p->p() < m_PTrackMin ) continue; - if ( p->pt() < m_PtTrackMin ) continue; - if ( ipl < m_IPmin ) continue; - if ( sqrt( chi2l ) < m_Signif ) continue; - - //----------------------------------------// - // save all the charge, with a good Pt // - // and a good IP/IPE as the input // - // particles for the Cone Jet Algo // - //----------------------------------------// - - PartIP.push_back( p ); - } - - //----------------------------------------// - // discard the trk coming any VTX // - //----------------------------------------// - - if ( m_PVveto ) - for ( iv = verts->begin(); iv != verts->end(); iv++ ) RemoveTracks( PartIP, **iv ); - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "nb of ghost (LiH) : " << testLiHi << endmsg; - debug() << "nb of clone : " << testClone << endmsg; - debug() << "nb of long : " << testLong << endmsg; - debug() << "nb of charged : " << testCharged << endmsg; - debug() << "nb of 'photon' : " << (int)all.size() - double( testCharged ) << endmsg; - debug() << "Particle ALL size : " << input_.size() << endmsg; - debug() << "Particle ALL size (JetCone) : " << all.size() << endmsg; - debug() << "Particle INPUT for JET : " << PartIP.size() << endmsg; - } - - int cntK0 = 0; - - //============================================================================= - // COMBINEPIPIMINK0: GET RIF OF K0S - //============================================================================= - for ( LHCb::Particle::ConstVector::const_iterator ip = PartIP.begin(); PartIP.end() != ip; ++ip ) { - const LHCb::Particle* p = *ip; - bool testK = false; - - for ( LHCb::Particle::ConstVector::const_iterator kp = KsDau.begin(); KsDau.end() != kp; ++kp ) - if ( *ip == *kp ) { - testK = true; - break; - } - - if ( !testK ) - for ( LHCb::Particle::ConstVector::const_iterator ip2 = ip + 1; PartIP.end() != ip2; ++ip2 ) { - const LHCb::Particle* p2 = *ip2; - - sc = m_dist->distance( ( p ), ( p2 ), dca12, dcae12, geometry ); - if ( !sc ) { warning() << "can't mesure the dist " << endmsg; } - if ( dca12 > m_DtrakMax ) continue; // dca too large bwt the tracks - - Gaudi::LorentzVector sum = ( p )->momentum() + ( p2 )->momentum(); - - //============================================================================= - // COMBINEPIPIMINK0: GET RIF OF K0S - //============================================================================= - if ( ( p )->particleID().abspid() == 211 && ( p2 )->particleID().abspid() == 211 && - ( ( p )->charge() ) * ( ( p2 )->charge() ) < 0 ) - if ( fabs( sum.M() - MK0 ) < m_DMK0 ) { - cntK0++; - testK = true; - - KsDau.push_back( p ); - KsDau.push_back( p2 ); - } - } - - if ( testK ) continue; - PartIPK0Sub.push_back( p ); - } - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "nb K0 : " << cntK0 << endmsg; - debug() << "Particle INPUT for JET no K0 : " << PartIPK0Sub.size() << endmsg; - } - - if ( PartIPK0Sub.size() < 2 ) { - Warning( "Not enough good part for seeding" ).ignore(); - return StatusCode::SUCCESS; - } - - Gaudi::XYZPoint BL_P = Gaudi::XYZPoint( 0, 0, 0 ); - Gaudi::LorentzVector BL_M = Gaudi::LorentzVector( 0, 0, 1, 0 ); - - if ( exist<LHCb::Particle::Range>( "/Event/BeamLine" ) ) { - const LHCb::Particle::Range BL = get<LHCb::Particle::Range>( "/Event/BeamLine" ); - const LHCb::Particle* tmp = *( BL.begin() ); - BL_P = Gaudi::XYZPoint( tmp->referencePoint() ); - BL_M = Gaudi::LorentzVector( tmp->momentum() ); - // m_BeamLine->setMomentum( tmp->momentum() ); - if ( msgLevel( MSG::DEBUG ) ) debug() << "Beam line position " << BL_P << " direction " << BL_M << endmsg; - } else { - if ( msgLevel( MSG::DEBUG ) ) - debug() << "No Beam line found at " - << "/Event/BeamLine" << endmsg; - - BL_P = Gaudi::XYZPoint( RecVert->position() ); - - if ( msgLevel( MSG::DEBUG ) ) debug() << "Beam line position " << BL_P << " direction " << BL_M << endmsg; - } - - // COMBIN PART - - LHCb::Particle::ConstVector::const_iterator jp, kp, lp; - for ( jp = PartIPK0Sub.begin(); jp != PartIPK0Sub.end(); jp++ ) { - for ( kp = jp + 1; kp != PartIPK0Sub.end(); kp++ ) { - - sc = m_dist->distance( ( *jp ), ( *kp ), dca12, dcae12, geometry ); - if ( !sc ) { warning() << "can't mesure the dist " << endmsg; } - if ( dca12 > m_DtrakMax ) continue; // dca too large bwt the tracks - - Gaudi::LorentzVector sum = ( *jp )->momentum() + ( *kp )->momentum(); - - if ( m_Triplets ) { - - // add Try 3-tracks Seed otpion - for ( lp = ( kp + 1 ); lp != PartIPK0Sub.end(); lp++ ) { - - double dca13, dcae13; - m_dist->distance( ( *jp ), ( *lp ), dca13, dcae13, geometry ).ignore(); - if ( dca13 > m_DtrakMax ) continue; // dca too large btw the trks - - double dca23, dcae23; - m_dist->distance( ( *lp ), ( *kp ), dca23, dcae23, geometry ).ignore(); - if ( dca23 > m_DtrakMax ) continue; // dca too large btw the trks - - StatusCode scvtx = m_fitter->fit( vtx, **jp, **kp, **lp, geometry ); - if ( scvtx.isFailure() ) { - warning() << "VTX Fit failed" << endmsg; - continue; - } - - if ( vtx.chi2PerDoF() > m_SeedsMaxChi2DoF ) continue; - // cut on chi2 - double mydz = -1; - if ( scvtx ) mydz = vtx.position().z() - RecVert->position().z(); - if ( mydz < 0 ) continue; - // the vtx is before the PV - - // intersection of the beam line with the XY plane, - // find the lambda parameter of the line. - double lambda = ( vtx.position().z() - BL_P.z() ) / BL_M.z(); - - // find x and y of intersection point - double x = BL_P.x() + lambda * BL_M.x(); - double y = BL_P.y() + lambda * BL_M.y(); - - x -= vtx.position().x(); - y -= vtx.position().y(); - - double dr = sqrt( x * x + y * y ); - if ( dr < m_DRmin || dr > m_DRmax ) continue; - - // the vtx is too close or to far from the beam direction - - //--------------------------- - // Seed level cuts ----------- - //--------------------------- - LHCb::VertexBase* PPvtx2 = vtx.clone(); - m_dist->distance( PPvtx2, RecVert, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMin || tof > m_TseedVtxMax ) continue; - // the vtx is too close or to far from the PV (~time of flight) - bool PVveto = false; - if ( m_PVveto ) - for ( std::vector<LHCb::VertexBase*>::iterator pv = PVs.begin(); pv != PVs.end(); pv++ ) { - m_dist->distance( PPvtx2, *pv, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMinAnyPV ) { - PVveto = true; - break; - } - } - if ( PVveto ) continue; - - Gaudi::LorentzVector sum2 = sum + ( *lp )->momentum(); - if ( sum2.Pt() < m_PtSeedsMin ) continue; - // Pt too soft - - LHCb::Particle::ConstVector daughters2; - LHCb::Particle pSeed2; - LHCb::Vertex vSeed2; - - daughters2.push_back( *lp ); - daughters2.push_back( *kp ); - daughters2.push_back( *jp ); - pSeed2.setParticleID( LHCb::ParticleID( m_seedID ) ); - pSeed2.setReferencePoint( RecVert->position() ); - pSeed2.setMomentum( sum2 ); - - StatusCode sc = m_combiner->combine( daughters2, pSeed2, vSeed2, geometry ); - - pSeed2.setEndVertex( vtx.clone() ); - // save the trks as the daugthers - // remove - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - Seeds.push_back( pSeed2.clone() ); - // remove lp to test another part - - } // end loop lp - } - - // 2-seed - LHCb::Particle pSeed; - LHCb::Particle::ConstVector daughters; - LHCb::Vertex vSeed; - - StatusCode scvtx = m_fitter->fit( vtx, **jp, **kp, geometry ); - - if ( scvtx.isFailure() ) { - warning() << "VTX Fit failed" << endmsg; - continue; - } - - if ( vtx.chi2PerDoF() > m_SeedsMaxChi2DoF ) continue; - - double mydz = -1; - if ( scvtx ) mydz = vtx.position().z() - RecVert->position().z(); - if ( mydz < 0 ) continue; - // the vtx is before the PV - - // Gaudi::XYZPoint dv = Gaudi::XYZPoint(vtx.position() - m_BeamLine->referencePoint().position()); - - // intersection of the beam line with the XY plane, - // find the lambda parameter of the line. - double lambda = ( vtx.position().z() - BL_P.z() ) / BL_M.z(); - - // find x and y of intersection point - double x = BL_P.x() + lambda * BL_M.x(); - double y = BL_P.y() + lambda * BL_M.y(); - - x -= vtx.position().x(); - y -= vtx.position().y(); - - double dr = sqrt( x * x + y * y ); - if ( dr < m_DRmin || dr > m_DRmax ) continue; - // the vtx is too close or to far from the beam direction - - //--------------------------- - // Seed level cuts ----------- - //--------------------------- - LHCb::VertexBase* PPvtx2 = vtx.clone(); - m_dist->distance( PPvtx2, RecVert, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMin || tof > m_TseedVtxMax ) continue; - // the vtx is too close or to far from the PV (~time of flight) - - bool PVveto = false; - if ( m_PVveto ) - for ( std::vector<LHCb::VertexBase*>::iterator pv = PVs.begin(); pv != PVs.end(); pv++ ) { - m_dist->distance( PPvtx2, *pv, tof, tofe ).ignore(); - if ( tof < m_TseedVtxMinAnyPV ) { - PVveto = true; - break; - } - } - if ( PVveto ) continue; - - if ( sum.Pt() < m_PtSeedsMin ) continue; - // Pt too soft - - daughters.push_back( *jp ); - daughters.push_back( *kp ); - - pSeed.setParticleID( LHCb::ParticleID( m_seedID ) ); - pSeed.setReferencePoint( RecVert->position() ); - pSeed.setMomentum( sum ); - StatusCode sc = m_combiner->combine( daughters, pSeed, vSeed, geometry ); - // save the trks as the daugthers - pSeed.setEndVertex( vtx.clone() ); - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - Seeds.push_back( pSeed.clone() ); - - } // end loop kp - } // end loop jp - - if ( Seeds.empty() ) { - Warning( "No Seeds from LoKiVVSeedFinder" ).ignore(); - return StatusCode::SUCCESS; - } - - //============================================================================= - // SORT THE SEEDS - //============================================================================= - if ( m_sort == 1 ) std::sort( Seeds.begin(), Seeds.end(), sortPt() ); - if ( m_sort == 2 ) std::sort( Seeds.begin(), Seeds.end(), sortE() ); - if ( m_sort == 3 ) std::sort( Seeds.begin(), Seeds.end(), sortEta() ); - if ( m_sort == 4 ) std::sort( Seeds.begin(), Seeds.end(), sortDauPt() ); - - IJetMaker::Jets ProtoJets, MergeProtoJets; - ProtoJets.reserve( Seeds.size() ); - MergeProtoJets.reserve( Seeds.size() ); - - int myj = 0; - for ( IJetMaker::Jets::iterator ijet = Seeds.begin(); Seeds.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "Seeds " << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - } - - //============================================================================= - // BUILDPROTOJETS: BUILDS JETS AROUND THE SEEDS FORMED WITH THE PRIMARY VERTEX AND THE VERTEXCOMBINATIONS - //============================================================================= - ProtoJets = JetCone( m_r, Seeds, input_, geometry ); - - myj = 0; - for ( IJetMaker::Jets::iterator ijet = ProtoJets.begin(); ProtoJets.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "BSeeds " << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - } - - //============================================================================= - // FILTERPROTOJETS: MERGES JETS WHICH ARE TOO CLOSE OR THE SAME - //============================================================================= - MergeProtoJets = FilterProtoJets( m_DeltaRSeeds, ProtoJets ); - - if ( msgLevel( MSG::INFO ) ) { m_nSeeds += Seeds.size(); } - - jets_ = MergeProtoJets; - - myj = 0; - for ( IJetMaker::Jets::iterator ijet = MergeProtoJets.begin(); MergeProtoJets.end() != ijet; ++ijet ) { - LHCb::Particle* pJet = *ijet; - - myj++; - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "CSeeds " << myj << "| nb of combin Seeds: " << pJet->weight() - << " nb part in the jet: " << pJet->daughtersVector().size() << " Pt: " << pJet->pt() << endmsg; - debug() << " | m : " << pJet->momentum() << endmsg; - debug() << " | eta: " << pJet->momentum().eta() << " / phi: " << pJet->momentum().phi() << endmsg; - } - } - - return StatusCode::SUCCESS; -} - -//============================================================================= -// JETCONE: THE HEART OF EVERYTHING, BUILDS A JET AROUND A SEED -//============================================================================= -IJetMaker::Jets LoKi::VVSeedFinder::JetCone( double const& Rmax, IJetMaker::Jets& seeds, IJetMaker::Input const& inputs, - IGeometryInfo const& geometry ) const { - - LHCb::Particle::Vector pJetVec; - - IJetMaker::Jets::const_iterator iseeds; - - for ( iseeds = seeds.begin(); iseeds != seeds.end(); ++iseeds ) { - - const LHCb::Particle* pSeed = *iseeds; - - LHCb::Particle pJet; - LHCb::Vertex vJet; - LHCb::Particle::ConstVector daughters; - - pJet.setParticleID( LHCb::ParticleID( 9999999 ) ); - - pJet.setReferencePoint( pSeed->referencePoint() ); - - const double PI = 3.14159265; - - double phiSeed = pSeed->momentum().phi(); - double etaSeed = pSeed->momentum().eta(); - Gaudi::LorentzVector ptot = Gaudi::LorentzVector( 0., 0., 0., 0. ); - - LHCb::Particle::ConstVector::const_iterator ip; - - for ( ip = inputs.begin(); ip != inputs.end(); ++ip ) { - - const LHCb::Particle* myPart = *ip; - Gaudi::LorentzVector p1 = myPart->momentum(); - - double phiPart = myPart->momentum().phi(); - double etaPart = myPart->momentum().eta(); - double Dphi = fabs( phiSeed - phiPart ); - if ( Dphi > PI ) { Dphi = 2 * PI - Dphi; }; - double R = sqrt( Dphi * Dphi + ( etaSeed - etaPart ) * ( etaSeed - etaPart ) ); - - if ( R > Rmax ) continue; - ptot = ptot + p1; - - daughters.push_back( myPart ); - } - - if ( daughters.empty() ) Warning( "Empty list of of daughter particles, skip it" ).ignore(); - - StatusCode sc = m_combiner->combine( daughters, pJet, vJet, geometry ); - if ( sc.isFailure() ) Warning( "Error from momentum combiner, skip", sc, 0 ).ignore(); - - pJet.setMomentum( Gaudi::LorentzVector( ptot ) ); - pJetVec.push_back( pJet.clone() ); - } - - return pJetVec; -} - -//============================================================================= -// JETCONE: END -//============================================================================= - -//============================================================================= -// FILTERPROTOJETS: MERGES JETS WHICH ARE TOO CLOSE OR THE SAME -//============================================================================= -IJetMaker::Jets LoKi::VVSeedFinder::FilterProtoJets( const double& DeltaRSeeds, IJetMaker::Jets& ProtoJets ) const { - //---------------------------------------------------------------------- - // filter protojets which are too close or the same-------- ------------ - //---------------------------------------------------------------------- - - // const double PI=3.14159265; - // erase the 6th element - // myvector.erase (myvector.begin()+5); - - if ( msgLevel( MSG::DEBUG ) ) debug() << "size of noFiltered ProtoJets = " << ProtoJets.size() << endmsg; - - IJetMaker::Jets NewProtoJets; - IJetMaker::Jets NewProtoJets2; - - bool testJets = true; - int mySize = 0; - int k = 0; - if ( m_preFilter ) { - for ( IJetMaker::Jets::iterator ip = ProtoJets.begin(); ip != ProtoJets.end(); ++ip ) { - LHCb::Particle* seed = ( *ip ); - if ( seed->momentum().E() > 0. ) { - int j = 0; - for ( IJetMaker::Jets::iterator ip2 = ip + 1; ip2 != ProtoJets.end(); ++ip2 ) { - j++; - - LHCb::Particle* seed2 = ( *ip2 ); - if ( seed2->momentum().E() > 0. ) { - - double testR = getDeltaR( seed, seed2 ); - - if ( testR < 0.00005 ) { - // if( testR < 0.05 || (fabs(v1.E() - v2.E())/v1.E() < 0.05 && testR < 0.1)){ - // SAME JET ! - ProtoJets.at( k )->setMomentum( Gaudi::LorentzVector( seed->momentum() ) ); - - double myWeight = seed->weight() + seed2->weight(); - - ProtoJets.at( k )->eraseInfo( LHCb::Particle::additionalInfo::Weight ); - ProtoJets.at( k )->addInfo( LHCb::Particle::additionalInfo::Weight, myWeight ); - - seed->eraseInfo( LHCb::Particle::additionalInfo::Weight ); - seed->addInfo( LHCb::Particle::additionalInfo::Weight, myWeight ); - - seed->setMomentum( Gaudi::LorentzVector( seed->momentum() ) ); - ProtoJets.at( k + j )->setMomentum( Gaudi::LorentzVector( 0., 0., 0., 0. ) ); - } - } - } - } - k++; - } - } - for ( IJetMaker::Jets::iterator ip = ProtoJets.begin(); ip != ProtoJets.end(); ++ip ) { - - LHCb::Particle* seed = ( *ip ); - if ( seed->momentum().E() > 0. ) { NewProtoJets.push_back( seed ); } - } - - ProtoJets.clear(); - ProtoJets = NewProtoJets; - - if ( msgLevel( MSG::DEBUG ) ) debug() << "size of preFiltered ProtoJets = " << ProtoJets.size() << endmsg; - - if ( m_jetFilter ) { - while ( testJets ) { - - NewProtoJets.clear(); - mySize = (int)ProtoJets.size(); - testJets = false; - int k = 0; - int J1test = 0; - int J2test = 0; - double testDR = 100.0; - - NewProtoJets = ProtoJets; - - for ( IJetMaker::Jets::iterator ip = ProtoJets.begin(); ip != ProtoJets.end(); ++ip ) { - LHCb::Particle* seed = ( *ip ); - - if ( seed->momentum().E() > 0. ) { - int j = 0; - for ( IJetMaker::Jets::iterator ip2 = ip + 1; ip2 != ProtoJets.end(); ++ip2 ) { - j++; - LHCb::Particle* seed2 = ( *ip2 ); - - if ( seed2->momentum().E() > 0. ) { - - double R = getDeltaR( seed, seed2 ); - - if ( R < DeltaRSeeds && R < testDR ) { - - testDR = R; - J1test = k; - J2test = j + k; - } - } - } - } - k++; - } - - if ( testDR < DeltaRSeeds ) { - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << " FilterProtoJet merge " << endmsg; - verbose() << "Nseed phi/eta= " << ( ProtoJets.at( J1test )->momentum() ).phi() << " " - << ( ProtoJets.at( J1test )->momentum() ).eta() << endmsg; - verbose() << "Nseed2phi/eta= " << ( ProtoJets.at( J2test )->momentum() ).phi() << " " - << ( ProtoJets.at( J2test )->momentum() ).eta() << endmsg; - verbose() << "E P1= " << ( ProtoJets.at( J1test )->momentum() ).E() - << " E P2= " << ( ProtoJets.at( J2test )->momentum() ).E() << endmsg; - verbose() << "E P1/P2= " - << fabs( ( ProtoJets.at( J1test )->momentum() ).E() - ( ProtoJets.at( J2test )->momentum() ).E() ) / - ( ProtoJets.at( J1test )->momentum() ).E() - << endmsg; - } - if ( testDR > LHCb::Math::looseTolerance || !m_preFilter ) { - - NewProtoJets.at( J1test )->setMomentum( - Gaudi::LorentzVector( ( ProtoJets.at( J1test )->momentum() + ProtoJets.at( J2test )->momentum() ) ) ); - - const LHCb::Particle::ConstVector DaughtersJ2 = ProtoJets.at( J2test )->daughtersVector(); - for ( LHCb::Particle::ConstVector::const_iterator ip3 = DaughtersJ2.begin(); ip3 != DaughtersJ2.end(); - ++ip3 ) { - const LHCb::Particle* d1 = *ip3; - NewProtoJets.at( J1test )->addToDaughters( d1 ); - } - - double myWeight = ProtoJets.at( J1test )->weight() + ProtoJets.at( J2test )->weight(); - double myCL = -1; - - if ( ProtoJets.at( J1test )->momentum().E() > ProtoJets.at( J2test )->momentum().E() ) { - NewProtoJets.at( J1test )->setReferencePoint( ProtoJets.at( J1test )->referencePoint() ); - myCL = ProtoJets.at( J1test )->confLevel(); - } else { - NewProtoJets.at( J1test )->setReferencePoint( ProtoJets.at( J2test )->referencePoint() ); - myCL = ProtoJets.at( J2test )->confLevel(); - } - NewProtoJets.at( J1test )->eraseInfo( LHCb::Particle::additionalInfo::Weight ); - NewProtoJets.at( J1test )->addInfo( LHCb::Particle::additionalInfo::Weight, myWeight ); - NewProtoJets.at( J1test )->eraseInfo( LHCb::Particle::additionalInfo::ConfLevel ); - NewProtoJets.at( J1test )->addInfo( LHCb::Particle::additionalInfo::ConfLevel, myCL ); - - } else { - - // Same jet... nothing to do.... just erase one. - } - - NewProtoJets.erase( NewProtoJets.begin() + J2test ); - } - - if ( msgLevel( MSG::VERBOSE ) ) { - verbose() << " ProtoJet Size : " << mySize << endmsg; - verbose() << " ProtoJet Size : " << ProtoJets.size() << endmsg; - verbose() << " NewFilterProtoJet Size: " << NewProtoJets.size() << endmsg; - } - - if ( (int)NewProtoJets.size() != mySize ) testJets = true; - - ProtoJets.clear(); - ProtoJets = NewProtoJets; - - if ( NewProtoJets.size() % 10 == 0 ) - if ( msgLevel( MSG::VERBOSE ) ) - verbose() << "Number of filtered ProtoJets is " << NewProtoJets.size() << endmsg; - } - - if ( msgLevel( MSG::DEBUG ) ) - debug() << "size of Filtered ProtoJets = " << NewProtoJets.size() << endmsg; - - NewProtoJets.clear(); - for ( IJetMaker::Jets::iterator ip2 = ProtoJets.begin(); ip2 != ProtoJets.end(); ++ip2 ) { - if ( ( *ip2 )->momentum().Pt() < m_PtMergedSeedsMin ) continue; - NewProtoJets.push_back( ( *ip2 ) ); - } - } - - return NewProtoJets; -} -//============================================================================= -// END OF FILTERPROTOJETS -//============================================================================= - -//============================================================================= -// GETDELTAR: CALCULATES DELTAR -//============================================================================= -double LoKi::VVSeedFinder::getDeltaR( LHCb::Particle* p1, LHCb::Particle* p2 ) const { - - const double PI = 3.14159265; - double R, Dphi, phi1, phi2, e1, e2; - phi1 = p1->momentum().phi(); - phi2 = p2->momentum().phi(); - e1 = p1->momentum().eta(); - e2 = p2->momentum().eta(); - Dphi = fabs( phi1 - phi2 ); - if ( Dphi > PI ) { Dphi = 2 * PI - Dphi; }; - R = sqrt( Dphi * Dphi + ( e1 - e2 ) * ( e1 - e2 ) ); - return R; -} - -//============================================================================= -// Remove trks form a PV -//============================================================================= - -void LoKi::VVSeedFinder::RemoveTracks( LHCb::Particle::ConstVector& particles, const LHCb::RecVertex PV ) const { - - // Remove all tracks from the PV - LHCb::Particle::ConstVector tmp; - SmartRefVector<LHCb::Track>::const_iterator iPV = PV.tracks().begin(); - int endkey = PV.tracks().back()->key(); - - for ( LHCb::Particle::ConstVector::const_iterator i = particles.begin(); i != particles.end(); ++i ) { - if ( ( *i )->proto()->track() == NULL ) continue; - const LHCb::Track* tk = ( *i )->proto()->track(); - - while ( ( ( *iPV )->key() < tk->key() ) && ( *iPV )->key() != endkey ) { iPV++; } - if ( ( *iPV )->key() == tk->key() ) { - if ( ( *iPV )->key() != endkey ) iPV++; - continue; - } - tmp.push_back( *i ); - } - // Copy back tmp to particles - particles = tmp; -} - -/// The factory -DECLARE_COMPONENT( LoKi::VVSeedFinder ) diff --git a/Phys/LoKiJets/src/LoKiVVSeedFinder.h b/Phys/LoKiJets/src/LoKiVVSeedFinder.h deleted file mode 100644 index e688d9eacc5729c8827bbde30663586905d26fd7..0000000000000000000000000000000000000000 --- a/Phys/LoKiJets/src/LoKiVVSeedFinder.h +++ /dev/null @@ -1,274 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -// ============================================================================ -#pragma once - -#include "GaudiAlg/GaudiTool.h" - -#include "Kernel/IJetMaker.h" -#include "Kernel/IParticleCombiner.h" - -#include "Event/Particle.h" -#include "Kernel/IDistanceCalculator.h" -#include "Kernel/IVertexFit.h" -#include "LoKi/Geometry.h" -#include "LoKi/Kinematics.h" - -#include "Event/RecVertex.h" -#include "GaudiAlg/ITupleTool.h" -#include "Kernel/IParticleTransporter.h" -#include "LHCbMath/LHCbMath.h" - -namespace LoKi { - - /** @class VVSeedFinder - * - * The VVSeedFinder (EPFL) implementaion of interface IJetMaker - * @see IJetMaker - * - * This file is a part of LoKi project - - * "C++ ToolKit for Smart and Friendly Physics Analysis" - * - * The package has been designed with the kind help from - * Galina PAKHLOVA and Sergey BARSUK. Many bright ideas, - * contributions and advices from G.Raven, J.van Tilburg, - * A.Golutvin, P.Koppenburg have been used in the design. - * - * @author Cedric POTTERAT cedric.potterat@cern.ch - * @date 2011-01-31 - */ - class VVSeedFinder : public virtual IJetMaker, public GaudiTool { - public: - /** The main method: seed-finding procedure - * - * @code - * - * // get the tool - * const IJetMaker* seedMaker = tool<IJetMaker> ( .... ) ; - * - * // input particles - * IJetMaker::Inputs input = ... - * // 1) - * // const Particles* particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles->begin() , particles->end() ) ; - * // 2) - * // LHCb::Particle::ConstVector particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * // 3) - * // LoKi::Range particles = .... ; - * // // create the input container - * // IJetMaker::Inputs input( particles.begin() , particles.end() ) ; - * - * // placeholder for "output" jets - * IJetMaker::Jets seeds ; - * - * // find the jets! - * StatusCode sc = seedMaker -> makeJets ( input , seeds ) ; - * - * // make a loop over jets: - * for ( IJetMaker::Jets::const_iterator iSeed = seeds.begin() ; - * seeds.end() != iSeed ; ++iSeed ) - * { - * // get the jet - * LHCb::Particle* seed = *iSeed ; - * } - * - * @endcode - * - * @attention It is a responsibility of users (e.g. the algorithm) - * to take care about the ownership of jets *AND* their - * vertices). The tool is not intended to do it! - * - * @param input contaainer of input particles - * @param seeds container of output seeds (of type Particle) - * @return status code - */ - StatusCode makeJets( const IJetMaker::Input& input, IJetMaker::Jets& jets ) const override; - StatusCode makeJets( const IJetMaker::Input& input, const LHCb::RecVertex& vtx, - IJetMaker::Jets& jets ) const override; - - /** the standard constructor - * - * @todo The default values for configuration parameters - * (especially for R-parameter) need to be adjusted - * according to EPFL studies. - * @link http://cdsweb.cern.ch/record/1211038/files/LHCb-INT-2009-023.pdf LHCb-INT-2009-023 - Bay/Potterat - * @endlink - * @link http://cdsweb.cern.ch/record/1266883/files/CERN-THESIS-2010-074.pdf CERN_THESIS-2010-074 - Potterat - * @endlink - * - */ - VVSeedFinder( const std::string& type, const std::string& name, const IInterface* parent ) - : GaudiTool( type, name, parent ) - // - , m_seedID( 20098 ) - // - , m_r( 0.15 ) - // - , m_sort( 4 ) - // - , m_combinerName( "MomentumCombiner" ) - , m_combiner( 0 ) - , m_dist( 0 ) - , m_fitter( 0 ) - , m_PtTrackMin( 600.0 ) - , m_PTrackMin( 1000.0 ) - , m_IPmin( 0.1 ) - , m_Signif( 2.5 ) - , m_DMK0( 10.0 ) - , m_TseedVtxMin( 1.0 ) - , m_TseedVtxMax( 200.0 ) - , m_TseedVtxMinAnyPV( 0.1 ) - , m_DtrakMax( 0.5 ) - , m_PtSeedsMin( 1000. ) - , m_PtMergedSeedsMin( 1000. ) - , m_SeedsMaxChi2DoF( 50. ) - , m_Triplets( false ) - , m_DRmin( 0.1 ) - , m_DRmax( 10. ) - , m_TrkChi2DoF( 2.5 ) - , m_DeltaRSeeds( 0.5 ) - , m_preFilter( true ) - , m_jetFilter( true ) - , m_PVveto( true ) { - // - declareInterface<IJetMaker>( this ); - // - declareProperty( "SeedID", m_seedID, "Particle ID for the Seed" ); - declareProperty( "SeedRParameter", m_r ); - declareProperty( "Sort", m_sort, - "Sorting Criteria for jets [0:none,1:pt,2:E,3:eta, 4:ProtoSeed the Dau then Pt, default:4]" ); - // define momentum combiner - declareProperty( "ParticleCombiner", m_combinerName ); - declareProperty( "SeedPtTrackMin", m_PtTrackMin, "pt of the track used for Vertexing" ); - declareProperty( "SeedPTrackMin", m_PTrackMin, "p of the track used for Vertexing" ); - declareProperty( "SeedIPmin", m_IPmin, "ip of the track used for Vertexing" ); - declareProperty( "SeedSignif", m_Signif, "signif oft he track used for Vertexing" ); - declareProperty( "SeedDMK0", m_DMK0, "mass window for Ks" ); - declareProperty( "SeedTseedVtxMin", m_TseedVtxMin, "min distance btw PV and the vtx" ); - declareProperty( "SeedTseedVtxMax", m_TseedVtxMax, "max distance btw PV and the vtx" ); - declareProperty( "SeedTseedVtxMinAnyPV", m_TseedVtxMinAnyPV, "min distance btw any PV and the vtx" ); - - declareProperty( "SeedDtrakMax", m_DtrakMax, "dca window for vtx" ); - declareProperty( "SeedPtSeedsMin", m_PtSeedsMin, "min pt of the seeds" ); - declareProperty( "SeedVtxMaxChi2PerDoF", m_SeedsMaxChi2DoF, "max chi2 per dof for the vtx fit of the seed" ); - declareProperty( "PtMergedSeedsMin", m_PtMergedSeedsMin, "min pt of the merged seeds" ); - - declareProperty( "SeedTriplets", m_Triplets, "built Vtx with 3 tracks" ); - declareProperty( "SeedDRmin", m_DRmin, "min positon in R of the vtx" ); - declareProperty( "SeedDRmax", m_DRmax, "max positon in R of the vtx" ); - declareProperty( "SeedTrkChi2PerDoF", m_TrkChi2DoF, "max chi2PerDoF for the track used for the vtx" ); - declareProperty( "SeedsDeltaR", m_DeltaRSeeds, "min DR btw 2 seeds otherwise:merge" ); - declareProperty( "SeedpreFilter", m_preFilter, "pre filter the protoseeds" ); - declareProperty( "SeedFilter", m_jetFilter, "filter the protoseeds" ); - declareProperty( "vetoPV", m_PVveto, - "exclude vertex near to any PV with distance fixe by 'SeedTseedVtxMinAnyPV' and exclude the trk " - "associated to any PV to construct a secondary vtx" ); - } - /** standard initialization of the tool - * @return status code - */ - StatusCode initialize() override; - - protected: - /// make the detailed check of all parameters - - inline StatusCode check() const { - if ( 0 > m_ptmin ) { Warning( "PtMin is negative " ).ignore(); } - return StatusCode::SUCCESS; - } - - int to_user_index( const int index ) const { return index + 10000; } - int from_user_index( const int index ) const { return index - 10000; } - // proposed jet ID - int m_seedID; ///< proposed jet ID - // R-parameter - double m_r; ///< R-parameters - // ptMin-parameter - double m_ptmin; ///< pt-min parameter - // jet sorting criteria - int m_sort; ///< jet sorting criteria - // combiner - std::string m_combinerName; - mutable IParticleCombiner* m_combiner; ///< combiner to be used - - IDistanceCalculator* m_dist; - IVertexFit* m_fitter; - - double m_Rmax; - double m_PtTrackMin; - double m_PTrackMin; - double m_IPmin; - double m_Signif; - double m_DMK0; - double m_TseedVtxMin; - double m_TseedVtxMax; - double m_TseedVtxMinAnyPV; - double m_DtrakMax; - double m_PtSeedsMin; - double m_PtMergedSeedsMin; - double m_SeedsMaxChi2DoF; - bool m_Triplets; - double m_DRmin; - double m_DRmax; - double m_TrkChi2DoF; - double m_DeltaRSeeds; - bool m_preFilter; - bool m_jetFilter; - bool m_PVveto; - - IJetMaker::Jets JetCone( const double& Rmax, IJetMaker::Jets& seeds, const IJetMaker::Input& inputs ) const; - - IJetMaker::Jets FilterProtoJets( const double& DeltaRSeeds, IJetMaker::Jets& ProtoJets ) const; - - double getDeltaR( LHCb::Particle* p1, LHCb::Particle* p2 ) const; - - void RemoveTracks( LHCb::Particle::ConstVector& particles, const LHCb::RecVertex PV ) const; - - private: - mutable Gaudi::Accumulators::StatCounter<> m_nSeeds{this, "#seeds"}; - }; - -} // namespace LoKi - -class sortDauPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - if ( obj1->weight() == obj2->weight() ) - if ( obj1->daughtersVector().size() == obj2->daughtersVector().size() ) - return obj1->pt() > obj2->pt(); - else - return obj1->daughtersVector().size() > obj2->daughtersVector().size(); - else - return obj1->weight() > obj2->weight(); - } -}; - -class sortE { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().E() > obj2->momentum().E(); - } -}; - -class sortEta { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { - return obj1->momentum().Eta() > obj2->momentum().Eta(); - } -}; - -class sortPt { -public: - inline bool operator()( LHCb::Particle* obj1, LHCb::Particle* obj2 ) { return obj1->pt() > obj2->pt(); } -};