Skip to content

Commit 5037c1c

Browse files
ffiondaalcaliva
ffionda
authored and
alcaliva
committed
update for generating a cocktail of X3872 + Psi2S to Jpsi pi pi (#1860)
(cherry picked from commit 198ebde)
1 parent d93a4c7 commit 5037c1c

File tree

5 files changed

+153
-0
lines changed

5 files changed

+153
-0
lines changed

‎MC/config/PWGDQ/EvtGen/DecayTablesEvtgen/X3872ANDPSI2STOJPSIPIPI.DEC‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
Decay X_1(3872)
2+
#### Breit-Wigner function for the pion distribution (S-Wave approximation)
3+
1.000 J/psi pi+ pi- XJPSIRHO0PIPI_S;
4+
Enddecay
5+
6+
Decay psi(2S)
7+
1.000 J/psi pi+ pi- VVPIPI;
8+
Enddecay
9+
10+
Decay J/psi
11+
### from DECAY.DEC
12+
1.000 e+ e- PHOTOS VLL;
13+
Enddecay
14+
15+
End

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

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -962,5 +962,30 @@ FairGenerator*
962962
}
963963

964964

965+
FairGenerator* GeneratorCocktailX3872AndPsi2StoJpsi_pp13TeV()
966+
{
967+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
968+
969+
auto genX3872 = new o2::eventgen::O2_GeneratorParamX3872MidY;
970+
genX3872->SetNSignalPerEvent(1); // number of jpsis per event
971+
auto genPsi2S = new o2::eventgen::O2_GeneratorParamPsiMidY;
972+
genPsi2S->SetNSignalPerEvent(1); // number of jpsis per event
973+
genCocktailEvtGen->AddGenerator(genX3872, 1); // add J/psi generator
974+
genCocktailEvtGen->AddGenerator(genPsi2S, 1); // add Psi(2S) generator
975+
976+
TString pdgs = "9920443;100443";
977+
std::string spdg;
978+
TObjArray* obj = pdgs.Tokenize(";");
979+
genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast());
980+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
981+
spdg = obj->At(i)->GetName();
982+
genCocktailEvtGen->AddPdg(std::stoi(spdg), i);
983+
printf("PDG %d \n", std::stoi(spdg));
984+
}
985+
TString pathO2 = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
986+
genCocktailEvtGen->SetDecayTable(Form("%s/X3872ANDPSI2STOJPSIPIPI.DEC", pathO2.Data()));
987+
988+
return genCocktailEvtGen;
989+
}
965990

966991

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,9 @@ public:
5151
case 9: // generate prompt psi2S to Jpsi pi pi at midrapidity
5252
mGeneratorParam = (Generator*)GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV("100443");
5353
break;
54+
case 10: // generate cocktail of prompt X_1(3872) and psi2S to Jpsi pi pi at midrapidity
55+
mGeneratorParam = (Generator*)GeneratorCocktailX3872AndPsi2StoJpsi_pp13TeV();
56+
break;
5457
}
5558
mGeneratorParam->Init();
5659
}

‎MC/config/PWGDQ/ini/Generator_InjectedPromptX3782andPsi2SMidy_TriggerGap.ini‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
### The external generator derives from GeneratorPythia8.
2+
[GeneratorExternal]
3+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C
4+
funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,10)
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg

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

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
int External()
2+
{
3+
int checkPdgSignal[] = {9920443,100443}; // pdg code X3872
4+
TString PdgSignalName[] = {"X(3872)", "Psi2S"};
5+
int checkPdgDecay[] = {443, 211, -211};
6+
int leptonPdg = 11;
7+
Double_t rapidityWindow = 1.0;
8+
std::string path{"o2sim_Kine.root"};
9+
for(int iSig =0; iSig < 2; iSig++) std::cout << "Check for\nsignal PDG " << checkPdgSignal[iSig] << "\n decay PDG " << checkPdgDecay[0] << ", " << checkPdgDecay[1] << ", " << checkPdgDecay[2] << "\n";
10+
TFile file(path.c_str(), "READ");
11+
if (file.IsZombie()) {
12+
std::cerr << "Cannot open ROOT file " << path << "\n";
13+
return 1;
14+
}
15+
auto tree = (TTree*)file.Get("o2sim");
16+
std::vector<o2::MCTrack>* tracks{};
17+
tree->SetBranchAddress("MCTrack", &tracks);
18+
19+
int nLeptons[]={0,0};
20+
int nAntileptons[]={0,0};
21+
int nLeptonPairs[]={0,0};
22+
int nLeptonPairsToBeDone[]={0,0};
23+
int nSignalX3872[]={0,0};
24+
int nSignalPionsPos[]={0,0};
25+
int nSignalPionsNeg[]={0,0};
26+
int nSignalPsi2S{};
27+
int nSignalX3872WithinAcc[]={0,0};
28+
int nSignalPionsPosWithinAcc[]={0,0};
29+
int nSignalPionsNegWithinAcc[]={0,0};
30+
auto nEvents = tree->GetEntries();
31+
o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
32+
Bool_t hasPsi2SMoth = kFALSE;
33+
34+
for (int i = 0; i < nEvents; i++) {
35+
tree->GetEntry(i);
36+
for (auto& track : *tracks) {
37+
auto pdg = track.GetPdgCode();
38+
auto rapidity = track.GetRapidity();
39+
auto idMoth = track.getMotherTrackId();
40+
int idX3872 = -1; int IdChild0 = -1; int IdChild1 = -1;
41+
for(int iSig=0; iSig<2; iSig++) {
42+
if (pdg == leptonPdg) {
43+
// count leptons
44+
nLeptons[iSig]++;
45+
} else if(pdg == -leptonPdg) {
46+
// count anti-leptons
47+
nAntileptons[iSig]++;
48+
} else if (pdg == checkPdgSignal[iSig]) {
49+
// check daughters
50+
std::cout << "Signal PDG: " << pdg << "\n";
51+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
52+
auto pdgDau = tracks->at(j).GetPdgCode();
53+
std::cout << "Daughter " << j << " is: " << pdgDau << "\n";
54+
if(TMath::Abs(pdgDau) == checkPdgDecay[0] ) { nSignalX3872[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalX3872WithinAcc[iSig]++; idX3872 = j; }
55+
if(pdgDau == checkPdgDecay[1] ) { nSignalPionsPos[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsPosWithinAcc[iSig]++; }
56+
if(pdgDau == checkPdgDecay[2] ) { nSignalPionsNeg[iSig]++; if( std::abs(track.GetRapidity()) < rapidityWindow) nSignalPionsNegWithinAcc[iSig]++; }
57+
}
58+
59+
auto trackX3872 = tracks->at(idX3872);
60+
for (int j{trackX3872.getFirstDaughterTrackId()}; j <= trackX3872.getLastDaughterTrackId(); ++j) {
61+
auto pdgDau = tracks->at(j).GetPdgCode();
62+
if(pdgDau == leptonPdg ) IdChild0 = j;
63+
if(pdgDau == -leptonPdg ) IdChild1 = j;
64+
}
65+
auto child0 = tracks->at(IdChild0);
66+
auto child1 = tracks->at(IdChild1);
67+
// check for parent-child relations
68+
auto pdg0 = child0.GetPdgCode();
69+
auto pdg1 = child1.GetPdgCode();
70+
std::cout << "Lepton daughter particles of mother " << trackX3872.GetPdgCode() << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
71+
if (std::abs(pdg0) == leptonPdg && std::abs(pdg1) == leptonPdg && pdg0 == -pdg1) {
72+
nLeptonPairs[iSig]++;
73+
if (child0.getToBeDone() && child1.getToBeDone()) {
74+
nLeptonPairsToBeDone[iSig]++;
75+
}
76+
}
77+
}
78+
}
79+
}
80+
}
81+
82+
std::cout << "#events: " << nEvents << "\n";
83+
for(int iSig=0; iSig < 2; iSig++){
84+
std::cout << "#leptons from " << PdgSignalName[iSig] << ": " << nLeptons[iSig] << "\n"
85+
<< "#antileptons from " << PdgSignalName[iSig] << ": " << nAntileptons[iSig] << "\n"
86+
<< "#signal (jpsi <-" << PdgSignalName[iSig] <<"): " << nSignalX3872[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalX3872WithinAcc[iSig] << "\n"
87+
<< "#signal (pi+ <-" << PdgSignalName[iSig] <<"): " << nSignalPionsPos[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsPosWithinAcc[iSig] << "\n"
88+
<< "#signal (pi- <-" << PdgSignalName[iSig] <<"): " << nSignalPionsNeg[iSig] << "; within acceptance (|y| < " << rapidityWindow << "): " << nSignalPionsNegWithinAcc[iSig] << "\n"
89+
<< "#lepton pairs from " << PdgSignalName[iSig] <<": " << nLeptonPairs[iSig] << "\n"
90+
<< "#lepton pairs to be done from " << PdgSignalName[iSig] <<": " << nLeptonPairs[iSig] << "\n";
91+
92+
93+
if (nLeptonPairs[iSig] == 0 || nLeptons[iSig] == 0 || nAntileptons[iSig] == 0) {
94+
std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
95+
return 1;
96+
}
97+
if (nLeptonPairs[iSig] != nLeptonPairsToBeDone[iSig]) {
98+
std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
99+
return 1;
100+
}
101+
102+
}
103+
return 0;
104+
}

0 commit comments

Comments
 (0)