Skip to content

Commit 262c869

Browse files
ffiondaffionda
authored
add X3872 and Psi2S to Jpsi pi pi (#1801)
* add X3872 and Psi2S to Jpsi pi pi
1 parent 82c2c51 commit 262c869

8 files changed

+351
-1
lines changed

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

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
Decay J/psi
2+
### from DECAY.DEC
3+
1.000 e+ e- PHOTOS VLL;
4+
Enddecay
5+
Decay psi(2S)
6+
1.000 J/psi pi+ pi- VVPIPI;
7+
Enddecay
8+
End

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

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
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+
Decay J/psi
6+
### from DECAY.DEC
7+
1.000 e+ e- PHOTOS VLL;
8+
Enddecay
9+
10+
End

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

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -638,6 +638,75 @@ class O2_GeneratorParamPsiPbPb5TeV : public GeneratorTGenerator
638638
GeneratorParam* paramPsi = nullptr;
639639
};
640640

641+
642+
class O2_GeneratorParamX3872MidY : public GeneratorTGenerator
643+
{
644+
645+
public:
646+
O2_GeneratorParamX3872MidY() : GeneratorTGenerator("ParamX3872MidY")
647+
{
648+
paramX3872 = new GeneratorParam(1, -1, PtX3872pp13TeV, YX3872pp13TeV, V2X3872pp13TeV, IpX3872pp13TeV);
649+
paramX3872->SetMomentumRange(0., 1.e6);
650+
paramX3872->SetPtRange(0., 1000.);
651+
paramX3872->SetYRange(-1.0, 1.0);
652+
paramX3872->SetPhiRange(0., 360.);
653+
paramX3872->SetDecayer(new TPythia6Decayer()); // Pythia
654+
paramX3872->SetForceDecay(kNoDecay); // particle left undecayed
655+
setTGenerator(paramX3872);
656+
};
657+
658+
~O2_GeneratorParamX3872MidY()
659+
{
660+
delete paramX3872;
661+
};
662+
663+
Bool_t Init() override
664+
{
665+
GeneratorTGenerator::Init();
666+
paramX3872->Init();
667+
return true;
668+
}
669+
670+
void SetNSignalPerEvent(Int_t nsig) { paramX3872->SetNumberParticles(nsig); }
671+
672+
//-------------------------------------------------------------------------//
673+
static Double_t PtX3872pp13TeV(const Double_t* px, const Double_t* /*dummy*/)
674+
{
675+
// prompt X3872 pT
676+
// pp, 13TeV (tuned LHCb pp 13 TeV)
677+
//
678+
const Double_t kC = 7.64519e+00 ;
679+
const Double_t kpt0 = 5.30628e+00;
680+
const Double_t kn = 3.30887e+00;
681+
Double_t pt = px[0];
682+
return kC * pt / TMath::Power((1. + (pt / kpt0) * (pt / kpt0)), kn);
683+
}
684+
685+
//-------------------------------------------------------------------------//
686+
static Double_t YX3872pp13TeV(const Double_t* py, const Double_t* /*dummy*/)
687+
{
688+
// flat rapidity distribution assumed at midrapidity
689+
return 1.;
690+
}
691+
692+
//-------------------------------------------------------------------------//
693+
static Double_t V2X3872pp13TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
694+
{
695+
// X3872 v2
696+
return 0.;
697+
}
698+
699+
//-------------------------------------------------------------------------//
700+
static Int_t IpX3872pp13TeV(TRandom*)
701+
{
702+
return 9920443;
703+
}
704+
705+
private:
706+
GeneratorParam* paramX3872 = nullptr;
707+
};
708+
709+
641710
} // namespace eventgen
642711
} // namespace o2
643712

@@ -741,6 +810,31 @@ FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV()
741810
return genCocktailEvtGen;
742811
}
743812

813+
814+
FairGenerator*
815+
GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV(TString pdgs = "100443")
816+
{
817+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamPsiMidY>();
818+
gen->SetNSignalPerEvent(1); // number of jpsis per event
819+
820+
std::string spdg;
821+
TObjArray* obj = pdgs.Tokenize(";");
822+
gen->SetSizePdg(obj->GetEntriesFast());
823+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
824+
spdg = obj->At(i)->GetName();
825+
gen->AddPdg(std::stoi(spdg), i);
826+
printf("PDG %d \n", std::stoi(spdg));
827+
}
828+
TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
829+
gen->SetDecayTable(Form("%s/PSITOJPSIPIPI.DEC", pathO2.Data()));
830+
831+
// print debug
832+
gen->PrintDebug();
833+
834+
return gen;
835+
}
836+
837+
744838
FairGenerator*
745839
GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV(TString pdgs = "443")
746840
{
@@ -844,6 +938,29 @@ FairGenerator*
844938
}
845939

846940

941+
FairGenerator*
942+
GeneratorParamX3872ToJpsiEvtGen_pp13TeV(TString pdgs = "9920443")
943+
{
944+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamX3872MidY>();
945+
gen->SetNSignalPerEvent(1); // number of jpsis per event
946+
947+
std::string spdg;
948+
TObjArray* obj = pdgs.Tokenize(";");
949+
gen->SetSizePdg(obj->GetEntriesFast());
950+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
951+
spdg = obj->At(i)->GetName();
952+
gen->AddPdg(std::stoi(spdg), i);
953+
printf("PDG %d \n", std::stoi(spdg));
954+
}
955+
TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
956+
gen->SetDecayTable(Form("%s/X3872TOJPSIPIPI.DEC", pathO2.Data()));
957+
958+
// print debug
959+
gen->PrintDebug();
960+
961+
return gen;
962+
}
963+
847964

848965

849966

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

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,12 @@ public:
4545
case 7: // generate prompt charmonia cocktail at forward rapidity at 5TeV
4646
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV();
4747
break;
48-
48+
case 8: // generate prompt X_1(3872) to Jpsi pi pi at midrapidity
49+
mGeneratorParam = (Generator*)GeneratorParamX3872ToJpsiEvtGen_pp13TeV("9920443");
50+
break;
51+
case 9: // generate prompt psi2S to Jpsi pi pi at midrapidity
52+
mGeneratorParam = (Generator*)GeneratorParamPromptPsiToJpsiPiPiEvtGen_pp13TeV("100443");
53+
break;
4954
}
5055
mGeneratorParam->Init();
5156
}

‎MC/config/PWGDQ/ini/Generator_InjectedPromptPsi2SToJpsiPiPiMidy_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,9)
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg

‎MC/config/PWGDQ/ini/Generator_InjectedPromptX3782Midy_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,8)
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg

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

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

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

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

0 commit comments

Comments
 (0)