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
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
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() )
00147 {
00148
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
00158
00159
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
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
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
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
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
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
00271
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" );