CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Alignment/ReferenceTrajectories/plugins/TwoBodyDecayTrajectoryFactory.cc

Go to the documentation of this file.
00001 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" 
00002 #include "DataFormats/GeometrySurface/interface/Surface.h" 
00003 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h" 
00004 #include "DataFormats/Math/interface/Vector.h" 
00005 #include "DataFormats/Math/interface/Error.h" 
00006 #include "TrackingTools/TrajectoryState/interface/CopyUsingClone.h" 
00007 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h" 
00008 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
00009 
00010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00011 
00012 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
00013 #include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
00014 
00015 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 
00016 
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00019 
00020 
00021 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
00022 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectory.h"
00023 
00024 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayFitter.h"
00025 
00031 class TwoBodyDecayTrajectoryFactory : public TrajectoryFactoryBase
00032 {
00033 public:
00034 
00035   typedef TwoBodyDecayVirtualMeasurement VirtualMeasurement;
00036   typedef TwoBodyDecayTrajectoryState::TsosContainer TsosContainer;
00037   typedef TwoBodyDecayTrajectory::ConstRecHitCollection ConstRecHitCollection;
00038 
00039   TwoBodyDecayTrajectoryFactory(const edm::ParameterSet &config);
00040   ~TwoBodyDecayTrajectoryFactory() {}
00041 
00043   virtual const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
00044                                                            const ConstTrajTrackPairCollection &tracks,
00045                                                            const reco::BeamSpot &beamSpot) const;
00046 
00047   virtual const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
00048                                                            const ConstTrajTrackPairCollection &tracks,
00049                                                            const ExternalPredictionCollection &external,
00050                                                            const reco::BeamSpot &beamSpot) const;
00051 
00052   virtual TwoBodyDecayTrajectoryFactory* clone() const { return new TwoBodyDecayTrajectoryFactory(*this); }
00053 
00054 protected:
00055   const ReferenceTrajectoryCollection constructTrajectories(const ConstTrajTrackPairCollection &tracks,
00056                                                             const TwoBodyDecay &tbd,
00057                                                             const MagneticField *magField,
00058                                                             const reco::BeamSpot &beamSpot,
00059                                                             bool setParameterErrors) const;
00060 
00061   bool match(const TrajectoryStateOnSurface &state,
00062              const TransientTrackingRecHit::ConstRecHitPointer &recHit) const;
00063     
00064   void produceVirtualMeasurement(const edm::ParameterSet &config);
00065 
00066   VirtualMeasurement theVM;
00067 
00068   TwoBodyDecayFitter* theFitter;
00069 
00070   double theNSigmaCutValue;
00071   bool theUseRefittedStateFlag;
00072   bool theConstructTsosWithErrorsFlag;
00073 
00074 };
00075 
00079 
00080 TwoBodyDecayTrajectoryFactory::TwoBodyDecayTrajectoryFactory( const edm::ParameterSet& config )
00081   : TrajectoryFactoryBase( config ),
00082     theFitter( new TwoBodyDecayFitter( config ) ),
00083     theNSigmaCutValue( config.getParameter< double >( "NSigmaCut" ) ),
00084     theUseRefittedStateFlag( config.getParameter< bool >( "UseRefittedState" ) ),
00085     theConstructTsosWithErrorsFlag( config.getParameter< bool >( "ConstructTsosWithErrors" ) )
00086 {
00087   produceVirtualMeasurement( config );
00088 }
00089 
00090 
00091 const TrajectoryFactoryBase::ReferenceTrajectoryCollection
00092 TwoBodyDecayTrajectoryFactory::trajectories(const edm::EventSetup &setup,
00093                                             const ConstTrajTrackPairCollection &tracks,
00094                                             const reco::BeamSpot &beamSpot) const
00095 {
00096   ReferenceTrajectoryCollection trajectories;
00097 
00098   edm::ESHandle< MagneticField > magneticField;
00099   setup.get< IdealMagneticFieldRecord >().get( magneticField );
00100 
00101   if ( tracks.size() == 2 )
00102   {
00103     // produce transient tracks from persistent tracks
00104     std::vector< reco::TransientTrack > transientTracks( 2 );
00105 
00106     transientTracks[0] = reco::TransientTrack( *tracks[0].second, magneticField.product() );
00107     transientTracks[0].setES( setup );
00108 
00109     transientTracks[1] = reco::TransientTrack( *tracks[1].second, magneticField.product() );
00110     transientTracks[1].setES( setup );
00111 
00112     // estimate the decay parameters
00113     TwoBodyDecay tbd = theFitter->estimate( transientTracks, theVM );
00114 
00115     if ( !tbd.isValid() )
00116     {
00117       trajectories.push_back( ReferenceTrajectoryPtr( new TwoBodyDecayTrajectory() ) );
00118       return trajectories;
00119     }
00120 
00121     return constructTrajectories( tracks, tbd, magneticField.product(), beamSpot, false );
00122   }
00123   else
00124   {
00125     edm::LogInfo( "ReferenceTrajectories" ) << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
00126                                             << "Need 2 tracks, got " << tracks.size() << ".\n";
00127   }
00128 
00129   return trajectories;
00130 }
00131 
00132 
00133 const TrajectoryFactoryBase::ReferenceTrajectoryCollection
00134 TwoBodyDecayTrajectoryFactory::trajectories(const edm::EventSetup &setup,
00135                                             const ConstTrajTrackPairCollection &tracks,
00136                                             const ExternalPredictionCollection &external,
00137                                             const reco::BeamSpot &beamSpot) const
00138 {
00139   ReferenceTrajectoryCollection trajectories;
00140 
00141   edm::ESHandle< MagneticField > magneticField;
00142   setup.get< IdealMagneticFieldRecord >().get( magneticField );
00143 
00144   if ( tracks.size() == 2 && external.size() == 2 )
00145   {
00146     if ( external[0].isValid() && external[1].isValid() ) // Include external estimates
00147     {
00148       // produce transient tracks from persistent tracks
00149       std::vector< reco::TransientTrack > transientTracks( 2 );
00150 
00151       transientTracks[0] = reco::TransientTrack( *tracks[0].second, magneticField.product() );
00152       transientTracks[0].setES( setup );
00153 
00154       transientTracks[1] = reco::TransientTrack( *tracks[1].second, magneticField.product() );
00155       transientTracks[1].setES( setup );
00156 
00157       // estimate the decay parameters. the transient tracks are not really associated to the
00158       // the external tsos, but this is o.k., because the only information retrieved from them
00159       // is the magnetic field.
00160       TwoBodyDecay tbd = theFitter->estimate( transientTracks, external, theVM );
00161 
00162       if ( !tbd.isValid() )
00163       {
00164         trajectories.push_back( ReferenceTrajectoryPtr( new TwoBodyDecayTrajectory() ) );
00165         return trajectories;
00166       }
00167 
00168       return constructTrajectories( tracks, tbd, magneticField.product(), beamSpot, true );
00169     }
00170     else
00171     {
00172       // Return without external estimate
00173       trajectories = this->trajectories(setup, tracks, beamSpot);
00174     }
00175   }
00176   else
00177   {
00178     edm::LogInfo( "ReferenceTrajectories" ) << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
00179                                             << "Need 2 tracks, got " << tracks.size() << ".\n";
00180   }
00181 
00182   return trajectories;
00183 }
00184 
00185 
00186 const TwoBodyDecayTrajectoryFactory::ReferenceTrajectoryCollection
00187 TwoBodyDecayTrajectoryFactory::constructTrajectories( const ConstTrajTrackPairCollection& tracks,
00188                                                       const TwoBodyDecay& tbd,
00189                                                       const MagneticField* magField,
00190                                                       const reco::BeamSpot &beamSpot,
00191                                                       bool setParameterErrors ) const
00192 {
00193   ReferenceTrajectoryCollection trajectories;
00194 
00195   // get innermost valid trajectory state and hits from the tracks
00196   TrajectoryInput input1 = this->innermostStateAndRecHits( tracks[0] );
00197   TrajectoryInput input2 = this->innermostStateAndRecHits( tracks[1] );
00198 
00199   if ( !( input1.first.isValid() && input2.first.isValid() ) ) return trajectories;
00200 
00201   // produce TwoBodyDecayTrajectoryState (input for TwoBodyDecayTrajectory)
00202   TsosContainer tsos( input1.first, input2.first );
00203   ConstRecHitCollection recHits( input1.second, input2.second );
00204   TwoBodyDecayTrajectoryState trajectoryState( tsos, tbd, theVM.secondaryMass(), magField );
00205 
00206   if ( !trajectoryState.isValid() )
00207   {
00208     trajectories.push_back( ReferenceTrajectoryPtr( new TwoBodyDecayTrajectory() ) );
00209     return trajectories;
00210   }
00211 
00212   // always use the refitted trajectory state for matching
00213   TransientTrackingRecHit::ConstRecHitPointer updatedRecHit1( recHits.first.front()->clone( tsos.first ) );
00214   bool valid1 = match( trajectoryState.trajectoryStates( true ).first,
00215                        updatedRecHit1 );
00216 
00217   TransientTrackingRecHit::ConstRecHitPointer updatedRecHit2( recHits.second.front()->clone( tsos.second ) );
00218   bool valid2 = match( trajectoryState.trajectoryStates( true ).second,
00219                        updatedRecHit2 );
00220 
00221   if ( !valid1 || !valid2 )
00222   {
00223     trajectories.push_back( ReferenceTrajectoryPtr( new TwoBodyDecayTrajectory() ) );
00224     return trajectories;
00225   }
00226 
00227   // set the flag for reversing the RecHits to false, since they are already in the correct order.
00228   TwoBodyDecayTrajectory* result = new TwoBodyDecayTrajectory( trajectoryState, recHits, magField,
00229                                                                materialEffects(), propagationDirection(),
00230                                                                false, beamSpot, theUseRefittedStateFlag,
00231                                                                theConstructTsosWithErrorsFlag );
00232   if ( setParameterErrors && tbd.hasError() ) result->setParameterErrors( tbd.covariance() );
00233   trajectories.push_back( ReferenceTrajectoryPtr( result ) );
00234   return trajectories;
00235 }
00236 
00237 
00238 bool TwoBodyDecayTrajectoryFactory::match( const TrajectoryStateOnSurface& state,
00239                                            const TransientTrackingRecHit::ConstRecHitPointer& recHit ) const
00240 {
00241   LocalPoint lp1 = state.localPosition();
00242   LocalPoint lp2 = recHit->localPosition();
00243 
00244   double deltaX = lp1.x() - lp2.x();
00245   double deltaY = lp1.y() - lp2.y();
00246 
00247   LocalError le = recHit->localPositionError();
00248 
00249   double varX = le.xx();
00250   double varY = le.yy();
00251 
00252   AlignmentPositionError* gape = recHit->det()->alignmentPositionError();
00253   if ( gape )
00254   {
00255     ErrorFrameTransformer eft;
00256     LocalError lape = eft.transform( gape->globalError(), recHit->det()->surface() );
00257 
00258     varX += lape.xx();
00259     varY += lape.yy();
00260   }
00261 
00262   return ( ( fabs(deltaX)/sqrt(varX) < theNSigmaCutValue ) && ( fabs(deltaY)/sqrt(varY) < theNSigmaCutValue ) );
00263 }
00264 
00265 
00266 void TwoBodyDecayTrajectoryFactory::produceVirtualMeasurement( const edm::ParameterSet& config )
00267 {
00268   const edm::ParameterSet bsc = config.getParameter< edm::ParameterSet >( "BeamSpot" );
00269   const edm::ParameterSet ppc = config.getParameter< edm::ParameterSet >( "ParticleProperties" );
00270   // FIXME: Should get 3D beamspot and errors from BeamSpot input from
00271   //        event with extrapolation of tracks to beam line?
00272 
00273   GlobalPoint theBeamSpot( bsc.getParameter< double >( "MeanX" ),
00274                            bsc.getParameter< double >( "MeanY" ),
00275                            bsc.getParameter< double >( "MeanZ" ) );
00276 
00277   GlobalError theBeamSpotError( bsc.getParameter< double >( "VarXX" ),
00278                                 bsc.getParameter< double >( "VarXY" ),
00279                                 bsc.getParameter< double >( "VarYY" ),
00280                                 bsc.getParameter< double >( "VarXZ" ),
00281                                 bsc.getParameter< double >( "VarYZ" ),
00282                                 bsc.getParameter< double >( "VarZZ" ) );
00283 
00284   theVM = VirtualMeasurement( ppc.getParameter< double >( "PrimaryMass" ),
00285                               ppc.getParameter< double >( "PrimaryWidth" ),
00286                               ppc.getParameter< double >( "SecondaryMass" ),
00287                               theBeamSpot, theBeamSpotError );
00288   return;
00289 }
00290 
00291 
00292 DEFINE_EDM_PLUGIN( TrajectoryFactoryPlugin, TwoBodyDecayTrajectoryFactory, "TwoBodyDecayTrajectoryFactory" );