PythiaFilterZgamma.cc

Go to the documentation of this file.

#include "GeneratorInterface/GenFilters/interface/PythiaFilterZgamma.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include <iostream>
#include<list>
#include<vector>
#include<cmath>

PythiaFilterZgamma::PythiaFilterZgamma(const edm::ParameterSet& iConfig) :
token_(consumes<edm::HepMCProduct>(edm::InputTag(iConfig.getUntrackedParameter("moduleLabel",std::string("generator")),"unsmeared"))),
selProc(iConfig.getUntrackedParameter<int>("SelectProcess")),
ptElMin(iConfig.getUntrackedParameter<double>("MinElPt", 5.0)),
ptMuMin(iConfig.getUntrackedParameter<double>("MinMuPt", 3.0)),
ptPhotonMin(iConfig.getUntrackedParameter<double>("MinPhotPt", 5.0)),
etaElMax(iConfig.getUntrackedParameter<double>("MaxElecEta", 2.7)),
etaMuMax(iConfig.getUntrackedParameter<double>("MaxMuonEta", 2.4)),
etaPhotonMax(iConfig.getUntrackedParameter<double>("MaxPhotEta", 2.7))
{  
  theNumberOfSelected = 0;
}


PythiaFilterZgamma::~PythiaFilterZgamma(){}


// ------------ method called to produce the data  ------------
bool PythiaFilterZgamma::filter(edm::Event& iEvent, const edm::EventSetup& iSetup){


  bool accepted = false;
  edm::Handle<edm::HepMCProduct> evt;
  iEvent.getByToken(token_, evt);

  const HepMC::GenEvent * myGenEvent = evt->GetEvent();

  std::vector<const HepMC::GenParticle *> el;
  std::vector<const HepMC::GenParticle *> mu;
  std::vector<const HepMC::GenParticle *> gam;

  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();  
               p != myGenEvent->particles_end(); ++p ) {
              
    if(selProc == 1 ) {
      if ( std::abs((*p)->pdg_id())==11 && (*p)->status()==1 )
        el.push_back(*p);
      if(el.size()>1) break;
    } else if(selProc == 2 ) {
      if ( std::abs((*p)->pdg_id())==13 && (*p)->status()==1 )
        mu.push_back(*p);
      if(mu.size()>1) break;        
    }
      
  } // end of first particle loop for finding Z0 decays

  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();  
               p != myGenEvent->particles_end(); ++p ) {

    if ( std::abs((*p)->pdg_id())==22 && (*p)->status()==1 )
      gam.push_back(*p);
    if(gam.size()==1) break;

  } // end of second particle loop for finding gamma

  if(selProc == 1 ) {
    if(el.size()>1){
      double ptEl1 = el[0]->momentum().perp();
      double ptEl2 = el[1]->momentum().perp();
      double etaEl1 = el[0]->momentum().eta();
      double etaEl2 = el[1]->momentum().eta();
      if (ptEl1 > ptElMin && ptEl2 > ptElMin  && 
    std::abs(etaEl1) < etaElMax &&
    std::abs(etaEl2) < etaElMax) {
        if(gam.size()==1) {
      double ptGam = gam[0]->momentum().perp();
          double etaGam = gam[0]->momentum().eta();
      if(ptGam > ptPhotonMin && std::abs(etaGam) < etaPhotonMax)
                 accepted=true;
    }     
      } 
    } 
  } else if(selProc == 2 ) {

    if(mu.size()>1){
      double ptMu1 = mu[0]->momentum().perp();
      double ptMu2 = mu[1]->momentum().perp();
      double etaMu1 = mu[0]->momentum().eta();
      double etaMu2 = mu[1]->momentum().eta();
      if (ptMu1 > ptMuMin && ptMu2 > ptMuMin  && 
    std::abs(etaMu1) < etaMuMax &&
    std::abs(etaMu2) < etaMuMax) {
        if(gam.size()==1) {
      double ptGam = gam[0]->momentum().perp();
          double etaGam = gam[0]->momentum().eta();
      if(ptGam > ptPhotonMin && std::abs(etaGam) < etaPhotonMax)
                 accepted=true;
    }     
      } 
    } 
    
  }
/*
  if(accepted) {
    std::cout << "Accepted event Number: " << theNumberOfSelected 
              << "  of category " << selProc << std::endl;
    if(selProc == 1 ) { 
      std::cout << "Electon 1 pt, eta = " << el[0]->momentum().perp() << ", " 
                <<  el[0]->momentum().eta() << std::endl;
      std::cout << "Electon 2 pt, eta = " << el[1]->momentum().perp() << ", " 
                <<  el[1]->momentum().eta() << std::endl;
      std::cout << "Photon pt, eta = " << gam[0]->momentum().perp() << ", " 
                <<  gam[0]->momentum().eta() << std::endl;
    } else if(selProc == 2 ) {
      std::cout << "Muon 1 pt, eta = " << mu[0]->momentum().perp() << ", " 
                <<  mu[0]->momentum().eta() << std::endl;
      std::cout << "Muon 2 pt, eta = " << mu[1]->momentum().perp() << ", " 
                <<  mu[1]->momentum().eta() << std::endl;
      std::cout << "Photon pt, eta = " << gam[0]->momentum().perp() << ", " 
                <<  gam[0]->momentum().eta() << std::endl;

    }  
  }
*/
  if (accepted) {
    theNumberOfSelected++;
    return true; 
  }
  else return false;

}