diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index a713a510a7..4421cc7e6d 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -5,8 +5,6 @@ on: [push, pull_request] jobs: format: runs-on: ubuntu-latest - strategy: - fail-fast: false steps: - uses: actions/checkout@v3 with: @@ -15,21 +13,6 @@ jobs: - name: Start container run: | docker run -it --name CI_container -v ${GITHUB_WORKSPACE}:/Package -v /cvmfs:/cvmfs:shared -d ghcr.io/aidasoft/centos7:latest /bin/bash - - name: Setup Spack - run: | - docker exec CI_container /bin/bash -c 'mkdir Spack - cd Spack - git clone https://github.com/key4hep/spack - git clone https://github.com/key4hep/key4hep-spack' - - name: Setup LLVM - run: | - docker exec CI_container /bin/bash -c 'cd ./Package - source ../Spack/spack/share/spack/setup-env.sh - source ../Spack/key4hep-spack/environments/key4hep-release-user/setup_clingo_centos7.sh - spack env activate ../Spack/key4hep-spack/environments/key4hep-release-user - spack add llvm - spack concretize -f' - # spack install' - name: Add upstream run: | docker exec CI_container /bin/bash -c 'cd ./Package @@ -38,12 +21,9 @@ jobs: - name: Run formatter run: | docker exec CI_container /bin/bash -c 'cd ./Package - source ../Spack/spack/share/spack/setup-env.sh - source ../Spack/key4hep-spack/environments/key4hep-release-user/setup_clingo_centos7.sh - spack env activate ../Spack/key4hep-spack/environments/key4hep-release-user - spack load llvm + source /cvmfs/sft.cern.ch/lcg/contrib/clang/14.0.6/x86_64-centos7/setup.sh git clang-format --style=file $(git merge-base upstream/master HEAD)' - name: Check cleanliness run: | docker exec CI_container /bin/bash -c 'cd ./Package - git diff --exit-code' + git diff' diff --git a/.zenodo.json b/.zenodo.json index 18ce474087..a4efdc202a 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,7 +1,7 @@ { "creators": [ { - "affiliation": "CERN", + "affiliation": "Ecole Polytechnique Fédérale de Lausanne", "name": "Helsens, Clement", "orcid": "0000-0002-9243-7554" }, diff --git a/CMakeLists.txt b/CMakeLists.txt index 190ed6c612..0df276bd99 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,7 @@ include(CTest) #--- options ------------------------------------------------------------------ -set(WITH_ACTS AUTO CACHE STRING "Build analyzers that need Acts") +set(WITH_ACTS OFF CACHE STRING "Build analyzers that need Acts") set_property(CACHE WITH_ACTS PROPERTY STRINGS AUTO ON OFF) set(WITH_DD4HEP AUTO CACHE STRING "Build analyzers that need DD4hep") @@ -62,6 +62,8 @@ find_package(podio) set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) find_package(FastJet) +find_package( Delphes REQUIRED ) + if(WITH_ACTS) find_package( Acts COMPONENTS Core ) @@ -88,6 +90,17 @@ if(WITH_DD4HEP) endif() endif() +if(WITH_ONNX) + find_package(ONNXRuntime) + if(ONNXRuntime_FOUND) + + elseif(WITH_ONNX STREQUAL AUTO) + message(WARNING "ONNXRuntime not found. Skipping ONNX-dependent analyzers.") + set(WITH_ONNX OFF) + else() + message(FATAL_ERROR "Failed to locate ONNXRuntime!") + endif() +endif() if(WITH_ONNX AND BUILD_TESTING) # currently these files are only needed by ONNX-parts # Grab the test files into a cached directory diff --git a/README.md b/README.md index 87b292525c..82bb375b1d 100644 --- a/README.md +++ b/README.md @@ -48,9 +48,18 @@ ROOT dataframe documentation is available ## Getting started -In order to use the FCC analysers within ROOT dataframe, a dictionary needs to -be built and put into `LD_LIBRARY_PATH` (this happens in `setup.sh`). The -following needs to be done when running local code and for developers. +In order to use the FCC analyzers within ROOT RDataFrame, a dictionary needs to +be built and put into `LD_LIBRARY_PATH`. In order to build and load FCCAnalyses +with default options one needs to run following two commands: + +```shell +source ./setup.sh +fccanalysis build +``` + +The FCCAnalyses is a CMake based project and any customizations can be provided +in classic CMake style, the following commands are equivalent to default version +of FCCAnalyses: ```shell source ./setup.sh @@ -65,6 +74,21 @@ cd .. > Each time changes are made in the C++ code, for example in > `analyzers/dataframe/` please do not forget to re-compile :) > +> To cleanly recompile the default version of FCCAnalyses one can use +> `fccanalysis build --clean-build`. + +In order to provide the possibility to keep developing an analysis with well +defined Key4hep stack, the sub-command `fccanalysis pin` is provided. One can +pin his/her analysis with +``` +source setup.sh +fccanalysis pin +``` + +To remove the pin run +``` +fccanalysis pin --clear +``` ## Generalities @@ -206,6 +230,16 @@ In an attempt to ease the development of new physics case studies, such as for t The preferred style of the C++ code in the FCCAnalyses is LLVM which is checked by CI job. -Currently `clang-format` is not available in the Key4hep stack, but one way to -obtain a recent version of it is through downloading of the +Currently `clang-format` is not available in the Key4hep stack, but one can +obtain a suitable version of it from CVMFS thanks to LCG: +``` +source /cvmfs/sft.cern.ch/lcg/contrib/clang/14.0.6/x86_64-centos7/setup.sh +``` + +Then to apply formatting to a given file: +``` +clang-format -i -style=file /path/to/file.cpp +``` + +Another way to obtain a recent version of `clang-format` is through downloading [Key4hep Spack instance](https://key4hep.github.io/key4hep-doc/spack-build-instructions-for-librarians/spack-setup.html#downloading-a-spack-instance). diff --git a/addons/FastJet/python/jetClusteringHelper.py b/addons/FastJet/python/jetClusteringHelper.py index 7b4753a6c8..5618974006 100644 --- a/addons/FastJet/python/jetClusteringHelper.py +++ b/addons/FastJet/python/jetClusteringHelper.py @@ -1,90 +1,90 @@ -import json -import ROOT - - -class ExclusiveJetClusteringHelper: - def __init__(self, coll, njets, tag=""): - - self.input_coll = coll - self.njets = njets - - self.tag = tag - if tag != "": - self.tag = "_{}".format(tag) - - part_px = "part{}_px".format(self.tag) - part_py = "part{}_py".format(self.tag) - part_pz = "part{}_pz".format(self.tag) - part_e = "part{}_e".format(self.tag) - part_m = "part{}_m".format(self.tag) - part_q = "part{}_q".format(self.tag) - - pjetc = "pjetc{}".format(self.tag) - - _jet = "_jet{}".format(self.tag) - jet = "jet{}".format(self.tag) - _jetc = "_jetc{}".format(self.tag) - jetc = "jetc{}".format(self.tag) - - # compute jet observables - - observables = ["p", "e", "mass", "phi", "theta", "nconst"] - - self.jet_obs = dict() - for obs in observables: - self.jet_obs[obs] = "jet_{}{}".format(obs, self.tag) - event_njet = "event_njet{}".format(self.tag) - - self.jets = jet - self.constituents = jetc - - self.definition = dict() - - # get single particle properties - self.definition[part_px] = "ReconstructedParticle::get_px({})".format(self.input_coll) - self.definition[part_py] = "ReconstructedParticle::get_py({})".format(self.input_coll) - self.definition[part_pz] = "ReconstructedParticle::get_pz({})".format(self.input_coll) - self.definition[part_e] = "ReconstructedParticle::get_e({})".format(self.input_coll) - self.definition[part_m] = "ReconstructedParticle::get_mass({})".format(self.input_coll) - self.definition[part_q] = "ReconstructedParticle::get_charge({})".format(self.input_coll) - - # form fastjet pseudo jets - self.definition[pjetc] = "JetClusteringUtils::set_pseudoJets({}, {}, {}, {})".format( - part_px, part_py, part_pz, part_e - ) - - # run jet clustering with all reconstructed particles. ee_kt_algorithm, R=1.5, inclusive clustering, E-scheme - self.definition[_jet] = "JetClustering::clustering_ee_kt(2, {}, 1, 0)({})".format(njets, pjetc) - - # get the jets out of the struct - self.definition[jet] = "JetClusteringUtils::get_pseudoJets({})".format(_jet) - - # get the jets constituents out of the struct - self.definition[_jetc] = "JetClusteringUtils::get_constituents({})".format(_jet) - - # get constituents - self.definition[jetc] = "JetConstituentsUtils::build_constituents_cluster({}, {})".format( - self.input_coll, _jetc - ) - - # compute jet observables - self.definition[self.jet_obs["p"]] = "JetClusteringUtils::get_p({})".format(self.jets) - self.definition[self.jet_obs["e"]] = "JetClusteringUtils::get_e({})".format(self.jets) - self.definition[self.jet_obs["mass"]] = "JetClusteringUtils::get_m({})".format(self.jets) - self.definition[self.jet_obs["phi"]] = "JetClusteringUtils::get_phi({})".format(self.jets) - self.definition[self.jet_obs["theta"]] = "JetClusteringUtils::get_theta({})".format(self.jets) - self.definition[self.jet_obs["nconst"]] = "JetConstituentsUtils::count_consts({})".format(self.constituents) - self.definition[event_njet] = "JetConstituentsUtils::count_jets({})".format(self.constituents) - - def define(self, df): - - for var, call in self.definition.items(): - df = df.Define(var, call) - - return df - - def outputBranches(self): - - out = list(self.jet_obs.values()) - out += [obs for obs in self.definition.keys() if "event_" in obs] - return out +import json +import ROOT + + +class ExclusiveJetClusteringHelper: + def __init__(self, coll, njets, tag=""): + + self.input_coll = coll + self.njets = njets + + self.tag = tag + if tag != "": + self.tag = "_{}".format(tag) + + part_px = "part{}_px".format(self.tag) + part_py = "part{}_py".format(self.tag) + part_pz = "part{}_pz".format(self.tag) + part_e = "part{}_e".format(self.tag) + part_m = "part{}_m".format(self.tag) + part_q = "part{}_q".format(self.tag) + + pjetc = "pjetc{}".format(self.tag) + + _jet = "_jet{}".format(self.tag) + jet = "jet{}".format(self.tag) + _jetc = "_jetc{}".format(self.tag) + jetc = "jetc{}".format(self.tag) + + # compute jet observables + + observables = ["p", "e", "mass", "phi", "theta", "nconst"] + + self.jet_obs = dict() + for obs in observables: + self.jet_obs[obs] = "jet_{}{}".format(obs, self.tag) + event_njet = "event_njet{}".format(self.tag) + + self.jets = jet + self.constituents = jetc + + self.definition = dict() + + # get single particle properties + self.definition[part_px] = "ReconstructedParticle::get_px({})".format(self.input_coll) + self.definition[part_py] = "ReconstructedParticle::get_py({})".format(self.input_coll) + self.definition[part_pz] = "ReconstructedParticle::get_pz({})".format(self.input_coll) + self.definition[part_e] = "ReconstructedParticle::get_e({})".format(self.input_coll) + self.definition[part_m] = "ReconstructedParticle::get_mass({})".format(self.input_coll) + self.definition[part_q] = "ReconstructedParticle::get_charge({})".format(self.input_coll) + + # form fastjet pseudo jets + self.definition[pjetc] = "JetClusteringUtils::set_pseudoJets({}, {}, {}, {})".format( + part_px, part_py, part_pz, part_e + ) + + # run jet clustering with all reconstructed particles. ee_kt_algorithm, R=1.5, inclusive clustering, E-scheme + self.definition[_jet] = "JetClustering::clustering_ee_kt(2, {}, 1, 0)({})".format(njets, pjetc) + + # get the jets out of the struct + self.definition[jet] = "JetClusteringUtils::get_pseudoJets({})".format(_jet) + + # get the jets constituents out of the struct + self.definition[_jetc] = "JetClusteringUtils::get_constituents({})".format(_jet) + + # get constituents + self.definition[jetc] = "JetConstituentsUtils::build_constituents_cluster({}, {})".format( + self.input_coll, _jetc + ) + + # compute jet observables + self.definition[self.jet_obs["p"]] = "JetClusteringUtils::get_p({})".format(self.jets) + self.definition[self.jet_obs["e"]] = "JetClusteringUtils::get_e({})".format(self.jets) + self.definition[self.jet_obs["mass"]] = "JetClusteringUtils::get_m({})".format(self.jets) + self.definition[self.jet_obs["phi"]] = "JetClusteringUtils::get_phi({})".format(self.jets) + self.definition[self.jet_obs["theta"]] = "JetClusteringUtils::get_theta({})".format(self.jets) + self.definition[self.jet_obs["nconst"]] = "JetConstituentsUtils::count_consts({})".format(self.constituents) + self.definition[event_njet] = "JetConstituentsUtils::count_jets({})".format(self.constituents) + + def define(self, df): + + for var, call in self.definition.items(): + df = df.Define(var, call) + + return df + + def outputBranches(self): + + out = list(self.jet_obs.values()) + out += [obs for obs in self.definition.keys() if "event_" in obs] + return out diff --git a/analyzers/dataframe/CMakeLists.txt b/analyzers/dataframe/CMakeLists.txt index a758321723..1c5cfb32d3 100644 --- a/analyzers/dataframe/CMakeLists.txt +++ b/analyzers/dataframe/CMakeLists.txt @@ -6,14 +6,19 @@ find_package(Vdt) message(STATUS "includes-------------------------- dataframe edm4hep: ${EDM4HEP_INCLUDE_DIRS}") message(STATUS "includes-------------------------- dataframe podio : ${podio_INCLUDE_DIR}") +message(STATUS "includes-------------------------- dataframe delphes: ${DELPHES_INCLUDE_DIR}") +message(STATUS "includes-------------------------- dataframe delphes EXt TrkCov: ${DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR}") +message(STATUS "includes-------------------------- dataframe delphes EXt: ${DELPHES_EXTERNALS_INCLUDE_DIR}") + +include_directories(${DELPHES_INCLUDE_DIR} + ${DELPHES_EXTERNALS_INCLUDE_DIR} + ${DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR} + ) file(GLOB sources src/*.cc) file(GLOB headers RELATIVE ${CMAKE_CURRENT_LIST_DIR} FCCAnalyses/*.h) -message(STATUS "includes headers ${headers}") -message(STATUS "includes sources ${sources}") - list(FILTER headers EXCLUDE REGEX "LinkDef.h") if(NOT WITH_DD4HEP) list(FILTER headers EXCLUDE REGEX "CaloNtupleizer.h") @@ -22,6 +27,8 @@ endif() if(NOT WITH_ONNX) list(FILTER headers EXCLUDE REGEX "JetFlavourUtils.h") list(FILTER sources EXCLUDE REGEX "JetFlavourUtils.cc") + list(FILTER headers EXCLUDE REGEX "WeaverUtils.h") + list(FILTER sources EXCLUDE REGEX "WeaverUtils.cc") endif() if(NOT WITH_ACTS) @@ -31,6 +38,8 @@ if(NOT WITH_ACTS) list(FILTER sources EXCLUDE REGEX "VertexFinderActs.cc") endif() +message(STATUS "includes headers ${headers}") +message(STATUS "includes sources ${sources}") message(STATUS "CMAKE_CURRENT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}") message(STATUS "CMAKE_INSTALL_INCLUDEDIR ${CMAKE_INSTALL_INCLUDEDIR}") @@ -41,8 +50,13 @@ target_include_directories(FCCAnalyses PUBLIC $ $ ${VDT_INCLUDE_DIR} + ${DELPHES_INCLUDE_DIR} + ${DELPHES_EXTERNALS_INCLUDE_DIR} + ${DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR} ) +message(STATUS " ====== DELPHES LIBRARY = " ${DELPHES_LIBRARY} ) +message(STATUS " ====== DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR = " ${DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR} ) target_link_libraries(FCCAnalyses @@ -53,6 +67,7 @@ target_link_libraries(FCCAnalyses EDM4HEP::edm4hep EDM4HEP::edm4hepDict podio::podio + ${DELPHES_LIBRARY} ${ADDONS_LIBRARIES} gfortran # todo: why necessary? ) diff --git a/analyzers/dataframe/FCCAnalyses/JetTaggingUtils.h b/analyzers/dataframe/FCCAnalyses/JetTaggingUtils.h index ba457cd43f..ac707c2739 100644 --- a/analyzers/dataframe/FCCAnalyses/JetTaggingUtils.h +++ b/analyzers/dataframe/FCCAnalyses/JetTaggingUtils.h @@ -1,45 +1,60 @@ -#ifndef JETTAGGINGUTILS_ANALYZERS_H -#define JETTAGGINGUTILS_ANALYZERS_H - -#include -#include "Math/Vector4D.h" -#include "ROOT/RVec.hxx" -#include "edm4hep/MCParticleData.h" -#include "fastjet/JetDefinition.hh" -#include "TRandom3.h" - -/** Jet tagging utilities interface. -This represents a set functions and utilities to perfom jet tagging from a list of jets. -*/ -namespace FCCAnalyses{ - -namespace JetTaggingUtils{ - - /** @name JetTaggingUtils - * Jet tagging interface utilities. - */ - - //Get flavour association of jet - ROOT::VecOps::RVec get_flavour(ROOT::VecOps::RVec in, ROOT::VecOps::RVec MCin); - //Get b-tags with an efficiency applied - ROOT::VecOps::RVec get_btag(ROOT::VecOps::RVec in, float efficiency, float mistag_c=0., float mistag_l=0., float mistag_g=0.); - //Get c-tags with an efficiency applied - ROOT::VecOps::RVec get_ctag(ROOT::VecOps::RVec in, float efficiency, float mistag_b=0., float mistag_l=0., float mistag_g=0.); - //Get l-tags with an efficiency applied - ROOT::VecOps::RVec get_ltag(ROOT::VecOps::RVec in, float efficiency, float mistag_b=0., float mistag_c=0., float mistag_g=0.); - //Get g-tags with an efficiency applied - ROOT::VecOps::RVec get_gtag(ROOT::VecOps::RVec in, float efficiency, float mistag_b=0., float mistag_c=0., float mistag_l=0.); - - /// select a list of jets depending on the status of a certain boolean flag (corresponding to its tagging state) - struct sel_tag { - bool m_pass; // if pass is true, select tagged jets. Otherwise select anti-tagged ones - sel_tag(bool arg_pass); - ROOT::VecOps::RVec operator() (ROOT::VecOps::RVec tags, ROOT::VecOps::RVec in); - }; - - ///@} -}//end NS JetTaggingUtils - -}//end NS FCCAnalyses - -#endif +#ifndef JETTAGGINGUTILS_ANALYZERS_H +#define JETTAGGINGUTILS_ANALYZERS_H + +#include "Math/Vector4D.h" +#include "ROOT/RVec.hxx" +#include "TRandom3.h" +#include "edm4hep/MCParticleData.h" +#include "fastjet/JetDefinition.hh" +#include + +/** Jet tagging utilities interface. +This represents a set functions and utilities to perfom jet tagging from a list +of jets. +*/ +namespace FCCAnalyses { + +namespace JetTaggingUtils { + +/** @name JetTaggingUtils + * Jet tagging interface utilities. + */ + +// Get flavour association of jet +ROOT::VecOps::RVec +get_flavour(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec MCin); +// Get b-tags with an efficiency applied +ROOT::VecOps::RVec get_btag(ROOT::VecOps::RVec in, float efficiency, + float mistag_c = 0., float mistag_l = 0., + float mistag_g = 0.); +// Get c-tags with an efficiency applied +ROOT::VecOps::RVec get_ctag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b = 0., float mistag_l = 0., + float mistag_g = 0.); +// Get l-tags with an efficiency applied +ROOT::VecOps::RVec get_ltag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b = 0., float mistag_c = 0., + float mistag_g = 0.); +// Get g-tags with an efficiency applied +ROOT::VecOps::RVec get_gtag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b = 0., float mistag_c = 0., + float mistag_l = 0.); + +/// select a list of jets depending on the status of a certain boolean flag +/// (corresponding to its tagging state) +struct sel_tag { + bool m_pass; // if pass is true, select tagged jets. Otherwise select + // anti-tagged ones + sel_tag(bool arg_pass); + ROOT::VecOps::RVec + operator()(ROOT::VecOps::RVec tags, + ROOT::VecOps::RVec in); +}; + +///@} +} // namespace JetTaggingUtils + +} // namespace FCCAnalyses + +#endif diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h index 637333d0cb..48f76453c0 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h @@ -6,8 +6,11 @@ #include #include "ROOT/RVec.hxx" +#include "edm4hep/Quantity.h" #include "edm4hep/ReconstructedParticleData.h" +#include "edm4hep/TrackData.h" #include "edm4hep/TrackState.h" +#include "edm4hep/TrackerHitData.h" #include #include #include diff --git a/analyzers/dataframe/FCCAnalyses/SmearObjects.h b/analyzers/dataframe/FCCAnalyses/SmearObjects.h index b414974374..eb2d2e498c 100644 --- a/analyzers/dataframe/FCCAnalyses/SmearObjects.h +++ b/analyzers/dataframe/FCCAnalyses/SmearObjects.h @@ -4,63 +4,98 @@ #include #include +#include "FCCAnalyses/ReconstructedParticle2Track.h" +#include "ROOT/RVec.hxx" #include "TLorentzVector.h" #include "TMatrixDSym.h" #include "TRandom.h" -#include -#include "ROOT/RVec.hxx" #include "edm4hep/MCParticleData.h" +#include -#include "FCCAnalyses/ReconstructedParticle2Track.h" +namespace FCCAnalyses { + +namespace SmearObjects { -namespace FCCAnalyses -{ +/// for a given MC particle, returns a "track state", i.e. a vector of 5 helix +/// parameters, in Delphes convention +TVectorD TrackParamFromMC_DelphesConv(edm4hep::MCParticleData aMCParticle); - namespace SmearObjects - { +/// generates new track states, by rescaling the covariance matrix of the tracks +struct SmearedTracks { + bool m_debug; + TRandom m_random; + float m_smear_parameters[5]; + SmearedTracks(float smear_d0, float smear_phi, float smear_omega, + float smear_z0, float smear_tlambda, bool debug); + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &alltracks, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles); +}; - /// for a given MC particle, returns a "track state", i.e. a vector of 5 helix parameters, in Delphes convention - TVectorD TrackParamFromMC_DelphesConv(edm4hep::MCParticleData aMCParticle); +/// used to validate the method above. Stores the "MC-truth track states", in a +/// collection that runs parallel to the full collection of tracks. +ROOT::VecOps::RVec mcTrackParameters( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &alltracks, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles); - /// generates new track states, by rescaling the covariance matrix of the tracks - struct SmearedTracks - { - bool m_debug; - TRandom m_random; - float m_smear_parameters[5]; - SmearedTracks(float smear_d0, float smear_phi, float smear_omega, float smear_z0, float smear_tlambda, bool debug); - ROOT::VecOps::RVec operator()( - const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &alltracks, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles); - }; +/// generates random values for a vector, given the covariance matrix of its +/// components, using a Choleski decomposition. Code from Franco Bedeschi +TVectorD CovSmear(TVectorD x, TMatrixDSym C, TRandom *ran, bool debug); - /// used to validate the method above. Stores the "MC-truth track states", in a collection that runs parallel to the full collection of tracks. - ROOT::VecOps::RVec mcTrackParameters(const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &alltracks, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles); +/// generates new track dNdx, by rescaling the poisson error of the cluster +/// count +struct SmearedTracksdNdx { + bool m_debug; + TRandom m_random; + float m_scale; + SmearedTracksdNdx(float m_scale, bool debug); + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &dNdx, + const ROOT::VecOps::RVec &length, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles); +}; - /// generates random values for a vector, given the covariance matrix of its components, using a Choleski decomposition. Code from Franco Bedeschi - TVectorD CovSmear(TVectorD x, TMatrixDSym C, TRandom *ran, bool debug); +/// generates new tracker hits, by rescaling the timing measurement +struct SmearedTracksTOF { + bool m_debug; + TRandom m_random; + float m_scale; + SmearedTracksTOF(float m_scale, bool debug); + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &trackdata, + const ROOT::VecOps::RVec &trackerhits, + const ROOT::VecOps::RVec &length, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles); +}; - /// generates new reco particles, smeared by given parameters - struct SmearedReconstructedParticle - { - bool m_debug; - float m_scale; - int m_type; - int m_mode; - SmearedReconstructedParticle(float scale, int type, int mode, bool debug); +/// generates new reco particles, smeared by given parameters +struct SmearedReconstructedParticle { + bool m_debug; + float m_scale; + int m_type; + int m_mode; + SmearedReconstructedParticle(float scale, int type, int mode, bool debug); - ROOT::VecOps::RVec operator()( - const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles); - }; + ROOT::VecOps::RVec + operator()(const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles); +}; - } -} +} // namespace SmearObjects +} // namespace FCCAnalyses #endif diff --git a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h index b397a07036..fa408a9ad3 100644 --- a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h +++ b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h @@ -22,6 +22,9 @@ #include "edm4hep/VertexData.h" #include "edm4hep/Vertex.h" +#include "VertexFit.h" // from Delphes - updates Franco, Jul 2022 +#include "VertexMore.h" + /** Vertex interface using Franco Bedeshi's code. This represents a set functions and utilities to perfom vertexing from a list of tracks. @@ -52,12 +55,10 @@ namespace VertexFitterSimple{ double bsc_x=0., double bsc_y=0., double bsc_z=0. ) ; /// Return the tracks that are flagged as coming from the primary vertex - ROOT::VecOps::RVec get_PrimaryTracks( VertexingUtils::FCCAnalysesVertex initialVertex, - ROOT::VecOps::RVec tracks, + ROOT::VecOps::RVec get_PrimaryTracks( ROOT::VecOps::RVec tracks, bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, - double bsc_x, double bsc_y, double bsc_z, - int ipass = 0 ) ; + double bsc_x, double bsc_y, double bsc_z ) ; /// Return the tracks that are NOT flagged as coming from the primary vertex @@ -69,7 +70,7 @@ namespace VertexFitterSimple{ ROOT::VecOps::RVec primaryTracks ) ; - +/* Double_t FastRv(TVectorD p1, TVectorD p2) ; TMatrixDSym RegInv3(TMatrixDSym &Smat0) ; TMatrixD Fill_A(TVectorD par, Double_t phi) ; @@ -80,6 +81,12 @@ namespace VertexFitterSimple{ TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q); TVector3 ParToP(TVectorD Par); + TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q); + TVector3 ParToP(TVectorD Par); +*/ + + + }//end NS VertexFitterSimple }//end NS FCCAnalyses diff --git a/analyzers/dataframe/FCCAnalyses/VertexingUtils.h b/analyzers/dataframe/FCCAnalyses/VertexingUtils.h index 9b6436db02..a7f3b4c300 100644 --- a/analyzers/dataframe/FCCAnalyses/VertexingUtils.h +++ b/analyzers/dataframe/FCCAnalyses/VertexingUtils.h @@ -19,12 +19,20 @@ #include "fastjet/JetDefinition.hh" + /** Vertexing utilities */ namespace FCCAnalyses{ namespace VertexingUtils{ + /// from delphes: returns track state parameters (delphes convention) for a given vertex (x), momentum (p) and charge + TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q); + + /// from delphes: returns the momentum corresponding to a given track state + TVector3 ParToP(TVectorD Par); + + /// Structure to keep useful track information that is related to the vertex struct FCCAnalysesVertex{ edm4hep::VertexData vertex; @@ -303,14 +311,30 @@ namespace VertexingUtils{ ROOT::VecOps::RVec> get_position_SV( ROOT::VecOps::RVec> vertices ); // --- for get_SV_jets --- // - // --- Internal methods needed by the code of Franco B : float get_trackMom( edm4hep::TrackState & atrack ); - TVectorD get_trackParam( edm4hep::TrackState & atrack) ; - TMatrixDSym get_trackCov( edm4hep::TrackState & atrack) ; + + +// --- Conversion methods between the Delphes and edm4hep conventions + +/// convert track parameters, from edm4hep to delphes conventions + TVectorD Edm4hep2Delphes_TrackParam( const TVectorD& param, bool Units_mm ); +/// convert track parameters, from delphes to edm4hep conventions + TVectorD Delphes2Edm4hep_TrackParam( const TVectorD& param, bool Units_mm ); +/// convert track covariance matrix, from edm4hep to delphes conventions + TMatrixDSym Edm4hep2Delphes_TrackCovMatrix( const std::array& covMatrix, bool Units_mm ); +/// convert track covariance matrix, from delphes to edm4hep conventions + std::array Delphes2Edm4hep_TrackCovMatrix( const TMatrixDSym& cov, bool Units_mm ) ; + + + /// --- Internal methods needed by the code of Franco B: + TVectorD get_trackParam( edm4hep::TrackState & atrack, bool Units_mm = false) ; + TMatrixDSym get_trackCov( edm4hep::TrackState & atrack, bool Units_mm = false) ; TVectorD ParToACTS(TVectorD Par); TMatrixDSym CovToACTS(TMatrixDSym Cov,TVectorD Par); + + }//end NS VertexingUtils }//end NS FCCAnalyses diff --git a/analyzers/dataframe/src/JetTaggingUtils.cc b/analyzers/dataframe/src/JetTaggingUtils.cc index 2600cd3d3f..7ed529c18a 100644 --- a/analyzers/dataframe/src/JetTaggingUtils.cc +++ b/analyzers/dataframe/src/JetTaggingUtils.cc @@ -1,148 +1,157 @@ -#include "FCCAnalyses/JetTaggingUtils.h" - -namespace FCCAnalyses{ - -namespace JetTaggingUtils{ - -ROOT::VecOps::RVec get_flavour(ROOT::VecOps::RVec in, - ROOT::VecOps::RVec MCin) -{ - ROOT::VecOps::RVec result(in.size(),0); - - int loopcount =0; - for (size_t i = 0; i < MCin.size(); ++i) { - auto & parton = MCin[i]; - //Select partons only (for pythia8 71-79, for pythia6 2): - if ((parton.generatorStatus>80 || - parton.generatorStatus<70) && - parton.generatorStatus != 2 ) continue; - if (std::abs(parton.PDG) > 5 && parton.PDG!=21) continue; - ROOT::Math::PxPyPzMVector lv(parton.momentum.x, parton.momentum.y, - parton.momentum.z, parton.mass); - - for (size_t j = 0; j < in.size(); ++j) { - auto & p = in[j]; - //float dEta = lv.Eta() - p.eta(); - //float dPhi = lv.Phi() - p.phi(); - //float deltaR = sqrt(dEta*dEta+dPhi*dPhi); - //if (deltaR <= 0.5 && gRandom->Uniform() <= efficiency) result[j] = true; - - Float_t dot = p.px() * parton.momentum.x - + p.py() * parton.momentum.y - + p.pz() * parton.momentum.z; - Float_t lenSq1 = p.px() * p.px() - + p.py() * p.py() - + p.pz() * p.pz(); - Float_t lenSq2 = parton.momentum.x * parton.momentum.x - + parton.momentum.y * parton.momentum.y - + parton.momentum.z * parton.momentum.z; - Float_t norm = sqrt(lenSq1*lenSq2); - Float_t angle = acos(dot/norm); - - if (angle <= 0.3) { - if (result[j]==21 or result[j]==0) { - // if no match before, or matched to gluon, match to - // this particle (favour quarks over gluons) - result[j] = std::abs ( parton.PDG ); - } - else if (parton.PDG!=21) { - // if matched to quark, and this is a quark, favour - // heavier flavours - result[j] = std::max(result[j], std::abs ( parton.PDG )); - } else { - // if matched to quark, and this is a gluon, keep - // previous result (favour quark) - ; - } - } - - - } - } - - return result; -} - -ROOT::VecOps::RVec -get_btag(ROOT::VecOps::RVec in, - float efficiency, float mistag_c, - float mistag_l, float mistag_g) { - - ROOT::VecOps::RVec result(in.size(),0); - - for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 5 && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) result[j] = 1; - } - return result; -} - -ROOT::VecOps::RVec -get_ctag(ROOT::VecOps::RVec in, - float efficiency, float mistag_b, - float mistag_l, float mistag_g) { - - ROOT::VecOps::RVec result(in.size(),0); - - for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 4 && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) result[j] = 1; - } - return result; -} - -ROOT::VecOps::RVec -get_ltag(ROOT::VecOps::RVec in, - float efficiency, float mistag_b, - float mistag_c, float mistag_g) { - - ROOT::VecOps::RVec result(in.size(),0); - - for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) < 4 && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) result[j] = 1; - } - return result; -} - -ROOT::VecOps::RVec -get_gtag(ROOT::VecOps::RVec in, - float efficiency, float mistag_b, - float mistag_c, float mistag_l) { - - ROOT::VecOps::RVec result(in.size(),0); - - for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 21 && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) result[j] = 1; - } - return result; -} - -sel_tag::sel_tag(bool arg_pass): m_pass(arg_pass) {}; -ROOT::VecOps::RVec -sel_tag::operator()(ROOT::VecOps::RVec tags, - ROOT::VecOps::RVec in){ - ROOT::VecOps::RVec result; - for (size_t i = 0; i < in.size(); ++i) { - if (m_pass) { - if (tags.at(i)) result.push_back(in.at(i)); - } - else { - if (!tags.at(i)) result.push_back(in.at(i)); - } - } - return result; -} - -}//end NS JetTaggingUtils - -}//end NS FCCAnalyses +#include "FCCAnalyses/JetTaggingUtils.h" + +namespace FCCAnalyses { + +namespace JetTaggingUtils { + +ROOT::VecOps::RVec +get_flavour(ROOT::VecOps::RVec in, + ROOT::VecOps::RVec MCin) { + ROOT::VecOps::RVec result(in.size(), 0); + + int loopcount = 0; + for (size_t i = 0; i < MCin.size(); ++i) { + auto &parton = MCin[i]; + // Select partons only (for pythia8 71-79, for pythia6 2): + if ((parton.generatorStatus > 80 || parton.generatorStatus < 70) && + parton.generatorStatus != 2) + continue; + if (std::abs(parton.PDG) > 5 && parton.PDG != 21) + continue; + ROOT::Math::PxPyPzMVector lv(parton.momentum.x, parton.momentum.y, + parton.momentum.z, parton.mass); + + for (size_t j = 0; j < in.size(); ++j) { + auto &p = in[j]; + // float dEta = lv.Eta() - p.eta(); + // float dPhi = lv.Phi() - p.phi(); + // float deltaR = sqrt(dEta*dEta+dPhi*dPhi); + // if (deltaR <= 0.5 && gRandom->Uniform() <= efficiency) result[j] = + // true; + + Float_t dot = p.px() * parton.momentum.x + p.py() * parton.momentum.y + + p.pz() * parton.momentum.z; + Float_t lenSq1 = p.px() * p.px() + p.py() * p.py() + p.pz() * p.pz(); + Float_t lenSq2 = parton.momentum.x * parton.momentum.x + + parton.momentum.y * parton.momentum.y + + parton.momentum.z * parton.momentum.z; + Float_t norm = sqrt(lenSq1 * lenSq2); + Float_t angle = acos(dot / norm); + + if (angle <= 0.3) { + if (result[j] == 21 or result[j] == 0) { + // if no match before, or matched to gluon, match to + // this particle (favour quarks over gluons) + result[j] = std::abs(parton.PDG); + } else if (parton.PDG != 21) { + // if matched to quark, and this is a quark, favour + // heavier flavours + result[j] = std::max(result[j], std::abs(parton.PDG)); + } else { + // if matched to quark, and this is a gluon, keep + // previous result (favour quark) + ; + } + } + } + } + + return result; +} + +ROOT::VecOps::RVec get_btag(ROOT::VecOps::RVec in, float efficiency, + float mistag_c, float mistag_l, + float mistag_g) { + + ROOT::VecOps::RVec result(in.size(), 0); + + for (size_t j = 0; j < in.size(); ++j) { + if (in.at(j) == 5 && gRandom->Uniform() <= efficiency) + result[j] = 1; + if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + result[j] = 1; + if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + result[j] = 1; + if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + result[j] = 1; + } + return result; +} + +ROOT::VecOps::RVec get_ctag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b, float mistag_l, + float mistag_g) { + + ROOT::VecOps::RVec result(in.size(), 0); + + for (size_t j = 0; j < in.size(); ++j) { + if (in.at(j) == 4 && gRandom->Uniform() <= efficiency) + result[j] = 1; + if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + result[j] = 1; + if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + result[j] = 1; + if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + result[j] = 1; + } + return result; +} + +ROOT::VecOps::RVec get_ltag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b, float mistag_c, + float mistag_g) { + + ROOT::VecOps::RVec result(in.size(), 0); + + for (size_t j = 0; j < in.size(); ++j) { + if (in.at(j) < 4 && gRandom->Uniform() <= efficiency) + result[j] = 1; + if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + result[j] = 1; + if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + result[j] = 1; + if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + result[j] = 1; + } + return result; +} + +ROOT::VecOps::RVec get_gtag(ROOT::VecOps::RVec in, float efficiency, + float mistag_b, float mistag_c, + float mistag_l) { + + ROOT::VecOps::RVec result(in.size(), 0); + + for (size_t j = 0; j < in.size(); ++j) { + if (in.at(j) == 21 && gRandom->Uniform() <= efficiency) + result[j] = 1; + if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + result[j] = 1; + if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + result[j] = 1; + if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + result[j] = 1; + } + return result; +} + +sel_tag::sel_tag(bool arg_pass) : m_pass(arg_pass){}; +ROOT::VecOps::RVec +sel_tag::operator()(ROOT::VecOps::RVec tags, + ROOT::VecOps::RVec in) { + ROOT::VecOps::RVec result; + for (size_t i = 0; i < in.size(); ++i) { + if (m_pass) { + if (tags.at(i)) + result.push_back(in.at(i)); + } else { + if (!tags.at(i)) + result.push_back(in.at(i)); + } + } + return result; +} + +} // namespace JetTaggingUtils + +} // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/SmearObjects.cc b/analyzers/dataframe/src/SmearObjects.cc index 63a2029564..b1d3e664bb 100644 --- a/analyzers/dataframe/src/SmearObjects.cc +++ b/analyzers/dataframe/src/SmearObjects.cc @@ -2,432 +2,640 @@ #include "FCCAnalyses/VertexFitterSimple.h" #include "FCCAnalyses/VertexingUtils.h" - #include "TDecompChol.h" #include -namespace FCCAnalyses -{ +namespace FCCAnalyses { - namespace SmearObjects - { +namespace SmearObjects { - // ------------------------------------------------------------------------------------------- +// ------------------------------------------------------------------------------------------- - TVectorD TrackParamFromMC_DelphesConv(edm4hep::MCParticleData aMCParticle) - { +TVectorD TrackParamFromMC_DelphesConv(edm4hep::MCParticleData aMCParticle) { - TVector3 p(aMCParticle.momentum.x, aMCParticle.momentum.y, aMCParticle.momentum.z); - TVector3 x(1e-3 * aMCParticle.vertex.x, 1e-3 * aMCParticle.vertex.y, 1e-3 * aMCParticle.vertex.z); // mm to m - float Q = aMCParticle.charge; - TVectorD param = VertexFitterSimple::XPtoPar(x, p, Q); // convention Franco + TVector3 p(aMCParticle.momentum.x, aMCParticle.momentum.y, + aMCParticle.momentum.z); + TVector3 x(1e-3 * aMCParticle.vertex.x, 1e-3 * aMCParticle.vertex.y, + 1e-3 * aMCParticle.vertex.z); // mm to m + float Q = aMCParticle.charge; + TVectorD param = VertexingUtils::XPtoPar(x, p, Q); // convention Franco - return param; - } + return param; +} + +// ------------------------------------------------------------------------------------------- + +SmearedTracks::SmearedTracks(float smear_d0, float smear_phi, float smear_omega, + float smear_z0, float smear_tlambda, + bool debug = false) { + + m_smear_parameters[0] = smear_d0; + m_smear_parameters[1] = smear_phi; + m_smear_parameters[2] = smear_omega; + m_smear_parameters[3] = smear_z0; + m_smear_parameters[4] = smear_tlambda; + + m_debug = debug; +} + +ROOT::VecOps::RVec SmearedTracks::operator()( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &alltracks, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles) { + + // returns a vector of TrackStates that is parallel to the collection of full + // tracks (alltracks), i.e. same number of entries, same order. The method + // retrieve the MC particle that is associated to a track, and builds a "track + // state" out of the MC particle (i.e. it determines the corresponding (d0, + // phi, omega, z0, tanlambda). From this vector, and from the covariance + // matrix of the track, which is scaled by the user, new track states are + // generated. + + int ntracks = alltracks.size(); + + ROOT::VecOps::RVec result; + result.resize(ntracks); - // ------------------------------------------------------------------------------------------- + edm4hep::TrackState dummy; - SmearedTracks::SmearedTracks(float smear_d0, float smear_phi, float smear_omega, float smear_z0, float smear_tlambda, bool debug = false) - { + TVectorD zero(5); + for (int k = 0; k < 5; k++) { + zero(k) = 0.; + } - m_smear_parameters[0] = smear_d0; - m_smear_parameters[1] = smear_phi; - m_smear_parameters[2] = smear_omega; - m_smear_parameters[3] = smear_z0; - m_smear_parameters[4] = smear_tlambda; + for (int itrack = 0; itrack < ntracks; itrack++) { + edm4hep::TrackState track = alltracks[itrack]; + edm4hep::TrackState smeared_track = track; - m_debug = debug; + // find the corresponding MC particle + int MCindex = -1; + for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) { + edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; + int track_index = rp.tracks_begin; + if (track_index == itrack) { + MCindex = RP2MC_indices[ireco]; + break; + } + } // end loop on RPs + + if (MCindex < 0 || + MCindex >= + mcParticles + .size()) { // in principle, this should not happen in delphes, + // each track should be matched to a MC particle. + result[itrack] = dummy; + continue; } - ROOT::VecOps::RVec SmearedTracks::operator()( - const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &alltracks, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles) - { - - // returns a vector of TrackStates that is parallel to the collection of full tracks (alltracks), - // i.e. same number of entries, same order. - // The method retrieve the MC particle that is associated to a track, and builds a "track state" - // out of the MC particle (i.e. it determines the corresponding (d0, phi, omega, z0, tanlambda). - // From this vector, and from the covariance matrix of the track, which is scaled by the user, - // new track states are generated. - - int ntracks = alltracks.size(); - - ROOT::VecOps::RVec result; - result.resize(ntracks); - - edm4hep::TrackState dummy; - - TVectorD zero(5); - for (int k=0; k<5; k++) { zero(k) = 0. ; } - - for (int itrack = 0; itrack < ntracks; itrack++) - { - edm4hep::TrackState track = alltracks[itrack]; - edm4hep::TrackState smeared_track = track; - - // find the corresponding MC particle - int MCindex = -1; - for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) - { - edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; - int track_index = rp.tracks_begin; - if (track_index == itrack) - { - MCindex = RP2MC_indices[ireco]; - break; - } - } // end loop on RPs - - if (MCindex < 0 || MCindex >= mcParticles.size()) - { // in principle, this should not happen in delphes, - // each track should be matched to a MC particle. - result[itrack] = dummy; - continue; - } - - edm4hep::MCParticleData MCpart = mcParticles[MCindex]; - - // the MC-truth track parameters, in Delphes's comvention - TVectorD mcTrackParam = TrackParamFromMC_DelphesConv(MCpart); - - // the covariance matrix of the track, in Delphes's convention - TMatrixDSym Cov = VertexingUtils::get_trackCov(track); - - // if the covMat of the track is pathological (numerical precision issue, fraction of tracks = 5e-6): return original track - if ( Cov.Determinant() <= 0 ) { - result[itrack] = smeared_track; - continue; - } - - // scale the covariance matrix - for (int j = 0; j < 5; j++) - { - for (int k = 0; k < 5; k++) - { - Cov[j][k] = Cov[j][k] * (m_smear_parameters[j] * m_smear_parameters[k]); - } - } - // if (m_debug) { - // Cov.Print(); - //} - - // generate a new track state (in Delphes's convention) - TVectorD smeared_param = CovSmear(mcTrackParam, Cov, &m_random, m_debug); - - if ( smeared_param == zero ) { // Cholesky decomposition failed - result[itrack] = smeared_track; - continue; - } - - // back to the edm4hep conventions.. - - double scale0 = 1e-3; // convert mm to m - double scale1 = 1; - double scale2 = 0.5 * 1e3; // C = rho/2, convert from mm-1 to m-1 - double scale3 = 1e-3; // convert mm to m - double scale4 = 1.; - scale2 = -scale2; // sign of omega - - smeared_track.D0 = smeared_param[0] / scale0; - smeared_track.phi = smeared_param[1] / scale1; - smeared_track.omega = smeared_param[2] / scale2; - smeared_track.Z0 = smeared_param[3] / scale3; - smeared_track.tanLambda = smeared_param[4] / scale4; - - // transform rescaled cov matrix from Delphes convention to EDM4hep convention - smeared_track.covMatrix = track.covMatrix; - - smeared_track.covMatrix[0] = Cov[0][0] / (scale0 * scale0); - smeared_track.covMatrix[1] = Cov[1][0] / (scale1 * scale0); - smeared_track.covMatrix[2] = Cov[1][1] / (scale1 * scale1); - smeared_track.covMatrix[3] = Cov[2][0] / (scale2 * scale0); - smeared_track.covMatrix[4] = Cov[2][1] / (scale2 * scale1); - smeared_track.covMatrix[5] = Cov[2][2] / (scale2 * scale2); - smeared_track.covMatrix[6] = Cov[3][0] / (scale3 * scale0); - smeared_track.covMatrix[7] = Cov[3][1] / (scale3 * scale1); - smeared_track.covMatrix[8] = Cov[3][2] / (scale3 * scale2); - smeared_track.covMatrix[9] = Cov[3][3] / (scale3 * scale3); - smeared_track.covMatrix[10] = Cov[4][0] / (scale4 * scale0); - smeared_track.covMatrix[11] = Cov[4][1] / (scale4 * scale1); - smeared_track.covMatrix[12] = Cov[4][2] / (scale4 * scale2); - smeared_track.covMatrix[13] = Cov[4][3] / (scale4 * scale3); - smeared_track.covMatrix[14] = Cov[4][4] / (scale4 * scale4); - - if (m_debug) - { - std::cout << std::endl - << "Original track " << track.D0 << " " << track.phi << " " << track.omega << " " << track.Z0 << " " << track.tanLambda << std::endl; - std::cout << "Smeared track " << smeared_track.D0 << " " << smeared_track.phi << " " << smeared_track.omega << " " << smeared_track.Z0 << " " << smeared_track.tanLambda << std::endl; - std::cout << "MC particle " << mcTrackParam[0] / scale0 << " " << mcTrackParam[1] / scale1 << " " << mcTrackParam[2] / scale2 << " " << mcTrackParam[3] / scale3 << " " << mcTrackParam[4] / scale4 << std::endl; - for (int j = 0; j < 15; j++) - std::cout << "smeared cov matrix(" << j << "): " << smeared_track.covMatrix[j] << ", scale factor: " << smeared_track.covMatrix[j] / track.covMatrix[j] << std::endl; - } - result[itrack] = smeared_track; - - } // end loop on tracks - - return result; + edm4hep::MCParticleData MCpart = mcParticles[MCindex]; + + // the MC-truth track parameters, in Delphes's comvention + TVectorD mcTrackParam = TrackParamFromMC_DelphesConv(MCpart); + // and in edm4hep convention + TVectorD mcTrackParam_edm4hep = + VertexingUtils::Delphes2Edm4hep_TrackParam(mcTrackParam, false); + + // the covariance matrix of the track, in Delphes's convention + TMatrixDSym Cov = VertexingUtils::get_trackCov(track); + + // if the covMat of the track is pathological (numerical precision issue, + // fraction of tracks = 5e-6): return original track + if (Cov.Determinant() <= 0) { + result[itrack] = smeared_track; + continue; } - // ------------------------------------------------------------------------------------------- - - // to validate the SmearedTracks method.. : retrieve the TrackStates of the MC particles - - ROOT::VecOps::RVec mcTrackParameters(const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &alltracks, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles) - { - - int ntracks = alltracks.size(); - ROOT::VecOps::RVec result; - - edm4hep::TrackState dummy; - - for (int itrack = 0; itrack < ntracks; itrack++) - { - edm4hep::TrackState track = alltracks[itrack]; - - // find the corresponding MC particle - int MCindex = -1; - for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) - { - edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; - int track_index = rp.tracks_begin; - if (track_index == itrack) - { - MCindex = RP2MC_indices[ireco]; - break; - } - } // end loop on RPs - - if (MCindex < 0 || MCindex >= mcParticles.size()) - { - result.push_back(dummy); - continue; - } - - edm4hep::MCParticleData MCpart = mcParticles[MCindex]; - - // the MC-truth track parameters - edm4hep::TrackState mcTrack; - TVectorD mcTrackParam = TrackParamFromMC_DelphesConv(MCpart); - // put the into an edm4hep TrackState : - double scale0 = 1e-3; // convert mm to m - double scale2 = 0.5 * 1e3; // C = rho/2, convert from mm-1 to m-1 - double scale3 = 1e-3; // convert mm to m - scale2 = -scale2; // sign of omega - mcTrack.D0 = mcTrackParam[0] / scale0; - mcTrack.phi = mcTrackParam[1]; - mcTrack.omega = mcTrackParam[2] / scale2; - mcTrack.Z0 = mcTrackParam[3] / scale3; - mcTrack.tanLambda = mcTrackParam[4]; - - result.push_back(mcTrack); + // scale the covariance matrix + for (int j = 0; j < 5; j++) { + for (int k = 0; k < 5; k++) { + Cov[j][k] = Cov[j][k] * (m_smear_parameters[j] * m_smear_parameters[k]); } - return result; + } + // if (m_debug) { + // Cov.Print(); + //} + + // generate a new track state (in Delphes's convention) + TVectorD smeared_param_delphes = + CovSmear(mcTrackParam, Cov, &m_random, m_debug); + + if (smeared_param_delphes == zero) { // Cholesky decomposition failed + result[itrack] = smeared_track; + continue; } - // -------------------------------------------------------- + // back to the edm4hep conventions.. + TVectorD smeared_param = VertexingUtils::Delphes2Edm4hep_TrackParam( + smeared_param_delphes, false); + + smeared_track.D0 = smeared_param[0]; + smeared_track.phi = smeared_param[1]; + smeared_track.omega = smeared_param[2]; + smeared_track.Z0 = smeared_param[3]; + smeared_track.tanLambda = smeared_param[4]; + + // transform rescaled cov matrix from Delphes convention to EDM4hep + // convention + std::array covMatrix_edm4hep = + VertexingUtils::Delphes2Edm4hep_TrackCovMatrix(Cov, false); + + smeared_track.covMatrix = covMatrix_edm4hep; + + if (m_debug) { + std::cout << std::endl + << "Original track " << track.D0 << " " << track.phi << " " + << track.omega << " " << track.Z0 << " " << track.tanLambda + << std::endl; + std::cout << "Smeared track " << smeared_track.D0 << " " + << smeared_track.phi << " " << smeared_track.omega << " " + << smeared_track.Z0 << " " << smeared_track.tanLambda + << std::endl; + std::cout << "MC particle " << mcTrackParam_edm4hep[0] << " " + << mcTrackParam_edm4hep[1] << " " << mcTrackParam_edm4hep[2] + << " " << mcTrackParam_edm4hep[3] << " " + << mcTrackParam_edm4hep[4] << std::endl; + for (int j = 0; j < 15; j++) + std::cout << "smeared cov matrix(" << j + << "): " << smeared_track.covMatrix[j] << ", scale factor: " + << smeared_track.covMatrix[j] / track.covMatrix[j] + << std::endl; + } + result[itrack] = smeared_track; - // code from FB + } // end loop on tracks - TVectorD CovSmear(TVectorD x, TMatrixDSym C, TRandom *ran, bool debug = false) - { + return result; +} - // - // Check arrays - // - // Consistency of dimensions - Int_t Nvec = x.GetNrows(); - Int_t Nmat = C.GetNrows(); - if (Nvec != Nmat || Nvec == 0) - { - std::cout << "TrkUtil::CovSmear: vector/matrix mismatch. Aborting." << std::endl; - exit(EXIT_FAILURE); - } - // Positive diagonal elements - for (Int_t i = 0; i < Nvec; i++) - { - if (C(i, i) <= 0.0) - { - std::cout << "TrkUtil::CovSmear: covariance matrix has negative diagonal elements. Aborting." << std::endl; - exit(EXIT_FAILURE); - } +// ------------------------------------------------------------------------------------------- + +// to validate the SmearedTracks method.. : retrieve the TrackStates of the MC +// particles + +ROOT::VecOps::RVec mcTrackParameters( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &alltracks, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles) { + + int ntracks = alltracks.size(); + ROOT::VecOps::RVec result; + + edm4hep::TrackState dummy; + + for (int itrack = 0; itrack < ntracks; itrack++) { + edm4hep::TrackState track = alltracks[itrack]; + + // find the corresponding MC particle + int MCindex = -1; + for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) { + edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; + int track_index = rp.tracks_begin; + if (track_index == itrack) { + MCindex = RP2MC_indices[ireco]; + break; } - // - // Do a Choleski decomposition and random number extraction, with appropriate stabilization - // - TMatrixDSym CvN = C; - TMatrixDSym DCv(Nvec); - DCv.Zero(); - TMatrixDSym DCvInv(Nvec); - DCvInv.Zero(); - for (Int_t id = 0; id < Nvec; id++) - { - Double_t dVal = TMath::Sqrt(C(id, id)); - DCv(id, id) = dVal; - DCvInv(id, id) = 1.0 / dVal; + } // end loop on RPs + + if (MCindex < 0 || MCindex >= mcParticles.size()) { + result.push_back(dummy); + continue; + } + + edm4hep::MCParticleData MCpart = mcParticles[MCindex]; + + // the MC-truth track parameters + edm4hep::TrackState mcTrack; + TVectorD mcTrackParam_delphes = + TrackParamFromMC_DelphesConv(MCpart); // delphes convention, units = m + TVectorD mcTrackParam = VertexingUtils::Delphes2Edm4hep_TrackParam( + mcTrackParam_delphes, false); // edm4hep convention + + mcTrack.D0 = mcTrackParam[0]; + mcTrack.phi = mcTrackParam[1]; + mcTrack.omega = mcTrackParam[2]; + mcTrack.Z0 = mcTrackParam[3]; + mcTrack.tanLambda = mcTrackParam[4]; + + result.push_back(mcTrack); + } + return result; +} + +// -------------------------------------------------------- + +// code from FB + +TVectorD CovSmear(TVectorD x, TMatrixDSym C, TRandom *ran, bool debug = false) { + + // + // Check arrays + // + // Consistency of dimensions + Int_t Nvec = x.GetNrows(); + Int_t Nmat = C.GetNrows(); + if (Nvec != Nmat || Nvec == 0) { + std::cout << "TrkUtil::CovSmear: vector/matrix mismatch. Aborting." + << std::endl; + exit(EXIT_FAILURE); + } + // Positive diagonal elements + for (Int_t i = 0; i < Nvec; i++) { + if (C(i, i) <= 0.0) { + std::cout << "TrkUtil::CovSmear: covariance matrix has negative diagonal " + "elements. Aborting." + << std::endl; + exit(EXIT_FAILURE); + } + } + // + // Do a Choleski decomposition and random number extraction, with appropriate + // stabilization + // + TMatrixDSym CvN = C; + TMatrixDSym DCv(Nvec); + DCv.Zero(); + TMatrixDSym DCvInv(Nvec); + DCvInv.Zero(); + for (Int_t id = 0; id < Nvec; id++) { + Double_t dVal = TMath::Sqrt(C(id, id)); + DCv(id, id) = dVal; + DCvInv(id, id) = 1.0 / dVal; + } + CvN.Similarity(DCvInv); // Normalize diagonal to 1 + TDecompChol Chl(CvN); + Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix + if (!OK) { + std::cout << "SmearingObjects::CovSmear: covariance matrix is not positive " + "definite. Will use the original track." + << std::endl; + // exit(EXIT_FAILURE); + TVectorD zero(5); + for (int k = 0; k < 5; k++) { + zero(k) = 0.; + } + return zero; + } + TMatrixD U = Chl.GetU(); // Get Upper triangular matrix + TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular) + TVectorD r(Nvec); + for (Int_t i = 0; i < Nvec; i++) + r(i) = ran->Gaus(0.0, 1.0); // Array of normal random numbers + if (debug) + std::cout << " random nb " << ran->Gaus(0.0, 1.0) << std::endl; + TVectorD xOut = x + DCv * (Ut * r); // Observed parameter vector + // + return xOut; +} + +// ------------------------------------------------------------------------------------------- + +SmearedTracksdNdx::SmearedTracksdNdx(float scale, bool debug = false) { + + // rescale resolution by this factor + m_scale = scale; + + // debug flag + m_debug = debug; +} + +ROOT::VecOps::RVec SmearedTracksdNdx::operator()( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &dNdx, + const ROOT::VecOps::RVec &length, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles) { + + // returns a vector of dNdx that is parallel to the collection of full + // tracks (alltracks), i.e. same number of entries, same order. The method + // retrieve the MC particle that is associated to a track, and builds a "track + // state" out of the MC particle and regenerates a new value of the dNdx + + ROOT::VecOps::RVec result; + edm4hep::Quantity dummy; + + int ntracks = dNdx.size(); + result.resize(ntracks); + + // for dNdx calculation + TVector3 mc_mom; + TrkUtil tu; + + for (int itrack = 0; itrack < ntracks; itrack++) { + edm4hep::Quantity dndx = dNdx[itrack]; + edm4hep::Quantity smeared_dndx = dndx; + + // find the corresponding MC particle + int MCindex = -1; + for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) { + edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; + int track_index = rp.tracks_begin; + if (track_index == itrack) { + MCindex = RP2MC_indices[ireco]; + break; } - CvN.Similarity(DCvInv); // Normalize diagonal to 1 - TDecompChol Chl(CvN); - Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix - if (!OK) - { - std::cout << "SmearingObjects::CovSmear: covariance matrix is not positive definite. Will use the original track." << std::endl; - // exit(EXIT_FAILURE); - TVectorD zero(5); - for (int k=0; k <5; k++) { zero(k) = 0. ; } - return zero ; + } // end loop on RPs + + if (MCindex < 0 || + MCindex >= + mcParticles + .size()) { // in principle, this should not happen in delphes, + // each track should be matched to a MC particle. + result[itrack] = smeared_dndx; + continue; + } + + edm4hep::MCParticleData mc_part = mcParticles[MCindex]; + + // mom and mass evaluated on gen particles + mc_mom.SetXYZ(mc_part.momentum.x, mc_part.momentum.y, mc_part.momentum.z); + + float bg = mc_mom.Mag() / mc_part.mass; // beta * gamma + float muClu = + tu.Nclusters(bg, 0) * length[itrack]; // avg. number of clusters + + float Ncl = dndx.value * length[itrack]; + Ncl = std::max(muClu + m_scale * (Ncl - muClu), float(0.)); + + result[itrack].type = 0; + result[itrack].value = Ncl / length[itrack]; + + if (m_debug) { + std::cout << std::endl + << "requested smearing dNdx factor: " << m_scale << std::endl + << "gen part (PID, p): " << mc_part.PDG << " " << mc_mom.Mag() + << std::endl + << "original dNdx: " << dNdx[itrack].value << std::endl; + std::cout << "smeared dNdx : " << result[itrack].value << std::endl; + } + + } // end loop on tracks + + return result; +} + +// ------------------------------------------------------------------------------------------- + +SmearedTracksTOF::SmearedTracksTOF(float scale, bool debug = false) { + + // rescale resolution by this factor + m_scale = scale; + + // debug flag + m_debug = debug; +} + +ROOT::VecOps::RVec SmearedTracksTOF::operator()( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &trackdata, + const ROOT::VecOps::RVec &trackerhits, + const ROOT::VecOps::RVec &length, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles) { + // returns a vector of dNdx that is parallel to the collection of full + // tracks (alltracks), i.e. same number of entries, same order. The method + // retrieve the MC particle that is associated to a track, and builds a "track + // state" out of the MC particle and regenerates a new value of the dNdx + + ROOT::VecOps::RVec result; + edm4hep::TrackerHitData dummy; + + int ntracks = length.size(); + int nhits = trackerhits.size(); // 3x size of tracks since 3 hits per track + result.resize(nhits); + + TLorentzVector gen_p4; + + float c_light = 2.99792458e+8; + float mm_to_sec = 1e-03 / c_light; + + edm4hep::TrackerHitData thits_0, thits_1, thits_2; + edm4hep::TrackerHitData smeared_thits_0, smeared_thits_1, smeared_thits_2; + + for (int itrack = 0; itrack < ntracks; itrack++) { + + int idx_tin = trackdata.at(itrack).trackerHits_begin; // at IP + int idx_tpix = + trackdata.at(itrack).trackerHits_begin + 1; // at 1st pixel layer + int idx_tout = trackdata.at(itrack).trackerHits_end - 1; // at calo + + smeared_thits_0 = trackerhits.at(idx_tin); + smeared_thits_1 = trackerhits.at(idx_tpix); + smeared_thits_2 = trackerhits.at(idx_tout); + + // find the corresponding MC particle + int MCindex = -1; + for (int ireco = 0; ireco < allRecoParticles.size(); ireco++) { + edm4hep::ReconstructedParticleData rp = allRecoParticles[ireco]; + int track_index = rp.tracks_begin; + if (track_index == itrack) { + MCindex = RP2MC_indices[ireco]; + break; } - TMatrixD U = Chl.GetU(); // Get Upper triangular matrix - TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular) - TVectorD r(Nvec); - for (Int_t i = 0; i < Nvec; i++) - r(i) = ran->Gaus(0.0, 1.0); // Array of normal random numbers - if (debug) - std::cout << " random nb " << ran->Gaus(0.0, 1.0) << std::endl; - TVectorD xOut = x + DCv * (Ut * r); // Observed parameter vector - // - return xOut; + } // end loop on RPs + + if (MCindex < 0 || + MCindex >= + mcParticles + .size()) { // in principle, this should not happen in delphes, + // each track should be matched to a MC particle. + result[idx_tin] = smeared_thits_0; + result[idx_tpix] = smeared_thits_1; + result[idx_tout] = smeared_thits_2; + continue; + } + + edm4hep::MCParticleData mc_part = mcParticles[MCindex]; + + gen_p4.SetXYZM(mc_part.momentum.x, mc_part.momentum.y, mc_part.momentum.z, + mc_part.mass); + + // everything in second + float mc_tin = mc_part.time * mm_to_sec; + float mc_tof = length[itrack] / gen_p4.Beta() * mm_to_sec; + float mc_tout = mc_tin + mc_tof; + float reco_tout = trackerhits.at(idx_tout).time; + float smeared_tout = mc_tout + m_scale * (reco_tout - mc_tout); + + smeared_thits_2.time = smeared_tout; + + result[idx_tin] = smeared_thits_0; + result[idx_tpix] = smeared_thits_1; + result[idx_tout] = smeared_thits_2; + + if (m_debug) { + std::cout << std::endl + << "requested smearing tof factor: " << m_scale << std::endl + << "gen part (PID, p , beta, t_in, L): " << mc_part.PDG << " " + << gen_p4.P() << " " << gen_p4.Beta() << " " << mc_tin * 1e12 + << " " << length[itrack] << std::endl + << "gen t_out (ps): " << mc_tout * 1e12 << std::endl; + std::cout << "reco t_out (ps) : " << reco_tout * 1e12 << std::endl; + std::cout << "smeared t_out (ps) : " << smeared_tout * 1e12 << std::endl; } - // ------------------------------------------------------------------------------------------- + } // end loop on tracks + + return result; +} + +// ------------------------------------------------------------------------------------------- + +SmearedReconstructedParticle::SmearedReconstructedParticle(float scale, + int type, int mode, + bool debug = false) { + + // rescale resolution by this factor + m_scale = scale; - SmearedReconstructedParticle::SmearedReconstructedParticle(float scale, int type, int mode, bool debug = false) - { + // apply rescaling only to particle of given type + // supported: 11 (electrons), 13 (muons), 130 (neutral hadrons), 22 (photon), + // 0 (charged hadrons), -1 (all) + m_type = type; - // rescale resolution by this factor - m_scale = scale; + // 0: energy, 1: momentum + m_mode = mode; - // apply rescaling only to particle of given type - // supported: 11 (electrons), 13 (muons), 130 (neutral hadrons), 22 (photon), 0 (charged hadrons), -1 (all) - m_type = type; + // debug flag + m_debug = debug; +} - // 0: energy, 1: momentum - m_mode = mode; +ROOT::VecOps::RVec +SmearedReconstructedParticle::operator()( + const ROOT::VecOps::RVec + &allRecoParticles, + const ROOT::VecOps::RVec &RP2MC_indices, + const ROOT::VecOps::RVec &mcParticles) { - // debug flag - m_debug = debug; + // returns a vector of ReconstructedParticleData + // The method retrieve the MC particle that is associated to the + // ReconstructedParticle, and creates a new ReconstructedParticle with smeared + // parameters. + + int npart = allRecoParticles.size(); + + ROOT::VecOps::RVec result; + result.resize(npart); + + TLorentzVector gen_p4, reco_p4, smeared_p4; + + for (int ipart = 0; ipart < npart; ipart++) { + + edm4hep::ReconstructedParticleData reco_part = allRecoParticles[ipart]; + edm4hep::ReconstructedParticleData smeared_part = reco_part; + + int reco_part_type = abs(reco_part.type); + + // have to manually infer pid of ele/mu from mass because type not stored in + // reco particles + if (abs(reco_part.charge) > 0 and + abs(reco_part.mass - 0.000510999) < 1.e-05) { + reco_part_type = 11; + } else if (abs(reco_part.charge) > 0 and + abs(reco_part.mass - 0.105658) < 1.e-03) { + reco_part_type = 13; } - ROOT::VecOps::RVec SmearedReconstructedParticle::operator()( - const ROOT::VecOps::RVec &allRecoParticles, - const ROOT::VecOps::RVec &RP2MC_indices, - const ROOT::VecOps::RVec &mcParticles) - { - - // returns a vector of ReconstructedParticleData - // The method retrieve the MC particle that is associated to the ReconstructedParticle, - // and creates a new ReconstructedParticle with smeared parameters. - - int npart = allRecoParticles.size(); - - ROOT::VecOps::RVec result; - result.resize(npart); - - TLorentzVector gen_p4, reco_p4, smeared_p4; - - for (int ipart = 0; ipart < npart; ipart++) - { - - edm4hep::ReconstructedParticleData reco_part = allRecoParticles[ipart]; - edm4hep::ReconstructedParticleData smeared_part = reco_part; - - int reco_part_type = abs(reco_part.type); - - // have to manually infer pid of ele/mu from mass because type not stored in reco particles - if (abs(reco_part.charge) > 0 and abs(reco_part.mass - 0.000510999) < 1.e-05) - { - reco_part_type = 11; - } - else if (abs(reco_part.charge) > 0 and abs(reco_part.mass - 0.105658) < 1.e-03) - { - reco_part_type = 13; - } - - // find the corresponding MC particle - int MCindex = -1; - MCindex = RP2MC_indices[ipart]; - - // smear particle only if MC particle found, else return original particle - // and if type == requested - if (MCindex >= 0 and MCindex < mcParticles.size() and reco_part_type == m_type) - { - edm4hep::MCParticleData mc_part = mcParticles[MCindex]; - - gen_p4.SetXYZM(mc_part.momentum.x, mc_part.momentum.y, mc_part.momentum.z, mc_part.mass); - reco_p4.SetXYZM(reco_part.momentum.x, reco_part.momentum.y, reco_part.momentum.z, reco_part.mass); - - float smeared_p = -1; - - if (m_mode == 0) - { - // rescale existing smearing of the energy - smeared_part.energy = gen_p4.E() + m_scale * (reco_p4.E() - gen_p4.E()); - - // recompute momentum magnitude - smeared_p = (smeared_part.energy - reco_p4.M()) > 0 ? std::sqrt(smeared_part.energy * smeared_part.energy - reco_p4.M() * reco_p4.M()) : smeared_part.energy; - - // recompute mom x, y, z using original reco particle direction - smeared_part.momentum.x = smeared_p * std::sin(reco_p4.Theta()) * std::cos(reco_p4.Phi()); - smeared_part.momentum.y = smeared_p * std::sin(reco_p4.Theta()) * std::sin(reco_p4.Phi()); - smeared_part.momentum.z = smeared_p * std::cos(reco_p4.Theta()); - - smeared_p4.SetXYZM(smeared_part.momentum.x, smeared_part.momentum.y, smeared_part.momentum.z, smeared_part.mass); - } - - // momentum resolution mode - else if (m_mode == 1) - { - // rescale existing momentum smearing - smeared_p = gen_p4.P() + m_scale * (reco_p4.P() - gen_p4.P()); - - // recompute energy - smeared_part.energy = std::sqrt(smeared_p * smeared_p + reco_p4.M() * reco_p4.M()); - // recompute mom x, y, z using original reco particle direction - smeared_part.momentum.x = smeared_p * std::sin(reco_p4.Theta()) * std::cos(reco_p4.Phi()); - smeared_part.momentum.y = smeared_p * std::sin(reco_p4.Theta()) * std::sin(reco_p4.Phi()); - smeared_part.momentum.z = smeared_p * std::cos(reco_p4.Theta()); - - } - - // return mc truth particle - else if (m_mode == -1) - { - smeared_part.energy = gen_p4.E(); - smeared_p = gen_p4.P(); - smeared_p4 = gen_p4; - - // recompute mom x, y, z using original reco particle direction - smeared_part.momentum.x = gen_p4.Px(); - smeared_part.momentum.y = gen_p4.Py(); - smeared_part.momentum.z = gen_p4.Pz(); - - // set type - smeared_part.type = mc_part.PDG; - } - - if (m_debug) - { - std::cout << std::endl - << "requested smearing energy factor: " << m_scale << std::endl - << "gen part (PID, p, theta, phi, m): " << mc_part.PDG << " " << gen_p4.P() << " " << gen_p4.Theta() << " " << gen_p4.Phi() << " " << gen_p4.M() << std::endl - << "reco part (PID, p, theta, phi, m): " << reco_part_type << " " << reco_p4.P() << " " << reco_p4.Theta() << " " << reco_p4.Phi() << " " << reco_p4.M() << std::endl; - std::cout << "smeared part (PID, p, theta, phi, m): " << smeared_part.type << " " << smeared_p4.P() << " " << smeared_p4.Theta() << " " << smeared_p4.Phi() << " " << smeared_p4.M() << std::endl; - } - } - - result[ipart] = smeared_part; - - } // end loop on particles - - return result; + // find the corresponding MC particle + int MCindex = -1; + MCindex = RP2MC_indices[ipart]; + + // smear particle only if MC particle found, else return original particle + // and if type == requested + if (MCindex >= 0 and MCindex < mcParticles.size() and + reco_part_type == m_type) { + edm4hep::MCParticleData mc_part = mcParticles[MCindex]; + + gen_p4.SetXYZM(mc_part.momentum.x, mc_part.momentum.y, mc_part.momentum.z, + mc_part.mass); + reco_p4.SetXYZM(reco_part.momentum.x, reco_part.momentum.y, + reco_part.momentum.z, reco_part.mass); + + float smeared_p = -1; + + if (m_mode == 0) { + // rescale existing smearing of the energy + smeared_part.energy = std::max( + gen_p4.E() + m_scale * (reco_p4.E() - gen_p4.E()), reco_p4.M()); + + // recompute momentum magnitude + smeared_p = std::sqrt(smeared_part.energy * smeared_part.energy - + reco_p4.M() * reco_p4.M()); + + // recompute mom x, y, z using original reco particle direction + smeared_part.momentum.x = + smeared_p * std::sin(reco_p4.Theta()) * std::cos(reco_p4.Phi()); + smeared_part.momentum.y = + smeared_p * std::sin(reco_p4.Theta()) * std::sin(reco_p4.Phi()); + smeared_part.momentum.z = smeared_p * std::cos(reco_p4.Theta()); + + smeared_p4.SetXYZM(smeared_part.momentum.x, smeared_part.momentum.y, + smeared_part.momentum.z, smeared_part.mass); + } + + // momentum resolution mode + else if (m_mode == 1) { + // rescale existing momentum smearing + smeared_p = + std::max(float(gen_p4.P() + m_scale * (reco_p4.P() - gen_p4.P())), + float(0.)); + + // recompute energy + smeared_part.energy = + std::sqrt(smeared_p * smeared_p + reco_p4.M() * reco_p4.M()); + // recompute mom x, y, z using original reco particle direction + smeared_part.momentum.x = + smeared_p * std::sin(reco_p4.Theta()) * std::cos(reco_p4.Phi()); + smeared_part.momentum.y = + smeared_p * std::sin(reco_p4.Theta()) * std::sin(reco_p4.Phi()); + smeared_part.momentum.z = smeared_p * std::cos(reco_p4.Theta()); + + } + + // return mc truth particle + else if (m_mode == -1) { + smeared_part.energy = gen_p4.E(); + smeared_p = gen_p4.P(); + smeared_p4 = gen_p4; + + // recompute mom x, y, z using original reco particle direction + smeared_part.momentum.x = gen_p4.Px(); + smeared_part.momentum.y = gen_p4.Py(); + smeared_part.momentum.z = gen_p4.Pz(); + + // set type + smeared_part.type = mc_part.PDG; + } + + if (m_debug) { + std::cout << std::endl + << "requested smearing energy factor: " << m_scale + << std::endl + << "gen part (PID, p, theta, phi, m): " << mc_part.PDG + << " " << gen_p4.P() << " " << gen_p4.Theta() << " " + << gen_p4.Phi() << " " << gen_p4.M() << std::endl + << "reco part (PID, p, theta, phi, m): " << reco_part_type + << " " << reco_p4.P() << " " << reco_p4.Theta() << " " + << reco_p4.Phi() << " " << reco_p4.M() << std::endl; + std::cout << "smeared part (PID, p, theta, phi, m): " + << smeared_part.type << " " << smeared_p4.P() << " " + << smeared_p4.Theta() << " " << smeared_p4.Phi() << " " + << smeared_p4.M() << std::endl; + } } - } // end NS -} // end NS FCCAnalyses + result[ipart] = smeared_part; + + } // end loop on particles + + return result; +} + +} // namespace SmearObjects +} // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/VertexFinderActs.cc b/analyzers/dataframe/src/VertexFinderActs.cc index ad8eac8726..75853a190c 100644 --- a/analyzers/dataframe/src/VertexFinderActs.cc +++ b/analyzers/dataframe/src/VertexFinderActs.cc @@ -68,7 +68,8 @@ VertexFinderAMVF(ROOT::VecOps::RVec tracks ){ // Set up deterministic annealing with user-defined temperatures std::vector temperatures{8.0, 4.0, 2.0, 1.4142136, 1.2247449, 1.0}; - Acts::AnnealingUtility::Config annealingConfig(temperatures); + Acts::AnnealingUtility::Config annealingConfig; + annealingConfig.setOfTemperatures = temperatures; Acts::AnnealingUtility annealingUtility(annealingConfig); @@ -101,7 +102,8 @@ VertexFinderAMVF(ROOT::VecOps::RVec tracks ){ using Finder = Acts::AdaptiveMultiVertexFinder; //using Finder = Acts::AdaptiveMultiVertexFinder; //Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator, linearizer); - Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator, linearizer, bField); + Finder::Config finderConfig = {std::move(fitter), seedFinder, ipEstimator, + std::move(linearizer), bField}; // We do not want to use a beamspot constraint here finderConfig.useBeamSpotConstraint = false; diff --git a/analyzers/dataframe/src/VertexFitterSimple.cc b/analyzers/dataframe/src/VertexFitterSimple.cc index e09b81225a..6d6ecd4edf 100644 --- a/analyzers/dataframe/src/VertexFitterSimple.cc +++ b/analyzers/dataframe/src/VertexFitterSimple.cc @@ -1,856 +1,364 @@ -#include "FCCAnalyses/VertexFitterSimple.h" -#include "FCCAnalyses/MCParticle.h" - -#include - -#include "TFile.h" -#include "TString.h" - -namespace FCCAnalyses{ - -namespace VertexFitterSimple{ - -TVector3 ParToP(TVectorD Par){ - double fB = 2; // 2 Tesla - - Double_t C = Par(2); - Double_t phi0 = Par(1); - Double_t ct = Par(4); - // - TVector3 Pval; - Double_t pt = fB*0.2998 / TMath::Abs(2 * C); - Pval(0) = pt*TMath::Cos(phi0); - Pval(1) = pt*TMath::Sin(phi0); - Pval(2) = pt*ct; - // - return Pval; -} - - -TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q){ - - double fB = 2; // 2 Tesla - - // - TVectorD Par(5); - // Transverse parameters - Double_t a = -Q*fB*0.2998; // Units are Tesla, GeV and meters - Double_t pt = p.Pt(); - Double_t C = a / (2 * pt); // Half curvature - //cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << endl; - Double_t r2 = x.Perp2(); - Double_t cross = x(0)*p(1) - x(1)*p(0); - Double_t T = TMath::Sqrt(pt*pt - 2 * a*cross + a*a*r2); - Double_t phi0 = TMath::ATan2((p(1) - a*x(0)) / T, (p(0) + a*x(1)) / T); // Phi0 - Double_t D; // Impact parameter D - if (pt < 10.0) D = (T - pt) / a; - else D = (-2 * cross + a*r2) / (T + pt); - // - Par(0) = D; // Store D - Par(1) = phi0; // Store phi0 - Par(2) = C; // Store C - //Longitudinal parameters - Double_t B = C*TMath::Sqrt(TMath::Max(r2 - D*D,0.0) / (1 + 2 * C*D)); - Double_t st = TMath::ASin(B) / C; - Double_t ct = p(2) / pt; - Double_t z0 = x(2) - ct*st; - // - Par(3) = z0; // Store z0 - Par(4) = ct; // Store cot(theta) - // - return Par; -} - - -// -//TH1F* hTry; -// -Double_t FastRv(TVectorD p1, TVectorD p2){ - // - // Find radius of intersection between two tracks in the transverse plane - // - // p = (D,phi, C) - // - // Solving matrix - TMatrixDSym H(2); - H(0, 0) = -TMath::Cos(p2(1)); - H(0, 1) = TMath::Cos(p1(1)); - H(1, 0) = -TMath::Sin(p2(1)); - H(1, 1) = TMath::Sin(p1(1)); - Double_t Det = TMath::Sin(p2(1) - p1(1)); - H *= 1.0 / Det; - // - // Convergence parameters - Int_t Ntry = 0; - Int_t NtryMax = 100; - Double_t eps = 1000.; - Double_t epsMin = 1.0e-6; - // - // Vertex finding loop - // - TVectorD cterm(2); - cterm(0) = p1(0); - cterm(1) = p2(0); - TVectorD xv(2); - Double_t R = 1000.; - while (eps > epsMin) - { - xv = H * cterm; - Ntry++; - if (Ntry > NtryMax) - { - std::cout << "FastRv: maximum number of iteration reached" << std::endl; - break; - } - Double_t Rnew = TMath::Sqrt(xv(0) * xv(0) + xv(1) * xv(1)); - eps = Rnew - R; - R = Rnew; - cterm(0) = p1(2) * R * R; - cterm(1) = p2(2) * R * R; - } - // - return R; -} -TMatrixDSym RegInv3(TMatrixDSym &Smat0){ - // - // Regularized inversion of symmetric 3x3 matrix with positive diagonal elements - // - TMatrixDSym Smat = Smat0; - Int_t N = Smat.GetNrows(); - if (N != 3){ - std::cout << "RegInv3 called with matrix size != 3. Abort & return standard inversion." << std::endl; - return Smat.Invert(); - } - TMatrixDSym D(N); D.Zero(); - Bool_t dZero = kTRUE; // No elements less or equal 0 on the diagonal - for (Int_t i = 0; i < N; i++) if (Smat(i, i) <= 0.0)dZero = kFALSE; - if (dZero){ - for (Int_t i = 0; i < N; i++) D(i, i) = 1.0 / TMath::Sqrt(Smat(i, i)); - TMatrixDSym RegMat = Smat.Similarity(D); - TMatrixDSym Q(2); - for (Int_t i = 0; i < 2; i++){ - for (Int_t j = 0; j < 2; j++)Q(i, j) = RegMat(i, j); - } - Double_t Det = 1 - Q(0, 1)*Q(1, 0); - TMatrixDSym H(2); - H = Q; - H(0, 1) = -Q(0, 1); - H(1, 0) = -Q(1, 0); - TVectorD p(2); - p(0) = RegMat(0, 2); - p(1) = RegMat(1, 2); - Double_t pHp = H.Similarity(p); - Double_t h = pHp-Det; - // - TMatrixDSym pp(2); pp.Rank1Update(p); - TMatrixDSym F = (h*H) - pp.Similarity(H); - F *= 1.0 / Det; - TVectorD b = H*p; - TMatrixDSym InvReg(3); - for (Int_t i = 0; i < 2; i++) - { - InvReg(i, 2) = b(i); - InvReg(2, i) = b(i); - for (Int_t j = 0; j < 2; j++) InvReg(i, j) = F(i, j); - } - InvReg(2, 2) = -Det; - // - InvReg *= 1.0 / h; - // - // - return InvReg.Similarity(D); - } - else - { - //std::cout << "RegInv3: found negative elements in diagonal. Return standard inversion." << std::endl; - return Smat.Invert(); - } -} -// -// -// -TMatrixD Fill_A(TVectorD par, Double_t phi){ - // - // Derivative of track 3D position vector with respect to track parameters at constant phase - // - // par = vector of track parameters - // phi = phase - // - TMatrixD A(3, 5); - // - // Decode input arrays - // - Double_t D = par(0); - Double_t p0 = par(1); - Double_t C = par(2); - Double_t z0 = par(3); - Double_t ct = par(4); - // - // Fill derivative matrix dx/d alpha - // D - A(0, 0) = -TMath::Sin(p0); - A(1, 0) = TMath::Cos(p0); - A(2, 0) = 0.0; - // phi0 - A(0, 1) = -D*TMath::Cos(p0) + (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C); - A(1, 1) = -D*TMath::Sin(p0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); - A(2, 1) = 0.0; - // C - A(0, 2) = -(TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C*C); - A(1, 2) = (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C*C); - A(2, 2) = -ct*phi / (2 * C*C); - // z0 - A(0, 3) = 0.0; - A(1, 3) = 0.0; - A(2, 3) = 1.0; - // ct = lambda - A(0, 4) = 0.0; - A(1, 4) = 0.0; - A(2, 4) = phi / (2 * C); - // - return A; -} - -// -TVectorD Fill_a(TVectorD par, Double_t phi){ - // - // Derivative of track 3D position vector with respect to phase at constant track parameters - // - // par = vector of track parameters - // phi = phase - // - TVectorD a(3); - // - // Decode input arrays - // - Double_t D = par(0); - Double_t p0 = par(1); - Double_t C = par(2); - Double_t z0 = par(3); - Double_t ct = par(4); - // - a(0) = TMath::Cos(phi + p0) / (2 * C); - a(1) = TMath::Sin(phi + p0) / (2 * C); - a(2) = ct / (2 * C); - // - return a; -} -// - -TVectorD Fill_x0(TVectorD par){ - // - // Calculate track 3D position at R = |D| (minimum approach to z-axis) - // - TVectorD x0(3); - // - // Decode input arrays - // - Double_t D = par(0); - Double_t p0 = par(1); - Double_t C = par(2); - Double_t z0 = par(3); - Double_t ct = par(4); - // - x0(0) = -D *TMath::Sin(p0); - x0(1) = D*TMath::Cos(p0); - x0(2) = z0; - // - return x0; -} - -// -TVectorD Fill_x(TVectorD par, Double_t phi){ - // - // Calculate track 3D position for a given phase, phi - // - TVectorD x(3); - // - // Decode input arrays - // - Double_t D = par(0); - Double_t p0 = par(1); - Double_t C = par(2); - Double_t z0 = par(3); - Double_t ct = par(4); - // - TVectorD x0 = Fill_x0(par); - x(0) = x0(0) + (TMath::Sin(phi + p0) - TMath::Sin(p0)) / (2 * C); - x(1) = x0(1) - (TMath::Cos(phi + p0) - TMath::Cos(p0)) / (2 * C); - x(2) = x0(2) + ct*phi / (2 * C); - // - return x; -} - - - -VertexingUtils::FCCAnalysesVertex VertexFitter( int Primary, - ROOT::VecOps::RVec recoparticles, - ROOT::VecOps::RVec thetracks, - bool BeamSpotConstraint, - double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, - double bsc_x, double bsc_y, double bsc_z ) { - - - - // input = a collection of recoparticles (in case one want to make associations to RecoParticles ?) - // and thetracks = the collection of all TrackState in the event - - VertexingUtils::FCCAnalysesVertex thevertex; - - // retrieve the tracks associated to the recoparticles - ROOT::VecOps::RVec tracks = ReconstructedParticle2Track::getRP2TRK( recoparticles, thetracks ); - - // and run the vertex fitter - - //FCCAnalysesVertex thevertex = VertexFitter_Tk( Primary, tracks, thetracks) ; - thevertex = VertexFitter_Tk( Primary, tracks, - thetracks, - BeamSpotConstraint, bsc_sigmax, bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z ); - - //fill the indices of the tracks - ROOT::VecOps::RVec reco_ind; - int Ntr = tracks.size(); - for (auto & p: recoparticles) { - //std::cout << " in VertexFitter: a recoparticle with charge = " << p.charge << std::endl; - if ( p.tracks_begin >=0 && p.tracks_begin tracks, - bool BeamSpotConstraint, - double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, - double bsc_x, double bsc_y, double bsc_z ) { - - ROOT::VecOps::RVec dummy; - return VertexFitter_Tk( Primary, tracks, dummy, BeamSpotConstraint, bsc_sigmax, bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z ); -} - - -VertexingUtils::FCCAnalysesVertex VertexFitter_Tk( int Primary, - ROOT::VecOps::RVec tracks, - const ROOT::VecOps::RVec& alltracks, - bool BeamSpotConstraint, - double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, - double bsc_x, double bsc_y, double bsc_z ) { - - // Units for the beam-spot : mum - // See https://github.com/HEP-FCC/FCCeePhysicsPerformance/tree/master/General#generating-events-under-realistic-fcc-ee-environment-conditions - - // final results : - VertexingUtils::FCCAnalysesVertex TheVertex; - - edm4hep::VertexData result; - ROOT::VecOps::RVec reco_chi2; - ROOT::VecOps::RVec< TVectorD > updated_track_parameters; - ROOT::VecOps::RVec reco_ind; - ROOT::VecOps::RVec final_track_phases; - ROOT::VecOps::RVec< TVector3 > updated_track_momentum_at_vertex; - - // if the collection of all tracks has been passed, keep trace of the indices of the tracks that are used to fit this vertex - if ( alltracks.size() > 0 ) { - for (int i=0; i < tracks.size(); i++) { // the fitted tracks - edm4hep::TrackState tr1 = tracks[i]; - for ( int j=0; j < alltracks.size(); j++) { // the collection of all tracks - edm4hep::TrackState tr2 = alltracks[j]; - if ( VertexingUtils::compare_Tracks( tr1, tr2 ) ) { - reco_ind.push_back( j ) ; - break; - } - } - } - } - - TheVertex.vertex = result; - TheVertex.reco_chi2 = reco_chi2; - TheVertex.updated_track_parameters = updated_track_parameters; - TheVertex.reco_ind = reco_ind; - TheVertex.final_track_phases = final_track_phases; - TheVertex.updated_track_momentum_at_vertex = updated_track_momentum_at_vertex; - - - int Ntr = tracks.size(); - TheVertex.ntracks = Ntr; - if ( Ntr <= 1) return TheVertex; // can not reconstruct a vertex with only one track... - - - bool debug = false; - if (debug) std::cout << " enter in VertexFitter_Tk for the Bs decay vertex " << std::endl; - - // if a beam-spot constraint is required : - TMatrixDSym BeamSpotCovI(3); - TVectorD BeamSpotPos(3); - if (BeamSpotConstraint) { // fill in the inverse of the covariance matrix. Convert the units into meters - BeamSpotCovI(0,0) = 1./pow( bsc_sigmax * 1e-6, 2) ; // mum to m - BeamSpotCovI(1,1) = 1./pow( bsc_sigmay * 1e-6, 2) ; - BeamSpotCovI(2,2) = 1./pow( bsc_sigmaz * 1e-6, 2) ; - BeamSpotPos(0) = bsc_x * 1e-6; - BeamSpotPos(1) = bsc_y * 1e-6 ; - BeamSpotPos(2) = bsc_z * 1e-6 ; - } - - Double_t *final_chi2 = new Double_t[Ntr]; - Double_t *final_phases = new Double_t[Ntr]; - std::vector< TVectorD > final_delta_alpha ; - TVectorD dummy(5); - for (int i=0; i < Ntr; i++) { - final_delta_alpha.push_back( dummy ); - } - - - // - // Vertex fit (units are meters) - // - // Initial variable definitions - TVectorD x0(3); for (Int_t v = 0; v < 3; v++)x0(v) = 100.; // set to large value - Double_t Chi2 = 0; - // - - TVectorD x(3); - TMatrixDSym covX(3); - - - // Stored quantities - Double_t *fi = new Double_t[Ntr]; // Phases - TVectorD **x0i = new TVectorD*[Ntr]; // Track expansion point - TVectorD **ai = new TVectorD*[Ntr]; // dx/dphi - Double_t *a2i = new Double_t[Ntr]; // a'Wa - TMatrixDSym **Di = new TMatrixDSym*[Ntr]; // W-WBW - TMatrixDSym **Wi = new TMatrixDSym*[Ntr]; // (ACA')^-1 - TMatrixDSym **Winvi = new TMatrixDSym*[Ntr]; // ACA' - TMatrixD **Ai = new TMatrixD*[Ntr]; // A - TMatrixDSym **Covi = new TMatrixDSym*[Ntr]; // Cov matrix of the track parameters - - // - // vertex radius approximation - // Maximum impact parameter - Double_t Rd = 0; - for (Int_t i = 0; i < Ntr; i++) - { - //ObsTrk* t = tracks[i]; - //TVectorD par = t->GetObsPar(); - edm4hep::TrackState t = tracks[i] ; - TVectorD par = VertexingUtils::get_trackParam( t ) ; - Double_t Dabs = TMath::Abs(par(0)); - if (Dabs > Rd)Rd = Dabs; - } - // - // Find track pair with largest phi difference - Int_t isel; Int_t jsel; // selected track indices - Double_t dphiMax = -9999.; // Max phi difference - for (Int_t i = 0; i < Ntr-1; i++) - { - //ObsTrk* ti = tracks[i]; - //TVectorD pari = ti->GetObsPar(); - edm4hep::TrackState ti = tracks[i] ; - TVectorD pari = VertexingUtils::get_trackParam( ti ); - Double_t phi1 = pari(1); - - for (Int_t j = i+1; j < Ntr; j++) - { - //ObsTrk* tj = tracks[j]; - //TVectorD parj = tj->GetObsPar(); - edm4hep::TrackState tj = tracks[j]; - TVectorD parj = VertexingUtils::get_trackParam( tj ); - Double_t phi2 = parj(1); - Double_t dphi = TMath::Abs(phi2 - phi1); - if (dphi > TMath::Pi())dphi = TMath::TwoPi() - dphi; - if (dphi > dphiMax) - { - isel = i; jsel = j; - dphiMax = dphi; - } - } - } - // - // - //ObsTrk* t1 = tracks[isel]; - //TVectorD p1 = t1->GetObsPar(); - edm4hep::TrackState t1 = tracks[isel]; - TVectorD p1 = VertexingUtils::get_trackParam( t1 ); - //ObsTrk* t2 = tracks[jsel]; - //TVectorD p2 = t2->GetObsPar(); - edm4hep::TrackState t2 = tracks[jsel]; - TVectorD p2 = VertexingUtils::get_trackParam( t2 ); - Double_t R = FastRv(p1, p2); - if (R > 1.0) R = Rd; - R = 0.5 * (R + Rd); - // - // Iteration properties - // - Int_t Ntry = 0; - Int_t TryMax = 100; - if (BeamSpotConstraint) TryMax = TryMax * 5; - Double_t eps = 1.0e-9; // vertex stability - Double_t epsi = 1000.; - // - while (epsi > eps && Ntry < TryMax) // Iterate until found vertex is stable - { - x.Zero(); - TVectorD cterm(3); TMatrixDSym H(3); TMatrixDSym DW1D(3); - covX.Zero(); // Reset vertex covariance - cterm.Zero(); // Reset constant term - H.Zero(); // Reset H matrix - DW1D.Zero(); - // - for (Int_t i = 0; i < Ntr; i++) - { - // Get track helix parameters and their covariance matrix - //ObsTrk *t = tracks[i]; - //TVectorD par = t->GetObsPar(); - //TMatrixDSym Cov = t->GetCov(); - edm4hep::TrackState t = tracks[i] ; - TVectorD par = VertexingUtils::get_trackParam( t ) ; - TMatrixDSym Cov = VertexingUtils::get_trackCov( t ); - Covi[i] = new TMatrixDSym(Cov); // Store matrix - Double_t fs; - if (Ntry <= 0) // Initialize all phases on first pass - { - Double_t D = par(0); - Double_t C = par(2); - Double_t arg = TMath::Max(1.0e-6, (R*R - D*D) / (1 + 2 * C*D)); - fs = 2 * TMath::ASin(C*TMath::Sqrt(arg)); - fi[i] = fs; - } - // - // Starting values - // - fs = fi[i]; // Get phase - TVectorD xs = Fill_x(par, fs); - x0i[i] = new TVectorD(xs); // Start helix position - // W matrix = (A*C*A')^-1; W^-1 = A*C*A' - TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters - Ai[i] = new TMatrixD(A); // Store matrix - TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' - Winvi[i] = new TMatrixDSym(Winv); // Store W^-1 matrix - TMatrixDSym W = RegInv3(Winv); // W = (A*C*A')^-1 - Wi[i] = new TMatrixDSym(W); // Store W matrix - TVectorD a = Fill_a(par, fs); // a = dx/ds = derivatives wrt phase - ai[i] = new TVectorD(a); // Store a - Double_t a2 = W.Similarity(a); - a2i[i] = a2; // Store a2 - // Build D matrix - TMatrixDSym B(3); - B.Rank1Update(a, 1.0); - B *= -1. / a2; - B.Similarity(W); - TMatrixDSym Ds = W+B; // D matrix - Di[i] = new TMatrixDSym(Ds); // Store D matrix - TMatrixDSym DsW1Ds = Winv.Similarity(Ds); // Service matrix to calculate covX - DW1D += DsW1Ds; - // Update hessian - H += Ds; - // update constant term - cterm += Ds * xs; - } // End loop on tracks - // - - TMatrixDSym H0 = H; - - if (BeamSpotConstraint) { - H += BeamSpotCovI ; - cterm += BeamSpotCovI * BeamSpotPos ; - DW1D += BeamSpotCovI ; - } - - // update vertex position - TMatrixDSym H1 = RegInv3(H); - x = H1*cterm; - - // Update vertex covariance - covX = DW1D.Similarity(H1); - - // Update phases and chi^2 - Chi2 = 0.0; - for (Int_t i = 0; i < Ntr; i++) - { - TVectorD lambda = (*Di[i])*(*x0i[i] - x); - TMatrixDSym Wm1 = *Winvi[i]; - Double_t addChi2 = Wm1.Similarity(lambda);; - //Chi2 += Wm1.Similarity(lambda); - Chi2 += addChi2; - final_chi2[i] = addChi2; - TVectorD a = *ai[i]; - TVectorD b = (*Wi[i])*(x - *x0i[i]); - for (Int_t j = 0; j < 3; j++)fi[i] += a(j)*b(j) / a2i[i]; - final_phases[i] = fi[i]; - - TMatrixD ta(TMatrixD::kTransposed, *Ai[i]); - TMatrixDSym kk(5); - kk = *Covi[i]; - final_delta_alpha[i] = kk * ta * lambda; // that's minus delta_alpha - } - // - - TVectorD dx = x - x0; - x0 = x; - // update vertex stability - TMatrixDSym Hess = RegInv3(covX); - - epsi = Hess.Similarity(dx); - Ntry++; - //if ( Ntry >= TryMax) std::cout << " ... in VertexFitterSimple, Ntry >= TryMax " << std::endl; - - if (BeamSpotConstraint) { - - // add the following term to the chi2 : - TVectorD dx_beamspot = x - BeamSpotPos ; - Double_t chi2_bsc = BeamSpotCovI.Similarity( dx_beamspot ); - //Chi2 += chi2_bsc -3; - Chi2 += chi2_bsc ; - - } - - - - // - // Cleanup - // - for (Int_t i = 0; i < Ntr; i++) - { - x0i[i]->Clear(); - Winvi[i]->Clear(); - Wi[i]->Clear(); - ai[i]->Clear(); - Di[i]->Clear(); - Ai[i]->Clear(); - Covi[i]->Clear(); - - delete x0i[i]; - delete Winvi[i]; - delete Wi[i]; - delete ai[i]; - delete Di[i]; - delete Ai[i]; - delete Covi[i]; - } - } - // - delete[] fi; // Phases - delete[] x0i; // Track expansion point - delete[] ai; // dx/dphi - delete[] a2i; // a'Wa - delete[] Di; // W-WBW - delete[] Wi; // (ACA')^-1 - delete[] Winvi; // ACA' - delete[] Ai ; // A - delete[] Covi; // Cov - - // - //return Chi2; - - // store the results in an edm4hep::VertexData object - // go back from meters to millimeters for the units - float conv = 1e3; - std::array covMatrix; // covMat in edm4hep is a LOWER-triangle matrix. - covMatrix[0] = covX(0,0) * pow(conv,2); - covMatrix[1] = covX(1,0) * pow(conv,2); - covMatrix[2] = covX(1,1) * pow(conv,2); - covMatrix[3] = covX(2,0) * pow(conv,2); - covMatrix[4] = covX(2,1) * pow(conv,2); - covMatrix[5] = covX(2,2) * pow(conv,2); - - float Ndof = 2.0 * Ntr - 3.0; ; - - result.primary = Primary; - result.chi2 = Chi2 /Ndof ; // I store the normalised chi2 here - result.position = edm4hep::Vector3f( x(0)*conv, x(1)*conv, x(2)*conv ) ; // store the vertex in mm - result.covMatrix = covMatrix; - result.algorithmType = 1; - - // Need to fill the associations ... - - double scale0 = 1e-3; //convert mm to m - double scale1 = 1; - double scale2 = 0.5*1e3; // C = rho/2, convert from mm-1 to m-1 - double scale3 = 1e-3 ; //convert mm to m - double scale4 = 1.; - - scale2 = -scale2 ; // sign of omega (sign convention) - - for (Int_t i = 0; i < Ntr; i++) { - - edm4hep::TrackState t = tracks[i] ; - TVectorD par = VertexingUtils::get_trackParam( t ) ; - - // initial momentum : - //TVector3 ptrack_ini = ParToP( par ); - //std::cout << "----- Track # " << i << " initial track momentum : " << std::endl; - //ptrack_ini.Print(); - - // uncomment below to get the post-fit track parameters : - par -= final_delta_alpha[i] ; - - //std::cout << " Track i = " << i << " --- delta_alpha : " << std::endl; - //final_delta_alpha[i].Print(); - - // ( px, py, pz) of the track - TVector3 ptrack = ParToP( par ); - //std::cout << " updates track param :" << std::endl; - //ptrack.Print(); - - // and (px, py) at the vertex instead of the dca : - double phi0 = par(1); - double phi = final_phases[i] ; - double px_at_vertex = ptrack.Pt() * TMath::Cos( phi0 + phi ); - double py_at_vertex = ptrack.Pt() * TMath::Sin( phi0 + phi ); - TVector3 ptrack_at_vertex( px_at_vertex, py_at_vertex, ptrack.Pz() ); - //std::cout << " momentum at the vertex : " << std::endl; - //std::cout << " phi0 at dca = " << phi0 << " phi at vertex = " << phi0+phi << " C = " << par(2) << " phase " << phi << std::endl; - //ptrack_at_vertex.Print(); - - updated_track_momentum_at_vertex.push_back( ptrack_at_vertex ); - - // back to EDM4HEP units... - par[0] = par[0] / scale0 ; - par[1] = par[1] / scale1 ; - par[2] = par[2] / scale2 ; - par[3] = par[3] / scale3 ; - par[4] = par[4] / scale4 ; - updated_track_parameters.push_back( par ); - - reco_chi2.push_back( final_chi2[i] ); - final_track_phases.push_back( final_phases[i] ); - - } - - TheVertex.vertex = result; - TheVertex.reco_chi2 = reco_chi2; - TheVertex.reco_ind = reco_ind; - TheVertex.updated_track_parameters = updated_track_parameters ; - TheVertex.updated_track_momentum_at_vertex = updated_track_momentum_at_vertex; - TheVertex.final_track_phases = final_track_phases; - - //std::cout << " end of VertexFitter " << std::endl; - /* - for ( Int_t i = 0; i < Ntr; i++) { - std::cout << " Track #" << i << " chi2 = " << reco_chi2[i] << std::endl; - std::cout << " Initial parameters: " << std::endl; - VertexingUtils::get_trackParam( tracks[i] ).Print(); - std::cout << " Updated parameters : " << std::endl; - updated_track_parameters[i].Print(); - } - */ - - delete[] final_chi2; - delete[] final_phases; - - return TheVertex; -} - - -//////////////////////////////////////////////////// - - - -ROOT::VecOps::RVec get_PrimaryTracks( VertexingUtils::FCCAnalysesVertex initialVertex, - ROOT::VecOps::RVec tracks, - bool BeamSpotConstraint, - double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, - double bsc_x, double bsc_y, double bsc_z, - int ipass ) { - - -// iterative procedure to determine the primary vertex - and the primary tracks -// Start from a vertex reconstructed from all tracks, remove the one with the highest chi2, fit again etc - -// tracks = the collection of tracks that was used in the first step - -//bool debug = true ; - bool debug = false; -float CHI2MAX = 25 ; - -if (debug) { - if (ipass == 0) std::cout << " \n --------------------------------------------------------\n" << std::endl; - std::cout << " ... enter in get_PrimaryTracks ipass = " << ipass << std::endl; - if (ipass == 0) std::cout << " initial number of tracks = " << tracks.size() << std::endl; -} - -ROOT::VecOps::RVec seltracks = tracks; -ROOT::VecOps::RVec reco_chi2 = initialVertex.reco_chi2; - -if ( seltracks.size() <= 1 ) return seltracks; - -int isPrimaryVertex = initialVertex.vertex.primary ; - -int maxElementIndex = std::max_element(reco_chi2.begin(),reco_chi2.end()) - reco_chi2.begin(); -auto minmax = std::minmax_element(reco_chi2.begin(), reco_chi2.end()); -float chi2max = *minmax.second ; - -if ( chi2max < CHI2MAX ) { - if (debug) { - std::cout << " --- DONE, all tracks have chi2 < CHI2MAX " << std::endl; - std::cout << " number of primary tracks selected = " << seltracks.size() << std::endl; - - } - return seltracks ; -} - - if (debug) { - std::cout << " remove a track that has chi2 = " << chi2max << std::endl; - } - -seltracks.erase( seltracks.begin() + maxElementIndex ); -ipass ++; - - VertexingUtils::FCCAnalysesVertex vtx = VertexFitter_Tk( isPrimaryVertex, - seltracks, - BeamSpotConstraint, - bsc_sigmax, bsc_sigmay, bsc_sigmaz, - bsc_x, bsc_y, bsc_z ) ; - - return get_PrimaryTracks( vtx, seltracks, BeamSpotConstraint, bsc_sigmax, bsc_sigmay, bsc_sigmaz, - bsc_x, bsc_y, bsc_z, ipass ) ; - - - -} - - -ROOT::VecOps::RVec get_NonPrimaryTracks( ROOT::VecOps::RVec allTracks, - ROOT::VecOps::RVec primaryTracks ) { - - ROOT::VecOps::RVec result; - for (auto & track: allTracks) { - bool isInPrimary = false; - for ( auto & primary: primaryTracks) { - if ( track.D0 == primary.D0 && track.Z0 == primary.Z0 && track.phi == primary.phi && track.omega == primary.omega && track.tanLambda == primary.tanLambda ) { - isInPrimary = true; - break; - } - } - if ( !isInPrimary) result.push_back( track ); - } - - return result; -} - - -ROOT::VecOps::RVec IsPrimary_forTracks( ROOT::VecOps::RVec allTracks, - ROOT::VecOps::RVec primaryTracks ) { - - ROOT::VecOps::RVec result; - for (auto & track: allTracks) { - bool isInPrimary = false; - for ( auto & primary: primaryTracks) { - if ( track.D0 == primary.D0 && track.Z0 == primary.Z0 && track.phi == primary.phi && track.omega == primary.omega && track.tanLambda == primary.tanLambda ) { - isInPrimary = true; - break; - } - } - result.push_back( isInPrimary ); - } - return result; -} - -}//end NS VertexFitterSimple - -}//end NS FCCAnalyses +#include "FCCAnalyses/VertexFitterSimple.h" +#include "FCCAnalyses/MCParticle.h" + +#include + +#include "TFile.h" +#include "TString.h" + +//#include "TrkUtil.h" // from delphes + +namespace FCCAnalyses { + +namespace VertexFitterSimple { + +// ----------------------------------------------------------------------------- + +VertexingUtils::FCCAnalysesVertex VertexFitter( + int Primary, + ROOT::VecOps::RVec recoparticles, + ROOT::VecOps::RVec thetracks, bool BeamSpotConstraint, + double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, double bsc_x, + double bsc_y, double bsc_z) { + + // input = a collection of recoparticles (in case one want to make + // associations to RecoParticles ?) and thetracks = the collection of all + // TrackState in the event + + VertexingUtils::FCCAnalysesVertex thevertex; + + // retrieve the tracks associated to the recoparticles + ROOT::VecOps::RVec tracks = + ReconstructedParticle2Track::getRP2TRK(recoparticles, thetracks); + + // and run the vertex fitter + + // FCCAnalysesVertex thevertex = VertexFitter_Tk( Primary, tracks, thetracks) + // ; + thevertex = + VertexFitter_Tk(Primary, tracks, thetracks, BeamSpotConstraint, + bsc_sigmax, bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z); + + // fill the indices of the tracks + ROOT::VecOps::RVec reco_ind; + int Ntr = tracks.size(); + for (auto &p : recoparticles) { + // std::cout << " in VertexFitter: a recoparticle with charge = " << + // p.charge << std::endl; + if (p.tracks_begin >= 0 && p.tracks_begin < thetracks.size()) { + reco_ind.push_back(p.tracks_begin); + } + } + if (reco_ind.size() != Ntr) + std::cout << " ... problem in Vertex, size of reco_ind != Ntr " + << std::endl; + + thevertex.reco_ind = reco_ind; + + return thevertex; +} + +// --------------------------------------------------------------------------------------------------------------------------- + +VertexingUtils::FCCAnalysesVertex +VertexFitter_Tk(int Primary, ROOT::VecOps::RVec tracks, + bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay, + double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z) { + + ROOT::VecOps::RVec dummy; + return VertexFitter_Tk(Primary, tracks, dummy, BeamSpotConstraint, bsc_sigmax, + bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z); +} + +// --------------------------------------------------------------------------------------------------------------------------- + +VertexingUtils::FCCAnalysesVertex +VertexFitter_Tk(int Primary, ROOT::VecOps::RVec tracks, + const ROOT::VecOps::RVec &alltracks, + bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay, + double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z) { + + // Units for the beam-spot : mum + // See + // https://github.com/HEP-FCC/FCCeePhysicsPerformance/tree/master/General#generating-events-under-realistic-fcc-ee-environment-conditions + + // final results : + VertexingUtils::FCCAnalysesVertex TheVertex; + + edm4hep::VertexData result; + ROOT::VecOps::RVec reco_chi2; + ROOT::VecOps::RVec updated_track_parameters; + ROOT::VecOps::RVec reco_ind; + ROOT::VecOps::RVec final_track_phases; + ROOT::VecOps::RVec updated_track_momentum_at_vertex; + + // if the collection of all tracks has been passed, keep trace of the indices + // of the tracks that are used to fit this vertex + if (alltracks.size() > 0) { + for (int i = 0; i < tracks.size(); i++) { // the fitted tracks + edm4hep::TrackState tr1 = tracks[i]; + for (int j = 0; j < alltracks.size(); + j++) { // the collection of all tracks + edm4hep::TrackState tr2 = alltracks[j]; + if (VertexingUtils::compare_Tracks(tr1, tr2)) { + reco_ind.push_back(j); + break; + } + } + } + } + + TheVertex.vertex = result; + TheVertex.reco_chi2 = reco_chi2; + TheVertex.updated_track_parameters = updated_track_parameters; + TheVertex.reco_ind = reco_ind; + TheVertex.final_track_phases = final_track_phases; + TheVertex.updated_track_momentum_at_vertex = updated_track_momentum_at_vertex; + + int Ntr = tracks.size(); + TheVertex.ntracks = Ntr; + if (Ntr <= 1) + return TheVertex; // can not reconstruct a vertex with only one track... + + TVectorD **trkPar = new TVectorD *[Ntr]; + TMatrixDSym **trkCov = new TMatrixDSym *[Ntr]; + + bool Units_mm = true; + + for (Int_t i = 0; i < Ntr; i++) { + edm4hep::TrackState t = tracks[i]; + TVectorD par = VertexingUtils::get_trackParam(t, Units_mm); + trkPar[i] = new TVectorD(par); + TMatrixDSym Cov = VertexingUtils::get_trackCov(t, Units_mm); + trkCov[i] = new TMatrixDSym(Cov); + } + + VertexFit theVertexFit(Ntr, trkPar, trkCov); + + if (BeamSpotConstraint) { + float conv_BSC = 1e-3; // convert mum to mm + TVectorD xv_BS(3); + xv_BS[0] = bsc_x * conv_BSC; + xv_BS[1] = bsc_y * conv_BSC; + xv_BS[2] = bsc_z * conv_BSC; + TMatrixDSym cov_BS(3); + cov_BS[0][0] = pow(bsc_sigmax * conv_BSC, 2); + cov_BS[1][1] = pow(bsc_sigmay * conv_BSC, 2); + cov_BS[2][2] = pow(bsc_sigmaz * conv_BSC, 2); + theVertexFit.AddVtxConstraint(xv_BS, cov_BS); + } + + TVectorD x = theVertexFit.GetVtx(); // this actually runs the fit + + result.position = + edm4hep::Vector3f(x(0), x(1), x(2)); // vertex position in mm + + // store the results in an edm4hep::VertexData object + + float Chi2 = theVertexFit.GetVtxChi2(); + float Ndof = 2.0 * Ntr - 3.0; + ; + result.chi2 = Chi2 / Ndof; + + // the chi2 of all the tracks : + TVectorD tracks_chi2 = theVertexFit.GetVtxChi2List(); + for (int it = 0; it < Ntr; it++) { + reco_chi2.push_back(tracks_chi2[it]); + } + + // std::cout << " Fitted vertex: " << x(0)*conv << " " << x(1)*conv << " " << + // x(2)*conv << std::endl; + TMatrixDSym covX = theVertexFit.GetVtxCov(); + std::array + covMatrix; // covMat in edm4hep is a LOWER-triangle matrix. + covMatrix[0] = covX(0, 0); + covMatrix[1] = covX(1, 0); + covMatrix[2] = covX(1, 1); + covMatrix[3] = covX(2, 0); + covMatrix[4] = covX(2, 1); + covMatrix[5] = covX(2, 2); + result.covMatrix = covMatrix; + + result.algorithmType = 1; + + result.primary = Primary; + + TheVertex.vertex = result; + + // Use VertexMore to retrieve more information : + VertexMore theVertexMore(&theVertexFit, Units_mm); + + for (Int_t i = 0; i < Ntr; i++) { + + TVectorD updated_par = + theVertexFit.GetNewPar(i); // updated track parameters + TVectorD updated_par_edm4hep = + VertexingUtils::Delphes2Edm4hep_TrackParam(updated_par, Units_mm); + updated_track_parameters.push_back(updated_par_edm4hep); + + // Momenta of the tracks at the vertex: + TVector3 ptrack_at_vertex = theVertexMore.GetMomentum(i); + updated_track_momentum_at_vertex.push_back(ptrack_at_vertex); + } + + TheVertex.updated_track_parameters = updated_track_parameters; + TheVertex.updated_track_momentum_at_vertex = updated_track_momentum_at_vertex; + TheVertex.final_track_phases = final_track_phases; + TheVertex.reco_chi2 = reco_chi2; + + // memory cleanup + for (Int_t i = 0; i < Ntr; i++) { + delete trkPar[i]; + delete trkCov[i]; + } + delete[] trkPar; + delete[] trkCov; + + return TheVertex; +} + +// --------------------------------------------------------------------------------------------------------------------------- + +ROOT::VecOps::RVec +get_PrimaryTracks(ROOT::VecOps::RVec tracks, + bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay, + double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z) { + + // iterative procedure to determine the primary vertex - and the primary + // tracks + + // Feb 2023: Avoid the recursive approach used before... else very very slow, + // with the new VertexFit objects + + // Units for the beam-spot : mum + // See + // https://github.com/HEP-FCC/FCCeePhysicsPerformance/tree/master/General#generating-events-under-realistic-fcc-ee-environment-conditions + + // bool debug = true ; + bool debug = false; + float CHI2MAX = 25; + // float CHI2MAX = 10; + + if (debug) { + std::cout << " ... enter in VertexFitterSimple::get_PrimaryTracks Ntr = " + << tracks.size() << std::endl; + } + + ROOT::VecOps::RVec seltracks = tracks; + + if (seltracks.size() <= 1) + return seltracks; + + int Ntr = tracks.size(); + + TVectorD **trkPar = new TVectorD *[Ntr]; + TMatrixDSym **trkCov = new TMatrixDSym *[Ntr]; + + for (Int_t i = 0; i < Ntr; i++) { + edm4hep::TrackState t = tracks[i]; + TVectorD par = VertexingUtils::get_trackParam(t); + trkPar[i] = new TVectorD(par); + TMatrixDSym Cov = VertexingUtils::get_trackCov(t); + trkCov[i] = new TMatrixDSym(Cov); + } + + VertexFit theVertexFit(Ntr, trkPar, trkCov); + + if (BeamSpotConstraint) { + TVectorD xv_BS(3); + xv_BS[0] = bsc_x * 1e-6; + xv_BS[1] = bsc_y * 1e-6; + xv_BS[2] = bsc_z * 1e-6; + TMatrixDSym cov_BS(3); + cov_BS[0][0] = pow(bsc_sigmax * 1e-6, 2); + cov_BS[1][1] = pow(bsc_sigmay * 1e-6, 2); + cov_BS[2][2] = pow(bsc_sigmaz * 1e-6, 2); + theVertexFit.AddVtxConstraint(xv_BS, cov_BS); + } + + TVectorD x = theVertexFit.GetVtx(); // this actually runs the fit + + float chi2_max = 1e30; + + while (chi2_max >= CHI2MAX) { + + TVectorD tracks_chi2 = theVertexFit.GetVtxChi2List(); + chi2_max = tracks_chi2.Max(); + + int n_removed = 0; + for (int i = 0; i < theVertexFit.GetNtrk(); i++) { + float track_chi2 = tracks_chi2[i]; + if (track_chi2 >= chi2_max) { + theVertexFit.RemoveTrk(i); + seltracks.erase(seltracks.begin() + i); + n_removed++; + } + } + if (n_removed > 0) { + if (theVertexFit.GetNtrk() > 1) { + // run the fit again: + x = theVertexFit.GetVtx(); + TVectorD new_tracks_chi2 = theVertexFit.GetVtxChi2List(); + chi2_max = new_tracks_chi2.Max(); + } else { + chi2_max = 0; // exit from the loop w/o crashing.. + } + } + } // end while + + // memory cleanup : + for (Int_t i = 0; i < Ntr; i++) { + delete trkPar[i]; + delete trkCov[i]; + } + delete[] trkPar; + delete[] trkCov; + + return seltracks; +} + +// --------------------------------------------------------------------------------------------------------------------------- + +ROOT::VecOps::RVec +get_NonPrimaryTracks(ROOT::VecOps::RVec allTracks, + ROOT::VecOps::RVec primaryTracks) { + + ROOT::VecOps::RVec result; + for (auto &track : allTracks) { + bool isInPrimary = false; + for (auto &primary : primaryTracks) { + if (VertexingUtils::compare_Tracks(track, primary)) { + isInPrimary = true; + break; + } + } + if (!isInPrimary) + result.push_back(track); + } + + return result; +} + +// --------------------------------------------------------------------------------------------------------------------------- + +ROOT::VecOps::RVec +IsPrimary_forTracks(ROOT::VecOps::RVec allTracks, + ROOT::VecOps::RVec primaryTracks) { + + ROOT::VecOps::RVec result; + for (auto &track : allTracks) { + bool isInPrimary = false; + for (auto &primary : primaryTracks) { + if (VertexingUtils::compare_Tracks(track, primary)) { + isInPrimary = true; + break; + } + } + result.push_back(isInPrimary); + } + return result; +} + +} // namespace VertexFitterSimple + +} // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/VertexingUtils.cc b/analyzers/dataframe/src/VertexingUtils.cc index b74a17bc76..9af606a7ab 100644 --- a/analyzers/dataframe/src/VertexingUtils.cc +++ b/analyzers/dataframe/src/VertexingUtils.cc @@ -1,1347 +1,1476 @@ -#include "FCCAnalyses/VertexingUtils.h" -#include "FCCAnalyses/VertexFitterSimple.h" -#include "FCCAnalyses/myUtils.h" - -namespace FCCAnalyses{ - -namespace VertexingUtils{ - -// -// Selection of particles based on the d0 / z0 significances of the associated track -// -selTracks::selTracks( float arg_d0sig_min, float arg_d0sig_max, float arg_z0sig_min, float arg_z0sig_max) : m_d0sig_min(arg_d0sig_min), - m_d0sig_max( arg_d0sig_max ), - m_z0sig_min( arg_z0sig_min ), - m_z0sig_max (arg_z0sig_max) { }; -ROOT::VecOps::RVec -selTracks::operator() (ROOT::VecOps::RVec recop, - ROOT::VecOps::RVec tracks ) { - - ROOT::VecOps::RVec result; - result.reserve(recop.size()); - - for (size_t i = 0; i < recop.size(); ++i) { - auto & p = recop[i]; - if (p.tracks_begin m_d0sig_max || fabs( d0sig ) < m_d0sig_min ) continue; - //double z0sig = fabs( tr.Z0 / sqrt( tr.covMatrix[12]) ); - double z0sig = fabs( tr.Z0 / sqrt( tr.covMatrix[9]) ); // covMat = lower-triangle - if ( fabs( z0sig ) > m_z0sig_max || fabs( z0sig ) < m_z0sig_min ) continue; - result.emplace_back(p); - } - } - return result; -} - - -// -// Selection of primary particles based on the matching of RecoParticles -// to MC particles -// -ROOT::VecOps::RVec -SelPrimaryTracks (ROOT::VecOps::RVec recind, ROOT::VecOps::RVec mcind, - ROOT::VecOps::RVec reco, - ROOT::VecOps::RVec mc, - TVector3 MC_EventPrimaryVertex) { - - ROOT::VecOps::RVec result; - result.reserve(reco.size()); - - // Event primary vertex: - double xvtx0 = MC_EventPrimaryVertex[0]; - double yvtx0 = MC_EventPrimaryVertex[1]; - double zvtx0 = MC_EventPrimaryVertex[2]; - - for (unsigned int i=0; i tracks) { - int nt = tracks.size(); - return nt; -} - -bool compare_Tracks( const edm4hep::TrackState& tr1, const edm4hep::TrackState& tr2 ) { - if ( tr1.D0 == tr2.D0 && - tr1.phi == tr2.phi && - tr1.omega == tr2.omega && - tr1.Z0 == tr2.Z0 && - tr1.tanLambda == tr2.tanLambda && - tr1.time == tr2.time && - tr1.referencePoint.x == tr2.referencePoint.x && - tr1.referencePoint.y == tr2.referencePoint.y && - tr1.referencePoint.z == tr2.referencePoint.z ) - return true; - return false; -} - - - -TVectorD -get_trackParam( edm4hep::TrackState & atrack) { - double d0 =atrack.D0 ; - double phi0 = atrack.phi ; - double omega = atrack.omega ; - double z0 = atrack.Z0 ; - double tanlambda = atrack.tanLambda ; - TVectorD res(5); - - double scale0 = 1e-3; //convert mm to m - double scale1 = 1; - double scale2 = 0.5*1e3; // C = rho/2, convert from mm-1 to m-1 - double scale3 = 1e-3 ; //convert mm to m - double scale4 = 1.; - - scale2 = -scale2 ; // sign of omega - - res[0] = d0 * scale0; - res[1] = phi0 * scale1 ; - res[2] = omega * scale2 ; - res[3] = z0 * scale3 ; - res[4] = tanlambda * scale4 ; - return res; -} - -float get_trackMom( edm4hep::TrackState & atrack ) { - double fB = 2; // 2 Tesla - - float C = -0.5*1e3 * atrack.omega; - float phi0 = atrack.phi; - float ct = atrack.tanLambda; - // - float pt = fB*0.2998 / TMath::Abs(2 * C); - TVector3 p(pt*TMath::Cos(phi0), pt*TMath::Sin(phi0), pt*ct); - float result = p.Mag(); - return result; -} - -TMatrixDSym -get_trackCov( edm4hep::TrackState & atrack) { - auto covMatrix = atrack.covMatrix; - TMatrixDSym covM(5); - - double scale0 = 1e-3; - double scale1 = 1.; - double scale2 = 0.5*1e3; - double scale3 = 1e-3 ; - double scale4 = 1.; - - scale2 = -scale2 ; // sign of omega - - // covMatrix = lower-triang;e - - covM[0][0] = covMatrix[0] *scale0 * scale0; - - covM[1][0] = covMatrix[1] *scale1 * scale0; - covM[1][1] = covMatrix[2] *scale1 * scale1; - - covM[0][1] = covM[1][0]; - - covM[2][0] = covMatrix[3] *scale2 * scale0; - covM[2][1] = covMatrix[4] *scale2 * scale1; - covM[2][2] = covMatrix[5] *scale2 * scale2; - - covM[0][2] = covM[2][0]; - covM[1][2] = covM[2][1]; - - covM[3][0] = covMatrix[6] *scale3 * scale0; - covM[3][1] = covMatrix[7] *scale3 * scale1; - covM[3][2] = covMatrix[8] *scale3 * scale2; - covM[3][3] = covMatrix[9] *scale3 * scale3; - - covM[0][3] = covM[3][0]; - covM[1][3] = covM[3][1]; - covM[2][3] = covM[3][2]; - - covM[4][0] = covMatrix[10] *scale4 * scale0; - covM[4][1] = covMatrix[11] *scale4 * scale1; - covM[4][2] = covMatrix[12] *scale4 * scale2; - covM[4][3] = covMatrix[13] *scale4 * scale3; - covM[4][4] = covMatrix[14] *scale4 * scale4; - - covM[0][4] = covM[4][0]; - covM[1][4] = covM[4][1]; - covM[2][4] = covM[4][2]; - covM[3][4] = covM[4][3]; - - return covM; -} - -// Following code is used in the Exotic Higgs Decays to LLPs analysis to find DVs with the SV finder of LCFI+ - -// function to merge vertices if the distance between them are less than 10*error on the position, or if they are within 1 mm from each other -ROOT::VecOps::RVec mergeVertices ( ROOT::VecOps::RVec vertices_in ) { - ROOT::VecOps::RVec result; - ROOT::VecOps::RVec merged_vertices; - int n = vertices_in.size(); - - // only try merging if there are 2 or more vertices - if (n == 0 || n == 1) { - result = vertices_in; - return result; - } - - // check distance between vertices pair-wise - for (int i=0; i < n-1; i++){ - FCCAnalysesVertex ivtx = vertices_in.at(i); - - for (int j=i+1; j < n; j++){ - FCCAnalysesVertex jvtx = vertices_in.at(j); - float distance = FCCAnalyses::myUtils::get_distanceVertex(ivtx.vertex, jvtx.vertex, -1); //input on edm4::vertexdata format, int comp = -1 to get distance in xyz - float error = FCCAnalyses::myUtils::get_distanceErrorVertex(ivtx.vertex, jvtx.vertex, -1); - float threshold = 10*error; - - if (distance < threshold || distance < 1){ - ROOT::VecOps::RVec tracks; - - for (edm4hep::TrackState track_i : ivtx.tracks) { - tracks.push_back(track_i); - } - - for (edm4hep::TrackState track_j : jvtx.tracks) { - tracks.push_back(track_j); - } - - // get the associated tracks to each vertex, combine and do a refit - FCCAnalysesVertex mergedVertex = VertexFitterSimple::VertexFitter_Tk(2, tracks); - - merged_vertices.push_back(mergedVertex); - - for (int k = 0; k < n; k++) { - if (k != i && k != j) { - merged_vertices.push_back(vertices_in.at(k)); - } - } - - //recursive process, repeat until vertices that all have distance > threshold are found - result = FCCAnalyses::VertexingUtils::mergeVertices(merged_vertices); - - return result; - } - } - - result.push_back(ivtx); - } - result.push_back(vertices_in.at(n-1)); - return result; -} - -// selection of tracks based on the transverse momentum pT -sel_pt_tracks::sel_pt_tracks(float arg_min_pt) : m_min_pt(arg_min_pt) {}; -ROOT::VecOps::RVec sel_pt_tracks::operator() (ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - auto & track = in[i]; - TVectorD track_param = FCCAnalyses::VertexingUtils::get_trackParam(track); //track parameters in edm4hep format - TVector3 tracks_p = FCCAnalyses::VertexFitterSimple::ParToP(track_param); //get the momentum of the tracks from the track_parameters - double tracks_pt = tracks_p.Pt(); //Pt: the transverse component - if (tracks_pt > m_min_pt) { - result.emplace_back(track); - } - } - return result; -} - -// selection of tracks based on the impact paramter d0 -sel_d0_tracks::sel_d0_tracks(float arg_min_d0) : m_min_d0(arg_min_d0) {}; -ROOT::VecOps::RVec sel_d0_tracks::operator() (ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - auto & track = in[i]; - double tr_d0 = fabs(track.D0); - if (tr_d0 > m_min_d0) { - result.emplace_back(track); - } - } - return result; -} - -// end of DV specific code - -FCCAnalysesVertex -get_FCCAnalysesVertex(ROOT::VecOps::RVec TheVertexColl, int index ){ - FCCAnalysesVertex result; - if (index TheVertexColl ){ - return TheVertexColl.size(); -} - - -edm4hep::VertexData get_VertexData( FCCAnalysesVertex TheVertex ) { - return TheVertex.vertex ; -} - -ROOT::VecOps::RVec get_VertexData( ROOT::VecOps::RVec TheVertexColl ) { - ROOT::VecOps::RVec result; - for (unsigned int i=0; i TheVertexColl, int index) { - edm4hep::VertexData result; - if (index get_VertexNtrk( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - for(auto & TheVertex: vertices){ - result.push_back(TheVertex.ntracks); - } - return result; -} - -ROOT::VecOps::RVec get_VertexRecoInd( FCCAnalysesVertex TheVertex ) { - return TheVertex.reco_ind; -} - -ROOT::VecOps::RVec get_VertexRecoParticlesInd( FCCAnalysesVertex TheVertex, - const ROOT::VecOps::RVec& reco ) { - - ROOT::VecOps::RVec result; - ROOT::VecOps::RVec indices_tracks = TheVertex.reco_ind; - for (int i=0; i < indices_tracks.size(); i++) { - int tk_index = indices_tracks[i]; - for ( int j=0; j < reco.size(); j++) { - auto & p = reco[j]; - if ( p.tracks_begin == p.tracks_end ) continue; - if ( p.tracks_begin == tk_index ) { - result.push_back( j ); - break; - } - } - } - return result; -} - - - - -TVectorD ParToACTS(TVectorD Par){ - - TVectorD pACTS(6); // Return vector - // - double fB=2.; - Double_t b = -0.29988*fB / 2.; - pACTS(0) = 1000*Par(0); // D from m to mm - pACTS(1) = 1000 * Par(3); // z0 from m to mm - pACTS(2) = Par(1); // Phi0 is unchanged - pACTS(3) = TMath::ATan2(1.0,Par(4)); // Theta in [0, pi] range - pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4))); // q/p in GeV - pACTS(5) = 0.0; // Time: currently undefined - // - return pACTS; -} - - -// Covariance conversion to ACTS format -TMatrixDSym CovToACTS(TMatrixDSym Cov, TVectorD Par){ - - double fB=2.; - TMatrixDSym cACTS(6); cACTS.Zero(); - Double_t b = -0.29988*fB / 2.; - // - // Fill derivative matrix - TMatrixD A(5, 5); A.Zero(); - Double_t ct = Par(4); // cot(theta) - Double_t C = Par(2); // half curvature - A(0, 0) = 1000.; // D-D conversion to mm - A(1, 2) = 1.0; // phi0-phi0 - A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b); // q/p-C - A(3, 1) = 1000.; // z0-z0 conversion to mm - A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta) - A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta) - // - TMatrixDSym Cv = Cov; - TMatrixD At(5, 5); - At.Transpose(A); - Cv.Similarity(At); - TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv; - cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes - // - return cACTS; -} - -//////////////////////////////////////////////////// - -// get all reconstructed vertices in a single vector -ROOT::VecOps::RVec get_all_vertices( FCCAnalysesVertex PV, - ROOT::VecOps::RVec SV ) { - // Returns a vector of all vertices (PV and SVs) - ROOT::VecOps::RVec result; - result.push_back(PV); - for (auto &p:SV){ - result.push_back(p); - } - return result; -} -// -ROOT::VecOps::RVec get_all_vertices( FCCAnalysesVertex PV, - ROOT::VecOps::RVec> SV ) { - // Returns a vector of all vertices (PV and SVs) - ROOT::VecOps::RVec result; - result.push_back(PV); - for (auto i_SV:SV){ - for (auto &p:i_SV){ - result.push_back(p); - } - } - return result; -} -// -// get all SVs in a single vector -ROOT::VecOps::RVec get_all_SVs( ROOT::VecOps::RVec> vertices ) { - // Returns a vector of all SVs - ROOT::VecOps::RVec result; - for (auto i_SV:vertices){ - for (auto &p:i_SV){ - result.push_back(p); - } - } - return result; -} - -// internal fns for SV finder - -// invariant mass of a two track vertex -double get_invM_pairs( FCCAnalysesVertex vertex, - double m1, - double m2 ) { - // CAUTION: m1 -> first track; m2 -> second track - - double result; - - ROOT::VecOps::RVec p_tracks = vertex.updated_track_momentum_at_vertex; - - TLorentzVector p4_vtx; - double m[2] = {m1, m2}; - int nTr = p_tracks.size(); - - for(unsigned int i=0; i get_invM_pairs( ROOT::VecOps::RVec vertices, - double m1, - double m2 ) { - // CAUTION: m1 -> first track; m2 -> second track - - ROOT::VecOps::RVec result; - for (auto & vertex: vertices) { - - double result_i; - ROOT::VecOps::RVec p_tracks = vertex.updated_track_momentum_at_vertex; - - TLorentzVector p4_vtx; - double m[2] = {m1, m2}; - int nTr = p_tracks.size(); - - for(unsigned int i=0; i p_tracks = vertex.updated_track_momentum_at_vertex; - - TLorentzVector p4_vtx; - const double m = 0.13957039; // pion mass - - for(TVector3 p_tr : p_tracks) { - TLorentzVector p4_tr; - p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); - p4_vtx += p4_tr; - } - - result = p4_vtx.M(); - return result; -} - -ROOT::VecOps::RVec get_invM( ROOT::VecOps::RVec vertices ){ - - ROOT::VecOps::RVec result; - for (auto & vertex: vertices) { - - double result_i; - ROOT::VecOps::RVec p_tracks = vertex.updated_track_momentum_at_vertex; - - TLorentzVector p4_vtx; - const double m = 0.13957039; // pion mass - - for(TVector3 p_tr : p_tracks) { - TLorentzVector p4_tr; - p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); - p4_vtx += p4_tr; - } - - result_i = p4_vtx.M(); - result.push_back(result_i); - } - return result; -} - -// cos(angle) b/n V0 candidate's (or any vtx) momentum & PV to V0 displacement vector -double get_PV2V0angle( FCCAnalysesVertex V0, - FCCAnalysesVertex PV ) { - double result; - - ROOT::VecOps::RVec p_tracks = V0.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - edm4hep::Vector3f r_V0 = V0.vertex.position; // in mm - edm4hep::Vector3f r_PV = PV.vertex.position; // in mm - - TVector3 r_V0_PV(r_V0[0] - r_PV[0], r_V0[1] - r_PV[1], r_V0[2] - r_PV[2]); - - double pDOTr = p_sum.Dot(r_V0_PV); - double p_mag = p_sum.Mag(); - double r_mag = r_V0_PV.Mag(); - - result = pDOTr / (p_mag * r_mag); - return result; -} - -// cos(angle) b/n track momentum sum & PV to vtx displacement vector -double get_PV2vtx_angle( ROOT::VecOps::RVec tracks, - FCCAnalysesVertex vtx, - FCCAnalysesVertex PV ) { - double result; - - TVector3 p_sum; - for(edm4hep::TrackState tr : tracks) { - TVectorD ipar = get_trackParam(tr); - TVector3 ip = VertexFitterSimple::ParToP(ipar); - p_sum += ip; - } - - edm4hep::Vector3f r_vtx = vtx.vertex.position; // in mm - edm4hep::Vector3f r_PV = PV.vertex.position; // in mm - - TVector3 r_vtx_PV(r_vtx[0] - r_PV[0], r_vtx[1] - r_PV[1], r_vtx[2] - r_PV[2]); - - double pDOTr = p_sum.Dot(r_vtx_PV); - double p_mag = p_sum.Mag(); - double r_mag = r_vtx_PV.Mag(); - - result = pDOTr / (p_mag * r_mag); - return result; -} - -// get track's energy assuming it to be a pion -double get_trackE( edm4hep::TrackState track ) { - - double result; - - const double m_pi = 0.13957039; - - TVectorD par = get_trackParam(track); - TVector3 p = VertexFitterSimple::ParToP(par); - TLorentzVector p4; - p4.SetXYZM(p[0], p[1], p[2], m_pi); - - result = p4.E(); - return result; -} - -//////////////////////////////////////////////// - -// no of reconstructed V0s -int get_n_SV( FCCAnalysesV0 SV ) { - int result = SV.vtx.size(); - return result; -} - -// vector of position of all reconstructed V0 (in mm) -ROOT::VecOps::RVec get_position_SV( FCCAnalysesV0 SV ) { - ROOT::VecOps::RVec result; - for(FCCAnalysesVertex ivtx : SV.vtx) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - result.push_back(xyz); - } - return result; -} - -// vector of PDG IDs of all reconstructed V0 -ROOT::VecOps::RVec get_pdg_V0( FCCAnalysesV0 V0 ) { - ROOT::VecOps::RVec result = V0.pdgAbs; - return result; -} - -// vector of invariant masses of all reconstructed V0 -ROOT::VecOps::RVec get_invM_V0( FCCAnalysesV0 V0 ) { - ROOT::VecOps::RVec result = V0.invM; - return result; -} - -// -// vector of momenta of all reconstructed V0 -ROOT::VecOps::RVec get_p_SV( FCCAnalysesV0 SV ) { - ROOT::VecOps::RVec result; - - for(FCCAnalysesVertex ivtx : SV.vtx) { - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - result.push_back(p_sum); - } - return result; -} - -// vector of chi2 of all reconstructed V0s -ROOT::VecOps::RVec get_chi2_SV( FCCAnalysesV0 SV ) { - ROOT::VecOps::RVec result; - - for(FCCAnalysesVertex ivtx : SV.vtx) { - int nDOF = 2*ivtx.ntracks - 3; - result.push_back(nDOF*ivtx.vertex.chi2); - } - return result; -} - -// passing a vector of FCCAnalysesVertex instead of new structs - -// no of reconstructed SVs -int get_n_SV( ROOT::VecOps::RVec vertices ) { - int result = vertices.size(); - return result; -} - -// vector of momenta of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_p_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - result.push_back(p_sum); - } - return result; -} - -// vector of position of all reconstructed SV (in mm) -ROOT::VecOps::RVec get_position_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - for(FCCAnalysesVertex ivtx : vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - result.push_back(xyz); - } - return result; -} - -// vector of momentum magnitude of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_pMag_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - result.push_back(p_sum.Mag()); - } - return result; -} - -// vector of chi2 of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_chi2_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - int nDOF = 2*ivtx.ntracks - 3; - result.push_back(nDOF*ivtx.vertex.chi2); - } - return result; -} - -// vector of chi2 (normalised) of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_norm_chi2_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) result.push_back(ivtx.vertex.chi2); - return result; -} - -// vector of nDOF of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_nDOF_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) result.push_back(2*ivtx.ntracks - 3); - return result; -} - -// vector of polar angle (theta) of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_theta_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - result.push_back(xyz.Theta()); - } - return result; -} - -// vector of azimuth angle (phi) of all reconstructed vertices (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_phi_SV( ROOT::VecOps::RVec vertices ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - result.push_back(xyz.Phi()); - } - return result; -} - -// vector of (cos of) angles b/n vtx momenta & PV to vtx displacement vectors -ROOT::VecOps::RVec get_pointingangle_SV( ROOT::VecOps::RVec vertices, - FCCAnalysesVertex PV ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - double iresult = 0.; - - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - edm4hep::Vector3f r_vtx = ivtx.vertex.position; // in mm - edm4hep::Vector3f r_PV = PV.vertex.position; // in mm - - TVector3 r_vtx_PV(r_vtx[0] - r_vtx[0], r_vtx[1] - r_PV[1], r_vtx[2] - r_PV[2]); - - double pDOTr = p_sum.Dot(r_vtx_PV); - double p_mag = p_sum.Mag(); - double r_mag = r_vtx_PV.Mag(); - - iresult = pDOTr / (p_mag * r_mag); - result.push_back(iresult); - } - return result; -} - -// vector of distances of all reconstructed SV from PV (in mm in xy plane) -ROOT::VecOps::RVec get_dxy_SV( ROOT::VecOps::RVec vertices, - FCCAnalysesVertex PV) { - ROOT::VecOps::RVec result; - TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], PV.vertex.position[2]); - for(auto & ivtx : vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 x_vtx_PV = x_vtx - x_PV; - - result.push_back(x_vtx_PV.Perp()); - } - return result; -} - -// vector of distances of all reconstructed SV from PV (in mm in 3D) -ROOT::VecOps::RVec get_d3d_SV( ROOT::VecOps::RVec vertices, - FCCAnalysesVertex PV) { - ROOT::VecOps::RVec result; - TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], PV.vertex.position[2]); - for(auto & ivtx : vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 x_vtx_PV = x_vtx - x_PV; - - result.push_back(x_vtx_PV.Mag()); - } - return result; -} - -// vector of distances of all reconstructed SV from given TVector3d (in mm in 3D) -ROOT::VecOps::RVec get_d3d_SV_obj( ROOT::VecOps::RVec vertices, - TVector3 object) { - ROOT::VecOps::RVec result; - for(auto & ivtx : vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - x_vtx = x_vtx - object; - - result.push_back(x_vtx.Mag()); - } - return result; -} - -// vector of distances of all reconstructed SV from given edm4hep::Vector3d (in mm in 3D) -ROOT::VecOps::RVec get_d3d_SV_obj( ROOT::VecOps::RVec vertices, - edm4hep::Vector3d object) { - ROOT::VecOps::RVec result; - for(auto & ivtx : vertices) { - double dx = ivtx.vertex.position[0] - object.x; - double dy = ivtx.vertex.position[1] - object.y; - double dz = ivtx.vertex.position[2] - object.z; - - TVector3 d3d(dx, dy, dz); - - result.push_back(d3d.Mag()); - } - return result; -} - -// vector of decay position distances of all reconstructed SV from given TVector3d (in mm in 3D) -ROOT::VecOps::RVec get_dR_SV_obj( ROOT::VecOps::RVec vertices, - TVector3 object) { - ROOT::VecOps::RVec result; - for(auto & ivtx : vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - result.push_back(1e3*TMath::Sqrt( pow(ivtx.vertex.position[0],2) + pow(ivtx.vertex.position[1],2) + pow(ivtx.vertex.position[2],2) ) - 1e3*TMath::Sqrt( pow(object.x(),2) + pow(object.y(),2) + pow(object.z(),2) )); - } - return result; -} - -// vector of decay position distances of all reconstructed SV from given edm4hep::Vector3d (in mm in 3D) -ROOT::VecOps::RVec get_dR_SV_obj( ROOT::VecOps::RVec vertices, - edm4hep::Vector3d object) { - ROOT::VecOps::RVec result; - for(auto & ivtx : vertices) { - result.push_back(1e3*TMath::Sqrt( pow(ivtx.vertex.position[0],2) + pow(ivtx.vertex.position[1],2) + pow(ivtx.vertex.position[2],2) ) - 1e3*TMath::Sqrt( pow(object.x,2) + pow(object.y,2) + pow(object.z,2) )); - } - return result; -} - -// vector of polar angle (theta) of all reconstructed vertices wrt jet axis (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_relTheta_SV( ROOT::VecOps::RVec vertices, - ROOT::VecOps::RVec nSV_jet, - ROOT::VecOps::RVec jets ) { - ROOT::VecOps::RVec result; - - unsigned int j = 0; - int nSV = nSV_jet[0]; - for(unsigned int i=0; i= nSV) { - j++; - nSV += nSV_jet[j]; - } - auto & ijet = jets[j]; - double jetTheta = ijet.theta(); - - result.push_back(xyz.Theta() - jetTheta); - } - return result; -} - -// vector of azimuthal angle (phi) of all reconstructed vertices wrt jet axis (SV.vtx or V0.vtx) -ROOT::VecOps::RVec get_relPhi_SV( ROOT::VecOps::RVec vertices, - ROOT::VecOps::RVec nSV_jet, - ROOT::VecOps::RVec jets ) { - ROOT::VecOps::RVec result; - - unsigned int j = 0; - int nSV = nSV_jet[0]; - for(unsigned int i=0; i= nSV) { - j++; - nSV += nSV_jet[j]; - } - auto & ijet = jets[j]; - TVector3 jetP(ijet.px(), ijet.py(), ijet.pz()); - - result.push_back(xyz.DeltaPhi(jetP)); - } - return result; -} - -// For get_SV_jets outputs - -// no of reconstructed SVs -int get_n_SV( ROOT::VecOps::RVec> vertices ) { - int result = 0; - if(vertices.size() != 0) { - for(auto SV_jets : vertices) result += SV_jets.size(); - } - return result; -} - -// Return the number of reconstructed SVs -ROOT::VecOps::RVec get_n_SV_jets( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec result; - if(vertices.size() != 0) { - for(auto SV_jets : vertices) result.push_back(SV_jets.size()); - } - return result; -} - -// -// separate V0s by jets -ROOT::VecOps::RVec> get_svInJets( ROOT::VecOps::RVec vertices, - ROOT::VecOps::RVec nSV_jet ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - int index=0; - for(unsigned int i : nSV_jet) { - for(unsigned int j=0; j> get_tracksInJets( ROOT::VecOps::RVec recoparticles, - ROOT::VecOps::RVec thetracks, - ROOT::VecOps::RVec jets, - std::vector> jet_consti ) { - std::vector> result; - std::vector iJet_tracks; - - int nJet = jets.size(); - // - for (unsigned int j=0; j i_jetconsti = jet_consti[j]; - - for(unsigned int ip : i_jetconsti) { - auto & p = recoparticles[ip]; - if(p.tracks_begin >= 0 && p.tracks_begin get_relTheta_SV( ROOT::VecOps::RVec vertices, fastjet::PseudoJet jet ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - // - result.push_back(xyz.Theta() - jet.theta()); - } - // - return result; -} - -// vector of azimuthal angle (phi) of all reconstructed vertices wrt jet axis [only for vertices from 1 jet] -ROOT::VecOps::RVec get_relPhi_SV( ROOT::VecOps::RVec vertices, fastjet::PseudoJet jet ) { - ROOT::VecOps::RVec result; - - for(auto & ivtx : vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 jetP(jet.px(), jet.py(), jet.pz()); - // - result.push_back(xyz.DeltaPhi(jetP)); - } - // - return result; -} - - -/////// vec of vec functions (for get_SV_jets) ///////// - -// SV invariant mass -ROOT::VecOps::RVec> get_invM( ROOT::VecOps::RVec> vertices ){ - - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - - for (auto & vertex: i_vertices) { - ROOT::VecOps::RVec p_tracks = vertex.updated_track_momentum_at_vertex; - // - TLorentzVector p4_vtx; - const double m = 0.13957039; // pion mass - // - for(TVector3 p_tr : p_tracks) { - TLorentzVector p4_tr; - p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); - p4_vtx += p4_tr; - } - i_result.push_back(p4_vtx.M()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV momentum -ROOT::VecOps::RVec> get_p_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - i_result.push_back(p_sum); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV momentum magnitude -ROOT::VecOps::RVec> get_pMag_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - i_result.push_back(p_sum.Mag()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV daughters multiplicity -ROOT::VecOps::RVec> get_VertexNtrk( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & TheVertex: i_vertices){ - i_result.push_back(TheVertex.ntracks); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV chi2 -ROOT::VecOps::RVec> get_chi2_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - int nDOF = 2*ivtx.ntracks - 3; - i_result.push_back(nDOF*ivtx.vertex.chi2); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV normalised chi2 -ROOT::VecOps::RVec> get_norm_chi2_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) i_result.push_back(ivtx.vertex.chi2); - result.push_back(i_result); - i_result.clear(); - } -return result; -} - -// SV no of DOF -ROOT::VecOps::RVec> get_nDOF_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) i_result.push_back(2*ivtx.ntracks - 3); - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV theta -ROOT::VecOps::RVec> get_theta_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - - for(auto & ivtx : i_vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - i_result.push_back(xyz.Theta()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV phi -ROOT::VecOps::RVec> get_phi_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - - for(auto & ivtx : i_vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - i_result.push_back(xyz.Phi()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV relative theta -ROOT::VecOps::RVec> get_relTheta_SV( ROOT::VecOps::RVec> vertices, - ROOT::VecOps::RVec jets ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - fastjet::PseudoJet i_jet = jets.at(i); - for(auto & ivtx : i_vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - // - i_result.push_back(xyz.Theta() - i_jet.theta()); - } - result.push_back(i_result); - i_result.clear(); - } - // - return result; -} - -// SV relative phi -ROOT::VecOps::RVec> get_relPhi_SV( ROOT::VecOps::RVec> vertices, - ROOT::VecOps::RVec jets ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_vertices = vertices.at(i); - fastjet::PseudoJet i_jet = jets.at(i); - for(auto & ivtx : i_vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 jetP(i_jet.px(), i_jet.py(), i_jet.pz()); - // - i_result.push_back(xyz.DeltaPhi(jetP)); - } - result.push_back(i_result); - i_result.clear(); - } - // - return result; -} - -// SV pointing angle wrt PV -ROOT::VecOps::RVec> get_pointingangle_SV( ROOT::VecOps::RVec> vertices, - FCCAnalysesVertex PV ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - edm4hep::Vector3f r_PV = PV.vertex.position; // in mm - - for(unsigned int i=0; i i_vertices = vertices.at(i); - for(auto & ivtx : i_vertices) { - double pointangle = 0.; - - ROOT::VecOps::RVec p_tracks = ivtx.updated_track_momentum_at_vertex; - TVector3 p_sum; - for(TVector3 p_tr : p_tracks) p_sum += p_tr; - - edm4hep::Vector3f r_vtx = ivtx.vertex.position; // in mm - - TVector3 r_vtx_PV(r_vtx[0] - r_vtx[0], r_vtx[1] - r_PV[1], r_vtx[2] - r_PV[2]); - - double pDOTr = p_sum.Dot(r_vtx_PV); - double p_mag = p_sum.Mag(); - double r_mag = r_vtx_PV.Mag(); - - pointangle = pDOTr / (p_mag * r_mag); - i_result.push_back(pointangle); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV distance from PV in xy -ROOT::VecOps::RVec> get_dxy_SV( ROOT::VecOps::RVec> vertices, - FCCAnalysesVertex PV) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], PV.vertex.position[2]); - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 x_vtx_PV = x_vtx - x_PV; - - i_result.push_back(x_vtx_PV.Perp()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV distance from PV in 3D -ROOT::VecOps::RVec> get_d3d_SV( ROOT::VecOps::RVec> vertices, - FCCAnalysesVertex PV) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], PV.vertex.position[2]); - - for(unsigned int i=0; i i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - TVector3 x_vtx_PV = x_vtx - x_PV; - - i_result.push_back(x_vtx_PV.Mag()); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// SV position in 3D -ROOT::VecOps::RVec> get_position_SV( ROOT::VecOps::RVec> vertices ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - for(unsigned int i=0; i i_result; - ROOT::VecOps::RVec i_vertices = vertices.at(i); - // - for(auto & ivtx : i_vertices) { - TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], ivtx.vertex.position[2]); - i_result.push_back(xyz); - } - result.push_back(i_result); - i_result.clear(); - } - return result; -} - -// V0 pdg -ROOT::VecOps::RVec> get_pdg_V0( ROOT::VecOps::RVec pdg, - ROOT::VecOps::RVec nSV_jet ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - int index=0; - for(unsigned int i : nSV_jet) { - for(unsigned int j=0; j> get_invM_V0( ROOT::VecOps::RVec invM, - ROOT::VecOps::RVec nSV_jet ) { - ROOT::VecOps::RVec> result; - ROOT::VecOps::RVec i_result; - - int index=0; - for(unsigned int i : nSV_jet) { - for(unsigned int j=0; j selTracks::operator()( + ROOT::VecOps::RVec recop, + ROOT::VecOps::RVec tracks) { + + ROOT::VecOps::RVec result; + result.reserve(recop.size()); + + for (size_t i = 0; i < recop.size(); ++i) { + auto &p = recop[i]; + if (p.tracks_begin < tracks.size()) { + auto &tr = tracks.at(p.tracks_begin); + double d0sig = fabs(tr.D0 / sqrt(tr.covMatrix[0])); + if (fabs(d0sig) > m_d0sig_max || fabs(d0sig) < m_d0sig_min) + continue; + // double z0sig = fabs( tr.Z0 / sqrt( tr.covMatrix[12]) ); + double z0sig = + fabs(tr.Z0 / sqrt(tr.covMatrix[9])); // covMat = lower-triangle + if (fabs(z0sig) > m_z0sig_max || fabs(z0sig) < m_z0sig_min) + continue; + result.emplace_back(p); + } + } + return result; +} + +// +// Selection of primary particles based on the matching of RecoParticles +// to MC particles +// +ROOT::VecOps::RVec +SelPrimaryTracks(ROOT::VecOps::RVec recind, ROOT::VecOps::RVec mcind, + ROOT::VecOps::RVec reco, + ROOT::VecOps::RVec mc, + TVector3 MC_EventPrimaryVertex) { + + ROOT::VecOps::RVec result; + result.reserve(reco.size()); + + // Event primary vertex: + double xvtx0 = MC_EventPrimaryVertex[0]; + double yvtx0 = MC_EventPrimaryVertex[1]; + double zvtx0 = MC_EventPrimaryVertex[2]; + + for (unsigned int i = 0; i < recind.size(); i++) { + double xvtx = mc.at(mcind.at(i)).vertex.x; + double yvtx = mc.at(mcind.at(i)).vertex.y; + double zvtx = mc.at(mcind.at(i)).vertex.z; + // primary particle ? + double zero = 1e-12; + if (fabs(xvtx - xvtx0) < zero && fabs(yvtx - yvtx0) < zero && + fabs(zvtx - zvtx0) < zero) { + int reco_idx = recind.at(i); + result.push_back(reco.at(reco_idx)); + } + } + return result; +} + +int get_nTracks(ROOT::VecOps::RVec tracks) { + int nt = tracks.size(); + return nt; +} + +bool compare_Tracks(const edm4hep::TrackState &tr1, + const edm4hep::TrackState &tr2) { + if (tr1.D0 == tr2.D0 && tr1.phi == tr2.phi && tr1.omega == tr2.omega && + tr1.Z0 == tr2.Z0 && tr1.tanLambda == tr2.tanLambda && + tr1.time == tr2.time && tr1.referencePoint.x == tr2.referencePoint.x && + tr1.referencePoint.y == tr2.referencePoint.y && + tr1.referencePoint.z == tr2.referencePoint.z) + return true; + return false; +} + +// ---------------------------------------------------------------------------------------- + +// --- Conversion methods between the Delphes and edm4hep conventions + +TVectorD Edm4hep2Delphes_TrackParam(const TVectorD ¶m, bool Units_mm) { + + double conv = 1e-3; // convert mm to m + if (Units_mm) + conv = 1.; + + double scale0 = conv; // convert mm to m if needed + double scale1 = 1; + double scale2 = 0.5 / conv; // C = rho/2, convert from mm-1 to m-1 + double scale3 = conv; // convert mm to m + double scale4 = 1.; + scale2 = -scale2; // sign of omega + + TVectorD result(5); + result[0] = param[0] * scale0; + result[1] = param[1] * scale1; + result[2] = param[2] * scale2; + result[3] = param[3] * scale3; + result[4] = param[4] * scale4; + + return result; +} + +TVectorD Delphes2Edm4hep_TrackParam(const TVectorD ¶m, bool Units_mm) { + + double conv = 1e-3; // convert mm to m + if (Units_mm) + conv = 1.; + + double scale0 = conv; // convert mm to m if needed + double scale1 = 1; + double scale2 = 0.5 / conv; // C = rho/2, convert from mm-1 to m-1 + double scale3 = conv; // convert mm to m + double scale4 = 1.; + scale2 = -scale2; // sign of omega + + TVectorD result(5); + result[0] = param[0] / scale0; + result[1] = param[1] / scale1; + result[2] = param[2] / scale2; + result[3] = param[3] / scale3; + result[4] = param[4] / scale4; + + return result; +} + +TMatrixDSym +Edm4hep2Delphes_TrackCovMatrix(const std::array &covMatrix, + bool Units_mm) { + + // careful: since summer 2022, the covariance matrix in edm4hep is an array of + // 21 floats because the time has been added as a 6th track parameter. But + // Delphes (and k4marlinWrapper) samples still fill it as an array of 15 + // floats. + + double conv = 1e-3; // convert mm to m + if (Units_mm) + conv = 1.; + + double scale0 = conv; // convert mm to m if needed + double scale1 = 1; + double scale2 = 0.5 / conv; // C = rho/2, convert from mm-1 to m-1 + double scale3 = conv; // convert mm to m + double scale4 = 1.; + scale2 = -scale2; // sign of omega + + TMatrixDSym covM(5); + + covM[0][0] = covMatrix[0] * scale0 * scale0; + + covM[1][0] = covMatrix[1] * scale1 * scale0; + covM[1][1] = covMatrix[2] * scale1 * scale1; + + covM[0][1] = covM[1][0]; + + covM[2][0] = covMatrix[3] * scale2 * scale0; + covM[2][1] = covMatrix[4] * scale2 * scale1; + covM[2][2] = covMatrix[5] * scale2 * scale2; + + covM[0][2] = covM[2][0]; + covM[1][2] = covM[2][1]; + + covM[3][0] = covMatrix[6] * scale3 * scale0; + covM[3][1] = covMatrix[7] * scale3 * scale1; + covM[3][2] = covMatrix[8] * scale3 * scale2; + covM[3][3] = covMatrix[9] * scale3 * scale3; + + covM[0][3] = covM[3][0]; + covM[1][3] = covM[3][1]; + covM[2][3] = covM[3][2]; + + covM[4][0] = covMatrix[10] * scale4 * scale0; + covM[4][1] = covMatrix[11] * scale4 * scale1; + covM[4][2] = covMatrix[12] * scale4 * scale2; + covM[4][3] = covMatrix[13] * scale4 * scale3; + covM[4][4] = covMatrix[14] * scale4 * scale4; + + covM[0][4] = covM[4][0]; + covM[1][4] = covM[4][1]; + covM[2][4] = covM[4][2]; + covM[3][4] = covM[4][3]; + + return covM; +} + +std::array Delphes2Edm4hep_TrackCovMatrix(const TMatrixDSym &cov, + bool Units_mm) { + + // careful: since summer 2022, the covariance matrix in edm4hep is an array of + // 21 floats because the time has been added as a 6th track parameter. But + // Delphes (and k4marlinWrapper) samples still fill it as an array of 15 + // floats. + + double conv = 1e-3; // convert mm to m + if (Units_mm) + conv = 1.; + + double scale0 = conv; // convert mm to m if needed + double scale1 = 1; + double scale2 = 0.5 / conv; // C = rho/2, convert from mm-1 to m-1 + double scale3 = conv; // convert mm to m + double scale4 = 1.; + scale2 = -scale2; // sign of omega + + std::array covMatrix; + covMatrix[0] = cov[0][0] / (scale0 * scale0); + covMatrix[1] = cov[1][0] / (scale1 * scale0); + covMatrix[2] = cov[1][1] / (scale1 * scale1); + covMatrix[3] = cov[2][0] / (scale0 * scale2); + covMatrix[4] = cov[2][1] / (scale1 * scale2); + covMatrix[5] = cov[2][2] / (scale2 * scale2); + covMatrix[6] = cov[3][0] / (scale3 * scale0); + covMatrix[7] = cov[3][1] / (scale3 * scale1); + covMatrix[8] = cov[3][2] / (scale3 * scale2); + covMatrix[9] = cov[3][3] / (scale3 * scale3); + covMatrix[10] = cov[4][0] / (scale4 * scale0); + covMatrix[11] = cov[4][1] / (scale4 * scale1); + covMatrix[12] = cov[4][2] / (scale4 * scale2); + covMatrix[13] = cov[4][3] / (scale4 * scale3); + covMatrix[14] = cov[4][4] / (scale4 * scale4); + + for (int i = 15; i < 21; i++) { + covMatrix[i] = 0.; + } + + return covMatrix; +} + +TVectorD get_trackParam(edm4hep::TrackState &atrack, bool Units_mm) { + double d0 = atrack.D0; + double phi0 = atrack.phi; + double omega = atrack.omega; + double z0 = atrack.Z0; + double tanlambda = atrack.tanLambda; + TVectorD param(5); + param[0] = d0; + param[1] = phi0; + param[2] = omega; + param[3] = z0; + param[4] = tanlambda; + TVectorD res = Edm4hep2Delphes_TrackParam(param, Units_mm); + + return res; +} + +TMatrixDSym get_trackCov(edm4hep::TrackState &atrack, bool Units_mm) { + auto covMatrix = atrack.covMatrix; + + TMatrixDSym covM = Edm4hep2Delphes_TrackCovMatrix(covMatrix, Units_mm); + + return covM; +} + +// ---------------------------------------------------------------------------------------- + +float get_trackMom(edm4hep::TrackState &atrack) { + double fB = 2; // 2 Tesla + + float C = -0.5 * 1e3 * atrack.omega; + float phi0 = atrack.phi; + float ct = atrack.tanLambda; + // + float pt = fB * 0.2998 / TMath::Abs(2 * C); + TVector3 p(pt * TMath::Cos(phi0), pt * TMath::Sin(phi0), pt * ct); + float result = p.Mag(); + return result; +} + +FCCAnalysesVertex +get_FCCAnalysesVertex(ROOT::VecOps::RVec TheVertexColl, + int index) { + FCCAnalysesVertex result; + if (index < TheVertexColl.size()) + result = TheVertexColl.at(index); + return result; +} + +int get_Nvertex(ROOT::VecOps::RVec TheVertexColl) { + return TheVertexColl.size(); +} + +edm4hep::VertexData get_VertexData(FCCAnalysesVertex TheVertex) { + return TheVertex.vertex; +} + +ROOT::VecOps::RVec +get_VertexData(ROOT::VecOps::RVec TheVertexColl) { + ROOT::VecOps::RVec result; + for (unsigned int i = 0; i < TheVertexColl.size(); i++) { + result.push_back(TheVertexColl.at(i).vertex); + } + return result; +} + +edm4hep::VertexData +get_VertexData(ROOT::VecOps::RVec TheVertexColl, int index) { + edm4hep::VertexData result; + if (index < TheVertexColl.size()) + result = TheVertexColl.at(index).vertex; + return result; +} + +int get_VertexNtrk(FCCAnalysesVertex TheVertex) { return TheVertex.ntracks; } + +ROOT::VecOps::RVec +get_VertexNtrk(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + for (auto &TheVertex : vertices) { + result.push_back(TheVertex.ntracks); + } + return result; +} + +ROOT::VecOps::RVec get_VertexRecoInd(FCCAnalysesVertex TheVertex) { + return TheVertex.reco_ind; +} + +ROOT::VecOps::RVec get_VertexRecoParticlesInd( + FCCAnalysesVertex TheVertex, + const ROOT::VecOps::RVec &reco) { + + ROOT::VecOps::RVec result; + ROOT::VecOps::RVec indices_tracks = TheVertex.reco_ind; + for (int i = 0; i < indices_tracks.size(); i++) { + int tk_index = indices_tracks[i]; + for (int j = 0; j < reco.size(); j++) { + auto &p = reco[j]; + if (p.tracks_begin == p.tracks_end) + continue; + if (p.tracks_begin == tk_index) { + result.push_back(j); + break; + } + } + } + return result; +} + +TVectorD ParToACTS(TVectorD Par) { + + TVectorD pACTS(6); // Return vector + // + double fB = 2.; + Double_t b = -0.29988 * fB / 2.; + pACTS(0) = 1000 * Par(0); // D from m to mm + pACTS(1) = 1000 * Par(3); // z0 from m to mm + pACTS(2) = Par(1); // Phi0 is unchanged + pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range + pACTS(4) = Par(2) / (b * TMath::Sqrt(1 + Par(4) * Par(4))); // q/p in GeV + pACTS(5) = 0.0; // Time: currently undefined + // + return pACTS; +} + +// Covariance conversion to ACTS format +TMatrixDSym CovToACTS(TMatrixDSym Cov, TVectorD Par) { + + double fB = 2.; + TMatrixDSym cACTS(6); + cACTS.Zero(); + Double_t b = -0.29988 * fB / 2.; + // + // Fill derivative matrix + TMatrixD A(5, 5); + A.Zero(); + Double_t ct = Par(4); // cot(theta) + Double_t C = Par(2); // half curvature + A(0, 0) = 1000.; // D-D conversion to mm + A(1, 2) = 1.0; // phi0-phi0 + A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b); // q/p-C + A(3, 1) = 1000.; // z0-z0 conversion to mm + A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta) + A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta) + // + TMatrixDSym Cv = Cov; + TMatrixD At(5, 5); + At.Transpose(A); + Cv.Similarity(At); + TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv; + cACTS(5, 5) = + 0.1; // Currently undefined: set to arbitrary value to avoid crashes + // + return cACTS; +} + +//////////////////////////////////////////////////// + +// get all reconstructed vertices in a single vector +ROOT::VecOps::RVec +get_all_vertices(FCCAnalysesVertex PV, + ROOT::VecOps::RVec SV) { + // Returns a vector of all vertices (PV and SVs) + ROOT::VecOps::RVec result; + result.push_back(PV); + for (auto &p : SV) { + result.push_back(p); + } + return result; +} +// +ROOT::VecOps::RVec +get_all_vertices(FCCAnalysesVertex PV, + ROOT::VecOps::RVec> SV) { + // Returns a vector of all vertices (PV and SVs) + ROOT::VecOps::RVec result; + result.push_back(PV); + for (auto i_SV : SV) { + for (auto &p : i_SV) { + result.push_back(p); + } + } + return result; +} +// +// get all SVs in a single vector +ROOT::VecOps::RVec get_all_SVs( + ROOT::VecOps::RVec> vertices) { + // Returns a vector of all SVs + ROOT::VecOps::RVec result; + for (auto i_SV : vertices) { + for (auto &p : i_SV) { + result.push_back(p); + } + } + return result; +} + +// internal fns for SV finder + +// invariant mass of a two track vertex +double get_invM_pairs(FCCAnalysesVertex vertex, double m1, double m2) { + // CAUTION: m1 -> first track; m2 -> second track + + double result; + + ROOT::VecOps::RVec p_tracks = + vertex.updated_track_momentum_at_vertex; + + TLorentzVector p4_vtx; + double m[2] = {m1, m2}; + int nTr = p_tracks.size(); + + for (unsigned int i = 0; i < nTr; i++) { + TLorentzVector p4_tr; + p4_tr.SetXYZM(p_tracks[i].X(), p_tracks[i].Y(), p_tracks[i].Z(), m[i]); + p4_vtx += p4_tr; + } + + result = p4_vtx.M(); + return result; +} + +ROOT::VecOps::RVec +get_invM_pairs(ROOT::VecOps::RVec vertices, double m1, + double m2) { + // CAUTION: m1 -> first track; m2 -> second track + + ROOT::VecOps::RVec result; + for (auto &vertex : vertices) { + + double result_i; + ROOT::VecOps::RVec p_tracks = + vertex.updated_track_momentum_at_vertex; + + TLorentzVector p4_vtx; + double m[2] = {m1, m2}; + int nTr = p_tracks.size(); + + for (unsigned int i = 0; i < nTr; i++) { + TLorentzVector p4_tr; + p4_tr.SetXYZM(p_tracks[i].X(), p_tracks[i].Y(), p_tracks[i].Z(), m[i]); + p4_vtx += p4_tr; + } + + result_i = p4_vtx.M(); + result.push_back(result_i); + } + return result; +} + +// invariant mass of a vertex (assuming all tracks to be pions) +double get_invM(FCCAnalysesVertex vertex) { + + double result; + + ROOT::VecOps::RVec p_tracks = + vertex.updated_track_momentum_at_vertex; + + TLorentzVector p4_vtx; + const double m = 0.13957039; // pion mass + + for (TVector3 p_tr : p_tracks) { + TLorentzVector p4_tr; + p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); + p4_vtx += p4_tr; + } + + result = p4_vtx.M(); + return result; +} + +ROOT::VecOps::RVec +get_invM(ROOT::VecOps::RVec vertices) { + + ROOT::VecOps::RVec result; + for (auto &vertex : vertices) { + + double result_i; + ROOT::VecOps::RVec p_tracks = + vertex.updated_track_momentum_at_vertex; + + TLorentzVector p4_vtx; + const double m = 0.13957039; // pion mass + + for (TVector3 p_tr : p_tracks) { + TLorentzVector p4_tr; + p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); + p4_vtx += p4_tr; + } + + result_i = p4_vtx.M(); + result.push_back(result_i); + } + return result; +} + +// cos(angle) b/n V0 candidate's (or any vtx) momentum & PV to V0 displacement +// vector +double get_PV2V0angle(FCCAnalysesVertex V0, FCCAnalysesVertex PV) { + double result; + + ROOT::VecOps::RVec p_tracks = V0.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + edm4hep::Vector3f r_V0 = V0.vertex.position; // in mm + edm4hep::Vector3f r_PV = PV.vertex.position; // in mm + + TVector3 r_V0_PV(r_V0[0] - r_PV[0], r_V0[1] - r_PV[1], r_V0[2] - r_PV[2]); + + double pDOTr = p_sum.Dot(r_V0_PV); + double p_mag = p_sum.Mag(); + double r_mag = r_V0_PV.Mag(); + + result = pDOTr / (p_mag * r_mag); + return result; +} + +// cos(angle) b/n track momentum sum & PV to vtx displacement vector +double get_PV2vtx_angle(ROOT::VecOps::RVec tracks, + FCCAnalysesVertex vtx, FCCAnalysesVertex PV) { + double result; + + TVector3 p_sum; + for (edm4hep::TrackState tr : tracks) { + TVectorD ipar = get_trackParam(tr); + TVector3 ip = ParToP(ipar); + p_sum += ip; + } + + edm4hep::Vector3f r_vtx = vtx.vertex.position; // in mm + edm4hep::Vector3f r_PV = PV.vertex.position; // in mm + + TVector3 r_vtx_PV(r_vtx[0] - r_PV[0], r_vtx[1] - r_PV[1], r_vtx[2] - r_PV[2]); + + double pDOTr = p_sum.Dot(r_vtx_PV); + double p_mag = p_sum.Mag(); + double r_mag = r_vtx_PV.Mag(); + + result = pDOTr / (p_mag * r_mag); + return result; +} + +// get track's energy assuming it to be a pion +double get_trackE(edm4hep::TrackState track) { + + double result; + + const double m_pi = 0.13957039; + + TVectorD par = get_trackParam(track); + TVector3 p = ParToP(par); + TLorentzVector p4; + p4.SetXYZM(p[0], p[1], p[2], m_pi); + + result = p4.E(); + return result; +} + +//////////////////////////////////////////////// + +// no of reconstructed V0s +int get_n_SV(FCCAnalysesV0 SV) { + int result = SV.vtx.size(); + return result; +} + +// vector of position of all reconstructed V0 (in mm) +ROOT::VecOps::RVec get_position_SV(FCCAnalysesV0 SV) { + ROOT::VecOps::RVec result; + for (FCCAnalysesVertex ivtx : SV.vtx) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + result.push_back(xyz); + } + return result; +} + +// vector of PDG IDs of all reconstructed V0 +ROOT::VecOps::RVec get_pdg_V0(FCCAnalysesV0 V0) { + ROOT::VecOps::RVec result = V0.pdgAbs; + return result; +} + +// vector of invariant masses of all reconstructed V0 +ROOT::VecOps::RVec get_invM_V0(FCCAnalysesV0 V0) { + ROOT::VecOps::RVec result = V0.invM; + return result; +} + +// +// vector of momenta of all reconstructed V0 +ROOT::VecOps::RVec get_p_SV(FCCAnalysesV0 SV) { + ROOT::VecOps::RVec result; + + for (FCCAnalysesVertex ivtx : SV.vtx) { + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + result.push_back(p_sum); + } + return result; +} + +// vector of chi2 of all reconstructed V0s +ROOT::VecOps::RVec get_chi2_SV(FCCAnalysesV0 SV) { + ROOT::VecOps::RVec result; + + for (FCCAnalysesVertex ivtx : SV.vtx) { + int nDOF = 2 * ivtx.ntracks - 3; + result.push_back(nDOF * ivtx.vertex.chi2); + } + return result; +} + +// passing a vector of FCCAnalysesVertex instead of new structs + +// no of reconstructed SVs +int get_n_SV(ROOT::VecOps::RVec vertices) { + int result = vertices.size(); + return result; +} + +// vector of momenta of all reconstructed vertices (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_p_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + result.push_back(p_sum); + } + return result; +} + +// vector of position of all reconstructed SV (in mm) +ROOT::VecOps::RVec +get_position_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + for (FCCAnalysesVertex ivtx : vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + result.push_back(xyz); + } + return result; +} + +// vector of momentum magnitude of all reconstructed vertices (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_pMag_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + result.push_back(p_sum.Mag()); + } + return result; +} + +// vector of chi2 of all reconstructed vertices (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_chi2_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + int nDOF = 2 * ivtx.ntracks - 3; + result.push_back(nDOF * ivtx.vertex.chi2); + } + return result; +} + +// vector of chi2 (normalised) of all reconstructed vertices (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_norm_chi2_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) + result.push_back(ivtx.vertex.chi2); + return result; +} + +// vector of nDOF of all reconstructed vertices (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_nDOF_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) + result.push_back(2 * ivtx.ntracks - 3); + return result; +} + +// vector of polar angle (theta) of all reconstructed vertices (SV.vtx or +// V0.vtx) +ROOT::VecOps::RVec +get_theta_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + result.push_back(xyz.Theta()); + } + return result; +} + +// vector of azimuth angle (phi) of all reconstructed vertices (SV.vtx or +// V0.vtx) +ROOT::VecOps::RVec +get_phi_SV(ROOT::VecOps::RVec vertices) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + result.push_back(xyz.Phi()); + } + return result; +} + +// vector of (cos of) angles b/n vtx momenta & PV to vtx displacement vectors +ROOT::VecOps::RVec +get_pointingangle_SV(ROOT::VecOps::RVec vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + double iresult = 0.; + + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + edm4hep::Vector3f r_vtx = ivtx.vertex.position; // in mm + edm4hep::Vector3f r_PV = PV.vertex.position; // in mm + + TVector3 r_vtx_PV(r_vtx[0] - r_vtx[0], r_vtx[1] - r_PV[1], + r_vtx[2] - r_PV[2]); + + double pDOTr = p_sum.Dot(r_vtx_PV); + double p_mag = p_sum.Mag(); + double r_mag = r_vtx_PV.Mag(); + + iresult = pDOTr / (p_mag * r_mag); + result.push_back(iresult); + } + return result; +} + +// vector of distances of all reconstructed SV from PV (in mm in xy plane) +ROOT::VecOps::RVec +get_dxy_SV(ROOT::VecOps::RVec vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec result; + TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], + PV.vertex.position[2]); + for (auto &ivtx : vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 x_vtx_PV = x_vtx - x_PV; + + result.push_back(x_vtx_PV.Perp()); + } + return result; +} + +// vector of distances of all reconstructed SV from PV (in mm in 3D) +ROOT::VecOps::RVec +get_d3d_SV(ROOT::VecOps::RVec vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec result; + TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], + PV.vertex.position[2]); + for (auto &ivtx : vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 x_vtx_PV = x_vtx - x_PV; + + result.push_back(x_vtx_PV.Mag()); + } + return result; +} + +// vector of distances of all reconstructed SV from given TVector3d (in mm in +// 3D) +ROOT::VecOps::RVec +get_d3d_SV_obj(ROOT::VecOps::RVec vertices, + TVector3 object) { + ROOT::VecOps::RVec result; + for (auto &ivtx : vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + x_vtx = x_vtx - object; + + result.push_back(x_vtx.Mag()); + } + return result; +} + +// vector of distances of all reconstructed SV from given edm4hep::Vector3d (in +// mm in 3D) +ROOT::VecOps::RVec +get_d3d_SV_obj(ROOT::VecOps::RVec vertices, + edm4hep::Vector3d object) { + ROOT::VecOps::RVec result; + for (auto &ivtx : vertices) { + double dx = ivtx.vertex.position[0] - object.x; + double dy = ivtx.vertex.position[1] - object.y; + double dz = ivtx.vertex.position[2] - object.z; + + TVector3 d3d(dx, dy, dz); + + result.push_back(d3d.Mag()); + } + return result; +} + +// vector of decay position distances of all reconstructed SV from given +// TVector3d (in mm in 3D) +ROOT::VecOps::RVec +get_dR_SV_obj(ROOT::VecOps::RVec vertices, TVector3 object) { + ROOT::VecOps::RVec result; + for (auto &ivtx : vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + result.push_back(1e3 * TMath::Sqrt(pow(ivtx.vertex.position[0], 2) + + pow(ivtx.vertex.position[1], 2) + + pow(ivtx.vertex.position[2], 2)) - + 1e3 * TMath::Sqrt(pow(object.x(), 2) + pow(object.y(), 2) + + pow(object.z(), 2))); + } + return result; +} + +// vector of decay position distances of all reconstructed SV from given +// edm4hep::Vector3d (in mm in 3D) +ROOT::VecOps::RVec +get_dR_SV_obj(ROOT::VecOps::RVec vertices, + edm4hep::Vector3d object) { + ROOT::VecOps::RVec result; + for (auto &ivtx : vertices) { + result.push_back(1e3 * TMath::Sqrt(pow(ivtx.vertex.position[0], 2) + + pow(ivtx.vertex.position[1], 2) + + pow(ivtx.vertex.position[2], 2)) - + 1e3 * TMath::Sqrt(pow(object.x, 2) + pow(object.y, 2) + + pow(object.z, 2))); + } + return result; +} + +// vector of polar angle (theta) of all reconstructed vertices wrt jet axis +// (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_relTheta_SV(ROOT::VecOps::RVec vertices, + ROOT::VecOps::RVec nSV_jet, + ROOT::VecOps::RVec jets) { + ROOT::VecOps::RVec result; + + unsigned int j = 0; + int nSV = nSV_jet[0]; + for (unsigned int i = 0; i < vertices.size(); i++) { + auto &ivtx = vertices[i]; + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + + if (i >= nSV) { + j++; + nSV += nSV_jet[j]; + } + auto &ijet = jets[j]; + double jetTheta = ijet.theta(); + + result.push_back(xyz.Theta() - jetTheta); + } + return result; +} + +// vector of azimuthal angle (phi) of all reconstructed vertices wrt jet axis +// (SV.vtx or V0.vtx) +ROOT::VecOps::RVec +get_relPhi_SV(ROOT::VecOps::RVec vertices, + ROOT::VecOps::RVec nSV_jet, + ROOT::VecOps::RVec jets) { + ROOT::VecOps::RVec result; + + unsigned int j = 0; + int nSV = nSV_jet[0]; + for (unsigned int i = 0; i < vertices.size(); i++) { + auto &ivtx = vertices[i]; + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + + if (i >= nSV) { + j++; + nSV += nSV_jet[j]; + } + auto &ijet = jets[j]; + TVector3 jetP(ijet.px(), ijet.py(), ijet.pz()); + + result.push_back(xyz.DeltaPhi(jetP)); + } + return result; +} + +// For get_SV_jets outputs + +// no of reconstructed SVs +int get_n_SV( + ROOT::VecOps::RVec> vertices) { + int result = 0; + if (vertices.size() != 0) { + for (auto SV_jets : vertices) + result += SV_jets.size(); + } + return result; +} + +// Return the number of reconstructed SVs +ROOT::VecOps::RVec get_n_SV_jets( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec result; + if (vertices.size() != 0) { + for (auto SV_jets : vertices) + result.push_back(SV_jets.size()); + } + return result; +} + +// +// separate V0s by jets +ROOT::VecOps::RVec> +get_svInJets(ROOT::VecOps::RVec vertices, + ROOT::VecOps::RVec nSV_jet) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + int index = 0; + for (unsigned int i : nSV_jet) { + for (unsigned int j = 0; j < i; j++) { + i_result.push_back(vertices[j + index]); + } + + result.push_back(i_result); + i_result.clear(); + index += i; + } + return result; +} + +// separate tracks by jet +std::vector> get_tracksInJets( + ROOT::VecOps::RVec recoparticles, + ROOT::VecOps::RVec thetracks, + ROOT::VecOps::RVec jets, + std::vector> jet_consti) { + std::vector> result; + std::vector iJet_tracks; + + int nJet = jets.size(); + // + for (unsigned int j = 0; j < nJet; j++) { + + std::vector i_jetconsti = jet_consti[j]; + + for (unsigned int ip : i_jetconsti) { + auto &p = recoparticles[ip]; + if (p.tracks_begin >= 0 && p.tracks_begin < thetracks.size()) + iJet_tracks.push_back(thetracks.at(p.tracks_begin)); + } + + result.push_back(iJet_tracks); + iJet_tracks.clear(); + } + return result; +} + +// vector of polar angle (theta) of reconstructed vertices of a jet wrt that jet +// axis [only for vertices from 1 jet] +ROOT::VecOps::RVec +get_relTheta_SV(ROOT::VecOps::RVec vertices, + fastjet::PseudoJet jet) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + // + result.push_back(xyz.Theta() - jet.theta()); + } + // + return result; +} + +// vector of azimuthal angle (phi) of all reconstructed vertices wrt jet axis +// [only for vertices from 1 jet] +ROOT::VecOps::RVec +get_relPhi_SV(ROOT::VecOps::RVec vertices, + fastjet::PseudoJet jet) { + ROOT::VecOps::RVec result; + + for (auto &ivtx : vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 jetP(jet.px(), jet.py(), jet.pz()); + // + result.push_back(xyz.DeltaPhi(jetP)); + } + // + return result; +} + +/////// vec of vec functions (for get_SV_jets) ///////// + +// SV invariant mass +ROOT::VecOps::RVec> +get_invM(ROOT::VecOps::RVec> vertices) { + + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + + for (auto &vertex : i_vertices) { + ROOT::VecOps::RVec p_tracks = + vertex.updated_track_momentum_at_vertex; + // + TLorentzVector p4_vtx; + const double m = 0.13957039; // pion mass + // + for (TVector3 p_tr : p_tracks) { + TLorentzVector p4_tr; + p4_tr.SetXYZM(p_tr.X(), p_tr.Y(), p_tr.Z(), m); + p4_vtx += p4_tr; + } + i_result.push_back(p4_vtx.M()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV momentum +ROOT::VecOps::RVec> +get_p_SV(ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + i_result.push_back(p_sum); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV momentum magnitude +ROOT::VecOps::RVec> get_pMag_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + i_result.push_back(p_sum.Mag()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV daughters multiplicity +ROOT::VecOps::RVec> get_VertexNtrk( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &TheVertex : i_vertices) { + i_result.push_back(TheVertex.ntracks); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV chi2 +ROOT::VecOps::RVec> get_chi2_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + int nDOF = 2 * ivtx.ntracks - 3; + i_result.push_back(nDOF * ivtx.vertex.chi2); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV normalised chi2 +ROOT::VecOps::RVec> get_norm_chi2_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) + i_result.push_back(ivtx.vertex.chi2); + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV no of DOF +ROOT::VecOps::RVec> get_nDOF_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) + i_result.push_back(2 * ivtx.ntracks - 3); + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV theta +ROOT::VecOps::RVec> get_theta_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + + for (auto &ivtx : i_vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + i_result.push_back(xyz.Theta()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV phi +ROOT::VecOps::RVec> +get_phi_SV(ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + + for (auto &ivtx : i_vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + i_result.push_back(xyz.Phi()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV relative theta +ROOT::VecOps::RVec> get_relTheta_SV( + ROOT::VecOps::RVec> vertices, + ROOT::VecOps::RVec jets) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < jets.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + fastjet::PseudoJet i_jet = jets.at(i); + for (auto &ivtx : i_vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + // + i_result.push_back(xyz.Theta() - i_jet.theta()); + } + result.push_back(i_result); + i_result.clear(); + } + // + return result; +} + +// SV relative phi +ROOT::VecOps::RVec> get_relPhi_SV( + ROOT::VecOps::RVec> vertices, + ROOT::VecOps::RVec jets) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < jets.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + fastjet::PseudoJet i_jet = jets.at(i); + for (auto &ivtx : i_vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 jetP(i_jet.px(), i_jet.py(), i_jet.pz()); + // + i_result.push_back(xyz.DeltaPhi(jetP)); + } + result.push_back(i_result); + i_result.clear(); + } + // + return result; +} + +// SV pointing angle wrt PV +ROOT::VecOps::RVec> get_pointingangle_SV( + ROOT::VecOps::RVec> vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + edm4hep::Vector3f r_PV = PV.vertex.position; // in mm + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + for (auto &ivtx : i_vertices) { + double pointangle = 0.; + + ROOT::VecOps::RVec p_tracks = + ivtx.updated_track_momentum_at_vertex; + TVector3 p_sum; + for (TVector3 p_tr : p_tracks) + p_sum += p_tr; + + edm4hep::Vector3f r_vtx = ivtx.vertex.position; // in mm + + TVector3 r_vtx_PV(r_vtx[0] - r_vtx[0], r_vtx[1] - r_PV[1], + r_vtx[2] - r_PV[2]); + + double pDOTr = p_sum.Dot(r_vtx_PV); + double p_mag = p_sum.Mag(); + double r_mag = r_vtx_PV.Mag(); + + pointangle = pDOTr / (p_mag * r_mag); + i_result.push_back(pointangle); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV distance from PV in xy +ROOT::VecOps::RVec> +get_dxy_SV(ROOT::VecOps::RVec> vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], + PV.vertex.position[2]); + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 x_vtx_PV = x_vtx - x_PV; + + i_result.push_back(x_vtx_PV.Perp()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV distance from PV in 3D +ROOT::VecOps::RVec> +get_d3d_SV(ROOT::VecOps::RVec> vertices, + FCCAnalysesVertex PV) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + TVector3 x_PV(PV.vertex.position[0], PV.vertex.position[1], + PV.vertex.position[2]); + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + TVector3 x_vtx(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + TVector3 x_vtx_PV = x_vtx - x_PV; + + i_result.push_back(x_vtx_PV.Mag()); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// SV position in 3D +ROOT::VecOps::RVec> get_position_SV( + ROOT::VecOps::RVec> vertices) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + for (unsigned int i = 0; i < vertices.size(); i++) { + ROOT::VecOps::RVec i_result; + ROOT::VecOps::RVec i_vertices = vertices.at(i); + // + for (auto &ivtx : i_vertices) { + TVector3 xyz(ivtx.vertex.position[0], ivtx.vertex.position[1], + ivtx.vertex.position[2]); + i_result.push_back(xyz); + } + result.push_back(i_result); + i_result.clear(); + } + return result; +} + +// V0 pdg +ROOT::VecOps::RVec> +get_pdg_V0(ROOT::VecOps::RVec pdg, ROOT::VecOps::RVec nSV_jet) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + int index = 0; + for (unsigned int i : nSV_jet) { + for (unsigned int j = 0; j < i; j++) { + i_result.push_back(pdg[j + index]); + } + + result.push_back(i_result); + i_result.clear(); + index += i; + } + return result; +} + +// V0 invariant mass +ROOT::VecOps::RVec> +get_invM_V0(ROOT::VecOps::RVec invM, ROOT::VecOps::RVec nSV_jet) { + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec i_result; + + int index = 0; + for (unsigned int i : nSV_jet) { + for (unsigned int j = 0; j < i; j++) { + i_result.push_back(invM[j + index]); + } + + result.push_back(i_result); + i_result.clear(); + index += i; + } + return result; +} + +} // namespace VertexingUtils + +} // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/myUtils.cc b/analyzers/dataframe/src/myUtils.cc index 4986bfd54a..fbc4a56541 100644 --- a/analyzers/dataframe/src/myUtils.cc +++ b/analyzers/dataframe/src/myUtils.cc @@ -1323,7 +1323,7 @@ ROOT::VecOps::RVec get_truetrack(ROOT::VecOps::RVec in TVector3 momentum ( tlv.Px(),tlv.Py(),tlv.Pz()); - TVectorD track_param = VertexFitterSimple::XPtoPar( vertexFB, momentum, charge ); + TVectorD track_param = VertexingUtils::XPtoPar( vertexFB, momentum, charge ); track.D0 = track_param[0] * 1e3 ; // from meters to mm @@ -1375,7 +1375,7 @@ ROOT::VecOps::RVec get_pseudotrack(ROOT::VecOps::RVec getFCCAnalysesComposite_track(ROOT::VecO vertex.at(p.vertex).vertex.position.z * norm); TVector3 momentum ( p.particle.Px(),p.particle.Py(),p.particle.Pz()); - TVectorD track_param = VertexFitterSimple::XPtoPar( vertexFB, momentum, p.charge ); + TVectorD track_param = VertexingUtils::XPtoPar( vertexFB, momentum, p.charge ); track.D0 = track_param[0] * 1e3 ; // from meters to mm diff --git a/bin/fccanalysis b/bin/fccanalysis index 6c859992f8..2d45603d14 100755 --- a/bin/fccanalysis +++ b/bin/fccanalysis @@ -1,37 +1,40 @@ #!/usr/bin/env python3 +import argparse -if __name__ == "__main__": - import argparse - import sys - parser = argparse.ArgumentParser(description='FCCAnalyses parser') - subparsers = parser.add_subparsers(help='types of running modes', dest='command') - parser_init = subparsers.add_parser('init', help="generate a RDataFrame based FCC analysis") - parser_build = subparsers.add_parser('build', help='build and install local analysis') - parser_pin = subparsers.add_parser('pin', help='pin fccanalyses to the current version of Key4hep stack') - parser_run = subparsers.add_parser('run', help="run a RDataFrame based FCC analysis") - parser_run_final = subparsers.add_parser('final', help="run a RDataFrame based FCC analysis final configuration") - parser_run_plots = subparsers.add_parser('plots', help="run a RDataFrame based FCC analysis plot configuration") +def main(): + parser = argparse.ArgumentParser(description='FCCAnalyses parser') + subparsers = parser.add_subparsers(help='types of running modes', dest='command') + parser_init = subparsers.add_parser('init', help="generate a RDataFrame based FCC analysis") + parser_build = subparsers.add_parser('build', help='build and install local analysis') + parser_pin = subparsers.add_parser('pin', help='pin fccanalyses to the current version of Key4hep stack') + parser_run = subparsers.add_parser('run', help="run a RDataFrame based FCC analysis") + parser_run_final = subparsers.add_parser('final', help="run a RDataFrame based FCC analysis final configuration") + parser_run_plots = subparsers.add_parser('plots', help="run a RDataFrame based FCC analysis plot configuration") + + import config.Parsers as fccpars + fccpars.setup_init_parser(parser_init) + fccpars.setup_build_parser(parser_build) + fccpars.setup_pin_parser(parser_pin) + fccpars.setup_run_parser(parser_run) + fccpars.setup_run_parser_final(parser_run_final) + fccpars.setup_run_parser_plots(parser_run_plots) - from config.Parsers import * - setup_init_parser(parser_init) - setup_build_parser(parser_build) - setup_pin_parser(parser_pin) - setup_run_parser(parser_run) - setup_run_parser_final(parser_run_final) - setup_run_parser_plots(parser_run_plots) + args = parser.parse_args() - args = parser.parse_args() + if args.command == 'init': + from config.FCCAnalysisSetup import setup + setup(parser) + elif args.command == 'build': + from config.build_analysis import build_analysis + build_analysis(parser) + elif args.command == 'pin': + from config.pin_analysis import PinAnalysis + PinAnalysis(parser) + else: + from config.FCCAnalysisRun import run + run(parser) - if args.command == 'init': - from config.FCCAnalysisSetup import setup - setup(parser) - elif args.command == 'build': - from config.build_analysis import build_analysis - build_analysis(parser) - elif args.command == 'pin': - from config.pin_analysis import PinAnalysis - PinAnalysis(parser) - else: - from config.FCCAnalysisRun import run - run(parser) + +if __name__ == "__main__": + main() diff --git a/cmake/FindDelphes.cmake b/cmake/FindDelphes.cmake new file mode 100644 index 0000000000..9c578d3368 --- /dev/null +++ b/cmake/FindDelphes.cmake @@ -0,0 +1,67 @@ +set(searchpath + $ENV{DELPHES_DIR} + $ENV{DELPHES_DIR}/external + $ENV{DELPHES_DIR}/lib + $ENV{DELPHES_DIR}/include + $ENV{DELPHES_DIR}/include/TrackCovariance +# $ENV{DELPHES} +# $ENV{DELPHES}/external +# $ENV{DELPHES}/lib +# $ENV{DELPHES}/include + ) + + +find_library(DELPHES_LIBRARY + NAMES Delphes delphes + HINTS ${searchpath} + PATH_SUFFIXES lib) + +find_path(DELPHES_INCLUDE_DIR + # NAMES DelphesClasses.h Delphes.h + NAMES classes/DelphesClasses.h modules/Delphes.h external/ExRootAnalysis + HINTS ${searchpath} + PATH_SUFFIXES include) + +find_path(DELPHES_EXTERNALS_INCLUDE_DIR + # NAMES DelphesClasses.h Delphes.h + NAMES ExRootAnalysis/ExRootConfReader.h + HINTS ${searchpath} + PATH_SUFFIXES include +) + +find_path(DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR + # NAMES DelphesClasses.h Delphes.h + NAMES TrkUtil.h + HINTS ${searchpath} + PATH_SUFFIXES include +) + +# Necessary to run the tests +find_path(DELPHES_BINARY_DIR + NAMES DelphesROOT + HINTS ${DELPHES_INCLUDE_DIR}/../bin +) + +find_path(DELPHES_CARDS_DIR + NAMES delphes_card_IDEA.tcl + HINTS ${searchpath} + PATH_SUFFIXES cards) + +unset(searchpath) + +set(DELPHES_INCLUDE_DIRS ${DELPHES_INCLUDE_DIR} ${DELPHES_EXTERNALS_INCLUDE_DIR}) +set(DELPHES_EXTERNALS_INCLUDE_DIRS ${DELPHES_EXTERNALS_INCLUDE_DIR}) +set(DELPHES_EXTERNALS_TKCOV_INCLUDE_DIRS ${DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR}) +set(DELPHES_LIBRARIES ${DELPHES_LIBRARY}) + +# Delphes does not offer an obvious version indicator, but we need to know +# whether the TrackCovariance module is available or not. So here we simply +# check whether the corresponding header is installed +find_file(DELPHES_TRACK_COV_HEADER modules/TrackCovariance.h PATHS ${DELPHES_INCLUDE_DIRS} NO_DEFAULT_PATHS) + +include(FindPackageHandleStandardArgs) +# handle the QUIETLY and REQUIRED arguments and set DELPHES_FOUND to TRUE +# if all listed variables are TRUE +find_package_handle_standard_args(Delphes DEFAULT_MSG DELPHES_INCLUDE_DIR DELPHES_EXTERNALS_INCLUDE_DIR DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR DELPHES_LIBRARY) + +mark_as_advanced(DELPHES_INCLUDE_DIR DELPHES_EXTERNALS_INCLUDE_DIR DELPHES_EXTERNALS_TKCOV_INCLUDE_DIR DELPHES_LIBRARY DELPHES_BINARY_DIR DELPHES_TRACK_COV_HEADER) diff --git a/config/FCCAnalysisRun.py b/config/FCCAnalysisRun.py index ec64d818d2..ddf47063cb 100644 --- a/config/FCCAnalysisRun.py +++ b/config/FCCAnalysisRun.py @@ -10,13 +10,6 @@ from config.common_defaults import deffccdicts import datetime -print ("----> Load cxx analyzers from libFCCAnalyses... ",) -ROOT.gSystem.Load("libFCCAnalyses") -ROOT.gErrorIgnoreLevel = ROOT.kFatal -#Is this still needed?? 01/04/2022 still to be the case -_fcc = ROOT.dummyLoader - - DATE = datetime.datetime.fromtimestamp(datetime.datetime.now().timestamp()).strftime('%Y-%m-%d_%H-%M-%S') #__________________________________________________________ @@ -433,7 +426,7 @@ def sendToBatch(rdfModule, chunkList, process, analysisFile): subprocess.getstatusoutput('chmod 777 %s'%(frunname)) frun.write('#!/bin/bash\n') - frun.write('source ${LOCAL_DIR}/setup.sh\n') + frun.write('source ' + localDir + '/setup.sh\n') #add userBatchConfig if any if userBatchConfig!="": @@ -448,9 +441,9 @@ def sendToBatch(rdfModule, chunkList, process, analysisFile): frun.write('cd job{}_chunk{}\n'.format(process,ch)) if not os.path.isabs(outputDir): - frun.write('$LOCAL_DIR/bin/fccanalysis run {} --batch --output {}chunk{}.root --files-list '.format(analysisFile, outputDir, ch)) + frun.write(localDir + '/bin/fccanalysis run {} --batch --output {}chunk{}.root --files-list '.format(analysisFile, outputDir, ch)) else: - frun.write('$LOCAL_DIR/bin/fccanalysis run {} --batch --output {}{}/chunk{}.root --files-list '.format(analysisFile, outputDir, process,ch)) + frun.write(localDir + '/bin/fccanalysis run {} --batch --output {}{}/chunk{}.root --files-list '.format(analysisFile, outputDir, process,ch)) for ff in range(len(chunkList[ch])): frun.write(' %s'%(chunkList[ch][ff])) @@ -482,7 +475,7 @@ def sendToBatch(rdfModule, chunkList, process, analysisFile): frun_condor.write('Log = {}/condor_job.{}.$(ClusterId).$(ProcId).log\n'.format(logDir,process)) frun_condor.write('Output = {}/condor_job.{}.$(ClusterId).$(ProcId).out\n'.format(logDir,process)) frun_condor.write('Error = {}/condor_job.{}.$(ClusterId).$(ProcId).error\n'.format(logDir,process)) - frun_condor.write('getenv = True\n') + frun_condor.write('getenv = False\n') frun_condor.write('environment = "LS_SUBCWD={}"\n'.format(logDir)) # not sure frun_condor.write('requirements = ( (OpSysAndVer =?= "CentOS7") && (Machine =!= LastRemoteHost) && (TARGET.has_avx2 =?= True) )\n') frun_condor.write('on_exit_remove = (ExitBySignal == False) && (ExitCode == 0)\n') @@ -500,9 +493,9 @@ def sendToBatch(rdfModule, chunkList, process, analysisFile): #__________________________________________________________ def addeosType(fileName): sfileName=fileName.split('/') - if sfileName[1]=='experiment': + if sfileName[2]=='experiment': fileName='root://eospublic.cern.ch/'+fileName - elif sfileName[1]=='user' or sfileName[1].contains('home-'): + elif sfileName[2]=='user' or sfileName[2].contains('home-'): fileName='root://eosuser.cern.ch/'+fileName else: print('unknown eos type, please check with developers as it might not run with best performances') @@ -517,7 +510,7 @@ def runLocal(rdfModule, fileList, args): nevents_local = 0 for fileName in fileList: - if fileName.split('/')[0]=='eos': + if fileName.split('/')[1]=='eos': fileName=addeosType(fileName) fileListRoot.push_back(fileName) @@ -1025,7 +1018,6 @@ def runFinal(rdfModule): #__________________________________________________________ def runPlots(analysisFile): - import config.doPlots as dp dp.run(analysisFile) @@ -1085,46 +1077,74 @@ def run(mainparser, subparser=None): print("specify a valid analysis script in the command line arguments") sys.exit(3) + print ("----> Info: Loading analyzers from libFCCAnalyses... ",) + ROOT.gSystem.Load("libFCCAnalyses") + ROOT.gErrorIgnoreLevel = ROOT.kFatal + #Is this still needed?? 01/04/2022 still to be the case + _fcc = ROOT.dummyLoader + #set the RDF ELogLevel try: verbosity = ROOT.Experimental.RLogScopedVerbosity(ROOT.Detail.RDF.RDFLogChannel(), getattr(ROOT.Experimental.ELogLevel,args.eloglevel)) except AttributeError: pass #load the analysis - analysisFile=os.path.abspath(analysisFile) - print ("--------------loading analysis file ",analysisFile) + analysisFile = os.path.abspath(analysisFile) + print('----> Info: Loading analysis file:') + print(' ' + analysisFile) rdfSpec = importlib.util.spec_from_file_location("rdfanalysis", analysisFile) rdfModule = importlib.util.module_from_spec(rdfSpec) rdfSpec.loader.exec_module(rdfModule) - try: - print(args.command) - args.command - if args.command == "run": runStages(args, rdfModule, args.preprocess, analysisFile) - elif args.command == "final": runFinal(rdfModule) - elif args.command == "plots": runPlots(analysisFile) + if hasattr(args, 'command'): + if args.command == "run": + try: + runStages(args, rdfModule, args.preprocess, analysisFile) + except Exception as excp: + print('----> Error: During the execution of the stage file:') + print(' ' + analysisFile) + print(' exception occurred:') + print(excp) + elif args.command == "final": + try: + runFinal(rdfModule) + except Exception as excp: + print('----> Error: During the execution of the final stage file:') + print(' ' + analysisFile) + print(' exception occurred:') + print(excp) + elif args.command == "plots": + try: + runPlots(analysisFile) + except Exception as excp: + print('----> Error: During the execution of the plots file:') + print(' ' + analysisFile) + print(' exception occurred:') + print(excp) return - except Exception as e: - print("============running the old way") + print('----> Info: Running the old way...') + print(' This way of running the analysis is deprecated and will') + print(' be removed in the next release!') - #below is legacy using the old way of runnig with options in "python config/FCCAnalysisRun.py analysis.py --options - #check if this is final analysis + # below is legacy using the old way of runnig with options in + # "python config/FCCAnalysisRun.py analysis.py --options check if this is + # final analysis if args.final: if args.plots: - print ('----> Can not have --plots with --final, exit') + print('----> Can not have --plots with --final, exit') sys.exit(3) if args.preprocess: - print ('----> Can not have --preprocess with --final, exit') + print('----> Can not have --preprocess with --final, exit') sys.exit(3) runFinal(rdfModule) elif args.plots: if args.final: - print ('----> Can not have --final with --plots, exit') + print('----> Can not have --final with --plots, exit') sys.exit(3) if args.preprocess: - print ('----> Can not have --preprocess with --plots, exit') + print('----> Can not have --preprocess with --plots, exit') sys.exit(3) runPlots(analysisFile) @@ -1134,10 +1154,10 @@ def run(mainparser, subparser=None): else: if args.preprocess: if args.plots: - print ('----> Can not have --plots with --preprocess, exit') + print('----> Can not have --plots with --preprocess, exit') sys.exit(3) if args.final: - print ('----> Can not have --final with --preprocess, exit') + print('----> Can not have --final with --preprocess, exit') sys.exit(3) runStages(args, rdfModule, args.preprocess, analysisFile) diff --git a/config/pin_analysis.py b/config/pin_analysis.py index d7d6a1ff2a..349e0b89e7 100644 --- a/config/pin_analysis.py +++ b/config/pin_analysis.py @@ -39,9 +39,8 @@ def show_pin(self): Show current pin ''' if not self.pin_path.is_file(): - print('----> Error: Analysis pin file not found!') - print(' Aborting...') - sys.exit(3) + print('----> Info: Analysis not pinned.') + sys.exit(0) with open(self.pin_path, 'r') as pinfile: lines = pinfile.readlines() diff --git a/examples/FCCee/fullSim/caloNtupleizer/analysis.py b/examples/FCCee/fullSim/caloNtupleizer/analysis.py index 7ef3e53cff..9d4e8aa659 100644 --- a/examples/FCCee/fullSim/caloNtupleizer/analysis.py +++ b/examples/FCCee/fullSim/caloNtupleizer/analysis.py @@ -149,7 +149,7 @@ def run(self): dict_outputBranchName_function["highestEnergyCluster_cells_transient"] = "myRange(PositionedCaloClusterCells, highestEnergyCluster_firstCell_index, highestEnergyCluster_lastCell_index)" dict_outputBranchName_function["highestEnergyCluster_cells_energy"] = "CaloNtupleizer::getCaloHit_energy(highestEnergyCluster_cells_transient)" - dict_outputBranchName_function["relative_highestEnergyCluster_cells_phi"] = "CaloNtupleizer::getCaloHit_phi(highestEnergyCluster_cells_transient) - highestEnergyCluster_phi" + dict_outputBranchName_function["highestEnergyCluster_cells_relative_phi"] = "CaloNtupleizer::getCaloHit_phi(highestEnergyCluster_cells_transient) - highestEnergyCluster_phi" dict_outputBranchName_function["highestEnergyCluster_cells_relative_theta"] = "CaloNtupleizer::getCaloHit_theta(highestEnergyCluster_cells_transient) - highestEnergyCluster_theta" dict_outputBranchName_function["highestEnergyCluster_cells_layer"] = "CaloNtupleizer::getCaloHit_layer(highestEnergyCluster_cells_transient)" dict_outputBranchName_function["highestEnergyCluster_cells_n"] = "highestEnergyCluster_cells_transient.size()" @@ -175,7 +175,7 @@ def run(self): weaver = WeaverUtils.setup_weaver(args.weaverFiles[0], args.weaverFiles[1], ('highestEnergyCluster_cells_energy', 'relative_highestEnergyCluster_cells_phi', 'highestEnergyCluster_cells_relative_theta', 'highestEnergyCluster_cells_layer')) df2 = (df2.Define("cells_e", "Utils::as_vector(highestEnergyCluster_cells_energy)") .Define("cells_theta", "Utils::as_vector(highestEnergyCluster_cells_relative_theta)") - .Define("cells_phi", "Utils::as_vector(relative_highestEnergyCluster_cells_phi)") + .Define("cells_phi", "Utils::as_vector(highestEnergyCluster_cells_relative_phi)") .Define("cells_layer", "Utils::as_vector(highestEnergyCluster_cells_layer)") .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_phi, cells_theta, cells_layer)") .Define("highestEnergyCluster_isPhoton_inferred", "WeaverUtils::get_weight(MVAVec, 0)") diff --git a/examples/FCCee/smearing/smear_jets.py b/examples/FCCee/smearing/smear_jets.py index b91db672b9..557e7cfd24 100644 --- a/examples/FCCee/smearing/smear_jets.py +++ b/examples/FCCee/smearing/smear_jets.py @@ -135,14 +135,14 @@ def analysers(df): ## produce smeared collections ### run same sequences but with smeared collection - scale_factors = [0.5, 1.0, 2.0] - collections_ip = deepcopy(collections) - collections_res = deepcopy(collections) + scale_factors = [1.0, 2.0, 5.0, 10.0] - - ## 1. do Impact parameter smearing first + for sf in scale_factors: + + ## 1. do Impact parameter smearing first + collections_ip = deepcopy(collections) ip_tag = "ip{}".format(sf).replace(".", "p") collections_ip["TrackState"] = "TrackState_{}".format(ip_tag) @@ -160,9 +160,9 @@ def analysers(df): df = jet_sequence(df, collections_ip, output_branches, ip_tag) - ## 2. do Neutral Hadron energy smearing - for sf in scale_factors: + ## 2. do Neutral Hadron energy smearing + collections_res = deepcopy(collections) res_tag = "res{}".format(sf).replace(".", "p") collections_res["PFParticles"] = "ReconstructedParticles_{}".format(res_tag) @@ -179,6 +179,49 @@ def analysers(df): ## run full sequence with energy nh smeared df = jet_sequence(df, collections_res, output_branches, res_tag) + ## 3. do dNdx smearing + + collections_dndx = deepcopy(collections) + dndx_tag = "dndx{}".format(sf).replace(".", "p") + collections_dndx["dNdx"] = "dNdx_{}".format(dndx_tag) + + df = df.Define( + collections_dndx["dNdx"], + ROOT.SmearObjects.SmearedTracksdNdx(sf, False), + [ + collections["PFParticles"], + collections["dNdx"], + collections["PathLength"], + "reco_mc_index", + collections["GenParticles"], + ], + ) + + ## run full sequence with energy nh smeared + df = jet_sequence(df, collections_dndx, output_branches, dndx_tag) + + + ## 4. do tof smearing + collections_tof = deepcopy(collections) + tof_tag = "tof{}".format(sf).replace(".", "p") + collections_tof["TrackerHits"] = "tof_{}".format(tof_tag) + + df = df.Define( + collections_tof["TrackerHits"], + ROOT.SmearObjects.SmearedTracksTOF(sf, False), + [ + collections["PFParticles"], + collections["PFTracks"], + collections["TrackerHits"], + collections["PathLength"], + "reco_mc_index", + collections["GenParticles"], + ], + ) + + ## run full sequence with energy nh smeared + df = jet_sequence(df, collections_tof, output_branches, tof_tag) + return df # __________________________________________________________ diff --git a/examples/FCCee/tutorials/vertexing/analysis_Bs2DsK_MCseeded.py b/examples/FCCee/tutorials/vertexing/analysis_Bs2DsK_MCseeded.py new file mode 100644 index 0000000000..075dedce5c --- /dev/null +++ b/examples/FCCee/tutorials/vertexing/analysis_Bs2DsK_MCseeded.py @@ -0,0 +1,355 @@ +#Optional test file +testFile="/eos/experiment/fcc/ee/generation/DelphesEvents/winter2023/IDEA/p8_ee_Zbb_ecm91_EvtGen_Bs2DsK/events_017659734.root" + + + +########################################################################################################## + +# +# Code dedicated to the Bs to Ds K analysis, not in central FCCAnalysis library +# + +########################################################################################################## + + +ReconstructedDs_code=''' + +// +// -- This method reconstructs the Ds pseudo-track using Franco's VertexMore. +// -- It returns a vector by convenience - but the vector only contains at most one TrackState. +// -- tracks in input = the DsTracks +// + +#include "FCCAnalyses/VertexFitterSimple.h" +#include "FCCAnalyses/VertexingUtils.h" + + +ROOT::VecOps::RVec ReconstructedDs( ROOT::VecOps::RVec tracks, + bool AddMassConstraints = false) { + + + ROOT::VecOps::RVec result; + + int Ntr = tracks.size(); + if ( Ntr != 3 ) return result; + + TVectorD** trkPar = new TVectorD*[Ntr]; + TMatrixDSym** trkCov = new TMatrixDSym*[Ntr]; + + bool Units_mm = true; + + for (Int_t i = 0; i < Ntr; i++) { + edm4hep::TrackState t = tracks[i] ; + TVectorD par = FCCAnalyses::VertexingUtils::get_trackParam( t, Units_mm ) ; + trkPar[i] = new TVectorD( par ); + TMatrixDSym Cov = FCCAnalyses::VertexingUtils::get_trackCov( t, Units_mm ); + trkCov[i] = new TMatrixDSym ( Cov ); + } + + VertexFit theVertexFit( Ntr, trkPar, trkCov ); + TVectorD x = theVertexFit.GetVtx() ; // this actually runs the fit + + VertexFit* vertexfit = &theVertexFit; + VertexMore vertexmore( vertexfit, Units_mm ); + + if ( AddMassConstraints ) { + + const double kaon_mass = 4.9367700e-01 ; + const double pion_mass = 0.140; + const double Ds_mass = 1.9683000e+00 ; + const double Phi_mass = 1.0194610e+00 ; + + double Ds_masses[3] = { kaon_mass, kaon_mass, pion_mass }; + int Ds_list[3] = { 0, 1, 2 }; + vertexmore.AddMassConstraint(Ds_mass, 3, Ds_masses, Ds_list); // Ds mass constraint + + /* + // Phi mass constraint: does not crash, but pulls of the resulting Bs vertex are bad. + double Phi_masses[2] = { kaon_mass, kaon_mass }; + int Phi_list[2] = { 0, 1 }; + vertexmore.AddMassConstraint(Phi_mass, 2, Phi_masses, Phi_list); // Phi mass constraint + */ + + vertexmore.MassConstrFit(); + } + + TVectorD Ds_track_param = vertexmore.GetVpar(); + TMatrixDSym cov = vertexmore.GetVcov(); + + TVectorD Ds_track_param_edm4hep = FCCAnalyses::VertexingUtils::Delphes2Edm4hep_TrackParam( Ds_track_param, Units_mm ); + edm4hep::TrackState track; + track.D0 = Ds_track_param_edm4hep[0] ; + track.phi = Ds_track_param_edm4hep[1]; + track.omega = Ds_track_param_edm4hep[2]; + track.Z0 = Ds_track_param_edm4hep[3] ; + track.tanLambda = Ds_track_param_edm4hep[4] ; + + // now the covariance matrix - lower-triangle : + + TMatrixDSym covM(5); + std::array covMatrix = FCCAnalyses::VertexingUtils::Delphes2Edm4hep_TrackCovMatrix( cov, Units_mm ) ; + + track.covMatrix = covMatrix ; + + result.push_back( track ); + + return result; +} +''' + + +########################################################################################################## + +Momentum_ReconstructedDs_code=''' + +#include "FCCAnalyses/VertexFitterSimple.h" +#include "FCCAnalyses/VertexingUtils.h" + +TVector3 Momentum_ReconstructedDs( edm4hep::TrackState Ds_pseudoTrack ) { + + TVectorD Param = FCCAnalyses::VertexingUtils::get_trackParam( Ds_pseudoTrack ); // track parameters, Franco's convention + TVector3 result = FCCAnalyses::VertexingUtils::ParToP( Param ); + + return result; + +} + +''' + + +########################################################################################################## + + +Tracks_for_the_Bs_vertex_code=''' +#include "FCCAnalyses/VertexingUtils.h" +#include "FCCAnalyses/ReconstructedParticle.h" + +ROOT::VecOps::RVec tracks_for_fitting_the_Bs_vertex( + ROOT::VecOps::RVec ReconstructedDs, + ROOT::VecOps::RVec BachelorKTrack) { + + ROOT::VecOps::RVec result; + if ( ReconstructedDs.size() != 1 ) return result; + if ( BachelorKTrack.size() != 1 ) return result; + + result.push_back( ReconstructedDs[0]) ; // the pseudo-Ds track + result.push_back( BachelorKTrack[0] ); // the bachelor K + + return result; +} +''' + + +########################################################################################################## + + +import ROOT +ROOT.gInterpreter.Declare(ReconstructedDs_code) +ROOT.gInterpreter.Declare(Momentum_ReconstructedDs_code) +ROOT.gInterpreter.Declare(Tracks_for_the_Bs_vertex_code) + + + +#Mandatory: RDFanalysis class where the use defines the operations on the TTree +class RDFanalysis(): + + #__________________________________________________________ + #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 + def analysers(df): + df2 = ( + df + + .Alias("Particle1", "Particle#1.index") + .Alias("MCRecoAssociations0", "MCRecoAssociations#0.index") + .Alias("MCRecoAssociations1", "MCRecoAssociations#1.index") + + # --------------------------------------------------------------------------------------- + # + # ----- Retrieve the indices of the MC Particles of interest + + # MC indices of the decay Bs_bar (PDG = -531) -> Ds+ (PDG = 431) K- (PDG = -321) + # Retrieves a vector of int's which correspond to indices in the Particle block + # vector[0] = the mother, and then the daughters in the order specified, i.e. here + # [1] = the Ds+, [2] = the K- + # Boolean arguments : + # 1st: stableDaughters. when set to true, the daughters specified in the list are looked + # for among the final, stable particles that come out from the mother, i.e. the decay tree is + # explored recursively if needed. + # 2nd: chargeConjugateMother + # 3rd: chargeConjugateDaughters + # 4th: inclusiveDecay: when set to false, if a mother is found, that decays + # into the particles specified in the list plus other particle(s), this decay is not selected. + # If the event contains more than one such decays,only the first one is kept. + .Define("Bs2DsK_indices", "MCParticle::get_indices( -531, {431, -321}, false, true, true, false) ( Particle, Particle1)" ) + + # select events for which the requested decay chain has been found: + .Filter("Bs2DsK_indices.size() > 0") + + .Define("Bs_MCindex", "return Bs2DsK_indices[0] ;" ) + .Define("Ds_MCindex", "return Bs2DsK_indices[1] ;" ) + .Define("BachelorK_MCindex", "return Bs2DsK_indices[2] ;" ) + + + # MC indices of (this) Ds+ -> K+ K- Pi+ + # Boolean arguments : + # 1st: stableDaughters. when set to true, the daughters specified in the list are looked + # for among the final, stable particles that come out from the mother, i.e. the decay tree is + # explored recursively if needed. + # 2nd: chargeConjugateDaughters + # 3rd: inclusiveDecay + .Define("Ds2KKPi_indices", "MCParticle::get_indices_MotherByIndex( Ds_MCindex, { 321, -321, 211 }, true, true, false, Particle, Particle1 )") + + .Define("Kplus_MCindex", "return Ds2KKPi_indices[1] ; ") + .Define("Kminus_MCindex", "return Ds2KKPi_indices[2] ; ") + .Define("Piplus_MCindex", "return Ds2KKPi_indices[3] ; ") + + # --------------------------------------------------------------------------------------- + + + # --------------------------------------------------------------------------------------- + # ----- The MC Particles : + + # the MC Bs : + .Define("Bs", "return Particle[ Bs_MCindex ]; ") + # the MC Ds : + .Define("Ds", "return Particle[ Ds_MCindex ]; ") + # the MC bachelor K- from the Bs decay : + .Define("BachelorK", "return Particle[ BachelorK_MCindex ]; " ) + # The MC legs from the Ds decay + .Define("Kplus", "return Particle[ Kplus_MCindex ] ; ") + .Define("Kminus", "return Particle[ Kminus_MCindex ] ; ") + .Define("Piplus", "return Particle[ Piplus_MCindex ] ; ") + + # some MC-truth kinematic quantities: + .Define("BachelorK_px", "return BachelorK.momentum.x ;") + .Define("BachelorK_py", "return BachelorK.momentum.y ;") + .Define("BachelorK_pz", "return BachelorK.momentum.z ;") + + .Define("Ds_px", "return Ds.momentum.x ;") + .Define("Ds_py", "return Ds.momentum.y ;") + .Define("Ds_pz", "return Ds.momentum.z ;") + .Define("Ds_pt", "ROOT::VecOps::RVec v; v.push_back( Ds ); return MCParticle::get_pt(v ) ;") + + + # --------------------------------------------------------------------------------------- + + + # --------------------------------------------------------------------------------------- + # ----- The MC-truth decay vertices of the Bs and of the Ds + + # Note: in case the Bs or Bsbar has oscillated before it decays, the vertex returned + # below is the decay vertex of the Bs after oscillation. + + # MC Decay vertex of the Bs = the production vertex of the Bachelor K + .Define("BsMCDecayVertex", "return BachelorK.vertex; " ) + # MC Decay vertex of the Ds = the production vertex of the Kplus + .Define("DsMCDecayVertex", "return Kplus.vertex ; " ) + + # --------------------------------------------------------------------------------------- + + + # --------------------------------------------------------------------------------------- + # ----- The RecoParticles that are MC-matched with the particles of the Ds decay + + # RecoParticles associated with the Ds decay + # Note: the size of DsRecoParticles below is always 3 provided that Ds2KKPi_indices is not empty. + # possibly including "dummy" particles in case one of the legs did not make a RecoParticle + # (e.g. because it is outside the tracker acceptance). + # This is done on purpose, in order to maintain the mapping with the indices - i.e. the 1st particle in + # the list BsRecoParticles is the Kminus, then the Kplus, then the Piplus. + # (selRP_matched_to_list ignores the unstable MC particles that are in the input list of indices + # hence the mother particle, which is the [0] element of the Ds2KKPi_indices vector). + # + # The matching between RecoParticles and MCParticles requires 4 collections. For more + # detail, see https://github.com/HEP-FCC/FCCAnalyses/tree/master/examples/basics + + .Define("DsRecoParticles", " ReconstructedParticle2MC::selRP_matched_to_list( Ds2KKPi_indices, MCRecoAssociations0,MCRecoAssociations1,ReconstructedParticles,Particle)") + + # the corresponding tracks - here, dummy particles, if any, are removed + .Define("DsTracks", "ReconstructedParticle2Track::getRP2TRK( DsRecoParticles, EFlowTrack_1)" ) + + # number of tracks used to reconstruct the Ds vertex + .Define("n_DsTracks", "ReconstructedParticle2Track::getTK_n( DsTracks )") + + # --------------------------------------------------------------------------------------- + # ------ Reco'ed vertex of the Ds + + # The index "3" below is just to indicate that this is a "tertiary" vertex, it has no influence + # on the vertex fitting. + # Note: no mass constraint is applied here. See the code of "ReconstructedDs" for an example + # of how mass constraints can be applied. + .Define("DsVertexObject", "VertexFitterSimple::VertexFitter_Tk( 3, DsTracks)" ) + .Define("DsVertex", "VertexingUtils::get_VertexData( DsVertexObject )") + + # ------------------------------------------------------------------------------------------------------- + + + # ------------------------------------------------------------------------------------------------------- + # ---------- Reconstruction of the Bs vertex + + # The Ds pseudoTrack (TrackState) - returned as a vector of TrackState, with one single element : + # boolean = true: add Ds mass constraint in the Ds vertex fit. + # Default = false (no mass constraint applied) + .Define("v_Ds_PseudoTrack", "ReconstructedDs( DsTracks, true )") + + # Number of Ds pseudoTracks. In principle it should always be equal to one at this stage. + .Define("n_pdt", "ReconstructedParticle2Track::getTK_n( v_Ds_PseudoTrack )") + # but it is explicitely required to be non zero, otherwise the code that comes next would crash: + .Filter("n_pdt > 0") + + # The momentum vector (TVector3) of the Ds : + .Define("Ds_momentum", "Momentum_ReconstructedDs( v_Ds_PseudoTrack[0] ) ") + .Define("RecoDs_px", "return Ds_momentum.x() ") + .Define("RecoDs_py", "return Ds_momentum.y() ") + .Define("RecoDs_pz", "return Ds_momentum.z() ") + + + # the RecoParticle associated with the bachelor K + .Define("BsRecoParticles", "ReconstructedParticle2MC::selRP_matched_to_list( Bs2DsK_indices, MCRecoAssociations0,MCRecoAssociations1,ReconstructedParticles,Particle)") + .Define("RecoBachelorK", " return BsRecoParticles[0] ; ") # only the bachelor K is stable, among the indices of Bs2DsK_indices + + .Define("RecoBachelorK_px", "return RecoBachelorK.momentum.x; ") + .Define("RecoBachelorK_py", "return RecoBachelorK.momentum.y; ") + .Define("RecoBachelorK_pz", "return RecoBachelorK.momentum.z; ") + + + # and the corresponding track + .Define("v_RecoBachelorK", "ROOT::VecOps::RVec v; v.push_back( RecoBachelorK ); return v; ") + .Define("v_BachelorKTrack", "ReconstructedParticle2Track::getRP2TRK( v_RecoBachelorK, EFlowTrack_1)" ) + + # Now we have the two tracks that we need for the Bs vertex : + .Define("BsTracks", "tracks_for_fitting_the_Bs_vertex( v_Ds_PseudoTrack, v_BachelorKTrack) ") + + + # --------------------------------------------------------------------------------------- + # ------ Reco'ed vertex of the Bs + + .Define("BsVertexObject", "VertexFitterSimple::VertexFitter_Tk( 2, BsTracks )" ) + .Define("n_BsTracks", "ReconstructedParticle2Track::getTK_n( BsTracks )") + + # This is the final Bs vertex + .Define("BsVertex", "VertexingUtils::get_VertexData( BsVertexObject )") + + + ) + return df2 + + + #__________________________________________________________ + #Mandatory: output function, please make sure you return the branchlist as a python list + def output(): + branchList = [ + "DsMCDecayVertex", + "BsMCDecayVertex", + "n_DsTracks", + "DsVertex", + "Ds_momentum", + "n_BsTracks", + "BsVertex", + "BachelorK_px", "BachelorK_py", "BachelorK_pz", + "RecoBachelorK_px","RecoBachelorK_py","RecoBachelorK_pz", + "Ds_px", "Ds_py", "Ds_pz","Ds_pt", + "RecoDs_px", "RecoDs_py", "RecoDs_pz", + + ] + return branchList diff --git a/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py b/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py index 3472492e2d..e95fcd58d3 100644 --- a/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py +++ b/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py @@ -40,7 +40,7 @@ def analysers(df): # This is not a good estimate of the primary vertex: even in a Z -> uds event, there are displaced tracks (e.g. Ks, Lambdas), which would bias the fit. # Below, we determine the "primary tracks" using an iterative algorithm - cf LCFI+. - .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( VertexObject_allTracks, EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0., 0)") + .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)") # Now we run again the vertex fit, but only on the primary tracks : .Define("PrimaryVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300) ") diff --git a/examples/FCCee/tutorials/vertexing/plots_Bs2DsK.x b/examples/FCCee/tutorials/vertexing/plots_Bs2DsK.x new file mode 100644 index 0000000000..23ec54c390 --- /dev/null +++ b/examples/FCCee/tutorials/vertexing/plots_Bs2DsK.x @@ -0,0 +1,55 @@ +{ + +// ------------------------------------------------------------------ + +// Ds vertex fit + +TString cut = "n_DsTracks==3"; + +// normalised chi2 of the fit : (Ndf = 2 x Ntracks - 3 ) +events->Draw("DsVertex.chi2", cut); + +// resolutions in x : +events->Draw("(DsVertex.position.x-DsMCDecayVertex.x)",cut) ; + +// pulls of the fitted vertex position : +events->Draw("(DsVertex.position.x-DsMCDecayVertex.x)/TMath::Sqrt( DsVertex.covMatrix[0] )",cut); +events->Draw("(DsVertex.position.y-DsMCDecayVertex.y)/TMath::Sqrt( DsVertex.covMatrix[2] )",cut) ; +events->Draw("(DsVertex.position.z-DsMCDecayVertex.z)/TMath::Sqrt( DsVertex.covMatrix[5] )",cut) ; + + +// ------------------------------------------------------------------ + + +// Bs vertex fit + +cut = "n_DsTracks==3 && n_BsTracks ==2" ; + + +// normalised chi2 : +events->Draw("BsVertex.chi2", cut); + +// resolutions: +events->Draw("(BsVertex.position.x-BsMCDecayVertex.x)",cut) ; + +// pulls in x : +events->Draw("(BsVertex.position.x-BsMCDecayVertex.x)/TMath::Sqrt( BsVertex.covMatrix[0] )",cut); + +// pulls of the flight distance : + +TString fld_mm = "TMath::Sqrt( pow( BsVertex.position.x, 2) + pow( BsVertex.position.y,2) + pow( BsVertex.position.z,2))"; +TString fld_gen_mm = "TMath::Sqrt( pow( BsMCDecayVertex.x[0], 2) + pow( BsMCDecayVertex.y[0],2) + pow( BsMCDecayVertex.z[0],2) )"; +TString fld_res_mm = fld_mm + " - " + fld_gen_mm; +TString term1 = " BsVertex.position.x * ( BsVertex.covMatrix[0] * BsVertex.position.x + BsVertex.covMatrix[1] * BsVertex.position.y + BsVertex.covMatrix[3] * BsVertex.position.z ) " ; +TString term2 = " BsVertex.position.y * ( BsVertex.covMatrix[1] * BsVertex.position.x + BsVertex.covMatrix[2] * BsVertex.position.y + BsVertex.covMatrix[4] * BsVertex.position.z ) " ; +TString term3 = " BsVertex.position.z * ( BsVertex.covMatrix[3] * BsVertex.position.x + BsVertex.covMatrix[4] * BsVertex.position.y + BsVertex.covMatrix[5] * BsVertex.position.z ) "; +TString tsum = term1 + " + " + term2 + " + " + term3; +TString fld_unc = " ( TMath::Sqrt( " + tsum + ") / " + fld_mm +" ) "; +TString fld_pull = "( " + fld_res_mm + " ) / " + fld_unc; +events->Draw(fld_pull , cut); + + + + +} + diff --git a/examples/FCCee/vertex/analysis.py b/examples/FCCee/vertex/analysis.py index 0efa96d9ae..1a3ebfcfce 100644 --- a/examples/FCCee/vertex/analysis.py +++ b/examples/FCCee/vertex/analysis.py @@ -84,9 +84,10 @@ def run(self): # --- now, determime the primary (and secondary) tracks without using the MC-matching: # First, reconstruct a vertex from all tracks - .Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300)") + #.Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300)") # Select the tracks that are reconstructed as primaries - .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( VertexObject_allTracks, EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0., 0)") + .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)") + .Define("n_RecoedPrimaryTracks", "ReconstructedParticle2Track::getTK_n( RecoedPrimaryTracks )") # the final primary vertex : .Define("FinalVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300) ") diff --git a/examples/FCCee/vertex_lcfiplus/analysis_SV.py b/examples/FCCee/vertex_lcfiplus/analysis_SV.py index d099827230..cf7a7b7930 100644 --- a/examples/FCCee/vertex_lcfiplus/analysis_SV.py +++ b/examples/FCCee/vertex_lcfiplus/analysis_SV.py @@ -1,3 +1,5 @@ +testFile="root://eospublic.cern.ch//eos/experiment/fcc/ee/generation/DelphesEvents/winter2023/IDEA/p8_ee_Zbb_ecm91/events_066726720.root" + #Mandatory: List of processes processList = { 'p8_ee_Zuds_ecm91':{'fraction':0.0001}, #Run 0.01% statistics in one output file named /p8_ee_Zuds_ecm91.root @@ -37,10 +39,9 @@ def analysers(df): ##### # determime the primary (and secondary) tracks without using the MC-matching: - # First, reconstruct a vertex from all tracks - .Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300)") # Select the tracks that are reconstructed as primaries - .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( VertexObject_allTracks, EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0., 0)") + .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)") + .Define("n_RecoedPrimaryTracks", "ReconstructedParticle2Track::getTK_n( RecoedPrimaryTracks )") # the final primary vertex : .Define("PrimaryVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300) ") @@ -72,7 +73,7 @@ def analysers(df): .Define("SV_jet", "VertexFinderLCFIPlus::get_SV_jets(ReconstructedParticles, EFlowTrack_1, PrimaryVertexObject, IsPrimary_based_on_reco, jets_ee_kt, jetconstituents_ee_kt)") # finding SVs in the event (two interfaces) #.Define("SV_evt1", "VertexFinderLCFIPlus::get_SV_event(ReconstructedParticles, EFlowTrack_1, PrimaryVertexObject, IsPrimary_based_on_reco)") - #.Define("SV_evt2", "VertexFinderLCFIPlus::get_SV_event(SecondaryTracks, PrimaryVertexObject)") + #.Define("SV_evt2", "VertexFinderLCFIPlus::get_SV_event(SecondaryTracks, EFlowTrack_1, PrimaryVertexObject)") # multiplicity #.Define("SV_evt2_n","VertexingUtils::get_n_SV(SV_evt2)") @@ -119,11 +120,6 @@ def output(): #'SV_evt1_position', #'SV_evt2_position', - # SV chi2 - 'SV_jet_chi2', - #'SV_evt1_chi2', - #'SV_evt2_chi2', - # more SV variables 'sv_mass', # 'sv_p', diff --git a/examples/FCCee/vertex_lcfiplus/analysis_V0.py b/examples/FCCee/vertex_lcfiplus/analysis_V0.py index ebd611685f..d6bcc35505 100644 --- a/examples/FCCee/vertex_lcfiplus/analysis_V0.py +++ b/examples/FCCee/vertex_lcfiplus/analysis_V0.py @@ -1,3 +1,5 @@ +testFile="root://eospublic.cern.ch//eos/experiment/fcc/ee/generation/DelphesEvents/winter2023/IDEA/p8_ee_Zbb_ecm91/events_066726720.root" + #Mandatory: List of processes processList = { 'p8_ee_Zuds_ecm91':{'fraction':0.0001}, #Run 0.01% statistics in one output file named /p8_ee_Zuds_ecm91.root @@ -44,10 +46,9 @@ def analysers(df): # determime the primary (and secondary) tracks without using the MC-matching: - # First, reconstruct a vertex from all tracks - .Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300)") # Select the tracks that are reconstructed as primaries - .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( VertexObject_allTracks, EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0., 0)") + .Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)") + .Define("n_RecoedPrimaryTracks", "ReconstructedParticle2Track::getTK_n( RecoedPrimaryTracks )") # the final primary vertex : .Define("PrimaryVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300) ") diff --git a/examples/basics/README.md b/examples/basics/README.md index b9db0c6c55..1cecc3d183 100644 --- a/examples/basics/README.md +++ b/examples/basics/README.md @@ -56,11 +56,36 @@ root[0] TBrowser b As shown in the screenshot above, there are two types of branches: - - Branches without a pound (#) in their name: Electron (1), Muon (2), AllMuon (3), EFlowNeutralHadron (4), Particle (5), Photon (6), ReconstructedParticles (7), EFlowPhoton (8), MCRecoAssociations (9), MissingET (10), ParticleIDs (11), Jet (12), EFlowTrack (13), EFlowTrack\_1 (14). They refer to collections of objects. + - Branches without a pound (#) in their name refer to collections of objects: Electron (1), Muon (2), AllMuon (3), EFlowNeutralHadron (4), Particle (5), Photon (6), ReconstructedParticles (7), EFlowPhoton (8), MCRecoAssociations (9), MissingET (10), ParticleIDs (11), Jet (12), EFlowTrack (13), EFlowTrack\_1 (14). - NB: "Particle" denotes the collection of Monte-Carlo particles. "Muon" contains the isolated muons, while "AllMuon" contains all muons, isolated or not. - Branches with a pound in their name: Each of the object collections listed above, e.g. "Collection", has up to six associated collections of references, i.e. indices that point to another or to the same object collection. They are labeled Collection#i, with i = 0 ... 5. For example, the Muon collection has one single - associated collection of references, Muon#0. + associated collection of references, Muon#0. + - NB2: With `winter2023` samples, the `CollectionID`s are different. The correct list can be obtained using `podio-dump` or `collInfo` script from [AuxTools](https://github.com/HEP-FCC/AuxTools), e.g.: + +``` +./collInfo -i /eos/experiment/fcc/ee/generation/DelphesEvents/winter2023/IDEA/p8_ee_ZZ_ecm240/events_051628351.root + +ID -> Collection +================== +1 -> MissingET +2 -> MCRecoAssociations +3 -> ParticleIDs +4 -> magFieldBz +5 -> TrackerHits +6 -> EFlowTrack +7 -> CalorimeterHits +8 -> Particle +9 -> Photon +10 -> EFlowTrack_L +11 -> Electron +12 -> EFlowPhoton +13 -> EFlowNeutralHadron +14 -> Jet +15 -> ReconstructedParticles +16 -> Muon +================== +``` To figure out which collection is pointed to by Muon#0 (or by any other collection of references), one can look at the value of Muon#0.collectionID (see screenshot below). The collectionID of Muon#0 is the collection number 7 (in the example file used here), which, in the list of "object collections" above, corresponds to the collection of ReconstructedParticles. diff --git a/setup.sh b/setup.sh index f8c09a38ed..d5dd2e4f28 100644 --- a/setup.sh +++ b/setup.sh @@ -21,8 +21,12 @@ if [ "${0}" != "${BASH_SOURCE}" ]; then export CMAKE_PREFIX_PATH=${LOCAL_DIR}/install:${CMAKE_PREFIX_PATH} export ROOT_INCLUDE_PATH=${LOCAL_DIR}/install/include:${ROOT_INCLUDE_PATH} - export ONNXRUNTIME_ROOT_DIR=`python -c "import onnxruntime; print(onnxruntime.__path__[0]+'/../../../..')"` - export LD_LIBRARY_PATH=$ONNXRUNTIME_ROOT_DIR/lib:$LD_LIBRARY_PATH + export ONNXRUNTIME_ROOT_DIR=`python -c "import onnxruntime; print(onnxruntime.__path__[0]+'/../../../..')" 2> /dev/null` + if [ -z "${ONNXRUNTIME_ROOT_DIR}" ]; then + echo "WARNING: ONNX Runtime not found! Related analyzers won't be build..." + else + export LD_LIBRARY_PATH=${ONNXRUNTIME_ROOT_DIR}/lib:${LD_LIBRARY_PATH} + fi else echo "ERROR: This script is meant to be sourced!" fi