double EM enrichment filter support for HepMCfilter

Configuration/Generator/python/doubleEMEnrichingHepMCfilter_cfi.py

@@ -0,0 +1,40 @@
+import FWCore.ParameterSet.Config as cms
+
+doubleEMenrichingHepMCfilterParams = cms.PSet(
+    # seed thresholds
+    PtSeedThr         = cms.untracked.double(5.0),
+    EtaSeedThr        = cms.untracked.double(2.8),
+
+    # photon thresholds
+    PtGammaThr        = cms.untracked.double(0.0),
+    EtaGammaThr       = cms.untracked.double(2.8),
+
+    # electron threshold
+    PtElThr           = cms.untracked.double(2.0),
+    EtaElThr          = cms.untracked.double(2.8),
+
+    dRSeedMax         = cms.untracked.double(0.0),
+    dPhiSeedMax       = cms.untracked.double(0.2),
+    dEtaSeedMax       = cms.untracked.double(0.12),
+    dRNarrowCone      = cms.untracked.double(0.02),
+
+    PtTkThr           = cms.untracked.double(1.6),
+    EtaTkThr          = cms.untracked.double(2.2),
+    dRTkMax           = cms.untracked.double(0.2),
+
+    PtMinCandidate1   = cms.untracked.double(15.0),
+    PtMinCandidate2   = cms.untracked.double(15.0),
+    EtaMaxCandidate   = cms.untracked.double(3.0),
+
+    NTkConeMax        = cms.untracked.int32(2),
+    NTkConeSum        = cms.untracked.int32(4),
+
+    # mass 80..Inf GeV
+    InvMassMin        = cms.untracked.double(80.0),
+    InvMassMax        = cms.untracked.double(14000.0),
+
+    EnergyCut         = cms.untracked.double(1.0),
+
+    AcceptPrompts     = cms.untracked.bool(True),
+    PromptPtThreshold = cms.untracked.double(15.0),
+    )

GeneratorInterface/Core/BuildFile.xml

@@ -9,6 +9,7 @@
 <use   name="clhep"/>
 <use   name="lhapdf"/>
 <use   name="f77compiler"/>
+<use   name="root"/>
 <export>
   <lib   name="1"/>
 </export>

GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h

@@ -0,0 +1,46 @@
+#ifndef PYTHIAHEPMCFILTERGAMMAGAMMA_h
+#define PYTHIAHEPMCFILTERGAMMAGAMMA_h
+
+//
+// Package:    GeneratorInterface/GenFilters
+// Class:      PythiaHepMCFilterGammaGamma
+// 
+// Original Author:  Matteo Sani
+//
+//
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "GeneratorInterface/Core/interface/BaseHepMCFilter.h"
+#include "FWCore/Framework/interface/EDFilter.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+namespace edm {
+	  class HepMCProduct;
+}
+
+class PythiaHepMCFilterGammaGamma : public BaseHepMCFilter {
+ public:
+  explicit PythiaHepMCFilterGammaGamma(const edm::ParameterSet&);
+  ~PythiaHepMCFilterGammaGamma() override;
+
+  /** @return true if this GenEvent passes the double EM enrichment
+      criterion */
+  bool filter(const HepMC::GenEvent* myGenEvent) override;
+ private:
+
+  const edm::EDGetTokenT<edm::HepMCProduct> token_;
+  const int maxEvents;
+
+  const double ptSeedThr, etaSeedThr, ptGammaThr, etaGammaThr, ptTkThr, etaTkThr;
+  const double ptElThr, etaElThr, dRTkMax, dRSeedMax, dPhiSeedMax, dEtaSeedMax, dRNarrowCone, pTMinCandidate1, pTMinCandidate2, etaMaxCandidate;
+  const double invMassMin, invMassMax;
+  const double energyCut;
+  const int nTkConeMax, nTkConeSum;
+  const bool acceptPrompts;
+  const double promptPtThreshold;
+
+};
+#endif

GeneratorInterface/Core/src/HepMCFilterDriver.cc

@@ -2,6 +2,7 @@
 #include "GeneratorInterface/Core/interface/GenericDauHepMCFilter.h"
 #include "GeneratorInterface/Core/interface/PartonShowerBsHepMCFilter.h"
 #include "GeneratorInterface/Core/interface/EmbeddingHepMCFilter.h"
+#include "GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h"
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
@@ -29,8 +30,11 @@ HepMCFilterDriver::HepMCFilterDriver(const edm::ParameterSet& pset) :
   else if (filterName=="EmbeddingHepMCFilter") {
     filter_ = new EmbeddingHepMCFilter(filterParameters);
   }
+  else if (filterName=="PythiaHepMCFilterGammaGamma") {
+    filter_ = new PythiaHepMCFilterGammaGamma(filterParameters);
+  }
   else {
-    edm::LogError("HepMCFilterDriver")<< "Invalid HepMCFilter name:" << filterName;
+    throw edm::Exception(edm::errors::Configuration,"HepMCFilterDriver") << "Invalid HepMCFilter name:" << filterName;
   }
 
 }

GeneratorInterface/Core/src/PythiaHepMCFilterGammaGamma.cc

@@ -0,0 +1,225 @@
+#include "GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+#include "Math/GenVector/VectorUtil.h"
+//#include "CLHEP/HepMC/GenParticle.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+
+#include "TLorentzVector.h"
+
+#include <iostream>
+
+using namespace edm;
+using namespace std;
+using namespace HepMC;
+
+
+PythiaHepMCFilterGammaGamma::PythiaHepMCFilterGammaGamma(const edm::ParameterSet& iConfig) :
+  maxEvents(iConfig.getUntrackedParameter<int>("maxEvents", 0)),
+  ptSeedThr(iConfig.getUntrackedParameter<double>("PtSeedThr")),
+  etaSeedThr(iConfig.getUntrackedParameter<double>("EtaSeedThr")),
+  ptGammaThr(iConfig.getUntrackedParameter<double>("PtGammaThr")),
+  etaGammaThr(iConfig.getUntrackedParameter<double>("EtaGammaThr")),
+  ptTkThr(iConfig.getUntrackedParameter<double>("PtTkThr")),
+  etaTkThr(iConfig.getUntrackedParameter<double>("EtaTkThr")),
+  ptElThr(iConfig.getUntrackedParameter<double>("PtElThr")),
+  etaElThr(iConfig.getUntrackedParameter<double>("EtaElThr")),
+  dRTkMax(iConfig.getUntrackedParameter<double>("dRTkMax")),
+  dRSeedMax(iConfig.getUntrackedParameter<double>("dRSeedMax")),
+  dPhiSeedMax(iConfig.getUntrackedParameter<double>("dPhiSeedMax")),
+  dEtaSeedMax(iConfig.getUntrackedParameter<double>("dEtaSeedMax")),
+  dRNarrowCone(iConfig.getUntrackedParameter<double>("dRNarrowCone")),
+  pTMinCandidate1(iConfig.getUntrackedParameter<double>("PtMinCandidate1")),
+  pTMinCandidate2(iConfig.getUntrackedParameter<double>("PtMinCandidate2")),
+  etaMaxCandidate(iConfig.getUntrackedParameter<double>("EtaMaxCandidate")),
+  invMassMin(iConfig.getUntrackedParameter<double>("InvMassMin")),
+  invMassMax(iConfig.getUntrackedParameter<double>("InvMassMax")),
+  energyCut(iConfig.getUntrackedParameter<double>("EnergyCut")),
+  nTkConeMax(iConfig.getUntrackedParameter<int>("NTkConeMax")),
+  nTkConeSum(iConfig.getUntrackedParameter<int>("NTkConeSum")),
+  acceptPrompts(iConfig.getUntrackedParameter<bool>("AcceptPrompts")), 
+  promptPtThreshold(iConfig.getUntrackedParameter<double>("PromptPtThreshold")) {
+
+  if (maxEvents != 0) edm::LogInfo("PythiaFilterGammaGamma::PythiaFilterGammaGamma") << "WARNING, ignoring unsuported option, maxEvents = " << maxEvents << endl;
+
+}
+
+PythiaHepMCFilterGammaGamma::~PythiaHepMCFilterGammaGamma() 
+{  
+}
+
+bool PythiaHepMCFilterGammaGamma::filter(const HepMC::GenEvent* myGenEvent) {
+
+  bool accepted = false;
+
+  std::vector<const GenParticle*> seeds, egamma, stable; 
+  std::vector<const GenParticle*>::const_iterator itPart, itStable, itEn;
+
+ // Loop on egamma
+  for(HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
+
+    if (
+	((*p)->status()==1&&(*p)->pdg_id() == 22) ||                   // gamma
+        ((*p)->status()==1&&abs((*p)->pdg_id()) == 11))                // electron
+
+      {       
+	// check for eta and pT threshold for seed in gamma, el
+	if ((*p)->momentum().perp() > ptSeedThr &&
+	    fabs((*p)->momentum().eta()) < etaSeedThr) {
+
+	  seeds.push_back(*p);
+	}
+	}   
+
+
+    if ((*p)->status() == 1) {
+
+      // save charged stable tracks
+      if (abs((*p)->pdg_id()) == 211 ||
+          abs((*p)->pdg_id()) == 321 ||
+          abs((*p)->pdg_id()) == 11 ||
+          abs((*p)->pdg_id()) == 13 ||
+          abs((*p)->pdg_id()) == 15) {
+        // check if it passes the cut
+        if ((*p)->momentum().perp() > ptTkThr &&
+            fabs((*p)->momentum().eta()) < etaTkThr) {
+          stable.push_back(*p);
+        }
+      }
+
+      // save egamma for candidate calculation
+      if ((*p)->pdg_id() == 22 &&
+	  (*p)->momentum().perp() > ptGammaThr &&
+	  fabs((*p)->momentum().eta()) < etaGammaThr) {
+	egamma.push_back(*p);
+      }
+      if (abs((*p)->pdg_id()) == 11 &&
+	  (*p)->momentum().perp() > ptElThr &&
+	  fabs((*p)->momentum().eta()) < etaElThr) {
+	egamma.push_back(*p); 
+      }
+    }
+  }
+
+  if (seeds.size() < 2) return accepted;
+
+  std::vector<int> nTracks;
+  std::vector<TLorentzVector> candidate, candidateNarrow, candidateSeed;
+  std::vector<const GenParticle*>::iterator itSeed;
+
+  const GenParticle* mom;
+  int this_id;
+  int first_different_id;
+
+  for(itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
+
+    TLorentzVector energy, narrowCone, temp1, temp2, tempseed;
+
+    tempseed.SetXYZM((*itSeed)->momentum().px(), (*itSeed)->momentum().py(), (*itSeed)->momentum().pz(), 0);   
+    for(itEn = egamma.begin(); itEn != egamma.end(); ++itEn) {
+      temp1.SetXYZM((*itEn)->momentum().px(), (*itEn)->momentum().py(), (*itEn)->momentum().pz(), 0);        
+
+      double DR = temp1.DeltaR(tempseed);
+      double DPhi = temp1.DeltaPhi(tempseed);
+      double DEta = (*itEn)->momentum().eta()-(*itSeed)->momentum().eta();
+      if(DPhi<0) DPhi=-DPhi;
+      if(DEta<0) DEta=-DEta;
+
+      if (DR < dRSeedMax || (DPhi<dPhiSeedMax&&DEta<dEtaSeedMax)) {
+	energy += temp1;
+      }
+      if (DR < dRNarrowCone) {
+	narrowCone += temp1;
+      }
+    }
+
+    int counter = 0;
+
+    if ( energy.Et() != 0. ) {
+      if (fabs(energy.Eta()) < etaMaxCandidate) {
+
+	temp2.SetXYZM(energy.Px(), energy.Py(), energy.Pz(), 0);        
+
+	for(itStable = stable.begin(); itStable != stable.end(); ++itStable) {  
+	  temp1.SetXYZM((*itStable)->momentum().px(), (*itStable)->momentum().py(), (*itStable)->momentum().pz(), 0);        
+	  double DR = temp1.DeltaR(temp2);
+	  if (DR < dRTkMax) counter++;        
+	}
+
+	if(acceptPrompts) {
+
+	    if ((*itSeed)->momentum().perp()>promptPtThreshold)
+	    {
+	      bool isPrompt=true;
+	      this_id = (*itSeed)->pdg_id();
+	      mom = (*itSeed);
+	      while (mom->pdg_id() == this_id) {
+
+		const GenParticle* mother = mom->production_vertex() ?       
+		  *(mom->production_vertex()->particles_in_const_begin()) : nullptr;
+
+		mom = mother;
+		if (mom == nullptr) {
+		  break;
+		}	  
+	      }
+
+	      first_different_id = mom->pdg_id();
+
+	      if (mom->status() == 2 && std::abs(first_different_id)>100) isPrompt=false;
+
+
+	      if(isPrompt) counter=0;
+	    }
+	}
+      }
+    }
+
+    candidate.push_back(energy);
+    candidateSeed.push_back(tempseed);
+    candidateNarrow.push_back(narrowCone);
+    nTracks.push_back(counter);
+  } 
+
+  if (candidate.size() <2) return accepted;
+
+  TLorentzVector minvMin, minvMax;
+
+  int i1, i2;
+  for(unsigned int i=0; i<candidate.size()-1; ++i) {
+
+    if (candidate[i].Energy() < energyCut) continue;
+    if(nTracks[i]>nTkConeMax) continue;
+    if (fabs(candidate[i].Eta()) > etaMaxCandidate) continue;
+
+    for(unsigned int j=i+1; j<candidate.size(); ++j) { 
+      if (candidate[j].Energy() < energyCut) continue;
+      if(nTracks[j]>nTkConeMax) continue;
+      if (fabs(candidate[j].Eta()) > etaMaxCandidate) continue;
+
+      if (nTracks[i] + nTracks[j] > nTkConeSum) continue;
+
+      if (candidate[i].Pt() > candidate[j].Pt()) {
+	i1 = i;
+	i2 = j;
+      } 
+      else {
+	i1 = j;
+	i2 = i;
+      }
+
+      if (candidate[i1].Pt() < pTMinCandidate1 || candidate[i2].Pt() < pTMinCandidate2) continue;
+
+      minvMin = candidate[i] + candidate[j];
+      if (minvMin.M() < invMassMin) continue;
+
+      minvMax = candidate[i] + candidate[j];
+      if (minvMax.M() > invMassMax) continue;
+
+      accepted = true;
+
+    }
+  }
+
+  return accepted;
+}
+

GeneratorInterface/GenFilters/interface/PythiaFilterGammaGamma.h

@@ -16,8 +16,10 @@
 #include "FWCore/Framework/interface/MakerMacros.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
+#include "GeneratorInterface/Core/interface/PythiaHepMCFilterGammaGamma.h"
+
 namespace edm {
-	  class HepMCProduct;
+  class HepMCProduct;
 }
 
 class PythiaFilterGammaGamma : public edm::global::EDFilter<> {
@@ -29,15 +31,14 @@ class PythiaFilterGammaGamma : public edm::global::EDFilter<> {
  private:
 
   const edm::EDGetTokenT<edm::HepMCProduct> token_;
-  const int maxEvents;
-
-  const double ptSeedThr, etaSeedThr, ptGammaThr, etaGammaThr, ptTkThr, etaTkThr;
-  const double ptElThr, etaElThr, dRTkMax, dRSeedMax, dPhiSeedMax, dEtaSeedMax, dRNarrowCone, pTMinCandidate1, pTMinCandidate2, etaMaxCandidate;
-  const double invMassMin, invMassMax;
-  const double energyCut;
-  const int nTkConeMax, nTkConeSum;
-  const bool acceptPrompts;
-  const double promptPtThreshold;
+
+  /** the actual implementation of the filter,
+      adapted to be used with HepMCFilterDriver.
+
+      We make this a pointer because EDFilter::filter() is const
+      while BaseHepMCFilter::filter() is not.
+ */
+  std::auto_ptr<PythiaHepMCFilterGammaGamma> hepMCFilter_;
 
 };
 #endif