Skip to content

Commit 56d5c8d

Browse files
lucamicheletti93lucamicheletti
authored and
alcaliva
committed
[PWGDQ] Adding scripts to run DQ simulations in Pb-Pb collisions (#1705)
* Adding script to run Pb-Pb sim with charmonia injected * Adding test file --------- Co-authored-by: Lucamicheletti93 <luca.mike93@gmail.com> (cherry picked from commit d13558d)
1 parent 825e9e2 commit 56d5c8d

File tree

5 files changed

+278
-0
lines changed

5 files changed

+278
-0
lines changed

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

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -492,6 +492,152 @@ class O2_GeneratorParamChiC2 : public GeneratorTGenerator
492492
GeneratorParam* paramChiC2 = nullptr;
493493
};
494494

495+
class O2_GeneratorParamJpsiPbPb5TeV : public GeneratorTGenerator
496+
{
497+
498+
public:
499+
O2_GeneratorParamJpsiPbPb5TeV() : GeneratorTGenerator("ParamJpsi")
500+
{
501+
paramJpsi = new GeneratorParam(1, -1, PtJPsiPbPb5TeV, YJPsiPbPb5TeV, V2JPsiPbPb5TeV, IpJPsiPbPb5TeV);
502+
paramJpsi->SetMomentumRange(0., 1.e6);
503+
paramJpsi->SetPtRange(0, 999.);
504+
paramJpsi->SetYRange(-4.2, -2.3);
505+
paramJpsi->SetPhiRange(0., 360.);
506+
paramJpsi->SetDecayer(new TPythia6Decayer());
507+
paramJpsi->SetForceDecay(kNoDecay); // particle left undecayed
508+
// - - paramJpsi->SetTrackingFlag(1); // (from AliGenParam) -> check this
509+
setTGenerator(paramJpsi);
510+
};
511+
512+
~O2_GeneratorParamJpsiPbPb5TeV()
513+
{
514+
delete paramJpsi;
515+
};
516+
517+
Bool_t Init() override
518+
{
519+
GeneratorTGenerator::Init();
520+
paramJpsi->Init();
521+
return true;
522+
}
523+
524+
void SetNSignalPerEvent(Int_t nsig) { paramJpsi->SetNumberParticles(nsig); }
525+
526+
//-------------------------------------------------------------------------//
527+
static Double_t PtJPsiPbPb5TeV(const Double_t* px, const Double_t* /*dummy*/)
528+
{
529+
// jpsi pT in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
530+
Double_t x = *px;
531+
Float_t p0, p1, p2, p3;
532+
p0 = 1.00715e6;
533+
p1 = 3.50274;
534+
p2 = 1.93403;
535+
p3 = 3.96363;
536+
return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
537+
}
538+
539+
//-------------------------------------------------------------------------//
540+
static Double_t YJPsiPbPb5TeV(const Double_t* py, const Double_t* /*dummy*/)
541+
{
542+
// jpsi y in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
543+
Double_t y = *py;
544+
Float_t p0, p1, p2;
545+
p0 = 1.09886e6;
546+
p1 = 0;
547+
p2 = 2.12568;
548+
return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
549+
}
550+
551+
//-------------------------------------------------------------------------//
552+
static Double_t V2JPsiPbPb5TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
553+
{
554+
// jpsi v2
555+
return 0.;
556+
}
557+
558+
//-------------------------------------------------------------------------//
559+
static Int_t IpJPsiPbPb5TeV(TRandom*)
560+
{
561+
return 443;
562+
}
563+
564+
private:
565+
GeneratorParam* paramJpsi = nullptr;
566+
};
567+
568+
class O2_GeneratorParamPsiPbPb5TeV : public GeneratorTGenerator
569+
{
570+
571+
public:
572+
O2_GeneratorParamPsiPbPb5TeV() : GeneratorTGenerator("ParamPsi")
573+
{
574+
paramPsi = new GeneratorParam(1, -1, PtPsiPbPb5TeV, YPsiPbPb5TeV, V2PsiPbPb5TeV, IpPsiPbPb5TeV);
575+
paramPsi->SetMomentumRange(0., 1.e6);
576+
paramPsi->SetPtRange(0, 999.);
577+
paramPsi->SetYRange(-4.2, -2.3);
578+
paramPsi->SetPhiRange(0., 360.);
579+
paramPsi->SetDecayer(new TPythia6Decayer());
580+
paramPsi->SetForceDecay(kNoDecay); // particle left undecayed
581+
// - - paramJpsi->SetTrackingFlag(1); // check this
582+
setTGenerator(paramPsi);
583+
};
584+
585+
~O2_GeneratorParamPsiPbPb5TeV()
586+
{
587+
delete paramPsi;
588+
};
589+
590+
Bool_t Init() override
591+
{
592+
GeneratorTGenerator::Init();
593+
paramPsi->Init();
594+
return true;
595+
}
596+
597+
void SetNSignalPerEvent(Int_t nsig) { paramPsi->SetNumberParticles(nsig); }
598+
599+
//-------------------------------------------------------------------------//
600+
static Double_t PtPsiPbPb5TeV(const Double_t* px, const Double_t* /*dummy*/)
601+
{
602+
// jpsi pT in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
603+
Double_t x = *px;
604+
Float_t p0, p1, p2, p3;
605+
p0 = 1.00715e6;
606+
p1 = 3.50274;
607+
p2 = 1.93403;
608+
p3 = 3.96363;
609+
return p0 * x / TMath::Power(1. + TMath::Power(x / p1, p2), p3);
610+
}
611+
612+
//-------------------------------------------------------------------------//
613+
static Double_t YPsiPbPb5TeV(const Double_t* py, const Double_t* /*dummy*/)
614+
{
615+
// jpsi y in PbPb, tuned on data (2015) -> Castillo embedding https://alice.its.cern.ch/jira/browse/ALIROOT-8174?jql=text%20~%20%22LHC19a2%22
616+
Double_t y = *py;
617+
Float_t p0, p1, p2;
618+
p0 = 1.09886e6;
619+
p1 = 0;
620+
p2 = 2.12568;
621+
return p0 * TMath::Exp(-(1. / 2.) * TMath::Power(((y - p1) / p2), 2));
622+
}
623+
624+
//-------------------------------------------------------------------------//
625+
static Double_t V2PsiPbPb5TeV(const Double_t* /*dummy*/, const Double_t* /*dummy*/)
626+
{
627+
// jpsi v2
628+
return 0.;
629+
}
630+
631+
//-------------------------------------------------------------------------//
632+
static Int_t IpPsiPbPb5TeV(TRandom*)
633+
{
634+
return 100443;
635+
}
636+
637+
private:
638+
GeneratorParam* paramPsi = nullptr;
639+
};
640+
495641
} // namespace eventgen
496642
} // namespace o2
497643

@@ -671,6 +817,32 @@ FairGenerator*
671817
return genCocktailEvtGen;
672818
}
673819

820+
FairGenerator*
821+
GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV()
822+
{
823+
auto genCocktailEvtGen = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
824+
825+
auto genJpsi = new o2::eventgen::O2_GeneratorParamJpsiPbPb5TeV;
826+
genJpsi->SetNSignalPerEvent(4); // 4 J/psi generated per event by GeneratorParam
827+
auto genPsi = new o2::eventgen::O2_GeneratorParamPsiPbPb5TeV;
828+
genPsi->SetNSignalPerEvent(2); // 2 Psi(2S) generated per event by GeneratorParam
829+
genCocktailEvtGen->AddGenerator(genJpsi, 1); // 2/3 J/psi
830+
genCocktailEvtGen->AddGenerator(genPsi, 1); // 1/3 Psi(2S)
831+
832+
TString pdgs = "443;100443";
833+
std::string spdg;
834+
TObjArray* obj = pdgs.Tokenize(";");
835+
genCocktailEvtGen->SetSizePdg(obj->GetEntriesFast());
836+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
837+
spdg = obj->At(i)->GetName();
838+
genCocktailEvtGen->AddPdg(std::stoi(spdg), i);
839+
printf("PDG %d \n", std::stoi(spdg));
840+
}
841+
genCocktailEvtGen->SetForceDecay(kEvtDiMuon);
842+
843+
return genCocktailEvtGen;
844+
}
845+
674846

675847

676848

‎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
@@ -42,6 +42,9 @@ public:
4242
case 6: // generate prompt ChiC1 + ChiC2 cocktail at midrapidity
4343
mGeneratorParam = (Generator*)GeneratorCocktailChiCToElectronEvtGen_pp13TeV();
4444
break;
45+
case 7: // generate prompt charmonia cocktail at forward rapidity at 5TeV
46+
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_PbPb5TeV();
47+
break;
4548

4649
}
4750
mGeneratorParam->Init();

‎MC/config/PWGDQ/ini/Generator_InjectedPromptCharmoniaFwdy_TriggerGap_PbPb5TeV.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_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_withInjectedPromptSignals_gaptriggered_dq.C
4+
funcName=GeneratorPythia8InjectedPromptCharmoniaGapTriggered(5,7)
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap_5TeV.cfg

‎MC/config/PWGDQ/ini/tests/Generator_InjectedPromptCharmoniaFwdy_TriggerGap_PbPb5TeV.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+
}

‎MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap_5TeV.cfg‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
### beams
2+
Beams:idA 2212 # proton
3+
Beams:idB 2212 # proton
4+
Beams:eCM 5360. # GeV
5+
6+
### processes
7+
SoftQCD:inelastic on # all inelastic processes
8+
9+
### decays
10+
ParticleDecays:limitTau0 on
11+
ParticleDecays:tau0Max 10.

0 commit comments

Comments
 (0)