From 27c290e54e53ecc9ce3900333991bc64a64d0c68 Mon Sep 17 00:00:00 2001
From: lintao <lintao51@gmail.com>
Date: Sun, 10 Nov 2019 22:50:36 +0800
Subject: [PATCH] WIP: start to implement the particle gun.
---
Generator/CMakeLists.txt | 35 +++++++----
Generator/src/GtGunTool.cpp | 119 ++++++++++++++++++++++++++++++++++++
Generator/src/GtGunTool.h | 15 +++++
3 files changed, 158 insertions(+), 11 deletions(-)
diff --git a/Generator/CMakeLists.txt b/Generator/CMakeLists.txt
index 5d0c6d94..16be6fa4 100644
--- a/Generator/CMakeLists.txt
+++ b/Generator/CMakeLists.txt
@@ -15,7 +15,10 @@ set(GenAlgo_srcs
)
set(GenAlgo_incs src)
-find_package(ROOT COMPONENTS RIO Tree TreePlayer MathCore Net Graf3d Graf Gpad REQUIRED)
+find_package(Geant4 REQUIRED)
+include(${Geant4_USE_FILE})
+
+find_package(ROOT COMPONENTS RIO Tree TreePlayer MathCore Net Graf3d Graf Gpad EG REQUIRED)
find_package(LCIO)
find_package(podio)
find_package(plcio)
@@ -39,16 +42,26 @@ endif(HepMC_FOUND)
INCLUDE_DIRECTORIES(${GenAlgo_incs})
gaudi_add_module(GenAlgo ${GenAlgo_srcs}
-INCLUDE_DIRS
-GaudiKernel
-FWCore
-${LCIO_INCLUDE_DIRS}
-${podio_INCLUDE_DIRS}
-${plcio_INCLUDE_DIRS}
-${ROOT_INCLUDE_DIRS}
-HepMC
-LINK_LIBRARIES GaudiKernel ${LCIO_LIBRARIES} ${podio_LIBRARIES} ROOT ${plcio_LIBRARY_DIR}/libplcio.so ${plcio_LIBRARY_DIR}/libplcioDict.so FWCore HepMC
-)
+ INCLUDE_DIRS
+ GaudiKernel
+ FWCore
+ Geant4
+ ${LCIO_INCLUDE_DIRS}
+ ${podio_INCLUDE_DIRS}
+ ${plcio_INCLUDE_DIRS}
+ ${ROOT_INCLUDE_DIRS}
+ HepMC
+ LINK_LIBRARIES
+ GaudiKernel
+ ${LCIO_LIBRARIES}
+ ${podio_LIBRARIES}
+ ROOT
+ ${plcio_LIBRARY_DIR}/libplcio.so
+ ${plcio_LIBRARY_DIR}/libplcioDict.so
+ FWCore
+ HepMC
+ Geant4
+ )
#gaudi_add_test(Reader FRAMEWORK options/read.py)
###########################
diff --git a/Generator/src/GtGunTool.cpp b/Generator/src/GtGunTool.cpp
index b4087cce..f9d7da9c 100644
--- a/Generator/src/GtGunTool.cpp
+++ b/Generator/src/GtGunTool.cpp
@@ -1,10 +1,47 @@
#include "GtGunTool.h"
+#include "TDatabasePDG.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Random/RandGauss.h"
+
DECLARE_COMPONENT(GtGunTool)
StatusCode
GtGunTool::initialize() {
StatusCode sc;
+ // the particles names/pdgs and energies should be specified.
+ if (m_particles.value().size()==0) {
+ error() << "Please specify the list of particle names/pdgs" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ if (m_energies.value().size() != m_particles.value().size()) {
+ error() << "Mismatched energies and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ // others should be empty or specify
+ if (m_thetamins.value().size()
+ && m_thetamins.value().size() != m_particles.value().size()) {
+ error() << "Mismatched thetamins and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ if (m_thetamaxs.value().size()
+ && m_thetamaxs.value().size() != m_particles.value().size()) {
+ error() << "Mismatched thetamaxs and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ if (m_phimins.value().size()
+ && m_phimins.value().size() != m_particles.value().size()) {
+ error() << "Mismatched phimins and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ if (m_phimaxs.value().size()
+ && m_phimaxs.value().size() != m_particles.value().size()) {
+ error() << "Mismatched phimaxs and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
return sc;
}
@@ -18,6 +55,88 @@ GtGunTool::finalize() {
bool
GtGunTool::mutate(MyHepMC::GenEvent& event) {
+ TDatabasePDG* db_pdg = TDatabasePDG::Instance();
+
+ // The default unit here is GeV.
+ // but we don't add the unit, because in geant4 it is multiplied.
+
+ for (int i = 0; i < m_particles.value().size(); ++i) {
+ const std::string& particle_name = m_particles.value()[i];
+ int pdgcode = 0;
+ double mass = 0;
+
+ TParticlePDG* particle = db_pdg->GetParticle(particle_name.c_str());
+ if (particle) {
+ pdgcode = particle->PdgCode();
+ mass = particle->Mass(); // GeV
+ } else {
+ // guess it is pdg code
+ pdgcode = atol(particle_name.c_str());
+ if (!pdgcode) {
+ error() << "Unsupported particle name/pdgcode " << particle_name << endmsg;
+ return false;
+ }
+ }
+
+ double energy = m_energies.value()[i];
+
+ // create the MC particle
+ plcio::MCParticle mcp = event.m_mc_vec.create();
+ mcp.setPDG(pdgcode);
+ mcp.setGeneratorStatus(1);
+ mcp.setSimulatorStatus(1);
+ // mcp.setCharge();
+ mcp.setTime(0.0);
+ // mcp.setVertex();
+ // mcp.setEndpoint();
+
+ // assume energy is momentum
+ double p = energy;
+
+ // direction
+ // by default, randomize the direction
+ double costheta = CLHEP::RandFlat::shoot(-1, 1);
+ double phi = 360*CLHEP::RandFlat::shoot()*CLHEP::degree;
+ double sintheta = sqrt(1.-costheta*costheta);
+
+ // check if theta min/max is set
+
+ if (i < m_thetamins.value().size()
+ && i < m_thetamaxs.value().size()) {
+ double thetamin = m_thetamins.value()[i];
+ double thetamax = m_thetamaxs.value()[i];
+
+ if (thetamin == thetamax) { // fixed theta
+ costheta = cos(thetamin);
+ sintheta = sin(thetamin);
+ info() << "theta is fixed: " << thetamin << endmsg;
+ }
+ }
+
+ if (i < m_phimins.value().size()
+ && i < m_phimaxs.value().size()) {
+ double phimin = m_phimins.value()[i];
+ double phimax = m_phimaxs.value()[i];
+
+ if (phimin == phimax) { // fixed phi
+ phi = phimin;
+ info() << "phi is fixed: " << phimin << endmsg;
+ }
+ }
+
+ double px = p*sintheta*cos(phi);
+ double py = p*sintheta*sin(phi);
+ double pz = p*costheta;
+
+ mcp.setMomentum(plcio::FloatThree(px,py,pz));
+ // mcp.setMomentumAtEndpoint();
+ // mcp.setSpin();
+ // mcp.setColorFlow();
+
+ }
+
+ // event.SetEventHeader( m_processed_event, -99, 9999, "Generator");
+
return true;
}
diff --git a/Generator/src/GtGunTool.h b/Generator/src/GtGunTool.h
index 4138a8ee..22b02a81 100644
--- a/Generator/src/GtGunTool.h
+++ b/Generator/src/GtGunTool.h
@@ -12,8 +12,11 @@
*/
#include <GaudiKernel/AlgTool.h>
+#include <GaudiKernel/Property.h>
#include "IGenTool.h"
+#include <vector>
+
class GtGunTool: public extends<AlgTool, IGenTool> {
public:
using extends::extends;
@@ -27,6 +30,18 @@ public:
bool finish() override;
bool configure_gentool() override;
+private:
+
+ Gaudi::Property<std::vector<std::string>> m_particles{this, "Particles"};
+
+ Gaudi::Property<std::vector<double>> m_energies{this, "Energies"};
+
+ Gaudi::Property<std::vector<double>> m_thetamins{this, "ThetaMins"};
+ Gaudi::Property<std::vector<double>> m_thetamaxs{this, "ThetaMaxs"};
+
+ Gaudi::Property<std::vector<double>> m_phimins{this, "PhiMins"};
+ Gaudi::Property<std::vector<double>> m_phimaxs{this, "PhiMaxs"};
+
};
--
GitLab