Included the betaboost in the MCSingleParticle filters

cms-sw · May 1, 2017 · 22ccc08 · 22ccc08
1 parent fb82d34
commit 22ccc08
Show file tree

Hide file tree

Showing 4 changed files with 46 additions and 6 deletions.
diff --git a/GeneratorInterface/GenFilters/interface/MCSingleParticleFilter.h b/GeneratorInterface/GenFilters/interface/MCSingleParticleFilter.h
@@ -38,6 +38,9 @@
 namespace edm {
   class HepMCProduct;
 }
+namespace HepMC {
+   class FourVector;
+}
 
 class MCSingleParticleFilter : public edm::EDFilter {
    public:
@@ -47,6 +50,9 @@ class MCSingleParticleFilter : public edm::EDFilter {
 
       virtual bool filter(edm::Event&, const edm::EventSetup&);
    private:
+      // ----------memeber function----------------------
+       HepMC::FourVector zboost(const HepMC::FourVector&);
+
       // ----------member data ---------------------------
 
        edm::EDGetTokenT<edm::HepMCProduct> token_;
@@ -55,5 +61,6 @@ class MCSingleParticleFilter : public edm::EDFilter {
        std::vector<double> etaMin;  
        std::vector<double> etaMax;
        std::vector<int> status;
+       double betaBoost;
 };
 #endif
diff --git a/GeneratorInterface/GenFilters/interface/MCSmartSingleParticleFilter.h b/GeneratorInterface/GenFilters/interface/MCSmartSingleParticleFilter.h
@@ -35,6 +35,9 @@
 namespace edm {
   class HepMCProduct;
 }
+namespace HepMC {
+   class FourVector;
+}
 
 class MCSmartSingleParticleFilter : public edm::EDFilter {
    public:
@@ -44,6 +47,9 @@ class MCSmartSingleParticleFilter : public edm::EDFilter {
 
       virtual bool filter(edm::Event&, const edm::EventSetup&);
    private:
+      // ----------memeber function----------------------
+       HepMC::FourVector zboost(const HepMC::FourVector&);
+
       // ----------member data ---------------------------
 
        edm::EDGetTokenT<edm::HepMCProduct> token_;
@@ -57,5 +63,6 @@ class MCSmartSingleParticleFilter : public edm::EDFilter {
        std::vector<double> decayRadiusMax;
        std::vector<double> decayZMin;  
        std::vector<double> decayZMax;
+       double betaBoost;
 };
 #endif
diff --git a/GeneratorInterface/GenFilters/src/MCSingleParticleFilter.cc b/GeneratorInterface/GenFilters/src/MCSingleParticleFilter.cc
@@ -9,7 +9,8 @@ using namespace std;
 
 
 MCSingleParticleFilter::MCSingleParticleFilter(const edm::ParameterSet& iConfig) :
-token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared")))
+token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared"))),
+betaBoost(iConfig.getUntrackedParameter("BetaBoost",0.))
 {
    //here do whatever other initialization is needed
    vector<int> defpid ;
@@ -95,8 +96,9 @@ bool MCSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSetup& i
      for (unsigned int i = 0; i < particleID.size(); i++){
        if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {
 
-	 if ( (*p)->momentum().perp() > ptMin[i] && (*p)->momentum().eta() > etaMin[i] 
-	      && (*p)->momentum().eta() < etaMax[i] && ((*p)->status() == status[i] || status[i] == 0)) { 
+         HepMC::FourVector mom = zboost((*p)->momentum());
+	 if ( mom.perp() > ptMin[i] && mom.eta() > etaMin[i] 
+	      && mom.eta() < etaMax[i] && ((*p)->status() == status[i] || status[i] == 0)) { 
           accepted = true; 
 	 }  
 
@@ -110,3 +112,14 @@ bool MCSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSetup& i
 
 }
 
+
+HepMC::FourVector MCSingleParticleFilter::zboost(const HepMC::FourVector& mom) {
+   //Boost this Lorentz vector (from TLorentzVector::Boost)
+   double b2 = betaBoost*betaBoost;
+   double gamma = 1.0 / sqrt(1.0 - b2);
+   double bp = betaBoost*mom.pz();
+   double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
+
+   return HepMC::FourVector(mom.px(), mom.py(), mom.pz() + gamma2*bp*betaBoost + gamma*betaBoost*mom.e(), gamma*(mom.e()+bp));
+}
+
diff --git a/GeneratorInterface/GenFilters/src/MCSmartSingleParticleFilter.cc b/GeneratorInterface/GenFilters/src/MCSmartSingleParticleFilter.cc
@@ -9,7 +9,8 @@ using namespace std;
 
 
 MCSmartSingleParticleFilter::MCSmartSingleParticleFilter(const edm::ParameterSet& iConfig) :
-token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared")))
+token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared"))),
+betaBoost(iConfig.getUntrackedParameter("BetaBoost",0.))
 {
    //here do whatever other initialization is needed
    vector<int> defpid ;
@@ -150,8 +151,9 @@ bool MCSmartSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSet
      for (unsigned int i = 0; i < particleID.size(); i++){
        if (particleID[i] == (*p)->pdg_id() || particleID[i] == 0) {
 
-	 if ( (*p)->momentum().rho() > pMin[i] && (*p)->momentum().perp() > ptMin[i]
-	      && (*p)->momentum().eta() > etaMin[i] && (*p)->momentum().eta() < etaMax[i]
+         HepMC::FourVector mom = zboost((*p)->momentum());
+	 if ( mom.rho() > pMin[i] && mom.perp() > ptMin[i]
+	      && mom.eta() > etaMin[i] && mom.eta() < etaMax[i]
 	      && ((*p)->status() == status[i] || status[i] == 0)) { 
 
 	   if (!((*p)->production_vertex())) continue;
@@ -179,3 +181,14 @@ bool MCSmartSingleParticleFilter::filter(edm::Event& iEvent, const edm::EventSet
 
 }
 
+
+HepMC::FourVector MCSmartSingleParticleFilter::zboost(const HepMC::FourVector& mom) {
+   //Boost this Lorentz vector (from TLorentzVector::Boost)
+   double b2 = betaBoost*betaBoost;
+   double gamma = 1.0 / sqrt(1.0 - b2);
+   double bp = betaBoost*mom.pz();
+   double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
+
+   return HepMC::FourVector(mom.px(), mom.py(), mom.pz() + gamma2*bp*betaBoost + gamma*betaBoost*mom.e(), gamma*(mom.e()+bp));
+}
+