Skip to content

Commit e0a2ea5

Browse files
njacazionjacazio
authored
Extend fast sim for nuclei + safeties (#97)
- Add safety when creating LUTs - Cleaning python files - Add triton for LUT making - Throw error if particle is available - Add config for nuclei
1 parent ca009d9 commit e0a2ea5

File tree

10 files changed

+98
-66
lines changed

10 files changed

+98
-66
lines changed

‎examples/aod/createO2tables.C‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,9 @@ int createO2tables(const char* inputFile = "delphes.root",
7171
TDatabasePDG::Instance()->AddParticle("deuteron", "deuteron", 1.8756134, kTRUE, 0.0, 3, "Nucleus", 1000010020);
7272
TDatabasePDG::Instance()->AddAntiParticle("anti-deuteron", -1000010020);
7373

74+
TDatabasePDG::Instance()->AddParticle("triton", "triton", 2.8089218, kTRUE, 0.0, 3, "Nucleus", 1000010030);
75+
TDatabasePDG::Instance()->AddAntiParticle("anti-triton", -1000010030);
76+
7477
TDatabasePDG::Instance()->AddParticle("helium3", "helium3", 2.80839160743, kTRUE, 0.0, 6, "Nucleus", 1000020030);
7578
TDatabasePDG::Instance()->AddAntiParticle("anti-helium3", -1000020030);
7679

@@ -102,7 +105,8 @@ int createO2tables(const char* inputFile = "delphes.root",
102105
mapPdgLut.insert(std::make_pair(2212, "lutCovm.pr.dat"));
103106
if (enable_nuclei) {
104107
mapPdgLut.insert(std::make_pair(1000010020, "lutCovm.de.dat"));
105-
mapPdgLut.insert(std::make_pair(1000020030, "lutCovm.he.dat"));
108+
mapPdgLut.insert(std::make_pair(1000010030, "lutCovm.tr.dat"));
109+
mapPdgLut.insert(std::make_pair(1000020030, "lutCovm.he3.dat"));
106110
}
107111
for (auto e : mapPdgLut) {
108112
if (!smearer.loadTable(e.first, e.second)) {

‎examples/pythia8/pythia_nuclei.cfg‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
#Config file to define the nuclei species that are not in vanilla pythia
2+
1000020030:all 3He 3Hebar 1 6 0 2.8094

‎examples/scripts/clean.sh‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,3 +17,4 @@ rm tmpscript*.sh
1717
rm lutWrite*.cc
1818
rm DetectorK.*
1919
rm -r fwdRes
20+
rm -r __pycache__

‎examples/scripts/createO2tables.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ def do_copy(in_file, out_file=".", in_path=None):
107107
lut_tag = opt("lut_tag")
108108
lut_particles = ["el", "mu", "pi", "ka", "pr"]
109109
if use_nuclei:
110-
lut_particles += ["de", "he3"]
110+
lut_particles += ["de", "tr", "he3"]
111111
if create_luts:
112112
# Creating LUTs
113113
minimum_track_radius = opt("minimum_track_radius")

‎examples/scripts/create_luts.sh‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -84,22 +84,32 @@ for i in $PARTICLES; do
8484
TDatabasePDG::Instance()->AddParticle("deuteron", "deuteron", 1.8756134, kTRUE, 0.0, 3, "Nucleus", 1000010020);
8585
TDatabasePDG::Instance()->AddAntiParticle("anti-deuteron", -1000010020);
8686
87+
TDatabasePDG::Instance()->AddParticle("triton", "triton", 2.8089218, kTRUE, 0.0, 3, "Nucleus", 1000010030);
88+
TDatabasePDG::Instance()->AddAntiParticle("anti-triton", -1000010030);
89+
8790
TDatabasePDG::Instance()->AddParticle("helium3", "helium3", 2.80839160743, kTRUE, 0.0, 6, "Nucleus", 1000020030);
8891
TDatabasePDG::Instance()->AddAntiParticle("anti-helium3", -1000020030);
8992
90-
const TString pn[7] = {"el", "mu", "pi", "ka", "pr", "de", "he3"};
91-
const int pc[7] = {11, 13, 211, 321, 2212, 1000010020, 1000020030 };
93+
const int N = 8;
94+
const TString pn[N] = {"el", "mu", "pi", "ka", "pr", "de", "tr", "he3"};
95+
const int pc[N] = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030 };
9296
const float field = ${FIELD}f;
9397
const float rmin = ${RMIN};
9498
const int i = ${i};
95-
lutWrite_${WHAT}("${OUT_PATH}/lutCovm." + pn[i] + "${OUT_TAG}.dat", pc[i], field, rmin);
99+
const TString out_file = "${OUT_PATH}/lutCovm." + pn[i] + "${OUT_TAG}.dat";
100+
Printf("Creating LUT for particle ID %i: %s with pdg code %i to %s", i, pn[i].Data(), pc[i], out_file.Data());
101+
if(i >= 0 && i < N){
102+
lutWrite_${WHAT}(out_file, pc[i], field, rmin);
103+
} else{
104+
Printf("Particle ID %i is too large or too small", i);
105+
}
96106
97107
EOF
98108
done
99109

100110
# Checking that the output LUTs are OK
101111
NullSize=""
102-
P=(el mu pi ka pr de he3)
112+
P=(el mu pi ka pr de tr he3)
103113
for i in $PARTICLES; do
104114
LUT_FILE=${OUT_PATH}/lutCovm.${P[$i]}${OUT_TAG}.dat
105115
if [[ ! -s ${LUT_FILE} ]]; then

‎examples/scripts/default_configfile.ini‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,17 @@ custom_gen = INPUT_FILES /tmp/AnalysisResults_*.hepmc
7676
[GUN]
7777
custom_gen = rpythia8-gun --pdg 421 --px 1. --py 0. --pz 0. --xProd 0. --yProd 0. --zProd 0. --config $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg --decay
7878

79+
[BOX_pion]
80+
custom_gen = rpythia8-box --pdg 211 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100
81+
82+
[BOX_proton]
83+
custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100
84+
7985
[BOX_deuteron]
80-
custom_gen = rpythia8-box --pdg 1000010020 --etamin -2. --etamax 2. --phimin 0. --phimax 3.14 --pmin 100 --pmax 1000 --xProd 0. --yProd 0. --zProd 0. --config $DELPHESO2_ROOT/examples/pythia8/decays/force_hadronic_D.cfg --decay --npart 100
86+
custom_gen = rpythia8-box --pdg 1000010020 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100
87+
88+
[BOX_helium3]
89+
custom_gen = rpythia8-box --pdg 1000020030 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg
8190

8291
[GUN_Lc_pKpi]
8392
custom_gen = rpythia8-gun --pdg 4122 --px 1. --py 0. --pz 0. --xProd 1. --yProd 0. --zProd 0. --config $O2DPG_ROOT/MC/config/PWGHF/pythia8/decayer/force_hadronic_D_forceLcChannel1.cfg --decay

‎rpythia8/.clang-format‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
../examples/aod/.clang-format

‎rpythia8/rpythia8-box.cc‎

Lines changed: 31 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ int main(int argc, char **argv)
1818
int npart;
1919
double xProd, yProd, zProd;
2020
bool verbose, decay;
21-
21+
2222
/** process arguments **/
2323
namespace po = boost::program_options;
2424
po::options_description desc("Options");
@@ -44,29 +44,40 @@ int main(int argc, char **argv)
4444
("verbose,V", po::bool_switch(&verbose)->default_value(false), "Verbose event listing")
4545
("seed", po::value<int>(&seed)->default_value(1), "initial seed")
4646
;
47-
47+
4848
po::variables_map vm;
4949
po::store(po::parse_command_line(argc, argv, desc), vm);
5050
po::notify(vm);
51-
51+
5252
if (vm.count("help")) {
5353
std::cout << desc << std::endl;
5454
return 1;
5555
}
56-
}
57-
catch(std::exception& e) {
56+
} catch (std::exception& e) {
5857
std::cerr << "Error: " << e.what() << std::endl;
5958
std::cout << desc << std::endl;
6059
return 1;
6160
}
62-
61+
6362
HepMC::Pythia8ToHepMC ToHepMC;
6463
HepMC::IO_GenEvent ascii_io(output, std::ios::out);
65-
64+
6665
// pythia
6766
Pythia pythia;
6867

68+
// configure pythia
69+
pythia.readString("ProcessLevel:all = off");
70+
// pythia.readString("SoftQCD:elastic on");
71+
if (!config.empty() && !pythia.readFile(config)) {
72+
std::cout << "Error: could not read config file \"" << config << "\"" << std::endl;
73+
return 1;
74+
}
75+
6976
// check valid pdg code
77+
if (!pythia.particleData.isParticle(pdg)) {
78+
std::cout << "Error: invalid PDG code \"" << pdg << "\" is not in the particle list" << std::endl;
79+
return 1;
80+
}
7081
if (!pythia.particleData.isLepton(pdg) &&
7182
!pythia.particleData.isHadron(pdg) &&
7283
!pythia.particleData.isResonance(pdg)) {
@@ -78,15 +89,7 @@ int main(int argc, char **argv)
7889
}
7990
}
8091

81-
// config
82-
pythia.readString("ProcessLevel:all = off");
83-
// pythia.readString("SoftQCD:elastic on");
84-
if (!config.empty() && !pythia.readFile(config)) {
85-
std::cout << "Error: could not read config file \"" << config << "\"" << std::endl;
86-
return 1;
87-
}
88-
89-
std::cout<<"Random:seed =" + std::to_string(seed)<<std::endl;
92+
std::cout << "Random:seed =" + std::to_string(seed) << std::endl;
9093
pythia.readString("Random:setSeed = on");
9194
pythia.readString("Random:seed =" + std::to_string(seed));
9295
// init
@@ -103,7 +106,7 @@ int main(int argc, char **argv)
103106
particle.zProd(zProd);
104107

105108
// background interface
106-
Pythia8::Pythia *pythia_bkg = nullptr;
109+
Pythia8::Pythia* pythia_bkg = nullptr;
107110
if (!background_config.empty()) {
108111
std::cout << "Background: configure from " << background_config << std::endl;
109112
pythia_bkg = new Pythia8::Pythia;
@@ -113,17 +116,15 @@ int main(int argc, char **argv)
113116
}
114117
pythia_bkg->init();
115118
}
116-
117119

118-
119120
// event loop
120121
double eta, phi, p, pt, pl, e;
121-
srand (time(NULL));
122+
srand(time(NULL));
122123
for (int iev = 0; iev < nevents; ++iev) {
123124

124125
// reset, add particle and decay
125126
pythia.event.reset();
126-
for(int ipart = 0; ipart < npart; ipart++){
127+
for (int ipart = 0; ipart < npart; ipart++) {
127128
eta = eta_min + static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / (eta_max - eta_min)));
128129
phi = phi_min + static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / (phi_max - phi_min)));
129130
p = p_min + static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / (p_max - p_min)));
@@ -135,21 +136,23 @@ int main(int argc, char **argv)
135136
particle.py(pt * sin(phi));
136137
particle.pz(pl);
137138
pythia.event.append(particle);
138-
if (decay) pythia.moreDecays();
139+
if (decay)
140+
pythia.moreDecays();
139141
pythia.next();
140142
}
141-
143+
142144
// print verbose
143-
if (verbose) pythia.event.list(1);
145+
if (verbose)
146+
pythia.event.list(1);
144147

145148
// background
146149
if (pythia_bkg) {
147150
pythia_bkg->next();
148151
pythia.event += pythia_bkg->event;
149-
}
150-
152+
}
153+
151154
// write HepMC
152-
HepMC::GenEvent *hepmcevt = new HepMC::GenEvent();
155+
HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
153156
ToHepMC.fill_next_event(pythia, hepmcevt);
154157
ascii_io << hepmcevt;
155158
delete hepmcevt;
@@ -161,7 +164,6 @@ int main(int argc, char **argv)
161164
pythia_bkg->stat();
162165
delete pythia_bkg;
163166
}
164-
167+
165168
return 0;
166-
167169
}

‎rpythia8/rpythia8-gun.cc‎

Lines changed: 30 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -6,15 +6,15 @@
66

77
using namespace Pythia8;
88

9-
int main(int argc, char **argv)
9+
int main(int argc, char** argv)
1010
{
1111

1212
int nevents, pdg, seed;
1313
std::string config, output, background_config;
1414
double px, py, pz;
1515
double xProd, yProd, zProd;
1616
bool verbose, decay;
17-
17+
1818
/** process arguments **/
1919
namespace po = boost::program_options;
2020
po::options_description desc("Options");
@@ -36,29 +36,40 @@ int main(int argc, char **argv)
3636
("verbose,V", po::bool_switch(&verbose)->default_value(false), "Verbose event listing")
3737
("seed", po::value<int>(&seed)->default_value(1), "initial seed")
3838
;
39-
39+
4040
po::variables_map vm;
4141
po::store(po::parse_command_line(argc, argv, desc), vm);
4242
po::notify(vm);
43-
43+
4444
if (vm.count("help")) {
4545
std::cout << desc << std::endl;
4646
return 1;
4747
}
48-
}
49-
catch(std::exception& e) {
48+
} catch (std::exception& e) {
5049
std::cerr << "Error: " << e.what() << std::endl;
5150
std::cout << desc << std::endl;
5251
return 1;
5352
}
54-
53+
5554
HepMC::Pythia8ToHepMC ToHepMC;
5655
HepMC::IO_GenEvent ascii_io(output, std::ios::out);
57-
56+
5857
// pythia
5958
Pythia pythia;
6059

60+
// configure pythia
61+
pythia.readString("ProcessLevel:all = off");
62+
// pythia.readString("SoftQCD:elastic on");
63+
if (!config.empty() && !pythia.readFile(config)) {
64+
std::cout << "Error: could not read config file \"" << config << "\"" << std::endl;
65+
return 1;
66+
}
67+
6168
// check valid pdg code
69+
if (!pythia.particleData.isParticle(pdg)) {
70+
std::cout << "Error: invalid PDG code \"" << pdg << "\" is not in the particle list" << std::endl;
71+
return 1;
72+
}
6273
if (!pythia.particleData.isLepton(pdg) &&
6374
!pythia.particleData.isHadron(pdg) &&
6475
!pythia.particleData.isResonance(pdg)) {
@@ -70,15 +81,7 @@ int main(int argc, char **argv)
7081
}
7182
}
7283

73-
// config
74-
pythia.readString("ProcessLevel:all = off");
75-
// pythia.readString("SoftQCD:elastic on");
76-
if (!config.empty() && !pythia.readFile(config)) {
77-
std::cout << "Error: could not read config file \"" << config << "\"" << std::endl;
78-
return 1;
79-
}
80-
81-
std::cout<<"Random:seed =" + std::to_string(seed)<<std::endl;
84+
std::cout << "Random:seed =" + std::to_string(seed) << std::endl;
8285
pythia.readString("Random:setSeed = on");
8386
pythia.readString("Random:seed =" + std::to_string(seed));
8487
// init
@@ -100,7 +103,7 @@ int main(int argc, char **argv)
100103
particle.zProd(zProd);
101104

102105
// background interface
103-
Pythia8::Pythia *pythia_bkg = nullptr;
106+
Pythia8::Pythia* pythia_bkg = nullptr;
104107
if (!background_config.empty()) {
105108
std::cout << "Background: configure from " << background_config << std::endl;
106109
pythia_bkg = new Pythia8::Pythia;
@@ -110,29 +113,29 @@ int main(int argc, char **argv)
110113
}
111114
pythia_bkg->init();
112115
}
113-
114116

115-
116117
// event loop
117118
for (int iev = 0; iev < nevents; ++iev) {
118119

119120
// reset, add particle and decay
120121
pythia.event.reset();
121122
pythia.event.append(particle);
122-
if (decay) pythia.moreDecays();
123+
if (decay)
124+
pythia.moreDecays();
123125
pythia.next();
124-
126+
125127
// print verbose
126-
if (verbose) pythia.event.list(1);
128+
if (verbose)
129+
pythia.event.list(1);
127130

128131
// background
129132
if (pythia_bkg) {
130133
pythia_bkg->next();
131134
pythia.event += pythia_bkg->event;
132-
}
133-
135+
}
136+
134137
// write HepMC
135-
HepMC::GenEvent *hepmcevt = new HepMC::GenEvent();
138+
HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
136139
ToHepMC.fill_next_event(pythia, hepmcevt);
137140
ascii_io << hepmcevt;
138141
delete hepmcevt;
@@ -144,7 +147,6 @@ int main(int argc, char **argv)
144147
pythia_bkg->stat();
145148
delete pythia_bkg;
146149
}
147-
150+
148151
return 0;
149-
150152
}

‎src/TrackSmearer.hh‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,15 +40,16 @@ public:
4040
case 321: return 3; // Kaon
4141
case 2212: return 4; // Proton
4242
case 1000010020: return 5; // Deuteron
43-
case 1000020030: return 6; // Helium3
43+
case 1000010030: return 6; // Triton
44+
case 1000020030: return 7; // Helium3
4445
default: return 2; // Default: pion
4546
};
4647
};
4748

4849
void setdNdEta(float val) { mdNdEta = val; };
4950

5051
protected:
51-
static constexpr unsigned int nLUTs = 7;
52+
static constexpr unsigned int nLUTs = 8; // Number of LUT available
5253
lutHeader_t *mLUTHeader[nLUTs] = {nullptr};
5354
lutEntry_t *****mLUTEntry[nLUTs] = {nullptr};
5455
bool mUseEfficiency = true;

0 commit comments

Comments
 (0)