CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc

Go to the documentation of this file.
00001 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateAssembler.h"
00002 
00003 #include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"
00004 #include "TrackingTools/GsfTools/src/TrajectoryStateLessWeight.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "FWCore/Utilities/interface/Exception.h"
00007 
00008 MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () :
00009   combinationDone(false),
00010   thePzError(false),
00011   theValidWeightSum(0.),
00012   theInvalidWeightSum(0.)
00013 {
00014   //
00015   // parameters (could be configurable)
00016   //
00017   sortStates = false;
00018   minValidFraction = 0.01;
00019   minFractionalWeight = 1.e-6;
00020   //   //
00021   //   // Timers
00022   //   //
00023   //   if ( theTimerAdd==0 ) {
00024   //     theTimerAdd = 
00025   //       &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::addState")]; 
00026   //     theTimerAdd->switchCPU(false);
00027   //     theTimerComb = 
00028   //       &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::combinedState")]; 
00029   //     theTimerComb->switchCPU(false);
00030   //   }
00031 }  
00032 
00033 void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) {
00034   //   // Timer
00035   //   TimeMe t(*theTimerAdd,false);
00036   //
00037   // refuse to add states after combination has been done
00038   //
00039   if ( combinationDone )
00040     throw cms::Exception("LogicError") 
00041       << "MultiTrajectoryStateAssembler: trying to add states after combination";
00042   //
00043   // Verify validity of state to be added
00044   //
00045   if ( !tsos.isValid() )
00046     throw cms::Exception("LogicError") << "MultiTrajectoryStateAssembler: trying to add invalid state";
00047   //
00048   // Add components (i.e. state to be added can be single or multi state)
00049   //
00050   MultiTSOS components(tsos.components());
00051   addStateVector(components);
00052 }
00053 
00054 void MultiTrajectoryStateAssembler::addStateVector (const MultiTSOS& states)
00055 {
00056   //
00057   // refuse to add states after combination has been done
00058   //
00059   if ( combinationDone )
00060     throw cms::Exception("LogicError") 
00061       << "MultiTrajectoryStateAssembler: trying to add states after combination";
00062   //
00063   // sum up weights (all components are supposed to be valid!!!) and
00064   // check for consistent pz
00065   //
00066   double sum(0.);
00067   double pzFirst = theStates.empty() ? 0. : theStates.front().localParameters().pzSign();
00068   for ( MultiTSOS::const_iterator i=states.begin();
00069         i!=states.end(); i++ ) {
00070     if ( !(i->isValid()) )
00071       throw cms::Exception("LogicError") 
00072         << "MultiTrajectoryStateAssembler: trying to add invalid state";
00073     // weights
00074     sum += i->weight();
00075     // check on p_z
00076     if ( !theStates.empty() && 
00077          pzFirst*i->localParameters().pzSign()<0. )  thePzError = true;
00078   }
00079   theValidWeightSum += sum;
00080   //
00081   // add to vector of states
00082   //
00083   theStates.insert(theStates.end(),states.begin(),states.end());
00084 }
00085 
00086 
00087 void MultiTrajectoryStateAssembler::addInvalidState (const double weight) {
00088   //
00089   // change status of combination (contains at least one invalid state)
00090   //
00091   theInvalidWeightSum += weight;
00092 }
00093 
00094 TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () {
00095   //   // Timer
00096   //   TimeMe t(*theTimerComb,false);
00097   //
00098   // Prepare resulting state vector
00099   //
00100   if ( !prepareCombinedState() )  return TSOS();
00101   //
00102   // If invalid states in input: use reweighting
00103   //
00104   if ( theInvalidWeightSum>0. )  
00105     return reweightedCombinedState(theValidWeightSum+theInvalidWeightSum);
00106   //
00107   // Return new multi state without reweighting
00108   //
00109   return TSOS(new BasicMultiTrajectoryState(theStates));
00110 }
00111 
00112 TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState (const float newWeight) {
00113   //   // Timer
00114   //   TimeMe t(*theTimerComb,false);
00115   //
00116   // Prepare resulting state vector
00117   //
00118   if ( !prepareCombinedState() )  return TSOS();
00119   //
00120   // return reweighted state
00121   //
00122   return reweightedCombinedState(newWeight);
00123 }
00124 
00125 bool
00126 MultiTrajectoryStateAssembler::prepareCombinedState () {
00127   //
00128   // Protect against empty combination (no valid input state)
00129   //
00130   if ( invalidCombinedState() )  return false;
00131   //
00132   // Check for states with wrong pz
00133   //
00134   if ( thePzError )  removeWrongPz();
00135   //
00136   // Check for minimum fraction of valid states
00137   //
00138   double allWeights(theValidWeightSum+theInvalidWeightSum);
00139   if ( theInvalidWeightSum>0. && (theValidWeightSum/allWeights)<minValidFraction )  return false;
00140   //
00141   // remaining part to be done only once
00142   //
00143   if ( combinationDone )  return true;
00144   else  combinationDone = true;
00145   //
00146   // Remove states with negligible weights
00147   //
00148   removeSmallWeights();
00149   if ( invalidCombinedState() )  return false;
00150   //
00151   // Sort output by weights?
00152   //
00153   if ( sortStates ) 
00154     sort(theStates.begin(),theStates.end(),TrajectoryStateLessWeight());
00155 
00156   return true;
00157 }
00158 
00159 TrajectoryStateOnSurface
00160 MultiTrajectoryStateAssembler::reweightedCombinedState (const double newWeight) const {
00161   //
00162   // check status
00163   //
00164   if ( invalidCombinedState() )  return TSOS();
00165   //
00166   // scaling factor
00167   //
00168   double factor = theValidWeightSum>0. ? newWeight/theValidWeightSum : 1;
00169   //
00170   // create new vector of states & combined state
00171   //
00172   MultiTSOS reweightedStates;
00173   reweightedStates.reserve(theStates.size());
00174   for ( MultiTSOS::const_iterator i=theStates.begin();
00175         i!=theStates.end(); i++ ) {
00176     double oldWeight = i->weight();
00177     reweightedStates.push_back(TrajectoryStateOnSurface(i->localParameters(),
00178                                                         i->localError(),
00179                                                         i->surface(),
00180                                                         &(i->globalParameters().magneticField()),
00181                                                         i->surfaceSide(),
00182                                                         factor*oldWeight));
00183   }
00184   return TSOS(new BasicMultiTrajectoryState(reweightedStates));
00185 }
00186 
00187 void
00188 MultiTrajectoryStateAssembler::removeSmallWeights()
00189 {
00190   //
00191   // check total weight
00192   //
00193   double totalWeight(theInvalidWeightSum+theValidWeightSum);
00194   if ( totalWeight == 0. ) {
00195     theStates.clear();
00196     return;
00197   }
00198   //
00199   // Loop until no more states are removed
00200   //
00201   bool redo;
00202   do {
00203     redo = false;
00204     for ( MultiTSOS::iterator i=theStates.begin();
00205           i!=theStates.end(); i++ ) {
00206       if ( (*i).weight()/totalWeight < minFractionalWeight ) {
00207         theStates.erase(i);
00208         redo = true;
00209         break;
00210       }
00211     }
00212   } while (redo);
00213 }
00214 
00215 void
00216 MultiTrajectoryStateAssembler::removeWrongPz () {
00217   //   edm::LogDebug("MultiTrajectoryStateAssembler") 
00218   //     << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n"
00219   //     << "  #state / weights before cleaning = " << theStates.size()
00220   //     << " / " << theValidWeightSum
00221   //     << " / " << theInvalidWeightSum;
00222   //
00223   // Calculate average pz
00224   //
00225   double meanPz(0.);
00226   for ( MultiTSOS::const_iterator is=theStates.begin();
00227         is!=theStates.end(); is++ ) {
00228     meanPz += is->weight()*is->localParameters().pzSign();
00229     //     edm::LogDebug("MultiTrajectoryStateAssembler") 
00230     //       << "  weight / pz / global position = " << is->weight() 
00231     //       << " " << is->localParameters().pzSign() 
00232     //       << " " << is->globalPosition();
00233   }
00234   meanPz /= theValidWeightSum;
00235   //
00236   // Now keep only states compatible with the average pz
00237   //
00238   //   double oldValidWeight(theValidWeightSum);
00239   theValidWeightSum = 0.;
00240   MultiTSOS oldStates(theStates);
00241   theStates.clear();
00242   for ( MultiTSOS::const_iterator is=oldStates.begin();
00243         is!=oldStates.end(); is++ ) {
00244     if ( meanPz*is->localParameters().pzSign()>=0. ) {
00245       theValidWeightSum += is->weight();
00246       theStates.push_back(*is);
00247     }
00248     else {
00249       theInvalidWeightSum += is->weight();
00250     }
00251   }
00252   //   edm::LogDebug("MultiTrajectoryStateAssembler") 
00253   //     << "  #state / weights after cleaning = " << theStates.size()
00254   //     << " / " << theValidWeightSum
00255   //     << " / " << theInvalidWeightSum;
00256 }
00257 
00258 // TimingReport::Item * MultiTrajectoryStateAssembler::theTimerAdd(0);
00259 // TimingReport::Item * MultiTrajectoryStateAssembler::theTimerComb(0);