CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PlotRecTracks.cc
Go to the documentation of this file.
5 
7 
11 
15 
21 
24 
28 
31 
34 
37 
39 
41 
43 
44 using namespace std;
45 
46 /*****************************************************************************/
48  (const edm::EventSetup& es_, std::string trackProducer_,
49  ofstream& file_) : es(es_), trackProducer(trackProducer_), file(file_)
50 {
51  // Get tracker geometry
52  edm::ESHandle<TrackerGeometry> trackerHandle;
53  es.get<TrackerDigiGeometryRecord>().get(trackerHandle);
54  theTracker = trackerHandle.product();
55 
56  // Get magnetic field
58  es.get<IdealMagneticFieldRecord>().get(magField);
59  theMagField = magField.product();
60 
61  // Get propagator
62  edm::ESHandle<Propagator> thePropagatorHandle;
63  es.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",
64  thePropagatorHandle);
65  thePropagator = thePropagatorHandle.product();
66 
67  // KFTrajectoryFitter
68 /*
69  edm::ESHandle<TrajectoryFitter> theFitterHandle;
70  es.get<TrackingComponentsRecord>().get("KFTrajectoryFitter", theFitterHandle);
71  theFitter = theFitterHandle.product();
72 */
73 }
74 
75 /*****************************************************************************/
77 {
78 }
79 
80 /*****************************************************************************/
82  (const TrackingRecHit* recHit, const ostringstream& o, const ostringstream& d)
83 {
84  const SiPixelRecHit* pixelRecHit =
85  dynamic_cast<const SiPixelRecHit *>(recHit);
86 
87  DetId id = recHit->geographicalId();
88  LocalPoint lpos = recHit->localPosition();
89  GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
90 
91  SiPixelRecHit::ClusterRef const& cluster = pixelRecHit->cluster();
92  std::vector<SiPixelCluster::Pixel> pixels = cluster->pixels();
93 
94  std::vector<PSimHit> simHits = theHitAssociator->associateHit(*recHit);
95 
96  std::string info = ", Text[StyleForm[\"" + o.str() + "\", URL->\"Track " + o.str()
97 + d.str();
98 
99  {
100  ostringstream o;
101  o << "simTrack (trackId=" << simHits[0].trackId()
102  << " pid=" << simHits[0].particleType()
103  << " proc=" << simHits[0].processType() << ")";
104 
105  info += " | " + o.str();
106  }
107 
108  {
109  ostringstream o;
110  o << theTracker->idToDet(id)->subDetector();
111 
112  info += " | " + o.str();
113  }
114 
115  info += HitInfo::getInfo(*recHit);
116 
117  info += "\"]";
118 
119  {
120  ostringstream o;
121  o << ", {" << p.x() << "," << p.y() << ",(" << p.z() << "-zs)*mz},"
122  << " {" << 0 << "," << -1 << "}]";
123 
124  info += o.str();
125  }
126 
127  return info;
128 }
129 
130 /*****************************************************************************/
132  (const TrackingRecHit* recHit, const ostringstream& o, const ostringstream& d)
133 {
134  DetId id = recHit->geographicalId();
135  LocalPoint lpos = recHit->localPosition();
136  GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
137 
138  std::vector<PSimHit> simHits = theHitAssociator->associateHit(*recHit);
139 
140  std::string info = ", Text[StyleForm[\"" + o.str() + "\", URL->\"Track " + o.str() + d.str();
141  const SiStripMatchedRecHit2D* stripMatchedRecHit =
142  dynamic_cast<const SiStripMatchedRecHit2D *>(recHit);
143  const ProjectedSiStripRecHit2D* stripProjectedRecHit =
144  dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit);
145  const SiStripRecHit2D* stripRecHit =
146  dynamic_cast<const SiStripRecHit2D *>(recHit);
147 
148  {
149  ostringstream o;
150  o << "simTrackId=" << simHits[0].trackId() ;
151 
152  info += " | " + o.str();
153  }
154 
155  {
156  ostringstream o;
157  o << theTracker->idToDet(id)->subDetector();
158 
159  info += " | " + o.str();
160  }
161 
162  info += HitInfo::getInfo(*recHit);
163 
164  if(stripMatchedRecHit != 0) info += " matched";
165  if(stripProjectedRecHit != 0) info += " projected";
166  if(stripRecHit != 0) info += " single";
167 
168 
169  info += "\"]";
170 
171  {
172  ostringstream o;
173  o << ", {" << p.x() << "," << p.y() << ",(" << p.z() << "-zs)*mz},"
174  << " {" << 0 << "," << -1 << "}]";
175 
176  info += o.str();
177  }
178 
179  return info;
180 }
181 
182 /*****************************************************************************/
184  (const reco::Track& track)
185 {
186  GlobalPoint position(track.outerX(),
187  track.outerY(),
188  track.outerZ());
189 
190  GlobalVector momentum(track.outerPx(),
191  track.outerPy(),
192  track.outerPz());
193 
195  track.charge(),theMagField);
196 
197  FreeTrajectoryState fts(gtp,track.outerStateCovariance());
198 
199  return fts;
200 }
201 
202 /*****************************************************************************/
204 {
205  theHitAssociator = new TrackerHitAssociator(ev);
206 
207  file << ", If[rt, {AbsolutePointSize[6]";
208 
209  edm::Handle<reco::TrackCollection> recTrackHandle;
210  ev.getByLabel(trackProducer, recTrackHandle);
211 
212  edm::Handle<std::vector<Trajectory> > trajectoryHandle;
213  ev.getByLabel(trackProducer, trajectoryHandle);
214 
215  const reco::TrackCollection* recTracks = recTrackHandle.product();
216  const std::vector<Trajectory>* trajectories = trajectoryHandle.product();
217 
218  edm::LogVerbatim("MinBiasTracking")
219  << " [EventPlotter] recTracks (" << trackProducer << ") "
220  << recTracks->size();
221 
222  PlotRecHits theRecHits(es,file);
223 
224  file << ", RGBColor[0,0,0.4]";
225  reco::TrackCollection::const_iterator recTrack = recTracks->begin();
226 
227  int i = 0;
228  for(std::vector<Trajectory>::const_iterator it = trajectories->begin();
229  it!= trajectories->end();
230  it++, i++, recTrack++)
231  {
232 /*
233 cerr << " track[" << i << "] " << recTrack->chi2() << " " << it->chiSquared() << std::endl;
234 
235 */
236 //theFitter->fit(*it);
237 
238  ostringstream o; o << i;
239 
240  ostringstream d; d << fixed << std::setprecision(2)
241  << " | d0=" << recTrack->d0() << " cm"
242  << " | z0=" << recTrack->dz() << " cm"
243  << " | pt=" << recTrack->pt() << " GeV/c";
244 
245  const Trajectory* trajectory = &(*it);
246 
247  for(std::vector<TrajectoryMeasurement>::const_reverse_iterator
248  meas = trajectory->measurements().rbegin();
249  meas!= trajectory->measurements().rend(); meas++)
250  {
251  const TrackingRecHit* recHit = meas->recHit()->hit();
252 
253  if(recHit->isValid())
254  {
255  DetId id = recHit->geographicalId();
256 
257  if(theTracker->idToDet(id)->subDetector() ==
259  theTracker->idToDet(id)->subDetector() ==
261  {
262  // Pixel
263  const SiPixelRecHit* pixelRecHit =
264  dynamic_cast<const SiPixelRecHit *>(recHit);
265 
266  if(pixelRecHit != 0)
267  {
268  theRecHits.printPixelRecHit(pixelRecHit);
269  file << getPixelInfo(recHit, o,d);
270  }
271  }
272  else
273  {
274  // Strip
275  const SiStripMatchedRecHit2D* stripMatchedRecHit =
276  dynamic_cast<const SiStripMatchedRecHit2D *>(recHit);
277  const ProjectedSiStripRecHit2D* stripProjectedRecHit =
278  dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit);
279  const SiStripRecHit2D* stripRecHit =
280  dynamic_cast<const SiStripRecHit2D *>(recHit);
281 
282  if(stripMatchedRecHit != 0)
283  {
284  auto m = stripMatchedRecHit->monoHit();
285  auto s = stripMatchedRecHit->stereoHit();
286  theRecHits.printStripRecHit(&m);
287  theRecHits.printStripRecHit(&s);
288 
289  DetId id = stripMatchedRecHit->geographicalId();
290  LocalPoint lpos = stripMatchedRecHit->localPosition();
291  GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
292 
293  file << ", Point[{"<< p.x()<<","<<p.y()<<",("<<p.z()<<"-zs)*mz}]" << std::endl;
294  }
295 
296  if(stripProjectedRecHit != 0)
297  theRecHits.printStripRecHit(&(stripProjectedRecHit->originalHit()));
298 
299  if(stripRecHit != 0)
300  theRecHits.printStripRecHit(stripRecHit);
301 
302  if(stripMatchedRecHit != 0 ||
303  stripProjectedRecHit != 0 ||
304  stripRecHit != 0)
305  file << getStripInfo(recHit, o,d);
306  }
307  }
308  }
309  }
310 
311  PlotUtils plotUtils;
312 
313  // Trajectory
314  recTrack = recTracks->begin();
315 
316  for(std::vector<Trajectory>::const_iterator it = trajectories->begin();
317  it!= trajectories->end();
318  it++, recTrack++)
319  {
320  int algo;
321  switch(recTrack->algo())
322  {
323  case reco::TrackBase::iter1: algo = 0; break;
324  case reco::TrackBase::iter2: algo = 1; break;
325  case reco::TrackBase::iter3: algo = 2; break;
326  default: algo = 0;
327  }
328 
329  if(algo == 0) file << ", RGBColor[1,0,0]";
330  if(algo == 1) file << ", RGBColor[0.2,0.6,0.2]";
331  if(algo == 2) file << ", RGBColor[0.2,0.2,0.6]";
332 
333  std::vector<TrajectoryMeasurement> meas = it->measurements();
334 
335  for(std::vector<TrajectoryMeasurement>::reverse_iterator im = meas.rbegin();
336  im!= meas.rend(); im++)
337  {
338  if(im == meas.rbegin())
339  {
340  GlobalPoint p2 = (*(im )).updatedState().globalPosition();
341  GlobalVector v2 = (*(im )).updatedState().globalDirection();
342  GlobalPoint p1 = GlobalPoint(recTrack->vertex().x(),
343  recTrack->vertex().y(),
344  recTrack->vertex().z());
345 
346  plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
347  }
348 
349  if(im+1 != meas.rend())
350  {
351  GlobalPoint p1 = (*(im )).updatedState().globalPosition();
352  GlobalPoint p2 = (*(im+1)).updatedState().globalPosition();
353  GlobalVector v2 = (*(im+1)).updatedState().globalDirection();
354 
355  plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
356  }
357  }
358 
359  // Ecal
360  GlobalPoint p1 = (*(meas.rend()-1)).forwardPredictedState().globalPosition();
361  FreeTrajectoryState fts = getTrajectoryAtOuterPoint(*recTrack);
364 
365  // Ecal Barrel
366  Cylinder::ConstCylinderPointer theCylinder =
367  Cylinder::build(Surface::PositionType(0.,0.,0), rot, 129.);
368  tsos = thePropagator->propagate(fts,*theCylinder);
369 
370  if(tsos.isValid() &&
371  fabs(tsos.globalPosition().z()) < 320.9)
372  {
373  GlobalPoint p2 = tsos.globalPosition();
374  GlobalVector v2 = tsos.globalDirection();
375  plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
376  }
377  else
378  {
379  // ECAL Endcap
380  Plane::ConstPlanePointer thePlanePos =
381  Plane::build(Surface::PositionType(0.,0.,320.9), rot);
382  tsos = thePropagator->propagate(fts,*thePlanePos);
383 
384  if(tsos.isValid())
385  {
386  GlobalPoint p2 = tsos.globalPosition();
387  GlobalVector v2 = tsos.globalDirection();
388  plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
389  }
390  else
391  {
392  // ECAL Endcap
393  Plane::ConstPlanePointer thePlaneNeg =
394  Plane::build(Surface::PositionType(0.,0.,-320.9), rot);
395  tsos = thePropagator->propagate(fts,*thePlaneNeg);
396 
397  if(tsos.isValid())
398  {
399  GlobalPoint p2 = tsos.globalPosition();
400  GlobalVector v2 = tsos.globalDirection();
401  plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
402  }
403  }
404  }
405  }
406 
407  file << "}]";
408 
409  delete theHitAssociator;
410 }
PlotRecTracks(const edm::EventSetup &es_, std::string trackProducer_, std::ofstream &file_)
static std::string getInfo(const DetId &id)
Definition: HitInfo.cc:24
int i
Definition: DBlmapReader.cc:9
void printHelix(const GlobalPoint &p1, const GlobalPoint &p2, const GlobalVector &n2, std::ofstream &outFile, int charge)
Definition: PlotUtils.cc:7
std::string getPixelInfo(const TrackingRecHit *recHit, const std::ostringstream &o, const std::ostringstream &d)
void printRecTracks(const edm::Event &ev)
double outerPy() const
y coordinate of momentum vector at the outermost hit position
Definition: Track.h:73
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:10
T y() const
Definition: PV3DBase.h:62
GlobalPoint globalPosition() const
static PlanePointer build(const PositionType &pos, const RotationType &rot, MediumProperties *mp=0)
Definition: Plane.h:25
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
double outerZ() const
z coordinate of the outermost hit position
Definition: Track.h:81
DataContainer const & measurements() const
Definition: Trajectory.h:203
void printPixelRecHit(const SiPixelRecHit *recHit)
Definition: PlotRecHits.cc:35
CovarianceMatrix outerStateCovariance() const
outermost trajectory state curvilinear errors
Definition: Track.h:51
T z() const
Definition: PV3DBase.h:63
double outerX() const
x coordinate of the outermost hit position
Definition: Track.h:77
double p2[4]
Definition: TauolaWrapper.h:90
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
double outerPz() const
z coordinate of momentum vector at the outermost hit position
Definition: Track.h:75
void printStripRecHit(const SiStripRecHit2D *recHit)
Definition: PlotRecHits.cc:86
Definition: DetId.h:20
T const * product() const
Definition: ESHandle.h:62
bool isValid() const
tuple simHits
Definition: trackerHits.py:16
static CylinderPointer build(const PositionType &pos, const RotationType &rot, Scalar radius, MediumProperties *mp=0)
Definition: Cylinder.h:29
T const * product() const
Definition: Handle.h:74
std::string getStripInfo(const TrackingRecHit *recHit, const std::ostringstream &o, const std::ostringstream &d)
double p1[4]
Definition: TauolaWrapper.h:89
double outerY() const
y coordinate of the outermost hit position
Definition: Track.h:79
int charge() const
track electric charge
Definition: TrackBase.h:113
FreeTrajectoryState getTrajectoryAtOuterPoint(const reco::Track &track)
DetId geographicalId() const
LimitAlgo * algo
Definition: Combine.cc:60
T x() const
Definition: PV3DBase.h:61
virtual LocalPoint localPosition() const =0
const SiStripRecHit2D & originalHit() const
double outerPx() const
x coordinate of momentum vector at the outermost hit position
Definition: Track.h:71
GlobalVector globalDirection() const
Pixel Reconstructed Hit.