Commit 1bd437c7 authored by Prabhakar Palni's avatar Prabhakar Palni
Browse files

Superchic Interface codesÃ

parent 1a494d70
......@@ -27,7 +27,7 @@ mainGenerators += ["Exhume", "Phojet", "Epos", "QGSJet"]
mainGenerators += ["ParticleGenerator", "ParticleGun"]
mainGenerators += ["CosmicGenerator", "BeamHaloGenerator"]
## Heavy ion generators
mainGenerators += ["Starlight", "Hijing", "Hydjet", "Reldis", "Pyquen"]
mainGenerators += ["Superchic","Starlight", "Hijing", "Hydjet", "Reldis", "Pyquen"]
## Misc generators
mainGenerators += ["AcerMC", "TopRex", "LPair"]
## Reading in fully-formed events
......@@ -42,7 +42,7 @@ afterburnerGenerators = ["Photos", "Photospp", "Tauola", "TauolaPP", "Tauolapp",
knownGenerators = inputGenerators + mainGenerators + afterburnerGenerators
## Note which generators should NOT be sanity tested by the TestHepMC alg
notesthepmcGenerators = ["ParticleDecayer", "ParticleGun", "CosmicGenerator", "BeamHaloGenerator", "FPMC",
notesthepmcGenerators = ["Superchic","ParticleDecayer", "ParticleGun", "CosmicGenerator", "BeamHaloGenerator", "FPMC",
"Hijing", "Hydjet", "Starlight", "PythiaRhad"]
## Generators with no flexibility/concept of a tune or PDF choice
......
# $Id: CMakeLists.txt 762209 2016-07-15 14:15:34Z krasznaa $
################################################################################
# Package: Superchic_i
################################################################################
# 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 )
## TODO Currently FindSuperchic.cmake do not find the superchic generator libraries as it is under atlasexternals
#include("FindSuperchic.cmake")
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_joboptions( share/*.py )
/*
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.
@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;
static HEPEVT* s_atlas_HEPEVT;
HEPEVT m_atlas_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_atlas_HEPEVT) s_atlas_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"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/GenEvent.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; // TODO check this ?
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
evgenConfig.description = "Superchic gamma + gamma UPC collisions at 5020 GeV"
evgenConfig.keywords = ["2photon","2lepton"]
#evgenConfig.weighting = 0
evgenConfig.contact = ["prabhakar.palni@cern.ch"]
if not os.path.exists('inputs'):
os.makedirs('inputs')
if not os.path.exists('evrecs'):
os.makedirs('evrecs')
from Superchic_i.Superchic_iConf import Superchic_i
genSeq += Superchic_i("Superchic")
genSeq.Superchic.McEventKey = "GEN_EVENT"
evgenConfig.generators += ["Superchic"]
from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
_evgenstream = AthenaPoolOutputStream("StreamEVGEN")
_evgenstream.ItemList = ["2101#*","133273#GEN_EVENT"]
del _evgenstream
# TODO: Sort out proper param setting based on runArgs.ecmEnergy
if int(runArgs.ecmEnergy) != 5020:
evgenLog.error(" Set beam energy in JO initialization with parameter rts ")
sys.exit(1)
genSeq.Superchic.Initialize = \
["rts 5.02d3", # set the COM collision energy (in fortran syntax)
"isurv 4", # Model of soft survival
"intag 'in5'", # for input files
"PDFname 'MMHT2014lo68cl'", # PDF set name
"PDFmember 0", # PDF member
"proc 59", # Process number (59 = gg->gg, 56 =gg->ee )
"beam 'ion'", # Beam type ('prot', 'ion')
"outtg 'out'", # for output file name
"sfaci .false.", # Include soft survival effects
"ncall 100000", # Number of calls for preconditioning
"itmx 1", # Number of iterations for preconditioning
"prec 1.0d0" # precision
]
****************************************************************************************
****** RE-RUN ./init IF FIRST FIVE PARAMETERS ARE CHANGED (and beam = 'prot'): *******
****************************************************************************************
5.02d3 ! [rts] : CMS collision energy (GeV)
4 ! [isurv] : Model of soft survival (from 1 -> 4, corresponding to arXiv:1306.2149)
'in5' ! [intag] for input files
****************************************************************************************
****************************************************************************************
'MMHT2014lo68cl' ! [PDFname] : PDF set
0 ! [PDFmember] : PDF member
****************************************************************************************
59 ! [proc] : Process number (see manual for labelling)
'ion' ! [beam] : Beam type ('prot', 'ion', 'ionp' or 'electron')
'out' ! [outtg] : for output file
.false. ! [sfaci] : Include soft survival effects
****************************************************************************************
208d0 ! [an] : Ion mass number
82d0 ! [az] : Ion atomic number
6.68d0 ! [rz] : Ion proton density - radius
0.447d0 ! [dz] : Ion proton density - skin thickness
6.7d0 ! [rn] : Ion neutron density - radius
0.55d0 ! [dn] : Ion neutron density - skin thickness
'coh' ! [ionqcd] : Coherent ('coh') or incoherent ('incoh') for QCD-induced processes
****************************************************************************************
************************* Integration parameters ***************************************
****************************************************************************************
10000 ! [ncall] : Number of calls for preconditioning
10 ! [itmx] : Number of iterations for preconditioning
1.0d0 ! [prec] : Relative accuracy (in %) in main run
10000 ! [ncall1] : Number of calls in first iteration
10000 ! [inccall] : Number of increase calls per iteration
1000 ! [itend] : Maximum number of iterations
23 ! [iseed] : Random number seed (integer > 0)
****************************************************************************************
8 ! [s2int] : Survival factor integration : see manual for recommendation
****************************************************************************************
******************************* Unweighted events **************************************
****************************************************************************************
.true. ! [genunw] : Generate unweighted events
1 ! [nev] : Number of events ( < 1000000 recommended)
'hepevt' ! [erec] : Event record format ('hepmc','lhe','hepevt')
.false. ! [readwt] : Set to true to read in pre-calculated maxium weight below
0d0 ! [wtmax] : Maximum weight
****************************************************************************************
******************************* general cuts **************************************
****************************************************************************************
-2.5d0 ! [ymin] : Minimum object rapidity Y_X
2.5d0 ! [ymax] : Maximum object rapidity Y_X
5d0 ! [mmin] : Minimum object mass M_X (redundant for resonance production)
50d0 ! [mmax] : Maximum object mass M_X (redundant for resonance production)
.true. ! [gencuts] : Generate cuts below
.true. ! [scorr] : Include spin correlations (for chi_c/b decays)
.true. ! [fwidth] : Include finite width (for chi_c decays)
****************************************************************************************
************************ See manual for momentum assignments ***************************
****************************************************************************************
*********************************** Proton Cuts ***************************************
****************************************************************************************
100d0 ! [ptxmax]
****************************************************************************************
************************* 2 body final states : p(a) + p(b) ****************************
****************************************************************************************
0d0 ! [ptamin]
0d0 ! [ptbmin]
-2.5d0 ! [etaamin]
2.5d0 ! [etaamax]
-2.5d0 ! [etabmin]
2.5d0 ! [etabmax]
10d0 ! [acoabmax]
****************************************************************************************
********************** 3 body final states : p(a) + p(b) + p(c) ************************
****************************************************************************************
20d0 ! [ptamin]
20d0 ! [ptbmin]
20d0 ! [ptcmin]
-2.4d0 ! [etaamin]
2.4d0 ! [etaamax]
-2.4d0 ! [etabmin]
2.4d0 ! [etabmax]
-2.4d0 ! [etacmin]
2.4d0 ! [etacmax]
****************************************************************************************
****************** 4 body final states : p(a) + p(b) + p(c)+ p(d) **********************
****************************************************************************************
20d0 ! [ptamin]
20d0 ! [ptbmin]
20d0 ! [ptcmin]
20d0 ! [ptdmin]
-2.4d0 ! [etaamin]
2.4d0 ! [etaamax]
-2.4d0 ! [etabmin]
2.4d0 ! [etabmax]
-2.4d0 ! [etacmin]
2.4d0 ! [etacmax]
-2.4d0 ! [etadmin]
2.4d0 ! [etadmax]
****************************************************************************************
*********** 6 body final states : p(a) + p(b) + p(c)+ p(d) + p(e) + p(f) ***************
****************************************************************************************
0d0 ! [ptamin]
0d0 ! [ptbmin]
0d0 ! [ptcmin]
0d0 ! [ptdmin]
0d0 ! [ptemin]
0d0 ! [ptfmin]
-2.5d1 ! [etaamin]
2.5d1 ! [etaamax]
-2.5d1 ! [etabmin]
2.5d1 ! [etabmax]
-2.5d1 ! [etacmin]
2.5d1 ! [etacmax]
-2.5d1 ! [etadmin]
2.5d1 ! [etadmax]
-2.5d1 ! [etaemin]
2.5d1 ! [etaemax]
-2.5d1 ! [etafmin]
2.5d1 ! [etafmax]
****************************************************************************************
******* Jet algorithm parameters
****************************************************************************************
0.6d0 ! [rjet] : Jet Radius
'antikt' ! [jalg] : Jet algorithm ('antikt','kt')
****************************************************************************************
****** chi_c/b two-body decays
****************************************************************************************
0.133d0 ! [m2b] : mass of decay particles
211 ! [pdgid1] : PDG number particle 1
-211 ! [pdgid2] : PDG number particle 2
****************************************************************************************
****** ALP parameters
****************************************************************************************
5d-2 ! [malp] : ALP mass (GeV)
1d-4 ! [gax] : ALP coupling (GeV^-2)
'ps' ! [alpt] : AlP type (scalar - 'sc', pseudoscalar - 'ps')
****************************************************************************************
****** Monopole parameters
****************************************************************************************
500d0 ! [mpol] : Monopole mass
933d0 ! [mmon] : Monopolium mass
10d0 ! [gamm] : Monopolium width
****************************************************************************************
****** SUSY parameters
****************************************************************************************
100d0 ! [mcharg] : Chargino/Slepton mass
80d0 ! [mneut] : Neutralino mass
****************************************************************************************
****************************************************************************************
//--------------------------------------------------------------------------
// File and Version Information:
// Atlas_HEPEVT_modified.cxx
//
// Description:
//
// Needed for interface of Superchic_i to Superchic3.03 generator
// to store the umodified HEPEVT common.
//
// Copied directly from the Tauolo_i and modified a bit
// Author List:
// Borut Paul Kersevan (BPK), June 2003
//
// Copyright Information:
//
// 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)
//------------------------------------------------------------------------
#include "Superchic_i/Atlas_HEPEVT_modified.h"
// set pointer to zero at start
Atlas_HEPEVT_modified::HEPEVT* Atlas_HEPEVT_modified::s_atlas_HEPEVT =0;
// Constructor
Atlas_HEPEVT_modified::Atlas_HEPEVT_modified()
{
m_atlas_HEPEVT.nevhep = 0;
m_atlas_HEPEVT.nhep = 0;
for(int i=0 ; i< s_lenNmxhep; i++ ){
m_atlas_HEPEVT.idhep[i] = 0;
m_atlas_HEPEVT.isthep[i] = 0;
}
m_dummy = 0;
m_realdummy = 0.;
}
// Destructor
Atlas_HEPEVT_modified::~Atlas_HEPEVT_modified()
{
}
void Atlas_HEPEVT_modified::fill()
{
init(); // check COMMON is initialized
m_atlas_HEPEVT= *(s_atlas_HEPEVT);
return;
}
void Atlas_HEPEVT_modified::spill()
{
s_atlas_HEPEVT =0; //re-init the pointer
init(); // check COMMON is initialized
*(s_atlas_HEPEVT)=m_atlas_HEPEVT;
return;
}
int& Atlas_HEPEVT_modified::nevhep() {
init(); // check COMMON is initialized
return s_atlas_HEPEVT->nevhep;
}
int& Atlas_HEPEVT_modified::nhep() {
init(); // check COMMON is initialized
return s_atlas_HEPEVT->nhep;
}
int& Atlas_HEPEVT_modified::isthep(int nh) {
init(); // check COMMON is initialized
if( nh < 1 || nh > lenNmxhep())
{
m_dummy = -999;
return m_dummy;
}
return s_atlas_HEPEVT->isthep[nh-1];
}
int& Atlas_HEPEVT_modified::idhep(int nh) {
init(); // check COMMON is initialized
if( nh < 1 || nh > lenNmxhep())
{
m_dummy = -999;
return m_dummy;