diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h
index 85cf7b30328ee21fa27d89e69491786703fbc302..3f947855d42c3b2c49b761dcbb504b8026c9aca7 100644
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetTrkInJetType.h
@@ -79,11 +79,14 @@ namespace InDet {
private:
TMVA::Reader* m_tmvaReader;
- double m_trkMinPtCut;
- float m_d0_limLow;
- float m_d0_limUpp;
- float m_Z0_limLow;
- float m_Z0_limUpp;
+ int m_trkSctHitsCut{};
+ int m_trkPixelHitsCut{};
+ float m_trkChi2Cut{};
+ float 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{};
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
index cfb8d0c68ede6af5645d378815a310a5d245551b..e6da7b7d5b417d5653faf99813c4d0a46da699f4 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/InDetVKalVxInJetTool/InDetVKalVxInJetTool.h
@@ -48,6 +48,7 @@
//
#include "Particle/TrackParticle.h"
#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTruth/TruthEventContainer.h"
//
//
#include "InDetRecToolInterfaces/ISecVertexInJetFinder.h"
@@ -120,17 +121,6 @@ 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
//
@@ -212,6 +202,8 @@ namespace InDet {
double m_zTrkErrorCut{};
double m_cutHFClass{};
double m_antiGarbageCut{};
+ double m_antiFragmentCut{};
+ double m_Vrt2TrMassLimit{};
bool m_fillHist{};
@@ -261,6 +253,10 @@ namespace InDet {
//-------------------------------------------
//For ntuples (only for development/tuning!)
+ int notFromBC(int PDGID) const;
+ const xAOD::TruthParticle * getPreviousParent(const xAOD::TruthParticle * child, int & ParentPDG) const;
+ int getIdHF(const Rec::TrackParticle* TP ) const;
+ int getIdHF(const xAOD::TrackParticle* TP ) const;
int getG4Inter( const Rec::TrackParticle* TP ) const;
int getG4Inter( const xAOD::TrackParticle* TP ) const;
int getMCPileup(const Rec::TrackParticle* TP ) const;
@@ -269,6 +265,7 @@ namespace InDet {
struct DevTuple
{
static const int maxNTrk=100;
+ static const int maxNVrt=100;
int nTrkInJet;
float ptjet;
float etajet;
@@ -288,19 +285,31 @@ namespace InDet {
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];
+ float TotM;
+ int nVrt;
+ float VrtDist2D[maxNVrt];
+ float VrtSig3D[maxNVrt];
+ float VrtSig2D[maxNVrt];
+ float VrtDR[maxNVrt];
+ float mass[maxNVrt];
+ float Chi2[maxNVrt];
+ int itrk[maxNVrt];
+ int jtrk[maxNVrt];
+ int badVrt[maxNVrt];
+ int ibl[maxNVrt];
+ int bl[maxNVrt];
int NTHF;
- int itHF[maxNTrk];
+ int itHF[maxNVrt];
+ //---
+ int nNVrt;
+ float NVrtDist2D[maxNVrt];
+ int NVrtNT[maxNVrt];
+ int NVrtTrkI[maxNVrt];
+ float NVrtM[maxNVrt];
+ float NVrtChi2[maxNVrt];
+ float NVrtMaxW[maxNVrt];
+ float NVrtAveW[maxNVrt];
+ float NVrtDR[maxNVrt];
};
DevTuple* m_curTup;
@@ -317,6 +326,7 @@ namespace InDet {
double Signif3DProj=0.;
double Signif2D=0.;
double Chi2=0.;
+ double dRSVPV=-1.;
};
@@ -326,7 +336,7 @@ namespace InDet {
float m_chiScale[11]{};
VKalVxInJetTemp* m_WorkArray{};
struct WrkVrt
- { bool Good{};
+ { bool Good=true;
std::deque<long int> SelTrk;
Amg::Vector3D vertex;
TLorentzVector vertexMom;
@@ -335,9 +345,9 @@ namespace InDet {
std::vector<double> Chi2PerTrk;
std::vector< std::vector<double> > TrkAtVrt;
double Chi2{};
- int nCloseVrt{};
- double dCloseVrt{};
- double ProjectedVrt{};
+ int nCloseVrt=0;
+ double dCloseVrt=1000000.;
+ double ProjectedVrt=0.;
int detachedTrack=-1;
};
@@ -411,7 +421,9 @@ namespace InDet {
double VrtRadiusError(const Amg::Vector3D & SecVrt, const std::vector<double> & VrtErr) const;
bool insideMatLayer(float ,float ) const;
- void fillVrtNTup( std::vector<Vrt2Tr> & all2TrVrt) const;
+ void fillVrtNTup( std::vector<Vrt2Tr> & all2TrVrt) const;
+ void fillNVrtNTup(std::vector<WrkVrt> & VrtSet, std::vector< std::vector<float> > & trkScore,
+ const xAOD::Vertex & PrimVrt, const TLorentzVector & JetDir)const;
TLorentzVector GetBDir( const xAOD::TrackParticle* trk1,
const xAOD::TrackParticle* trk2,
@@ -482,6 +494,9 @@ namespace InDet {
template <class Particle>
double mergeAndRefitVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2, WrkVrt & newvrt,
std::vector<const Particle*> & AllTrackList) const;
+ template <class Particle>
+ void mergeAndRefitOverlapVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2,
+ std::vector<const Particle*> & AllTrackList) const;
template <class Particle>
double improveVertexChi2( std::vector<WrkVrt> *WrkVrtSet, int V, std::vector<const Particle*> & AllTracks)const;
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/README b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/README
new file mode 100644
index 0000000000000000000000000000000000000000..73266aa7a361a4a6c8b15ee9f63027a099356eaf
--- /dev/null
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/Root/README
@@ -0,0 +1,11 @@
+Tool calibration data can be obtained by running the InDetVKalVxInJetTool with the following jobO:
+
+FillHist = 1
+AntiGarbageCut = 1.0
+
+Produce calibation files using different MC samples (ttbar, Zp, dijet,...)
+Edit Train.C to provide correct pathes to all this calibration files
+Run corrected Train.C using root 6.12 or above.
+
+Use HF tracks from b-jets
+Use fragmentation tracks from L+HF jets
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
index ecf1ccae1066ee6bad93e1b4659568b50b1b87ad..104fb338a7223343bd1d62803a82f80142df04a5 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSec.cxx
@@ -36,10 +36,20 @@ namespace InDet{
float median(std::vector<float> &Vec){
int N=Vec.size();
if(N==1) return Vec[0];
- if(N>1){
+ if(N==2) return (Vec[0]+Vec[1])/2.;
+ if(N==3) return Vec[1];
+ if(N>3){
std::vector<float> tmp(Vec);
- std::sort(tmp.begin(),tmp.end()); //can use nth_element instead of completely sorting, it's quicker
- return (tmp[(N-1)/2]+tmp[N/2])/2.; //only true if the number of elements is even?
+ int N05m=(N-1)/2, N05p=N/2;
+ //std::sort(tmp.begin(),tmp.end()); //can use nth_element instead of completely sorting, it's quicker
+ if(N05m==N05p){
+ std::nth_element(tmp.begin(),tmp.begin()+N05m,tmp.end());
+ return tmp[N05m];
+ } else {
+ std::nth_element(tmp.begin(),tmp.begin()+N05m,tmp.end());
+ std::nth_element(tmp.begin(),tmp.begin()+N05p,tmp.end());
+ return (tmp[N05m]+tmp[N05p])/2.;
+ }
}
return 0.;
}
@@ -508,6 +518,7 @@ namespace InDet{
float R=JetDir.DeltaR(TLorentzVector(FitVertex.x()-PrimVrt.x(),FitVertex.y()-PrimVrt.y(),
FitVertex.z()-PrimVrt.z(), 1.e4));
m_hb_deltaRSVPV->Fill( R, m_w_1);
+ if(m_curTup)m_curTup->TotM=Momentum.M();
}
//-------------------------------------------------------------------------------------
@@ -713,7 +724,6 @@ namespace InDet{
//
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);
@@ -731,7 +741,6 @@ namespace InDet{
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);
@@ -743,7 +752,7 @@ namespace InDet{
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(getIdHF(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];
@@ -759,10 +768,8 @@ namespace InDet{
}
}
}
- if (not m_curTup) return; //something very wrong
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
ListSecondTracks.reserve(2*NTracks); // Reserve memory for single vertex
@@ -776,8 +783,13 @@ namespace InDet{
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
+ if( std::max(trkScore[i][1],trkScore[j][1]) > m_antiFragmentCut ) continue; // Remove definite fragmentation tracks
+ bool goodPair=false;
+ float minWgtB=std::min(trkScore[i][0],trkScore[j][0]);
+ float maxWgtB=std::max(trkScore[i][0],trkScore[j][0]);
+ if( minWgtB>m_cutHFClass ) goodPair=true;
+ if( maxWgtB>0.5 && minWgtB>m_cutHFClass/2. )goodPair=true;
+ if( !goodPair ) continue;
}
int BadTracks = 0; //Bad tracks identification
TracksForFit.resize(2);
@@ -805,12 +817,15 @@ namespace InDet{
Dist2D=tmpVrt.FitVertex.perp();
if(Dist2D > 180. ) continue; // can't be from B decay
double mass_PiPi = tmpVrt.Momentum.M();
- if(mass_PiPi > c_vrtBCMassLimit) continue; // can't be from B decay
+ if(mass_PiPi > m_Vrt2TrMassLimit) 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);
//---
TVector3 SVmPV(tmpVrt.FitVertex.x()-PrimVrt.x(),tmpVrt.FitVertex.y()-PrimVrt.y(),tmpVrt.FitVertex.z()-PrimVrt.z());
+ tmpVrt.dRSVPV=JetDir.DeltaR(TLorentzVector(SVmPV, 1.)); //DeltaR SV-PV vs jet
+ if(tmpVrt.dRSVPV > m_coneForTag ) continue; // SV is outside of the jet cone
+//---
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*/
@@ -972,18 +987,27 @@ namespace InDet{
// { if( trkTypeSV[all2TrVrt[iv].i]>0 && trkTypeSV[all2TrVrt[iv].j]>0) all2TrVrt.erase(all2TrVrt.begin()+iv); else iv++; }
//
*/
+//-- Cleaning. Remove small wgtB tracks attached to one vertex only.
+// std::vector<int> inVrt(NTracks,0);
+// for( auto vv : all2TrVrt){ inVrt[vv.i]++; inVrt[vv.j]++; }
+// ////std::map<float,int> trkInOneV; for( int kt=0; kt<NTracks; kt++) if(inVrt[kt]==1) trkInOneV[trkScore[kt][0]]=kt;
+// for(int kk=0; kk<NTracks; kk++){
+// if( inVrt[kk]==1 && trkScore[kk][0]<2.5*m_cutHFClass ) {
+// int iv=0; while (iv<(int)all2TrVrt.size()){ if( all2TrVrt[iv].i==kk || all2TrVrt[iv].j==kk ) { all2TrVrt.erase(all2TrVrt.begin()+iv); break;} else iv++; }
+// } }
+//============================================================================
//-- Save results
ListSecondTracks.clear();
std::map<float,int> trkHF;
- for( auto vv : all2TrVrt){
+ 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]++; }
+ 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);
}
//
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
index 0a73eb71e87e7eb1a6e7b96e77c5e3ea70c2ea91..368e41ff99e85bf41f806575c921a8315be57905 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/BTagVrtSecMulti.cxx
@@ -56,8 +56,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
const
{
- const double probVrtMergeLimit=0.04;
- const int useMaterialRejection=1;
+ const double probVrtMergeLimit=0.01;
m_NRefPVTrk=0;
int inpNPart=0;
@@ -74,7 +73,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
if(msgLvl(MSG::DEBUG))msg(MSG::DEBUG) << "InDet GetVrtSecMulti() called with NPart=" <<inpNPart<< endmsg;
//std::vector<const Rec::TrackParticle*> listJetTracks, tmpListTracks, listSecondTracks, TracksForFit;
- std::vector<xAOD::Vertex*> finalVertices;
+ std::vector<xAOD::Vertex*> finalVertices(0);
if( inpNPart < 2 ) { return finalVertices;} // 0,1 track => nothing to do!
@@ -205,7 +204,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
newvrt.Good = true;
newvrt.nCloseVrt = 0;
newvrt.dCloseVrt = 1000000.;
- newvrt.ProjectedVrt=JetProjDist(newvrt.vertex, PrimVrt, JetDir);
+ newvrt.ProjectedVrt=JetProjDist(newvrt.vertex, PrimVrt, JetDir); //3D SV-PV distance
WrkVrtSet->push_back(newvrt);
}
//==================================================================================
@@ -220,45 +219,23 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
// 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));}}
// RemoveDoubleEntries(tempTrk);
-//--- Initial cleaning of solutions
//
-//- Try to merge vertices with 3 and more tracks which have at least 2 common tracks
- int NSoluI=(*WrkVrtSet).size();
- for(int iv=0; iv<NSoluI-1; iv++ ){
- if(!(*WrkVrtSet)[iv].Good ) continue;
- int nTrk1=(*WrkVrtSet)[iv].SelTrk.size();
- for(int jv=iv+1; jv<NSoluI; jv++){
- if(!(*WrkVrtSet)[jv].Good ) continue;
- int nTrk2=(*WrkVrtSet)[jv].SelTrk.size();
- if( nTrk1<3 || nTrk2<3 || nTrk1+nTrk2<7) continue;
- int nTCom=nTrkCommon( WrkVrtSet, iv, jv);
- if( nTrk1-nTCom==1 || nTrk2-nTCom==1){
- double probV=-1.;
- if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, xAODwrk->listJetTracks);
- else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, RECwork->listJetTracks);
- if(probV<probVrtMergeLimit)continue;
- newvrt.Good = true;
- newvrt.ProjectedVrt=JetProjDist(newvrt.vertex, PrimVrt, JetDir);
- (*WrkVrtSet).push_back(newvrt);
- (*WrkVrtSet)[iv].Good=false;
- (*WrkVrtSet)[jv].Good=false;
- }
- }
- }
+//========= Initial cleaning of solutions
//-Remove worst track from vertices with very bad Chi2
- NSoluI=(*WrkVrtSet).size();
- for(int iv=0; iv<NSoluI; iv++){
- if(!(*WrkVrtSet)[iv].Good ) continue;
- if( (*WrkVrtSet)[iv].SelTrk.size() == 2 ) continue;
- //if( (*WrkVrtSet)[iv].Chi2 > (5.*(*WrkVrtSet)[iv].SelTrk.size()) ){
- if( TMath::Prob( (*WrkVrtSet)[iv].Chi2, 2*(*WrkVrtSet)[iv].SelTrk.size()-3) <1.e-2){
+ bool disassembled=false;
+ int NSoluI=0;
+ do{
+ disassembled=false;
+ NSoluI=(*WrkVrtSet).size();
+ for(int iv=0; iv<NSoluI; iv++){
+ if(!(*WrkVrtSet)[iv].Good || (*WrkVrtSet)[iv].SelTrk.size() == 2 ) continue;
+ if( TMath::Prob( (*WrkVrtSet)[iv].Chi2, 2*(*WrkVrtSet)[iv].SelTrk.size()-3) <1.e-3){
+ //printWrkSet(WrkVrtSet,"BigChi2Vertex present");
if (xAODwrk)DisassembleVertex(WrkVrtSet,iv,xAODwrk->listJetTracks);
else if(RECwork)DisassembleVertex(WrkVrtSet,iv,RECwork->listJetTracks);
- (*WrkVrtSet)[iv].ProjectedVrt=JetProjDist((*WrkVrtSet)[iv].vertex, PrimVrt, JetDir);
- int kpp=(*WrkVrtSet).size()-1;
- (*WrkVrtSet)[kpp].ProjectedVrt=JetProjDist((*WrkVrtSet)[kpp].vertex, PrimVrt, JetDir);
- }
- }
+ disassembled=true;
+ } }
+ }while(disassembled);
//-Remove vertices fully contained in other vertices
while( (*WrkVrtSet).size()>1 ){
int tmpN=(*WrkVrtSet).size(); int iv=0;
@@ -272,32 +249,50 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
} 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++;
- 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++ ){
- if(!(*WrkVrtSet)[iv].Good ) continue;
- for(int jv=iv+1; jv<(int)(*WrkVrtSet).size(); jv++){
- 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) {
- double probV=0.;
- if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, xAODwrk->listJetTracks);
- else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, RECwork->listJetTracks);
- if(probV>probVrtMergeLimit){ // Good merged vertex found
- double tstDst=JetProjDist(newvrt.vertex, PrimVrt, JetDir);
- if(tstDst>0.){ // only positive vertex directions are accepted as merging result
- newvrt.ProjectedVrt=tstDst;
- (*WrkVrtSet).push_back(newvrt);
- (*WrkVrtSet)[iv].Good=false;
- (*WrkVrtSet)[jv].Good=false;
- break;
- }
+//
+//- Try to merge all vertices with common tracks
+ std::multimap<int,std::pair<int,int>> vrtWithCommonTrk;
+ std::multimap<int,std::pair<int,int>>::reverse_iterator icvrt;
+ do{
+ NSoluI=(*WrkVrtSet).size();
+ vrtWithCommonTrk.clear();
+ for(int iv=0; iv<NSoluI-1; iv++ ){ for(int jv=iv+1; jv<NSoluI; jv++){
+ if(!(*WrkVrtSet)[iv].Good || !(*WrkVrtSet)[jv].Good) continue;
+ int nTCom=nTrkCommon( WrkVrtSet, iv, jv); if(!nTCom)continue;
+ vrtWithCommonTrk.emplace(nTCom,std::make_pair(iv,jv));
+ } }
+ //============================== DEBUG output
+ //printWrkSet(WrkVrtSet,"InitialVrts");
+ //for(auto ku : vrtWithCommonTrk)std::cout<<" nCom="<<ku.first<<" v1="<<ku.second.first<<" v2="<<ku.second.second<<'\n';
+ //===========================================
+ for(icvrt=vrtWithCommonTrk.rbegin(); icvrt!=vrtWithCommonTrk.rend(); icvrt++){
+ int nTCom=(*icvrt).first;
+ int iv=(*icvrt).second.first;
+ int jv=(*icvrt).second.second;
+ int nTrkI=(*WrkVrtSet)[iv].SelTrk.size();
+ int nTrkJ=(*WrkVrtSet)[jv].SelTrk.size();
+ double probV=-1.;
+ if (xAODwrk)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, xAODwrk->listJetTracks);
+ else if(RECwork)probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, RECwork->listJetTracks);
+ ////std::cout<<" ntcommon="<<(*icvrt).first<<" prb="<<probV<<" itrk="<<nTrkI<<" jtrk="<<nTrkJ<<'\n';
+ if(probV<probVrtMergeLimit){
+ if(nTrkI==2 || nTrkJ==2 || nTCom<2) continue;
+ if((nTCom>nTrkI-nTCom || nTCom>nTrkJ-nTCom)){ //2 and more common tracks for NTr>=3 vertices. Merge anyway.
+ if (xAODwrk)mergeAndRefitOverlapVertices( WrkVrtSet, iv, jv, xAODwrk->listJetTracks);
+ else if(RECwork)mergeAndRefitOverlapVertices( WrkVrtSet, iv, jv, RECwork->listJetTracks);
+ break; //Vertex list is changed. Restart merging from scratch.
}
- }
+ continue; //Continue merging loop
+ }
+ newvrt.Good = true;
+ (*WrkVrtSet).push_back(newvrt);
+ (*WrkVrtSet)[iv].Good=false;
+ (*WrkVrtSet)[jv].Good=false;
+ break; //Merging successful. Break merging loop and remake list of connected vertices
}
- }
+ } while( icvrt != vrtWithCommonTrk.rend() );
+ 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); }
//-Identify/remove vertices behind the PV wrt jet direction
//-Identify remaining 2-track vertices with very bad Chi2 and mass (b-tagging)
// for(int iv=0; iv<(int)(*WrkVrtSet).size(); iv++ ){
@@ -308,9 +303,24 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
// if( (*WrkVrtSet)[iv].Chi2 > 10.) (*WrkVrtSet)[iv].Good=false;
// if( (*WrkVrtSet)[iv].vertexMom.M()>c_vrtBCMassLimit) (*WrkVrtSet)[iv].Good=false; //VK B-tagging specific requirement
// }
+//
//-Remove all bad vertices from the working set
+ for( auto &tmpV : (*WrkVrtSet) ) {
+ if(tmpV.vertex.perp()>m_rLayer3+10.)tmpV.Good=false; //Vertices outside Pixel detector
+ TLorentzVector SVPV(tmpV.vertex.x()-PrimVrt.x(),tmpV.vertex.y()-PrimVrt.y(),tmpV.vertex.z()-PrimVrt.z(),1.);
+ if(JetDir.DeltaR(SVPV)>m_coneForTag)tmpV.Good=false; // SV is outside of the jet cone
+ }
int tmpV=0; while( tmpV<(int)(*WrkVrtSet).size() )if( !(*WrkVrtSet)[tmpV].Good ) { (*WrkVrtSet).erase((*WrkVrtSet).begin()+tmpV);} else {tmpV++;}
+ if((*WrkVrtSet).size()==0){ // No vertices at all
+ delete WrkVrtSet;
+ return finalVertices;
+ }
+ //------
//printWrkSet(WrkVrtSet,"Interm");
+ //------
+ std::vector< std::vector<float> > trkScore(0);
+ if(xAODwrk){ for(auto &trk : xAODwrk->listJetTracks) trkScore.push_back(m_trackClassificator->trkTypeWgts(trk, PrimVrt, JetDir)); }
+ for( auto &tmpV : (*WrkVrtSet) ) tmpV.ProjectedVrt=JetProjDist(tmpV.vertex, PrimVrt, JetDir); //Setup ProjectedVrt
//----------------------------------------------------------------------------
// Add primary vertex
//
@@ -342,11 +352,13 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
// }
//
//----------------------------------------------------------------------------
-// Here we have the overlapping solutions
+// Here we have the overlapping solutions.
+// Vertices may have only 1 common track.
// Now solution cleaning
int nGoodVertices=0; // Final number of good vertices
- int n2trVrt=0; // N vertices with 2 and more tracks
+ int n2trVrt=0; // N of vertices with 2 tracks
+ int nNtrVrt=0; // N vertices with 3 and more tracks
std::vector< std::deque<long int> > *TrkInVrt =new std::vector< std::deque<long int> >(NTracks);
double foundMaxT; long int SelectedTrack, SelectedVertex;
int foundV1, foundV2;
@@ -425,7 +437,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
minDistVV=minVrtVrtDistNext( WrkVrtSet, foundV1, foundV2);
}
//
-// Try to improve vertices with big Chi2
+// Try to improve vertices with big Chi2 if something went wrong. Just precaution.
for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
if(!(*WrkVrtSet)[iv].Good ) continue; //don't work on vertex which is already bad
if( (*WrkVrtSet)[iv].SelTrk.size()<3 ) continue;
@@ -541,23 +553,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
}
//
//--- Check interactions on pixel layers
-//
- //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());
- 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(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(Signif3D<20.) continue;
- }
+ if(m_fillHist && nth==2){ m_hb_r2d->Fill( curVrt.vertex.perp(), m_w_1); }
//
//--- Check V0s and conversions
if(nth==2 && curVrt.vertexCharge==0 && curVrt.detachedTrack<0){
@@ -574,19 +570,9 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
//---
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(tmpProb>0.01)continue; // Vertex can be associated with PV
-// }
//---
curVrt.Good = true; /* Vertex is absolutely good */
nGoodVertices++;
- if(nth>=2)n2trVrt++;
}
//
//--Final cleaning of the 1-track vertices set. Must be behind all other cleanings.
@@ -599,13 +585,9 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
//Checks
std::vector<WrkVrt> GoodVertices(0);
nGoodVertices=0; // Final number of good vertices
- n2trVrt=0; // N vertices with 2 and more tracks
- int nNtrVrt=0; // N vertices with 3 and more tracks
+ n2trVrt=nNtrVrt=0; // N of vertices with different amount of tracks
for(auto & iv : (*WrkVrtSet) ) {
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( iv.Good) {
nGoodVertices++;
GoodVertices.emplace_back(iv); /* add it */
@@ -657,8 +639,10 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
double Signif=0.; std::vector<double> Impact,ImpactError;
for(int it=0; it<(int)nonusedTrk.size(); it++){
MatchedSV tmpV = {0, 1.e9};
- for(int iv=0; iv<(int)GoodVertices.size(); iv++){
- if(GoodVertices[iv].SelTrk.size()<2) continue;
+ for(int iv=0; iv<(int)GoodVertices.size(); iv++){ //--Find vertex closest to the given track
+ if(!GoodVertices[iv].Good) continue;
+ if(GoodVertices[iv].SelTrk.size()<2) continue;
+ if( VrtVrtDist(PrimVrt, GoodVertices[iv].vertex, GoodVertices[iv].vertexCov, JetDir) < 10.) continue; //--Too close to PV
if (RECwork)Signif = m_fitSvc->VKalGetImpact(RECwork->listJetTracks[nonusedTrk[it]], GoodVertices[iv].vertex, 1, Impact, ImpactError);
else if(xAODwrk)Signif = m_fitSvc->VKalGetImpact(xAODwrk->listJetTracks[nonusedTrk[it]], GoodVertices[iv].vertex, 1, Impact, ImpactError);
if(Signif<tmpV.Signif3D){tmpV.Signif3D=Signif; tmpV.indVrt=iv;}
@@ -674,7 +658,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
if(addTrk.size()>1){atrk++;if(atrk->first<4.){newV.SelTrk.push_back(nonusedTrk[atrk->second]);}}
if(addedT){ if (xAODwrk)vProb = RefitVertex( newV, xAODwrk->listJetTracks);
else if(RECwork)vProb = RefitVertex( newV, RECwork->listJetTracks);
- if(vProb>0.01)GoodVertices[iv]=newV;
+ if(vProb>0.01)GoodVertices[iv]=newV;
else{ std::vector<WrkVrt> TestVertices(1,newV);
if (xAODwrk)vProb=improveVertexChi2( &TestVertices, 0, xAODwrk->listJetTracks);
else if(RECwork)vProb=improveVertexChi2( &TestVertices, 0, RECwork->listJetTracks);
@@ -682,44 +666,24 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
}
}
}
-//-------------------------------------------------------------
-// Check if very close vertices are present and merge them. Not big effect in ttbar
-//
- bool vrtMerged=false; //to check whether something is really merged
- for(int iv=0; iv<(int)GoodVertices.size()-1; iv++ ){ if(!GoodVertices[iv].Good) continue;
- for(int jv=iv+1; jv<(int)GoodVertices.size(); jv++){ if(!GoodVertices[jv].Good) continue;
- if(GoodVertices[iv].SelTrk.size()>2 && GoodVertices[jv].SelTrk.size()>2)continue;
- double probV=-1.;
- if (xAODwrk)probV=mergeAndRefitVertices( &GoodVertices, iv, jv, newvrt, xAODwrk->listJetTracks);
- else if(RECwork)probV=mergeAndRefitVertices( &GoodVertices, iv, jv, newvrt, RECwork->listJetTracks);
- if(probV>0.1){ // Merged vertex is very good
- std::swap(GoodVertices[iv],newvrt);
- GoodVertices[iv].ProjectedVrt=JetProjDist(GoodVertices[iv].vertex, PrimVrt, JetDir);
- GoodVertices[jv].Good=false; //Drop vertex
- vrtMerged=true;
- } } }
- if(vrtMerged){ //-Remove all bad vertices from the good vertex set
- int iV=0; while( iV<(int)GoodVertices.size() )
- { if( !GoodVertices[iV].Good ) { GoodVertices.erase(GoodVertices.begin()+iV);} else {iV++;} }
- }
-//--
-// printWrkSet(&GoodVertices,"FinalVrtSet");
+/////
+ //if(GoodVertices.size()>=4)
+ //printWrkSet(&GoodVertices,"FinalVrtSet");
//-------------------------------------------
// Saving and final cleaning
//
nGoodVertices=0; // Final number of good vertices
int n1trVrt=0; // Final number of good 1-track vertices
TLorentzVector VertexMom;
- double trackPt, trackPtMax=0.;
for(int iv=0; iv<(int)GoodVertices.size(); iv++) {
nth=GoodVertices[iv].SelTrk.size();
if(xAODwrk)xAODwrk->tmpListTracks.clear(); else if(RECwork)RECwork->tmpListTracks.clear();
+ float vrtSumW=0.;
for(i=0;i<nth;i++) {
j=GoodVertices[iv].SelTrk[i]; /*Track number*/
if (xAODwrk)xAODwrk->tmpListTracks.push_back( xAODwrk->listJetTracks[j] );
else if(RECwork)RECwork->tmpListTracks.push_back( RECwork->listJetTracks[j] );
- trackPt=pTvsDir(Amg::Vector3D(JetDir.Vect().X(),JetDir.Vect().X(),JetDir.Vect().Z()), GoodVertices[iv].TrkAtVrt[i]);
- if(trackPt>trackPtMax)trackPtMax=trackPt;
+ if(xAODwrk)vrtSumW+=trkScore[j][0];
}
if( m_fillHist ){
if(nth==1)m_hb_r1dc->Fill( GoodVertices[iv].vertex.perp(), m_w_1);
@@ -732,7 +696,8 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
else { m_hb_totmass2T2->Fill( GoodVertices[iv].vertexMom.M(), m_w_1);}
m_hb_sig3D2tr->Fill( Signif3D , m_w_1);
if(GoodVertices[iv].vertexCharge==0)m_hb_totmassEE->Fill(massV0(GoodVertices[iv].TrkAtVrt,m_massE,m_massE),m_w_1);
- } else if( nth==1){
+ //} else if( nth==1){
+ } else if( GoodVertices[iv].vertexMom.M()>6000.){
m_hb_sig3D1tr->Fill( Signif3D, m_w_1);
} else {
m_hb_sig3DNtr->Fill( Signif3D, m_w_1);
@@ -785,7 +750,9 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
//==============================================
if(m_fillHist){m_hb_goodvrtN->Fill( nGoodVertices+0.1, m_w_1);
- m_hb_goodvrtN->Fill( n1trVrt+15., m_w_1);}
+ if(n1trVrt)m_hb_goodvrtN->Fill( n1trVrt+15., m_w_1);
+ fillNVrtNTup( GoodVertices, trkScore, PrimVrt, JetDir);
+ }
if(nGoodVertices == 0){
delete WrkVrtSet;
delete TrkInVrt;
@@ -854,12 +821,18 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
WrkVrt newvrt; newvrt.Good=true;
std::vector<const Particle*> ListBaseTracks;
int NTrk=(*WrkVrtSet)[iv].SelTrk.size(), SelT=-1;
+ if(NTrk<3)return;
+ StatusCode sc; sc.setChecked();
+//=== To get robust definition of most bad outlier
+ m_fitSvc->setRobustness(5);
+ sc = RefitVertex( WrkVrtSet, iv, AllTracks);
+ if(sc.isFailure()){ (*WrkVrtSet)[iv].Good=false; return; }
+ m_fitSvc->setRobustness(0);
+//--------------------------------------------------
double Chi2Max=0.;
for(int i=0; i<NTrk; i++){
if( (*WrkVrtSet)[iv].Chi2PerTrk[i]>Chi2Max) { Chi2Max=(*WrkVrtSet)[iv].Chi2PerTrk[i]; SelT=i;}
}
- if(SelT==-1)return;
- StatusCode sc;
int NVrtCur=WrkVrtSet->size();
for(int i=0; i<NTrk; i++){
if(i==SelT)continue;
@@ -925,24 +898,19 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
for(int mtv=0; mtv<(int)WrkVrtSet->size(); mtv++) {
if( (*WrkVrtSet)[mtv].SelTrk.size()<2 ) continue;
if(!(*WrkVrtSet)[mtv].Good ) continue;
- if( countVT[mtv] < 2 ) continue;
- double distM=1000000.;
+ if( countVT[mtv] < 1 ) continue;
+ double distM=1.e9;
int best1TVrt=-1;
for(int i1tv=0; i1tv<(int)WrkVrtSet->size(); i1tv++) {
if( (*WrkVrtSet)[i1tv].SelTrk.size()!=1 ) continue;
if(!(*WrkVrtSet)[i1tv].Good ) continue;
if( linkedVrt[i1tv] != mtv ) continue;
- double dist=((*WrkVrtSet)[mtv].vertex - (*WrkVrtSet)[i1tv].vertex).mag();
+ //double dist=((*WrkVrtSet)[mtv].vertex - (*WrkVrtSet)[i1tv].vertex).mag();
+ double dist=((*WrkVrtSet)[mtv].vertexMom+(*WrkVrtSet)[i1tv].vertexMom).M(); //Use
if( dist < distM ) {distM=dist; best1TVrt=i1tv;}
(*WrkVrtSet)[i1tv].Good=false;
}
- if(best1TVrt>-1)(*WrkVrtSet)[best1TVrt].Good=true;
- }
-//
-//--- Final Chi2 (last 2track vertex with the given track)
- for(int i1tv=0; i1tv<(int)WrkVrtSet->size(); i1tv++) {
- if( (*WrkVrtSet)[i1tv].SelTrk.size()!=1 || !(*WrkVrtSet)[i1tv].Good ) continue;
- if( (*WrkVrtSet)[i1tv].Chi2>8.) (*WrkVrtSet)[i1tv].Good=false;
+ if(best1TVrt>-1 && distM<c_vrtBCMassLimit)(*WrkVrtSet)[best1TVrt].Good=true;
}
}
@@ -951,12 +919,13 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
std::vector< std::deque<long int> > *TrkInVrt)
const
{
- int NSet=WrkVrtSet->size(); int iv, jt, NTrkAtVrt, tracknum;
- for( iv=0; iv<NSet; iv++) {
- NTrkAtVrt=(*WrkVrtSet)[iv].SelTrk.size();
+ int NSet=WrkVrtSet->size();
+ for(int iv=0; iv<NSet; iv++) {
+ if(!(*WrkVrtSet)[iv].Good) continue;
+ int NTrkAtVrt=(*WrkVrtSet)[iv].SelTrk.size();
if(NTrkAtVrt<2) continue;
- for( jt=0; jt<NTrkAtVrt; jt++){
- tracknum=(*WrkVrtSet)[iv].SelTrk[jt];
+ for(int jt=0; jt<NTrkAtVrt; jt++){
+ int tracknum=(*WrkVrtSet)[iv].SelTrk[jt];
(*TrkInVrt).at(tracknum).push_back(iv);
}
}
@@ -969,29 +938,34 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
long int & SelectedVertex)
const
{
- long int NTrack=TrkInVrt->size(); int it, jv, itmp, NVrt, VertexNumber;
+ long int NTrack=TrkInVrt->size();
+ int it, jv, itmp, NVrt, VertexNumber;
- double MaxOf=-999999999; double Chi2Red; double SelectedProb=-1.;
+ double MaxOf=-999999999, Chi2Red=0., SelectedProb=-1.;
//
- long int NShMax=0;
+ int NShareMax=0;
for( it=0; it<NTrack; it++) {
NVrt=(*TrkInVrt)[it].size(); /*Number of vertices for this track*/
- if( NVrt > NShMax ) NShMax=NVrt;
+ if( NVrt > NShareMax ) NShareMax=NVrt;
}
- if(NShMax<=1)return MaxOf; /* No shared tracks */
+ if(NShareMax<=1)return MaxOf; /* No shared tracks */
//
for( it=0; it<NTrack; it++) {
NVrt=(*TrkInVrt)[it].size(); /*Number of vertices for this track*/
if( NVrt <= 1 ) continue; /*Should ALWAYS be checked for safety*/
- if( NVrt < NShMax ) continue; /*Not a shared track with maximal sharing*/
+ if( NVrt < NShareMax ) continue; /*Not a shared track with maximal sharing*/
+ int N2trVrt=0;
+ for(auto &vrtn :(*TrkInVrt)[it] ){ if((*WrkVrtSet).at(vrtn).SelTrk.size()==2)N2trVrt++; } //count number of 2-track vertices
for( jv=0; jv<NVrt; jv++ ){
VertexNumber=(*TrkInVrt)[it][jv];
- int NTrkInVrt=(*WrkVrtSet).at(VertexNumber).SelTrk.size();
- if( NTrkInVrt <= 1) continue; // one track vertex - nothing to do
+ if(!(*WrkVrtSet).at(VertexNumber).Good)continue;
+ int NTrkInVrt=(*WrkVrtSet)[VertexNumber].SelTrk.size();
+ if( NTrkInVrt <= 1) continue; // one track vertex - nothing to do
+ if( N2trVrt>0 && N2trVrt<NShareMax && NTrkInVrt>2) continue; // Mixture of multi-track and 2-track vrts. Skip multi-track then.
for( itmp=0; itmp<NTrkInVrt; itmp++){
if( (*WrkVrtSet)[VertexNumber].SelTrk[itmp] == it ) { /* Track found*/
//Chi2Red=(*WrkVrtSet)[VertexNumber].Chi2PerTrk.at(itmp)/m_chiScale[(NTrkInVrt<10?NTrkInVrt:10)]; // Reduced Chi2
- //if(NTrkInVrt==2){ Chi2Red += 30./((*WrkVrtSet)[VertexNumber].vertex.perp()+5.);} //VK Reduce vrt multiplicity. May decrease fake rate
+ //if(NTrkInVrt==2){ Chi2Red += 30./((*WrkVrtSet)[VertexNumber].vertex.perp()+5.);} //VK Reduce vrt multiplicity.
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
@@ -999,7 +973,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
}
double prob_vrt = TMath::Prob( (*WrkVrtSet)[VertexNumber].Chi2, 2*(*WrkVrtSet)[VertexNumber].SelTrk.size()-3);
if( MaxOf < Chi2Red ){
- if(MaxOf>0 && prob_vrt>0.01 && SelectedProb<0.01 ) continue; // Don't disassemble good vertices is a bad one is present
+ if(MaxOf>0 && prob_vrt>0.01 && SelectedProb<0.01 ) continue; // Don't disassemble good vertices if a bad one is present
MaxOf = Chi2Red;
SelectedTrack=it; SelectedVertex=VertexNumber;
SelectedProb = prob_vrt;
@@ -1064,6 +1038,7 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
(*TrkInVrt)[LeftTrack].erase(it); break;
}
}
+ if( TMath::Prob( (*WrkVrtSet)[SelectedVertex].Chi2, 1) < 0.05 ) (*WrkVrtSet)[SelectedVertex].Good=false; // Not really good Chi2 for one-track vertex
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
@@ -1188,6 +1163,95 @@ const double c_vrtBCMassLimit=5500.; // Mass limit to consider a vertex not com
if( newvrt.Chi2PerTrk.size()==2) newvrt.Chi2PerTrk[0]=newvrt.Chi2PerTrk[1]=newvrt.Chi2/2.;
return TMath::Prob( newvrt.Chi2, 2*newvrt.SelTrk.size()-3);
}
+ //================== Intelligent merge of multitrack vertices with 2 and more common tracks
+ template <class Particle>
+ void InDetVKalVxInJetTool::mergeAndRefitOverlapVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2,
+ std::vector<const Particle*> & AllTrackList)
+ const
+ {
+ if(!(*WrkVrtSet).at(V1).Good)return; //bad vertex
+ if(!(*WrkVrtSet).at(V2).Good)return; //bad vertex
+ WrkVrt newvrt;
+ newvrt.Good=true;
+ if( nTrkCommon( WrkVrtSet, V1, V2)<2 )return; //No overlap
+ //- V1 should become ref. vertex. Another Vrt tracks will be added to it.
+ if( (*WrkVrtSet)[V1].SelTrk.size() < (*WrkVrtSet)[V2].SelTrk.size() )
+ {int itmp=V2; V2=V1; V1=itmp;} //Vertex with NTrk=max is chosen
+ else if( (*WrkVrtSet)[V1].SelTrk.size() == (*WrkVrtSet)[V2].SelTrk.size() )
+ { if( (*WrkVrtSet)[V1].Chi2 > (*WrkVrtSet)[V2].Chi2 ) {int itmp=V2; V2=V1; V1=itmp;} } // Vertex with minimal Chi2 is chosen
+ //- Fill base track list for new vertex
+ newvrt.SelTrk.resize( (*WrkVrtSet)[V1].SelTrk.size() );
+ std::copy((*WrkVrtSet)[V1].SelTrk.begin(),(*WrkVrtSet)[V1].SelTrk.end(), newvrt.SelTrk.begin());
+ //- Identify non-common tracks in second vertex
+ std::vector<int> noncommonTrk(0);
+ for(auto &it : (*WrkVrtSet)[V2].SelTrk){
+ if( std::find((*WrkVrtSet)[V1].SelTrk.begin(), (*WrkVrtSet)[V1].SelTrk.end(), it) == (*WrkVrtSet)[V1].SelTrk.end() )
+ noncommonTrk.push_back(it);
+ }
+ //
+ // Try to add non-common tracks one by one
+ std::vector<const Particle*> fitTrackList(0);
+ std::vector<int> detachedTrk(0);
+ StatusCode sc; sc.setChecked();
+ WrkVrt bestVrt;
+ bool foundMerged=false;
+ for( auto nct : noncommonTrk){
+ fitTrackList.clear();
+ for(int it=0; it<(int)newvrt.SelTrk.size(); it++)fitTrackList.push_back( AllTrackList[newvrt.SelTrk[it]] );
+ fitTrackList.push_back( AllTrackList.at(nct) );
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V1].vertex[0],(*WrkVrtSet)[V1].vertex[1],(*WrkVrtSet)[V1].vertex[2]);
+ sc=VKalVrtFitBase(fitTrackList, newvrt.vertex, newvrt.vertexMom, newvrt.vertexCharge, newvrt.vertexCov,
+ newvrt.Chi2PerTrk, newvrt.TrkAtVrt, newvrt.Chi2);
+ if( sc.isFailure() || TMath::Prob( newvrt.Chi2, 2*newvrt.SelTrk.size()+2-3)<0.001 ) {
+ detachedTrk.push_back(nct);
+ continue;
+ }
+ newvrt.SelTrk.push_back(nct); // Compatible track. Add to common vertex.
+ bestVrt=newvrt;
+ foundMerged=true;
+ }
+ if(foundMerged)(*WrkVrtSet)[V1]=bestVrt;
+ (*WrkVrtSet)[V2].Good=false;
+ //
+ // Now detached tracks
+ if(detachedTrk.size()>1){
+ WrkVrt nVrt;
+ fitTrackList.clear(); nVrt.SelTrk.clear();
+ for(auto nct : detachedTrk){ fitTrackList.push_back( AllTrackList[nct] ); nVrt.SelTrk.push_back(nct); }
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V2].vertex[0],(*WrkVrtSet)[V2].vertex[1],(*WrkVrtSet)[V2].vertex[2]);
+ sc=VKalVrtFitBase(fitTrackList, nVrt.vertex, nVrt.vertexMom, nVrt.vertexCharge, nVrt.vertexCov,
+ nVrt.Chi2PerTrk, nVrt.TrkAtVrt, nVrt.Chi2);
+ if(sc.isSuccess()) (*WrkVrtSet).push_back(nVrt);
+ } else if( detachedTrk.size()==1 ){
+ bool tFound=false;
+ for( auto &vrt : (*WrkVrtSet) ){
+ if( !vrt.Good || vrt.SelTrk.size()<2 )continue;
+ if( std::find(vrt.SelTrk.begin(), vrt.SelTrk.end(), detachedTrk[0]) != vrt.SelTrk.end() ) { tFound=true; break; }
+ }
+ if( !tFound ) { //Track is not present in other vertices.
+ double Chi2min=1.e9; int selectedTrk=-1;
+ WrkVrt saveVrt;
+ fitTrackList.resize(2);
+ fitTrackList[0]= AllTrackList[detachedTrk[0]];
+ for(auto trk : (*WrkVrtSet)[V2].SelTrk){
+ if(trk==detachedTrk[0])continue;
+ WrkVrt nVrt; nVrt.Good=true;
+ fitTrackList[1]=AllTrackList[trk];
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V2].vertex[0],(*WrkVrtSet)[V2].vertex[1],(*WrkVrtSet)[V2].vertex[2]);
+ sc=VKalVrtFitBase(fitTrackList, nVrt.vertex, nVrt.vertexMom, nVrt.vertexCharge, nVrt.vertexCov,
+ nVrt.Chi2PerTrk, nVrt.TrkAtVrt, nVrt.Chi2);
+ if(sc.isSuccess() && nVrt.Chi2<Chi2min) {Chi2min=nVrt.Chi2; saveVrt=nVrt; selectedTrk=trk; }
+ }
+ if(Chi2min<1.e9){
+ saveVrt.SelTrk.resize(1); saveVrt.SelTrk[0]=detachedTrk[0];
+ saveVrt.detachedTrack=selectedTrk;
+ saveVrt.vertexMom=MomAtVrt(saveVrt.TrkAtVrt[0]); //redefine vertex momentum
+ (*WrkVrtSet).push_back(saveVrt);
+ }
+ }
+ }
+ return ;
+ }
//
// Iterate track removal until vertex get good Chi2
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
index f813ae45f756a2fc42c43b5efde7ea39e64e7b25..935a4abc3710720144ff3abd6252f09ab771d79c 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/CutTrk.cxx
@@ -127,7 +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;
+ std::multimap<double,const xAOD::TrackParticle*> orderedTrk;
int NPrimTrk=0;
for (i_ntrk = InpTrk.begin(); i_ntrk < InpTrk.end(); ++i_ntrk) {
//
@@ -215,7 +215,7 @@ namespace InDet{
if( sc.isFailure() ) continue;
//double rankBTrk=RankBTrk((*i_ntrk)->pt(),JetDir.Perp(),ImpactSignif);
if(trkRank[1]>0.5)NPrimTrk += 1;
- orderedTrk[trkRank[0]]= *i_ntrk;
+ orderedTrk.emplace(trkRank[0],*i_ntrk);
}
//---- Order tracks according to ranks
std::map<double,const xAOD::TrackParticle*>::reverse_iterator rt=orderedTrk.rbegin();
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx
index c7b07e5d15b497b2c384617740ec6c7199e15ada..3da8d5e80ab0c21570b7ea84ed1c6b15c0aa0c29 100644
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetTrkInJetType.cxx
@@ -20,6 +20,9 @@ InDetTrkInJetType::InDetTrkInJetType(const std::string& type,
const IInterface* parent):
AthAlgTool(type,name,parent),
m_tmvaReader(0),
+ m_trkSctHitsCut(4),
+ m_trkPixelHitsCut(1),
+ m_trkChi2Cut(5.),
m_trkMinPtCut(700.),
m_d0_limLow(-3.),
m_d0_limUpp( 5.),
@@ -29,6 +32,9 @@ InDetTrkInJetType::InDetTrkInJetType(const std::string& type,
m_fitterSvc("Trk::TrkVKalVrtFitter/VertexFitterTool",this)
{
declareInterface<IInDetTrkInJetType>(this);
+ declareProperty("trkSctHits", m_trkSctHitsCut , "Cut on track SCT hits number" );
+ declareProperty("trkPixelHits", m_trkPixelHitsCut , "Cut on track Pixel hits number" );
+ declareProperty("trkChi2", m_trkChi2Cut , "Cut on track Chi2/Ndf" );
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" );
@@ -97,9 +103,21 @@ InDetTrkInJetType::InDetTrkInJetType(const std::string& type,
std::vector<float> InDetTrkInJetType::trkTypeWgts(const xAOD::TrackParticle * Trk, const xAOD::Vertex & PV, const TLorentzVector & Jet)
{
+//-- Track quality checks
std::vector<float> safeReturn(3,0.);
- if( !m_initialised ) return safeReturn;
- if(Jet.Perp()>2500000.)return safeReturn;
+ if( !m_initialised ) return safeReturn;
+ if(Jet.Perp() > 2500000.) return safeReturn;
+ if(Trk->pt() < m_trkMinPtCut) return safeReturn;
+ if(Trk->pt() > Jet.Pt()) return safeReturn;
+ if(Trk->numberDoF() == 0) return safeReturn; //Safety
+ if(Trk->chiSquared()/Trk->numberDoF() > m_trkChi2Cut) return safeReturn;
+ uint8_t PixelHits,SctHits;
+ if( !(Trk->summaryValue(PixelHits,xAOD::numberOfPixelHits)) ) return safeReturn; // No Pixel hits. Bad.
+ if( !(Trk->summaryValue( SctHits,xAOD::numberOfSCTHits)) ) return safeReturn; // No SCT hits. Bad.
+ if( PixelHits < m_trkPixelHitsCut ) return safeReturn;
+ if( SctHits < m_trkSctHitsCut ) return safeReturn;
+
+
std::vector<double> Impact,ImpactError;
m_Sig3D=m_fitSvc->VKalGetImpact(Trk, PV.position(), 1, Impact, ImpactError);
AmgVector(5) tmpPerigee = Trk->perigeeParameters().parameters();
@@ -110,7 +128,7 @@ InDetTrkInJetType::InDetTrkInJetType(const std::string& type,
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) );
-//---
+//---Calibrated range selection
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;
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
index b947561ad208c494f40304a845dba9f2b627ec26..4fd96f2179c2338dc6b008409aed5fc4ff822431 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/InDetVKalVxInJetTool.cxx
@@ -51,7 +51,9 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
m_a0TrkErrorCut(1.0),
m_zTrkErrorCut(5.0),
m_cutHFClass(0.1),
- m_antiGarbageCut(0.85),
+ m_antiGarbageCut(0.80),
+ m_antiFragmentCut(0.80),
+ m_Vrt2TrMassLimit(4000.),
m_fillHist(false),
m_existIBL(true),
m_RobustFit(1),
@@ -105,7 +107,9 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
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("AntiGarbageCut", m_antiGarbageCut, "Cut on Garbage classification weight for track removal" );
+ declareProperty("AntiFragmentCut",m_antiFragmentCut, "Cut on Fragmentation classification weight for track removal" );
+ declareProperty("Vrt2TrMassLimit",m_Vrt2TrMassLimit, "Maximal allowed mass for 2-track vertices" );
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" );
@@ -339,7 +343,9 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
//-------------------------------------------------------
m_curTup=new DevTuple();
m_tuple = new TTree("Tracks","Tracks");
- std::string TreeDir("/file1/stat/SVrtInJet"+m_instanceName+"/");
+ std::string TreeDir;
+ if(m_multiVertex) TreeDir="/file1/stat/MSVrtInJet"+m_instanceName+"/";
+ else 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");
@@ -363,17 +369,32 @@ InDetVKalVxInJetTool::InDetVKalVxInJetTool(const std::string& type,
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("TotM", &m_curTup->TotM, "TotM/F");
+ //-----
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("VrtDR", &m_curTup->VrtDR, "VrtDR[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");
+ //-----
+ m_tuple->Branch("nNVrt", &m_curTup->nNVrt, "nNVrt/I");
+ m_tuple->Branch("NVrtDist2D", &m_curTup->NVrtDist2D, "NVrtDist2D[nNVrt]/F");
+ m_tuple->Branch("NVrtNT", &m_curTup->NVrtNT, "NVrtNT[nNVrt]/I");
+ m_tuple->Branch("NVrtTrkI", &m_curTup->NVrtTrkI, "NVrttrkI[nNVrt]/I");
+ m_tuple->Branch("NVrtM", &m_curTup->NVrtM, "NVrtM[nNVrt]/F");
+ m_tuple->Branch("NVrtChi2", &m_curTup->NVrtChi2, "NVrtChi2[nNVrt]/F");
+ m_tuple->Branch("NVrtMaxW", &m_curTup->NVrtMaxW, "NVrtMaxW[nNVrt]/F");
+ m_tuple->Branch("NVrtAveW", &m_curTup->NVrtAveW, "NVrtAveW[nNVrt]/F");
+ m_tuple->Branch("NVrtDR", &m_curTup->NVrtDR, "NVrtDR[nNVrt]/F");
}
}
@@ -425,7 +446,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; }
+ if(m_curTup){ m_curTup->nVrt=0; m_curTup->nTrkInJet=0; m_curTup->NTHF=0; m_curTup->nNVrt=0;}
int pseudoVrt = 0;
@@ -498,7 +519,7 @@ 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; }
+ if(m_curTup){ m_curTup->nVrt=0; m_curTup->nTrkInJet=0; m_curTup->nNVrt=0; m_curTup->NTHF=0; }
xAOD::Vertex xaodPrimVrt;
xaodPrimVrt.setPosition(PrimVrt.position());
diff --git a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
index 5c322d61c54f76eaeb1226dc026a40559d98e2aa..bd65a75f02d475ccfa2d3274e9d30846848c2fe0 100755
--- a/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
+++ b/InnerDetector/InDetRecTools/InDetVKalVxInJetTool/src/Utilities.cxx
@@ -7,7 +7,6 @@
#include "TrkNeutralParameters/NeutralParameters.h"
#include "TrkTrackSummary/TrackSummary.h"
#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
-#include "xAODTruth/TruthParticleContainer.h"
//-------------------------------------------------
// Other stuff
#include <cmath>
@@ -742,8 +741,8 @@ namespace InDet{
const
{ if(!m_curTup)return;
int ipnt=0;
- for(auto vrt : all2TrVrt) {
- if(ipnt==100)break;
+ for(auto & vrt : all2TrVrt) {
+ if(ipnt==DevTuple::maxNTrk)break;
m_curTup->VrtDist2D[ipnt]=vrt.FitVertex.perp();
m_curTup->VrtSig3D[ipnt]=vrt.Signif3D;
m_curTup->VrtSig2D[ipnt]=vrt.Signif2D;
@@ -752,10 +751,83 @@ namespace InDet{
m_curTup->mass[ipnt]=vrt.Momentum.M();
m_curTup->Chi2[ipnt]=vrt.Chi2;
m_curTup->badVrt[ipnt]=vrt.badVrt;
+ m_curTup->VrtDR[ipnt]=vrt.dRSVPV;
ipnt++; m_curTup->nVrt=ipnt;
}
}
+ void InDetVKalVxInJetTool::fillNVrtNTup(std::vector<WrkVrt> & VrtSet, std::vector< std::vector<float> > & trkScore,
+ const xAOD::Vertex & PV, const TLorentzVector & JetDir)
+ const
+ { if(!m_curTup)return;
+ int ipnt=0;
+ TLorentzVector VertexMom;
+ for(auto & vrt : VrtSet) {
+ if(ipnt==DevTuple::maxNVrt)break;
+ m_curTup->NVrtDist2D[ipnt]=vrt.vertex.perp();
+ m_curTup->NVrtNT[ipnt]=vrt.SelTrk.size();
+ m_curTup->NVrtTrkI[ipnt]=vrt.SelTrk[0];
+ m_curTup->NVrtM[ipnt]=vrt.vertexMom.M();
+ m_curTup->NVrtChi2[ipnt]=vrt.Chi2;
+ float maxW=0., sumW=0.;
+ for(auto trk : vrt.SelTrk){ sumW+=trkScore[trk][0]; maxW=std::max(trkScore[trk][0], maxW);}
+ m_curTup->NVrtMaxW[ipnt]=maxW;
+ m_curTup->NVrtAveW[ipnt]=sumW/vrt.SelTrk.size();
+ TLorentzVector SVPV(vrt.vertex.x()-PV.x(),vrt.vertex.y()-PV.y(),vrt.vertex.z()-PV.z(),1.);
+ m_curTup->NVrtDR[ipnt]=JetDir.DeltaR(SVPV);
+ VertexMom += vrt.vertexMom;
+ ipnt++; m_curTup->nNVrt=ipnt;
+ }
+ m_curTup->TotM=VertexMom.M();
+ }
+
+
+ int InDetVKalVxInJetTool::getIdHF(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 0;
+ if( TP->auxdata< float >( "truthMatchProbability" ) < 0.75 ) return 0;
+ if( (*tplink)->barcode() > 200000) return 0;
+ if( (*tplink)->hasProdVtx()){
+ if( (*tplink)->prodVtx()->nIncomingParticles()==1){
+ int PDGID1=0, PDGID2=0, PDGID3=0;
+ const xAOD::TruthParticle * parTP1=getPreviousParent(*tplink, PDGID1);
+ const xAOD::TruthParticle * parTP2=0;
+ int noBC1=notFromBC(PDGID1);
+ if(noBC1) parTP2 = getPreviousParent(parTP1, PDGID2);
+ int noBC2=notFromBC(PDGID2);
+ if(noBC2 && parTP2) getPreviousParent(parTP2, PDGID3);
+ int noBC3=notFromBC(PDGID3);
+ if(noBC1 && noBC2 && noBC3)return 0;
+ return 1; //This is a reconstructed track from B/C decays
+ } } }
+ return 0;
+ }
+
+ int InDetVKalVxInJetTool::notFromBC(int PDGID) const {
+ int noBC=0;
+ if(PDGID<=0)return 1;
+ if(PDGID>600 && PDGID<4000)noBC=1;
+ if(PDGID<400 || PDGID>5600)noBC=1;
+ if(PDGID==513 || PDGID==523 || PDGID==533 || PDGID==543)noBC=1; //Remove tracks from B* (they are in PV)
+ if(PDGID==5114 || PDGID==5214 || PDGID==5224 || PDGID==5314 || PDGID==5324)noBC=1; //Remove tracks from B_Barions* (they are in PV)
+ //if(PDGID==413 || PDGID==423 || PDGID==433 )continue; //Keep tracks from D* (they are from B vertex)
+ //if(PDGID==4114 || PDGID==4214 || PDGID==4224 || PDGID==4314 || PDGID==4324)continue;
+ return noBC;
+ }
+ const xAOD::TruthParticle * InDetVKalVxInJetTool::getPreviousParent(const xAOD::TruthParticle * child, int & ParentPDG) const {
+ ParentPDG=0;
+ if( child->hasProdVtx() ){
+ if( child->prodVtx()->nIncomingParticles()==1 ){
+ ParentPDG = abs((*(child->prodVtx()->incomingParticleLinks())[0])->pdgId());
+ return *(child->prodVtx()->incomingParticleLinks())[0];
+ }
+ }
+ return 0;
+ }
+
+
int InDetVKalVxInJetTool::getG4Inter(const xAOD::TrackParticle* TP ) const {
if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
const ElementLink<xAOD::TruthParticleContainer>& tplink =
@@ -772,6 +844,7 @@ namespace InDet{
} else { return 1; }
return 0;
}
+ int InDetVKalVxInJetTool::getIdHF(const Rec::TrackParticle*) const { return 0; }
int InDetVKalVxInJetTool::getG4Inter(const Rec::TrackParticle* ) const { return 0; }
int InDetVKalVxInJetTool::getMCPileup(const Rec::TrackParticle* ) const { return 0; }