TwoBodyDecayTrajectoryFactory.cc

Go to the documentation of this file.

#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
#include "DataFormats/Math/interface/Vector.h"
#include "DataFormats/Math/interface/Error.h"
#include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h"
#include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
 
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
 
#include "Geometry/CommonDetUnit/interface/TrackerGeomDet.h"
 
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
 
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
#include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectory.h"
 
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayFitter.h"
 
class TwoBodyDecayTrajectoryFactory : public TrajectoryFactoryBase {
public:
  typedef TwoBodyDecayVirtualMeasurement VirtualMeasurement;
  typedef TwoBodyDecayTrajectoryState::TsosContainer TsosContainer;
  typedef TwoBodyDecayTrajectory::ConstRecHitCollection ConstRecHitCollection;
 
  TwoBodyDecayTrajectoryFactory(const edm::ParameterSet &config);
  ~TwoBodyDecayTrajectoryFactory() override;
 
  const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
                                                   const ConstTrajTrackPairCollection &tracks,
                                                   const reco::BeamSpot &beamSpot) const override;
 
  const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
                                                   const ConstTrajTrackPairCollection &tracks,
                                                   const ExternalPredictionCollection &external,
                                                   const reco::BeamSpot &beamSpot) const override;
 
  TwoBodyDecayTrajectoryFactory *clone() const override { return new TwoBodyDecayTrajectoryFactory(*this); }
 
protected:
  const ReferenceTrajectoryCollection constructTrajectories(const ConstTrajTrackPairCollection &tracks,
                                                            const TwoBodyDecay &tbd,
                                                            const MagneticField *magField,
                                                            const reco::BeamSpot &beamSpot,
                                                            bool setParameterErrors) const;
 
  bool match(const TrajectoryStateOnSurface &state, const TransientTrackingRecHit::ConstRecHitPointer &recHit) const;
 
  TwoBodyDecayFitter theFitter;
 
  double thePrimaryMass;
  double thePrimaryWidth;
  double theSecondaryMass;
 
  double theNSigmaCutValue;
  double theChi2CutValue;
  bool theUseRefittedStateFlag;
  bool theConstructTsosWithErrorsFlag;
};
 
 
TwoBodyDecayTrajectoryFactory::TwoBodyDecayTrajectoryFactory(const edm::ParameterSet &config)
    : TrajectoryFactoryBase(config, 2), theFitter(config) {
  const edm::ParameterSet ppc = config.getParameter<edm::ParameterSet>("ParticleProperties");
  thePrimaryMass = ppc.getParameter<double>("PrimaryMass");
  thePrimaryWidth = ppc.getParameter<double>("PrimaryWidth");
  theSecondaryMass = ppc.getParameter<double>("SecondaryMass");
 
  theNSigmaCutValue = config.getParameter<double>("NSigmaCut");
  theChi2CutValue = config.getParameter<double>("Chi2Cut");
  theUseRefittedStateFlag = config.getParameter<bool>("UseRefittedState");
  theConstructTsosWithErrorsFlag = config.getParameter<bool>("ConstructTsosWithErrors");
}
 
TwoBodyDecayTrajectoryFactory::~TwoBodyDecayTrajectoryFactory() {}
 
const TrajectoryFactoryBase::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::trajectories(
    const edm::EventSetup &setup, const ConstTrajTrackPairCollection &tracks, const reco::BeamSpot &beamSpot) const {
  ReferenceTrajectoryCollection trajectories;
 
  edm::ESHandle<MagneticField> magneticField;
  setup.get<IdealMagneticFieldRecord>().get(magneticField);
 
  if (tracks.size() == 2) {
    // produce transient tracks from persistent tracks
    std::vector<reco::TransientTrack> transientTracks(2);
 
    transientTracks[0] = reco::TransientTrack(*tracks[0].second, magneticField.product());
    transientTracks[0].setES(setup);
 
    transientTracks[1] = reco::TransientTrack(*tracks[1].second, magneticField.product());
    transientTracks[1].setES(setup);
 
    // estimate the decay parameters
    VirtualMeasurement vm(thePrimaryMass, thePrimaryWidth, theSecondaryMass, beamSpot);
    TwoBodyDecay tbd = theFitter.estimate(transientTracks, vm);
 
    if (!tbd.isValid() || (tbd.chi2() > theChi2CutValue)) {
      trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
      return trajectories;
    }
 
    return constructTrajectories(tracks, tbd, magneticField.product(), beamSpot, false);
  } else {
    edm::LogInfo("ReferenceTrajectories") << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
                                          << "Need 2 tracks, got " << tracks.size() << ".\n";
  }
 
  return trajectories;
}
 
const TrajectoryFactoryBase::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::trajectories(
    const edm::EventSetup &setup,
    const ConstTrajTrackPairCollection &tracks,
    const ExternalPredictionCollection &external,
    const reco::BeamSpot &beamSpot) const {
  ReferenceTrajectoryCollection trajectories;
 
  edm::ESHandle<MagneticField> magneticField;
  setup.get<IdealMagneticFieldRecord>().get(magneticField);
 
  if (tracks.size() == 2 && external.size() == 2) {
    if (external[0].isValid() && external[1].isValid())  // Include external estimates
    {
      // produce transient tracks from persistent tracks
      std::vector<reco::TransientTrack> transientTracks(2);
 
      transientTracks[0] = reco::TransientTrack(*tracks[0].second, magneticField.product());
      transientTracks[0].setES(setup);
 
      transientTracks[1] = reco::TransientTrack(*tracks[1].second, magneticField.product());
      transientTracks[1].setES(setup);
 
      // estimate the decay parameters. the transient tracks are not really associated to the
      // the external tsos, but this is o.k., because the only information retrieved from them
      // is the magnetic field.
      VirtualMeasurement vm(thePrimaryMass, thePrimaryWidth, theSecondaryMass, beamSpot);
      TwoBodyDecay tbd = theFitter.estimate(transientTracks, external, vm);
 
      if (!tbd.isValid() || (tbd.chi2() > theChi2CutValue)) {
        trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
        return trajectories;
      }
 
      return constructTrajectories(tracks, tbd, magneticField.product(), beamSpot, true);
    } else {
      // Return without external estimate
      trajectories = this->trajectories(setup, tracks, beamSpot);
    }
  } else {
    edm::LogInfo("ReferenceTrajectories") << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
                                          << "Need 2 tracks, got " << tracks.size() << ".\n";
  }
 
  return trajectories;
}
 
const TwoBodyDecayTrajectoryFactory::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::constructTrajectories(
    const ConstTrajTrackPairCollection &tracks,
    const TwoBodyDecay &tbd,
    const MagneticField *magField,
    const reco::BeamSpot &beamSpot,
    bool setParameterErrors) const {
  ReferenceTrajectoryCollection trajectories;
 
  // get innermost valid trajectory state and hits from the tracks
  TrajectoryInput input1 = this->innermostStateAndRecHits(tracks[0]);
  TrajectoryInput input2 = this->innermostStateAndRecHits(tracks[1]);
 
  if (!(input1.first.isValid() && input2.first.isValid()))
    return trajectories;
 
  // produce TwoBodyDecayTrajectoryState (input for TwoBodyDecayTrajectory)
  TsosContainer tsos(input1.first, input2.first);
  ConstRecHitCollection recHits(input1.second, input2.second);
  TwoBodyDecayTrajectoryState trajectoryState(tsos, tbd, theSecondaryMass, magField);
 
  if (!trajectoryState.isValid()) {
    trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
    return trajectories;
  }
 
  // always use the refitted trajectory state for matching
  // FIXME FIXME CLONE
  //TrackingRecHit::ConstRecHitPointer updatedRecHit1( recHits.first.front()->clone( tsos.first ) );
  // TrackingRecHit::ConstRecHitPointer updatedRecHit2( recHits.second.front()->clone( tsos.second ) );
 
  TrackingRecHit::ConstRecHitPointer updatedRecHit1(recHits.first.front());
  TrackingRecHit::ConstRecHitPointer updatedRecHit2(recHits.second.front());
 
  bool valid1 = match(trajectoryState.trajectoryStates(true).first, updatedRecHit1);
 
  bool valid2 = match(trajectoryState.trajectoryStates(true).second, updatedRecHit2);
 
  if (!valid1 || !valid2) {
    trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
    return trajectories;
  }
 
  ReferenceTrajectoryBase::Config config(materialEffects(), propagationDirection());
  config.useBeamSpot = useBeamSpot_;
  config.includeAPEs = includeAPEs_;
  config.allowZeroMaterial = allowZeroMaterial_;
  config.useRefittedState = theUseRefittedStateFlag;
  config.constructTsosWithErrors = theConstructTsosWithErrorsFlag;
  // set the flag for reversing the RecHits to false, since they are already in the correct order.
  config.hitsAreReverse = false;
  TwoBodyDecayTrajectory *result = new TwoBodyDecayTrajectory(trajectoryState, recHits, magField, beamSpot, config);
  if (setParameterErrors && tbd.hasError())
    result->setParameterErrors(tbd.covariance());
  trajectories.push_back(ReferenceTrajectoryPtr(result));
  return trajectories;
}
 
bool TwoBodyDecayTrajectoryFactory::match(const TrajectoryStateOnSurface &state,
                                          const TransientTrackingRecHit::ConstRecHitPointer &recHit) const {
  LocalPoint lp1 = state.localPosition();
  LocalPoint lp2 = recHit->localPosition();
 
  double deltaX = lp1.x() - lp2.x();
  double deltaY = lp1.y() - lp2.y();
 
  LocalError le = recHit->localPositionError();
 
  double varX = le.xx();
  double varY = le.yy();
 
  return ((fabs(deltaX) / sqrt(varX) < theNSigmaCutValue) && (fabs(deltaY) / sqrt(varY) < theNSigmaCutValue));
}
 
DEFINE_EDM_PLUGIN(TrajectoryFactoryPlugin, TwoBodyDecayTrajectoryFactory, "TwoBodyDecayTrajectoryFactory");