Skip to content

Commit 20ee6e4

Browse files
nburmasonburmaso
committed
Simplify parameterized single mu generator
1 parent 247bd31 commit 20ee6e4

File tree

1 file changed

+120
-95
lines changed

1 file changed

+120
-95
lines changed

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

Lines changed: 120 additions & 95 deletions
Original file line numberDiff line numberDiff line change
@@ -1,59 +1,29 @@
11
// Parameterized generator for muons
2-
// Port of the Run 2 generator by P. Pillot:
2+
// Adaptation of the Run 2 generator by P. Pillot:
33
// https://github.com/alisw/AliDPG/blob/master/MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C
44

5-
// clang-format off
6-
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen)
7-
// clang-format on
8-
R__LOAD_LIBRARY(libpythia6)
5+
#include "FairGenerator.h"
6+
#include "TF1.h"
7+
#include "TRandom.h"
8+
#include "TDatabasePDG.h"
9+
#include "TParticlePDG.h"
910

10-
#include "GeneratorEvtGen.C"
11-
12-
namespace o2
13-
{
14-
namespace eventgen
15-
{
16-
17-
class O2_GeneratorParamMuon : public GeneratorTGenerator
11+
class O2_GeneratorParamMuon : public FairGenerator
1812
{
1913
public:
20-
// parameters tuned to Pb-Pb @ 5.02 TeV
21-
inline static Double_t fPtP0{797.446};
22-
inline static Double_t fPtP1{0.830278};
23-
inline static Double_t fPtP2{0.632177};
24-
inline static Double_t fPtP3{10.2202};
25-
inline static Double_t fPtP4{-0.000614809};
26-
inline static Double_t fPtP5{-1.70993};
27-
28-
inline static Double_t fYP0{1.87732};
29-
inline static Double_t fYP1{0.00658212};
30-
inline static Double_t fYP2{-0.0988071};
31-
inline static Double_t fYP3{-0.000452746};
32-
inline static Double_t fYP4{0.00269782};
33-
34-
O2_GeneratorParamMuon() : GeneratorTGenerator("ParamMuon")
14+
O2_GeneratorParamMuon(int npart = 2, int pdg = 13, double ymin = -4.3, double ymax = -2.3, double ptmin = 0.1, double ptmax = 999.)
15+
: FairGenerator("ParaMuon"), fNParticles(npart), fPDGCode(pdg), fYMin(ymin), fYMax(ymax), fPtMin(ptmin), fPtMax(ptmax)
3516
{
36-
fParamMuon = new GeneratorParam(2, -1, PtMuon, YMuon, V2Muon, IpMuon);
37-
fParamMuon->SetPtRange(0.1, 999.);
38-
fParamMuon->SetYRange(-4.3, -2.3);
39-
fParamMuon->SetPhiRange(0., 360.);
40-
fParamMuon->SetDecayer(new TPythia6Decayer()); // "no decayer" error otherwise
41-
fParamMuon->SetForceDecay(kNoDecay);
42-
setTGenerator(fParamMuon);
17+
TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(fPDGCode);
18+
fMass = particle->Mass();
19+
fMass2 = fMass * fMass;
4320
}
4421

45-
~O2_GeneratorParamMuon() { delete fParamMuon; };
22+
~O2_GeneratorParamMuon() = default;
4623

47-
Bool_t Init() override
48-
{
49-
GeneratorTGenerator::Init();
50-
fParamMuon->Init();
51-
return true;
52-
}
24+
void SetRandomCharge(bool flag) { fRandomizeCharge = flag; }
5325

54-
// for tuning steps
55-
static void SetPtPars(Double_t p0, Double_t p1, Double_t p2, Double_t p3,
56-
Double_t p4, Double_t p5)
26+
void SetPtPars(double p0, double p1, double p2, double p3, double p4, double p5)
5727
{
5828
fPtP0 = p0;
5929
fPtP1 = p1;
@@ -63,9 +33,7 @@ class O2_GeneratorParamMuon : public GeneratorTGenerator
6333
fPtP5 = p5;
6434
}
6535

66-
// for tuning steps
67-
static void SetYPars(Double_t p0, Double_t p1, Double_t p2, Double_t p3,
68-
Double_t p4)
36+
void SetYPars(double p0, double p1, double p2, double p3, double p4)
6937
{
7038
fYP0 = p0;
7139
fYP1 = p1;
@@ -74,68 +42,125 @@ class O2_GeneratorParamMuon : public GeneratorTGenerator
7442
fYP4 = p4;
7543
}
7644

77-
void SetNSignalPerEvent(Int_t nsig) { fParamMuon->SetNumberParticles(nsig); }
78-
79-
// muon composition
80-
static Int_t IpMuon(TRandom* ran)
45+
void InitParaFuncs()
8146
{
82-
if (ran->Rndm() < 0.5) {
83-
return 13;
84-
} else {
85-
return -13;
86-
}
47+
fPtPara = new TF1("pt-para", PtMuon, fPtMin, fPtMax, 6);
48+
fPtPara->SetParameter(0, fPtP0);
49+
fPtPara->SetParameter(1, fPtP1);
50+
fPtPara->SetParameter(2, fPtP2);
51+
fPtPara->SetParameter(3, fPtP3);
52+
fPtPara->SetParameter(4, fPtP4);
53+
fPtPara->SetParameter(5, fPtP5);
54+
fYPara = new TF1("y-para", YMuon, fYMin, fYMax, 5);
55+
fYPara->SetParameter(0, fYP0);
56+
fYPara->SetParameter(1, fYP1);
57+
fYPara->SetParameter(2, fYP2);
58+
fYPara->SetParameter(3, fYP3);
59+
fYPara->SetParameter(4, fYP4);
8760
}
8861

89-
// muon pT
90-
static Double_t PtMuon(const Double_t* px, const Double_t*)
62+
static double PtMuon(const double* xx, const double* par)
9163
{
92-
Double_t x = px[0];
93-
return fPtP0 * (1. / TMath::Power(fPtP1 + TMath::Power(x, fPtP2), fPtP3) +
94-
fPtP4 * TMath::Exp(fPtP5 * x));
64+
double x = xx[0];
65+
double p0 = par[0];
66+
double p1 = par[1];
67+
double p2 = par[2];
68+
double p3 = par[3];
69+
double p4 = par[4];
70+
double p5 = par[5];
71+
return p0 * (1. / std::pow(p1 + std::pow(x, p2), p3) + p4 * std::exp(p5 * x));
9572
}
9673

97-
// muon y
98-
static Double_t YMuon(const Double_t* py, const Double_t*)
74+
static double YMuon(const double* xx, const double* par)
9975
{
100-
Double_t y = py[0];
101-
return fYP0 * (1. + fYP1 * y + fYP2 * y * y + fYP3 * y * y * y +
102-
fYP4 * y * y * y * y);
76+
double x = xx[0];
77+
double p0 = par[0];
78+
double p1 = par[1];
79+
double p2 = par[2];
80+
double p3 = par[3];
81+
double p4 = par[4];
82+
double x2 = x * x;
83+
return p0 * (1. + p1 * x + p2 * x2 + p3 * x2 * x + p4 * x2 * x2);
10384
}
10485

105-
static Double_t V2Muon(const Double_t*, const Double_t*)
86+
bool ReadEvent(FairPrimaryGenerator* primGen) override
10687
{
107-
// muon v2
108-
return 0.;
88+
int iPart = fNParticles;
89+
// no kinematic cuts -> accepting all
90+
for (int i = 0; i < fNParticles; ++i) {
91+
double pt = fPtPara->GetRandom();
92+
double ty = std::tanh(fYPara->GetRandom());
93+
double xmt = std::sqrt(pt * pt + fMass2);
94+
double pl = xmt * ty / std::sqrt((1. - ty * ty));
95+
double phi = gRandom->Uniform(0., 2. * M_PI);
96+
double px = pt * std::cos(phi);
97+
double py = pt * std::sin(phi);
98+
int pdg = fPDGCode;
99+
if (fRandomizeCharge) {
100+
int charge;
101+
if (gRandom->Rndm() < 0.5) {
102+
charge = 1;
103+
} else {
104+
charge = -1;
105+
}
106+
pdg = std::abs(pdg) * charge;
107+
}
108+
primGen->AddTrack(pdg, px, py, pl, 0, 0, 0);
109+
printf("Add track %d %.2f %.2f %.2f \n", pdg, px, py, pl);
110+
}
111+
return kTRUE;
109112
}
110113

111114
private:
112-
GeneratorParam* fParamMuon{nullptr};
113-
};
115+
TF1* fPtPara{nullptr};
116+
TF1* fYPara{nullptr};
117+
TF1* fdNdPhi{nullptr};
114118

115-
} // namespace eventgen
116-
} // namespace o2
119+
// parameters tuned to Pb-Pb @ 5.02 TeV
120+
double fPtP0{797.446};
121+
double fPtP1{0.830278};
122+
double fPtP2{0.632177};
123+
double fPtP3{10.2202};
124+
double fPtP4{-0.000614809};
125+
double fPtP5{-1.70993};
126+
127+
double fYP0{1.87732};
128+
double fYP1{0.00658212};
129+
double fYP2{-0.0988071};
130+
double fYP3{-0.000452746};
131+
double fYP4{0.00269782};
132+
133+
const double fDeltaPt{0.01};
134+
135+
// configuration
136+
int fPDGCode{13};
137+
int fNParticles{2};
138+
double fYMin{-4.3};
139+
double fYMax{-2.3};
140+
double fPtMin{0.1};
141+
double fPtMax{999.};
142+
const double fPhiMin{0.};
143+
const double fPhiMax{2. * M_PI};
144+
bool fRandomizeCharge{true};
145+
double fMass{0.10566};
146+
double fMass2{0};
147+
};
117148

118-
FairGenerator* paramMuGen(Double_t ptP0 = 797.446, Double_t ptP1 = 0.830278,
119-
Double_t ptP2 = 0.632177, Double_t ptP3 = 10.2202,
120-
Double_t ptP4 = -0.000614809, Double_t ptP5 = -1.70993,
121-
Double_t yP0 = 1.87732, Double_t yP1 = 0.00658212,
122-
Double_t yP2 = -0.0988071, Double_t yP3 = -0.000452746,
123-
Double_t yP4 = 0.00269782, Int_t nMuons = 2, TString pdgs = "13")
149+
FairGenerator* paramMuGen(double ptP0 = 797.446, double ptP1 = 0.830278,
150+
double ptP2 = 0.632177, double ptP3 = 10.2202,
151+
double ptP4 = -0.000614809, double ptP5 = -1.70993,
152+
double yP0 = 1.87732, double yP1 = 0.00658212,
153+
double yP2 = -0.0988071, double yP3 = -0.000452746,
154+
double yP4 = 0.00269782,
155+
int nPart = 2, int pdg = 13,
156+
double ymin = -4.3, double ymax = -2.3,
157+
double ptmin = 0.1, float ptmax = 999.,
158+
int randCharge = 1)
124159
{
125-
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamMuon>();
126-
o2::eventgen::O2_GeneratorParamMuon::SetPtPars(ptP0, ptP1, ptP2, ptP3, ptP4, ptP5);
127-
o2::eventgen::O2_GeneratorParamMuon::SetYPars(yP0, yP1, yP2, yP3, yP4);
128-
gen->SetNSignalPerEvent(nMuons); // number of muons per event
129-
130-
std::string spdg;
131-
TObjArray* obj = pdgs.Tokenize(";");
132-
gen->SetSizePdg(obj->GetEntriesFast());
133-
for (int i = 0; i < obj->GetEntriesFast(); i++) {
134-
spdg = obj->At(i)->GetName();
135-
gen->AddPdg(std::stoi(spdg), i);
136-
printf("PDG %d \n", std::stoi(spdg));
137-
}
138-
139-
gen->PrintDebug();
160+
auto* gen = new O2_GeneratorParamMuon(nPart, pdg, ymin, ymax, ptmin, ptmax);
161+
gen->SetPtPars(ptP0, ptP1, ptP2, ptP3, ptP4, ptP5);
162+
gen->SetYPars(yP0, yP1, yP2, yP3, yP4);
163+
gen->InitParaFuncs();
164+
gen->SetRandomCharge(randCharge);
140165
return gen;
141166
}

0 commit comments

Comments
 (0)