GsfTools/src/MultiTrajectoryStateAssembler.cc

0001 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateAssembler.h"
0002 #include "TrackingTools/GsfTools/interface/GetComponents.h"
0003 #include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"
0004 #include "TrackingTools/GsfTools/src/TrajectoryStateLessWeight.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler()
0009     : combinationDone(false), thePzError(false), theValidWeightSum(0.), theInvalidWeightSum(0.) {
0010   //
0011   // parameters (could be configurable)
0012   //
0013   sortStates = false;
0014   minValidFraction = 0.01;
0015   minFractionalWeight = 1.e-6;  // 4;
0016 }
0017
0018 void MultiTrajectoryStateAssembler::addState(const TrajectoryStateOnSurface tsos) {
0019   //
0020   // refuse to add states after combination has been done
0021   //
0022   if (combinationDone)
0023     throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add states after combination";
0024   //
0025   // Verify validity of state to be added
0026   //
0027   if (!tsos.isValid())
0028     throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add invalid state";
0029   //
0030   // Add components (i.e. state to be added can be single or multi state)
0031   //
0032   GetComponents comps(tsos);
0033   const MultiTSOS &components(comps());
0034   addStateVector(components);
0035 }
0036
0037 void MultiTrajectoryStateAssembler::addStateVector(const MultiTSOS &states) {
0038   //
0039   // refuse to add states after combination has been done
0040   //
0041   if (combinationDone)
0042     throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add states after combination";
0043   //
0044   // sum up weights (all components are supposed to be valid!!!) and
0045   // check for consistent pz
0046   //
0047   double sum(0.);
0048   double pzFirst = theStates.empty() ? 0. : theStates.front().localParameters().pzSign();
0049   for (MultiTSOS::const_iterator i = states.begin(); i != states.end(); i++) {
0050     if (!(i->isValid()))
0051       throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add invalid state";
0052     // weights
0053     sum += i->weight();
0054     // check on p_z
0055     if (!theStates.empty() && pzFirst * i->localParameters().pzSign() < 0.)
0056       thePzError = true;
0057   }
0058   theValidWeightSum += sum;
0059   //
0060   // add to vector of states
0061   //
0062   theStates.insert(theStates.end(), states.begin(), states.end());
0063 }
0064
0065 void MultiTrajectoryStateAssembler::addInvalidState(const double weight) {
0066   //
0067   // change status of combination (contains at least one invalid state)
0068   //
0069   theInvalidWeightSum += weight;
0070 }
0071
0072 TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState() {
0073   //
0074   // Prepare resulting state vector
0075   //
0076   if (!prepareCombinedState())
0077     return TSOS();
0078   //
0079   // If invalid states in input: use reweighting
0080   //
0081   if (theInvalidWeightSum > 0.)
0082     return reweightedCombinedState(theValidWeightSum + theInvalidWeightSum);
0083   //
0084   // Return new multi state without reweighting
0085   //
0086   return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(theStates)));
0087 }
0088
0089 TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState(const float newWeight) {
0090   //
0091   // Prepare resulting state vector
0092   //
0093   if (!prepareCombinedState())
0094     return TSOS();
0095   //
0096   // return reweighted state
0097   //
0098   return reweightedCombinedState(newWeight);
0099 }
0100
0101 bool MultiTrajectoryStateAssembler::prepareCombinedState() {
0102   //
0103   // Protect against empty combination (no valid input state)
0104   //
0105   if (invalidCombinedState())
0106     return false;
0107   //
0108   // Check for states with wrong pz
0109   //
0110   if (thePzError)
0111     removeWrongPz();
0112   //
0113   // Check for minimum fraction of valid states
0114   //
0115   double allWeights(theValidWeightSum + theInvalidWeightSum);
0116   if (theInvalidWeightSum > 0. && theValidWeightSum < minValidFraction * allWeights)
0117     return false;
0118   //
0119   // remaining part to be done only once
0120   //
0121   if (combinationDone)
0122     return true;
0123   combinationDone = true;
0124   //
0125   // Remove states with negligible weights
0126   //
0127   removeSmallWeights();
0128   if (invalidCombinedState())
0129     return false;
0130   //
0131   // Sort output by weights?
0132   //
0133   if (sortStates)
0134     sort(theStates.begin(), theStates.end(), TrajectoryStateLessWeight());
0135
0136   return true;
0137 }
0138
0139 TrajectoryStateOnSurface MultiTrajectoryStateAssembler::reweightedCombinedState(const double newWeight) const {
0140   //
0141   // check status
0142   //
0143   if (invalidCombinedState())
0144     return TSOS();
0145   //
0146   // scaling factor
0147   //
0148   double factor = theValidWeightSum > 0. ? newWeight / theValidWeightSum : 1;
0149   //
0150   // create new vector of states & combined state
0151   //
0152   MultiTSOS reweightedStates;
0153   reweightedStates.reserve(theStates.size());
0154   for (auto const &is : theStates) {
0155     auto oldWeight = is.weight();
0156     reweightedStates.emplace_back(factor * oldWeight,
0157                                   is.localParameters(),
0158                                   is.localError(),
0159                                   is.surface(),
0160                                   &(is.globalParameters().magneticField()),
0161                                   is.surfaceSide());
0162   }
0163   return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(reweightedStates)));
0164 }
0165
0166 void MultiTrajectoryStateAssembler::removeSmallWeights() {
0167   //
0168   // check total weight
0169   //
0170   auto totalWeight(theInvalidWeightSum + theValidWeightSum);
0171   if (totalWeight == 0.) {
0172     theStates.clear();
0173     return;
0174   }
0175   theStates.erase(
0176       std::remove_if(theStates.begin(),
0177                      theStates.end(),
0178                      [&](MultiTSOS::value_type const &s) { return s.weight() < minFractionalWeight * totalWeight; }),
0179       theStates.end());
0180 }
0181
0182 void MultiTrajectoryStateAssembler::removeWrongPz() {
0183   LogDebug("GsfTrackFitters") << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n"
0184                               << "  #state / weights before cleaning = " << theStates.size() << " / "
0185                               << theValidWeightSum << " / " << theInvalidWeightSum;
0186   //
0187   // Calculate average pz
0188   //
0189   double meanPz(0.);
0190   for (auto const &is : theStates)
0191     meanPz += is.weight() * is.localParameters().pzSign();
0192   meanPz /= theValidWeightSum;
0193
0194   // When meanPz=0 it is not possible to check which are the states non agreeing with the average pZ sign
0195   if (meanPz == 0.) {
0196     edm::LogError("MultiTrajectoryStateAssembler")
0197         << " input multistate has average pZ sign == 0. Rejecting!" << std::endl;
0198     theStates.clear();
0199     return;
0200   }
0201
0202   //
0203   // Now keep only states compatible with the average pz
0204   //
0205   theValidWeightSum = 0.;
0206   MultiTSOS oldStates(theStates);
0207   theStates.clear();
0208   for (auto const &is : oldStates) {
0209     if (meanPz * is.localParameters().pzSign() > 0.) {
0210       theValidWeightSum += is.weight();
0211       theStates.push_back(is);
0212     } else {
0213       theInvalidWeightSum += is.weight();
0214       LogDebug("GsfTrackFitters") << "removing  weight / pz / global position = " << is.weight() << " "
0215                                   << is.localParameters().pzSign() << " " << is.globalPosition();
0216     }
0217   }
0218   LogDebug("GsfTrackFitters") << "  #state / weights after cleaning = " << theStates.size() << " / "
0219                               << theValidWeightSum << " / " << theInvalidWeightSum;
0220 }