Skip to content

Commit fe8c5cf

Browse files
njacazionjacazio
authored
Extend simulation with nuclei + general improvements (#98)
- Update drawing script - Saving pngs for all particles - Add triton to nuclei pythia config - Add several clang format - Update LUT reader to work with large sized ones - Merging with the same number of jobs - In LUT making: using negative mass when charge is > 2 - Add triton in config file
1 parent e0a2ea5 commit fe8c5cf

File tree

11 files changed

+128
-87
lines changed

11 files changed

+128
-87
lines changed

‎examples/pythia8/pythia_nuclei.cfg‎

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

‎examples/scripts/default_configfile.ini‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,9 @@ custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phim
8585
[BOX_deuteron]
8686
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
8787

88+
[BOX_triton]
89+
custom_gen = rpythia8-box --pdg 1000010030 --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
90+
8891
[BOX_helium3]
8992
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
9093

‎examples/scripts/diagnostic_tools/doanalysis.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -241,7 +241,7 @@ def build_list_of_files(file_list):
241241
"is already found, remove it before merging, you can use the --mergeonly flag to avoid running the analysis again")
242242
continue
243243
merged_files.append(merged_file)
244-
run_cmd(f"hadd -f {merged_file} " +
244+
run_cmd(f"hadd -j {njobs} -f {merged_file} " +
245245
" ".join(files_per_type[i]))
246246
if len(merged_files) == 0:
247247
warning_msg("Merged no files")

‎examples/smearing/.clang-format‎

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

‎examples/smearing/draw.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -207,10 +207,11 @@ def draw_alice3_label(x=0.5, y=0.9):
207207
if tag_name is not None:
208208
drawn.append(latex.DrawLatexNDC(0.5, 0.80, tag_name[0]))
209209
drawn.append(latex.DrawLatexNDC(0.5, 0.75, f"#eta = {int(eta)}" +
210-
" R_{min} = " + rmin))
210+
(" R_{min} = " + rmin if rmin is not None else "") ))
211211

212212
adjust_pad()
213213
canvas_all_species.Update()
214+
canvas_all_species.SaveAs(f"/tmp/{canvas_all_species.GetName()}.png")
214215
if not background:
215216
input("Done, press enter to continue")
216217

‎examples/smearing/lutRead_dca.C‎

Lines changed: 37 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,58 +1,66 @@
11
#include "lutCovm.hh"
22

3-
TGraph *
4-
lutRead_dca(const char *filename = "lutCovm.dat", double eta = 0.)
3+
TGraph* lutRead_dca(const char* filename = "lutCovm.dat",
4+
double eta = 0.,
5+
double nch = 100.)
56
{
6-
7+
78
// input file
89
ifstream lutFile(filename, std::ofstream::binary);
910

1011
// read header
1112
lutHeader_t lutHeader;
12-
lutFile.read(reinterpret_cast<char *>(&lutHeader), sizeof(lutHeader));
13+
lutFile.read(reinterpret_cast<char*>(&lutHeader), sizeof(lutHeader));
1314
lutHeader.print();
1415

1516
// entries
1617
const int nnch = lutHeader.nchmap.nbins;
1718
const int nrad = lutHeader.radmap.nbins;
1819
const int neta = lutHeader.etamap.nbins;
1920
const int npt = lutHeader.ptmap.nbins;
20-
lutEntry_t lutTable[nnch][nrad][neta][npt];
21-
21+
lutEntry_t lutTable;
22+
23+
const int nch_bin = lutHeader.nchmap.find(nch);
24+
const int rad_bin = lutHeader.radmap.find(0.);
25+
const int eta_bin = lutHeader.etamap.find(eta);
26+
27+
// create graph of pt resolution at eta = 0
28+
auto gpt = new TGraph();
29+
gpt->SetName(filename);
30+
gpt->SetTitle(filename);
31+
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
32+
gpt->GetYaxis()->SetTitle("pointing resolution (#mum)");
33+
2234
// read entries
2335
for (int inch = 0; inch < nnch; ++inch) {
2436
for (int irad = 0; irad < nrad; ++irad) {
2537
for (int ieta = 0; ieta < neta; ++ieta) {
2638
for (int ipt = 0; ipt < npt; ++ipt) {
27-
lutFile.read(reinterpret_cast<char *>(&lutTable[inch][irad][ieta][ipt]), sizeof(lutEntry_t));
28-
// lutTable[inch][irad][ieta][ipt].print();
39+
lutFile.read(reinterpret_cast<char*>(&lutTable), sizeof(lutEntry_t));
40+
// lutTable.print();
41+
auto lutEntry = &lutTable;
42+
if (!lutEntry->valid) {
43+
std::cout << " ipt = " << ipt << " is not valid " << std::endl;
44+
continue;
45+
}
46+
if (nch_bin != inch) {
47+
continue;
48+
}
49+
if (eta_bin != ieta) {
50+
continue;
51+
}
52+
if (rad_bin != irad) {
53+
continue;
54+
}
55+
auto cen = lutEntry->pt;
56+
auto val = sqrt(lutEntry->covm[0]) * 1.e4;
57+
gpt->SetPoint(gpt->GetN(), cen, val);
2958
}
3059
}
3160
}
3261
}
3362

3463
lutFile.close();
3564

36-
// create graph of pt resolution at eta = 0
37-
auto inch = lutHeader.nchmap.find(0.);
38-
auto irad = lutHeader.radmap.find(0.);
39-
auto ieta = lutHeader.etamap.find(eta);
40-
auto gpt = new TGraph();
41-
gpt->SetName(filename);
42-
gpt->SetTitle(filename);
43-
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
44-
gpt->GetYaxis()->SetTitle("pointing resolution (#mum)");
45-
for (int ipt = 0; ipt < npt; ++ipt) {
46-
auto lutEntry = &lutTable[inch][irad][ieta][ipt];
47-
if (!lutEntry->valid) {
48-
std::cout << " ipt = " << ipt << " is not valid " << std::endl;
49-
continue;
50-
}
51-
auto cen = lutEntry->pt;
52-
auto val = sqrt(lutEntry->covm[0]) * 1.e4;
53-
gpt->SetPoint(gpt->GetN(), cen, val);
54-
}
55-
5665
return gpt;
57-
5866
}

‎examples/smearing/lutRead_eff.C‎

Lines changed: 37 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,58 +1,66 @@
11
#include "lutCovm.hh"
22

3-
TGraph *
4-
lutRead_eff(const char *filename = "lutCovm.dat", double eta = 0.)
3+
TGraph* lutRead_eff(const char* filename = "lutCovm.dat",
4+
double eta = 0.,
5+
double nch = 100.)
56
{
6-
7+
78
// input file
89
ifstream lutFile(filename, std::ofstream::binary);
910

1011
// read header
1112
lutHeader_t lutHeader;
12-
lutFile.read(reinterpret_cast<char *>(&lutHeader), sizeof(lutHeader));
13+
lutFile.read(reinterpret_cast<char*>(&lutHeader), sizeof(lutHeader));
1314
lutHeader.print();
1415

1516
// entries
1617
const int nnch = lutHeader.nchmap.nbins;
1718
const int nrad = lutHeader.radmap.nbins;
1819
const int neta = lutHeader.etamap.nbins;
1920
const int npt = lutHeader.ptmap.nbins;
20-
lutEntry_t lutTable[nnch][nrad][neta][npt];
21-
21+
lutEntry_t lutTable;
22+
23+
const int nch_bin = lutHeader.nchmap.find(nch);
24+
const int rad_bin = lutHeader.radmap.find(0.);
25+
const int eta_bin = lutHeader.etamap.find(eta);
26+
27+
// create graph of pt resolution at eta = 0
28+
auto gpt = new TGraph();
29+
gpt->SetName(filename);
30+
gpt->SetTitle(filename);
31+
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
32+
gpt->GetYaxis()->SetTitle("efficiency (%)");
33+
2234
// read entries
2335
for (int inch = 0; inch < nnch; ++inch) {
2436
for (int irad = 0; irad < nrad; ++irad) {
2537
for (int ieta = 0; ieta < neta; ++ieta) {
2638
for (int ipt = 0; ipt < npt; ++ipt) {
27-
lutFile.read(reinterpret_cast<char *>(&lutTable[inch][irad][ieta][ipt]), sizeof(lutEntry_t));
28-
// lutTable[inch][irad][ieta][ipt].print();
39+
lutFile.read(reinterpret_cast<char*>(&lutTable), sizeof(lutEntry_t));
40+
// lutTable.print();
41+
auto lutEntry = &lutTable;
42+
if (!lutEntry->valid) {
43+
std::cout << " ipt = " << ipt << " is not valid " << std::endl;
44+
continue;
45+
}
46+
if (nch_bin != inch) {
47+
continue;
48+
}
49+
if (eta_bin != ieta) {
50+
continue;
51+
}
52+
if (rad_bin != irad) {
53+
continue;
54+
}
55+
auto cen = lutEntry->pt;
56+
auto val = lutEntry->eff * 1.e2;
57+
gpt->SetPoint(gpt->GetN(), cen, val);
2958
}
3059
}
3160
}
3261
}
3362

3463
lutFile.close();
3564

36-
// create graph of pt resolution at eta = 0
37-
auto inch = lutHeader.nchmap.find(0.);
38-
auto irad = lutHeader.radmap.find(0.);
39-
auto ieta = lutHeader.etamap.find(eta);
40-
auto gpt = new TGraph();
41-
gpt->SetName(filename);
42-
gpt->SetTitle(filename);
43-
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
44-
gpt->GetYaxis()->SetTitle("efficiency (%)");
45-
for (int ipt = 0; ipt < npt; ++ipt) {
46-
auto lutEntry = &lutTable[inch][irad][ieta][ipt];
47-
if (!lutEntry->valid) {
48-
std::cout << " ipt = " << ipt << " is not valid " << std::endl;
49-
continue;
50-
}
51-
auto cen = lutEntry->pt;
52-
auto val = lutEntry->eff * 1.e2;
53-
gpt->SetPoint(gpt->GetN(), cen, val);
54-
}
55-
5665
return gpt;
57-
5866
}

‎examples/smearing/lutRead_pt.C‎

Lines changed: 37 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -1,55 +1,66 @@
11
#include "lutCovm.hh"
22

3-
TGraph *
4-
lutRead_pt(const char *filename = "lutCovm.dat", double eta = 0.)
3+
TGraph* lutRead_pt(const char* filename = "lutCovm.dat",
4+
double eta = 0.,
5+
double nch = 100.)
56
{
6-
7+
78
// input file
89
ifstream lutFile(filename, std::ofstream::binary);
910

1011
// read header
1112
lutHeader_t lutHeader;
12-
lutFile.read(reinterpret_cast<char *>(&lutHeader), sizeof(lutHeader));
13+
lutFile.read(reinterpret_cast<char*>(&lutHeader), sizeof(lutHeader));
1314
lutHeader.print();
1415

1516
// entries
1617
const int nnch = lutHeader.nchmap.nbins;
1718
const int nrad = lutHeader.radmap.nbins;
1819
const int neta = lutHeader.etamap.nbins;
1920
const int npt = lutHeader.ptmap.nbins;
20-
lutEntry_t lutTable[nnch][nrad][neta][npt];
21-
21+
lutEntry_t lutTable;
22+
23+
const int nch_bin = lutHeader.nchmap.find(nch);
24+
const int rad_bin = lutHeader.radmap.find(0.);
25+
const int eta_bin = lutHeader.etamap.find(eta);
26+
27+
// create graph of pt resolution at eta = 0
28+
auto gpt = new TGraph();
29+
gpt->SetName(filename);
30+
gpt->SetTitle(filename);
31+
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
32+
gpt->GetYaxis()->SetTitle("momentum resolution (%)");
33+
2234
// read entries
2335
for (int inch = 0; inch < nnch; ++inch) {
2436
for (int irad = 0; irad < nrad; ++irad) {
2537
for (int ieta = 0; ieta < neta; ++ieta) {
2638
for (int ipt = 0; ipt < npt; ++ipt) {
27-
lutFile.read(reinterpret_cast<char *>(&lutTable[inch][irad][ieta][ipt]), sizeof(lutEntry_t));
28-
// lutTable[inch][irad][ieta][ipt].print();
39+
lutFile.read(reinterpret_cast<char*>(&lutTable), sizeof(lutEntry_t));
40+
// lutTable.print();
41+
auto lutEntry = &lutTable;
42+
if (!lutEntry->valid) {
43+
std::cout << " ipt = " << ipt << " is not valid " << std::endl;
44+
continue;
45+
}
46+
if (nch_bin != inch) {
47+
continue;
48+
}
49+
if (eta_bin != ieta) {
50+
continue;
51+
}
52+
if (rad_bin != irad) {
53+
continue;
54+
}
55+
auto cen = lutEntry->pt;
56+
auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.;
57+
gpt->SetPoint(gpt->GetN(), cen, val);
2958
}
3059
}
3160
}
3261
}
3362

3463
lutFile.close();
3564

36-
// create graph of pt resolution at eta = 0
37-
auto inch = lutHeader.nchmap.find(0.);
38-
auto irad = lutHeader.radmap.find(0.);
39-
auto ieta = lutHeader.etamap.find(eta);
40-
auto gpt = new TGraph();
41-
gpt->SetName(filename);
42-
gpt->SetTitle(filename);
43-
gpt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
44-
gpt->GetYaxis()->SetTitle("momentum resolution (%)");
45-
for (int ipt = 0; ipt < npt; ++ipt) {
46-
auto lutEntry = &lutTable[inch][irad][ieta][ipt];
47-
if (!lutEntry->valid) continue;
48-
auto cen = lutEntry->pt;
49-
auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.;
50-
gpt->SetPoint(gpt->GetN(), cen, val);
51-
}
52-
5365
return gpt;
54-
5566
}

‎src/.clang-format‎

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

‎src/TrackSmearer.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ TrackSmearer::loadTable(int pdg, const char *filename, bool forceReload)
1919
{
2020
auto ipdg = getIndexPDG(pdg);
2121
if (mLUTHeader[ipdg] && !forceReload) {
22-
std::cout << " --- LUT table for PDG " << pdg << " has been already loaded " << std::endl;
22+
std::cout << " --- LUT table for PDG " << pdg << " has been already loaded with index " << ipdg << std::endl;
2323
return false;
2424
}
2525
mLUTHeader[ipdg] = new lutHeader_t;

0 commit comments

Comments
 (0)