diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/CMakeLists.txt b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/CMakeLists.txt
index 7635b18fc8136c9a70691a1e3e7d39b26356c4e7..971985c8edd16b085a70c8579241a351aee1bedc 100644
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/CMakeLists.txt
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/CMakeLists.txt
@@ -28,15 +28,25 @@ atlas_depends_on_subdirs( PUBLIC
# External dependencies:
find_package( Boost COMPONENTS filesystem thread system )
find_package( Eigen )
-find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
+find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread TMVA )
+
+# Library to link against:
+#atlas_add_library( InDetVKalVxInJetToolLib
+# src/*.cxx src/components/*.cxx
+# PUBLIC_HEADERS InDetVKalVxInJetTool InDetTrkInJetType
+# INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}
+# LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps xAODTracking GaudiKernel
+# InDetRecToolInterfaces Particle TrkVKalVrtFitterLib GeoPrimitives AnalysisUtilsLib TrkNeutralParameters
+# TrkParticleBase TrkTrackSummary VxSecVertex VxVertex TrkToolInterfaces TrkVertexFitterInterfaces PathResolver )
# Component(s) in the package:
-atlas_add_component( InDetVKalVxInJetTool
- src/*.cxx
- src/components/*.cxx
+atlas_add_component( InDetVKalVxInJetTool
+ src/*.cxx src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps xAODTracking GaudiKernel InDetRecToolInterfaces Particle TrkVKalVrtFitterLib GeoPrimitives AnalysisUtilsLib TrkNeutralParameters TrkParticleBase TrkTrackSummary VxSecVertex VxVertex TrkToolInterfaces TrkVertexFitterInterfaces )
-
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps xAODTracking GaudiKernel
+ InDetRecToolInterfaces Particle TrkVKalVrtFitterLib GeoPrimitives AnalysisUtilsLib TrkNeutralParameters
+ TrkParticleBase TrkTrackSummary VxSecVertex VxVertex TrkToolInterfaces TrkVertexFitterInterfaces PathResolver )
+
# Install files from the package:
atlas_install_headers( InDetVKalVxInJetTool )
atlas_install_python_modules( python/*.py )
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h
new file mode 100644
index 0000000000000000000000000000000000000000..f192f9ef81402a3383385b03a531d4490b338b8b
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h
@@ -0,0 +1,106 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+//
+// InDetTrkInJetType.h - Description
+//
+/*
+ Tool to classify track origins in a jet.
+ Track types:
+ 0 - Heavy Flavour (Signal)
+ 1 - Fragmentation tracks (Fragment)
+ 2 - Garbage (Interactions+V0s+Pileup)
+
+ The tool works for JetPt<2.5TeV.
+ Tool is trained using ttbar+Z'(1.5,3,5TeV)+JZ4,5,6 samples
+
+ Author: Vadim Kostyukhin
+ e-mail: vadim.kostyukhin@cern.ch
+*/
+#ifndef InDet_InDetTrkInJetType_H
+#define InDet_InDetTrkInJetType_H
+
+#include <vector>
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "TLorentzVector.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+#include "Particle/TrackParticle.h"
+//
+
+namespace TMVA { class Reader; }
+
+namespace InDet {
+
+//------------------------------------------------------------------------
+ static const InterfaceID IID_IInDetTrkInJetType("IInDetTrkInJetType", 1, 0);
+
+ class IInDetTrkInJetType : virtual public IAlgTool {
+ public:
+ static const InterfaceID& interfaceID() { return IID_IInDetTrkInJetType;}
+//---------------------------------------------------------------------------
+//Interface itself
+
+ virtual std::vector<float> trkTypeWgts( const xAOD::TrackParticle *, const xAOD::Vertex &, const TLorentzVector &) =0;
+ virtual std::vector<float> trkTypeWgts( const Rec::TrackParticle *, const xAOD::Vertex &, const TLorentzVector &) =0;
+
+ };
+
+
+
+
+ class InDetTrkInJetType : public AthAlgTool, virtual public IInDetTrkInJetType
+ {
+
+ public:
+ /* Constructor */
+ InDetTrkInJetType(const std::string& type, const std::string& name, const IInterface* parent);
+ /* Destructor */
+ virtual ~InDetTrkInJetType();
+
+
+ StatusCode initialize();
+ StatusCode finalize();
+
+ std::vector<float> trkTypeWgts(const xAOD::TrackParticle *, const xAOD::Vertex &, const TLorentzVector &);
+ std::vector<float> trkTypeWgts(const Rec::TrackParticle *, const xAOD::Vertex &, const TLorentzVector &);
+
+//------------------------------------------------------------------------------------------------------------------
+// Private data and functions
+//
+
+ private:
+
+ TMVA::Reader* m_tmvaReader;
+ double m_trkMinPtCut;
+ float m_d0_limLow;
+ float m_d0_limUpp;
+ float m_Z0_limLow;
+ float m_Z0_limUpp;
+ std::string m_calibFileName;
+ ToolHandle < Trk::IVertexFitter > m_fitterSvc;
+ Trk::TrkVKalVrtFitter* m_fitSvc;
+
+ int m_initialised;
+
+ float m_prbS;
+ float m_Sig3D;
+ float m_prbP;
+ float m_d0;
+ float m_vChi2;
+ float m_pTvsJet;
+ float m_prodTJ;
+ float m_SigZ;
+ float m_SigR;
+ float m_ptjet;
+ float m_etajet;
+ float m_ibl;
+ float m_bl;
+ };
+
+
+
+
+} //end namespace
+#endif
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
index 42e82086d62958057bec1d4e562acaeecedb5a56..cfb8d0c68ede6af5645d378815a310a5d245551b 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
//
@@ -50,13 +50,15 @@
#include "xAODTracking/TrackParticleContainer.h"
//
//
-
#include "InDetRecToolInterfaces/ISecVertexInJetFinder.h"
+//#include "InDetMaterialRejTool/InDetMaterialRejTool.h"
+#include "InDetVKalVxInJetTool/InDetTrkInJetType.h"
class TH1D;
class TH2D;
class TH1F;
class TProfile;
+class TTree;
namespace Trk{
class TrkVKalVrtFitter;
@@ -70,12 +72,10 @@ typedef std::vector<double> dvect;
//------------------------------------------------------------------------
namespace InDet {
-
class VKalVxInJetTemp{
public:
std::vector<int> m_Incomp;
std::vector<int> m_Prmtrack;
-
};
struct workVectorArrxAOD{
@@ -120,6 +120,17 @@ namespace InDet {
const TLorentzVector & jetMomentum,
const std::vector<const xAOD::IParticle*> & inputTracks) const;
+//--------------------------------------------- For package debugging
+ public:
+ void setVBTrk(std::vector<const xAOD::TrackParticle*> * BTrk) const { m_dbgVBTrk=BTrk; };
+ private:
+ int isBTrk(const xAOD::TrackParticle* prt) const
+ { if(!m_dbgVBTrk) return 0;
+ if(std::find((*m_dbgVBTrk).begin(),(*m_dbgVBTrk).end(), prt) != (*m_dbgVBTrk).end()) return 1;
+ return 0;
+ }
+ int isBTrk(const Rec::TrackParticle*) const { return 0; }
+ mutable std::vector<const xAOD::TrackParticle*> *m_dbgVBTrk=nullptr;
//------------------------------------------------------------------------------------------------------------------
// Private data and functions
//
@@ -127,6 +138,7 @@ namespace InDet {
private:
double m_w_1{};
+ TTree* m_tuple{};
TH1D* m_hb_massPiPi{};
TH1D* m_hb_massPiPi1{};
TH1D* m_hb_massPPi{};
@@ -141,6 +153,7 @@ namespace InDet {
TH1D* m_hb_impact{};
TH1D* m_hb_impactR{};
TH2D* m_hb_impactRZ{};
+ TH1D* m_hb_trkD0{};
TH1D* m_hb_ntrkjet{};
TH1D* m_hb_impactZ{};
TH1D* m_hb_r2d{};
@@ -161,7 +174,6 @@ namespace InDet {
TH1D* m_hb_sig3D2tr{};
TH1D* m_hb_sig3DNtr{};
TH1D* m_hb_trkPtMax{};
- TH1D* m_hb_tr2SelVar{};
TH1F* m_hb_blshared{};
TH1F* m_hb_pxshared{};
TH1F* m_hb_rawVrtN{};
@@ -169,12 +181,7 @@ namespace InDet {
TH1F* m_hb_trkPErr{};
TH1F* m_hb_deltaRSVPV{};
//--
- TH1D* m_hb_massJetTrkSV{};
- TH1D* m_hb_ratioJetTrkSV{};
- TH1D* m_hb_DST_JetTrkSV{};
- TH1D* m_hb_NImpJetTrkSV{};
- TH1D* m_hb_nHImpTrkCnt{};
-//--
+ TProfile * m_pr_NSelTrkMean{};
TProfile * m_pr_effVrt2tr{};
TProfile * m_pr_effVrt2trEta{};
TProfile * m_pr_effVrt{};
@@ -187,69 +194,60 @@ namespace InDet {
- long int m_CutSctHits{};
- long int m_CutPixelHits{};
- long int m_CutSiHits{};
- long int m_CutBLayHits{};
- long int m_CutSharedHits{};
- double m_CutPt{};
- double m_CutZVrt{};
- double m_CutA0{};
- double m_CutChi2{};
- double m_SecTrkChi2Cut{};
- double m_ConeForTag{};
- double m_Sel2VrtChi2Cut{};
- double m_Sel2VrtSigCut{};
- double m_TrkSigCut{};
- double m_TrkSigNTrkDep{};
- double m_TrkSigSumCut{};
- double m_A0TrkErrorCut{};
- double m_ZTrkErrorCut{};
- double m_AntiPileupSigRCut{};
- double m_AntiPileupSigZCut{};
- double m_AntiFake2trVrtCut{};
- double m_JetPtFractionCut{};
- int m_TrackInJetNumberLimit{};
- double m_pseudoSigCut{};
- double m_hadronIntPtCut{};
-
- bool m_FillHist{};
+ long int m_cutSctHits{};
+ long int m_cutPixelHits{};
+ long int m_cutSiHits{};
+ long int m_cutBLayHits{};
+ long int m_cutSharedHits{};
+ double m_cutPt{};
+ double m_cutZVrt{};
+ double m_cutA0{};
+ double m_cutChi2{};
+ double m_secTrkChi2Cut{};
+ double m_coneForTag{};
+ double m_sel2VrtChi2Cut{};
+ double m_sel2VrtSigCut{};
+ double m_trkSigCut{};
+ double m_a0TrkErrorCut{};
+ double m_zTrkErrorCut{};
+ double m_cutHFClass{};
+ double m_antiGarbageCut{};
+
+ bool m_fillHist{};
bool m_existIBL{};
long int m_RobustFit{};
- double m_Xbeampipe{};
- double m_Ybeampipe{};
- double m_XlayerB{};
- double m_YlayerB{};
- double m_Xlayer1{};
- double m_Ylayer1{};
- double m_Xlayer2{};
- double m_Ylayer2{};
- double m_Rbeampipe{};
- double m_RlayerB{};
- double m_Rlayer1{};
- double m_Rlayer2{};
- double m_Rlayer3{};
- double m_SVResolutionR{};
-
- bool m_useMaterialRejection{};
+ double m_beampipeX{};
+ double m_beampipeY{};
+ double m_xLayerB{};
+ double m_yLayerB{};
+ double m_xLayer1{};
+ double m_yLayer1{};
+ double m_xLayer2{};
+ double m_yLayer2{};
+ double m_beampipeR{};
+ double m_rLayerB{};
+ double m_rLayer1{};
+ double m_rLayer2{};
+ double m_rLayer3{};
+
bool m_useVertexCleaning{};
- int m_MassType{};
- bool m_MultiVertex{};
- bool m_MultiWithPrimary{};
+ bool m_multiVertex{};
+ bool m_multiWithPrimary{};
bool m_getNegativeTail{};
bool m_getNegativeTag{};
- bool m_MultiWithOneTrkVrt{};
+ bool m_multiWithOneTrkVrt{};
- double m_VertexMergeCut{};
- double m_TrackDetachCut{};
+ double m_vertexMergeCut{};
+ double m_trackDetachCut{};
ToolHandle < Trk::IVertexFitter > m_fitterSvc;
Trk::TrkVKalVrtFitter* m_fitSvc{};
+ ToolHandle < IInDetTrkInJetType > m_trackClassificator;
double m_massPi {};
double m_massP {};
@@ -257,9 +255,70 @@ namespace InDet {
double m_massK0{};
double m_massLam{};
double m_massB{};
- mutable int m_NRefTrk{}; //Measure of track in jet multiplicity
+ mutable int m_NRefPVTrk{}; //Measure of track in jet multiplicity
std::string m_instanceName;
+//-------------------------------------------
+//For ntuples (only for development/tuning!)
+
+ int getG4Inter( const Rec::TrackParticle* TP ) const;
+ int getG4Inter( const xAOD::TrackParticle* TP ) const;
+ int getMCPileup(const Rec::TrackParticle* TP ) const;
+ int getMCPileup(const xAOD::TrackParticle* TP ) const;
+
+ struct DevTuple
+ {
+ static const int maxNTrk=100;
+ int nTrkInJet;
+ float ptjet;
+ float etajet;
+ float phijet;
+ float p_prob[maxNTrk];
+ float s_prob[maxNTrk];
+ int idMC[maxNTrk];
+ float SigR[maxNTrk];
+ float SigZ[maxNTrk];
+ float d0[maxNTrk];
+ float Z0[maxNTrk];
+ float pTvsJet[maxNTrk];
+ float prodTJ[maxNTrk];
+ float wgtB[maxNTrk];
+ float wgtL[maxNTrk];
+ float wgtG[maxNTrk];
+ float Sig3D[maxNTrk];
+ int chg[maxNTrk];
+ int nVrtT[maxNTrk];
+ int nVrt;
+ float VrtDist2D[maxNTrk];
+ float VrtSig3D[maxNTrk];
+ float VrtSig2D[maxNTrk];
+ float mass[maxNTrk];
+ float Chi2[maxNTrk];
+ int itrk[maxNTrk];
+ int jtrk[maxNTrk];
+ int badVrt[maxNTrk];
+ int ibl[maxNTrk];
+ int bl[maxNTrk];
+ int NTHF;
+ int itHF[maxNTrk];
+ };
+ DevTuple* m_curTup;
+
+ struct Vrt2Tr
+ { int i=0, j=0;
+ int badVrt=0;
+ Amg::Vector3D FitVertex;
+ TLorentzVector Momentum;
+ long int vertexCharge;
+ std::vector<double> ErrorMatrix;
+ std::vector<double> Chi2PerTrk;
+ std::vector< std::vector<double> > TrkAtVrt;
+ double Signif3D=0.;
+ double Signif3DProj=0.;
+ double Signif2D=0.;
+ double Chi2=0.;
+ };
+
// For multivertex version only
@@ -280,10 +339,21 @@ namespace InDet {
double dCloseVrt{};
double ProjectedVrt{};
int detachedTrack=-1;
- };
+ };
//
+// Private struct for track oredering (happened to be not needed)
+//
+// struct compareTrk {
+// compareTrk(const TLorentzVector J):JetDir(J.Vect()){};
+// bool operator() (const xAOD::TrackParticle* p1, const xAOD::TrackParticle*p2) {
+// double cmp1=p1->p4().Rho()*(1.-sin(p1->p4().Angle(JetDir)));
+// double cmp2=p2->p4().Rho()*(1.-sin(p2->p4().Angle(JetDir)));
+// return (cmp1>cmp2); }
+// TVector3 JetDir;
+// };
+
// Private technical functions
//
//
@@ -306,13 +376,6 @@ namespace InDet {
const TLorentzVector & JetDir,
std::vector<double> & Results) const;
- xAOD::Vertex* tryPseudoVertex(const std::vector<const xAOD::TrackParticle*>& Tracks,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- const TLorentzVector & TrkJet,
- const int & nTrkLead,
- std::vector<double> & Results) const;
-
void TrackClassification(std::vector<WrkVrt> *WrkVrtSet,
std::vector< std::deque<long int> > *TrkInVrt) const;
@@ -336,7 +399,7 @@ namespace InDet {
//
// Gives correct mass assignment in case of nonequal masses
double massV0( std::vector< std::vector<double> >& TrkAtVrt, double massP, double massPi ) const;
- int FindMax( std::vector<double>& Chi2PerTrk, std::vector<int>& countT) const;
+ int FindMax( std::vector<double>& Chi2PerTrk, std::vector<float>& rank) const;
TLorentzVector TotalMom(const std::vector<const Trk::Perigee*>& InpTrk) const;
@@ -345,8 +408,10 @@ namespace InDet {
double TotalTMom(const std::vector<const Trk::Perigee*> & InpTrk) const;
double TotalTVMom(const Amg::Vector3D &Dir, const std::vector<const Trk::Perigee*>& InpTrk) const;
double pTvsDir(const Amg::Vector3D &Dir, const std::vector<double>& InpTrk) const;
+ double VrtRadiusError(const Amg::Vector3D & SecVrt, const std::vector<double> & VrtErr) const;
- bool insideMatLayer(float ,float ) const;
+ bool insideMatLayer(float ,float ) const;
+ void fillVrtNTup( std::vector<Vrt2Tr> & all2TrVrt) const;
TLorentzVector GetBDir( const xAOD::TrackParticle* trk1,
const xAOD::TrackParticle* trk2,
@@ -364,6 +429,8 @@ namespace InDet {
const std::vector<double> VrtErr,double& Signif ) const;
double VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
const std::vector<double> VrtErr,double& Signif ) const;
+ double VrtVrtDist2D(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> VrtErr,double& Signif ) const;
double VrtVrtDist(const Trk::RecVertex & PrimVrt, const Amg::Vector3D & SecVrt,
const std::vector<double> SecVrtErr, const TLorentzVector & JetDir) const;
double VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
@@ -375,8 +442,10 @@ namespace InDet {
void DisassembleVertex(std::vector<WrkVrt> *WrkVrtSet, int iv,
std::vector<const Particle*> AllTracks) const;
- double ProjPos(const Amg::Vector3D & Vrt, const TLorentzVector & JetDir) const;
- double ProjPosT(const Amg::Vector3D & Vrt, const TLorentzVector & JetDir) const;
+ double ProjSV_PV(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Jet) const;
+
+ double RankBTrk(double TrkPt, double JetPt, double Signif) const;
+
const Trk::Perigee* GetPerigee( const xAOD::TrackParticle* ) const;
const Trk::Perigee* GetPerigee( const Rec::TrackParticle* ) const;
@@ -386,7 +455,7 @@ namespace InDet {
template <class Trk>
double FitCommonVrt(std::vector<const Trk*>& ListSecondTracks,
- std::vector<int> & cntComb,
+ std::vector<float> & trkRank,
const xAOD::Vertex & PrimVrt,
const TLorentzVector & JetDir,
std::vector<double> & InpMass,
@@ -396,13 +465,10 @@ namespace InDet {
std::vector< std::vector<double> > & TrkAtVrt) const;
template <class Trk>
- void RemoveEntryInList(std::vector<const Trk*>& , std::vector<int>&, int) const;
+ void RemoveEntryInList(std::vector<const Trk*>& , std::vector<float>&, int) const;
template <class Trk>
void RemoveDoubleEntries(std::vector<const Trk*>& ) const;
- template <class Trk>
- double RemoveNegImpact(std::vector<const Trk*>& , const xAOD::Vertex &, const TLorentzVector&, const double ) const;
-
template <class Particle>
StatusCode RefitVertex( std::vector<WrkVrt> *WrkVrtSet, int SelectedVertex,
std::vector<const Particle*> & SelectedTracks) const;
@@ -428,14 +494,6 @@ namespace InDet {
const xAOD::Vertex & PrimVrt,
const TLorentzVector & JetDir,
std::vector<const xAOD::TrackParticle*>& SelPart) const;
- int SelGoodTrkParticleRelax( const std::vector<const Rec::TrackParticle*>& InpPart,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- std::vector<const Rec::TrackParticle*>& SelPart) const;
- int SelGoodTrkParticleRelax( const std::vector<const xAOD::TrackParticle*>& InpPart,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- std::vector<const xAOD::TrackParticle*>& SelPart) const;
template <class Trk>
@@ -449,7 +507,6 @@ namespace InDet {
Amg::MatrixX makeVrtCovMatrix( std::vector<double> & ErrorMatrix ) const;
- double trkPtCorr(double pT) const;
Amg::MatrixX SetFullMatrix(int NTrk, std::vector<double> & Matrix) const;
@@ -462,11 +519,9 @@ namespace InDet {
StatusCode VKalVrtFitFastBase(const std::vector<const xAOD::TrackParticle*>& listPart,Amg::Vector3D& Vertex) const;
template <class Track>
- bool Check2TrVertexInPixel( const Track* p1, const Track* p2, Amg::Vector3D &, double ) const;
- template <class Track>
- bool Check1TrVertexInPixel( const Track* p1, Amg::Vector3D &) const;
+ bool Check2TrVertexInPixel( const Track* p1, const Track* p2, Amg::Vector3D &, std::vector<double> &) const;
template <class Track>
- double medianPtF( std::vector<const Track*> & Trks) const;
+ bool Check1TrVertexInPixel( const Track* p1, Amg::Vector3D &, std::vector<double> & ) const;
void getPixelLayers(const Rec::TrackParticle* Part, int &blHit, int &l1Hit, int &l2Hit, int &nLay) const;
void getPixelLayers(const xAOD::TrackParticle* Part, int &blHit, int &l1Hit, int &l2Hit, int &nLay) const;
@@ -496,16 +551,17 @@ namespace InDet {
PartToBase(const std::vector<const Rec::TrackParticle*> & listPart) const;
const std::vector<const Rec::TrackParticle*>
BaseToPart(const std::vector<const Trk::TrackParticleBase*> & listBase) const;
+ StatusCode GetTrkFitWeights(std::vector<double> & wgt) const;
};
template <class Trk>
- void InDetVKalVxInJetTool::RemoveEntryInList(std::vector<const Trk*>& ListTracks, std::vector<int> &cnt, int Outlier) const
+ void InDetVKalVxInJetTool::RemoveEntryInList(std::vector<const Trk*>& ListTracks, std::vector<float> & rank, int Outlier) const
{
if(Outlier < 0 ) return;
if(Outlier >= (int)ListTracks.size() ) return;
ListTracks.erase( ListTracks.begin()+Outlier);
- cnt.erase( cnt.begin()+Outlier);
+ rank.erase( rank.begin()+Outlier);
}
template <class Trk>
@@ -534,39 +590,39 @@ namespace InDet {
};
template <class Track>
- bool InDetVKalVxInJetTool::Check1TrVertexInPixel( const Track* p1, Amg::Vector3D &FitVertex)
+ bool InDetVKalVxInJetTool::Check1TrVertexInPixel( const Track* p1, Amg::Vector3D &FitVertex, std::vector<double> &VrtCov)
const
{
int blTrk=0, blP=0, l1Trk=0, l1P=0, l2Trk=0, nLays=0;
getPixelLayers( p1, blTrk , l1Trk, l2Trk, nLays );
getPixelProblems(p1, blP, l1P );
- double xDif=FitVertex.x()-m_XlayerB, yDif=FitVertex.y()-m_YlayerB ;
+ double xDif=FitVertex.x()-m_xLayerB, yDif=FitVertex.y()-m_yLayerB ;
double Dist2DBL=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2DBL < m_RlayerB-m_SVResolutionR){ //----------------------------------------- Inside B-layer
+ if (Dist2DBL < m_rLayerB-VrtRadiusError(FitVertex, VrtCov)){ //----------------------------------------- Inside B-layer
if( blTrk<1 && l1Trk<1 ) return false;
if( nLays <2 ) return false; // Less than 2 layers on track 0
return true;
- }else if(Dist2DBL > m_RlayerB+m_SVResolutionR){ //----------------------------------------- Outside b-layer
+ }else if(Dist2DBL > m_rLayerB+VrtRadiusError(FitVertex, VrtCov)){ //----------------------------------------- Outside b-layer
if( blTrk>0 && blP==0 ) return false; // Good hit in b-layer is present
}
//
// L1 and L2 are considered only if vertex is in acceptance
//
if(fabs(FitVertex.z())<400.){
- xDif=FitVertex.x()-m_Xlayer1, yDif=FitVertex.y()-m_Ylayer1 ;
+ xDif=FitVertex.x()-m_xLayer1, yDif=FitVertex.y()-m_yLayer1 ;
double Dist2DL1=sqrt(xDif*xDif+yDif*yDif);
- xDif=FitVertex.x()-m_Xlayer2, yDif=FitVertex.y()-m_Ylayer2 ;
+ xDif=FitVertex.x()-m_xLayer2, yDif=FitVertex.y()-m_yLayer2 ;
double Dist2DL2=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2DL1 < m_Rlayer1-m_SVResolutionR) { //------------------------------------------ Inside 1st-layer
+ if (Dist2DL1 < m_rLayer1-VrtRadiusError(FitVertex, VrtCov)) { //------------------------------------------ Inside 1st-layer
if( l1Trk<1 && l2Trk<1 ) return false; // Less than 1 hits on track 0
return true;
- }else if(Dist2DL1 > m_Rlayer1+m_SVResolutionR) { //------------------------------------------- Outside 1st-layer
+ }else if(Dist2DL1 > m_rLayer1+VrtRadiusError(FitVertex, VrtCov)) { //------------------------------------------- Outside 1st-layer
if( l1Trk>0 && l1P==0 ) return false; // Good L1 hit is present
}
- if (Dist2DL2 < m_Rlayer2-m_SVResolutionR) { //------------------------------------------- Inside 2nd-layer
+ if (Dist2DL2 < m_rLayer2-VrtRadiusError(FitVertex, VrtCov)) { //------------------------------------------- Inside 2nd-layer
if( l2Trk==0 ) return false; // At least one L2 hit must be present
- }else if(Dist2DL2 > m_Rlayer2+m_SVResolutionR) {
+ }else if(Dist2DL2 > m_rLayer2+VrtRadiusError(FitVertex, VrtCov)) {
// if( l2Trk>0 ) return false; // L2 hits are present
}
} else {
@@ -577,11 +633,6 @@ namespace InDet {
return true;
}
- template <class Track>
- double InDetVKalVxInJetTool::medianPtF( std::vector<const Track*> & Trks) const{
- int NT=Trks.size(); if(NT==0)return 0.;
- return (Trks[(NT-1)/2]->pt()+Trks[NT/2]->pt())/2.;
- }
} //end namespace
#endif
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/Train.C b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/Train.C
new file mode 100644
index 0000000000000000000000000000000000000000..62242ed720c222fbab0a0168d0a7b0e669cec9b3
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/Train.C
@@ -0,0 +1,137 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+/*-------------------------------------------------------------------------
+ Script for multiclass TMVA training using track origin information.
+ Don't forget to set AntiGarbageCut=1 in SSVF to create the training data.
+ Author: V.Kostyukhin
+*/
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <string>
+
+#include "TChain.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TString.h"
+#include "TObjString.h"
+#include "TSystem.h"
+#include "TROOT.h"
+
+#include "TMVA/Factory.h"
+#include "TMVA/Tools.h"
+#include "TMVA/TMVAGui.h"
+
+
+void TrainMulti( )
+{
+ TString outfileName( "TMVAMulti.root" );
+ TFile* outputFile = TFile::Open( outfileName, "RECREATE" );
+ TMVA::Factory *factory = new TMVA::Factory( "TMVAMulticlass", outputFile,
+ "!V:!Silent:Color:DrawProgressBar:Transformations=I;P;G,D:AnalysisType=multiclass" );
+ TMVA::DataLoader *dataloader=new TMVA::DataLoader("dataset");
+
+ //dataloader->AddVariable( "prbS", 'F' );
+ dataloader->AddVariable( "Sig3D", 'F' );
+ dataloader->AddVariable( "prbP", 'F' );
+ dataloader->AddVariable( "pTvsJet", 'F' );
+ //dataloader->AddVariable( "prodTJ", 'F' );
+ dataloader->AddVariable( "d0", 'F' );
+ dataloader->AddVariable( "SigR", 'F' );
+ dataloader->AddVariable( "SigZ", 'F' );
+ dataloader->AddVariable( "ptjet", 'F' );
+ dataloader->AddVariable( "ibl", 'I' );
+ dataloader->AddVariable( "bl", 'I' );
+ dataloader->AddVariable( "etajet", 'F' );
+ //dataloader->AddVariable( "vChi2", 'F' );
+ //dataloader->AddSpectator( "idMC", 'I' );
+
+ TChain *chainB =new TChain("/stat/SVrtInJetTrackExample.BJetSVFinder/Tracks");
+ chainB->Add("../run/allcalib.root");
+ chainB->Add("../runZp1500/allcalib.root");
+ chainB->Add("../runZp3000/allcalib.root");
+ chainB->Add("../runZp5000/allcalib.root");
+ chainB->Add("../runLJ4/allcalib.root");
+ chainB->Add("../runLJ6/allcalib.root");
+ chainB->Add("../runLJ8/allcalib.root");
+ TChain *chainL =new TChain("/stat/SVrtInJetTrackExample.LJetSVFinder/Tracks");
+ chainL->Add("../run/allcalib.root");
+ chainL->Add("../runZp1500/allcalib.root");
+ chainL->Add("../runZp3000/allcalib.root");
+ chainL->Add("../runZp5000/allcalib.root");
+ chainL->Add("../runLJ4/allcalib.root");
+ chainL->Add("../runLJ6/allcalib.root");
+ chainL->Add("../runLJ8/allcalib.root");
+
+ dataloader->AddTree(chainB,"Signal");
+ dataloader->AddTree(chainL,"Signal");
+ dataloader->AddTree(chainB,"Fragment");
+ dataloader->AddTree(chainL,"Fragment");
+ dataloader->AddTree(chainB,"Garbage");
+ dataloader->AddTree(chainL,"Garbage");
+ //dataloader->AddTree(chainB,"Pileup");
+ //dataloader->AddTree(chainL,"Pileup");
+
+ //dataloader->SetCut(TCut("idMC==2&&prbS>0&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12."),"Signal");
+ //dataloader->SetCut(TCut("idMC==0&&prbS>0&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12."),"Fragment");
+ //dataloader->SetCut(TCut("idMC==1&&prbS>0&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12."),"Garbage");
+ //dataloader->SetCut(TCut("idMC==3&&prbS>0&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12."),"Pileup");
+ dataloader->SetCut(TCut("idMC==2&&(SigR*SigR+SigZ*SigZ)>1.&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12.&&ptjet<2500000"),"Signal");
+ dataloader->SetCut(TCut("idMC==0&&(SigR*SigR+SigZ*SigZ)>1.&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12.&&ptjet<2500000"),"Fragment");
+ dataloader->SetCut(TCut("(idMC==1||idMC==3)&&(SigR*SigR+SigZ*SigZ)>1.&&(-3)<d0&&d0<5.&&(-8)<Z0&&Z0<12.&&ptjet<2500000"),"Garbage");
+
+ dataloader->PrepareTrainingAndTestTree( "", "SplitMode=Random:NormMode=NumEvents:!V" );
+ //dataloader->PrepareTrainingAndTestTree( "", "SplitMode=Random:NormMode=NumEvents:!V:nTrain_Signal=5000:nTrain_Fragment=5000:nTrain_Garbage=10000:nTest_Signal=5000:nTest_Fragment=5000:nTest_Garbage=10000" );
+ //dataloader->PrepareTrainingAndTestTree( "", "SplitMode=Random:NormMode=NumEvents:!V:nTrain_Signal=5000:nTrain_Fragment=5000:nTrain_Garbage=5000:nTrain_Pileup=5000:nTest_Signal=5000:nTest_Fragment=5000:nTest_Garbage=5000:nTest_Pileup=5000" );
+
+ factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", "!H:!V:NTrees=1000:BoostType=Grad:Shrinkage=0.25:UseBaggedBoost:BaggedSampleFraction=0.50:nCuts=20:MaxDepth=5");
+ //factory->BookMethod( TMVA::Types::kPDEFoam, "PDEFoam", "!H:!V:TailCut=0.001:VolFrac=0.0666:nActiveCells=500:nSampl=2000:nBin=5:Nmin=100:Kernel=None:Compress=T" );
+ // Train MVAs using the set of training events
+ factory->TrainAllMethods();
+
+ // ---- Evaluate all MVAs using the set of test events
+ factory->TestAllMethods();
+
+ // ----- Evaluate and compare performance of all configured MVAs
+ factory->EvaluateAllMethods();
+
+ // --------------------------------------------------------------
+
+ // Save the output
+ outputFile->Close();
+
+ std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;
+ std::cout << "==> TMVAClassification is done!" << std::endl;
+
+ delete factory;
+}
+
+//
+
+void TMVAMultiApply( )
+{
+ // create the Reader object
+ TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );
+ Float_t prbS, prbP, ptjet, etajet;
+ Int_t ntrk;
+ reader->AddVariable( "prbS", &prbS );
+ reader->AddVariable( "prbP", &prbP );
+ reader->AddVariable( "ptjet", &ptjet );
+ reader->AddVariable( "etajet", &etajet );
+ //reader->AddSpectator( "ntrk", &ntrk );
+
+ TString methodName = "BDTG";
+ TString weightfile = "weights/TMVAMulticlass_BDTG.weights.xml";
+ reader->BookMVA( methodName, weightfile );
+
+ prbS=0.6;
+ prbP=0.9;
+ ptjet=250000.;
+ etajet=1.;
+ const std::vector<Float_t> &result=reader->EvaluateMulticlass("BDTG");
+ std::cout<<" TEST="<<result.size()<<" wgt="<<result[0]<<","<<result[1]<<","<<result[2]<<'\n';
+
+ delete reader;
+}
+
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
index 780819de28b06e29dccf38edaa67eab32bd4dd2c..e18a662ac0508e9aaa2caf0cdf34e603ef6b097a 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// Header include
@@ -15,6 +15,7 @@
#include "TH2D.h"
#include "TMath.h"
//
+//#include<iostream>
//----------------------------------------------------------------------------------------
// GetVrtSec resurns the vector Results with the following
@@ -32,8 +33,20 @@
//----------------------------------------------------------------------------------------
namespace InDet{
- //inline double cutNDepNorm(int N, double Prb){return 1.41421356237*TMath::ErfInverse(1.-sqrt(2.*Prb/(N*N-N)));}
- inline double cutNDepNorm(int N, double Prb){return TMath::ChisquareQuantile(1.-sqrt(2.*Prb/(N*N-N)),2.);}
+ float median(std::vector<float> &Vec){
+ int N=Vec.size();
+ if(N==1) return Vec[0];
+ if(N>1){
+ std::vector<float> tmp(Vec);
+ std::sort(tmp.begin(),tmp.end());
+ return (tmp[(N-1)/2]+tmp[N/2])/2.;
+ }
+ return 0.;
+ }
+
+ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not coming from B,C-decays
+
+//inline double cutNDepNorm(int N, double Prb){return TMath::ChisquareQuantile(1.-sqrt(2.*Prb/(N*N-N)),2.);}
xAOD::Vertex* InDetVKalVxInJetTool::GetVrtSec(const std::vector<const Rec::TrackParticle*>& InpTrk,
const xAOD::Vertex & PrimVrt,
@@ -46,26 +59,25 @@ namespace InDet{
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "GetVrtSec() called with Rec::TrackParticle=" <<InpTrk.size()<< endmsg;
- std::vector<const Rec::TrackParticle*> SelectedTracks,SelectedTracksRelax;
- SelectedTracks.clear(); SelectedTracksRelax.clear();
+ std::vector<const Rec::TrackParticle*> SelectedTracks;
+ SelectedTracks.clear();
ListSecondTracks.clear();
Results.clear();
- m_NRefTrk=0;
+ m_NRefPVTrk=0;
if( InpTrk.size() < 2 ) { return 0;} // 0,1 track => nothing to do!
- int NPVParticle = SelGoodTrkParticle( InpTrk, PrimVrt, JetDir, SelectedTracks);
+ SelGoodTrkParticle( InpTrk, PrimVrt, JetDir, SelectedTracks);
long int NTracks = (int) (SelectedTracks.size());
- if(m_FillHist){m_hb_ntrkjet->Fill( (double) NTracks, m_w_1); }
+ if(m_fillHist){m_hb_ntrkjet->Fill( (double) NTracks, m_w_1); }
if( NTracks < 2 ) { return 0;} // 0,1 selected track => nothing to do!
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "Number of selected tracks inside jet= " <<NTracks << endmsg;
- SelGoodTrkParticleRelax( InpTrk, PrimVrt, JetDir, SelectedTracksRelax);
- TLorentzVector MomentumJet = TotalMom(GetPerigeeVector(SelectedTracksRelax));
- if(m_FillHist){m_hb_jmom->Fill( MomentumJet.E(), m_w_1); }
+ TLorentzVector MomentumJet = TotalMom(GetPerigeeVector(SelectedTracks));
+ if(m_fillHist){m_hb_jmom->Fill( MomentumJet.E(), m_w_1); }
//--------------------------------------------------------------------------------------------
@@ -93,11 +105,9 @@ namespace InDet{
double Vrt2TrackNumber = (double) ListSecondTracks.size()/2.;
RemoveDoubleEntries(ListSecondTracks);
AnalysisUtils::Sort::pT (&ListSecondTracks);
-//--Number of 2tr vertices where each track is used
- std::vector<int> combCount(ListSecondTracks.size());
- for(int tk=0;tk<(int)ListSecondTracks.size(); tk++){
- combCount[tk]=std::count(saveSecondTracks.begin(),saveSecondTracks.end(),ListSecondTracks[tk]);
- }
+//--Ranking of selected tracks
+ std::vector<float> trkRank(0);
+ for(auto tk : ListSecondTracks) trkRank.push_back( m_trackClassificator->trkTypeWgts(tk, PrimVrt, JetDir)[0] );
//---
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Found different tracks in pairs="<< ListSecondTracks.size()<<endmsg;
if(ListSecondTracks.size() < 2 ) { return 0;} // Less than 2 tracks left
@@ -114,7 +124,7 @@ namespace InDet{
double Signif3D, Signif3DP=0, Signif3DS=0;
std::vector<double> Impact,ImpactError;
- double Chi2 = FitCommonVrt( ListSecondTracks,combCount, xaodPrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
+ double Chi2 = FitCommonVrt( ListSecondTracks, trkRank, xaodPrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
//
if( Chi2 < 0) { return 0; } // Vertex not reconstructed
//
@@ -129,7 +139,7 @@ namespace InDet{
Signif3DS = m_fitSvc->VKalGetImpact(i_ntrk, FitVertex , 1, Impact, ImpactError);
if( Signif3DS > 10.) continue;
Signif3DP = m_fitSvc->VKalGetImpact(i_ntrk, PrimVrt.position(), 1, Impact, ImpactError);
- if(m_FillHist){ m_hb_diffPS->Fill( Signif3DP-Signif3DS, m_w_1); }
+ if(m_fillHist){ m_hb_diffPS->Fill( Signif3DP-Signif3DS, m_w_1); }
if(Signif3DP-Signif3DS>1.0) AdditionalTracks.push_back(i_ntrk);
}
}
@@ -137,8 +147,9 @@ namespace InDet{
// Add found tracks and refit
if( AdditionalTracks.size() > 0){
for (auto i_ntrk : AdditionalTracks) ListSecondTracks.push_back(i_ntrk);
- std::vector<int> tmpCount(ListSecondTracks.size(),1);
- Chi2 = FitCommonVrt( ListSecondTracks, tmpCount, xaodPrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
+ trkRank.clear();
+ for(auto tk : ListSecondTracks) trkRank.push_back( m_trackClassificator->trkTypeWgts(tk, PrimVrt, JetDir)[0] );
+ Chi2 = FitCommonVrt( ListSecondTracks, trkRank, xaodPrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
if( Chi2 < 0) { return 0; } // Vertex not reconstructed
}
//
@@ -158,16 +169,16 @@ namespace InDet{
tCnt++;
}
if(m_useVertexCleaning){
- if(!Check2TrVertexInPixel(ListSecondTracks[0],ListSecondTracks[1],FitVertex,Signif3D)) return 0;
- double xDif=FitVertex.x()-m_XlayerB, yDif=FitVertex.y()-m_YlayerB ;
+ if(!Check2TrVertexInPixel(ListSecondTracks[0],ListSecondTracks[1],FitVertex,ErrorMatrix)) return 0;
+ double xDif=FitVertex.x()-m_xLayerB, yDif=FitVertex.y()-m_yLayerB ;
double Dist2D=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2D < m_RlayerB-m_SVResolutionR) { // Inside B-layer
- if(m_FillHist){ if(Charge){m_hb_totmass2T1->Fill(Momentum.M(),m_w_1);}else{m_hb_totmass2T0->Fill(Momentum.M(),m_w_1);} }
- if(m_FillHist){ m_hb_blshared->Fill((float)BLshared,m_w_1); }
- //if(BLshared>m_CutSharedHits) return 0;
- } else if(Dist2D > m_RlayerB+m_SVResolutionR) { //Outside b-layer
- if(m_FillHist){ m_hb_pxshared->Fill((float)PXshared,m_w_1); }
- //if(PXshared>m_CutSharedHits) return 0;
+ if (Dist2D < m_rLayerB) { // Inside B-layer
+ if(m_fillHist){ if(Charge){m_hb_totmass2T1->Fill(Momentum.M(),m_w_1);}else{m_hb_totmass2T0->Fill(Momentum.M(),m_w_1);} }
+ if(m_fillHist){ m_hb_blshared->Fill((float)BLshared,m_w_1); }
+ //if(BLshared>m_cutSharedHits) return 0;
+ } else if(Dist2D > m_rLayerB) { //Outside b-layer
+ if(m_fillHist){ m_hb_pxshared->Fill((float)PXshared,m_w_1); }
+ //if(PXshared>m_cutSharedHits) return 0;
}
} // end of 2tr vertex cleaning code
//
@@ -180,15 +191,7 @@ namespace InDet{
}
return 0;
}
-//-- Protection against fake vertices far from interaction point
- if(NPVParticle<1)NPVParticle=1;
- double vvdist3D=VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
- double t3Dimp1= m_fitSvc->VKalGetImpact(ListSecondTracks[0], PrimVrt.position(), 1, Impact, ImpactError)/fabs(TrkAtVrt[0][2]);
- double t3Dimp2= m_fitSvc->VKalGetImpact(ListSecondTracks[1], PrimVrt.position(), 1, Impact, ImpactError)/fabs(TrkAtVrt[1][2]);
- double selVar=(t3Dimp1<t3Dimp2 ? t3Dimp1 : t3Dimp2)/sqrt((double)NPVParticle)/vvdist3D/500.;
- if(m_FillHist){ m_hb_tr2SelVar->Fill( selVar , m_w_1); }
- if(selVar<m_AntiFake2trVrtCut)return 0;
- if(m_FillHist){ m_hb_totmass2T2->Fill(Momentum.M(),m_w_1); }
+ if(m_fillHist){ m_hb_totmass2T2->Fill(Momentum.M(),m_w_1); }
}
@@ -205,16 +208,16 @@ namespace InDet{
double xvt=FitVertex.x(); double yvt=FitVertex.y();
- double Dist2DBP=sqrt( (xvt-m_Xbeampipe)*(xvt-m_Xbeampipe) + (yvt-m_Ybeampipe)*(yvt-m_Ybeampipe) );
- double Dist2DBL=sqrt( (xvt-m_XlayerB)*(xvt-m_XlayerB) + (yvt-m_YlayerB)*(yvt-m_YlayerB) );
- double Dist2DL1=sqrt( (xvt-m_Xlayer1)*(xvt-m_Xlayer1) + (yvt-m_Ylayer1)*(yvt-m_Ylayer1) );
- double Dist2DL2=sqrt( (xvt-m_Xlayer2)*(xvt-m_Xlayer2) + (yvt-m_Ylayer2)*(yvt-m_Ylayer2) );
+ double Dist2DBP=sqrt( (xvt-m_beampipeX)*(xvt-m_beampipeX) + (yvt-m_beampipeY)*(yvt-m_beampipeY) );
+ double Dist2DBL=sqrt( (xvt-m_xLayerB)*(xvt-m_xLayerB) + (yvt-m_yLayerB)*(yvt-m_yLayerB) );
+ double Dist2DL1=sqrt( (xvt-m_xLayer1)*(xvt-m_xLayer1) + (yvt-m_yLayer1)*(yvt-m_yLayer1) );
+ double Dist2DL2=sqrt( (xvt-m_xLayer2)*(xvt-m_xLayer2) + (yvt-m_yLayer2)*(yvt-m_yLayer2) );
double minDstMat=39.9;
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DBP-m_Rbeampipe));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DBL-m_RlayerB));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DL1-m_Rlayer1));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_Rlayer2));
- if(m_existIBL) minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_Rlayer3)); // 4-layer pixel detector
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DBP-m_beampipeR));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DBL-m_rLayerB));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DL1-m_rLayer1));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_rLayer2));
+ if(m_existIBL) minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_rLayer3)); // 4-layer pixel detector
VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
if(JetVrtDir < 0) Signif3D = -Signif3D;
@@ -222,10 +225,7 @@ namespace InDet{
Amg::Vector3D DirForPt( FitVertex.x()-PrimVrt.x(),
FitVertex.y()-PrimVrt.y(),
FitVertex.z()-PrimVrt.z());
- if( m_MassType == 3 ) Results.push_back( TotalTVMom(DirForPt, GetPerigeeVector(ListSecondTracks)));
- if( m_MassType == 2 ) Results.push_back(TotalTMom(GetPerigeeVector(ListSecondTracks))*1.15);
- if( m_MassType == 1 ) Results.push_back(Momentum.M()); //1st
-
+ Results.push_back(Momentum.M()); //1st
double eRatio = Momentum.E()/MomentumJet.E();
Results.push_back( eRatio<0.99999 ? eRatio : 0.99999); //2nd
Results.push_back(Vrt2TrackNumber); //3rd
@@ -239,7 +239,7 @@ namespace InDet{
Results.push_back((Momentum.M()-2.*m_massPi)*eRatio/m_massB); //10th "Product" variable
Results.push_back((Momentum.Pt()/Momentum.M())*(m_massB/JetDir.Pt()) ); //11th "Boost" variable
- if(m_FillHist){
+ if(m_fillHist){
if(ListSecondTracks.size()==2) m_hb_r2dc->Fill( FitVertex.perp(), m_w_1);
else m_hb_rNdc->Fill( FitVertex.perp(), m_w_1);
m_hb_mom->Fill( MomentumJet.E(), m_w_1);
@@ -254,7 +254,7 @@ namespace InDet{
if(trackPt>trackPtMax)trackPtMax=trackPt;
}
m_hb_trkPtMax->Fill( trackPtMax, m_w_1);
- m_pr_effVrt->Fill((float)m_NRefTrk,1.);
+ m_pr_effVrt->Fill((float)m_NRefPVTrk,1.);
m_pr_effVrtEta->Fill( JetDir.Eta(),1.);
}
@@ -281,133 +281,6 @@ namespace InDet{
- xAOD::Vertex* InDetVKalVxInJetTool::tryPseudoVertex(const std::vector<const xAOD::TrackParticle*>& SelectedTracks,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- const TLorentzVector & TrkJet,
- const int & nTrkLead,
- std::vector<double> & Results)
- const
- {
-//---------------First try jet axis+track
- if((int)SelectedTracks.size()<nTrkLead)return 0;
- //-----------------------------------------
- TLorentzVector sumSelTrk(0.,0.,0.,0.); //Tracks coming from any SV
-
- Amg::Vector3D FitVertex;
- std::vector<double> ErrorMatrix;
- std::vector< std::vector<double> > TrkAtVrt;
- TLorentzVector Momentum;
- std::vector<double> Chi2PerTrk;
- long int tChrg=0;
- double Chi2=0.;
- std::vector<const xAOD::TrackParticle*> tTrkForFit(2,0);
- std::vector<float> tmpCov(15,0.); tmpCov[0]=1e-4; tmpCov[2]=4.e-4; tmpCov[5]=4.e-4; tmpCov[9]=4.e-4; tmpCov[14]=1.e-10;
- StatusCode scode; scode.setChecked();
- std::vector<const xAOD::TrackParticle*> reducedTrkSet(SelectedTracks.begin(),SelectedTracks.begin()+nTrkLead);
- double maxImp=RemoveNegImpact(reducedTrkSet,PrimVrt,JetDir,m_pseudoSigCut);
- if(reducedTrkSet.size()==0) return 0;
- if(reducedTrkSet.size()==1 && maxImp<m_pseudoSigCut+1.)return 0;
- if(m_FillHist)m_hb_rawVrtN->Fill(reducedTrkSet.size(),1.);
-//
-//------
- int sel1T=-1; double selDST=0.;
- if(reducedTrkSet.size()>0){
- m_fitSvc->setApproximateVertex(PrimVrt.x(), PrimVrt.y(), PrimVrt.z());
- xAOD::TrackParticle *tmpBTP=new xAOD::TrackParticle(); tmpBTP->makePrivateStore();
- tmpBTP->setDefiningParameters(0., 0., JetDir.Phi(), JetDir.Theta(), sin(JetDir.Theta())/2.e5); // Pt=200GeV track
- tmpBTP->setParametersOrigin(PrimVrt.x(),PrimVrt.y(),PrimVrt.z());
- tmpBTP->setDefiningParametersCovMatrixVec(tmpCov);
- tTrkForFit[1]=tmpBTP;
- TLorentzVector sumB(0.,0.,0.,0.);
- int nvFitted=0;
- int nPRT=reducedTrkSet.size();
- for(int it=0; it<nPRT; it++){
- tTrkForFit[0]=reducedTrkSet[it];
- if(tTrkForFit[0]->pt()<2000.)continue;
- if(nPRT==1 && tTrkForFit[0]->pt()<300.*log(JetDir.Pt()))continue;
- m_fitSvc->setApproximateVertex(0., 0., 0.);
- scode=VKalVrtFitBase( tTrkForFit, FitVertex, Momentum, tChrg, ErrorMatrix, Chi2PerTrk, TrkAtVrt, Chi2); nvFitted++;
- if(scode.isFailure() || Chi2>6.)continue;
- if(FitVertex.perp()>reducedTrkSet[it]->radiusOfFirstHit())continue;
- //double dSVPV=ProjPosT(FitVertex-PrimVrt.position(),JetDir);
- double Signif3D; VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
- if(FitVertex.perp()>m_Rbeampipe && Signif3D<2.) continue; // Cleaning of material interactions
- if(m_FillHist)m_hb_DST_JetTrkSV->Fill(Signif3D,1.);
- if(tTrkForFit[0]->pt()<2.e4){
- if(!Check2TrVertexInPixel(tTrkForFit[0],tTrkForFit[0],FitVertex,Signif3D)) continue;
- }
- //sumSelTrk += MomAtVrt(TrkAtVrt[0]) ; sumB += MomAtVrt(TrkAtVrt[1]) ;
- if(Signif3D>selDST){ sel1T=it; selDST=Signif3D; sumSelTrk=MomAtVrt(TrkAtVrt[0]); sumB=MomAtVrt(TrkAtVrt[1]);}
- }
- if(sumSelTrk.Pt()>0. && sel1T>=0 ){
- Results.resize(4);
- double pt1=sumSelTrk.Pt(sumB.Vect()); double E1=sqrt(pt1*pt1+sumSelTrk.M2());
- Results[0]=pt1+E1; //Invariant mass
- Results[1]=sumSelTrk.Pt()/TrkJet.Pt(); //Ratio
- Results[2]=0.; //Should be
- Results[3]=reducedTrkSet.size(); //Found leading tracks with high positive impact
- tTrkForFit[0]=reducedTrkSet[sel1T]; //------------------------- Refit selected vertex for xAOD::Vertex
- if(nvFitted>1){ //If only one fit attempt made - don't need to refit
- scode=VKalVrtFitBase( tTrkForFit, FitVertex, Momentum, tChrg, ErrorMatrix, Chi2PerTrk, TrkAtVrt, Chi2);
- if(scode.isFailure()){sumSelTrk.SetXYZT(0.,0.,0.,0.); Results.clear();}
- }
- }
- delete tmpBTP;
- }
-//
-//------ Plane-Plane crossing doesn't provide good vertices for the moment
-// int sel2TI=-1,sel2TJ=-1;
-// if(reducedTrkSet.size()<1 && sel1T<0){
-// int nPRT=reducedTrkSet.size(); Amg::Vector3D VB1,VB2; float distMax=1.e10;
-// for(int it=0; it<nPRT-1; it++){ for(int jt=it+1; jt<nPRT; jt++){
-// TLorentzVector pseudoB=GetBDir( reducedTrkSet[it], reducedTrkSet[jt], PrimVrt, VB1, VB2);
-// if( VB1.dot(VB2)==0. || pseudoB.DeltaR(JetDir)>0.3 ) continue;
-// if(Amg::distance(VB1,VB2)<distMax){ sel2TI=it; sel2TJ=jt; distMax=Amg::distance(VB1,VB2);}
-// } }
-// if(sel2TI>=0){
-// tTrkForFit[0]=reducedTrkSet[sel2TI]; tTrkForFit[1]=reducedTrkSet[sel2TJ];
-// //float RMHIT=TMath::Min(tTrkForFit[0]->radiusOfFirstHit(),tTrkForFit[1]->radiusOfFirstHit()); // Closest hit on both tracks
-// m_fitSvc->setApproximateVertex(0., 0., 0.);
-// scode=VKalVrtFitBase( tTrkForFit, FitVertex, Momentum, tChrg, ErrorMatrix, Chi2PerTrk, TrkAtVrt, Chi2);
-// if(scode.isSuccess() && ProjPosT(FitVertex-PrimVrt.position(),JetDir)>0 && FitVertex.perp()<180.){
-// sumSelTrk += MomAtVrt(TrkAtVrt[0]) ; sumSelTrk += MomAtVrt(TrkAtVrt[1]) ;
-// Results.resize(4);
-// Results[0]=Momentum.M(); //Invariant mass
-// Results[1]=Momentum.Pt()/TrkJet.Pt(); //Ratio
-// Results[2]=0.; //Should be
-// Results[3]=nPRT; //Found leading tracks with high positive impact
-// } }
-// }
-//
-//------
- if(sumSelTrk.Pt()==0)return 0; //---------- Nothing found. Damn it! Else - return xAOD::Vertex
- Results.resize(7);
- Results[4]=0.; Results[5]=0.;
- Results[6]=TrkJet.E();
- xAOD::Vertex * tmpVertex=new xAOD::Vertex();
- tmpVertex->makePrivateStore();
- tmpVertex->setPosition(FitVertex);
- std::vector<float> floatErrMtx; floatErrMtx.resize(ErrorMatrix.size());
- for(int i=0; i<(int)ErrorMatrix.size(); i++) floatErrMtx[i]=ErrorMatrix[i];
- tmpVertex->setCovariance(floatErrMtx);
- tmpVertex->setFitQuality(Chi2, (float)(tTrkForFit.size()*2.-3.));
- std::vector<Trk::VxTrackAtVertex> & tmpVTAV=tmpVertex->vxTrackAtVertex(); tmpVTAV.clear();
- AmgSymMatrix(5) *CovMtxP=new AmgSymMatrix(5); (*CovMtxP).setIdentity();
- Trk::Perigee * tmpMeasPer = new Trk::Perigee( 0.,0., TrkAtVrt[0][0], TrkAtVrt[0][1], TrkAtVrt[0][2],
- Trk::PerigeeSurface(FitVertex), CovMtxP );
- tmpVTAV.push_back( Trk::VxTrackAtVertex( 1., tmpMeasPer) );
- ElementLink<xAOD::TrackParticleContainer> TEL; TEL.setElement( tTrkForFit[0] );
- const xAOD::TrackParticleContainer* cont = (const xAOD::TrackParticleContainer* ) (tTrkForFit[0]->container() );
- TEL.setStorableObject(*cont);
- tmpVertex->addTrackAtVertex(TEL,1.);
- if(m_FillHist){m_hb_massJetTrkSV ->Fill(Results[0],1.);
- m_hb_ratioJetTrkSV->Fill(Results[1],1.);
- m_hb_NImpJetTrkSV ->Fill(Results[3],1.);}
- return tmpVertex;
- }
-
-
@@ -428,40 +301,33 @@ namespace InDet{
TLorentzVector Momentum;
double Signif3D=0., Chi2=0.;
- std::vector<const xAOD::TrackParticle*> SelectedTracks(0),SelectedTracksRelax(0);
+ std::vector<const xAOD::TrackParticle*> SelectedTracks(0);
Results.clear();
ListSecondTracks.clear();
- m_NRefTrk=0;
+ m_NRefPVTrk=0;
if( InpTrk.size() < 2 ) { return 0;} // 0,1 track => nothing to do!
- int NPVParticle = SelGoodTrkParticle( InpTrk, PrimVrt, JetDir, SelectedTracks);
- while(SelectedTracks.size() && SelectedTracks[0]->pt()/JetDir.Pt()>1.)SelectedTracks.erase(SelectedTracks.begin());
- if((int)SelectedTracks.size()>m_TrackInJetNumberLimit){
- SelectedTracks.resize(m_TrackInJetNumberLimit); // SelectedTracks are ordered in pT
- }
- while( SelectedTracks.size()>4 && medianPtF(SelectedTracks)/JetDir.Pt()<0.01) SelectedTracks.pop_back();
-
+ SelGoodTrkParticle( InpTrk, PrimVrt, JetDir, SelectedTracks);
+ if(m_fillHist){m_hb_ntrkjet->Fill( (double)SelectedTracks.size(), m_w_1);
+ m_pr_NSelTrkMean->Fill(JetDir.Pt(),(double)SelectedTracks.size()); }
long int NTracks = (int) (SelectedTracks.size());
- if(m_FillHist){m_hb_ntrkjet->Fill( (double) NTracks, m_w_1); }
if( NTracks < 2 ) { return 0;} // 0,1 selected track => nothing to do!
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "Number of selected tracks inside jet= " <<NTracks << endmsg;
- SelGoodTrkParticleRelax( InpTrk, PrimVrt, JetDir, SelectedTracksRelax);
- TLorentzVector MomentumJet = TotalMom(SelectedTracksRelax);
- if(m_FillHist){m_hb_jmom->Fill( MomentumJet.E(), m_w_1); }
+ TLorentzVector MomentumJet = TotalMom(SelectedTracks);
+ if(m_fillHist){m_hb_jmom->Fill( MomentumJet.E(), m_w_1); }
//--------------------------------------------------------------------------------------------
// Initial xAOD::TrackParticle list ready
float Vrt2TrackNumber =0;
- const int nTrkLead=5;
-
std::vector<const xAOD::TrackParticle*> TracksForFit;
//std::vector<double> InpMass; for(int i=0; i<NTracks; i++){InpMass.push_back(m_massPi);}
std::vector<double> InpMass(NTracks,m_massPi);
- Select2TrVrt(SelectedTracks, TracksForFit, PrimVrt, JetDir, InpMass, TrkFromV0, ListSecondTracks);
+ Select2TrVrt(SelectedTracks, TracksForFit, PrimVrt, JetDir, InpMass, TrkFromV0,
+ ListSecondTracks);
m_WorkArray->m_Incomp.clear(); // Not needed for single vertex version
//
//--- Cleaning
@@ -478,15 +344,14 @@ namespace InDet{
if(itf!=SelectedTracks.end()) SelectedTracks.erase(itf);}
//---
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Found different xAOD tracks in pairs="<< ListSecondTracks.size()<<endmsg;
- if(ListSecondTracks.size() < 2 ) return tryPseudoVertex( SelectedTracks, PrimVrt, JetDir, MomentumJet, nTrkLead, Results);
-
-//--- At least 2 secondary tracks are found.
-
- AnalysisUtils::Sort::pT (&ListSecondTracks);
-//--Number of 2tr vertices where each track is used
- std::vector<int> combCount(ListSecondTracks.size(),0);
- for(int tk=0;tk<(int)ListSecondTracks.size(); tk++){
- combCount[tk]=std::count(saveSecondTracks.begin(),saveSecondTracks.end(),ListSecondTracks[tk]);
+ if(ListSecondTracks.size() < 2 ) return 0;
+
+//--Ranking of selected tracks
+ std::vector<float> trkRank(0);
+ for(auto tk : ListSecondTracks) trkRank.push_back( m_trackClassificator->trkTypeWgts(tk, PrimVrt, JetDir)[0] );
+ while( median(trkRank)<0.3 && trkRank.size()>3 ) {
+ int Smallest= std::min_element(trkRank.begin(),trkRank.end()) - trkRank.begin();
+ RemoveEntryInList(ListSecondTracks,trkRank,Smallest);
}
//
//-----------------------------------------------------------------------------------------------------
@@ -495,57 +360,55 @@ namespace InDet{
//
double Signif3DP=0, Signif3DS=0;
- Chi2 = FitCommonVrt( ListSecondTracks, combCount, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
+ Chi2 = FitCommonVrt( ListSecondTracks, trkRank, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
if( Chi2 < 0 && ListSecondTracks.size()>2 ) { // Vertex not reconstructed. Try to remove one track with biggest pt.
double tpmax=0; int ipmax=-1;
for(int it=0; it<(int)ListSecondTracks.size(); it++) if(tpmax<ListSecondTracks[it]->pt()){tpmax=ListSecondTracks[it]->pt(); ipmax=it;}
- if(ipmax>=0)RemoveEntryInList(ListSecondTracks,combCount,ipmax);
- Chi2 = FitCommonVrt( ListSecondTracks, combCount, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
- if( Chi2 < 0 && ListSecondTracks.size()>2 ) { // Vertex not reconstructed. Try to remove another track with biggest pt.
- tpmax=0.; ipmax=-1;
- for(int it=0; it<(int)ListSecondTracks.size(); it++) if(tpmax<ListSecondTracks[it]->pt()){tpmax=ListSecondTracks[it]->pt(); ipmax=it;}
- if(ipmax>=0)RemoveEntryInList(ListSecondTracks,combCount,ipmax);
- Chi2 = FitCommonVrt( ListSecondTracks, combCount, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
- }
+ if(ipmax>=0)RemoveEntryInList(ListSecondTracks,trkRank,ipmax);
+ Chi2 = FitCommonVrt( ListSecondTracks, trkRank, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Second FitCommonVrt try="<< Chi2<<" Ntrk="<<ListSecondTracks.size()<<endmsg;
}
- if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" FitCommonVrt result="<< Chi2<<endmsg;
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" FitCommonVrt result="<< Chi2<<" Ntrk="<<ListSecondTracks.size()<<endmsg;
//
- if( Chi2 < 0) return tryPseudoVertex( SelectedTracks, PrimVrt, JetDir, MomentumJet, nTrkLead, Results);
+ if( Chi2 < 0) return 0;
//
// Check jet tracks not in secondary vertex
std::map<double,const xAOD::TrackParticle*> AdditionalTracks;
VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
if(Signif3D>8.){
- for (auto i_ntrk : SelectedTracks) {
- std::vector<const xAOD::TrackParticle*>::const_iterator i_found =
- find( ListSecondTracks.begin(), ListSecondTracks.end(), i_ntrk);
- if( i_found != ListSecondTracks.end() ) continue;
- if(i_ntrk->pt()<m_JetPtFractionCut*JetDir.Perp())continue;
- if(!Check1TrVertexInPixel(i_ntrk,FitVertex)) continue;
- Signif3DS = m_fitSvc->VKalGetImpact(i_ntrk, FitVertex , 1, Impact, ImpactError);
- if( Signif3DS > 10.) continue;
- if(i_ntrk->radiusOfFirstHit()>60 && FitVertex.perp()<m_RlayerB-m_SVResolutionR)continue; //VK no hit in IBL and BL
- Signif3DP = m_fitSvc->VKalGetImpact(i_ntrk, PrimVrt.position(), 1, Impact, ImpactError);
- if(m_FillHist){ m_hb_diffPS->Fill( Signif3DP-Signif3DS, m_w_1); }
- if(Signif3DP-Signif3DS>-1.0) AdditionalTracks[Signif3DP-Signif3DS]=i_ntrk;
- }
+ int hitL1=0, nLays=0, hitIBL=0, hitBL=0;
+ for (auto i_ntrk : SelectedTracks) {
+ if( find( ListSecondTracks.begin(), ListSecondTracks.end(), i_ntrk) != ListSecondTracks.end() ) continue; // Track is used already
+ std::vector<float> trkScore=m_trackClassificator->trkTypeWgts(i_ntrk, PrimVrt, JetDir);
+ if( trkScore[0] < m_cutHFClass/2.) continue;
+ Signif3DS = m_fitSvc->VKalGetImpact(i_ntrk, FitVertex , 1, Impact, ImpactError);
+ if( Signif3DS > 10.) continue;
+ getPixelLayers(i_ntrk , hitIBL , hitBL, hitL1, nLays );
+ if( hitIBL<=0 && hitBL<=0 ) continue; // No IBL and BL pixel hits => non-precise track
+ Signif3DP = m_fitSvc->VKalGetImpact(i_ntrk, PrimVrt.position(), 1, Impact, ImpactError);
+ if(Signif3DP<1.)continue;
+ if(m_fillHist){ m_hb_diffPS->Fill( Signif3DP-Signif3DS, m_w_1); }
+ if(Signif3DP-Signif3DS>4.0) AdditionalTracks[Signif3DP-Signif3DS]=i_ntrk;
+ }
}
//
// Add found tracks and refit
+//
if( AdditionalTracks.size() > 0){
-while (AdditionalTracks.size()>3) AdditionalTracks.erase(AdditionalTracks.begin());//Tracks are in increasing DIFF order.
-for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second); //3tracks with max DIFF are selected
- std::vector<int> tmpCount(ListSecondTracks.size(),1);
- Chi2 = FitCommonVrt( ListSecondTracks, tmpCount, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
+ while (AdditionalTracks.size()>3) AdditionalTracks.erase(AdditionalTracks.begin());//Tracks are in increasing DIFF order.
+ for (auto atrk : AdditionalTracks) ListSecondTracks.push_back(atrk.second); //3tracks with max DIFF are selected
+ trkRank.clear();
+ for(auto tk : ListSecondTracks) trkRank.push_back( m_trackClassificator->trkTypeWgts(tk, PrimVrt, JetDir)[0] );
+ Chi2 = FitCommonVrt( ListSecondTracks, trkRank, PrimVrt, JetDir, InpMass, FitVertex, ErrorMatrix, Momentum, TrkAtVrt);
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Added track FitCommonVrt output="<< Chi2<<endmsg;
- if( Chi2 < 0) return tryPseudoVertex( SelectedTracks, PrimVrt, JetDir, MomentumJet, nTrkLead, Results);
+ if( Chi2 < 0) return 0;
}
//
// Saving of results
//
-//
//
if( ListSecondTracks.size()==2 ){ // If there are 2 only tracks
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Start Ntr=2 vertex check"<<endmsg;
int Charge=0;
uint8_t BLshared=0;
uint8_t PXshared=0;
@@ -555,22 +418,18 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
if( i_ntrk->summaryValue( retval, xAOD::numberOfPixelSharedHits) ) PXshared += retval;
if( i_ntrk->summaryValue( retval, xAOD::numberOfInnermostPixelLayerSharedHits) ) BLshared += retval;
}
-
- double vvdist3D=VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
+ VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
if(m_useVertexCleaning){
- if(!Check2TrVertexInPixel(ListSecondTracks[0],ListSecondTracks[1],FitVertex,Signif3D)) return 0;
- double xDif=FitVertex.x()-m_XlayerB, yDif=FitVertex.y()-m_YlayerB ;
- double Dist2D=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2D < m_RlayerB-m_SVResolutionR) { // Inside B-layer
- if(m_FillHist){ if(Charge){m_hb_totmass2T1->Fill(Momentum.M(),m_w_1);}else{m_hb_totmass2T0->Fill(Momentum.M(),m_w_1);} }
- if(m_FillHist){ m_hb_blshared->Fill((float)BLshared,m_w_1); }
- //if(BLshared>m_CutSharedHits) return 0; //VK Kills more b-jets events
- }else if(Dist2D > m_RlayerB+m_SVResolutionR) { //Outside b-layer
- if(m_FillHist){ m_hb_pxshared->Fill((float)PXshared,m_w_1); }
- //if(PXshared>m_CutSharedHits) return 0;
- }
+ if(!Check2TrVertexInPixel(ListSecondTracks[0],ListSecondTracks[1],FitVertex,ErrorMatrix)) return 0;
+ if(m_fillHist){
+ double xDif=FitVertex.x()-m_xLayerB, yDif=FitVertex.y()-m_yLayerB ;
+ double Dist2D=sqrt(xDif*xDif+yDif*yDif);
+ if (Dist2D < m_rLayerB-VrtRadiusError(FitVertex,ErrorMatrix)) m_hb_blshared->Fill((float)BLshared,m_w_1);
+ else if(Dist2D > m_rLayerB+VrtRadiusError(FitVertex,ErrorMatrix)) m_hb_pxshared->Fill((float)PXshared,m_w_1);
+ }
} //end 2tr vertex cleaning code
//
+ if(m_fillHist){ if(Charge){m_hb_totmass2T1->Fill(Momentum.M(),m_w_1);}else{m_hb_totmass2T0->Fill(Momentum.M(),m_w_1);} }
if( !Charge && fabs(Momentum.M()-m_massK0)<15. ) { // Final rejection of K0
TrkFromV0.push_back(ListSecondTracks[0]);
TrkFromV0.push_back(ListSecondTracks[1]);
@@ -580,21 +439,12 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
}
return 0;
}
-//-- Protection against fake vertices far from interaction point
- if(NPVParticle<1)NPVParticle=1;
- double t3Dimp1= m_fitSvc->VKalGetImpact(ListSecondTracks[0], PrimVrt.position(), 1, Impact, ImpactError)/fabs(TrkAtVrt[0][2]);
- double t3Dimp2= m_fitSvc->VKalGetImpact(ListSecondTracks[1], PrimVrt.position(), 1, Impact, ImpactError)/fabs(TrkAtVrt[1][2]);
- double selVar=(t3Dimp1<t3Dimp2 ? t3Dimp1 : t3Dimp2)/sqrt((double)NPVParticle)/vvdist3D/500.;
- if(m_FillHist){ m_hb_tr2SelVar->Fill( selVar , m_w_1); }
- if(selVar<m_AntiFake2trVrtCut)return 0;
- if(m_FillHist){ m_hb_totmass2T2->Fill(Momentum.M(),m_w_1); }
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Ntr=2 vertex check passed"<<endmsg;
}
- double JetVrtDir =
- JetDir.Px()*(FitVertex.x()-PrimVrt.x())
- + JetDir.Py()*(FitVertex.y()-PrimVrt.y())
- + JetDir.Pz()*(FitVertex.z()-PrimVrt.z());
+ double JetVrtDir = ProjSV_PV(FitVertex,PrimVrt,JetDir);
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<"Combined SV neg.dir="<<JetVrtDir<<endmsg;
if( m_getNegativeTag )
{ if( JetVrtDir>0. ) return 0; }
else if( m_getNegativeTail )
@@ -603,30 +453,26 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
{ if( JetVrtDir<0. ) return 0; }
double xvt=FitVertex.x(); double yvt=FitVertex.y();
- double Dist2DBP=sqrt( (xvt-m_Xbeampipe)*(xvt-m_Xbeampipe) + (yvt-m_Ybeampipe)*(yvt-m_Ybeampipe) );
- double Dist2DBL=sqrt( (xvt-m_XlayerB)*(xvt-m_XlayerB) + (yvt-m_YlayerB)*(yvt-m_YlayerB) );
- double Dist2DL1=sqrt( (xvt-m_Xlayer1)*(xvt-m_Xlayer1) + (yvt-m_Ylayer1)*(yvt-m_Ylayer1) );
- double Dist2DL2=sqrt( (xvt-m_Xlayer2)*(xvt-m_Xlayer2) + (yvt-m_Ylayer2)*(yvt-m_Ylayer2) );
+ double Dist2DBP=sqrt( (xvt-m_beampipeX)*(xvt-m_beampipeX) + (yvt-m_beampipeY)*(yvt-m_beampipeY) );
+ double Dist2DBL=sqrt( (xvt-m_xLayerB)*(xvt-m_xLayerB) + (yvt-m_yLayerB)*(yvt-m_yLayerB) );
+ double Dist2DL1=sqrt( (xvt-m_xLayer1)*(xvt-m_xLayer1) + (yvt-m_yLayer1)*(yvt-m_yLayer1) );
+ double Dist2DL2=sqrt( (xvt-m_xLayer2)*(xvt-m_xLayer2) + (yvt-m_yLayer2)*(yvt-m_yLayer2) );
double minDstMat=39.9;
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DBP-m_Rbeampipe));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DBL-m_RlayerB));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DL1-m_Rlayer1));
- minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_Rlayer2));
- if(m_existIBL) minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_Rlayer3)); // 4-layer pixel detector
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DBP-m_beampipeR));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DBL-m_rLayerB));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DL1-m_rLayer1));
+ minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_rLayer2));
+ if(m_existIBL) minDstMat=TMath::Min(minDstMat,fabs(Dist2DL2-m_rLayer3)); // 4-layer pixel detector
VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
if(JetVrtDir < 0) Signif3D = -Signif3D;
- if(FitVertex.perp()>m_Rbeampipe && Signif3D<20.) return 0; // Final cleaning of material interactions
Amg::Vector3D DirForPt( FitVertex.x()-PrimVrt.x(),
FitVertex.y()-PrimVrt.y(),
FitVertex.z()-PrimVrt.z());
- //if( m_MassType == 3 ) Results.push_back( TotalTVMom(DirForPt, GetPerigeeVector(ListSecondTracks)));
- //if( m_MassType == 2 ) Results.push_back(TotalTMom(GetPerigeeVector(ListSecondTracks))*1.15);
- if( m_MassType == 1 ) Results.push_back(Momentum.M()); //1st
-
+ Results.push_back(Momentum.M()); //1st
double eRatio = Momentum.E()/MomentumJet.E();
Results.push_back( eRatio<0.99999 ? eRatio : 0.99999); //2nd
Results.push_back(Vrt2TrackNumber); //3rd
@@ -640,18 +486,17 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
Results.push_back((Momentum.M()-2.*m_massPi)*eRatio/m_massB); //10th "Product" variable
Results.push_back((Momentum.Pt()/Momentum.M())*(m_massB/JetDir.Pt()) ); //11th "Boost" variable
- if(m_FillHist){
+ if(m_fillHist){
// Find highest track Pt with respect to jet direction
double trackPt, trackPtMax=0.;
for (int tr=0; tr<(int)ListSecondTracks.size(); tr++) {
+ if(ListSecondTracks[tr]->pt()/JetDir.Pt()>0.5)continue;
trackPt=pTvsDir(Amg::Vector3D(JetDir.X(),JetDir.Y(),JetDir.Z()) , TrkAtVrt[tr]);
if(trackPt>trackPtMax)trackPtMax=trackPt;
}
- if(ListSecondTracks.size()==2) m_hb_r2dc->Fill( FitVertex.perp(), m_w_1);
- else if(ListSecondTracks.size()==3) m_hb_r3dc->Fill( FitVertex.perp(), m_w_1);
- else m_hb_rNdc->Fill( FitVertex.perp(), m_w_1);
+ m_hb_rNdc->Fill( FitVertex.perp(), m_w_1);
m_hb_trkPtMax->Fill( trackPtMax, m_w_1);
- m_pr_effVrt->Fill((float)m_NRefTrk,1.);
+ m_pr_effVrt->Fill((float)m_NRefPVTrk,1.);
m_pr_effVrtEta->Fill( JetDir.Eta(),1.);
m_hb_mom->Fill( MomentumJet.E(), m_w_1);
m_hb_ratio->Fill( Results[1], m_w_1);
@@ -693,14 +538,6 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
-
-
-
-
-
-
-
-
//
//--------------------------------------------------------
// Template routine for global secondary vertex fitting
@@ -708,7 +545,7 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
template <class Track>
double InDetVKalVxInJetTool::FitCommonVrt(std::vector<const Track*>& ListSecondTracks,
- std::vector<int> & cntComb,
+ std::vector<float> & trkRank,
const xAOD::Vertex & PrimVrt,
const TLorentzVector & JetDir,
std::vector<double> & InpMass,
@@ -722,7 +559,6 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
//preparation
StatusCode sc;
std::vector<double> Chi2PerTrk;
- const double maxRecMASS=6000.;
long int Charge;
double Chi2 = 0.;
Amg::Vector3D tmpVertex;
@@ -735,37 +571,29 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
m_fitSvc->setMassInputParticles( InpMass ); // Use pions masses
m_fitSvc->setMomCovCalc(1); /* Total momentum and its covariance matrix are calculated*/
sc=VKalVrtFitFastBase(ListSecondTracks,FitVertex); /* Fast crude estimation */
- if(sc.isFailure() || FitVertex.perp() > m_Rlayer2*2. ) { /* No initial estimation */
+ if(sc.isFailure() || FitVertex.perp() > m_rLayer2*2. ) { /* No initial estimation */
m_fitSvc->setApproximateVertex(PrimVrt.x(), PrimVrt.y(), PrimVrt.z()); /* Use as starting point */
} else {
m_fitSvc->setApproximateVertex(FitVertex.x(),FitVertex.y(),FitVertex.z()); /*Use as starting point*/
}
-// m_fitSvc-> setVertexForConstraint(PrimVrt.x(),PrimVrt.y(),PrimVrt.z());
-// m_fitSvc->setCovVrtForConstraint(PrimVrt.errorPosition().covValue(Trk::x,Trk::x),
-// PrimVrt.errorPosition().covValue(Trk::y,Trk::x),
-// PrimVrt.errorPosition().covValue(Trk::y,Trk::y),
-// PrimVrt.errorPosition().covValue(Trk::z,Trk::x),
-// PrimVrt.errorPosition().covValue(Trk::z,Trk::y),
-// PrimVrt.errorPosition().covValue(Trk::z,Trk::z) );
-// m_fitSvc->setCnstType(7);
- if(m_RobustFit)m_fitSvc->setRobustness(m_RobustFit);
- else m_fitSvc->setRobustness(0);
//fit itself
int NTracksVrt = ListSecondTracks.size(); double FitProb=0.;
-
+ std::vector<double> trkFitWgt(0);
for (i=0; i < NTracksVrt-1; i++) {
-// sc=m_fitSvc->VKalVrtFit(ListSecondTracks,FitVertex, Momentum,Charge,
-// ErrorMatrix,Chi2PerTrk,TrkAtVrt,Chi2);
- sc=VKalVrtFitBase(ListSecondTracks,FitVertex, Momentum,Charge,
+ if(m_RobustFit)m_fitSvc->setRobustness(m_RobustFit);
+ else m_fitSvc->setRobustness(0);
+ sc=VKalVrtFitBase(ListSecondTracks,FitVertex,Momentum,Charge,
ErrorMatrix,Chi2PerTrk,TrkAtVrt,Chi2);
if(sc.isFailure() || Chi2 > 1000000. ) { return -10000.;} // No fit
- //for( int mmt=0; mmt<(int) Chi2PerTrk.size(); mmt++) Chi2PerTrk[mmt] -= TMath::Min(trkSigPV[mmt],9.);
- Outlier = FindMax( Chi2PerTrk, cntComb );
+ sc=GetTrkFitWeights(trkFitWgt);
+ if(sc.isFailure()){ return -10000.;} // No weights
+ Outlier=std::min_element(trkFitWgt.begin(),trkFitWgt.end())-trkFitWgt.begin();
+ //////Outlier = FindMax( Chi2PerTrk, trkRank );
FitProb=TMath::Prob( Chi2, 2*ListSecondTracks.size()-3);
if(ListSecondTracks.size() == 2 ) break; // Only 2 tracks left
//////////////////////////////
double Signif3Dproj=VrtVrtDist( PrimVrt, FitVertex, ErrorMatrix, JetDir);
- if(Signif3Dproj<0 && (!m_getNegativeTail) && (!m_getNegativeTag) && FitProb < 0.1){
+ if(Signif3Dproj<0 && (!m_getNegativeTail) && (!m_getNegativeTag)){
double maxDst=-1.e12; int maxT=-1; double minChi2=1.e12;
for(int it=0; it<(int)ListSecondTracks.size(); it++){
std::vector<const Track*> tmpList(ListSecondTracks);
@@ -776,17 +604,19 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
if(Signif3Dproj>maxDst && maxDst<10. ){maxDst=Signif3Dproj; maxT=it; minChi2=Chi2;}
else if(Signif3Dproj>0. && maxDst>10. && Chi2<minChi2) {minChi2=Chi2; maxT=it;}
}
- if(maxT>=0){ Outlier=maxT; RemoveEntryInList(ListSecondTracks,cntComb,Outlier);
+ if(maxT>=0){ Outlier=maxT; RemoveEntryInList(ListSecondTracks,trkRank,Outlier);
m_fitSvc->setApproximateVertex(FitVertex.x(),FitVertex.y(),FitVertex.z());
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Remove negative outlier="<< maxT<<" from "
+ <<ListSecondTracks.size()+1<<" tracks"<<endmsg;
continue;}
}
////////////////////////////////////////
-// if( Momentum.m() <6000.) {
-// if( Chi2PerTrk[Outlier] < m_SecTrkChi2Cut && FitProb > 0.001) break; // Solution found
+// if( Momentum.m() < c_vrtBCMassLimit) {
+// if( Chi2PerTrk[Outlier] < m_secTrkChi2Cut && FitProb > 0.001) break; // Solution found
// }
if( FitProb > 0.001) {
- if( Momentum.M() <maxRecMASS) {
- if( Chi2PerTrk[Outlier] < m_SecTrkChi2Cut*m_chiScale[ListSecondTracks.size()<10?ListSecondTracks.size():10]) break; // Solution found
+ if( Momentum.M() <c_vrtBCMassLimit) {
+ if( Chi2PerTrk[Outlier] < m_secTrkChi2Cut*m_chiScale[ListSecondTracks.size()<10?ListSecondTracks.size():10]) break; // Solution found
} else {
double minM=1.e12; int minT=-1; double minChi2=1.e12;
for(int it=0; it<(int)ListSecondTracks.size(); it++){
@@ -794,28 +624,33 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
tmpList.erase(tmpList.begin()+it);
sc=VKalVrtFitBase(tmpList,tmpVertex,Momentum,Charge,ErrorMatrix,Chi2PerTrk,TrkAtVrt,Chi2);
if(sc.isFailure())continue;
- if(Momentum.M()<minM && minM>maxRecMASS){minM=Momentum.M(); minT=it; minChi2=Chi2;}
- else if(Momentum.M()<maxRecMASS && minM<maxRecMASS && Chi2<minChi2){minChi2=Chi2; minT=it;}
+ if(ProjSV_PV(tmpVertex,PrimVrt,JetDir)<0.)continue; // Drop negative direction
+ Chi2 += trkRank[it]; // Remove preferably non-HF-tracks
+ if(Momentum.M()<minM && minM>c_vrtBCMassLimit){minM=Momentum.M(); minT=it; minChi2=Chi2;}
+ else if(Momentum.M()<c_vrtBCMassLimit && minM<c_vrtBCMassLimit && Chi2<minChi2){minChi2=Chi2; minT=it;}
}
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" Big mass. Remove trk="<<minT<<" New mass="<<minM<<" New Chi2="<<minChi2<<endmsg;
if(minT>=0)Outlier=minT;
}
}
- RemoveEntryInList(ListSecondTracks,cntComb,Outlier);
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<"SecVrt remove trk="<<Outlier<<" from "<< ListSecondTracks.size()<<" tracks"<<endmsg;
+ RemoveEntryInList(ListSecondTracks,trkRank,Outlier);
m_fitSvc->setApproximateVertex(FitVertex.x(),FitVertex.y(),FitVertex.z()); /*Use as starting point*/
}
//--
- if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" SecVrt fit converged="<< ListSecondTracks.size()<<", "
- <<Chi2<<", "<<Chi2PerTrk[Outlier]<<" Mass="<<Momentum.M()<<endmsg;
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" SecVrt fit converged. Ntr="<< ListSecondTracks.size()<<" Chi2="<<Chi2
+ <<" Chi2_trk="<<Chi2PerTrk[Outlier]<<" Prob="<<FitProb<<" M="<<Momentum.M()<<" Dir="<<ProjSV_PV(FitVertex,PrimVrt,JetDir)<<endmsg;
//--
if( ListSecondTracks.size()==2 ){
- if( Momentum.M() > 6000. || FitProb < 0.001 || Chi2PerTrk[Outlier] > m_SecTrkChi2Cut) { return -10000.; }
+ if( Momentum.M() > c_vrtBCMassLimit || FitProb < 0.001 || Chi2PerTrk[Outlier] > m_secTrkChi2Cut) return -10000.;
+ if(std::max(trkRank[0],trkRank[1])<m_cutHFClass*2.) return -10000.;
}
//
//-- To kill remnants of conversion
double Dist2D=sqrt(FitVertex.x()*FitVertex.x()+FitVertex.y()*FitVertex.y());
if( ListSecondTracks.size()==2 && (Dist2D > 20.) && Charge==0 ) {
double mass_EE = massV0( TrkAtVrt,m_massE,m_massE);
- if(m_FillHist){m_hb_totmassEE->Fill( mass_EE, m_w_1); }
+ if(m_fillHist){m_hb_totmassEE->Fill( mass_EE, m_w_1); }
if( mass_EE < 40. ) return -40.;
}
//-- Test creation of Trk::Track
@@ -856,15 +691,17 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
std::vector<const Track*> & ListSecondTracks)
const
{
- Amg::Vector3D FitVertex, vDist;
- std::vector<double> ErrorMatrix,Chi2PerTrk,VKPerigee,CovPerigee,closeVrtSig(0);
- std::vector< std::vector<double> > TrkAtVrt;
- TLorentzVector Momentum;
+ std::vector<double> Chi2PerTrk,VKPerigee,CovPerigee,closeVrtSig(0),closeVrtCh2(0);
+ //TLorentzVector Momentum;
std::vector<double> Impact,ImpactError;
double ImpactSignif=0;
double Chi2, Signif3D, Dist2D, JetVrtDir;
long int Charge;
int i,j;
+ StatusCode sc; sc.setChecked();
+ Vrt2Tr tmpVrt;
+ std::vector<Vrt2Tr> all2TrVrt(0);
+ TLorentzVector TLV;
//
int NTracks = (int) (SelectedTracks.size());
@@ -873,14 +710,15 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
//
// Impact parameters with sign calculations
//
- double SignifR,SignifZ;
- std::vector<double> TrackSignif(NTracks),TrackPt(NTracks),TrackSignifBase(NTracks),adpt(NTracks),TrkSigZ(NTracks);
+ std::vector<float> covPV=PrimVrt.covariance();
+ double SignifR=0.,SignifZ=0.;
+ int nTrkHF=0;
+ std::vector<int> hitIBL(NTracks,0), hitBL(NTracks,0);
+ std::vector<double> TrackSignif(NTracks),TrkSig3D(NTracks);
+ std::vector< std::vector<float> > trkScore(NTracks);
AmgVector(5) tmpPerigee; tmpPerigee<<0.,0.,0.,0.,0.;
- int NPrimTrk=0, NSecTrk=0;
- m_NRefTrk=0;
for (i=0; i<NTracks; i++) {
- ImpactSignif = m_fitSvc->VKalGetImpact(SelectedTracks[i], PrimVrt.position(), 1, Impact, ImpactError);
- TrackSignifBase[i]=ImpactSignif;
+ TrkSig3D[i] = m_fitSvc->VKalGetImpact(SelectedTracks[i], PrimVrt.position(), 1, Impact, ImpactError);
tmpPerigee = GetPerigee(SelectedTracks[i])->parameters();
if( sin(tmpPerigee[2]-JetDir.Phi())*Impact[0] < 0 ){ Impact[0] = -fabs(Impact[0]);}
else{ Impact[0] = fabs(Impact[0]);}
@@ -888,48 +726,44 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
else{ Impact[1] = fabs(Impact[1]);}
SignifR = Impact[0]/ sqrt(ImpactError[0]);
SignifZ = Impact[1]/ sqrt(ImpactError[2]);
- if(m_FillHist){
+ TrackSignif[i] = sqrt( SignifR*SignifR + SignifZ*SignifZ);
+ int hL1=0, nLays=0; getPixelLayers(SelectedTracks[i] , hitIBL[i] , hitBL[i], hL1, nLays );
+ //----
+ trkScore[i]=m_trackClassificator->trkTypeWgts(SelectedTracks[i], PrimVrt, JetDir);
+ if( trkScore[i][0]>trkScore[i][1] && trkScore[i][0]>trkScore[i][2] ) nTrkHF++; //Good HF track
+ if(m_fillHist){
m_hb_impactR->Fill( SignifR, m_w_1);
m_hb_impactZ->Fill( SignifZ, m_w_1);
m_hb_impactRZ->Fill(SignifR, SignifZ, m_w_1);
- }
- TrackPt[i] = sin(tmpPerigee[3])/fabs(tmpPerigee[4]);
- if(ImpactSignif < 3.) { NPrimTrk += 1;}
- else{NSecTrk += 1;}
- if( SignifR<0 || SignifZ<0 ) m_NRefTrk++;
- if(m_getNegativeTail){
- ImpactSignif = sqrt( SignifR*SignifR + SignifZ*SignifZ);
- }else if(m_getNegativeTag){
- ImpactSignif = sqrt( (SignifR-0.6)*(SignifR-0.6)
- + (SignifZ-0.5)*(SignifZ-0.5) );
- }else{
- ImpactSignif = sqrt( (SignifR+0.6)*(SignifR+0.6)
- + (SignifZ+0.5)*(SignifZ+0.5) );
- }
- if(fabs(SignifR) < m_AntiPileupSigRCut) { // cut against tracks from pileup vertices
- if(SignifZ > 1.+m_AntiPileupSigZCut ) ImpactSignif=0.;
- if(SignifZ < 1.-m_AntiPileupSigZCut ) ImpactSignif=0.;
- }
- TrackSignif[i]=ImpactSignif;
- TrkSigZ[i]=SignifZ;
- if(m_FillHist){m_hb_impact->Fill( ImpactSignif, m_w_1);}
- adpt[i]=pow(TrackPt[i]/JetDir.Perp(),0.5);
- //adpt[i]=trkPtCorr(TrackPt[i]);
+ m_hb_impact->Fill( ImpactSignif, m_w_1);
+ if(i<DevTuple::maxNTrk && m_curTup){
+ m_curTup->p_prob[i]=RankBTrk(SelectedTracks[i]->pt(),JetDir.Pt(),0.);
+ m_curTup->s_prob[i]=RankBTrk(0.,0.,ImpactSignif);
+ m_curTup->SigR[i]=SignifR; m_curTup->SigZ[i]=SignifZ;
+ m_curTup->d0[i]=Impact[0]; m_curTup->Z0[i]=Impact[1];
+ m_curTup->idMC[i]=getG4Inter(SelectedTracks[i]);
+ if(isBTrk(SelectedTracks[i]))m_curTup->idMC[i]=2;
+ if(getMCPileup(SelectedTracks[i]))m_curTup->idMC[i]=3;
+ m_curTup->wgtB[i]=trkScore[i][0]; m_curTup->wgtL[i]=trkScore[i][1]; m_curTup->wgtG[i]=trkScore[i][2];
+ m_curTup->Sig3D[i]=TrkSig3D[i];
+ m_curTup->chg[i]=tmpPerigee[4]<0. ? 1: -1;
+ m_curTup->ibl[i]=hitIBL[i];
+ m_curTup->bl[i]=hitBL[i];
+ TLorentzVector TLV;
+ TLV.SetPtEtaPhiE(SelectedTracks[i]->pt(),SelectedTracks[i]->eta(),SelectedTracks[i]->phi(),SelectedTracks[i]->e());
+ m_curTup->pTvsJet[i]=TLV.Perp(JetDir.Vect());
+ TLorentzVector normJ; normJ.SetPtEtaPhiM(1.,JetDir.Eta(),JetDir.Phi(),0.);
+ m_curTup->prodTJ[i]=sqrt(TLV.Dot(normJ));
+ m_curTup->nVrtT[i]=0;
+ }
+ }
}
+ if(m_fillHist){ m_curTup->ptjet=JetDir.Perp(); m_curTup->etajet=fabs(JetDir.Eta()); m_curTup->phijet=JetDir.Phi();
+ m_curTup->nTrkInJet=std::min(NTracks,DevTuple::maxNTrk); };
+ if(nTrkHF==0) return; //====== No at all good HF tracks
- m_NRefTrk=TMath::Max(NPrimTrk,TMath::Max(2,(int)(0.3*NTracks)));
- //double addNTrkDep=TMath::Min((JetDir.Perp()/1.e6),m_TrkSigNTrkDep)*m_NRefTrk; // NTrk dependence
- double addNTrkDep=TMath::Min((JetDir.Perp()/150.e3),1.)*m_TrkSigNTrkDep*m_NRefTrk; // NTrk dependence
- double SelLim = m_TrkSigCut;
+ ListSecondTracks.reserve(2*NTracks); // Reserve memory for single vertex
- if(m_FillHist){ int nSelTPairs=0;
- for(i=0; i<NTracks; i++){ if(TrackSignif[i] < SelLim+m_TrkSigSumCut/2+adpt[i])continue; nSelTPairs++;}
- m_hb_nHImpTrkCnt->Fill((double)nSelTPairs,m_w_1);
- }
-
- StatusCode sc; if(sc.isSuccess())ImpactSignif=0; //Safety !
- //if(m_MultiVertex || m_MultiWithPrimary) m_WorkArray->m_Incomp.reserve(NTracks*(NTracks-1)); // Reserve memory for PGRAPH multivertex
- ListSecondTracks.reserve(2*NTracks); // Reserve memory for sigle vertex
Amg::Vector3D iniVrt(0.,0.,0.);
m_fitSvc->setDefault();
@@ -937,87 +771,76 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
m_fitSvc->setMomCovCalc(1); // Total momentum and its covariance matrix are calculated
for (i=0; i<NTracks-1; i++) {
for (j=i+1; j<NTracks; j++) {
- //--if(m_MultiVertex || m_MultiWithPrimary){m_WorkArray->m_Incomp.push_back(i);m_WorkArray->m_Incomp.push_back(j);} // For PGRAPH multivertex !!!ALL SELECTIONS MUST BE AFTER THIS LINE!!!
- if(TrackSignif[i]==0.)continue; // Pileup and other problems
- if(TrackSignif[j]==0.)continue; // Pileup and other problems
- double cutTrI= SelLim +adpt[i] +addNTrkDep/2.;
- double cutTrJ= SelLim +adpt[j] +addNTrkDep/2.;
- if(!m_MultiWithPrimary) { // Not used for multi-vertex with primary one search
- if(TrackSignif[i] < cutTrI) continue;
- if(TrackSignif[j] < cutTrJ) continue;
- if(TrackSignif[i]+TrackSignif[j] < cutTrI+cutTrJ +m_TrkSigSumCut +addNTrkDep) continue;
+ if(trkScore[i][0]==0.)continue;
+ if(trkScore[j][0]==0.)continue;
+ if(!m_multiWithPrimary) { // Not used for multi-vertex with primary one search
+ if( trkScore[i][0]<m_cutHFClass )continue; //Not classified HF tracks
+ if( trkScore[j][0]<m_cutHFClass )continue; //Not classified HF tracks
}
int BadTracks = 0; //Bad tracks identification
- TracksForFit.clear();
+ TracksForFit.resize(2);
m_fitSvc->setDefault(); //Reset VKalVrt settings
m_fitSvc->setMomCovCalc(1); // Total momentum and its covariance matrix are calculated
- TracksForFit.push_back( SelectedTracks[i] );
- TracksForFit.push_back( SelectedTracks[j] );
- sc=VKalVrtFitFastBase(TracksForFit,FitVertex); /* Fast crude estimation*/
- if( sc.isFailure() || FitVertex.perp() > m_Rlayer2*2. ) { /* No initial estimation */
+ TracksForFit[0]=SelectedTracks[i];
+ TracksForFit[1]=SelectedTracks[j];
+ sc=VKalVrtFitFastBase(TracksForFit,tmpVrt.FitVertex); /* Fast crude estimation*/
+ if( sc.isFailure() || tmpVrt.FitVertex.perp() > m_rLayer2*2. ) { /* No initial estimation */
iniVrt=PrimVrt.position();
- if( m_MultiWithPrimary ) iniVrt.setZero();
+ if( m_multiWithPrimary ) iniVrt.setZero();
} else {
- vDist=FitVertex-PrimVrt.position();
- JetVrtDir = JetDir.Px()*vDist.x() + JetDir.Py()*vDist.y() + JetDir.Pz()*vDist.z();
- if( m_MultiWithPrimary ) JetVrtDir=fabs(JetVrtDir); /* Always positive when primary vertex is seeked for*/
- if( JetVrtDir>0. ) iniVrt=FitVertex; /* Good initial estimation */
+ JetVrtDir = ProjSV_PV(tmpVrt.FitVertex,PrimVrt,JetDir);
+ if( m_multiWithPrimary ) JetVrtDir=fabs(JetVrtDir); /* Always positive when primary vertex is seeked for*/
+ if( JetVrtDir>0. ) iniVrt=tmpVrt.FitVertex; /* Good initial estimation */
else iniVrt=PrimVrt.position();
}
m_fitSvc->setApproximateVertex(iniVrt.x(), iniVrt.y(), iniVrt.z());
- sc=VKalVrtFitBase(TracksForFit,FitVertex, Momentum,Charge,
- ErrorMatrix,Chi2PerTrk,TrkAtVrt,Chi2);
- if(sc.isFailure()) continue; /* No fit */
- if(Chi2 > m_Sel2VrtChi2Cut+4.) continue; /* Bad Chi2 */
- if(fabs(FitVertex.z())> 650.) continue; // definitely outside of Pixel detector
- Dist2D=FitVertex.perp();
+ tmpVrt.i=i; tmpVrt.j=j;
+ sc=VKalVrtFitBase(TracksForFit,tmpVrt.FitVertex, tmpVrt.Momentum, Charge,
+ tmpVrt.ErrorMatrix, tmpVrt.Chi2PerTrk, tmpVrt.TrkAtVrt, tmpVrt.Chi2);
+ if(sc.isFailure()) continue; /* No fit */
+ if(tmpVrt.Chi2 > m_sel2VrtChi2Cut) continue; /* Bad Chi2 */
+ if(fabs(tmpVrt.FitVertex.z())> 650.) continue; // definitely outside of Pixel detector
+ Dist2D=tmpVrt.FitVertex.perp();
if(Dist2D > 180. ) continue; // can't be from B decay
- if(m_useMaterialRejection && Dist2D>m_Rbeampipe-2.)
- { if( TrkSigZ[i]>25. || TrkSigZ[j]>25. || TrkSigZ[i]<-10. || TrkSigZ[j]<-10.) continue; }
- VrtVrtDist(PrimVrt, FitVertex, ErrorMatrix, Signif3D);
+ double mass_PiPi = tmpVrt.Momentum.M();
+ if(mass_PiPi > c_vrtBCMassLimit) continue; // can't be from B decay
+ VrtVrtDist(PrimVrt, tmpVrt.FitVertex, tmpVrt.ErrorMatrix, Signif3D);
+ tmpVrt.Signif3D=Signif3D;
+ VrtVrtDist2D(PrimVrt, tmpVrt.FitVertex, tmpVrt.ErrorMatrix, tmpVrt.Signif2D);
//---
- vDist=FitVertex-PrimVrt.position();
- JetVrtDir = JetDir.Px()*vDist.x() + JetDir.Py()*vDist.y() + JetDir.Pz()*vDist.z();
- double vPos=(vDist.x()*Momentum.Px()+vDist.y()*Momentum.Py()+vDist.z()*Momentum.Pz())/Momentum.Rho();
- if((!m_MultiWithPrimary) &&(!m_getNegativeTail) && (!m_getNegativeTag) && JetVrtDir<0. ) continue; /* secondary vertex behind primary*/
- if(vPos<-150.) continue; /* Secondary vertex is too far behind primary*/
+ TVector3 SVmPV(tmpVrt.FitVertex.x()-PrimVrt.x(),tmpVrt.FitVertex.y()-PrimVrt.y(),tmpVrt.FitVertex.z()-PrimVrt.z());
+ JetVrtDir = SVmPV.Dot(JetDir.Vect());
+ double vPos=SVmPV.Dot(tmpVrt.Momentum.Vect())/tmpVrt.Momentum.Rho();
+ if((!m_multiWithPrimary) &&(!m_getNegativeTail) && (!m_getNegativeTag) && JetVrtDir<0. ) continue; /* secondary vertex behind primary*/
+ if(vPos<-100.) continue; /* Secondary vertex is too far behind primary*/
//
// Check pixel hits vs vertex positions.
- if(m_useVertexCleaning){ if(!Check2TrVertexInPixel(SelectedTracks[i],SelectedTracks[j],FitVertex,Signif3D)) continue; }
+ if(m_useVertexCleaning){ if(!Check2TrVertexInPixel(SelectedTracks[i],SelectedTracks[j],tmpVrt.FitVertex,tmpVrt.ErrorMatrix)) continue; }
//--------
//
- double Signif3Dproj=VrtVrtDist( PrimVrt, FitVertex, ErrorMatrix, JetDir);
+ double Signif3Dproj=VrtVrtDist( PrimVrt, tmpVrt.FitVertex, tmpVrt.ErrorMatrix, JetDir);
+ tmpVrt.Signif3DProj=Signif3Dproj;
double Signif3DSign=Signif3D; if(JetVrtDir<0) Signif3DSign=-Signif3D;
- if(m_FillHist)m_hb_signif3D->Fill( Signif3DSign, m_w_1);
- if(Signif3DSign<12. && Chi2>m_Sel2VrtChi2Cut) continue; /* Bad Chi2 */
-
- if( m_MultiWithPrimary || m_MultiVertex) { // For multivertex
- double limPtTr=m_JetPtFractionCut*JetDir.Perp();
- if(TrackPt[i]<limPtTr && Signif3D<9. && m_MultiWithPrimary) continue;
- if(TrackPt[j]<limPtTr && Signif3D<9. && m_MultiWithPrimary) continue;
- //m_WorkArray->m_Incomp.pop_back();m_WorkArray->m_Incomp.pop_back(); //For PGRAPH
- add_edge(i,j,*m_compatibilityGraph);
- }
-
- if( m_MultiWithPrimary ) continue; /* Multivertex with primary one. All below is not needed */
- double mass_PiPi = Momentum.M();
- if(mass_PiPi > 6000.) continue; // can't be from B decay
+ if(m_fillHist){ m_hb_signif3D->Fill( Signif3DSign, m_w_1); m_hb_sig3DNtr->Fill(Signif3Dproj, 1.);}
+
+ //if( m_multiWithPrimary || m_multiVertex) { // For multivertex
+ // add_edge(i,j,*m_compatibilityGraph);
+ //}
+ if( m_multiWithPrimary ) continue; /* Multivertex with primary one. All below is not needed */
//
// Check if V0 or material interaction on Pixel layer is present
//
if( Charge == 0 && Signif3D>8. && mass_PiPi<900.) {
- double mass_PPi = massV0( TrkAtVrt,m_massP,m_massPi);
- double mass_EE = massV0( TrkAtVrt,m_massE,m_massE);
- if(m_FillHist && !m_MultiVertex){m_hb_massEE->Fill( mass_EE, m_w_1);}
+ double mass_PPi = massV0( tmpVrt.TrkAtVrt,m_massP,m_massPi);
+ double mass_EE = massV0( tmpVrt.TrkAtVrt,m_massE,m_massE);
+ if(m_fillHist && !m_multiVertex){m_hb_massEE->Fill( mass_EE, m_w_1);}
if( mass_EE < 40.) {
BadTracks = 3;
}else{
- if(m_FillHist && !m_MultiVertex){m_hb_massPiPi->Fill( mass_PiPi, m_w_1);} /* Total mass with input particles masses*/
- if(m_FillHist && !m_MultiVertex){m_hb_massPPi->Fill( mass_PPi, m_w_1);}
+ if(m_fillHist && !m_multiVertex){m_hb_massPiPi->Fill( mass_PiPi, m_w_1);} /* Total mass with input particles masses*/
+ if(m_fillHist && !m_multiVertex){m_hb_massPPi->Fill( mass_PPi, m_w_1);}
if( fabs(mass_PiPi-m_massK0) < 22. ) BadTracks = 1;
if( fabs(mass_PPi-m_massLam) < 8. ) BadTracks = 2;
- //double TransMass=TotalTMom(GetPerigeeVector(TracksForFit));
- //if( TransMass<400. && m_FillHist)m_hb_massPiPi1->Fill( mass_PiPi, m_w_1);
}
//
// Creation of V0 tracks
@@ -1034,7 +857,7 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
m_fitSvc->setCnstType(1); // Set mass constraint
}
if( BadTracks == 2 ) { // Lambda case
- if( fabs(1./TrkAtVrt[0][2]) > fabs(1./TrkAtVrt[1][2]) ) {
+ if( fabs(1./tmpVrt.TrkAtVrt[0][2]) > fabs(1./tmpVrt.TrkAtVrt[1][2]) ) {
InpMassV0.push_back(m_massP);InpMassV0.push_back(m_massPi);
}else{ InpMassV0.push_back(m_massPi);InpMassV0.push_back(m_massP); }
m_fitSvc->setMassInputParticles( InpMassV0 );
@@ -1045,41 +868,45 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
InpMassV0.push_back(m_massE);InpMassV0.push_back(m_massE);
m_fitSvc->setMassInputParticles( InpMassV0 );
m_fitSvc->setCnstType(12); // Set 3d angular constraint
-// m_fitSvc->setVertexForConstraint(PrimVrt.x(),PrimVrt.y(),PrimVrt.z());
-// m_fitSvc->setCovVrtForConstraint(0.015*0.015,0.,0.015*0.015,0.,0.,56.*56);
-// m_fitSvc->setCnstType(7); // Set pointing to primary vertex constraint
}
- m_fitSvc->setApproximateVertex(FitVertex.x(),FitVertex.y(),FitVertex.z());
+ m_fitSvc->setApproximateVertex(tmpVrt.FitVertex.x(),tmpVrt.FitVertex.y(),tmpVrt.FitVertex.z());
TLorentzVector MomentumV0;
Amg::Vector3D FitVertexV0;
- std::vector< std::vector<double> > TrkAtVrt0;
- double Chi2_0;
+ std::vector< std::vector<double> > TrkAtVrtV0;
+ std::vector<double> ErrorMatrixV0;
+ double Chi2V0;
sc=VKalVrtFitBase(TracksForFit, FitVertexV0, MomentumV0, Charge,
- ErrorMatrix,Chi2PerTrk,TrkAtVrt0,Chi2_0);
+ ErrorMatrixV0,Chi2PerTrk,TrkAtVrtV0,Chi2V0);
if(sc.isSuccess()) {
-
- sc=m_fitSvc->VKalVrtCvtTool(FitVertexV0,MomentumV0,ErrorMatrix,0,VKPerigee,CovPerigee);
+ sc=m_fitSvc->VKalVrtCvtTool(FitVertexV0,MomentumV0,ErrorMatrixV0,0,VKPerigee,CovPerigee);
if(sc.isSuccess()) {
const Trk::Track* TT = m_fitSvc->CreateTrkTrack(VKPerigee,CovPerigee);
- ImpactSignif=m_fitSvc->VKalGetImpact(TT, PrimVrt.position(), 0, Impact, ImpactError);
- if(m_FillHist){m_hb_impV0->Fill( ImpactSignif, m_w_1); }
- if(ImpactSignif>3.5) BadTracks=0;
+ double ImpactSignifV0=m_fitSvc->VKalGetImpact(TT, PrimVrt.position(), 0, Impact, ImpactError);
+ if(m_fillHist){m_hb_impV0->Fill( ImpactSignifV0, m_w_1);}
+ if(ImpactSignifV0>3.0 ) BadTracks=0;
delete TT;
} else { BadTracks=0;}
} // else { BadTracks=0;}
}
}
//
-// Check interactions on pixel layers
+// Check interactions on material layers
//
- if(m_FillHist){ m_hb_r2d->Fill( FitVertex.perp(), m_w_1); }
- if(m_useMaterialRejection && Dist2D>m_Rbeampipe-2.){
- float ptLim=TMath::Max(m_hadronIntPtCut,m_JetPtFractionCut*JetDir.Perp());
-
- if( TMath::Min(TrackPt[i],TrackPt[j])<ptLim ){
- if( insideMatLayer(FitVertex.x(), FitVertex.y()) ) BadTracks = 4;
- }
- }
+ float minWgtI = std::min(trkScore[i][2],trkScore[j][2]);
+ if( minWgtI >0.50 && Dist2D > m_beampipeR-VrtRadiusError(tmpVrt.FitVertex, tmpVrt.ErrorMatrix) ) BadTracks = 4;
+ //if( (trkScore[i][2]>0.4 || trkScore[j][2]>0.4)
+ // && insideMatLayer(tmpVrt.FitVertex.x(),tmpVrt.FitVertex.y()) ) BadTracks=5;
+//
+//-----------------------------------------------
+ tmpVrt.badVrt=BadTracks; //
+ all2TrVrt.push_back(tmpVrt); //
+//-----------------------------------------------
+ if(m_fillHist){ m_hb_r2d->Fill( tmpVrt.FitVertex.perp(), m_w_1); }
+ //if(m_useMaterialRejection && Dist2D>m_beampipeR-2.){
+ //if(m_materialMap){
+ // if(m_materialMap->inMaterial(tmpVrt.FitVertex)) BadTracks=4;
+ // if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" MaterialMap test="<< BadTracks<<endreq;
+ //}
//
// Creation of tracks from V0 list
//
@@ -1087,78 +914,86 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
TrkFromV0.push_back(SelectedTracks[i]);
TrkFromV0.push_back(SelectedTracks[j]);
}else{
- double limPtTr=m_JetPtFractionCut*JetDir.Perp();
- if(TrackPt[i]<limPtTr && Signif3D<15.)continue;
- if(TrackPt[j]<limPtTr && Signif3D<15.)continue;
//
double JetVrtDir =
- JetDir.Px()*(FitVertex.x()-PrimVrt.x())
- + JetDir.Py()*(FitVertex.y()-PrimVrt.y())
- + JetDir.Pz()*(FitVertex.z()-PrimVrt.z());
+ JetDir.Px()*(tmpVrt.FitVertex.x()-PrimVrt.x())
+ + JetDir.Py()*(tmpVrt.FitVertex.y()-PrimVrt.y())
+ + JetDir.Pz()*(tmpVrt.FitVertex.z()-PrimVrt.z());
if(m_getNegativeTail) JetVrtDir=fabs(JetVrtDir); // For negative TAIL
// accepts also negative
// track pairs
if(m_getNegativeTag) JetVrtDir=-JetVrtDir; // For negative TAG
// accepts only negative
// track pairs
- if( (Signif3D>m_Sel2VrtSigCut) && (JetVrtDir>0) ) {
- if(fabs(Signif3Dproj)<m_Sel2VrtSigCut)continue;
+ if( (Signif3D>m_sel2VrtSigCut) ) {
+ if(fabs(Signif3Dproj)<m_sel2VrtSigCut)continue;
ListSecondTracks.push_back(SelectedTracks[i]);
ListSecondTracks.push_back(SelectedTracks[j]);
closeVrtSig.push_back(Signif3D);
+ closeVrtCh2.push_back(Chi2);
}
}
}
}
+// if(m_fillHist) fillVrtNTup(all2TrVrt);
+
//
-//-- Post-selection cleaning
+//-- Start vertex analysis
+/*
+ std::vector<int> trkTypeSV(NTracks,-1); // Track identification: -1 -unknown, 2-fragmentation, 1-Interaction, 0-HF
+ std::vector<int> nFrVrtT(NTracks,0); // Fragmentation vertex per track counter
+ int foundHF=0;
+ for( auto vv : all2TrVrt){
+ if( 1.-1./11.*vv.Signif2D > std::min(trkScore[vv.i][0],trkScore[vv.j][0]) ){
+ if( std::min(trkScore[vv.i][1],trkScore[vv.j][1]) > 0.5 ){ nFrVrtT[vv.i]++; nFrVrtT[vv.j]++; }
+ continue;
+ }
+ if( trkScore[vv.i][0]+trkScore[vv.j][0] < trkScore[vv.i][1]+trkScore[vv.j][1]) continue;
+ if( trkScore[vv.i][0]+trkScore[vv.j][0] < trkScore[vv.i][2]+trkScore[vv.j][2]) continue;
+ trkTypeSV[vv.i]=0; trkTypeSV[vv.j]=0; foundHF++;
+ //if( trkScore[vv.i][0]>trkScore[vv.i][1] && trkScore[vv.j][0]>trkScore[vv.j][1]
+ // && trkScore[vv.i][0]>trkScore[vv.i][2] && trkScore[vv.j][0]>trkScore[vv.j][2] ){ trkTypeSV[vv.i]=0; trkTypeSV[vv.j]=0; foundHF++; }
+ }
+ for( auto vv : all2TrVrt){ //Now interactions+V0s
+ if( foundHF==1 && (trkTypeSV[vv.i]==0 || trkTypeSV[vv.j]==0 )) continue; //preserve single HF-vertex
+ if( vv.badVrt ){ trkTypeSV[vv.i]=1; trkTypeSV[vv.j]=1;}
+ }
+ for( auto vv : all2TrVrt){ //Now fragmentation
+ if( trkTypeSV[vv.i]>=0 || trkTypeSV[vv.j]>=0 ) continue; //Skip identified tracks
+ if( trkScore[vv.i][1]>trkScore[vv.i][0] && trkScore[vv.j][1]>trkScore[vv.j][0]
+ && trkScore[vv.i][1]>trkScore[vv.i][2] && trkScore[vv.j][1]>trkScore[vv.j][2] ){ trkTypeSV[vv.i]=2; trkTypeSV[vv.j]=2;}
+ }
+ for (i=0; i<NTracks; i++) if( trkTypeSV[i]==0 && nFrVrtT[i]>0 && trkScore[i][1]>0.5 ) trkTypeSV[i]=2;
+//
+//-- Remove FF, IF(some) and II vertices
+// iv=0; while ( iv < (int)all2TrVrt.size() )
+// { if( trkTypeSV[all2TrVrt[iv].i]>0 && trkTypeSV[all2TrVrt[iv].j]>0) all2TrVrt.erase(all2TrVrt.begin()+iv); else iv++; }
//
- if(ListSecondTracks.size()==1 && closeVrtSig[0]<9.){ //Turned off for the moment (ListSecondTracks.size() can't be 1!)
- auto it0=std::find(SelectedTracks.begin(),SelectedTracks.end(),ListSecondTracks[0]);
- auto it1=std::find(SelectedTracks.begin(),SelectedTracks.end(),ListSecondTracks[1]);
- int eGood=1; double cutNTrkDep=cutNDepNorm(m_NRefTrk,0.05); // NTrk dependence
- if(it0!=SelectedTracks.end() && pow(TrackSignifBase[it0-SelectedTracks.begin()],2.) < cutNTrkDep ) eGood=0;
- if(it1!=SelectedTracks.end() && pow(TrackSignifBase[it1-SelectedTracks.begin()],2.) < cutNTrkDep ) eGood=0;
- if(!eGood)ListSecondTracks.clear();
+*/
+//-- Save results
+ ListSecondTracks.clear();
+ std::map<float,int> trkHF;
+ for( auto vv : all2TrVrt){
+ if( m_multiWithPrimary || m_multiVertex) add_edge(vv.i,vv.j,*m_compatibilityGraph);
+ trkHF[trkScore[vv.i][0]]=vv.i; trkHF[trkScore[vv.j][0]]=vv.j;
+ }
+ for( auto it : trkHF) { ListSecondTracks.push_back(SelectedTracks[it.second]); }
+//-Debug
+ if( m_fillHist && m_curTup ){
+ for( auto it : trkHF) { m_curTup->itHF[m_curTup->NTHF++]=it.second; }
+ for( auto vv : all2TrVrt){ m_curTup->nVrtT[vv.i]++; m_curTup->nVrtT[vv.j]++; }
+ fillVrtNTup(all2TrVrt);
}
- //--------------------------------------------------------------------
+//
+//--------------------------------------------------------------------
+//-- Post-selection checks
+//--------------------------------------------------------------------
if(ListSecondTracks.size()>0 ){
- if(m_FillHist){ m_pr_effVrt2tr->Fill((float)m_NRefTrk,1.);
+ if(m_fillHist){ m_pr_effVrt2tr->Fill((float)m_NRefPVTrk,1.);
//m_pr_effVrt2trEta->Fill( JetDir.Eta(),1.);}
m_pr_effVrt2trEta->Fill( JetDir.Eta(),(double)ListSecondTracks.size()/2.);}
- //----- Only vertices with unique tracks are considered as bad
- if(TrkFromV0.size()){
- std::vector<const Track*> tmpVec(0);
- for(int tk=0;tk<(int)TrkFromV0.size(); tk+=2){
- int nCheck1=std::count(TrkFromV0.begin(),TrkFromV0.end(),TrkFromV0[tk]);
- int nCheck2=std::count(TrkFromV0.begin(),TrkFromV0.end(),TrkFromV0[tk+1]);
- int nFound1=std::count(ListSecondTracks.begin(),ListSecondTracks.end(),TrkFromV0[tk]);
- int nFound2=std::count(ListSecondTracks.begin(),ListSecondTracks.end(),TrkFromV0[tk+1]);
- if(nFound1+nFound2 && nCheck1==1 && nCheck2==1 ){ //Unique in TrkFromV0 but duplicated in ListSecondTracks
- ListSecondTracks.push_back(TrkFromV0[tk]);
- ListSecondTracks.push_back(TrkFromV0[tk+1]);
- }else{
- tmpVec.push_back(TrkFromV0[tk]);
- tmpVec.push_back(TrkFromV0[tk+1]);
- }
- }
- TrkFromV0.swap(tmpVec);
- }
- } else if(ListSecondTracks.size()==0) { if(m_FillHist){m_pr_effVrt2tr->Fill((float)m_NRefTrk,0.);
+ } else if(ListSecondTracks.size()==0) { if(m_fillHist){m_pr_effVrt2tr->Fill((float)m_NRefPVTrk,0.);
m_pr_effVrt2trEta->Fill(JetDir.Eta(),0.); }}
-/////////// Attempt to find iasolated tracks with high impact parameter. They are RARE!!! No worth to use them!
-// TLorentzVector psum,tmp; int nprt=0;
-// for (i=0; i<NTracks; i++) {
-// if( TrackSignif[i]<2.*SelLim || signTrk[i]<0) continue;
-// auto it0=std::find(ListSecondTracks.begin(),ListSecondTracks.end(),SelectedTracks[i]); if( it0!=ListSecondTracks.end() )continue;
-// it0=std::find(TrkFromV0.begin(),TrkFromV0.end(),SelectedTracks[i]); if( it0!=TrkFromV0.end() )continue;
-// it0=std::find(ListCloseTracks.begin(),ListCloseTracks.end(),SelectedTracks[i]); if( it0!=ListCloseTracks.end() )continue;
-// int ibl,bl,l1,nlay; getPixelLayers(SelectedTracks[i], ibl, bl, l1, nlay); if(ibl+bl<2)continue;
-// nprt++; tmp.SetPtEtaPhiM(SelectedTracks[i]->pt(),SelectedTracks[i]->eta(),SelectedTracks[i]->phi(),m_massPi); psum += tmp;
-// } if(nprt<1)return;
-// if(nprt==1) { tmp.SetPtEtaPhiM(psum.Pt(),JetDir.Eta(),JetDir.Phi(),m_massPi); psum += tmp; }
-// if(m_FillHist)m_hb_massPiPi1->Fill( psum.M(), m_w_1);
-/////////////////////////////////////
return;
}
@@ -1167,8 +1002,7 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
template <class Track>
bool InDetVKalVxInJetTool::Check2TrVertexInPixel( const Track* p1, const Track* p2,
- //Amg::Vector3D &FitVertex, double Signif3D)
- Amg::Vector3D &FitVertex, double)
+ Amg::Vector3D &FitVertex, std::vector<double> & VrtErr)
const
{
int blTrk[2]={0,0};
@@ -1181,46 +1015,40 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
getPixelLayers( p2, blTrk[1] , l1Trk[1], l2Trk[1], nLays[1] ); // Very close to PV. Both b-layer hits are mandatory.
getPixelProblems(p1, blP[0], l1P[0] );
getPixelProblems(p2, blP[1], l1P[1] );
- //if( Signif3D<15. && FitVertex.perp()<15. ){
- // if( blTrk[0]<1 && l1Trk[0]<1 ) return false;
- // if( blTrk[1]<1 && l1Trk[1]<1 ) return false;
- //}
- double xDif=FitVertex.x()-m_XlayerB, yDif=FitVertex.y()-m_YlayerB ;
+ double xDif=FitVertex.x()-m_xLayerB, yDif=FitVertex.y()-m_yLayerB ;
double Dist2DBL=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2DBL < m_RlayerB-m_SVResolutionR){ //----------------------------------------- Inside B-layer
+ if (Dist2DBL < m_rLayerB-VrtRadiusError(FitVertex, VrtErr)) { //----------------------------------------- Inside B-layer
if( blTrk[0]==0 && blTrk[1]==0) return false; // No b-layer hits at all, but all expected
if( blTrk[0]<1 && l1Trk[0]<1 ) return false;
if( blTrk[1]<1 && l1Trk[1]<1 ) return false;
if( nLays[0] <2 ) return false; // Less than 2 layers on track 0
if( nLays[1] <2 ) return false; // Less than 2 layers on track 1
return true;
- }else if(Dist2DBL > m_RlayerB+m_SVResolutionR){ //----------------------------------------- Outside b-layer
- if( blTrk[0]>0 && blP[0]==0 ) return false; // Good hit in b-layer is present
- if( blTrk[1]>0 && blP[1]==0 ) return false; // Good hit in b-layer is present
- }
+ }else if(Dist2DBL > m_rLayerB+VrtRadiusError(FitVertex, VrtErr)){ //----------------------------------------- Outside b-layer
+ if( blTrk[0]>0 && blP[0]==0 && blTrk[1]>0 && blP[1]==0 ) return false; // Good hit in b-layer is present
+ }
//
// L1 and L2 are considered only if vertex is in acceptance
//
if(fabs(FitVertex.z())<400.){
- xDif=FitVertex.x()-m_Xlayer1, yDif=FitVertex.y()-m_Ylayer1 ;
+ xDif=FitVertex.x()-m_xLayer1, yDif=FitVertex.y()-m_yLayer1 ;
double Dist2DL1=sqrt(xDif*xDif+yDif*yDif);
- xDif=FitVertex.x()-m_Xlayer2, yDif=FitVertex.y()-m_Ylayer2 ;
+ xDif=FitVertex.x()-m_xLayer2, yDif=FitVertex.y()-m_yLayer2 ;
double Dist2DL2=sqrt(xDif*xDif+yDif*yDif);
- if (Dist2DL1 < m_Rlayer1-m_SVResolutionR) { //------------------------------------------ Inside 1st-layer
+ if (Dist2DL1 < m_rLayer1-VrtRadiusError(FitVertex, VrtErr)) { //------------------------------------------ Inside 1st-layer
if( l1Trk[0]==0 && l1Trk[1]==0 ) return false; // No L1 hits at all
if( l1Trk[0]<1 && l2Trk[0]<1 ) return false; // Less than 1 hits on track 0
if( l1Trk[1]<1 && l2Trk[1]<1 ) return false; // Less than 1 hits on track 1
return true;
- }else if(Dist2DL1 > m_Rlayer1+m_SVResolutionR) { //------------------------------------------- Outside 1st-layer
- if( l1Trk[0]>0 && l1P[0]==0 ) return false; // Good L1 hit is present
- if( l1Trk[1]>0 && l1P[1]==0 ) return false; // Good L1 hit is present
+ }else if(Dist2DL1 > m_rLayer1+VrtRadiusError(FitVertex, VrtErr)) { //------------------------------------------- Outside 1st-layer
+ if( l1Trk[0]>0 && l1P[0]==0 && l1Trk[1]>0 && l1P[1]==0 ) return false; // Good L1 hit is present
}
- if (Dist2DL2 < m_Rlayer2-m_SVResolutionR) { //------------------------------------------- Inside 2nd-layer
+ if (Dist2DL2 < m_rLayer2-VrtRadiusError(FitVertex, VrtErr)) { //------------------------------------------- Inside 2nd-layer
if( (l2Trk[0]+l2Trk[1])==0 ) return false; // At least one L2 hit must be present
- }else if(Dist2DL2 > m_Rlayer2+m_SVResolutionR) {
+ }else if(Dist2DL2 > m_rLayer2+VrtRadiusError(FitVertex, VrtErr)) {
// if( (l2Trk[0]+l2Trk[1])>0 ) return false; // L2 hits are present
- }
+ }
} else {
int d0Trk[2]={0,0};
int d1Trk[2]={0,0};
@@ -1233,60 +1061,4 @@ for (auto atrk : AdditionalTracks)ListSecondTracks.push_back(atrk.second);
return true;
}
-
-
-
- template <class Track>
- double InDetVKalVxInJetTool::RemoveNegImpact(std::vector<const Track*> & inTracks,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- double Limit)
- const
- {
- int blTrk=0, l1Trk=0, l2Trk=0, nLays=0;
- int i=0;
- while( i<(int)inTracks.size()){
- getPixelLayers( inTracks[i], blTrk , l1Trk, l2Trk, nLays );
- //if(nLays<1 || inTracks[i]->pt()>JetDir.Pt()) inTracks.erase(inTracks.begin()+i); //bad track: no pixel hit OR trk_pt>jet_pt
- if(nLays<1) inTracks.erase(inTracks.begin()+i); //bad track: no pixel hit
- else i++;
- }
-//----
- int NTracks = inTracks.size();
- std::vector<double> ImpR(NTracks);
- std::vector<double> Impact,ImpactError;
- AmgVector(5) tmpPerigee; tmpPerigee<<0.,0.,0.,0.,0.;
- double maxImp=-1.e10, zImp=0.;
- for (i=0; i<NTracks; i++) {
- m_fitSvc->VKalGetImpact(inTracks[i], PrimVrt.position(), 1, Impact, ImpactError);
- tmpPerigee = GetPerigee(inTracks[i])->parameters();
- if( sin(tmpPerigee[2]-JetDir.Phi())*Impact[0] < 0 ){ Impact[0] = -fabs(Impact[0]);}
- else{ Impact[0] = fabs(Impact[0]);}
- if( (tmpPerigee[3]-JetDir.Theta())*Impact[1] < 0 ){ Impact[1] = -fabs(Impact[1]);}
- else{ Impact[1] = fabs(Impact[1]);}
- double SignifR = Impact[0]/ sqrt(ImpactError[0]); ImpR[i]=SignifR;
- if(fabs(SignifR) < m_AntiPileupSigRCut) { // cut against tracks from pileup vertices
- if( fabs(Impact[1])/sqrt(ImpactError[2]) > m_AntiPileupSigZCut ) ImpR[i]=-9999.;
- }
- if(fabs(Impact[1])/sqrt(ImpactError[2])>fabs(ImpR[i]) && Impact[1]<0 && ImpR[i]>0) ImpR[i]*=-1.;
- if(ImpR[i]>maxImp){maxImp=ImpR[i]; zImp=Impact[1]/sqrt(ImpactError[2]);}
- }
- if(maxImp<Limit){ inTracks.clear(); return maxImp;}
- double rmin=1.e6;
- do{ rmin=1.e6; int jpm=0;
- for(i=0; i<(int)ImpR.size(); i++){ if(rmin>ImpR[i]){rmin=ImpR[i]; jpm=i;}}; if(rmin>Limit)continue;
- ImpR.erase(ImpR.begin()+jpm); inTracks.erase(inTracks.begin()+jpm);
- }while(rmin<=Limit);
- if(inTracks.size()==1 && zImp<1.){ inTracks.clear(); maxImp=0.;}
- return maxImp;
- }
-
-
-template
-bool InDetVKalVxInJetTool::Check2TrVertexInPixel( const xAOD::TrackParticle* p1,
- const xAOD::TrackParticle* p2,
- Amg::Vector3D &FitVertex, double)
- const;
-
-
} //end of namespace
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
index 460fc48ec0f01f9fd075384591f6e9266fe8d81a..0a73eb71e87e7eb1a6e7b96e77c5e3ea70c2ea91 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// Header include
@@ -30,19 +30,15 @@
// 7) Jet energy used in (2) calculation
//----------------------------------------------------------------------------------------
-namespace Trk {
-extern int pgraphm_(
- long int *weit, long int *edges, long int *nodes,
- long int *set, long int *nptr, long int *nth);
-}
-//extern "C" {
-// float prob_(const float & Chi2, const long int& ndf);
-//}
+//namespace Trk {
+// extern int pgraphm_(
+// long int *weit, long int *edges, long int *nodes,
+// long int *set, long int *nptr, long int *nth); }
namespace InDet{
-const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comoming from B,C-decays
+const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not coming from B,C-decays
// std::vector<xAOD::Vertex*> InDetVKalVxInJetTool::GetVrtSecMulti(
@@ -61,8 +57,9 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
{
const double probVrtMergeLimit=0.04;
+ const int useMaterialRejection=1;
- m_NRefTrk=0;
+ m_NRefPVTrk=0;
int inpNPart=0;
int i,j;
if(xAODwrk){
@@ -95,21 +92,11 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
NTracks = RECwork->listJetTracks.size();
MomentumJet = TotalMom(GetPerigeeVector(RECwork->listJetTracks));}
- if(NTracks>m_TrackInJetNumberLimit){
- if (xAODwrk ) xAODwrk->listJetTracks.resize(m_TrackInJetNumberLimit);
- else if(RECwork ) RECwork->listJetTracks.resize(m_TrackInJetNumberLimit);
- NTracks=m_TrackInJetNumberLimit;
- }
if( NTracks < 2 ) { return finalVertices;} // 0,1 selected track => nothing to do!
- if(xAODwrk){
- while( xAODwrk->listJetTracks.size()>4 && medianPtF(xAODwrk->listJetTracks)/JetDir.Pt()<0.01)
- xAODwrk->listJetTracks.pop_back();
- NTracks = xAODwrk->listJetTracks.size();
- }
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "Number of selected tracks inside jet= " <<NTracks << endmsg;
- if(m_FillHist){m_hb_jmom->Fill( MomentumJet.Perp(), m_w_1);
+ if(m_fillHist){m_hb_jmom->Fill( MomentumJet.Perp(), m_w_1);
m_hb_ntrkjet->Fill( (double) NTracks, m_w_1); }
//
@@ -162,8 +149,10 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
//================================================ PGRAPH version
// int iRet=0;
-// long int* weit = new long int[m_WorkArray->m_Incomp.size()];
-// long int* Solution = new long int[NTracks];
+// //long int* weit = new long int[m_WorkArray->m_Incomp.size()];
+// //long int* Solution = new long int[NTracks];
+// std::unique_ptr<long int[]> weit = std::unique_ptr<long int[]>(new long int[m_WorkArray->m_Incomp.size()])
+// std::unique_ptr<long int[]> Solution = std::unique_ptr<long int[]>(new long int[NTracks])
// for(i=0; i<(int)m_WorkArray->m_Incomp.size(); i++) *(weit+i)=(long int) (m_WorkArray->m_Incomp[i]+1); /* +1 is needed for PGRAPH*/
// long int edges = m_WorkArray->m_Incomp.size()/2;
// while(true) {
@@ -174,8 +163,6 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
// newvrt.SelTrk.clear();
// for(i=0;i<NPTR;i++) { newvrt.SelTrk.push_back(Solution[i]-1);}//std::cout<<"Solution="<<Solution[i]<<'\n';
//================================================== Boost version (don't forget to uncomment addEdge in Select2TrVrt()
- const long int* const weit=0; //coverity 105578
- const long int* const Solution=0; //coverity 105577
std::vector< std::vector<int> > allCliques;
bron_kerbosch_all_cliques(*m_compatibilityGraph, clique_visitor(allCliques));
for(int cq=0; cq<(int)allCliques.size();cq++){
@@ -190,13 +177,13 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
}
if (xAODwrk)sc = VKalVrtFitFastBase(xAODwrk->tmpListTracks,newvrt.vertex);
else if(RECwork)sc = VKalVrtFitFastBase(RECwork->tmpListTracks,newvrt.vertex);
- if( sc.isFailure() || newvrt.vertex.perp() > m_Rlayer2*2. ) { /* No initial estimation */
+ if( sc.isFailure() || newvrt.vertex.perp() > m_rLayer2*2. ) { /* No initial estimation */
m_fitSvc->setApproximateVertex(PrimVrt.x(), PrimVrt.y(), PrimVrt.z()); /*Use as starting point*/
- if( m_MultiWithPrimary ) m_fitSvc->setApproximateVertex(0., 0., 0.);
+ if( m_multiWithPrimary ) m_fitSvc->setApproximateVertex(0., 0., 0.);
} else {
Amg::Vector3D vDist=newvrt.vertex-PrimVrt.position();
double JetVrtDir = JetDir.Px()*vDist.x() + JetDir.Py()*vDist.y() + JetDir.Pz()*vDist.z();
- if( m_MultiWithPrimary ) JetVrtDir=fabs(JetVrtDir); /* Always positive when primary vertex is seeked for*/
+ if( m_multiWithPrimary ) JetVrtDir=fabs(JetVrtDir); /* Always positive when primary vertex is seeked for*/
if( JetVrtDir>0. ) { /* Good initial estimation */
m_fitSvc->setApproximateVertex(newvrt.vertex.x(),newvrt.vertex.y(),newvrt.vertex.z()); /*Use as starting point*/
}else{
@@ -285,7 +272,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
} if(jv!=tmpN) break; // One vertex is erased. Restart check
} if(iv==tmpN-1) break; // No vertex deleted
}
- if(m_FillHist){ int cvgood=0; for(int iv=0; iv<(int)(*WrkVrtSet).size(); iv++) if((*WrkVrtSet)[iv].Good)cvgood++;
+ if(m_fillHist){ int cvgood=0; for(int iv=0; iv<(int)(*WrkVrtSet).size(); iv++) if((*WrkVrtSet)[iv].Good)cvgood++;
m_hb_rawVrtN->Fill( (float)cvgood, m_w_1); }
//- Try to merge all vertices which have at least 2 common track
for(int iv=0; iv<(int)(*WrkVrtSet).size()-1; iv++ ){
@@ -294,7 +281,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if(!(*WrkVrtSet)[jv].Good ) continue;
if(nTrkCommon( WrkVrtSet, iv, jv)<2) continue;
if( VrtVrtDist((*WrkVrtSet)[iv].vertex,(*WrkVrtSet)[iv].vertexCov,
- (*WrkVrtSet)[jv].vertex,(*WrkVrtSet)[jv].vertexCov) < m_VertexMergeCut) {
+ (*WrkVrtSet)[jv].vertex,(*WrkVrtSet)[jv].vertexCov) < m_vertexMergeCut) {
double probV=0.;
if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, xAODwrk->listJetTracks);
else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, RECwork->listJetTracks);
@@ -319,7 +306,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
// if( (*WrkVrtSet)[iv].SelTrk.size()!=2)continue;
// if( (*WrkVrtSet)[iv].ProjectedVrt <0.) (*WrkVrtSet)[iv].Good =false;
// if( (*WrkVrtSet)[iv].Chi2 > 10.) (*WrkVrtSet)[iv].Good=false;
-// if( (*WrkVrtSet)[iv].vertexMom.M()>VrtBCMassLimit) (*WrkVrtSet)[iv].Good=false; //VK B-tagging specific requirement
+// if( (*WrkVrtSet)[iv].vertexMom.M()>c_vrtBCMassLimit) (*WrkVrtSet)[iv].Good=false; //VK B-tagging specific requirement
// }
//-Remove all bad vertices from the working set
int tmpV=0; while( tmpV<(int)(*WrkVrtSet).size() )if( !(*WrkVrtSet)[tmpV].Good ) { (*WrkVrtSet).erase((*WrkVrtSet).begin()+tmpV);} else {tmpV++;}
@@ -373,18 +360,18 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
//printWrkSet(WrkVrtSet,"Iterat");
double foundMinVrtDst = 1000000.;
- if(foundMaxT<m_TrackDetachCut) foundMinVrtDst = minVrtVrtDist( WrkVrtSet, foundV1, foundV2);
+ if(foundMaxT<m_trackDetachCut) foundMinVrtDst = minVrtVrtDist( WrkVrtSet, foundV1, foundV2);
//Choice of action
- if( foundMaxT<m_TrackDetachCut && foundMinVrtDst<m_VertexMergeCut && nTrkCommon( WrkVrtSet, foundV1, foundV2)){
- //if( foundMaxT<m_TrackDetachCut && foundMinVrtDst<m_VertexMergeCut){
+ if( foundMaxT<m_trackDetachCut && foundMinVrtDst<m_vertexMergeCut && nTrkCommon( WrkVrtSet, foundV1, foundV2)){
+ //if( foundMaxT<m_trackDetachCut && foundMinVrtDst<m_vertexMergeCut){
bool vrtMerged=false; //to check whether something is really merged
- while(foundMinVrtDst<m_VertexMergeCut){
+ while(foundMinVrtDst<m_vertexMergeCut){
if(foundV1<foundV2) { int tmp=foundV1; foundV1=foundV2; foundV2=tmp;} /*Always drop vertex with smallest number*/
double probV=0.;
if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, xAODwrk->listJetTracks);
else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, RECwork->listJetTracks);
- if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<VrtBCMassLimit){ // Good merged vertex found
+ if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){ // Good merged vertex found
double tstDst=JetProjDist(newvrt.vertex, PrimVrt, JetDir);
if(tstDst>0.){ // only positive vertex directions are accepted as merging result
std::swap((*WrkVrtSet)[foundV1],newvrt);
@@ -419,12 +406,12 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
// Final check/merge for close vertices
//
double minDistVV=minVrtVrtDist( WrkVrtSet, foundV1, foundV2); //recalculate VV distances
- while ( minDistVV < m_VertexMergeCut) {
+ while ( minDistVV < m_vertexMergeCut) {
if(foundV1<foundV2) { int tmp=foundV1; foundV1=foundV2; foundV2=tmp;}
double probV=0.;
if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, xAODwrk->listJetTracks);
else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, RECwork->listJetTracks);
- if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<VrtBCMassLimit){ // Good merged vertex found
+ if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){ // Good merged vertex found
double tstDst=JetProjDist(newvrt.vertex, PrimVrt, JetDir);
if(tstDst>0.){ // only positive vertex directions are accepted as merging result
std::swap((*WrkVrtSet)[foundV1],newvrt);
@@ -498,20 +485,19 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if( fabs(curVrt.vertex.z())>650. ){curVrt.Good=false; continue;} //vertex outside Pixel det. For ALL vertices
if(curVrt.SelTrk.size() != 1) continue;
curVrt.Good=false; // Make them bad by default
- if(MomAtVrt(curVrt.TrkAtVrt[0]).Pt() < m_JetPtFractionCut*JetDir.Perp() ) continue; //Low Pt
- if(m_MultiWithOneTrkVrt){ /* 1track vertices left unassigned from good 2tr vertices */
+ if(m_multiWithOneTrkVrt){ /* 1track vertices left unassigned from good 2tr vertices */
VrtVrtDist(PrimVrt,curVrt.vertex, curVrt.vertexCov, Signif3D); //VK non-projected Signif3D is worse
double tmpProb=TMath::Prob( curVrt.Chi2, 1); //Chi2 of the original 2tr vertex
bool trkGood=false;
- if (RECwork)trkGood=Check1TrVertexInPixel(RECwork->listJetTracks[curVrt.SelTrk[0]],curVrt.vertex);
- else if(xAODwrk)trkGood=Check1TrVertexInPixel(xAODwrk->listJetTracks[curVrt.SelTrk[0]],curVrt.vertex);
+ if (RECwork)trkGood=Check1TrVertexInPixel(RECwork->listJetTracks[curVrt.SelTrk[0]],curVrt.vertex,curVrt.vertexCov);
+ else if(xAODwrk)trkGood=Check1TrVertexInPixel(xAODwrk->listJetTracks[curVrt.SelTrk[0]],curVrt.vertex,curVrt.vertexCov);
if(trkGood && tmpProb>0.01){ /* accept only good tracks coming from good 2tr vertex*/
- //if( m_useMaterialRejection && insideMatLayer(curVrt.vertex.x(),curVrt.vertex.y()) ) continue;
+ //if( useMaterialRejection && insideMatLayer(curVrt.vertex.x(),curVrt.vertex.y()) ) continue;
std::vector<double> Impact,ImpactError; double Signif3DP = 0;
if (xAODwrk) Signif3DP=m_fitSvc->VKalGetImpact(xAODwrk->listJetTracks[curVrt.SelTrk[0]],PrimVrt.position(), 1, Impact, ImpactError);
else if(RECwork) Signif3DP=m_fitSvc->VKalGetImpact(RECwork->listJetTracks[curVrt.SelTrk[0]],PrimVrt.position(), 1, Impact, ImpactError);
- if(m_FillHist&&curVrt.vertex.perp()>20.){m_hb_diffPS->Fill( Signif3DP, m_w_1); }
- if( Signif3DP>2.*m_TrkSigCut && Signif3D>m_Sel2VrtSigCut) curVrt.Good=true; // accept only tracks which are far from primary vertex
+ if(m_fillHist&&curVrt.vertex.perp()>20.){m_hb_diffPS->Fill( Signif3DP, m_w_1); }
+ if( Signif3DP>2.*m_trkSigCut && Signif3D>m_sel2VrtSigCut) curVrt.Good=true; // accept only tracks which are far from primary vertex
}
}
}
@@ -548,29 +534,28 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
//-----------------------------------------------------------------------------------------
if(nth==2 && m_useVertexCleaning){
if(RECwork){
- if(!Check2TrVertexInPixel(RECwork->tmpListTracks[0],RECwork->tmpListTracks[1],curVrt.vertex,Signif3D))continue;
+ if(!Check2TrVertexInPixel(RECwork->tmpListTracks[0],RECwork->tmpListTracks[1],curVrt.vertex,curVrt.vertexCov))continue;
}else if(xAODwrk){
- if(!Check2TrVertexInPixel(xAODwrk->tmpListTracks[0],xAODwrk->tmpListTracks[1],curVrt.vertex,Signif3D))continue;
+ if(!Check2TrVertexInPixel(xAODwrk->tmpListTracks[0],xAODwrk->tmpListTracks[1],curVrt.vertex,curVrt.vertexCov))continue;
}
}
//
//--- Check interactions on pixel layers
//
- //if( curVrt.vertex.perp()>m_Rbeampipe-2. && curVrt.detachedTrack<0) {
- if( curVrt.vertex.perp()>m_Rbeampipe-2.) {
+ //if( curVrt.vertex.perp()>m_beampipeR-2. && curVrt.detachedTrack<0) {
+ if( curVrt.vertex.perp()>m_beampipeR-2.) {
double minPt=1.e9; for(int ti=0; ti<nth; ti++) minPt=TMath::Min(minPt,MomAtVrt(curVrt.TrkAtVrt[ti]).Pt());
- double ptLim=TMath::Max(m_hadronIntPtCut,m_JetPtFractionCut*JetDir.Perp());
- if(m_FillHist){
+ if(m_fillHist){
m_hb_totmass2T0->Fill( curVrt.vertexMom.M()-nth*m_massPi, m_w_1);
m_hb_r2d->Fill( curVrt.vertex.perp(), m_w_1);
}
- if(m_useMaterialRejection && minPt<ptLim){
+ if(useMaterialRejection){
if( insideMatLayer(curVrt.vertex.x(),curVrt.vertex.y()) ) continue;
}
//double dR=0; for(int mi=0;mi<nth-1;mi++)for(int mj=mi+1;mj<nth;mj++)
// dR=TMath::Max(dR,MomAtVrt(curVrt.TrkAtVrt[mi]).DeltaR(MomAtVrt(curVrt.TrkAtVrt[mj])));
- //if(m_FillHist)m_hb_deltaRSVPV->Fill(dR,m_w_1);
- //if( m_killHighPtIBLFakes && curVrt.vertex.perp()<m_Rlayer1 && dR<0.015)continue;
+ //if(m_fillHist)m_hb_deltaRSVPV->Fill(dR,m_w_1);
+ //if( m_killHighPtIBLFakes && curVrt.vertex.perp()<m_rLayer1 && dR<0.015)continue;
if(Signif3D<20.) continue;
}
//
@@ -579,7 +564,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
double mass_PiPi = curVrt.vertexMom.M();
double mass_PPi = massV0(curVrt.TrkAtVrt,m_massP,m_massPi);
double mass_EE = massV0(curVrt.TrkAtVrt,m_massE,m_massE);
- if(m_FillHist){ m_hb_massPiPi->Fill( mass_PiPi, m_w_1);
+ if(m_fillHist){ m_hb_massPiPi->Fill( mass_PiPi, m_w_1);
m_hb_massPPi ->Fill( mass_PPi, m_w_1);
if( curVrt.vertex.perp()>20.)m_hb_massEE ->Fill( mass_EE, m_w_1); }
if( fabs(mass_PiPi-m_massK0) < 22.) continue;
@@ -587,15 +572,15 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if( mass_EE < 60. && curVrt.vertex.perp() > 20.) continue;
}
//---
- if(m_FillHist){m_hb_sig3DTot->Fill( Signif3D, m_w_1); }
- if(Signif3D<m_Sel2VrtSigCut)continue; //Main PV-SV distance quality cut
+ if(m_fillHist){m_hb_sig3DTot->Fill( Signif3D, m_w_1); }
+ if(Signif3D<m_sel2VrtSigCut)continue; //Main PV-SV distance quality cut
//-----
// float Dist2D= (*WrkVrtSet)[iv].vertex.perp();
// if(Dist2D<2.){
// double tmpProb=0.;
// if (xAODwrk)tmpProb=FitVertexWithPV( WrkVrtSet, iv, PrimVrt, xAODwrk->listJetTracks);
// else if(RECwork)tmpProb=FitVertexWithPV( WrkVrtSet, iv, PrimVrt, RECwork->listJetTracks);
-// if(m_FillHist){m_hb_trkPtMax->Fill( tmpProb*1.e5, m_w_1); }
+// if(m_fillHist){m_hb_trkPtMax->Fill( tmpProb*1.e5, m_w_1); }
// if(tmpProb>0.01)continue; // Vertex can be associated with PV
// }
//---
@@ -605,7 +590,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
}
//
//--Final cleaning of the 1-track vertices set. Must be behind all other cleanings.
- if(m_MultiWithOneTrkVrt) Clean1TrVertexSet(WrkVrtSet);
+ if(m_multiWithOneTrkVrt) Clean1TrVertexSet(WrkVrtSet);
//------------ Debug
// std::vector<const xAOD::TrackParticle*> tempTrk(0);
// for(auto & iv : (*WrkVrtSet)){ if(iv.Good){for(auto & it : iv.SelTrk)tempTrk.push_back(xAODwrk->listJetTracks.at(it));}}
@@ -620,7 +605,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
nth=iv.SelTrk.size(); if(nth == 0) continue; /* Definitely bad vertices */
Amg::Vector3D tmpVec=iv.vertex-PrimVrt.position();
TLorentzVector Momentum(tmpVec.x(),tmpVec.y(),tmpVec.z(),m_massPi);
- //if(Momentum.DeltaR(JetDir)>m_ConeForTag) iv.Good=false; /* Vertex outside jet cone??? Very bad cut*/
+ //if(Momentum.DeltaR(JetDir)>m_coneForTag) iv.Good=false; /* Vertex outside jet cone??? Very bad cut*/
if( iv.Good) {
nGoodVertices++;
GoodVertices.emplace_back(iv); /* add it */
@@ -629,7 +614,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
}
}
if(nGoodVertices == 0 || (n2trVrt+nNtrVrt)==0){ // No good vertices at all
- delete WrkVrtSet; delete TrkInVrt; if(weit)delete[] weit; if(Solution)delete[] Solution;
+ delete WrkVrtSet; delete TrkInVrt;
return finalVertices;
}
//
@@ -654,7 +639,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if(nGoodVertices>1){
if( GoodVertices[1].vertexMom.M()-GoodVertices[0].vertexMom.M() > 5000.) std::swap( GoodVertices[0], GoodVertices[1] );
}
- if(m_FillHist){m_hb_distVV->Fill( minVrtVrtDist( WrkVrtSet, foundV1, foundV2), m_w_1); }
+ if(m_fillHist){m_hb_distVV->Fill( minVrtVrtDist( WrkVrtSet, foundV1, foundV2), m_w_1); }
//----------------------------------------------------------------------------------
// Nonused tracks for one-track-vertex search
//
@@ -736,7 +721,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
trackPt=pTvsDir(Amg::Vector3D(JetDir.Vect().X(),JetDir.Vect().X(),JetDir.Vect().Z()), GoodVertices[iv].TrkAtVrt[i]);
if(trackPt>trackPtMax)trackPtMax=trackPt;
}
- if( m_FillHist ){
+ if( m_fillHist ){
if(nth==1)m_hb_r1dc->Fill( GoodVertices[iv].vertex.perp(), m_w_1);
if(nth==2)m_hb_r2dc->Fill( GoodVertices[iv].vertex.perp(), m_w_1);
if(nth==3)m_hb_r3dc->Fill( GoodVertices[iv].vertex.perp(), m_w_1);
@@ -785,7 +770,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
nGoodVertices++;
if(nth==1)n1trVrt++;
//-----
- if( iv==0 && m_MultiWithPrimary ) continue; //skip primary vertex if present
+ if( iv==0 && m_multiWithPrimary ) continue; //skip primary vertex if present
VertexMom += GoodVertices[iv].vertexMom;
}
//===================Fake vertex with list of all tracks for optimisation
@@ -799,11 +784,11 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
// } finalVertices.push_back(tmpVertex);
//==============================================
- if(m_FillHist){m_hb_goodvrtN->Fill( nGoodVertices+0.1, m_w_1);
+ if(m_fillHist){m_hb_goodvrtN->Fill( nGoodVertices+0.1, m_w_1);
m_hb_goodvrtN->Fill( n1trVrt+15., m_w_1);}
if(nGoodVertices == 0){
delete WrkVrtSet;
- delete TrkInVrt; if(weit)delete[] weit; if(Solution)delete[] Solution;
+ delete TrkInVrt;
return finalVertices;
}
@@ -822,12 +807,12 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
Results.push_back(0.); //6th - not clear what to use here -> return 0.
Results.push_back(MomentumJet.E()); //7th
- if(m_FillHist){m_hb_ratio->Fill( Results[1], m_w_1); }
- if(m_FillHist){m_hb_totmass->Fill( Results[0], m_w_1); }
- if(m_FillHist){m_hb_nvrt2->Fill( Results[2], m_w_1); }
- if(m_FillHist){m_hb_mom->Fill( MomentumJet.Perp(), m_w_1);}
+ if(m_fillHist){m_hb_ratio->Fill( Results[1], m_w_1); }
+ if(m_fillHist){m_hb_totmass->Fill( Results[0], m_w_1); }
+ if(m_fillHist){m_hb_nvrt2->Fill( Results[2], m_w_1); }
+ if(m_fillHist){m_hb_mom->Fill( MomentumJet.Perp(), m_w_1);}
- delete WrkVrtSet; delete TrkInVrt; if(weit)delete[] weit; if(Solution)delete[] Solution;
+ delete WrkVrtSet; delete TrkInVrt;
return finalVertices;
@@ -886,7 +871,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
newvrt.SelTrk[1]=(*WrkVrtSet)[iv].SelTrk[SelT];
sc = VKalVrtFitFastBase(ListBaseTracks,newvrt.vertex);
if( sc.isFailure() ) continue;
- if( newvrt.vertex.perp() > m_Rlayer2*2. ) newvrt.vertex=Amg::Vector3D(0.,0.,0.);
+ if( newvrt.vertex.perp() > m_rLayer2*2. ) newvrt.vertex=Amg::Vector3D(0.,0.,0.);
m_fitSvc->setApproximateVertex(newvrt.vertex[0],newvrt.vertex[1],newvrt.vertex[2]);
sc=VKalVrtFitBase(ListBaseTracks,
newvrt.vertex,
@@ -929,7 +914,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if(!(*WrkVrtSet)[mtv].Good) continue;
if( std::find((*WrkVrtSet)[mtv].SelTrk.begin(),(*WrkVrtSet)[mtv].SelTrk.end(), Trk1) != (*WrkVrtSet)[mtv].SelTrk.end()){
//double m2Vrt=((*WrkVrtSet)[mtv].vertexMom+(*WrkVrtSet)[i1tv].vertexMom).M(); //VK Commented. M cut in other places
- //if(m2Vrt>VrtBCMassLimit){ (*WrkVrtSet)[i1tv].Good=false; break; } //1Tr + manyTr system is too heavy
+ //if(m2Vrt>c_vrtBCMassLimit){ (*WrkVrtSet)[i1tv].Good=false; break; } //1Tr + manyTr system is too heavy
foundInGoodVrt++; countVT[mtv]++; linkedVrt[i1tv]=mtv; //Linked vertex found
}
}
@@ -1010,7 +995,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
Chi2Red=(*WrkVrtSet)[VertexNumber].Chi2PerTrk.at(itmp); // Normal Chi2 seems the best
if(NTrkInVrt==2){
Chi2Red=(*WrkVrtSet)[VertexNumber].Chi2/2.; //VK 2track vertices with Normal Chi2Red
- if((*WrkVrtSet)[VertexNumber].vertexMom.M()>VrtBCMassLimit)Chi2Red=100.; //VK break immediately very heavy 2tr vertices
+ if((*WrkVrtSet)[VertexNumber].vertexMom.M()>c_vrtBCMassLimit)Chi2Red=100.; //VK break immediately very heavy 2tr vertices
}
double prob_vrt = TMath::Prob( (*WrkVrtSet)[VertexNumber].Chi2, 2*(*WrkVrtSet)[VertexNumber].SelTrk.size()-3);
if( MaxOf < Chi2Red ){
@@ -1034,8 +1019,8 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if (SelectedVertex==v1 && dst2<dst1) SelectedVertex=v2; // Swap to remove the closest to PV vertex
else if(SelectedVertex==v2 && dst1<dst2) SelectedVertex=v1; // Swap to remove the closest to PV vertex
double M1=(*WrkVrtSet)[v1].vertexMom.M(); double M2=(*WrkVrtSet)[v2].vertexMom.M();
- if( M1>VrtBCMassLimit && M2<VrtBCMassLimit ) SelectedVertex=v1;
- if( M1<VrtBCMassLimit && M2>VrtBCMassLimit ) SelectedVertex=v2;
+ if( M1>c_vrtBCMassLimit && M2<c_vrtBCMassLimit ) SelectedVertex=v1;
+ if( M1<c_vrtBCMassLimit && M2>c_vrtBCMassLimit ) SelectedVertex=v2;
}
if( (*WrkVrtSet)[v1].SelTrk.size()+(*WrkVrtSet)[v2].SelTrk.size() > 4){
if( (*WrkVrtSet)[v1].SelTrk.size()==2 && dst1>dst2) SelectedVertex=v2;
@@ -1079,7 +1064,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
(*TrkInVrt)[LeftTrack].erase(it); break;
}
}
- if( (*WrkVrtSet)[SelectedVertex].vertexMom.M()>VrtBCMassLimit)(*WrkVrtSet)[SelectedVertex].Good=false; // Vertex is too heavy
+ if( (*WrkVrtSet)[SelectedVertex].vertexMom.M()>c_vrtBCMassLimit)(*WrkVrtSet)[SelectedVertex].Good=false; // Vertex is too heavy
int ipos=0; if(posInVrtFit==0)ipos=1; // Position of remaining track in previous 2track vertex fit
(*WrkVrtSet)[SelectedVertex].vertexMom=MomAtVrt((*WrkVrtSet)[SelectedVertex].TrkAtVrt[ipos]); //Redefine vertexMom using remaining track
if((*TrkInVrt)[LeftTrack].size()>0)(*WrkVrtSet)[SelectedVertex].Good=false; //Vertex is declared false only if remaining track
@@ -1173,14 +1158,6 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
{
if(!(*WrkVrtSet).at(V1).Good)return -1.; //bad vertex
if(!(*WrkVrtSet).at(V2).Good)return -1.; //bad vertex
- double RMin=TMath::Min((*WrkVrtSet)[V1].vertex.perp(),(*WrkVrtSet)[V2].vertex.perp());
- double RMax=TMath::Max((*WrkVrtSet)[V1].vertex.perp(),(*WrkVrtSet)[V2].vertex.perp());
- if(RMax-RMin>m_SVResolutionR){
- if(RMin<m_RlayerB && m_RlayerB<RMax) return -1.;
- if(RMin<m_Rlayer1 && m_Rlayer1<RMax) return -1.;
- if(RMin<m_Rlayer2 && m_Rlayer2<RMax) return -1.;
- }
-
newvrt.Good=true;
int NTrk_V1=(*WrkVrtSet)[V1].SelTrk.size();
int NTrk_V2=(*WrkVrtSet)[V2].SelTrk.size();
@@ -1369,7 +1346,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
//For possble future use
-/* if( m_MultiWithOneTrkVrt && (!m_MultiWithPrimary) && nGoodVertices<5){ // Addition of one-track vertices is allowed
+/* if( m_multiWithOneTrkVrt && (!m_multiWithPrimary) && nGoodVertices<5){ // Addition of one-track vertices is allowed
double addVrtChi2Cut =5.0;
double addSig3DCut =5.0;
int tmpNTrk=0; if(xAODwrk)tmpNTrk=xAODwrk->listJetTracks.size(); else if(RECwork)tmpNTrk=RECwork->listJetTracks.size();
@@ -1390,7 +1367,7 @@ const double VrtBCMassLimit=6000.; // Mass limit to consider a vertex not comom
if(tmpNBLHits>0){ // accept only tracks with b-layer hit
if (RECwork)Signif3DP = m_fitSvc->VKalGetImpact(RECwork->listJetTracks[atr], PrimVrt.position(), 1, Impact, ImpactError);
else if(xAODwrk)Signif3DP = m_fitSvc->VKalGetImpact(xAODwrk->listJetTracks[atr], PrimVrt.position(), 1, Impact, ImpactError);
- if( Signif3DP > m_TrkSigCut ) nonusedTracks.push_back(atr);
+ if( Signif3DP > m_trkSigCut ) nonusedTracks.push_back(atr);
}
}
}
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
index eec26c5391299373febf0a2696ff94244828f641..f813ae45f756a2fc42c43b5efde7ea39e64e7b25 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// Header include
@@ -21,20 +21,22 @@ namespace InDet{
{
double Pt = sin(ThetaVert)/fabs(PInvVert);
//- Track quality
- if(Pt < m_CutPt) return StatusCode::FAILURE;
-
- if(!m_MultiWithPrimary){ //Must not be used for primary vertex search
- if(fabs(ZVert) > m_CutZVrt/sin(ThetaVert)) return StatusCode::FAILURE;
+ if(Pt < m_cutPt) return StatusCode::FAILURE;
+//std::cout<<" ZVert="<<ZVert<<", "<<sin(ThetaVert)<<'\n';
+//std::cout<<" Chi2="<<Chi2<<'\n';
+//std::cout<<" A0Vert="<<A0Vert<<'\n';
+ if(!m_multiWithPrimary){ //Must not be used for primary vertex search
+ if(fabs(ZVert) > m_cutZVrt) return StatusCode::FAILURE;
}
- if(Chi2 > m_CutChi2) return StatusCode::FAILURE;
- if(fabs(A0Vert) > m_CutA0) return StatusCode::FAILURE;
+ if(Chi2 > m_cutChi2) return StatusCode::FAILURE;
+ if(fabs(A0Vert) > m_cutA0) return StatusCode::FAILURE;
- if(PixelHits < m_CutPixelHits) return StatusCode::FAILURE;
- if(SctHits < m_CutSctHits) return StatusCode::FAILURE;
- if((PixelHits+SctHits) < m_CutSiHits) return StatusCode::FAILURE;
- if(BLayHits < m_CutBLayHits) return StatusCode::FAILURE;
- if(SharedHits > m_CutSharedHits) return StatusCode::FAILURE;
+ if(PixelHits < m_cutPixelHits) return StatusCode::FAILURE;
+ if(SctHits < m_cutSctHits) return StatusCode::FAILURE;
+ if((PixelHits+SctHits) < m_cutSiHits) return StatusCode::FAILURE;
+ if(BLayHits < m_cutBLayHits) return StatusCode::FAILURE;
+ if(SharedHits > m_cutSharedHits) return StatusCode::FAILURE;
return StatusCode::SUCCESS;
}
@@ -69,14 +71,15 @@ namespace InDet{
double CovTrkMtx11 = (*(mPer->covariance()))(0,0);
double CovTrkMtx22 = (*(mPer->covariance()))(1,1);
- if ( CovTrkMtx11 > m_A0TrkErrorCut*m_A0TrkErrorCut ) continue;
- if ( CovTrkMtx22 > m_ZTrkErrorCut*m_ZTrkErrorCut ) continue;
- if( ConeDist(VectPerig,JetDir) > m_ConeForTag ) continue;
+ if ( CovTrkMtx11 > m_a0TrkErrorCut*m_a0TrkErrorCut ) continue;
+ if ( CovTrkMtx22 > m_zTrkErrorCut*m_zTrkErrorCut ) continue;
+ if( ConeDist(VectPerig,JetDir) > m_coneForTag ) continue;
+ if( (*i_ntrk)->pt() > JetDir.Pt() ) continue;
double trkP=1./fabs(VectPerig[4]);
double CovTrkMtx55 = (*(mPer->covariance()))(4,4);
if(trkP>10000.){ double trkPErr=sqrt(CovTrkMtx55)*trkP;
- if(m_FillHist)m_hb_trkPErr->Fill( trkPErr , m_w_1);
+ if(m_fillHist)m_hb_trkPErr->Fill( trkPErr , m_w_1);
if(trkPErr>0.5) continue; }
long int PixelHits = 3;
@@ -110,49 +113,8 @@ namespace InDet{
AnalysisUtils::Sort::pT (&SelectedTracks); // no equivalent for TrkTrack yet...
return NPrimTrk;
}
-//
-// Simplified TrackParticle selection to get reliably reconstructed tracks for jet momentum estimation
-//
- int InDetVKalVxInJetTool::SelGoodTrkParticleRelax( const std::vector<const Rec::TrackParticle*>& InpTrk,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- std::vector<const Rec::TrackParticle*> & SelectedTracks)
- const
- {
-
- std::vector<const Rec::TrackParticle*>::const_iterator i_ntrk;
- AmgVector(5) VectPerig; VectPerig<<0.,0.,0.,0.,0.;
- const Trk::Perigee* mPer;
- std::vector<double> Impact,ImpactError;
- int NPrimTrk=0;
- for (i_ntrk = InpTrk.begin(); i_ntrk < InpTrk.end(); ++i_ntrk) {
-//
-//-- MeasuredPerigee in TrackParticle
-//
- mPer=GetPerigee( (*i_ntrk) ) ;
- if( mPer == NULL ){ continue; }
- VectPerig = mPer->parameters();
- if( ConeDist(VectPerig,JetDir) > m_ConeForTag ) continue;
-//----------------------------------- Summary tools
- const Trk::TrackSummary* testSum = (*i_ntrk)->trackSummary();
- long int PixelHits = (long int) testSum->get(Trk::numberOfPixelHits);
- long int SctHits = (long int) testSum->get(Trk::numberOfSCTHits);
- if(PixelHits < 0 ) PixelHits=0;
- if(SctHits < 0 ) SctHits=0;
- double ImpactSignif = m_fitSvc->VKalGetImpact((*i_ntrk), PrimVrt.position(), 1, Impact, ImpactError);
- if(fabs(Impact[0]) > m_CutA0) continue;
- if(fabs(Impact[1]) > m_CutZVrt/sin(VectPerig[3])) continue;
- if(PixelHits < m_CutPixelHits) continue;
- if(SctHits < m_CutSctHits) continue;
- if((PixelHits+SctHits) < m_CutSiHits) continue;
- if(ImpactSignif < 3.)NPrimTrk += 1;
- SelectedTracks.push_back(*i_ntrk);
- }
- AnalysisUtils::Sort::pT (&SelectedTracks); // no equivalent for TrkTrack yet...
- return NPrimTrk;
- }
//==============================================================================================================
// xAOD based stuff
//
@@ -165,6 +127,7 @@ namespace InDet{
std::vector<const xAOD::TrackParticle*>::const_iterator i_ntrk;
std::vector<double> Impact,ImpactError;
+ std::map<double,const xAOD::TrackParticle*> orderedTrk;
int NPrimTrk=0;
for (i_ntrk = InpTrk.begin(); i_ntrk < InpTrk.end(); ++i_ntrk) {
//
@@ -183,13 +146,14 @@ namespace InDet{
- if ( CovTrkMtx11 > m_A0TrkErrorCut*m_A0TrkErrorCut ) continue;
- if ( CovTrkMtx22 > m_ZTrkErrorCut*m_ZTrkErrorCut ) continue;
- if ( ConeDist(VectPerig,JetDir) > m_ConeForTag ) continue;
+ if ( CovTrkMtx11 > m_a0TrkErrorCut*m_a0TrkErrorCut ) continue;
+ if ( CovTrkMtx22 > m_zTrkErrorCut*m_zTrkErrorCut ) continue;
+ if ( ConeDist(VectPerig,JetDir) > m_coneForTag ) continue;
+ if( (*i_ntrk)->pt() > JetDir.Pt() ) continue;
double trkP=1./fabs(VectPerig[4]);
if(trkP>10000.){ double trkPErr=sqrt(CovTrkMtx55)*trkP;
- if(m_FillHist)m_hb_trkPErr->Fill( trkPErr , m_w_1);
+ if(m_fillHist)m_hb_trkPErr->Fill( trkPErr , m_w_1);
if(trkPErr>0.5) continue; }
uint8_t PixelHits,SctHits,BLayHits;
@@ -220,11 +184,12 @@ namespace InDet{
//std::cout<<"NwInnerM="<<(long int)InmHits<<", "<<(long int)InmSharedH<<", "<<(long int)InmSplitH<<", "<<(long int)InmOutlier<<'\n';
- double ImpactSignif = m_fitSvc->VKalGetImpact((*i_ntrk), PrimVrt.position(), 1, Impact, ImpactError);
- double ImpactA0=VectPerig[0]; // Temporary
- double ImpactZ=VectPerig[1]-PrimVrt.position().z(); // Temporary
- ImpactA0=Impact[0];
- ImpactZ=Impact[1];
+ m_fitSvc->VKalGetImpact((*i_ntrk), PrimVrt.position(), 1, Impact, ImpactError);
+ //double ImpactA0=VectPerig[0]; // Temporary
+ //double ImpactZ=VectPerig[1]-PrimVrt.position().z(); // Temporary
+ double ImpactA0=Impact[0];
+ double ImpactZ=Impact[1];
+ if(m_fillHist){ m_hb_trkD0->Fill( ImpactA0, m_w_1); }
//---- Improved cleaning
///// if(PixelHits<=2 && ( outPixHits || splPixHits )) continue; //VK Bad idea at high Pt!
if(fabs((*i_ntrk)->eta())>2. ) {
@@ -233,59 +198,31 @@ namespace InDet{
else {if(PixelHits)PixelHits -=1;} // 3-layer pixel detector
}
if(fabs((*i_ntrk)->eta())>1.65) if(SctHits)SctHits -=1;
+//----Anti-pileup cut
+// double SignifR = Impact[0]/ sqrt(ImpactError[0]);
+// if(fabs(SignifR) < m_AntiPileupSigRCut) { // cut against tracks from pileup vertices
+// double SignifZ = Impact[1]/ sqrt(ImpactError[2]);
+// if(SignifZ > 1.+m_AntiPileupSigZCut ) continue;
+// if(SignifZ < 1.-m_AntiPileupSigZCut ) continue;
+// }
+//---- Use classificator to remove Pileup+Interactions
+ std::vector<float> trkRank=m_trackClassificator->trkTypeWgts(*i_ntrk, PrimVrt, JetDir);
+ if(trkRank[2] > m_antiGarbageCut)continue;
//----
StatusCode sc = CutTrk( VectPerig[4] , VectPerig[3],
ImpactA0 , ImpactZ, trkChi2,
- PixelHits, SctHits, SharedHits, BLayHits); //
+ PixelHits, SctHits, SharedHits, BLayHits);
if( sc.isFailure() ) continue;
- if(ImpactSignif < 3.)NPrimTrk += 1;
- SelectedTracks.push_back(*i_ntrk);
+ //double rankBTrk=RankBTrk((*i_ntrk)->pt(),JetDir.Perp(),ImpactSignif);
+ if(trkRank[1]>0.5)NPrimTrk += 1;
+ orderedTrk[trkRank[0]]= *i_ntrk;
}
- AnalysisUtils::Sort::pT (&SelectedTracks); // no equivalent for TrkTrack yet...
+//---- Order tracks according to ranks
+ std::map<double,const xAOD::TrackParticle*>::reverse_iterator rt=orderedTrk.rbegin();
+ SelectedTracks.resize(orderedTrk.size());
+ for ( int cntt=0; rt!=orderedTrk.rend(); ++rt,++cntt) {SelectedTracks[cntt]=(*rt).second;}
+ m_NRefPVTrk=std::max(NPrimTrk,1); // VK set reference multiplicity here
return NPrimTrk;
}
-//
-// Simplified xAOD::TrackParticle selection to get reliably reconstructed tracks for jet momentum estimation
-//
- int InDetVKalVxInJetTool::SelGoodTrkParticleRelax( const std::vector<const xAOD::TrackParticle*>& InpTrk,
- const xAOD::Vertex & PrimVrt,
- const TLorentzVector & JetDir,
- std::vector<const xAOD::TrackParticle*>& SelectedTracks)
- const
- {
- std::vector<const xAOD::TrackParticle*>::const_iterator i_ntrk;
- std::vector<double> Impact,ImpactError;
- int NPrimTrk=0;
- for (i_ntrk = InpTrk.begin(); i_ntrk < InpTrk.end(); ++i_ntrk) {
-//
-//-- MeasuredPerigee in xAOD::TrackParticle
-//
- const Trk::Perigee mPer=(*i_ntrk)->perigeeParameters() ;
- AmgVector(5) VectPerig = mPer.parameters();
- if( ConeDist(VectPerig,JetDir) > m_ConeForTag ) continue;
-//----------------------------------- Summary tools
- uint8_t PixelHits,SctHits;
- if( !((*i_ntrk)->summaryValue(PixelHits,xAOD::numberOfPixelHits)) ) PixelHits=0;
- if( !((*i_ntrk)->summaryValue( SctHits,xAOD::numberOfSCTHits)) ) SctHits=0;
-
-
- double ImpactSignif = m_fitSvc->VKalGetImpact((*i_ntrk), PrimVrt.position(), 1, Impact, ImpactError);
- if(fabs(Impact[0]) > m_CutA0) continue;
- if(fabs(Impact[1]) > m_CutZVrt/sin(VectPerig[3])) continue;
-
- int currCutPixelHits=m_CutPixelHits; if(fabs((*i_ntrk)->eta())>2. )currCutPixelHits +=1;
- int currCutSctHits =m_CutSctHits; if(fabs((*i_ntrk)->eta())>1.65)currCutSctHits +=1;
-
- if(PixelHits < currCutPixelHits) continue;
- if(SctHits < currCutSctHits) continue;
- if((PixelHits+SctHits) < m_CutSiHits) continue;
- if(ImpactSignif < 3.)NPrimTrk += 1;
- SelectedTracks.push_back(*i_ntrk);
- }
- AnalysisUtils::Sort::pT (&SelectedTracks); // no equivalent for TrkTrack yet...
- return NPrimTrk;
- }
-
-
}//end namespace
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..49a0c844966a441ef511188f7690496bb4749876
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx
@@ -0,0 +1,170 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "InDetVKalVxInJetTool/InDetTrkInJetType.h"
+#include "TMVA/MethodBase.h"
+#include "TMVA/Reader.h"
+#include "PathResolver/PathResolver.h"
+//
+//-------------------------------------------------
+namespace InDet {
+//
+//Constructor--------------------------------------------------------------
+InDetTrkInJetType::InDetTrkInJetType(const std::string& type,
+ const std::string& name,
+ const IInterface* parent):
+ AthAlgTool(type,name,parent),
+ m_tmvaReader(0),
+ m_trkMinPtCut(700.),
+ m_d0_limLow(-3.),
+ m_d0_limUpp( 5.),
+ m_Z0_limLow(-8.),
+ m_Z0_limUpp(12.),
+ m_calibFileName("TrackClassif_3cl.v01.xml"),
+ m_fitterSvc("Trk::TrkVKalVrtFitter/VertexFitterTool",this)
+ {
+ declareInterface<IInDetTrkInJetType>(this);
+ declareProperty("trkMinPt", m_trkMinPtCut , "Minimal track Pt cut" );
+ declareProperty("d0_limLow", m_d0_limLow , "Low d0 impact cut" );
+ declareProperty("d0_limUpp", m_d0_limUpp , "Upper d0 impact cut" );
+ declareProperty("Z0_limLow", m_Z0_limLow , "Low Z0 impact cut" );
+ declareProperty("Z0_limUpp", m_Z0_limUpp , "Upper Z0 impact cut" );
+ }
+
+//Destructor---------------------------------------------------------------
+ InDetTrkInJetType::~InDetTrkInJetType(){
+ delete m_tmvaReader;
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<< "InDetTrkInJetType destructor called" << endmsg;
+ }
+
+//Initialize---------------------------------------------------------------
+ StatusCode InDetTrkInJetType::initialize(){
+ m_initialised = 0;
+ m_tmvaReader = new TMVA::Reader();
+ //m_tmvaReader->AddVariable( "prbS", &m_prbS );
+ m_tmvaReader->AddVariable( "Sig3D", &m_Sig3D );
+ m_tmvaReader->AddVariable( "prbP", &m_prbP );
+ m_tmvaReader->AddVariable( "pTvsJet", &m_pTvsJet );
+ //m_tmvaReader->AddVariable( "prodTJ", &m_prodTJ );
+ m_tmvaReader->AddVariable( "d0", &m_d0 );
+ m_tmvaReader->AddVariable( "SigR", &m_SigR );
+ m_tmvaReader->AddVariable( "SigZ", &m_SigZ );
+ m_tmvaReader->AddVariable( "ptjet", &m_ptjet );
+ m_tmvaReader->AddVariable( "ibl" , &m_ibl );
+ m_tmvaReader->AddVariable( "bl" , &m_bl );
+ m_tmvaReader->AddVariable( "etajet", &m_etajet );
+ //m_tmvaReader->AddVariable( "vChi2", &m_vChi2 );
+//
+//-- Calibration file
+//
+// std::string fullPathToFile = PathResolverFindCalibFile("InDetVKalVxInJetTool/TrackClassif_3cl.v01.xml");
+ std::string fullPathToFile = PathResolverFindCalibFile("InDetVKalVxInJetTool/"+m_calibFileName);
+ if(fullPathToFile != ""){
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) <<"TrackClassification calibration file" << fullPathToFile << endmsg;
+ m_tmvaReader->BookMVA("BDTG", fullPathToFile);
+ }else{
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) <<"Error! No calibration for TrackClassification found." << endmsg;
+ return StatusCode::SUCCESS;
+ }
+ //-------
+ if (m_fitterSvc.retrieve().isFailure()) {
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "Could not find Trk::TrkVKalVrtFitter" << endmsg;
+ return StatusCode::SUCCESS;
+ } else {
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "InDetTrkInJetTool TrkVKalVrtFitter found" << endmsg;
+ }
+ m_fitSvc = dynamic_cast<Trk::TrkVKalVrtFitter*>(&(*m_fitterSvc));
+ if(!m_fitSvc){
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<<" No implemented Trk::ITrkVKalVrtFitter interface" << endmsg;
+ return StatusCode::SUCCESS;
+ }
+ m_initialised = 1; // Tool is initialised successfully.
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode InDetTrkInJetType::finalize()
+ {
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) <<"InDetTrkInJetType finalize()" << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ std::vector<float> InDetTrkInJetType::trkTypeWgts(const Rec::TrackParticle *, const xAOD::Vertex &, const TLorentzVector &)
+ { return std::vector<float>(3,0.); }
+
+ std::vector<float> InDetTrkInJetType::trkTypeWgts(const xAOD::TrackParticle * Trk, const xAOD::Vertex & PV, const TLorentzVector & Jet)
+ {
+ std::vector<float> safeReturn(3,0.);
+ if( !m_initialised ) return safeReturn;
+ if(Jet.Perp()>2500000.)return safeReturn;
+ std::vector<double> Impact,ImpactError;
+ m_Sig3D=m_fitSvc->VKalGetImpact(Trk, PV.position(), 1, Impact, ImpactError);
+ AmgVector(5) tmpPerigee = Trk->perigeeParameters().parameters();
+ if( sin(tmpPerigee[2]-Jet.Phi())*Impact[0] < 0 ){ Impact[0] = -fabs(Impact[0]);}
+ else{ Impact[0] = fabs(Impact[0]);}
+ if( (tmpPerigee[3]-Jet.Theta())*Impact[1] < 0 ){ Impact[1] = -fabs(Impact[1]);}
+ else{ Impact[1] = fabs(Impact[1]);}
+ double SignifR = Impact[0]/ sqrt(ImpactError[0]);
+ double SignifZ = Impact[1]/ sqrt(ImpactError[2]);
+ double trkSignif = sqrt( (SignifR+0.6)*(SignifR+0.6) + (SignifZ+0.0)*(SignifZ+0.0) );
+//---
+ if(Impact[0]<m_d0_limLow || Impact[0]>m_d0_limUpp) return safeReturn;
+ if(Impact[0]<m_Z0_limLow || Impact[0]>m_Z0_limUpp) return safeReturn;
+ if( sqrt(SignifR*SignifR +SignifZ*SignifZ) < 1.) return safeReturn;
+//---IBL/BL hits
+ int hitIBL=0, hitBL=0;
+ uint8_t IBLhit,BLhit,IBLexp,BLexp;
+ if(!Trk->summaryValue( IBLhit, xAOD::numberOfInnermostPixelLayerHits) ) IBLhit = 0;
+ if(!Trk->summaryValue( BLhit, xAOD::numberOfNextToInnermostPixelLayerHits) ) BLhit = 0;
+ if(!Trk->summaryValue( IBLexp, xAOD::expectInnermostPixelLayerHit) ) IBLexp = 0;
+ if(!Trk->summaryValue( BLexp, xAOD::expectNextToInnermostPixelLayerHit) ) BLexp = 0;
+ hitIBL=IBLhit; if( IBLexp==0 ) hitIBL=-1;
+ hitBL = BLhit; if( BLexp==0 ) hitBL =-1;
+/*---PV constraint (doesn't improve rejection in first try)
+ Amg::Vector3D FitVrt;
+ TLorentzVector Momentum;
+ long int Charge=0;
+ std::vector<double> ErrorMatrix, Chi2PerTrk;
+ std::vector< std::vector<double> > TrkAtVrt;
+ std::vector<const xAOD::TrackParticle *> TrkForFit(1,Trk);
+ std::vector<const xAOD::NeutralParticle*> netralDummy(0);
+ m_fitSvc->setDefault(); //Reset VKalVrt settings
+ std::vector<float> covPV=PV.covariance();
+ m_fitSvc->setVertexForConstraint(PV.x(),PV.y(),PV.z());
+ m_fitSvc->setCovVrtForConstraint(covPV[0],covPV[1],covPV[2],covPV[3],covPV[4],covPV[5]);
+ m_fitSvc->setCnstType(6); // Set primary vertex constraint
+ StatusCode sc=m_fitSvc->VKalVrtFit( TrkForFit, netralDummy, FitVrt, Momentum, Charge, ErrorMatrix, Chi2PerTrk, TrkAtVrt, Chi2);
+ if(sc.isFailure())Chi2=exp(11.);
+ if(Chi2>exp(11.))Chi2=exp(11.);
+*/
+//====================== BDT weights
+ double coeffPt=10.;
+ double pfrac=(Trk->pt()-m_trkMinPtCut)/sqrt(Jet.Perp());
+ m_prbP= pfrac/(coeffPt+pfrac);
+//---
+ double coeffSig=1.0;
+ if(trkSignif<coeffSig) return safeReturn;
+ m_prbS=(trkSignif-coeffSig)/trkSignif;
+ if(m_prbS<0.) return safeReturn;
+//---
+ m_d0=Impact[0];
+ m_SigZ=SignifZ;
+ m_SigR=SignifR;
+//---
+ m_ptjet=Jet.Perp();
+ m_etajet=fabs(Jet.Eta());
+//---
+ m_ibl = (float)hitIBL;
+ m_bl = (float)hitBL;
+//---
+ TLorentzVector TLV;
+ TLV.SetPtEtaPhiE(Trk->pt(),Trk->eta(),Trk->phi(),Trk->e());
+ m_pTvsJet=TLV.Perp(Jet.Vect());
+//---
+ TLorentzVector normJ; normJ.SetPtEtaPhiM(1.,Jet.Eta(),Jet.Phi(),0.);
+ m_prodTJ=sqrt(TLV.Dot(normJ));
+ return m_tmvaReader->EvaluateMulticlass("BDTG");
+
+ }
+
+}// close namespace
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
index 522aed13051fc19a3a723baf82db9991a51f1220..f606e53685d0deb6d435bf48d68c248a142604e9 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// Header include
@@ -11,14 +11,14 @@
#include "GaudiKernel/ITHistSvc.h"
#include "TH1D.h"
#include "TH2D.h"
+#include "TTree.h"
#include "TProfile.h"
#include "TMath.h"
//
//-------------------------------------------------
// Other stuff
-//
-#include<iostream>
+//#include<iostream>
//extern "C" {
// float chisin_(const float & prob, const long int& ndf);
@@ -31,61 +31,52 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
const std::string& name,
const IInterface* parent):
AthAlgTool(type,name,parent),
- m_CutSctHits(4),
- m_CutPixelHits(1),
- m_CutSiHits(7),
- m_CutBLayHits(0),
- m_CutSharedHits(1),
- m_CutPt(700.),
- m_CutZVrt(25.),
- m_CutA0(5.),
- m_CutChi2(5.),
- m_SecTrkChi2Cut(10.),
- m_ConeForTag(0.4),
- m_Sel2VrtChi2Cut(4.5),
- m_Sel2VrtSigCut(3.0),
- m_TrkSigCut(2.0),
- m_TrkSigNTrkDep(0.20),
- m_TrkSigSumCut(2.),
- m_A0TrkErrorCut(1.0),
- m_ZTrkErrorCut(5.0),
- m_AntiPileupSigRCut(2.0),
- m_AntiPileupSigZCut(6.0),
- m_AntiFake2trVrtCut(0.5),
- m_JetPtFractionCut(0.01),
- m_TrackInJetNumberLimit(25),
- m_pseudoSigCut(3.),
- m_hadronIntPtCut(5000.),
- m_FillHist(false),
+ m_cutSctHits(4),
+ m_cutPixelHits(1),
+ m_cutSiHits(7),
+ m_cutBLayHits(0),
+ m_cutSharedHits(1),
+ m_cutPt(700.),
+ m_cutZVrt(15.),
+ m_cutA0(5.),
+ m_cutChi2(5.),
+ m_secTrkChi2Cut(10.),
+ m_coneForTag(0.4),
+ m_sel2VrtChi2Cut(10.0),
+ m_sel2VrtSigCut(4.0),
+ m_trkSigCut(2.0),
+ m_a0TrkErrorCut(1.0),
+ m_zTrkErrorCut(5.0),
+ m_cutHFClass(0.1),
+ m_antiGarbageCut(0.85),
+ m_fillHist(false),
m_existIBL(true),
m_RobustFit(1),
- m_Xbeampipe (0.),
- m_Ybeampipe (0.),
- m_XlayerB (0.),
- m_YlayerB (0.),
- m_Xlayer1 (0.),
- m_Ylayer1 (0.),
- m_Xlayer2 (0.),
- m_Ylayer2 (0.),
- m_Rbeampipe (0.), //Correct values are filled
- m_RlayerB (0.), // in jobO or initialize()
- m_Rlayer1 (0.),
- m_Rlayer2 (0.),
- m_SVResolutionR(3.),
- m_useMaterialRejection(true),
- m_useVertexCleaning(true),
- m_MassType (1),
- m_MultiVertex(false),
- m_MultiWithPrimary(false),
+ m_beampipeX (0.),
+ m_beampipeY (0.),
+ m_xLayerB (0.),
+ m_yLayerB (0.),
+ m_xLayer1 (0.),
+ m_yLayer1 (0.),
+ m_xLayer2 (0.),
+ m_yLayer2 (0.),
+ m_beampipeR (0.), //Correct values are filled
+ m_rLayerB (0.), // in jobO or initialize()
+ m_rLayer1 (0.),
+ m_rLayer2 (0.),
+ m_useVertexCleaning(false),
+ m_multiVertex(false),
+ m_multiWithPrimary(false),
m_getNegativeTail(false),
m_getNegativeTag(false),
- m_MultiWithOneTrkVrt(true),
- //m_killHighPtIBLFakes(false),
- m_VertexMergeCut(3.),
- m_TrackDetachCut(6.),
- m_fitterSvc("Trk::TrkVKalVrtFitter/VertexFitterTool",this)
+ m_multiWithOneTrkVrt(true),
+ m_vertexMergeCut(3.),
+ m_trackDetachCut(6.),
+ m_fitterSvc("Trk::TrkVKalVrtFitter/VertexFitterTool",this),
+// m_useMaterialRejection(true),
// m_materialMap ("InDet::InDetMaterialRejTool", this)
// m_fitSvc("Trk::TrkVKalVrtFitter/VKalVrtFitter",this)
+ m_trackClassificator("InDet::InDetTrkInJetType",this)
{
//
// Declare additional interface
@@ -94,67 +85,57 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
// Properties
//
//
- declareProperty("CutSctHits", m_CutSctHits , "Remove track is it has less SCT hits" );
- declareProperty("CutPixelHits", m_CutPixelHits, "Remove track is it has less Pixel hits");
- declareProperty("CutSiHits", m_CutSiHits, "Remove track is it has less Pixel+SCT hits" );
- declareProperty("CutBLayHits", m_CutBLayHits, "Remove track is it has less B-layer hits" );
- declareProperty("CutSharedHits", m_CutSharedHits,"Reject final 2tr vertices if tracks have shared hits" );
-
- declareProperty("CutPt", m_CutPt, "Track Pt selection cut" );
- declareProperty("CutA0", m_CutA0, "Track A0 selection cut" );
- declareProperty("CutZVrt", m_CutZVrt, "Track Z impact selection cut");
- declareProperty("ConeForTag", m_ConeForTag,"Cone around jet direction for track selection");
- declareProperty("CutChi2", m_CutChi2, "Track Chi2 selection cut" );
- declareProperty("TrkSigCut", m_TrkSigCut, "Track 3D impact significance w/r primary vertex" );
- declareProperty("TrkSigSumCut", m_TrkSigSumCut, "Sum of 3D track significances cut for 2tr vertex search");
- declareProperty("TrkSigNTrkDep", m_TrkSigNTrkDep, "NTrack in jet dependent increase of TrkSigCut and TrkSigSumCut");
- declareProperty("SecTrkChi2Cut", m_SecTrkChi2Cut,"Track - common secondary vertex association cut. Single Vertex Finder only");
-
- declareProperty("A0TrkErrorCut", m_A0TrkErrorCut, "Track A0 error cut" );
- declareProperty("ZTrkErrorCut", m_ZTrkErrorCut, "Track Z impact error cut" );
-
- declareProperty("Sel2VrtChi2Cut", m_Sel2VrtChi2Cut, "Cut on Chi2 of 2-track vertex for initial selection" );
- declareProperty("Sel2VrtSigCut", m_Sel2VrtSigCut, "Cut on significance of 3D distance between initial 2-track vertex and PV" );
- declareProperty("AntiPileupSigRCut", m_AntiPileupSigRCut, "Remove tracks with low Rphi and big Z impacts presumably coming from pileup" );
- declareProperty("AntiPileupSigZCut", m_AntiPileupSigZCut, "Remove tracks with low Rphi and big Z impacts presumably coming from pileup" );
- declareProperty("AntiFake2trVrtCut", m_AntiFake2trVrtCut, "Cut to reduce fake 2-track vertices contribution.Single Vertex Finder only" );
- declareProperty("JetPtFractionCut", m_JetPtFractionCut, "Reduce high Pt fakes. Jet HLV input is mandatory, direction is not enough. Multi and single vertex versions are affected" );
- declareProperty("TrackInJetNumberLimit", m_TrackInJetNumberLimit, " Use only limited number of highest pT tracks in jet for vertex search" );
- declareProperty("PseudoSigCut", m_pseudoSigCut, " Cut on track impact significance for pseudo-vertex search" );
- declareProperty("HadronIntPtCut", m_hadronIntPtCut, "Pt cut to select hadronic interactions" );
-
- declareProperty("FillHist", m_FillHist, "Fill technical histograms" );
+ declareProperty("CutSctHits", m_cutSctHits , "Remove track is it has less SCT hits" );
+ declareProperty("CutPixelHits", m_cutPixelHits, "Remove track is it has less Pixel hits");
+ declareProperty("CutSiHits", m_cutSiHits, "Remove track is it has less Pixel+SCT hits" );
+ declareProperty("CutBLayHits", m_cutBLayHits, "Remove track is it has less B-layer hits" );
+ declareProperty("CutSharedHits", m_cutSharedHits,"Reject final 2tr vertices if tracks have shared hits" );
+
+ declareProperty("CutPt", m_cutPt, "Track Pt selection cut" );
+ declareProperty("CutA0", m_cutA0, "Track A0 selection cut" );
+ declareProperty("CutZVrt", m_cutZVrt, "Track Z impact selection cut");
+ declareProperty("ConeForTag", m_coneForTag,"Cone around jet direction for track selection");
+ declareProperty("CutChi2", m_cutChi2, "Track Chi2 selection cut" );
+ declareProperty("TrkSigCut", m_trkSigCut, "Track 3D impact significance w/r primary vertex" );
+ declareProperty("SecTrkChi2Cut", m_secTrkChi2Cut,"Track - common secondary vertex association cut. Single Vertex Finder only");
+
+ declareProperty("A0TrkErrorCut", m_a0TrkErrorCut, "Track A0 error cut" );
+ declareProperty("ZTrkErrorCut", m_zTrkErrorCut, "Track Z impact error cut" );
+ declareProperty("CutHFClass", m_cutHFClass, "Cut on HF classification weight" );
+ declareProperty("AntiGarbageCut", m_antiGarbageCut, "Cut on Garbage classification weight for removal" );
+
+ declareProperty("Sel2VrtChi2Cut", m_sel2VrtChi2Cut, "Cut on Chi2 of 2-track vertex for initial selection" );
+ declareProperty("Sel2VrtSigCut", m_sel2VrtSigCut, "Cut on significance of 3D distance between initial 2-track vertex and PV" );
+
+ declareProperty("FillHist", m_fillHist, "Fill technical histograms" );
declareProperty("ExistIBL", m_existIBL, "Inform whether 3-layer or 4-layer detector is used " );
declareProperty("RobustFit", m_RobustFit, "Use vertex fit with RobustFit functional(VKalVrt) for common secondary vertex fit" );
- declareProperty("Xbeampipe", m_Xbeampipe);
- declareProperty("Ybeampipe", m_Ybeampipe);
- declareProperty("XlayerB", m_XlayerB );
- declareProperty("YlayerB", m_YlayerB );
- declareProperty("Xlayer1", m_Xlayer1 );
- declareProperty("Ylayer1", m_Ylayer1 );
- declareProperty("Xlayer2", m_Xlayer2 );
- declareProperty("Ylayer2", m_Ylayer2 );
- declareProperty("Rbeampipe", m_Rbeampipe);
- declareProperty("RlayerB", m_RlayerB );
- declareProperty("Rlayer1", m_Rlayer1 );
- declareProperty("Rlayer2", m_Rlayer2 );
-
- declareProperty("SVResolutionR", m_SVResolutionR, "Radial resolution of SVs. Needed for correct pixel layers crossing checks" );
- declareProperty("useMaterialRejection", m_useMaterialRejection, "Reject vertices from hadronic interactions in detector material" );
+ declareProperty("Xbeampipe", m_beampipeX);
+ declareProperty("Ybeampipe", m_beampipeY);
+ declareProperty("XlayerB", m_xLayerB );
+ declareProperty("YlayerB", m_yLayerB );
+ declareProperty("Xlayer1", m_xLayer1 );
+ declareProperty("Ylayer1", m_yLayer1 );
+ declareProperty("Xlayer2", m_xLayer2 );
+ declareProperty("Ylayer2", m_yLayer2 );
+ declareProperty("Rbeampipe", m_beampipeR);
+ declareProperty("RlayerB", m_rLayerB );
+ declareProperty("Rlayer1", m_rLayer1 );
+ declareProperty("Rlayer2", m_rLayer2 );
+
+// declareProperty("useMaterialRejection", m_useMaterialRejection, "Reject vertices from hadronic interactions in detector material" );
declareProperty("useVertexCleaning", m_useVertexCleaning, "Clean vertices by requiring pixel hit presence according to vertex position" );
- declareProperty("MassType", m_MassType, "Type of vertex mass returned by finder. Single Vertex Finder only!" );
- declareProperty("MultiVertex", m_MultiVertex, "Run Multiple Secondary Vertices in jet finder" );
- declareProperty("MultiWithPrimary", m_MultiWithPrimary, "Find Multiple Secondary Vertices + primary vertex in jet. MultiVertex Finder only!" );
- declareProperty("MultiWithOneTrkVrt", m_MultiWithOneTrkVrt,"Allow one-track-vertex addition to already found secondary vertices. MultiVertex Finder only! ");
- //declareProperty("KillHighPtIBLFakes", m_killHighPtIBLFakes,"Remove fake vertices produced by tracking. MultiVertex Finder only! ");
+ declareProperty("MultiVertex", m_multiVertex, "Run Multiple Secondary Vertices in jet finder" );
+ declareProperty("MultiWithPrimary", m_multiWithPrimary, "Find Multiple Secondary Vertices + primary vertex in jet. MultiVertex Finder only!" );
+ declareProperty("MultiWithOneTrkVrt", m_multiWithOneTrkVrt,"Allow one-track-vertex addition to already found secondary vertices. MultiVertex Finder only! ");
declareProperty("getNegativeTail", m_getNegativeTail, "Allow secondary vertex behind the primary one (negative) w/r jet direction (not for multivertex!)" );
declareProperty("getNegativeTag", m_getNegativeTag, "Return ONLY negative secondary vertices (not for multivertex!)" );
- declareProperty("VertexMergeCut", m_VertexMergeCut, "To allow vertex merging for MultiVertex Finder" );
- declareProperty("TrackDetachCut", m_TrackDetachCut, "To allow track from vertex detachment for MultiVertex Finder" );
+ declareProperty("VertexMergeCut", m_vertexMergeCut, "To allow vertex merging for MultiVertex Finder" );
+ declareProperty("TrackDetachCut", m_trackDetachCut, "To allow track from vertex detachment for MultiVertex Finder" );
declareProperty("VertexFitterTool", m_fitterSvc);
// declareProperty("MaterialMap", m_materialMap);
@@ -170,6 +151,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
m_WorkArray = 0;
m_compatibilityGraph = nullptr;
m_instanceName=name;
+ m_curTup = 0;
}
@@ -177,9 +159,9 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
InDetVKalVxInJetTool::~InDetVKalVxInJetTool(){
//MsgStream log( msgSvc(), name() ) ;
//log << MSG::DEBUG << "InDetVKalVxInJetTool destructor called" << endmsg;
- if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<< "InDetVKalVxInJetTool destructor called" << endmsg;
if(m_WorkArray) delete m_WorkArray;
if(m_compatibilityGraph)delete m_compatibilityGraph;
+ if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG)<< "InDetVKalVxInJetTool destructor called" << endmsg;
}
//Initialize---------------------------------------------------------------
@@ -228,22 +210,22 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
//------------------------------------------
// Chose whether IBL is installed
if(m_existIBL){ // 4-layer pixel detector
- if( m_Rbeampipe==0.) m_Rbeampipe=24.0;
- if( m_RlayerB ==0.) m_RlayerB =34.0;
- if( m_Rlayer1 ==0.) m_Rlayer1 =51.6;
- if( m_Rlayer2 ==0.) m_Rlayer2 =90.0;
- m_Rlayer3 =122.5;
+ if( m_beampipeR==0.) m_beampipeR=24.0;
+ if( m_rLayerB ==0.) m_rLayerB =34.0;
+ if( m_rLayer1 ==0.) m_rLayer1 =51.6;
+ if( m_rLayer2 ==0.) m_rLayer2 =90.0;
+ m_rLayer3 =122.5;
} else { // 3-layer pixel detector
- if( m_Rbeampipe==0.) m_Rbeampipe=29.4;
- if( m_RlayerB ==0.) m_RlayerB =51.5;
- if( m_Rlayer1 ==0.) m_Rlayer1 =90.0;
- if( m_Rlayer2 ==0.) m_Rlayer2 =122.5;
+ if( m_beampipeR==0.) m_beampipeR=29.4;
+ if( m_rLayerB ==0.) m_rLayerB =51.5;
+ if( m_rLayer1 ==0.) m_rLayer1 =90.0;
+ if( m_rLayer2 ==0.) m_rLayer2 =122.5;
}
//
//
ITHistSvc* hist_root=0;
- if(m_FillHist){
+ if(m_fillHist){
StatusCode sc = service( "THistSvc", hist_root);
if( sc.isFailure() ) {
@@ -251,7 +233,8 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
}
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "InDetVKalVxInJetTool Histograms found" << endmsg;
- m_hb_massPiPi = new TH1D("massPiPi"," massPiPi",160,200., 1000.);
+ m_hb_massPiPi = new TH1D("massPiPi"," mass PiPi",160,200., 1000.);
+ m_hb_massPiPi1 = new TH1D("massPiPi1"," mass PiPi",100,200., 2000.);
m_hb_massPPi = new TH1D("massPPi"," massPPi", 100,1000., 1250.);
m_hb_massEE = new TH1D("massEE"," massEE", 100,0., 200.);
m_hb_nvrt2 = new TH1D("nvrt2"," vertices2", 50,0., 50.);
@@ -262,9 +245,10 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
m_hb_totmass2T1 = new TH1D("mass2trcomvrt1"," totmass 2tr common vertex", 200,0., 10000.);
m_hb_totmass2T2 = new TH1D("mass2trcomvrt2"," totmass 2tr common vertex", 200,0., 10000.);
m_hb_impact = new TH1D("impact", " impact", 100,0., 20.);
- m_hb_impactR = new TH1D("impactR"," impactR", 100,-30., 70.);
+ m_hb_impactR = new TH1D("impactR"," impactR", 400,-30., 70.);
m_hb_impactZ = new TH1D("impactZ"," impactZ", 100,-30., 70.);
m_hb_impactRZ = new TH2D("impactRZ"," impactRZ", 40,-10., 10., 60, -30.,30. );
+ m_hb_trkD0 = new TH1D("trkD0"," d0 of tracks", 100,-20., 20.);
m_hb_r2d = new TH1D("r2interact","Interaction radius 2tr selected", 150,0., 150.);
m_hb_r1dc = new TH1D("r1interactCommon","Interaction 1tr radius common", 150,0., 150.);
m_hb_r2dc = new TH1D("r2interactCommon","Interaction 2tr radius common", 150,0., 150.);
@@ -276,7 +260,6 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
m_hb_jmom = new TH1D("jetmom"," Jet mom", 200,0., 2000000.);
m_hb_mom = new TH1D("jetmomvrt"," Jet mom with sec. vertex", 200,0., 2000000.);
m_hb_signif3D = new TH1D("signif3D"," Signif3D for initial 2tr vertex", 140,-20., 50.);
- m_hb_massPiPi1 = new TH1D("massPiPi1"," massPiPi",100,0., 4000.);
m_hb_sig3DTot = new TH1D("sig3dcommon"," Signif3D for common vertex", 140,-20., 50.);
m_hb_sig3D1tr = new TH1D("sig3D1tr","Signif3D for 1tr vertices", 140,-20., 50.);
m_hb_sig3D2tr = new TH1D("sig3D2tr","Signif3D for 2tr single vertex", 140,-20., 50.);
@@ -284,30 +267,24 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
m_hb_goodvrtN = new TH1F("goodvrtN","Number of good vertices", 20,0., 20.);
m_hb_distVV = new TH1D("distvv","Vertex-Vertex dist", 100,0., 20.);
m_hb_diffPS = new TH1D("diffPS","Primary-Secondary assoc", 200,-20., 20.);
- m_hb_trkPtMax = new TH1D("trkPtMax","Maximal track Pt to jet", 100, 0., 10000.);
- m_hb_tr2SelVar = new TH1D("tr2SelVar","New 2tr variable", 200, 0., 10.);
+ m_hb_trkPtMax = new TH1D("trkPtMax","Maximal track Pt to jet", 100, 0., 5000.);
m_hb_blshared = new TH1F("blshared","Number of shared hits in B-layer for R<BL", 5, 0., 5.);
m_hb_pxshared = new TH1F("pxshared","Number of shared hits in pixel for R>BL", 5, 0., 5.);
- m_hb_rawVrtN = new TH1F("rawVrtN","Number of raw vertices multivertex case", 10, 0., 10.);
+ m_hb_rawVrtN = new TH1F("rawVrtN","Number of raw vertices multivertex case", 20, 0., 20.);
m_hb_lifetime = new TH1F("lifetime","Distance/momentum", 100, 0., 5.);
m_hb_trkPErr = new TH1F("trkPErr","Track momentum error for P>10 GeV", 100, 0., 0.5);
m_hb_deltaRSVPV = new TH1F("deltaRSVPV","SV-PV vs jet dR ", 200, 0., 1.);
//---
- m_hb_massJetTrkSV = new TH1D("PSEUmassJetTrkSV","SV mass for jet+track case", 250, 0., 10000.);
- m_hb_ratioJetTrkSV = new TH1D("PSEUratioJetTrkSV","SV ratio for jet+track case", 51,0., 1.02);
- m_hb_DST_JetTrkSV = new TH1D("PSEUDST_JetTrkSV", "DST PV-SV for jet+track case", 100,0., 20.);
- m_hb_NImpJetTrkSV = new TH1D("PSEUnTrkJetTrkSV", "N Track selected for jet+track case", 10,0., 10.);
-//---
- m_hb_nHImpTrkCnt = new TH1D("NHImpTrkCnt", "N Big Impact Track selected in jet", 30,0., 30.);
-//---
+ m_pr_NSelTrkMean = new TProfile("NSelTrkMean"," NTracks selected vs jet pt", 80, 0., 1600000.);
m_pr_effVrt2tr = new TProfile("effVrt2tr"," 2tr vertex efficiency vs Ntrack", 50, 0., 50.);
m_pr_effVrt2trEta= new TProfile("effVrt2trEta"," 2tr vertex efficiency vs eta", 50, -3., 3.);
m_pr_effVrt = new TProfile("effVrt","Full vertex efficiency vs Ntrack", 50, 0., 50.);
m_pr_effVrtEta= new TProfile("effVrtEta","Full vertex efficiency vs eta", 50, -3., 3.);
std::string histDir;
- if(m_MultiVertex) histDir="/file1/stat/MSVrtInJet"+m_instanceName+"/";
+ if(m_multiVertex) histDir="/file1/stat/MSVrtInJet"+m_instanceName+"/";
else histDir="/file1/stat/SVrtInJet"+m_instanceName+"/";
sc = hist_root->regHist(histDir+"massPiPi", m_hb_massPiPi);
+ sc = hist_root->regHist(histDir+"massPiPi1", m_hb_massPiPi1);
sc = hist_root->regHist(histDir+"massPPi", m_hb_massPPi);
sc = hist_root->regHist(histDir+"massEE", m_hb_massEE );
sc = hist_root->regHist(histDir+"nvrt2", m_hb_nvrt2);
@@ -321,6 +298,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
sc = hist_root->regHist(histDir+"impactR", m_hb_impactR);
sc = hist_root->regHist(histDir+"impactZ", m_hb_impactZ);
sc = hist_root->regHist(histDir+"impactRZ", m_hb_impactRZ);
+ sc = hist_root->regHist(histDir+"trkD0", m_hb_trkD0);
sc = hist_root->regHist(histDir+"r2interact", m_hb_r2d);
sc = hist_root->regHist(histDir+"r1interactCommon", m_hb_r1dc);
sc = hist_root->regHist(histDir+"r2interactCommon", m_hb_r2dc);
@@ -332,7 +310,6 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
sc = hist_root->regHist(histDir+"jetmom", m_hb_jmom);
sc = hist_root->regHist(histDir+"jetmomvrt", m_hb_mom);
sc = hist_root->regHist(histDir+"signif3D", m_hb_signif3D);
- sc = hist_root->regHist(histDir+"massPiPi1", m_hb_massPiPi1);
sc = hist_root->regHist(histDir+"sig3dcommon", m_hb_sig3DTot);
sc = hist_root->regHist(histDir+"sig3D1tr", m_hb_sig3D1tr);
sc = hist_root->regHist(histDir+"sig3D2tr", m_hb_sig3D2tr);
@@ -341,30 +318,63 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
sc = hist_root->regHist(histDir+"distVV", m_hb_distVV);
sc = hist_root->regHist(histDir+"diffPS", m_hb_diffPS);
sc = hist_root->regHist(histDir+"trkPtMax", m_hb_trkPtMax);
- sc = hist_root->regHist(histDir+"tr2SelVar", m_hb_tr2SelVar);
sc = hist_root->regHist(histDir+"blshared", m_hb_blshared);
sc = hist_root->regHist(histDir+"pxshared", m_hb_pxshared);
sc = hist_root->regHist(histDir+"rawVrtN", m_hb_rawVrtN);
sc = hist_root->regHist(histDir+"lifetime", m_hb_lifetime);
sc = hist_root->regHist(histDir+"trkPErr", m_hb_trkPErr);
sc = hist_root->regHist(histDir+"deltaRSVPV", m_hb_deltaRSVPV);
+ sc = hist_root->regHist(histDir+"NSelTrkMean", m_pr_NSelTrkMean);
sc = hist_root->regHist(histDir+"effVrt2tr", m_pr_effVrt2tr);
sc = hist_root->regHist(histDir+"effVrt2trEta", m_pr_effVrt2trEta);
sc = hist_root->regHist(histDir+"effVrt", m_pr_effVrt);
sc = hist_root->regHist(histDir+"effVrtEta", m_pr_effVrtEta);
- sc = hist_root->regHist(histDir+"PSEUmassJetTrkSV", m_hb_massJetTrkSV);
- sc = hist_root->regHist(histDir+"PSEUratioJetTrkSV",m_hb_ratioJetTrkSV);
- sc = hist_root->regHist(histDir+"PSEUDST_JetTrkSV", m_hb_DST_JetTrkSV);
- sc = hist_root->regHist(histDir+"PSEUnTrkJetTrkSV", m_hb_NImpJetTrkSV);
- sc = hist_root->regHist(histDir+"NHImpTrkCnt", m_hb_nHImpTrkCnt);
if( sc.isFailure() ) { // Check of StatusCode
if(msgLvl(MSG::INFO))msg(MSG::INFO) << "BTagVrtSec Histogram registration failure!!!" << endmsg;
}
m_w_1 = 1.;
-
+//-------------------------------------------------------
+ m_curTup=new DevTuple();
+ m_tuple = new TTree("Tracks","Tracks");
+ std::string TreeDir("/file1/stat/SVrtInJet"+m_instanceName+"/");
+ sc = hist_root->regTree(TreeDir,m_tuple);
+ if (sc.isSuccess()) {
+ m_tuple->Branch("ptjet", &m_curTup->ptjet, "ptjet/F");
+ m_tuple->Branch("etajet", &m_curTup->etajet, "etajet/F");
+ m_tuple->Branch("phijet", &m_curTup->phijet, "phijet/F");
+ m_tuple->Branch("ntrk", &m_curTup->nTrkInJet, "ntrk/I");
+ m_tuple->Branch("prbS", &m_curTup->s_prob, "prbS[ntrk]/F");
+ m_tuple->Branch("prbP", &m_curTup->p_prob, "prbP[ntrk]/F");
+ m_tuple->Branch("wgtB", &m_curTup->wgtB, "wgtB[ntrk]/F");
+ m_tuple->Branch("wgtL", &m_curTup->wgtL, "wgtL[ntrk]/F");
+ m_tuple->Branch("wgtG", &m_curTup->wgtG, "wgtG[ntrk]/F");
+ m_tuple->Branch("Sig3D", &m_curTup->Sig3D, "Sig3D[ntrk]/F");
+ m_tuple->Branch("idMC", &m_curTup->idMC, "idMC[ntrk]/I");
+ m_tuple->Branch("ibl", &m_curTup->ibl, "ibl[ntrk]/I");
+ m_tuple->Branch("bl", &m_curTup->bl, "bl[ntrk]/I");
+ m_tuple->Branch("SigR", &m_curTup->SigR, "SigR[ntrk]/F");
+ m_tuple->Branch("SigZ", &m_curTup->SigZ, "SigZ[ntrk]/F");
+ m_tuple->Branch("d0", &m_curTup->d0, "d0[ntrk]/F");
+ m_tuple->Branch("Z0", &m_curTup->Z0, "Z0[ntrk]/F");
+ m_tuple->Branch("pTvsJet", &m_curTup->pTvsJet, "pTvsJet[ntrk]/F");
+ m_tuple->Branch("prodTJ", &m_curTup->prodTJ, "prodTJ[ntrk]/F");
+ m_tuple->Branch("nVrtT", &m_curTup->nVrtT, "nVrtT[ntrk]/I");
+ m_tuple->Branch("chg", &m_curTup->chg, "chg[ntrk]/I");
+ m_tuple->Branch("nvrt", &m_curTup->nVrt, "nvrt/I");
+ m_tuple->Branch("VrtDist2D", &m_curTup->VrtDist2D, "VrtDist2D[nvrt]/F");
+ m_tuple->Branch("VrtSig3D", &m_curTup->VrtSig3D, "VrtSig3D[nvrt]/F");
+ m_tuple->Branch("VrtSig2D", &m_curTup->VrtSig2D, "VrtSig2D[nvrt]/F");
+ m_tuple->Branch("itrk", &m_curTup->itrk, "itrk[nvrt]/I");
+ m_tuple->Branch("jtrk", &m_curTup->jtrk, "jtrk[nvrt]/I");
+ m_tuple->Branch("badV", &m_curTup->badVrt, "badV[nvrt]/I");
+ m_tuple->Branch("mass", &m_curTup->mass, "mass[nvrt]/F");
+ m_tuple->Branch("Chi2", &m_curTup->Chi2, "Chi2[nvrt]/F");
+ m_tuple->Branch("ntHF", &m_curTup->NTHF, "ntHF/I");
+ m_tuple->Branch("itHF", &m_curTup->itHF, "itHF[ntHF]/I");
+ }
}
- if(!m_MultiVertex)m_MultiWithPrimary = false;
+ if(!m_multiVertex)m_multiWithPrimary = false;
if(m_getNegativeTag){
if(msgLvl(MSG::INFO))msg(MSG::INFO) << " Negative TAG is requested! " << endmsg;
@@ -412,6 +422,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
double EnergyJet = 0.;
int N2trVertices = 0 ;
int NBigImpTrk = 0 ;
+ if(m_curTup){ m_curTup->nVrt=0; m_curTup->nTrkInJet=0; m_curTup->NTHF=0; }
int pseudoVrt = 0;
@@ -422,7 +433,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
if(tmp)InpTrk.push_back(tmp);
}
- if(m_MultiVertex){
+ if(m_multiVertex){
workVectorArrxAOD * tmpVectxAOD=new workVectorArrxAOD();
tmpVectxAOD->InpTrk.resize(InpTrk.size());
std::copy(InpTrk.begin(),InpTrk.end(), tmpVectxAOD->InpTrk.begin());
@@ -434,7 +445,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
}else{
xAOD::Vertex* secVrt = GetVrtSec( InpTrk,PrimVrt,JetDir,Results,SelSecTrk,xaodTrkFromV0);
if(secVrt != 0) listVrtSec.push_back(secVrt);
- else if(m_FillHist){ m_pr_effVrt->Fill((float)m_NRefTrk,0.);
+ else if(m_fillHist){ m_pr_effVrt->Fill((float)m_NRefPVTrk,0.);
m_pr_effVrtEta->Fill( JetDir.Eta(),0.);}
}
if(Results.size()<3) {
@@ -460,6 +471,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
}
+ if(m_fillHist){ m_tuple->Fill(); };
m_fitSvc->clearMemory();
m_compatibilityGraph->clear();
std::vector<int> zytmp(1000); m_WorkArray->m_Incomp.swap(zytmp); // Deallocate memory
@@ -483,12 +495,13 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
double RatioE = 0.;
double EnergyJet = 0.;
int N2trVertices = 0 ;
+ if(m_curTup){ m_curTup->nVrt=0; m_curTup->nTrkInJet=0; }
xAOD::Vertex xaodPrimVrt;
xaodPrimVrt.setPosition(PrimVrt.position());
xaodPrimVrt.setCovariancePosition(PrimVrt.covariancePosition());
- if(m_MultiVertex){
+ if(m_multiVertex){
workVectorArrREC * tmpVectREC=new workVectorArrREC();
tmpVectREC->InpTrk.resize(InpTrk.size());
std::copy(InpTrk.begin(),InpTrk.end(), tmpVectREC->InpTrk.begin());
@@ -500,7 +513,7 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
}else{
xAOD::Vertex* secVrt = GetVrtSec( InpTrk,xaodPrimVrt,JetDir,Results,SelSecTrk,TrkFromV0);
if(secVrt != 0) listVrtSec.push_back(secVrt);
- else if(m_FillHist){ m_pr_effVrt->Fill((float)m_NRefTrk,0.);
+ else if(m_fillHist){ m_pr_effVrt->Fill((float)m_NRefPVTrk,0.);
m_pr_effVrtEta->Fill( JetDir.Eta(),0.);}
}
if(Results.size()<3) {
@@ -514,6 +527,8 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
const Trk::VxSecVKalVertexInfo* res =
new Trk::VxSecVKalVertexInfo(listVrtSec, SecVtxMass, RatioE, N2trVertices, EnergyJet, PartToBase(TrkFromV0) );
if(Results.size()>8)res->setDstToMatLay(Results[7]);
+
+ if(m_fillHist){ m_tuple->Fill(); };
m_fitSvc->clearMemory();
m_compatibilityGraph->clear();
std::vector<int> zytmp(1000); m_WorkArray->m_Incomp.swap(zytmp); // Deallocate memory
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
index ddd1cff33020f0d00f67c9258918b21f2e902fab..b1fd5d2caa6a4d85dae19f01904fce3d61d47ff8 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
// Header include
@@ -7,14 +7,30 @@
#include "TrkNeutralParameters/NeutralParameters.h"
#include "TrkTrackSummary/TrackSummary.h"
#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
-
+#include "xAODTruth/TruthParticleContainer.h"
//-------------------------------------------------
// Other stuff
#include <cmath>
+//#include<iostream>
namespace InDet{
+ double InDetVKalVxInJetTool::RankBTrk(double TrkPt, double JetPt, double Signif) const
+ {
+ double coeffSig=1.0;
+ double s_prob=(Signif-coeffSig)/Signif; // Old probability to be b-track
+ double coeffPt=10.;
+ double pfrac=(TrkPt-m_cutPt)/sqrt(JetPt);
+ double p_prob= pfrac/(coeffPt+pfrac); // Old probability to be b-track
+ if(TrkPt + JetPt == 0.) return s_prob;
+ else if(Signif == 0.) return p_prob;
+ //----------------------------------Initial definition of selective variable
+ double contrib=0.4;
+ return (1.+contrib)*std::max(s_prob,0.)+(1.-contrib)*p_prob;
+ }
+
+
TLorentzVector InDetVKalVxInJetTool::GetBDir( const xAOD::TrackParticle* trk1,
const xAOD::TrackParticle* trk2,
const xAOD::Vertex & PrimVrt,
@@ -22,7 +38,6 @@ namespace InDet{
const
{ // B hadron flight direction based on 2 separate tracks and PV. Calculated via plane-plane crossing
Amg::Vector3D PVRT(PrimVrt.x(),PrimVrt.y(),PrimVrt.z());
-
//----------------------------------------------------------------------------
Amg::Vector3D pnt1=trk1->perigeeParameters().position()-PVRT;
Amg::Vector3D mom1((trk1->p4()).Px(),(trk1->p4()).Py(),(trk1->p4()).Pz());
@@ -55,40 +70,54 @@ namespace InDet{
}
void InDetVKalVxInJetTool::printWrkSet(const std::vector<WrkVrt> *, const std::string ) const {
-
+/* void InDetVKalVxInJetTool::printWrkSet(const std::vector<WrkVrt> *WrkVrtSet, const std::string name) const {
+ int nGoodV=0;
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
+ std::cout<<name
+ <<"= "<<(*WrkVrtSet)[iv].vertex[0]
+ <<", "<<(*WrkVrtSet)[iv].vertex[1]
+ <<", "<<(*WrkVrtSet)[iv].vertex[2]
+ <<" NTrk="<<(*WrkVrtSet)[iv].SelTrk.size()
+ <<" is good="<<std::boolalpha<<(*WrkVrtSet)[iv].Good<<std::noboolalpha
+ <<" Chi2="<<(*WrkVrtSet)[iv].Chi2
+ <<" Mass="<<(*WrkVrtSet)[iv].vertexMom.M()
+ <<" detached="<<(*WrkVrtSet)[iv].detachedTrack
+ <<" proj.dist="<<(*WrkVrtSet)[iv].ProjectedVrt
+ <<" trk=";
+ for(int kk=0; kk<(int)(*WrkVrtSet)[iv].SelTrk.size(); kk++) {
+ std::cout<<", "<<(*WrkVrtSet)[iv].SelTrk[kk];}
+ //for(int kk=0; kk<(int)(*WrkVrtSet)[iv].SelTrk.size(); kk++) {
+ // std::cout<<", "<<MomAtVrt((*WrkVrtSet)[iv].TrkAtVrt[kk]).Perp();}
+ std::cout<<'\n';
+ if((*WrkVrtSet)[iv].Good)nGoodV++;
+ }
+ std::cout<<name<<" N="<<nGoodV<<'\n';*/
}
/* Technicalities */
- double InDetVKalVxInJetTool::ProjPos(const Amg::Vector3D & Vrt, const TLorentzVector & JetDir)
- const
- {
- //Amg::Vector3D Vrt=SV-PV;
- return (Vrt.x()*JetDir.Px() + Vrt.y()*JetDir.Py() + Vrt.z()*JetDir.Pz())/JetDir.P();
- }
-
- double InDetVKalVxInJetTool::ProjPosT(const Amg::Vector3D & Vrt, const TLorentzVector & JetDir)
- const
- {
- return (Vrt.x()*JetDir.Px() + Vrt.y()*JetDir.Py())/JetDir.Pt();
+ double InDetVKalVxInJetTool::ProjSV_PV(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Jet) const
+ {
+ TVector3 SV_PV( SV.x()-PV.x(), SV.y()-PV.y(), SV.z()-PV.z() );
+ return Jet.Vect().Unit()*SV_PV.Unit();
}
bool InDetVKalVxInJetTool::insideMatLayer(float xvt,float yvt) const
{
- float Dist2DBP=sqrt( (xvt-m_Xbeampipe)*(xvt-m_Xbeampipe) + (yvt-m_Ybeampipe)*(yvt-m_Ybeampipe) );
- float Dist2DBL=sqrt( (xvt-m_XlayerB)*(xvt-m_XlayerB) + (yvt-m_YlayerB)*(yvt-m_YlayerB) );
- float Dist2DL1=sqrt( (xvt-m_Xlayer1)*(xvt-m_Xlayer1) + (yvt-m_Ylayer1)*(yvt-m_Ylayer1) );
- float Dist2DL2=sqrt( (xvt-m_Xlayer2)*(xvt-m_Xlayer2) + (yvt-m_Ylayer2)*(yvt-m_Ylayer2) );
+ float Dist2DBP=sqrt( (xvt-m_beampipeX)*(xvt-m_beampipeX) + (yvt-m_beampipeY)*(yvt-m_beampipeY) );
+ float Dist2DBL=sqrt( (xvt-m_xLayerB)*(xvt-m_xLayerB) + (yvt-m_yLayerB)*(yvt-m_yLayerB) );
+ float Dist2DL1=sqrt( (xvt-m_xLayer1)*(xvt-m_xLayer1) + (yvt-m_yLayer1)*(yvt-m_yLayer1) );
+ float Dist2DL2=sqrt( (xvt-m_xLayer2)*(xvt-m_xLayer2) + (yvt-m_yLayer2)*(yvt-m_yLayer2) );
if(m_existIBL){ // 4-layer pixel detector
- if( fabs(Dist2DBP-m_Rbeampipe)< 1.0) return true; // Beam Pipe removal
- if( fabs(Dist2DBL-m_RlayerB) < 2.5) return true;
- if( fabs(Dist2DL1-m_Rlayer1) < 3.0) return true;
- if( fabs(Dist2DL2-m_Rlayer2) < 3.0) return true;
- //if( fabs(Dist2DL2-m_Rlayer3) < 4.0) return true;
+ if( fabs(Dist2DBP-m_beampipeR)< 1.0) return true; // Beam Pipe removal
+ if( fabs(Dist2DBL-m_rLayerB) < 2.5) return true;
+ if( fabs(Dist2DL1-m_rLayer1) < 3.0) return true;
+ if( fabs(Dist2DL2-m_rLayer2) < 3.0) return true;
+ //if( fabs(Dist2DL2-m_rLayer3) < 4.0) return true;
}else{ // 3-layer pixel detector
- if( fabs(Dist2DBP-m_Rbeampipe)< 1.5) return true; // Beam Pipe removal
- if( fabs(Dist2DBL-m_RlayerB) < 3.5) return true;
- if( fabs(Dist2DL1-m_Rlayer1) < 4.0) return true;
- if( fabs(Dist2DL2-m_Rlayer2) < 5.0) return true;
+ if( fabs(Dist2DBP-m_beampipeR)< 1.5) return true; // Beam Pipe removal
+ if( fabs(Dist2DBL-m_rLayerB) < 3.5) return true;
+ if( fabs(Dist2DL1-m_rLayer1) < 4.0) return true;
+ if( fabs(Dist2DL2-m_rLayer2) < 5.0) return true;
}
return false;
}
@@ -158,6 +187,31 @@ namespace InDet{
return sqrt(distx*distx+disty*disty+distz*distz);
}
+ double InDetVKalVxInJetTool::VrtVrtDist2D(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> SecVrtErr, double& Signif)
+ const
+ {
+ double distx = PrimVrt.x()- SecVrt.x();
+ double disty = PrimVrt.y()- SecVrt.y();
+
+
+ AmgSymMatrix(3) PrimCovMtx=PrimVrt.covariancePosition(); //Create
+ AmgSymMatrix(2) CovMtx;
+ CovMtx(0,0) = PrimCovMtx(0,0) + SecVrtErr[0];
+ CovMtx(0,1) = PrimCovMtx(0,1) + SecVrtErr[1];
+ CovMtx(1,0) = PrimCovMtx(1,0) + SecVrtErr[1];
+ CovMtx(1,1) = PrimCovMtx(1,1) + SecVrtErr[2];
+
+ AmgSymMatrix(2) WgtMtx = CovMtx.inverse();
+
+ Signif = distx*WgtMtx(0,0)*distx
+ +disty*WgtMtx(1,1)*disty
+ +2.*distx*WgtMtx(0,1)*disty;
+ Signif=sqrt(Signif);
+ if( Signif!=Signif ) Signif = 0.;
+ return sqrt(distx*distx+disty*disty);
+ }
+
//--------------------------------------------------
// Significance along jet direction
//--------------------------------------------------
@@ -267,6 +321,22 @@ namespace InDet{
//
+//----------------------------
+// Vertex error along radius
+//----------------------------
+ double InDetVKalVxInJetTool::VrtRadiusError(const Amg::Vector3D & SecVrt, const std::vector<double> & VrtErr) const
+ {
+ double DirX=SecVrt.x(), DirY=SecVrt.y();
+ double Covar = DirX*VrtErr[0]*DirX
+ +2.*DirX*VrtErr[1]*DirY
+ +DirY*VrtErr[2]*DirY;
+ Covar /= DirX*DirX + DirY*DirY;
+ Covar=sqrt(Covar);
+ if(Covar != Covar) Covar = 0.;
+ return Covar;
+ }
+
+
double InDetVKalVxInJetTool::ConeDist(const AmgVector(5) & VectPerig, const TLorentzVector & JetDir)
const
@@ -304,16 +374,15 @@ namespace InDet{
- int InDetVKalVxInJetTool::FindMax( std::vector<double>& Chi2PerTrk, std::vector<int> & cntTrk)
+ int InDetVKalVxInJetTool::FindMax( std::vector<double>& Chi2PerTrk, std::vector<float> & Rank)
const
{
double Chi2Ref=0.;
- int Position=0;
+ int Position=-1;
if( Chi2PerTrk.size() < 1 ) return Position ;
for (int i=0; i< (int)Chi2PerTrk.size(); i++){
- if( Chi2PerTrk[i]/cntTrk[i] > Chi2Ref) { Chi2Ref=Chi2PerTrk[i]/cntTrk[i]; Position=i;}
+ if(Chi2PerTrk[i]/std::max(Rank[i],(float)0.1) > Chi2Ref) { Chi2Ref=Chi2PerTrk[i]/std::max(Rank[i],(float)0.1); Position=i;}
}
-
return Position;
}
@@ -519,8 +588,10 @@ namespace InDet{
}
-
-
+ StatusCode InDetVKalVxInJetTool::GetTrkFitWeights(std::vector<double> & wgt) const
+ {
+ return m_fitSvc->VKalGetTrkWeights(wgt);
+ }
/*************************************************************************************************************/
void InDetVKalVxInJetTool::getPixelLayers(const Rec::TrackParticle* Part, int &blHit, int &l1Hit, int &l2Hit, int &nLays ) const
{
@@ -558,7 +629,7 @@ namespace InDet{
// bitH=HitPattern&((int)pow(2,Trk::pixelBarrel1));
} else { // 3-layer pixel detector
uint8_t BLhit,NPlay,NHoles,IBLhit;
- if(!Part->summaryValue( BLhit, xAOD::numberOfInnermostPixelLayerHits) ) BLhit = 0;
+ if(!Part->summaryValue( BLhit, xAOD::numberOfBLayerHits) ) BLhit = 0;
if(!Part->summaryValue(IBLhit, xAOD::numberOfInnermostPixelLayerHits) ) IBLhit = 0; // Some safety
BLhit=BLhit>IBLhit ? BLhit : IBLhit; // Some safety
if(!Part->summaryValue( NPlay, xAOD::numberOfContribPixelLayers) ) NPlay = 0;
@@ -666,10 +737,41 @@ namespace InDet{
return VrtCovMtx;
}
- double InDetVKalVxInJetTool::trkPtCorr(double pT) const
- {
- double adp=pT/64000.; if(adp<0.)adp=0; if(adp>1.)adp=1.; adp=sqrt(adp)/2.;
- return adp;
+ void InDetVKalVxInJetTool::fillVrtNTup( std::vector<Vrt2Tr> & all2TrVrt)
+ const
+ { if(!m_curTup)return;
+ int ipnt=0;
+ for(auto vrt : all2TrVrt) {
+ if(ipnt==100)break;
+ m_curTup->VrtDist2D[ipnt]=vrt.FitVertex.perp();
+ m_curTup->VrtSig3D[ipnt]=vrt.Signif3D;
+ m_curTup->VrtSig2D[ipnt]=vrt.Signif2D;
+ m_curTup->itrk[ipnt]=vrt.i;
+ m_curTup->jtrk[ipnt]=vrt.j;
+ m_curTup->mass[ipnt]=vrt.Momentum.M();
+ m_curTup->Chi2[ipnt]=vrt.Chi2;
+ m_curTup->badVrt[ipnt]=vrt.badVrt;
+ ipnt++; m_curTup->nVrt=ipnt;
+ }
+ }
+
+ int InDetVKalVxInJetTool::getG4Inter(const xAOD::TrackParticle* TP ) const {
+ if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
+ const ElementLink<xAOD::TruthParticleContainer>& tplink =
+ TP->auxdata< ElementLink< xAOD::TruthParticleContainer > >("truthParticleLink");
+ if( tplink.isValid() && (*tplink)->barcode()>200000) return 1;
+ }
+ return 0;
+ }
+ int InDetVKalVxInJetTool::getMCPileup(const xAOD::TrackParticle* TP ) const {
+ if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
+ const ElementLink<xAOD::TruthParticleContainer>& tplink =
+ TP->auxdata< ElementLink< xAOD::TruthParticleContainer > >("truthParticleLink");
+ if( !tplink.isValid() ) return 1;
+ } else { return 1; }
+ return 0;
}
+ int InDetVKalVxInJetTool::getG4Inter(const Rec::TrackParticle* ) const { return 0; }
+ int InDetVKalVxInJetTool::getMCPileup(const Rec::TrackParticle* ) const { return 0; }
} //end namespace
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/components/InDetVKalVxInJetTool_entries.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/components/InDetVKalVxInJetTool_entries.cxx
index 24809f0610bb5ed8bdb3564e872b67f2681936e6..6262d64189554b94732996c72e8528a14ac90e6b 100644
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/components/InDetVKalVxInJetTool_entries.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/components/InDetVKalVxInJetTool_entries.cxx
@@ -1,6 +1,8 @@
#include "InDetVKalVxInJetTool/InDetVKalVxInJetTool.h"
+#include "InDetVKalVxInJetTool/InDetTrkInJetType.h"
using namespace InDet;
DECLARE_COMPONENT( InDetVKalVxInJetTool )
+DECLARE_COMPONENT( InDetTrkInJetType )