MultiTrajectoryStateAssembler.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTools/interface/MultiTrajectoryStateAssembler.h"
#include "TrackingTools/GsfTools/interface/GetComponents.h"
#include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"
#include "TrackingTools/GsfTools/src/TrajectoryStateLessWeight.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () :
  combinationDone(false),
  thePzError(false),
  theValidWeightSum(0.),
  theInvalidWeightSum(0.)
{
  //
  // parameters (could be configurable)
  //
  sortStates = false;
  minValidFraction = 0.01;
  minFractionalWeight = 1.e-6; // 4;
}  

void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) {
  //
  // refuse to add states after combination has been done
  //
  if ( combinationDone )
    throw cms::Exception("LogicError") 
      << "MultiTrajectoryStateAssembler: trying to add states after combination";
  //
  // Verify validity of state to be added
  //
  if ( !tsos.isValid() )
    throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add invalid state";
  //
  // Add components (i.e. state to be added can be single or multi state)
  //
  GetComponents comps(tsos);
  MultiTSOS components(comps());
  addStateVector(components);
}

void MultiTrajectoryStateAssembler::addStateVector (const MultiTSOS& states)
{
  //
  // refuse to add states after combination has been done
  //
  if ( combinationDone )
    throw cms::Exception("LogicError") 
      << "MultiTrajectoryStateAssembler: trying to add states after combination";
  //
  // sum up weights (all components are supposed to be valid!!!) and
  // check for consistent pz
  //
  double sum(0.);
  double pzFirst = theStates.empty() ? 0. : theStates.front().localParameters().pzSign();
  for ( MultiTSOS::const_iterator i=states.begin();
    i!=states.end(); i++ ) {
    if ( !(i->isValid()) )
      throw cms::Exception("LogicError") 
    << "MultiTrajectoryStateAssembler: trying to add invalid state";
    // weights
    sum += i->weight();
    // check on p_z
    if ( !theStates.empty() && 
     pzFirst*i->localParameters().pzSign()<0. )  thePzError = true;
  }
  theValidWeightSum += sum;
  //
  // add to vector of states
  //
  theStates.insert(theStates.end(),states.begin(),states.end());
}


void MultiTrajectoryStateAssembler::addInvalidState (const double weight) {
  //
  // change status of combination (contains at least one invalid state)
  //
  theInvalidWeightSum += weight;
}

TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () {
  //
  // Prepare resulting state vector
  //
  if ( !prepareCombinedState() )  return TSOS();
  //
  // If invalid states in input: use reweighting
  //
  if ( theInvalidWeightSum>0. )  
    return reweightedCombinedState(theValidWeightSum+theInvalidWeightSum);
  //
  // Return new multi state without reweighting
  //
  return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(theStates)));
}

TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState (const float newWeight) {
  //
  // Prepare resulting state vector
  //
  if ( !prepareCombinedState() )  return TSOS();
  //
  // return reweighted state
  //
  return reweightedCombinedState(newWeight);
}

bool
MultiTrajectoryStateAssembler::prepareCombinedState () {
  //
  // Protect against empty combination (no valid input state)
  //
  if ( invalidCombinedState() )  return false;
  //
  // Check for states with wrong pz
  //
  if ( thePzError )  removeWrongPz();
  //
  // Check for minimum fraction of valid states
  //
  double allWeights(theValidWeightSum+theInvalidWeightSum);
  if ( theInvalidWeightSum>0. && theValidWeightSum <minValidFraction*allWeights )  return false;
  //
  // remaining part to be done only once
  //
  if ( combinationDone )  return true;
  combinationDone = true;
  //
  // Remove states with negligible weights
  //
  removeSmallWeights();
  if ( invalidCombinedState() )  return false;
  //
  // Sort output by weights?
  //
  if ( sortStates ) 
    sort(theStates.begin(),theStates.end(),TrajectoryStateLessWeight());

  return true;
}

TrajectoryStateOnSurface
MultiTrajectoryStateAssembler::reweightedCombinedState (const double newWeight) const {
  //
  // check status
  //
  if ( invalidCombinedState() )  return TSOS();
  //
  // scaling factor
  //
  double factor = theValidWeightSum>0. ? newWeight/theValidWeightSum : 1;
  //
  // create new vector of states & combined state
  //
  MultiTSOS reweightedStates;
  reweightedStates.reserve(theStates.size());
  for (auto const & is : theStates) {
    auto oldWeight = is.weight();
    reweightedStates.emplace_back(factor*oldWeight,
                                                        is.localParameters(),
                            is.localError(),
                            is.surface(),
                            &(is.globalParameters().magneticField()),
                            is.surfaceSide()
                              );
  }
  return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(reweightedStates)));
}

void
MultiTrajectoryStateAssembler::removeSmallWeights()
{
  //
  // check total weight
  //
  auto totalWeight(theInvalidWeightSum+theValidWeightSum);
  if ( totalWeight == 0. ) {
    theStates.clear();
    return;
  }
  theStates.erase(std::remove_if(theStates.begin(),theStates.end(),
                  [&](MultiTSOS::value_type const & s){ return s.weight() < minFractionalWeight*totalWeight;}),
                 theStates.end());
}

void
MultiTrajectoryStateAssembler::removeWrongPz () {
   LogDebug("GsfTrackFitters")
       << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n"
       << "  #state / weights before cleaning = " << theStates.size()
       << " / " << theValidWeightSum
       << " / " << theInvalidWeightSum;
  //
  // Calculate average pz
  //
  double meanPz(0.);
  for (auto const & is : theStates)
    meanPz += is.weight()*is.localParameters().pzSign();
  meanPz /= theValidWeightSum;
  //
  // Now keep only states compatible with the average pz
  //
  theValidWeightSum = 0.;
  MultiTSOS oldStates(theStates);
  theStates.clear();
  for (auto const & is :oldStates) {
    if ( meanPz*is.localParameters().pzSign()>=0. ) {
      theValidWeightSum += is.weight();
      theStates.push_back(std::move(is));
    }
    else {
      theInvalidWeightSum += is.weight();
        LogDebug("GsfTrackFitters")
          << "removing  weight / pz / global position = " << is.weight()
          << " " << is.localParameters().pzSign()
          << " " << is.globalPosition();

    }
  }
  LogDebug("GsfTrackFitters") 
       << "  #state / weights after cleaning = " << theStates.size()
       << " / " << theValidWeightSum
       << " / " << theInvalidWeightSum;
}