Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-12 07:47:35

0001 //////////////////////////
0002 //  Producer by Anders  //
0003 //     and Emmanuele    //
0004 //    july 2012 @ CU    //
0005 //////////////////////////
0006 
0007 ////////////////////
0008 // FRAMEWORK HEADERS
0009 #include "FWCore/PluginManager/interface/ModuleDef.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 //
0012 #include "FWCore/Framework/interface/one/EDProducer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 //
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 
0023 ///////////////////////
0024 // DATA FORMATS HEADERS
0025 #include "DataFormats/Common/interface/Handle.h"
0026 #include "DataFormats/Common/interface/Ref.h"
0027 //
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0030 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0031 #include "DataFormats/Common/interface/DetSetVector.h"
0032 //
0033 #include "L1Trigger/TrackFindingTracklet/interface/SLHCEvent.h"
0034 #include "L1Trigger/TrackFindingTracklet/interface/L1TStub.h"
0035 
0036 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0037 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0038 #include "SimDataFormats/Track/interface/SimTrack.h"
0039 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0040 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0041 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0042 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0043 #include "SimTracker/TrackTriggerAssociation/interface/TTStubAssociationMap.h"
0044 #include "SimTracker/TrackTriggerAssociation/interface/TTClusterAssociationMap.h"
0045 //
0046 #include "DataFormats/Math/interface/LorentzVector.h"
0047 #include "DataFormats/Math/interface/Vector3D.h"
0048 //
0049 #include "DataFormats/L1TrackTrigger/interface/TTCluster.h"
0050 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0051 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0052 #include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h"
0053 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0054 #include "DataFormats/L1TrackTrigger/interface/TTDTC.h"
0055 #include "L1Trigger/TrackerDTC/interface/Setup.h"
0056 //
0057 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0058 #include "DataFormats/Candidate/interface/Candidate.h"
0059 //
0060 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0061 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0062 #include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
0063 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0064 
0065 ////////////////////////////
0066 // DETECTOR GEOMETRY HEADERS
0067 #include "MagneticField/Engine/interface/MagneticField.h"
0068 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0069 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0070 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0071 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0072 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0073 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0074 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0075 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0076 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0077 //
0078 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0079 
0080 ///////////////
0081 // Tracklet emulation
0082 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0083 #include "L1Trigger/TrackFindingTracklet/interface/Sector.h"
0084 #include "L1Trigger/TrackFindingTracklet/interface/Track.h"
0085 #include "L1Trigger/TrackFindingTracklet/interface/TrackletEventProcessor.h"
0086 
0087 ////////////////
0088 // PHYSICS TOOLS
0089 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0090 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0091 
0092 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0093 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0094 
0095 #include "L1Trigger/TrackTrigger/interface/StubPtConsistency.h"
0096 #include "L1Trigger/TrackTrigger/interface/L1TrackQuality.h"
0097 
0098 //////////////
0099 // STD HEADERS
0100 #include <memory>
0101 #include <string>
0102 #include <iostream>
0103 #include <fstream>
0104 
0105 //////////////
0106 // NAMESPACES
0107 using namespace edm;
0108 using namespace std;
0109 
0110 //////////////////////////////
0111 //                          //
0112 //     CLASS DEFINITION     //
0113 //                          //
0114 //////////////////////////////
0115 
0116 /////////////////////////////////////
0117 // this class is needed to make a map
0118 // between different types of stubs
0119 struct L1TStubCompare {
0120 public:
0121   bool operator()(const trklet::L1TStub& a, const trklet::L1TStub& b) const {
0122     if (a.x() != b.x())
0123       return (b.x() > a.x());
0124     else {
0125       if (a.y() != b.y())
0126         return (b.y() > a.y());
0127       else
0128         return (a.z() > b.z());
0129     }
0130   }
0131 };
0132 
0133 class L1FPGATrackProducer : public edm::one::EDProducer<edm::one::WatchRuns> {
0134 public:
0135   /// Constructor/destructor
0136   explicit L1FPGATrackProducer(const edm::ParameterSet& iConfig);
0137   ~L1FPGATrackProducer() override;
0138 
0139 private:
0140   int eventnum;
0141 
0142   /// Containers of parameters passed by python configuration file
0143   edm::ParameterSet config;
0144 
0145   bool readMoreMcTruth_;
0146 
0147   /// File path for configuration files
0148   edm::FileInPath fitPatternFile;
0149   edm::FileInPath memoryModulesFile;
0150   edm::FileInPath processingModulesFile;
0151   edm::FileInPath wiresFile;
0152 
0153   edm::FileInPath tableTEDFile;
0154   edm::FileInPath tableTREFile;
0155 
0156   string asciiEventOutName_;
0157   std::ofstream asciiEventOut_;
0158 
0159   // settings containing various constants for the tracklet processing
0160   trklet::Settings settings;
0161 
0162   // event processor for the tracklet track finding
0163   trklet::TrackletEventProcessor eventProcessor;
0164 
0165   unsigned int nHelixPar_;
0166   bool extended_;
0167 
0168   bool trackQuality_;
0169   std::unique_ptr<L1TrackQuality> trackQualityModel_;
0170 
0171   std::map<string, vector<int>> dtclayerdisk;
0172 
0173   edm::InputTag MCTruthClusterInputTag;
0174   edm::InputTag MCTruthStubInputTag;
0175   edm::InputTag TrackingParticleInputTag;
0176 
0177   const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0178 
0179   edm::EDGetTokenT<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> ttClusterMCTruthToken_;
0180   edm::EDGetTokenT<std::vector<TrackingParticle>> TrackingParticleToken_;
0181   edm::EDGetTokenT<TTDTC> tokenDTC_;
0182 
0183   // helper class to store DTC configuration
0184   trackerDTC::Setup setup_;
0185 
0186   // Setup token
0187   edm::ESGetToken<trackerDTC::Setup, trackerDTC::SetupRcd> esGetToken_;
0188   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0189 
0190   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tGeomToken_;
0191   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0192 
0193   /// ///////////////// ///
0194   /// MANDATORY METHODS ///
0195   void beginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
0196   void endRun(edm::Run const&, edm::EventSetup const&) override;
0197   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0198 };
0199 
0200 //////////////
0201 // CONSTRUCTOR
0202 L1FPGATrackProducer::L1FPGATrackProducer(edm::ParameterSet const& iConfig)
0203     : config(iConfig),
0204       readMoreMcTruth_(iConfig.getParameter<bool>("readMoreMcTruth")),
0205       MCTruthClusterInputTag(readMoreMcTruth_ ? config.getParameter<edm::InputTag>("MCTruthClusterInputTag")
0206                                               : edm::InputTag()),
0207       MCTruthStubInputTag(readMoreMcTruth_ ? config.getParameter<edm::InputTag>("MCTruthStubInputTag")
0208                                            : edm::InputTag()),
0209       TrackingParticleInputTag(readMoreMcTruth_ ? iConfig.getParameter<edm::InputTag>("TrackingParticleInputTag")
0210                                                 : edm::InputTag()),
0211       bsToken_(consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("BeamSpotSource"))),
0212       tokenDTC_(consumes<TTDTC>(edm::InputTag(iConfig.getParameter<edm::InputTag>("InputTagTTDTC")))),
0213       magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0214       tGeomToken_(esConsumes()),
0215       tTopoToken_(esConsumes()) {
0216   if (readMoreMcTruth_) {
0217     ttClusterMCTruthToken_ = consumes<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>>(MCTruthClusterInputTag);
0218     TrackingParticleToken_ = consumes<std::vector<TrackingParticle>>(TrackingParticleInputTag);
0219   }
0220 
0221   produces<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>("Level1TTTracks").setBranchAlias("Level1TTTracks");
0222 
0223   asciiEventOutName_ = iConfig.getUntrackedParameter<string>("asciiFileName", "");
0224 
0225   fitPatternFile = iConfig.getParameter<edm::FileInPath>("fitPatternFile");
0226   processingModulesFile = iConfig.getParameter<edm::FileInPath>("processingModulesFile");
0227   memoryModulesFile = iConfig.getParameter<edm::FileInPath>("memoryModulesFile");
0228   wiresFile = iConfig.getParameter<edm::FileInPath>("wiresFile");
0229 
0230   extended_ = iConfig.getParameter<bool>("Extended");
0231   nHelixPar_ = iConfig.getParameter<unsigned int>("Hnpar");
0232 
0233   if (extended_) {
0234     tableTEDFile = iConfig.getParameter<edm::FileInPath>("tableTEDFile");
0235     tableTREFile = iConfig.getParameter<edm::FileInPath>("tableTREFile");
0236   }
0237 
0238   // book ES product
0239   esGetToken_ = esConsumes<trackerDTC::Setup, trackerDTC::SetupRcd, edm::Transition::BeginRun>();
0240 
0241   // --------------------------------------------------------------------------------
0242   // set options in Settings based on inputs from configuration files
0243   // --------------------------------------------------------------------------------
0244 
0245   settings.setExtended(extended_);
0246   settings.setNHelixPar(nHelixPar_);
0247 
0248   settings.setFitPatternFile(fitPatternFile.fullPath());
0249   settings.setProcessingModulesFile(processingModulesFile.fullPath());
0250   settings.setMemoryModulesFile(memoryModulesFile.fullPath());
0251   settings.setWiresFile(wiresFile.fullPath());
0252 
0253   if (extended_) {
0254     settings.setTableTEDFile(tableTEDFile.fullPath());
0255     settings.setTableTREFile(tableTREFile.fullPath());
0256 
0257     //FIXME: The TED and TRE tables are currently disabled by default, so we
0258     //need to allow for the additional tracklets that will eventually be
0259     //removed by these tables, once they are finalized
0260     settings.setNbitstrackletindex(10);
0261   }
0262 
0263   eventnum = 0;
0264   if (not asciiEventOutName_.empty()) {
0265     asciiEventOut_.open(asciiEventOutName_.c_str());
0266   }
0267 
0268   if (settings.debugTracklet()) {
0269     edm::LogVerbatim("Tracklet") << "fit pattern :     " << fitPatternFile.fullPath()
0270                                  << "\n process modules : " << processingModulesFile.fullPath()
0271                                  << "\n memory modules :  " << memoryModulesFile.fullPath()
0272                                  << "\n wires          :  " << wiresFile.fullPath();
0273     if (extended_) {
0274       edm::LogVerbatim("Tracklet") << "table_TED    :  " << tableTEDFile.fullPath()
0275                                    << "\n table_TRE    :  " << tableTREFile.fullPath();
0276     }
0277   }
0278 
0279   trackQuality_ = iConfig.getParameter<bool>("TrackQuality");
0280   if (trackQuality_) {
0281     trackQualityModel_ = std::make_unique<L1TrackQuality>(iConfig.getParameter<edm::ParameterSet>("TrackQualityPSet"));
0282   }
0283 }
0284 
0285 /////////////
0286 // DESTRUCTOR
0287 L1FPGATrackProducer::~L1FPGATrackProducer() {
0288   if (asciiEventOut_.is_open()) {
0289     asciiEventOut_.close();
0290   }
0291 }
0292 
0293 ///////END RUN
0294 //
0295 void L1FPGATrackProducer::endRun(const edm::Run& run, const edm::EventSetup& iSetup) {}
0296 
0297 ////////////
0298 // BEGIN JOB
0299 void L1FPGATrackProducer::beginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0300   ////////////////////////
0301   // GET MAGNETIC FIELD //
0302   const MagneticField* theMagneticField = &iSetup.getData(magneticFieldToken_);
0303   double mMagneticFieldStrength = theMagneticField->inTesla(GlobalPoint(0, 0, 0)).z();
0304   settings.setBfield(mMagneticFieldStrength);
0305 
0306   setup_ = iSetup.getData(esGetToken_);
0307 
0308   // initialize the tracklet event processing (this sets all the processing & memory modules, wiring, etc)
0309   eventProcessor.init(settings);
0310 }
0311 
0312 //////////
0313 // PRODUCE
0314 void L1FPGATrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0315   typedef std::map<trklet::L1TStub,
0316                    edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>,
0317                    L1TStubCompare>
0318       stubMapType;
0319   typedef edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0320       TTClusterRef;
0321 
0322   /// Prepare output
0323   auto L1TkTracksForOutput = std::make_unique<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>();
0324 
0325   stubMapType stubMap;
0326 
0327   ////////////
0328   // GET BS //
0329   edm::Handle<reco::BeamSpot> beamSpotHandle;
0330   iEvent.getByToken(bsToken_, beamSpotHandle);
0331   math::XYZPoint bsPosition = beamSpotHandle->position();
0332 
0333   eventnum++;
0334   trklet::SLHCEvent ev;
0335   ev.setEventNum(eventnum);
0336   ev.setIP(bsPosition.x(), bsPosition.y());
0337 
0338   // tracking particles
0339   edm::Handle<std::vector<TrackingParticle>> TrackingParticleHandle;
0340   if (readMoreMcTruth_)
0341     iEvent.getByToken(TrackingParticleToken_, TrackingParticleHandle);
0342 
0343   // tracker topology
0344   const TrackerTopology* const tTopo = &iSetup.getData(tTopoToken_);
0345   const TrackerGeometry* const theTrackerGeom = &iSetup.getData(tGeomToken_);
0346 
0347   ////////////////////////
0348   // GET THE PRIMITIVES //
0349   edm::Handle<TTDTC> handleDTC;
0350   iEvent.getByToken<TTDTC>(tokenDTC_, handleDTC);
0351 
0352   // must be defined for code to compile, even if it's not used unless readMoreMcTruth_ is true
0353   map<edm::Ptr<TrackingParticle>, int> translateTP;
0354 
0355   // MC truth association maps
0356   edm::Handle<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTClusterHandle;
0357   if (readMoreMcTruth_) {
0358     iEvent.getByToken(ttClusterMCTruthToken_, MCTruthTTClusterHandle);
0359 
0360     ////////////////////////////////////////////////
0361     /// LOOP OVER TRACKING PARTICLES & GET SIMTRACKS
0362 
0363     int ntps = 1;  //count from 1 ; 0 will mean invalid
0364 
0365     int this_tp = 0;
0366     for (const auto& iterTP : *TrackingParticleHandle) {
0367       edm::Ptr<TrackingParticle> tp_ptr(TrackingParticleHandle, this_tp);
0368       this_tp++;
0369 
0370       // only keep TPs producing a cluster
0371       if (MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).empty())
0372         continue;
0373 
0374       if (iterTP.g4Tracks().empty()) {
0375         continue;
0376       }
0377 
0378       int sim_eventid = iterTP.g4Tracks().at(0).eventId().event();
0379       int sim_type = iterTP.pdgId();
0380       float sim_pt = iterTP.pt();
0381       float sim_eta = iterTP.eta();
0382       float sim_phi = iterTP.phi();
0383 
0384       float vx = iterTP.vertex().x();
0385       float vy = iterTP.vertex().y();
0386       float vz = iterTP.vertex().z();
0387 
0388       if (sim_pt < 1.0 || std::abs(vz) > 100.0 || hypot(vx, vy) > 50.0)
0389         continue;
0390 
0391       ev.addL1SimTrack(sim_eventid, ntps, sim_type, sim_pt, sim_eta, sim_phi, vx, vy, vz);
0392 
0393       translateTP[tp_ptr] = ntps;
0394       ntps++;
0395 
0396     }  //end loop over TPs
0397 
0398   }  // end if (readMoreMcTruth_)
0399 
0400   /////////////////////////////////
0401   /// READ DTC STUB INFORMATION ///
0402   /////////////////////////////////
0403 
0404   // Process stubs in each region and channel within that region
0405   for (const int& region : handleDTC->tfpRegions()) {
0406     for (const int& channel : handleDTC->tfpChannels()) {
0407       // Get the DTC name form the channel
0408 
0409       static string dtcbasenames[12] = {
0410           "PS10G_1", "PS10G_2", "PS10G_3", "PS10G_4", "PS_1", "PS_2", "2S_1", "2S_2", "2S_3", "2S_4", "2S_5", "2S_6"};
0411 
0412       string dtcname = dtcbasenames[channel % 12];
0413 
0414       if (channel % 24 >= 12)
0415         dtcname = "neg" + dtcname;
0416 
0417       dtcname += (channel < 24) ? "_A" : "_B";
0418 
0419       // Get the stubs from the DTC
0420       const TTDTC::Stream& streamFromDTC{handleDTC->stream(region, channel)};
0421 
0422       // Prepare the DTC stubs for the IR
0423       for (size_t stubIndex = 0; stubIndex < streamFromDTC.size(); ++stubIndex) {
0424         const TTDTC::Frame& stub{streamFromDTC[stubIndex]};
0425 
0426         if (stub.first.isNull()) {
0427           continue;
0428         }
0429 
0430         const GlobalPoint& ttPos = setup_.stubPos(stub.first);
0431 
0432         //Get the 2 bits for the layercode
0433         string layerword = stub.second.to_string().substr(61, 2);
0434         unsigned int layercode = 2 * (layerword[0] - '0') + layerword[1] - '0';
0435         assert(layercode < 4);
0436 
0437         //translation from the two bit layercode to the layer/disk number of each of the
0438         //12 channels (dtcs)
0439         static int layerdisktab[12][4] = {{0, 6, 8, 10},
0440                                           {0, 7, 9, -1},
0441                                           {1, 7, -1, -1},
0442                                           {6, 8, 10, -1},
0443                                           {2, 7, -1, -1},
0444                                           {2, 9, -1, -1},
0445                                           {3, 4, -1, -1},
0446                                           {4, -1, -1, -1},
0447                                           {5, -1, -1, -1},
0448                                           {5, 8, -1, -1},
0449                                           {6, 9, -1, -1},
0450                                           {7, 10, -1, -1}};
0451 
0452         int layerdisk = layerdisktab[channel % 12][layercode];
0453         assert(layerdisk != -1);
0454 
0455         //Get the 36 bit word - skip the lowest 3 buts (status and layer code)
0456         constexpr int DTCLinkWordSize = 64;
0457         constexpr int StubWordSize = 36;
0458         constexpr int LayerandStatusCodeSize = 3;
0459         string stubword =
0460             stub.second.to_string().substr(DTCLinkWordSize - StubWordSize - LayerandStatusCodeSize, StubWordSize);
0461         string stubwordhex = "";
0462 
0463         //Loop over the 9 words in the 36 bit stub word
0464         for (unsigned int i = 0; i < 9; i++) {
0465           bitset<4> bits(stubword.substr(i * 4, 4));
0466           ulong val = bits.to_ulong();
0467           stubwordhex += ((val < 10) ? ('0' + val) : ('A' + val - 10));
0468         }
0469 
0470         /// Get the Inner and Outer TTCluster
0471         edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0472             innerCluster = stub.first->clusterRef(0);
0473         edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0474             outerCluster = stub.first->clusterRef(1);
0475 
0476         // -----------------------------------------------------
0477         // check module orientation, if flipped, need to store that information for track fit
0478         // -----------------------------------------------------
0479 
0480         const DetId innerDetId = innerCluster->getDetId();
0481         const GeomDetUnit* det_inner = theTrackerGeom->idToDetUnit(innerDetId);
0482         const auto* theGeomDet_inner = dynamic_cast<const PixelGeomDetUnit*>(det_inner);
0483         const PixelTopology* topol_inner = dynamic_cast<const PixelTopology*>(&(theGeomDet_inner->specificTopology()));
0484 
0485         MeasurementPoint coords_inner = innerCluster->findAverageLocalCoordinatesCentered();
0486         LocalPoint clustlp_inner = topol_inner->localPosition(coords_inner);
0487         GlobalPoint posStub_inner = theGeomDet_inner->surface().toGlobal(clustlp_inner);
0488 
0489         const DetId outerDetId = outerCluster->getDetId();
0490         const GeomDetUnit* det_outer = theTrackerGeom->idToDetUnit(outerDetId);
0491         const auto* theGeomDet_outer = dynamic_cast<const PixelGeomDetUnit*>(det_outer);
0492         const PixelTopology* topol_outer = dynamic_cast<const PixelTopology*>(&(theGeomDet_outer->specificTopology()));
0493 
0494         MeasurementPoint coords_outer = outerCluster->findAverageLocalCoordinatesCentered();
0495         LocalPoint clustlp_outer = topol_outer->localPosition(coords_outer);
0496         GlobalPoint posStub_outer = theGeomDet_outer->surface().toGlobal(clustlp_outer);
0497 
0498         bool isFlipped = (posStub_outer.mag() < posStub_inner.mag());
0499 
0500         vector<int> assocTPs;
0501 
0502         for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0503 
0504           const TTClusterRef& ttClusterRef = stub.first->clusterRef(iClus);
0505 
0506           // Now identify all TP's contributing to either cluster in stub.
0507           vector<edm::Ptr<TrackingParticle>> vecTpPtr = MCTruthTTClusterHandle->findTrackingParticlePtrs(ttClusterRef);
0508 
0509           for (const edm::Ptr<TrackingParticle>& tpPtr : vecTpPtr) {
0510             if (translateTP.find(tpPtr) != translateTP.end()) {
0511               if (iClus == 0) {
0512                 assocTPs.push_back(translateTP.at(tpPtr));
0513               } else {
0514                 assocTPs.push_back(-translateTP.at(tpPtr));
0515               }
0516               // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0517             } else {
0518               assocTPs.push_back(0);
0519             }
0520           }
0521         }
0522 
0523         double stubbend = stub.first->bendFE();  //stub.first->rawBend()
0524         if (ttPos.z() < -120) {
0525           stubbend = -stubbend;
0526         }
0527 
0528         ev.addStub(dtcname,
0529                    region,
0530                    layerdisk,
0531                    stubwordhex,
0532                    setup_.psModule(setup_.dtcId(region, channel)),
0533                    isFlipped,
0534                    ttPos.x(),
0535                    ttPos.y(),
0536                    ttPos.z(),
0537                    stubbend,
0538                    stub.first->innerClusterPosition(),
0539                    assocTPs);
0540 
0541         const trklet::L1TStub& lastStub = ev.lastStub();
0542         stubMap[lastStub] = stub.first;
0543       }
0544     }
0545   }
0546 
0547   //////////////////////////
0548   // NOW RUN THE L1 tracking
0549 
0550   if (!asciiEventOutName_.empty()) {
0551     ev.write(asciiEventOut_);
0552   }
0553 
0554   const std::vector<trklet::Track>& tracks = eventProcessor.tracks();
0555 
0556   // this performs the actual tracklet event processing
0557   eventProcessor.event(ev);
0558 
0559   int ntracks = 0;
0560 
0561   for (const auto& track : tracks) {
0562     if (track.duplicate())
0563       continue;
0564 
0565     ntracks++;
0566 
0567     // this is where we create the TTTrack object
0568     double tmp_rinv = track.rinv(settings);
0569     double tmp_phi = track.phi0(settings);
0570     double tmp_tanL = track.tanL(settings);
0571     double tmp_z0 = track.z0(settings);
0572     double tmp_d0 = track.d0(settings);
0573     double tmp_chi2rphi = track.chisqrphi();
0574     double tmp_chi2rz = track.chisqrz();
0575     unsigned int tmp_hit = track.hitpattern();
0576 
0577     TTTrack<Ref_Phase2TrackerDigi_> aTrack(tmp_rinv,
0578                                            tmp_phi,
0579                                            tmp_tanL,
0580                                            tmp_z0,
0581                                            tmp_d0,
0582                                            tmp_chi2rphi,
0583                                            tmp_chi2rz,
0584                                            0,
0585                                            0,
0586                                            0,
0587                                            tmp_hit,
0588                                            settings.nHelixPar(),
0589                                            settings.bfield());
0590 
0591     unsigned int trksector = track.sector();
0592     unsigned int trkseed = (unsigned int)abs(track.seed());
0593 
0594     aTrack.setPhiSector(trksector);
0595     aTrack.setTrackSeedType(trkseed);
0596 
0597     const vector<trklet::L1TStub>& stubptrs = track.stubs();
0598     vector<trklet::L1TStub> stubs;
0599 
0600     stubs.reserve(stubptrs.size());
0601     for (const auto& stubptr : stubptrs) {
0602       stubs.push_back(stubptr);
0603     }
0604 
0605     stubMapType::const_iterator it;
0606     for (const auto& itstubs : stubs) {
0607       it = stubMap.find(itstubs);
0608       if (it != stubMap.end()) {
0609         aTrack.addStubRef(it->second);
0610       } else {
0611         // could not find stub in stub map
0612       }
0613     }
0614 
0615     // pt consistency
0616     aTrack.setStubPtConsistency(
0617         StubPtConsistency::getConsistency(aTrack, theTrackerGeom, tTopo, settings.bfield(), settings.nHelixPar()));
0618 
0619     // set TTTrack word
0620     aTrack.setTrackWordBits();
0621 
0622     if (trackQuality_) {
0623       trackQualityModel_->setL1TrackQuality(aTrack);
0624     }
0625 
0626     // test track word
0627     //aTrack.testTrackWordBits();
0628 
0629     L1TkTracksForOutput->push_back(aTrack);
0630   }
0631 
0632   iEvent.put(std::move(L1TkTracksForOutput), "Level1TTTracks");
0633 
0634 }  /// End of produce()
0635 
0636 // ///////////////////////////
0637 // // DEFINE THIS AS A PLUG-IN
0638 DEFINE_FWK_MODULE(L1FPGATrackProducer);