Skip to content

Commit

Permalink
Simplify parameterized single mu generator
Browse files Browse the repository at this point in the history
  • Loading branch information
@nburmaso
nburmaso committed Feb 10, 2025
1 parent 247bd31 commit 20ee6e4
Showing 1 changed file with 120 additions and 95 deletions.
215 changes: 120 additions & 95 deletions MC/config/PWGDQ/external/generator/GeneratorParamSingleMuon_PbPb_5TeV.C
View file
Original file line number Diff line number Diff line change
@@ -1,59 +1,29 @@
// Parameterized generator for muons
// Port of the Run 2 generator by P. Pillot:
// Adaptation of the Run 2 generator by P. Pillot:
// https://github.com/alisw/AliDPG/blob/master/MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C

// clang-format off
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen)
// clang-format on
R__LOAD_LIBRARY(libpythia6)
#include "FairGenerator.h"
#include "TF1.h"
#include "TRandom.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"

#include "GeneratorEvtGen.C"

namespace o2
{
namespace eventgen
{

class O2_GeneratorParamMuon : public GeneratorTGenerator
class O2_GeneratorParamMuon : public FairGenerator
{
public:
// parameters tuned to Pb-Pb @ 5.02 TeV
inline static Double_t fPtP0{797.446};
inline static Double_t fPtP1{0.830278};
inline static Double_t fPtP2{0.632177};
inline static Double_t fPtP3{10.2202};
inline static Double_t fPtP4{-0.000614809};
inline static Double_t fPtP5{-1.70993};

inline static Double_t fYP0{1.87732};
inline static Double_t fYP1{0.00658212};
inline static Double_t fYP2{-0.0988071};
inline static Double_t fYP3{-0.000452746};
inline static Double_t fYP4{0.00269782};

O2_GeneratorParamMuon() : GeneratorTGenerator("ParamMuon")
O2_GeneratorParamMuon(int npart = 2, int pdg = 13, double ymin = -4.3, double ymax = -2.3, double ptmin = 0.1, double ptmax = 999.)
: FairGenerator("ParaMuon"), fNParticles(npart), fPDGCode(pdg), fYMin(ymin), fYMax(ymax), fPtMin(ptmin), fPtMax(ptmax)
{
fParamMuon = new GeneratorParam(2, -1, PtMuon, YMuon, V2Muon, IpMuon);
fParamMuon->SetPtRange(0.1, 999.);
fParamMuon->SetYRange(-4.3, -2.3);
fParamMuon->SetPhiRange(0., 360.);
fParamMuon->SetDecayer(new TPythia6Decayer()); // "no decayer" error otherwise
fParamMuon->SetForceDecay(kNoDecay);
setTGenerator(fParamMuon);
TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(fPDGCode);
fMass = particle->Mass();
fMass2 = fMass * fMass;
}

~O2_GeneratorParamMuon() { delete fParamMuon; };
~O2_GeneratorParamMuon() = default;

Bool_t Init() override
{
GeneratorTGenerator::Init();
fParamMuon->Init();
return true;
}
void SetRandomCharge(bool flag) { fRandomizeCharge = flag; }

// for tuning steps
static void SetPtPars(Double_t p0, Double_t p1, Double_t p2, Double_t p3,
Double_t p4, Double_t p5)
void SetPtPars(double p0, double p1, double p2, double p3, double p4, double p5)
{
fPtP0 = p0;
fPtP1 = p1;
Expand All @@ -63,9 +33,7 @@ class O2_GeneratorParamMuon : public GeneratorTGenerator
fPtP5 = p5;
}

// for tuning steps
static void SetYPars(Double_t p0, Double_t p1, Double_t p2, Double_t p3,
Double_t p4)
void SetYPars(double p0, double p1, double p2, double p3, double p4)
{
fYP0 = p0;
fYP1 = p1;
Expand All @@ -74,68 +42,125 @@ class O2_GeneratorParamMuon : public GeneratorTGenerator
fYP4 = p4;
}

void SetNSignalPerEvent(Int_t nsig) { fParamMuon->SetNumberParticles(nsig); }

// muon composition
static Int_t IpMuon(TRandom* ran)
void InitParaFuncs()
{
if (ran->Rndm() < 0.5) {
return 13;
} else {
return -13;
}
fPtPara = new TF1("pt-para", PtMuon, fPtMin, fPtMax, 6);
fPtPara->SetParameter(0, fPtP0);
fPtPara->SetParameter(1, fPtP1);
fPtPara->SetParameter(2, fPtP2);
fPtPara->SetParameter(3, fPtP3);
fPtPara->SetParameter(4, fPtP4);
fPtPara->SetParameter(5, fPtP5);
fYPara = new TF1("y-para", YMuon, fYMin, fYMax, 5);
fYPara->SetParameter(0, fYP0);
fYPara->SetParameter(1, fYP1);
fYPara->SetParameter(2, fYP2);
fYPara->SetParameter(3, fYP3);
fYPara->SetParameter(4, fYP4);
}

// muon pT
static Double_t PtMuon(const Double_t* px, const Double_t*)
static double PtMuon(const double* xx, const double* par)
{
Double_t x = px[0];
return fPtP0 * (1. / TMath::Power(fPtP1 + TMath::Power(x, fPtP2), fPtP3) +
fPtP4 * TMath::Exp(fPtP5 * x));
double x = xx[0];
double p0 = par[0];
double p1 = par[1];
double p2 = par[2];
double p3 = par[3];
double p4 = par[4];
double p5 = par[5];
return p0 * (1. / std::pow(p1 + std::pow(x, p2), p3) + p4 * std::exp(p5 * x));
}

// muon y
static Double_t YMuon(const Double_t* py, const Double_t*)
static double YMuon(const double* xx, const double* par)
{
Double_t y = py[0];
return fYP0 * (1. + fYP1 * y + fYP2 * y * y + fYP3 * y * y * y +
fYP4 * y * y * y * y);
double x = xx[0];
double p0 = par[0];
double p1 = par[1];
double p2 = par[2];
double p3 = par[3];
double p4 = par[4];
double x2 = x * x;
return p0 * (1. + p1 * x + p2 * x2 + p3 * x2 * x + p4 * x2 * x2);
}

static Double_t V2Muon(const Double_t*, const Double_t*)
bool ReadEvent(FairPrimaryGenerator* primGen) override
{
// muon v2
return 0.;
int iPart = fNParticles;
// no kinematic cuts -> accepting all
for (int i = 0; i < fNParticles; ++i) {
double pt = fPtPara->GetRandom();
double ty = std::tanh(fYPara->GetRandom());
double xmt = std::sqrt(pt * pt + fMass2);
double pl = xmt * ty / std::sqrt((1. - ty * ty));
double phi = gRandom->Uniform(0., 2. * M_PI);
double px = pt * std::cos(phi);
double py = pt * std::sin(phi);
int pdg = fPDGCode;
if (fRandomizeCharge) {
int charge;
if (gRandom->Rndm() < 0.5) {
charge = 1;
} else {
charge = -1;
}
pdg = std::abs(pdg) * charge;
}
primGen->AddTrack(pdg, px, py, pl, 0, 0, 0);
printf("Add track %d %.2f %.2f %.2f \n", pdg, px, py, pl);
}
return kTRUE;
}

private:
GeneratorParam* fParamMuon{nullptr};
};
TF1* fPtPara{nullptr};
TF1* fYPara{nullptr};
TF1* fdNdPhi{nullptr};

} // namespace eventgen
} // namespace o2
// parameters tuned to Pb-Pb @ 5.02 TeV
double fPtP0{797.446};
double fPtP1{0.830278};
double fPtP2{0.632177};
double fPtP3{10.2202};
double fPtP4{-0.000614809};
double fPtP5{-1.70993};

double fYP0{1.87732};
double fYP1{0.00658212};
double fYP2{-0.0988071};
double fYP3{-0.000452746};
double fYP4{0.00269782};

const double fDeltaPt{0.01};

// configuration
int fPDGCode{13};
int fNParticles{2};
double fYMin{-4.3};
double fYMax{-2.3};
double fPtMin{0.1};
double fPtMax{999.};
const double fPhiMin{0.};
const double fPhiMax{2. * M_PI};
bool fRandomizeCharge{true};
double fMass{0.10566};
double fMass2{0};
};

FairGenerator* paramMuGen(Double_t ptP0 = 797.446, Double_t ptP1 = 0.830278,
Double_t ptP2 = 0.632177, Double_t ptP3 = 10.2202,
Double_t ptP4 = -0.000614809, Double_t ptP5 = -1.70993,
Double_t yP0 = 1.87732, Double_t yP1 = 0.00658212,
Double_t yP2 = -0.0988071, Double_t yP3 = -0.000452746,
Double_t yP4 = 0.00269782, Int_t nMuons = 2, TString pdgs = "13")
FairGenerator* paramMuGen(double ptP0 = 797.446, double ptP1 = 0.830278,
double ptP2 = 0.632177, double ptP3 = 10.2202,
double ptP4 = -0.000614809, double ptP5 = -1.70993,
double yP0 = 1.87732, double yP1 = 0.00658212,
double yP2 = -0.0988071, double yP3 = -0.000452746,
double yP4 = 0.00269782,
int nPart = 2, int pdg = 13,
double ymin = -4.3, double ymax = -2.3,
double ptmin = 0.1, float ptmax = 999.,
int randCharge = 1)
{
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::O2_GeneratorParamMuon>();
o2::eventgen::O2_GeneratorParamMuon::SetPtPars(ptP0, ptP1, ptP2, ptP3, ptP4, ptP5);
o2::eventgen::O2_GeneratorParamMuon::SetYPars(yP0, yP1, yP2, yP3, yP4);
gen->SetNSignalPerEvent(nMuons); // number of muons per event

std::string spdg;
TObjArray* obj = pdgs.Tokenize(";");
gen->SetSizePdg(obj->GetEntriesFast());
for (int i = 0; i < obj->GetEntriesFast(); i++) {
spdg = obj->At(i)->GetName();
gen->AddPdg(std::stoi(spdg), i);
printf("PDG %d \n", std::stoi(spdg));
}

gen->PrintDebug();
auto* gen = new O2_GeneratorParamMuon(nPart, pdg, ymin, ymax, ptmin, ptmax);
gen->SetPtPars(ptP0, ptP1, ptP2, ptP3, ptP4, ptP5);
gen->SetYPars(yP0, yP1, yP2, yP3, yP4);
gen->InitParaFuncs();
gen->SetRandomCharge(randCharge);
return gen;
}

0 comments on commit 20ee6e4

Please sign in to comment.