From 7ec380a9caf2e7ce451d232ac4316fc5009dffa3 Mon Sep 17 00:00:00 2001
From: "guofangyi@ihep.ac.cn" <guofangyi@ihep.ac.cn>
Date: Wed, 14 Aug 2024 03:35:56 +0000
Subject: [PATCH] beam induced bkg simulation with time info
---
Examples/options/beambkgsim.py | 127 ------------
Examples/options/sim_beambkg.py | 197 +++++++++++++++++++
Generator/src/BeamBackgroundFileParserV1.cpp | 31 ++-
Generator/src/BeamBackgroundFileParserV1.h | 3 +-
Generator/src/BeamBackgroundFileParserV2.cpp | 52 +++--
Generator/src/BeamBackgroundFileParserV2.h | 8 +-
Generator/src/GtBeamBackgroundTool.cpp | 20 +-
Generator/src/GtBeamBackgroundTool.h | 6 +
Generator/src/GtGunTool.cpp | 13 +-
Generator/src/GtGunTool.h | 3 +
Generator/src/StdHepRdr.cpp | 2 +-
Generator/src/StdHepRdr.h | 1 +
12 files changed, 289 insertions(+), 174 deletions(-)
delete mode 100644 Examples/options/beambkgsim.py
create mode 100755 Examples/options/sim_beambkg.py
diff --git a/Examples/options/beambkgsim.py b/Examples/options/beambkgsim.py
deleted file mode 100644
index a786b3c7..00000000
--- a/Examples/options/beambkgsim.py
+++ /dev/null
@@ -1,127 +0,0 @@
-#!/usr/bin/env python
-#Author: Zhan Li <lizhan@ihep.ac.cn>
-#Created [2024-03-07 Thu 14:53]
-
-import os
-import sys
-
-import Gaudi.Configuration
-from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_, k4DataSvc, GeomSvc
-from Configurables import TimeProjectionChamberSensDetTool
-from Configurables import GenAlgo
-from Configurables import GtGunTool
-from Configurables import StdHepRdr
-from Configurables import SLCIORdr
-from Configurables import HepMCRdr
-from Configurables import GenPrinter
-from Configurables import GtBeamBackgroundTool
-from Configurables import DetSimSvc
-from Configurables import DetSimAlg
-from Configurables import AnExampleDetElemTool
-from Configurables import PodioOutput
-from Configurables import ApplicationMgr
-
-seed = [42]
-
-rndmengine = Gaudi.Configuration.HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine")
-rndmengine.SetSingleton = True
-rndmengine.Seeds = seed
-
-rndmgensvc = RndmGenSvc("RndmGenSvc")
-rndmgensvc.Engine = rndmengine.name()
-
-
-dsvc = k4DataSvc("EventDataSvc")
-#geometry_option = "CepC_v4-onlyVXD.xml"
-#geometry_option = "CepC_v4_onlyTracker.xml"
-geometry_option = "CepC_v4.xml"
-
-geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option)
-geosvc = GeomSvc("GeomSvc")
-geosvc.compact = geometry_path
-
-#Previously I do not have these 2 lines
-tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
-tpc_sensdettool.TypeOption = 1
-
-
-# Physics Generator
-# ------------------ WARNING ------------------
-# Files for code test only.
-# Contact @Haoyu Shi (shihy@ihep.ac.cn) for reliable simulation samples and rates.
-# ---------------------------------------------
-bg = GtBeamBackgroundTool("GtBeamBackgroundTool")
-bg.InputFileMap = {
- "BGB":"/cefs/higgs/guofy/CEPCSW_BeamBkgSim/run/BGB_Higgs.root",
- "BGC":"/cefs/higgs/guofy/CEPCSW_BeamBkgSim/run/BGC_Higgs.root",
- "BTH":"/cefs/higgs/guofy/CEPCSW_BeamBkgSim/run/BTH_Higgs.root",
-# "TSC":"",
-# "Pair":"",
- "SR":"/cefs/higgs/guofy/CEPCSW_BeamBkgSim/run/SRBkg/"
-}
-bg.InputFormatMap = {
- "BGB":"BeamBackgroundFileParserV1",
- "BGC":"BeamBackgroundFileParserV1",
- "BTH":"BeamBackgroundFileParserV1",
-# "TSC":"BeamBackgroundFileParserV1",
-# "Pair":"BeamBackgroundFileParserV1",
- "SR":"BeamBackgroundFileParserV2"
-}
-bg.InputRateMap = { # in Hz
- "BGB":1e7,
- "BGC":5e6,
- "BTH":5e5,
-# "TSC":"5e5",
-# "Pair":"",
- "SR":1
-}
-bg.TimeWindow = 1e-6
-bg.InputBeamEnergy = 120.
-bg.RotationAlongYMap = {
- "BGB":16.5e-3,
- "BGC":16.5e-3,
- "BTH":16.5e-3,
-# "TSC":"16.5e-3",
-# "Pair":"",
- "SR":0
-}
-
-gun = GtGunTool("GtGunTool")
-gun.Particles = ["gamma"]
-gun.EnergyMins = [10]
-gun.EnergyMaxs = [10]
-gun.ThetaMins = [90]
-gun.ThetaMaxs = [90]
-gun.PhiMins = [0]
-gun.PhiMaxs = [0]
-
-genprinter = GenPrinter("GenPrinter")
-
-genalg = GenAlgo("GenAlgo")
-genalg.GenTools = ["GtBeamBackgroundTool", "GtGunTool"]
-
-detsimsvc = DetSimSvc("DetSimSvc")
-
-detsimalg = DetSimAlg("DetSimAlg")
-detsimalg.RandomSeeds = seed
-
-
-detsimalg.RunCmds = []
-detsimalg.AnaElems = [
- "Edm4hepWriterAnaElemTool"
-]
-detsimalg.RootDetElem = "WorldDetElemTool"
-
-example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
-
-
-# POD I/O
-out = PodioOutput("outputalg")
-out.filename = "beam-SETv0.root"
-out.outputCommands = ["keep *"]
-
-ApplicationMgr( TopAlg = [genalg, detsimalg, out],
- EvtSel = 'NONE',
- EvtMax = 3,
- ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
-)
diff --git a/Examples/options/sim_beambkg.py b/Examples/options/sim_beambkg.py
new file mode 100755
index 00000000..e53127f8
--- /dev/null
+++ b/Examples/options/sim_beambkg.py
@@ -0,0 +1,197 @@
+#!/usr/bin/env python
+import os
+from Gaudi.Configuration import *
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+
+seed = [556]
+Nevt = 1
+yyy_input = "/cefs/higgs/songwz/summer24/bkgGenerator/pairHiggs/pair4Higgs/"
+yyy_output = '../simdir/sim_bkg_556.root'
+
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = seed
+
+rndmgensvc = RndmGenSvc("RndmGenSvc")
+rndmgensvc.Engine = rndmengine.name()
+
+geometry_option = "TDR_o1_v01/TDR_o1_v01.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+#ntsvc.Output = ["MyTuples DATAFILE='result.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+from Configurables import GtBeamBackgroundTool
+
+# gun = GtGunTool("GtGunTool")
+# gun.Particles = ["mu-"]
+# gun.PositionXs = [0.]
+# gun.PositionYs = [0.]
+# gun.PositionZs = [0.]
+# gun.EnergyMins = [5.0] # GeV
+# gun.EnergyMaxs = [5.0] # GeV
+# gun.ThetaMins = [90] # deg
+# gun.ThetaMaxs = [90] # deg
+# gun.PhiMins = [0] # deg
+# gun.PhiMaxs = [0] # deg
+# gun.Times = [0.75e3] # ns
+
+bg = GtBeamBackgroundTool("GtBeamBackgroundTool")
+Nbunch = 10
+TbunchSpacing = 355. # 50MW, higgs mode, and in ns
+bg.TimeWindow = Nbunch*TbunchSpacing
+bg.InputFileMap = {
+ "BGB":"/cefs/higgs/songwz/summer24/bkgGenerator/singleBeamHiggs/BGB.root",
+ "BGC":"/cefs/higgs/songwz/summer24/bkgGenerator/singleBeamHiggs/BGC.root",
+ "BTH":"/cefs/higgs/songwz/summer24/bkgGenerator/singleBeamHiggs/BTH.root",
+}
+bg.InputFileMap.update(dict([(f'pair{i}', yyy_input) for i in range(Nbunch)]))
+bg.InputFormatMap = {
+ "BGB":"BeamBackgroundFileParserV1",
+ "BGC":"BeamBackgroundFileParserV1",
+ "BTH":"BeamBackgroundFileParserV1",
+}
+bg.InputFormatMap.update(dict([(f'pair{i}', "BeamBackgroundFileParserV2") for i in range(Nbunch)]))
+bg.InputRateMap = { # in Hz
+ "BGB":83280.65,
+ "BGC":884002.12,
+ "BTH":623520.09,
+}
+bg.InputBeamEnergy = 120.
+bg.RotationAlongYMap = {
+ "BGB":16.5e-3,
+ "BGC":16.5e-3,
+ "BTH":16.5e-3
+}
+bg.RotationAlongYMap.update(dict([(f'pair{i}', 0) for i in range(Nbunch)]))
+bg.TimeBkgMap = dict([(f'pair{i}', TbunchSpacing*i) for i in range(Nbunch)]) #ns
+bg.NumberMcParticle = { #-1: pair use one file and single beam use rate*time; >=0: fixed number
+ "BGB":-1,
+ "BGC":-1,
+ "BTH":-1
+}
+bg.NumberMcParticle.update(dict([(f'pair{i}', -1) for i in range(Nbunch)]))
+
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC240/higgs/Higgs_10M/data/E240.Pnnh_gg.e0.p0.whizard195/nnh_gg.e0.p0.00100.stdhep"
+# stdheprdr.StartTime = 2.4e3 #ns
+
+# # lciordr = SLCIORdr("SLCIORdr")
+# # lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+# # hepmcrdr = HepMCRdr("HepMCRdr")
+# # hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtBeamBackgroundTool"]
+# genalg.GenTools = ["GtGunTool"]
+# genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+detsimsvc = DetSimSvc("DetSimSvc")
+
+from Configurables import DetSimAlg
+detsimalg = DetSimAlg("DetSimAlg")
+detsimalg.RandomSeeds = seed
+# detsimalg.VisMacs = ["vis.mac"]
+detsimalg.RunCmds = [
+# "/tracking/verbose 2",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaDoseElemTool",
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+detsimalg.PhysicsList = "QGSP_BERT_EMV"
+
+# from Configurables import ExampleAnaDoseElemTool
+# dosesimtool = ExampleAnaDoseElemTool("ExampleAnaDoseElemTool")
+# dosesimtool.Gridnbins = [220, 220, 60]
+# dosesimtool.Coormin = [0, 0, 0]
+# dosesimtool.Coormax = [2200, 2200, 600]
+# # dosesimtool.Gridnbins = [2200,1,4400]
+# # dosesimtool.Coormin = [0,-0.5,-2200]
+# # dosesimtool.Coormax = [2200,0.5,2200]
+# dosesimtool.Regionhist = [20,1.e-12,10]
+# dosesimtool.filename = "output/testdose_yyy_n_"
+# dosesimtool.Dosecofffilename = "dosecoffe/"
+# dosesimtool.HEHadroncut = 0.02
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+#gearsvc.GearXMLFile = "../../Detector/DetCEPCv4/compact/FullDetGear.xml"
+
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import TimeProjectionChamberSensDetTool
+tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
+tpc_sensdettool.TypeOption = 1
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+cal_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+cal_sensdettool.CalNamesMergeDisable = ["CaloDetector"]
+# cal_sensdettool.CalNamesApplyBirks = ["HcalBarrel"]
+
+
+# output
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = yyy_output
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr(
+ TopAlg = [genalg, detsimalg, out], #digiVXD, digiSIT, digiSET, digiFTD, spSET, digiTPC, tracking, forward, subset, full,
+ EvtSel = 'NONE',
+ EvtMax = Nevt,
+ ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc],
+ #ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc],
+ OutputLevel=INFO
+)
diff --git a/Generator/src/BeamBackgroundFileParserV1.cpp b/Generator/src/BeamBackgroundFileParserV1.cpp
index 7f1e5408..a39d2cd3 100644
--- a/Generator/src/BeamBackgroundFileParserV1.cpp
+++ b/Generator/src/BeamBackgroundFileParserV1.cpp
@@ -10,7 +10,8 @@ BeamBackgroundFileParserV1::BeamBackgroundFileParserV1(const std::string& filena
const std::string& treename,
double beam_energy,
double rate,
- double timewindow ) {
+ double timewindow,
+ int Nmcp) {
m_inputFile.reset(TFile::Open(filename.c_str(), "read"));
@@ -28,12 +29,12 @@ BeamBackgroundFileParserV1::BeamBackgroundFileParserV1(const std::string& filena
m_readTree->SetBranchAddress("dp", &dp);
m_readTree->SetBranchAddress("dz", &dz);
m_readTree->SetBranchAddress("pid", &pid);
- if(m_readTree->GetBranch("cosz")){
- m_readTree->SetBranchAddress("cosz", &cosz);
- }
+ m_readTree->SetBranchAddress("charge", &charge);
+ m_readTree->SetBranchAddress("cosz", &cosz);
m_rate = rate;
m_timewindow = timewindow;
+ m_Nmcp = Nmcp;
}
bool BeamBackgroundFileParserV1::load(IBeamBackgroundFileParser::BeamBackgroundData& result, int iEntry) {
@@ -47,24 +48,17 @@ bool BeamBackgroundFileParserV1::load(IBeamBackgroundFileParser::BeamBackgroundD
m_readTree->GetEntry(iEntry);
double p = m_beam_energy*(1+dp);
-
- // Now, we get a almost valid data
const double m2mm = 1e3; // convert from m to mm
- result.pdgid = pid;
- result.charge = (pid == 11) ? -1 : (pid == -11) ? 1 : -1;
result.x = x * m2mm;
result.y = y * m2mm;
result.z = (z+dz) * m2mm;
-
result.px = p * cosx;
result.py = p * cosy;
- if(m_readTree->GetBranch("cosz")){
- result.pz = p * cosz;}
- else{
- result.pz = p * std::sqrt(1-cosx*cosx-cosy*cosy);}
-
+ result.pz = p * cosz;
result.mass = 0.000511; // assume e-/e+, mass is 0.511 MeV
-
+ result.t = CLHEP::RandFlat::shoot(0., m_timewindow); // uniformly distributed within a time window, and ns
+ result.pdgid = pid;
+ result.charge = charge;
return true;
}
@@ -74,8 +68,11 @@ bool BeamBackgroundFileParserV1::load(IBeamBackgroundFileParser::BeamBackgroundD
bool BeamBackgroundFileParserV1::SampleParticleNum(int& npart, int& start ){
- npart = int(m_rate * m_timewindow);
- npart = CLHEP::RandPoisson::shoot(npart);
+ if(m_Nmcp==-1){
+ npart = int(m_rate * m_timewindow);
+ npart = CLHEP::RandPoisson::shoot(npart);
+ }
+ else npart = m_Nmcp;
int nTotPartInFile = m_readTree->GetEntries();
if(npart>nTotPartInFile){
diff --git a/Generator/src/BeamBackgroundFileParserV1.h b/Generator/src/BeamBackgroundFileParserV1.h
index 00527a2d..bba73ab1 100644
--- a/Generator/src/BeamBackgroundFileParserV1.h
+++ b/Generator/src/BeamBackgroundFileParserV1.h
@@ -8,7 +8,7 @@
class BeamBackgroundFileParserV1: public IBeamBackgroundFileParser {
public:
- BeamBackgroundFileParserV1(const std::string& filename, const std::string& treename, double beam_energy, double rate, double timewindow);
+ BeamBackgroundFileParserV1(const std::string& filename, const std::string& treename, double beam_energy, double rate, double timewindow, int Nmcp);
bool load(IBeamBackgroundFileParser::BeamBackgroundData&) { return 0; }
bool load(IBeamBackgroundFileParser::BeamBackgroundData&, int iEntry);
@@ -21,6 +21,7 @@ private:
double m_beam_energy;
double m_rate;
double m_timewindow;
+ int m_Nmcp;
double x, y, z, cosx, cosy, dz, dp, cosz;
int pid, charge;
diff --git a/Generator/src/BeamBackgroundFileParserV2.cpp b/Generator/src/BeamBackgroundFileParserV2.cpp
index 748ed9b8..f9de4bf5 100644
--- a/Generator/src/BeamBackgroundFileParserV2.cpp
+++ b/Generator/src/BeamBackgroundFileParserV2.cpp
@@ -6,7 +6,9 @@
BeamBackgroundFileParserV2::BeamBackgroundFileParserV2(const std::string& dirpath,
const std::string& treename,
- double beam_energy) {
+ double beam_energy,
+ double beam_time,
+ int Nmcp) {
//Choose one root file from the dir
std::vector<std::string> rootFiles;
@@ -20,18 +22,23 @@ BeamBackgroundFileParserV2::BeamBackgroundFileParserV2(const std::string& dirpat
int randomIndex = CLHEP::RandFlat::shoot(0., rootFiles.size()-1.);
std::cout << " Chosen file index "<<randomIndex<<", file name "<<rootFiles[randomIndex]<<std::endl;
m_inputFile.reset(TFile::Open(rootFiles[randomIndex].c_str(), "read"));
- std::cout << "BeamBackgroundFileParserV1: readin file name "<<m_inputFile->GetName()<<", status " << m_inputFile->IsOpen() << std::endl;
+ std::cout << "BeamBackgroundFileParserV2: readin file name "<<m_inputFile->GetName()<<", status " << m_inputFile->IsOpen() << std::endl;
m_readTree = (TTree*)m_inputFile->Get(treename.c_str());
std::cout << " Get Tree: entries " << m_readTree->GetEntries() << std::endl;
m_beam_energy = beam_energy;
-
- m_readTree->SetBranchAddress("X", &x);
- m_readTree->SetBranchAddress("Y", &y);
- m_readTree->SetBranchAddress("Z", &z);
- m_readTree->SetBranchAddress("PX", &px);
- m_readTree->SetBranchAddress("PY", &py);
- m_readTree->SetBranchAddress("DP", &dp);
- m_readTree->SetBranchAddress("PID", &pid);
+ m_beam_time = beam_time;
+ m_Nmcp = Nmcp;
+
+ m_readTree->SetBranchAddress("x", &x);
+ m_readTree->SetBranchAddress("y", &y);
+ m_readTree->SetBranchAddress("z", &z);
+ m_readTree->SetBranchAddress("cosx", &cosx);
+ m_readTree->SetBranchAddress("cosy", &cosy);
+ m_readTree->SetBranchAddress("dp", &dp);
+ m_readTree->SetBranchAddress("dz", &dz);
+ m_readTree->SetBranchAddress("pid", &pid);
+ m_readTree->SetBranchAddress("charge", &charge);
+ m_readTree->SetBranchAddress("cosz", &cosz);
}
@@ -47,17 +54,16 @@ bool BeamBackgroundFileParserV2::load(IBeamBackgroundFileParser::BeamBackgroundD
double p = m_beam_energy*(1+dp);
const double m2mm = 1e3; // convert from m to mm
- result.pdgid = (int)pid;
result.x = x * m2mm;
result.y = y * m2mm;
- result.z = z * m2mm;
-
- result.px = px;
- result.py = py;
- result.pz = sqrt(p*p - px*px - py*py);
-
- if(pid==22) result.mass = 0.;
- else result.mass = 0.000511; // assume e-/e+, mass is 0.511 MeV
+ result.z = (z+dz) * m2mm;
+ result.px = p * cosx;
+ result.py = p * cosy;
+ result.pz = p * cosz;
+ result.mass = 0.000511; // assume e-/e+, mass is 0.511 MeV
+ result.t = m_beam_time; // bunch time
+ result.pdgid = pid;
+ result.charge = charge;
return true;
@@ -67,8 +73,12 @@ bool BeamBackgroundFileParserV2::load(IBeamBackgroundFileParser::BeamBackgroundD
}
bool BeamBackgroundFileParserV2::SampleParticleNum(int& npart, int& start ){
- npart = m_readTree->GetEntries();
- start = 0;
+ if(m_Nmcp==-1){
+ npart = m_readTree->GetEntries();
+ }
+ else npart = m_Nmcp;
+
+ start = 0;
return true;
}
diff --git a/Generator/src/BeamBackgroundFileParserV2.h b/Generator/src/BeamBackgroundFileParserV2.h
index 0f9cb260..7799d94e 100644
--- a/Generator/src/BeamBackgroundFileParserV2.h
+++ b/Generator/src/BeamBackgroundFileParserV2.h
@@ -10,7 +10,7 @@
class BeamBackgroundFileParserV2: public IBeamBackgroundFileParser {
public:
- BeamBackgroundFileParserV2(const std::string& dirname, const std::string& treename, double beam_energy);
+ BeamBackgroundFileParserV2(const std::string& dirname, const std::string& treename, double beam_energy, double beam_time, int Nmcp);
bool load(IBeamBackgroundFileParser::BeamBackgroundData&) { return 0; }
bool load(IBeamBackgroundFileParser::BeamBackgroundData&, int iEntry);
@@ -20,9 +20,11 @@ public:
private:
std::unique_ptr<TFile> m_inputFile;
TTree* m_readTree;
- double m_beam_energy;
+ double m_beam_energy, m_beam_time;;
+ int m_Nmcp;
- float x, y, z, px, py, dp, pid;
+ double x, y, z, cosx, cosy, dz, dp, cosz;
+ int pid, charge;
};
diff --git a/Generator/src/GtBeamBackgroundTool.cpp b/Generator/src/GtBeamBackgroundTool.cpp
index 82847a23..cd8a3727 100644
--- a/Generator/src/GtBeamBackgroundTool.cpp
+++ b/Generator/src/GtBeamBackgroundTool.cpp
@@ -136,8 +136,12 @@ bool GtBeamBackgroundTool::init_BeamBackgroundFileParserV1(const std::string& la
error() << "Did not find the beam process rate for: " << label << endmsg;
return false;
}
-
- m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV1>(inputfn, "BeamTree", m_Ebeam, itRate->second, m_timewindow);
+ auto itNmcp = m_Nmcpmaps.find(label);
+ if(itNmcp == m_Nmcpmaps.end()){
+ error() << "Did not find the beam process Nmcp for: " << label << endmsg;
+ return false;
+ }
+ m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV1>(inputfn, "BeamTree", m_Ebeam, itRate->second, m_timewindow, itNmcp->second);
return true;
}
@@ -146,7 +150,17 @@ bool GtBeamBackgroundTool::init_BeamBackgroundFileParserV2(const std::string& la
const std::string& inputfn) {
info() << "Initializing beam background ... " << label << inputfn <<endmsg;
- m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV2>(inputfn, "RBBG", m_Ebeam);
+ auto itTime = m_timebkgmaps.find(label);
+ if(itTime == m_timebkgmaps.end()){
+ error() << "Did not find the beam process time for: " << label << endmsg;
+ return false;
+ }
+ auto itNmcp = m_Nmcpmaps.find(label);
+ if(itNmcp == m_Nmcpmaps.end()){
+ error() << "Did not find the beam process Nmcp for: " << label << endmsg;
+ return false;
+ }
+ m_beaminputs[label] = std::make_shared<BeamBackgroundFileParserV2>(inputfn, "BeamTree", m_Ebeam, itTime->second, itNmcp->second);
return true;
}
diff --git a/Generator/src/GtBeamBackgroundTool.h b/Generator/src/GtBeamBackgroundTool.h
index 6affcb85..40a6ce66 100644
--- a/Generator/src/GtBeamBackgroundTool.h
+++ b/Generator/src/GtBeamBackgroundTool.h
@@ -67,6 +67,12 @@ private:
// Rotation along Y for single beam background. Unit: rad
Gaudi::Property<std::map<std::string, double>> m_rotYMap {this, "RotationAlongYMap"};
+ // Time of bunch crossing. Unit: ns, consistent with simulation
+ Gaudi::Property<std::map<std::string, double>> m_timebkgmaps{this, "TimeBkgMap"};
+
+ // Number of McParticles in different beambkg. -1: pair use one file and single beam use rate*time; >=0: fixed number
+ Gaudi::Property<std::map<std::string, int>> m_Nmcpmaps{this, "NumberMcParticle"};
+
private:
std::map<std::string, std::shared_ptr<IBeamBackgroundFileParser>> m_beaminputs;
diff --git a/Generator/src/GtGunTool.cpp b/Generator/src/GtGunTool.cpp
index 7f9b56a3..48697c0b 100644
--- a/Generator/src/GtGunTool.cpp
+++ b/Generator/src/GtGunTool.cpp
@@ -95,6 +95,15 @@ GtGunTool::initialize() {
return StatusCode::FAILURE;
}
+ // Time
+ if (m_times.value().size()==0){
+ for(int i=0; i<m_particles.value().size(); i++) m_times.value().push_back(0);
+ }
+ else if (m_times.value().size() != m_particles.value().size()) {
+ error() << "Mismatched times and particles." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
return sc;
}
@@ -163,7 +172,9 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
mcp.setGeneratorStatus(1);
mcp.setSimulatorStatus(1);
mcp.setCharge(static_cast<float>(charge));
- mcp.setTime(0.0);
+ mcp.setTime(m_times.value()[i]);
+ //if (i<m_times.value().size()) { mcp.setTime(m_times.value()[i]); }
+ //else { mcp.setTime(0.0); }
mcp.setMass(mass);
// Unit is mm
diff --git a/Generator/src/GtGunTool.h b/Generator/src/GtGunTool.h
index d4b8d9f8..0f3f7716 100644
--- a/Generator/src/GtGunTool.h
+++ b/Generator/src/GtGunTool.h
@@ -57,6 +57,9 @@ private:
Gaudi::Property<std::vector<double>> m_phimins{this, "PhiMins"};
Gaudi::Property<std::vector<double>> m_phimaxs{this, "PhiMaxs"};
+ // For time
+ Gaudi::Property<std::vector<double>> m_times{this, "Times"};
+
};
diff --git a/Generator/src/StdHepRdr.cpp b/Generator/src/StdHepRdr.cpp
index 4255e4bd..2eccb9f1 100644
--- a/Generator/src/StdHepRdr.cpp
+++ b/Generator/src/StdHepRdr.cpp
@@ -48,7 +48,7 @@ bool StdHepRdr::mutate(MyHepMC::GenEvent& event){
mcp.setGeneratorStatus (mc->getGeneratorStatus());
mcp.setSimulatorStatus (mc->getSimulatorStatus());
mcp.setCharge (mc->getCharge());
- mcp.setTime (mc->getTime());
+ mcp.setTime (m_starttime + mc->getTime());
mcp.setMass (mc->getMass());
mcp.setVertex (mc->getVertex());
mcp.setEndpoint (mc->getEndpoint());
diff --git a/Generator/src/StdHepRdr.h b/Generator/src/StdHepRdr.h
index 7517430b..4796eeda 100644
--- a/Generator/src/StdHepRdr.h
+++ b/Generator/src/StdHepRdr.h
@@ -34,6 +34,7 @@ private:
// input file name
Gaudi::Property<std::string> m_filename{this, "Input"};
+ Gaudi::Property<double> m_starttime{this, "StartTime", 0, "Default start time"};
};
--
GitLab