CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/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/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   // Get tracker geometry
00052   edm::ESHandle<TrackerGeometry> trackerHandle;
00053   es.get<TrackerDigiGeometryRecord>().get(trackerHandle);
00054   theTracker = trackerHandle.product();
00055 
00056   // Get magnetic field
00057   edm::ESHandle<MagneticField> magField;
00058   es.get<IdealMagneticFieldRecord>().get(magField);
00059   theMagField = magField.product();
00060 
00061   // Get propagator
00062   edm::ESHandle<Propagator> thePropagatorHandle;
00063   es.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",
00064                                           thePropagatorHandle);
00065   thePropagator = thePropagatorHandle.product();
00066 
00067   // KFTrajectoryFitter
00068 /*
00069   edm::ESHandle<TrajectoryFitter> theFitterHandle;
00070   es.get<TrackingComponentsRecord>().get("KFTrajectoryFitter", theFitterHandle);
00071   theFitter = theFitterHandle.product();
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 cerr << " track[" << i << "] " << recTrack->chi2() << " " << it->chiSquared() << std::endl;
00234 
00235 */
00236 //theFitter->fit(*it);
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         // Pixel
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         // Strip
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   // Trajectory
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     // Ecal
00369     GlobalPoint p1 = (*(meas.rend()-1)).forwardPredictedState().globalPosition();
00370     FreeTrajectoryState fts = getTrajectoryAtOuterPoint(*recTrack);
00371     Surface::RotationType rot;
00372     TrajectoryStateOnSurface tsos;
00373 
00374     // Ecal Barrel
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       // ECAL Endcap
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         // ECAL Endcap
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 }