Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-26 00:19:29

0001 ///////////////////////////////////////////////////////////////////////////
0002 //                                                                       //
0003 // Producer of TkJet,                                                    //
0004 // Cluster L1 tracks using fastjet                                       //
0005 //                                                                       //
0006 ///////////////////////////////////////////////////////////////////////////
0007 
0008 // system include files
0009 #include <memory>
0010 
0011 // user include files
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "DataFormats/Common/interface/Handle.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0023 #include "DataFormats/Math/interface/LorentzVector.h"
0024 
0025 // L1 objects
0026 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0027 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0028 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0029 #include "DataFormats/L1Trigger/interface/Vertex.h"
0030 
0031 // geometry
0032 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0033 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0034 
0035 #include <fastjet/JetDefinition.hh>
0036 
0037 #include <string>
0038 #include "TMath.h"
0039 #include "TH1.h"
0040 
0041 using namespace l1t;
0042 using namespace edm;
0043 using namespace std;
0044 
0045 //////////////////////////////
0046 //                          //
0047 //     CLASS DEFINITION     //
0048 //                          //
0049 //////////////////////////////
0050 
0051 class L1TrackFastJetProducer : public edm::stream::EDProducer<> {
0052 public:
0053   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0054   typedef std::vector<L1TTTrackType> L1TTTrackCollectionType;
0055 
0056   explicit L1TrackFastJetProducer(const edm::ParameterSet&);
0057   ~L1TrackFastJetProducer() override;
0058 
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060 
0061 private:
0062   //virtual void beginJob();
0063   void produce(edm::Event&, const edm::EventSetup&) override;
0064   //virtual void endJob();
0065 
0066   // track selection criteria
0067   const float trkZMax_;          // in [cm]
0068   const float trkChi2dofMax_;    // maximum track chi2dof
0069   const double trkBendChi2Max_;  // maximum track bendchi2
0070   const float trkPtMin_;         // in [GeV]
0071   const float trkEtaMax_;        // in [rad]
0072   const int trkNStubMin_;        // minimum number of stubs
0073   const int trkNPSStubMin_;      // minimum number of PS stubs
0074   const double deltaZ0Cut_;      // save with |L1z-z0| < maxZ0
0075   const double coneSize_;        // Use anti-kt with this cone size
0076   const bool doTightChi2_;
0077   const float trkPtTightChi2_;
0078   const float trkChi2dofTightChi2_;
0079   const bool displaced_;  //use prompt/displaced tracks
0080 
0081   const edm::EDGetTokenT<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > trackToken_;
0082   edm::EDGetTokenT<VertexCollection> pvToken_;
0083   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0084 };
0085 
0086 // constructor
0087 L1TrackFastJetProducer::L1TrackFastJetProducer(const edm::ParameterSet& iConfig)
0088     : trkZMax_((float)iConfig.getParameter<double>("trk_zMax")),
0089       trkChi2dofMax_((float)iConfig.getParameter<double>("trk_chi2dofMax")),
0090       trkBendChi2Max_(iConfig.getParameter<double>("trk_bendChi2Max")),
0091       trkPtMin_((float)iConfig.getParameter<double>("trk_ptMin")),
0092       trkEtaMax_((float)iConfig.getParameter<double>("trk_etaMax")),
0093       trkNStubMin_((int)iConfig.getParameter<int>("trk_nStubMin")),
0094       trkNPSStubMin_((int)iConfig.getParameter<int>("trk_nPSStubMin")),
0095       deltaZ0Cut_((float)iConfig.getParameter<double>("deltaZ0Cut")),
0096       coneSize_((float)iConfig.getParameter<double>("coneSize")),
0097       doTightChi2_(iConfig.getParameter<bool>("doTightChi2")),
0098       trkPtTightChi2_((float)iConfig.getParameter<double>("trk_ptTightChi2")),
0099       trkChi2dofTightChi2_((float)iConfig.getParameter<double>("trk_chi2dofTightChi2")),
0100       displaced_(iConfig.getParameter<bool>("displaced")),
0101       trackToken_(consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > >(
0102           iConfig.getParameter<edm::InputTag>("L1TrackInputTag"))),
0103       pvToken_(consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("L1PrimaryVertexTag"))),
0104       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))) {
0105   if (displaced_)
0106     produces<TkJetCollection>("L1TrackFastJetsExtended");
0107   else
0108     produces<TkJetCollection>("L1TrackFastJets");
0109 }
0110 
0111 // destructor
0112 L1TrackFastJetProducer::~L1TrackFastJetProducer() {}
0113 
0114 // producer
0115 void L1TrackFastJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0116   std::unique_ptr<TkJetCollection> L1TrackFastJets(new TkJetCollection);
0117 
0118   // L1 tracks
0119   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > TTTrackHandle;
0120   iEvent.getByToken(trackToken_, TTTrackHandle);
0121   std::vector<TTTrack<Ref_Phase2TrackerDigi_> >::const_iterator iterL1Track;
0122 
0123   // Tracker Topology
0124   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
0125 
0126   edm::Handle<l1t::VertexCollection> L1PrimaryVertexHandle;
0127   iEvent.getByToken(pvToken_, L1PrimaryVertexHandle);
0128 
0129   fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, coneSize_);
0130   std::vector<fastjet::PseudoJet> JetInputs;
0131 
0132   float recoVtx = L1PrimaryVertexHandle->begin()->z0();
0133   unsigned int this_l1track = 0;
0134   for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
0135     this_l1track++;
0136     float trk_pt = iterL1Track->momentum().perp();
0137     float trk_z0 = iterL1Track->z0();
0138     float trk_chi2dof = iterL1Track->chi2Red();
0139     float trk_bendchi2 = iterL1Track->stubPtConsistency();
0140     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
0141         theStubs = iterL1Track->getStubRefs();
0142     int trk_nstub = (int)theStubs.size();
0143 
0144     if (std::abs(trk_z0) > trkZMax_)
0145       continue;
0146     if (std::abs(iterL1Track->momentum().eta()) > trkEtaMax_)
0147       continue;
0148     if (trk_pt < trkPtMin_)
0149       continue;
0150     if (trk_nstub < trkNStubMin_)
0151       continue;
0152     if (trk_chi2dof > trkChi2dofMax_)
0153       continue;
0154     if (trk_bendchi2 > trkBendChi2Max_)
0155       continue;
0156     if (doTightChi2_ && (trk_pt > trkPtTightChi2_ && trk_chi2dof > trkChi2dofTightChi2_))
0157       continue;
0158 
0159     int trk_nPS = 0;
0160     for (int istub = 0; istub < trk_nstub; istub++) {
0161       DetId detId(theStubs.at(istub)->getDetId());
0162       bool tmp_isPS = false;
0163       if (detId.det() == DetId::Detector::Tracker) {
0164         if (detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3)
0165           tmp_isPS = true;
0166         else if (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)
0167           tmp_isPS = true;
0168       }
0169       if (tmp_isPS)
0170         trk_nPS++;
0171     }
0172     if (trk_nPS < trkNPSStubMin_)
0173       continue;
0174     if (std::abs(recoVtx - trk_z0) > deltaZ0Cut_)
0175       continue;
0176 
0177     fastjet::PseudoJet psuedoJet(iterL1Track->momentum().x(),
0178                                  iterL1Track->momentum().y(),
0179                                  iterL1Track->momentum().z(),
0180                                  iterL1Track->momentum().mag());
0181     JetInputs.push_back(psuedoJet);                     // input tracks for clustering
0182     JetInputs.back().set_user_index(this_l1track - 1);  // save track index in the collection
0183   }                                                     // end loop over tracks
0184 
0185   fastjet::ClusterSequence cs(JetInputs, jet_def);  // define the output jet collection
0186   std::vector<fastjet::PseudoJet> JetOutputs =
0187       fastjet::sorted_by_pt(cs.inclusive_jets(0));  // output jet collection, pT-ordered
0188 
0189   for (unsigned int ijet = 0; ijet < JetOutputs.size(); ++ijet) {
0190     math::XYZTLorentzVector jetP4(
0191         JetOutputs[ijet].px(), JetOutputs[ijet].py(), JetOutputs[ijet].pz(), JetOutputs[ijet].modp());
0192     float sumpt = 0;
0193     float avgZ = 0;
0194     std::vector<edm::Ptr<L1TTTrackType> > L1TrackPtrs;
0195     std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(cs.constituents(JetOutputs[ijet]));
0196 
0197     for (unsigned int i = 0; i < fjConstituents.size(); ++i) {
0198       auto index = fjConstituents[i].user_index();
0199       edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, index);
0200       L1TrackPtrs.push_back(trkPtr);  // L1Tracks in the jet
0201       sumpt = sumpt + trkPtr->momentum().perp();
0202       avgZ = avgZ + trkPtr->momentum().perp() * trkPtr->z0();
0203     }
0204     avgZ = avgZ / sumpt;
0205     edm::Ref<JetBxCollection> jetRef;
0206     TkJet trkJet(jetP4, jetRef, L1TrackPtrs, avgZ);
0207     L1TrackFastJets->push_back(trkJet);
0208   }  //end loop over Jet Outputs
0209 
0210   if (displaced_)
0211     iEvent.put(std::move(L1TrackFastJets), "L1TrackFastJetsExtended");
0212   else
0213     iEvent.put(std::move(L1TrackFastJets), "L1TrackFastJets");
0214 }
0215 
0216 //void L1TrackFastJetProducer::beginJob() {}
0217 
0218 //void L1TrackFastJetProducer::endJob() {}
0219 
0220 void L1TrackFastJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0221   //The following says we do not know what parameters are allowed so do no validation
0222   // Please change this to state exactly what you do use, even if it is no parameters
0223   edm::ParameterSetDescription desc;
0224   desc.setUnknown();
0225   descriptions.addDefault(desc);
0226 }
0227 
0228 //define this as a plug-in
0229 DEFINE_FWK_MODULE(L1TrackFastJetProducer);