File indexing completed on 2024-04-06 12:13:41
0001 #include "GeneratorInterface/GenFilters/plugins/PythiaFilterIsolatedTrack.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003
0004 #include "CLHEP/Random/RandomEngine.h"
0005 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0006
0007 #include <iostream>
0008 #include <list>
0009 #include <vector>
0010 #include <cmath>
0011
0012 std::pair<double, double> PythiaFilterIsolatedTrack::GetEtaPhiAtEcal(
0013 double etaIP, double phiIP, double pT, int charge, double vtxZ) {
0014 double deltaPhi;
0015 double etaEC = 100;
0016 double phiEC = 100;
0017 double Rcurv = pT * 33.3 * 100 / 38;
0018 double theta = 2 * atan(exp(-etaIP));
0019 double zNew;
0020 if (theta > CLHEP::halfpi)
0021 theta = CLHEP::pi - theta;
0022 if (fabs(etaIP) < 1.479) {
0023 deltaPhi = -charge * asin(0.5 * ecRad_ / Rcurv);
0024 double alpha1 = 2 * asin(0.5 * ecRad_ / Rcurv);
0025 double z = ecRad_ / tan(theta);
0026 if (etaIP > 0)
0027 zNew = z * (Rcurv * alpha1) / ecRad_ + vtxZ;
0028 else
0029 zNew = -z * (Rcurv * alpha1) / ecRad_ + vtxZ;
0030 double zAbs = fabs(zNew);
0031 if (zAbs < ecDist_) {
0032 etaEC = -log(tan(0.5 * atan(ecRad_ / zAbs)));
0033 deltaPhi = -charge * asin(0.5 * ecRad_ / Rcurv);
0034 }
0035 if (zAbs > ecDist_) {
0036 zAbs = (fabs(etaIP) / etaIP) * ecDist_;
0037 double Zflight = fabs(zAbs - vtxZ);
0038 double alpha = (Zflight * ecRad_) / (z * Rcurv);
0039 double Rec = 2 * Rcurv * sin(alpha / 2);
0040 deltaPhi = -charge * alpha / 2;
0041 etaEC = -log(tan(0.5 * atan(Rec / ecDist_)));
0042 }
0043 } else {
0044 zNew = (fabs(etaIP) / etaIP) * ecDist_;
0045 double Zflight = fabs(zNew - vtxZ);
0046 double Rvirt = fabs(Zflight * tan(theta));
0047 double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
0048 deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
0049 etaEC = -log(tan(0.5 * atan(Rec / ecDist_)));
0050 }
0051
0052 if (zNew < 0)
0053 etaEC = -etaEC;
0054 phiEC = phiIP + deltaPhi;
0055
0056 if (phiEC < -CLHEP::pi)
0057 phiEC = 2 * CLHEP::pi + phiEC;
0058 if (phiEC > CLHEP::pi)
0059 phiEC = -2 * CLHEP::pi + phiEC;
0060
0061 std::pair<double, double> retVal(etaEC, phiEC);
0062 return retVal;
0063 }
0064
0065 double PythiaFilterIsolatedTrack::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0066 double dR, Rec;
0067 if (fabs(eta1) < 1.479)
0068 Rec = ecRad_;
0069 else
0070 Rec = ecDist_;
0071 double ce1 = cosh(eta1);
0072 double ce2 = cosh(eta2);
0073 double te1 = tanh(eta1);
0074 double te2 = tanh(eta2);
0075
0076 double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
0077 if (z != 0)
0078 dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
0079 else
0080 dR = 999999.;
0081 return dR;
0082 }
0083
0084 PythiaFilterIsolatedTrack::PythiaFilterIsolatedTrack(const edm::ParameterSet &iConfig,
0085 const PythiaFilterIsoTracks::Counters *counters)
0086 : token_(consumes<edm::HepMCProduct>(
0087 iConfig.getUntrackedParameter("moduleLabel", edm::InputTag("generator", "unsmeared")))),
0088 maxSeedEta_(iConfig.getUntrackedParameter<double>("maxSeedEta", 2.3)),
0089 minSeedEta_(iConfig.getUntrackedParameter<double>("minSeedEta", 0.0)),
0090 minSeedMom_(iConfig.getUntrackedParameter<double>("minSeedMom", 20.)),
0091 minIsolTrackMom_(iConfig.getUntrackedParameter<double>("minIsolTrackMom", 2.0)),
0092 isolCone_(iConfig.getUntrackedParameter<double>("isolCone", 40.0)),
0093 onlyHadrons_(iConfig.getUntrackedParameter<bool>("onlyHadrons", true)),
0094 nAll_(0),
0095 nGood_(0),
0096 ecDist_(317.0),
0097 ecRad_(129.0),
0098 pdt_(esConsumes<HepPDT::ParticleDataTable, PDTRecord>()) {
0099 edm::LogVerbatim("PythiaFilter") << "PythiaFilterIsolatedTrack: Eta " << minSeedEta_ << ":" << maxSeedEta_ << " p > "
0100 << minSeedMom_ << " Isolation Cone " << isolCone_ << " with p > " << minIsolTrackMom_
0101 << " OnlyHadron " << onlyHadrons_;
0102 }
0103
0104 PythiaFilterIsolatedTrack::~PythiaFilterIsolatedTrack() {}
0105
0106
0107 bool PythiaFilterIsolatedTrack::filter(edm::Event &iEvent, edm::EventSetup const &iSetup) {
0108 ++nAll_;
0109
0110 auto const &pdt = iSetup.getData(pdt_);
0111
0112 edm::Handle<edm::HepMCProduct> evt;
0113 iEvent.getByToken(token_, evt);
0114
0115 const HepMC::GenEvent *myGenEvent = evt->GetEvent();
0116
0117
0118 std::vector<const HepMC::GenParticle *> chargedParticles;
0119
0120
0121 std::vector<const HepMC::GenParticle *> seeds;
0122
0123
0124 for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0125 iter != myGenEvent->particles_end();
0126 ++iter) {
0127 const HepMC::GenParticle *p = *iter;
0128 if (!(pdt.particle(p->pdg_id())))
0129 continue;
0130 int charge3 = pdt.particle(p->pdg_id())->ID().threeCharge();
0131 int status = p->status();
0132 double momentum = p->momentum().rho();
0133 double abseta = fabs(p->momentum().eta());
0134
0135
0136 if (abs(charge3) == 3 && status == 1 && momentum > minIsolTrackMom_ && abseta < maxSeedEta_ + 0.5) {
0137 chargedParticles.push_back(p);
0138 if (momentum > minSeedMom_ && abseta < maxSeedEta_ && abseta >= minSeedEta_) {
0139 seeds.push_back(p);
0140 }
0141 }
0142 }
0143
0144
0145 unsigned int ntrk(0);
0146 for (std::vector<const HepMC::GenParticle *>::const_iterator it1 = seeds.begin(); it1 != seeds.end(); ++it1) {
0147 const HepMC::GenParticle *p1 = *it1;
0148 if (!(pdt.particle(p1->pdg_id())))
0149 continue;
0150 if (p1->pdg_id() < -100 || p1->pdg_id() > 100 || (!onlyHadrons_)) {
0151 std::pair<double, double> EtaPhi1 = GetEtaPhiAtEcal(p1->momentum().eta(),
0152 p1->momentum().phi(),
0153 p1->momentum().perp(),
0154 (pdt.particle(p1->pdg_id()))->ID().threeCharge() / 3,
0155 0.0);
0156
0157
0158 bool failsIso = false;
0159 for (std::vector<const HepMC::GenParticle *>::const_iterator it2 = chargedParticles.begin();
0160 it2 != chargedParticles.end();
0161 ++it2) {
0162 const HepMC::GenParticle *p2 = *it2;
0163
0164
0165 if (p1 != p2) {
0166 std::pair<double, double> EtaPhi2 = GetEtaPhiAtEcal(p2->momentum().eta(),
0167 p2->momentum().phi(),
0168 p2->momentum().perp(),
0169 (pdt.particle(p2->pdg_id()))->ID().threeCharge() / 3,
0170 0.0);
0171
0172
0173
0174
0175 if (getDistInCM(EtaPhi1.first, EtaPhi1.second, EtaPhi2.first, EtaPhi2.second) < isolCone_) {
0176 failsIso = true;
0177 break;
0178 }
0179 }
0180 }
0181
0182 if (!failsIso)
0183 ++ntrk;
0184 }
0185 }
0186 if (ntrk > 0) {
0187 ++nGood_;
0188 return true;
0189 } else {
0190 return false;
0191 }
0192 }
0193
0194 void PythiaFilterIsolatedTrack::endStream() {
0195 globalCache()->nAll_ += nAll_;
0196 globalCache()->nGood_ += nGood_;
0197 }
0198
0199 void PythiaFilterIsolatedTrack::globalEndJob(const PythiaFilterIsoTracks::Counters *count) {
0200 edm::LogVerbatim("PythiaFilter") << "PythiaFilterIsolatedTrack::Accepts " << count->nGood_ << " events out of "
0201 << count->nAll_ << std::endl;
0202 }