CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LowPtGsfElectronSCProducer.cc
Go to the documentation of this file.
17 #include <iostream>
18 
20 //
22  : gsfPfRecTracks_{consumes<reco::GsfPFRecTrackCollection>(cfg.getParameter<edm::InputTag>("gsfPfRecTracks"))},
23  ecalClusters_{consumes<reco::PFClusterCollection>(cfg.getParameter<edm::InputTag>("ecalClusters"))},
24  dr2_{cfg.getParameter<double>("MaxDeltaR2")} {
25  produces<reco::CaloClusterCollection>();
26  produces<reco::SuperClusterCollection>();
27  produces<edm::ValueMap<reco::SuperClusterRef> >();
28 }
29 
31 //
33 
35 //
37  // Input GsfPFRecTracks collection
39 
40  // Input EcalClusters collection
42 
43  // Output SuperClusters collection and getRefBeforePut
44  auto superClusters = std::make_unique<reco::SuperClusterCollection>(reco::SuperClusterCollection());
45  superClusters->reserve(gsfPfRecTracks->size());
46  const reco::SuperClusterRefProd superClustersRefProd = event.getRefBeforePut<reco::SuperClusterCollection>();
47 
48  // Output ValueMap container of GsfPFRecTrackRef index to SuperClusterRef
49  std::vector<reco::SuperClusterRef> superClustersValueMap;
50 
51  // Output CaloClusters collection
52  auto caloClusters = std::make_unique<reco::CaloClusterCollection>(reco::CaloClusterCollection());
53  caloClusters->reserve(ecalClusters.size());
54 
55  // Index[GSF track][trajectory point] for "best" CaloCluster
56  std::vector<std::vector<int> > cluster_idx;
57  cluster_idx.reserve(gsfPfRecTracks->size());
58 
59  // Index[GSF track][trajectory point] for "best" PFCluster
60  std::vector<std::vector<int> > pfcluster_idx;
61  pfcluster_idx.reserve(gsfPfRecTracks->size());
62 
63  // dr2min[GSF track][trajectory point] for "best" CaloCluster
64  std::vector<std::vector<float> > cluster_dr2min;
65  cluster_dr2min.reserve(gsfPfRecTracks->size());
66 
67  // Construct list of trajectory points from the GSF track and electron brems
68  std::vector<std::vector<const reco::PFTrajectoryPoint*> > points;
69  points.reserve(gsfPfRecTracks->size());
70  for (auto const& trk : *gsfPfRecTracks) {
71  // Extrapolated track
72  std::vector<const reco::PFTrajectoryPoint*> traj;
73  traj.reserve(trk.PFRecBrem().size() + 1);
74  traj.push_back(&trk.extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::ECALShowerMax));
75  // Extrapolated brem trajectories
76  for (auto const& brem : trk.PFRecBrem()) {
77  traj.push_back(&brem.extrapolatedPoint(reco::PFTrajectoryPoint::LayerType::ECALShowerMax));
78  }
79  auto size = traj.size();
80  points.push_back(std::move(traj));
81  // Size containers
82  cluster_idx.emplace_back(size, -1);
83  pfcluster_idx.emplace_back(size, -1);
84  cluster_dr2min.emplace_back(size, 1.e6);
85  }
86 
87  // For each cluster, find closest trajectory point ("global" arbitration at event level)
88  for (size_t iclu = 0; iclu < ecalClusters.size(); ++iclu) { // Cluster loop
89  std::pair<int, int> point = std::make_pair(-1, -1);
90  float dr2min = 1.e6;
91  for (size_t ipoint = 0; ipoint < points.size(); ++ipoint) { // GSF track loop
92  for (size_t jpoint = 0; jpoint < points[ipoint].size(); ++jpoint) { // Traj point loop
93  if (points[ipoint][jpoint]->isValid()) {
94  float dr2 = reco::deltaR2(ecalClusters[iclu], points[ipoint][jpoint]->positionREP());
95  if (dr2 < dr2min) {
96  // Store nearest point to this cluster
97  dr2min = dr2;
98  point = std::make_pair(ipoint, jpoint);
99  }
100  }
101  }
102  }
103  if (point.first >= 0 && point.second >= 0 && // if this cluster is matched to a point ...
104  dr2min < cluster_dr2min[point.first][point.second]) { // ... and cluster is closest to the same point
105  // Copy CaloCluster to new collection
106  caloClusters->push_back(ecalClusters[iclu]);
107  // Store cluster index for creation of SC later
108  cluster_idx[point.first][point.second] = caloClusters->size() - 1;
109  pfcluster_idx[point.first][point.second] = iclu;
110  cluster_dr2min[point.first][point.second] = dr2min;
111  }
112  }
113 
114  // Put CaloClusters in event and get orphan handle
115  const edm::OrphanHandle<reco::CaloClusterCollection>& caloClustersH = event.put(std::move(caloClusters));
116 
117  // Loop through GSF tracks
118  for (size_t itrk = 0; itrk < gsfPfRecTracks->size(); ++itrk) {
119  // Used to create SC
120  float energy = 0.;
121  float X = 0., Y = 0., Z = 0.;
124  std::vector<const reco::PFCluster*> barePtrs;
125 
126  // Find closest match in dr2 from points associated to given track
127  int index = -1;
128  float dr2 = 1.e6;
129  for (size_t ipoint = 0; ipoint < cluster_idx[itrk].size(); ++ipoint) {
130  if (cluster_idx[itrk][ipoint] < 0) {
131  continue;
132  }
133  if (cluster_dr2min[itrk][ipoint] < dr2) {
134  dr2 = cluster_dr2min[itrk][ipoint];
135  index = cluster_idx[itrk][ipoint];
136  }
137  }
138 
139  // For each track, loop through points and use associated cluster
140  for (size_t ipoint = 0; ipoint < cluster_idx[itrk].size(); ++ipoint) {
141  if (cluster_idx[itrk][ipoint] < 0) {
142  continue;
143  }
144  reco::CaloClusterPtr clu(caloClustersH, cluster_idx[itrk][ipoint]);
145  if (clu.isNull()) {
146  continue;
147  }
148  if (!(cluster_dr2min[itrk][ipoint] < dr2_ || // Require cluster to be closer than dr2_ ...
149  index == cluster_idx[itrk][ipoint])) {
150  continue;
151  } // ... unless it is the closest one ...
152  if (seed.isNull()) {
153  seed = clu;
154  }
155  clusters.push_back(clu);
156  energy += clu->correctedEnergy();
157  X += clu->position().X() * clu->correctedEnergy();
158  Y += clu->position().Y() * clu->correctedEnergy();
159  Z += clu->position().Z() * clu->correctedEnergy();
160  auto index = pfcluster_idx[itrk][ipoint];
161  if (index < static_cast<decltype(index)>(ecalClusters.size())) {
162  barePtrs.push_back(&(ecalClusters[index]));
163  }
164  }
165  X /= energy;
166  Y /= energy;
167  Z /= energy;
168 
169  // Create SC
170  if (seed.isNonnull()) {
171  reco::SuperCluster sc(energy, math::XYZPoint(X, Y, Z));
172  sc.setCorrectedEnergy(energy);
173  sc.setSeed(seed);
174  sc.setClusters(clusters);
175  PFClusterWidthAlgo pfwidth(barePtrs);
176  sc.setEtaWidth(pfwidth.pflowEtaWidth());
177  sc.setPhiWidth(pfwidth.pflowPhiWidth());
178  sc.rawEnergy(); // Cache the value of raw energy
179  superClusters->push_back(sc);
180 
181  // Populate ValueMap container
182  superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd, superClusters->size() - 1));
183  } else {
184  superClustersValueMap.push_back(reco::SuperClusterRef(superClustersRefProd.id()));
185  }
186 
187  } // GSF tracks
188 
189  // Put SuperClusters in event
190  event.put(std::move(superClusters));
191 
192  auto ptr = std::make_unique<edm::ValueMap<reco::SuperClusterRef> >(edm::ValueMap<reco::SuperClusterRef>());
194  filler.insert(gsfPfRecTracks, superClustersValueMap.begin(), superClustersValueMap.end());
195  filler.fill();
196  event.put(std::move(ptr));
197 }
198 
200 //
203  std::vector<int>& matched) {
204  reco::PFClusterRef closest;
205  if (point.isValid()) {
206  float dr2min = dr2_;
207  for (size_t ii = 0; ii < clusters->size(); ++ii) {
208  if (std::find(matched.begin(), matched.end(), ii) == matched.end()) {
209  float dr2 = reco::deltaR2(clusters->at(ii), point.positionREP());
210  if (dr2 < dr2min) {
211  closest = reco::PFClusterRef(clusters, ii);
212  dr2min = dr2;
213  }
214  }
215  }
216  if (dr2min < (dr2_ - 1.e-6)) {
217  matched.push_back(closest.index());
218  }
219  }
220  return closest;
221 }
222 
224 //
227  desc.add<edm::InputTag>("gsfPfRecTracks", edm::InputTag("lowPtGsfElePfGsfTracks"));
228  desc.add<edm::InputTag>("ecalClusters", edm::InputTag("particleFlowClusterECAL"));
229  desc.add<edm::InputTag>("hcalClusters", edm::InputTag("particleFlowClusterHCAL"));
230  desc.add<double>("MaxDeltaR2", 0.5);
231  descriptions.add("lowPtGsfElectronSuperClusters", desc);
232 }
233 
235 //
const REPPoint & positionREP() const
trajectory position in (rho, eta, phi) base
size
Write out results.
T getParameter(std::string const &) const
const edm::EDGetTokenT< reco::PFClusterCollection > ecalClusters_
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:149
double pflowPhiWidth() const
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:53
#define X(str)
Definition: MuonsGrabber.cc:38
const edm::EDGetTokenT< reco::GsfPFRecTrackCollection > gsfPfRecTracks_
key_type index() const
Definition: Ref.h:253
void setSeed(const CaloClusterPtr &r)
list of used xtals by DetId // now inherited by CaloCluster
Definition: SuperCluster.h:107
static void fillDescriptions(edm::ConfigurationDescriptions &)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
void setPhiWidth(double pw)
Definition: SuperCluster.h:73
double pflowEtaWidth() const
void setClusters(const CaloClusterPtrVector &clusters)
Definition: SuperCluster.h:110
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:547
void setEtaWidth(double ew)
Definition: SuperCluster.h:74
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void produce(edm::Event &, const edm::EventSetup &) override
LowPtGsfElectronSCProducer(const edm::ParameterSet &)
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
void setCorrectedEnergy(double cenergy)
Definition: CaloCluster.h:136
edm::Ref< PFClusterCollection > PFClusterRef
persistent reference to PFCluster objects
Definition: PFClusterFwd.h:15
bool isNull() const
Checks for null.
Definition: Ptr.h:142
std::vector< CaloCluster > CaloClusterCollection
collection of CaloCluster objects
reco::PFClusterRef closestCluster(const reco::PFTrajectoryPoint &point, const edm::Handle< reco::PFClusterCollection > &clusters, std::vector< int > &matched)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
ii
Definition: cuy.py:590
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:124
double rawEnergy() const
raw uncorrected energy (sum of energies of component BasicClusters)
Definition: SuperCluster.h:58
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
bool isValid() const
is this point valid ?
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
void add(std::string const &label, ParameterSetDescription const &psetDescription)
std::vector< PFCluster > PFClusterCollection
collection of PFCluster objects
Definition: PFClusterFwd.h:9
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
auto makeValid(const U &iOtherHandleType) noexcept(false)
Definition: ValidHandle.h:52
def move(src, dest)
Definition: eostools.py:511
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
Definition: event.py:1