Skip to content

Commit e84384a

Browse files
mcoquet642mcoquet642
committed
Charmonium generator for pO
1 parent 58146a1 commit e84384a

File tree

4 files changed

+270
-1
lines changed

4 files changed

+270
-1
lines changed

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

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -787,6 +787,153 @@ class O2_GeneratorParamPsipp5TeV : public GeneratorTGenerator
787787
GeneratorParam* paramPsi = nullptr;
788788
};
789789

790+
class O2_GeneratorParamJpsipp96TeV : public GeneratorTGenerator
791+
{
792+
793+
public:
794+
O2_GeneratorParamJpsipp96TeV() : GeneratorTGenerator("ParamJpsi")
795+
{
796+
paramJpsi = new GeneratorParam(1, -1, PtJPsipp96TeV, YJPsipp96TeV, V2JPsipp96TeV, IpJPsipp96TeV);
797+
paramJpsi->SetMomentumRange(0., 1.e6);
798+
paramJpsi->SetPtRange(0, 999.);
799+
paramJpsi->SetYRange(-4.2, -2.3);
800+
paramJpsi->SetPhiRange(0., 360.);
801+
paramJpsi->SetDecayer(new TPythia6Decayer());
802+
paramJpsi->SetForceDecay(kNoDecay); // particle left undecayed
803+
// - - paramJpsi->SetTrackingFlag(1); // (from AliGenParam) -> check this
804+
setTGenerator(paramJpsi);
805+
};
806+
807+
~O2_GeneratorParamJpsipp96TeV()
808+
{
809+
delete paramJpsi;
810+
};
811+
812+
Bool_t Init() override
813+
{
814+
GeneratorTGenerator::Init();
815+
paramJpsi->Init();
816+
return true;
817+
}
818+
819+
void SetNSignalPerEvent(Int_t nsig) { paramJpsi->SetNumberParticles(nsig); }
820+
821+
//-------------------------------------------------------------------------//
822+
static Double_t PtJPsipp96TeV(const Double_t* px, const Double_t* /*dummy*/)
823+
{
824+
// Parameters extrapolated linearly between 5 TeV and 13 TeV as a function of log(sqrt(s))
825+
Double_t x = *px;
826+
Float_t p0, p1, p2, p3;
827+
p0 = 1;
828+
p1 = 4.61097;
829+
p2 = 1.7333;
830+
p3 = 4.45508;
831+
return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
832+
}
833+
834+
//-------------------------------------------------------------------------//
835+
static Double_t YJPsipp96TeV(const Double_t* py, const Double_t* /*dummy*/)
836+
{
837+
// Parameters extrapolated linearly between 5 TeV and 13 TeV as a function of log(sqrt(s))
838+
Double_t y = *py;
839+
Float_t p0, p1, p2;
840+
p0 = 1;
841+
p1 = 0.0107769;
842+
p2 = 2.98205;
843+
return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
844+
}
845+
846+
//-------------------------------------------------------------------------//
847+
static Double_t V2JPsipp96TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
848+
{
849+
// jpsi v2
850+
return 0.;
851+
}
852+
853+
//-------------------------------------------------------------------------//
854+
static Int_t IpJPsipp96TeV(TRandom*)
855+
{
856+
return 443;
857+
}
858+
859+
private:
860+
GeneratorParam* paramJpsi = nullptr;
861+
};
862+
863+
class O2_GeneratorParamPsipp96TeV : public GeneratorTGenerator
864+
{
865+
866+
public:
867+
O2_GeneratorParamPsipp96TeV() : GeneratorTGenerator("ParamPsi")
868+
{
869+
paramPsi = new GeneratorParam(1, -1, PtPsipp96TeV, YPsipp96TeV, V2Psipp96TeV, IpPsipp96TeV);
870+
paramPsi->SetMomentumRange(0., 1.e6);
871+
paramPsi->SetPtRange(0, 999.);
872+
paramPsi->SetYRange(-4.2, -2.3);
873+
paramPsi->SetPhiRange(0., 360.);
874+
paramPsi->SetDecayer(new TPythia6Decayer());
875+
paramPsi->SetForceDecay(kNoDecay); // particle left undecayed
876+
// - - paramJpsi->SetTrackingFlag(1); // check this
877+
setTGenerator(paramPsi);
878+
};
879+
880+
~O2_GeneratorParamPsipp96TeV()
881+
{
882+
delete paramPsi;
883+
};
884+
885+
Bool_t Init() override
886+
{
887+
GeneratorTGenerator::Init();
888+
paramPsi->Init();
889+
return true;
890+
}
891+
892+
void SetNSignalPerEvent(Int_t nsig) { paramPsi->SetNumberParticles(nsig); }
893+
894+
//-------------------------------------------------------------------------//
895+
static Double_t PtPsipp96TeV(const Double_t* px, const Double_t* /*dummy*/)
896+
{
897+
// Taking same parameters as Psi(2S) at 13 TeV
898+
Double_t x = *px;
899+
Float_t p0, p1, p2, p3;
900+
p0 = 1;
901+
p1 = 4.75208;
902+
p2 = 1.69247;
903+
p3 = 4.49224;
904+
return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
905+
}
906+
907+
//-------------------------------------------------------------------------//
908+
static Double_t YPsipp96TeV(const Double_t* py, const Double_t* /*dummy*/)
909+
{
910+
// Taking same parameters as Psi(2S) at 13 TeV
911+
Double_t y = *py;
912+
Float_t p0, p1, p2;
913+
p0 = 1;
914+
p1 = 0;
915+
p2 = 2.98887;
916+
return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
917+
}
918+
919+
//-------------------------------------------------------------------------//
920+
static Double_t V2Psipp96TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
921+
{
922+
// jpsi v2
923+
return 0.;
924+
}
925+
926+
//-------------------------------------------------------------------------//
927+
static Int_t IpPsipp96TeV(TRandom*)
928+
{
929+
return 100443;
930+
}
931+
932+
private:
933+
GeneratorParam* paramPsi = nullptr;
934+
};
935+
936+
790937
class O2_GeneratorParamJpsiPbPb5TeV : public GeneratorTGenerator
791938
{
792939

@@ -1131,6 +1278,33 @@ FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp5TeV()
11311278
return genCocktailEvtGen;
11321279
}
11331280

1281+
FairGenerator* GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp96TeV()
1282+
{
1283+
1284+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
1285+
1286+
auto genJpsi = new o2::eventgen::O2_GeneratorParamJpsipp96TeV;
1287+
genJpsi->SetNSignalPerEvent(1); // 1 J/psi generated per event by GeneratorParam
1288+
auto genPsi = new o2::eventgen::O2_GeneratorParamPsipp96TeV;
1289+
genPsi->SetNSignalPerEvent(1); // 1 Psi(2S) generated per event by GeneratorParam
1290+
genCocktailEvtGen->AddGenerator(genJpsi, 1); // add J/psi generator
1291+
genCocktailEvtGen->AddGenerator(genPsi, 1); // add Psi(2S) generator
1292+
1293+
TString pdgs = "443;100443";
1294+
std::string spdg;
1295+
TObjArray* obj = pdgs.Tokenize(";");
1296+
genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast());
1297+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
1298+
spdg = obj->At(i)->GetName();
1299+
genCocktailEvtGen->AddPdg(std::stoi(spdg), i);
1300+
printf("PDG %d \n", std::stoi(spdg));
1301+
}
1302+
genCocktailEvtGen->SetForceDecay(kEvtDiMuon);
1303+
1304+
return genCocktailEvtGen;
1305+
}
1306+
1307+
11341308
FairGenerator* GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp5TeV()
11351309
{
11361310

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

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -54,12 +54,15 @@ public:
5454
case 10: // generate cocktail of prompt X_1(3872) and psi2S to Jpsi pi pi at midrapidity
5555
mGeneratorParam = (Generator*)GeneratorCocktailX3872AndPsi2StoJpsi_pp13TeV();
5656
break;
57-
case 11: // generate cocktail of prompt X_1(3872) and psi2S to Jpsi pi pi at midrapidity
57+
case 11: // generate prompt charmonium at forward rapidity at 5TeV
5858
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp5TeV();
5959
break;
6060
case 12: // generate prompt charmonia cocktail at mid rapidity at 5TeV
6161
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp5TeV();
6262
break;
63+
case 13: // generate prompt charmonia cocktail at mid rapidity at 9.6TeV
64+
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp96TeV();
65+
break;
6366
}
6467
mGeneratorParam->Init();
6568
}

‎MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaFwdy_TriggerGap_pO96TeV.ini‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
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,13)
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_pO_961.cfg

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

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
int External()
2+
{
3+
int checkPdgSignal[] = {443,100443};
4+
int checkPdgDecay = 13;
5+
double rapiditymin = -4.3; double rapiditymax = -2.3;
6+
std::string path{"o2sim_Kine.root"};
7+
std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n";
8+
TFile file(path.c_str(), "READ");
9+
if (file.IsZombie()) {
10+
std::cerr << "Cannot open ROOT file " << path << "\n";
11+
return 1;
12+
}
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 nSignalPsi2S{};
24+
int nSignalJpsiWithinAcc{};
25+
int nSignalPsi2SWithinAcc{};
26+
auto nEvents = tree->GetEntries();
27+
o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
28+
Bool_t isInjected = kFALSE;
29+
30+
for (int i = 0; i < nEvents; i++) {
31+
tree->GetEntry(i);
32+
for (auto& track : *tracks) {
33+
auto pdg = track.GetPdgCode();
34+
auto rapidity = track.GetRapidity();
35+
auto idMoth = track.getMotherTrackId();
36+
if (pdg == checkPdgDecay) {
37+
// count leptons
38+
nLeptons++;
39+
} else if(pdg == -checkPdgDecay) {
40+
// count anti-leptons
41+
nAntileptons++;
42+
} else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1]) {
43+
if(idMoth < 0) {
44+
// count signal PDG
45+
pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
46+
// count signal PDG within acceptance
47+
if(rapidity > rapiditymin && rapidity < rapiditymax) { pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;}
48+
}
49+
auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
50+
auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
51+
if (child0 != nullptr && child1 != nullptr) {
52+
// check for parent-child relations
53+
auto pdg0 = child0->GetPdgCode();
54+
auto pdg1 = child1->GetPdgCode();
55+
std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
56+
if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
57+
nLeptonPairs++;
58+
if (child0->getToBeDone() && child1->getToBeDone()) {
59+
nLeptonPairsToBeDone++;
60+
}
61+
}
62+
}
63+
}
64+
}
65+
}
66+
std::cout << "#events: " << nEvents << "\n"
67+
<< "#leptons: " << nLeptons << "\n"
68+
<< "#antileptons: " << nAntileptons << "\n"
69+
<< "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
70+
<< "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
71+
<< "#lepton pairs: " << nLeptonPairs << "\n"
72+
<< "#lepton pairs to be done: " << nLeptonPairs << "\n";
73+
74+
75+
if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) {
76+
std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
77+
return 1;
78+
}
79+
if (nLeptonPairs != nLeptonPairsToBeDone) {
80+
std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
81+
return 1;
82+
}
83+
84+
return 0;
85+
}

0 commit comments

Comments
 (0)