Merge branch 'antonioForCommit' into '21.2'

QGTaggerTool nTracks based in DerivationFramework

See merge request atlas/athena!16380

Former-commit-id: 065de8785632fc11a88c8432bd5ed73f2eb8a51c

PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/CMakeLists.txt

@@ -25,7 +25,9 @@ atlas_depends_on_subdirs( PUBLIC
                           PhysicsAnalysis/Interfaces/FTagAnalysisInterfaces
                           PhysicsAnalysis/Interfaces/JetAnalysisInterfaces
                           Reconstruction/PFlow/PFlowUtils
-                          Tools/PathResolver )
+                          Tools/PathResolver 
+			  InnerDetector/InDetRecTools/InDetTrackSelectionTool # QGTaggerTool for nTracks ###
+			  )
 
 # External dependencies:
 find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
@@ -36,7 +38,9 @@ atlas_add_library( DerivationFrameworkJetEtMissLib
                    NO_PUBLIC_HEADERS
                    INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
                    LINK_LIBRARIES ${ROOT_LIBRARIES} GaudiKernel
-                   PRIVATE_LINK_LIBRARIES AthenaBaseComps xAODCore xAODJet xAODPFlow xAODTracking xAODTrigger xAODTruth JetInterface FTagAnalysisInterfacesLib JetJvtEfficiency JetAnalysisInterfacesLib PathResolver PFlowUtils AssociationUtilsLib )
+                   PRIVATE_LINK_LIBRARIES AthenaBaseComps xAODCore xAODJet xAODPFlow xAODTracking xAODTrigger xAODTruth JetInterface FTagAnalysisInterfacesLib JetJvtEfficiency JetAnalysisInterfacesLib PathResolver PFlowUtils AssociationUtilsLib 
+InDetTrackSelectionToolLib # QGTaggerTool for nTracks ### 
+			   )
 
 atlas_add_component( DerivationFrameworkJetEtMiss
                      src/components/*.cxx

PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/python/ExtendedJetCommon.py

@@ -407,6 +407,26 @@ def addOriginCorrection(jetalg, sequence, algname,vertexPrefix):
     extjetlog.info('ExtendedJetCommon: Adding OriginCorrection for jet collection: '+jetalg)
     applyJetAugmentation(jetalg,algname,sequence,jetaugtool)
 
+### Schedule Q/G-tagging decorations ### QGTaggerTool #####
+def addQGTaggerTool(jetalg, sequence, algname ):
+    jetaugtool = getJetAugmentationTool(jetalg)
+    if(jetaugtool==None):
+        extjetlog.warning('*** addQGTaggerTool called but corresponding augmentation tool does not exist! ***')
+
+    QGTaggerToolName = 'DFQGTaggerTool' + '_InDetTrackSelectionTool_' + jetalg
+    from AthenaCommon.AppMgr import ToolSvc
+    if hasattr(ToolSvc, QGTaggerToolName):
+        jetaugtool.TrackSelectionTool = getattr(ToolSvc, QGTaggerToolName)
+    else:
+        trackselectiontool = CfgMgr.InDet__InDetTrackSelectionTool( QGTaggerToolName )
+        trackselectiontool.CutLevel = "Loose"
+        ToolSvc += trackselectiontool
+        jetaugtool.TrackSelectionTool = trackselectiontool
+
+    extjetlog.info('ExtendedJetCommon: Adding QGTaggerTool for jet collection: '+jetalg)
+    applyJetAugmentation(jetalg, algname, sequence, jetaugtool)
+#############################################################
+
 def applyOverlapRemoval(sequence=DerivationFrameworkJob):
     from AssociationUtils.config import recommended_tools
     from AssociationUtils.AssociationUtilsConf import OverlapRemovalGenUseAlg

PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/share/JETM3.py

@@ -210,6 +210,9 @@ if DerivationFrameworkIsMonteCarlo:
 from DerivationFrameworkJetEtMiss.PFlowCommon import applyPFOAugmentation
 applyPFOAugmentation(DerivationFrameworkJob)
 
+# QGTaggerTool ###
+addQGTaggerTool(jetalg="AntiKt4EMTopo", sequence=jetm3Seq, algname="QGTaggerToolAlg")
+
 #====================================================================
 # Add the containers to the output stream - slimming done here
 #====================================================================

PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/src/JetAugmentationTool.cxx

@@ -40,7 +40,12 @@ namespace DerivationFramework {
     dec_originphi(0),
     dec_originm(0),
     m_jetPtAssociationTool(""),
-    m_decorateptassociation(false)
+    m_decorateptassociation(false),
+    m_decorateQGVariables(false), // QGTaggerTool ---
+    dec_AssociatedNTracks(0),
+    dec_AssociatedTracksWidth(0),
+    dec_AssociatedTracksC1(0),
+    m_trkSelectionTool("")
   {
     declareInterface<DerivationFramework::IAugmentationTool>(this);
     declareProperty("MomentPrefix",   m_momentPrefix = "DFCommonJets_");
@@ -56,6 +61,7 @@ namespace DerivationFramework {
     declareProperty("JetTrackSumMomentsTool", m_jetTrackSumMomentsTool);
     declareProperty("JetPtAssociationTool", m_jetPtAssociationTool);
     declareProperty("JetOriginCorrectionTool",m_jetOriginCorrectionTool);
+    declareProperty("TrackSelectionTool", m_trkSelectionTool); // QGTaggerTool ---
   }
 
   StatusCode JetAugmentationTool::initialize()
@@ -107,13 +113,25 @@ namespace DerivationFramework {
     }
 
     // This tool creates the GhostTruthAssociation decorations recommended for truth matching //
-    if(!m_jetPtAssociationTool.empty()) {
-      CHECK(m_jetPtAssociationTool.retrieve());
-      ATH_MSG_INFO("Augmenting jets with GhostTruthAssociation moments Link and Fraction");
-      m_decorateptassociation = true;
-      dec_GhostTruthAssociationFraction = new SG::AuxElement::Decorator<float>("GhostTruthAssociationFraction");
-      dec_GhostTruthAssociationLink     = new SG::AuxElement::Decorator< ElementLink<xAOD::JetContainer> >("GhostTruthAssociationLink");
+    if(!m_jetPtAssociationTool.empty()) 
+      {
+	CHECK(m_jetPtAssociationTool.retrieve());
+	ATH_MSG_INFO("Augmenting jets with GhostTruthAssociation moments Link and Fraction");
+	m_decorateptassociation = true;
+	dec_GhostTruthAssociationFraction = new SG::AuxElement::Decorator<float>("GhostTruthAssociationFraction");
+	dec_GhostTruthAssociationLink     = new SG::AuxElement::Decorator< ElementLink<xAOD::JetContainer> >("GhostTruthAssociationLink");
     }
+    
+    // Here it for ntracks decoration --- QGTaggerTool ---
+    // set up InDet selection tool
+    if(!m_trkSelectionTool.empty()) {
+      CHECK( m_trkSelectionTool.retrieve() );
+      m_decorateQGVariables = true; 
+      dec_AssociatedNTracks     = new SG::AuxElement::Decorator<int>(m_momentPrefix + "QGTagger_NTracks");
+      dec_AssociatedTracksWidth = new SG::AuxElement::Decorator<float>(m_momentPrefix + "QGTagger_TracksWidth");
+      dec_AssociatedTracksC1    = new SG::AuxElement::Decorator<float>(m_momentPrefix + "QGTagger_TracksC1");
+    } // now works
+
 
     if(!m_jetOriginCorrectionTool.empty()) {
       CHECK(m_jetOriginCorrectionTool.retrieve());
@@ -153,11 +171,20 @@ namespace DerivationFramework {
       delete dec_tracksumpt;
     }
 
-    if(m_decorateptassociation){
+    if(m_decorateptassociation)
+      {
       delete dec_GhostTruthAssociationFraction;
       delete dec_GhostTruthAssociationLink;
     }
     
+    // QGTaggerTool ---
+    if(m_decorateQGVariables)
+      {
+	delete dec_AssociatedNTracks;
+	delete dec_AssociatedTracksWidth;
+        delete dec_AssociatedTracksC1;
+      }
+    
     if(m_decorateorigincorrection){
       delete dec_origincorrection;
       delete dec_originpt;
@@ -221,6 +248,7 @@ namespace DerivationFramework {
       }
     }
 
+
     // loop over the copies
     for(const auto& jet : *jets_copy) {
       // get the original jet so we can decorate it
@@ -274,8 +302,108 @@ namespace DerivationFramework {
         (*dec_GhostTruthAssociationLink)(jet_orig) = jet->getAttribute< ElementLink<xAOD::JetContainer> >("GhostTruthAssociationLink");
         ATH_MSG_VERBOSE("GhostTruthAssociationLink: " << (*dec_GhostTruthAssociationLink)(jet_orig) );
       }
-    }
+
+      // QGTaggerTool ---
+      if(m_decorateQGVariables)
+	{
+	  ATH_MSG_DEBUG("Test Decorate QG ");
+	  std::vector<const xAOD::IParticle*> jettracks;
+	  jet->getAssociatedObjects<xAOD::IParticle>(xAOD::JetAttribute::GhostTrack,jettracks);
+	  
+	  int nTracksCount = 0;
+	  double TracksWidth = 0., SumTracks_pTs = 0., TracksC1 = 0., beta = 0.2;
+	  bool invalidJet = false;
+	  
+	  const xAOD::Vertex *pv = 0;
+	  const xAOD::VertexContainer* vxCont = 0;
+	  if(evtStore()->retrieve( vxCont, "PrimaryVertices" ).isFailure()){
+	    ATH_MSG_WARNING("Unable to retrieve primary vertex container PrimaryVertices");
+	    nTracksCount = -1;
+	    TracksWidth = -1.;
+            TracksC1 = -1.;
+	    invalidJet = true;
+	  }
+	  else if(vxCont->empty()){
+	    ATH_MSG_WARNING("Event has no primary vertices!");
+	    nTracksCount = -1;
+	    TracksWidth = -1.;
+            TracksC1 = -1.;
+	    invalidJet = true;
+	  }
+	  else{
+	    for(const auto& vx : *vxCont){
+	      // take the first vertex in the list that is a primary vertex
+	      if(vx->vertexType()==xAOD::VxType::PriVtx){
+		pv = vx;
+		break;
+	      }
+	    }
+	  }
+
+	  std::vector<bool> IsGoodTrack;
+          TLorentzVector tracki_TLV, trackj_TLV;
+          TLorentzVector jet_TLV = jet -> p4();
+	  for (size_t i = 0; i < jettracks.size(); i++) {
+	    
+	    if(invalidJet) continue;
+	    
+	    const xAOD::TrackParticle* trk = static_cast<const xAOD::TrackParticle*>(jettracks[i]);
+
+	    // only count tracks with selections
+	    // 1) pt>500 MeV
+	    // 2) accepted track from InDetTrackSelectionTool with CutLevel==Loose
+	    // 3) associated to primary vertex OR within 3mm of the primary vertex
+	    bool accept = (trk->pt()>500 &&
+			   m_trkSelectionTool->accept(*trk) &&
+			   (trk->vertex()==pv || (!trk->vertex() && fabs((trk->z0()+trk->vz()-pv->z())*sin(trk->theta()))<3.))
+			   );	    
+	    
+	    ATH_MSG_DEBUG("Test Decorate QG: trkSelTool output " << m_trkSelectionTool->accept(*trk) );
+
+	    IsGoodTrack.push_back(accept);
+	    if (!accept) continue;
+	    
+	    nTracksCount++;
+
+	    tracki_TLV = trk -> p4();
+	    double DR_tracki_jet = tracki_TLV.DeltaR(jet_TLV);
+            TracksWidth += trk -> pt() * DR_tracki_jet;
+            SumTracks_pTs += trk -> pt();
+	    
+	  }// end loop over jettracks
+	  
+	  if(SumTracks_pTs>0.) TracksWidth = TracksWidth / SumTracks_pTs;
+	  else TracksWidth = -1.;
+
+	  for(size_t i = 0; i < jettracks.size(); i++) {
+            if(invalidJet) continue;
+            const xAOD::TrackParticle* trki = static_cast<const xAOD::TrackParticle*>(jettracks[i]);
+	    if( !( IsGoodTrack.at(i) ) ) continue;
+	    
+            for(size_t j = i+1; j < jettracks.size(); j++) {
+              const xAOD::TrackParticle* trkj = static_cast<const xAOD::TrackParticle*>(jettracks[j]);
+	      if( !( IsGoodTrack.at(j) ) ) continue;
+
+              tracki_TLV = trki -> p4();
+              trackj_TLV = trkj -> p4();
+              double DR_tracki_trackj = tracki_TLV.DeltaR(trackj_TLV);
+              TracksC1 += trki -> pt() * trkj -> pt() * pow( DR_tracki_trackj, beta) ;
+
+            }//end loop over j
+          }//end double loop over ij
+
+          if(SumTracks_pTs>0.) TracksC1 = TracksC1 / ( pow(SumTracks_pTs, 2.) );
+	  else TracksC1 = -1.;
+
+	  (*dec_AssociatedNTracks)(jet_orig)     = nTracksCount;
+	  (*dec_AssociatedTracksWidth)(jet_orig) = TracksWidth;
+          (*dec_AssociatedTracksC1)(jet_orig)    = TracksC1;
+
+	}// end if m_decorateQGVariables
+
+    }//end loop on jets copies
 
     return StatusCode::SUCCESS;
-  }
+  }//end addBranches
+
 }

PhysicsAnalysis/DerivationFramework/DerivationFrameworkJetEtMiss/src/JetAugmentationTool.h

@@ -23,6 +23,8 @@
 #include "FTagAnalysisInterfaces/IBTaggingSelectionTool.h"
 #include "xAODJet/JetContainer.h"
 
+#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h" // QGTaggerTool
+
 namespace DerivationFramework {
 
   class JetAugmentationTool : public AthAlgTool, public IAugmentationTool {
@@ -39,6 +41,7 @@ namespace DerivationFramework {
     //
     // implement augmentations explicitly to avoid need to parse lists of moments to copy
     //
+    
     // calibration
     SG::AuxElement::Decorator<float>* dec_calibpt;
     SG::AuxElement::Decorator<float>* dec_calibeta;
@@ -84,6 +87,14 @@ namespace DerivationFramework {
     bool m_decorateptassociation;
     SG::AuxElement::Decorator<float>* dec_GhostTruthAssociationFraction;
     SG::AuxElement::Decorator< ElementLink<xAOD::JetContainer> >* dec_GhostTruthAssociationLink;
+
+    // Ntracks for QGTaggerTool ---
+    bool m_decorateQGVariables;
+    SG::AuxElement::Decorator<int>* dec_AssociatedNTracks;
+    SG::AuxElement::Decorator<float>* dec_AssociatedTracksWidth;
+    SG::AuxElement::Decorator<float>* dec_AssociatedTracksC1;
+    ToolHandle<InDet::IInDetTrackSelectionTool> m_trkSelectionTool;
+
   };
 }
 

PhysicsAnalysis/DerivationFramework/DerivationFrameworkSM/share/STDM6.py

@@ -86,6 +86,9 @@ addHadRecoilMETMap(STDM6Sequence, STDM6Stream, "STDM6")
 #svcMgr += createThinningSvc( svcName="STDM6ThinningSvc", outStreams=[evtStream] )
 
 
+# QGTaggerTool ### 
+addQGTaggerTool(jetalg="AntiKt4EMTopo", sequence=STDM6Sequence, algname="QGTaggerToolAlg")
+
 #====================================================================
 # Jet reconstruction/retagging
 #====================================================================
@@ -93,10 +96,10 @@ addHadRecoilMETMap(STDM6Sequence, STDM6Stream, "STDM6")
 #re-tag PFlow jets so they have b-tagging info.
 FlavorTagInit(JetCollections = ['AntiKt4EMPFlowJets'], Sequencer = STDM6Sequence)
 
-
 #====================================================================
 # Add the containers to the output stream - slimming done here
 #====================================================================
+
 from DerivationFrameworkCore.SlimmingHelper import SlimmingHelper
 from DerivationFrameworkSM.STDMExtraContent import *
 
@@ -140,5 +143,3 @@ addJetOutputs(STDM6SlimmingHelper,["STDM6","STDM6Jets"])
 addMETOutputs(STDM6SlimmingHelper,["AntiKt4EMPFlow"])
 
 STDM6SlimmingHelper.AppendContentToStream(STDM6Stream)
-
-

PhysicsAnalysis/DerivationFramework/DerivationFrameworkSM/share/STDM9.py

@@ -20,6 +20,7 @@ from DerivationFrameworkCore.LHE3WeightMetadata import *
 
 stdm9Seq = CfgMgr.AthSequencer("STDM9Sequence")
 
+
 #====================================================================
 # SET UP STREAM   
 #====================================================================
@@ -174,6 +175,9 @@ stdm9Seq += CfgMgr.DerivationFramework__DerivationKernel("STDM9Kernel", Thinning
 
 addRscanJets("AntiKt",0.6,"LCTopo",stdm9Seq,"STDM9")
 
+# QGTaggerTool ###
+addQGTaggerTool(jetalg="AntiKt4EMTopo", sequence=stdm9Seq, algname="QGTaggerToolAlg")
+
 #====================================================================
 # Add the containers to the output stream - slimming done here
 #====================================================================