Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:56

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);