CMS 3D CMS Logo

HGCalTrackCollectionProducer.cc
Go to the documentation of this file.
1 // S. Zenz, 12 February 2015
2 //
3 // Splits a track collection into two, based on whether they propagate to the HGCal or not
4 // Tracks with bad pt resolution (suspected fakes) are dropped and not in either collection
5 
7 
13 
17 
18 //#include "DataFormats/TrackReco/interface/Track.h"
22 
23 
24 // for track propagation through HGC
25 // N.B. we are only propogating to first layer, so check these later
33 
34 //geometry records
38 
41 
42 #include<unordered_map>
43 
45 
46 public:
48 private:
49  void produce( edm::Event &, const edm::EventSetup & ) override;
51  const edm::EventSetup&) override;
52 
54 
55  // variables needed for copied goodPtResolution function
56  // need to go back and figure out sensible values
57  bool debug_;
58  const std::vector<double> DPtovPtCut_;
59  const std::vector<unsigned> NHitCut_;
60  const bool useIterTracking_;
61  // const bool _useFirstLayerOnly; // always true now
62 
63  // variables needed for copied extrapolation
66  std::array<std::string,1> hgc_names_; // 3 --> 1; extrapolate to hgcee only
67  std::array<edm::ESHandle<HGCalGeometry>,1> hgcGeometries_; // 3 --> 1; extrapolate to hgcee only
68  std::array<std::vector<ReferenceCountingPointer<BoundDisk> >,1> plusSurface_,minusSurface_; // 3 --> 1; extrapolate to hgcee only
69  std::unique_ptr<PropagatorWithMaterial> mat_prop_;
70 
73 };
74 
76  src_(consumes<edm::View<reco::PFRecTrack> >(iConfig.getParameter<edm::InputTag> ("src"))),
77  debug_(iConfig.getParameter<bool>("debug")),
78  DPtovPtCut_(iConfig.getParameter<std::vector<double> >("DPtOverPtCuts_byTrackAlgo")),
79  NHitCut_(iConfig.getParameter<std::vector<unsigned> >("NHitCuts_byTrackAlgo")),
80  useIterTracking_(iConfig.getParameter<bool>("useIterativeTracking"))
81 {
82 
83  LogDebug("HGCalTrackCollectionProducer") << " HGCalTrackCollectionProducer::HGCalTrackCollectionProducer " << std::endl;
84 
85  const edm::ParameterSet& geoconf = iConfig.getParameterSet("hgcalGeometryNames");
86  hgc_names_[0] = geoconf.getParameter<std::string>("HGC_ECAL");
87 
88  produces<reco::PFRecTrackCollection>("TracksInHGCal");
89  produces<reco::PFRecTrackCollection>("TracksNotInHGCal");
90 
91 }
92 
93 // From https://github.com/cms-sw/cmssw/blob/CMSSW_6_2_X_SLHC/RecoParticleFlow/PFClusterProducer/src/HGCClusterizer.cc#L441-L447 and beyond
94 // TODO: we only need the front of the calorimeter, so modify this
96  constexpr float m_pion = 0.1396;
97  // get dependencies for setting up propagator
100  // get HGC geometries (assume that layers are ordered in Z!)
101  for( unsigned i = 0; i < hgcGeometries_.size(); ++i ) {
103  }
104 
105  // make propagator
107  // setup HGC layers for track propagation
108  Surface::RotationType rot; //unit rotation matrix
109  for( unsigned i = 0; i < hgcGeometries_.size(); ++i ) {
110  minusSurface_[i].clear();
111  plusSurface_[i].clear();
112  const HGCalDDDConstants &dddCons=hgcGeometries_[i]->topology().dddConstants();
113  std::map<float,float> zrhoCoord;
114  std::map<float,float> innerRadiusCoord;
115  auto theTrForms = dddCons.getTrForms();
116  const auto& firstLayerIt = theTrForms.back();
117  float Z(std::abs(firstLayerIt.h3v.z()));
118  // use hardcoded radii for now (FIX ME)
119  diskInnerRadius_ = 31.5;
120  diskOuterRadius_ = 161.0f;
121  LogDebug("HGCalTrackCollectionProducer") << "O HAI I'm making a bound disk with Outer R=" << diskOuterRadius_ << " Inner R=" << diskInnerRadius_ << " and Z=" << Z << std::endl;
123  new SimpleDiskBounds( diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
125  new SimpleDiskBounds( diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
126  }
127 }
128 
130 
132  evt.getByToken(src_,trackHandle);
133  const auto& tracks = *trackHandle;
134 
135  auto outputInHGCal = std::make_unique<reco::PFRecTrackCollection>();
136  auto outputNotInHGCal = std::make_unique<reco::PFRecTrackCollection>();
137 
138  for ( unsigned int i = 0 ; i < tracks.size() ; i++) {
139  const auto track = tracks.ptrAt(i);
141  LogDebug("HGCalTrackCollectionProducer") << "HGCalTrackCollectionProducer Track number " << i << " has a goodPtResolution result of " << isGood << std::endl;
142  if (!isGood) continue;
143  bool found = false;
145  auto detbegin = myTSOS.globalPosition().z() > 0 ? plusSurface_.begin() : minusSurface_.begin();
146  auto detend = myTSOS.globalPosition().z() > 0 ? plusSurface_.end() : minusSurface_.end();
147  for( auto det = detbegin; det != detend; ++det ) {
148  LogDebug("HGCalTrackCollectionProducer") << "at HGC detector: " << std::distance(detbegin,det) << std::endl;
149  unsigned layer_count = 1;
150  for( const auto& layer : *det ) {
151  LogDebug("HGCalTrackCollectionProducer") << " at DET layer: " << layer_count++ << std::endl;
152  TrajectoryStateOnSurface piStateAtSurface = mat_prop_->propagate(myTSOS, *layer);
153  if( piStateAtSurface.isValid() ) {
154  LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is valid!" << std::endl;
155  GlobalPoint pt = piStateAtSurface.globalPosition();
156  if (pt.perp() < diskOuterRadius_) {
157  if (pt.perp() > diskInnerRadius_) {
158  LogDebug("HGCalTrackCollectionProducer") << "(x,y,z,r)=(" << pt.x() << ", " << pt.y() << ", " << pt.z() << ", " << sqrt(pt.x()*pt.x() + pt.y()*pt.y()) << ")" << std::endl;
159  if (std::abs(track->trackRef()->eta()) < 1.47) LogDebug("HGCalTrackCollectionProducer") << " ETA IN BARREL REGION: " << track->trackRef()->eta()
160  << " (PT: " << track->trackRef()->pt() << ")" << std::endl;
161  found = true;
162  } else {
163  LogDebug("HGCalTrackCollectionProducer") << " but r=" << pt.perp() << " < diskInnerRadius=" << diskInnerRadius_ << " so skipping " << std::endl;
164  }
165  } else {
166  LogDebug("HGCalTrackCollectionProducer") << " but r=" << pt.perp() << " > diskOuterRadius=" << diskOuterRadius_ << " so skipping " << std::endl;
167  }
168  } else {
169  LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is NOT valid!" << std::endl;
170  }
171  }
172  }
173  if (found) {
174  LogDebug("HGCalTrackCollectionProducer") << " Track going to outputInHGCal pt eta " << track->trackRef()->pt() << " " << track->trackRef()->eta() << std::endl;
175  outputInHGCal->push_back(*track);
176  } else {
177  outputNotInHGCal->push_back(*track);
178  }
179  } // track loop
180 
181  evt.put(std::move(outputInHGCal),"TracksInHGCal");
182  evt.put(std::move(outputNotInHGCal),"TracksNotInHGCal");
183 }
184 
#define LogDebug(id)
T getParameter(std::string const &) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:137
T perp() const
Definition: PV3DBase.h:72
std::array< std::vector< ReferenceCountingPointer< BoundDisk > >, 1 > minusSurface_
TrajectoryStateOnSurface outerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field, bool withErr=true)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:579
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
const std::vector< unsigned > NHitCut_
std::array< std::string, 1 > hgc_names_
#define constexpr
edm::EDGetTokenT< edm::View< reco::PFRecTrack > > src_
std::array< std::vector< ReferenceCountingPointer< BoundDisk > >, 1 > plusSurface_
std::vector< HGCalParameters::hgtrform > getTrForms() const
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
std::unique_ptr< PropagatorWithMaterial > mat_prop_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void beginLuminosityBlock(const edm::LuminosityBlock &, const edm::EventSetup &) override
edm::ESHandle< MagneticField > bField_
HGCalTrackCollectionProducer(const edm::ParameterSet &)
ParameterSet const & getParameterSet(std::string const &) const
Disk BoundDisk
Definition: BoundDisk.h:62
const std::vector< double > DPtovPtCut_
fixed size matrix
HLT enums.
T get() const
Definition: EventSetup.h:63
std::array< edm::ESHandle< HGCalGeometry >, 1 > hgcGeometries_
edm::ESHandle< TrackerGeometry > tkGeom_
T x() const
Definition: PV3DBase.h:62
T const * product() const
Definition: ESHandle.h:86
bool goodPtResolution(const reco::TrackRef &, const std::vector< double > &DPtovPtCut, const std::vector< unsigned > &NHitCut, bool useIterTracking, bool debug=false)
void produce(edm::Event &, const edm::EventSetup &) override
def move(src, dest)
Definition: eostools.py:510