GenFilters/plugins/HerwigMaxPtPartonFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HerwigMaxPtPartonFilter
0004 // Class:      HerwigMaxPtPartonFilter
0005 //
0006 /**\class HerwigMaxPtPartonFilter HerwigMaxPtPartonFilter.cc IOMC/HerwigMaxPtPartonFilter/src/HerwigMaxPtPartonFilter.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Brian Dorney
0015 //         Created:  July 27th 2010
0016 // $Id: HerwigMaxPtPartonFilter.h v1.0
0017 //
0018 // Modified From: PythiaFilter.cc
0019 //
0020 // Special Thanks to Filip Moortgat
0021 //
0022 /*
0023 Date: July 29th 2010
0024 Version: 2.2
0025 First Release In: CMSSW_3_8_X
0026
0027 PURPOSE: This Filter is designed to run on Herwig Monte Carlo Event Files
0028 (Pythia status codes are assumed to NOT BE EMULATED!!!!)
0029
0030 For a description of Herwig Status Codes, See:
0031 https://arxiv.org/abs/hep-ph/0011363
0032 (Section 8.3.1)
0033
0034 This Filter Finds all final state quarks (pdg_id=1,2,3,4 or 5, status=158 or 159) with Pt>1 GeV/c
0035 that occur before the first cluster (pdg_id=91) appears in the event cascade. This is done per event.
0036
0037 Then a histogram (which is RESET EACH EVENT) 2D histogram is formed, the Final State quarks
0038 are then plotted in eta-phi space.  These histogram entries are weighted by the quark Pt.
0039
0040 This forms a 2D eta-phi "jet" topology for each event, and acts as a very rudimentary jet algorithm
0041
0042 The maximum bin entry (i.e. "jet") in this histogram is the highest Pt "Jet" in the event.
0043
0044 This is then used for filtering.
0045
0046 The size of each bin in this 2D histogram corresponds roughly to a cone radius of 0.5
0047
0048 i.e. This Filter Checks:
0049 minptcut <= Highest Pt "Jet" < maxptcut
0050
0051 If this is true, the event is accepted.
0052 */
0053
0054 #include "FWCore/Framework/interface/Event.h"
0055 #include "FWCore/Framework/interface/Frameworkfwd.h"
0056 #include "FWCore/Framework/interface/stream/EDFilter.h"
0057 #include "FWCore/Framework/interface/MakerMacros.h"
0058 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0059 #include "FWCore/Utilities/interface/EDGetToken.h"
0060 #include "FWCore/Utilities/interface/InputTag.h"
0061 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0062
0063 #include "TH2.h"
0064
0065 #include <cmath>
0066 #include <cstdlib>
0067
0068 class HerwigMaxPtPartonFilter : public edm::stream::EDFilter<> {
0069 public:
0070   explicit HerwigMaxPtPartonFilter(const edm::ParameterSet&);
0071
0072   bool filter(edm::Event&, const edm::EventSetup&) override;
0073
0074 private:
0075   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0076
0077   const double minptcut;
0078   const double maxptcut;
0079   const int processID;
0080
0081   TH2D hFSPartons_JS_PtWgting;
0082 };
0083
0084 using namespace edm;
0085 using namespace std;
0086
0087 HerwigMaxPtPartonFilter::HerwigMaxPtPartonFilter(const edm::ParameterSet& iConfig)
0088     : token_(consumes<edm::HepMCProduct>(
0089           iConfig.getUntrackedParameter("moduleLabel", edm::InputTag("generator", "unsmeared")))),
0090       minptcut(iConfig.getUntrackedParameter("MinPt", 0.)),
0091       maxptcut(iConfig.getUntrackedParameter("MaxPt", 10000.)),
0092       processID(iConfig.getUntrackedParameter("ProcessID", 0)),
0093       hFSPartons_JS_PtWgting(
0094           "hFSPartons_JS_PtWgting", "#phi-#eta Space of FS Partons (p_{T} wgt'ing)", 20, -5.205, 5.205, 32, -M_PI, M_PI) {
0095 }
0096
0097 bool HerwigMaxPtPartonFilter::filter(edm::Event& iEvent, const edm::EventSetup&) {
0098   //Histogram, reset each event
0099   hFSPartons_JS_PtWgting.Reset();
0100
0101   bool accepted = false;     //The Accept/Reject Variable
0102   bool isFSQuark = false;    //Keeps track of whether a particle is a Final State Quark
0103   double maxPartonPt = 0.0;  //Self Explanatory
0104
0105   //int ChosenPartonId=0, ChosenPartonSt=0;
0106
0107   int pos1stCluster = 0;  //keeps track of the position of the first herwig cluster in the event
0108
0109   //This is when Hadronization within the event occurs.
0110   long counter = 0;  //keeps track of the particle index in the event
0111
0112   Handle<HepMCProduct> evt;
0113   iEvent.getByToken(token_, evt);
0114
0115   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0116
0117   if (processID == 0 || processID == myGenEvent->signal_process_id()) {  //processId if statement
0118
0119     //Find the Position of the 1st Herwig Cluster
0120     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0121          ++p) {
0122       if (abs((*p)->pdg_id()) == 91) {
0123         break;
0124       }
0125       pos1stCluster++;  //Starts at Zero, like the Collection
0126     }
0127
0128     //Loop through the all particles in the event
0129     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0130          ++p) {
0131       //"Garbage" Cut, 1 GeV/c Pt Cut on All Particles Considered
0132       if ((*p)->momentum().perp() > 1.0) {
0133         //Final State Quark criterion
0134         if (abs((*p)->pdg_id()) == 1 || abs((*p)->pdg_id()) == 2 || abs((*p)->pdg_id()) == 3 ||
0135             abs((*p)->pdg_id()) == 4 || abs((*p)->pdg_id()) == 5) {
0136           if (counter < pos1stCluster && ((*p)->status() == 158 || (*p)->status() == 159)) {
0137             isFSQuark = true;
0138           }
0139         }  //end if FS Quark criterion
0140       }    //End "Garbage" Cut
0141
0142       if (isFSQuark) {
0143         hFSPartons_JS_PtWgting.Fill(
0144             (*p)->momentum().eta(), (*p)->momentum().phi(), (*p)->momentum().perp());  //weighted by Particle Pt
0145       }
0146
0147       counter++;
0148       isFSQuark = false;
0149     }  //end all particles loop
0150
0151     maxPartonPt = hFSPartons_JS_PtWgting.GetMaximum();
0152
0153     //The Actual Filtering Process
0154     if (maxPartonPt >= minptcut && maxPartonPt < maxptcut) {
0155       accepted = true;  //Accept the Event
0156
0157     }  //End Filtering
0158   }    //end processId if statement
0159
0160   else {
0161     accepted = true;
0162   }
0163   return accepted;
0164 }
0165
0166 DEFINE_FWK_MODULE(HerwigMaxPtPartonFilter);