TrackFindingTracklet/plugins/ProducerKFout.cc

0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002 #include "FWCore/Framework/interface/Run.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Utilities/interface/EDGetToken.h"
0007 #include "FWCore/Utilities/interface/EDPutToken.h"
0008 #include "FWCore/Utilities/interface/ESGetToken.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012
0013 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0014 #include "L1Trigger/TrackerTFP/interface/DataFormats.h"
0015 #include "L1Trigger/TrackTrigger/interface/L1TrackQuality.h"
0016
0017 #include <string>
0018 #include <numeric>
0019
0020 using namespace std;
0021 using namespace edm;
0022 using namespace trackerTFP;
0023 using namespace tt;
0024
0025 namespace trklet {
0026
0027   /*! \class  trklet::ProducerKFout
0028    *  \brief  Converts KF output into tttrack collection and TFP output
0029    *  A bit accurate emulation of the track transformation, the
0030    *  eta routing and splitting of the 96-bit track words into 64-bit
0031    *  packets. Also run is a bit accurate emulation of the track quality
0032    *  BDT, whose output is also added to the track word.
0033    *  \author Christopher Brown
0034    *  \date   2021, Aug
0035    *  \update 2024, June by Claire Savard
0036    */
0037   class ProducerKFout : public stream::EDProducer<> {
0038   public:
0039     explicit ProducerKFout(const ParameterSet&);
0040     ~ProducerKFout() override {}
0041
0042   private:
0043     void beginRun(const Run&, const EventSetup&) override;
0044     void produce(Event&, const EventSetup&) override;
0045     void endJob() {}
0046
0047     // ED input token of kf stubs
0048     EDGetTokenT<StreamsStub> edGetTokenStubs_;
0049     // ED input token of kf tracks
0050     EDGetTokenT<StreamsTrack> edGetTokenTracks_;
0051     // ED output token for accepted kfout tracks
0052     EDPutTokenT<StreamsTrack> edPutTokenAccepted_;
0053     // ED output token for TTTracks
0054     EDPutTokenT<TTTracks> edPutTokenTTTracks_;
0055     // ED output token for truncated kfout tracks
0056     EDPutTokenT<StreamsTrack> edPutTokenLost_;
0057     // Setup token
0058     ESGetToken<Setup, SetupRcd> esGetTokenSetup_;
0059     // DataFormats token
0060     ESGetToken<DataFormats, DataFormatsRcd> esGetTokenDataFormats_;
0061     // configuration
0062     ParameterSet iConfig_;
0063     // helper class to store configurations
0064     const Setup* setup_;
0065     // helper class to extract structured data from tt::Frames
0066     const DataFormats* dataFormats_;
0067     // Bins for dPhi/dZ use to create weight LUT
0068     vector<double> dPhiBins_;
0069     vector<double> dZBins_;
0070
0071     std::unique_ptr<L1TrackQuality> trackQualityModel_;
0072     double tqTanlScale_;
0073     double tqZ0Scale_;
0074     static constexpr double ap_fixed_rescale = 32.0;
0075
0076     // For convenience and keeping readable code, accessed many times
0077     int numWorkers_;
0078     int partialTrackWordBits_;
0079
0080     // Helper function to convert floating value to bin
0081     template <typename T>
0082     unsigned int digitise(const T& bins, double value, double factor) {
0083       unsigned int bin = 0;
0084       for (unsigned int i = 0; i < bins.size() - 1; i++) {
0085         if (value * factor > bins[i] && value * factor <= bins[i + 1])
0086           break;
0087         bin++;
0088       }
0089       return bin;
0090     }
0091   };
0092
0093   ProducerKFout::ProducerKFout(const ParameterSet& iConfig) : iConfig_(iConfig) {
0094     const string& labelKF = iConfig.getParameter<string>("LabelKF");
0095     const string& branchStubs = iConfig.getParameter<string>("BranchAcceptedStubs");
0096     const string& branchTracks = iConfig.getParameter<string>("BranchAcceptedTracks");
0097     const string& branchTTTracks = iConfig.getParameter<string>("BranchAcceptedTTTracks");
0098     const string& branchLost = iConfig.getParameter<string>("BranchLostTracks");
0099     // book in- and output ED products
0100     edGetTokenStubs_ = consumes<StreamsStub>(InputTag(labelKF, branchStubs));
0101     edGetTokenTracks_ = consumes<StreamsTrack>(InputTag(labelKF, branchTracks));
0102     edPutTokenAccepted_ = produces<StreamsTrack>(branchTracks);
0103     edPutTokenTTTracks_ = produces<TTTracks>(branchTTTracks);
0104     edPutTokenLost_ = produces<StreamsTrack>(branchLost);
0105     // book ES products
0106     esGetTokenSetup_ = esConsumes<Setup, SetupRcd, Transition::BeginRun>();
0107     esGetTokenDataFormats_ = esConsumes<DataFormats, DataFormatsRcd, Transition::BeginRun>();
0108     // initial ES products
0109     setup_ = nullptr;
0110     dataFormats_ = nullptr;
0111
0112     trackQualityModel_ = std::make_unique<L1TrackQuality>(iConfig.getParameter<edm::ParameterSet>("TrackQualityPSet"));
0113     edm::ParameterSet trackQualityPSset = iConfig.getParameter<edm::ParameterSet>("TrackQualityPSet");
0114     tqTanlScale_ = trackQualityPSset.getParameter<double>("tqemu_TanlScale");
0115     tqZ0Scale_ = trackQualityPSset.getParameter<double>("tqemu_Z0Scale");
0116   }
0117
0118   void ProducerKFout::beginRun(const Run& iRun, const EventSetup& iSetup) {
0119     // helper class to store configurations
0120     setup_ = &iSetup.getData(esGetTokenSetup_);
0121     if (!setup_->configurationSupported())
0122       return;
0123     // check process history if desired
0124     if (iConfig_.getParameter<bool>("CheckHistory"))
0125       setup_->checkHistory(iRun.processHistory());
0126     // helper class to extract structured data from tt::Frames
0127     dataFormats_ = &iSetup.getData(esGetTokenDataFormats_);
0128
0129     // Calculate 1/dz**2 and 1/dphi**2 bins for v0 and v1 weightings
0130     float temp_dphi = 0.0;
0131     float temp_dz = 0.0;
0132     for (int i = 0;
0133          i < pow(2, dataFormats_->width(Variable::dPhi, Process::kfin)) / pow(2, setup_->weightBinFraction());
0134          i++) {
0135       temp_dphi =
0136           pow(dataFormats_->base(Variable::dPhi, Process::kfin) * (i + 1) * pow(2, setup_->weightBinFraction()), -2);
0137       temp_dphi = temp_dphi / setup_->dphiTruncation();
0138       temp_dphi = std::floor(temp_dphi);
0139       dPhiBins_.push_back(temp_dphi * setup_->dphiTruncation());
0140     }
0141     for (int i = 0; i < pow(2, dataFormats_->width(Variable::dZ, Process::kfin)) / pow(2, setup_->weightBinFraction());
0142          i++) {
0143       temp_dz =
0144           pow(dataFormats_->base(Variable::dZ, Process::kfin) * (i + 1) * pow(2, setup_->weightBinFraction()), -2);
0145       temp_dz = temp_dz * setup_->dzTruncation();
0146       temp_dz = std::ceil(temp_dz);
0147       dZBins_.push_back(temp_dz / setup_->dzTruncation());
0148     }
0149     numWorkers_ = setup_->kfNumWorker();
0150     partialTrackWordBits_ = TTBV::S_ / 2;
0151   }
0152
0153   void ProducerKFout::produce(Event& iEvent, const EventSetup& iSetup) {
0154     // empty KFout product
0155     StreamsTrack accepted(setup_->numRegions() * setup_->tfpNumChannel());
0156     StreamsTrack lost(setup_->numRegions() * setup_->tfpNumChannel());
0157     // read in KF Product and produce KFout product
0158     if (setup_->configurationSupported()) {
0159       Handle<StreamsStub> handleStubs;
0160       iEvent.getByToken<StreamsStub>(edGetTokenStubs_, handleStubs);
0161       const StreamsStub& streamsStubs = *handleStubs.product();
0162       Handle<StreamsTrack> handleTracks;
0163       iEvent.getByToken<StreamsTrack>(edGetTokenTracks_, handleTracks);
0164       const StreamsTrack& streamsTracks = *handleTracks.product();
0165
0166       // Setup KFout track collection
0167       TrackKFOutSAPtrCollection KFoutTracks;
0168
0169       // Setup containers for track quality
0170       float tempTQMVAPreSig = 0.0;
0171       // Due to ap_fixed implementation in CMSSW this 10,5 must be specified at compile time, TODO make this a changeable parameter
0172       std::vector<ap_fixed<10, 5>> trackQuality_inputs = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
0173
0174       // calculate track quality and fill TTTracks
0175       TTTracks ttTracks;
0176       int nTracks(0);
0177       for (const StreamTrack& stream : streamsTracks)
0178         nTracks += accumulate(stream.begin(), stream.end(), 0, [](int sum, const FrameTrack& frame) {
0179           return sum + (frame.first.isNonnull() ? 1 : 0);
0180         });
0181       ttTracks.reserve(nTracks);
0182       for (int iLink = 0; iLink < (int)streamsTracks.size(); iLink++) {
0183         for (int iTrack = 0; iTrack < (int)streamsTracks[iLink].size(); iTrack++) {
0184           const auto& track = streamsTracks[iLink].at(iTrack);
0185           TrackKF inTrack(track, dataFormats_);
0186
0187           double temp_z0 = inTrack.zT() - ((inTrack.cot() * setup_->chosenRofZ()));
0188           // Correction to Phi calcuation depending if +ve/-ve phi sector
0189           const double baseSectorCorr = inTrack.sectorPhi() ? -setup_->baseSector() : setup_->baseSector();
0190           double temp_phi0 = inTrack.phiT() - ((inTrack.inv2R()) * setup_->hybridChosenRofPhi()) + baseSectorCorr;
0191           double temp_tanL = inTrack.cotGlobal();
0192
0193           TTBV hitPattern(0, setup_->numLayers());
0194           double tempchi2rphi = 0;
0195           double tempchi2rz = 0;
0196           int temp_nstub = 0;
0197           int temp_ninterior = 0;
0198           bool counter = false;
0199
0200           vector<StubKF> stubs;
0201           stubs.reserve(setup_->numLayers());
0202           for (int iStub = 0; iStub < setup_->numLayers(); iStub++) {
0203             const auto& stub = streamsStubs[setup_->numLayers() * iLink + iStub].at(iTrack);
0204             StubKF inStub(stub, dataFormats_, iStub);
0205             if (stub.first.isNonnull())
0206               stubs.emplace_back(stub, dataFormats_, iStub);
0207
0208             if (!stub.first.isNonnull()) {
0209               if (counter)
0210                 temp_ninterior += 1;
0211               continue;
0212             }
0213
0214             counter = true;
0215
0216             hitPattern.set(iStub);
0217             temp_nstub += 1;
0218             double phiSquared = pow(inStub.phi(), 2);
0219             double zSquared = pow(inStub.z(), 2);
0220
0221             double tempv0 = dPhiBins_[(inStub.dPhi() / (dataFormats_->base(Variable::dPhi, Process::kfin) *
0222                                                         pow(2, setup_->weightBinFraction())))];
0223             double tempv1 = dZBins_[(
0224                 inStub.dZ() / (dataFormats_->base(Variable::dZ, Process::kfin) * pow(2, setup_->weightBinFraction())))];
0225
0226             double tempRphi = phiSquared * tempv0;
0227             double tempRz = zSquared * tempv1;
0228
0229             tempchi2rphi += tempRphi;
0230             tempchi2rz += tempRz;
0231           }  // Iterate over track stubs
0232
0233           // Create bit vectors for each output, including digitisation of chi2
0234           // TODO implement extraMVA, bendChi2, d0
0235           TTBV trackValid(1, TTTrack_TrackWord::TrackBitWidths::kValidSize, false);
0236           TTBV extraMVA(0, TTTrack_TrackWord::TrackBitWidths::kMVAOtherSize, false);
0237           TTBV bendChi2(0, TTTrack_TrackWord::TrackBitWidths::kBendChi2Size, false);
0238           TTBV chi2rphi(digitise(TTTrack_TrackWord::chi2RPhiBins, tempchi2rphi, (double)setup_->kfoutchi2rphiConv()),
0239                         TTTrack_TrackWord::TrackBitWidths::kChi2RPhiSize,
0240                         false);
0241           TTBV chi2rz(digitise(TTTrack_TrackWord::chi2RZBins, tempchi2rz, (double)setup_->kfoutchi2rzConv()),
0242                       TTTrack_TrackWord::TrackBitWidths::kChi2RZSize,
0243                       false);
0244           TTBV d0(0, TTTrack_TrackWord::TrackBitWidths::kD0Size, false);
0245           TTBV z0(
0246               temp_z0, dataFormats_->base(Variable::zT, Process::kf), TTTrack_TrackWord::TrackBitWidths::kZ0Size, true);
0247           TTBV tanL(temp_tanL,
0248                     dataFormats_->base(Variable::cot, Process::kf),
0249                     TTTrack_TrackWord::TrackBitWidths::kTanlSize,
0250                     true);
0251           TTBV phi0(temp_phi0,
0252                     dataFormats_->base(Variable::phiT, Process::kf),
0253                     TTTrack_TrackWord::TrackBitWidths::kPhiSize,
0254                     true);
0255           TTBV invR(-inTrack.inv2R(),
0256                     dataFormats_->base(Variable::inv2R, Process::kf),
0257                     TTTrack_TrackWord::TrackBitWidths::kRinvSize + 1,
0258                     true);
0259           invR.resize(TTTrack_TrackWord::TrackBitWidths::kRinvSize);
0260
0261           // conversion to tttrack to calculate bendchi2
0262           // temporary fix for MVA1 while bendchi2 not implemented
0263           TTTrack temp_tttrack = inTrack.ttTrack(stubs);
0264           double tempbendchi2 = temp_tttrack.chi2BendRed();
0265
0266           // Create input vector for BDT
0267           trackQuality_inputs = {
0268               (std::trunc(tanL.val() / tqTanlScale_)) / ap_fixed_rescale,
0269               (std::trunc(z0.val() / tqZ0Scale_)) / ap_fixed_rescale,
0270               digitise(TTTrack_TrackWord::bendChi2Bins, tempbendchi2, 1.),
0271               temp_nstub,
0272               temp_ninterior,
0273               digitise(TTTrack_TrackWord::chi2RPhiBins, tempchi2rphi, (double)setup_->kfoutchi2rphiConv()),
0274               digitise(TTTrack_TrackWord::chi2RZBins, tempchi2rz, (double)setup_->kfoutchi2rzConv())};
0275
0276           // Run BDT emulation and package output into 3 bits
0277           // output needs sigmoid transformation applied
0278           tempTQMVAPreSig = trackQualityModel_->runEmulatedTQ(trackQuality_inputs);
0279           TTBV tqMVA(digitise(L1TrackQuality::getTqMVAPreSigBins(), tempTQMVAPreSig, 1.0),
0280                      TTTrack_TrackWord::TrackBitWidths::kMVAQualitySize,
0281                      false);
0282
0283           // Build 32 bit partial tracks for outputting in 64 bit packets
0284           //                  12 +  3       +  7         +  3    +  6
0285           TTBV partialTrack3((d0 + bendChi2 + hitPattern + tqMVA + extraMVA), partialTrackWordBits_, false);
0286           //                  16   + 12    + 4
0287           TTBV partialTrack2((tanL + z0 + chi2rz), partialTrackWordBits_, false);
0288           //                    1        + 15   +  12 +    4
0289           TTBV partialTrack1((trackValid + invR + phi0 + chi2rphi), partialTrackWordBits_, false);
0290
0291           int sortKey = (inTrack.sectorEta() < (int)(setup_->numSectorsEta() / 2)) ? 0 : 1;
0292           int nonantId = iLink / setup_->kfNumWorker();
0293           // Set correct bit to valid for track valid
0294           TrackKFOut temp_track(partialTrack1.set((partialTrackWordBits_ - 1)),
0295                                 partialTrack2,
0296                                 partialTrack3,
0297                                 sortKey,
0298                                 nonantId,
0299                                 track,
0300                                 iTrack,
0301                                 iLink,
0302                                 true);
0303           KFoutTracks.push_back(std::make_shared<TrackKFOut>(temp_track));
0304
0305           // add MVA to tttrack and add tttrack to collection
0306           temp_tttrack.settrkMVA1(1. / (1. + exp(tempTQMVAPreSig)));
0307           temp_tttrack.setTrackWordBits();
0308           ttTracks.emplace_back(temp_tttrack);
0309         }  // Iterate over Tracks
0310       }  // Iterate over Links
0311       const OrphanHandle<tt::TTTracks> orphanHandleTTTracks = iEvent.emplace(edPutTokenTTTracks_, std::move(ttTracks));
0312
0313       // sort partial KFout tracks into 18 separate links (nonant idx * eta idx) with tttrack ref info
0314       // 0th index order: [nonant 0 + negative eta, nonant 0 + positive eta, nonant 1 + negative eta, ...]
0315       struct kfoTrack_info {
0316         TTBV partialBits;
0317         TTTrackRef trackRef;
0318       };
0319       vector<vector<kfoTrack_info>> sortedPartialTracks(setup_->numRegions() * setup_->tfpNumChannel(),
0320                                                         vector<kfoTrack_info>(0));
0321       for (int i = 0; i < (int)KFoutTracks.size(); i++) {
0322         auto& kfoTrack = KFoutTracks.at(i);
0323         if (kfoTrack->dataValid()) {
0324           sortedPartialTracks[kfoTrack->nonantId() * setup_->tfpNumChannel() + kfoTrack->sortKey()].push_back(
0325               {kfoTrack->PartialTrack1(), TTTrackRef(orphanHandleTTTracks, i)});
0326           sortedPartialTracks[kfoTrack->nonantId() * setup_->tfpNumChannel() + kfoTrack->sortKey()].push_back(
0327               {kfoTrack->PartialTrack2(), TTTrackRef(orphanHandleTTTracks, i)});
0328           sortedPartialTracks[kfoTrack->nonantId() * setup_->tfpNumChannel() + kfoTrack->sortKey()].push_back(
0329               {kfoTrack->PartialTrack3(), TTTrackRef(orphanHandleTTTracks, i)});
0330         }
0331       }
0332       // fill remaining tracks allowed on each link (setup_->numFramesIO()) with null info
0333       kfoTrack_info nullTrack_info;
0334       for (int i = 0; i < (int)sortedPartialTracks.size(); i++) {
0335         // will not fill if any additional tracks if already above limit
0336         while ((int)sortedPartialTracks.at(i).size() < setup_->numFramesIO() * 2)
0337           sortedPartialTracks.at(i).push_back(nullTrack_info);
0338       }
0339
0340       // combine sorted partial tracks into proper format:
0341       // < TTTrackRef A, first 64 A bits >
0342       // < TTTrackRef B, last 32 A bits + first 32 B bits >
0343       // < TTTrackRef null, last 64 B bits >
0344       // ... repeat for next tracks
0345       const TTBV nullPartialBits(0, partialTrackWordBits_, false);
0346       const TTTrackRef nullTrackRef;
0347       int partialFactor = TTBV::S_ / partialTrackWordBits_;  //how many partial track words to combine in an output
0348       for (int iLink = 0; iLink < (int)sortedPartialTracks.size(); iLink++) {
0349         for (int iTrack = 0; iTrack < (int)sortedPartialTracks[iLink].size(); iTrack += partialFactor) {
0350           // if a partial track has no pair, pair it with null partial track
0351           if (iTrack + 1 == (int)sortedPartialTracks[iLink].size())
0352             sortedPartialTracks[iLink].push_back({nullPartialBits, nullTrackRef});
0353           // keep TTTrackRef null every third (96 bits / 32 partial bits) output packet
0354           TTTrackRef fillTrackRef;
0355           if ((iTrack / partialFactor + 1) % (TTTrack_TrackWord::kTrackWordSize / partialTrackWordBits_) != 0)
0356             fillTrackRef = sortedPartialTracks[iLink][iTrack + 1].trackRef;
0357
0358           // if there are too many output packets, truncate and put remaining outputs in lost collection
0359           if (iTrack / partialFactor < setup_->numFramesIO())
0360             accepted[iLink].emplace_back(
0361                 std::make_pair(fillTrackRef,
0362                                (sortedPartialTracks[iLink][iTrack].partialBits.slice(partialTrackWordBits_) +
0363                                 sortedPartialTracks[iLink][iTrack + 1].partialBits.slice(partialTrackWordBits_))
0364                                    .bs()));
0365           else
0366             lost[iLink].emplace_back(
0367                 std::make_pair(fillTrackRef,
0368                                (sortedPartialTracks[iLink][iTrack].partialBits.slice(partialTrackWordBits_) +
0369                                 sortedPartialTracks[iLink][iTrack + 1].partialBits.slice(partialTrackWordBits_))
0370                                    .bs()));
0371         }
0372       }
0373     }  // Config Supported
0374
0375     // store products
0376     iEvent.emplace(edPutTokenAccepted_, std::move(accepted));
0377     iEvent.emplace(edPutTokenLost_, std::move(lost));
0378   }
0379 }  // namespace trklet
0380
0381 DEFINE_FWK_MODULE(trklet::ProducerKFout);