TrackFindingTracklet/test/L1TrackNtupleMaker.cc

0001 //////////////////////////////////////////////////////////////////////
0002 //                                                                  //
0003 //  Analyzer for making mini-ntuple for L1 track performance plots  //
0004 //                                                                  //
0005 //////////////////////////////////////////////////////////////////////
0006
0007 ////////////////////
0008 // FRAMEWORK HEADERS
0009 #include "FWCore/PluginManager/interface/ModuleDef.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Run.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/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021 ///////////////////////
0022 // DATA FORMATS HEADERS
0023 #include "DataFormats/Common/interface/Handle.h"
0024 #include "DataFormats/Common/interface/Ref.h"
0025
0026 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0027 #include "DataFormats/L1TrackTrigger/interface/TTCluster.h"
0028 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0029 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0030 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0031 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0032 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0033 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0034 #include "SimDataFormats/Associations/interface/TTClusterAssociationMap.h"
0035 #include "SimDataFormats/Associations/interface/TTStubAssociationMap.h"
0036 #include "SimDataFormats/Associations/interface/TTTrackAssociationMap.h"
0037 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0038
0039 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0040 #include "DataFormats/JetReco/interface/GenJet.h"
0041
0042 ////////////////////////////
0043 // DETECTOR GEOMETRY HEADERS
0044 #include "MagneticField/Engine/interface/MagneticField.h"
0045 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0046 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0047 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0049 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0050
0051 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0052 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0053 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0054 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0055
0056 ////////////////
0057 // PHYSICS TOOLS
0058 #include "L1Trigger/TrackFindingTracklet/interface/HitPatternHelper.h"
0059 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0060 #include "CLHEP/Units/PhysicalConstants.h"
0061
0062 ///////////////
0063 // ROOT HEADERS
0064 #include <TROOT.h>
0065 #include <TCanvas.h>
0066 #include <TTree.h>
0067 #include <TFile.h>
0068 #include <TF1.h>
0069 #include <TH2F.h>
0070 #include <TH1F.h>
0071
0072 //////////////
0073 // STD HEADERS
0074 #include <memory>
0075 #include <string>
0076 #include <iostream>
0077
0078 //////////////
0079 // NAMESPACES
0080 using namespace std;
0081 using namespace edm;
0082
0083 //////////////////////////////
0084 //                          //
0085 //     CLASS DEFINITION     //
0086 //                          //
0087 //////////////////////////////
0088
0089 class L1TrackNtupleMaker : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0090 public:
0091   // Constructor/destructor
0092   explicit L1TrackNtupleMaker(const edm::ParameterSet& iConfig);
0093   ~L1TrackNtupleMaker() override;
0094
0095   // Mandatory methods
0096   void beginJob() override;
0097   void endJob() override;
0098   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0099   void beginRun(const Run& iEvent, const EventSetup& iSetup) override {}
0100   void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0101
0102 protected:
0103 private:
0104   //-----------------------------------------------------------------------------------------------
0105   // Containers of parameters passed by python configuration file
0106   edm::ParameterSet config;
0107
0108   int MyProcess;       // 11/13/211 for single electrons/muons/pions, 6/15 for pions from ttbar/taus, 1 for inclusive
0109   bool DebugMode;      // lots of debug printout statements
0110   bool SaveAllTracks;  // store in ntuples not only truth-matched tracks but ALL tracks
0111   bool SaveStubs;      // option to save also stubs in the ntuples (makes them large...)
0112   int L1Tk_nPar;       // use 4 or 5 parameter track fit?
0113   int TP_minNStub;  // require TPs to have >= minNStub (defining efficiency denominator) (==0 means to only require >= 1 cluster)
0114   int TP_minNStubLayer;  // require TPs to have stubs in >= minNStubLayer layers/disks (defining efficiency denominator)
0115   double TP_minPt;       // save TPs with pt > minPt
0116   double TP_maxEta;      // save TPs with |eta| < maxEta
0117   double TP_maxZ0;       // save TPs with |z0| < maxZ0
0118   int L1Tk_minNStub;     // require L1 tracks to have >= minNStub (this is mostly for tracklet purposes)
0119
0120   bool TrackingInJets;  // do tracking in jets?
0121
0122   edm::InputTag L1TrackInputTag;       // L1 track collection
0123   edm::InputTag MCTruthTrackInputTag;  // MC truth collection
0124   edm::InputTag MCTruthClusterInputTag;
0125   edm::InputTag L1StubInputTag;
0126   edm::InputTag MCTruthStubInputTag;
0127   edm::InputTag TrackingParticleInputTag;
0128   edm::InputTag TrackingVertexInputTag;
0129   edm::InputTag GenJetInputTag;
0130
0131   edm::EDGetTokenT<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_> > > ttClusterToken_;
0132   edm::EDGetTokenT<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> > > ttStubToken_;
0133   edm::EDGetTokenT<TTClusterAssociationMap<Ref_Phase2TrackerDigi_> > ttClusterMCTruthToken_;
0134   edm::EDGetTokenT<TTStubAssociationMap<Ref_Phase2TrackerDigi_> > ttStubMCTruthToken_;
0135
0136   edm::EDGetTokenT<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > ttTrackToken_;
0137   edm::EDGetTokenT<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> > ttTrackMCTruthToken_;
0138
0139   edm::EDGetTokenT<std::vector<TrackingParticle> > TrackingParticleToken_;
0140   edm::EDGetTokenT<std::vector<TrackingVertex> > TrackingVertexToken_;
0141
0142   edm::EDGetTokenT<std::vector<reco::GenJet> > GenJetToken_;
0143
0144   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> getTokenTrackerGeom_;
0145   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> getTokenTrackerTopo_;
0146   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> getTokenBField_;
0147   edm::ESGetToken<hph::Setup, hph::SetupRcd> getTokenHPHSetup_;
0148   //-----------------------------------------------------------------------------------------------
0149   // tree & branches for mini-ntuple
0150
0151   bool available_;  // ROOT file for histograms is open.
0152
0153   TTree* eventTree;
0154
0155   // all L1 tracks
0156   std::vector<float>* m_trk_pt;
0157   std::vector<float>* m_trk_eta;
0158   std::vector<float>* m_trk_phi;
0159   std::vector<float>* m_trk_d0;  // (filled if L1Tk_nPar==5, else 999)
0160   std::vector<float>* m_trk_z0;
0161   std::vector<float>* m_trk_chi2;
0162   std::vector<float>* m_trk_chi2_dof;
0163   std::vector<float>* m_trk_chi2rphi;
0164   std::vector<float>* m_trk_chi2rphi_dof;
0165   std::vector<float>* m_trk_chi2rz;
0166   std::vector<float>* m_trk_chi2rz_dof;
0167   std::vector<float>* m_trk_bendchi2;
0168   std::vector<int>* m_trk_nstub;
0169   std::vector<int>* m_trk_lhits;
0170   std::vector<int>* m_trk_dhits;
0171   std::vector<int>* m_trk_seed;
0172   std::vector<int>* m_trk_hitpattern;
0173   std::vector<int>* m_trk_lhits_hitpattern;  // 6-digit hit mask (barrel layer only) dervied from hitpattern
0174   std::vector<int>* m_trk_dhits_hitpattern;  // disk only
0175   std::vector<int>* m_trk_nPSstub_hitpattern;
0176   std::vector<int>* m_trk_n2Sstub_hitpattern;
0177   std::vector<int>* m_trk_nLostPSstub_hitpattern;
0178   std::vector<int>* m_trk_nLost2Sstub_hitpattern;
0179   std::vector<int>* m_trk_nLoststub_V1_hitpattern;  // Same as the definiton of "nlaymiss_interior" in TrackQuality.cc
0180   std::vector<int>* m_trk_nLoststub_V2_hitpattern;  // A tighter version of "nlaymiss_interior"
0181   std::vector<unsigned int>* m_trk_phiSector;
0182   std::vector<int>* m_trk_etaSector;
0183   std::vector<int>* m_trk_genuine;
0184   std::vector<int>* m_trk_loose;
0185   std::vector<int>* m_trk_unknown;
0186   std::vector<int>* m_trk_combinatoric;
0187   std::vector<int>* m_trk_fake;  //0 fake, 1 track from primary interaction, 2 secondary track
0188   std::vector<float>* m_trk_MVA1;
0189   std::vector<int>* m_trk_matchtp_pdgid;
0190   std::vector<float>* m_trk_matchtp_pt;
0191   std::vector<float>* m_trk_matchtp_eta;
0192   std::vector<float>* m_trk_matchtp_phi;
0193   std::vector<float>* m_trk_matchtp_z0;
0194   std::vector<float>* m_trk_matchtp_dxy;
0195   std::vector<float>* m_trk_matchtp_d0;
0196   std::vector<int>* m_trk_injet;          //is the track within dR<0.4 of a genjet with pt > 30 GeV?
0197   std::vector<int>* m_trk_injet_highpt;   //is the track within dR<0.4 of a genjet with pt > 100 GeV?
0198   std::vector<int>* m_trk_injet_vhighpt;  //is the track within dR<0.4 of a genjet with pt > 200 GeV?
0199
0200   // all tracking particles
0201   std::vector<float>* m_tp_pt;
0202   std::vector<float>* m_tp_eta;
0203   std::vector<float>* m_tp_phi;
0204   std::vector<float>* m_tp_dxy;
0205   std::vector<float>* m_tp_d0;
0206   std::vector<float>* m_tp_z0;
0207   std::vector<float>* m_tp_d0_prod;
0208   std::vector<float>* m_tp_z0_prod;
0209   std::vector<int>* m_tp_pdgid;
0210   std::vector<int>* m_tp_nmatch;
0211   std::vector<int>* m_tp_nstub;
0212   std::vector<int>* m_tp_eventid;
0213   std::vector<int>* m_tp_charge;
0214   std::vector<int>* m_tp_injet;
0215   std::vector<int>* m_tp_injet_highpt;
0216   std::vector<int>* m_tp_injet_vhighpt;
0217
0218   // *L1 track* properties if m_tp_nmatch > 0
0219   std::vector<float>* m_matchtrk_pt;
0220   std::vector<float>* m_matchtrk_eta;
0221   std::vector<float>* m_matchtrk_phi;
0222   std::vector<float>* m_matchtrk_d0;  //this variable is only filled if L1Tk_nPar==5
0223   std::vector<float>* m_matchtrk_z0;
0224   std::vector<float>* m_matchtrk_chi2;
0225   std::vector<float>* m_matchtrk_chi2_dof;
0226   std::vector<float>* m_matchtrk_chi2rphi;
0227   std::vector<float>* m_matchtrk_chi2rphi_dof;
0228   std::vector<float>* m_matchtrk_chi2rz;
0229   std::vector<float>* m_matchtrk_chi2rz_dof;
0230   std::vector<float>* m_matchtrk_bendchi2;
0231   std::vector<float>* m_matchtrk_MVA1;
0232   std::vector<int>* m_matchtrk_nstub;
0233   std::vector<int>* m_matchtrk_lhits;
0234   std::vector<int>* m_matchtrk_dhits;
0235   std::vector<int>* m_matchtrk_seed;
0236   std::vector<int>* m_matchtrk_hitpattern;
0237   std::vector<int>* m_matchtrk_injet;
0238   std::vector<int>* m_matchtrk_injet_highpt;
0239   std::vector<int>* m_matchtrk_injet_vhighpt;
0240
0241   // ALL stubs
0242   std::vector<float>* m_allstub_x;
0243   std::vector<float>* m_allstub_y;
0244   std::vector<float>* m_allstub_z;
0245
0246   std::vector<int>* m_allstub_isBarrel;  // stub is in barrel (1) or in disk (0)
0247   std::vector<int>* m_allstub_layer;
0248   std::vector<int>* m_allstub_isPSmodule;
0249   std::vector<int>* m_allstub_isTiltedBarrel;
0250
0251   std::vector<float>* m_allstub_trigDisplace;
0252   std::vector<float>* m_allstub_trigOffset;
0253   std::vector<float>* m_allstub_trigPos;
0254   std::vector<float>* m_allstub_trigBend;
0255
0256   // stub associated with tracking particle ?
0257   std::vector<int>* m_allstub_matchTP_pdgid;  // -999 if not matched
0258   std::vector<float>* m_allstub_matchTP_pt;   // -999 if not matched
0259   std::vector<float>* m_allstub_matchTP_eta;  // -999 if not matched
0260   std::vector<float>* m_allstub_matchTP_phi;  // -999 if not matched
0261
0262   std::vector<int>* m_allstub_genuine;
0263
0264   // track jet variables (for each gen jet, store the sum of pt of TPs / tracks inside jet cone)
0265   std::vector<float>* m_jet_eta;
0266   std::vector<float>* m_jet_phi;
0267   std::vector<float>* m_jet_pt;
0268   std::vector<float>* m_jet_tp_sumpt;
0269   std::vector<float>* m_jet_trk_sumpt;
0270   std::vector<float>* m_jet_matchtrk_sumpt;
0271 };
0272
0273 //////////////////////////////////
0274 //                              //
0275 //     CLASS IMPLEMENTATION     //
0276 //                              //
0277 //////////////////////////////////
0278
0279 //////////////
0280 // CONSTRUCTOR
0281 L1TrackNtupleMaker::L1TrackNtupleMaker(edm::ParameterSet const& iConfig) : config(iConfig) {
0282   usesResource("TFileService");
0283   MyProcess = iConfig.getParameter<int>("MyProcess");
0284   DebugMode = iConfig.getParameter<bool>("DebugMode");
0285   SaveAllTracks = iConfig.getParameter<bool>("SaveAllTracks");
0286   SaveStubs = iConfig.getParameter<bool>("SaveStubs");
0287   L1Tk_nPar = iConfig.getParameter<int>("L1Tk_nPar");
0288   TP_minNStub = iConfig.getParameter<int>("TP_minNStub");
0289   TP_minNStubLayer = iConfig.getParameter<int>("TP_minNStubLayer");
0290   TP_minPt = iConfig.getParameter<double>("TP_minPt");
0291   TP_maxEta = iConfig.getParameter<double>("TP_maxEta");
0292   TP_maxZ0 = iConfig.getParameter<double>("TP_maxZ0");
0293   L1TrackInputTag = iConfig.getParameter<edm::InputTag>("L1TrackInputTag");
0294   MCTruthTrackInputTag = iConfig.getParameter<edm::InputTag>("MCTruthTrackInputTag");
0295   L1Tk_minNStub = iConfig.getParameter<int>("L1Tk_minNStub");
0296
0297   TrackingInJets = iConfig.getParameter<bool>("TrackingInJets");
0298
0299   L1StubInputTag = iConfig.getParameter<edm::InputTag>("L1StubInputTag");
0300   MCTruthClusterInputTag = iConfig.getParameter<edm::InputTag>("MCTruthClusterInputTag");
0301   MCTruthStubInputTag = iConfig.getParameter<edm::InputTag>("MCTruthStubInputTag");
0302   TrackingParticleInputTag = iConfig.getParameter<edm::InputTag>("TrackingParticleInputTag");
0303   TrackingVertexInputTag = iConfig.getParameter<edm::InputTag>("TrackingVertexInputTag");
0304   GenJetInputTag = iConfig.getParameter<edm::InputTag>("GenJetInputTag");
0305
0306   ttTrackToken_ = consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > >(L1TrackInputTag);
0307   ttTrackMCTruthToken_ = consumes<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> >(MCTruthTrackInputTag);
0308   ttStubToken_ = consumes<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> > >(L1StubInputTag);
0309   ttClusterMCTruthToken_ = consumes<TTClusterAssociationMap<Ref_Phase2TrackerDigi_> >(MCTruthClusterInputTag);
0310   ttStubMCTruthToken_ = consumes<TTStubAssociationMap<Ref_Phase2TrackerDigi_> >(MCTruthStubInputTag);
0311
0312   TrackingParticleToken_ = consumes<std::vector<TrackingParticle> >(TrackingParticleInputTag);
0313   TrackingVertexToken_ = consumes<std::vector<TrackingVertex> >(TrackingVertexInputTag);
0314   GenJetToken_ = consumes<std::vector<reco::GenJet> >(GenJetInputTag);
0315
0316   getTokenTrackerGeom_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0317   getTokenTrackerTopo_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0318   getTokenBField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0319   getTokenHPHSetup_ = esConsumes<hph::Setup, hph::SetupRcd>();
0320 }
0321
0322 /////////////
0323 // DESTRUCTOR
0324 L1TrackNtupleMaker::~L1TrackNtupleMaker() {}
0325
0326 //////////
0327 // END JOB
0328 void L1TrackNtupleMaker::endJob() {
0329   // things to be done at the exit of the event Loop
0330   edm::LogVerbatim("Tracklet") << "L1TrackNtupleMaker::endJob";
0331
0332   // clean up
0333   delete m_trk_pt;
0334   delete m_trk_eta;
0335   delete m_trk_phi;
0336   delete m_trk_z0;
0337   delete m_trk_d0;
0338   delete m_trk_chi2;
0339   delete m_trk_chi2_dof;
0340   delete m_trk_chi2rphi;
0341   delete m_trk_chi2rphi_dof;
0342   delete m_trk_chi2rz;
0343   delete m_trk_chi2rz_dof;
0344   delete m_trk_bendchi2;
0345   delete m_trk_nstub;
0346   delete m_trk_lhits;
0347   delete m_trk_dhits;
0348   delete m_trk_seed;
0349   delete m_trk_hitpattern;
0350   delete m_trk_lhits_hitpattern;
0351   delete m_trk_dhits_hitpattern;
0352   delete m_trk_nPSstub_hitpattern;
0353   delete m_trk_n2Sstub_hitpattern;
0354   delete m_trk_nLostPSstub_hitpattern;
0355   delete m_trk_nLost2Sstub_hitpattern;
0356   delete m_trk_nLoststub_V1_hitpattern;
0357   delete m_trk_nLoststub_V2_hitpattern;
0358   delete m_trk_phiSector;
0359   delete m_trk_etaSector;
0360   delete m_trk_genuine;
0361   delete m_trk_loose;
0362   delete m_trk_unknown;
0363   delete m_trk_combinatoric;
0364   delete m_trk_fake;
0365   delete m_trk_MVA1;
0366   delete m_trk_matchtp_pdgid;
0367   delete m_trk_matchtp_pt;
0368   delete m_trk_matchtp_eta;
0369   delete m_trk_matchtp_phi;
0370   delete m_trk_matchtp_z0;
0371   delete m_trk_matchtp_dxy;
0372   delete m_trk_matchtp_d0;
0373   delete m_trk_injet;
0374   delete m_trk_injet_highpt;
0375   delete m_trk_injet_vhighpt;
0376
0377   delete m_tp_pt;
0378   delete m_tp_eta;
0379   delete m_tp_phi;
0380   delete m_tp_dxy;
0381   delete m_tp_d0;
0382   delete m_tp_z0;
0383   delete m_tp_d0_prod;
0384   delete m_tp_z0_prod;
0385   delete m_tp_pdgid;
0386   delete m_tp_nmatch;
0387   delete m_tp_nstub;
0388   delete m_tp_eventid;
0389   delete m_tp_charge;
0390   delete m_tp_injet;
0391   delete m_tp_injet_highpt;
0392   delete m_tp_injet_vhighpt;
0393
0394   delete m_matchtrk_pt;
0395   delete m_matchtrk_eta;
0396   delete m_matchtrk_phi;
0397   delete m_matchtrk_z0;
0398   delete m_matchtrk_d0;
0399   delete m_matchtrk_chi2;
0400   delete m_matchtrk_chi2_dof;
0401   delete m_matchtrk_chi2rphi;
0402   delete m_matchtrk_chi2rphi_dof;
0403   delete m_matchtrk_chi2rz;
0404   delete m_matchtrk_chi2rz_dof;
0405   delete m_matchtrk_bendchi2;
0406   delete m_matchtrk_MVA1;
0407   delete m_matchtrk_nstub;
0408   delete m_matchtrk_dhits;
0409   delete m_matchtrk_lhits;
0410   delete m_matchtrk_seed;
0411   delete m_matchtrk_hitpattern;
0412   delete m_matchtrk_injet;
0413   delete m_matchtrk_injet_highpt;
0414   delete m_matchtrk_injet_vhighpt;
0415
0416   delete m_allstub_x;
0417   delete m_allstub_y;
0418   delete m_allstub_z;
0419   delete m_allstub_isBarrel;
0420   delete m_allstub_layer;
0421   delete m_allstub_isPSmodule;
0422   delete m_allstub_trigDisplace;
0423   delete m_allstub_trigOffset;
0424   delete m_allstub_trigPos;
0425   delete m_allstub_trigBend;
0426   delete m_allstub_matchTP_pdgid;
0427   delete m_allstub_matchTP_pt;
0428   delete m_allstub_matchTP_eta;
0429   delete m_allstub_matchTP_phi;
0430   delete m_allstub_genuine;
0431
0432   delete m_jet_eta;
0433   delete m_jet_phi;
0434   delete m_jet_pt;
0435   delete m_jet_tp_sumpt;
0436   delete m_jet_trk_sumpt;
0437   delete m_jet_matchtrk_sumpt;
0438 }
0439
0440 ////////////
0441 // BEGIN JOB
0442 void L1TrackNtupleMaker::beginJob() {
0443   // things to be done before entering the event Loop
0444   edm::LogVerbatim("Tracklet") << "L1TrackNtupleMaker::beginJob";
0445
0446   //-----------------------------------------------------------------------------------------------
0447   // book histograms / make ntuple
0448   edm::Service<TFileService> fs;
0449   available_ = fs.isAvailable();
0450   if (not available_)
0451     return;  // No ROOT file open.
0452
0453   // initilize
0454   m_trk_pt = new std::vector<float>;
0455   m_trk_eta = new std::vector<float>;
0456   m_trk_phi = new std::vector<float>;
0457   m_trk_z0 = new std::vector<float>;
0458   m_trk_d0 = new std::vector<float>;
0459   m_trk_chi2 = new std::vector<float>;
0460   m_trk_chi2_dof = new std::vector<float>;
0461   m_trk_chi2rphi = new std::vector<float>;
0462   m_trk_chi2rphi_dof = new std::vector<float>;
0463   m_trk_chi2rz = new std::vector<float>;
0464   m_trk_chi2rz_dof = new std::vector<float>;
0465   m_trk_bendchi2 = new std::vector<float>;
0466   m_trk_nstub = new std::vector<int>;
0467   m_trk_lhits = new std::vector<int>;
0468   m_trk_dhits = new std::vector<int>;
0469   m_trk_seed = new std::vector<int>;
0470   m_trk_hitpattern = new std::vector<int>;
0471   m_trk_lhits_hitpattern = new std::vector<int>;
0472   m_trk_dhits_hitpattern = new std::vector<int>;
0473   m_trk_nPSstub_hitpattern = new std::vector<int>;
0474   m_trk_n2Sstub_hitpattern = new std::vector<int>;
0475   m_trk_nLostPSstub_hitpattern = new std::vector<int>;
0476   m_trk_nLost2Sstub_hitpattern = new std::vector<int>;
0477   m_trk_nLoststub_V1_hitpattern = new std::vector<int>;
0478   m_trk_nLoststub_V2_hitpattern = new std::vector<int>;
0479   m_trk_phiSector = new std::vector<unsigned int>;
0480   m_trk_etaSector = new std::vector<int>;
0481   m_trk_genuine = new std::vector<int>;
0482   m_trk_loose = new std::vector<int>;
0483   m_trk_unknown = new std::vector<int>;
0484   m_trk_combinatoric = new std::vector<int>;
0485   m_trk_fake = new std::vector<int>;
0486   m_trk_MVA1 = new std::vector<float>;
0487   m_trk_matchtp_pdgid = new std::vector<int>;
0488   m_trk_matchtp_pt = new std::vector<float>;
0489   m_trk_matchtp_eta = new std::vector<float>;
0490   m_trk_matchtp_phi = new std::vector<float>;
0491   m_trk_matchtp_z0 = new std::vector<float>;
0492   m_trk_matchtp_dxy = new std::vector<float>;
0493   m_trk_matchtp_d0 = new std::vector<float>;
0494   m_trk_injet = new std::vector<int>;
0495   m_trk_injet_highpt = new std::vector<int>;
0496   m_trk_injet_vhighpt = new std::vector<int>;
0497
0498   m_tp_pt = new std::vector<float>;
0499   m_tp_eta = new std::vector<float>;
0500   m_tp_phi = new std::vector<float>;
0501   m_tp_dxy = new std::vector<float>;
0502   m_tp_d0 = new std::vector<float>;
0503   m_tp_z0 = new std::vector<float>;
0504   m_tp_d0_prod = new std::vector<float>;
0505   m_tp_z0_prod = new std::vector<float>;
0506   m_tp_pdgid = new std::vector<int>;
0507   m_tp_nmatch = new std::vector<int>;
0508   m_tp_nstub = new std::vector<int>;
0509   m_tp_eventid = new std::vector<int>;
0510   m_tp_charge = new std::vector<int>;
0511   m_tp_injet = new std::vector<int>;
0512   m_tp_injet_highpt = new std::vector<int>;
0513   m_tp_injet_vhighpt = new std::vector<int>;
0514
0515   m_matchtrk_pt = new std::vector<float>;
0516   m_matchtrk_eta = new std::vector<float>;
0517   m_matchtrk_phi = new std::vector<float>;
0518   m_matchtrk_z0 = new std::vector<float>;
0519   m_matchtrk_d0 = new std::vector<float>;
0520   m_matchtrk_chi2 = new std::vector<float>;
0521   m_matchtrk_chi2_dof = new std::vector<float>;
0522   m_matchtrk_chi2rphi = new std::vector<float>;
0523   m_matchtrk_chi2rphi_dof = new std::vector<float>;
0524   m_matchtrk_chi2rz = new std::vector<float>;
0525   m_matchtrk_chi2rz_dof = new std::vector<float>;
0526   m_matchtrk_bendchi2 = new std::vector<float>;
0527   m_matchtrk_MVA1 = new std::vector<float>;
0528   m_matchtrk_nstub = new std::vector<int>;
0529   m_matchtrk_dhits = new std::vector<int>;
0530   m_matchtrk_lhits = new std::vector<int>;
0531   m_matchtrk_seed = new std::vector<int>;
0532   m_matchtrk_hitpattern = new std::vector<int>;
0533   m_matchtrk_injet = new std::vector<int>;
0534   m_matchtrk_injet_highpt = new std::vector<int>;
0535   m_matchtrk_injet_vhighpt = new std::vector<int>;
0536
0537   m_allstub_x = new std::vector<float>;
0538   m_allstub_y = new std::vector<float>;
0539   m_allstub_z = new std::vector<float>;
0540
0541   m_allstub_isBarrel = new std::vector<int>;
0542   m_allstub_layer = new std::vector<int>;
0543   m_allstub_isPSmodule = new std::vector<int>;
0544   m_allstub_isTiltedBarrel = new std::vector<int>;
0545   m_allstub_trigDisplace = new std::vector<float>;
0546   m_allstub_trigOffset = new std::vector<float>;
0547   m_allstub_trigPos = new std::vector<float>;
0548   m_allstub_trigBend = new std::vector<float>;
0549
0550   m_allstub_matchTP_pdgid = new std::vector<int>;
0551   m_allstub_matchTP_pt = new std::vector<float>;
0552   m_allstub_matchTP_eta = new std::vector<float>;
0553   m_allstub_matchTP_phi = new std::vector<float>;
0554
0555   m_allstub_genuine = new std::vector<int>;
0556
0557   m_jet_eta = new std::vector<float>;
0558   m_jet_phi = new std::vector<float>;
0559   m_jet_pt = new std::vector<float>;
0560   m_jet_tp_sumpt = new std::vector<float>;
0561   m_jet_trk_sumpt = new std::vector<float>;
0562   m_jet_matchtrk_sumpt = new std::vector<float>;
0563
0564   // ntuple
0565   eventTree = fs->make<TTree>("eventTree", "Event tree");
0566
0567   if (SaveAllTracks) {
0568     eventTree->Branch("trk_pt", &m_trk_pt);
0569     eventTree->Branch("trk_eta", &m_trk_eta);
0570     eventTree->Branch("trk_phi", &m_trk_phi);
0571     eventTree->Branch("trk_d0", &m_trk_d0);
0572     eventTree->Branch("trk_z0", &m_trk_z0);
0573     eventTree->Branch("trk_chi2", &m_trk_chi2);
0574     eventTree->Branch("trk_chi2_dof", &m_trk_chi2_dof);
0575     eventTree->Branch("trk_chi2rphi", &m_trk_chi2rphi);
0576     eventTree->Branch("trk_chi2rphi_dof", &m_trk_chi2rphi_dof);
0577     eventTree->Branch("trk_chi2rz", &m_trk_chi2rz);
0578     eventTree->Branch("trk_chi2rz_dof", &m_trk_chi2rz_dof);
0579     eventTree->Branch("trk_bendchi2", &m_trk_bendchi2);
0580     eventTree->Branch("trk_nstub", &m_trk_nstub);
0581     eventTree->Branch("trk_lhits", &m_trk_lhits);
0582     eventTree->Branch("trk_dhits", &m_trk_dhits);
0583     eventTree->Branch("trk_seed", &m_trk_seed);
0584     eventTree->Branch("trk_hitpattern", &m_trk_hitpattern);
0585     eventTree->Branch("trk_lhits_hitpattern", &m_trk_lhits_hitpattern);
0586     eventTree->Branch("trk_dhits_hitpattern", &m_trk_dhits_hitpattern);
0587     eventTree->Branch("trk_nPSstub_hitpattern", &m_trk_nPSstub_hitpattern);
0588     eventTree->Branch("trk_n2Sstub_hitpattern", &m_trk_n2Sstub_hitpattern);
0589     eventTree->Branch("trk_nLostPSstub_hitpattern", &m_trk_nLostPSstub_hitpattern);
0590     eventTree->Branch("trk_nLost2Sstub_hitpattern", &m_trk_nLost2Sstub_hitpattern);
0591     eventTree->Branch("trk_nLoststub_V1_hitpattern", &m_trk_nLoststub_V1_hitpattern);
0592     eventTree->Branch("trk_nLoststub_V2_hitpattern", &m_trk_nLoststub_V2_hitpattern);
0593     eventTree->Branch("trk_phiSector", &m_trk_phiSector);
0594     eventTree->Branch("trk_etaSector", &m_trk_etaSector);
0595     eventTree->Branch("trk_genuine", &m_trk_genuine);
0596     eventTree->Branch("trk_loose", &m_trk_loose);
0597     eventTree->Branch("trk_unknown", &m_trk_unknown);
0598     eventTree->Branch("trk_combinatoric", &m_trk_combinatoric);
0599     eventTree->Branch("trk_fake", &m_trk_fake);
0600     eventTree->Branch("trk_MVA1", &m_trk_MVA1);
0601     eventTree->Branch("trk_matchtp_pdgid", &m_trk_matchtp_pdgid);
0602     eventTree->Branch("trk_matchtp_pt", &m_trk_matchtp_pt);
0603     eventTree->Branch("trk_matchtp_eta", &m_trk_matchtp_eta);
0604     eventTree->Branch("trk_matchtp_phi", &m_trk_matchtp_phi);
0605     eventTree->Branch("trk_matchtp_z0", &m_trk_matchtp_z0);
0606     eventTree->Branch("trk_matchtp_dxy", &m_trk_matchtp_dxy);
0607     eventTree->Branch("trk_matchtp_d0", &m_trk_matchtp_d0);
0608     if (TrackingInJets) {
0609       eventTree->Branch("trk_injet", &m_trk_injet);
0610       eventTree->Branch("trk_injet_highpt", &m_trk_injet_highpt);
0611       eventTree->Branch("trk_injet_vhighpt", &m_trk_injet_vhighpt);
0612     }
0613   }
0614
0615   eventTree->Branch("tp_pt", &m_tp_pt);
0616   eventTree->Branch("tp_eta", &m_tp_eta);
0617   eventTree->Branch("tp_phi", &m_tp_phi);
0618   eventTree->Branch("tp_dxy", &m_tp_dxy);
0619   eventTree->Branch("tp_d0", &m_tp_d0);
0620   eventTree->Branch("tp_z0", &m_tp_z0);
0621   eventTree->Branch("tp_d0_prod", &m_tp_d0_prod);
0622   eventTree->Branch("tp_z0_prod", &m_tp_z0_prod);
0623   eventTree->Branch("tp_pdgid", &m_tp_pdgid);
0624   eventTree->Branch("tp_nmatch", &m_tp_nmatch);
0625   eventTree->Branch("tp_nstub", &m_tp_nstub);
0626   eventTree->Branch("tp_eventid", &m_tp_eventid);
0627   eventTree->Branch("tp_charge", &m_tp_charge);
0628   if (TrackingInJets) {
0629     eventTree->Branch("tp_injet", &m_tp_injet);
0630     eventTree->Branch("tp_injet_highpt", &m_tp_injet_highpt);
0631     eventTree->Branch("tp_injet_vhighpt", &m_tp_injet_vhighpt);
0632   }
0633
0634   eventTree->Branch("matchtrk_pt", &m_matchtrk_pt);
0635   eventTree->Branch("matchtrk_eta", &m_matchtrk_eta);
0636   eventTree->Branch("matchtrk_phi", &m_matchtrk_phi);
0637   eventTree->Branch("matchtrk_z0", &m_matchtrk_z0);
0638   eventTree->Branch("matchtrk_d0", &m_matchtrk_d0);
0639   eventTree->Branch("matchtrk_chi2", &m_matchtrk_chi2);
0640   eventTree->Branch("matchtrk_chi2_dof", &m_matchtrk_chi2_dof);
0641   eventTree->Branch("matchtrk_chi2rphi", &m_matchtrk_chi2rphi);
0642   eventTree->Branch("matchtrk_chi2rphi_dof", &m_matchtrk_chi2rphi_dof);
0643   eventTree->Branch("matchtrk_chi2rz", &m_matchtrk_chi2rz);
0644   eventTree->Branch("matchtrk_chi2rz_dof", &m_matchtrk_chi2rz_dof);
0645   eventTree->Branch("matchtrk_bendchi2", &m_matchtrk_bendchi2);
0646   eventTree->Branch("matchtrk_MVA1", &m_matchtrk_MVA1);
0647   eventTree->Branch("matchtrk_nstub", &m_matchtrk_nstub);
0648   eventTree->Branch("matchtrk_lhits", &m_matchtrk_lhits);
0649   eventTree->Branch("matchtrk_dhits", &m_matchtrk_dhits);
0650   eventTree->Branch("matchtrk_seed", &m_matchtrk_seed);
0651   eventTree->Branch("matchtrk_hitpattern", &m_matchtrk_hitpattern);
0652   if (TrackingInJets) {
0653     eventTree->Branch("matchtrk_injet", &m_matchtrk_injet);
0654     eventTree->Branch("matchtrk_injet_highpt", &m_matchtrk_injet_highpt);
0655     eventTree->Branch("matchtrk_injet_vhighpt", &m_matchtrk_injet_vhighpt);
0656   }
0657
0658   if (SaveStubs) {
0659     eventTree->Branch("allstub_x", &m_allstub_x);
0660     eventTree->Branch("allstub_y", &m_allstub_y);
0661     eventTree->Branch("allstub_z", &m_allstub_z);
0662
0663     eventTree->Branch("allstub_isBarrel", &m_allstub_isBarrel);
0664     eventTree->Branch("allstub_layer", &m_allstub_layer);
0665     eventTree->Branch("allstub_isPSmodule", &m_allstub_isPSmodule);
0666     eventTree->Branch("allstub_isTiltedBarrel", &m_allstub_isTiltedBarrel);
0667
0668     eventTree->Branch("allstub_trigDisplace", &m_allstub_trigDisplace);
0669     eventTree->Branch("allstub_trigOffset", &m_allstub_trigOffset);
0670     eventTree->Branch("allstub_trigPos", &m_allstub_trigPos);
0671     eventTree->Branch("allstub_trigBend", &m_allstub_trigBend);
0672
0673     eventTree->Branch("allstub_matchTP_pdgid", &m_allstub_matchTP_pdgid);
0674     eventTree->Branch("allstub_matchTP_pt", &m_allstub_matchTP_pt);
0675     eventTree->Branch("allstub_matchTP_eta", &m_allstub_matchTP_eta);
0676     eventTree->Branch("allstub_matchTP_phi", &m_allstub_matchTP_phi);
0677
0678     eventTree->Branch("allstub_genuine", &m_allstub_genuine);
0679   }
0680
0681   if (TrackingInJets) {
0682     eventTree->Branch("jet_eta", &m_jet_eta);
0683     eventTree->Branch("jet_phi", &m_jet_phi);
0684     eventTree->Branch("jet_pt", &m_jet_pt);
0685     eventTree->Branch("jet_tp_sumpt", &m_jet_tp_sumpt);
0686     eventTree->Branch("jet_trk_sumpt", &m_jet_trk_sumpt);
0687     eventTree->Branch("jet_matchtrk_sumpt", &m_jet_matchtrk_sumpt);
0688   }
0689 }
0690
0691 //////////
0692 // ANALYZE
0693 void L1TrackNtupleMaker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0694   if (not available_)
0695     return;  // No ROOT file open.
0696
0697   if (!(MyProcess == 13 || MyProcess == 11 || MyProcess == 211 || MyProcess == 6 || MyProcess == 15 ||
0698         MyProcess == 1)) {
0699     edm::LogVerbatim("Tracklet") << "The specified MyProcess is invalid! Exiting...";
0700     return;
0701   }
0702
0703   if (!(L1Tk_nPar == 4 || L1Tk_nPar == 5)) {
0704     edm::LogVerbatim("Tracklet") << "Invalid number of track parameters, specified L1Tk_nPar == " << L1Tk_nPar
0705                                  << " but only 4/5 are valid options! Exiting...";
0706     return;
0707   }
0708
0709   // clear variables
0710   if (SaveAllTracks) {
0711     m_trk_pt->clear();
0712     m_trk_eta->clear();
0713     m_trk_phi->clear();
0714     m_trk_d0->clear();
0715     m_trk_z0->clear();
0716     m_trk_chi2->clear();
0717     m_trk_chi2_dof->clear();
0718     m_trk_chi2rphi->clear();
0719     m_trk_chi2rphi_dof->clear();
0720     m_trk_chi2rz->clear();
0721     m_trk_chi2rz_dof->clear();
0722     m_trk_bendchi2->clear();
0723     m_trk_nstub->clear();
0724     m_trk_lhits->clear();
0725     m_trk_dhits->clear();
0726     m_trk_seed->clear();
0727     m_trk_hitpattern->clear();
0728     m_trk_lhits_hitpattern->clear();
0729     m_trk_dhits_hitpattern->clear();
0730     m_trk_nPSstub_hitpattern->clear();
0731     m_trk_n2Sstub_hitpattern->clear();
0732     m_trk_nLostPSstub_hitpattern->clear();
0733     m_trk_nLost2Sstub_hitpattern->clear();
0734     m_trk_nLoststub_V1_hitpattern->clear();
0735     m_trk_nLoststub_V2_hitpattern->clear();
0736     m_trk_phiSector->clear();
0737     m_trk_etaSector->clear();
0738     m_trk_genuine->clear();
0739     m_trk_loose->clear();
0740     m_trk_unknown->clear();
0741     m_trk_combinatoric->clear();
0742     m_trk_fake->clear();
0743     m_trk_MVA1->clear();
0744     m_trk_matchtp_pdgid->clear();
0745     m_trk_matchtp_pt->clear();
0746     m_trk_matchtp_eta->clear();
0747     m_trk_matchtp_phi->clear();
0748     m_trk_matchtp_z0->clear();
0749     m_trk_matchtp_dxy->clear();
0750     m_trk_matchtp_d0->clear();
0751     m_trk_injet->clear();
0752     m_trk_injet_highpt->clear();
0753     m_trk_injet_vhighpt->clear();
0754   }
0755
0756   m_tp_pt->clear();
0757   m_tp_eta->clear();
0758   m_tp_phi->clear();
0759   m_tp_dxy->clear();
0760   m_tp_d0->clear();
0761   m_tp_z0->clear();
0762   m_tp_d0_prod->clear();
0763   m_tp_z0_prod->clear();
0764   m_tp_pdgid->clear();
0765   m_tp_nmatch->clear();
0766   m_tp_nstub->clear();
0767   m_tp_eventid->clear();
0768   m_tp_charge->clear();
0769   m_tp_injet->clear();
0770   m_tp_injet_highpt->clear();
0771   m_tp_injet_vhighpt->clear();
0772
0773   m_matchtrk_pt->clear();
0774   m_matchtrk_eta->clear();
0775   m_matchtrk_phi->clear();
0776   m_matchtrk_z0->clear();
0777   m_matchtrk_d0->clear();
0778   m_matchtrk_chi2->clear();
0779   m_matchtrk_chi2_dof->clear();
0780   m_matchtrk_chi2rphi->clear();
0781   m_matchtrk_chi2rphi_dof->clear();
0782   m_matchtrk_chi2rz->clear();
0783   m_matchtrk_chi2rz_dof->clear();
0784   m_matchtrk_bendchi2->clear();
0785   m_matchtrk_MVA1->clear();
0786   m_matchtrk_nstub->clear();
0787   m_matchtrk_lhits->clear();
0788   m_matchtrk_dhits->clear();
0789   m_matchtrk_seed->clear();
0790   m_matchtrk_hitpattern->clear();
0791   m_matchtrk_injet->clear();
0792   m_matchtrk_injet_highpt->clear();
0793   m_matchtrk_injet_vhighpt->clear();
0794
0795   if (SaveStubs) {
0796     m_allstub_x->clear();
0797     m_allstub_y->clear();
0798     m_allstub_z->clear();
0799
0800     m_allstub_isBarrel->clear();
0801     m_allstub_layer->clear();
0802     m_allstub_isPSmodule->clear();
0803     m_allstub_isTiltedBarrel->clear();
0804
0805     m_allstub_trigDisplace->clear();
0806     m_allstub_trigOffset->clear();
0807     m_allstub_trigPos->clear();
0808     m_allstub_trigBend->clear();
0809
0810     m_allstub_matchTP_pdgid->clear();
0811     m_allstub_matchTP_pt->clear();
0812     m_allstub_matchTP_eta->clear();
0813     m_allstub_matchTP_phi->clear();
0814
0815     m_allstub_genuine->clear();
0816   }
0817
0818   m_jet_eta->clear();
0819   m_jet_phi->clear();
0820   m_jet_pt->clear();
0821   m_jet_tp_sumpt->clear();
0822   m_jet_trk_sumpt->clear();
0823   m_jet_matchtrk_sumpt->clear();
0824
0825   // -----------------------------------------------------------------------------------------------
0826   // retrieve various containers
0827   // -----------------------------------------------------------------------------------------------
0828
0829   // L1 tracks
0830   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > TTTrackHandle;
0831   iEvent.getByToken(ttTrackToken_, TTTrackHandle);
0832
0833   // L1 stubs
0834   edm::Handle<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> > > TTStubHandle;
0835   if (SaveStubs)
0836     iEvent.getByToken(ttStubToken_, TTStubHandle);
0837
0838   // MC truth association maps
0839   edm::Handle<TTClusterAssociationMap<Ref_Phase2TrackerDigi_> > MCTruthTTClusterHandle;
0840   iEvent.getByToken(ttClusterMCTruthToken_, MCTruthTTClusterHandle);
0841   edm::Handle<TTStubAssociationMap<Ref_Phase2TrackerDigi_> > MCTruthTTStubHandle;
0842   iEvent.getByToken(ttStubMCTruthToken_, MCTruthTTStubHandle);
0843   edm::Handle<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> > MCTruthTTTrackHandle;
0844   iEvent.getByToken(ttTrackMCTruthToken_, MCTruthTTTrackHandle);
0845
0846   // tracking particles
0847   edm::Handle<std::vector<TrackingParticle> > TrackingParticleHandle;
0848   edm::Handle<std::vector<TrackingVertex> > TrackingVertexHandle;
0849   iEvent.getByToken(TrackingParticleToken_, TrackingParticleHandle);
0850   //iEvent.getByToken(TrackingVertexToken_, TrackingVertexHandle);
0851
0852   // -----------------------------------------------------------------------------------------------
0853   // more for TTStubs
0854   edm::ESHandle<TrackerGeometry> tGeomHandle = iSetup.getHandle(getTokenTrackerGeom_);
0855
0856   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(getTokenTrackerTopo_);
0857
0858   edm::ESHandle<MagneticField> bFieldHandle = iSetup.getHandle(getTokenBField_);
0859
0860   edm::ESHandle<hph::Setup> hphHandle = iSetup.getHandle(getTokenHPHSetup_);
0861
0862   const TrackerTopology* const tTopo = tTopoHandle.product();
0863   const TrackerGeometry* const theTrackerGeom = tGeomHandle.product();
0864   const hph::Setup* hphSetup = hphHandle.product();
0865
0866   // ----------------------------------------------------------------------------------------------
0867   // loop over L1 stubs
0868   // ----------------------------------------------------------------------------------------------
0869
0870   if (SaveStubs) {
0871     for (auto gd = theTrackerGeom->dets().begin(); gd != theTrackerGeom->dets().end(); gd++) {
0872       DetId detid = (*gd)->geographicalId();
0873       if (detid.subdetId() != StripSubdetector::TOB && detid.subdetId() != StripSubdetector::TID)
0874         continue;
0875       if (!tTopo->isLower(detid))
0876         continue;                              // loop on the stacks: choose the lower arbitrarily
0877       DetId stackDetid = tTopo->stack(detid);  // Stub module detid
0878
0879       if (TTStubHandle->find(stackDetid) == TTStubHandle->end())
0880         continue;
0881
0882       // Get the DetSets of the Clusters
0883       edmNew::DetSet<TTStub<Ref_Phase2TrackerDigi_> > stubs = (*TTStubHandle)[stackDetid];
0884       const GeomDetUnit* det0 = theTrackerGeom->idToDetUnit(detid);
0885       const auto* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(det0);
0886       const PixelTopology* topol = dynamic_cast<const PixelTopology*>(&(theGeomDet->specificTopology()));
0887
0888       // loop over stubs
0889       for (auto stubIter = stubs.begin(); stubIter != stubs.end(); ++stubIter) {
0890         edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > tempStubPtr =
0891             edmNew::makeRefTo(TTStubHandle, stubIter);
0892
0893         int isBarrel = 0;
0894         int layer = -999999;
0895         if (detid.subdetId() == StripSubdetector::TOB) {
0896           isBarrel = 1;
0897           layer = static_cast<int>(tTopo->layer(detid));
0898         } else if (detid.subdetId() == StripSubdetector::TID) {
0899           isBarrel = 0;
0900           layer = static_cast<int>(tTopo->layer(detid));
0901         } else {
0902           edm::LogVerbatim("Tracklet") << "WARNING -- neither TOB or TID stub, shouldn't happen...";
0903           layer = -1;
0904         }
0905
0906         int isPSmodule = 0;
0907         if (topol->nrows() == 960)
0908           isPSmodule = 1;
0909
0910         const unsigned int tobSide = tTopo->tobSide(detid);  // nonBarrel = 0, tiltedMinus = 1, tiltedPlus = 2, flat = 3
0911         int isTiltedBarrel = 0;
0912         if (isBarrel == 1 && (tobSide == 1 || tobSide == 2))
0913           isTiltedBarrel = 1;
0914
0915         MeasurementPoint coords = tempStubPtr->clusterRef(0)->findAverageLocalCoordinatesCentered();
0916         LocalPoint clustlp = topol->localPosition(coords);
0917         GlobalPoint posStub = theGeomDet->surface().toGlobal(clustlp);
0918
0919         double tmp_stub_x = posStub.x();
0920         double tmp_stub_y = posStub.y();
0921         double tmp_stub_z = posStub.z();
0922
0923         float trigDisplace = tempStubPtr->rawBend();
0924         float trigOffset = tempStubPtr->bendOffset();
0925         float trigPos = tempStubPtr->innerClusterPosition();
0926         float trigBend = tempStubPtr->bendFE();
0927
0928         m_allstub_x->push_back(tmp_stub_x);
0929         m_allstub_y->push_back(tmp_stub_y);
0930         m_allstub_z->push_back(tmp_stub_z);
0931
0932         m_allstub_isBarrel->push_back(isBarrel);
0933         m_allstub_layer->push_back(layer);
0934         m_allstub_isPSmodule->push_back(isPSmodule);
0935         m_allstub_isTiltedBarrel->push_back(isTiltedBarrel);
0936
0937         m_allstub_trigDisplace->push_back(trigDisplace);
0938         m_allstub_trigOffset->push_back(trigOffset);
0939         m_allstub_trigPos->push_back(trigPos);
0940         m_allstub_trigBend->push_back(trigBend);
0941
0942         // matched to tracking particle?
0943         edm::Ptr<TrackingParticle> my_tp = MCTruthTTStubHandle->findTrackingParticlePtr(tempStubPtr);
0944
0945         int myTP_pdgid = -999;
0946         float myTP_pt = -999;
0947         float myTP_eta = -999;
0948         float myTP_phi = -999;
0949
0950         if (my_tp.isNull() == false) {
0951           int tmp_eventid = my_tp->eventId().event();
0952
0953           if (tmp_eventid > 0)
0954             continue;  // this means stub from pileup track
0955
0956           myTP_pdgid = my_tp->pdgId();
0957           myTP_pt = my_tp->p4().pt();
0958           myTP_eta = my_tp->p4().eta();
0959           myTP_phi = my_tp->p4().phi();
0960         }
0961
0962         m_allstub_matchTP_pdgid->push_back(myTP_pdgid);
0963         m_allstub_matchTP_pt->push_back(myTP_pt);
0964         m_allstub_matchTP_eta->push_back(myTP_eta);
0965         m_allstub_matchTP_phi->push_back(myTP_phi);
0966
0967         int tmp_stub_genuine = 0;
0968         if (MCTruthTTStubHandle->isGenuine(tempStubPtr))
0969           tmp_stub_genuine = 1;
0970
0971         m_allstub_genuine->push_back(tmp_stub_genuine);
0972       }
0973     }
0974   }
0975
0976   // ----------------------------------------------------------------------------------------------
0977   // tracking in jets
0978   // ----------------------------------------------------------------------------------------------
0979
0980   std::vector<math::XYZTLorentzVector> v_jets;
0981   std::vector<int> v_jets_highpt;
0982   std::vector<int> v_jets_vhighpt;
0983
0984   if (TrackingInJets) {
0985     // gen jets
0986     if (DebugMode)
0987       edm::LogVerbatim("Tracklet") << "get genjets";
0988     edm::Handle<std::vector<reco::GenJet> > GenJetHandle;
0989     iEvent.getByToken(GenJetToken_, GenJetHandle);
0990
0991     if (GenJetHandle.isValid()) {
0992       if (DebugMode)
0993         edm::LogVerbatim("Tracklet") << "loop over genjets";
0994       std::vector<reco::GenJet>::const_iterator iterGenJet;
0995       for (iterGenJet = GenJetHandle->begin(); iterGenJet != GenJetHandle->end(); ++iterGenJet) {
0996         reco::GenJet myJet = reco::GenJet(*iterGenJet);
0997
0998         if (myJet.pt() < 30.0)
0999           continue;
1000         if (std::abs(myJet.eta()) > 2.5)
1001           continue;
1002
1003         if (DebugMode)
1004           edm::LogVerbatim("Tracklet") << "genjet pt = " << myJet.pt() << ", eta = " << myJet.eta();
1005
1006         bool ishighpt = false;
1007         bool isveryhighpt = false;
1008         if (myJet.pt() > 100.0)
1009           ishighpt = true;
1010         if (myJet.pt() > 200.0)
1011           isveryhighpt = true;
1012
1013         const math::XYZTLorentzVector& jetP4 = myJet.p4();
1014         v_jets.push_back(jetP4);
1015         if (ishighpt)
1016           v_jets_highpt.push_back(1);
1017         else
1018           v_jets_highpt.push_back(0);
1019         if (isveryhighpt)
1020           v_jets_vhighpt.push_back(1);
1021         else
1022           v_jets_vhighpt.push_back(0);
1023
1024       }  // end loop over genjets
1025     }  // end isValid
1026
1027   }  // end TrackingInJets
1028
1029   const int NJETS = 10;
1030   float jets_tp_sumpt[NJETS] = {0};        //sum pt of TPs with dR<0.4 of jet
1031   float jets_matchtrk_sumpt[NJETS] = {0};  //sum pt of tracks matched to TP with dR<0.4 of jet
1032   float jets_trk_sumpt[NJETS] = {0};       //sum pt of all tracks with dR<0.4 of jet
1033
1034   // ----------------------------------------------------------------------------------------------
1035   // loop over L1 tracks
1036   // ----------------------------------------------------------------------------------------------
1037
1038   if (SaveAllTracks) {
1039     if (DebugMode) {
1040       edm::LogVerbatim("Tracklet") << "\n Loop over L1 tracks!";
1041       edm::LogVerbatim("Tracklet") << "\n Looking at " << L1Tk_nPar << "-parameter tracks!";
1042     }
1043
1044     int this_l1track = 0;
1045     std::vector<TTTrack<Ref_Phase2TrackerDigi_> >::const_iterator iterL1Track;
1046     for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
1047       edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_> > l1track_ptr(TTTrackHandle, this_l1track);
1048       this_l1track++;
1049
1050       float tmp_trk_pt = iterL1Track->momentum().perp();
1051       float tmp_trk_eta = iterL1Track->momentum().eta();
1052       float tmp_trk_phi = iterL1Track->momentum().phi();
1053       float tmp_trk_z0 = iterL1Track->z0();  //cm
1054       float tmp_trk_tanL = iterL1Track->tanL();
1055       bool usingNewKF = hphSetup->useNewKF();
1056       if (usingNewKF) {
1057         // Skip crazy tracks to avoid crash (as NewKF applies no cuts to kill them).
1058         constexpr float crazy_z0_cut = 30.;  // Cut to kill any crazy tracks found by New KF (which applies no cuts)
1059         if (fabs(tmp_trk_z0) > crazy_z0_cut)
1060           continue;
1061       }
1062
1063       int tmp_trk_hitpattern = 0;
1064       tmp_trk_hitpattern = (int)iterL1Track->hitPattern();
1065       hph::HitPatternHelper hph(hphSetup, tmp_trk_hitpattern, tmp_trk_tanL, tmp_trk_z0);
1066       std::vector<int> hitpattern_expanded_binary = hph.binary();
1067       int tmp_trk_lhits_hitpattern = 0;
1068       int tmp_trk_dhits_hitpattern = 0;
1069       for (int i = 0; i < (int)hitpattern_expanded_binary.size(); i++) {
1070         if (hitpattern_expanded_binary[i]) {
1071           if (i < 6) {
1072             tmp_trk_lhits_hitpattern += pow(10, i);
1073           } else {
1074             tmp_trk_dhits_hitpattern += pow(10, i - 6);
1075           }
1076         }
1077       }
1078       int tmp_trk_nPSstub_hitpattern = hph.numPS();
1079       int tmp_trk_n2Sstub_hitpattern = hph.num2S();
1080       int tmp_trk_nLostPSstub_hitpattern = hph.numMissingPS();
1081       int tmp_trk_nLost2Sstub_hitpattern = hph.numMissing2S();
1082       int tmp_trk_nLoststub_V1_hitpattern = hph.numMissingInterior1();
1083       int tmp_trk_nLoststub_V2_hitpattern = hph.numMissingInterior2();
1084
1085       float tmp_trk_d0 = -999;
1086       if (L1Tk_nPar == 5) {
1087         float tmp_trk_x0 = iterL1Track->POCA().x();
1088         float tmp_trk_y0 = iterL1Track->POCA().y();
1089         tmp_trk_d0 = tmp_trk_x0 * sin(tmp_trk_phi) - tmp_trk_y0 * cos(tmp_trk_phi);
1090       }
1091
1092       float tmp_trk_chi2 = iterL1Track->chi2();
1093       float tmp_trk_chi2rphi = iterL1Track->chi2XY();
1094       float tmp_trk_chi2rz = iterL1Track->chi2Z();
1095       float tmp_trk_bendchi2 = iterL1Track->stubPtConsistency();
1096       float tmp_trk_MVA1 = iterL1Track->trkMVA1();
1097
1098       std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
1099           stubRefs = iterL1Track->getStubRefs();
1100       int tmp_trk_nstub = (int)stubRefs.size();
1101       int ndof = 2 * tmp_trk_nstub - L1Tk_nPar;
1102       int ndofrphi = tmp_trk_nstub - L1Tk_nPar + 2;
1103       int ndofrz = tmp_trk_nstub - 2;
1104       float tmp_trk_chi2_dof = (float)tmp_trk_chi2 / ndof;
1105       float tmp_trk_chi2rphi_dof = (float)tmp_trk_chi2rphi / ndofrphi;
1106       float tmp_trk_chi2rz_dof = (float)tmp_trk_chi2rz / ndofrz;
1107
1108       int tmp_trk_seed = 0;
1109       tmp_trk_seed = (int)iterL1Track->trackSeedType();
1110
1111       unsigned int tmp_trk_phiSector = iterL1Track->phiSector();
1112       int tmp_trk_etaSector = hph.etaSector();
1113
1114       // ----------------------------------------------------------------------------------------------
1115       // loop over stubs on tracks
1116
1117       //float tmp_trk_bend_chi2 = 0;
1118       int tmp_trk_dhits = 0;
1119       int tmp_trk_lhits = 0;
1120
1121       if (true) {
1122         // loop over stubs
1123         for (int is = 0; is < tmp_trk_nstub; is++) {
1124           //detID of stub
1125           DetId detIdStub = theTrackerGeom->idToDet((stubRefs.at(is)->clusterRef(0))->getDetId())->geographicalId();
1126
1127           MeasurementPoint coords = stubRefs.at(is)->clusterRef(0)->findAverageLocalCoordinatesCentered();
1128           const GeomDet* theGeomDet = theTrackerGeom->idToDet(detIdStub);
1129           Global3DPoint posStub = theGeomDet->surface().toGlobal(theGeomDet->topology().localPosition(coords));
1130
1131           double x = posStub.x();
1132           double y = posStub.y();
1133           double z = posStub.z();
1134
1135           int layer = -999999;
1136           if (detIdStub.subdetId() == StripSubdetector::TOB) {
1137             layer = static_cast<int>(tTopo->layer(detIdStub));
1138             if (DebugMode)
1139               edm::LogVerbatim("Tracklet")
1140                   << "   stub in layer " << layer << " at position x y z = " << x << " " << y << " " << z;
1141             tmp_trk_lhits += pow(10, layer - 1);
1142           } else if (detIdStub.subdetId() == StripSubdetector::TID) {
1143             layer = static_cast<int>(tTopo->layer(detIdStub));
1144             if (DebugMode)
1145               edm::LogVerbatim("Tracklet")
1146                   << "   stub in disk " << layer << " at position x y z = " << x << " " << y << " " << z;
1147             tmp_trk_dhits += pow(10, layer - 1);
1148           }
1149
1150         }  //end loop over stubs
1151       }
1152       // ----------------------------------------------------------------------------------------------
1153
1154       int tmp_trk_genuine = 0;
1155       int tmp_trk_loose = 0;
1156       int tmp_trk_unknown = 0;
1157       int tmp_trk_combinatoric = 0;
1158       if (MCTruthTTTrackHandle->isLooselyGenuine(l1track_ptr))
1159         tmp_trk_loose = 1;
1160       if (MCTruthTTTrackHandle->isGenuine(l1track_ptr))
1161         tmp_trk_genuine = 1;
1162       if (MCTruthTTTrackHandle->isUnknown(l1track_ptr))
1163         tmp_trk_unknown = 1;
1164       if (MCTruthTTTrackHandle->isCombinatoric(l1track_ptr))
1165         tmp_trk_combinatoric = 1;
1166
1167       if (DebugMode) {
1168         edm::LogVerbatim("Tracklet") << "L1 track,"
1169                                      << " pt: " << tmp_trk_pt << " eta: " << tmp_trk_eta << " phi: " << tmp_trk_phi
1170                                      << " z0: " << tmp_trk_z0 << " chi2: " << tmp_trk_chi2
1171                                      << " chi2rphi: " << tmp_trk_chi2rphi << " chi2rz: " << tmp_trk_chi2rz
1172                                      << " nstub: " << tmp_trk_nstub;
1173         if (tmp_trk_genuine)
1174           edm::LogVerbatim("Tracklet") << "    (is genuine)";
1175         if (tmp_trk_unknown)
1176           edm::LogVerbatim("Tracklet") << "    (is unknown)";
1177         if (tmp_trk_combinatoric)
1178           edm::LogVerbatim("Tracklet") << "    (is combinatoric)";
1179       }
1180
1181       m_trk_pt->push_back(tmp_trk_pt);
1182       m_trk_eta->push_back(tmp_trk_eta);
1183       m_trk_phi->push_back(tmp_trk_phi);
1184       m_trk_z0->push_back(tmp_trk_z0);
1185       if (L1Tk_nPar == 5)
1186         m_trk_d0->push_back(tmp_trk_d0);
1187       else
1188         m_trk_d0->push_back(999.);
1189       m_trk_chi2->push_back(tmp_trk_chi2);
1190       m_trk_chi2_dof->push_back(tmp_trk_chi2_dof);
1191       m_trk_chi2rphi->push_back(tmp_trk_chi2rphi);
1192       m_trk_chi2rphi_dof->push_back(tmp_trk_chi2rphi_dof);
1193       m_trk_chi2rz->push_back(tmp_trk_chi2rz);
1194       m_trk_chi2rz_dof->push_back(tmp_trk_chi2rz_dof);
1195       m_trk_bendchi2->push_back(tmp_trk_bendchi2);
1196       m_trk_MVA1->push_back(tmp_trk_MVA1);
1197       m_trk_nstub->push_back(tmp_trk_nstub);
1198       m_trk_dhits->push_back(tmp_trk_dhits);
1199       m_trk_lhits->push_back(tmp_trk_lhits);
1200       m_trk_seed->push_back(tmp_trk_seed);
1201       m_trk_hitpattern->push_back(tmp_trk_hitpattern);
1202       m_trk_lhits_hitpattern->push_back(tmp_trk_lhits_hitpattern);
1203       m_trk_dhits_hitpattern->push_back(tmp_trk_dhits_hitpattern);
1204       m_trk_nPSstub_hitpattern->push_back(tmp_trk_nPSstub_hitpattern);
1205       m_trk_n2Sstub_hitpattern->push_back(tmp_trk_n2Sstub_hitpattern);
1206       m_trk_nLostPSstub_hitpattern->push_back(tmp_trk_nLostPSstub_hitpattern);
1207       m_trk_nLost2Sstub_hitpattern->push_back(tmp_trk_nLost2Sstub_hitpattern);
1208       m_trk_nLoststub_V1_hitpattern->push_back(tmp_trk_nLoststub_V1_hitpattern);
1209       m_trk_nLoststub_V2_hitpattern->push_back(tmp_trk_nLoststub_V2_hitpattern);
1210       m_trk_phiSector->push_back(tmp_trk_phiSector);
1211       m_trk_etaSector->push_back(tmp_trk_etaSector);
1212       m_trk_genuine->push_back(tmp_trk_genuine);
1213       m_trk_loose->push_back(tmp_trk_loose);
1214       m_trk_unknown->push_back(tmp_trk_unknown);
1215       m_trk_combinatoric->push_back(tmp_trk_combinatoric);
1216
1217       // ----------------------------------------------------------------------------------------------
1218       // for studying the fake rate
1219       // ----------------------------------------------------------------------------------------------
1220
1221       edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(l1track_ptr);
1222
1223       int myFake = 0;
1224
1225       int tmp_matchtp_pdgid = -999;
1226       float tmp_matchtp_pt = -999;
1227       float tmp_matchtp_eta = -999;
1228       float tmp_matchtp_phi = -999;
1229       float tmp_matchtp_z0 = -999;
1230       float tmp_matchtp_dxy = -999;
1231       float tmp_matchtp_d0 = -999;
1232
1233       if (my_tp.isNull())
1234         myFake = 0;
1235       else {
1236         int tmp_eventid = my_tp->eventId().event();
1237
1238         if (tmp_eventid > 0)
1239           myFake = 2;
1240         else
1241           myFake = 1;
1242
1243         tmp_matchtp_pdgid = my_tp->pdgId();
1244         tmp_matchtp_pt = my_tp->pt();
1245         tmp_matchtp_eta = my_tp->eta();
1246         tmp_matchtp_phi = my_tp->phi();
1247
1248         float tmp_matchtp_vz = my_tp->vz();
1249         float tmp_matchtp_vx = my_tp->vx();
1250         float tmp_matchtp_vy = my_tp->vy();
1251         tmp_matchtp_dxy = sqrt(tmp_matchtp_vx * tmp_matchtp_vx + tmp_matchtp_vy * tmp_matchtp_vy);
1252
1253         // ----------------------------------------------------------------------------------------------
1254         // get d0/z0 propagated back to the IP
1255
1256         float tmp_matchtp_t = 1.0 / tan(2.0 * atan(exp(-tmp_matchtp_eta)));
1257
1258         float delx = -tmp_matchtp_vx;
1259         float dely = -tmp_matchtp_vy;
1260
1261         float b_field = bFieldHandle.product()->inTesla(GlobalPoint(0, 0, 0)).z();
1262         float c_converted = CLHEP::c_light / 1.0E5;
1263         float r2_inv = my_tp->charge() * c_converted * b_field / tmp_matchtp_pt / 2.0;
1264
1265         float tmp_matchtp_x0p = delx - (1. / (2. * r2_inv) * sin(tmp_matchtp_phi));
1266         float tmp_matchtp_y0p = dely + (1. / (2. * r2_inv) * cos(tmp_matchtp_phi));
1267         float tmp_matchtp_rp = sqrt(tmp_matchtp_x0p * tmp_matchtp_x0p + tmp_matchtp_y0p * tmp_matchtp_y0p);
1268         tmp_matchtp_d0 = my_tp->charge() * tmp_matchtp_rp - (1. / (2. * r2_inv));
1269
1270         static double pi = M_PI;
1271         float delphi = tmp_matchtp_phi - atan2(-r2_inv * tmp_matchtp_x0p, r2_inv * tmp_matchtp_y0p);
1272         if (delphi < -pi)
1273           delphi += 2.0 * pi;
1274         if (delphi > pi)
1275           delphi -= 2.0 * pi;
1276         tmp_matchtp_z0 = tmp_matchtp_vz + tmp_matchtp_t * delphi / (2.0 * r2_inv);
1277         // ----------------------------------------------------------------------------------------------
1278
1279         if (DebugMode) {
1280           edm::LogVerbatim("Tracklet") << "TP matched to track has pt = " << my_tp->p4().pt()
1281                                        << " eta = " << my_tp->momentum().eta() << " phi = " << my_tp->momentum().phi()
1282                                        << " z0 = " << my_tp->vertex().z() << " pdgid = " << my_tp->pdgId()
1283                                        << " dxy = " << tmp_matchtp_dxy;
1284         }
1285       }
1286
1287       m_trk_fake->push_back(myFake);
1288
1289       m_trk_matchtp_pdgid->push_back(tmp_matchtp_pdgid);
1290       m_trk_matchtp_pt->push_back(tmp_matchtp_pt);
1291       m_trk_matchtp_eta->push_back(tmp_matchtp_eta);
1292       m_trk_matchtp_phi->push_back(tmp_matchtp_phi);
1293       m_trk_matchtp_z0->push_back(tmp_matchtp_z0);
1294       m_trk_matchtp_dxy->push_back(tmp_matchtp_dxy);
1295       m_trk_matchtp_d0->push_back(tmp_matchtp_d0);
1296
1297       // ----------------------------------------------------------------------------------------------
1298       // for tracking in jets
1299       // ----------------------------------------------------------------------------------------------
1300
1301       if (TrackingInJets) {
1302         if (DebugMode)
1303           edm::LogVerbatim("Tracklet") << "doing tracking in jets now";
1304
1305         int InJet = 0;
1306         int InJetHighpt = 0;
1307         int InJetVeryHighpt = 0;
1308
1309         for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1310           float deta = tmp_trk_eta - (v_jets.at(ij)).eta();
1311           float dphi = tmp_trk_phi - (v_jets.at(ij)).phi();
1312           while (dphi > 3.14159)
1313             dphi = std::abs(2 * 3.14159 - dphi);
1314           float dR = sqrt(deta * deta + dphi * dphi);
1315
1316           if (dR < 0.4) {
1317             InJet = 1;
1318             if (v_jets_highpt.at(ij) == 1)
1319               InJetHighpt = 1;
1320             if (v_jets_vhighpt.at(ij) == 1)
1321               InJetVeryHighpt = 1;
1322             if (ij < NJETS)
1323               jets_trk_sumpt[ij] += tmp_trk_pt;
1324           }
1325         }
1326
1327         m_trk_injet->push_back(InJet);
1328         m_trk_injet_highpt->push_back(InJetHighpt);
1329         m_trk_injet_vhighpt->push_back(InJetVeryHighpt);
1330
1331       }  //end tracking in jets
1332
1333     }  //end track loop
1334
1335   }  //end if SaveAllTracks
1336
1337   // ----------------------------------------------------------------------------------------------
1338   // loop over tracking particles
1339   // ----------------------------------------------------------------------------------------------
1340
1341   if (DebugMode)
1342     edm::LogVerbatim("Tracklet") << "\n Loop over tracking particles!";
1343
1344   int this_tp = 0;
1345   std::vector<TrackingParticle>::const_iterator iterTP;
1346   for (iterTP = TrackingParticleHandle->begin(); iterTP != TrackingParticleHandle->end(); ++iterTP) {
1347     edm::Ptr<TrackingParticle> tp_ptr(TrackingParticleHandle, this_tp);
1348     this_tp++;
1349
1350     int tmp_eventid = iterTP->eventId().event();
1351     if (MyProcess != 1 && tmp_eventid > 0)
1352       continue;  //only care about tracking particles from the primary interaction (except for MyProcess==1, i.e. looking at all TPs)
1353
1354     float tmp_tp_pt = iterTP->pt();
1355     float tmp_tp_charge = iterTP->charge();
1356     float tmp_tp_eta = iterTP->eta();
1357
1358     if (tmp_tp_pt < TP_minPt)  // Save CPU by applying these cuts here.
1359       continue;
1360     if (tmp_tp_charge == 0.)
1361       continue;
1362     if (std::abs(tmp_tp_eta) > TP_maxEta)
1363       continue;
1364
1365     float tmp_tp_phi = iterTP->phi();
1366     float tmp_tp_vz = iterTP->vz();
1367     float tmp_tp_vx = iterTP->vx();
1368     float tmp_tp_vy = iterTP->vy();
1369     int tmp_tp_pdgid = iterTP->pdgId();
1370     float tmp_tp_z0_prod = tmp_tp_vz;
1371     float tmp_tp_d0_prod = tmp_tp_vx * sin(tmp_tp_phi) - tmp_tp_vy * cos(tmp_tp_phi);
1372
1373     // ----------------------------------------------------------------------------------------------
1374     // get d0/z0 propagated back to the IP
1375
1376     float tmp_tp_t = 1.0 / tan(2.0 * atan(exp(-tmp_tp_eta)));
1377
1378     float delx = -tmp_tp_vx;
1379     float dely = -tmp_tp_vy;
1380
1381     float b_field = bFieldHandle.product()->inTesla(GlobalPoint(0, 0, 0)).z();
1382     float c_converted = CLHEP::c_light / 1.0E5;
1383     float r2_inv = tmp_tp_charge * c_converted * b_field / tmp_tp_pt / 2.0;
1384
1385     float tmp_tp_x0p = delx - (1. / (2. * r2_inv) * sin(tmp_tp_phi));
1386     float tmp_tp_y0p = dely + (1. / (2. * r2_inv) * cos(tmp_tp_phi));
1387     float tmp_tp_rp = sqrt(tmp_tp_x0p * tmp_tp_x0p + tmp_tp_y0p * tmp_tp_y0p);
1388     float tmp_tp_d0 = tmp_tp_charge * tmp_tp_rp - (1. / (2. * r2_inv));
1389
1390     static double pi = M_PI;
1391     float delphi = tmp_tp_phi - atan2(-r2_inv * tmp_tp_x0p, r2_inv * tmp_tp_y0p);
1392     if (delphi < -pi)
1393       delphi += 2.0 * pi;
1394     if (delphi > pi)
1395       delphi -= 2.0 * pi;
1396     float tmp_tp_z0 = tmp_tp_vz + tmp_tp_t * delphi / (2.0 * r2_inv);
1397     // ----------------------------------------------------------------------------------------------
1398
1399     if (MyProcess == 13 && abs(tmp_tp_pdgid) != 13)
1400       continue;
1401     if (MyProcess == 11 && abs(tmp_tp_pdgid) != 11)
1402       continue;
1403     if ((MyProcess == 6 || MyProcess == 15 || MyProcess == 211) && abs(tmp_tp_pdgid) != 211)
1404       continue;
1405
1406     if (std::abs(tmp_tp_z0) > TP_maxZ0)
1407       continue;
1408
1409     // for pions in ttbar, only consider TPs coming from near the IP!
1410     float dxy = sqrt(tmp_tp_vx * tmp_tp_vx + tmp_tp_vy * tmp_tp_vy);
1411     float tmp_tp_dxy = dxy;
1412     if (MyProcess == 6 && (dxy > 1.0))
1413       continue;
1414
1415     if (DebugMode)
1416       edm::LogVerbatim("Tracklet") << "Tracking particle, pt: " << tmp_tp_pt << " eta: " << tmp_tp_eta
1417                                    << " phi: " << tmp_tp_phi << " z0: " << tmp_tp_z0 << " d0: " << tmp_tp_d0
1418                                    << " z_prod: " << tmp_tp_z0_prod << " d_prod: " << tmp_tp_d0_prod
1419                                    << " pdgid: " << tmp_tp_pdgid << " eventID: " << iterTP->eventId().event()
1420                                    << " ttclusters " << MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).size()
1421                                    << " ttstubs " << MCTruthTTStubHandle->findTTStubRefs(tp_ptr).size() << " tttracks "
1422                                    << MCTruthTTTrackHandle->findTTTrackPtrs(tp_ptr).size();
1423
1424     // ----------------------------------------------------------------------------------------------
1425     // only consider TPs associated with >= 1 cluster, or >= X stubs, or have stubs in >= X layers (configurable options)
1426
1427     if (MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).empty()) {
1428       if (DebugMode)
1429         edm::LogVerbatim("Tracklet") << "No matching TTClusters for TP, continuing...";
1430       continue;
1431     }
1432
1433     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
1434         theStubRefs = MCTruthTTStubHandle->findTTStubRefs(tp_ptr);
1435     int nStubTP = (int)theStubRefs.size();
1436
1437     // how many layers/disks have stubs?
1438     int hasStubInLayer[11] = {0};
1439     for (auto& theStubRef : theStubRefs) {
1440       DetId detid(theStubRef->getDetId());
1441
1442       int layer = -1;
1443       if (detid.subdetId() == StripSubdetector::TOB) {
1444         layer = static_cast<int>(tTopo->layer(detid)) - 1;  //fill in array as entries 0-5
1445       } else if (detid.subdetId() == StripSubdetector::TID) {
1446         layer = static_cast<int>(tTopo->layer(detid)) + 5;  //fill in array as entries 6-10
1447       }
1448
1449       //bool isPS = (theTrackerGeom->getDetectorType(detid)==TrackerGeometry::ModuleType::Ph2PSP);
1450
1451       //treat genuine stubs separately (==2 is genuine, ==1 is not)
1452       if (MCTruthTTStubHandle->findTrackingParticlePtr(theStubRef).isNull() && hasStubInLayer[layer] < 2)
1453         hasStubInLayer[layer] = 1;
1454       else
1455         hasStubInLayer[layer] = 2;
1456     }
1457
1458     int nStubLayerTP = 0;
1459     int nStubLayerTP_g = 0;
1460     for (int isum : hasStubInLayer) {
1461       if (isum >= 1)
1462         nStubLayerTP += 1;
1463       if (isum == 2)
1464         nStubLayerTP_g += 1;
1465     }
1466
1467     if (DebugMode)
1468       edm::LogVerbatim("Tracklet") << "TP is associated with " << nStubTP << " stubs, and has stubs in " << nStubLayerTP
1469                                    << " different layers/disks, and has GENUINE stubs in " << nStubLayerTP_g
1470                                    << " layers ";
1471
1472     if (TP_minNStub > 0) {
1473       if (DebugMode)
1474         edm::LogVerbatim("Tracklet") << "Only consider TPs with >= " << TP_minNStub << " stubs";
1475       if (nStubTP < TP_minNStub) {
1476         if (DebugMode)
1477           edm::LogVerbatim("Tracklet") << "TP fails minimum nbr stubs requirement! Continuing...";
1478         continue;
1479       }
1480     }
1481     if (TP_minNStubLayer > 0) {
1482       if (DebugMode)
1483         edm::LogVerbatim("Tracklet") << "Only consider TPs with stubs in >= " << TP_minNStubLayer << " layers/disks";
1484       if (nStubLayerTP < TP_minNStubLayer) {
1485         if (DebugMode)
1486           edm::LogVerbatim("Tracklet") << "TP fails stubs in minimum nbr of layers/disks requirement! Continuing...";
1487         continue;
1488       }
1489     }
1490
1491     // ----------------------------------------------------------------------------------------------
1492     // look for L1 tracks matched to the tracking particle
1493
1494     std::vector<edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_> > > matchedTracks =
1495         MCTruthTTTrackHandle->findTTTrackPtrs(tp_ptr);
1496
1497     int nMatch = 0;
1498     int i_track = -1;
1499     float i_chi2dof = 99999;
1500
1501     if (!matchedTracks.empty()) {
1502       if (DebugMode && (matchedTracks.size() > 1))
1503         edm::LogVerbatim("Tracklet") << "TrackingParticle has more than one matched L1 track!";
1504
1505       // ----------------------------------------------------------------------------------------------
1506       // loop over matched L1 tracks
1507       // here, "match" means tracks that can be associated to a TrackingParticle with at least one hit of at least one of its clusters
1508       // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SLHCTrackerTriggerSWTools#MC_truth_for_TTTrack
1509
1510       for (int it = 0; it < (int)matchedTracks.size(); it++) {
1511         bool tmp_trk_genuine = false;
1512         bool tmp_trk_loosegenuine = false;
1513         if (MCTruthTTTrackHandle->isGenuine(matchedTracks.at(it)))
1514           tmp_trk_genuine = true;
1515         if (MCTruthTTTrackHandle->isLooselyGenuine(matchedTracks.at(it)))
1516           tmp_trk_loosegenuine = true;
1517         if (!tmp_trk_loosegenuine)
1518           continue;
1519
1520         if (DebugMode) {
1521           if (MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it)).isNull()) {
1522             edm::LogVerbatim("Tracklet") << "track matched to TP is NOT uniquely matched to a TP";
1523           } else {
1524             edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it));
1525             edm::LogVerbatim("Tracklet") << "TP matched to track matched to TP ... tp pt = " << my_tp->p4().pt()
1526                                          << " eta = " << my_tp->momentum().eta() << " phi = " << my_tp->momentum().phi()
1527                                          << " z0 = " << my_tp->vertex().z();
1528           }
1529           edm::LogVerbatim("Tracklet") << "   ... matched L1 track has pt = " << matchedTracks.at(it)->momentum().perp()
1530                                        << " eta = " << matchedTracks.at(it)->momentum().eta()
1531                                        << " phi = " << matchedTracks.at(it)->momentum().phi()
1532                                        << " chi2 = " << matchedTracks.at(it)->chi2()
1533                                        << " consistency = " << matchedTracks.at(it)->stubPtConsistency()
1534                                        << " z0 = " << matchedTracks.at(it)->z0()
1535                                        << " nstub = " << matchedTracks.at(it)->getStubRefs().size();
1536           if (tmp_trk_genuine)
1537             edm::LogVerbatim("Tracklet") << "    (genuine!) ";
1538           if (tmp_trk_loosegenuine)
1539             edm::LogVerbatim("Tracklet") << "    (loose genuine!) ";
1540         }
1541
1542         // ----------------------------------------------------------------------------------------------
1543         // further require L1 track to be (loosely) genuine, that there is only one TP matched to the track
1544         // + have >= L1Tk_minNStub stubs for it to be a valid match (only relevant is your track collection
1545         // e.g. stores 3-stub tracks but at plot level you require >= 4 stubs (--> tracklet case)
1546
1547         std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
1548             stubRefs = matchedTracks.at(it)->getStubRefs();
1549         int tmp_trk_nstub = stubRefs.size();
1550
1551         if (tmp_trk_nstub < L1Tk_minNStub)
1552           continue;
1553
1554         /*
1555     // PS stubs
1556     int tmp_trk_nPSstub = 0;
1557     for (int is=0; is<tmp_trk_nstub; is++) {
1558       DetId detIdStub = theTrackerGeom->idToDet( (stubRefs.at(is)->clusterRef(0))->getDetId() )->geographicalId();
1559       DetId stackDetid = tTopo->stack(detIdStub);
1560       bool isPS = (theTrackerGeom->getDetectorType(stackDetid)==TrackerGeometry::ModuleType::Ph2PSP);
1561       if (isPS) tmp_trk_nPSstub++;
1562     }
1563     */
1564
1565         float dmatch_pt = 999;
1566         float dmatch_eta = 999;
1567         float dmatch_phi = 999;
1568         int match_id = 999;
1569
1570         edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it));
1571         dmatch_pt = std::abs(my_tp->p4().pt() - tmp_tp_pt);
1572         dmatch_eta = std::abs(my_tp->p4().eta() - tmp_tp_eta);
1573         dmatch_phi = std::abs(my_tp->p4().phi() - tmp_tp_phi);
1574         match_id = my_tp->pdgId();
1575
1576         float tmp_trk_chi2dof = (matchedTracks.at(it)->chi2()) / (2 * tmp_trk_nstub - L1Tk_nPar);
1577
1578         // ensure that track is uniquely matched to the TP we are looking at!
1579         if (dmatch_pt < 0.1 && dmatch_eta < 0.1 && dmatch_phi < 0.1 && tmp_tp_pdgid == match_id && tmp_trk_genuine) {
1580           nMatch++;
1581           if (i_track < 0 || tmp_trk_chi2dof < i_chi2dof) {
1582             i_track = it;
1583             i_chi2dof = tmp_trk_chi2dof;
1584           }
1585         }
1586
1587       }  // end loop over matched L1 tracks
1588
1589     }  // end has at least 1 matched L1 track
1590     // ----------------------------------------------------------------------------------------------
1591
1592     float tmp_matchtrk_pt = -999;
1593     float tmp_matchtrk_eta = -999;
1594     float tmp_matchtrk_phi = -999;
1595     float tmp_matchtrk_z0 = -999;
1596     float tmp_matchtrk_d0 = -999;
1597     float tmp_matchtrk_chi2 = -999;
1598     float tmp_matchtrk_chi2_dof = -999;
1599     float tmp_matchtrk_chi2rphi = -999;
1600     float tmp_matchtrk_chi2rphi_dof = -999;
1601     float tmp_matchtrk_chi2rz = -999;
1602     float tmp_matchtrk_chi2rz_dof = -999;
1603     float tmp_matchtrk_bendchi2 = -999;
1604     float tmp_matchtrk_MVA1 = -999;
1605     int tmp_matchtrk_nstub = -999;
1606     int tmp_matchtrk_dhits = -999;
1607     int tmp_matchtrk_lhits = -999;
1608     int tmp_matchtrk_seed = -999;
1609     int tmp_matchtrk_hitpattern = -999;
1610
1611     if (nMatch > 1 && DebugMode)
1612       edm::LogVerbatim("Tracklet") << "WARNING *** 2 or more matches to genuine L1 tracks ***";
1613
1614     if (nMatch > 0) {
1615       tmp_matchtrk_pt = matchedTracks.at(i_track)->momentum().perp();
1616       tmp_matchtrk_eta = matchedTracks.at(i_track)->momentum().eta();
1617       tmp_matchtrk_phi = matchedTracks.at(i_track)->momentum().phi();
1618       tmp_matchtrk_z0 = matchedTracks.at(i_track)->z0();
1619
1620       if (L1Tk_nPar == 5) {
1621         float tmp_matchtrk_x0 = matchedTracks.at(i_track)->POCA().x();
1622         float tmp_matchtrk_y0 = matchedTracks.at(i_track)->POCA().y();
1623         tmp_matchtrk_d0 = tmp_matchtrk_x0 * sin(tmp_matchtrk_phi) - tmp_matchtrk_y0 * cos(tmp_matchtrk_phi);
1624       }
1625
1626       tmp_matchtrk_chi2 = matchedTracks.at(i_track)->chi2();
1627       tmp_matchtrk_chi2rphi = matchedTracks.at(i_track)->chi2XY();
1628       tmp_matchtrk_chi2rz = matchedTracks.at(i_track)->chi2Z();
1629       tmp_matchtrk_bendchi2 = matchedTracks.at(i_track)->stubPtConsistency();
1630       tmp_matchtrk_MVA1 = matchedTracks.at(i_track)->trkMVA1();
1631       tmp_matchtrk_nstub = (int)matchedTracks.at(i_track)->getStubRefs().size();
1632       tmp_matchtrk_seed = (int)matchedTracks.at(i_track)->trackSeedType();
1633       tmp_matchtrk_hitpattern = (int)matchedTracks.at(i_track)->hitPattern();
1634
1635       int ndof = 2 * tmp_matchtrk_nstub - L1Tk_nPar;
1636       int ndofrphi = tmp_matchtrk_nstub - L1Tk_nPar + 2;
1637       int ndofrz = tmp_matchtrk_nstub - 2;
1638       tmp_matchtrk_chi2_dof = (float)tmp_matchtrk_chi2 / ndof;
1639       tmp_matchtrk_chi2rphi_dof = (float)tmp_matchtrk_chi2rphi / ndofrphi;
1640       tmp_matchtrk_chi2rz_dof = (float)tmp_matchtrk_chi2rz / ndofrz;
1641
1642       // ------------------------------------------------------------------------------------------
1643
1644       //float tmp_matchtrk_bend_chi2 = 0;
1645
1646       tmp_matchtrk_dhits = 0;
1647       tmp_matchtrk_lhits = 0;
1648
1649       std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
1650           stubRefs = matchedTracks.at(i_track)->getStubRefs();
1651       int tmp_nstub = stubRefs.size();
1652
1653       for (int is = 0; is < tmp_nstub; is++) {
1654         DetId detIdStub = theTrackerGeom->idToDet((stubRefs.at(is)->clusterRef(0))->getDetId())->geographicalId();
1655         /*
1656     MeasurementPoint coords = stubRefs.at(is)->clusterRef(0)->findAverageLocalCoordinatesCentered();
1657     const GeomDet* theGeomDet = theTrackerGeom->idToDet(detIdStub);
1658     Global3DPoint posStub = theGeomDet->surface().toGlobal( theGeomDet->topology().localPosition(coords) );
1659     */
1660
1661         int layer = -999999;
1662         if (detIdStub.subdetId() == StripSubdetector::TOB) {
1663           layer = static_cast<int>(tTopo->layer(detIdStub));
1664           tmp_matchtrk_lhits += pow(10, layer - 1);
1665         } else if (detIdStub.subdetId() == StripSubdetector::TID) {
1666           layer = static_cast<int>(tTopo->layer(detIdStub));
1667           tmp_matchtrk_dhits += pow(10, layer - 1);
1668         }
1669
1670         // ------------------------------------------------------------------------------------------
1671       }
1672     }
1673
1674     m_tp_pt->push_back(tmp_tp_pt);
1675     m_tp_eta->push_back(tmp_tp_eta);
1676     m_tp_phi->push_back(tmp_tp_phi);
1677     m_tp_dxy->push_back(tmp_tp_dxy);
1678     m_tp_z0->push_back(tmp_tp_z0);
1679     m_tp_d0->push_back(tmp_tp_d0);
1680     m_tp_z0_prod->push_back(tmp_tp_z0_prod);
1681     m_tp_d0_prod->push_back(tmp_tp_d0_prod);
1682     m_tp_pdgid->push_back(tmp_tp_pdgid);
1683     m_tp_nmatch->push_back(nMatch);
1684     m_tp_nstub->push_back(nStubTP);
1685     m_tp_eventid->push_back(tmp_eventid);
1686     m_tp_charge->push_back(tmp_tp_charge);
1687
1688     m_matchtrk_pt->push_back(tmp_matchtrk_pt);
1689     m_matchtrk_eta->push_back(tmp_matchtrk_eta);
1690     m_matchtrk_phi->push_back(tmp_matchtrk_phi);
1691     m_matchtrk_z0->push_back(tmp_matchtrk_z0);
1692     m_matchtrk_d0->push_back(tmp_matchtrk_d0);
1693     m_matchtrk_chi2->push_back(tmp_matchtrk_chi2);
1694     m_matchtrk_chi2rphi->push_back(tmp_matchtrk_chi2rphi);
1695     m_matchtrk_chi2rz->push_back(tmp_matchtrk_chi2rz);
1696     m_matchtrk_bendchi2->push_back(tmp_matchtrk_bendchi2);
1697     m_matchtrk_MVA1->push_back(tmp_matchtrk_MVA1);
1698     m_matchtrk_nstub->push_back(tmp_matchtrk_nstub);
1699     m_matchtrk_dhits->push_back(tmp_matchtrk_dhits);
1700     m_matchtrk_lhits->push_back(tmp_matchtrk_lhits);
1701     m_matchtrk_seed->push_back(tmp_matchtrk_seed);
1702     m_matchtrk_hitpattern->push_back(tmp_matchtrk_hitpattern);
1703     m_matchtrk_chi2_dof->push_back(tmp_matchtrk_chi2_dof);
1704     m_matchtrk_chi2rphi_dof->push_back(tmp_matchtrk_chi2rphi_dof);
1705     m_matchtrk_chi2rz_dof->push_back(tmp_matchtrk_chi2rz_dof);
1706
1707     // ----------------------------------------------------------------------------------------------
1708     // for tracking in jets
1709     // ----------------------------------------------------------------------------------------------
1710
1711     if (TrackingInJets) {
1712       if (DebugMode)
1713         edm::LogVerbatim("Tracklet") << "check if TP/matched track is within jet";
1714
1715       int tp_InJet = 0;
1716       int matchtrk_InJet = 0;
1717       int tp_InJetHighpt = 0;
1718       int matchtrk_InJetHighpt = 0;
1719       int tp_InJetVeryHighpt = 0;
1720       int matchtrk_InJetVeryHighpt = 0;
1721
1722       for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1723         float deta = tmp_tp_eta - (v_jets.at(ij)).eta();
1724         float dphi = tmp_tp_phi - (v_jets.at(ij)).phi();
1725         while (dphi > 3.14159)
1726           dphi = std::abs(2 * 3.14159 - dphi);
1727         float dR = sqrt(deta * deta + dphi * dphi);
1728         if (dR < 0.4) {
1729           tp_InJet = 1;
1730           if (v_jets_highpt.at(ij) == 1)
1731             tp_InJetHighpt = 1;
1732           if (v_jets_vhighpt.at(ij) == 1)
1733             tp_InJetVeryHighpt = 1;
1734           if (ij < NJETS)
1735             jets_tp_sumpt[ij] += tmp_tp_pt;
1736         }
1737
1738         if (nMatch > 0) {
1739           deta = tmp_matchtrk_eta - (v_jets.at(ij)).eta();
1740           dphi = tmp_matchtrk_phi - (v_jets.at(ij)).phi();
1741           while (dphi > 3.14159)
1742             dphi = std::abs(2 * 3.14159 - dphi);
1743           dR = sqrt(deta * deta + dphi * dphi);
1744           if (dR < 0.4) {
1745             matchtrk_InJet = 1;
1746             if (v_jets_highpt.at(ij) == 1)
1747               matchtrk_InJetHighpt = 1;
1748             if (v_jets_vhighpt.at(ij) == 1)
1749               matchtrk_InJetVeryHighpt = 1;
1750             if (ij < NJETS)
1751               jets_matchtrk_sumpt[ij] += tmp_matchtrk_pt;
1752           }
1753         }
1754       }
1755
1756       m_tp_injet->push_back(tp_InJet);
1757       m_tp_injet_highpt->push_back(tp_InJetHighpt);
1758       m_tp_injet_vhighpt->push_back(tp_InJetVeryHighpt);
1759       m_matchtrk_injet->push_back(matchtrk_InJet);
1760       m_matchtrk_injet_highpt->push_back(matchtrk_InJetHighpt);
1761       m_matchtrk_injet_vhighpt->push_back(matchtrk_InJetVeryHighpt);
1762
1763     }  //end TrackingInJets
1764
1765   }  //end loop tracking particles
1766
1767   if (TrackingInJets) {
1768     for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1769       if (ij < NJETS) {
1770         m_jet_eta->push_back((v_jets.at(ij)).eta());
1771         m_jet_phi->push_back((v_jets.at(ij)).phi());
1772         m_jet_pt->push_back((v_jets.at(ij)).pt());
1773         m_jet_tp_sumpt->push_back(jets_tp_sumpt[ij]);
1774         m_jet_trk_sumpt->push_back(jets_trk_sumpt[ij]);
1775         m_jet_matchtrk_sumpt->push_back(jets_matchtrk_sumpt[ij]);
1776       }
1777     }
1778   }
1779
1780   eventTree->Fill();
1781
1782 }  // end of analyze()
1783
1784 ///////////////////////////
1785 // DEFINE THIS AS A PLUG-IN
1786 DEFINE_FWK_MODULE(L1TrackNtupleMaker);