TrackingNtuple.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    NtupleDump/TrackingNtuple
// Class:      TrackingNtuple
//
//
// Original Author:  Giuseppe Cerati
//         Created:  Tue, 25 Aug 2015 13:22:49 GMT
//
//


// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CommonTools/Utils/interface/DynArray.h"
#include "DataFormats/Provenance/interface/ProductID.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/transform.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"

#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/RZLine.h"
#include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
#include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h"

#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"

#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"

#include "SimTracker/Records/interface/TrackAssociatorRecord.h"
#include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
#include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
#include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
#include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
#include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
#include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"

#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"

#include "Validation/RecoTrack/interface/trackFromSeedFitFailed.h"

#include <set>
#include <map>
#include <unordered_set>
#include <unordered_map>

#include "TTree.h"

/*
todo: 
add refitted hit position after track/seed fit
add local angle, path length!
*/

namespace {
  std::string subdetstring(int subdet) {
    switch(subdet) {
    case StripSubdetector::TIB:         return "- TIB";
    case StripSubdetector::TOB:         return "- TOB";
    case StripSubdetector::TEC:         return "- TEC";
    case StripSubdetector::TID:         return "- TID";
    case PixelSubdetector::PixelBarrel: return "- PixBar";
    case PixelSubdetector::PixelEndcap: return "- PixFwd";
    default:                            return "UNKNOWN TRACKER HIT TYPE";
    }
  }

  struct ProductIDSetPrinter {
    ProductIDSetPrinter(const std::set<edm::ProductID>& set): set_(set) {}

    void print(std::ostream& os) const {
      for(const auto& item: set_) {
        os << item << " ";
      }
    }

    const std::set<edm::ProductID>& set_;
  };
  std::ostream& operator<<(std::ostream& os, const ProductIDSetPrinter& o) {
    o.print(os);
    return os;
  }
  template <typename T>
  struct ProductIDMapPrinter {
    ProductIDMapPrinter(const std::map<edm::ProductID, T>& map): map_(map) {}

    void print(std::ostream& os) const {
      for(const auto& item: map_) {
        os << item.first << " ";
      }
    }

    const std::map<edm::ProductID, T>& map_;
  };
  template <typename T>
  auto make_ProductIDMapPrinter(const std::map<edm::ProductID, T>& map) {
    return ProductIDMapPrinter<T>(map);
  }
  template <typename T>
  std::ostream& operator<<(std::ostream& os, const ProductIDMapPrinter<T>& o) {
    o.print(os);
    return os;
  }

  template <typename T>
  struct VectorPrinter {
    VectorPrinter(const std::vector<T>& vec): vec_(vec) {}

    void print(std::ostream& os) const {
      for(const auto& item: vec_) {
        os << item << " ";
      }
    }

    const std::vector<T>& vec_;
  };
  template <typename T>
  auto make_VectorPrinter(const std::vector<T>& vec) {
    return VectorPrinter<T>(vec);
  }
  template <typename T>
  std::ostream& operator<<(std::ostream& os, const VectorPrinter<T>& o) {
    o.print(os);
    return os;
  }

  void checkProductID(const std::set<edm::ProductID>& set, const edm::ProductID& id, const char *name) {
    if(set.find(id) == set.end())
      throw cms::Exception("Configuration") << "Got " << name << " with a hit with ProductID " << id
                                            << " which does not match to the set of ProductID's for the hits: "
                                            << ProductIDSetPrinter(set)
                                            << ". Usually this is caused by a wrong hit collection in the configuration.";
  }

  template <typename SimLink, typename Func>
  void forEachMatchedSimLink(const edm::DetSet<SimLink>& digiSimLinks, uint32_t channel, Func func) {
    for(const auto& link: digiSimLinks) {
      if(link.channel() == channel) {
        func(link);
      }
    }
  }

  std::map<unsigned int, double>  chargeFraction(const SiPixelCluster& cluster, const DetId& detId,
                                                 const edm::DetSetVector<PixelDigiSimLink>& digiSimLink) {
    std::map<unsigned int, double> simTrackIdToAdc;

    auto idetset = digiSimLink.find(detId);
    if(idetset == digiSimLink.end())
      return simTrackIdToAdc;

    double adcSum = 0;
    PixelDigiSimLink found;
    for(int iPix=0; iPix != cluster.size(); ++iPix) {
      const SiPixelCluster::Pixel& pixel = cluster.pixel(iPix);
      adcSum += pixel.adc;
      uint32_t channel = PixelChannelIdentifier::pixelToChannel(pixel.x, pixel.y);
      forEachMatchedSimLink(*idetset, channel, [&](const PixelDigiSimLink& simLink){
          double& adc = simTrackIdToAdc[simLink.SimTrackId()];
          adc += pixel.adc*simLink.fraction();
        });
    }

    for(auto& pair: simTrackIdToAdc) {
      if(adcSum == 0.)
        pair.second = 0.;
      else
        pair.second /= adcSum;
    }

    return simTrackIdToAdc;
  }

  std::map<unsigned int, double> chargeFraction(const SiStripCluster& cluster, const DetId& detId,
                                                const edm::DetSetVector<StripDigiSimLink>& digiSimLink) {
    std::map<unsigned int, double> simTrackIdToAdc;

    auto idetset = digiSimLink.find(detId);
    if(idetset == digiSimLink.end())
      return simTrackIdToAdc;

    double adcSum = 0;
    StripDigiSimLink found;
    int first  = cluster.firstStrip();
    for(size_t i=0; i<cluster.amplitudes().size(); ++i) {
      adcSum += cluster.amplitudes()[i];
      forEachMatchedSimLink(*idetset, first+i, [&](const StripDigiSimLink& simLink){
          double& adc = simTrackIdToAdc[simLink.SimTrackId()];
          adc += cluster.amplitudes()[i]*simLink.fraction();
        });

      for(const auto& pair: simTrackIdToAdc) {
        simTrackIdToAdc[pair.first] = (adcSum != 0. ? pair.second/adcSum : 0.);
      }

    }
    return simTrackIdToAdc;
  }


  std::map<unsigned int, double> chargeFraction(const Phase2TrackerCluster1D& cluster, const DetId& detId,
                                                const edm::DetSetVector<StripDigiSimLink>& digiSimLink) {
    std::map<unsigned int, double> simTrackIdToAdc;
    throw cms::Exception("LogicError") << "Not possible to use StripDigiSimLink with Phase2TrackerCluster1D! ";
    return simTrackIdToAdc;
  }

  //In the OT, there is no measurement of the charge, so no ADC value.
  //Only in the SSA chip (so in PSs) you have one "threshold" flag that tells you if the charge of at least one strip in the cluster exceeded 1.2 MIPs.
  std::map<unsigned int, double> chargeFraction(const Phase2TrackerCluster1D& cluster, const DetId& detId,
                                                const edm::DetSetVector<PixelDigiSimLink>& digiSimLink) {
    std::map<unsigned int, double> simTrackIdToAdc;
    return simTrackIdToAdc;
  }

}

//
// class declaration
//

class TrackingNtuple : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit TrackingNtuple(const edm::ParameterSet&);
  ~TrackingNtuple();

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);


private:
  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;

  void clearVariables();

  // This pattern is copied from QuickTrackAssociatorByHitsImpl. If
  // further needs arises, it wouldn't hurt to abstract it somehow.
  typedef std::unordered_set<reco::RecoToSimCollection::index_type> TrackingParticleRefKeySet;
  typedef std::unordered_map<reco::RecoToSimCollection::index_type, size_t> TrackingParticleRefKeyToIndex;
  typedef TrackingParticleRefKeyToIndex TrackingVertexRefKeyToIndex;
  typedef std::pair<TrackPSimHitRef::key_type, edm::ProductID> SimHitFullKey;
  typedef std::map<SimHitFullKey, size_t> SimHitRefKeyToIndex;

  enum class HitType {
    Pixel = 0,
    Strip = 1,
    Glued = 2,
    Invalid = 3,
    Phase2OT = 4,
    Unknown = 99
  };

  // This gives the "best" classification of a reco hit
  // In case of reco hit mathing to multiple sim, smaller number is
  // considered better
  // To be kept in synch with class HitSimType in ntuple.py
  enum class HitSimType {
    Signal = 0,
    ITPileup = 1,
    OOTPileup = 2,
    Noise = 3,
    Unknown = 99
  };

  struct TPHitIndex {
    TPHitIndex(unsigned int tp=0, unsigned int simHit=0, float to=0, unsigned int id=0): tpKey(tp), simHitIdx(simHit), tof(to), detId(id) {}
    unsigned int tpKey;
    unsigned int simHitIdx;
    float tof;
    unsigned int detId;
  };
  static bool tpHitIndexListLess(const TPHitIndex& i, const TPHitIndex& j) { return (i.tpKey < j.tpKey); }
  static bool tpHitIndexListLessSort(const TPHitIndex& i, const TPHitIndex& j) {
    if(i.tpKey == j.tpKey) {
      if(edm::isNotFinite(i.tof) && edm::isNotFinite(j.tof)) {
        return i.detId < j.detId;
      }
      return i.tof < j.tof; // works as intended if either one is NaN
    }
    return i.tpKey < j.tpKey;
  }

  void fillBeamSpot(const reco::BeamSpot& bs);
  void fillPixelHits(const edm::Event& iEvent,
                     const ClusterTPAssociation& clusterToTPMap,
                     const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                     const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                     const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
                     const TransientTrackingRecHitBuilder& theTTRHBuilder,
                     const TrackerTopology& tTopo,
                     const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                     std::set<edm::ProductID>& hitProductIds
                     );

  void fillStripRphiStereoHits(const edm::Event& iEvent,
                               const ClusterTPAssociation& clusterToTPMap,
                               const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                               const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                               const edm::DetSetVector<StripDigiSimLink>& digiSimLink,
                               const TransientTrackingRecHitBuilder& theTTRHBuilder,
                               const TrackerTopology& tTopo,
                               const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                               std::set<edm::ProductID>& hitProductIds
                               );

  void fillStripMatchedHits(const edm::Event& iEvent,
                            const TransientTrackingRecHitBuilder& theTTRHBuilder,
                            const TrackerTopology& tTopo,
                            std::vector<std::pair<int, int> >& monoStereoClusterList
                            );

  void fillPhase2OTHits(const edm::Event& iEvent,
                        const ClusterTPAssociation& clusterToTPMap,
                        const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                        const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                        const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
                        const TransientTrackingRecHitBuilder& theTTRHBuilder,
                        const TrackerTopology& tTopo,
                        const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                        std::set<edm::ProductID>& hitProductIds
                        );

  void fillSeeds(const edm::Event& iEvent,
                 const TrackingParticleRefVector& tpCollection,
                 const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                 const reco::BeamSpot& bs,
                 const reco::TrackToTrackingParticleAssociator& associatorByHits,
                 const TransientTrackingRecHitBuilder& theTTRHBuilder,
                 const MagneticField *theMF,
                 const std::vector<std::pair<int, int> >& monoStereoClusterList,
                 const std::set<edm::ProductID>& hitProductIds,
                 std::map<edm::ProductID, size_t>& seedToCollIndex
                 );

  void fillTracks(const edm::RefToBaseVector<reco::Track>& tracks,
                  const TrackingParticleRefVector& tpCollection,
                  const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                  const reco::BeamSpot& bs,
                  const reco::TrackToTrackingParticleAssociator& associatorByHits,
                  const TransientTrackingRecHitBuilder& theTTRHBuilder,
                  const TrackerTopology& tTopo,
                  const std::set<edm::ProductID>& hitProductIds,
                  const std::map<edm::ProductID, size_t>& seedToCollIndex
                  );

  void fillSimHits(const TrackerGeometry& tracker,
                   const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                   const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                   const TrackerTopology& tTopo,
                   SimHitRefKeyToIndex& simHitRefKeyToIndex,
                   std::vector<TPHitIndex>& tpHitList);

  void fillTrackingParticles(const edm::Event& iEvent, const edm::EventSetup& iSetup,
                             const edm::RefToBaseVector<reco::Track>& tracks,
                             const TrackingParticleRefVector& tpCollection,
                             const TrackingVertexRefKeyToIndex& tvKeyToIndex,
                             const reco::TrackToTrackingParticleAssociator& associatorByHits,
                             const std::vector<TPHitIndex>& tpHitList
                             );

  void fillVertices(const reco::VertexCollection& vertices,
                    const edm::RefToBaseVector<reco::Track>& tracks);

  void fillTrackingVertices(const TrackingVertexRefVector& trackingVertices,
                            const TrackingParticleRefKeyToIndex& tpKeyToIndex
                            );



  struct SimHitData {
    std::vector<int> matchingSimHit;
    std::vector<float> chargeFraction;
    std::vector<int> bunchCrossing;
    std::vector<int> event;
    HitSimType type = HitSimType::Unknown;
  };

  template <typename SimLink>
  SimHitData matchCluster(const OmniClusterRef& cluster,
                          DetId hitId, int clusterKey,
                          const TransientTrackingRecHit::RecHitPointer& ttrh,
                          const ClusterTPAssociation& clusterToTPMap,
                          const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                          const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                          const edm::DetSetVector<SimLink>& digiSimLinks,
                          const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                          HitType hitType
                          );

  // ----------member data ---------------------------
  std::vector<edm::EDGetTokenT<edm::View<reco::Track> > > seedTokens_;
  edm::EDGetTokenT<edm::View<reco::Track> > trackToken_;
  edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
  edm::EDGetTokenT<TrackingParticleRefVector> trackingParticleRefToken_;
  edm::EDGetTokenT<ClusterTPAssociation> clusterTPMapToken_;
  edm::EDGetTokenT<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitTPMapToken_;
  edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> trackAssociatorToken_;
  edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > pixelSimLinkToken_;
  edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink> > stripSimLinkToken_;
  edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinksToken_;
  bool includeStripHits_, includePhase2OTHits_;
  edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
  edm::EDGetTokenT<SiPixelRecHitCollection> pixelRecHitToken_;
  edm::EDGetTokenT<SiStripRecHit2DCollection> stripRphiRecHitToken_;
  edm::EDGetTokenT<SiStripRecHit2DCollection> stripStereoRecHitToken_;
  edm::EDGetTokenT<SiStripMatchedRecHit2DCollection> stripMatchedRecHitToken_;
  edm::EDGetTokenT<Phase2TrackerRecHit1DCollectionNew> phase2OTRecHitToken_;
  edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
  edm::EDGetTokenT<TrackingVertexCollection> trackingVertexToken_;
  edm::EDGetTokenT<edm::ValueMap<unsigned int> > tpNLayersToken_;
  edm::EDGetTokenT<edm::ValueMap<unsigned int> > tpNPixelLayersToken_;
  edm::EDGetTokenT<edm::ValueMap<unsigned int> > tpNStripStereoLayersToken_;
  std::string builderName_;
  std::string parametersDefinerName_;
  const bool includeSeeds_;
  const bool includeAllHits_;

  TTree* t;
  // event
  edm::RunNumber_t ev_run;
  edm::LuminosityBlockNumber_t ev_lumi;
  edm::EventNumber_t ev_event;

  // tracks
  // (first) index runs through tracks
  std::vector<float> trk_px       ;
  std::vector<float> trk_py       ;
  std::vector<float> trk_pz       ;
  std::vector<float> trk_pt       ;
  std::vector<float> trk_inner_px ;
  std::vector<float> trk_inner_py ;
  std::vector<float> trk_inner_pz ;
  std::vector<float> trk_inner_pt ;
  std::vector<float> trk_outer_px ;
  std::vector<float> trk_outer_py ;
  std::vector<float> trk_outer_pz ;
  std::vector<float> trk_outer_pt ;
  std::vector<float> trk_eta      ;
  std::vector<float> trk_lambda   ;
  std::vector<float> trk_cotTheta ;
  std::vector<float> trk_phi      ;
  std::vector<float> trk_dxy      ;
  std::vector<float> trk_dz       ;
  std::vector<float> trk_ptErr    ;
  std::vector<float> trk_etaErr   ;
  std::vector<float> trk_lambdaErr;
  std::vector<float> trk_phiErr   ;
  std::vector<float> trk_dxyErr   ;
  std::vector<float> trk_dzErr    ;
  std::vector<float> trk_refpoint_x;
  std::vector<float> trk_refpoint_y;
  std::vector<float> trk_refpoint_z;
  std::vector<float> trk_nChi2    ;
  std::vector<int> trk_q       ;
  std::vector<unsigned int> trk_nValid  ;
  std::vector<unsigned int> trk_nInvalid;
  std::vector<unsigned int> trk_nPixel  ;
  std::vector<unsigned int> trk_nStrip  ;
  std::vector<unsigned int> trk_nPixelLay;
  std::vector<unsigned int> trk_nStripLay;
  std::vector<unsigned int> trk_n3DLay  ;
  std::vector<unsigned int> trk_nOuterLost;
  std::vector<unsigned int> trk_nInnerLost;
  std::vector<unsigned int> trk_algo    ;
  std::vector<unsigned int> trk_originalAlgo;
  std::vector<decltype(reco::TrackBase().algoMaskUL())> trk_algoMask;
  std::vector<unsigned int> trk_stopReason;
  std::vector<short> trk_isHP    ;
  std::vector<int> trk_seedIdx ;
  std::vector<int> trk_vtxIdx;
  std::vector<std::vector<float> > trk_shareFrac; // second index runs through matched TrackingParticles
  std::vector<std::vector<int> > trk_simTrkIdx;   // second index runs through matched TrackingParticles
  std::vector<std::vector<int> > trk_hitIdx;      // second index runs through hits
  std::vector<std::vector<int> > trk_hitType;     // second index runs through hits
  // sim tracks
  // (first) index runs through TrackingParticles
  std::vector<int>   sim_event    ;
  std::vector<int>   sim_bunchCrossing;
  std::vector<int>   sim_pdgId    ;
  std::vector<float> sim_px       ;
  std::vector<float> sim_py       ;
  std::vector<float> sim_pz       ;
  std::vector<float> sim_pt       ;
  std::vector<float> sim_eta      ;
  std::vector<float> sim_phi      ;
  std::vector<float> sim_pca_pt   ;
  std::vector<float> sim_pca_eta  ;
  std::vector<float> sim_pca_lambda;
  std::vector<float> sim_pca_cotTheta;
  std::vector<float> sim_pca_phi  ;
  std::vector<float> sim_pca_dxy  ;
  std::vector<float> sim_pca_dz   ;
  std::vector<int> sim_q       ;
  std::vector<unsigned int> sim_nValid  ;
  std::vector<unsigned int> sim_nPixel  ;
  std::vector<unsigned int> sim_nStrip  ;
  std::vector<unsigned int> sim_nLay;
  std::vector<unsigned int> sim_nPixelLay;
  std::vector<unsigned int> sim_n3DLay  ;
  std::vector<std::vector<int> > sim_trkIdx;      // second index runs through matched tracks
  std::vector<std::vector<float> > sim_shareFrac; // second index runs through matched tracks
  std::vector<int> sim_parentVtxIdx;
  std::vector<std::vector<int> > sim_decayVtxIdx; // second index runs through decay vertices
  std::vector<std::vector<int> > sim_simHitIdx;   // second index runs through SimHits
  // pixel hits
  // (first) index runs through hits
  std::vector<short> pix_isBarrel ;
  std::vector<unsigned short> pix_det    ;
  std::vector<unsigned short> pix_lay    ;
  std::vector<unsigned int> pix_detId    ;
  std::vector<std::vector<int> > pix_trkIdx;    // second index runs through tracks containing this hit
  std::vector<std::vector<int> > pix_seeIdx;    // second index runs through seeds containing this hit
  std::vector<std::vector<int> > pix_simHitIdx; // second index runs through SimHits inducing this hit
  std::vector<std::vector<float> > pix_chargeFraction; // second index runs through SimHits inducing this hit
  std::vector<unsigned short> pix_simType;
  std::vector<float> pix_x    ;
  std::vector<float> pix_y    ;
  std::vector<float> pix_z    ;
  std::vector<float> pix_xx   ;
  std::vector<float> pix_xy   ;
  std::vector<float> pix_yy   ;
  std::vector<float> pix_yz   ;
  std::vector<float> pix_zz   ;
  std::vector<float> pix_zx   ;
  std::vector<float> pix_radL ;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
  std::vector<float> pix_bbxi ;
  // strip hits
  // (first) index runs through hits
  std::vector<short> str_isBarrel ;
  std::vector<short> str_isStereo ;
  std::vector<unsigned short> str_det    ;
  std::vector<unsigned short> str_lay    ;
  std::vector<unsigned int> str_detId    ;
  std::vector<std::vector<int> > str_trkIdx;    // second index runs through tracks containing this hit
  std::vector<std::vector<int> > str_seeIdx;    // second index runs through seeds containing this hitw
  std::vector<std::vector<int> > str_simHitIdx; // second index runs through SimHits inducing this hit
  std::vector<std::vector<float> > str_chargeFraction; // second index runs through SimHits inducing this hit
  std::vector<unsigned short> str_simType;
  std::vector<float> str_x    ;
  std::vector<float> str_y    ;
  std::vector<float> str_z    ;
  std::vector<float> str_xx   ;
  std::vector<float> str_xy   ;
  std::vector<float> str_yy   ;
  std::vector<float> str_yz   ;
  std::vector<float> str_zz   ;
  std::vector<float> str_zx   ;
  std::vector<float> str_radL ;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
  std::vector<float> str_bbxi ;
  // strip matched hits
  // (first) index runs through hits
  std::vector<short> glu_isBarrel ;
  std::vector<unsigned int> glu_det      ;
  std::vector<unsigned int> glu_lay      ;
  std::vector<unsigned int> glu_detId    ;
  std::vector<int> glu_monoIdx  ;
  std::vector<int> glu_stereoIdx;
  std::vector<std::vector<int> > glu_seeIdx; // second index runs through seeds containing this hit
  std::vector<float> glu_x    ;
  std::vector<float> glu_y    ;
  std::vector<float> glu_z    ;
  std::vector<float> glu_xx   ;
  std::vector<float> glu_xy   ;
  std::vector<float> glu_yy   ;
  std::vector<float> glu_yz   ;
  std::vector<float> glu_zz   ;
  std::vector<float> glu_zx   ;
  std::vector<float> glu_radL ;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
  std::vector<float> glu_bbxi ;
  // phase2 Outer Tracker hits
  // (first) index runs through hits
  std::vector<short> ph2_isBarrel ;
  std::vector<unsigned short> ph2_det    ;
  std::vector<unsigned short> ph2_lay    ;
  std::vector<unsigned int> ph2_detId    ;
  std::vector<std::vector<int> > ph2_trkIdx;    // second index runs through tracks containing this hit
  std::vector<std::vector<int> > ph2_seeIdx;    // second index runs through seeds containing this hit
  std::vector<std::vector<int> > ph2_simHitIdx; // second index runs through SimHits inducing this hit
  //std::vector<std::vector<float> > ph2_chargeFraction; // Not supported at the moment for Phase2
  std::vector<unsigned short> ph2_simType;
  std::vector<float> ph2_x    ;
  std::vector<float> ph2_y    ;
  std::vector<float> ph2_z    ;
  std::vector<float> ph2_xx   ;
  std::vector<float> ph2_xy   ;
  std::vector<float> ph2_yy   ;
  std::vector<float> ph2_yz   ;
  std::vector<float> ph2_zz   ;
  std::vector<float> ph2_zx   ;
  std::vector<float> ph2_radL ;  //http://cmslxr.fnal.gov/lxr/source/DataFormats/GeometrySurface/interface/MediumProperties.h
  std::vector<float> ph2_bbxi ;

  // invalid (missing/inactive/etc) hits
  // (first) index runs through hits
  std::vector<short> inv_isBarrel;
  std::vector<unsigned short> inv_det;
  std::vector<unsigned short> inv_lay;
  std::vector<unsigned int> inv_detId;
  std::vector<unsigned short> inv_type;
  // sim hits
  // (first) index runs through hits
  std::vector<unsigned short> simhit_det;
  std::vector<unsigned short> simhit_lay;
  std::vector<unsigned int> simhit_detId;
  std::vector<float> simhit_x;
  std::vector<float> simhit_y;
  std::vector<float> simhit_z;
  std::vector<int> simhit_particle;
  std::vector<short> simhit_process;
  std::vector<float> simhit_eloss;
  std::vector<float> simhit_tof;
  //std::vector<unsigned int> simhit_simTrackId; // can be useful for debugging, but not much of general interest
  std::vector<int> simhit_simTrkIdx;
  std::vector<std::vector<int> > simhit_hitIdx;  // second index runs through induced reco hits
  std::vector<std::vector<int> > simhit_hitType; // second index runs through induced reco hits
  // beam spot
  float bsp_x;
  float bsp_y;
  float bsp_z;
  float bsp_sigmax;
  float bsp_sigmay;
  float bsp_sigmaz;
  // seeds
  // (first) index runs through seeds
  std::vector<short> see_fitok     ;
  std::vector<float> see_px       ;
  std::vector<float> see_py       ;
  std::vector<float> see_pz       ;
  std::vector<float> see_pt       ;
  std::vector<float> see_eta      ;
  std::vector<float> see_phi      ;
  std::vector<float> see_dxy      ;
  std::vector<float> see_dz       ;
  std::vector<float> see_ptErr    ;
  std::vector<float> see_etaErr   ;
  std::vector<float> see_phiErr   ;
  std::vector<float> see_dxyErr   ;
  std::vector<float> see_dzErr    ;
  std::vector<float> see_chi2     ;
  std::vector<int> see_q       ;
  std::vector<unsigned int> see_nValid  ;
  std::vector<unsigned int> see_nPixel  ;
  std::vector<unsigned int> see_nGlued  ;
  std::vector<unsigned int> see_nStrip  ;
  std::vector<unsigned int> see_nPhase2OT;
  std::vector<unsigned int> see_algo    ;
  std::vector<int> see_trkIdx;
  std::vector<std::vector<float> > see_shareFrac; // second index runs through matched TrackingParticles
  std::vector<std::vector<int> > see_simTrkIdx;   // second index runs through matched TrackingParticles
  std::vector<std::vector<int> > see_hitIdx;      // second index runs through hits
  std::vector<std::vector<int> > see_hitType;     // second index runs through hits
  //seed algo offset, index runs through iterations
  std::vector<unsigned int> see_offset  ;


  // Vertices
  // (first) index runs through vertices
  std::vector<float> vtx_x;
  std::vector<float> vtx_y;
  std::vector<float> vtx_z;
  std::vector<float> vtx_xErr;
  std::vector<float> vtx_yErr;
  std::vector<float> vtx_zErr;
  std::vector<float> vtx_ndof;
  std::vector<float> vtx_chi2;
  std::vector<short> vtx_fake;
  std::vector<short> vtx_valid;
  std::vector<std::vector<int> > vtx_trkIdx; // second index runs through tracks used in the vertex fit

  // Tracking vertices
  // (first) index runs through TrackingVertices
  std::vector<int>   simvtx_event;
  std::vector<int>   simvtx_bunchCrossing;
  std::vector<unsigned int> simvtx_processType; // only from first SimVertex of TrackingVertex
  std::vector<float> simvtx_x;
  std::vector<float> simvtx_y;
  std::vector<float> simvtx_z;
  std::vector<std::vector<int> > simvtx_sourceSimIdx;   // second index runs through source TrackingParticles
  std::vector<std::vector<int> > simvtx_daughterSimIdx; // second index runs through daughter TrackingParticles
  std::vector<int> simpv_idx;
};

//
// constructors and destructor
//
TrackingNtuple::TrackingNtuple(const edm::ParameterSet& iConfig):
  seedTokens_(edm::vector_transform(iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("seedTracks"), [&](const edm::InputTag& tag) {
        return consumes<edm::View<reco::Track> >(tag);
      })),
  trackToken_(consumes<edm::View<reco::Track> >(iConfig.getUntrackedParameter<edm::InputTag>("tracks"))),
  clusterTPMapToken_(consumes<ClusterTPAssociation>(iConfig.getUntrackedParameter<edm::InputTag>("clusterTPMap"))),
  simHitTPMapToken_(consumes<SimHitTPAssociationProducer::SimHitTPAssociationList>(iConfig.getUntrackedParameter<edm::InputTag>("simHitTPMap"))),
  trackAssociatorToken_(consumes<reco::TrackToTrackingParticleAssociator>(iConfig.getUntrackedParameter<edm::InputTag>("trackAssociator"))),
  pixelSimLinkToken_(consumes<edm::DetSetVector<PixelDigiSimLink> >(iConfig.getUntrackedParameter<edm::InputTag>("pixelDigiSimLink"))),
  stripSimLinkToken_(consumes<edm::DetSetVector<StripDigiSimLink> >(iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink"))),
  siphase2OTSimLinksToken_(consumes<edm::DetSetVector<PixelDigiSimLink> >(iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink"))),
  includeStripHits_(iConfig.getUntrackedParameter<edm::InputTag>("stripDigiSimLink").label() != ""),
  includePhase2OTHits_(iConfig.getUntrackedParameter<edm::InputTag>("phase2OTSimLink").label() != ""),
  beamSpotToken_(consumes<reco::BeamSpot>(iConfig.getUntrackedParameter<edm::InputTag>("beamSpot"))),
  pixelRecHitToken_(consumes<SiPixelRecHitCollection>(iConfig.getUntrackedParameter<edm::InputTag>("pixelRecHits"))),
  stripRphiRecHitToken_(consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripRphiRecHits"))),
  stripStereoRecHitToken_(consumes<SiStripRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripStereoRecHits"))),
  stripMatchedRecHitToken_(consumes<SiStripMatchedRecHit2DCollection>(iConfig.getUntrackedParameter<edm::InputTag>("stripMatchedRecHits"))),
  phase2OTRecHitToken_(consumes<Phase2TrackerRecHit1DCollectionNew>(iConfig.getUntrackedParameter<edm::InputTag>("phase2OTRecHits"))),
  vertexToken_(consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertices"))),
  trackingVertexToken_(consumes<TrackingVertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackingVertices"))),
  tpNLayersToken_(consumes<edm::ValueMap<unsigned int> >(iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNlayers"))),
  tpNPixelLayersToken_(consumes<edm::ValueMap<unsigned int> >(iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNpixellayers"))),
  tpNStripStereoLayersToken_(consumes<edm::ValueMap<unsigned int> >(iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleNstripstereolayers"))),
  builderName_(iConfig.getUntrackedParameter<std::string>("TTRHBuilder")),
  parametersDefinerName_(iConfig.getUntrackedParameter<std::string>("parametersDefiner")),
  includeSeeds_(iConfig.getUntrackedParameter<bool>("includeSeeds")),
  includeAllHits_(iConfig.getUntrackedParameter<bool>("includeAllHits"))
{
  const bool tpRef = iConfig.getUntrackedParameter<bool>("trackingParticlesRef");
  const auto tpTag = iConfig.getUntrackedParameter<edm::InputTag>("trackingParticles");
  if(tpRef) {
    trackingParticleRefToken_ = consumes<TrackingParticleRefVector>(tpTag);
  }
  else {
    trackingParticleToken_ = consumes<TrackingParticleCollection>(tpTag);
  }

  usesResource(TFileService::kSharedResource);
  edm::Service<TFileService> fs;
  t = fs->make<TTree>("tree","tree");

  t->Branch("event"        , &ev_event);
  t->Branch("lumi"         , &ev_lumi);
  t->Branch("run"          , &ev_run);

  //tracks
  t->Branch("trk_px"       , &trk_px);
  t->Branch("trk_py"       , &trk_py);
  t->Branch("trk_pz"       , &trk_pz);
  t->Branch("trk_pt"       , &trk_pt);
  t->Branch("trk_inner_px" , &trk_inner_px);
  t->Branch("trk_inner_py" , &trk_inner_py);
  t->Branch("trk_inner_pz" , &trk_inner_pz);
  t->Branch("trk_inner_pt" , &trk_inner_pt);
  t->Branch("trk_outer_px" , &trk_outer_px);
  t->Branch("trk_outer_py" , &trk_outer_py);
  t->Branch("trk_outer_pz" , &trk_outer_pz);
  t->Branch("trk_outer_pt" , &trk_outer_pt);
  t->Branch("trk_eta"      , &trk_eta);
  t->Branch("trk_lambda"   , &trk_lambda);
  t->Branch("trk_cotTheta" , &trk_cotTheta);
  t->Branch("trk_phi"      , &trk_phi);
  t->Branch("trk_dxy"      , &trk_dxy      );
  t->Branch("trk_dz"       , &trk_dz       );
  t->Branch("trk_ptErr"    , &trk_ptErr    );
  t->Branch("trk_etaErr"   , &trk_etaErr   );
  t->Branch("trk_lambdaErr", &trk_lambdaErr);
  t->Branch("trk_phiErr"   , &trk_phiErr   );
  t->Branch("trk_dxyErr"   , &trk_dxyErr   );
  t->Branch("trk_dzErr"    , &trk_dzErr    );
  t->Branch("trk_refpoint_x", &trk_refpoint_x);
  t->Branch("trk_refpoint_y", &trk_refpoint_y);
  t->Branch("trk_refpoint_z", &trk_refpoint_z);
  t->Branch("trk_nChi2"    , &trk_nChi2);
  t->Branch("trk_q"        , &trk_q);
  t->Branch("trk_nValid"   , &trk_nValid  );
  t->Branch("trk_nInvalid" , &trk_nInvalid);
  t->Branch("trk_nPixel"   , &trk_nPixel  );
  t->Branch("trk_nStrip"   , &trk_nStrip  );
  t->Branch("trk_nPixelLay", &trk_nPixelLay);
  t->Branch("trk_nStripLay", &trk_nStripLay);
  t->Branch("trk_n3DLay"   , &trk_n3DLay  );
  t->Branch("trk_nOuterLost", &trk_nOuterLost  );
  t->Branch("trk_nInnerLost", &trk_nInnerLost  );
  t->Branch("trk_algo"     , &trk_algo    );
  t->Branch("trk_originalAlgo", &trk_originalAlgo);
  t->Branch("trk_algoMask" , &trk_algoMask);
  t->Branch("trk_stopReason", &trk_stopReason);
  t->Branch("trk_isHP"     , &trk_isHP    );
  if(includeSeeds_) {
    t->Branch("trk_seedIdx"  , &trk_seedIdx );
  }
  t->Branch("trk_vtxIdx"   , &trk_vtxIdx);
  t->Branch("trk_shareFrac", &trk_shareFrac);
  t->Branch("trk_simTrkIdx", &trk_simTrkIdx  );
  if(includeAllHits_) {
    t->Branch("trk_hitIdx" , &trk_hitIdx);
    t->Branch("trk_hitType", &trk_hitType);
  }
  //sim tracks
  t->Branch("sim_event"    , &sim_event    );
  t->Branch("sim_bunchCrossing", &sim_bunchCrossing);
  t->Branch("sim_pdgId"    , &sim_pdgId    );
  t->Branch("sim_px"       , &sim_px       );
  t->Branch("sim_py"       , &sim_py       );
  t->Branch("sim_pz"       , &sim_pz       );
  t->Branch("sim_pt"       , &sim_pt       );
  t->Branch("sim_eta"      , &sim_eta      );
  t->Branch("sim_phi"      , &sim_phi      );
  t->Branch("sim_pca_pt"   , &sim_pca_pt  );
  t->Branch("sim_pca_eta"  , &sim_pca_eta  );
  t->Branch("sim_pca_lambda", &sim_pca_lambda);
  t->Branch("sim_pca_cotTheta", &sim_pca_cotTheta);
  t->Branch("sim_pca_phi"  , &sim_pca_phi  );
  t->Branch("sim_pca_dxy"  , &sim_pca_dxy  );
  t->Branch("sim_pca_dz"   , &sim_pca_dz   );
  t->Branch("sim_q"        , &sim_q        );
  t->Branch("sim_nValid"   , &sim_nValid   );
  t->Branch("sim_nPixel"   , &sim_nPixel   );
  t->Branch("sim_nStrip"   , &sim_nStrip   );
  t->Branch("sim_nLay"     , &sim_nLay     );
  t->Branch("sim_nPixelLay", &sim_nPixelLay);
  t->Branch("sim_n3DLay"   , &sim_n3DLay   );
  t->Branch("sim_trkIdx"   , &sim_trkIdx   );
  t->Branch("sim_shareFrac", &sim_shareFrac);
  t->Branch("sim_parentVtxIdx", &sim_parentVtxIdx);
  t->Branch("sim_decayVtxIdx", &sim_decayVtxIdx);
  if(includeAllHits_) {
    t->Branch("sim_simHitIdx" , &sim_simHitIdx );
  }
  if(includeAllHits_) {
    //pixels
    t->Branch("pix_isBarrel"  , &pix_isBarrel );
    t->Branch("pix_det"       , &pix_det      );
    t->Branch("pix_lay"       , &pix_lay      );
    t->Branch("pix_detId"     , &pix_detId    );
    t->Branch("pix_trkIdx"    , &pix_trkIdx   );
    if(includeSeeds_) {
      t->Branch("pix_seeIdx"    , &pix_seeIdx   );
    }
    t->Branch("pix_simHitIdx" , &pix_simHitIdx);
    t->Branch("pix_chargeFraction", &pix_chargeFraction);
    t->Branch("pix_simType", &pix_simType);
    t->Branch("pix_x"     , &pix_x    );
    t->Branch("pix_y"     , &pix_y    );
    t->Branch("pix_z"     , &pix_z    );
    t->Branch("pix_xx"    , &pix_xx   );
    t->Branch("pix_xy"    , &pix_xy   );
    t->Branch("pix_yy"    , &pix_yy   );
    t->Branch("pix_yz"    , &pix_yz   );
    t->Branch("pix_zz"    , &pix_zz   );
    t->Branch("pix_zx"    , &pix_zx   );
    t->Branch("pix_radL"  , &pix_radL );
    t->Branch("pix_bbxi"  , &pix_bbxi );
    t->Branch("pix_bbxi"  , &pix_bbxi );
    //strips
    if(includeStripHits_){
      t->Branch("str_isBarrel"  , &str_isBarrel );
      t->Branch("str_isStereo"  , &str_isStereo );
      t->Branch("str_det"       , &str_det      );
      t->Branch("str_lay"       , &str_lay      );
      t->Branch("str_detId"     , &str_detId    );
      t->Branch("str_trkIdx"    , &str_trkIdx   );
      if(includeSeeds_) {
        t->Branch("str_seeIdx"    , &str_seeIdx   );
      }
      t->Branch("str_simHitIdx" , &str_simHitIdx);
      t->Branch("str_chargeFraction", &str_chargeFraction);
      t->Branch("str_simType", &str_simType);
      t->Branch("str_x"     , &str_x    );
      t->Branch("str_y"     , &str_y    );
      t->Branch("str_z"     , &str_z    );
      t->Branch("str_xx"    , &str_xx   );
      t->Branch("str_xy"    , &str_xy   );
      t->Branch("str_yy"    , &str_yy   );
      t->Branch("str_yz"    , &str_yz   );
      t->Branch("str_zz"    , &str_zz   );
      t->Branch("str_zx"    , &str_zx   );
      t->Branch("str_radL"  , &str_radL );
      t->Branch("str_bbxi"  , &str_bbxi );
      //matched hits
      t->Branch("glu_isBarrel"  , &glu_isBarrel );
      t->Branch("glu_det"       , &glu_det      );
      t->Branch("glu_lay"       , &glu_lay      );
      t->Branch("glu_detId"     , &glu_detId    );
      t->Branch("glu_monoIdx"   , &glu_monoIdx  );
      t->Branch("glu_stereoIdx" , &glu_stereoIdx);
      if(includeSeeds_) {
        t->Branch("glu_seeIdx"    , &glu_seeIdx   );
      }
      t->Branch("glu_x"         , &glu_x        );
      t->Branch("glu_y"         , &glu_y        );
      t->Branch("glu_z"         , &glu_z        );
      t->Branch("glu_xx"        , &glu_xx       );
      t->Branch("glu_xy"        , &glu_xy       );
      t->Branch("glu_yy"        , &glu_yy       );
      t->Branch("glu_yz"        , &glu_yz       );
      t->Branch("glu_zz"        , &glu_zz       );
      t->Branch("glu_zx"        , &glu_zx       );
      t->Branch("glu_radL"      , &glu_radL     );
      t->Branch("glu_bbxi"      , &glu_bbxi     );
    }
    //phase2 OT
    if(includePhase2OTHits_){
      t->Branch("ph2_isBarrel"  , &ph2_isBarrel );
      t->Branch("ph2_det"       , &ph2_det      );
      t->Branch("ph2_lay"       , &ph2_lay      );
      t->Branch("ph2_detId"     , &ph2_detId    );
      t->Branch("ph2_trkIdx"    , &ph2_trkIdx   );
      if(includeSeeds_) {
        t->Branch("ph2_seeIdx"    , &ph2_seeIdx   );
      }
      t->Branch("ph2_simHitIdx" , &ph2_simHitIdx);
      t->Branch("ph2_simType", &ph2_simType);
      t->Branch("ph2_x"     , &ph2_x    );
      t->Branch("ph2_y"     , &ph2_y    );
      t->Branch("ph2_z"     , &ph2_z    );
      t->Branch("ph2_xx"    , &ph2_xx   );
      t->Branch("ph2_xy"    , &ph2_xy   );
      t->Branch("ph2_yy"    , &ph2_yy   );
      t->Branch("ph2_yz"    , &ph2_yz   );
      t->Branch("ph2_zz"    , &ph2_zz   );
      t->Branch("ph2_zx"    , &ph2_zx   );
      t->Branch("ph2_radL"  , &ph2_radL );
      t->Branch("ph2_bbxi"  , &ph2_bbxi );
      t->Branch("ph2_bbxi"  , &ph2_bbxi );
    }
    //invalid hits
    t->Branch("inv_isBarrel"  , &inv_isBarrel );
    t->Branch("inv_det"       , &inv_det      );
    t->Branch("inv_lay"       , &inv_lay      );
    t->Branch("inv_detId"     , &inv_detId    );
    t->Branch("inv_type"      , &inv_type    );
    //simhits
    t->Branch("simhit_det"     ,  &simhit_det);
    t->Branch("simhit_lay"    ,  &simhit_lay);
    t->Branch("simhit_detId"   ,  &simhit_detId);
    t->Branch("simhit_x"       ,  &simhit_x);
    t->Branch("simhit_y"       ,  &simhit_y);
    t->Branch("simhit_z"       ,  &simhit_z);
    t->Branch("simhit_particle",  &simhit_particle);
    t->Branch("simhit_process" ,  &simhit_process);
    t->Branch("simhit_eloss"   ,  &simhit_eloss);
    t->Branch("simhit_tof"     ,  &simhit_tof);
    //t->Branch("simhit_simTrackId", &simhit_simTrackId);
    t->Branch("simhit_simTrkIdx", &simhit_simTrkIdx);
    t->Branch("simhit_hitIdx"  ,  &simhit_hitIdx);
    t->Branch("simhit_hitType" ,  &simhit_hitType);
  }
  //beam spot
  t->Branch("bsp_x" , &bsp_x , "bsp_x/F");
  t->Branch("bsp_y" , &bsp_y , "bsp_y/F");
  t->Branch("bsp_z" , &bsp_z , "bsp_z/F");
  t->Branch("bsp_sigmax" , &bsp_sigmax , "bsp_sigmax/F");
  t->Branch("bsp_sigmay" , &bsp_sigmay , "bsp_sigmay/F");
  t->Branch("bsp_sigmaz" , &bsp_sigmaz , "bsp_sigmaz/F");
  if(includeSeeds_) {
    //seeds
    t->Branch("see_fitok"    , &see_fitok   );
    t->Branch("see_px"       , &see_px      );
    t->Branch("see_py"       , &see_py      );
    t->Branch("see_pz"       , &see_pz      );
    t->Branch("see_pt"       , &see_pt      );
    t->Branch("see_eta"      , &see_eta     );
    t->Branch("see_phi"      , &see_phi     );
    t->Branch("see_dxy"      , &see_dxy     );
    t->Branch("see_dz"       , &see_dz      );
    t->Branch("see_ptErr"    , &see_ptErr   );
    t->Branch("see_etaErr"   , &see_etaErr  );
    t->Branch("see_phiErr"   , &see_phiErr  );
    t->Branch("see_dxyErr"   , &see_dxyErr  );
    t->Branch("see_dzErr"    , &see_dzErr   );
    t->Branch("see_chi2"     , &see_chi2    );
    t->Branch("see_q"        , &see_q       );
    t->Branch("see_nValid"   , &see_nValid  );
    t->Branch("see_nPixel"   , &see_nPixel  );
    t->Branch("see_nGlued"   , &see_nGlued  );
    t->Branch("see_nStrip"   , &see_nStrip  );
    t->Branch("see_nPhase2OT", &see_nPhase2OT);
    t->Branch("see_algo"     , &see_algo    );
    t->Branch("see_trkIdx"   , &see_trkIdx  );
    t->Branch("see_shareFrac", &see_shareFrac);
    t->Branch("see_simTrkIdx", &see_simTrkIdx  );
    if(includeAllHits_) {
      t->Branch("see_hitIdx" , &see_hitIdx  );
      t->Branch("see_hitType", &see_hitType );
    }
    //seed algo offset
    t->Branch("see_offset"  , &see_offset );
  }

  //vertices
  t->Branch("vtx_x"       , &vtx_x);
  t->Branch("vtx_y"       , &vtx_y);
  t->Branch("vtx_z"       , &vtx_z);
  t->Branch("vtx_xErr"    , &vtx_xErr);
  t->Branch("vtx_yErr"    , &vtx_yErr);
  t->Branch("vtx_zErr"    , &vtx_zErr);
  t->Branch("vtx_ndof"    , &vtx_ndof);
  t->Branch("vtx_chi2"    , &vtx_chi2);
  t->Branch("vtx_fake"    , &vtx_fake);
  t->Branch("vtx_valid"   , &vtx_valid);
  t->Branch("vtx_trkIdx"  , &vtx_trkIdx);

  // tracking vertices
  t->Branch("simvtx_event"   , &simvtx_event    );
  t->Branch("simvtx_bunchCrossing", &simvtx_bunchCrossing);
  t->Branch("simvtx_processType", &simvtx_processType);
  t->Branch("simvtx_x"       , &simvtx_x);
  t->Branch("simvtx_y"       , &simvtx_y);
  t->Branch("simvtx_z"       , &simvtx_z);
  t->Branch("simvtx_sourceSimIdx", &simvtx_sourceSimIdx);
  t->Branch("simvtx_daughterSimIdx", &simvtx_daughterSimIdx);

  t->Branch("simpv_idx"      , &simpv_idx);

  //t->Branch("" , &);
}


TrackingNtuple::~TrackingNtuple() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}


//
// member functions
//
void TrackingNtuple::clearVariables() {

  ev_run = 0;
  ev_lumi = 0;
  ev_event = 0;

  //tracks
  trk_px       .clear();
  trk_py       .clear();
  trk_pz       .clear();
  trk_pt       .clear();
  trk_inner_px .clear();
  trk_inner_py .clear();
  trk_inner_pz .clear();
  trk_inner_pt .clear();
  trk_outer_px .clear();
  trk_outer_py .clear();
  trk_outer_pz .clear();
  trk_outer_pt .clear();
  trk_eta      .clear();
  trk_lambda   .clear();
  trk_cotTheta .clear();
  trk_phi      .clear();
  trk_dxy      .clear();
  trk_dz       .clear();
  trk_ptErr    .clear();
  trk_etaErr   .clear();
  trk_lambdaErr.clear();
  trk_phiErr   .clear();
  trk_dxyErr   .clear();
  trk_dzErr    .clear();
  trk_refpoint_x.clear();
  trk_refpoint_y.clear();
  trk_refpoint_z.clear();
  trk_nChi2    .clear();
  trk_q        .clear();
  trk_nValid   .clear();
  trk_nInvalid .clear();
  trk_nPixel   .clear();
  trk_nStrip   .clear();
  trk_nPixelLay.clear();
  trk_nStripLay.clear();
  trk_n3DLay   .clear();
  trk_nOuterLost.clear();
  trk_nInnerLost.clear();
  trk_algo     .clear();
  trk_originalAlgo.clear();
  trk_algoMask .clear();
  trk_stopReason.clear();
  trk_isHP     .clear();
  trk_seedIdx  .clear();
  trk_vtxIdx   .clear();
  trk_shareFrac.clear();
  trk_simTrkIdx.clear();
  trk_hitIdx   .clear();
  trk_hitType  .clear();
  //sim tracks
  sim_event    .clear();
  sim_bunchCrossing.clear();
  sim_pdgId    .clear();
  sim_px       .clear();
  sim_py       .clear();
  sim_pz       .clear();
  sim_pt       .clear();
  sim_eta      .clear();
  sim_phi      .clear();
  sim_pca_pt   .clear();
  sim_pca_eta  .clear();
  sim_pca_lambda.clear();
  sim_pca_cotTheta.clear();
  sim_pca_phi  .clear();
  sim_pca_dxy  .clear();
  sim_pca_dz   .clear();
  sim_q        .clear();
  sim_nValid   .clear();
  sim_nPixel   .clear();
  sim_nStrip   .clear();
  sim_nLay     .clear();
  sim_nPixelLay.clear();
  sim_n3DLay   .clear();
  sim_trkIdx   .clear();
  sim_shareFrac.clear();
  sim_parentVtxIdx.clear();
  sim_decayVtxIdx.clear();
  sim_simHitIdx   .clear();
  //pixels
  pix_isBarrel .clear();
  pix_det      .clear();
  pix_lay      .clear();
  pix_detId    .clear();
  pix_trkIdx   .clear();
  pix_seeIdx   .clear();
  pix_simHitIdx.clear();
  pix_chargeFraction.clear();
  pix_simType.clear();
  pix_x    .clear();
  pix_y    .clear();
  pix_z    .clear();
  pix_xx   .clear();
  pix_xy   .clear();
  pix_yy   .clear();
  pix_yz   .clear();
  pix_zz   .clear();
  pix_zx   .clear();
  pix_radL .clear();
  pix_bbxi .clear();
  //strips
  str_isBarrel .clear();
  str_isStereo .clear();
  str_det      .clear();
  str_lay      .clear();
  str_detId    .clear();
  str_trkIdx   .clear();
  str_seeIdx   .clear();
  str_simHitIdx.clear();
  str_chargeFraction.clear();
  str_simType.clear();
  str_x    .clear();
  str_y    .clear();
  str_z    .clear();
  str_xx   .clear();
  str_xy   .clear();
  str_yy   .clear();
  str_yz   .clear();
  str_zz   .clear();
  str_zx   .clear();
  str_radL .clear();
  str_bbxi .clear();
  //matched hits
  glu_isBarrel .clear();
  glu_det      .clear();
  glu_lay      .clear();
  glu_detId    .clear();
  glu_monoIdx  .clear();
  glu_stereoIdx.clear();
  glu_seeIdx   .clear();
  glu_x        .clear();
  glu_y        .clear();
  glu_z        .clear();
  glu_xx       .clear();
  glu_xy       .clear();
  glu_yy       .clear();
  glu_yz       .clear();
  glu_zz       .clear();
  glu_zx       .clear();
  glu_radL     .clear();
  glu_bbxi     .clear();
  //phase2 OT
  ph2_isBarrel .clear();
  ph2_det      .clear();
  ph2_lay      .clear();
  ph2_detId    .clear();
  ph2_trkIdx   .clear();
  ph2_seeIdx   .clear();
  ph2_simHitIdx.clear();
  ph2_simType.clear();
  ph2_x    .clear();
  ph2_y    .clear();
  ph2_z    .clear();
  ph2_xx   .clear();
  ph2_xy   .clear();
  ph2_yy   .clear();
  ph2_yz   .clear();
  ph2_zz   .clear();
  ph2_zx   .clear();
  ph2_radL .clear();
  ph2_bbxi .clear();
  //invalid hits
  inv_isBarrel .clear();
  inv_det      .clear();
  inv_lay      .clear();
  inv_detId    .clear();
  inv_type     .clear();
  // simhits
  simhit_det.clear();
  simhit_lay.clear();
  simhit_detId.clear();
  simhit_x.clear();
  simhit_y.clear();
  simhit_z.clear();
  simhit_particle.clear();
  simhit_process.clear();
  simhit_eloss.clear();
  simhit_tof.clear();
  //simhit_simTrackId.clear();
  simhit_simTrkIdx.clear();
  simhit_hitIdx.clear();
  simhit_hitType.clear();
  //beamspot
  bsp_x = -9999.;
  bsp_y = -9999.;
  bsp_z = -9999.;
  bsp_sigmax = -9999.;
  bsp_sigmay = -9999.;
  bsp_sigmaz = -9999.;
  //seeds
  see_fitok   .clear();
  see_px      .clear();
  see_py      .clear();
  see_pz      .clear();
  see_pt      .clear();
  see_eta     .clear();
  see_phi     .clear();
  see_dxy     .clear();
  see_dz      .clear();
  see_ptErr   .clear();
  see_etaErr  .clear();
  see_phiErr  .clear();
  see_dxyErr  .clear();
  see_dzErr   .clear();
  see_chi2    .clear();
  see_q       .clear();
  see_nValid  .clear();
  see_nPixel  .clear();
  see_nGlued  .clear();
  see_nStrip  .clear();
  see_nPhase2OT.clear();
  see_algo    .clear();
  see_trkIdx  .clear();
  see_shareFrac.clear();
  see_simTrkIdx.clear();
  see_hitIdx  .clear();
  see_hitType .clear();
  //seed algo offset
  see_offset.clear();

  // vertices
  vtx_x.clear();
  vtx_y.clear();
  vtx_z.clear();
  vtx_xErr.clear();
  vtx_yErr.clear();
  vtx_zErr.clear();
  vtx_ndof.clear();
  vtx_chi2.clear();
  vtx_fake.clear();
  vtx_valid.clear();
  vtx_trkIdx.clear();

  // Tracking vertices
  simvtx_event.clear();
  simvtx_bunchCrossing.clear();
  simvtx_x.clear();
  simvtx_y.clear();
  simvtx_z.clear();
  simvtx_sourceSimIdx.clear();
  simvtx_daughterSimIdx.clear();
  simpv_idx.clear();
}


// ------------ method called for each event  ------------
void TrackingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {

  using namespace edm;
  using namespace reco;
  using namespace std;

  edm::ESHandle<MagneticField> theMF;
  iSetup.get<IdealMagneticFieldRecord>().get(theMF);

  edm::ESHandle<TransientTrackingRecHitBuilder> theTTRHBuilder;
  iSetup.get<TransientRecHitRecord>().get(builderName_,theTTRHBuilder);

  edm::ESHandle<TrackerTopology> tTopoHandle;
  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology& tTopo = *tTopoHandle;

  edm::ESHandle<TrackerGeometry> geometryHandle;
  iSetup.get<TrackerDigiGeometryRecord>().get(geometryHandle);
  const TrackerGeometry &tracker = *geometryHandle;

  edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
  iEvent.getByToken(trackAssociatorToken_, theAssociator);
  const reco::TrackToTrackingParticleAssociator& associatorByHits = *theAssociator;

  LogDebug("TrackingNtuple") << "Analyzing new event";

  //initialize tree variables
  clearVariables();

  // FIXME: we really need to move to edm::View for reading the
  // TrackingParticles... Unfortunately it has non-trivial
  // consequences on the associator/association interfaces etc.
  TrackingParticleRefVector tmpTP;
  TrackingParticleRefKeySet tmpTPkeys;
  const TrackingParticleRefVector *tmpTPptr = nullptr;
  edm::Handle<TrackingParticleCollection>  TPCollectionH;
  edm::Handle<TrackingParticleRefVector>  TPCollectionHRefVector;

  if(!trackingParticleToken_.isUninitialized()) {
    iEvent.getByToken(trackingParticleToken_, TPCollectionH);
    for(size_t i=0, size=TPCollectionH->size(); i<size; ++i) {
      tmpTP.push_back(TrackingParticleRef(TPCollectionH, i));
    }
    tmpTPptr = &tmpTP;
  }
  else {
    iEvent.getByToken(trackingParticleRefToken_, TPCollectionHRefVector);
    tmpTPptr = TPCollectionHRefVector.product();
    for(const auto& ref: *tmpTPptr) {
      tmpTPkeys.insert(ref.key());
    }
  }
  const TrackingParticleRefVector& tpCollection = *tmpTPptr;

  // Fill mapping from Ref::key() to index
  TrackingParticleRefKeyToIndex tpKeyToIndex;
  for(size_t i=0; i<tpCollection.size(); ++i) {
    tpKeyToIndex[tpCollection[i].key()] = i;
  }

  // tracking vertices
  edm::Handle<TrackingVertexCollection> htv;
  iEvent.getByToken(trackingVertexToken_, htv);
  const TrackingVertexCollection& tvs = *htv;

  // Fill mapping from Ref::key() to index
  TrackingVertexRefVector tvRefs;
  TrackingVertexRefKeyToIndex tvKeyToIndex;
  for(size_t i=0; i<tvs.size(); ++i) {
    const TrackingVertex& v = tvs[i];
    if(v.eventId().bunchCrossing() != 0) // Ignore OOTPU; would be better to not to hardcode?
      continue;
    tvKeyToIndex[i] = tvRefs.size();
    tvRefs.push_back(TrackingVertexRef(htv, i));
  }

  //get association maps, etc.
  Handle<ClusterTPAssociation> pCluster2TPListH;
  iEvent.getByToken(clusterTPMapToken_, pCluster2TPListH);
  const ClusterTPAssociation& clusterToTPMap = *pCluster2TPListH;
  edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
  iEvent.getByToken(simHitTPMapToken_, simHitsTPAssoc);

  // SimHit key -> index mapping
  SimHitRefKeyToIndex simHitRefKeyToIndex;

  //make a list to link TrackingParticles to its simhits
  std::vector<TPHitIndex> tpHitList;

  std::set<edm::ProductID> hitProductIds;
  std::map<edm::ProductID, size_t> seedCollToOffset;

  ev_run = iEvent.id().run();
  ev_lumi = iEvent.id().luminosityBlock();
  ev_event = iEvent.id().event();

  // Digi->Sim links for pixels and strips
  edm::Handle<edm::DetSetVector<PixelDigiSimLink> > pixelDigiSimLinksHandle;
  iEvent.getByToken(pixelSimLinkToken_, pixelDigiSimLinksHandle);
  const auto& pixelDigiSimLinks = *pixelDigiSimLinksHandle;

  edm::Handle<edm::DetSetVector<StripDigiSimLink> > stripDigiSimLinksHandle;
  iEvent.getByToken(stripSimLinkToken_, stripDigiSimLinksHandle);

  // Phase2 OT DigiSimLink
  edm::Handle<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinksHandle;
  iEvent.getByToken(siphase2OTSimLinksToken_, siphase2OTSimLinksHandle);

  //beamspot
  Handle<reco::BeamSpot> recoBeamSpotHandle;
  iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
  BeamSpot const & bs = *recoBeamSpotHandle;
  fillBeamSpot(bs);


  //prapare list to link matched hits to collection
  vector<pair<int,int> > monoStereoClusterList;
  if(includeAllHits_) {
    // simhits
    fillSimHits(tracker, tpKeyToIndex, *simHitsTPAssoc, tTopo, simHitRefKeyToIndex, tpHitList);

    //pixel hits
    fillPixelHits(iEvent, clusterToTPMap, tpKeyToIndex, *simHitsTPAssoc, pixelDigiSimLinks, *theTTRHBuilder, tTopo, simHitRefKeyToIndex, hitProductIds);

    //strip hits
    if(includeStripHits_){
      LogDebug("TrackingNtuple") << "foundStripSimLink" ;
      const auto& stripDigiSimLinks = *stripDigiSimLinksHandle;
      fillStripRphiStereoHits(iEvent, clusterToTPMap, tpKeyToIndex, *simHitsTPAssoc, stripDigiSimLinks, *theTTRHBuilder, tTopo, simHitRefKeyToIndex, hitProductIds);

      //matched hits
      fillStripMatchedHits(iEvent, *theTTRHBuilder, tTopo, monoStereoClusterList);
    }

    if(includePhase2OTHits_){
      LogDebug("TrackingNtuple") << "foundPhase2OTSimLinks" ;
      const auto& phase2OTSimLinks = *siphase2OTSimLinksHandle;
      fillPhase2OTHits(iEvent, clusterToTPMap, tpKeyToIndex, *simHitsTPAssoc, phase2OTSimLinks, *theTTRHBuilder, tTopo, simHitRefKeyToIndex, hitProductIds);
    }
  }

  //seeds
  if(includeSeeds_) {
    fillSeeds(iEvent, tpCollection, tpKeyToIndex, bs, associatorByHits, *theTTRHBuilder, theMF.product(), monoStereoClusterList, hitProductIds, seedCollToOffset);
  }

  //tracks
  edm::Handle<edm::View<reco::Track> > tracksHandle;
  iEvent.getByToken(trackToken_, tracksHandle);
  const edm::View<reco::Track>& tracks = *tracksHandle;
  // The associator interfaces really need to be fixed...
  edm::RefToBaseVector<reco::Track> trackRefs;
  for(edm::View<Track>::size_type i=0; i<tracks.size(); ++i) {
    trackRefs.push_back(tracks.refAt(i));
  }
  fillTracks(trackRefs, tpCollection, tpKeyToIndex, bs, associatorByHits, *theTTRHBuilder, tTopo, hitProductIds, seedCollToOffset);

  //tracking particles
  //sort association maps with simHits
  std::sort( tpHitList.begin(), tpHitList.end(), tpHitIndexListLessSort );
  fillTrackingParticles(iEvent, iSetup, trackRefs, tpCollection, tvKeyToIndex, associatorByHits, tpHitList);

  // vertices
  edm::Handle<reco::VertexCollection> vertices;
  iEvent.getByToken(vertexToken_, vertices);
  fillVertices(*vertices, trackRefs);

  // tracking vertices
  fillTrackingVertices(tvRefs, tpKeyToIndex);

  t->Fill();

}

void TrackingNtuple::fillBeamSpot(const reco::BeamSpot& bs) {
  bsp_x = bs.x0();
  bsp_y = bs.y0();
  bsp_z = bs.x0();
  bsp_sigmax = bs.BeamWidthX();
  bsp_sigmay = bs.BeamWidthY();
  bsp_sigmaz = bs.sigmaZ();
}

namespace {
  template <typename SimLink> struct GetCluster;
  template <>
  struct GetCluster<PixelDigiSimLink> {
    static const SiPixelCluster& call(const OmniClusterRef& cluster) { return cluster.pixelCluster(); }
  };
  template <>
  struct GetCluster<StripDigiSimLink> {
    static const SiStripCluster& call(const OmniClusterRef& cluster) { return cluster.stripCluster(); }
  };
}

template <typename SimLink>
TrackingNtuple::SimHitData TrackingNtuple::matchCluster(const OmniClusterRef& cluster,
                                                        DetId hitId, int clusterKey,
                                                        const TransientTrackingRecHit::RecHitPointer& ttrh,
                                                        const ClusterTPAssociation& clusterToTPMap,
                                                        const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                                        const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                                                        const edm::DetSetVector<SimLink>& digiSimLinks,
                                                        const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                                                        HitType hitType
                                                        ) {
  SimHitData ret;

  std::map<unsigned int, double> simTrackIdToChargeFraction;
  if(hitType == HitType::Phase2OT) simTrackIdToChargeFraction = chargeFraction(cluster.phase2OTCluster(), hitId, digiSimLinks);
  else simTrackIdToChargeFraction = chargeFraction(GetCluster<SimLink>::call(cluster), hitId, digiSimLinks);

  ret.type = HitSimType::Noise;
  auto range = clusterToTPMap.equal_range( cluster );
  if( range.first != range.second ) {
    for( auto ip=range.first; ip != range.second; ++ip ) {
      const TrackingParticleRef& trackingParticle = ip->second;

      // Find out if the cluster is from signal/ITPU/OOTPU
      const auto event = trackingParticle->eventId().event();
      const auto bx = trackingParticle->eventId().bunchCrossing();
      HitSimType type = HitSimType::OOTPileup;
      if(bx == 0) {
        type = (event == 0 ? HitSimType::Signal : HitSimType::ITPileup);
      }
      ret.type = static_cast<HitSimType>(std::min(static_cast<int>(ret.type), static_cast<int>(type)));

      // Limit to only input TrackingParticles (usually signal+ITPU)
      auto tpIndex = tpKeyToIndex.find(trackingParticle.key());
      if( tpIndex == tpKeyToIndex.end())
        continue;

      //now get the corresponding sim hit
      std::pair<TrackingParticleRef, TrackPSimHitRef> simHitTPpairWithDummyTP(trackingParticle,TrackPSimHitRef());
      //SimHit is dummy: for simHitTPAssociationListGreater sorting only the TP is needed
      auto range = std::equal_range(simHitsTPAssoc.begin(), simHitsTPAssoc.end(),
                                    simHitTPpairWithDummyTP, SimHitTPAssociationProducer::simHitTPAssociationListGreater);
      int simHitKey = -1;
      bool foundElectron = false;
      edm::ProductID simHitID;
      for(auto ip = range.first; ip != range.second; ++ip) {
        TrackPSimHitRef TPhit = ip->second;
        DetId dId = DetId(TPhit->detUnitId());
        if (dId.rawId()==hitId.rawId()) {
          // skip electron SimHits for non-electron TPs also here
          if(std::abs(TPhit->particleType()) == 11 && std::abs(trackingParticle->pdgId()) != 11) {
            foundElectron = true;
            continue;
          }

          simHitKey = TPhit.key();
          simHitID = TPhit.id();
          break;
        }
      }
      if(simHitKey < 0) {
        // In case we didn't find a simhit because of filtered-out
        // electron SimHit, just ignore the missing SimHit.
        if(foundElectron)
          continue;

        auto ex = cms::Exception("LogicError") << "Did not find SimHit for reco hit DetId " << hitId.rawId()
                                               << " for TP " << trackingParticle.key() << " bx:event " << bx << ":" << event
                                               << ".\nFound SimHits from detectors ";
        for(auto ip = range.first; ip != range.second; ++ip) {
          TrackPSimHitRef TPhit = ip->second;
          DetId dId = DetId(TPhit->detUnitId());
          ex << dId.rawId() << " ";
        }
        if(trackingParticle->eventId().event() != 0) {
          ex << "\nSince this is a TrackingParticle from pileup, check that you're running the pileup mixing in playback mode.";
        }
        throw ex;
      }
      auto simHitIndex = simHitRefKeyToIndex.at(std::make_pair(simHitKey, simHitID));
      ret.matchingSimHit.push_back(simHitIndex);

      double chargeFraction = 0.;
      for(const SimTrack& simtrk: trackingParticle->g4Tracks()) {
        auto found = simTrackIdToChargeFraction.find(simtrk.trackId());
        if(found != simTrackIdToChargeFraction.end()) {
          chargeFraction += found->second;
        }
      }
      ret.chargeFraction.push_back(chargeFraction);

      // only for debug prints
      ret.bunchCrossing.push_back(bx);
      ret.event.push_back(event);

      simhit_hitIdx[simHitIndex].push_back(clusterKey);
      simhit_hitType[simHitIndex].push_back(static_cast<int>(hitType));
    }
  }

  return ret;
}

void TrackingNtuple::fillSimHits(const TrackerGeometry& tracker,
                                 const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                 const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                                 const TrackerTopology& tTopo,
                                 SimHitRefKeyToIndex& simHitRefKeyToIndex,
                                 std::vector<TPHitIndex>& tpHitList) {

  for(const auto& assoc: simHitsTPAssoc) {
    auto tpKey = assoc.first.key();

    // SimHitTPAssociationList can contain more TrackingParticles than
    // what are given to this EDAnalyzer, so we can filter those out here.
    auto found = tpKeyToIndex.find(tpKey);
    if(found == tpKeyToIndex.end())
      continue;
    const auto tpIndex = found->second;

    // skip non-tracker simhits (mostly muons)
    const auto& simhit = *(assoc.second);
    auto detId = DetId(simhit.detUnitId());
    if(detId.det() != DetId::Tracker) continue;

    // Skip electron SimHits for non-electron TrackingParticles to
    // filter out delta rays. The delta ray hits just confuse. If we
    // need them later, let's add them as a separate "collection" of
    // hits of a TP
    const TrackingParticle& tp = *(assoc.first);
    if(std::abs(simhit.particleType()) == 11 && std::abs(tp.pdgId()) != 11) continue;

    auto simHitKey = std::make_pair(assoc.second.key(), assoc.second.id());

    if(simHitRefKeyToIndex.find(simHitKey) != simHitRefKeyToIndex.end()) {
      for(const auto& assoc2: simHitsTPAssoc) {
        if(std::make_pair(assoc2.second.key(), assoc2.second.id()) == simHitKey) {

#ifdef EDM_ML_DEBUG
          auto range1 = std::equal_range(simHitsTPAssoc.begin(), simHitsTPAssoc.end(),
                                         std::make_pair(assoc.first, TrackPSimHitRef()),
                                         SimHitTPAssociationProducer::simHitTPAssociationListGreater);
          auto range2 = std::equal_range(simHitsTPAssoc.begin(), simHitsTPAssoc.end(),
                                         std::make_pair(assoc2.first, TrackPSimHitRef()),
                                         SimHitTPAssociationProducer::simHitTPAssociationListGreater);

          LogTrace("TrackingNtuple") << "Earlier TP " << assoc2.first.key() << " SimTrack Ids";
          for(const auto& simTrack: assoc2.first->g4Tracks()) {
            edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event " << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
          }
          for(auto iHit = range2.first; iHit != range2.second; ++iHit) {
            LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof " << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event " << iHit->second->eventId().bunchCrossing() << ":" << iHit->second->eventId().event();
          }
          LogTrace("TrackingNtuple") << "Current TP " << assoc.first.key() << " SimTrack Ids";
          for(const auto& simTrack: assoc.first->g4Tracks()) {
            edm::LogPrint("TrackingNtuple") << " SimTrack " << simTrack.trackId() << " BX:event " << simTrack.eventId().bunchCrossing() << ":" << simTrack.eventId().event();
          }
          for(auto iHit = range1.first; iHit != range1.second; ++iHit) {
            LogTrace("TrackingNtuple") << " SimHit " << iHit->second.key() << " " << iHit->second.id() << " tof " << iHit->second->tof() << " trackId " << iHit->second->trackId() << " BX:event " << iHit->second->eventId().bunchCrossing() << ":" << iHit->second->eventId().event();
          }
#endif

          throw cms::Exception("LogicError") << "Got second time the SimHit " << simHitKey.first << " of " << simHitKey.second << ", first time with TrackingParticle " << assoc2.first.key() << ", now with " << tpKey;
        }
      }
      throw cms::Exception("LogicError") << "Got second time the SimHit " << simHitKey.first << " of " << simHitKey.second << ", now with TrackingParticle " << tpKey << ", but I didn't find the first occurrance!";
    }

    auto det = tracker.idToDetUnit(detId);
    if(!det)
      throw cms::Exception("LogicError") << "Did not find a det unit for DetId " << simhit.detUnitId() << " from tracker geometry";

    const auto pos = det->surface().toGlobal(simhit.localPosition());
    const float tof = simhit.timeOfFlight();

    const auto simHitIndex = simhit_detId.size();
    simHitRefKeyToIndex[simHitKey] = simHitIndex;

    simhit_det.push_back(detId.subdetId());
    simhit_lay.push_back(tTopo.layer(detId));
    simhit_detId.push_back(detId.rawId());
    simhit_x.push_back(pos.x());
    simhit_y.push_back(pos.y());
    simhit_z.push_back(pos.z());
    simhit_particle.push_back(simhit.particleType());
    simhit_process.push_back(simhit.processType());
    simhit_eloss.push_back(simhit.energyLoss());
    simhit_tof.push_back(tof);
    //simhit_simTrackId.push_back(simhit.trackId());

    simhit_simTrkIdx.push_back(tpIndex);

    simhit_hitIdx.emplace_back();  // filled in matchCluster
    simhit_hitType.emplace_back(); // filled in matchCluster

    tpHitList.emplace_back(tpKey, simHitIndex, tof, simhit.detUnitId());
  }
}

void TrackingNtuple::fillPixelHits(const edm::Event& iEvent,
                                   const ClusterTPAssociation& clusterToTPMap,
                                   const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                   const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                                   const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
                                   const TransientTrackingRecHitBuilder& theTTRHBuilder,
                                   const TrackerTopology& tTopo,
                                   const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                                   std::set<edm::ProductID>& hitProductIds
                                   ) {
  edm::Handle<SiPixelRecHitCollection> pixelHits;
  iEvent.getByToken(pixelRecHitToken_, pixelHits);
  for (auto it = pixelHits->begin(); it!=pixelHits->end(); it++ ) {
    const DetId hitId = it->detId();
    for (auto hit = it->begin(); hit!=it->end(); hit++ ) {
      TransientTrackingRecHit::RecHitPointer ttrh = theTTRHBuilder.build(&*hit);

      hitProductIds.insert(hit->cluster().id());

      const int key = hit->cluster().key();
      const int lay = tTopo.layer(hitId);
      SimHitData simHitData = matchCluster(hit->firstClusterRef(), hitId, key, ttrh,
                                           clusterToTPMap, tpKeyToIndex, simHitsTPAssoc, digiSimLink, simHitRefKeyToIndex, HitType::Pixel);

      pix_isBarrel .push_back( hitId.subdetId()==1 );
      pix_det      .push_back( hitId.subdetId() );
      pix_lay      .push_back( lay );
      pix_detId    .push_back( hitId.rawId() );
      pix_trkIdx   .emplace_back(); // filled in fillTracks
      pix_seeIdx   .emplace_back(); // filled in fillSeeds
      pix_simHitIdx.push_back( simHitData.matchingSimHit );
      pix_simType.push_back( static_cast<int>(simHitData.type) );
      pix_x    .push_back( ttrh->globalPosition().x() );
      pix_y    .push_back( ttrh->globalPosition().y() );
      pix_z    .push_back( ttrh->globalPosition().z() );
      pix_xx   .push_back( ttrh->globalPositionError().cxx() );
      pix_xy   .push_back( ttrh->globalPositionError().cyx() );
      pix_yy   .push_back( ttrh->globalPositionError().cyy() );
      pix_yz   .push_back( ttrh->globalPositionError().czy() );
      pix_zz   .push_back( ttrh->globalPositionError().czz() );
      pix_zx   .push_back( ttrh->globalPositionError().czx() );
      pix_chargeFraction.push_back( simHitData.chargeFraction );
      pix_radL .push_back( ttrh->surface()->mediumProperties().radLen() );
      pix_bbxi .push_back( ttrh->surface()->mediumProperties().xi() );
      LogTrace("TrackingNtuple") << "pixHit cluster=" << key
                                 << " subdId=" << hitId.subdetId()
                                 << " lay=" << lay
                                 << " rawId=" << hitId.rawId()
                                 << " pos =" << ttrh->globalPosition()
                                 << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
      if(!simHitData.matchingSimHit.empty()) {
        const auto simHitIdx = simHitData.matchingSimHit[0];
        LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx
                                   << " simHitPos=" << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
                                   << " energyLoss=" << simhit_eloss[simHitIdx]
                                   << " particleType=" << simhit_particle[simHitIdx]
                                   << " processType=" << simhit_process[simHitIdx]
                                   << " bunchCrossing=" << simHitData.bunchCrossing[0]
                                   << " event=" << simHitData.event[0];
      }
    }
  }
}


void TrackingNtuple::fillStripRphiStereoHits(const edm::Event& iEvent,
                                             const ClusterTPAssociation& clusterToTPMap,
                                             const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                             const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                                             const edm::DetSetVector<StripDigiSimLink>& digiSimLink,
                                             const TransientTrackingRecHitBuilder& theTTRHBuilder,
                                             const TrackerTopology& tTopo,
                                             const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                                             std::set<edm::ProductID>& hitProductIds
                                             ) {
  //index strip hit branches by cluster index
  edm::Handle<SiStripRecHit2DCollection> rphiHits;
  iEvent.getByToken(stripRphiRecHitToken_, rphiHits);
  edm::Handle<SiStripRecHit2DCollection> stereoHits;
  iEvent.getByToken(stripStereoRecHitToken_, stereoHits);
  int totalStripHits = rphiHits->dataSize()+stereoHits->dataSize();
  str_isBarrel .resize(totalStripHits);
  str_isStereo .resize(totalStripHits);
  str_det      .resize(totalStripHits);
  str_lay      .resize(totalStripHits);
  str_detId    .resize(totalStripHits);
  str_trkIdx   .resize(totalStripHits); // filled in fillTracks
  str_seeIdx   .resize(totalStripHits); // filled in fillSeeds
  str_simHitIdx.resize(totalStripHits);
  str_simType  .resize(totalStripHits);
  str_x    .resize(totalStripHits);
  str_y    .resize(totalStripHits);
  str_z    .resize(totalStripHits);
  str_xx   .resize(totalStripHits);
  str_xy   .resize(totalStripHits);
  str_yy   .resize(totalStripHits);
  str_yz   .resize(totalStripHits);
  str_zz   .resize(totalStripHits);
  str_zx   .resize(totalStripHits);
  str_chargeFraction.resize(totalStripHits);
  str_radL .resize(totalStripHits);
  str_bbxi .resize(totalStripHits);

  auto fill = [&](const SiStripRecHit2DCollection& hits, const char *name, bool isStereo) {
    for(const auto& detset: hits) {
      const DetId hitId = detset.detId();
      for(const auto& hit: detset) {
        TransientTrackingRecHit::RecHitPointer ttrh = theTTRHBuilder.build(&hit);

        hitProductIds.insert(hit.cluster().id());

        const int key = hit.cluster().key();
        const int lay = tTopo.layer(hitId);
        SimHitData simHitData = matchCluster(hit.firstClusterRef(), hitId, key, ttrh,
                                             clusterToTPMap, tpKeyToIndex, simHitsTPAssoc, digiSimLink, simHitRefKeyToIndex, HitType::Strip);
        str_isBarrel [key] = (hitId.subdetId()==StripSubdetector::TIB || hitId.subdetId()==StripSubdetector::TOB);
        str_isStereo [key] = isStereo;
        str_det      [key] = hitId.subdetId();
        str_lay      [key] = lay;
        str_detId    [key] = hitId.rawId();
        str_simHitIdx[key] = simHitData.matchingSimHit;
        str_simType  [key] = static_cast<int>(simHitData.type);
        str_x    [key] = ttrh->globalPosition().x();
        str_y    [key] = ttrh->globalPosition().y();
        str_z    [key] = ttrh->globalPosition().z();
        str_xx   [key] = ttrh->globalPositionError().cxx();
        str_xy   [key] = ttrh->globalPositionError().cyx();
        str_yy   [key] = ttrh->globalPositionError().cyy();
        str_yz   [key] = ttrh->globalPositionError().czy();
        str_zz   [key] = ttrh->globalPositionError().czz();
        str_zx   [key] = ttrh->globalPositionError().czx();
        str_chargeFraction[key] = simHitData.chargeFraction;
        str_radL [key] = ttrh->surface()->mediumProperties().radLen();
        str_bbxi [key] = ttrh->surface()->mediumProperties().xi();
        LogTrace("TrackingNtuple") << name << " cluster=" << key
                                   << " subdId=" << hitId.subdetId()
                                   << " lay=" << lay
                                   << " rawId=" << hitId.rawId()
                                   << " pos =" << ttrh->globalPosition()
                                   << " nMatchingSimHit=" << simHitData.matchingSimHit.size();
        if(!simHitData.matchingSimHit.empty()) {
          const auto simHitIdx = simHitData.matchingSimHit[0];
          LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx
                                     << " simHitPos=" << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
                                     << " simHitPos=" << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
                                     << " energyLoss=" << simhit_eloss[simHitIdx]
                                     << " particleType=" << simhit_particle[simHitIdx]
                                     << " processType=" << simhit_process[simHitIdx]
                                     << " bunchCrossing=" << simHitData.bunchCrossing[0]
                                     << " event=" << simHitData.event[0];
        }
      }
    }
  };

  fill(*rphiHits, "stripRPhiHit", false);
  fill(*stereoHits, "stripStereoHit", true);
}

void TrackingNtuple::fillStripMatchedHits(const edm::Event& iEvent,
                                          const TransientTrackingRecHitBuilder& theTTRHBuilder,
                                          const TrackerTopology& tTopo,
                                          std::vector<std::pair<int, int> >& monoStereoClusterList
                                          ) {
  edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
  iEvent.getByToken(stripMatchedRecHitToken_, matchedHits);
  for (auto it = matchedHits->begin(); it!=matchedHits->end(); it++ ) {
    const DetId hitId = it->detId();
    for (auto hit = it->begin(); hit!=it->end(); hit++ ) {
      TransientTrackingRecHit::RecHitPointer ttrh = theTTRHBuilder.build(&*hit);
      const int lay = tTopo.layer(hitId);
      monoStereoClusterList.emplace_back(hit->monoHit().cluster().key(),hit->stereoHit().cluster().key());
      glu_isBarrel .push_back( (hitId.subdetId()==StripSubdetector::TIB || hitId.subdetId()==StripSubdetector::TOB) );
      glu_det      .push_back( hitId.subdetId() );
      glu_lay      .push_back( tTopo.layer(hitId) );
      glu_detId    .push_back( hitId.rawId() );
      glu_monoIdx  .push_back( hit->monoHit().cluster().key() );
      glu_stereoIdx.push_back( hit->stereoHit().cluster().key() );
      glu_seeIdx   .emplace_back(); // filled in fillSeeds
      glu_x        .push_back( ttrh->globalPosition().x() );
      glu_y        .push_back( ttrh->globalPosition().y() );
      glu_z        .push_back( ttrh->globalPosition().z() );
      glu_xx       .push_back( ttrh->globalPositionError().cxx() );
      glu_xy       .push_back( ttrh->globalPositionError().cyx() );
      glu_yy       .push_back( ttrh->globalPositionError().cyy() );
      glu_yz       .push_back( ttrh->globalPositionError().czy() );
      glu_zz       .push_back( ttrh->globalPositionError().czz() );
      glu_zx       .push_back( ttrh->globalPositionError().czx() );
      glu_radL     .push_back( ttrh->surface()->mediumProperties().radLen() );
      glu_bbxi     .push_back( ttrh->surface()->mediumProperties().xi() );
      LogTrace("TrackingNtuple") << "stripMatchedHit"
                                 << " cluster0=" << hit->stereoHit().cluster().key()
                                 << " cluster1=" << hit->monoHit().cluster().key()
                                 << " subdId=" << hitId.subdetId()
                                 << " lay=" << lay
                                 << " rawId=" << hitId.rawId()
                                 << " pos =" << ttrh->globalPosition();
    }
  }
}

void TrackingNtuple::fillPhase2OTHits(const edm::Event& iEvent,
                                   const ClusterTPAssociation& clusterToTPMap,
                                   const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                   const SimHitTPAssociationProducer::SimHitTPAssociationList& simHitsTPAssoc,
                                   const edm::DetSetVector<PixelDigiSimLink>& digiSimLink,
                                   const TransientTrackingRecHitBuilder& theTTRHBuilder,
                                   const TrackerTopology& tTopo,
                                   const SimHitRefKeyToIndex& simHitRefKeyToIndex,
                                   std::set<edm::ProductID>& hitProductIds
                                   ) {
  edm::Handle<Phase2TrackerRecHit1DCollectionNew> phase2OTHits;
  iEvent.getByToken(phase2OTRecHitToken_, phase2OTHits);
  for (auto it = phase2OTHits->begin(); it!=phase2OTHits->end(); it++ ) {
    const DetId hitId = it->detId();
    for (auto hit = it->begin(); hit!=it->end(); hit++ ) {
      TransientTrackingRecHit::RecHitPointer ttrh = theTTRHBuilder.build(&*hit);

      hitProductIds.insert(hit->cluster().id());

      const int key = hit->cluster().key();
      const int lay = tTopo.layer(hitId);
      SimHitData simHitData = matchCluster(hit->firstClusterRef(), hitId, key, ttrh,
                                           clusterToTPMap, tpKeyToIndex, simHitsTPAssoc, digiSimLink, simHitRefKeyToIndex, HitType::Phase2OT);

      ph2_isBarrel .push_back( hitId.subdetId()==1 );
      ph2_det      .push_back( hitId.subdetId() );
      ph2_lay      .push_back( lay );
      ph2_detId    .push_back( hitId.rawId() );
      ph2_trkIdx   .emplace_back(); // filled in fillTracks
      ph2_seeIdx   .emplace_back(); // filled in fillSeeds
      ph2_simHitIdx.push_back( simHitData.matchingSimHit );
      ph2_simType.push_back( static_cast<int>(simHitData.type) );
      ph2_x    .push_back( ttrh->globalPosition().x() );
      ph2_y    .push_back( ttrh->globalPosition().y() );
      ph2_z    .push_back( ttrh->globalPosition().z() );
      ph2_xx   .push_back( ttrh->globalPositionError().cxx() );
      ph2_xy   .push_back( ttrh->globalPositionError().cyx() );
      ph2_yy   .push_back( ttrh->globalPositionError().cyy() );
      ph2_yz   .push_back( ttrh->globalPositionError().czy() );
      ph2_zz   .push_back( ttrh->globalPositionError().czz() );
      ph2_zx   .push_back( ttrh->globalPositionError().czx() );
      ph2_radL .push_back( ttrh->surface()->mediumProperties().radLen() );
      ph2_bbxi .push_back( ttrh->surface()->mediumProperties().xi() );

      LogTrace("TrackingNtuple") << "phase2 OT cluster=" << key
                                 << " subdId=" << hitId.subdetId()
                                 << " lay=" << lay
                                 << " rawId=" << hitId.rawId()
                                 << " pos =" << ttrh->globalPosition()
                                 << " nMatchingSimHit=" << simHitData.matchingSimHit.size();

      if(!simHitData.matchingSimHit.empty()) {
        const auto simHitIdx = simHitData.matchingSimHit[0];
        LogTrace("TrackingNtuple") << " firstMatchingSimHit=" << simHitIdx
                                   << " simHitPos=" << GlobalPoint(simhit_x[simHitIdx], simhit_y[simHitIdx], simhit_z[simHitIdx])
                                   << " energyLoss=" << simhit_eloss[simHitIdx]
                                   << " particleType=" << simhit_particle[simHitIdx]
                                   << " processType=" << simhit_process[simHitIdx]
                                   << " bunchCrossing=" << simHitData.bunchCrossing[0]
                                   << " event=" << simHitData.event[0];
      }
    }
  }
}


void TrackingNtuple::fillSeeds(const edm::Event& iEvent,
                               const TrackingParticleRefVector& tpCollection,
                               const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                               const reco::BeamSpot& bs,
                               const reco::TrackToTrackingParticleAssociator& associatorByHits,
                               const TransientTrackingRecHitBuilder& theTTRHBuilder,
                               const MagneticField *theMF,
                               const std::vector<std::pair<int, int> >& monoStereoClusterList,
                               const std::set<edm::ProductID>& hitProductIds,
                               std::map<edm::ProductID, size_t>& seedCollToOffset
                               ) {
  TSCBLBuilderNoMaterial tscblBuilder;
  for(const auto& seedToken: seedTokens_) {
    edm::Handle<edm::View<reco::Track> > seedTracksHandle;
    iEvent.getByToken(seedToken, seedTracksHandle);
    const auto& seedTracks = *seedTracksHandle;

    if(seedTracks.empty())
      continue;

    // The associator interfaces really need to be fixed...
    edm::RefToBaseVector<reco::Track> seedTrackRefs;
    for(edm::View<reco::Track>::size_type i=0; i<seedTracks.size(); ++i) {
      seedTrackRefs.push_back(seedTracks.refAt(i));
    }
    reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(seedTrackRefs, tpCollection);

    edm::EDConsumerBase::Labels labels;
    labelsForToken(seedToken, labels);
    TString label = labels.module;
    //format label to match algoName
    label.ReplaceAll("seedTracks", "");
    label.ReplaceAll("Seeds","");
    label.ReplaceAll("muonSeeded","muonSeededStep");
    int algo = reco::TrackBase::algoByName(label.Data());

    edm::ProductID id = seedTracks[0].seedRef().id();
    auto inserted = seedCollToOffset.emplace(id, see_fitok.size());
    if(!inserted.second)
      throw cms::Exception("Configuration") << "Trying to add seeds with ProductID " << id << " for a second time from collection " << labels.module << ", seed algo " << label << ". Typically this is caused by a configuration problem.";
    see_offset.push_back(see_fitok.size());

    LogTrace("TrackingNtuple") << "NEW SEED LABEL: " << label << " size: " << seedTracks.size() << " algo=" << algo
                               << " ProductID " << id;

    for(const auto& seedTrackRef: seedTrackRefs) {

      const auto& seedTrack = *seedTrackRef;
      const auto& seedRef = seedTrack.seedRef();
      const auto& seed = *seedRef;

      if(seedRef.id() != id)
        throw cms::Exception("LogicError") << "All tracks in 'TracksFromSeeds' collection should point to seeds in the same collection. Now the element 0 had ProductID " << id << " while the element " << seedTrackRef.key() << " had " << seedTrackRef.id() << ". The source collection is " << labels.module << ".";

      std::vector<float> sharedFraction;
      std::vector<int> tpIdx;
      auto foundTPs = recSimColl.find(seedTrackRef);
      if (foundTPs != recSimColl.end()) {
        for(const auto tpQuality: foundTPs->val) {
          sharedFraction.push_back(tpQuality.second);
          tpIdx.push_back( tpKeyToIndex.at( tpQuality.first.key() ) );
        }
      }


      const bool seedFitOk = !trackFromSeedFitFailed(seedTrack);
      const int charge = seedTrack.charge();
      const float pt  = seedFitOk ? seedTrack.pt()  : 0;
      const float eta = seedFitOk ? seedTrack.eta() : 0;
      const float phi = seedFitOk ? seedTrack.phi() : 0;
      const int nHits = seedTrack.numberOfValidHits();

      const auto seedIndex = see_fitok.size();

      see_fitok   .push_back(seedFitOk);

      see_px      .push_back( seedFitOk ? seedTrack.px() : 0 );
      see_py      .push_back( seedFitOk ? seedTrack.py() : 0 );
      see_pz      .push_back( seedFitOk ? seedTrack.pz() : 0 );
      see_pt      .push_back( pt );
      see_eta     .push_back( eta );
      see_phi     .push_back( phi );
      see_q       .push_back( charge );
      see_nValid  .push_back( nHits );

      see_dxy     .push_back( seedFitOk ? seedTrack.dxy(bs.position()) : 0);
      see_dz      .push_back( seedFitOk ? seedTrack.dz(bs.position()) : 0);
      see_ptErr   .push_back( seedFitOk ? seedTrack.ptError() : 0);
      see_etaErr  .push_back( seedFitOk ? seedTrack.etaError() : 0);
      see_phiErr  .push_back( seedFitOk ? seedTrack.phiError() : 0);
      see_dxyErr  .push_back( seedFitOk ? seedTrack.dxyError() : 0);
      see_dzErr   .push_back( seedFitOk ? seedTrack.dzError() : 0);
      see_algo    .push_back( algo );

      see_trkIdx  .push_back(-1); // to be set correctly in fillTracks
      see_shareFrac.push_back( sharedFraction );
      see_simTrkIdx.push_back( tpIdx );

      /*
      TransientTrackingRecHit::RecHitPointer lastRecHit = theTTRHBuilder.build(&*(seed.recHits().second-1));
      TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( itSeed->startingState(), lastRecHit->surface(), theMF);
      float pt  = state.globalParameters().momentum().perp();
      float eta = state.globalParameters().momentum().eta();
      float phi = state.globalParameters().momentum().phi();
      see_px      .push_back( state.globalParameters().momentum().x() );
      see_py      .push_back( state.globalParameters().momentum().y() );
      see_pz      .push_back( state.globalParameters().momentum().z() );
      */

      std::vector<int> hitIdx;
      std::vector<int> hitType;

      for (auto hit=seed.recHits().first; hit!=seed.recHits().second; ++hit) {
    TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder.build(&*hit);
    int subid = recHit->geographicalId().subdetId();
    if (subid == (int) PixelSubdetector::PixelBarrel || subid == (int) PixelSubdetector::PixelEndcap) {
      const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*recHit);
          const auto& clusterRef = bhit->firstClusterRef();
          const auto clusterKey = clusterRef.cluster_pixel().key();
          if(includeAllHits_) {
            checkProductID(hitProductIds, clusterRef.id(), "seed");
            pix_seeIdx[clusterKey].push_back(seedIndex);
          }
          hitIdx.push_back( clusterKey );
          hitType.push_back( static_cast<int>(HitType::Pixel) );
    } else if (subid == (int) StripSubdetector::TOB || subid == (int) StripSubdetector::TID ||
               subid == (int) StripSubdetector::TIB || subid == (int) StripSubdetector::TEC) {
      if (trackerHitRTTI::isMatched(*recHit)) {
        const SiStripMatchedRecHit2D * matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(&*recHit);
            if(includeAllHits_) {
              checkProductID(hitProductIds, matchedHit->monoClusterRef().id(), "seed");
              checkProductID(hitProductIds, matchedHit->stereoClusterRef().id(), "seed");
            }
        int monoIdx = matchedHit->monoClusterRef().key();
        int stereoIdx = matchedHit->stereoClusterRef().key();

            std::vector<std::pair<int,int> >::const_iterator pos = find( monoStereoClusterList.begin(), monoStereoClusterList.end(), std::make_pair(monoIdx,stereoIdx) );
            const auto gluedIndex = std::distance(monoStereoClusterList.begin(), pos);
            if(includeAllHits_) glu_seeIdx[gluedIndex].push_back(seedIndex);
            hitIdx.push_back( gluedIndex );
            hitType.push_back( static_cast<int>(HitType::Glued) );
      } else {
        const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*recHit);
            const auto& clusterRef = bhit->firstClusterRef();
            unsigned int clusterKey;
            if(clusterRef.isPhase2()){
              clusterKey = clusterRef.cluster_phase2OT().key();
            } else {
              clusterKey = clusterRef.cluster_strip().key();
            }

            if(includeAllHits_) {
              checkProductID(hitProductIds, clusterRef.id(), "seed");
              if(clusterRef.isPhase2()){
                ph2_seeIdx[clusterKey].push_back(seedIndex);
              } else {
                str_seeIdx[clusterKey].push_back(seedIndex);
              }
            }

            hitIdx.push_back( clusterKey );
            if(clusterRef.isPhase2()){
              hitType.push_back( static_cast<int>(HitType::Phase2OT) );
            } else {
              hitType.push_back( static_cast<int>(HitType::Strip) );
            }
      }
    } else {
          LogTrace("TrackingNtuple") << " not pixel and not Strip detector";
        }
      }
      see_hitIdx   .push_back( hitIdx  );
      see_hitType  .push_back( hitType );
      see_nPixel   .push_back( std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Pixel)) );
      see_nGlued   .push_back( std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Glued)) );
      see_nStrip   .push_back( std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Strip)) );
      see_nPhase2OT.push_back( std::count(hitType.begin(), hitType.end(), static_cast<int>(HitType::Phase2OT)) );
      //the part below is not strictly needed
      float chi2 = -1;
      if (nHits==2) {
    TransientTrackingRecHit::RecHitPointer recHit0 = theTTRHBuilder.build(&*(seed.recHits().first));
    TransientTrackingRecHit::RecHitPointer recHit1 = theTTRHBuilder.build(&*(seed.recHits().first+1));
        std::vector<GlobalPoint> gp(2);
        std::vector<GlobalError> ge(2);
    gp[0] = recHit0->globalPosition();
    ge[0] = recHit0->globalPositionError();
    gp[1] = recHit1->globalPosition();
    ge[1] = recHit1->globalPositionError();
        LogTrace("TrackingNtuple") << "seed " << seedTrackRef.key()
                                   << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
                                   << " - PAIR - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId()
                                   << " hitpos: " << gp[0] << " " << gp[1]
                                   << " trans0: " << (recHit0->transientHits().size()>1 ? recHit0->transientHits()[0]->globalPosition() : GlobalPoint(0,0,0))
                                   << " " << (recHit0->transientHits().size()>1 ? recHit0->transientHits()[1]->globalPosition() : GlobalPoint(0,0,0))
                                   << " trans1: " << (recHit1->transientHits().size()>1 ? recHit1->transientHits()[0]->globalPosition() : GlobalPoint(0,0,0))
                                   << " " << (recHit1->transientHits().size()>1 ? recHit1->transientHits()[1]->globalPosition() : GlobalPoint(0,0,0))
                                   << " eta,phi: " << gp[0].eta() << "," << gp[0].phi();
      } else if (nHits==3) {
    TransientTrackingRecHit::RecHitPointer recHit0 = theTTRHBuilder.build(&*(seed.recHits().first));
    TransientTrackingRecHit::RecHitPointer recHit1 = theTTRHBuilder.build(&*(seed.recHits().first+1));
    TransientTrackingRecHit::RecHitPointer recHit2 = theTTRHBuilder.build(&*(seed.recHits().first+2));
    declareDynArray(GlobalPoint,4, gp);
    declareDynArray(GlobalError,4, ge);
    declareDynArray(bool,4, bl);
    gp[0] = recHit0->globalPosition();
    ge[0] = recHit0->globalPositionError();
    int subid0 = recHit0->geographicalId().subdetId();
    bl[0] = (subid0 == StripSubdetector::TIB || subid0 == StripSubdetector::TOB || subid0 == (int) PixelSubdetector::PixelBarrel);
    gp[1] = recHit1->globalPosition();
    ge[1] = recHit1->globalPositionError();
    int subid1 = recHit1->geographicalId().subdetId();
    bl[1] = (subid1 == StripSubdetector::TIB || subid1 == StripSubdetector::TOB || subid1 == (int) PixelSubdetector::PixelBarrel);
    gp[2] = recHit2->globalPosition();
    ge[2] = recHit2->globalPositionError();
    int subid2 = recHit2->geographicalId().subdetId();
    bl[2] = (subid2 == StripSubdetector::TIB || subid2 == StripSubdetector::TOB || subid2 == (int) PixelSubdetector::PixelBarrel);
    RZLine rzLine(gp,ge,bl);
    float seed_chi2 = rzLine.chi2();
    //float seed_pt = state.globalParameters().momentum().perp();
        float seed_pt = pt;
    LogTrace("TrackingNtuple") << "seed " << seedTrackRef.key()
                                   << " pt=" << pt << " eta=" << eta << " phi=" << phi << " q=" << charge
                                   << " - TRIPLET - ids: " << recHit0->geographicalId().rawId() << " " << recHit1->geographicalId().rawId() << " " << recHit2->geographicalId().rawId()
                                   << " hitpos: " << gp[0] << " " << gp[1] << " " << gp[2]
                                   << " trans0: " << (recHit0->transientHits().size()>1 ? recHit0->transientHits()[0]->globalPosition() : GlobalPoint(0,0,0))
                                   << " " << (recHit0->transientHits().size()>1 ? recHit0->transientHits()[1]->globalPosition() : GlobalPoint(0,0,0))
                                   << " trans1: " << (recHit1->transientHits().size()>1 ? recHit1->transientHits()[0]->globalPosition() : GlobalPoint(0,0,0))
                                   << " " << (recHit1->transientHits().size()>1 ? recHit1->transientHits()[1]->globalPosition() : GlobalPoint(0,0,0))
                                   << " trans2: " << (recHit2->transientHits().size()>1 ? recHit2->transientHits()[0]->globalPosition() : GlobalPoint(0,0,0))
                                   << " " << (recHit2->transientHits().size()>1 ? recHit2->transientHits()[1]->globalPosition() : GlobalPoint(0,0,0))
                                   << " local: " << recHit2->localPosition()
          //<< " tsos pos, mom: " << state.globalPosition()<<" "<<state.globalMomentum()
                                   << " eta,phi: " << gp[0].eta() << "," << gp[0].phi()
                                   << " pt,chi2: " << seed_pt << "," << seed_chi2;
    chi2 = seed_chi2;
      }
      see_chi2   .push_back( chi2 );
    }
  }
}

void TrackingNtuple::fillTracks(const edm::RefToBaseVector<reco::Track>& tracks,
                                const TrackingParticleRefVector& tpCollection,
                                const TrackingParticleRefKeyToIndex& tpKeyToIndex,
                                const reco::BeamSpot& bs,
                                const reco::TrackToTrackingParticleAssociator& associatorByHits,
                                const TransientTrackingRecHitBuilder& theTTRHBuilder,
                                const TrackerTopology& tTopo,
                                const std::set<edm::ProductID>& hitProductIds,
                                const std::map<edm::ProductID, size_t>& seedCollToOffset
                                ) {
  reco::RecoToSimCollection recSimColl = associatorByHits.associateRecoToSim(tracks, tpCollection);
  edm::EDConsumerBase::Labels labels;
  labelsForToken(trackToken_, labels);
  LogTrace("TrackingNtuple") << "NEW TRACK LABEL: " << labels.module;
  for(size_t iTrack = 0; iTrack<tracks.size(); ++iTrack) {
    const auto& itTrack = tracks[iTrack];
    int nSimHits = 0;
    bool isSimMatched = false;
    std::vector<float> sharedFraction;
    std::vector<int> tpIdx;
    auto foundTPs = recSimColl.find(itTrack);
    if (foundTPs != recSimColl.end()) {
      if (!foundTPs->val.empty()) {
    nSimHits = foundTPs->val[0].first->numberOfTrackerHits();
    isSimMatched = true;
      }
      for(const auto tpQuality: foundTPs->val) {
        sharedFraction.push_back(tpQuality.second);
    tpIdx.push_back( tpKeyToIndex.at( tpQuality.first.key() ) );
      }
    }
    int charge = itTrack->charge();
    float pt = itTrack->pt();
    float eta = itTrack->eta();
    const double lambda = itTrack->lambda();
    float chi2 = itTrack->normalizedChi2();
    float phi = itTrack->phi();
    int nHits = itTrack->numberOfValidHits();
    const reco::HitPattern& hp = itTrack->hitPattern();
    trk_px       .push_back(itTrack->px());
    trk_py       .push_back(itTrack->py());
    trk_pz       .push_back(itTrack->pz());
    trk_pt       .push_back(pt);
    trk_inner_px.push_back(itTrack->innerMomentum().x());
    trk_inner_py.push_back(itTrack->innerMomentum().y());
    trk_inner_pz.push_back(itTrack->innerMomentum().z());
    trk_inner_pt.push_back(itTrack->innerMomentum().rho());
    trk_outer_px.push_back(itTrack->outerMomentum().x());
    trk_outer_py.push_back(itTrack->outerMomentum().y());
    trk_outer_pz.push_back(itTrack->outerMomentum().z());
    trk_outer_pt.push_back(itTrack->outerMomentum().rho());
    trk_eta      .push_back(eta);
    trk_lambda   .push_back(lambda);
    trk_cotTheta .push_back(1/tan(M_PI*0.5-lambda));
    trk_phi      .push_back(phi);
    trk_dxy      .push_back(itTrack->dxy(bs.position()));
    trk_dz       .push_back(itTrack->dz(bs.position()));
    trk_ptErr    .push_back(itTrack->ptError());
    trk_etaErr   .push_back(itTrack->etaError());
    trk_lambdaErr.push_back(itTrack->lambdaError());
    trk_phiErr   .push_back(itTrack->phiError());
    trk_dxyErr   .push_back(itTrack->dxyError());
    trk_dzErr    .push_back(itTrack->dzError());
    trk_refpoint_x.push_back(itTrack->vx());
    trk_refpoint_y.push_back(itTrack->vy());
    trk_refpoint_z.push_back(itTrack->vz());
    trk_nChi2    .push_back( itTrack->normalizedChi2());
    trk_shareFrac.push_back(sharedFraction);
    trk_q        .push_back(charge);
    trk_nValid   .push_back(hp.numberOfValidHits());
    trk_nInvalid .push_back(hp.numberOfLostHits(reco::HitPattern::TRACK_HITS));
    trk_nPixel   .push_back(hp.numberOfValidPixelHits());
    trk_nStrip   .push_back(hp.numberOfValidStripHits());
    trk_nPixelLay.push_back(hp.pixelLayersWithMeasurement());
    trk_nStripLay.push_back(hp.stripLayersWithMeasurement());
    trk_n3DLay   .push_back(hp.numberOfValidStripLayersWithMonoAndStereo()+hp.pixelLayersWithMeasurement());
    trk_nOuterLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
    trk_nInnerLost.push_back(hp.numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
    trk_algo     .push_back(itTrack->algo());
    trk_originalAlgo.push_back(itTrack->originalAlgo());
    trk_algoMask .push_back(itTrack->algoMaskUL());
    trk_stopReason.push_back(itTrack->stopReason());
    trk_isHP     .push_back(itTrack->quality(reco::TrackBase::highPurity));
    if(includeSeeds_) {
      auto offset = seedCollToOffset.find(itTrack->seedRef().id());
      if(offset == seedCollToOffset.end()) {
        throw cms::Exception("Configuration") << "Track algo '" << reco::TrackBase::algoName(itTrack->algo())
                                              << "' originalAlgo '" << reco::TrackBase::algoName(itTrack->originalAlgo())
                                              << "' refers to seed collection " << itTrack->seedRef().id()
                                              << ", but that seed collection is not given as an input. The following collections were given as an input " << make_ProductIDMapPrinter(seedCollToOffset);
      }

      const auto seedIndex = offset->second + itTrack->seedRef().key();
      trk_seedIdx  .push_back(seedIndex);
      if(see_trkIdx[seedIndex] != -1) {
        throw cms::Exception("LogicError") << "Track index has already been set for seed " << seedIndex << " to " << see_trkIdx[seedIndex] << "; was trying to set it to " << iTrack;
      }
      see_trkIdx[seedIndex] = iTrack;
    }
    trk_vtxIdx   .push_back(-1); // to be set correctly in fillVertices
    trk_simTrkIdx.push_back(tpIdx);
    LogTrace("TrackingNtuple") << "Track #" << itTrack.key() << " with q=" << charge
                               << ", pT=" << pt << " GeV, eta: " << eta << ", phi: " << phi
                               << ", chi2=" << chi2
                               << ", Nhits=" << nHits
                               << ", algo=" << itTrack->algoName(itTrack->algo()).c_str()
                               << " hp=" << itTrack->quality(reco::TrackBase::highPurity)
                               << " seed#=" << itTrack->seedRef().key()
                               << " simMatch=" << isSimMatched
                               << " nSimHits=" << nSimHits
                               << " sharedFraction=" << (sharedFraction.empty()?-1:sharedFraction[0])
                               << " tpIdx=" << (tpIdx.empty()?-1:tpIdx[0]);
    std::vector<int> hitIdx;
    std::vector<int> hitType;

    for(auto i=itTrack->recHitsBegin(); i!=itTrack->recHitsEnd(); i++) {
      TransientTrackingRecHit::RecHitPointer hit=theTTRHBuilder.build(&**i );
      DetId hitId = hit->geographicalId();
      LogTrace("TrackingNtuple") << "hit #" << std::distance(itTrack->recHitsBegin(), i) << " subdet=" << hitId.subdetId();
      if(hitId.det() != DetId::Tracker)
        continue;

      LogTrace("TrackingNtuple") << " " << subdetstring(hitId.subdetId()) << " " << tTopo.layer(hitId);

      if (hit->isValid()) {
        //ugly... but works
        const BaseTrackerRecHit* bhit = dynamic_cast<const BaseTrackerRecHit*>(&*hit);
        const auto& clusterRef = bhit->firstClusterRef();
        unsigned int clusterKey;
        if(clusterRef.isPixel()){
          clusterKey = clusterRef.cluster_pixel().key();
        } else if(clusterRef.isPhase2()){
          clusterKey = clusterRef.cluster_phase2OT().key();
        } else {
          clusterKey = clusterRef.cluster_strip().key();
        }

        LogTrace("TrackingNtuple") << " id: " << hitId.rawId() << " - globalPos =" << hit->globalPosition()
                                   << " cluster=" << clusterKey
                                   << " clusterRef ID=" << clusterRef.id()
                                   << " eta,phi: " << hit->globalPosition().eta() << "," << hit->globalPosition().phi();
        if(includeAllHits_) {
          checkProductID(hitProductIds, clusterRef.id(), "track");
          if(clusterRef.isPixel()){
            pix_trkIdx[clusterKey].push_back(iTrack);
          } else if(clusterRef.isPhase2()){
            ph2_trkIdx[clusterKey].push_back(iTrack);
          } else {
            str_trkIdx[clusterKey].push_back(iTrack);
          }

        }

        hitIdx.push_back(clusterKey);
        if(clusterRef.isPixel()){
          hitType.push_back( static_cast<int>(HitType::Pixel));
        } else if(clusterRef.isPhase2()){
          hitType.push_back( static_cast<int>(HitType::Phase2OT));
        } else {
          hitType.push_back( static_cast<int>(HitType::Strip));
        }
      } else  {
        LogTrace("TrackingNtuple") << " - invalid hit";

        hitIdx.push_back( inv_isBarrel.size() );
        hitType.push_back( static_cast<int>(HitType::Invalid) );

        inv_isBarrel.push_back( hitId.subdetId()==1 );
        inv_det     .push_back( hitId.subdetId() );
        inv_lay     .push_back( tTopo.layer(hitId) );
        inv_detId   .push_back( hitId.rawId() );
        inv_type    .push_back( hit->getType() );

      }
    }

    trk_hitIdx.push_back(hitIdx);
    trk_hitType.push_back(hitType);
  }
}


void TrackingNtuple::fillTrackingParticles(const edm::Event& iEvent, const edm::EventSetup& iSetup,
                                           const edm::RefToBaseVector<reco::Track>& tracks,
                                           const TrackingParticleRefVector& tpCollection,
                                           const TrackingVertexRefKeyToIndex& tvKeyToIndex,
                                           const reco::TrackToTrackingParticleAssociator& associatorByHits,
                                           const std::vector<TPHitIndex>& tpHitList
                                           ) {
  edm::ESHandle<ParametersDefinerForTP> parametersDefinerH;
  iSetup.get<TrackAssociatorRecord>().get(parametersDefinerName_, parametersDefinerH);
  const ParametersDefinerForTP *parametersDefiner = parametersDefinerH.product();

  // Number of 3D layers for TPs
  edm::Handle<edm::ValueMap<unsigned int>> tpNLayersH;
  iEvent.getByToken(tpNLayersToken_, tpNLayersH);
  const auto& nLayers_tPCeff = *tpNLayersH;

  iEvent.getByToken(tpNPixelLayersToken_, tpNLayersH);
  const auto& nPixelLayers_tPCeff = *tpNLayersH;

  iEvent.getByToken(tpNStripStereoLayersToken_, tpNLayersH);
  const auto& nStripMonoAndStereoLayers_tPCeff = *tpNLayersH;

  reco::SimToRecoCollection simRecColl = associatorByHits.associateSimToReco(tracks, tpCollection);

  for(const TrackingParticleRef& tp: tpCollection) {
    LogTrace("TrackingNtuple") << "tracking particle pt=" << tp->pt() << " eta=" << tp->eta() << " phi=" << tp->phi();
    bool isRecoMatched = false;
    std::vector<int> tkIdx;
    std::vector<float> sharedFraction;
    auto foundTracks = simRecColl.find(tp);
    if(foundTracks != simRecColl.end()) {
      isRecoMatched = true;
      for(const auto trackQuality: foundTracks->val) {
        sharedFraction.push_back(trackQuality.second);
        tkIdx.push_back(trackQuality.first.key());
      }
    }
    LogTrace("TrackingNtuple") << "matched to tracks = " << make_VectorPrinter(tkIdx) << " isRecoMatched=" << isRecoMatched;
    sim_event    .push_back(tp->eventId().event());
    sim_bunchCrossing.push_back(tp->eventId().bunchCrossing());
    sim_pdgId    .push_back(tp->pdgId());
    sim_px       .push_back(tp->px());
    sim_py       .push_back(tp->py());
    sim_pz       .push_back(tp->pz());
    sim_pt       .push_back(tp->pt());
    sim_eta      .push_back(tp->eta());
    sim_phi      .push_back(tp->phi());
    sim_q        .push_back(tp->charge());
    sim_trkIdx   .push_back(tkIdx);
    sim_shareFrac.push_back(sharedFraction);
    sim_parentVtxIdx.push_back( tvKeyToIndex.at(tp->parentVertex().key()) );
    std::vector<int> decayIdx;
    for(const auto& v: tp->decayVertices())
      decayIdx.push_back( tvKeyToIndex.at(v.key()) );
    sim_decayVtxIdx.push_back(decayIdx);

    //Calcualte the impact parameters w.r.t. PCA
    TrackingParticle::Vector momentum = parametersDefiner->momentum(iEvent,iSetup,tp);
    TrackingParticle::Point vertex = parametersDefiner->vertex(iEvent,iSetup,tp);
    float dxySim = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
    float dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y())/sqrt(momentum.perp2())
      * momentum.z()/sqrt(momentum.perp2());
    const double lambdaSim = M_PI/2 - momentum.theta();
    sim_pca_pt       .push_back(std::sqrt(momentum.perp2()));
    sim_pca_eta      .push_back(momentum.Eta());
    sim_pca_lambda   .push_back(lambdaSim);
    sim_pca_cotTheta .push_back(1/tan(M_PI*0.5-lambdaSim));
    sim_pca_phi      .push_back(momentum.phi());
    sim_pca_dxy      .push_back(dxySim);
    sim_pca_dz       .push_back(dzSim);

    std::vector<int> hitIdx;
    int nPixel=0, nStrip=0;
    auto rangeHit = std::equal_range(tpHitList.begin(), tpHitList.end(), TPHitIndex(tp.key()), tpHitIndexListLess);
    for(auto ip = rangeHit.first; ip != rangeHit.second; ++ip) {
      auto type = HitType::Unknown;
      if(!simhit_hitType[ip->simHitIdx].empty())
        type = static_cast<HitType>(simhit_hitType[ip->simHitIdx][0]);
      LogTrace("TrackingNtuple") << "simhit=" << ip->simHitIdx << " type=" << static_cast<int>(type);
      hitIdx.push_back(ip->simHitIdx);
      const auto detid = DetId(simhit_detId[ip->simHitIdx]);
      if(detid.det() != DetId::Tracker) {
        throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId " << detid.rawId() << " whose det() is not Tracker but " << detid.det();
      }
      const auto subdet = detid.subdetId();
      switch(subdet) {
      case PixelSubdetector::PixelBarrel:
      case PixelSubdetector::PixelEndcap:
        ++nPixel;
        break;
      case StripSubdetector::TIB:
      case StripSubdetector::TID:
      case StripSubdetector::TOB:
      case StripSubdetector::TEC:
        ++nStrip;
        break;
      default:
        throw cms::Exception("LogicError") << "Encountered SimHit for TP " << tp.key() << " with DetId " << detid.rawId() << " whose subdet is not recognized, is " << subdet;
      };
    }
    sim_nValid   .push_back( hitIdx.size() );
    sim_nPixel   .push_back( nPixel );
    sim_nStrip   .push_back( nStrip );

    const auto nSimLayers = nLayers_tPCeff[tp];
    const auto nSimPixelLayers = nPixelLayers_tPCeff[tp];
    const auto nSimStripMonoAndStereoLayers = nStripMonoAndStereoLayers_tPCeff[tp];
    sim_nLay     .push_back( nSimLayers );
    sim_nPixelLay.push_back( nSimPixelLayers );
    sim_n3DLay   .push_back( nSimPixelLayers+nSimStripMonoAndStereoLayers );

    sim_simHitIdx.push_back(hitIdx);
  }
}

void TrackingNtuple::fillVertices(const reco::VertexCollection& vertices,
                                  const edm::RefToBaseVector<reco::Track>& tracks) {
  for(size_t iVertex=0, size=vertices.size(); iVertex<size; ++iVertex) {
    const reco::Vertex& vertex = vertices[iVertex];
    vtx_x.push_back(vertex.x());
    vtx_y.push_back(vertex.y());
    vtx_z.push_back(vertex.z());
    vtx_xErr.push_back(vertex.xError());
    vtx_yErr.push_back(vertex.yError());
    vtx_zErr.push_back(vertex.zError());
    vtx_chi2.push_back(vertex.chi2());
    vtx_ndof.push_back(vertex.ndof());
    vtx_fake.push_back(vertex.isFake());
    vtx_valid.push_back(vertex.isValid());

    std::vector<int> trkIdx;
    for(auto iTrack = vertex.tracks_begin(); iTrack != vertex.tracks_end(); ++iTrack) {
      // Ignore link if vertex was fit from a track collection different from the input
      if(iTrack->id() != tracks.id())
        continue;

      trkIdx.push_back(iTrack->key());

      if(trk_vtxIdx[iTrack->key()] != -1) {
        throw cms::Exception("LogicError") << "Vertex index has already been set for track " << iTrack->key() << " to " << trk_vtxIdx[iTrack->key()] << "; was trying to set it to " << iVertex;
      }
      trk_vtxIdx[iTrack->key()] = iVertex;
    }
    vtx_trkIdx.push_back(trkIdx);
  }
}

void TrackingNtuple::fillTrackingVertices(const TrackingVertexRefVector& trackingVertices,
                                          const TrackingParticleRefKeyToIndex& tpKeyToIndex
                                          ) {
  int current_event = -1;
  for(const auto& ref: trackingVertices) {
    const TrackingVertex v = *ref;
    if(v.eventId().event() != current_event) {
      // next PV
      current_event = v.eventId().event();
      simpv_idx.push_back(simvtx_x.size());
    }

    unsigned int processType = std::numeric_limits<unsigned int>::max();
    if(!v.g4Vertices().empty()) {
      processType = v.g4Vertices()[0].processType();
    }

    simvtx_event.push_back(v.eventId().event());
    simvtx_bunchCrossing.push_back(v.eventId().bunchCrossing());
    simvtx_processType.push_back(processType);

    simvtx_x.push_back(v.position().x());
    simvtx_y.push_back(v.position().y());
    simvtx_z.push_back(v.position().z());

    auto fill = [&](const TrackingParticleRefVector& tps, std::vector<int>& idx) {
      for(const auto& tpRef: tps) {
        auto found = tpKeyToIndex.find(tpRef.key());
        if(found != tpKeyToIndex.end()) {
          idx.push_back(tpRef.key());
        }
      }
    };

    std::vector<int> sourceIdx;
    std::vector<int> daughterIdx;
    fill(v.sourceTracks(), sourceIdx);
    fill(v.daughterTracks(), daughterIdx);

    simvtx_sourceSimIdx.push_back(sourceIdx);
    simvtx_daughterSimIdx.push_back(daughterIdx);
  }
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void TrackingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.addUntracked<std::vector<edm::InputTag> >("seedTracks", std::vector<edm::InputTag>{
      edm::InputTag("seedTracksinitialStepSeeds"),
      edm::InputTag("seedTracksdetachedTripletStepSeeds"),
      edm::InputTag("seedTrackspixelPairStepSeeds"),
      edm::InputTag("seedTrackslowPtTripletStepSeeds"),
      edm::InputTag("seedTracksmixedTripletStepSeeds"),
      edm::InputTag("seedTrackspixelLessStepSeeds"),
      edm::InputTag("seedTrackstobTecStepSeeds"),
      edm::InputTag("seedTracksjetCoreRegionalStepSeeds"),
      edm::InputTag("seedTracksmuonSeededSeedsInOut"),
      edm::InputTag("seedTracksmuonSeededSeedsOutIn")
  });
  desc.addUntracked<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
  desc.addUntracked<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
  desc.addUntracked<bool>("trackingParticlesRef", false);
  desc.addUntracked<edm::InputTag>("clusterTPMap", edm::InputTag("tpClusterProducer"));
  desc.addUntracked<edm::InputTag>("simHitTPMap", edm::InputTag("simHitTPAssocProducer"));
  desc.addUntracked<edm::InputTag>("trackAssociator", edm::InputTag("quickTrackAssociatorByHits"));
  desc.addUntracked<edm::InputTag>("pixelDigiSimLink", edm::InputTag("simSiPixelDigis"));
  desc.addUntracked<edm::InputTag>("stripDigiSimLink", edm::InputTag("simSiStripDigis"));
  desc.addUntracked<edm::InputTag>("phase2OTSimLink", edm::InputTag(""));
  desc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
  desc.addUntracked<edm::InputTag>("pixelRecHits", edm::InputTag("siPixelRecHits"));
  desc.addUntracked<edm::InputTag>("stripRphiRecHits", edm::InputTag("siStripMatchedRecHits", "rphiRecHit"));
  desc.addUntracked<edm::InputTag>("stripStereoRecHits", edm::InputTag("siStripMatchedRecHits", "stereoRecHit"));
  desc.addUntracked<edm::InputTag>("stripMatchedRecHits", edm::InputTag("siStripMatchedRecHits", "matchedRecHit"));
  desc.addUntracked<edm::InputTag>("phase2OTRecHits", edm::InputTag("siPhase2RecHits"));
  desc.addUntracked<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
  desc.addUntracked<edm::InputTag>("trackingVertices", edm::InputTag("mix", "MergedTrackTruth"));
  desc.addUntracked<edm::InputTag>("trackingParticleNlayers", edm::InputTag("trackingParticleNumberOfLayersProducer", "trackerLayers"));
  desc.addUntracked<edm::InputTag>("trackingParticleNpixellayers", edm::InputTag("trackingParticleNumberOfLayersProducer", "pixelLayers"));
  desc.addUntracked<edm::InputTag>("trackingParticleNstripstereolayers", edm::InputTag("trackingParticleNumberOfLayersProducer", "stripStereoLayers"));
  desc.addUntracked<std::string>("TTRHBuilder", "WithTrackAngle");
  desc.addUntracked<std::string>("parametersDefiner", "LhcParametersDefinerForTP");
  desc.addUntracked<bool>("includeSeeds", false);
  desc.addUntracked<bool>("includeAllHits", false);
  descriptions.add("trackingNtuple",desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(TrackingNtuple);