File indexing completed on 2022-04-15 00:38:11
0001 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include <iostream>
0006 #include <limits>
0007 #include "TLorentzVector.h"
0008
0009
0010
0011 namespace HepMCValidationHelper {
0012 void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
0013 const HepMC::GenEvent* all,
0014 double deltaRcut,
0015 std::vector<const HepMC::GenParticle*>& fsrphotons) {
0016 std::vector<const HepMC::GenParticle*> status1;
0017 allStatus1(all, status1);
0018 findFSRPhotons(leptons, status1, deltaRcut, fsrphotons);
0019 }
0020
0021 void findFSRPhotons(const std::vector<const HepMC::GenParticle*>& leptons,
0022 const std::vector<const HepMC::GenParticle*>& all,
0023 double deltaRcut,
0024 std::vector<const HepMC::GenParticle*>& fsrphotons) {
0025
0026 std::vector<const HepMC::GenParticle*> allphotons;
0027 for (unsigned int i = 0; i < all.size(); ++i) {
0028 if (all[i]->status() == 1 && all[i]->pdg_id() == 22)
0029 allphotons.push_back(all[i]);
0030 }
0031
0032
0033 for (unsigned int ipho = 0; ipho < allphotons.size(); ++ipho) {
0034 bool close = false;
0035 for (unsigned int ilep = 0; ilep < leptons.size(); ++ilep) {
0036 if (deltaR(allphotons[ipho]->momentum(), leptons[ilep]->momentum()) < deltaRcut) {
0037 close = true;
0038 break;
0039 }
0040 }
0041 if (close)
0042 fsrphotons.push_back(allphotons[ipho]);
0043 }
0044 }
0045
0046
0047 bool isChargedLepton(const HepMC::GenParticle* part) {
0048 int status = part->status();
0049 unsigned int pdg_id = abs(part->pdg_id());
0050 if (status == 2)
0051 return pdg_id == 15;
0052 else
0053 return status == 1 && (pdg_id == 11 || pdg_id == 13);
0054 }
0055
0056
0057 bool isNeutrino(const HepMC::GenParticle* part) {
0058 int status = part->status();
0059 unsigned int pdg_id = abs(part->pdg_id());
0060 return status == 1 && (pdg_id == 12 || pdg_id == 14 || pdg_id == 16);
0061 }
0062
0063
0064 bool isTau(const HepMC::GenParticle* part) { return part->status() == 2 && abs(part->pdg_id()) == 15; }
0065
0066
0067
0068
0069
0070
0071
0072
0073 void allStatus1(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0074 for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0075 if ((*iter)->status() == 1)
0076 status1.push_back(*iter);
0077 }
0078 }
0079
0080 void allStatus2(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0081 for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0082 if ((*iter)->status() == 2)
0083 status1.push_back(*iter);
0084 }
0085 }
0086
0087 void allStatus3(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& status1) {
0088 for (HepMC::GenEvent::particle_const_iterator iter = all->particles_begin(); iter != all->particles_end(); ++iter) {
0089 if ((*iter)->status() == 3)
0090 status1.push_back(*iter);
0091 }
0092 }
0093
0094 void findDescendents(const HepMC::GenParticle* a, std::vector<const HepMC::GenParticle*>& descendents) {
0095 HepMC::GenVertex* decayVertex = a->end_vertex();
0096 if (!decayVertex)
0097 return;
0098 HepMC::GenVertex::particles_out_const_iterator ipart;
0099 for (ipart = decayVertex->particles_out_const_begin(); ipart != decayVertex->particles_out_const_end(); ++ipart) {
0100 if ((*ipart)->status() == 1)
0101 descendents.push_back(*ipart);
0102 else
0103 findDescendents(*ipart, descendents);
0104 }
0105 }
0106
0107 void removeIsolatedLeptons(const HepMC::GenEvent* all,
0108 double deltaRcut,
0109 double sumPtCut,
0110 std::vector<const HepMC::GenParticle*>& pruned) {
0111
0112 std::vector<const HepMC::GenParticle*> status1;
0113 allStatus1(all, status1);
0114 std::vector<const HepMC::GenParticle*> toRemove;
0115
0116 for (unsigned int i = 0; i < status1.size(); ++i) {
0117
0118 if (isNeutrino(status1[i]) || (isChargedLepton(status1[i]) && !isTau(status1[i]))) {
0119
0120
0121 std::vector<const HepMC::GenParticle*> leptons;
0122 leptons.push_back(status1[i]);
0123 std::vector<const HepMC::GenParticle*> removedForIsolation;
0124 removedForIsolation.push_back(status1[i]);
0125 if (isChargedLepton(status1[i]))
0126 findFSRPhotons(leptons, status1, deltaRcut, removedForIsolation);
0127 #ifdef DEBUG_HepMCValidationHelper
0128
0129 #endif
0130
0131 std::vector<const HepMC::GenParticle*> forIsolation;
0132 std::vector<const HepMC::GenParticle*>::iterator iiso;
0133 for (iiso = status1.begin(); iiso != status1.end(); ++iiso) {
0134 std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0135 bool marked = false;
0136 for (iremove = removedForIsolation.begin(); iremove != removedForIsolation.end(); ++iremove) {
0137 if ((*iiso)->barcode() == (*iremove)->barcode()) {
0138 #ifdef DEBUG_HepMCValidationHelper
0139
0140 #endif
0141 marked = true;
0142 break;
0143 }
0144 }
0145 if (!marked)
0146 forIsolation.push_back(*iiso);
0147 }
0148
0149 double sumIso = 0;
0150 for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso) {
0151 if (deltaR(leptons.front()->momentum(), (*iiso)->momentum()) < deltaRcut) {
0152 sumIso += (*iiso)->momentum().perp();
0153 }
0154 }
0155
0156 if (sumIso < sumPtCut) {
0157 #ifdef DEBUG_HepMCValidationHelper
0158 std::cout << "particle " << *status1[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
0159 #endif
0160 toRemove.insert(toRemove.end(), removedForIsolation.begin(), removedForIsolation.end());
0161 }
0162 #ifdef DEBUG_HepMCValidationHelper
0163 else {
0164 std::cout << "NOT isolated! " << *status1[i] << " is considered not isolated, with sumPt " << sumIso
0165 << std::endl;
0166 }
0167 #endif
0168 }
0169 }
0170
0171
0172 std::vector<const HepMC::GenParticle*> status2;
0173 allStatus2(all, status2);
0174 std::vector<const HepMC::GenParticle*> taus;
0175
0176 for (unsigned int i = 0; i < status2.size(); ++i) {
0177 if (isTau(status2[i])) {
0178
0179
0180 bool duplicate = false;
0181 TLorentzVector taumomentum(status2[i]->momentum().x(),
0182 status2[i]->momentum().y(),
0183 status2[i]->momentum().z(),
0184 status2[i]->momentum().t());
0185 for (unsigned int j = 0; j < taus.size(); ++j) {
0186
0187 TLorentzVector othermomentum(
0188 taus[j]->momentum().x(), taus[j]->momentum().y(), taus[j]->momentum().z(), taus[j]->momentum().t());
0189 othermomentum -= taumomentum;
0190 if (status2[i]->pdg_id() == taus[j]->pdg_id() &&
0191 othermomentum.E() < 0.1 &&
0192 othermomentum.P() < 0.1) {
0193 duplicate = true;
0194 break;
0195 }
0196 }
0197 if (!duplicate)
0198 taus.push_back(status2[i]);
0199 }
0200 }
0201
0202 for (unsigned int i = 0; i < taus.size(); ++i) {
0203 std::vector<const HepMC::GenParticle*> taudaughters;
0204 findDescendents(taus[i], taudaughters);
0205 if (taudaughters.empty()) {
0206 throw cms::Exception("NoTauDaugters")
0207 << " HepMCValidationHelper found no daughters for Tau within index " << i << " and info \n"
0208 << *taus[i] << "\n This should not be able to happen and needs to be fixed.";
0209 }
0210 const HepMC::FourVector& taumom = taus[i]->momentum();
0211
0212 std::vector<const HepMC::GenParticle*> forIsolation;
0213 std::vector<const HepMC::GenParticle*>::iterator iiso;
0214 for (iiso = status1.begin(); iiso < status1.end(); ++iiso) {
0215 bool marked = false;
0216 std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0217 for (iremove = taudaughters.begin(); iremove != taudaughters.end(); ++iremove) {
0218 if ((*iiso)->barcode() == (*iremove)->barcode()) {
0219 #ifdef DEBUG_HepMCValidationHelper
0220
0221 #endif
0222 marked = true;
0223 break;
0224 }
0225 }
0226 if (!marked)
0227 forIsolation.push_back(*iiso);
0228 }
0229
0230 double sumIso = 0;
0231 for (iiso = forIsolation.begin(); iiso < forIsolation.end(); ++iiso) {
0232 if (deltaR(taumom, (*iiso)->momentum()) < deltaRcut) {
0233 sumIso += (*iiso)->momentum().perp();
0234 }
0235 }
0236
0237 if (sumIso < sumPtCut) {
0238 #ifdef DEBUG_HepMCValidationHelper
0239 std::cout << "particle " << *taus[i] << " is considered isolated, with sumPt " << sumIso << std::endl;
0240 #endif
0241 toRemove.insert(toRemove.end(), taudaughters.begin(), taudaughters.end());
0242 }
0243 }
0244
0245
0246 pruned.clear();
0247 for (unsigned int i = 0; i < status1.size(); ++i) {
0248 bool marked = false;
0249 std::vector<const HepMC::GenParticle*>::const_iterator iremove;
0250 for (iremove = toRemove.begin(); iremove != toRemove.end(); ++iremove) {
0251 if (status1[i]->barcode() == (*iremove)->barcode()) {
0252 marked = true;
0253 break;
0254 }
0255 }
0256 if (!marked)
0257 pruned.push_back(status1[i]);
0258 }
0259
0260 #ifdef DEBUG_HepMCValidationHelper
0261 std::cout << "list of remaining particles:" << std::endl;
0262 for (unsigned int i = 0; i < pruned.size(); ++i) {
0263 std::cout << *pruned[i] << std::endl;
0264 }
0265 #endif
0266 }
0267
0268
0269 void allVisibleParticles(const HepMC::GenEvent* all, std::vector<const HepMC::GenParticle*>& visible) {
0270 std::vector<const HepMC::GenParticle*> status1;
0271 visible.clear();
0272 allStatus1(all, status1);
0273 for (unsigned int i = 0; i < status1.size(); ++i) {
0274 if (!isNeutrino(status1[i]))
0275 visible.push_back(status1[i]);
0276 }
0277 }
0278
0279
0280 TLorentzVector genMet(const HepMC::GenEvent* all, double etamin, double etamax) {
0281 std::vector<const HepMC::GenParticle*> visible;
0282 allVisibleParticles(all, visible);
0283 TLorentzVector momsum(0., 0., 0., 0.);
0284 for (unsigned int i = 0; i < visible.size(); ++i) {
0285 if (visible[i]->momentum().eta() > etamin && visible[i]->momentum().eta() < etamax) {
0286 TLorentzVector mom(visible[i]->momentum().x(),
0287 visible[i]->momentum().y(),
0288 visible[i]->momentum().z(),
0289 visible[i]->momentum().t());
0290 momsum += mom;
0291 }
0292 }
0293 TLorentzVector met(-momsum.Px(), -momsum.Py(), 0., momsum.Pt());
0294 return met;
0295 }
0296 }