L1TTrackMatch/plugins/TPFastJetProducer.cc

0001 ///////////////////////////////////////////////////////////////////////////
0002 //                                                                       //
0003 // Producer of TPFastJets,                                               //
0004 // Cluster tracking particles using fastjet                              //
0005 //                                                                       //
0006 // Created by: Claire Savard (Oct. 2023)                                 //
0007 //                                                                       //
0008 ///////////////////////////////////////////////////////////////////////////
0009
0010 // system include files
0011 #include <memory>
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0025 #include "DataFormats/Math/interface/LorentzVector.h"
0026
0027 // L1 objects
0028 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0029 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0030 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0031 #include "DataFormats/L1Trigger/interface/Vertex.h"
0032 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0033
0034 // truth object
0035 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0036 #include "SimDataFormats/Associations/interface/TTStubAssociationMap.h"
0037
0038 // geometry
0039 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0040 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0041
0042 #include <fastjet/JetDefinition.hh>
0043
0044 #include <string>
0045 #include "TMath.h"
0046 #include "TH1.h"
0047
0048 using namespace l1t;
0049 using namespace edm;
0050 using namespace std;
0051
0052 //////////////////////////////
0053 //                          //
0054 //     CLASS DEFINITION     //
0055 //                          //
0056 //////////////////////////////
0057
0058 class TPFastJetProducer : public edm::stream::EDProducer<> {
0059 public:
0060   explicit TPFastJetProducer(const edm::ParameterSet&);
0061   ~TPFastJetProducer() override = default;
0062
0063   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0064
0065 private:
0066   void produce(edm::Event&, const edm::EventSetup&) override;
0067
0068   // track selection criteria
0069   const float tpPtMin_;
0070   const float tpEtaMax_;
0071   const float tpZMax_;
0072   const int tpNStubMin_;
0073   const int tpNStubLayerMin_;
0074   const float coneSize_;  // Use anti-kt with this cone size
0075
0076   edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticleToken_;
0077   edm::EDGetTokenT<TTStubAssociationMap<Ref_Phase2TrackerDigi_>> ttStubMCTruthToken_;
0078   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0079 };
0080
0081 // constructor
0082 TPFastJetProducer::TPFastJetProducer(const edm::ParameterSet& iConfig)
0083     : tpPtMin_((float)iConfig.getParameter<double>("tp_ptMin")),
0084       tpEtaMax_((float)iConfig.getParameter<double>("tp_etaMax")),
0085       tpZMax_((float)iConfig.getParameter<double>("tp_zMax")),
0086       tpNStubMin_((int)iConfig.getParameter<int>("tp_nStubMin")),
0087       tpNStubLayerMin_((int)iConfig.getParameter<int>("tp_nStubLayerMin")),
0088       coneSize_((float)iConfig.getParameter<double>("coneSize")),
0089       trackingParticleToken_(
0090           consumes<std::vector<TrackingParticle>>(iConfig.getParameter<edm::InputTag>("TrackingParticleInputTag"))),
0091       ttStubMCTruthToken_(consumes<TTStubAssociationMap<Ref_Phase2TrackerDigi_>>(
0092           iConfig.getParameter<edm::InputTag>("MCTruthStubInputTag"))),
0093       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))) {
0094   produces<TkJetCollection>("TPFastJets");
0095 }
0096
0097 // producer
0098 void TPFastJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0099   std::unique_ptr<TkJetCollection> TPFastJets(new TkJetCollection);
0100
0101   // Tracking particles
0102   edm::Handle<std::vector<TrackingParticle>> TrackingParticleHandle;
0103   iEvent.getByToken(trackingParticleToken_, TrackingParticleHandle);
0104   std::vector<TrackingParticle>::const_iterator iterTP;
0105
0106   // MC truth association maps
0107   edm::Handle<TTStubAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTStubHandle;
0108   iEvent.getByToken(ttStubMCTruthToken_, MCTruthTTStubHandle);
0109
0110   // Tracker Topology
0111   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
0112
0113   fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, coneSize_);
0114   std::vector<fastjet::PseudoJet> JetInputs;
0115
0116   // loop over tps
0117   unsigned int this_tp = 0;
0118   for (iterTP = TrackingParticleHandle->begin(); iterTP != TrackingParticleHandle->end(); iterTP++) {
0119     edm::Ptr<TrackingParticle> tp_ptr(TrackingParticleHandle, this_tp);
0120     this_tp++;
0121
0122     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
0123         theStubRefs = MCTruthTTStubHandle->findTTStubRefs(tp_ptr);
0124     int nStubTP = (int)theStubRefs.size();
0125
0126     // how many layers/disks have stubs?
0127     int hasStubInLayer[11] = {0};
0128     for (auto& theStubRef : theStubRefs) {
0129       DetId detid(theStubRef->getDetId());
0130
0131       int layer = -1;
0132       if (detid.subdetId() == StripSubdetector::TOB) {
0133         layer = static_cast<int>(tTopo.layer(detid)) - 1;  //fill in array as entries 0-5
0134       } else if (detid.subdetId() == StripSubdetector::TID) {
0135         layer = static_cast<int>(tTopo.layer(detid)) + 5;  //fill in array as entries 6-10
0136       }
0137
0138       //treat genuine stubs separately (==2 is genuine, ==1 is not)
0139       if (MCTruthTTStubHandle->findTrackingParticlePtr(theStubRef).isNull() && hasStubInLayer[layer] < 2)
0140         hasStubInLayer[layer] = 1;
0141       else
0142         hasStubInLayer[layer] = 2;
0143     }
0144
0145     int nStubLayerTP = 0;
0146     for (int isum : hasStubInLayer) {
0147       if (isum >= 1)
0148         nStubLayerTP += 1;
0149     }
0150
0151     // tp quality cuts to match L1 tracks
0152     if (iterTP->pt() < tpPtMin_)
0153       continue;
0154     if (fabs(iterTP->eta()) > tpEtaMax_)
0155       continue;
0156     if (nStubTP < tpNStubMin_)
0157       continue;
0158     if (nStubLayerTP < tpNStubLayerMin_)
0159       continue;
0160     if (fabs(iterTP->z0()) > tpZMax_)
0161       continue;
0162     if (iterTP->charge() == 0.)  // extra check that all tps are charged
0163       continue;
0164     if (iterTP->eventId().event() > 0)  // only select hard interaction tps
0165       continue;
0166
0167     fastjet::PseudoJet psuedoJet(iterTP->px(), iterTP->py(), iterTP->pz(), iterTP->energy());
0168     JetInputs.push_back(psuedoJet);                // input tps for clustering
0169     JetInputs.back().set_user_index(this_tp - 1);  // save tp index in the collection
0170   }  // end loop over tps
0171
0172   fastjet::ClusterSequence cs(JetInputs, jet_def);  // define the output jet collection
0173   std::vector<fastjet::PseudoJet> JetOutputs =
0174       fastjet::sorted_by_pt(cs.inclusive_jets(0));  // output jet collection, pT-ordered
0175
0176   for (unsigned int ijet = 0; ijet < JetOutputs.size(); ++ijet) {
0177     math::XYZTLorentzVector jetP4(
0178         JetOutputs[ijet].px(), JetOutputs[ijet].py(), JetOutputs[ijet].pz(), JetOutputs[ijet].modp());
0179     float sumpt = 0;
0180     float avgZ = 0;
0181     std::vector<edm::Ptr<TrackingParticle>> tpPtrs;
0182     std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(cs.constituents(JetOutputs[ijet]));
0183
0184     for (unsigned int i = 0; i < fjConstituents.size(); ++i) {
0185       auto index = fjConstituents[i].user_index();
0186       edm::Ptr<TrackingParticle> tpPtr(TrackingParticleHandle, index);
0187       tpPtrs.push_back(tpPtr);  // tracking particles in the jet
0188       sumpt = sumpt + tpPtr->pt();
0189       avgZ = avgZ + tpPtr->pt() * tpPtr->z0();
0190     }
0191     avgZ = avgZ / sumpt;
0192     edm::Ref<JetBxCollection> jetRef;
0193     std::vector<edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>>> dummyL1TrackPtrs;  // can't create TkJet with tp references
0194     TkJet tpJet(jetP4, dummyL1TrackPtrs, avgZ, fjConstituents.size(), 0, 0, 0, false);
0195     TPFastJets->push_back(tpJet);
0196   }  //end loop over Jet Outputs
0197
0198   iEvent.put(std::move(TPFastJets), "TPFastJets");
0199 }
0200
0201 void TPFastJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0202   // The following says we do not know what parameters are allowed so do no validation
0203   edm::ParameterSetDescription desc;
0204   desc.add<edm::InputTag>("TrackingParticleInputTag", edm::InputTag("mix", "MergedTrackTruth"));
0205   desc.add<edm::InputTag>("MCTruthStubInputTag", edm::InputTag("TTStubAssociatorFromPixelDigis", "StubAccepted"));
0206   desc.add<double>("tp_ptMin", 2.0);
0207   desc.add<double>("tp_etaMax", 2.4);
0208   desc.add<double>("tp_zMax", 15.);
0209   desc.add<int>("tp_nStubMin", 4);
0210   desc.add<int>("tp_nStubLayerMin", 4);
0211   desc.add<double>("coneSize", 0.4);
0212   descriptions.add("tpFastJets", desc);
0213 }
0214
0215 //define this as a plug-in
0216 DEFINE_FWK_MODULE(TPFastJetProducer);