From 68c5ba4baa871d59b238f3265483cb9cb18ebf72 Mon Sep 17 00:00:00 2001 From: Ivan Razumov Date: Tue, 19 Apr 2022 11:53:43 +0200 Subject: [PATCH] Drop PhotosInterface (legacy Fortran Photos) --- .../ExternalDecays/src/ExternalDecayDriver.cc | 10 - .../test/Py6GenFilter_Photos_cfg.py | 190 -------- .../test/Py8GenFilter_Photos_cfg.py | 87 ---- .../interface/PhotosInterface.h | 53 --- .../PhotosInterface/plugins/BuildFile.xml | 10 - .../plugins/PhotosInterface.cc | 432 ------------------ 6 files changed, 782 deletions(-) delete mode 100644 GeneratorInterface/ExternalDecays/test/Py6GenFilter_Photos_cfg.py delete mode 100644 GeneratorInterface/ExternalDecays/test/Py8GenFilter_Photos_cfg.py delete mode 100644 GeneratorInterface/PhotosInterface/interface/PhotosInterface.h delete mode 100644 GeneratorInterface/PhotosInterface/plugins/PhotosInterface.cc diff --git a/GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc b/GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc index f9a14fb12baff..428afc193e059 100644 --- a/GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc +++ b/GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc @@ -43,18 +43,8 @@ ExternalDecayDriver::ExternalDecayDriver(const ParameterSet& pset, edm::Consumes } else if (curSet == "Tauola" || curSet == "Tauolapp" || curSet == "Tauolapp114") { fTauolaInterface = std::unique_ptr( TauolaFactory::get()->create("Tauolapp114", pset.getUntrackedParameter(curSet), iCollector)); - fPhotosInterface = std::unique_ptr( - PhotosFactory::get()->create("Photos2155", pset.getUntrackedParameter(curSet))); - fPhotosInterface->configureOnlyFor(15); - fPhotosInterface->avoidTauLeptonicDecays(); exSharedResources.emplace_back(edm::SharedResourceNames::kTauola); exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - } else if (curSet == "Photos" || curSet == "Photos2155") { - if (!fPhotosInterface) { - fPhotosInterface = std::unique_ptr( - PhotosFactory::get()->create("Photos2155", pset.getUntrackedParameter(curSet))); - exSharedResources.emplace_back(edm::SharedResourceNames::kPhotos); - } } else if (curSet == "Photospp" || curSet == "Photospp356") { if (!fPhotosInterface) { fPhotosInterface = std::unique_ptr( diff --git a/GeneratorInterface/ExternalDecays/test/Py6GenFilter_Photos_cfg.py b/GeneratorInterface/ExternalDecays/test/Py6GenFilter_Photos_cfg.py deleted file mode 100644 index 8d171f018c170..0000000000000 --- a/GeneratorInterface/ExternalDecays/test/Py6GenFilter_Photos_cfg.py +++ /dev/null @@ -1,190 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -# -# WARNING: This is NOT an example for users - -# it's my private (JY) "development" cfg, for testing -# newly implemented PhotosInterface - which is NOT yet -# released via ExternalDecayDeriver -# - -process = cms.Process("TEST") -process.load("FWCore.Framework.test.cmsExceptionsFatal_cff") -process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") -#process.load("SimGeneral.HepPDTESSource.pdt_cfi") - -process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", - generator = cms.PSet( - initialSeed = cms.untracked.uint32(123456789), - engineName = cms.untracked.string('HepJamesRandom') - ) -) - - -# The following three lines reduce the clutter of repeated printouts -# of the same exception message. -process.load("FWCore.MessageLogger.MessageLogger_cfi") - -process.MessageLogger.cerr.enableStatistics = False - - -process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(20)) - -process.source = cms.Source("EmptySource") - -process.generator = cms.EDFilter("Pythia6GeneratorFilter", - pythiaHepMCVerbosity = cms.untracked.bool(True), - maxEventsToPrint = cms.untracked.int32(20), - pythiaPylistVerbosity = cms.untracked.int32(1), - # this shows how to turn ON some of the general Py6 printouts, like banner... - ## --> displayPythiaBanner = cms.untracked.bool(True), - ## --> displayPythiaCards = cms.untracked.bool(True), - comEnergy = cms.double(7000.0), - - ExternalDecays = cms.PSet( - Photos = cms.untracked.PSet(), - parameterSets = cms.vstring( "Photos" ) - ), - - PythiaParameters = cms.PSet( - - pythiaSimpleSettings = cms.vstring( - 'PMAS(5,1)=4.8 ! b quark mass', - 'PMAS(6,1)=172.3 ! t quark mass', - 'MSTP(61)=0 ! turn off initial state radiation', - 'mstj(41)=1' # per Steve M., instead of mstp(71)... btw, shoult it be 0 or 1 ? - #'MSTP(71)=0 ! final state radiation (or final-state showers ?)' - - ), - pythiaSpecialSettings = cms.vstring( - 'PMAS(25,1)=190.0 !mass of Higgs', - 'MSEL=0 !(D=1) to select between full user control (0, then use MSUB) and some preprogrammed alternative: QCD hight pT processes (1, then ISUB=11, 12, 13, 28, 53, 68), QCD low pT processes (2, then ISUB=11, 12, 13, 28, 53, 68, 91, 92, 94, 95)', - 'MSTJ(11)=3 !Choice of the fragmentation function', - 'MSTJ(41)=1 !Switch off Pythia QED bremsshtrahlung', - 'MSTP(51)=7 !structure function chosen', - 'MSTP(61)=0 ! no initial-state showers', - 'MSTP(71)=0 ! no final-state showers', - 'MSTP(81)=0 ! no multiple interactions', - 'MSTP(111)=0 ! no hadronization', - 'MSTU(21)=1 !Check on possible errors during program execution', - 'MSUB(102)=1 !ggH', - 'MSUB(123)=1 !ZZ fusion to H', - 'MSUB(124)=1 !WW fusion to H', - # these 4 below are irrelevant if the hadronization is off (mstp(111)=0) - 'PARP(82)=1.9 !pt cutoff for multiparton interactions', - 'PARP(83)=0.5 !Multiple interactions: matter distrbn parameter Registered by Chris.Seez@cern.ch', - 'PARP(84)=0.4 !Multiple interactions: matter distribution parameter Registered by Chris.Seez@cern.ch', - 'PARP(90)=0.16 !Multiple interactions: rescaling power Registered by Chris.Seez@cern.ch', - #........................................ - 'CKIN(45)=5. !high mass cut on m2 in 2 to 2 process Registered by Chris.Seez@cern.ch', - 'CKIN(46)=150. !high mass cut on secondary resonance m1 in 2->1->2 process Registered by Alexandre.Nikitenko@cern.ch', - 'CKIN(47)=5. !low mass cut on secondary resonance m2 in 2->1->2 process Registered by Alexandre.Nikitenko@cern.ch', - 'CKIN(48)=150. !high mass cut on secondary resonance m2 in 2->1->2 process Registered by Alexandre.Nikitenko@cern.ch' - ), - pythiaHZZDecays = cms.vstring( - 'MDME(174,1)=0 !Z decay into d dbar', - 'MDME(175,1)=0 !Z decay into u ubar', - 'MDME(176,1)=0 !Z decay into s sbar', - 'MDME(177,1)=0 !Z decay into c cbar', - # set it to 4 for the 1st Z to go to b bbar... - 'MDME(178,1)=0 !Z decay into b bbar', - # ............................................ - 'MDME(179,1)=0 !Z decay into t tbar', - 'MDME(182,1)=0 !Z decay into e- e+', - 'MDME(183,1)=0 !Z decay into nu_e nu_ebar', - 'MDME(184,1)=4 !Z decay into mu- mu+', - 'MDME(185,1)=0 !Z decay into nu_mu nu_mubar', - # ... and this one to 5, for the 2nd Z to goe to tau- tau+ - 'MDME(186,1)=5 !Z decay into tau- tau+', - # ............................................ - 'MDME(187,1)=0 !Z decay into nu_tau nu_taubar', - 'MDME(210,1)=0 !Higgs decay into dd', - 'MDME(211,1)=0 !Higgs decay into uu', - 'MDME(212,1)=0 !Higgs decay into ss', - 'MDME(213,1)=0 !Higgs decay into cc', - 'MDME(214,1)=0 !Higgs decay into bb', - 'MDME(215,1)=0 !Higgs decay into tt', - 'MDME(216,1)=0 !Higgs decay into', - 'MDME(217,1)=0 !Higgs decay into Higgs decay', - 'MDME(218,1)=0 !Higgs decay into e nu e', - 'MDME(219,1)=0 !Higgs decay into mu nu mu', - 'MDME(220,1)=0 !Higgs decay into tau nu tau', - 'MDME(221,1)=0 !Higgs decay into Higgs decay', - 'MDME(222,1)=0 !Higgs decay into g g', - 'MDME(223,1)=0 !Higgs decay into gam gam', - 'MDME(224,1)=0 !Higgs decay into gam Z', - 'MDME(225,1)=1 !Higgs decay into Z Z', - 'MDME(226,1)=0 !Higgs decay into W W' - ), - pythiaTauL = cms.vstring( - "mdme(89,1)=1", # tau -> e - "mdme(90,1)=1", # tau -> mu - # all other tau decays OFF - "mdme(91,1)=0", - "mdme(92,1)=0", - "mdme(93,1)=0", - "mdme(94,1)=0", - "mdme(95,1)=0", - "mdme(96,1)=0", - "mdme(97,1)=0", - "mdme(98,1)=0", - "mdme(99,1)=0", - "mdme(100,1)=0", - "mdme(101,1)=0", - "mdme(102,1)=0", - "mdme(103,1)=0", - "mdme(104,1)=0", - "mdme(105,1)=0", - "mdme(106,1)=0", - "mdme(107,1)=0", - "mdme(108,1)=0", - "mdme(109,1)=0", - "mdme(110,1)=0", - "mdme(111,1)=0", - "mdme(112,1)=0", - "mdme(113,1)=0", - "mdme(114,1)=0", - "mdme(115,1)=0", - "mdme(116,1)=0", - "mdme(117,1)=0", - "mdme(118,1)=0", - "mdme(119,1)=0", - "mdme(120,1)=0", - "mdme(121,1)=0", - "mdme(122,1)=0", - "mdme(123,1)=0", - "mdme(124,1)=0", - "mdme(125,1)=0", - "mdme(126,1)=0", - "mdme(127,1)=0", - "mdme(128,1)=0", - "mdme(129,1)=0", - "mdme(130,1)=0", - "mdme(131,1)=0", - "mdme(132,1)=0", - "mdme(133,1)=0", - "mdme(134,1)=0", - "mdme(135,1)=0", - "mdme(136,1)=0", - "mdme(137,1)=0", - "mdme(138,1)=0", - "mdme(139,1)=0", - "mdme(140,1)=0", - "mdme(141,1)=0", - "mdme(142,1)=0" - ), - # This is a vector of ParameterSet names to be read, in this order - parameterSets = cms.vstring( - 'pythiaSpecialSettings', - ##'pythiaSimpleSettings', - 'pythiaHZZDecays', 'pythiaTauL' ) - ) -) - -process.GEN = cms.OutputModule("PoolOutputModule", - fileName = cms.untracked.string('TestHZZleptons.root') -) - -process.p = cms.Path(process.generator) -process.outpath = cms.EndPath(process.GEN) - -process.schedule = cms.Schedule(process.p, process.outpath) diff --git a/GeneratorInterface/ExternalDecays/test/Py8GenFilter_Photos_cfg.py b/GeneratorInterface/ExternalDecays/test/Py8GenFilter_Photos_cfg.py deleted file mode 100644 index 4aab273830836..0000000000000 --- a/GeneratorInterface/ExternalDecays/test/Py8GenFilter_Photos_cfg.py +++ /dev/null @@ -1,87 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -process = cms.Process("TEST") - -process.load("FWCore.Framework.test.cmsExceptionsFatal_cff") -process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") - -process.source = cms.Source("EmptySource") - -process.generator = cms.EDFilter("Pythia8GeneratorFilter", - maxEventsToPrint = cms.untracked.int32(10), - pythiaPylistVerbosity = cms.untracked.int32(1), - filterEfficiency = cms.untracked.double(1.0), - pythiaHepMCVerbosity = cms.untracked.bool(True), - comEnergy = cms.double(7000.), - - ExternalDecays = cms.PSet( - Photos = cms.untracked.PSet(), - parameterSets = cms.vstring( "Photos" ) - ), - - PythiaParameters = cms.PSet( - py8SpecialSettings = cms.vstring( 'HadronLevel:Hadronize = off', # mstp(111)=0 - 'ParticleDecays:limitTau0 = on', # mstj(22)=2 - decay unstable particles - 'ParticleDecays:tau0Max = 10.', # parj(71)=10.- for which ctau < 10 mm - 'PartonLevel:ISR = off', # mstp(61)=0 - 'PartonLevel:FSR = off', # mstp(71)=0 ; what about mstj(41) ??? - 'PartonLevel:MI = off' # mstp(81)=0 - # for pp intractions, SpaceShower things are not important, - # but TimeShower settings are relevat - 'TimeShower:QEDshowerByL = off', - 'TimeShower:QEDshowerByQ = off' - ), - py8ProcessSettings = cms.vstring( 'StringZ:usePetersonB = on', # these 2 together == - 'StringZ:usePetersonC = on', # mstj(11)=3 - 'HiggsSM:gg2H = on', # mstp(102)=1 - 'HiggsSM:ff2Hff(t:ZZ) = on', # msub(123)=1 - ZZ fusion to H - 'HiggsSM:ff2Hff(t:WW) = on', # msub(124)=1 - WW fusion to H - '25:m0 = 190.0', # pmas(25,1)=190.0 - mass of H - 'PhaseSpace:pTHatMin = 5.', # ckin(45)=5. - 'PhaseSpace:pTHatMax = 150.' # ckin(46)=150. - # what about choice of structure function ??? (mstp(51)=7) - ), - py8HiggsDecaySettings = cms.vstring('25:onMode = off', # turn OFF all Higgs decays - '25:onIfAny = 23' # now turn ON H->ZZ - ), - py8ZDecaySettings = cms.vstring( '23:onMode = off', # turn OFF all Z decays - '23:onIfAny = 13', # turn ON Z->mumu - '23:onIfAny = 15' # turn ON Z->tautau - ), - py8TauDecaySettings = cms.vstring('15:onMode = off', # turn OFF all tau decays - '15:onIfAny = 11', # turn ON tau -> e - '15:onIfAny = 13' # turn ON tau -> mu - ), - parameterSets = cms.vstring('py8SpecialSettings','py8ProcessSettings', - 'py8HiggsDecaySettings', - 'py8ZDecaySettings','py8TauDecaySettings') - ) -) - -# The following three lines reduce the clutter of repeated printouts -# of the same exception message. -process.load("FWCore.MessageLogger.MessageLogger_cfi") - -process.MessageLogger.cerr.enableStatistics = False - - -process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", - generator = cms.PSet( - initialSeed = cms.untracked.uint32(123456789), - engineName = cms.untracked.string('HepJamesRandom') - ) -) - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(10) -) - -process.GEN = cms.OutputModule("PoolOutputModule", - fileName = cms.untracked.string('Py8ZHH_Photos.root') -) - -process.p = cms.Path(process.generator) -process.outpath = cms.EndPath(process.GEN) - -process.schedule = cms.Schedule(process.p, process.outpath) - diff --git a/GeneratorInterface/PhotosInterface/interface/PhotosInterface.h b/GeneratorInterface/PhotosInterface/interface/PhotosInterface.h deleted file mode 100644 index 0528e89245552..0000000000000 --- a/GeneratorInterface/PhotosInterface/interface/PhotosInterface.h +++ /dev/null @@ -1,53 +0,0 @@ -#ifndef gen_PhotosInterface_PhotosInterface_h -#define gen_PhotosInterface_PhotosInterface_h - -// #include "HepPDT/ParticleDataTable.hh" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" - -#include "HepMC/SimpleVector.h" -#include "GeneratorInterface/PhotosInterface/interface/PhotosInterfaceBase.h" - -namespace HepMC { - class GenEvent; - class GenVertex; -} // namespace HepMC - -namespace gen { - class PhotosInterface : public PhotosInterfaceBase { - public: - // ctor & dtor - PhotosInterface(); - PhotosInterface(const edm::ParameterSet&); - ~PhotosInterface() override; - - void init() override; - const std::vector& specialSettings() override { return fSpecialSettings; } - HepMC::GenEvent* apply(HepMC::GenEvent*) override; - void configureOnlyFor(int) override; - void avoidTauLeptonicDecays() override { - fAvoidTauLeptonicDecays = true; - return; - } - bool isTauLeptonicDecay(HepMC::GenVertex*); - void setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine) override; - static double flat(); - - private: - int fOnlyPDG; - bool fAvoidTauLeptonicDecays; - std::vector fBarcodes; - std::vector fSecVtxStore; - bool fIsInitialized; - - void applyToVertex(HepMC::GenEvent*, int); - void applyToBranch(HepMC::GenEvent*, int); - void attachParticles(HepMC::GenEvent*, HepMC::GenVertex*, int); - - static CLHEP::HepRandomEngine* fRandomEngine; - }; -} // namespace gen - -#endif diff --git a/GeneratorInterface/PhotosInterface/plugins/BuildFile.xml b/GeneratorInterface/PhotosInterface/plugins/BuildFile.xml index d45c12e79e311..7dcc4ceae0995 100644 --- a/GeneratorInterface/PhotosInterface/plugins/BuildFile.xml +++ b/GeneratorInterface/PhotosInterface/plugins/BuildFile.xml @@ -1,13 +1,3 @@ - - - - - - - - - - diff --git a/GeneratorInterface/PhotosInterface/plugins/PhotosInterface.cc b/GeneratorInterface/PhotosInterface/plugins/PhotosInterface.cc deleted file mode 100644 index 42f3429c9162b..0000000000000 --- a/GeneratorInterface/PhotosInterface/plugins/PhotosInterface.cc +++ /dev/null @@ -1,432 +0,0 @@ -#include -#include "GeneratorInterface/PhotosInterface/interface/PhotosInterface.h" -#include "FWCore/PluginManager/interface/ModuleDef.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h" - -#include "FWCore/ServiceRegistry/interface/Service.h" - -#include "HepMC/GenEvent.h" -#include "HepMC/IO_HEPEVT.h" -#include "HepMC/HEPEVT_Wrapper.h" - -using namespace gen; -using namespace edm; -using namespace std; - -CLHEP::HepRandomEngine* PhotosInterface::fRandomEngine = nullptr; - -extern "C" { - -void phoini_(void); -void photos_(int&); - -double phoran_(int* idummy) { return PhotosInterface::flat(); } - -extern struct { - // bool qedrad[NMXHEP]; - bool qedrad[4000]; // hardcoded for now... -} phoqed_; -} - -PhotosInterface::PhotosInterface() : fOnlyPDG(-1) { - fSpecialSettings.push_back("QED-brem-off:all"); - fAvoidTauLeptonicDecays = false; - fIsInitialized = false; -} - -PhotosInterface::PhotosInterface(const edm::ParameterSet&) : fOnlyPDG(-1) { - fSpecialSettings.push_back("QED-brem-off:all"); - fIsInitialized = false; -} - -PhotosInterface::~PhotosInterface() {} - -void PhotosInterface::setRandomEngine(CLHEP::HepRandomEngine* decayRandomEngine) { fRandomEngine = decayRandomEngine; } - -void PhotosInterface::configureOnlyFor(int ipdg) { - fOnlyPDG = ipdg; - // std::ostringstream command; - // command << "QED-brem-off:" << fOnlyPDG ; - fSpecialSettings.clear(); - // fSpecialSettings.push_back( command.str() ); - - return; -} - -void PhotosInterface::init() { - if (fIsInitialized) - return; // do init only once - - phoini_(); - - fIsInitialized = true; - - return; -} - -HepMC::GenEvent* PhotosInterface::apply(HepMC::GenEvent* evt) { - if (!fIsInitialized) - return evt; // conv.read_next_event(); - - // loop over HepMC::GenEvent, find vertices - - // for ( int ip=0; ipparticles_size(); ip++ ) - for (int ip = 0; ip < 4000; ip++) // 4000 is the max size of the array - { - phoqed_.qedrad[ip] = true; - } - - // - // now do actual job - // - - for (int iv = 1; iv <= evt->vertices_size(); iv++) { - bool legalVtx = false; - - fSecVtxStore.clear(); - - HepMC::GenVertex* vtx = evt->barcode_to_vertex(-iv); - - if (vtx->particles_in_size() != 1) - continue; // more complex than we need - if (vtx->particles_out_size() <= 1) - continue; // no outcoming particles - - if ((*(vtx->particles_in_const_begin()))->pdg_id() == 111) - continue; // pi0 decay vtx - no point to try - - if (fOnlyPDG != 1 && (*(vtx->particles_in_const_begin()))->pdg_id() != fOnlyPDG) { - continue; - } else { - // requested for specific PDG ID only, typically tau (15) - // - // first check if a brem vertex, where outcoming are the same pdg id and a photon - // - bool same = false; - for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children); - pitr != vtx->particles_end(HepMC::children); - ++pitr) { - if ((*pitr)->pdg_id() == fOnlyPDG) { - same = true; - break; - } - } - if (same) - continue; - - // OK, we get here if incoming fOnlyPDG and something else outcoming - // call it for the whole branch starting at vtx barcode iv, and go on - // NOTE: theoretically, it has a danger of double counting in vertices - // down the decay branch originating from fOnlyPDG, but in practice - // it's unlikely that down the branchg there'll be more fOnlyPDG's - - // cross-check printout - // vtx->print(); - - // overprotection... - // - if (fOnlyPDG == 15 && fAvoidTauLeptonicDecays && isTauLeptonicDecay(vtx)) - continue; - - applyToBranch(evt, -iv); - continue; - } - - // configured for all types of particles - // - for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children); - pitr != vtx->particles_end(HepMC::children); - ++pitr) { - // quark or gluon out of this vertex - no good ! - if (abs((*pitr)->pdg_id()) >= 1 && abs((*pitr)->pdg_id()) <= 8) - break; - if (abs((*pitr)->pdg_id()) == 21) - break; - - if ((*pitr)->status() == 1 || (*pitr)->end_vertex()) { - // OK, legal already ! - legalVtx = true; - break; - } - } - - if (!legalVtx) - continue; - - applyToVertex(evt, -iv); - - } // end of master loop - - // restore event number in HEPEVT (safety measure, somehow needed by Hw6) - HepMC::HEPEVT_Wrapper::set_event_number(evt->event_number()); - - return evt; -} - -void PhotosInterface::applyToVertex(HepMC::GenEvent* evt, int vtxbcode) { - HepMC::GenVertex* vtx = evt->barcode_to_vertex(vtxbcode); - - if (fAvoidTauLeptonicDecays && isTauLeptonicDecay(vtx)) - return; - - // cross-check printout - // - // vtx->print(); - - // first, flush out HEPEVT & tmp barcode storage - // - HepMC::HEPEVT_Wrapper::zero_everything(); - fBarcodes.clear(); - - // add incoming particle - // - int index = 1; - HepMC::HEPEVT_Wrapper::set_id(index, (*(vtx->particles_in_const_begin()))->pdg_id()); - HepMC::FourVector vec4; - vec4 = (*(vtx->particles_in_const_begin()))->momentum(); - HepMC::HEPEVT_Wrapper::set_momentum(index, vec4.x(), vec4.y(), vec4.z(), vec4.e()); - HepMC::HEPEVT_Wrapper::set_mass(index, (*(vtx->particles_in_const_begin()))->generated_mass()); - HepMC::HEPEVT_Wrapper::set_position( - index, vtx->position().x(), vtx->position().y(), vtx->position().z(), vtx->position().t()); - HepMC::HEPEVT_Wrapper::set_status(index, (*(vtx->particles_in_const_begin()))->status()); - HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0); - fBarcodes.push_back((*(vtx->particles_in_const_begin()))->barcode()); - - // add outcoming particles (decay products) - // - int lastDau = 1; - for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children); - pitr != vtx->particles_end(HepMC::children); - ++pitr) { - if ((*pitr)->status() == 1 || (*pitr)->end_vertex()) { - index++; - vec4 = (*pitr)->momentum(); - HepMC::HEPEVT_Wrapper::set_id(index, (*pitr)->pdg_id()); - HepMC::HEPEVT_Wrapper::set_momentum(index, vec4.x(), vec4.y(), vec4.z(), vec4.e()); - HepMC::HEPEVT_Wrapper::set_mass(index, (*pitr)->generated_mass()); - vec4 = (*pitr)->production_vertex()->position(); - HepMC::HEPEVT_Wrapper::set_position(index, vec4.x(), vec4.y(), vec4.z(), vec4.t()); - HepMC::HEPEVT_Wrapper::set_status(index, (*pitr)->status()); - HepMC::HEPEVT_Wrapper::set_parents(index, 1, 1); - fBarcodes.push_back((*pitr)->barcode()); - lastDau++; - } - if ((*pitr)->end_vertex()) { - fSecVtxStore.push_back((*pitr)->end_vertex()->barcode()); - } - } - - // store, further to set NHEP in HEPEVT - // - int nentries = index; - - // reset master pointer to mother - index = 1; - HepMC::HEPEVT_Wrapper::set_children(index, 2, lastDau); // FIXME: need to check - // if last daughter>=2 !!! - - // finally, set number of entries (NHEP) in HEPEVT - // - HepMC::HEPEVT_Wrapper::set_number_entries(nentries); - - // cross-check printout HEPEVT - // HepMC::HEPEVT_Wrapper::print_hepevt(); - - // OK, 1-level vertex is formed - now, call PHOTOS - // - photos_(index); - - // another cross-check printout HEPEVT - after photos - // HepMC::HEPEVT_Wrapper::print_hepevt(); - - // now check if something has been generated - // and make all adjustments to underlying vtx/parts - // - attachParticles(evt, vtx, nentries); - - // ugh, done with this vertex ! - - return; -} - -void PhotosInterface::applyToBranch(HepMC::GenEvent* evt, int vtxbcode) { - fSecVtxStore.clear(); - - // 1st level vertex - // - applyToVertex(evt, vtxbcode); - - // now look down the branch for more vertices, if any - // - // Note: fSecVtxStore gets filled up in applyToVertex, if necessary - // - unsigned int vcounter = 0; - - while (vcounter < fSecVtxStore.size()) { - applyToVertex(evt, fSecVtxStore[vcounter]); - vcounter++; - } - - fSecVtxStore.clear(); - - return; -} - -void PhotosInterface::attachParticles(HepMC::GenEvent* evt, HepMC::GenVertex* vtx, int nentries) { - if (HepMC::HEPEVT_Wrapper::number_entries() > nentries) { - // yes, need all corrections and adjustments - - // figure out how many photons and what particles in - // the decay branch have changes; - // also, follow up each one and correct accordingly; - // at the same time, add photon(s) to the GenVertex - // - - // vtx->print(); - - int largestBarcode = -1; - int Nbcodes = fBarcodes.size(); - - for (int ip = 1; ip < Nbcodes; ip++) { - int bcode = fBarcodes[ip]; - HepMC::GenParticle* prt = evt->barcode_to_particle(bcode); - if (bcode > largestBarcode) - largestBarcode = bcode; - double px = HepMC::HEPEVT_Wrapper::px(ip + 1); - double py = HepMC::HEPEVT_Wrapper::py(ip + 1); - double pz = HepMC::HEPEVT_Wrapper::pz(ip + 1); - double e = HepMC::HEPEVT_Wrapper::e(ip + 1); - double m = HepMC::HEPEVT_Wrapper::m(ip + 1); - - if (prt->end_vertex()) { - HepMC::GenVertex* endVtx = prt->end_vertex(); - - std::vector secVtxStorage; - secVtxStorage.clear(); - - secVtxStorage.push_back(endVtx->barcode()); - - const HepMC::FourVector& mom4 = prt->momentum(); - - // now rescale all descendants - double bet1[3], bet2[3], gam1, gam2, pb; - double mass = mom4.m(); - bet1[0] = -(mom4.px() / mass); - bet1[1] = -(mom4.py() / mass); - bet1[2] = -(mom4.pz() / mass); - bet2[0] = px / m; - bet2[1] = py / m; - bet2[2] = pz / m; - gam1 = mom4.e() / mass; - gam2 = e / m; - - unsigned int vcounter = 0; - - while (vcounter < secVtxStorage.size()) { - HepMC::GenVertex* theVtx = evt->barcode_to_vertex(secVtxStorage[vcounter]); - - for (HepMC::GenVertex::particle_iterator pitr = theVtx->particles_begin(HepMC::children); - pitr != theVtx->particles_end(HepMC::children); - ++pitr) { - if ((*pitr)->end_vertex()) { - secVtxStorage.push_back((*pitr)->end_vertex()->barcode()); - } - - if (theVtx->particles_out_size() == 1 && (*pitr)->pdg_id() == prt->pdg_id()) { - // carbon copy - (*pitr)->set_momentum(HepMC::FourVector(px, py, pz, e)); - continue; - } - - HepMC::FourVector dmom4 = (*pitr)->momentum(); - - // Boost vector to parent rest frame... - pb = bet1[0] * dmom4.px() + bet1[1] * dmom4.py() + bet1[2] * dmom4.pz(); - double dpx = dmom4.px() + bet1[0] * (dmom4.e() + pb / (gam1 + 1.)); - double dpy = dmom4.py() + bet1[1] * (dmom4.e() + pb / (gam1 + 1.)); - double dpz = dmom4.pz() + bet1[2] * (dmom4.e() + pb / (gam1 + 1.)); - double de = gam1 * dmom4.e() + pb; - // ...and boost back to modified parent frame - pb = bet2[0] * dpx + bet2[1] * dpy + bet2[2] * dpz; - dpx += bet2[0] * (de + pb / (gam2 + 1.)); - dpy += bet2[1] * (de + pb / (gam2 + 1.)); - dpz += bet2[2] * (de + pb / (gam2 + 1.)); - de *= gam2; - de += pb; - - (*pitr)->set_momentum(HepMC::FourVector(dpx, dpy, dpz, de)); - } - vcounter++; - } - - secVtxStorage.clear(); - } - - prt->set_momentum(HepMC::FourVector(px, py, pz, e)); - - } // ok, all affected particles update, but the photon(s) still not inserted - - int newlyGen = HepMC::HEPEVT_Wrapper::number_entries() - nentries; - - if (largestBarcode < evt->particles_size()) { - // need to adjust barcodes down from the affected vertex/particles - // such that we "free up" barcodes for newly generated photons - // in the middle of the event record - // - for (int ipp = evt->particles_size(); ipp > largestBarcode; ipp--) { - (evt->barcode_to_particle(ipp))->suggest_barcode(ipp + newlyGen); - } - } - - // now attach new generated photons to the vertex - // - for (int ipnw = 1; ipnw <= newlyGen; ipnw++) { - int nbcode = largestBarcode + ipnw; - int pdg_id = HepMC::HEPEVT_Wrapper::id(nentries + ipnw); - int status = HepMC::HEPEVT_Wrapper::status(nentries + ipnw); - double px = HepMC::HEPEVT_Wrapper::px(nentries + ipnw); - double py = HepMC::HEPEVT_Wrapper::py(nentries + ipnw); - double pz = HepMC::HEPEVT_Wrapper::pz(nentries + ipnw); - double e = HepMC::HEPEVT_Wrapper::e(nentries + ipnw); - double m = HepMC::HEPEVT_Wrapper::m(nentries + ipnw); - - HepMC::GenParticle* NewPart = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, e), pdg_id, status); - NewPart->set_generated_mass(m); - NewPart->suggest_barcode(nbcode); - vtx->add_particle_out(NewPart); - } - - //vtx->print(); - //std::cout << " leaving attachParticles() " << std::endl; - - } // end of global if-statement - - return; -} - -bool PhotosInterface::isTauLeptonicDecay(HepMC::GenVertex* vtx) { - if (abs((*(vtx->particles_in_const_begin()))->pdg_id()) != 15) - return false; - - for (HepMC::GenVertex::particle_iterator pitr = vtx->particles_begin(HepMC::children); - pitr != vtx->particles_end(HepMC::children); - ++pitr) { - if (abs((*pitr)->pdg_id()) == 11 || abs((*pitr)->pdg_id()) == 13) { - return true; - } - } - - return false; -} - -double PhotosInterface::flat() { - if (!fRandomEngine) { - throw cms::Exception("LogicError") - << "PhotosInterface::flat: Attempt to generate random number when engine pointer is null\n" - << "This might mean that the code was modified to generate a random number outside the\n" - << "event and beginLuminosityBlock methods, which is not allowed.\n"; - } - return fRandomEngine->flat(); -} - -DEFINE_EDM_PLUGIN(PhotosFactory, gen::PhotosInterface, "Photos2155");