00001 #include "QCDAnalysis/ChargedHadronSpectra/interface/PlotRecTracks.h"
00002 #include "QCDAnalysis/ChargedHadronSpectra/interface/PlotUtils.h"
00003 #include "QCDAnalysis/ChargedHadronSpectra/interface/PlotRecHits.h"
00004 #include "QCDAnalysis/ChargedHadronSpectra/interface/HitInfo.h"
00005
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "FWCore/Framework/interface/ESHandle.h"
00010 #include "FWCore/Framework/interface/Event.h"
00011
00012 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00013 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00015
00016 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00017 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00018 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00019 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00020 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00021
00022 #include "DataFormats/TrackReco/interface/Track.h"
00023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00024 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00025
00026 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00027 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
00028
00029 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00030 #include "MagneticField/Engine/interface/MagneticField.h"
00031
00032 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00033 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00034
00035 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00036
00037 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
00038
00039 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00040
00041 using namespace std;
00042
00043
00044 PlotRecTracks::PlotRecTracks
00045 (const edm::EventSetup& es_, std::string trackProducer_,
00046 ofstream& file_) : es(es_), trackProducer(trackProducer_), file(file_)
00047 {
00048
00049 edm::ESHandle<TrackerGeometry> trackerHandle;
00050 es.get<TrackerDigiGeometryRecord>().get(trackerHandle);
00051 theTracker = trackerHandle.product();
00052
00053
00054 edm::ESHandle<MagneticField> magField;
00055 es.get<IdealMagneticFieldRecord>().get(magField);
00056 theMagField = magField.product();
00057
00058
00059 edm::ESHandle<Propagator> thePropagatorHandle;
00060 es.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",
00061 thePropagatorHandle);
00062 thePropagator = thePropagatorHandle.product();
00063
00064
00065
00066
00067
00068
00069
00070 }
00071
00072
00073 PlotRecTracks::~PlotRecTracks()
00074 {
00075 }
00076
00077
00078 string PlotRecTracks::getPixelInfo(const TrackingRecHit* recHit, const TrackerTopology* tTopo, const ostringstream& o, const ostringstream& d)
00079 {
00080 const SiPixelRecHit* pixelRecHit =
00081 dynamic_cast<const SiPixelRecHit *>(recHit);
00082
00083 DetId id = recHit->geographicalId();
00084 LocalPoint lpos = recHit->localPosition();
00085 GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
00086
00087 SiPixelRecHit::ClusterRef const& cluster = pixelRecHit->cluster();
00088 std::vector<SiPixelCluster::Pixel> pixels = cluster->pixels();
00089
00090 std::vector<PSimHit> simHits = theHitAssociator->associateHit(*recHit);
00091
00092 std::string info = ", Text[StyleForm[\"" + o.str() + "\", URL->\"Track " + o.str()
00093 + d.str();
00094
00095 {
00096 ostringstream o;
00097 o << "simTrack (trackId=" << simHits[0].trackId()
00098 << " pid=" << simHits[0].particleType()
00099 << " proc=" << simHits[0].processType() << ")";
00100
00101 info += " | " + o.str();
00102 }
00103
00104 {
00105 ostringstream o;
00106 o << theTracker->idToDet(id)->subDetector();
00107
00108 info += " | " + o.str();
00109 }
00110
00111 info += HitInfo::getInfo(*recHit, tTopo);
00112
00113 info += "\"]";
00114
00115 {
00116 ostringstream o;
00117 o << ", {" << p.x() << "," << p.y() << ",(" << p.z() << "-zs)*mz},"
00118 << " {" << 0 << "," << -1 << "}]";
00119
00120 info += o.str();
00121 }
00122
00123 return info;
00124 }
00125
00126
00127 string PlotRecTracks::getStripInfo
00128 (const TrackingRecHit* recHit, const TrackerTopology* tTopo, const ostringstream& o, const ostringstream& d)
00129 {
00130 DetId id = recHit->geographicalId();
00131 LocalPoint lpos = recHit->localPosition();
00132 GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
00133
00134 std::vector<PSimHit> simHits = theHitAssociator->associateHit(*recHit);
00135
00136 std::string info = ", Text[StyleForm[\"" + o.str() + "\", URL->\"Track " + o.str() + d.str();
00137 const SiStripMatchedRecHit2D* stripMatchedRecHit =
00138 dynamic_cast<const SiStripMatchedRecHit2D *>(recHit);
00139 const ProjectedSiStripRecHit2D* stripProjectedRecHit =
00140 dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit);
00141 const SiStripRecHit2D* stripRecHit =
00142 dynamic_cast<const SiStripRecHit2D *>(recHit);
00143
00144 {
00145 ostringstream o;
00146 o << "simTrackId=" << simHits[0].trackId() ;
00147
00148 info += " | " + o.str();
00149 }
00150
00151 {
00152 ostringstream o;
00153 o << theTracker->idToDet(id)->subDetector();
00154
00155 info += " | " + o.str();
00156 }
00157
00158 info += HitInfo::getInfo(*recHit, tTopo);
00159
00160 if(stripMatchedRecHit != 0) info += " matched";
00161 if(stripProjectedRecHit != 0) info += " projected";
00162 if(stripRecHit != 0) info += " single";
00163
00164
00165 info += "\"]";
00166
00167 {
00168 ostringstream o;
00169 o << ", {" << p.x() << "," << p.y() << ",(" << p.z() << "-zs)*mz},"
00170 << " {" << 0 << "," << -1 << "}]";
00171
00172 info += o.str();
00173 }
00174
00175 return info;
00176 }
00177
00178
00179 FreeTrajectoryState PlotRecTracks::getTrajectoryAtOuterPoint
00180 (const reco::Track& track)
00181 {
00182 GlobalPoint position(track.outerX(),
00183 track.outerY(),
00184 track.outerZ());
00185
00186 GlobalVector momentum(track.outerPx(),
00187 track.outerPy(),
00188 track.outerPz());
00189
00190 GlobalTrajectoryParameters gtp(position,momentum,
00191 track.charge(),theMagField);
00192
00193 FreeTrajectoryState fts(gtp,track.outerStateCovariance());
00194
00195 return fts;
00196 }
00197
00198
00199 void PlotRecTracks::printRecTracks(const edm::Event& ev, const edm::EventSetup& es)
00200 {
00201
00202 edm::ESHandle<TrackerTopology> tTopoHandle;
00203 es.get<IdealGeometryRecord>().get(tTopoHandle);
00204 const TrackerTopology* const tTopo = tTopoHandle.product();
00205
00206 theHitAssociator = new TrackerHitAssociator(ev);
00207
00208 file << ", If[rt, {AbsolutePointSize[6]";
00209
00210 edm::Handle<reco::TrackCollection> recTrackHandle;
00211 ev.getByLabel(trackProducer, recTrackHandle);
00212
00213 edm::Handle<std::vector<Trajectory> > trajectoryHandle;
00214 ev.getByLabel(trackProducer, trajectoryHandle);
00215
00216 const reco::TrackCollection* recTracks = recTrackHandle.product();
00217 const std::vector<Trajectory>* trajectories = trajectoryHandle.product();
00218
00219 edm::LogVerbatim("MinBiasTracking")
00220 << " [EventPlotter] recTracks (" << trackProducer << ") "
00221 << recTracks->size();
00222
00223 PlotRecHits theRecHits(es,file);
00224
00225 file << ", RGBColor[0,0,0.4]";
00226 reco::TrackCollection::const_iterator recTrack = recTracks->begin();
00227
00228 int i = 0;
00229 for(std::vector<Trajectory>::const_iterator it = trajectories->begin();
00230 it!= trajectories->end();
00231 it++, i++, recTrack++)
00232 {
00233
00234
00235
00236
00237
00238
00239 ostringstream o; o << i;
00240
00241 ostringstream d; d << fixed << std::setprecision(2)
00242 << " | d0=" << recTrack->d0() << " cm"
00243 << " | z0=" << recTrack->dz() << " cm"
00244 << " | pt=" << recTrack->pt() << " GeV/c";
00245
00246 const Trajectory* trajectory = &(*it);
00247
00248 for(std::vector<TrajectoryMeasurement>::const_reverse_iterator
00249 meas = trajectory->measurements().rbegin();
00250 meas!= trajectory->measurements().rend(); meas++)
00251 {
00252 const TrackingRecHit* recHit = meas->recHit()->hit();
00253
00254 if(recHit->isValid())
00255 {
00256 DetId id = recHit->geographicalId();
00257
00258 if(theTracker->idToDet(id)->subDetector() ==
00259 GeomDetEnumerators::PixelBarrel ||
00260 theTracker->idToDet(id)->subDetector() ==
00261 GeomDetEnumerators::PixelEndcap)
00262 {
00263
00264 const SiPixelRecHit* pixelRecHit =
00265 dynamic_cast<const SiPixelRecHit *>(recHit);
00266
00267 if(pixelRecHit != 0)
00268 {
00269 theRecHits.printPixelRecHit(pixelRecHit);
00270 file << getPixelInfo(recHit, tTopo, o, d);
00271 }
00272 }
00273 else
00274 {
00275
00276 const SiStripMatchedRecHit2D* stripMatchedRecHit =
00277 dynamic_cast<const SiStripMatchedRecHit2D *>(recHit);
00278 const ProjectedSiStripRecHit2D* stripProjectedRecHit =
00279 dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit);
00280 const SiStripRecHit2D* stripRecHit =
00281 dynamic_cast<const SiStripRecHit2D *>(recHit);
00282
00283 if(stripMatchedRecHit != 0)
00284 {
00285 auto m = stripMatchedRecHit->monoHit();
00286 auto s = stripMatchedRecHit->stereoHit();
00287 theRecHits.printStripRecHit(&m);
00288 theRecHits.printStripRecHit(&s);
00289
00290 DetId id = stripMatchedRecHit->geographicalId();
00291 LocalPoint lpos = stripMatchedRecHit->localPosition();
00292 GlobalPoint p = theTracker->idToDet(id)->toGlobal(lpos);
00293
00294 file << ", Point[{"<< p.x()<<","<<p.y()<<",("<<p.z()<<"-zs)*mz}]" << std::endl;
00295 }
00296
00297 if(stripProjectedRecHit != 0)
00298 theRecHits.printStripRecHit(&(stripProjectedRecHit->originalHit()));
00299
00300 if(stripRecHit != 0)
00301 theRecHits.printStripRecHit(stripRecHit);
00302
00303 if(stripMatchedRecHit != 0 ||
00304 stripProjectedRecHit != 0 ||
00305 stripRecHit != 0)
00306 file << getStripInfo(recHit, tTopo, o,d);
00307 }
00308 }
00309 }
00310 }
00311
00312 PlotUtils plotUtils;
00313
00314
00315 recTrack = recTracks->begin();
00316
00317 for(std::vector<Trajectory>::const_iterator it = trajectories->begin();
00318 it!= trajectories->end();
00319 it++, recTrack++)
00320 {
00321 int algo;
00322 switch(recTrack->algo())
00323 {
00324 case reco::TrackBase::iter1: algo = 0; break;
00325 case reco::TrackBase::iter2: algo = 1; break;
00326 case reco::TrackBase::iter3: algo = 2; break;
00327 default: algo = 0;
00328 }
00329
00330 if(algo == 0) file << ", RGBColor[1,0,0]";
00331 if(algo == 1) file << ", RGBColor[0.2,0.6,0.2]";
00332 if(algo == 2) file << ", RGBColor[0.2,0.2,0.6]";
00333
00334 std::vector<TrajectoryMeasurement> meas = it->measurements();
00335
00336 for(std::vector<TrajectoryMeasurement>::reverse_iterator im = meas.rbegin();
00337 im!= meas.rend(); im++)
00338 {
00339 if(im == meas.rbegin())
00340 {
00341 GlobalPoint p2 = (*(im )).updatedState().globalPosition();
00342 GlobalVector v2 = (*(im )).updatedState().globalDirection();
00343 GlobalPoint p1 = GlobalPoint(recTrack->vertex().x(),
00344 recTrack->vertex().y(),
00345 recTrack->vertex().z());
00346
00347 plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
00348 }
00349
00350 if(im+1 != meas.rend())
00351 {
00352 GlobalPoint p1 = (*(im )).updatedState().globalPosition();
00353 GlobalPoint p2 = (*(im+1)).updatedState().globalPosition();
00354 GlobalVector v2 = (*(im+1)).updatedState().globalDirection();
00355
00356 plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
00357 }
00358 }
00359
00360
00361 GlobalPoint p1 = (*(meas.rend()-1)).forwardPredictedState().globalPosition();
00362 FreeTrajectoryState fts = getTrajectoryAtOuterPoint(*recTrack);
00363 Surface::RotationType rot;
00364 TrajectoryStateOnSurface tsos;
00365
00366
00367 Cylinder::ConstCylinderPointer theCylinder =
00368 Cylinder::build(129.f, Surface::PositionType(0.,0.,0), rot);
00369 tsos = thePropagator->propagate(fts,*theCylinder);
00370
00371 if(tsos.isValid() &&
00372 fabs(tsos.globalPosition().z()) < 320.9)
00373 {
00374 GlobalPoint p2 = tsos.globalPosition();
00375 GlobalVector v2 = tsos.globalDirection();
00376 plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
00377 }
00378 else
00379 {
00380
00381 Plane::ConstPlanePointer thePlanePos =
00382 Plane::build(Surface::PositionType(0.,0.,320.9), rot);
00383 tsos = thePropagator->propagate(fts,*thePlanePos);
00384
00385 if(tsos.isValid())
00386 {
00387 GlobalPoint p2 = tsos.globalPosition();
00388 GlobalVector v2 = tsos.globalDirection();
00389 plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
00390 }
00391 else
00392 {
00393
00394 Plane::ConstPlanePointer thePlaneNeg =
00395 Plane::build(Surface::PositionType(0.,0.,-320.9), rot);
00396 tsos = thePropagator->propagate(fts,*thePlaneNeg);
00397
00398 if(tsos.isValid())
00399 {
00400 GlobalPoint p2 = tsos.globalPosition();
00401 GlobalVector v2 = tsos.globalDirection();
00402 plotUtils.printHelix(p1,p2,v2, file, recTrack->charge());
00403 }
00404 }
00405 }
00406 }
00407
00408 file << "}]";
00409
00410 delete theHitAssociator;
00411 }