00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018
00019
00020 #include "RecoTracker/NuclearSeedGenerator/interface/NuclearTrackCorrector.h"
00021 #include "FWCore/Framework/interface/EventSetup.h"
00022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00023 #include "RecoTracker/CkfPattern/interface/TransientInitialStateEstimator.h"
00024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00025 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00026
00027 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00028
00029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00030
00031 using namespace edm;
00032 using namespace std;
00033 using namespace reco;
00034
00035
00036 NuclearTrackCorrector::NuclearTrackCorrector(const edm::ParameterSet& iConfig) :
00037 conf_(iConfig),
00038 theInitialState(0)
00039 {
00040 str_Input_Trajectory = iConfig.getParameter<std::string> ("InputTrajectory");
00041 str_Input_NuclearInteraction = iConfig.getParameter<std::string> ("InputNuclearInteraction");
00042 verbosity = iConfig.getParameter<int> ("Verbosity");
00043 KeepOnlyCorrectedTracks = iConfig.getParameter<bool> ("KeepOnlyCorrectedTracks");
00044
00045
00046 theAlgo = new TrackProducerAlgorithm<reco::Track>(iConfig);
00047
00048 produces< TrajectoryCollection >();
00049 produces< TrajectoryToTrajectoryMap >();
00050
00051 produces< reco::TrackExtraCollection >();
00052 produces< reco::TrackCollection >();
00053 produces< TrackToTrajectoryMap >();
00054
00055 produces< TrackToTrackMap >();
00056 }
00057
00058
00059 NuclearTrackCorrector::~NuclearTrackCorrector()
00060 {
00061 }
00062
00063 void
00064 NuclearTrackCorrector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00065 {
00066
00067
00068 std::auto_ptr<TrajectoryCollection> Output_traj ( new TrajectoryCollection );
00069 std::auto_ptr<TrajectoryToTrajectoryMap> Output_trajmap ( new TrajectoryToTrajectoryMap );
00070
00071 std::auto_ptr<reco::TrackExtraCollection> Output_trackextra ( new reco::TrackExtraCollection );
00072 std::auto_ptr<reco::TrackCollection> Output_track ( new reco::TrackCollection );
00073 std::auto_ptr<TrackToTrajectoryMap> Output_trackmap ( new TrackToTrajectoryMap );
00074
00075 std::auto_ptr<TrackToTrackMap> Output_tracktrackmap ( new TrackToTrackMap );
00076
00077
00078
00079
00080
00081
00082
00083 std::string fitterName = conf_.getParameter<std::string>("Fitter");
00084 iSetup.get<TrajectoryFitter::Record>().get(fitterName,theFitter);
00085
00086 std::string propagatorName = conf_.getParameter<std::string>("Propagator");
00087 iSetup.get<TrackingComponentsRecord>().get(propagatorName,thePropagator);
00088
00089 iSetup.get<TrackerDigiGeometryRecord>().get(theG);
00090
00091 reco::TrackExtraRefProd rTrackExtras = iEvent.getRefBeforePut<reco::TrackExtraCollection>();
00092
00093 iSetup.get<IdealMagneticFieldRecord>().get(theMF);
00094
00095
00096
00097 edm::Handle< TrajectoryCollection > temp_m_TrajectoryCollection;
00098 iEvent.getByLabel( str_Input_Trajectory.c_str(), temp_m_TrajectoryCollection );
00099 const TrajectoryCollection m_TrajectoryCollection = *(temp_m_TrajectoryCollection.product());
00100
00101 edm::Handle< NuclearInteractionCollection > temp_m_NuclearInteractionCollection;
00102 iEvent.getByLabel( str_Input_NuclearInteraction.c_str(), temp_m_NuclearInteractionCollection );
00103 const NuclearInteractionCollection m_NuclearInteractionCollection = *(temp_m_NuclearInteractionCollection.product());
00104
00105 edm::Handle< TrajTrackAssociationCollection > h_TrajToTrackCollection;
00106 iEvent.getByLabel( str_Input_Trajectory.c_str(), h_TrajToTrackCollection );
00107 m_TrajToTrackCollection = h_TrajToTrackCollection.product();
00108
00109
00110
00111
00112 if(verbosity>=1){
00113 LogDebug("NuclearTrackCorrector")
00114 <<"Number of trajectories = "<<m_TrajectoryCollection.size() <<std::endl
00115 <<"Number of nuclear interactions = "<<m_NuclearInteractionCollection.size();
00116 }
00117
00118 std::map<reco::TrackRef,TrajectoryRef> m_TrackToTrajMap;
00119 swap_map(temp_m_TrajectoryCollection, m_TrackToTrajMap);
00120
00121 for(unsigned int i = 0 ; i < m_NuclearInteractionCollection.size() ; i++)
00122 {
00123 reco::NuclearInteraction ni = m_NuclearInteractionCollection[i];
00124 if( ni.likelihood()<0.4) continue;
00125
00126 reco::TrackRef primTrackRef = ni.primaryTrack().castTo<reco::TrackRef>();
00127
00128 TrajectoryRef trajRef = m_TrackToTrajMap[primTrackRef];
00129
00130 Trajectory newTraj;
00131 if( newTrajNeeded(newTraj, trajRef, ni) ) {
00132
00133 AlgoProductCollection algoResults;
00134 bool isOK = getTrackFromTrajectory( newTraj , trajRef, algoResults);
00135
00136 if( isOK ) {
00137
00138 pair<unsigned int, unsigned int> tempory_pair;
00139 tempory_pair.first = Output_track->size();
00140 tempory_pair.second = i;
00141 Indice_Map.push_back(tempory_pair);
00142
00143 reco::TrackExtraRef teref= reco::TrackExtraRef ( rTrackExtras, i );
00144 reco::TrackExtra newTrackExtra = getNewTrackExtra(algoResults);
00145 (algoResults[0].second.first)->setExtra( teref );
00146
00147 Output_track->push_back(*algoResults[0].second.first);
00148 Output_trackextra->push_back( newTrackExtra );
00149 Output_traj->push_back(newTraj);
00150
00151 }
00152 }
00153 else {
00154 if(!KeepOnlyCorrectedTracks) {
00155 Output_track->push_back(*primTrackRef);
00156 Output_trackextra->push_back( *primTrackRef->extra() );
00157 Output_traj->push_back(*trajRef);
00158 }
00159 }
00160
00161 }
00162 const edm::OrphanHandle<TrajectoryCollection> Handle_traj = iEvent.put(Output_traj);
00163 const edm::OrphanHandle<reco::TrackCollection> Handle_tracks = iEvent.put(Output_track);
00164 iEvent.put(Output_trackextra);
00165
00166
00167
00168 if(Handle_tracks->size() != Handle_traj->size() )
00169 {
00170 printf("ERROR Handle_tracks->size() != Handle_traj->size() \n");
00171 return;
00172 }
00173
00174
00175
00176 for(unsigned int i = 0 ; i < Indice_Map.size() ; i++)
00177 {
00178 TrajectoryRef InTrajRef ( temp_m_TrajectoryCollection, Indice_Map[i].second );
00179 TrajectoryRef OutTrajRef ( Handle_traj, Indice_Map[i].first );
00180 reco::TrackRef TrackRef ( Handle_tracks, Indice_Map[i].first );
00181
00182 Output_trajmap ->insert(OutTrajRef,InTrajRef);
00183 Output_trackmap->insert(TrackRef,InTrajRef);
00184
00185 try{
00186 reco::TrackRef PrimaryTrackRef = m_TrajToTrackCollection->operator[]( InTrajRef );
00187 Output_tracktrackmap->insert(TrackRef,PrimaryTrackRef);
00188 }catch(edm::Exception event){}
00189
00190 }
00191 iEvent.put(Output_trajmap);
00192 iEvent.put(Output_trackmap);
00193 iEvent.put(Output_tracktrackmap);
00194
00195
00196 if(verbosity>=3)printf("-----------------------\n");
00197 }
00198
00199
00200 void
00201 NuclearTrackCorrector::beginRun(edm::Run & run, const edm::EventSetup& iSetup)
00202 {
00203 }
00204
00205
00206 void
00207 NuclearTrackCorrector::endJob() {
00208 }
00209
00210 bool NuclearTrackCorrector::newTrajNeeded(Trajectory& newtrajectory, const TrajectoryRef& trajRef, const NuclearInteraction& ni) {
00211
00212 bool needNewTraj=false;
00213 reco::Vertex::Point vtx_pos = ni.vertex().position();
00214 double vtx_pos_mag = sqrt (vtx_pos.X()*vtx_pos.X()+vtx_pos.Y()*vtx_pos.Y()+vtx_pos.Z()*vtx_pos.Z());
00215 if(verbosity>=2) printf("Nuclear Interaction pos = %f\n",vtx_pos_mag );
00216
00217
00218 newtrajectory = Trajectory(trajRef->seed(), alongMomentum);
00219
00220
00221 Trajectory::DataContainer Measurements = trajRef->measurements();
00222 if(verbosity>=2)LogDebug("NuclearTrackCorrector")<<"Size of Measurements = "<<Measurements.size();
00223
00224 for(unsigned int m=Measurements.size()-1 ;m!=(unsigned int)-1 ; m--){
00225
00226 if(!Measurements[m].recHit()->isValid() )continue;
00227 GlobalPoint hit_pos = theG->idToDet(Measurements[m].recHit()->geographicalId())->surface().toGlobal(Measurements[m].recHit()->localPosition());
00228
00229 if(verbosity>=2)printf("Hit pos = %f",hit_pos.mag() );
00230
00231 if(hit_pos.mag()>vtx_pos_mag){
00232 if(verbosity>=2)printf(" X ");
00233 needNewTraj=true;
00234 }else{
00235 newtrajectory.push(Measurements[m]);
00236 }
00237 if(verbosity>=2)printf("\n");
00238 }
00239
00240 return needNewTraj;
00241 }
00242
00243
00244 bool NuclearTrackCorrector::getTrackFromTrajectory(const Trajectory& newTraj , const TrajectoryRef& initialTrajRef, AlgoProductCollection& algoResults) {
00245
00246 const Trajectory* it = &newTraj;
00247
00248 TransientTrackingRecHit::RecHitContainer hits;
00249 it->validRecHits( hits );
00250
00251
00252 float ndof=0;
00253 for(unsigned int h=0 ; h<hits.size() ; h++)
00254 {
00255 if( hits[h]->isValid() )
00256 {
00257 ndof = ndof + hits[h]->dimension() * hits[h]->weight();
00258 }
00259 else {
00260 LogDebug("NuclearSeedGenerator") << " HIT IS INVALID ???";
00261 }
00262 }
00263
00264
00265 ndof = ndof - 5;
00266 reco::TrackRef theT = m_TrajToTrackCollection->operator[]( initialTrajRef );
00267 LogDebug("NuclearSeedGenerator") << " TrackCorrector - number of valid hits" << hits.size() << "\n"
00268 << " - number of hits from Track " << theT->recHitsSize() << "\n"
00269 << " - number of valid hits from initial track " << theT->numberOfValidHits();
00270
00271
00272 if( hits.size() > 1){
00273
00274 TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(&(*theT),hits,theG.product(),theMF.product()
00275 );
00276
00277 reco::BeamSpot bs;
00278 return theAlgo->buildTrack(theFitter.product(), thePropagator.product(), algoResults, hits, theInitialStateForRefitting ,it->seed(), ndof, bs, theT->seedRef());
00279 }
00280
00281 return false;
00282 }
00283
00284 reco::TrackExtra NuclearTrackCorrector::getNewTrackExtra(const AlgoProductCollection& algoResults) {
00285 Trajectory* theTraj = algoResults[0].first;
00286 PropagationDirection seedDir = algoResults[0].second.second;
00287
00288 TrajectoryStateOnSurface outertsos;
00289 TrajectoryStateOnSurface innertsos;
00290 unsigned int innerId, outerId;
00291 if (theTraj->direction() == alongMomentum) {
00292 outertsos = theTraj->lastMeasurement().updatedState();
00293 innertsos = theTraj->firstMeasurement().updatedState();
00294 outerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
00295 innerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
00296 } else {
00297 outertsos = theTraj->firstMeasurement().updatedState();
00298 innertsos = theTraj->lastMeasurement().updatedState();
00299 outerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
00300 innerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
00301 }
00302
00303 GlobalPoint v = outertsos.globalParameters().position();
00304 GlobalVector p = outertsos.globalParameters().momentum();
00305 math::XYZVector outmom( p.x(), p.y(), p.z() );
00306 math::XYZPoint outpos( v.x(), v.y(), v.z() );
00307 v = innertsos.globalParameters().position();
00308 p = innertsos.globalParameters().momentum();
00309 math::XYZVector inmom( p.x(), p.y(), p.z() );
00310 math::XYZPoint inpos( v.x(), v.y(), v.z() );
00311
00312 return reco::TrackExtra (outpos, outmom, true, inpos, inmom, true,
00313 outertsos.curvilinearError(), outerId,
00314 innertsos.curvilinearError(), innerId, seedDir);
00315
00316 }
00317
00318 TrajectoryStateOnSurface NuclearTrackCorrector::getInitialState(const reco::Track * theT,
00319 TransientTrackingRecHit::RecHitContainer& hits,
00320 const TrackingGeometry * theG,
00321 const MagneticField * theMF){
00322
00323 TrajectoryStateOnSurface theInitialStateForRefitting;
00324
00325 TrajectoryStateTransform transformer;
00326
00327 TrajectoryStateOnSurface innerStateFromTrack=transformer.innerStateOnSurface(*theT,*theG,theMF);
00328 TrajectoryStateOnSurface outerStateFromTrack=transformer.outerStateOnSurface(*theT,*theG,theMF);
00329 TrajectoryStateOnSurface initialStateFromTrack =
00330 ( (innerStateFromTrack.globalPosition()-hits.front()->globalPosition()).mag2() <
00331 (outerStateFromTrack.globalPosition()-hits.front()->globalPosition()).mag2() ) ?
00332 innerStateFromTrack: outerStateFromTrack;
00333
00334
00335 initialStateFromTrack.rescaleError(100);
00336 theInitialStateForRefitting = TrajectoryStateOnSurface(initialStateFromTrack.localParameters(),
00337 initialStateFromTrack.localError(),
00338 initialStateFromTrack.surface(),
00339 theMF);
00340 return theInitialStateForRefitting;
00341 }
00342
00343 void NuclearTrackCorrector::swap_map(const edm::Handle< TrajectoryCollection >& trajColl , std::map< reco::TrackRef, edm::Ref<TrajectoryCollection> >& result) {
00344 for(unsigned int i = 0 ; i < trajColl->size() ; i++)
00345 {
00346 TrajectoryRef InTrajRef ( trajColl, i);
00347 reco::TrackRef PrimaryTrackRef = m_TrajToTrackCollection->operator[]( InTrajRef );
00348 result[ PrimaryTrackRef ] = InTrajRef;
00349 }
00350 }
1.7.3