Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:33:31 

0001 // -*- C++ -*-
0002 //
0003 // Package:    trackJet/JetCoreMCtruthSeedGenerator
0004 // Class:      JetCoreMCtruthSeedGenerator
0005 //
0006 /**\class JetCoreMCtruthSeedGenerator JetCoreMCtruthSeedGenerator.cc trackJet/JetCoreMCtruthSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.cc
0007  Description: [one line class summary]
0008  Implementation:
0009      [Notes on implementation]
0010 */
0011 //
0012 // Original Author:  Valerio Bertacchi
0013 //         Created:  Mon, 18 Dec 2017 16:35:04 GMT
0014 //
0015 //
0016 
0017 // system include files
0018 
0019 #include <memory>
0020 
0021 // user include files
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/one/EDProducer.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/ESHandle.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 
0030 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0031 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0032 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0033 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0034 
0035 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0036 #include "DataFormats/Common/interface/Handle.h"
0037 #include "DataFormats/Common/interface/DetSetVector.h"
0038 #include "DataFormats/Common/interface/DetSet.h"
0039 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0040 #include "DataFormats/VertexReco/interface/Vertex.h"
0041 #include "DataFormats/TrackReco/interface/Track.h"
0042 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0043 #include "DataFormats/JetReco/interface/Jet.h"
0044 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0045 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0046 #include "DataFormats/Math/interface/Point3D.h"
0047 #include "DataFormats/Math/interface/Vector3D.h"
0048 #include "DataFormats/Candidate/interface/Candidate.h"
0049 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0050 
0051 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0052 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0053 
0054 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0055 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0056 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0057 
0058 #include "MagneticField/Engine/interface/MagneticField.h"
0059 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0060 
0061 #include "FWCore/ServiceRegistry/interface/Service.h"
0062 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0063 
0064 #include "SimDataFormats/Track/interface/SimTrack.h"
0065 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0066 
0067 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0068 
0069 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0070 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0071 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0072 
0073 //
0074 // class declaration
0075 //
0076 
0077 class JetCoreMCtruthSeedGenerator : public edm::one::EDProducer<edm::one::SharedResources> {
0078 public:
0079   explicit JetCoreMCtruthSeedGenerator(const edm::ParameterSet&);
0080   ~JetCoreMCtruthSeedGenerator() override;
0081 
0082   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0083 
0084   double jetPt_;
0085   double jetEta_;
0086   double pitchX_ = 0.01;              //100 um (pixel pitch in X)
0087   double pitchY_ = 0.015;             //150 um (pixel pitch in Y)
0088   static constexpr int jetDimX = 30;  //pixel dimension of NN window on layer2
0089   static constexpr int jetDimY = 30;  //pixel dimension of NN window on layer2
0090   bool inclusiveConeSeed_ =
0091       true;  //true= fill tracks in a cone of deltaR_, false=fill tracks which have SimHit on globDet
0092 
0093 private:
0094   void beginJob() override;
0095   void produce(edm::Event&, const edm::EventSetup&) override;
0096   void endJob() override;
0097 
0098   // ----------member data ---------------------------
0099   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> geomEsToken_;
0100   const edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> pixelCPEEsToken_;
0101   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoEsToken_;
0102   edm::ESHandle<GlobalTrackingGeometry> geometry_;
0103 
0104   edm::EDGetTokenT<std::vector<reco::Vertex>> vertices_;
0105   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> pixelClusters_;
0106   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> inputPixelClusters_;
0107   edm::EDGetTokenT<edm::View<reco::Candidate>> cores_;
0108   edm::EDGetTokenT<std::vector<SimTrack>> simtracksToken_;
0109   edm::EDGetTokenT<std::vector<SimVertex>> simvertexToken_;
0110   edm::EDGetTokenT<std::vector<PSimHit>> PSimHitToken_;
0111   edm::Handle<std::vector<PSimHit>> simhits_;
0112 
0113   double ptMin_;
0114   double deltaR_;
0115   double chargeFracMin_;
0116   double centralMIPCharge_;
0117 
0118   std::pair<bool, Basic3DVector<float>> findIntersection(const GlobalVector&,
0119                                                          const reco::Candidate::Point&,
0120                                                          const GeomDet*);
0121 
0122   const GeomDet* DetectorSelector(int,
0123                                   const reco::Candidate&,
0124                                   GlobalVector,
0125                                   const reco::Vertex&,
0126                                   const TrackerTopology* const,
0127                                   const edmNew::DetSetVector<SiPixelCluster>&);
0128 
0129   std::vector<GlobalVector> splittedClusterDirections(
0130       const reco::Candidate&,
0131       const TrackerTopology* const,
0132       const PixelClusterParameterEstimator*,
0133       const reco::Vertex&,
0134       int,
0135       const edmNew::DetSetVector<SiPixelCluster>&);  //if not working,: args=2 auto
0136 
0137   std::vector<PSimHit> coreHitsFilling(std::vector<PSimHit>,
0138                                        const GeomDet*,
0139                                        GlobalVector,
0140                                        const reco::Vertex&);  //if not working,: args=0 auto
0141 
0142   std::pair<std::vector<SimTrack>, std::vector<SimVertex>> coreTracksFilling(
0143       std::vector<PSimHit>,
0144       const std::vector<SimTrack>,
0145       const std::vector<SimVertex>);  //if not working,: args=1,2 auto
0146 
0147   std::vector<std::array<double, 5>> seedParFilling(std::pair<std::vector<SimTrack>, std::vector<SimVertex>>,
0148                                                     const GeomDet*,
0149                                                     const reco::Candidate&);
0150 
0151   std::pair<std::vector<SimTrack>, std::vector<SimVertex>> coreTracksFillingDeltaR(
0152       const std::vector<SimTrack>,
0153       const std::vector<SimVertex>,
0154       const GeomDet*,
0155       const reco::Candidate&,
0156       const reco::Vertex&);  //if not working,: args=0,1 auto
0157 };
0158 
0159 JetCoreMCtruthSeedGenerator::JetCoreMCtruthSeedGenerator(const edm::ParameterSet& iConfig)
0160     : geomEsToken_(esConsumes()),
0161       pixelCPEEsToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCPE")))),
0162       tTopoEsToken_(esConsumes()),
0163       vertices_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0164       pixelClusters_(
0165           consumes<edmNew::DetSetVector<SiPixelCluster>>(iConfig.getParameter<edm::InputTag>("pixelClusters"))),
0166       cores_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("cores"))),
0167       simtracksToken_(consumes<std::vector<SimTrack>>(iConfig.getParameter<edm::InputTag>("simTracks"))),
0168       simvertexToken_(consumes<std::vector<SimVertex>>(iConfig.getParameter<edm::InputTag>("simVertex"))),
0169       PSimHitToken_(consumes<std::vector<PSimHit>>(iConfig.getParameter<edm::InputTag>("simHit"))),
0170       ptMin_(iConfig.getParameter<double>("ptMin")),
0171       deltaR_(iConfig.getParameter<double>("deltaR")),
0172       chargeFracMin_(iConfig.getParameter<double>("chargeFractionMin")),
0173       centralMIPCharge_(iConfig.getParameter<double>("centralMIPCharge")) {
0174   produces<TrajectorySeedCollection>();
0175   produces<reco::TrackCollection>();
0176 }
0177 
0178 JetCoreMCtruthSeedGenerator::~JetCoreMCtruthSeedGenerator() {}
0179 
0180 void JetCoreMCtruthSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0181   auto result = std::make_unique<TrajectorySeedCollection>();
0182   auto resultTracks = std::make_unique<reco::TrackCollection>();
0183 
0184   using namespace edm;
0185   using namespace reco;
0186 
0187   geometry_ = iSetup.getHandle(geomEsToken_);
0188 
0189   const auto& inputPixelClusters_ = iEvent.get(pixelClusters_);
0190   const auto& simtracksVector = iEvent.get(simtracksToken_);
0191   const auto& simvertexVector = iEvent.get(simvertexToken_);
0192   const auto& simhits_ = iEvent.get(PSimHitToken_);
0193   const auto& vertices = iEvent.get(vertices_);
0194   const auto& cores = iEvent.get(cores_);
0195 
0196   const PixelClusterParameterEstimator* pixelCPE = &iSetup.getData(pixelCPEEsToken_);
0197   const TrackerTopology* const tTopo = &iSetup.getData(tTopoEsToken_);
0198 
0199   auto output = std::make_unique<edmNew::DetSetVector<SiPixelCluster>>();
0200 
0201   for (const auto& jet : cores) {  //jet loop
0202 
0203     if (jet.pt() > ptMin_) {
0204       std::set<long long int> ids;
0205       const reco::Vertex& jetVertex = vertices[0];
0206 
0207       std::vector<GlobalVector> splitClustDirSet =
0208           splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_);
0209       if (splitClustDirSet.empty()) {  //if layer 1 is broken find direcitons on layer 2
0210         splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_);
0211       }
0212       if (inclusiveConeSeed_)
0213         splitClustDirSet.clear();
0214       splitClustDirSet.emplace_back(GlobalVector(jet.px(), jet.py(), jet.pz()));
0215 
0216       for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) {
0217         GlobalVector bigClustDir = splitClustDirSet[cc];
0218 
0219         jetEta_ = jet.eta();
0220         jetPt_ = jet.pt();
0221 
0222         const auto& jetVert = jetVertex;  //trackInfo filling
0223 
0224         std::vector<PSimHit> goodSimHit;
0225 
0226         const GeomDet* globDet = DetectorSelector(
0227             2, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_);  //select detector mostly hitten by the jet
0228 
0229         if (globDet == nullptr)
0230           continue;
0231 
0232         std::pair<std::vector<SimTrack>, std::vector<SimVertex>> goodSimTkVx;
0233 
0234         if (inclusiveConeSeed_) {
0235           goodSimTkVx = JetCoreMCtruthSeedGenerator::coreTracksFillingDeltaR(
0236               simtracksVector, simvertexVector, globDet, jet, jetVert);
0237         } else {
0238           std::vector<PSimHit> goodSimHit =
0239               JetCoreMCtruthSeedGenerator::coreHitsFilling(simhits_, globDet, bigClustDir, jetVertex);
0240           goodSimTkVx = JetCoreMCtruthSeedGenerator::coreTracksFilling(goodSimHit, simtracksVector, simvertexVector);
0241         }
0242         edm::LogInfo("PerfectSeeder") << "seed number in deltaR cone =" << goodSimTkVx.first.size();
0243 
0244         std::vector<std::array<double, 5>> seedVector =
0245             JetCoreMCtruthSeedGenerator::seedParFilling(goodSimTkVx, globDet, jet);
0246         edm::LogInfo("PerfectSeeder") << "seedVector.size()=" << seedVector.size();
0247 
0248         for (uint tk = 0; tk < seedVector.size(); tk++) {
0249           for (int pp = 0; pp < 5; pp++) {
0250             edm::LogInfo("PerfectSeeder")
0251                 << "seed " << tk << ", int par " << pp << "=" << seedVector[tk][pp] << std::endl;
0252           }
0253           LocalPoint localSeedPoint = LocalPoint(seedVector[tk][0], seedVector[tk][1], 0);
0254           double track_theta = 2 * std::atan(std::exp(-seedVector[tk][2]));
0255           double track_phi = seedVector[tk][3];
0256           double pt = 1. / seedVector[tk][4];
0257 
0258           double normdirR = pt / sin(track_theta);
0259           const GlobalVector globSeedDir(
0260               GlobalVector::Polar(Geom::Theta<double>(track_theta), Geom::Phi<double>(track_phi), normdirR));
0261           LocalVector localSeedDir = globDet->surface().toLocal(globSeedDir);
0262 
0263           int64_t seedid = (int64_t(localSeedPoint.x() * 200.) << 0) + (int64_t(localSeedPoint.y() * 200.) << 16) +
0264                            (int64_t(seedVector[tk][2] * 400.) << 32) + (int64_t(track_phi * 400.) << 48);
0265           if (ids.count(seedid) != 0) {
0266             edm::LogInfo("PerfectSeeder") << "seed not removed with DeepCore cleaner";
0267           }
0268           ids.insert(seedid);
0269 
0270           //Covariance matrix, currently the hadrcoded variances = NN residuals width (see documentation of DeepCore)
0271           //in general: if are not compared with DeepCore but another algo-->to state-of-the art errors
0272           //The "perfect seeds" has no intrinsic error, but the CTF needs errors to propagate...
0273           float em[15] = {
0274               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};  // (see LocalTrajectoryError for details), order as follow:
0275           em[0] = 0.15 * 0.15;                               // q/pt
0276           em[2] = 0.5e-5;                                    // dxdz
0277           em[5] = 0.5e-5;                                    // dydz
0278           em[9] = 2e-5;                                      // x
0279           em[14] = 2e-5;                                     // y
0280           long int detId = globDet->geographicalId();
0281           LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1));
0282           result->emplace_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0),
0283                                               edm::OwnVector<TrackingRecHit>(),
0284                                               PropagationDirection::alongMomentum));
0285 
0286           GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint);
0287           reco::Track::CovarianceMatrix mm;
0288           resultTracks->emplace_back(
0289               reco::Track(1,
0290                           1,
0291                           reco::Track::Point(globalSeedPoint.x(), globalSeedPoint.y(), globalSeedPoint.z()),
0292                           reco::Track::Vector(globSeedDir.x(), globSeedDir.y(), globSeedDir.z()),
0293                           1,
0294                           mm));
0295           edm::LogInfo("PerfectSeeder") << "seed " << tk << ", out,  pt=" << pt << ", eta=" << globSeedDir.eta()
0296                                         << ", phi=" << globSeedDir.phi() << std::endl;
0297         }
0298 
0299       }  //bigcluster
0300     }  //jet > pt
0301   }  //jet
0302   iEvent.put(std::move(result));
0303   iEvent.put(std::move(resultTracks));
0304 }
0305 
0306 std::pair<bool, Basic3DVector<float>> JetCoreMCtruthSeedGenerator::findIntersection(
0307     const GlobalVector& dir, const reco::Candidate::Point& vertex, const GeomDet* det) {
0308   StraightLinePlaneCrossing vertexPlane(Basic3DVector<float>(vertex.x(), vertex.y(), vertex.z()),
0309                                         Basic3DVector<float>(dir.x(), dir.y(), dir.z()));
0310 
0311   std::pair<bool, Basic3DVector<float>> pos = vertexPlane.position(det->specificSurface());
0312 
0313   return pos;
0314 }
0315 
0316 const GeomDet* JetCoreMCtruthSeedGenerator::DetectorSelector(int llay,
0317                                                              const reco::Candidate& jet,
0318                                                              GlobalVector jetDir,
0319                                                              const reco::Vertex& jetVertex,
0320                                                              const TrackerTopology* const tTopo,
0321                                                              const edmNew::DetSetVector<SiPixelCluster>& clusters) {
0322   struct trkNumCompare {
0323     bool operator()(std::pair<int, const GeomDet*> x, std::pair<int, const GeomDet*> y) const {
0324       return x.first > y.first;
0325     }
0326   };
0327   std::set<std::pair<int, const GeomDet*>, trkNumCompare> track4detSet;
0328 
0329   double minDist = 0.0;
0330   const GeomDet* output = nullptr;
0331   for (const auto& detset : clusters) {
0332     auto aClusterID = detset.id();
0333     if (DetId(aClusterID).subdetId() != 1)
0334       continue;
0335     const GeomDet* det = geometry_->idToDet(aClusterID);
0336     int lay = tTopo->layer(det->geographicalId());
0337     if (lay != llay)
0338       continue;
0339     std::pair<bool, Basic3DVector<float>> interPair =
0340         findIntersection(jetDir, (reco::Candidate::Point)jetVertex.position(), det);
0341     if (interPair.first == false)
0342       continue;
0343     Basic3DVector<float> inter = interPair.second;
0344     auto localInter = det->specificSurface().toLocal((GlobalPoint)inter);
0345     if ((minDist == 0.0 || std::abs(localInter.x()) < minDist) && std::abs(localInter.y()) < 3.35) {
0346       minDist = std::abs(localInter.x());
0347       output = det;
0348     }
0349   }  //detset
0350   return output;
0351 }
0352 
0353 std::vector<GlobalVector> JetCoreMCtruthSeedGenerator::splittedClusterDirections(
0354     const reco::Candidate& jet,
0355     const TrackerTopology* const tTopo,
0356     const PixelClusterParameterEstimator* pixelCPE,
0357     const reco::Vertex& jetVertex,
0358     int layer,
0359     const edmNew::DetSetVector<SiPixelCluster>& clusters) {
0360   std::vector<GlobalVector> clustDirs;
0361   for (const auto& detset_int : clusters) {
0362     const GeomDet* det_int = geometry_->idToDet(detset_int.id());
0363     int lay = tTopo->layer(det_int->geographicalId());
0364     if (lay != layer)
0365       continue;  //NB: saved bigclusetr on all the layers!!
0366     auto detUnit = *geometry_->idToDetUnit(detset_int.id());
0367     for (const auto& aCluster : detset_int) {
0368       GlobalPoint clustPos = det_int->surface().toGlobal(pixelCPE->localParametersV(aCluster, detUnit)[0].first);
0369       GlobalPoint vertexPos(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z());
0370       GlobalVector clusterDir = clustPos - vertexPos;
0371       GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
0372       if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
0373         clustDirs.emplace_back(clusterDir);
0374       }
0375     }
0376   }
0377   return clustDirs;
0378 }
0379 
0380 std::vector<PSimHit> JetCoreMCtruthSeedGenerator::coreHitsFilling(std::vector<PSimHit> simhits,
0381                                                                   const GeomDet* globDet,
0382                                                                   GlobalVector bigClustDir,
0383                                                                   const reco::Vertex& jetVertex) {
0384   std::vector<PSimHit> goodSimHit;
0385   for (const auto& sh : simhits) {
0386     const GeomDet* det = geometry_->idToDet(sh.detUnitId());
0387     if (det != globDet)
0388       continue;
0389     std::pair<bool, Basic3DVector<float>> interPair =
0390         findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), det);
0391     if (interPair.first == false)
0392       continue;
0393     Basic3DVector<float> inter = interPair.second;
0394     auto localInter = det->specificSurface().toLocal((GlobalPoint)inter);
0395 
0396     if (std::abs((sh.localPosition()).x() - localInter.x()) / pitchX_ <= jetDimX / 2 &&
0397         std::abs((sh.localPosition()).y() - localInter.y()) / pitchY_ <= jetDimY / 2) {
0398       goodSimHit.emplace_back(sh);
0399     }
0400   }
0401   return goodSimHit;
0402 }
0403 
0404 std::pair<std::vector<SimTrack>, std::vector<SimVertex>> JetCoreMCtruthSeedGenerator::coreTracksFilling(
0405     std::vector<PSimHit> goodSimHit,
0406     const std::vector<SimTrack> simtracksVector,
0407     const std::vector<SimVertex> simvertexVector) {
0408   std::vector<SimTrack> goodSimTrk;
0409   std::vector<SimVertex> goodSimVtx;
0410 
0411   for (uint j = 0; j < simtracksVector.size(); j++) {
0412     for (std::vector<PSimHit>::const_iterator it = goodSimHit.begin(); it != goodSimHit.end(); ++it) {
0413       SimTrack st = simtracksVector[j];
0414       if (st.trackId() == (*it).trackId()) {
0415         for (uint v = 0; v < simvertexVector.size(); v++) {
0416           SimVertex sv = simvertexVector[v];
0417           if ((int)sv.vertexId() == (int)st.vertIndex()) {
0418             goodSimTrk.emplace_back(st);
0419             goodSimVtx.emplace_back(sv);
0420           }
0421         }
0422       }
0423     }
0424   }
0425   std::pair<std::vector<SimTrack>, std::vector<SimVertex>> output(goodSimTrk, goodSimVtx);
0426   return output;
0427 }
0428 
0429 std::pair<std::vector<SimTrack>, std::vector<SimVertex>> JetCoreMCtruthSeedGenerator::coreTracksFillingDeltaR(
0430     const std::vector<SimTrack> simtracksVector,
0431     const std::vector<SimVertex> simvertexVector,
0432     const GeomDet* globDet,
0433     const reco::Candidate& jet,
0434     const reco::Vertex& jetVertex) {
0435   std::vector<SimTrack> goodSimTrk;
0436   std::vector<SimVertex> goodSimVtx;
0437 
0438   GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
0439 
0440   for (uint j = 0; j < simtracksVector.size(); j++) {
0441     SimTrack st = simtracksVector[j];
0442     GlobalVector trkDir(st.momentum().Px(), st.momentum().Py(), st.momentum().Pz());
0443     if (st.charge() == 0)
0444       continue;
0445     if (Geom::deltaR(jetDir, trkDir) < deltaR_) {
0446       for (uint v = 0; v < simvertexVector.size(); v++) {
0447         SimVertex sv = simvertexVector[v];
0448         if ((int)sv.vertexId() == (int)st.vertIndex()) {
0449           goodSimTrk.emplace_back(st);
0450           goodSimVtx.emplace_back(sv);
0451         }
0452       }
0453     }
0454   }
0455   std::pair<std::vector<SimTrack>, std::vector<SimVertex>> output(goodSimTrk, goodSimVtx);
0456   return output;
0457 }
0458 
0459 std::vector<std::array<double, 5>> JetCoreMCtruthSeedGenerator::seedParFilling(
0460     std::pair<std::vector<SimTrack>, std::vector<SimVertex>> goodSimTkVx,
0461     const GeomDet* globDet,
0462     const reco::Candidate& jet) {
0463   std::vector<std::array<double, 5>> output;
0464   std::vector<SimTrack> goodSimTrk = goodSimTkVx.first;
0465   std::vector<SimVertex> goodSimVtx = goodSimTkVx.second;
0466 
0467   edm::LogInfo("PerfectSeeder") << "goodSimTrk size" << goodSimTrk.size();
0468   for (uint j = 0; j < goodSimTrk.size(); j++) {
0469     SimTrack st = goodSimTrk[j];
0470     SimVertex sv = goodSimVtx[j];
0471     GlobalVector trkMom(st.momentum().x(), st.momentum().y(), st.momentum().z());
0472     GlobalPoint trkPos(sv.position().x(), sv.position().y(), sv.position().z());
0473     edm::LogInfo("PerfectSeeder") << "seed " << j << ", very int pt" << st.momentum().Pt()
0474                                   << ", eta=" << st.momentum().Eta() << ", phi=" << st.momentum().Phi()
0475                                   << "------ internal point=" << trkMom.x() << "," << trkMom.y() << "," << trkMom.z()
0476                                   << "," << trkPos.x() << "," << trkPos.y() << "," << trkPos.z() << std::endl;
0477 
0478     std::pair<bool, Basic3DVector<float>> trkInterPair;
0479     trkInterPair = findIntersection(trkMom, (reco::Candidate::Point)trkPos, globDet);
0480     if (trkInterPair.first == false) {
0481       GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
0482       continue;
0483     }
0484     Basic3DVector<float> trkInter = trkInterPair.second;
0485 
0486     auto localTrkInter = globDet->specificSurface().toLocal((GlobalPoint)trkInter);  //trkInter->trkPos if par at vertex
0487     std::array<double, 5> tkPar{
0488         {localTrkInter.x(), localTrkInter.y(), st.momentum().Eta(), st.momentum().Phi(), 1 / st.momentum().Pt()}};
0489     output.emplace_back(tkPar);
0490   }
0491   return output;
0492 }
0493 
0494 // ------------ method called once each job just before starting event loop  ------------
0495 void JetCoreMCtruthSeedGenerator::beginJob() {}
0496 
0497 // ------------ method called once each job just after ending the event loop  ------------
0498 void JetCoreMCtruthSeedGenerator::endJob() {}
0499 
0500 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0501 void JetCoreMCtruthSeedGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0502   edm::ParameterSetDescription desc;
0503   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0504   desc.add<edm::InputTag>("pixelClusters", edm::InputTag("siPixelClustersPreSplitting"));
0505   desc.add<edm::InputTag>("cores", edm::InputTag("jetsForCoreTracking"));
0506   desc.add<double>("ptMin", 300);
0507   desc.add<double>("deltaR", 0.3);
0508   desc.add<double>("chargeFractionMin", 18000.0);
0509   desc.add<edm::InputTag>("simTracks", edm::InputTag("g4SimHits"));
0510   desc.add<edm::InputTag>("simVertex", edm::InputTag("g4SimHits"));
0511   desc.add<edm::InputTag>("simHit", edm::InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof"));
0512   desc.add<double>("centralMIPCharge", 2.);
0513   desc.add<std::string>("pixelCPE", "PixelCPEGeneric");
0514   descriptions.add("JetCoreMCtruthSeedGenerator", desc);
0515 }
0516 
0517 DEFINE_FWK_MODULE(JetCoreMCtruthSeedGenerator);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!