CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/QCDAnalysis/ChargedHadronSpectra/src/PlotRecTracks.cc

Go to the documentation of this file.
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   // Get tracker geometry
00049   edm::ESHandle<TrackerGeometry> trackerHandle;
00050   es.get<TrackerDigiGeometryRecord>().get(trackerHandle);
00051   theTracker = trackerHandle.product();
00052 
00053   // Get magnetic field
00054   edm::ESHandle<MagneticField> magField;
00055   es.get<IdealMagneticFieldRecord>().get(magField);
00056   theMagField = magField.product();
00057 
00058   // Get propagator
00059   edm::ESHandle<Propagator> thePropagatorHandle;
00060   es.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",
00061                                           thePropagatorHandle);
00062   thePropagator = thePropagatorHandle.product();
00063 
00064   // KFTrajectoryFitter
00065 /*
00066   edm::ESHandle<TrajectoryFitter> theFitterHandle;
00067   es.get<TrackingComponentsRecord>().get("KFTrajectoryFitter", theFitterHandle);
00068   theFitter = theFitterHandle.product();
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   //Retrieve tracker topology from geometry
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 cerr << " track[" << i << "] " << recTrack->chi2() << " " << it->chiSquared() << std::endl;
00235 
00236 */
00237 //theFitter->fit(*it);
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         // Pixel
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         // Strip
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   // Trajectory
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     // Ecal
00361     GlobalPoint p1 = (*(meas.rend()-1)).forwardPredictedState().globalPosition();
00362     FreeTrajectoryState fts = getTrajectoryAtOuterPoint(*recTrack);
00363     Surface::RotationType rot;
00364     TrajectoryStateOnSurface tsos;
00365 
00366     // Ecal Barrel
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       // ECAL Endcap
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         // ECAL Endcap
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 }