added cut on outliers

source/Refitting/include/RefitFinal.h

@@ -91,6 +91,7 @@ class RefitFinal : public marlin::Processor {
 
   bool _extrapolateForward = true;
   int _minClustersOnTrackAfterFit = 0;
+  int _maxOutliersAllowed = 99;
 
   float _initialTrackError_d0{};
   float _initialTrackError_phi0{};

source/Refitting/src/RefitFinal.cc

@@ -8,6 +8,8 @@
 #include <MarlinTrk/IMarlinTrack.h>
 #include <MarlinTrk/MarlinTrkUtils.h>
 
+#include <DD4hep/Detector.h>
+
 #include <EVENT/TrackerHit.h>
 #include <IMPL/LCCollectionVec.h>
 #include <IMPL/LCFlagImpl.h>
@@ -85,6 +87,9 @@ RefitFinal::RefitFinal() : Processor("RefitFinal") {
   registerProcessorParameter("MinClustersOnTrackAfterFit", "Final minimum number of track clusters",
                              _minClustersOnTrackAfterFit, int(4));
 
+  registerProcessorParameter("MaxOutliersAllowed", "Maximum number of outliers allowed on the refitted track",
+                             _maxOutliersAllowed, int(99));
+
   registerProcessorParameter("ReducedChi2Cut", "Cut on maximum allowed reduced chi2",
                              _ReducedChi2Cut, double(-1.0));
 }
@@ -177,24 +182,6 @@ void RefitFinal::processEvent(LCEvent* evt) {
       ++hitInSubDet[_encoder->operator[](UTIL::LCTrackerCellID::subdet())];
     }
 
-    //int init_status = FitInit2(track, marlin_trk.get());
-
-    //if (init_status != 0) {
-    //  continue;
-    //}
-
-    //streamlog_out(DEBUG4) << "Refit: Trackstate after initialisation\n" << marlin_trk->toString() << std::endl;
-
-    //streamlog_out(DEBUG5) << "track initialised " << std::endl;
-
-    //int fit_status = marlin_trk->fit();
-
-    //streamlog_out(DEBUG4) << "RefitHit: Trackstate after fit()\n" << marlin_trk->toString() << std::endl;
-
-    //if (fit_status != 0) {
-    //  continue;
-    //}
-
     auto lcio_trk = std::unique_ptr<IMPL::TrackImpl>(new IMPL::TrackImpl());
 
     // setup initial dummy covariance matrix
@@ -211,6 +198,13 @@ void RefitFinal::processEvent(LCEvent* evt) {
     covMatrix[9] = (_initialTrackError_z0);    // sigma_z0^2
     covMatrix[14] = (_initialTrackError_tanL); // sigma_tanl^2
 
+    // get magnetic field at the origin
+    dd4hep::Detector& lcdd = dd4hep::Detector::getInstance();
+    const double position[3] = {0, 0, 0};      // position to calculate magnetic field at (the origin in this case)
+    double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field
+    lcdd.field().magneticField(position, magneticFieldVector); // get the magnetic field vector from DD4hep
+    _bField = magneticFieldVector[2] / dd4hep::tesla;
+
     const bool fit_direction = _extrapolateForward ? MarlinTrk::IMarlinTrack::forward : MarlinTrk::IMarlinTrack::backward;
     int return_code = MarlinTrk::createFinalisedLCIOTrack(marlin_trk.get(), trkHits, lcio_trk.get(), fit_direction, covMatrix, _bField,
                                                     _Max_Chi2_Incr);
@@ -242,6 +236,14 @@ void RefitFinal::processEvent(LCEvent* evt) {
 
     marlin_trk->getOutliers(outliers);
 
+    if (int(outliers.size()) > _maxOutliersAllowed) {
+      streamlog_out(DEBUG3) << "More than " << _maxOutliersAllowed
+                            << " outliers: Track "
+                               "Discarded. Number of outliers =  "
+                            << outliers.size() << std::endl;
+      continue;
+    }
+
     std::vector<TrackerHit*> all_hits;
     all_hits.reserve(hits_in_fit.size() + outliers.size());
 
@@ -316,24 +318,3 @@ LCCollection* RefitFinal::GetCollection(LCEvent* evt, std::string colName) {
 
   return col;
 }
-
-int RefitFinal::FitInit2(Track* track, MarlinTrk::IMarlinTrack* marlinTrk) {
-  TrackStateImpl trackState;
-
-  if (_refPoint == -1) {
-    trackState = TrackStateImpl(TrackState::AtOther, track->getD0(), track->getPhi(), track->getOmega(), track->getZ0(),
-                                track->getTanLambda(), track->getCovMatrix(), track->getReferencePoint());
-  } else {
-    const TrackState* trackAtHit = track->getTrackState(_refPoint);
-    if (not trackAtHit) {
-      streamlog_out(ERROR) << "Cannot find trackstate for " << _refPoint << std::endl;
-      return MarlinTrk::IMarlinTrack::error;
-    }
-    trackState = TrackStateImpl(*trackAtHit);
-  }
-
-  const bool direction = _extrapolateForward ? MarlinTrk::IMarlinTrack::forward : MarlinTrk::IMarlinTrack::backward;
-  marlinTrk->initialise(trackState, _bField, direction);
-
-  return MarlinTrk::IMarlinTrack::success;
-}

source/Utils/include/FilterTracks.h

@@ -52,9 +52,11 @@ class FilterTracks : public marlin::Processor {
 private:
   //! Input track collection
   std::string _InputTrackCollection{};
+  std::string _InputTrackRelationCollection{};
 
   //! Output track collection
   std::string _OutputTrackCollection{};
+  std::string _OutputTrackRelationCollection{};
 
   bool _BarrelOnly = false;
   bool _HasCaloState = false;
@@ -71,15 +73,16 @@ class FilterTracks : public marlin::Processor {
   //! Cut off for maximum number of holes on track
   int _MaxHoles = 0;
 
-  //! Cut off for maximum number of outliers on track
-  int _MaxOutliers = 20;
-
   //! Cut off for momentum (GeV)
   float _MinPt = 1.0; // units GeV
 
   //! Cut off for spatial and temporal chi squared values
   float _Chi2Spatial = 0;
 
+  //! Cut off for impact parameters [mm]
+  float _MaxD0 = 999.;
+  float _MaxZ0 = 999.;
+
   //! Default magnetic field value (Tesla)
   float _Bz = 5.0;
 };

source/Utils/src/FilterTracks.cc

@@ -9,6 +9,10 @@
 #include <IMPL/LCCollectionVec.h>
 #include <UTIL/CellIDDecoder.h>
 #include <UTIL/LCTrackerConf.h>
+#include <EVENT/MCParticle.h>
+#include <IMPL/LCRelationImpl.h>
+
+#include <UTIL/LCRelationNavigator.h>
 
 FilterTracks aFilterTracks;
 
@@ -18,6 +22,20 @@ FilterTracks::FilterTracks() : Processor("FilterTracks") {
                  "filtered collection";
 
   // register steering parameters: name, description, class-variable, default value
+  registerInputCollection(LCIO::TRACK, "InputTrackCollectionName", "Name of the input collection",
+                          _InputTrackCollection, _InputTrackCollection);
+
+  registerOutputCollection(LCIO::TRACK, "OutputTrackCollectionName", "Name of output collection",
+                           _OutputTrackCollection, std::string("FilteredTracks"));
+
+  registerInputCollection(LCIO::LCRELATION, "InputRelationCollectionName",
+                          "Name of the input track to MCParticle relation collection", _InputTrackRelationCollection,
+                          _InputTrackRelationCollection);
+
+  registerOutputCollection(LCIO::LCRELATION, "OutputRelationCollectionName",
+                           "Refit Track to MCParticle relation collection Name", _OutputTrackRelationCollection,
+                           _OutputTrackRelationCollection);
+
   registerProcessorParameter("BarrelOnly", "If true, just keep tracks with only barrel hits", _BarrelOnly, _BarrelOnly);
 
   registerProcessorParameter("HasCaloState",
@@ -34,17 +52,13 @@ FilterTracks::FilterTracks() : Processor("FilterTracks") {
 
   registerProcessorParameter("MaxHoles", "Maximum number of holes on track", _MaxHoles, _MaxHoles);
 
-  registerProcessorParameter("MaxOutliers", "Max number of outliers", _MaxOutliers, _MaxOutliers);
-
   registerProcessorParameter("MinPt", "Minimum transverse momentum", _MinPt, _MinPt);
 
-  registerProcessorParameter("Chi2Spatial", "Spatial chi squared", _Chi2Spatial, _Chi2Spatial);
+  registerProcessorParameter("MaxD0", "Max d0", _MaxD0, _MaxD0);
 
-  registerInputCollection(LCIO::TRACK, "InputTrackCollectionName", "Name of the input collection",
-                          _InputTrackCollection, _InputTrackCollection);
+  registerProcessorParameter("MaxZ0", "Max z0", _MaxZ0, _MaxZ0);
 
-  registerOutputCollection(LCIO::TRACK, "OutputTrackCollectionName", "Name of output collection",
-                           _OutputTrackCollection, std::string("FilteredTracks"));
+  registerProcessorParameter("Chi2Spatial", "Spatial chi squared", _Chi2Spatial, _Chi2Spatial);
 }
 
 void FilterTracks::init() {
@@ -70,13 +84,20 @@ void FilterTracks::processEvent(LCEvent* evt) {
   LCCollectionVec* OutputTrackCollection = new LCCollectionVec(LCIO::TRACK);
   OutputTrackCollection->setSubset(true);
 
+  // track-mcpart relation output collections
+  UTIL::LCRelationNavigator myNav = UTIL::LCRelationNavigator(LCIO::TRACK, LCIO::MCPARTICLE); 
+
   // Get input collection
   LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection);
 
   if (InputTrackCollection->getTypeName() != lcio::LCIO::TRACK) {
     throw EVENT::Exception("Invalid collection type: " + InputTrackCollection->getTypeName());
   }
 
+  // Get input relations
+  LCCollection* InputTrackRelationCollection = evt->getCollection(_InputTrackRelationCollection);
+  std::shared_ptr<LCRelationNavigator> relation = std::make_shared<LCRelationNavigator>(InputTrackRelationCollection);
+
   // Filter
   std::string encoderString = lcio::LCTrackerCellID::encoding_string();
   UTIL::CellIDDecoder<lcio::TrackerHit> decoder(encoderString);
@@ -122,14 +143,24 @@ void FilterTracks::processEvent(LCEvent* evt) {
       }
     } else { // track property cuts
       if (nhittotal > _NHitsTotal && nhitvertex > _NHitsVertex && nhitinner > _NHitsInner && nhitouter > _NHitsOuter &&
-          pt > _MinPt && chi2spatial > _Chi2Spatial && nholes <= _MaxHoles) {
+          pt > _MinPt && chi2spatial > _Chi2Spatial && nholes <= _MaxHoles &&
+          fabs(trk->getD0()) < _MaxD0 && fabs(trk->getZ0()) < _MaxZ0) {
         OutputTrackCollection->addElement(trk);
+
+        auto mcParticleVec = relation->getRelatedToObjects(trk);
+        auto weightVec = relation->getRelatedToWeights(trk);
+        for (size_t irel = 0; irel < mcParticleVec.size(); ++irel) {
+          myNav.addRelation(trk, mcParticleVec[irel], weightVec[irel]);
+        }
+
       }
     }
   }
 
   // Save output track collection
   evt->addCollection(OutputTrackCollection, _OutputTrackCollection);
+  LCCollection* outputTrackHitRel = myNav.createLCCollection();
+  evt->addCollection(outputTrackHitRel, _OutputTrackRelationCollection);
 }
 
 void FilterTracks::end() {}