Core/src/PythiaHepMCFilterGammaGamma.cc

0001 #include "GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h"
0002 #include "DataFormats/Math/interface/LorentzVector.h"
0003 #include "Math/GenVector/VectorUtil.h"
0004 //#include "CLHEP/HepMC/GenParticle.h"
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006
0007 #include "TLorentzVector.h"
0008
0009 #include <iostream>
0010
0011 using namespace edm;
0012 using namespace std;
0013 using namespace HepMC;
0014
0015 PythiaHepMCFilterGammaGamma::PythiaHepMCFilterGammaGamma(const edm::ParameterSet& iConfig)
0016     : ptSeedThr(iConfig.getParameter<double>("PtSeedThr")),
0017       etaSeedThr(iConfig.getParameter<double>("EtaSeedThr")),
0018       ptGammaThr(iConfig.getParameter<double>("PtGammaThr")),
0019       etaGammaThr(iConfig.getParameter<double>("EtaGammaThr")),
0020       ptTkThr(iConfig.getParameter<double>("PtTkThr")),
0021       etaTkThr(iConfig.getParameter<double>("EtaTkThr")),
0022       ptElThr(iConfig.getParameter<double>("PtElThr")),
0023       etaElThr(iConfig.getParameter<double>("EtaElThr")),
0024       dRTkMax(iConfig.getParameter<double>("dRTkMax")),
0025       dRSeedMax(iConfig.getParameter<double>("dRSeedMax")),
0026       dPhiSeedMax(iConfig.getParameter<double>("dPhiSeedMax")),
0027       dEtaSeedMax(iConfig.getParameter<double>("dEtaSeedMax")),
0028       dRNarrowCone(iConfig.getParameter<double>("dRNarrowCone")),
0029       pTMinCandidate1(iConfig.getParameter<double>("PtMinCandidate1")),
0030       pTMinCandidate2(iConfig.getParameter<double>("PtMinCandidate2")),
0031       etaMaxCandidate(iConfig.getParameter<double>("EtaMaxCandidate")),
0032       invMassMin(iConfig.getParameter<double>("InvMassMin")),
0033       invMassMax(iConfig.getParameter<double>("InvMassMax")),
0034       energyCut(iConfig.getParameter<double>("EnergyCut")),
0035       nTkConeMax(iConfig.getParameter<int>("NTkConeMax")),
0036       nTkConeSum(iConfig.getParameter<int>("NTkConeSum")),
0037       acceptPrompts(iConfig.getParameter<bool>("AcceptPrompts")),
0038       promptPtThreshold(iConfig.getParameter<double>("PromptPtThreshold")) {}
0039
0040 PythiaHepMCFilterGammaGamma::~PythiaHepMCFilterGammaGamma() {}
0041
0042 bool PythiaHepMCFilterGammaGamma::filter(const HepMC::GenEvent* myGenEvent) {
0043   bool accepted = false;
0044
0045   // electron/photon seeds
0046   std::vector<const GenParticle*> seeds;
0047
0048   // other electrons/photons to be added to seeds
0049   // to form candidates
0050   std::vector<const GenParticle*> egamma;
0051
0052   // charged tracks to be taken into account in the isolation cones
0053   // around candidates
0054   std::vector<const GenParticle*> stable;
0055
0056   //----------
0057   // 1. find electron/photon seeds
0058   //----------
0059   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0060        ++p) {
0061     if (((*p)->status() == 1 && (*p)->pdg_id() == 22) ||     // gamma
0062         ((*p)->status() == 1 && abs((*p)->pdg_id()) == 11))  // electron
0063
0064     {
0065       // check for eta and pT threshold for seed in gamma, el
0066       if ((*p)->momentum().perp() > ptSeedThr && fabs((*p)->momentum().eta()) < etaSeedThr) {
0067         seeds.push_back(*p);
0068       }
0069     }
0070
0071     if ((*p)->status() == 1) {
0072       // save charged stable tracks
0073       if (abs((*p)->pdg_id()) == 211 ||  // charged pion
0074           abs((*p)->pdg_id()) == 321 ||  // charged kaon
0075           abs((*p)->pdg_id()) == 11 ||   // electron
0076           abs((*p)->pdg_id()) == 13 ||   // muon
0077           abs((*p)->pdg_id()) == 15) {   // tau
0078         // check if it passes the cut
0079         if ((*p)->momentum().perp() > ptTkThr && fabs((*p)->momentum().eta()) < etaTkThr) {
0080           stable.push_back(*p);
0081         }
0082       }
0083
0084       // save egamma for candidate calculation
0085       if ((*p)->pdg_id() == 22 && (*p)->momentum().perp() > ptGammaThr && fabs((*p)->momentum().eta()) < etaGammaThr) {
0086         egamma.push_back(*p);
0087       }
0088       if (abs((*p)->pdg_id()) == 11 && (*p)->momentum().perp() > ptElThr && fabs((*p)->momentum().eta()) < etaElThr) {
0089         egamma.push_back(*p);
0090       }
0091     }
0092   }
0093
0094   if (seeds.size() < 2)
0095     return accepted;
0096
0097   //----------
0098   // 2. loop over seeds to build candidates
0099   //
0100   //    (adding nearby electrons/photons
0101   //     to the seed electrons/photons to obtain the total
0102   //     electromagnetic energy)
0103   //----------
0104
0105   // number of tracks around each of the candidates
0106   std::vector<int> nTracks;
0107
0108   // the candidates (four momenta) formed from the
0109   // seed electrons/photons and nearby electrons/photons
0110   std::vector<TLorentzVector> candidate;
0111
0112   // these are filled but then not used afterwards (could be removed)
0113   std::vector<TLorentzVector> candidateNarrow, candidateSeed;
0114
0115   std::vector<const GenParticle*>::iterator itSeed;
0116
0117   const GenParticle* mom;
0118   int this_id;
0119   int first_different_id;
0120
0121   for (itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
0122     TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
0123
0124     tempseed.SetXYZM((*itSeed)->momentum().px(), (*itSeed)->momentum().py(), (*itSeed)->momentum().pz(), 0);
0125     for (auto itEn = egamma.begin(); itEn != egamma.end(); ++itEn) {
0126       temp1.SetXYZM((*itEn)->momentum().px(), (*itEn)->momentum().py(), (*itEn)->momentum().pz(), 0);
0127
0128       double DR = temp1.DeltaR(tempseed);
0129       double DPhi = temp1.DeltaPhi(tempseed);
0130       double DEta = (*itEn)->momentum().eta() - (*itSeed)->momentum().eta();
0131       if (DPhi < 0)
0132         DPhi = -DPhi;
0133       if (DEta < 0)
0134         DEta = -DEta;
0135
0136       // accept if within cone or within rectangular region around seed
0137       if (DR < dRSeedMax || (DPhi < dPhiSeedMax && DEta < dEtaSeedMax)) {
0138         energy += temp1;
0139       }
0140       if (DR < dRNarrowCone) {
0141         narrowCone += temp1;
0142       }
0143     }
0144
0145     // number of stable charged particles found within dRTkMax
0146     // around candidate
0147     int counter = 0;
0148
0149     if (energy.Et() != 0.) {
0150       if (fabs(energy.Eta()) < etaMaxCandidate) {
0151         temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0);
0152
0153         // count number of stable particles within cone around candidate
0154         for (auto itStable = stable.begin(); itStable != stable.end(); ++itStable) {
0155           temp1.SetXYZM((*itStable)->momentum().px(), (*itStable)->momentum().py(), (*itStable)->momentum().pz(), 0);
0156           double DR = temp1.DeltaR(temp2);
0157           if (DR < dRTkMax)
0158             counter++;
0159         }
0160
0161         if (acceptPrompts) {
0162           if ((*itSeed)->momentum().perp() > promptPtThreshold) {
0163             // check if *itSeed is a prompt particle
0164
0165             bool isPrompt = true;
0166             this_id = (*itSeed)->pdg_id();
0167             mom = (*itSeed);
0168             while (mom->pdg_id() == this_id) {
0169               const GenParticle* mother =
0170                   mom->production_vertex() ? *(mom->production_vertex()->particles_in_const_begin()) : nullptr;
0171
0172               mom = mother;
0173               if (mom == nullptr) {
0174                 break;
0175               }
0176             }
0177
0178             first_different_id = mom->pdg_id();
0179
0180             if (mom->status() == 2 && std::abs(first_different_id) > 100)
0181               isPrompt = false;
0182
0183             // ignore charged particles around prompt particles
0184             if (isPrompt)
0185               counter = 0;
0186           }
0187         }
0188       }
0189     }
0190
0191     candidate.push_back(energy);
0192     candidateSeed.push_back(tempseed);
0193     candidateNarrow.push_back(narrowCone);
0194     nTracks.push_back(counter);
0195   }
0196
0197   if (candidate.size() < 2)
0198     return accepted;
0199
0200   TLorentzVector minvMin, minvMax;
0201
0202   //----------
0203   // 3. perform further checks on candidates
0204   //
0205   //    (energy, charged isolation requirements etc.)
0206   //----------
0207
0208   int i1, i2;
0209   for (unsigned int i = 0; i < candidate.size() - 1; ++i) {
0210     if (candidate[i].Energy() < energyCut)
0211       continue;
0212     if (nTracks[i] > nTkConeMax)
0213       continue;
0214     if (fabs(candidate[i].Eta()) > etaMaxCandidate)
0215       continue;
0216
0217     for (unsigned int j = i + 1; j < candidate.size(); ++j) {
0218       // check features of second candidate alone
0219       if (candidate[j].Energy() < energyCut)
0220         continue;
0221       if (nTracks[j] > nTkConeMax)
0222         continue;
0223       if (fabs(candidate[j].Eta()) > etaMaxCandidate)
0224         continue;
0225
0226       // check requirement on sum of tracks in both isolation cones
0227       if (nTracks[i] + nTracks[j] > nTkConeSum)
0228         continue;
0229
0230       // swap candidates to have pt[i1] >= pt[i2]
0231       if (candidate[i].Pt() > candidate[j].Pt()) {
0232         i1 = i;
0233         i2 = j;
0234       } else {
0235         i1 = j;
0236         i2 = i;
0237       }
0238
0239       // require minimum pt on leading and subleading candidate
0240       if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2)
0241         continue;
0242
0243       // apply requirements on candidate pair mass
0244       minvMin = candidate[i] + candidate[j];
0245       if (minvMin.M() < invMassMin)
0246         continue;
0247
0248       minvMax = candidate[i] + candidate[j];
0249       if (minvMax.M() > invMassMax)
0250         continue;
0251
0252       accepted = true;
0253     }
0254   }
0255
0256   return accepted;
0257 }