CMS 3D CMS Logo

GenParticlePruner.cc
Go to the documentation of this file.
12 
13 namespace helper {
14  struct SelectCode {
15  enum KeepOrDrop { kKeep, kDrop };
16  enum FlagDepth { kNone, kFirst, kAll };
19  bool all_;
20  };
21 } // namespace helper
22 
24 public:
26 
27 private:
28  void produce(edm::Event &, const edm::EventSetup &) override;
33  std::vector<std::string> selection_;
34  std::vector<std::pair<StringCutObjectSelector<reco::GenParticle>, helper::SelectCode>> select_;
35  std::vector<int> flags_;
36  std::vector<size_t> indices_;
37  void parse(const std::string &selection, helper::SelectCode &code, std::string &cut) const;
38  void flagDaughters(const reco::GenParticle &, int);
39  void flagMothers(const reco::GenParticle &, int);
40  void recursiveFlagDaughters(size_t, const reco::GenParticleCollection &, int, std::vector<size_t> &);
41  void recursiveFlagMothers(size_t, const reco::GenParticleCollection &, int, std::vector<size_t> &);
42  void getDaughterKeys(std::vector<size_t> &, std::vector<size_t> &, const reco::GenParticleRefVector &) const;
43  void getMotherKeys(std::vector<size_t> &, std::vector<size_t> &, const reco::GenParticleRefVector &) const;
44 };
45 
46 using namespace edm;
47 using namespace std;
48 using namespace reco;
49 
50 const int keep = 1, drop = -1;
51 
53  using namespace ::helper;
54  size_t f = selection.find_first_not_of(' ');
55  size_t n = selection.size();
56  string command;
57  char c;
58  for (; (c = selection[f]) != ' ' && f < n; ++f) {
59  command.push_back(c);
60  }
61  if (command[0] == '+') {
62  command.erase(0, 1);
63  if (command[0] == '+') {
64  command.erase(0, 1);
65  code.mothersDepth_ = SelectCode::kAll;
66  } else {
67  code.mothersDepth_ = SelectCode::kFirst;
68  }
69  } else
71 
72  if (command[command.size() - 1] == '+') {
73  command.erase(command.size() - 1);
74  if (command[command.size() - 1] == '+') {
75  command.erase(command.size() - 1);
76  code.daughtersDepth_ = SelectCode::kAll;
77  } else {
78  code.daughtersDepth_ = SelectCode::kFirst;
79  }
80  } else
82 
83  if (command == "keep")
84  code.keepOrDrop_ = SelectCode::kKeep;
85  else if (command == "drop")
86  code.keepOrDrop_ = SelectCode::kDrop;
87  else {
88  throw Exception(errors::Configuration) << "invalid selection command: " << command << "\n" << endl;
89  }
90  for (; f < n; ++f) {
91  if (selection[f] != ' ')
92  break;
93  }
94  cut = string(selection, f);
95  if (cut[0] == '*')
96  cut = string(cut, 0, cut.find_first_of(' '));
97  code.all_ = cut == "*";
98 }
99 
101  : firstEvent_(true),
102  srcToken_(consumes<GenParticleCollection>(cfg.getParameter<InputTag>("src"))),
103  tableToken_(esConsumes()),
104  keepOrDropAll_(drop),
105  selection_(cfg.getParameter<vector<string>>("select")) {
106  using namespace ::helper;
107  produces<GenParticleCollection>();
108  produces<edm::Association<reco::GenParticleCollection>>();
109 }
110 
112  const GenParticleRefVector &daughters = gen.daughterRefVector();
113  for (GenParticleRefVector::const_iterator i = daughters.begin(); i != daughters.end(); ++i)
114  flags_[i->key()] = keepOrDrop;
115 }
116 
118  const GenParticleRefVector &mothers = gen.motherRefVector();
119  for (GenParticleRefVector::const_iterator i = mothers.begin(); i != mothers.end(); ++i)
120  flags_[i->key()] = keepOrDrop;
121 }
122 
125  int keepOrDrop,
126  std::vector<size_t> &allIndices) {
127  GenParticleRefVector daughters = src[index].daughterRefVector();
128  // avoid infinite recursion if the daughters are set to "this" particle.
129  size_t cachedIndex = index;
130  for (GenParticleRefVector::const_iterator i = daughters.begin(); i != daughters.end(); ++i) {
131  index = i->key();
132  // To also avoid infinite recursion if a "loop" is found in the daughter list,
133  // check to make sure the index hasn't already been added.
134  if (find(allIndices.begin(), allIndices.end(), index) == allIndices.end()) {
135  allIndices.push_back(index);
136  if (cachedIndex != index) {
137  flags_[index] = keepOrDrop;
138  recursiveFlagDaughters(index, src, keepOrDrop, allIndices);
139  }
140  }
141  }
142 }
143 
146  int keepOrDrop,
147  std::vector<size_t> &allIndices) {
148  GenParticleRefVector mothers = src[index].motherRefVector();
149  // avoid infinite recursion if the mothers are set to "this" particle.
150  size_t cachedIndex = index;
151  for (GenParticleRefVector::const_iterator i = mothers.begin(); i != mothers.end(); ++i) {
152  index = i->key();
153  // To also avoid infinite recursion if a "loop" is found in the daughter list,
154  // check to make sure the index hasn't already been added.
155  if (find(allIndices.begin(), allIndices.end(), index) == allIndices.end()) {
156  allIndices.push_back(index);
157  if (cachedIndex != index) {
158  flags_[index] = keepOrDrop;
159  recursiveFlagMothers(index, src, keepOrDrop, allIndices);
160  }
161  }
162  }
163 }
164 
166  if (firstEvent_) {
167  auto const &pdt = es.getData(tableToken_);
168 
169  for (vector<string>::const_iterator i = selection_.begin(); i != selection_.end(); ++i) {
170  string cut;
171  ::helper::SelectCode code;
172  parse(*i, code, cut);
173  if (code.all_) {
174  if (i != selection_.begin())
176  << "selections \"keep *\" and \"drop *\" can be used only as first options. Here used in position # "
177  << (i - selection_.begin()) + 1 << "\n"
178  << endl;
179  switch (code.keepOrDrop_) {
180  case ::helper::SelectCode::kDrop:
182  break;
183  case ::helper::SelectCode::kKeep:
185  };
186  } else {
187  cut = pdgEntryReplace(cut, pdt);
188  select_.push_back(make_pair(StringCutObjectSelector<GenParticle>(cut), code));
189  }
190  }
191  firstEvent_ = false;
192  }
193 
194  using namespace ::helper;
196  evt.getByToken(srcToken_, src);
197  const size_t n = src->size();
198  flags_.clear();
199  flags_.resize(n, keepOrDropAll_);
200  for (size_t j = 0; j < select_.size(); ++j) {
201  const pair<StringCutObjectSelector<GenParticle>, SelectCode> &sel = select_[j];
202  SelectCode code = sel.second;
204  for (size_t i = 0; i < n; ++i) {
205  const GenParticle &p = (*src)[i];
206  if (cut(p)) {
207  int keepOrDrop = keep;
208  switch (code.keepOrDrop_) {
209  case SelectCode::kKeep:
210  keepOrDrop = keep;
211  break;
212  case SelectCode::kDrop:
213  keepOrDrop = drop;
214  };
215  flags_[i] = keepOrDrop;
216  std::vector<size_t> allIndicesDa;
217  std::vector<size_t> allIndicesMo;
218  switch (code.daughtersDepth_) {
219  case SelectCode::kAll:
220  recursiveFlagDaughters(i, *src, keepOrDrop, allIndicesDa);
221  break;
222  case SelectCode::kFirst:
223  flagDaughters(p, keepOrDrop);
224  break;
225  case SelectCode::kNone:;
226  };
227  switch (code.mothersDepth_) {
228  case SelectCode::kAll:
229  recursiveFlagMothers(i, *src, keepOrDrop, allIndicesMo);
230  break;
231  case SelectCode::kFirst:
232  flagMothers(p, keepOrDrop);
233  break;
234  case SelectCode::kNone:;
235  };
236  }
237  }
238  }
239  indices_.clear();
240  int counter = 0;
241  for (size_t i = 0; i < n; ++i) {
242  if (flags_[i] == keep) {
243  indices_.push_back(i);
244  flags_[i] = counter++;
245  } else {
246  flags_[i] = -1; //set to invalid ref
247  }
248  }
249 
250  auto out = std::make_unique<GenParticleCollection>();
252  out->reserve(counter);
253 
254  for (vector<size_t>::const_iterator i = indices_.begin(); i != indices_.end(); ++i) {
255  size_t index = *i;
256  const GenParticle &gen = (*src)[index];
257  const LeafCandidate &part = gen;
258  out->push_back(GenParticle(part));
259  GenParticle &newGen = out->back();
260  //fill status flags
261  newGen.statusFlags() = gen.statusFlags();
262  // The "daIndxs" and "moIndxs" keep a list of the keys for the mother/daughter
263  // parentage/descendency. In some cases, a circular referencing is encountered,
264  // which would result in an infinite loop. The list is checked to
265  // avoid this.
266  vector<size_t> daIndxs, daNewIndxs;
267  getDaughterKeys(daIndxs, daNewIndxs, gen.daughterRefVector());
268  std::sort(daNewIndxs.begin(), daNewIndxs.end());
269  for (size_t i = 0; i < daNewIndxs.size(); ++i)
270  newGen.addDaughter(GenParticleRef(outRef, daNewIndxs[i]));
271 
272  vector<size_t> moIndxs, moNewIndxs;
273  getMotherKeys(moIndxs, moNewIndxs, gen.motherRefVector());
274  std::sort(moNewIndxs.begin(), moNewIndxs.end());
275  for (size_t i = 0; i < moNewIndxs.size(); ++i)
276  newGen.addMother(GenParticleRef(outRef, moNewIndxs[i]));
277  }
278 
280  auto orig2new = std::make_unique<edm::Association<reco::GenParticleCollection>>(oh);
282  orig2newFiller.insert(src, flags_.begin(), flags_.end());
283  orig2newFiller.fill();
284  evt.put(std::move(orig2new));
285 }
286 
287 void GenParticlePruner::getDaughterKeys(vector<size_t> &daIndxs,
288  vector<size_t> &daNewIndxs,
289  const GenParticleRefVector &daughters) const {
290  for (GenParticleRefVector::const_iterator j = daughters.begin(); j != daughters.end(); ++j) {
291  GenParticleRef dau = *j;
292  if (find(daIndxs.begin(), daIndxs.end(), dau.key()) == daIndxs.end()) {
293  daIndxs.push_back(dau.key());
294  int idx = flags_[dau.key()];
295  if (idx > 0) {
296  daNewIndxs.push_back(idx);
297  } else {
298  const GenParticleRefVector &daus = dau->daughterRefVector();
299  if (!daus.empty())
300  getDaughterKeys(daIndxs, daNewIndxs, daus);
301  }
302  }
303  }
304 }
305 
306 void GenParticlePruner::getMotherKeys(vector<size_t> &moIndxs,
307  vector<size_t> &moNewIndxs,
308  const GenParticleRefVector &mothers) const {
309  for (GenParticleRefVector::const_iterator j = mothers.begin(); j != mothers.end(); ++j) {
310  GenParticleRef mom = *j;
311  if (find(moIndxs.begin(), moIndxs.end(), mom.key()) == moIndxs.end()) {
312  moIndxs.push_back(mom.key());
313  int idx = flags_[mom.key()];
314  if (idx >= 0) {
315  moNewIndxs.push_back(idx);
316  } else {
317  const GenParticleRefVector &moms = mom->motherRefVector();
318  if (!moms.empty())
319  getMotherKeys(moIndxs, moNewIndxs, moms);
320  }
321  }
322  }
323 }
324 
326 
edm::RefProd< GenParticleCollection >
counter
Definition: counter.py:1
edm::Event::getRefBeforePut
RefProd< PROD > getRefBeforePut()
Definition: Event.h:158
GenParticlePruner::selection_
std::vector< std::string > selection_
Definition: GenParticlePruner.cc:33
Handle.h
mps_fire.i
i
Definition: mps_fire.py:428
GenParticlePruner::flagDaughters
void flagDaughters(const reco::GenParticle &, int)
Definition: GenParticlePruner.cc:111
edm::Association::Filler
Definition: Association.h:78
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
helper::SelectCode::keepOrDrop_
KeepOrDrop keepOrDrop_
Definition: GenParticlePruner.cc:17
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
reco::GenParticle
Definition: GenParticle.h:21
edm::EDGetTokenT< reco::GenParticleCollection >
edm
HLT enums.
Definition: AlignableModifier.h:19
GenParticlePruner::select_
std::vector< std::pair< StringCutObjectSelector< reco::GenParticle >, helper::SelectCode > > select_
Definition: GenParticlePruner.cc:34
edm::RefVector::begin
const_iterator begin() const
Initialize an iterator over the RefVector.
Definition: RefVector.h:223
reco::CompositeRefCandidateT::addMother
void addMother(const typename mothers::value_type &)
add a daughter via a reference
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:13
GenParticlePruner::indices_
std::vector< size_t > indices_
Definition: GenParticlePruner.cc:36
GenParticlePruner::flags_
std::vector< int > flags_
Definition: GenParticlePruner.cc:35
GenParticlePruner::GenParticlePruner
GenParticlePruner(const edm::ParameterSet &)
Definition: GenParticlePruner.cc:100
EDProducer.h
GenParticlePruner::parse
void parse(const std::string &selection, helper::SelectCode &code, std::string &cut) const
Definition: GenParticlePruner.cc:52
GenParticlePruner::flagMothers
void flagMothers(const reco::GenParticle &, int)
Definition: GenParticlePruner.cc:117
edm::RefVector< GenParticleCollection >
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
Association.h
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
mps_check.command
list command
Definition: mps_check.py:25
GenParticlePruner::recursiveFlagMothers
void recursiveFlagMothers(size_t, const reco::GenParticleCollection &, int, std::vector< size_t > &)
Definition: GenParticlePruner.cc:144
edm::Handle
Definition: AssociativeIterator.h:50
GenParticle
Definition: GenParticle.py:1
edm::Ref< GenParticleCollection >
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
GenParticle.h
EDMException.h
MakerMacros.h
edm::RefVector::end
const_iterator end() const
Termination of iteration.
Definition: RefVector.h:228
part
part
Definition: HCALResponse.h:20
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::RefVector::empty
bool empty() const
Is the RefVector empty.
Definition: RefVector.h:99
BPhysicsValidation_cfi.daughters
daughters
Definition: BPhysicsValidation_cfi.py:11
helper::SelectCode::all_
bool all_
Definition: GenParticlePruner.cc:19
corrVsCorr.selection
selection
main part
Definition: corrVsCorr.py:100
helper::SelectCode::kKeep
Definition: GenParticlePruner.cc:15
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
GenParticlePruner::recursiveFlagDaughters
void recursiveFlagDaughters(size_t, const reco::GenParticleCollection &, int, std::vector< size_t > &)
Definition: GenParticlePruner.cc:123
gen
Definition: PythiaDecays.h:13
GenParticlePruner::keepOrDropAll_
int keepOrDropAll_
Definition: GenParticlePruner.cc:32
funct::true
true
Definition: Factorize.h:173
helper::SelectCode::FlagDepth
FlagDepth
Definition: GenParticlePruner.cc:16
edm::ParameterSet
Definition: ParameterSet.h:47
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
TrackRefitter_38T_cff.src
src
Definition: TrackRefitter_38T_cff.py:24
Event.h
ParticleDataTable.h
pdgEntryReplace.h
GenParticlePruner::srcToken_
edm::EDGetTokenT< reco::GenParticleCollection > srcToken_
Definition: GenParticlePruner.cc:30
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
reco::GenParticle::statusFlags
const GenStatusFlags & statusFlags() const
Definition: GenParticle.h:38
helper::SelectCode::kFirst
Definition: GenParticlePruner.cc:16
helper::SelectCode
Definition: GenParticlePruner.cc:14
helper
Definition: helper.py:1
GenParticlePruner::firstEvent_
bool firstEvent_
Definition: GenParticlePruner.cc:29
edm::Event::put
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:133
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
edm::stream::EDProducer
Definition: EDProducer.h:36
pdgEntryReplace
std::string pdgEntryReplace(const std::string &, HepPDT::ParticleDataTable const &)
edm::EventSetup
Definition: EventSetup.h:58
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord >
InputTag.h
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:127
looper.cfg
cfg
Definition: looper.py:296
WDecay::kNone
Definition: TopGenEvent.h:27
helper::SelectCode::kNone
Definition: GenParticlePruner.cc:16
GenParticle.GenParticle
GenParticle
Definition: GenParticle.py:18
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
edm::OrphanHandle
Definition: EDProductfwd.h:39
StringCutObjectSelector.h
helper::SelectCode::daughtersDepth_
FlagDepth daughtersDepth_
Definition: GenParticlePruner.cc:18
PA_MinBiasSkim_cff.cut
cut
Definition: PA_MinBiasSkim_cff.py:13
GenParticlePruner::getMotherKeys
void getMotherKeys(std::vector< size_t > &, std::vector< size_t > &, const reco::GenParticleRefVector &) const
Definition: GenParticlePruner.cc:306
StringCutObjectSelector
Definition: StringCutObjectSelector.h:16
Exception
Definition: hltDiff.cc:245
edm::RefVectorIterator
Definition: EDProductfwd.h:33
GenParticlePruner::getDaughterKeys
void getDaughterKeys(std::vector< size_t > &, std::vector< size_t > &, const reco::GenParticleRefVector &) const
Definition: GenParticlePruner.cc:287
GenParticlePruner::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: GenParticlePruner.cc:165
relval_steps.gen
def gen(fragment, howMuch)
Production test section ####.
Definition: relval_steps.py:506
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
helper::SelectCode::mothersDepth_
FlagDepth mothersDepth_
Definition: GenParticlePruner.cc:18
edm::Ref::key
key_type key() const
Accessor for product key.
Definition: Ref.h:250
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
reco::CompositeRefCandidateT::addDaughter
void addDaughter(const typename daughters::value_type &)
add a daughter via a reference
drop
const int drop
Definition: GenParticlePruner.cc:50
helper::SelectCode::KeepOrDrop
KeepOrDrop
Definition: GenParticlePruner.cc:15
reco::LeafCandidate
Definition: LeafCandidate.h:16
GenParticlePruner
Definition: GenParticlePruner.cc:23
helper::SelectCode::kAll
Definition: GenParticlePruner.cc:16
ParameterSet.h
EgammaValidation_Wenu_cff.sel
sel
Definition: EgammaValidation_Wenu_cff.py:33
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:56
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
DeDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::Event
Definition: Event.h:73
GenParticlePruner::tableToken_
edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord > tableToken_
Definition: GenParticlePruner.cc:31
keep
const int keep
Definition: GenParticlePruner.cc:50
edm::errors::Configuration
Definition: EDMException.h:36
edm::InputTag
Definition: InputTag.h:15
helper::SelectCode::kDrop
Definition: GenParticlePruner.cc:15