Commit 1151bc36 authored by William Patrick Mccormack's avatar William Patrick Mccormack Committed by Ewelina Maria Lobodzinska
Browse files

New Superchic Interface

parent a19e6c61
......@@ -16,7 +16,7 @@ genseeds = {
'Pythia8B' : ["PYTHIA8 OFFSET {rnd} 4789899 {rnd}", "PYTHIA8_INIT 820021 2347532"],
'Herwig' : ["HERWIG OFFSET {rnd} 30450275 {rnd}", "HERWIG_INIT 620021 5347532"],
'Herwigpp' : ["Herwigpp OFFSET {rnd} {rnd} 31122001"],
'Herwig7' : ["Herwig7 OFFSET {rnd} {rnd} 31122001"],
'Herwig7' : ["Herwig7 OFFSET {rnd} {rnd} 31122001"],
'Jimmy' : ["JIMMY OFFSET {rnd} 39002061 {rnd}", "JIMMY_INIT 720021 4347532"],
'Cascade' : ["CASCADE OFFSET {rnd} 4789899 {rnd}", "CASCADE_INIT 889223465 78782321"],
'Tauola' : ["TAUOLA OFFSET {rnd} 10480275 {rnd}", "TAUOLA_INIT 920021 3347532"],
......@@ -40,7 +40,8 @@ genseeds = {
'Exhume' : ["ExhumeRand OFFSET {rnd} 4475757 {rnd}"],
'Pomwig' : ["POMWIG OFFSET {rnd} 37489241 {rnd}", "POMWIG_INIT 21219421 1984121"],
'Starlight' : ["STARLIGHT OFFSET {rnd} {rnd} 31122001"],
'BeamHaloGenerator' : ["BeamHalo OFFSET {rnd} 3524752 {rnd}"]
'BeamHaloGenerator' : ["BeamHalo OFFSET {rnd} 3524752 {rnd}"],
'Superchic' : ["SUPERCHIC OFFSET {rnd} {rnd} 10480275"]
}
## Decide whether to use the RanLux or standard random number service
......
......@@ -6,50 +6,15 @@
# Declare the package name:
atlas_subdir( Superchic_i )
# Declare the package's dependencies:
atlas_depends_on_subdirs(
PUBLIC
Generators/GeneratorFortranCommon
Generators/GeneratorModules
PRIVATE
Control/AthenaKernel
Control/StoreGate
GaudiKernel
Generators/GeneratorUtils
Generators/TruthUtils
Projects/AthGeneration
Projects/Athena
)
# External dependencies:
find_package( CLHEP )
find_package( HepMC COMPONENTS HepMC HepMCfio )
find_package( Apfel )
find_package( Lhapdf )
find_package( Superchic )
atlas_disable_as_needed()
# Component(s) in the package:
atlas_add_library( Superchic_iLib
Superchic_i/*.h src/*.cxx src/*.F
PUBLIC_HEADERS Superchic_i
INCLUDE_DIRS ${SUPERCHIC_INCLUDE_DIRS}
PRIVATE_INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
PRIVATE_DEFINITIONS ${CLHEP_DEFINITIONS}
LINK_LIBRARIES ${SUPERCHIC_LIBRARIES} GeneratorModulesLib
PRIVATE_LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaKernel
GaudiKernel TruthUtils )
atlas_add_component( Superchic_i
src/components/*.cxx
INCLUDE_DIRS ${SUPERCHIC_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS}
LINK_LIBRARIES ${SUPERCHIC_LIBRARIES} ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} GeneratorFortranCommonLib GeneratorModulesLib AthenaKernel StoreGateLib SGtests GaudiKernel GeneratorObjects TruthUtils Superchic_iLib )
# Install files from the package:
atlas_install_python_modules( python/*.py )
atlas_install_joboptions( share/*.py )
atlas_add_component( Superchic_i
src/components/*.cxx
)
# Set Superchic specific environment variable(s).
set( SuperchicEnvironment_DIR ${CMAKE_CURRENT_SOURCE_DIR}
CACHE PATH "Location of SuperchicEnvironment.cmake" )
find_package( SuperchicEnvironment )
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
Copied directly from the Tauola_i/Atlas_HEPEVT with minor changes
*/
#ifndef Atlas_HEPEVT_modified_h
#define Atlas_HEPEVT_modified_h
#include "CLIDSvc/CLASS_DEF.h"
extern "C" {
void* atlas_hepevt_modified_address_(void);
}
/**
@class Atlas_HEPEVT_modified.h
@brief Needed fot interface of Superchic_i to Superchic3.03 generator
to store the umodified HEPEVT common.
This code was directly copied from Tailola_i package https://gitlab.cern.ch/atlas/athena/blob/21.6/Generators/Tauola_i/Tauola_i/Atlas_HEPEVT.h
with minor changes, original author list is following. This modified class takes into account the length of the array in hepevt parameteri nmxhep which is set to (nmxhep=4000) instead of HEPEVT standard event common which uses nmxhep=10000
@author Borut Paul Kersevan (BPK), June 2003
HEPEVT standard event common
PARAMETER (NMXHEP=10000)
COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
& JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
*/
class Atlas_HEPEVT_modified {
public:
Atlas_HEPEVT_modified();
~Atlas_HEPEVT_modified();
void fill();
void spill();
int& nevhep();
int& nhep();
int& isthep(int nh);
int& idhep(int nh);
int& jmohep(int i, int nh);
int& jdahep(int i, int nh);
double& phep(int j, int nh);
double& vhep(int k, int nh);
inline void init(); // inlined for speed of access (small function)
// return common array lengths
int depthRel() const {return s_depthRel;}
int depthPhep() const {return s_depthPhep;}
int depthVhep() const {return s_depthVhep;}
int lenNmxhep() const {return s_lenNmxhep;}
private:
// Lengths of array in HEPEVT common
static const int s_depthRel = 2 ; // Relatives -Mother & Daughter
static const int s_depthPhep = 5 ;
static const int s_depthVhep = 4 ;
static const int s_lenNmxhep = 4000 ;
struct HEPEVT;
friend struct HEPEVT;
struct HEPEVT {
int nevhep;
int nhep;
int isthep[s_lenNmxhep];
int idhep[s_lenNmxhep];
int jmohep[s_lenNmxhep][s_depthRel];
int jdahep[s_lenNmxhep][s_depthRel];
double phep[s_lenNmxhep][s_depthPhep];
double vhep[s_lenNmxhep][s_depthVhep];
};
int m_dummy;
double m_realdummy;
// s_HEPEVT is needed access the umodified HEPEVT common block information
static HEPEVT* s_HEPEVT;
HEPEVT m_HEPEVT;
};
CLASS_DEF( Atlas_HEPEVT_modified, 9966, 1)
#include "Superchic_i/Atlas_HEPEVT_modified.icc"
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// Inline implementations for Atlas_HEPEVT_modified
// initialise pointer
void Atlas_HEPEVT_modified::init(void) {
if (!s_HEPEVT) s_HEPEVT = static_cast<HEPEVT*>(atlas_hepevt_modified_address_());
}
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// --------------------------------------------------
//
// File: Generators/Superchic_i.h
// Description:
// This code is used to get a Starlight Monte Carlo event.
// genInitialize() is used to read parameters
// callGenerator() makes the event
// genFinalize() writes log files etc
// fillEvt(GeneratorEvent* evt) passes the event to HepMC
//
// The output will be stored in the transient event store so it can be
// passed to the simulation.
//
// AuthorList:
// Prabhakar Palni, 29 August 2018
#ifndef GENERATORMODULESSUPERCHIC_H
#define GENERATORMODULESSUPERCHIC_H
#include "GeneratorModules/GenModule.h"
#include "Superchic_i/Atlas_HEPEVT_modified.h"
using std::string;
typedef std::vector<std::string> CommandVector;
class StoreGateSvc;
// ---------------------
class Superchic_i:public GenModule {
public:
Superchic_i(const std::string& name, ISvcLocator* pSvcLocator);
virtual ~Superchic_i();
virtual StatusCode genInitialize ();
virtual StatusCode callGenerator ();
virtual StatusCode genFinalize ();
virtual StatusCode fillEvt (HepMC::GenEvent* evt);
protected:
// event counter
int m_events ;
// inputs parameter (SuperChic initialization)
std::string m_configFileName;
std::string m_rts;
int m_isurv;
std::string m_intag;
std::string m_PDFname;
unsigned int m_PDFmember;
//pdfinfo
HepMC::PdfInfo m_pdfinfo;
// inputs parameter (events generation)
unsigned int m_proc;
std::string m_beam;
std::string m_outtg;
std::string m_sfaci;
std::string m_diff;
std::string m_an;
std::string m_az;
std::string m_rz;
std::string m_dz;
std::string m_rn;
std::string m_dn;
std::string m_ionqcd;
unsigned int m_ncall;
unsigned int m_itmx;
std::string m_prec;
unsigned int m_ncall1;
unsigned int m_inccall;
unsigned int m_itend;
int m_iseed;
//int m_s2int;
std::string m_genunw;
int m_nev;
std::string m_erec;
std::string m_readwt;
std::string m_wtmax;
std::string m_ymin;
std::string m_ymax;
std::string m_mmin;
std::string m_mmax;
std::string m_gencuts;
std::string m_scorr;
std::string m_fwidth;
std::string m_ptxmax;
std::string m_acoabmax;
std::string m_ptamin;
std::string m_ptbmin;
std::string m_ptcmin;
std::string m_ptdmin;
std::string m_ptemin;
std::string m_ptfmin;
std::string m_etaamin;
std::string m_etabmin;
std::string m_etacmin;
std::string m_etadmin;
std::string m_etaemin;
std::string m_etafmin;
std::string m_etaamax;
std::string m_etabmax;
std::string m_etacmax;
std::string m_etadmax;
std::string m_etaemax;
std::string m_etafmax;
std::string m_rjet;
std::string m_jalg;
std::string m_m2b;
int m_pdgid1;
int m_pdgid2;
std::string m_malp;
std::string m_gax;
std::string m_alpt;
std::string m_mpol;
std::string m_mmon;
std::string m_gamm;
std::string m_mcharg;
std::string m_mneut;
// Random numbers seed
std::vector<long int> m_seeds;
// Commands to setup superchic
CommandVector m_InitializeVector;
bool set_user_params();
bool prepare_params_file();
// I/O to HEPEVT
static Atlas_HEPEVT_modified* s_atlas_HEPEVT;
};
#endif
from AthenaCommon.AppMgr import ServiceMgr as svcMgr
from GeneratorModules.EvgenAlg import EvgenAlg
from AthenaPython.PyAthena import StatusCode
from AthenaPython.PyAthena import HepMC, StatusCode
import McParticleEvent.Pythonizations
import os
import ROOT
import random
class LheEVNTFiller(EvgenAlg):
def __init__(self, name="LheEVNTFiller"):
super(LheEVNTFiller, self).__init__(name=name)
fileName = "evrecs/evrecout.dat"
eventsProcessed = 0
def initialize(self):
seed = None
if(os.path.isfile(self.fileName)):
print(self.fileName)
return StatusCode.Success
else:
return StatusCode.Failure
def fillEvent(self, evt):
eventsSeen = 0
firstLine = True
particlesSeen = 0
with open(self.fileName,'r') as inputfile:
event = False
for line in inputfile:
if not event and not '<event>' in line:
continue
# Check if we are just starting an event
if not event and '<event>' in line and eventsSeen == self.eventsProcessed:
event = True
continue
# Check if we have finished and are doing something else
if '<' in line or '>' in line:
event = False
eventsSeen += 1
firstLine = True
continue
if event and firstLine:
firstLine = False
evt.weights().push_back(float(line.split()[2]))
continue
if event:
pos = HepMC.FourVector(0.0, 0.0, 0.0, 0.0)
gv = HepMC.GenVertex(pos)
ROOT.SetOwnership(gv, False)
evt.add_vertex(gv)
mom = HepMC.FourVector( float(line.split()[6])*1000. , float(line.split()[7])*1000. , float(line.split()[8])*1000. , float(line.split()[9])*1000. )
gp = HepMC.GenParticle()
gp.set_status(int(line.split()[1]) )
gp.set_pdg_id(int(line.split()[0]) )
gp.set_momentum(mom)
print(float(line.split()[10]) * 1000.)
gp.set_generated_mass(float(line.split()[10]) * 1000.)
ROOT.SetOwnership(gp, False)
gv.add_particle_out(gp)
self.eventsProcessed += 1
return StatusCode.Success
import subprocess, os
class SuperChicConfig:
superchicpath = os.environ['SUPERCHICPATH']
#SuperChic specific variables for input.DAT, see writeInputDAT function for more elaboration
rts = "13d3" #collision energy (GeV)
isurv = "4" #Model of soft survival (from 1 -> 4, corresponding to arXiv:1306.2149)
intag = "'in13'" #tag for input files
PDFname = "'MMHT2014lo68cl'"
PDFmember = "0"
proc = "57" #Please consult Superchic Manual https://superchic.hepforge.org/
beam = "'prot'"
outtg = "'out'"
sfaci = ".true."
diff = "'el'" #interaction: elastic ('el'), single ('sd','sda','sdb') and double ('dd') dissociation.
an = "208d0"
az = "82d0"
rz = "6.68d0"
dz = "0.447d0"
rn = "6.7d0"
dn = "0.55d0"
ionqcd = "'coh'"
ncall = "10000"
itmx = "10"
prec = "0.5d0"
ncall1 = "10000"
inccall = "10000"
itend = "1000"
iseed = "34"
genunw = ".true"
nev = "500"
erec = "'lhe'" #output file type
readwt = ".false."
wtmax = "0d0"
ymin = "-2.4d0" # Minimum object rapidity Y_X
ymax = "2.4d0" # Maximum object rapidity Y_X
mmin = "6d0" # Minimum object mass M_X
mmax = "500d0" # Maximum object mass M_X
gencuts = ".true." # Generate cuts below
scorr = ".true."
fwidth = ".true."
ptxmax = "100d0"
ptamin = "3.0d0" # Minimum pT of outgoing object a (gamma)
ptbmin = "3.0d0" # Minimum pT of outgoing object b (gamma)
etaamin = "-2.4d0" # Minimum eta of outgoing object a
etaamax = "2.4d0" # Maximum eta of outgoing object a
etabmin = "-2.4d0" # Minimum eta of outgoing object b
etabmax = "2.4d0" # Maximum eta of outgoing object b
acoabmax = "100d0"
ptcmin = "0d0"
etacmin = "-2.5d0"
etacmax = "2.5d0"
ptdmin = "0d0"
etadmin = "-2.5d0"
etadmax = "2.5d0"
ptemin = "0d0"
etaemin = "-2.5d0"
etaemax = "2.5d0"
ptfmin = "0d0"
etafmin = "-2.5d0"
etafmax = "2.5d0"
rjet = "0.6d0"
jalg = "'antikt'"
m2b = "0.133d0"
pdgid1 = "211"
pdgid2 = "-211"
malp = "1000d0"
gax = "0.001d0"
alpt = "'ps'"
mpol = "500d0"
mmon = "933d0"
gamm = "10d0"
mcharg = "100d0"
mneut = "80d0"
def writeInputDAT(Init):
outF = open("input.DAT", "w")
outF.write("***********************************************************************************\n")
outF.write("****** Initialize afain if FIRST FIVE PARAMETERS ARE CHANGED (and beam = 'prot'):******\n")
outF.write("***********************************************************************************\n")
outF.write(Init.rts + " ! [rts] : CMS collision energy (GeV) \n")
outF.write(Init.isurv + " ! [isurv] : Model of soft survival (from 1 -> 4)\n")
outF.write(Init.intag + " ! [intag] for input files \n")
outF.write("***********************************************************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.PDFname + " ! [PDFname] : PDF set \n")
outF.write(Init.PDFmember + " ! [PDFmember] : PDF member \n")
outF.write("***********************************************************************************\n")
outF.write(Init.proc + " ! [proc] : Process number \n")
outF.write(Init.beam + " ! [beam] : Beam type ('prot', 'ion', 'ionp' or 'el') \n")
outF.write(Init.outtg + " ! [outtg] : for output file \n")
outF.write(Init.sfaci + " ! [sfaci] : Include soft survival effects \n")
outF.write("***********************************************************************************\n")
outF.write("***********************************************************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.diff + " ! [diff] : dissociation flag \n")
outF.write("***********************************************************************************\n")
outF.write(Init.an + " ! [an] : Ion mass number \n")
outF.write(Init.az + " ! [az] : Ion atomic number \n")
outF.write(Init.rz + " ! [rz] : Ion proton density - radius \n")
outF.write(Init.dz + " ! [dz] : Ion proton density - skin thickness \n")
outF.write(Init.rn + " ! [rn] : Ion neutron density - radius \n")
outF.write(Init.dn + " ! [dn] : Ion neutron density - skin thickness \n")
outF.write(Init.ionqcd + " ! [ionqcd] : Coherent ('coh') or incoherent ('incoh') for QCD-induced processes \n")
outF.write("***********************************************************************************\n")
outF.write("*************Integration parameters************************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.ncall + " ! [ncall] : Number of calls for preconditioning \n")
outF.write(Init.itmx + " ! [itmx] : Number of iterations for preconditioning \n")
outF.write(Init.prec + " ! [prec] : Relative accuracy (in %) in main run \n")
outF.write(Init.ncall1 + " ! [ncall1] : Number of calls in first iteration \n")
outF.write(Init.inccall + " ! [inccall] : Number of increase calls per iteration \n")
outF.write(Init.itend + " ! [itend] : Maximum number of iterations \n")
outF.write(Init.iseed + " ! [iseed] : Random number seed (integer > 0) \n")
outF.write("***********************************************************************************\n")
outF.write("********************Unweighted events**********************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.genunw + " ! [genunw] : Generate unweighted events \n")
outF.write(Init.nev + " ! [nev] : Number of events (1 or a small value is recommended) ( the number of produced total events are controlled by athena ) \n")
outF.write(Init.erec + " ! [erec] : Event record format ('hepmc','lhe','hepevt') \n")
outF.write(Init.readwt + " ! [readwt] : Set to true to read in pre-calculated maxium weight below \n")
outF.write(Init.wtmax + " ! [wtmax] : Maximum weight \n")
outF.write("***********************************************************************************\n")
outF.write("***********************************************************************************\n")
outF.write("******************* general cuts ************************************************\n")
outF.write(Init.ymin + " ! [ymin] : Minimum object rapidity Y_X \n")
outF.write(Init.ymax + " ! [ymax] : Maximum object rapidity Y_X \n")
outF.write(Init.mmin + " ! [mmin] : Minimum object mass M_X (redundant for resonance production) \n")
outF.write(Init.mmax + " ! [mmax] : Maximum object mass M_X (redundant for resonance production) \n")
outF.write(Init.gencuts + " ! [gencuts] : Generate cuts below \n")
outF.write(Init.scorr + " ! [scorr] : Include spin correlations (for chi_c/b decays) \n")
outF.write(Init.fwidth + " ! [fwidth] : Include finite width (for chi_c decays) \n")
outF.write("***********************************************************************************\n")
outF.write("************ See manual for momentum assignments***********************************\n")
outF.write("***********************************************************************************\n")
outF.write("******************* Proton Cuts ***************************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.ptxmax + " ! [ptxmax] : max pT of the system \n")
outF.write("***********************************************************************************\n")
outF.write("**********2 body final states : p(a) + p(b) ***************************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init. ptamin + " ! [ptamin] \n")
outF.write(Init. ptbmin + " ! [ptbmin] \n")
outF.write(Init. etaamin + " ! [etaamin] \n")
outF.write(Init. etaamax + " ! [etaamax] \n")
outF.write(Init. etabmin + " ! [etabmin] \n")
outF.write(Init. etabmax + " ! [etabmax] \n")
outF.write(Init. acoabmax + " ! [acoabmax] \n")
outF.write("***********************************************************************************\n")
outF.write("****** 3 body final states : p(a) + p(b) + p(c) ***********************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.ptamin + " ! [ptamin] \n")
outF.write(Init.ptbmin + " ! [ptbmin] \n")
outF.write(Init.ptcmin + " ! [ptcmin] \n")
outF.write(Init.etaamin + " ! [etaamin] \n")
outF.write(Init.etaamax + " ! [etaamax] \n")
outF.write(Init.etabmin + " ! [etabmin] \n")
outF.write(Init.etabmax + " ! [etabmax] \n")
outF.write(Init.etacmin + " ! [etacmin] \n")
outF.write(Init.etacmax + " ! [etacmax] \n")
outF.write("***********************************************************************************\n")
outF.write("****** 4 body final states : p(a) + p(b) + p(c) + p(d) ****************************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.ptamin + " ! [ptamin] \n")
outF.write(Init.ptbmin + " ! [ptbmin] \n")
outF.write(Init.ptcmin + " ! [ptcmin] \n")
outF.write(Init.ptdmin + " ! [ptdmin] \n")
outF.write(Init.etaamin + " ! [etaamin] \n")
outF.write(Init.etaamax + " ! [etaamax] \n")
outF.write(Init.etabmin + " ! [etabmin] \n")
outF.write(Init.etabmax + " ! [etabmax] \n")
outF.write(Init.etacmin + " ! [etacmin] \n")
outF.write(Init.etacmax + " ! [etacmax] \n")
outF.write(Init.etadmin + " ! [etacmin] \n")
outF.write(Init.etadmax + " ! [etadmax] \n")
outF.write("***********************************************************************************\n")
outF.write("****** 6 body final states : p(a) + p(b) + p(c) + p(d) + p(e) + p(f) **************\n")
outF.write("***********************************************************************************\n")
outF.write(Init.ptamin + " ! [ptamin] \n")
outF.write(Init.ptbmin + " ! [ptbmin] \n")