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