Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-07 22:33:42 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!