diff --git a/macro/run/run_transport_beam.C b/macro/run/run_transport_beam.C
new file mode 100644
index 0000000000000000000000000000000000000000..8a8a5c772dc218e30543a44982baf9b20ceb3e88
--- /dev/null
+++ b/macro/run/run_transport_beam.C
@@ -0,0 +1,134 @@
+/** @file run_transport_beam.C
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 9 September 2020
+ **/
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include "CbmBeamGenerator";
+#include "CbmTransport.h"
+#include "FairSystemInfo.h"
+#include "TStopwatch.h"
+#endif
+
+
+/** @brief run_transport_beam
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @since 9 September 2020
+ ** @param nEvents Number of events from input to transport
+ ** @param setupName Name of predefined geometry setup
+ ** @param output Name of output data set
+ ** @param inputFile Name of input file
+ **
+ ** Macro for transport simulation of beam particles in the CBM setup.
+ ** Properties of the beam have to be specified in the macro body.
+ **
+ ** For the setup (geometry and field), predefined setups can be chosen
+ ** by the second argument. Available setups are in geometry/setup.
+ **
+ ** The output file will be named [output].tra.root.
+ ** A parameter file [output].par.root will be created.
+ ** The geometry (TGeoManager) will be written into [output].geo.root.
+ **
+ ** For the options for the transport simulation, see the documentation of
+ ** CbmTransport.
+ **/
+void run_transport_beam(Int_t nEvents = 1000,
+ const char* setupName = "sis100_electron",
+ const char* output = "beam") {
+
+ // --- Logger settings ----------------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "run_transport_beam"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+
+ // ----- File names ---------------------------------------------------
+ TString dataset(output);
+ TString outFile = dataset + ".tra.root";
+ TString parFile = dataset + ".par.root";
+ TString geoFile = dataset + ".geo.root";
+ // ------------------------------------------------------------------------
+
+
+ // ---- Define the beam -----------------------------------------------
+ // Beam species here is full stripped Gold
+ UInt_t beamZ = 79; // atomic number
+ UInt_t beamA = 197; // mass number
+ UInt_t beamQ = 79; // charge (fully stripped)
+ Double_t beamP = 12.; // Momentum 12 GeV per nucleon
+ Double_t beamStartZ = -1.; // Starts at z = -1 cm
+ // ------------------------------------------------------------------------
+
+
+ // ---- Define the beam profile ---------------------------------------
+ Double_t beamFocusZ = 0.; // z coordinate of beam focal plane
+ Double_t beamMeanX = 0.; // mean x position of beam [cm]
+ Double_t beamSigmaX = 0.1; // Gaussian sigma of beam x position [cm]
+ Double_t beamMeanY = 0.; // mean y position of beam [cm]
+ Double_t beamSigmaY = 0.1; // Gaussian sigma of beam y position [cm]
+ Double_t beamMeanTx = 0.; // mean x-z angle of beam [rad]
+ Double_t beamSigmaTx = 1.; // Gaussian sigma of beam x-z angle [rad]
+ Double_t beamMeanTy = 0.; // mean x-y angle of beam [rad]
+ Double_t beamSigmaTy = 1.; // Gaussian sigma of beam x-y angle [rad]
+ // ------------------------------------------------------------------------
+
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+
+ // --- Transport run ----------------------------------------------------
+ CbmTransport run;
+ run.SetOutFileName(outFile);
+ run.SetParFileName(parFile);
+ run.SetGeoFileName(geoFile);
+ run.LoadSetup(setupName);
+ run.SetTarget("Gold", 0.025, 2.5);
+ run.SetBeamPosition(beamMeanX, beamMeanY, beamSigmaX, beamSigmaY, beamFocusZ);
+ run.SetBeamAngle(beamMeanTx, beamMeanTy, beamSigmaTx, beamSigmaTy);
+ run.AddInput(new CbmBeamGenerator(beamZ, beamA, beamQ, beamP, beamStartZ));
+ run.Run(nEvents);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << outFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << "s"
+ << std::endl
+ << std::endl;
+ // ------------------------------------------------------------------------
+
+
+ // ----- Resource monitoring ------------------------------------------
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+ // ------------------------------------------------------------------------
+
+} // End of macro
diff --git a/sim/transport/base/CbmEventGenerator.cxx b/sim/transport/base/CbmEventGenerator.cxx
index dce2d1f351200faf6a70b507c21b24cc738236ba..87a7b9a27e53dbe906af386058c02a4ff069e326 100644
--- a/sim/transport/base/CbmEventGenerator.cxx
+++ b/sim/transport/base/CbmEventGenerator.cxx
@@ -6,12 +6,12 @@
#include "CbmEventGenerator.h"
#include "CbmBeam.h"
-#include "FairGenericStack.h"
-#include "FairLogger.h"
-#include "FairMCEventHeader.h"
-#include "TMath.h"
-#include "TRandom.h"
-#include "TVector3.h"
+#include <FairGenericStack.h>
+#include <FairLogger.h>
+#include <FairMCEventHeader.h>
+#include <TMath.h>
+#include <TRandom.h>
+#include <TVector3.h>
#include <cassert>
@@ -20,7 +20,9 @@ CbmEventGenerator::CbmEventGenerator()
: FairPrimaryGenerator()
, fBeamProfile()
, fTarget()
- , fForceVertexInTarget(kTRUE) {
+ , fForceVertexInTarget(kTRUE)
+ , fForceVertexAtZ(kFALSE)
+ , fVertexZ(0.) {
// Mother class members
fName = "EventGenerator";
fBeamAngle = kTRUE;
@@ -33,6 +35,15 @@ CbmEventGenerator::~CbmEventGenerator() {}
// -------------------------------------------------------------------------
+// ----- Force vertex at a given z position ----------------------------
+void CbmEventGenerator::ForceVertexAtZ(Double_t zVertex) {
+ fForceVertexAtZ = kTRUE;
+ fForceVertexInTarget = kFALSE;
+ fVertexZ = zVertex;
+}
+// -------------------------------------------------------------------------
+
+
// ----- Generate the input event --------------------------------------
Bool_t CbmEventGenerator::GenerateEvent(FairGenericStack* stack) {
@@ -101,7 +112,9 @@ Bool_t CbmEventGenerator::GenerateEvent(FairGenericStack* stack) {
// ----- Generate event vertex -----------------------------------------
void CbmEventGenerator::MakeVertex() {
- if (fForceVertexInTarget && fTarget)
+ if (fForceVertexAtZ)
+ MakeVertexAtZ();
+ else if (fForceVertexInTarget && fTarget)
MakeVertexInTarget();
else
MakeVertexInFocalPlane();
@@ -109,6 +122,21 @@ void CbmEventGenerator::MakeVertex() {
// -------------------------------------------------------------------------
+// ----- Generate event vertex in the beam focal plane -----------------
+void CbmEventGenerator::MakeVertexAtZ() {
+
+ std::unique_ptr<CbmBeam> beam;
+ beam = fBeamProfile.GenerateBeam();
+ TVector3 point(0., 0., fVertexZ); /// A point in the plane z = zVertex
+ TVector3 norm(0, 0., 1.); /// Normal on the plane z = const.
+ fVertex = beam->ExtrapolateToPlane(point, norm);
+ fBeamAngleX = beam->GetThetaX();
+ fBeamAngleY = beam->GetThetaY();
+ fBeamDirection = beam->GetDirection();
+}
+// -------------------------------------------------------------------------
+
+
// ----- Generate event vertex in the beam focal plane -----------------
void CbmEventGenerator::MakeVertexInFocalPlane() {
@@ -191,7 +219,10 @@ void CbmEventGenerator::Print(Option_t*) const {
<< fPhiMax << " rad";
else
LOG(info) << "Fixed event plane angle = 0";
- LOG(info) << "Number of generators = " << fGenList->GetEntries();
+ LOG(info) << "Number of generators " << fGenList->GetEntries();
+ for (Int_t iGen = 0; iGen < fGenList->GetEntries(); iGen++) {
+ fGenList->At(iGen)->Print();
+ }
}
// -------------------------------------------------------------------------
diff --git a/sim/transport/base/CbmEventGenerator.h b/sim/transport/base/CbmEventGenerator.h
index e73f907a09388da952d27d858497b69098760fb3..44d92b570a2068389c188d57e832eec4b15f7083 100644
--- a/sim/transport/base/CbmEventGenerator.h
+++ b/sim/transport/base/CbmEventGenerator.h
@@ -10,7 +10,7 @@
#include "CbmBeamProfile.h"
#include "CbmTarget.h"
-#include "FairPrimaryGenerator.h"
+#include <FairPrimaryGenerator.h>
#include <memory>
class FairGenericStack;
@@ -42,6 +42,15 @@ public:
virtual ~CbmEventGenerator();
+ /** @brief Force event vertex to be at a given z
+ ** @param zVertex z coordinate of event vertex
+ **
+ ** The event vertex will be sampled in x and y from the specified
+ ** beam properties (profile in focal plane and angular distribution),
+ **/
+ void ForceVertexAtZ(Double_t zVertex);
+
+
/** @brief Enable or disable forcing the vertex to be in the target
** @param choice If true, the vertex will be generated in the target
**/
@@ -130,6 +139,8 @@ private:
CbmBeamProfile fBeamProfile; ///< Beam properties
std::shared_ptr<const CbmTarget> fTarget; //! Target properties
Bool_t fForceVertexInTarget; ///< If set, vertex must be in target
+ Bool_t fForceVertexAtZ; ///< If set, vertex must be at given z
+ Double_t fVertexZ; ///< forced z coordinate of event vertex
/** @brief Generate beam angle
@@ -146,19 +157,29 @@ private:
**
** Will be called from FairPrimaryGenerator::GenerateEvent().
** Beam position and angles in the focal plane are sampled
- ** from the specified distributions. If a target is specified,
- ** the event vertex is sampled from the beam trajectory within the
- ** target - flat distributions between entry and exit point, or
- ** always in the target centre plane, depending on the choice set
- ** with SetSmearVertexZ(). If no target is specified, the event vertex
- ** is the beam position in the focal plane.
+ ** from the specified distributions. There are three options to generate
+ ** the event vertex:
+ ** 1. If a vertex z position is specified by a call to ForceVertexAtZ,
+ ** the event vertex is the beam extrapolation to the specified z.
+ ** 2. Else, if a target is specified, the event vertex is sampled from
+ ** the beam trajectory within the target - flat distributions between
+ ** entry and exit point, or always in the target centre plane, depending
+ ** on the choice set with SetSmearVertexZ().
+ ** 3. Else, the event vertex is the beam position in the focal plane.
**/
virtual void MakeVertex();
+ /** @brief Generate event vertex position at a given z
+ **
+ ** Will be used if ForceVertexAtZ was called.
+ **/
+ void MakeVertexAtZ();
+
+
/** @brief Generate event vertex position in the beam focal plane
**
- ** Will be called if no target was specified.
+ ** Will be used if no target was specified.
**/
virtual void MakeVertexInFocalPlane();
diff --git a/sim/transport/generators/CMakeLists.txt b/sim/transport/generators/CMakeLists.txt
index 684e3a8a4faa98240c5e5b84453abfa3cc77d679..f03069cc8189244821877afb841a926c60cd038c 100644
--- a/sim/transport/generators/CMakeLists.txt
+++ b/sim/transport/generators/CMakeLists.txt
@@ -8,6 +8,7 @@ Set(LIBRARY_NAME CbmSimGenerators)
# ----- Compilation sources ----------------------------
set(SRCS
+CbmBeamGenerator.cxx
CbmPhsdGenerator.cxx
CbmPlutoGenerator.cxx
CbmShieldGenerator.cxx
@@ -58,6 +59,7 @@ ${Boost_LIBRARY_DIRS}
Set(DEPENDENCIES
CbmBase
CbmData
+CbmSimBase
)
# ---------------------------------------------------------
diff --git a/sim/transport/generators/CbmBeamGenerator.cxx b/sim/transport/generators/CbmBeamGenerator.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4f5fb4095c6cc37b868b552c648efb9dbc4a67f0
--- /dev/null
+++ b/sim/transport/generators/CbmBeamGenerator.cxx
@@ -0,0 +1,107 @@
+/** @file CbmBeamGenerator.h
+ ** @author V. Friese <v.friese@gsi.de>
+ ** @since 8 September 2020
+ **/
+
+
+#include "CbmBeamGenerator.h"
+
+#include "CbmEventGenerator.h"
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include "FairLogger.h"
+#include "FairMCEventHeader.h"
+#include "FairPrimaryGenerator.h"
+#include "FairRunSim.h"
+#include "TFile.h"
+#include "TRandom.h"
+#include "TTree.h"
+#include "TVector3.h"
+#include "UEvent.h"
+#include "UParticle.h"
+#include "URun.h"
+#include <cassert>
+#include <sstream>
+
+
+// ----- Default constructor --------------------------------------------
+CbmBeamGenerator::CbmBeamGenerator()
+ : FairGenerator("BeamGenerator", "CBM generator")
+ , fP(0.)
+ , fStartZ(0.)
+ , fIon(nullptr) {}
+// --------------------------------------------------------------------------
+
+
+// ----- Standard constructor --------------------------------------------
+CbmBeamGenerator::CbmBeamGenerator(UInt_t beamZ,
+ UInt_t beamA,
+ UInt_t beamQ,
+ Double_t momentum,
+ Double_t startZ)
+ : FairGenerator("BeamGenerator", "CBM generator")
+ , fP(momentum * Double_t(beamA))
+ , fStartZ(startZ)
+ , fIon(nullptr) {
+
+ // --- Create the ion species and add it to the particle list
+ char name[20];
+ sprintf(name, "Beam_%d_%d_%d", beamZ, beamA, beamQ);
+ fIon = new FairIon(name, beamZ, beamA, beamQ);
+ FairRunSim* run = FairRunSim::Instance();
+ assert(run);
+ run->AddNewIon(fIon);
+
+ // --- Tell the event generator to force the vertex at the given z
+ FairPrimaryGenerator* primGen = run->GetPrimaryGenerator();
+ CbmEventGenerator* evGen = dynamic_cast<CbmEventGenerator*>(primGen);
+ assert(evGen);
+ evGen->ForceVertexAtZ(startZ);
+ evGen->SmearVertexZ(kFALSE);
+}
+// --------------------------------------------------------------------------
+
+
+// ----- Destructor -----------------------------------------------------
+CbmBeamGenerator::~CbmBeamGenerator() {}
+// --------------------------------------------------------------------------
+
+
+// ----- Print info -----------------------------------------------------
+void CbmBeamGenerator::Print(Option_t*) const { LOG(info) << ToString(); }
+// --------------------------------------------------------------------------
+
+
+// ----- Generate input event -------------------------------------------
+Bool_t CbmBeamGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
+
+ TParticlePDG* ion = TDatabasePDG::Instance()->GetParticle(fIon->GetName());
+ assert(ion);
+
+ // Note that the beam position in x and y and the beam angle are generated
+ // by CbmEventGenerator. Here, we generate the ion thus in z direction
+ // and at zero coordinates. It will be properly placed and rotated
+ // according to the beam profile specified to CbmEventGenerator.
+ primGen->AddTrack(ion->PdgCode(), 0., 0., fP, 0., 0., 0.);
+
+ return kTRUE;
+}
+// --------------------------------------------------------------------------
+
+
+// ----- Info -----------------------------------------------------------
+std::string CbmBeamGenerator::ToString() const {
+
+ std::stringstream ss;
+ ss << GetName() << " ion: " << fIon->GetName() << " Z " << fIon->GetZ()
+ << " A " << fIon->GetA() << " Q " << fIon->GetQ() << " mass "
+ << fIon->GetMass() << ", momentum " << fP
+ << " GeV/u, start Z = " << fStartZ;
+
+ return ss.str();
+}
+// --------------------------------------------------------------------------
+
+
+ClassImp(CbmBeamGenerator);
diff --git a/sim/transport/generators/CbmBeamGenerator.h b/sim/transport/generators/CbmBeamGenerator.h
new file mode 100644
index 0000000000000000000000000000000000000000..d7f7e9ab5cd6f6161e25a7e6744150031ffb120d
--- /dev/null
+++ b/sim/transport/generators/CbmBeamGenerator.h
@@ -0,0 +1,86 @@
+/** @file CbmBeamGenerator.h
+ ** @author V.Friese <v.friese@gsi.de>
+ ** @since 8 September 2020
+ **/
+
+
+#ifndef CBMBEAMGENERATOR_H
+#define CBMBEAMGENERATOR_H 1
+
+
+#include <FairGenerator.h>
+#include <FairIon.h>
+#include <Rtypes.h>
+#include <string>
+#include <vector>
+
+class FairPrimaryGenerator;
+
+
+/** @class CbmBeamGenerator.h
+ ** @brief Generates beam ions for transport simulation
+ ** @author V.Friese <v.friese@gsi.de>
+ ** @since 8 September 2020
+ **
+ ** The BeamGenerator is intended for the simulation of (non-interacting)
+ ** beam particles through the setup, to study the beam-related background
+ ** e.g. from delta electrons created by the beam passing the target or
+ ** other materials in the setup. One beam particle (ion) per event is created.
+ ** The user has to specify the ion species, the momentum or kinetic energy,
+ ** and the starting point of the beam. The beam properties (profile and angular
+ ** distribution) have to be specified to the CbmEventGenerator.
+ **
+ ** It is not recommended to use the BeamGenerator together with other
+ ** generators in the same transport run, because it will force the event
+ ** vertex to be at the specified z position and will deactivate vertex smearing
+ ** in z as well as event plane rotation.
+ **/
+class CbmBeamGenerator : public FairGenerator {
+
+public:
+ /** @brief Default constructor (should not be used) **/
+ CbmBeamGenerator();
+
+
+ /** @brief Standard constructor
+ ** @param beamZ Atomic number (number of protons)
+ ** @param beamA Atomic mass number (number of nucleons)
+ ** @param beamQ Electric charge
+ ** @param momentum Momentum per nucleon [GeV]
+ ** @param zStart z coordinate of beam start position
+ **/
+ CbmBeamGenerator(UInt_t beamZ,
+ UInt_t beamA,
+ UInt_t beamQ,
+ Double_t momentum,
+ Double_t zStart);
+
+
+ /** @brief Destructor **/
+ virtual ~CbmBeamGenerator();
+
+
+ /** @brief Print info to logger **/
+ virtual void Print(Option_t* opt = "") const;
+
+
+ /** @brief Generate one event (abstract in base class)
+ ** @param primGen Pointer to the FairPrimaryGenerator
+ **/
+ virtual Bool_t ReadEvent(FairPrimaryGenerator* primGen);
+
+
+ /** @brief Info to string **/
+ std::string ToString() const;
+
+
+private:
+ Double_t fP; ///< Total momentum [GeV]
+ Double_t fStartZ; ///< z coordinate of start point
+ FairIon* fIon; ///< Ion type
+
+
+ ClassDef(CbmBeamGenerator, 1);
+};
+
+#endif
diff --git a/sim/transport/generators/CbmSimGeneratorsLinkDef.h b/sim/transport/generators/CbmSimGeneratorsLinkDef.h
index 7ebdb36cf339618d63c615898bebc6b322db9b33..953051b44eec1f785c2dd85f5348f848bcb9422f 100644
--- a/sim/transport/generators/CbmSimGeneratorsLinkDef.h
+++ b/sim/transport/generators/CbmSimGeneratorsLinkDef.h
@@ -5,6 +5,7 @@
#pragma link off all functions;
// --- transport/base
+#pragma link C++ class CbmBeamGenerator + ;
#pragma link C++ class CbmPhsdGenerator + ;
#pragma link C++ class CbmPlutoGenerator + ;
#pragma link C++ class CbmShieldGenerator + ;
diff --git a/sim/transport/steer/CbmTransport.cxx b/sim/transport/steer/CbmTransport.cxx
index 5e3994a81a59b8b3221831c452776dd043129e8f..e56866d543a125d7a6eda7a6a047a483f8015624 100644
--- a/sim/transport/steer/CbmTransport.cxx
+++ b/sim/transport/steer/CbmTransport.cxx
@@ -89,6 +89,7 @@ CbmTransport::CbmTransport()
// By default, vertex smearing along the beam is activated
fEventGen->SmearVertexZ(kTRUE);
+ fRun->SetGenerator(fEventGen); // has to be available for some generators
}
// --------------------------------------------------------------------------
@@ -171,6 +172,14 @@ void CbmTransport::ConfigureVMC() {
// --------------------------------------------------------------------------
+// ----- Force creation of event vertex at a given z position -----------
+void CbmTransport::ForceVertexAtZ(Double_t zVertex) {
+ assert(fEventGen);
+ fEventGen->ForceVertexAtZ(zVertex);
+}
+// --------------------------------------------------------------------------
+
+
// ----- Force creation of event vertex in the target -------------------
void CbmTransport::ForceVertexInTarget(Bool_t choice) {
assert(fEventGen);
@@ -276,10 +285,6 @@ void CbmTransport::InitEventGenerator() {
<< ": Beam profile is not consistent with target!";
} //? Target not consistent with beam
} //? Target specified
-
-
- // --- Register event generator to run
- fRun->SetGenerator(fEventGen);
}
// --------------------------------------------------------------------------
diff --git a/sim/transport/steer/CbmTransport.h b/sim/transport/steer/CbmTransport.h
index aaf670902ce4681d2d9543e98fb353c199ccc5f6..5e927912d953a5f7616f8230f04fa66ae8ff6618 100644
--- a/sim/transport/steer/CbmTransport.h
+++ b/sim/transport/steer/CbmTransport.h
@@ -65,6 +65,16 @@ public:
void ConfigureVMC();
+ /** @brief Force the event vertex to be at a given z position
+ ** @param zVertex z position of event vertex
+ **
+ ** The event vertex will be determined by the beam intersection with
+ ** the specified z plane. The beam properties (position in focal plane
+ ** and direction) are sampled from the specified beam profile.
+ **/
+ void ForceVertexAtZ(Double_t zVertex);
+
+
/** @brief Enable or disable forcing the vertex to be in the target
** @param choice If true, the vertex will be generated in the target
**/