Page MenuHomeHEPForge
Files F19244024

EvtPhotosEngine.cpp
No OneTemporary
Actions

Size
6 KB
Referenced Files
None
Subscribers
None

EvtPhotosEngine.cpp
View Options

//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 2011 University of Warwick, UK
//
// Module: EvtPhotosEngine
//
// Description: Interface to the PHOTOS external generator
//
// Modification history:
//
// John Back May 2011 Module created
//
//------------------------------------------------------------------------
#include "EvtGenModels/EvtPhotosEngine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtPhotonParticle.hh"
#include "EvtGenBase/EvtReport.hh"
#include "Photos.h"
#include "PhotosHepMCEvent.h"
#include "PhotosHepMCParticle.h"
#include "PhotosParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/SimpleVector.h"
#include "HepMC/Units.h"
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
using std::endl;
EvtPhotosEngine::EvtPhotosEngine(std::string photonType) {
_gammaId = EvtPDL::getId(photonType);
if (_gammaId == EvtId(-1, -1)) {
report(INFO,"EvtGen")<<"Error in EvtPhotosEngine. Do not recognise the photon type "
<<photonType<<". Setting this to \"gamma\". "<<endl;
_gammaId = EvtPDL::getId("gamma");
}
_mPhoton = EvtPDL::getMeanMass(_gammaId);
_initialised = false;
}
EvtPhotosEngine::~EvtPhotosEngine() {
}
void EvtPhotosEngine::initialise() {
if (_initialised == false) {
report(INFO,"EvtGen")<<"Initialising PHOTOS."<<endl;
Photos::initialize();
// Set minimum photon energy (50keV at 1 GeV scale)
Photos::setInfraredCutOff(50.0e-6);
// Increase the maximum possible value of the interference weight
Photos::maxWtInterference(4.0); // 2^n, where n = number of charges (+,-)
_initialised = true;
}
}
bool EvtPhotosEngine::doDecay(EvtParticle* theMother) {
if (theMother == 0) {return false;}
if (_initialised == false) {this->initialise();}
// Create a dummy HepMC GenEvent containing a single vertex, with the mother
// assigned as the incoming particle and its daughters as outgoing particles.
// We then pass this event to Photos for processing.
// It will return a modified version of the event, updating the momentum of
// the original particles and will contain any new photon particles.
// We add these extra photons to the mother particle daughter list.
// Skip running Photos if the particle has no daughters, since we can't add FSR.
int nDaug(theMother->getNDaug());
if (nDaug == 0) {return false;}
// Create the dummy event.
HepMC::GenEvent* theEvent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
// Create the decay "vertex".
HepMC::GenVertex* theVertex = new HepMC::GenVertex();
theEvent->add_vertex(theVertex);
// Add the mother particle as the incoming particle to the vertex.
HepMC::GenParticle* hepMCMother = this->createGenParticle(theMother, true);
theVertex->add_particle_in(hepMCMother);
// Find all daughter particles and assign them as outgoing particles to the vertex.
int iDaug(0);
for (iDaug = 0; iDaug < nDaug; iDaug++) {
EvtParticle* theDaughter = theMother->getDaug(iDaug);
HepMC::GenParticle* hepMCDaughter = this->createGenParticle(theDaughter, false);
theVertex->add_particle_out(hepMCDaughter);
}
// Now pass the event to Photos for processing
// Create a Photos event object
PhotosHepMCEvent photosEvent(theEvent);
// Run the Photos algorithm
photosEvent.process();
// See if Photos has created new particles. If not, do nothing extra.
int nDecayPart = theVertex->particles_out_size();
int iLoop(0);
if (nDecayPart > nDaug) {
// We have extra particles from Photos
// Get the iterator of outgoing particles for this vertex
HepMC::GenVertex::particles_out_const_iterator outIter;
for (outIter = theVertex->particles_out_const_begin();
outIter != theVertex->particles_out_const_end(); ++outIter) {
// Get the next HepMC GenParticle
HepMC::GenParticle *outParticle = *outIter;
// Get the three-momentum Photos result for this particle
double px(0.0), py(0.0), pz(0.0);
if (outParticle != 0) {
HepMC::FourVector HepMCP4 = outParticle->momentum();
px = HepMCP4.px();
py = HepMCP4.py();
pz = HepMCP4.pz();
}
// Create an empty 4-momentum vector for the new/modified daughters
EvtVector4R newP4;
if (iLoop < nDaug) {
// Original daughters
EvtParticle* daugParticle = theMother->getDaug(iLoop);
if (daugParticle != 0) {
// Keep the original particle mass, but set the three-momentum
// according to what Photos has modified. However, this will
// violate energy conservation (from what Photos has provided).
double mass = daugParticle->mass();
double energy = sqrt(mass*mass + px*px + py*py + pz*pz);
newP4.set(energy, px, py, pz);
// Set the new four-momentum (FSR applied)
daugParticle->setP4WithFSR(newP4);
}
} else {
// Extra photon particle. Setup the four-momentum object.
double energy = sqrt(_mPhoton*_mPhoton + px*px + py*py + pz*pz);
newP4.set(energy, px, py, pz);
// Create a new photon particle and add it to the list of daughters
EvtPhotonParticle* gamma = new EvtPhotonParticle();
gamma->init(_gammaId, newP4);
gamma->setFSRP4toZero();
gamma->addDaug(theMother); // Let the mother know about this new particle
}
// Increment the loop counter for detecting additional photon particles
iLoop++;
}
}
// Cleanup
theEvent->clear();
delete theEvent;
return true;
}
HepMC::GenParticle* EvtPhotosEngine::createGenParticle(EvtParticle* theParticle, bool incoming) {
// Method to create an HepMC::GenParticle version of the given EvtParticle.
if (theParticle == 0) {return 0;}
// Get the 4-momentum (E, px, py, pz) for the EvtParticle
EvtVector4R p4(0.0, 0.0, 0.0, 0.0);
if (incoming == true) {
p4 = theParticle->getP4Restframe();
} else {
p4 = theParticle->getP4();
}
// Convert this to the HepMC 4-momentum
double E = p4.get(0);
double px = p4.get(1);
double py = p4.get(2);
double pz = p4.get(3);
HepMC::FourVector hepMC_p4(px, py, pz, E);
int PDGInt = EvtPDL::getStdHep(theParticle->getId());
// Set the status flag for the particle. This is required, otherwise Photos++
// will crash from out-of-bounds array index problems.
int status = PhotosParticle::HISTORY;
if (incoming == false) {status = PhotosParticle::STABLE;}
HepMC::GenParticle* genParticle = new HepMC::GenParticle(hepMC_p4, PDGInt, status);
return genParticle;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:38 AM (4 h, 14 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6564497
Default Alt Text
EvtPhotosEngine.cpp (6 KB)

Event Timeline

HEPForge.org · .