GenMuonRadiationAlgorithm Class Reference

#include <GenMuonRadiationAlgorithm.h>

Public Member Functions
reco::Candidate::LorentzVector	compFSR (const edm::StreamID &streamID, const reco::Candidate::LorentzVector &, int, const reco::Candidate::LorentzVector &, int &)
	GenMuonRadiationAlgorithm (const edm::ParameterSet &)
	~GenMuonRadiationAlgorithm ()

Private Types
enum	{ kPYTHIA, kPHOTOS }

Private Attributes
double	beamEnergy_
int	mode_
gen::PhotosInterfaceBase *	photos_
myPythia6ServiceWithCallback *	pythia_
int	verbosity_

Static Private Attributes
static bool	photos_isInitialized_ = false
static bool	pythia_isInitialized_ = false

Detailed Description

Auxiliary class to correct for muon –> muon + photon final state radiation (FSR) in selected Z –> mu+ mu- events. The FSR is estimated on event-by-event basis via a Monte Carlo technique: for each reconstructed muon, PHOTOS is used to obtain a random amount of radiated photon energy. The radiated photon energy is then added to the energy of the reconstructed muon and the energy of generator level tau leptons is set to the sum.

NOTE: FSR increases with the energy of the muon. So, in principle the energy of muon + photon (energy of muon before FSR), not the energy of the reconstructed muon (energy of muon after FSR) would need to be used as input to PHOTOS. As the amount of radiated photon energy is typically small (< 1 GeV on average), taking the energy of the reconstructed muon is a good approximation.

Author

Christian Veelken, LLR

Version

Revision:

1.6

Id:

GenMuonRadiationAlgorithm.h,v 1.6 2013/01/27 13:53:44 veelken Exp

Definition at line 36 of file GenMuonRadiationAlgorithm.h.

Member Enumeration Documentation

anonymous enum

private

Enumerator
kPYTHIA
kPHOTOS

Definition at line 47 of file GenMuonRadiationAlgorithm.h.

{ kPYTHIA, kPHOTOS };

Constructor & Destructor Documentation

GenMuonRadiationAlgorithm::GenMuonRadiationAlgorithm ( const edm::ParameterSet &  cfg )

explicit

Definition at line 163 of file GenMuonRadiationAlgorithm.cc.

References edm::hlt::Exception, edm::ParameterSet::exists(), reco::get(), edm::ParameterSet::getParameter(), kPHOTOS, kPYTHIA, mode_, photos_, pythia_, AlCaHLTBitMon_QueryRunRegistry::string, and verbosity_.

  : beamEnergy_(cfg.getParameter<double>("beamEnergy")),
    photos_(0),
    pythia_(0)
{
  std::string mode_string = cfg.getParameter<std::string>("mode");
  if      ( mode_string == "pythia" ) mode_ = kPYTHIA;
  else if ( mode_string == "photos" ) mode_ = kPHOTOS;
  else throw cms::Exception("Configuration")
    << " Invalid Configuration Parameter 'mode' = " << mode_string << " !!\n";

  if ( mode_ == kPYTHIA ) pythia_ = new myPythia6ServiceWithCallback(cfg);
  if ( mode_ == kPHOTOS ){
    photos_ =  (gen::PhotosInterfaceBase*)(PhotosFactory::get()->create("Photos2155", cfg.getParameter<edm::ParameterSet>("PhotosOptions")));
    //settings?
    //usesResource(edm::SharedResourceNames::kPhotos);
  }

  verbosity_ = ( cfg.exists("verbosity") ) ? 
    cfg.getParameter<int>("verbosity") : 0;
}

GenMuonRadiationAlgorithm::~GenMuonRadiationAlgorithm

(

)

Definition at line 185 of file GenMuonRadiationAlgorithm.cc.

References photos_, and pythia_.

{
  delete photos_;
  call_pystat(1);
  delete pythia_;
}

Member Function Documentation

reco::Candidate::LorentzVector GenMuonRadiationAlgorithm::compFSR	(	const edm::StreamID &	streamID,
		const reco::Candidate::LorentzVector &	muonP4,
		int	muonCharge,
		const reco::Candidate::LorentzVector &	otherP4,
		int &	errorFlag
	)

Definition at line 227 of file GenMuonRadiationAlgorithm.cc.

References gen::PhotosInterfaceBase::apply(), beamEnergy_, call_pyedit(), gen::call_pygive(), gen::call_pylist(), conv, gather_cfg::cout, decayRandomEngine, GenParticle::GenParticle, edm::RandomNumberGenerator::getEngine(), gen::PhotosInterfaceBase::init(), edm::Service< T >::isAvailable(), kPHOTOS, kPYTHIA, mode_, photos_, photos_isInitialized_, protonMass, pythia_, pythia_isInitialized_, repairBarcodes(), myPythia6ServiceWithCallback::setEvent(), gen::Pythia6Service::setGeneralParams(), gen::PhotosInterfaceBase::setRandomEngine(), mathSSE::sqrt(), Clusterizer1DCommons::square(), and verbosity_.

Referenced by GenMuonRadCorrAnalyzer::analyze(), and ParticleReplacerZtautau::produce().

{
  edm::Service<edm::RandomNumberGenerator> rng;
  if ( !rng.isAvailable() ) 
    throw cms::Exception("Configuration")
      << "The GenMuonRadiationAlgorithm module requires the RandomNumberGeneratorService\n"
      << "which appears to be absent. Please add that service to your configuration\n"
      << "or remove the modules that require it.\n";
  
  // this is a global variable defined in GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc
  CLHEP::HepRandomEngine& decayRandomEngine = rng->getEngine(streamID);
  photos_->setRandomEngine(&decayRandomEngine);

  if ( verbosity_ ) {
    std::cout << "<GenMuonRadiationAlgorithm::compMuonFSR>:" << std::endl;
    std::cout << " muon: En = " << muonP4.E() << ", Pt = " << muonP4.pt() << ", eta = " << muonP4.eta() << ", phi = " << muonP4.phi() << ", charge = " << muonCharge << std::endl;
    std::cout << " other: En = " << otherP4.E() << ", Pt = " << otherP4.pt() << ", eta = " << otherP4.eta() << ", phi = " << otherP4.phi() << std::endl;
  }

  reco::Candidate::LorentzVector muonP4_limited = getP4_limited(muonP4, muonP4.mass());  
  int muonPdgId = -13*muonCharge;
  reco::Candidate::LorentzVector otherP4_limited = getP4_limited(otherP4, otherP4.mass());
  int otherPdgId = +13*muonCharge;
  reco::Candidate::LorentzVector zP4 = muonP4_limited + otherP4_limited;
  int zPdgId = 23;
  if ( verbosity_ ) {
    std::cout << "muon(limited): En = " << muonP4_limited.E() << ", Pt = " << muonP4_limited.pt() << ", eta = " << muonP4_limited.eta() << ", phi = " << muonP4_limited.phi() << std::endl;
    std::cout << "other: En = " << otherP4_limited.E() << ", Pt = " << otherP4_limited.pt() << ", eta = " << otherP4_limited.eta() << ", phi = " << otherP4_limited.phi() << std::endl;
    std::cout << "Z: En = " << zP4.E() << ", Pt = " << zP4.pt() << ", eta = " << zP4.eta() << ", phi = " << zP4.phi() << std::endl;
  }

  HepMC::GenEvent* genEvt_beforeFSR = new HepMC::GenEvent();
  HepMC::GenVertex* genVtx_in = new HepMC::GenVertex(HepMC::FourVector(0.,0.,0.,0.));
  double protonEn = beamEnergy_;
  double protonPz = sqrt(square(protonEn) - square(protonMass));
  genVtx_in->add_particle_in(new HepMC::GenParticle(HepMC::FourVector(0., 0., +protonPz, protonEn), 2212, 3));
  genVtx_in->add_particle_in(new HepMC::GenParticle(HepMC::FourVector(0., 0., -protonPz, protonEn), 2212, 3));
  genEvt_beforeFSR->add_vertex(genVtx_in);  
  reco::Candidate::LorentzVector ppP4_scattered(-zP4.px(), -zP4.py(), -zP4.pz(), 2.*beamEnergy_ - zP4.E());
  double protonEn_rf = 0.5*ppP4_scattered.mass();
  double protonP_rf  = sqrt(square(protonEn_rf) - square(protonMass));
  reco::Candidate::LorentzVector proton1P4_scattered_rf(0., 0., +protonP_rf, protonEn_rf);
  reco::Candidate::LorentzVector proton2P4_scattered_rf(0., 0., -protonP_rf, protonEn_rf);
  reco::Candidate::LorentzVector proton1P4_scattered = boostToLab(ppP4_scattered, proton1P4_scattered_rf);
  reco::Candidate::LorentzVector proton2P4_scattered = boostToLab(ppP4_scattered, proton2P4_scattered_rf);
  genVtx_in->add_particle_out(new HepMC::GenParticle(HepMC::FourVector(proton1P4_scattered.px(), proton1P4_scattered.py(), proton1P4_scattered.pz(), proton1P4_scattered.E()), 2212, 1));
  genVtx_in->add_particle_out(new HepMC::GenParticle(HepMC::FourVector(proton2P4_scattered.px(), proton2P4_scattered.py(), proton2P4_scattered.pz(), proton2P4_scattered.E()), 2212, 1));  
  HepMC::GenVertex* genVtx_out = new HepMC::GenVertex(HepMC::FourVector(0.,0.,0.,0.));
  HepMC::GenParticle* genZ = new HepMC::GenParticle((HepMC::FourVector)zP4, zPdgId, 2, HepMC::Flow(), HepMC::Polarization(0,0));
  genVtx_in->add_particle_out(genZ);
  genVtx_out->add_particle_in(genZ);
  HepMC::GenParticle* genMuon = new HepMC::GenParticle((HepMC::FourVector)muonP4_limited, muonPdgId, 1, HepMC::Flow(), HepMC::Polarization(0,0));
  genVtx_out->add_particle_out(genMuon);
  HepMC::GenParticle* genOther = new HepMC::GenParticle((HepMC::FourVector)otherP4_limited, otherPdgId, 1, HepMC::Flow(), HepMC::Polarization(0,0));
  genVtx_out->add_particle_out(genOther);
  genEvt_beforeFSR->add_vertex(genVtx_out);
  if ( verbosity_ ) {
    std::cout << "genEvt (beforeFSR):" << std::endl;
    genEvt_beforeFSR->print(std::cout);
  }

  repairBarcodes(genEvt_beforeFSR);

  HepMC::GenEvent* genEvt_afterFSR = 0;

  if ( mode_ == kPYTHIA ) {
    gen::Pythia6Service::InstanceWrapper pythia6InstanceGuard(pythia_);
    
    if ( !pythia_isInitialized_ ) {
      pythia_->setGeneralParams();
      //pythia_->setCSAParams();
      //pythia_->setSLHAParams();   
      gen::call_pygive("MSTU(12) = 12345"); 
      call_pyinit("USER", "", "", 0.);
      pythia_isInitialized_ = true;
    }

    pythia_->setEvent(genEvt_beforeFSR);
    if ( verbosity_ ) {
      std::cout << "HEPEUP (beforeFSR):" << std::endl;
      call_pylist(7);
      std::cout.flush();
    }
    call_pyevnt();
    call_pyedit(0); // CV: remove comment lines (confuses conversion to HepMC::Event)
    if ( verbosity_ ) {
      std::cout << "PYJETS (afterFSR):" << std::endl;
      call_pylist(1);
      std::cout.flush();
    }
    call_pyhepc(1); 
    if ( verbosity_ ) {
      std::cout << "HepMC (afterFSR):" << std::endl;
      call_pylist(5);
      std::cout.flush();
    }

    HepMC::IO_HEPEVT conv;
    conv.set_print_inconsistency_errors(false);
    genEvt_afterFSR = conv.read_next_event();
  } else if ( mode_ == kPHOTOS ) {    

    if ( !photos_isInitialized_ ) {
      photos_->init();
      photos_isInitialized_ = true;
    }
    
    HepMC::IO_HEPEVT conv;
    conv.write_event(genEvt_beforeFSR); 
    genEvt_afterFSR = photos_->apply(genEvt_beforeFSR);
  } else assert(0);

  if ( verbosity_ ) {
    std::cout << "genEvt (afterFSR):" << std::endl;
    genEvt_afterFSR->print(std::cout);
  }

  // find lowest energy muon
  // (= muon after FSR) 
  reco::Candidate::LorentzVector muonP4afterFSR = muonP4;
  for ( HepMC::GenEvent::particle_iterator genParticle = genEvt_afterFSR->particles_begin();
    genParticle != genEvt_afterFSR->particles_end(); ++genParticle ) {    
    if ( (*genParticle)->pdg_id() == muonPdgId && (*genParticle)->momentum().e() < muonP4afterFSR.E() ) {
      muonP4afterFSR.SetPxPyPzE(
    (*genParticle)->momentum().px(), 
    (*genParticle)->momentum().py(), 
    (*genParticle)->momentum().pz(), 
    (*genParticle)->momentum().e());
    }    
  }
  reco::Candidate::LorentzVector sumPhotonP4 = muonP4 - muonP4afterFSR;
  if ( verbosity_ ) {
    std::cout << "sumPhotons: En = " << sumPhotonP4.E() << ", Pt = " << sumPhotonP4.pt() << ", eta = " << sumPhotonP4.eta() << ", phi = " << sumPhotonP4.phi() << std::endl;
  }

  delete genEvt_beforeFSR;
  if ( genEvt_afterFSR != genEvt_beforeFSR ) delete genEvt_afterFSR;
    
  return sumPhotonP4;
}

Member Data Documentation

double GenMuonRadiationAlgorithm::beamEnergy_

private

Definition at line 45 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR().

int GenMuonRadiationAlgorithm::mode_

private

Definition at line 48 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR(), and GenMuonRadiationAlgorithm().

gen::PhotosInterfaceBase* GenMuonRadiationAlgorithm::photos_

private

Definition at line 50 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR(), GenMuonRadiationAlgorithm(), and ~GenMuonRadiationAlgorithm().

bool GenMuonRadiationAlgorithm::photos_isInitialized_ = false

staticprivate

Definition at line 51 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR().

myPythia6ServiceWithCallback* GenMuonRadiationAlgorithm::pythia_

private

Definition at line 53 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR(), GenMuonRadiationAlgorithm(), and ~GenMuonRadiationAlgorithm().

bool GenMuonRadiationAlgorithm::pythia_isInitialized_ = false

staticprivate

Definition at line 54 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR().

int GenMuonRadiationAlgorithm::verbosity_

private

Definition at line 56 of file GenMuonRadiationAlgorithm.h.

Referenced by compFSR(), and GenMuonRadiationAlgorithm().