Skip to content

Commit d1293a0

Browse files
ypwanggypwangg
committed
add test file and update centrality cut
1 parent c97b94a commit d1293a0

5 files changed

Lines changed: 422 additions & 0 deletions

‎MC/config/PWGDQ/external/generator/generator_pythia8_StarlightTriggered_withGap.C‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -146,6 +146,12 @@ protected:
146146

147147
int getNumSignalsForImpactParameter(double impactParameter) const
148148
{
149+
const double photoproductionImpactParameterCut = 4.0; // tmp cut to embed signal > 10% centrality
150+
if (impactParameter < photoproductionImpactParameterCut) {
151+
std::cout<<"impact parameter: "<<impactParameter<<" fm, below photoproduction cut "
152+
<<photoproductionImpactParameterCut<<" fm; embedding no signal events"<<std::endl;
153+
return 0;
154+
}
149155
return std::max(0, static_cast<int>(std::lround(5.0 + 0.886202881 * std::pow(std::max(0.0, 17.5 - impactParameter), 1.7))));
150156
}
151157

‎MC/config/PWGDQ/ini/tests/Generator_InjectedCoherentCharmoniaFwdy_PbPb5TeV.C‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
TFile file(path.c_str(), "READ");
4+
if (file.IsZombie()) {
5+
std::cerr << "Cannot open ROOT file " << path << "\n";
6+
return 1;
7+
}
8+
auto tree = (TTree *)file.Get("o2sim");
9+
if (!tree)
10+
{
11+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
12+
return 1;
13+
}
14+
std::vector<o2::MCTrack>* tracks{};
15+
tree->SetBranchAddress("MCTrack", &tracks);
16+
17+
int checkPdgSignal[2] = {443, 100443};
18+
int checkPdgDecay = 11;
19+
double rapiditymin = 2.0; double rapiditymax = 4.5;
20+
21+
int nLeptons{};
22+
int nAntileptons{};
23+
int nLeptonPairs{};
24+
int nPrimaryLepton{};
25+
int nPrimaryAntilepton{};
26+
int nLeptonPairsToBeDone{};
27+
int nSignalJpsi{};
28+
int nSignalPsi2S{};
29+
int nPions{};
30+
int nPhotonsFromPsi{};
31+
int nPhotonsFromJPsi{};
32+
int nSignalJpsiWithinAcc{};
33+
int nSignalPsi2SWithinAcc{};
34+
35+
auto nEvents = tree->GetEntries();
36+
for (int i = 0; i < nEvents; i++) {
37+
tree->GetEntry(i);
38+
for (auto& track : *tracks) {
39+
auto pdg = track.GetPdgCode();
40+
auto rapidity = track.GetRapidity();
41+
auto idMoth = track.getMotherTrackId();
42+
if (pdg == checkPdgDecay) {
43+
// count leptons
44+
nLeptons++;
45+
if (idMoth == -1) nPrimaryLepton++;
46+
} else if (pdg == -checkPdgDecay) {
47+
// count anti-leptons
48+
nAntileptons++;
49+
if (idMoth == -1) nPrimaryAntilepton++;
50+
} else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) {
51+
if (idMoth < 0) {
52+
// count signal PDG
53+
pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
54+
// count signal PDG within acceptance
55+
if (rapidity > rapiditymin && rapidity < rapiditymax) {
56+
pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;
57+
}
58+
}
59+
int child0Id = track.getFirstDaughterTrackId();
60+
int child1Id = track.getLastDaughterTrackId();
61+
for (int iDaugh = child0Id; iDaugh <= child1Id; iDaugh++) {
62+
auto daugh = tracks->at(iDaugh);
63+
if (daugh.GetPdgCode() == 22) {
64+
if (pdg == checkPdgSignal[0]) {
65+
nPhotonsFromJPsi++;
66+
} else if (pdg == checkPdgSignal[1]) {
67+
nPhotonsFromPsi++;
68+
}
69+
}
70+
if (TMath::Abs(daugh.GetPdgCode()) == 211) nPions++;
71+
}
72+
auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
73+
auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
74+
if (child0 != nullptr && child1 != nullptr) {
75+
// check for parent-child relations
76+
auto pdg0 = child0->GetPdgCode();
77+
auto pdg1 = child1->GetPdgCode();
78+
// std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
79+
if (pdg == checkPdgSignal[0]) std::cout << "First and last children of parent Jpsi are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
80+
if (pdg == checkPdgSignal[1]) std::cout << "First and last children of parent Psi(2S) are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
81+
if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
82+
nLeptonPairs++;
83+
if (child0->getToBeDone() && child1->getToBeDone()) {
84+
nLeptonPairsToBeDone++;
85+
}
86+
}
87+
}
88+
}
89+
}
90+
}
91+
92+
std::cout << "#events: " << nEvents << "\n"
93+
<< "#leptons: " << nLeptons << "\n"
94+
<< "#antileptons: " << nAntileptons << "\n"
95+
<< "#primary leptons: " << nPrimaryLepton << "\n"
96+
<< "#primary antileptons: " << nPrimaryAntilepton << "\n"
97+
<< "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
98+
<< "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
99+
<< "#photons from Jpsi: " << nPhotonsFromJPsi <<" #photons from Psi(2S): "<< nPhotonsFromPsi<<"\n"
100+
<< "#pions from Psi(2S): "<< nPions<<"\n"
101+
<< "#lepton pairs: " << nLeptonPairs <<" #lepton pairs to be done: "<< nLeptonPairsToBeDone<<"\n";
102+
103+
return 0;
104+
}

‎MC/config/PWGDQ/ini/tests/Generator_InjectedCoherentCharmoniaMidy_PbPb5TeV.C‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
TFile file(path.c_str(), "READ");
4+
if (file.IsZombie()) {
5+
std::cerr << "Cannot open ROOT file " << path << "\n";
6+
return 1;
7+
}
8+
auto tree = (TTree *)file.Get("o2sim");
9+
if (!tree)
10+
{
11+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
12+
return 1;
13+
}
14+
std::vector<o2::MCTrack>* tracks{};
15+
tree->SetBranchAddress("MCTrack", &tracks);
16+
17+
int checkPdgSignal[2] = {443, 100443};
18+
int checkPdgDecay = 11;
19+
double rapiditymin = -1.5; double rapiditymax = 1.5;
20+
21+
int nLeptons{};
22+
int nAntileptons{};
23+
int nLeptonPairs{};
24+
int nPrimaryLepton{};
25+
int nPrimaryAntilepton{};
26+
int nLeptonPairsToBeDone{};
27+
int nSignalJpsi{};
28+
int nSignalPsi2S{};
29+
int nPions{};
30+
int nPhotonsFromPsi{};
31+
int nPhotonsFromJPsi{};
32+
int nSignalJpsiWithinAcc{};
33+
int nSignalPsi2SWithinAcc{};
34+
35+
auto nEvents = tree->GetEntries();
36+
for (int i = 0; i < nEvents; i++) {
37+
tree->GetEntry(i);
38+
for (auto& track : *tracks) {
39+
auto pdg = track.GetPdgCode();
40+
auto rapidity = track.GetRapidity();
41+
auto idMoth = track.getMotherTrackId();
42+
if (pdg == checkPdgDecay) {
43+
// count leptons
44+
nLeptons++;
45+
if (idMoth == -1) nPrimaryLepton++;
46+
} else if (pdg == -checkPdgDecay) {
47+
// count anti-leptons
48+
nAntileptons++;
49+
if (idMoth == -1) nPrimaryAntilepton++;
50+
} else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) {
51+
if (idMoth < 0) {
52+
// count signal PDG
53+
pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
54+
// count signal PDG within acceptance
55+
if (rapidity > rapiditymin && rapidity < rapiditymax) {
56+
pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;
57+
}
58+
}
59+
int child0Id = track.getFirstDaughterTrackId();
60+
int child1Id = track.getLastDaughterTrackId();
61+
for (int iDaugh = child0Id; iDaugh <= child1Id; iDaugh++) {
62+
auto daugh = tracks->at(iDaugh);
63+
if (daugh.GetPdgCode() == 22) {
64+
if (pdg == checkPdgSignal[0]) {
65+
nPhotonsFromJPsi++;
66+
} else if (pdg == checkPdgSignal[1]) {
67+
nPhotonsFromPsi++;
68+
}
69+
}
70+
if (TMath::Abs(daugh.GetPdgCode()) == 211) nPions++;
71+
}
72+
auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
73+
auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
74+
if (child0 != nullptr && child1 != nullptr) {
75+
// check for parent-child relations
76+
auto pdg0 = child0->GetPdgCode();
77+
auto pdg1 = child1->GetPdgCode();
78+
// std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
79+
if (pdg == checkPdgSignal[0]) std::cout << "First and last children of parent Jpsi are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
80+
if (pdg == checkPdgSignal[1]) std::cout << "First and last children of parent Psi(2S) are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
81+
if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
82+
nLeptonPairs++;
83+
if (child0->getToBeDone() && child1->getToBeDone()) {
84+
nLeptonPairsToBeDone++;
85+
}
86+
}
87+
}
88+
}
89+
}
90+
}
91+
92+
std::cout << "#events: " << nEvents << "\n"
93+
<< "#leptons: " << nLeptons << "\n"
94+
<< "#antileptons: " << nAntileptons << "\n"
95+
<< "#primary leptons: " << nPrimaryLepton << "\n"
96+
<< "#primary antileptons: " << nPrimaryAntilepton << "\n"
97+
<< "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
98+
<< "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
99+
<< "#photons from Jpsi: " << nPhotonsFromJPsi <<" #photons from Psi(2S): "<< nPhotonsFromPsi<<"\n"
100+
<< "#pions from Psi(2S): "<< nPions<<"\n"
101+
<< "#lepton pairs: " << nLeptonPairs <<" #lepton pairs to be done: "<< nLeptonPairsToBeDone<<"\n";
102+
103+
return 0;
104+
}

‎MC/config/PWGDQ/ini/tests/Generator_InjectedIncoherentCharmoniaFwdy_PbPb5TeV.C‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
TFile file(path.c_str(), "READ");
4+
if (file.IsZombie()) {
5+
std::cerr << "Cannot open ROOT file " << path << "\n";
6+
return 1;
7+
}
8+
auto tree = (TTree *)file.Get("o2sim");
9+
if (!tree)
10+
{
11+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
12+
return 1;
13+
}
14+
std::vector<o2::MCTrack>* tracks{};
15+
tree->SetBranchAddress("MCTrack", &tracks);
16+
17+
int checkPdgSignal[2] = {443, 100443};
18+
int checkPdgDecay = 11;
19+
double rapiditymin = 2.0; double rapiditymax = 4.5;
20+
21+
int nLeptons{};
22+
int nAntileptons{};
23+
int nLeptonPairs{};
24+
int nPrimaryLepton{};
25+
int nPrimaryAntilepton{};
26+
int nLeptonPairsToBeDone{};
27+
int nSignalJpsi{};
28+
int nSignalPsi2S{};
29+
int nPions{};
30+
int nPhotonsFromPsi{};
31+
int nPhotonsFromJPsi{};
32+
int nSignalJpsiWithinAcc{};
33+
int nSignalPsi2SWithinAcc{};
34+
35+
auto nEvents = tree->GetEntries();
36+
for (int i = 0; i < nEvents; i++) {
37+
tree->GetEntry(i);
38+
for (auto& track : *tracks) {
39+
auto pdg = track.GetPdgCode();
40+
auto rapidity = track.GetRapidity();
41+
auto idMoth = track.getMotherTrackId();
42+
if (pdg == checkPdgDecay) {
43+
// count leptons
44+
nLeptons++;
45+
if (idMoth == -1) nPrimaryLepton++;
46+
} else if (pdg == -checkPdgDecay) {
47+
// count anti-leptons
48+
nAntileptons++;
49+
if (idMoth == -1) nPrimaryAntilepton++;
50+
} else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) {
51+
if (idMoth < 0) {
52+
// count signal PDG
53+
pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
54+
// count signal PDG within acceptance
55+
if (rapidity > rapiditymin && rapidity < rapiditymax) {
56+
pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;
57+
}
58+
}
59+
int child0Id = track.getFirstDaughterTrackId();
60+
int child1Id = track.getLastDaughterTrackId();
61+
for (int iDaugh = child0Id; iDaugh <= child1Id; iDaugh++) {
62+
auto daugh = tracks->at(iDaugh);
63+
if (daugh.GetPdgCode() == 22) {
64+
if (pdg == checkPdgSignal[0]) {
65+
nPhotonsFromJPsi++;
66+
} else if (pdg == checkPdgSignal[1]) {
67+
nPhotonsFromPsi++;
68+
}
69+
}
70+
if (TMath::Abs(daugh.GetPdgCode()) == 211) nPions++;
71+
}
72+
auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
73+
auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
74+
if (child0 != nullptr && child1 != nullptr) {
75+
// check for parent-child relations
76+
auto pdg0 = child0->GetPdgCode();
77+
auto pdg1 = child1->GetPdgCode();
78+
// std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
79+
if (pdg == checkPdgSignal[0]) std::cout << "First and last children of parent Jpsi are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
80+
if (pdg == checkPdgSignal[1]) std::cout << "First and last children of parent Psi(2S) are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
81+
if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
82+
nLeptonPairs++;
83+
if (child0->getToBeDone() && child1->getToBeDone()) {
84+
nLeptonPairsToBeDone++;
85+
}
86+
}
87+
}
88+
}
89+
}
90+
}
91+
92+
std::cout << "#events: " << nEvents << "\n"
93+
<< "#leptons: " << nLeptons << "\n"
94+
<< "#antileptons: " << nAntileptons << "\n"
95+
<< "#primary leptons: " << nPrimaryLepton << "\n"
96+
<< "#primary antileptons: " << nPrimaryAntilepton << "\n"
97+
<< "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
98+
<< "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
99+
<< "#photons from Jpsi: " << nPhotonsFromJPsi <<" #photons from Psi(2S): "<< nPhotonsFromPsi<<"\n"
100+
<< "#pions from Psi(2S): "<< nPions<<"\n"
101+
<< "#lepton pairs: " << nLeptonPairs <<" #lepton pairs to be done: "<< nLeptonPairsToBeDone<<"\n";
102+
103+
return 0;
104+
}

0 commit comments

Comments
 (0)