CMS 3D CMS Logo

NuclearTrackCorrector.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: NuclearTrackCorrector
4 // Class: NuclearTrackCorrector
5 //
13 //
14 // Original Author: Loic QUERTENMONT, Vincent ROBERFROID
15 // Created: Tue Sep 18 14:22:48 CEST 2007
16 //
17 //
18 
25 
27 
29 
30 using namespace edm;
31 using namespace std;
32 using namespace reco;
33 
34 
36 conf_(iConfig),
37 theInitialState(nullptr)
38 {
39  str_Input_Trajectory = iConfig.getParameter<std::string> ("InputTrajectory");
40  str_Input_NuclearInteraction = iConfig.getParameter<std::string> ("InputNuclearInteraction");
41  verbosity = iConfig.getParameter<int> ("Verbosity");
42  KeepOnlyCorrectedTracks = iConfig.getParameter<bool> ("KeepOnlyCorrectedTracks");
43 
44 
46 
47  produces< TrajectoryCollection >();
48  produces< TrajectoryToTrajectoryMap >();
49 
50  produces< reco::TrackExtraCollection >();
51  produces< reco::TrackCollection >();
52  produces< TrackToTrajectoryMap >();
53 
54  produces< TrackToTrackMap >();
55 }
56 
57 
59 {
60 }
61 
62 void
64 {
65  // Create Output Collections
66  // --------------------------------------------------------------------------------------------------
67  auto Output_traj = std::make_unique<TrajectoryCollection>();
68  auto Output_trajmap = std::make_unique<TrajectoryToTrajectoryMap>();
69 
70  auto Output_trackextra = std::make_unique<reco::TrackExtraCollection>();
71  auto Output_track = std::make_unique<reco::TrackCollection>();
72  auto Output_trackmap = std::make_unique<TrackToTrajectoryMap>();
73 
74 
75 
76 
77 
78 
79  // Load Reccord
80  // --------------------------------------------------------------------------------------------------
81  std::string fitterName = conf_.getParameter<std::string>("Fitter");
82  iSetup.get<TrajectoryFitter::Record>().get(fitterName,theFitter);
83 
85  iSetup.get<TrackingComponentsRecord>().get(propagatorName,thePropagator);
86 
87  iSetup.get<TrackerDigiGeometryRecord>().get(theG);
88 
90 
91  iSetup.get<IdealMagneticFieldRecord>().get(theMF);
92 
93  // Load Inputs
94  // --------------------------------------------------------------------------------------------------
95  edm::Handle< TrajectoryCollection > temp_m_TrajectoryCollection;
96  iEvent.getByLabel( str_Input_Trajectory, temp_m_TrajectoryCollection );
97  const TrajectoryCollection m_TrajectoryCollection = *(temp_m_TrajectoryCollection.product());
98 
99  edm::Handle< NuclearInteractionCollection > temp_m_NuclearInteractionCollection;
100  iEvent.getByLabel( str_Input_NuclearInteraction, temp_m_NuclearInteractionCollection );
101  const NuclearInteractionCollection m_NuclearInteractionCollection = *(temp_m_NuclearInteractionCollection.product());
102 
103  edm::Handle< TrajTrackAssociationCollection > h_TrajToTrackCollection;
104  iEvent.getByLabel( str_Input_Trajectory, h_TrajToTrackCollection );
105  m_TrajToTrackCollection = h_TrajToTrackCollection.product();
106 
107 
108  // Correct the trajectories (Remove trajectory's hits that are located after the nuclear interacion)
109  // --------------------------------------------------------------------------------------------------
110  if(verbosity>=1){
111  LogDebug("NuclearTrackCorrector")
112  <<"Number of trajectories = "<<m_TrajectoryCollection.size() <<std::endl
113  <<"Number of nuclear interactions = "<<m_NuclearInteractionCollection.size();
114  }
115 
116  std::map<reco::TrackRef,TrajectoryRef> m_TrackToTrajMap;
117  swap_map(temp_m_TrajectoryCollection, m_TrackToTrajMap);
118 
119  for(unsigned int i = 0 ; i < m_NuclearInteractionCollection.size() ; i++)
120  {
121  reco::NuclearInteraction ni = m_NuclearInteractionCollection[i];
122  if( ni.likelihood()<0.4) continue;
123 
124  reco::TrackRef primTrackRef = ni.primaryTrack().castTo<reco::TrackRef>();
125 
126  TrajectoryRef trajRef = m_TrackToTrajMap[primTrackRef];
127 
128  Trajectory newTraj;
129  if( newTrajNeeded(newTraj, trajRef, ni) ) {
130 
131  AlgoProductCollection algoResults;
132  bool isOK = getTrackFromTrajectory( newTraj , trajRef, algoResults);
133 
134  if( isOK ) {
135 
136  pair<unsigned int, unsigned int> tempory_pair;
137  tempory_pair.first = Output_track->size();
138  tempory_pair.second = i;
139  Indice_Map.push_back(tempory_pair);
140 
141  reco::TrackExtraRef teref= reco::TrackExtraRef ( rTrackExtras, i );
142  reco::TrackExtra newTrackExtra = getNewTrackExtra(algoResults);
143  (algoResults[0].track)->setExtra( teref );
144 
145  Output_track->push_back(*algoResults[0].track);
146  Output_trackextra->push_back( newTrackExtra );
147  Output_traj->push_back(newTraj);
148 
149  }
150  }
151  else {
153  Output_track->push_back(*primTrackRef);
154  Output_trackextra->push_back( *primTrackRef->extra() );
155  Output_traj->push_back(*trajRef);
156  }
157  }
158 
159  }
160  const edm::OrphanHandle<TrajectoryCollection> Handle_traj = iEvent.put(std::move(Output_traj));
161  const edm::OrphanHandle<reco::TrackCollection> Handle_tracks = iEvent.put(std::move(Output_track));
162  iEvent.put(std::move(Output_trackextra));
163 
164  // Make Maps between elements
165  // --------------------------------------------------------------------------------------------------
166  if(Handle_tracks->size() != Handle_traj->size() )
167  {
168  printf("ERROR Handle_tracks->size() != Handle_traj->size() \n");
169  return;
170  }
171 
172  auto Output_tracktrackmap = std::make_unique<TrackToTrackMap>(Handle_tracks, m_TrajToTrackCollection->refProd().val);
173 
174  for(unsigned int i = 0 ; i < Indice_Map.size() ; i++)
175  {
176  TrajectoryRef InTrajRef ( temp_m_TrajectoryCollection, Indice_Map[i].second );
177  TrajectoryRef OutTrajRef ( Handle_traj, Indice_Map[i].first );
178  reco::TrackRef TrackRef ( Handle_tracks, Indice_Map[i].first );
179 
180  Output_trajmap ->insert(OutTrajRef,InTrajRef);
181  Output_trackmap->insert(TrackRef,InTrajRef);
182 
183  try{
184  reco::TrackRef PrimaryTrackRef = m_TrajToTrackCollection->operator[]( InTrajRef );
185  Output_tracktrackmap->insert(TrackRef,PrimaryTrackRef);
186  }catch(edm::Exception const&){}
187 
188  }
189  iEvent.put(std::move(Output_trajmap));
190  iEvent.put(std::move(Output_trackmap));
191  iEvent.put(std::move(Output_tracktrackmap));
192 
193 
194  if(verbosity>=3)printf("-----------------------\n");
195 }
196 
197 // ------------ method called once each job just after ending the event loop ------------
198 void
200 }
201 //----------------------------------------------------------------------------------------
202 bool NuclearTrackCorrector::newTrajNeeded(Trajectory& newtrajectory, const TrajectoryRef& trajRef, const NuclearInteraction& ni) {
203 
204  bool needNewTraj=false;
205  reco::Vertex::Point vtx_pos = ni.vertex().position();
206  double vtx_pos_mag = sqrt (vtx_pos.X()*vtx_pos.X()+vtx_pos.Y()*vtx_pos.Y()+vtx_pos.Z()*vtx_pos.Z());
207  if(verbosity>=2) printf("Nuclear Interaction pos = %f\n",vtx_pos_mag );
208 
209 
210  newtrajectory = Trajectory(trajRef->seed(), alongMomentum);
211 
212  // Look all the Hits of the trajectory and keep only Hits before seeds
213  Trajectory::DataContainer Measurements = trajRef->measurements();
214  if(verbosity>=2)LogDebug("NuclearTrackCorrector")<<"Size of Measurements = "<<Measurements.size();
215 
216  for(unsigned int m=Measurements.size()-1 ;m!=(unsigned int)-1 ; m--){
217 
218  if(!Measurements[m].recHit()->isValid() )continue;
219  GlobalPoint hit_pos = theG->idToDet(Measurements[m].recHit()->geographicalId())->surface().toGlobal(Measurements[m].recHit()->localPosition());
220 
221  if(verbosity>=2)printf("Hit pos = %f",hit_pos.mag() );
222 
223  if(hit_pos.mag()>vtx_pos_mag){
224  if(verbosity>=2)printf(" X ");
225  needNewTraj=true;
226  }else{
227  newtrajectory.push(Measurements[m]);
228  }
229  if(verbosity>=2)printf("\n");
230  }
231 
232  return needNewTraj;
233 }
234 
235 //----------------------------------------------------------------------------------------
236 bool NuclearTrackCorrector::getTrackFromTrajectory(const Trajectory& newTraj , const TrajectoryRef& initialTrajRef, AlgoProductCollection& algoResults) {
237 
238  const Trajectory* it = &newTraj;
239 
241  it->validRecHits( hits );
242 
243 
244  float ndof=0;
245  for(unsigned int h=0 ; h<hits.size() ; h++)
246  {
247  if( hits[h]->isValid() )
248  {
249  ndof = ndof + hits[h]->dimension() * hits[h]->weight();
250  }
251  else {
252  LogDebug("NuclearSeedGenerator") << " HIT IS INVALID ???";
253  }
254  }
255 
256 
257  ndof = ndof - 5;
258  reco::TrackRef theT = m_TrajToTrackCollection->operator[]( initialTrajRef );
259  LogDebug("NuclearSeedGenerator") << " TrackCorrector - number of valid hits" << hits.size() << "\n"
260  << " - number of hits from Track " << theT->recHitsSize() << "\n"
261  << " - number of valid hits from initial track " << theT->numberOfValidHits();
262 
263 
264  if( hits.size() > 1){
265 
266  TrajectoryStateOnSurface theInitialStateForRefitting = getInitialState(&(*theT),hits,theG.product(),theMF.product()
267 );
268 
269  reco::BeamSpot bs;
270  return theAlgo->buildTrack(theFitter.product(), thePropagator.product(), algoResults, hits, theInitialStateForRefitting ,it->seed(), ndof, bs, theT->seedRef());
271  }
272 
273  return false;
274 }
275 //----------------------------------------------------------------------------------------
277  Trajectory* theTraj = algoResults[0].trajectory;
278  PropagationDirection seedDir = algoResults[0].pDir;
279 
280  TrajectoryStateOnSurface outertsos;
281  TrajectoryStateOnSurface innertsos;
282  unsigned int innerId, outerId;
283  if (theTraj->direction() == alongMomentum) {
284  outertsos = theTraj->lastMeasurement().updatedState();
285  innertsos = theTraj->firstMeasurement().updatedState();
286  outerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
287  innerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
288  } else {
289  outertsos = theTraj->firstMeasurement().updatedState();
290  innertsos = theTraj->lastMeasurement().updatedState();
291  outerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
292  innerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
293  }
294 
295  GlobalPoint v = outertsos.globalParameters().position();
296  GlobalVector p = outertsos.globalParameters().momentum();
297  math::XYZVector outmom( p.x(), p.y(), p.z() );
298  math::XYZPoint outpos( v.x(), v.y(), v.z() );
299  v = innertsos.globalParameters().position();
300  p = innertsos.globalParameters().momentum();
301  math::XYZVector inmom( p.x(), p.y(), p.z() );
302  math::XYZPoint inpos( v.x(), v.y(), v.z() );
303 
304  return reco::TrackExtra (outpos, outmom, true, inpos, inmom, true,
305  outertsos.curvilinearError(), outerId,
306  innertsos.curvilinearError(), innerId, seedDir);
307 
308 }
309 //----------------------------------------------------------------------------------------
312  const TrackingGeometry * theG,
313  const MagneticField * theMF){
314 
315  TrajectoryStateOnSurface theInitialStateForRefitting;
316  //the starting state is the state closest to the first hit along seedDirection.
317 
318  //avoiding to use transientTrack, it should be faster;
319  TrajectoryStateOnSurface innerStateFromTrack=trajectoryStateTransform::innerStateOnSurface(*theT,*theG,theMF);
320  TrajectoryStateOnSurface outerStateFromTrack=trajectoryStateTransform::outerStateOnSurface(*theT,*theG,theMF);
321  TrajectoryStateOnSurface initialStateFromTrack =
322  ( (innerStateFromTrack.globalPosition()-hits.front()->globalPosition()).mag2() <
323  (outerStateFromTrack.globalPosition()-hits.front()->globalPosition()).mag2() ) ?
324  innerStateFromTrack: outerStateFromTrack;
325 
326  // error is rescaled, but correlation are kept.
327  initialStateFromTrack.rescaleError(100);
328  theInitialStateForRefitting = TrajectoryStateOnSurface(initialStateFromTrack.localParameters(),
329  initialStateFromTrack.localError(),
330  initialStateFromTrack.surface(),
331  theMF);
332  return theInitialStateForRefitting;
333 }
334 //----------------------------------------------------------------------------------------
336  for(unsigned int i = 0 ; i < trajColl->size() ; i++)
337  {
338  TrajectoryRef InTrajRef ( trajColl, i);
339  reco::TrackRef PrimaryTrackRef = m_TrajToTrackCollection->operator[]( InTrajRef );
340  result[ PrimaryTrackRef ] = InTrajRef;
341  }
342 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
T getParameter(std::string const &) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
ConstRecHitPointer const & recHit() const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
TrajectorySeed const & seed() const
Access to the seed used to reconstruct the Trajectory.
Definition: Trajectory.h:285
const LocalTrajectoryParameters & localParameters() const
TrajectoryStateOnSurface outerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field, bool withErr=true)
NuclearTrackCorrector(const edm::ParameterSet &)
edm::ESHandle< MagneticField > theMF
const edm::RefToBase< reco::Track > & primaryTrack() const
return the base reference to the primary track
void produce(edm::Event &, const edm::EventSetup &) override
const CurvilinearTrajectoryError & curvilinearError() const
std::vector< std::pair< unsigned int, unsigned int > > Indice_Map
#define nullptr
T y() const
Definition: PV3DBase.h:63
std::vector< ConstRecHitPointer > RecHitContainer
GlobalPoint globalPosition() const
bool getTrackFromTrajectory(const Trajectory &newTraj, const TrajectoryRef &initialTrajRef, AlgoProductCollection &algoResults)
Get the refitted track from the Trajectory.
edm::Ref< TrackExtraCollection > TrackExtraRef
persistent reference to a TrackExtra
Definition: TrackExtraFwd.h:17
PropagationDirection
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
const Point & position() const
position
Definition: Vertex.h:109
bool buildTrack(const TrajectoryFitter *, const Propagator *, AlgoProductCollection &, TransientTrackingRecHit::RecHitContainer &, TrajectoryStateOnSurface &, const TrajectorySeed &, float, const reco::BeamSpot &, SeedRef seedRef=SeedRef(), int qualityMask=0, signed char nLoops=0)
Construct Tracks to be put in the event.
PropagationDirection const & direction() const
Definition: Trajectory.cc:140
U second(std::pair< T, U > const &p)
int iEvent
Definition: GenABIO.cc:224
const SurfaceType & surface() const
std::vector< TrajectoryMeasurement > DataContainer
Definition: Trajectory.h:44
T mag() const
Definition: PV3DBase.h:67
const ref_type & refProd() const
return ref-prod structure
edm::ESHandle< TrackerGeometry > theG
T sqrt(T t)
Definition: SSEVec.h:18
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:174
T z() const
Definition: PV3DBase.h:64
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
void validRecHits(ConstRecHitContainer &cont) const
Definition: Trajectory.cc:133
std::vector< NuclearInteraction > NuclearInteractionCollection
collection of NuclearInteractions
TrackProducerAlgorithm< reco::Track > * theAlgo
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
const LocalTrajectoryError & localError() const
double likelihood() const
return the likelihood ~ probability that the vertex is a real nuclear interaction ...
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:480
RefProd< PROD > getRefBeforePut()
Definition: Event.h:150
edm::ESHandle< Propagator > thePropagator
edm::ESHandle< TrajectoryFitter > theFitter
std::vector< TrackExtra > TrackExtraCollection
collection of TrackExtra objects
Definition: TrackExtraFwd.h:11
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:187
T const * product() const
Definition: Handle.h:74
void insert(const key_type &k, const data_type &v)
insert an association
const GlobalTrajectoryParameters & globalParameters() const
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
reco::TrackExtra getNewTrackExtra(const AlgoProductCollection &algoresults)
get a new TrackExtra from an AlgoProductCollection
REF castTo() const
Definition: RefToBase.h:289
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
Definition: TrackFwd.h:21
bool newTrajNeeded(Trajectory &newtrajectory, const TrajectoryRef &trajRef, const reco::NuclearInteraction &ni)
check if the trajectory has to be refitted and get the new trajectory
std::vector< Trajectory > TrajectoryCollection
const TrajTrackAssociationCollection * m_TrajToTrackCollection
const reco::Vertex & vertex() const
return the vertex
fixed size matrix
HLT enums.
T get() const
Definition: EventSetup.h:71
const TrackerGeomDet * idToDet(DetId) const override
TrajectoryStateOnSurface const & updatedState() const
T x() const
Definition: PV3DBase.h:62
void swap_map(const edm::Handle< TrajectoryCollection > &trajColl, std::map< reco::TrackRef, edm::Ref< TrajectoryCollection > > &result)
T const * product() const
Definition: ESHandle.h:86
void push(const TrajectoryMeasurement &tm)
Definition: Trajectory.cc:50
TrackProducerAlgorithm< reco::Track >::AlgoProductCollection AlgoProductCollection
TrajectoryStateOnSurface getInitialState(const reco::Track *theT, TransientTrackingRecHit::RecHitContainer &hits, const TrackingGeometry *theG, const MagneticField *theMF)
Calculate the inital state to be used to buil the track.
def move(src, dest)
Definition: eostools.py:511
TrajectoryStateOnSurface innerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field, bool withErr=true)
std::string str_Input_NuclearInteraction