From 97a973d5989e9e58bd32fa682f5c47c7c2413600 Mon Sep 17 00:00:00 2001
From: gianelle <gianelle@pd.infn.it>
Date: Fri, 11 Oct 2024 15:31:24 +0200
Subject: [PATCH 01/12] Add new filter parameters:  - MaxHoles  - MinNdf  -
 Neural Network (BDT like) parameters

---
 source/Utils/include/FilterTracks.h |  13 ++-
 source/Utils/src/FilterTracks.cc    | 127 +++++++++++++++++++++++-----
 2 files changed, 117 insertions(+), 23 deletions(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index 3e5cfc9..3c94d6f 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -70,14 +70,25 @@ class FilterTracks : public marlin::Processor
    int _NHitsInner = 2;
    //! Cut off for number of hits in outer tracker (barrel and endcap combined)
    int _NHitsOuter = 1;
+   //! Cut off for number of holes
+   int _MaxHoles = 1;
 
    //! Cut off for momentum (GeV)
    float _MinPt = 1.0;   //units GeV
 
+   //! Cut off for the value ndf
+   int _MinNdf = 1;  
+
    //! Cut off for spatial and temporal chi squared values
    float _Chi2Spatial = 0;
 
+   //! NN parameters
+   std::string _NNmethod = ""; // if defined apply the NN
+   std::string _NNweights = ""; // xml file with weights
+   std::vector<std::string> _NNvars; // sorted list of variables used by NN
+   float _NNthr = 0; // NN threshold   
+
    //! Default magnetic field value (Tesla)
-   float _Bz = 3.57;
+   float _Bz = 5.0;
 };
 
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 4e28536..3a8b3d8 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -10,8 +10,13 @@
 #include <UTIL/CellIDDecoder.h>
 #include <UTIL/LCTrackerConf.h>
 
+#include "TMVA/Reader.h"
+
+#include <iostream> 
+
 FilterTracks aFilterTracks ;
 
+
 FilterTracks::FilterTracks()
   : Processor("FilterTracks")
 {
@@ -56,11 +61,54 @@ FilterTracks::FilterTracks()
 		 	      );
 
   registerProcessorParameter("Chi2Spatial",
-		  	     "Spatial chi squared",
+		  	     "iMinimum value for Spatial chi squared",
 			     _Chi2Spatial,
 			     _Chi2Spatial
 		 	      );
 
+  registerProcessorParameter("MinNdf",
+             "Minimum value for ndf",
+           _MinNdf,
+           _MinNdf
+            );
+
+  registerProcessorParameter("MaxHoles",
+             "Max number of holes",
+           _MaxHoles,
+           _MaxHoles
+            );
+
+  registerProcessorParameter("Bz",
+             "Magnetic field in Tesla (default: 5)",
+           _Bz,
+           _Bz
+            );
+
+  registerProcessorParameter("NNmethod",
+           "Name of the NN method, if empty uses standard cuts",
+           _NNmethod,
+           std::string("")
+            );
+
+  registerProcessorParameter("NNweights",
+             "File xml with the weights of the NN",
+           _NNweights,
+           std::string("")
+            );
+
+  registerProcessorParameter("NNvars",
+             "Sorted list with the names of the variables used in the training"
+              "Supported variables are: trtvhn, trtihn, trtohn, trthn, trtnh, trch2, trndf",
+           _NNvars,
+           _NNvars
+            );
+
+  registerProcessorParameter("NNthr",
+             "NN threshold (Takes tracks with prediction > NNthr)",
+           _NNthr,
+           _NNthr
+            );
+  
   registerInputCollection( LCIO::TRACK,
 		  	   "InputTrackCollectionName" ,
 			   "Name of the input collection",
@@ -81,7 +129,7 @@ void FilterTracks::init()
 {
   // Print the initial parameters
   printParameters() ;
-  buildBfield() ;
+  //buildBfield() ;
 }
 
 void FilterTracks::processRunHeader( LCRunHeader* /*run*/)
@@ -108,6 +156,8 @@ void FilterTracks::processEvent( LCEvent * evt )
   LCCollectionVec *OutputTrackCollection = new LCCollectionVec(LCIO::TRACK);
   OutputTrackCollection->setSubset(true);
 
+  TMVA::Reader* reader = new TMVA::Reader();
+
   // Get input collection
   LCCollection* InputTrackCollection  =evt->getCollection(_InputTrackCollection);
 
@@ -118,45 +168,78 @@ void FilterTracks::processEvent( LCEvent * evt )
   std::string encoderString = lcio::LCTrackerCellID::encoding_string();
   UTIL::CellIDDecoder<lcio::TrackerHit> decoder(encoderString);
 
+  std::unordered_map<std::string, float> vars = {
+    {"trtvhn", 0},
+    {"trtihn", 0},
+    {"trtohn", 0},
+    {"trthn", 0},
+    {"trtnh", 0},
+    {"trch2", 0},
+    {"trndf", 0},
+  };
+
+  if ( not _NNmethod.empty() ) {
+
+    for (unsigned i=0,N=_NNvars.size(); i<N; ++i) {
+      std::string name = _NNvars[i];
+      if (vars.find(name) == vars.end())
+        throw EVENT::Exception( "NN variable not supported: " + name ) ;  
+      reader->AddVariable( name, &vars[name] );
+    }
+    reader->BookMVA(_NNmethod, _NNweights);
+  }
+
   for(int i=0; i<InputTrackCollection->getNumberOfElements(); i++) {
     EVENT::Track *trk=dynamic_cast<EVENT::Track*>(InputTrackCollection->getElementAt(i));
 
-    int nhittotal  = trk->getTrackerHits().size();
+    vars["trtnh"]  = trk->getTrackerHits().size();
 
     const EVENT::IntVec& subdetectorHits = trk->getSubdetectorHitNumbers();
-    int nhitvertex = subdetectorHits[1]+subdetectorHits[2];
-    int nhitinner = subdetectorHits[3]+subdetectorHits[4];
-    int nhitouter = subdetectorHits[5]+subdetectorHits[6];
+    vars["trtvhn"] = subdetectorHits[1]+subdetectorHits[2];
+    vars["trtihn"] = subdetectorHits[3]+subdetectorHits[4];
+    vars["trtohn"] = subdetectorHits[5]+subdetectorHits[6];
 
+    vars["trthn"] = trk->getNholes();
     float pt = fabs(0.3*_Bz/trk->getOmega()/1000);
 
-    float chi2spatial = trk->getChi2();
+    vars["trch2"] = trk->getChi2();
+
+    vars["trndf"] = trk->getNdf();
 
     if(_BarrelOnly == true) {
       bool endcaphits = false;
-      for(int j=0; j<nhittotal; ++j) {
-	//Find what subdetector the hit is on 
-	uint32_t systemID = decoder(trk->getTrackerHits()[j])["system"];
-	if(systemID == 2 || systemID == 4 || systemID == 6) {
-	  endcaphits = true;
-	  break;
-	}
+      for(int j=0; j<vars["trtnh"]; ++j) {
+	      //Find what subdetector the hit is on 
+	      uint32_t systemID = decoder(trk->getTrackerHits()[j])["system"];
+	      if(systemID == 2 || systemID == 4 || systemID == 6) {
+	        endcaphits = true;
+	        break;
+	      }
       }
       if(endcaphits == false) { OutputTrackCollection->addElement(trk); }
     } else { // track property cuts
-      if(nhittotal    > _NHitsTotal  &&
-	 nhitvertex   > _NHitsVertex &&
-	 nhitinner    > _NHitsInner  &&
-	 nhitouter    > _NHitsOuter  &&
-	 pt           > _MinPt       &&
-	 chi2spatial  > _Chi2Spatial)
-	{ OutputTrackCollection->addElement(trk); }
+      if ( not _NNmethod.empty() ) { // NN cuts
+        auto prediction = reader->EvaluateMVA(_NNmethod);
+        if ( prediction > _NNthr )
+          OutputTrackCollection->addElement(trk); 
+      } else { // user cuts
+        if (vars["trtnh"]  > _NHitsTotal  &&
+	          vars["trtvhn"] > _NHitsVertex &&
+	          vars["trtihn"] > _NHitsInner  &&
+	          vars["trtohn"] > _NHitsOuter  &&
+	          pt             > _MinPt       &&
+	          vars["trch2"]  > _Chi2Spatial &&
+            vars["trndf"]  > _MinNdf      &&
+            vars["trthn"]  < _MaxHoles)
+	         { OutputTrackCollection->addElement(trk); }
+      }
     }
   }
 
   // Save output track collection
   evt->addCollection(OutputTrackCollection, _OutputTrackCollection);  
+  delete reader;
 }
 
 void FilterTracks::end()
-{ }
+{  }

From 05f4a0a64216761044ac4a91f19c2044b26518a6 Mon Sep 17 00:00:00 2001
From: gianelle <gianelle@pd.infn.it>
Date: Mon, 14 Oct 2024 16:44:12 +0200
Subject: [PATCH 02/12] Add TMVA

---
 CMakeLists.txt | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index bd8823e..5acca25 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -68,8 +68,8 @@ find_package( DD4hep REQUIRED COMPONENTS DDRec)
 set(CMAKE_MODULE_PATH  ${CMAKE_MODULE_PATH}  ${DD4hep_ROOT}/cmake )
 include( DD4hep )
 
-find_package( ROOT REQUIRED )
-set( ROOT_COMPONENT_LIBRARIES Geom Reflex)
+#find_package( ROOT REQUIRED )
+#set( ROOT_COMPONENT_LIBRARIES Geom Reflex TMVA)
 
 if(DD4HEP_USE_XERCESC)
   find_package(XercesC)
@@ -78,8 +78,7 @@ include(DD4hep_XML_setup)
 INCLUDE_DIRECTORIES( BEFORE SYSTEM ${DD4hep_INCLUDE_DIRS} )
 LINK_LIBRARIES( ${DD4hep_LIBRARIES} ${DD4hep_COMPONENT_LIBRARIES}  )
 
-#FIND_PACKAGE( ROOT REQUIRED ) 
-
+FIND_PACKAGE( ROOT REQUIRED COMPONENTS Geom TMVA) 
 INCLUDE_DIRECTORIES( SYSTEM ${ROOT_INCLUDE_DIRS} )
 LINK_LIBRARIES( ${ROOT_LIBRARIES} )
 ADD_DEFINITIONS( ${ROOT_DEFINITIONS} )

From 692194c722fd4be5d6335974100e29e85695561b Mon Sep 17 00:00:00 2001
From: gianelle <gianelle@pd.infn.it>
Date: Mon, 14 Oct 2024 17:04:54 +0200
Subject: [PATCH 03/12] Add new parameter in the filter: MaxPt

---
 source/Utils/include/FilterTracks.h |  3 ++-
 source/Utils/src/FilterTracks.cc    | 10 +++++++++-
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index 3c94d6f..564ec0f 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -74,7 +74,8 @@ class FilterTracks : public marlin::Processor
    int _MaxHoles = 1;
 
    //! Cut off for momentum (GeV)
-   float _MinPt = 1.0;   //units GeV
+   float _MinPt = 0.5;    //units GeV
+   float _MaxPt = 1000.0; //units GeV
 
    //! Cut off for the value ndf
    int _MinNdf = 1;  
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 3a8b3d8..cab4077 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -60,6 +60,12 @@ FilterTracks::FilterTracks()
 			     _MinPt
 		 	      );
 
+  registerProcessorParameter("MaxPt",
+             "Max transverse momentum",
+           _MaxPt,
+           _MaxPt
+            );
+
   registerProcessorParameter("Chi2Spatial",
 		  	     "iMinimum value for Spatial chi squared",
 			     _Chi2Spatial,
@@ -129,7 +135,8 @@ void FilterTracks::init()
 {
   // Print the initial parameters
   printParameters() ;
-  //buildBfield() ;
+  // it require that geometry is initialized
+  // buildBfield() ;
 }
 
 void FilterTracks::processRunHeader( LCRunHeader* /*run*/)
@@ -228,6 +235,7 @@ void FilterTracks::processEvent( LCEvent * evt )
 	          vars["trtihn"] > _NHitsInner  &&
 	          vars["trtohn"] > _NHitsOuter  &&
 	          pt             > _MinPt       &&
+            pt             < _MaxPt       &&
 	          vars["trch2"]  > _Chi2Spatial &&
             vars["trndf"]  > _MinNdf      &&
             vars["trthn"]  < _MaxHoles)

From 61b7e439f767f8fe76ea62dbc8da9141ca6f6f72 Mon Sep 17 00:00:00 2001
From: gianelle <gianelle@pd.infn.it>
Date: Fri, 29 Nov 2024 10:53:09 +0100
Subject: [PATCH 04/12] Add possibility to cut on theta and on hits outliers

---
 source/Utils/include/FilterTracks.h | 15 ++++++++++++++-
 source/Utils/src/FilterTracks.cc    | 30 +++++++++++++++++++++++++----
 2 files changed, 40 insertions(+), 5 deletions(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index 564ec0f..5d55b7c 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -16,7 +16,13 @@ namespace TrackPerf
  *  @parameter NHitsVertex Minimum number of hits on vertex detector
  *  @parameter NHitsInner Minimum number of hits on inner tracker
  *  @parameter NHitsOuter Minimum number of hits on outer tracker
+ *  @parameter MaxOutliers Maximum number of outliers hits on track
+ *  @parameter MaxHoles Maximum number of holes on track
+ *  @parameter MinNdf Minimum value for ndf
  *  @parameter MinPt Minimum transverse momentum
+ *  @parameter MaxPt Max transverse momentum
+ *  @parameter MinTheta Minimum theta
+ *  @parameter MaxTheta Max theta
  *  @parameter Chi2Spatial Spatial chi squared
  *
  * @author N. Bruhwiler
@@ -77,8 +83,15 @@ class FilterTracks : public marlin::Processor
    float _MinPt = 0.5;    //units GeV
    float _MaxPt = 1000.0; //units GeV
 
+   //! Cut off for theta (rad)
+   float _MinTheta = 0;    
+   float _MaxTheta = 3.14; 
+
    //! Cut off for the value ndf
-   int _MinNdf = 1;  
+   int _MinNdf = 1;
+
+   //! Cut off for outliers number   
+   int _MaxOutl = 0;  
 
    //! Cut off for spatial and temporal chi squared values
    float _Chi2Spatial = 0;
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index cab4077..c014d24 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -7,6 +7,7 @@
 #include <EVENT/Track.h>
 #include <EVENT/TrackerHit.h>
 #include <IMPL/LCCollectionVec.h>
+#include <IMPL/LCFlagImpl.h>
 #include <UTIL/CellIDDecoder.h>
 #include <UTIL/LCTrackerConf.h>
 
@@ -66,6 +67,18 @@ FilterTracks::FilterTracks()
            _MaxPt
             );
 
+  registerProcessorParameter("MinTheta",
+             "Minimum theta",
+           _MinTheta,
+           _MinTheta
+            );
+
+  registerProcessorParameter("MaxTheta",
+             "Max theta",
+           _MaxTheta,
+           _MaxTheta
+            );
+
   registerProcessorParameter("Chi2Spatial",
 		  	     "iMinimum value for Spatial chi squared",
 			     _Chi2Spatial,
@@ -84,6 +97,12 @@ FilterTracks::FilterTracks()
            _MaxHoles
             );
 
+  registerProcessorParameter("MaxOutliers",
+             "Max number of outliers",
+           _MaxOutl,
+           _MaxOutl
+          );
+
   registerProcessorParameter("Bz",
              "Magnetic field in Tesla (default: 5)",
            _Bz,
@@ -135,7 +154,7 @@ void FilterTracks::init()
 {
   // Print the initial parameters
   printParameters() ;
-  // it require that geometry is initialized
+  // it requires that geometry has been already initialized
   // buildBfield() ;
 }
 
@@ -166,8 +185,8 @@ void FilterTracks::processEvent( LCEvent * evt )
   TMVA::Reader* reader = new TMVA::Reader();
 
   // Get input collection
-  LCCollection* InputTrackCollection  =evt->getCollection(_InputTrackCollection);
-
+  LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection);
+  
   if( InputTrackCollection->getTypeName() != lcio::LCIO::TRACK )
     { throw EVENT::Exception( "Invalid collection type: " + InputTrackCollection->getTypeName() ) ; }
 
@@ -208,7 +227,7 @@ void FilterTracks::processEvent( LCEvent * evt )
 
     vars["trthn"] = trk->getNholes();
     float pt = fabs(0.3*_Bz/trk->getOmega()/1000);
-
+    float theta = M_PI_2-atan(trk->getTanLambda());
     vars["trch2"] = trk->getChi2();
 
     vars["trndf"] = trk->getNdf();
@@ -236,8 +255,11 @@ void FilterTracks::processEvent( LCEvent * evt )
 	          vars["trtohn"] > _NHitsOuter  &&
 	          pt             > _MinPt       &&
             pt             < _MaxPt       &&
+            theta          > _MinTheta    &&
+            theta          < _MaxTheta    &&
 	          vars["trch2"]  > _Chi2Spatial &&
             vars["trndf"]  > _MinNdf      &&
+            vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
             vars["trthn"]  < _MaxHoles)
 	         { OutputTrackCollection->addElement(trk); }
       }

From 74100ae6c19cf8f3b6bd568dcf64babd1f842a06 Mon Sep 17 00:00:00 2001
From: gianelle <gianelle@pd.infn.it>
Date: Wed, 11 Dec 2024 14:59:32 +0100
Subject: [PATCH 05/12] It saves real tracks, not pointers. It saves also the
 list of hits. Change two default values.

---
 source/Utils/include/FilterTracks.h |  4 ++--
 source/Utils/src/FilterTracks.cc    | 25 ++++++++++++++++++++-----
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index 5d55b7c..e11fa65 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -77,7 +77,7 @@ class FilterTracks : public marlin::Processor
    //! Cut off for number of hits in outer tracker (barrel and endcap combined)
    int _NHitsOuter = 1;
    //! Cut off for number of holes
-   int _MaxHoles = 1;
+   int _MaxHoles = 10;
 
    //! Cut off for momentum (GeV)
    float _MinPt = 0.5;    //units GeV
@@ -91,7 +91,7 @@ class FilterTracks : public marlin::Processor
    int _MinNdf = 1;
 
    //! Cut off for outliers number   
-   int _MaxOutl = 0;  
+   int _MaxOutl = 10;  
 
    //! Cut off for spatial and temporal chi squared values
    float _Chi2Spatial = 0;
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index c014d24..811599a 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -8,6 +8,7 @@
 #include <EVENT/TrackerHit.h>
 #include <IMPL/LCCollectionVec.h>
 #include <IMPL/LCFlagImpl.h>
+#include <IMPL/TrackImpl.h>
 #include <UTIL/CellIDDecoder.h>
 #include <UTIL/LCTrackerConf.h>
 
@@ -180,7 +181,13 @@ void FilterTracks::processEvent( LCEvent * evt )
 {
   // Make the output track collection
   LCCollectionVec *OutputTrackCollection = new LCCollectionVec(LCIO::TRACK);
-  OutputTrackCollection->setSubset(true);
+  // do not store pointers but real tracks
+  OutputTrackCollection->setSubset(false);
+
+  // if we want to point back to the hits we need to set the flag
+  LCFlagImpl trkFlag(0);
+  trkFlag.setBit(LCIO::TRBIT_HITS);
+  OutputTrackCollection->setFlag(trkFlag.getFlag());
 
   TMVA::Reader* reader = new TMVA::Reader();
 
@@ -242,12 +249,17 @@ void FilterTracks::processEvent( LCEvent * evt )
 	        break;
 	      }
       }
-      if(endcaphits == false) { OutputTrackCollection->addElement(trk); }
+      if (endcaphits == false) { 
+        auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
+        OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
+      }
     } else { // track property cuts
       if ( not _NNmethod.empty() ) { // NN cuts
         auto prediction = reader->EvaluateMVA(_NNmethod);
-        if ( prediction > _NNthr )
-          OutputTrackCollection->addElement(trk); 
+        if ( prediction > _NNthr ) {
+          auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
+          OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
+        }
       } else { // user cuts
         if (vars["trtnh"]  > _NHitsTotal  &&
 	          vars["trtvhn"] > _NHitsVertex &&
@@ -261,7 +273,10 @@ void FilterTracks::processEvent( LCEvent * evt )
             vars["trndf"]  > _MinNdf      &&
             vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
             vars["trthn"]  < _MaxHoles)
-	         { OutputTrackCollection->addElement(trk); }
+	          { 
+              auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
+              OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
+            }
       }
     }
   }

From 130c4ba4744090f02b6dc8f3bc0ba626203e84cf Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Wed, 5 Feb 2025 12:21:14 +0100
Subject: [PATCH 06/12] LP: fake jet removal in regional tracking filter

---
 source/Utils/include/FilterJetConeHits.h |  92 +++++
 source/Utils/src/FilterJetConeHits.cc    | 434 +++++++++++++++++++++++
 2 files changed, 526 insertions(+)
 create mode 100644 source/Utils/include/FilterJetConeHits.h
 create mode 100644 source/Utils/src/FilterJetConeHits.cc

diff --git a/source/Utils/include/FilterJetConeHits.h b/source/Utils/include/FilterJetConeHits.h
new file mode 100644
index 0000000..1c8f05f
--- /dev/null
+++ b/source/Utils/include/FilterJetConeHits.h
@@ -0,0 +1,92 @@
+#ifndef FilterJetConeHits_h
+#define FilterJetConeHits_h 1
+
+#include "marlin/Processor.h"
+#include "lcio.h"
+#include <string>
+#include <vector>
+
+#include <TH1F.h>
+
+#include <EVENT/MCParticle.h>
+#include <EVENT/ReconstructedParticle.h>
+#include "IMPL/ReconstructedParticleImpl.h"
+
+using namespace lcio ;
+using namespace marlin ;
+
+
+/** Utility processor that selects and saves the tracker hits that are included in
+ *  a DeltaR cone around the jet direction direction along with the corresponding
+ *  sim hits and the reco-sim relations.
+ *
+ *  @parameter MCParticleCollection name of the MCParticle collection
+ *  @parameter TrackerHitInputCollections name of the tracker hit input collections
+ *  @parameter TrackerSimHitInputCollections name of the tracker simhit input collections
+ *  @parameter TrackerHitInputRelations name of the tracker hit relation input collections
+ *  @parameter TrackerHitOutputCollections name of the tracker hit output collections
+ *  @parameter TrackerSimHitOutputCollections name of the tracker simhit output collections
+ *  @parameter TrackerHitOutputRelations name of the tracker hit relation output collections
+ *  @parameter DeltaRCut maximum angular distance between the hits and the particle direction
+ *  @parameter FillHistograms flag to fill the diagnostic histograms
+ *
+ * @author L. Sestini, INFN Padova
+ * @date  18 March 2021
+ * @version $Id: $
+ */
+
+class FilterJetConeHits : public Processor {
+  
+ public:
+  
+  virtual Processor*  newProcessor() { return new FilterJetConeHits ; }
+  
+  
+  FilterJetConeHits() ;
+  
+  virtual void init() ;
+  
+  virtual void processRunHeader( LCRunHeader* run ) ;
+  
+  virtual void processEvent( LCEvent * evt ) ; 
+  
+  virtual void check( LCEvent * evt ) ; 
+  
+  virtual void end() ;
+
+  bool filterJetBib(ReconstructedParticle* jet);
+  
+  
+ protected:
+
+  // --- Input/output collection names:
+  std::string m_inputJetCaloCollName{} ;
+  std::vector<std::string> m_inputTrackerHitsCollNames{} ;
+  std::vector<std::string> m_inputTrackerSimHitsCollNames{} ;
+  std::vector<std::string> m_inputTrackerHitRelNames{} ;
+  std::vector<std::string> m_outputTrackerHitsCollNames{} ;
+  std::vector<std::string> m_outputTrackerSimHitsCollNames{} ;
+  std::vector<std::string> m_outputTrackerHitRelNames{} ;
+
+
+  // --- Processor parameters:
+  bool m_fillHistos{} ;
+  double m_deltaRCut{} ;
+
+  // Jet filter parameters with BIB:
+  double m_minDaughterMaxPt{} ;
+  int m_minNTracks{} ;
+
+  // --- Diagnostic histograms:
+  TH1F* m_dist = nullptr ;
+  
+  // --- Run and event counters:
+  int _nRun{} ;
+  int _nEvt{} ;
+
+} ;
+
+#endif
+
+
+
diff --git a/source/Utils/src/FilterJetConeHits.cc b/source/Utils/src/FilterJetConeHits.cc
new file mode 100644
index 0000000..bc54335
--- /dev/null
+++ b/source/Utils/src/FilterJetConeHits.cc
@@ -0,0 +1,434 @@
+#include "FilterJetConeHits.h"
+#include <iostream>
+#include <cmath>
+#include <set>
+
+
+#include <math.h>
+
+#include <IMPL/LCCollectionVec.h>
+#include <IMPL/TrackerHitPlaneImpl.h>
+#include <IMPL/SimTrackerHitImpl.h>
+#include <IMPL/LCRelationImpl.h>
+#include <IMPL/LCFlagImpl.h>
+
+#include <marlin/AIDAProcessor.h>
+#include <marlinutil/GeometryUtil.h>
+
+#include "HelixClass_double.h"
+
+using namespace lcio ;
+using namespace marlin ;
+
+
+FilterJetConeHits aFilterJetConeHits ;
+
+
+FilterJetConeHits::FilterJetConeHits() : Processor("FilterJetConeHits") {
+
+  // --- Processor description:
+
+  _description = "FilterJetConeHits selects tracker hits in a cone opened around a jet cone direction";
+
+    
+  // --- Processor parameters:
+  
+  registerProcessorParameter("JetCaloCollection",
+			     "Name of the JetCalo collection",
+			     m_inputJetCaloCollName,
+			     std::string("JetCalo") );
+
+  registerProcessorParameter("TrackerHitInputCollections",
+			     "Name of the tracker hit input collections",
+			     m_inputTrackerHitsCollNames,
+			     {} );
+
+  registerProcessorParameter("TrackerSimHitInputCollections",
+			     "Name of the tracker simhit input collections",
+			     m_inputTrackerSimHitsCollNames,
+			     {} );
+
+  registerProcessorParameter("TrackerHitInputRelations",
+			     "Name of the tracker hit relation collections",
+			     m_inputTrackerHitRelNames,
+			     {} );
+
+  registerProcessorParameter("TrackerHitOutputCollections",
+			     "Name of the tracker hit output collections",
+			     m_outputTrackerHitsCollNames,
+			     {} );
+
+  registerProcessorParameter("TrackerSimHitOutputCollections",
+			     "Name of the tracker simhit output collections",
+			     m_outputTrackerSimHitsCollNames,
+			     {} );
+
+  registerProcessorParameter("TrackerHitOutputRelations",
+			     "Name of the tracker hit relation collections",
+			     m_outputTrackerHitRelNames,
+			     {} );
+
+  registerProcessorParameter( "DeltaRCut" ,
+			      "Maximum angular distance between the hits and the particle direction" ,
+			      m_deltaRCut,
+			      double(1.) );
+
+  registerProcessorParameter( "FillHistograms",
+			      "Flag to fill the diagnostic histograms",
+			      m_fillHistos,
+			      false );
+
+  registerProcessorParameter( "MinJetTracks",
+			      "Min number of tracks in jet to use it as filter",
+			      m_minNTracks,
+			      int(1) );  
+  
+  registerProcessorParameter( "MinDaughterMaxPt",
+			      "Min pT of the highest-pT track in jet to use it as filter",
+			      m_minDaughterMaxPt,
+			      double(2.) );  
+
+    
+}
+
+
+
+void FilterJetConeHits::init() { 
+
+  streamlog_out(DEBUG) << "   init called  " << std::endl ;
+  
+  // --- Print the processor parameters:
+
+  printParameters() ;
+
+  
+  // --- Initialize the run and event counters:
+
+  _nRun = 0 ;
+  _nEvt = 0 ;
+
+
+  // --- Initialize the AIDAProcessor and book the diagnostic histograms: 
+
+  AIDAProcessor::histogramFactory(this);
+
+  m_dist = new TH1F("m_dist", "deltaR distance between hit and jet axis", 1000, 0., 6.);
+
+}
+
+
+void FilterJetConeHits::processRunHeader( LCRunHeader* run) { 
+
+  _nRun++ ;
+
+} 
+
+
+
+void FilterJetConeHits::processEvent( LCEvent * evt ) { 
+
+  // --- Check whether the number of input and output collections match
+  bool noSimColl = false, noRelColl = false;
+  if ( m_inputTrackerSimHitsCollNames.size() == 0 ) noSimColl = true;
+  if ( m_inputTrackerHitRelNames.size() == 0 ) noRelColl = true;
+  if ( noSimColl != noRelColl ) {
+    std::stringstream err_msg;
+    err_msg << "Error: RelColl are necessary to save SimColl!"
+	    << std::endl ;
+
+    throw EVENT::Exception( err_msg.str() ) ;
+  }
+
+  if ( m_inputTrackerHitsCollNames.size() != m_outputTrackerHitsCollNames.size() ||
+      ( !noSimColl && m_inputTrackerSimHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() ) ||
+      ( !noRelColl && m_inputTrackerHitRelNames.size() != m_outputTrackerHitRelNames.size() )  ) {
+
+    std::stringstream err_msg;
+    err_msg << "Mismatch between input and output collections"
+	    << std::endl ;
+
+    throw EVENT::Exception( err_msg.str() ) ;
+  }
+
+  if ( !noSimColl && ( m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() ||
+                      m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size() ) ){
+
+    std::stringstream err_msg;
+    err_msg << "Mismatch between the reco and sim hits input collections"
+	    << std::endl ;
+
+    throw EVENT::Exception( err_msg.str() ) ;
+
+  }
+
+  if ( !noSimColl && ( m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() ||
+                      m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size() ) ){
+
+    std::stringstream err_msg;
+    err_msg << "Mismatch between the reco and sim hits output collections"
+	    << std::endl ;
+
+    throw EVENT::Exception( err_msg.str() ) ;
+
+  }
+  
+
+  // --- Get the JetCalo collection:
+
+  LCCollection* m_inputJetCalo = NULL;
+  try {
+    m_inputJetCalo = evt->getCollection( m_inputJetCaloCollName );
+  }
+  catch( lcio::DataNotAvailableException& e ) {
+    streamlog_out(WARNING) << m_inputJetCaloCollName << " collection not available" << std::endl;
+    return;
+  }
+
+
+  // --- Get the input hit collections and create the corresponding output collections:
+
+  const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size();
+  std::vector<LCCollection*> inputHitColls(nTrackerHitCol);
+  std::vector<LCCollection*> inputSimHitColls(nTrackerHitCol);
+  std::vector<LCCollection*> inputHitRels(nTrackerHitCol);
+
+  std::vector<LCCollectionVec*> outputTrackerHitColls(nTrackerHitCol);
+  std::vector<LCCollectionVec*> outputTrackerSimHitColls(nTrackerHitCol);
+  std::vector<LCCollectionVec*> outputTrackerHitRels(nTrackerHitCol);
+
+  for (unsigned int icol=0; icol<nTrackerHitCol ; ++icol) {
+
+    // get the reco hits
+    try {
+      inputHitColls[icol] = evt->getCollection(m_inputTrackerHitsCollNames[icol]);
+    }
+    catch( lcio::DataNotAvailableException& e ) {
+      streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol]
+			     << " collection not available" << std::endl;
+      continue;
+    }
+
+    // get the sim hits
+    if(!noSimColl){
+
+      try {
+        inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]);
+      }
+      catch( lcio::DataNotAvailableException& e ) {
+        streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol]
+			     << " collection not available" << std::endl;
+        continue;
+      }
+
+      try {
+        inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]);
+      }
+      catch( lcio::DataNotAvailableException& e ) {
+        streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol]
+			     << " collection not available" << std::endl;
+        continue;
+      }
+    } 
+
+    // reco hit output collections
+    std::string encoderString = inputHitColls[icol]->getParameters().getStringVal( "CellIDEncoding" );
+    outputTrackerHitColls[icol] = new LCCollectionVec( inputHitColls[icol]->getTypeName() );
+    outputTrackerHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString );
+    LCFlagImpl lcFlag(inputHitColls[icol]->getFlag());
+    outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag());
+    
+    // sim hit output collections
+    if(!noSimColl){
+      outputTrackerSimHitColls[icol] = new LCCollectionVec( inputSimHitColls[icol]->getTypeName() );
+      outputTrackerSimHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString );
+      LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag());
+      outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag());
+    
+      outputTrackerHitRels[icol] = new LCCollectionVec( inputHitRels[icol]->getTypeName() );
+      LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()) ;
+      outputTrackerHitRels[icol]->setFlag( lcFlag_rel.getFlag() ) ;
+    }
+    
+  }
+
+
+  // --- Loop over the JetCalo:
+  
+  std::vector<std::set<int> > hits_to_save(nTrackerHitCol);
+  
+  for (int ipart=0; ipart<m_inputJetCalo->getNumberOfElements(); ++ipart){
+
+    ReconstructedParticle* part = dynamic_cast<ReconstructedParticle*>( m_inputJetCalo->getElementAt(ipart) );
+
+    if( !FilterJetConeHits::filterJetBib(part) ) continue; //quality selection to reject fake jets
+
+   // double part_p = sqrt( part->getMomentum()[0]*part->getMomentum()[0] +part->getMomentum()[1]*part->getMomentum()[1] + part->getMomentum()[2]*part->getMomentum()[2] );
+
+    // --- Loop over the tracker hits and select hits inside a cone around the jet axis:
+
+    for (unsigned int icol=0; icol<inputHitColls.size(); ++icol){
+
+      LCCollection* hit_col  =  inputHitColls[icol];
+      if( !hit_col ) continue ;
+      
+      for (int ihit=0; ihit<hit_col->getNumberOfElements(); ++ihit){
+
+	TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));
+
+	// --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis:
+	if ( ( hit->getPosition()[0]*part->getMomentum()[0] +
+	       hit->getPosition()[1]*part->getMomentum()[1] +
+	       hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue;
+
+
+	// --- Get the distance between the hit and the jet axis
+	
+	double jet_p = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) + pow(part->getMomentum()[2],2)  );
+	double jet_theta = acos(part->getMomentum()[2]/jet_p);
+	double jet_eta = -std::log(tan(jet_theta/2.));
+	
+	double hit_d = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) + pow(hit->getPosition()[2],2)  );
+	double hit_theta = acos(hit->getPosition()[2]/hit_d);
+        double hit_eta = -std::log(tan(hit_theta/2.));
+
+	double jet_pxy = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) );
+	double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) );
+	double deltaPhi = acos( (part->getMomentum()[0]*hit->getPosition()[0] + part->getMomentum()[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
+
+	double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2)  );
+
+	if ( m_fillHistos ){
+
+	  m_dist->Fill(deltaR);
+
+	}
+	
+	if ( deltaR < m_deltaRCut )
+	  hits_to_save[icol].insert(ihit);
+	
+      } // ihit loop
+       
+    } // icol loop
+
+  } // ipart loop
+
+
+  // --- Add the filtered hits to the output collections:
+  
+  for (unsigned int icol=0; icol<inputHitColls.size(); ++icol){
+
+    for ( auto& ihit: hits_to_save[icol] ){
+
+      TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(inputHitColls[icol]->getElementAt(ihit));
+      TrackerHitPlaneImpl* hit_new = new TrackerHitPlaneImpl();
+
+      hit_new->setCellID0(hit->getCellID0());
+      hit_new->setCellID1(hit->getCellID1());
+      hit_new->setType(hit->getType());
+      hit_new->setPosition(hit->getPosition());
+      hit_new->setU(hit->getU());         
+      hit_new->setV(hit->getV());     
+      hit_new->setdU(hit->getdU());
+      hit_new->setdV(hit->getdV());
+      hit_new->setEDep(hit->getEDep());      
+      hit_new->setEDepError(hit->getEDepError()); 
+      hit_new->setTime(hit->getTime());
+      hit_new->setQuality(hit->getQuality());   
+
+      outputTrackerHitColls[icol]->addElement( hit_new );
+
+      if( noSimColl ) continue;
+
+      LCRelation* rel = dynamic_cast<LCRelation*>(inputHitRels[icol]->getElementAt(ihit));
+
+      
+      SimTrackerHit* simhit = dynamic_cast<SimTrackerHit*>(rel->getTo());
+      SimTrackerHitImpl* simhit_new = new SimTrackerHitImpl();
+
+      simhit_new->setCellID0(simhit->getCellID0());
+      simhit_new->setCellID1(simhit->getCellID1());
+      simhit_new->setPosition(simhit->getPosition());
+      simhit_new->setEDep(simhit->getEDep());
+      simhit_new->setTime(simhit->getTime());
+      simhit_new->setMCParticle(simhit->getMCParticle());
+      simhit_new->setMomentum(simhit->getMomentum());
+      simhit_new->setPathLength(simhit->getPathLength());
+      simhit_new->setQuality(simhit->getQuality());
+      simhit_new->setOverlay(simhit->isOverlay());
+      simhit_new->setProducedBySecondary(simhit->isProducedBySecondary());
+
+      outputTrackerSimHitColls[icol]->addElement( simhit_new );
+
+      
+      LCRelationImpl* rel_new = new LCRelationImpl();
+      
+      rel_new->setFrom(hit_new);
+      rel_new->setTo(simhit_new);
+      rel_new->setWeight(rel->getWeight());
+
+      outputTrackerHitRels[icol]->addElement( rel_new );
+
+    } // ihit loop
+
+    streamlog_out( MESSAGE ) << " " << hits_to_save[icol].size() << " hits added to the collections: "
+			     << m_outputTrackerHitsCollNames[icol] << std::endl;
+
+    evt->addCollection( outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol] ) ;
+    if(!noSimColl){
+      evt->addCollection( outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol] ) ;
+      evt->addCollection( outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol] ) ;
+    }
+
+    streamlog_out( DEBUG5 ) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type "
+			    << outputTrackerHitColls[icol]->getTypeName() << " added to the event" << std::endl;
+		if(!noSimColl) streamlog_out( DEBUG5 )  << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type "
+			    << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n"
+			    << " output collection " << m_outputTrackerHitRelNames[icol] << " of type "
+			    << outputTrackerHitRels[icol]->getTypeName() << " added to the event  "
+			    << std::endl ;
+
+  } // icol loop
+
+  streamlog_out(DEBUG) << "   processing event: " << evt->getEventNumber() 
+		       << "   in run:  " << evt->getRunNumber() << std::endl ;
+  
+  _nEvt ++ ;
+
+}
+
+
+
+void FilterJetConeHits::check( LCEvent * evt ) { 
+}
+
+
+
+void FilterJetConeHits::end(){ 
+
+  std::cout << "FilterJetConeHits::end()  " << name() 
+   	    << " processed " << _nEvt << " events in " << _nRun << " runs "
+     	    << std::endl ;
+
+}
+
+// jet quality selection
+bool FilterJetConeHits::filterJetBib(ReconstructedParticle* jet){
+
+  auto daughters = jet->getParticles();
+
+  double maxPt = 0., pt = 0.;
+  int ntracks = 0;
+
+  for(unsigned int i = 0; i < daughters.size(); i++){
+    
+    if(daughters[i]->getCharge() == 0) continue;  //only care about tracks
+    ntracks++;
+    pt = sqrt( pow( daughters[i]->getMomentum()[0] ,2) + pow( daughters[i]->getMomentum()[1] ,2) );
+    if(pt > maxPt) maxPt = pt;
+  }
+
+  if(ntracks < m_minNTracks || maxPt < m_minDaughterMaxPt) return false; //check quality of jet, if bad reject
+
+  streamlog_out(DEBUG) << "Accepted jet with ntracks = " << ntracks << " and max daughter pT = " << maxPt << std::endl;
+  return true;
+}

From 73e7f5a36028228b0134cd7b6d7c7f1872a372df Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Thu, 13 Feb 2025 15:44:02 +0100
Subject: [PATCH 07/12] LP: added new filter param in TrackFilter

---
 source/Utils/include/FilterTracks.h | 3 +++
 source/Utils/src/FilterTracks.cc    | 7 +++++++
 2 files changed, 10 insertions(+)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index e11fa65..ff1e038 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -93,6 +93,9 @@ class FilterTracks : public marlin::Processor
    //! Cut off for outliers number   
    int _MaxOutl = 10;  
 
+   //! Cut off for ratio outliers/tot hits
+   float _MaxOutlOverHits = 1.;  
+
    //! Cut off for spatial and temporal chi squared values
    float _Chi2Spatial = 0;
 
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 811599a..435ff8a 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -104,6 +104,12 @@ FilterTracks::FilterTracks()
            _MaxOutl
           );
 
+  registerProcessorParameter("MaxOutliersOverHits",
+            "Max ratio of outliers/hits",
+            _MaxOutlOverHits,
+            _MaxOutlOverHits
+         );
+
   registerProcessorParameter("Bz",
              "Magnetic field in Tesla (default: 5)",
            _Bz,
@@ -272,6 +278,7 @@ void FilterTracks::processEvent( LCEvent * evt )
 	          vars["trch2"]  > _Chi2Spatial &&
             vars["trndf"]  > _MinNdf      &&
             vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
+            (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"] < _MaxOutlOverHits &&
             vars["trthn"]  < _MaxHoles)
 	          { 
               auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);

From accf181444fcd325de52558078a3ea418ca53296 Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Fri, 21 Feb 2025 09:56:41 +0100
Subject: [PATCH 08/12] LP: added FakeRemovedJets storing

---
 source/Utils/include/FilterJetConeHits.h |   4 +
 source/Utils/src/FilterJetConeHits.cc    | 105 +++++++++++++++++------
 2 files changed, 85 insertions(+), 24 deletions(-)

diff --git a/source/Utils/include/FilterJetConeHits.h b/source/Utils/include/FilterJetConeHits.h
index 1c8f05f..5a35c73 100644
--- a/source/Utils/include/FilterJetConeHits.h
+++ b/source/Utils/include/FilterJetConeHits.h
@@ -55,12 +55,15 @@ class FilterJetConeHits : public Processor {
   virtual void end() ;
 
   bool filterJetBib(ReconstructedParticle* jet);
+
+  void saveJet( ReconstructedParticle* jet, LCCollectionVec* jetsColl );
   
   
  protected:
 
   // --- Input/output collection names:
   std::string m_inputJetCaloCollName{} ;
+  std::string m_filteredJetCaloCollName{} ;
   std::vector<std::string> m_inputTrackerHitsCollNames{} ;
   std::vector<std::string> m_inputTrackerSimHitsCollNames{} ;
   std::vector<std::string> m_inputTrackerHitRelNames{} ;
@@ -76,6 +79,7 @@ class FilterJetConeHits : public Processor {
   // Jet filter parameters with BIB:
   double m_minDaughterMaxPt{} ;
   int m_minNTracks{} ;
+  bool m_createFilteredJets{} ;
 
   // --- Diagnostic histograms:
   TH1F* m_dist = nullptr ;
diff --git a/source/Utils/src/FilterJetConeHits.cc b/source/Utils/src/FilterJetConeHits.cc
index bc54335..b171340 100644
--- a/source/Utils/src/FilterJetConeHits.cc
+++ b/source/Utils/src/FilterJetConeHits.cc
@@ -86,8 +86,17 @@ FilterJetConeHits::FilterJetConeHits() : Processor("FilterJetConeHits") {
   registerProcessorParameter( "MinDaughterMaxPt",
 			      "Min pT of the highest-pT track in jet to use it as filter",
 			      m_minDaughterMaxPt,
-			      double(2.) );  
+			      double(2.) );   
 
+  registerProcessorParameter( "FilteredJetsCollectionName",
+            "Name of the (optional) output filtered jets collection",
+            m_filteredJetCaloCollName,
+            std::string("FilteredJets") ); 
+
+  registerProcessorParameter( "MakeFilteredJetsCollection",
+            "Flag for producing collection of filtered jets",
+            m_createFilteredJets,
+            bool(false) );
     
 }
 
@@ -255,14 +264,21 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
   // --- Loop over the JetCalo:
   
   std::vector<std::set<int> > hits_to_save(nTrackerHitCol);
+
+  LCCollectionVec* filterJetsColl = new LCCollectionVec( LCIO::RECONSTRUCTEDPARTICLE );
+  LCFlagImpl flagJets( m_inputJetCalo->getFlag() );
+  filterJetsColl->setFlag( flagJets.getFlag() );
+
+  unsigned int nGoodJets = 0;
   
   for (int ipart=0; ipart<m_inputJetCalo->getNumberOfElements(); ++ipart){
 
     ReconstructedParticle* part = dynamic_cast<ReconstructedParticle*>( m_inputJetCalo->getElementAt(ipart) );
 
     if( !FilterJetConeHits::filterJetBib(part) ) continue; //quality selection to reject fake jets
+    nGoodJets++;
 
-   // double part_p = sqrt( part->getMomentum()[0]*part->getMomentum()[0] +part->getMomentum()[1]*part->getMomentum()[1] + part->getMomentum()[2]*part->getMomentum()[2] );
+    if( m_createFilteredJets ) saveJet( part , filterJetsColl); //save filtered jet
 
     // --- Loop over the tracker hits and select hits inside a cone around the jet axis:
 
@@ -273,38 +289,38 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
       
       for (int ihit=0; ihit<hit_col->getNumberOfElements(); ++ihit){
 
-	TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));
+	      TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));
 
-	// --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis:
-	if ( ( hit->getPosition()[0]*part->getMomentum()[0] +
-	       hit->getPosition()[1]*part->getMomentum()[1] +
-	       hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue;
+	      // --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis:
+        if ( ( hit->getPosition()[0]*part->getMomentum()[0] +
+              hit->getPosition()[1]*part->getMomentum()[1] +
+              hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue;
 
 
-	// --- Get the distance between the hit and the jet axis
-	
-	double jet_p = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) + pow(part->getMomentum()[2],2)  );
-	double jet_theta = acos(part->getMomentum()[2]/jet_p);
-	double jet_eta = -std::log(tan(jet_theta/2.));
-	
-	double hit_d = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) + pow(hit->getPosition()[2],2)  );
-	double hit_theta = acos(hit->getPosition()[2]/hit_d);
+        // --- Get the distance between the hit and the jet axis
+        
+        double jet_p = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) + pow(part->getMomentum()[2],2)  );
+        double jet_theta = acos(part->getMomentum()[2]/jet_p);
+        double jet_eta = -std::log(tan(jet_theta/2.));
+        
+        double hit_d = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) + pow(hit->getPosition()[2],2)  );
+        double hit_theta = acos(hit->getPosition()[2]/hit_d);
         double hit_eta = -std::log(tan(hit_theta/2.));
 
-	double jet_pxy = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) );
-	double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) );
-	double deltaPhi = acos( (part->getMomentum()[0]*hit->getPosition()[0] + part->getMomentum()[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
+        double jet_pxy = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) );
+        double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) );
+        double deltaPhi = acos( (part->getMomentum()[0]*hit->getPosition()[0] + part->getMomentum()[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
 
-	double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2)  );
+        double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2)  );
 
-	if ( m_fillHistos ){
+        if ( m_fillHistos ){
 
-	  m_dist->Fill(deltaR);
+          m_dist->Fill(deltaR);
 
-	}
+        }
 	
-	if ( deltaR < m_deltaRCut )
-	  hits_to_save[icol].insert(ihit);
+        if ( deltaR < m_deltaRCut )
+          hits_to_save[icol].insert(ihit);
 	
       } // ihit loop
        
@@ -389,6 +405,14 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
 
   } // icol loop
 
+  if( m_createFilteredJets ){
+    evt->addCollection(filterJetsColl, m_filteredJetCaloCollName);
+    streamlog_out(MESSAGE) << "Saved Filtered Jets collection, with " << filterJetsColl->getNumberOfElements() << " jets." << std::endl;
+  }
+  else{
+    delete filterJetsColl;
+  }
+
   streamlog_out(DEBUG) << "   processing event: " << evt->getEventNumber() 
 		       << "   in run:  " << evt->getRunNumber() << std::endl ;
   
@@ -432,3 +456,36 @@ bool FilterJetConeHits::filterJetBib(ReconstructedParticle* jet){
   streamlog_out(DEBUG) << "Accepted jet with ntracks = " << ntracks << " and max daughter pT = " << maxPt << std::endl;
   return true;
 }
+
+// save jet in collection
+void FilterJetConeHits::saveJet( ReconstructedParticle* jet, LCCollectionVec* jetsColl ){
+
+  // this way only a pointer to the original jet is saved
+  // if necessary implement here the phys copy save
+  //jetsColl->addElement(jet);
+  ReconstructedParticleImpl* j = new ReconstructedParticleImpl();
+  j->setType( jet->getType() );
+  j->setMomentum( jet->getMomentum() );
+  j->setEnergy( jet->getEnergy() );
+  j->setMass( jet->getMass() );
+  j->setCharge( jet->getCharge() );
+  j->setReferencePoint( jet->getReferencePoint() );
+
+  const EVENT::ReconstructedParticleVec rpvec = jet->getParticles();
+  for(unsigned int i = 0; i < rpvec.size(); i++) j->addParticle(rpvec[i]);
+
+  const EVENT::ClusterVec clvec = jet->getClusters();
+  for(unsigned int i = 0; i < clvec.size(); i++) j->addCluster(clvec[i]);
+
+  const EVENT::TrackVec trvec = jet->getTracks();
+  for(unsigned int i = 0; i < trvec.size(); i++) j->addTrack(trvec[i]);
+
+  jetsColl->addElement(j);
+
+  streamlog_out( MESSAGE ) << "Saving Jet: p = ( " << j->getMomentum()[0] << " , "
+                                                    << j->getMomentum()[1] << " , "
+                                                    << j->getMomentum()[2] << " ) " << std::endl;
+
+  return;
+}
+  
\ No newline at end of file

From 13d0d095895c5835405db959aa9846b10fec8a4d Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Mon, 24 Feb 2025 18:09:34 +0100
Subject: [PATCH 09/12] LP: add new uncertainty-based track filters

---
 source/Utils/include/FilterTracks.h |  7 +++++
 source/Utils/src/FilterTracks.cc    | 48 ++++++++++++++++++++++++++++-
 2 files changed, 54 insertions(+), 1 deletion(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index ff1e038..5060a82 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -99,6 +99,13 @@ class FilterTracks : public marlin::Processor
    //! Cut off for spatial and temporal chi squared values
    float _Chi2Spatial = 0;
 
+   //! Cut off for all track parameters uncertainties
+   float _MaxSigD0 = 999;
+   float _MaxSigZ0 = 999;
+   float _MaxSigTheta = 999;
+   float _MaxSigPhi = 999;
+   float _MaxSigQoverP = 999;
+
    //! NN parameters
    std::string _NNmethod = ""; // if defined apply the NN
    std::string _NNweights = ""; // xml file with weights
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 435ff8a..9497cdf 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -110,6 +110,36 @@ FilterTracks::FilterTracks()
             _MaxOutlOverHits
          );
 
+  registerProcessorParameter("MaxSigmaD0",
+          "Max sigma d0",
+          _MaxSigD0,
+          _MaxSigD0
+       );
+
+  registerProcessorParameter("MaxSigmaZ0",
+          "Max sigma z0",
+          _MaxSigZ0,
+          _MaxSigZ0
+      );
+
+  registerProcessorParameter("MaxSigmaTheta",
+          "Max sigma theta",
+          _MaxSigTheta,
+          _MaxSigTheta
+        );
+
+  registerProcessorParameter("MaxSigmaPhi",
+          "Max sigma phi",
+          _MaxSigPhi,
+          _MaxSigPhi
+        );
+
+  registerProcessorParameter("MaxSigmaQoverP",
+          "Max sigma q/p",
+          _MaxSigQoverP,
+          _MaxSigQoverP
+        );
+
   registerProcessorParameter("Bz",
              "Magnetic field in Tesla (default: 5)",
            _Bz,
@@ -215,6 +245,11 @@ void FilterTracks::processEvent( LCEvent * evt )
     {"trtnh", 0},
     {"trch2", 0},
     {"trndf", 0},
+    {"tr_sigd0",0},
+    {"tr_sigz0",0},
+    {"tr_sigtheta",0},
+    {"tr_sigphi",0},
+    {"tr_sigqoverp",0}
   };
 
   if ( not _NNmethod.empty() ) {
@@ -245,6 +280,12 @@ void FilterTracks::processEvent( LCEvent * evt )
 
     vars["trndf"] = trk->getNdf();
 
+    vars["tr_sigd0"] = sqrt(trk->getCovMatrix()[0]);
+    vars["tr_sigz0"] = sqrt(trk->getCovMatrix()[2]);
+    vars["tr_sigtheta"] = sqrt(trk->getCovMatrix()[9]);
+    vars["tr_sigphi"] = sqrt(trk->getCovMatrix()[5]);
+    vars["tr_sigqoverp"] = sqrt(trk->getCovMatrix()[14]);
+
     if(_BarrelOnly == true) {
       bool endcaphits = false;
       for(int j=0; j<vars["trtnh"]; ++j) {
@@ -279,7 +320,12 @@ void FilterTracks::processEvent( LCEvent * evt )
             vars["trndf"]  > _MinNdf      &&
             vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
             (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"] < _MaxOutlOverHits &&
-            vars["trthn"]  < _MaxHoles)
+            vars["trthn"]  < _MaxHoles    &&
+            vars["tr_sigd0"] < _MaxSigD0  &&
+            vars["tr_sigz0"] < _MaxSigZ0  &&
+            vars["tr_sigtheta"] < _MaxSigTheta  &&
+            vars["tr_sigphi"] < _MaxSigPhi  &&
+            vars["tr_sigqoverp"] < _MaxSigQoverP)
 	          { 
               auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
               OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 

From 97ff149f7aa2c609dccdbb5e5164c5a9e68bcc0c Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Wed, 26 Feb 2025 10:03:15 +0100
Subject: [PATCH 10/12] LP: added jet direction correction in FilterHitsJets

---
 source/Utils/include/FilterJetConeHits.h |  8 +++
 source/Utils/src/FilterJetConeHits.cc    | 88 +++++++++++++++++++++---
 2 files changed, 87 insertions(+), 9 deletions(-)

diff --git a/source/Utils/include/FilterJetConeHits.h b/source/Utils/include/FilterJetConeHits.h
index 5a35c73..ad39914 100644
--- a/source/Utils/include/FilterJetConeHits.h
+++ b/source/Utils/include/FilterJetConeHits.h
@@ -58,6 +58,7 @@ class FilterJetConeHits : public Processor {
 
   void saveJet( ReconstructedParticle* jet, LCCollectionVec* jetsColl );
   
+  void directionCorrection( const double* p, double* pcorr ) ;
   
  protected:
 
@@ -81,6 +82,13 @@ class FilterJetConeHits : public Processor {
   int m_minNTracks{} ;
   bool m_createFilteredJets{} ;
 
+  // Jet direction correction params
+  bool m_makeDirCorrection{} ;
+  double m_corrConst{};
+  double m_corrLin{};
+  double m_corrQuad{};
+  double m_corrCub{};
+
   // --- Diagnostic histograms:
   TH1F* m_dist = nullptr ;
   
diff --git a/source/Utils/src/FilterJetConeHits.cc b/source/Utils/src/FilterJetConeHits.cc
index b171340..8c0b9db 100644
--- a/source/Utils/src/FilterJetConeHits.cc
+++ b/source/Utils/src/FilterJetConeHits.cc
@@ -97,7 +97,31 @@ FilterJetConeHits::FilterJetConeHits() : Processor("FilterJetConeHits") {
             "Flag for producing collection of filtered jets",
             m_createFilteredJets,
             bool(false) );
-    
+
+  registerProcessorParameter( "MakeDirectionCorrection",
+            "Flag for correcting direction of filtered jets",
+            m_makeDirCorrection,
+             bool(false) );
+ 
+  registerProcessorParameter( "CorrectionConstant",
+            "Direction correction constant term",
+            m_corrConst,
+            double(0.) );
+
+  registerProcessorParameter( "CorrectionLinear",
+            "Direction correction linear term",
+            m_corrLin,
+            double(1.) );
+
+  registerProcessorParameter( "CorrectionQuadratic",
+            "Direction correction quadratic term",
+            m_corrQuad,
+            double(0.) );
+
+  registerProcessorParameter( "CorrectionCubic",
+            "Direction correction cubic term",
+            m_corrCub,
+            double(0.) );
 }
 
 
@@ -280,6 +304,14 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
 
     if( m_createFilteredJets ) saveJet( part , filterJetsColl); //save filtered jet
 
+    double mom[3];
+    if( !m_makeDirCorrection ){
+      for(int n = 0; n < 3; n++) mom[n] = part->getMomentum()[n]; //> to avoid const cast
+    }
+    else{
+      directionCorrection( part->getMomentum(), mom );
+    }
+
     // --- Loop over the tracker hits and select hits inside a cone around the jet axis:
 
     for (unsigned int icol=0; icol<inputHitColls.size(); ++icol){
@@ -292,24 +324,24 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
 	      TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));
 
 	      // --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis:
-        if ( ( hit->getPosition()[0]*part->getMomentum()[0] +
-              hit->getPosition()[1]*part->getMomentum()[1] +
-              hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue;
+        if ( ( hit->getPosition()[0]*mom[0] +
+              hit->getPosition()[1]*mom[1] +
+              hit->getPosition()[2]*mom[2] ) < 0. ) continue;
 
 
         // --- Get the distance between the hit and the jet axis
         
-        double jet_p = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) + pow(part->getMomentum()[2],2)  );
-        double jet_theta = acos(part->getMomentum()[2]/jet_p);
+        double jet_p = sqrt( pow(mom[0],2) + pow(mom[1],2) + pow(mom[2],2)  );
+        double jet_theta = acos(mom[2]/jet_p);
         double jet_eta = -std::log(tan(jet_theta/2.));
         
         double hit_d = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) + pow(hit->getPosition()[2],2)  );
         double hit_theta = acos(hit->getPosition()[2]/hit_d);
         double hit_eta = -std::log(tan(hit_theta/2.));
 
-        double jet_pxy = sqrt( pow(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) );
+        double jet_pxy = sqrt( pow(mom[0],2) + pow(mom[1],2) );
         double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) );
-        double deltaPhi = acos( (part->getMomentum()[0]*hit->getPosition()[0] + part->getMomentum()[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
+        double deltaPhi = acos( (mom[0]*hit->getPosition()[0] + mom[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
 
         double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2)  );
 
@@ -465,7 +497,14 @@ void FilterJetConeHits::saveJet( ReconstructedParticle* jet, LCCollectionVec* je
   //jetsColl->addElement(jet);
   ReconstructedParticleImpl* j = new ReconstructedParticleImpl();
   j->setType( jet->getType() );
-  j->setMomentum( jet->getMomentum() );
+
+  double pcorr[3];
+  if( !m_makeDirCorrection ) j->setMomentum( jet->getMomentum() );
+  else{
+    directionCorrection(jet->getMomentum(), pcorr);
+    j->setMomentum( pcorr );
+  }
+
   j->setEnergy( jet->getEnergy() );
   j->setMass( jet->getMass() );
   j->setCharge( jet->getCharge() );
@@ -486,6 +525,37 @@ void FilterJetConeHits::saveJet( ReconstructedParticle* jet, LCCollectionVec* je
                                                     << j->getMomentum()[1] << " , "
                                                     << j->getMomentum()[2] << " ) " << std::endl;
 
+  if( m_makeDirCorrection ){
+    const double* pold = jet->getMomentum();
+    double oldtheta = acos( pold[2] / sqrt(pold[0]*pold[0] + pold[1]*pold[1] + pold[2]*pold[2]) ) * 180. / 3.14159265;
+    const double* pnew = j->getMomentum();
+    double newtheta = acos( pnew[2] / sqrt(pnew[0]*pnew[0] + pnew[1]*pnew[1] + pnew[2]*pnew[2]) ) * 180. / 3.14159265;
+    streamlog_out( MESSAGE ) << "Corrected: old THETA = " << oldtheta << " => new THETA = " << newtheta << std::endl;
+  }
+
+  return;
+}
+
+void FilterJetConeHits::directionCorrection( const double* p, double* pcorr ){
+
+  double ptot = sqrt( p[0]*p[0] + p[1]*p[1] + p[2]*p[2] );
+  double theta = acos( p[2] / ptot ) * 180. / 3.14159265;
+  double theta_pre = theta;
+
+  //> only in transition region, hardcoded for now...
+  if(theta > 30. && theta < 60.){
+    theta = m_corrConst + m_corrLin * theta + m_corrQuad * pow(theta, 2.) + m_corrCub * pow(theta, 3.);
+  }
+  else if(theta > 120. && theta < 150.){
+    theta = 180. - theta;
+    theta = m_corrConst + m_corrLin * theta + m_corrQuad * pow(theta, 2.) + m_corrCub * pow(theta, 3.);
+    theta = 180. - theta;
+  }
+
+  pcorr[0] = p[0] * sin( theta * 3.14159265 / 180. ) / sin( theta_pre * 3.14159265 / 180. );
+  pcorr[1] = p[1] * sin( theta * 3.14159265 / 180. ) / sin( theta_pre * 3.14159265 / 180. );
+  pcorr[2] = p[2] * cos( theta * 3.14159265 / 180. ) / cos( theta_pre * 3.14159265 / 180. );
+
   return;
 }
   
\ No newline at end of file

From 9bc8c848dea545ed9e77506e4078bf55b6b8233d Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Wed, 17 Sep 2025 18:02:37 +0200
Subject: [PATCH 11/12] LP: new track filters and BDT variables update

---
 source/Utils/include/FilterTracks.h |  3 ++
 source/Utils/src/FilterTracks.cc    | 53 ++++++++++++++++++++---------
 2 files changed, 39 insertions(+), 17 deletions(-)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index 5060a82..fe61cda 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -106,6 +106,9 @@ class FilterTracks : public marlin::Processor
    float _MaxSigPhi = 999;
    float _MaxSigQoverP = 999;
 
+   float _MaxD0 = 999.;
+   float _MaxZ0 = 999.;
+
    //! NN parameters
    std::string _NNmethod = ""; // if defined apply the NN
    std::string _NNweights = ""; // xml file with weights
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 9497cdf..3f3c1ea 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -140,6 +140,18 @@ FilterTracks::FilterTracks()
           _MaxSigQoverP
         );
 
+  registerProcessorParameter("MaxPromptD0",
+          "Max d0",
+          _MaxD0,
+          _MaxD0
+        );
+
+  registerProcessorParameter("MaxPromptZ0",
+          "Max z0",
+          _MaxZ0,
+          _MaxZ0
+        );
+
   registerProcessorParameter("Bz",
              "Magnetic field in Tesla (default: 5)",
            _Bz,
@@ -245,11 +257,14 @@ void FilterTracks::processEvent( LCEvent * evt )
     {"trtnh", 0},
     {"trch2", 0},
     {"trndf", 0},
+    {"trnoh", 0},
     {"tr_sigd0",0},
     {"tr_sigz0",0},
     {"tr_sigtheta",0},
     {"tr_sigphi",0},
-    {"tr_sigqoverp",0}
+    {"tr_sigqoverp",0},
+    {"tr_d0",0},
+    {"tr_z0",0}
   };
 
   if ( not _NNmethod.empty() ) {
@@ -279,13 +294,16 @@ void FilterTracks::processEvent( LCEvent * evt )
     vars["trch2"] = trk->getChi2();
 
     vars["trndf"] = trk->getNdf();
-
+    vars["trnoh"] = (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"];
     vars["tr_sigd0"] = sqrt(trk->getCovMatrix()[0]);
     vars["tr_sigz0"] = sqrt(trk->getCovMatrix()[2]);
     vars["tr_sigtheta"] = sqrt(trk->getCovMatrix()[9]);
     vars["tr_sigphi"] = sqrt(trk->getCovMatrix()[5]);
     vars["tr_sigqoverp"] = sqrt(trk->getCovMatrix()[14]);
 
+    vars["tr_d0"] = trk->getD0();
+    vars["tr_z0"] = trk->getZ0();
+
     if(_BarrelOnly == true) {
       bool endcaphits = false;
       for(int j=0; j<vars["trtnh"]; ++j) {
@@ -301,14 +319,8 @@ void FilterTracks::processEvent( LCEvent * evt )
         OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
       }
     } else { // track property cuts
-      if ( not _NNmethod.empty() ) { // NN cuts
-        auto prediction = reader->EvaluateMVA(_NNmethod);
-        if ( prediction > _NNthr ) {
-          auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
-          OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
-        }
-      } else { // user cuts
-        if (vars["trtnh"]  > _NHitsTotal  &&
+
+      if (!(vars["trtnh"]  > _NHitsTotal  &&
 	          vars["trtvhn"] > _NHitsVertex &&
 	          vars["trtihn"] > _NHitsInner  &&
 	          vars["trtohn"] > _NHitsOuter  &&
@@ -319,18 +331,25 @@ void FilterTracks::processEvent( LCEvent * evt )
 	          vars["trch2"]  > _Chi2Spatial &&
             vars["trndf"]  > _MinNdf      &&
             vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
-            (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"] < _MaxOutlOverHits &&
+            vars["trnoh"] < _MaxOutlOverHits &&
             vars["trthn"]  < _MaxHoles    &&
             vars["tr_sigd0"] < _MaxSigD0  &&
             vars["tr_sigz0"] < _MaxSigZ0  &&
             vars["tr_sigtheta"] < _MaxSigTheta  &&
             vars["tr_sigphi"] < _MaxSigPhi  &&
-            vars["tr_sigqoverp"] < _MaxSigQoverP)
-	          { 
-              auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
-              OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
-            }
-      }
+            vars["tr_sigqoverp"] < _MaxSigQoverP &&
+            std::fabs(vars["tr_d0"]) < _MaxD0 &&
+            std::fabs(vars["tr_z0"]) < _MaxZ0)) continue;
+
+      if ( not _NNmethod.empty() ) { // NN cuts
+        
+        auto prediction = reader->EvaluateMVA(_NNmethod);
+        if ( prediction <= _NNthr ) continue;
+      } 
+
+      auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
+      OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); 
+      
     }
   }
 

From 24015f233a47c79b97949e3ca1f05ce273e344dc Mon Sep 17 00:00:00 2001
From: LeonardoPalombini <palombinileonardo@gmail.com>
Date: Fri, 3 Oct 2025 14:57:28 +0200
Subject: [PATCH 12/12] LP: add reduced chi2 track filter

---
 source/Utils/include/FilterTracks.h | 1 +
 source/Utils/src/FilterTracks.cc    | 9 +++++++++
 2 files changed, 10 insertions(+)

diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h
index fe61cda..50ef28e 100644
--- a/source/Utils/include/FilterTracks.h
+++ b/source/Utils/include/FilterTracks.h
@@ -89,6 +89,7 @@ class FilterTracks : public marlin::Processor
 
    //! Cut off for the value ndf
    int _MinNdf = 1;
+   float _MaxReducedChi2 = 1.e20;
 
    //! Cut off for outliers number   
    int _MaxOutl = 10;  
diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc
index 3f3c1ea..c4ff1c4 100644
--- a/source/Utils/src/FilterTracks.cc
+++ b/source/Utils/src/FilterTracks.cc
@@ -140,6 +140,12 @@ FilterTracks::FilterTracks()
           _MaxSigQoverP
         );
 
+  registerProcessorParameter("MaxReducedChi2",
+          "Max reduced Chi2",
+          _MaxReducedChi2,
+          _MaxReducedChi2
+        );
+
   registerProcessorParameter("MaxPromptD0",
           "Max d0",
           _MaxD0,
@@ -263,6 +269,7 @@ void FilterTracks::processEvent( LCEvent * evt )
     {"tr_sigtheta",0},
     {"tr_sigphi",0},
     {"tr_sigqoverp",0},
+    {"tr_redchi2",0},
     {"tr_d0",0},
     {"tr_z0",0}
   };
@@ -294,6 +301,7 @@ void FilterTracks::processEvent( LCEvent * evt )
     vars["trch2"] = trk->getChi2();
 
     vars["trndf"] = trk->getNdf();
+    vars["tr_redchi2"] = vars["trch2"] / vars["trndf"];
     vars["trnoh"] = (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"];
     vars["tr_sigd0"] = sqrt(trk->getCovMatrix()[0]);
     vars["tr_sigz0"] = sqrt(trk->getCovMatrix()[2]);
@@ -330,6 +338,7 @@ void FilterTracks::processEvent( LCEvent * evt )
             theta          < _MaxTheta    &&
 	          vars["trch2"]  > _Chi2Spatial &&
             vars["trndf"]  > _MinNdf      &&
+            vars["tr_redchi2"]  < _MaxReducedChi2  &&
             vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
             vars["trnoh"] < _MaxOutlOverHits &&
             vars["trthn"]  < _MaxHoles    &&