CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/RecoMuon/src/MuonSeedTrack.cc

Go to the documentation of this file.
00001 
00013 #include "Validation/RecoMuon/src/MuonSeedTrack.h"
00014 
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "DataFormats/TrackReco/interface/TrackExtra.h"
00019 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00020 
00021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00022 
00023 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00025 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00026 #include "TrackingTools/GeomPropagators/interface/TrackerBounds.h"
00027 
00028 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00029 #include "RecoMuon/TrackingTools/interface/MuonUpdatorAtVertex.h"
00030 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00031 
00032 using namespace reco;
00033 using namespace edm;
00034 using namespace std;
00035 
00036 typedef TrajectoryStateOnSurface TSOS;
00037 
00038 
00039 //
00040 // constructors
00041 //
00042 MuonSeedTrack::MuonSeedTrack(const edm::ParameterSet& pset)
00043 {
00044   // service parameters
00045   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
00046   theService = new MuonServiceProxy(serviceParameters);
00047   
00048   ParameterSet updatorPar = pset.getParameter<ParameterSet>("MuonUpdatorAtVertexParameters");
00049   //theSeedPropagatorName = updatorPar.getParameter<string>("Propagator");
00050   theSeedsLabel = pset.getParameter<InputTag>("MuonSeed");
00051   theUpdatorAtVtx = new MuonUpdatorAtVertex(updatorPar,theService);
00052 
00053   theAllowNoVtxFlag = pset.getUntrackedParameter<bool>("AllowNoVertex",false);
00054 
00055   //register products
00056   setAlias(pset.getParameter<std::string>("@module_label"));
00057   produces<reco::TrackCollection>().setBranchAlias(theAlias + "Tracks");
00058 
00059 }
00060 
00061 //
00062 // destructor
00063 //
00064 MuonSeedTrack::~MuonSeedTrack()
00065 {
00066   if (theService) delete theService;
00067   if (theUpdatorAtVtx) delete theUpdatorAtVtx;
00068 }
00069 
00070 //
00071 // member functions
00072 //
00073 
00077 void
00078 MuonSeedTrack::produce(edm::Event& event, const edm::EventSetup& eventSetup)
00079 {
00080    using namespace edm;
00081 
00082   // Update the services
00083   theService->update(eventSetup);
00084 
00085   // the track collectios; they will be loaded in the event  
00086   auto_ptr<reco::TrackCollection> trackCollection(new reco::TrackCollection());
00087   // ... and its reference into the event
00088   reco::TrackRefProd trackCollectionRefProd = event.getRefBeforePut<reco::TrackCollection>();
00089 
00090 
00091   Handle<TrajectorySeedCollection> seeds;
00092   event.getByLabel(theSeedsLabel, seeds);
00093 
00094   for ( TrajectorySeedCollection::const_iterator iSeed = seeds->begin();
00095         iSeed != seeds->end(); iSeed++ ) {
00096     pair<bool,reco::Track> resultOfTrackExtrapAtPCA = buildTrackAtPCA(*iSeed);
00097     if(!resultOfTrackExtrapAtPCA.first) continue;
00098     // take the "bare" track at PCA
00099     reco::Track &track = resultOfTrackExtrapAtPCA.second;
00100     // fill the TrackCollection
00101     trackCollection->push_back(track);
00102   }
00103   
00104   event.put(trackCollection);
00105 
00106 }
00107 
00111 void
00112 MuonSeedTrack::beginJob()
00113 {
00114 }
00115 
00116 
00120 void
00121 MuonSeedTrack::endJob() {
00122 }
00123 
00124 
00128 TrajectoryStateOnSurface MuonSeedTrack::getSeedTSOS(const TrajectorySeed& seed) const{
00129 
00130   // Get the Trajectory State on Det (persistent version of a TSOS) from the seed
00131   PTrajectoryStateOnDet pTSOD = seed.startingState();
00132 
00133   // Transform it in a TrajectoryStateOnSurface
00134   
00135 
00136   DetId seedDetId(pTSOD.detId());
00137 
00138   const GeomDet* gdet = theService->trackingGeometry()->idToDet( seedDetId );
00139 
00140   TrajectoryStateOnSurface initialState = trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*theService->magneticField());
00141 
00142   /*
00143   // Get the layer on which the seed relies
00144   const DetLayer *initialLayer = theService->detLayerGeometry()->idToLayer( seedDetId );
00145 
00146   PropagationDirection detLayerOrder = oppositeToMomentum;
00147 
00148   // ask for compatible layers
00149   vector<const DetLayer*> detLayers;
00150   detLayers = initialLayer->compatibleLayers( *initialState.freeState(),detLayerOrder);
00151   
00152   TrajectoryStateOnSurface result = initialState;
00153   if(detLayers.size()){
00154     const DetLayer* finalLayer = detLayers.back();
00155     const TrajectoryStateOnSurface propagatedState = theService->propagator(theSeedPropagatorName)->propagate(initialState, finalLayer->surface());
00156     if(propagatedState.isValid())
00157       result = propagatedState;
00158   }
00159   
00160   return result;
00161   */
00162 
00163   return initialState;
00164 
00165 }
00166 
00167 
00172 pair<bool,reco::Track> MuonSeedTrack::buildTrackAtPCA(const TrajectorySeed& seed) const {
00173 
00174   const string metname = "MuonSeedTrack";
00175 
00176   MuonPatternRecoDumper debug;
00177 
00178   TSOS seedTSOS = getSeedTSOS(seed);
00179   // This is needed to extrapolate the tsos at vertex
00180   LogTrace(metname) << "Propagate to PCA...";
00181   pair<bool,FreeTrajectoryState> 
00182     extrapolationResult = theUpdatorAtVtx->propagateToNominalLine(seedTSOS);
00183   FreeTrajectoryState ftsAtVtx;
00184   
00185   if(extrapolationResult.first) {
00186     ftsAtVtx = extrapolationResult.second;
00187   }  else {    
00188     if(TrackerBounds::isInside(seedTSOS.globalPosition())){
00189       LogWarning(metname) << "Track in the Tracker: taking the innermost state instead of the state at PCA";
00190       ftsAtVtx = *seedTSOS.freeState();
00191     }
00192     else{
00193       if ( theAllowNoVtxFlag ) {
00194         LogWarning(metname) << "Propagation to PCA failed, taking the innermost state instead of the state at PCA";
00195         ftsAtVtx = *seedTSOS.freeState();
00196       } else {
00197         LogWarning(metname) << "Stand Alone track: this track will be rejected";
00198         return pair<bool,reco::Track>(false,reco::Track());
00199       }
00200     }
00201   }
00202     
00203   LogTrace(metname) << "TSOS after the extrapolation at vtx";
00204   LogTrace(metname) << debug.dumpFTS(ftsAtVtx);
00205   
00206   GlobalPoint pca = ftsAtVtx.position();
00207   math::XYZPoint persistentPCA(pca.x(),pca.y(),pca.z());
00208   GlobalVector p = ftsAtVtx.momentum();
00209   math::XYZVector persistentMomentum(p.x(),p.y(),p.z());
00210   
00211   double dummyNDOF = 1.0;
00212   //double ndof = computeNDOF(seed);
00213   double dummyChi2 = 1.0;
00214   
00215   reco::Track track(dummyChi2, 
00216                     dummyNDOF,
00217                     persistentPCA,
00218                     persistentMomentum,
00219                     ftsAtVtx.charge(),
00220                     ftsAtVtx.curvilinearError());
00221   
00222   return pair<bool,reco::Track>(true,track);
00223 }
00224 
00228 double MuonSeedTrack::computeNDOF(const TrajectorySeed& trajectory) const {
00229   const string metname = "MuonSeedTrack";
00230 
00231   BasicTrajectorySeed::const_iterator recHits1 = (trajectory.recHits().first);
00232   BasicTrajectorySeed::const_iterator recHits2 = (trajectory.recHits().second);
00233   
00234   double ndof = 0.;
00235 
00236   if ((*recHits1).isValid()) ndof += (*recHits1).dimension();
00237   if ((*recHits2).isValid()) ndof += (*recHits2).dimension();
00238 
00239   //const Trajectory::RecHitContainer transRecHits = trajectory.recHits();
00240   //for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin();
00241   //  rechit != transRecHits.end(); ++rechit)
00242   //if ((*rechit)->isValid()) ndof += (*rechit)->dimension();
00243   
00244   // FIXME! in case of Boff is dof - 4
00245   return max(ndof - 5., 0.);
00246 }
00247