GenFilters/plugins/PythiaFilterIsolatedTrack.cc

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;  //r(m)=pT(GeV)*33.3/B(kG)
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;  //new z-coordinate of track
0028     else
0029       zNew = -z * (Rcurv * alpha1) / ecRad_ + vtxZ;  //new z-coordinate of track
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 // ------------ method called to produce the data  ------------
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   // all of the stable, charged particles with momentum>minIsolTrackMom_ and |eta|<maxSeedEta_+0.5
0118   std::vector<const HepMC::GenParticle *> chargedParticles;
0119
0120   // all of the stable, charged particles with momentum>minSeedMom_ and minSeedEta_<=|eta|<maxSeedEta_
0121   std::vector<const HepMC::GenParticle *> seeds;
0122
0123   // fill the vector of charged particles and seeds in the event
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     // only consider stable, charged particles
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   // loop over all the seeds and see if any of them are isolated
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_)) {  // Select hadrons only
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       // loop over all of the other charged particles in the event, and see if any are close by
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         // don't consider the seed particle among the other charge particles
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           // find out how far apart the particles are
0173           // if the seed fails the isolation requirement, try a different seed
0174           // occasionally allow a seed to pass to isolation requirement
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   }  //loop over seeds
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 }