SteppingHelixPropagator/test/SteppingHelixPropagatorAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackPropagation/SteppingHelixPropagator
0004 // Class:      SteppingHelixPropagatorAnalyzer
0005 //
0006 /**\class SteppingHelixPropagatorAnalyzer
0007    Description: Analyzer of SteppingHelixPropagator performance
0008
0009    Implementation:
0010    Use simTracks and simVertices as initial points. For all  muon PSimHits in the event
0011    extrapolate/propagate from the previous point (starting from a vertex)
0012    to the hit position (detector surface).
0013    Fill an nTuple (could've been an EventProduct) with expected (given by the propagator)
0014    and actual (PSimHits)
0015    positions of a muon in the detector.
0016 */
0017 //
0018 // Original Author:  Vyacheslav Krutelyov
0019 //         Created:  Fri Mar  3 16:01:24 CST 2006
0020 //
0021 //
0022
0023 // system include files
0024 #include <memory>
0025
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034
0035 #include "DataFormats/Common/interface/Handle.h"
0036
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039
0040 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0041 #include "DataFormats/GeometrySurface/interface/Plane.h"
0042
0043 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0044 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0045
0046 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0047 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0048 //#include "Geometry/RPCGeometry/interface/RPCRoll.h"
0049 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0050
0051 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0052
0053 #include "MagneticField/Engine/interface/MagneticField.h"
0054 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0055
0056 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0057
0058 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0059 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0060 #include "SimDataFormats/Track/interface/SimTrack.h"
0061 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0062 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0063 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0064
0065 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0066 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0067 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0068
0069 #include "CLHEP/Random/RandFlat.h"
0070 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0071 //#include "CLHEP/Matrix/DiagMatrix.h"
0072
0073 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0074 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixStateInfo.h"
0075
0076 #include "TFile.h"
0077 #include "TTree.h"
0078
0079 #include <map>
0080
0081 //
0082 // class decleration
0083 //
0084
0085 class SteppingHelixPropagatorAnalyzer : public edm::one::EDAnalyzer<> {
0086 public:
0087   explicit SteppingHelixPropagatorAnalyzer(const edm::ParameterSet&);
0088   ~SteppingHelixPropagatorAnalyzer();
0089
0090   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0091   virtual void endJob();
0092
0093 protected:
0094   struct GlobalSimHit {
0095     const PSimHit* hit;
0096     const Surface* surf;
0097     DetId id;
0098     CLHEP::Hep3Vector r3;
0099     CLHEP::Hep3Vector p3;
0100   };
0101
0102   void loadNtVars(int ind,
0103                   int eType,
0104                   int pStatus,
0105                   int id,  //defs offset: 0 for R, 1*3 for Z and, 2*3 for P
0106                   const CLHEP::Hep3Vector& p3,
0107                   const CLHEP::Hep3Vector& r3,
0108                   const CLHEP::Hep3Vector& p3R,
0109                   const CLHEP::Hep3Vector& r3R,
0110                   int charge,
0111                   const AlgebraicSymMatrix66& cov);
0112
0113   FreeTrajectoryState getFromCLHEP(const CLHEP::Hep3Vector& p3,
0114                                    const CLHEP::Hep3Vector& r3,
0115                                    int charge,
0116                                    const AlgebraicSymMatrix66& cov,
0117                                    const MagneticField* field);
0118   void getFromFTS(const FreeTrajectoryState& fts,
0119                   CLHEP::Hep3Vector& p3,
0120                   CLHEP::Hep3Vector& r3,
0121                   int& charge,
0122                   AlgebraicSymMatrix66& cov);
0123
0124   void addPSimHits(const edm::Event& iEvent,
0125                    const std::string instanceName,
0126                    const edm::ESHandle<GlobalTrackingGeometry>& geom,
0127                    std::vector<SteppingHelixPropagatorAnalyzer::GlobalSimHit>& hits) const;
0128
0129 private:
0130   // ----------member data ---------------------------
0131   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
0132   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propToken;
0133   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> globTrkGeomToken;
0134
0135   TFile* ntFile_;
0136   TTree* tr_;
0137
0138   int nPoints_;
0139   int q_[1000];
0140   int eType_[1000];
0141   int pStatus_[1000][3];
0142   float p3_[1000][9];
0143   float r3_[1000][9];
0144   int id_[1000];
0145   float p3R_[1000][3];
0146   float r3R_[1000][3];
0147   float covFlat_[1000][21];
0148
0149   bool debug_;
0150   int run_;
0151   int event_;
0152
0153   int trkIndOffset_;
0154
0155   bool doneMapping_;
0156
0157   bool radX0CorrectionMode_;
0158
0159   bool testPCAPropagation_;
0160
0161   bool ntupleTkHits_;
0162
0163   bool startFromPrevHit_;
0164
0165   std::string g4SimName_;
0166   edm::InputTag simTracksTag_;
0167   edm::InputTag simVertexesTag_;
0168 };
0169
0170 //
0171 // constructors and destructor
0172 //
0173 SteppingHelixPropagatorAnalyzer::SteppingHelixPropagatorAnalyzer(const edm::ParameterSet& iConfig)
0174     : magFieldToken(esConsumes()),
0175       propToken(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0176       globTrkGeomToken(esConsumes()),
0177       simTracksTag_(iConfig.getParameter<edm::InputTag>("simTracksTag")),
0178       simVertexesTag_(iConfig.getParameter<edm::InputTag>("simVertexesTag")) {
0179   //now do what ever initialization is needed
0180
0181   ntFile_ = new TFile(iConfig.getParameter<std::string>("NtFile").c_str(), "recreate");
0182   tr_ = new TTree("MuProp", "MuProp");
0183   tr_->Branch("nPoints", &nPoints_, "nPoints/I");
0184   tr_->Branch("q", q_, "q[nPoints]/I");
0185   tr_->Branch("pStatus", pStatus_, "pStatus[nPoints][3]/I");
0186   tr_->Branch("p3", p3_, "p3[nPoints][9]/F");
0187   tr_->Branch("r3", r3_, "r3[nPoints][9]/F");
0188   tr_->Branch("id", id_, "id[nPoints]/I");
0189   tr_->Branch("p3R", p3R_, "p3R[nPoints][3]/F");
0190   tr_->Branch("r3R", r3R_, "r3R[nPoints][3]/F");
0191   tr_->Branch("covFlat", covFlat_, "covFlat[nPoints][21]/F");
0192   tr_->Branch("run", &run_, "run/I");
0193   tr_->Branch("event_", &event_, "event/I");
0194
0195   trkIndOffset_ = iConfig.getParameter<int>("trkIndOffset");
0196   debug_ = iConfig.getParameter<bool>("debug");
0197   radX0CorrectionMode_ = iConfig.getParameter<bool>("radX0CorrectionMode");
0198
0199   testPCAPropagation_ = iConfig.getParameter<bool>("testPCAPropagation");
0200
0201   ntupleTkHits_ = iConfig.getParameter<bool>("ntupleTkHits");
0202
0203   startFromPrevHit_ = iConfig.getParameter<bool>("startFromPrevHit");
0204
0205   g4SimName_ = iConfig.getParameter<std::string>("g4SimName");
0206
0207   //initialize the variables
0208   for (int i = 0; i < 1000; ++i) {
0209     q_[i] = 0;
0210     eType_[i] = 0;
0211     for (int j = 0; j < 3; ++j)
0212       pStatus_[i][j] = 0;
0213     for (int j = 0; j < 9; ++j)
0214       p3_[i][j] = 0;
0215     for (int j = 0; j < 9; ++j)
0216       r3_[i][j] = 0;
0217     id_[i] = 0;
0218     for (int j = 0; j < 3; ++j)
0219       p3R_[i][j] = 0;
0220     for (int j = 0; j < 3; ++j)
0221       r3R_[i][j] = 0;
0222     for (int j = 0; j < 21; ++j)
0223       covFlat_[i][j] = 0;
0224   }
0225 }
0226
0227 SteppingHelixPropagatorAnalyzer::~SteppingHelixPropagatorAnalyzer() = default;
0228
0229 //
0230 // member functions
0231 //
0232
0233 // ------------ method called to produce the data  ------------
0234 void SteppingHelixPropagatorAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0235   constexpr char const* metname = "SteppingHelixPropagatorAnalyzer";
0236   (void)metname;
0237   using namespace edm;
0238   const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
0239   const ESHandle<Propagator> shProp = iSetup.getHandle(propToken);
0240   const ESHandle<Propagator> shPropAny = iSetup.getHandle(propToken);
0241   const ESHandle<GlobalTrackingGeometry> geomESH = iSetup.getHandle(globTrkGeomToken);
0242   if (debug_) {
0243     LogTrace(metname) << "Got GlobalTrackingGeometry " << std::endl;
0244   }
0245
0246   run_ = (int)iEvent.id().run();
0247   event_ = (int)iEvent.id().event();
0248   if (debug_) {
0249     LogTrace(metname) << "Begin for run:event ==" << run_ << ":" << event_ << std::endl;
0250   }
0251
0252   const double FPRP_MISMATCH = 150.;
0253   int pStatus = 0;  //1 will be bad
0254
0255   Handle<SimTrackContainer> simTracks;
0256   iEvent.getByLabel<SimTrackContainer>(simTracksTag_, simTracks);
0257   if (!simTracks.isValid()) {
0258     std::cout << "No tracks found" << std::endl;
0259     return;
0260   }
0261   if (debug_) {
0262     LogTrace(metname) << "Got simTracks of size " << simTracks->size() << std::endl;
0263   }
0264
0265   Handle<SimVertexContainer> simVertices;
0266   iEvent.getByLabel<SimVertexContainer>(simVertexesTag_, simVertices);
0267   if (!simVertices.isValid()) {
0268     std::cout << "No tracks found" << std::endl;
0269     return;
0270   }
0271   if (debug_) {
0272     LogTrace(metname) << "Got simVertices of size " << simVertices->size() << std::endl;
0273   }
0274
0275   std::vector<GlobalSimHit> allSimHits;
0276   allSimHits.clear();
0277
0278   addPSimHits(iEvent, "MuonDTHits", geomESH, allSimHits);
0279   addPSimHits(iEvent, "MuonCSCHits", geomESH, allSimHits);
0280   addPSimHits(iEvent, "MuonRPCHits", geomESH, allSimHits);
0281   addPSimHits(iEvent, "TrackerHitsPixelBarrelLowTof", geomESH, allSimHits);
0282   addPSimHits(iEvent, "TrackerHitsPixelEndcapLowTof", geomESH, allSimHits);
0283   addPSimHits(iEvent, "TrackerHitsTIBLowTof", geomESH, allSimHits);
0284   addPSimHits(iEvent, "TrackerHitsTIDLowTof", geomESH, allSimHits);
0285   addPSimHits(iEvent, "TrackerHitsTOBLowTof", geomESH, allSimHits);
0286   addPSimHits(iEvent, "TrackerHitsTECLowTof", geomESH, allSimHits);
0287
0288   SimTrackContainer::const_iterator tracksCI = simTracks->begin();
0289   for (; tracksCI != simTracks->end(); tracksCI++) {
0290     int trkPDG = tracksCI->type();
0291     if (abs(trkPDG) != 13) {
0292       if (debug_) {
0293         LogTrace(metname) << "Skip " << trkPDG << std::endl;
0294       }
0295       continue;
0296     }
0297     CLHEP::Hep3Vector p3T(tracksCI->momentum().x(), tracksCI->momentum().y(), tracksCI->momentum().z());
0298     if (p3T.mag() < 2.)
0299       continue;
0300
0301     int vtxInd = tracksCI->vertIndex();
0302     unsigned int trkInd = tracksCI->genpartIndex() - trkIndOffset_;
0303     CLHEP::Hep3Vector r3T(0., 0., 0.);
0304     if (vtxInd < 0) {
0305       std::cout << "Track with no vertex, defaulting to (0,0,0)" << std::endl;
0306     } else {
0307       r3T = CLHEP::Hep3Vector((*simVertices)[vtxInd].position().x(),
0308                               (*simVertices)[vtxInd].position().y(),
0309                               (*simVertices)[vtxInd].position().z());
0310     }
0311     AlgebraicSymMatrix66 covT = AlgebraicMatrixID();
0312     covT *= 1e-20;  // initialize to sigma=1e-10 .. should get overwhelmed by MULS
0313
0314     CLHEP::Hep3Vector p3F, r3F;  //propagated state
0315     AlgebraicSymMatrix66 covF;
0316     int charge = trkPDG > 0 ? -1 : 1;  //works for muons
0317
0318     nPoints_ = 0;
0319     pStatus = 0;
0320     loadNtVars(nPoints_, 0, pStatus, 0, p3T, r3T, p3T, r3T, charge, covT);
0321     nPoints_++;
0322     FreeTrajectoryState ftsTrack = getFromCLHEP(p3T, r3T, charge, covT, &*bField);
0323     FreeTrajectoryState ftsStart = ftsTrack;
0324     SteppingHelixStateInfo siStart(ftsStart);
0325     SteppingHelixStateInfo siDest;
0326     TrajectoryStateOnSurface tSOSDest;
0327
0328     std::map<double, const GlobalSimHit*> simHitsByDistance;
0329     std::map<double, const GlobalSimHit*> simHitsByTof;
0330     for (std::vector<GlobalSimHit>::const_iterator allHitsCI = allSimHits.begin(); allHitsCI != allSimHits.end();
0331          allHitsCI++) {
0332       if (allHitsCI->hit->trackId() != trkInd)
0333         continue;
0334       if (abs(allHitsCI->hit->particleType()) != 13)
0335         continue;
0336       if (allHitsCI->p3.mag() < 0.5)
0337         continue;
0338
0339       double distance = (allHitsCI->r3 - r3T).mag();
0340       double tof = allHitsCI->hit->timeOfFlight();
0341       simHitsByDistance[distance] = &*allHitsCI;
0342       simHitsByTof[tof] = &*allHitsCI;
0343     }
0344
0345     {  //new scope for timing purposes only
0346       if (testPCAPropagation_) {
0347         FreeTrajectoryState ftsDest;
0348         GlobalPoint pDest1(10., 10., 0.);
0349         GlobalPoint pDest2(10., 10., 10.);
0350         const SteppingHelixPropagator* shPropAnyCPtr = dynamic_cast<const SteppingHelixPropagator*>(&*shPropAny);
0351
0352         ftsDest = shPropAnyCPtr->propagate(ftsStart, pDest1);
0353         std::cout << "----------------------------------------------" << std::endl;
0354         ftsDest = shPropAnyCPtr->propagate(ftsStart, pDest1, pDest2);
0355         std::cout << "----------------------------------------------" << std::endl;
0356       }
0357
0358       //now we are supposed to have a sorted list of hits
0359       std::map<double, const GlobalSimHit*>::const_iterator simHitsCI = simHitsByDistance.begin();
0360       for (; simHitsCI != simHitsByDistance.end(); simHitsCI++) {
0361         const GlobalSimHit* igHit = simHitsCI->second;
0362
0363         if (debug_) {
0364           [[clang::suppress]]
0365           const PSimHit* iHit = simHitsCI->second->hit;
0366           LogTrace(metname) << igHit->id.rawId() << " r3L:" << iHit->localPosition() << " r3G:" << igHit->r3
0367                             << " p3L:" << iHit->momentumAtEntry() << " p3G:" << igHit->p3
0368                             << " pId:" << iHit->particleType() << " tId:" << iHit->trackId() << std::endl;
0369         }
0370
0371         if (debug_) {
0372           LogTrace(metname) << "Will propagate to surface: " << igHit->surf->position() << " "
0373                             << igHit->surf->rotation() << std::endl;
0374         }
0375         if (startFromPrevHit_) {
0376           if (simHitsCI != simHitsByDistance.begin()) {
0377             std::map<double, const GlobalSimHit*>::const_iterator simHitPrevCI = simHitsCI;
0378             simHitPrevCI--;
0379             const GlobalSimHit* gpHit = simHitPrevCI->second;
0380             ftsStart = getFromCLHEP(gpHit->p3, gpHit->r3, charge, covT, &*bField);
0381             siStart = SteppingHelixStateInfo(ftsTrack);
0382           }
0383         }
0384
0385         if (radX0CorrectionMode_) {
0386           const SteppingHelixPropagator* shPropCPtr = dynamic_cast<const SteppingHelixPropagator*>(&*shProp);
0387           shPropCPtr->propagate(siStart, *igHit->surf, siDest);
0388           if (siDest.isValid()) {
0389             siStart = siDest;
0390             siStart.getFreeState(ftsStart);
0391             getFromFTS(ftsStart, p3F, r3F, charge, covF);
0392             pStatus = 0;
0393           } else
0394             pStatus = 1;
0395           if (pStatus == 1 || (r3F - igHit->r3).mag() > FPRP_MISMATCH) {
0396             //start from the beginning if failed with previous
0397             siStart = SteppingHelixStateInfo(ftsTrack);
0398             pStatus = 1;
0399           }
0400         } else {
0401           tSOSDest = shProp->propagate(ftsStart, *igHit->surf);
0402           if (tSOSDest.isValid()) {
0403             ftsStart = *tSOSDest.freeState();
0404             getFromFTS(ftsStart, p3F, r3F, charge, covF);
0405             pStatus = 0;
0406           } else
0407             pStatus = 1;
0408           if (pStatus == 1 || (r3F - igHit->r3).mag() > FPRP_MISMATCH) {
0409             //start from the beginning if failed with previous
0410             ftsStart = ftsTrack;
0411             pStatus = 1;
0412           }
0413         }
0414
0415         if (debug_) {
0416           LogTrace(metname) << "Got to "
0417                             << " r3Prp:" << r3F << " r3Hit:" << igHit->r3 << " p3Prp:" << p3F << " p3Hit:" << igHit->p3
0418                             << " pPrp:" << p3F.mag() << " pHit:" << igHit->p3.mag() << std::endl;
0419         }
0420         loadNtVars(nPoints_, 0, pStatus, igHit->id.rawId(), p3F, r3F, igHit->p3, igHit->r3, charge, covF);
0421         nPoints_++;
0422       }
0423     }
0424     if (tr_)
0425       tr_->Fill();  //fill this track prop info
0426   }
0427 }
0428
0429 // "endJob" is an inherited method that you may implement to do post-EOF processing
0430 // and produce final output.
0431 //
0432 void SteppingHelixPropagatorAnalyzer::endJob() {
0433   ntFile_->cd();
0434   tr_->Write();
0435   delete ntFile_;
0436   ntFile_ = 0;
0437 }
0438
0439 void SteppingHelixPropagatorAnalyzer::loadNtVars(int ind,
0440                                                  int eType,
0441                                                  int pStatus,
0442                                                  int id,
0443                                                  const CLHEP::Hep3Vector& p3,
0444                                                  const CLHEP::Hep3Vector& r3,
0445                                                  const CLHEP::Hep3Vector& p3R,
0446                                                  const CLHEP::Hep3Vector& r3R,
0447                                                  int charge,
0448                                                  const AlgebraicSymMatrix66& cov) {
0449   p3_[ind][eType * 3 + 0] = p3.x();
0450   p3_[ind][eType * 3 + 1] = p3.y();
0451   p3_[ind][eType * 3 + 2] = p3.z();
0452   r3_[ind][eType * 3 + 0] = r3.x();
0453   r3_[ind][eType * 3 + 1] = r3.y();
0454   r3_[ind][eType * 3 + 2] = r3.z();
0455   id_[ind] = id;
0456   p3R_[ind][0] = p3R.x();
0457   p3R_[ind][1] = p3R.y();
0458   p3R_[ind][2] = p3R.z();
0459   r3R_[ind][0] = r3R.x();
0460   r3R_[ind][1] = r3R.y();
0461   r3R_[ind][2] = r3R.z();
0462   int flatInd = 0;
0463   for (int i = 1; i <= cov.kRows; i++)
0464     for (int j = 1; j <= i; j++) {
0465       covFlat_[ind][flatInd] = cov(i - 1, j - 1);
0466       flatInd++;
0467     }
0468   q_[ind] = charge;
0469   eType_[ind] = eType;
0470   pStatus_[ind][eType] = pStatus;
0471 }
0472
0473 FreeTrajectoryState SteppingHelixPropagatorAnalyzer::getFromCLHEP(const CLHEP::Hep3Vector& p3,
0474                                                                   const CLHEP::Hep3Vector& r3,
0475                                                                   int charge,
0476                                                                   const AlgebraicSymMatrix66& cov,
0477                                                                   const MagneticField* field) {
0478   GlobalVector p3GV(p3.x(), p3.y(), p3.z());
0479   GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
0480   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
0481
0482   CartesianTrajectoryError tCov(cov);
0483
0484   return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars);
0485 }
0486
0487 void SteppingHelixPropagatorAnalyzer::getFromFTS(const FreeTrajectoryState& fts,
0488                                                  CLHEP::Hep3Vector& p3,
0489                                                  CLHEP::Hep3Vector& r3,
0490                                                  int& charge,
0491                                                  AlgebraicSymMatrix66& cov) {
0492   GlobalVector p3GV = fts.momentum();
0493   GlobalPoint r3GP = fts.position();
0494
0495   p3.set(p3GV.x(), p3GV.y(), p3GV.z());
0496   r3.set(r3GP.x(), r3GP.y(), r3GP.z());
0497
0498   charge = fts.charge();
0499   cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
0500 }
0501
0502 void SteppingHelixPropagatorAnalyzer ::addPSimHits(
0503     const edm::Event& iEvent,
0504     const std::string instanceName,
0505     const edm::ESHandle<GlobalTrackingGeometry>& geom,
0506     std::vector<SteppingHelixPropagatorAnalyzer::GlobalSimHit>& hits) const {
0507   static std::string metname = "SteppingHelixPropagatorAnalyzer";
0508   edm::Handle<edm::PSimHitContainer> handle;
0509   iEvent.getByLabel(g4SimName_, instanceName, handle);
0510   if (!handle.isValid()) {
0511     std::cout << "No hits found" << std::endl;
0512     return;
0513   }
0514   if (debug_) {
0515     LogTrace(metname) << "Got " << instanceName << " of size " << handle->size() << std::endl;
0516   }
0517
0518   edm::PSimHitContainer::const_iterator pHits_CI = handle->begin();
0519   for (; pHits_CI != handle->end(); pHits_CI++) {
0520     int dtId = pHits_CI->detUnitId();
0521     DetId wId(dtId);
0522
0523     const GeomDet* layer = geom->idToDet(wId);
0524
0525     GlobalSimHit gHit;
0526     gHit.hit = &*pHits_CI;
0527     gHit.surf = &layer->surface();
0528     gHit.id = wId;
0529
0530     //place the hit onto the surface
0531     float dzFrac = gHit.hit->entryPoint().z() / (gHit.hit->exitPoint().z() - gHit.hit->entryPoint().z());
0532     GlobalPoint r3Hit =
0533         gHit.surf->toGlobal(gHit.hit->entryPoint() - dzFrac * (gHit.hit->exitPoint() - gHit.hit->entryPoint()));
0534     gHit.r3.set(r3Hit.x(), r3Hit.y(), r3Hit.z());
0535     GlobalVector p3Hit = gHit.surf->toGlobal(gHit.hit->momentumAtEntry());
0536     gHit.p3.set(p3Hit.x(), p3Hit.y(), p3Hit.z());
0537     hits.push_back(gHit);
0538   }
0539 }
0540
0541 //define this as a plug-in
0542 DEFINE_FWK_MODULE(SteppingHelixPropagatorAnalyzer);