Skip to content

GeneratorSpectators: add proton generation, divergence and crossing-angle sampling + fixes #31

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fcatalan92 wants to merge 8 commits into
base: master
Choose a base branch
from
+156 −44
Open

GeneratorSpectators: add proton generation, divergence and crossing-angle sampling + fixes #31

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
a21c3d5
Add functions to generate protons
fcatalan92 Dec 5, 2025
0d7ccc7
Fully expose configuration in function
fcatalan92 May 6, 2026
f4f7d79
Add possibility to sample beam divergence and beam crossing angle fro…
fcatalan92 May 6, 2026
5d270fc
Fix sampling of beam-crossing angle
fcatalan92 Jul 2, 2026
3d8ce01
Fix crossing-angle calculation and use plane information
fcatalan92 Jul 2, 2026
a44de8b
Fix minor issues
fcatalan92 Jul 2, 2026
2dd4c82
Make crossing plane configurable
fcatalan92 Jul 2, 2026
398f029
Minor fix
fcatalan92 Jul 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 91 additions & 9 deletions GeneratorSpectators/GeneratorSpectators.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,19 +6,22 @@
class GeneratorSpectatorsO2 : public o2::eventgen::GeneratorTGenerator {

public:
GeneratorSpectatorsO2() : GeneratorTGenerator("GeneratorSpectators", "GeneratorSpectators") {
GeneratorSpectatorsO2() : GeneratorTGenerator("GeneratorSpectators", "GeneratorSpectators")
{
spec = new GeneratorSpectators();
setTGenerator(spec);
};

~GeneratorSpectatorsO2() { delete spec; };

Bool_t Init() override {
Bool_t Init() override
{
spec->Init();
return true;
}

void updateHeader(o2::dataformats::MCEventHeader *eventHeader) final {
void updateHeader(o2::dataformats::MCEventHeader *eventHeader) final
{
// not updating the header, just printing particles for the time being
LOG(info) << Form("--- Number of particles: %zu ---", mParticles.size());
for (int i = 0; i < mParticles.size(); ++i) {
Expand All @@ -36,25 +39,104 @@ private:
GeneratorSpectators *spec = nullptr;
};

FairGenerator *GeneratorSingleNeutron() {

FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
float beamDiv = 0.000032,
float beamDivMin = -999.,
float beamDivMax = -999.,
float beamCrossAngle = 0.,
float beamCrossAngleMin = -999.,
float beamCrossAngleMax = -999.,
int beamCrossPlane = 2)
{
auto wrap = new GeneratorSpectatorsO2();
auto spec = wrap->getGenerator();
spec->SetNpart(1);
spec->SetParticle(2112);
spec->SetMomentum(mom);
spec->SetDirection(0, 0., 0., 1.);
spec->SetDivergence(beamDiv);
spec->SetSampleDivergence(beamDivMin, beamDivMax);
spec->SetCrossing(beamCrossAngle, beamCrossPlane);
spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
return wrap;
}

FairGenerator *GeneratorSingleProton(float mom = 2680.,
float beamDiv = 0.000032,
float beamDivMin = -999.,
float beamDivMax = -999.,
float beamCrossAngle = 0.,
float beamCrossAngleMin = -999.,
float beamCrossAngleMax = -999.,
int beamCrossPlane = 2)
{
auto wrap = new GeneratorSpectatorsO2();
auto spec = wrap->getGenerator();
spec->SetNpart(1);
spec->SetMomentum(5360. / 2.);
spec->SetParticle(2212);
spec->SetMomentum(mom);
spec->SetDirection(0, 0., 0., 1.);
spec->SetDivergence(0.000032); // beam divergence in murad
spec->SetDivergence(beamDiv);
spec->SetSampleDivergence(beamDivMin, beamDivMax);
spec->SetCrossing(beamCrossAngle, beamCrossPlane);
spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
return wrap;
}

FairGenerator *GeneratorNeutrons(int nNeutrons = -1, float b = -1., bool useFluctuation = false) {
FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
float b = -1.,
bool useFluctuation = false,
float mom = 2680.,
float beamDiv = 0.000032,
float beamDivMin = -999.,
float beamDivMax = -999.,
float beamCrossAngle = 0.,
float beamCrossAngleMin = -999.,
float beamCrossAngleMax = -999.,
int beamCrossPlane = 2)
{
auto wrap = new GeneratorSpectatorsO2();
auto spec = wrap->getGenerator();
spec->SetParticle(2112);
spec->SetImpactParameter(b);
spec->SetNpartFluctuation(useFluctuation);
spec->SetNpart(nNeutrons);
spec->SetMomentum(5360. / 2.);
spec->SetMomentum(mom);
spec->SetDirection(0, 0., 0., 1.);
spec->SetDivergence(beamDiv);
spec->SetSampleDivergence(beamDivMin, beamDivMax);
spec->SetCrossing(beamCrossAngle, beamCrossPlane);
spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
return wrap;
}

FairGenerator *GeneratorProtons(int nProtons = -1,
float b = -1.,
bool useFluctuation = false,
float mom = 2680.,
float beamDiv = 0.000032,
float beamDivMin = -999.,
float beamDivMax = -999.,
float beamCrossAngle = 0.,
float beamCrossAngleMin = -999.,
float beamCrossAngleMax = -999.,
int beamCrossPlane = 2)
{
auto wrap = new GeneratorSpectatorsO2();
auto spec = wrap->getGenerator();
spec->SetParticle(2212);
spec->SetImpactParameter(b);
// TODO: implement parameterizations for the number of protons produced
// At the moment, those for the neutrons are used
spec->SetNpartFluctuation(useFluctuation);
spec->SetNpart(nProtons);
spec->SetMomentum(mom);
spec->SetDirection(0, 0., 0., 1.);
spec->SetDivergence(0.000032); // beam divergence in murad
spec->SetDivergence(beamDiv);
spec->SetSampleDivergence(beamDivMin, beamDivMax);
spec->SetCrossing(beamCrossAngle, beamCrossPlane);
spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
return wrap;
}

82 changes: 50 additions & 32 deletions GeneratorSpectators/GeneratorSpectators.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,9 @@ GeneratorSpectators::GeneratorSpectators()
SetDirection();
SetFermi();
SetDivergence();
SetSampleDivergence();
SetCrossing();
SetSampleCrossing();

for(Int_t i=0; i<201; i++){
fProbintp[i] = 0;
Expand All @@ -49,14 +51,15 @@ GeneratorSpectators::GeneratorSpectators()

void GeneratorSpectators::Init() {
printf("\n **** GeneratorSpectators initialization:\n");
printf(" Impact parameter: %f fm\n", fImpactParameter);
printf(" Number of particles to be generated (overwrites estimation from impact parameter): %d\n", fNpart);
printf(" Particle PDG: %d, Track: cosx = %f cosy = %f cosz = %f \n", fPDGcode, fCosx, fCosy, fCosz);
printf(" Maximum momentum: %f MeV/c", fPtot);
printf(" Fermi flag: %d, Beam divergence: %f, Crossing angle: %f, plane: %d\n\n",
fFermiflag, fBeamDiv, fBeamCrossAngle, fBeamCrossPlane);
// Initialize Fermi momentum distributions for Pb-Pb
FermiTwoGaussian(208.);
printf(" Impact parameter: %f fm\n", fImpactParameter);
printf(" Number of particles to be generated (overwrites estimation from impact parameter): %d\n", fNpart);
printf(" Particle PDG: %d, Track: cosx = %f cosy = %f cosz = %f\n", fPDGcode, fCosx, fCosy, fCosz);
printf(" Maximum momentum: %f GeV/c, Fermi flag: %d\n", fPtot, fFermiflag);
printf(" Beam divergence: %f, Crossing angle: %f, plane: %d\n", fBeamDiv, fBeamCrossAngle, fBeamCrossPlane);
printf(" Sample beam divergence: %f-%f (overwrites beam divergence)\n", fBeamDivMin, fBeamDivMax);
printf(" Sample crossing angle: %f-%f (overwrites crossing angle)\n", fBeamCrossAngleMin, fBeamCrossAngleMax);

FermiTwoGaussian(208.); // Initialize Fermi momentum distributions for Pb-Pb
}

void GeneratorSpectators::GenerateEvent() {
Expand Down Expand Up @@ -88,24 +91,24 @@ void GeneratorSpectators::GenerateEvent() {
pLab[2] = ptot * TMath::Cos(scang);
}

for (int i = 0; i < 3; i++)
fP[i] = pLab[i];
for (int j = 0; j < 3; j++)
fP[j] = pLab[j];

if (fDebug) {
printf("\n Particle momentum before divergence and crossing: ");
printf(" pLab = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
printf("pLab = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
}

// Beam divergence and crossing angle
if (TMath::Abs(fBeamCrossAngle) > 0.) {
if (TMath::Abs(fBeamCrossAngle) > 0. || (fBeamCrossAngleMin > -999. && fBeamCrossAngleMax > -999.)) {
BeamCrossing(pLab);
for (int i = 0; i < 3; i++)
fP[i] = pLab[i];
for (int j = 0; j < 3; j++)
fP[j] = pLab[j];
}
if (TMath::Abs(fBeamDiv) > 0.) {
if (TMath::Abs(fBeamDiv) > 0. || (fBeamDivMin >= 0. && fBeamDivMax > 0.)) {
BeamDivergence(pLab);
for (int i = 0; i < 3; i++)
fP[i] = pLab[i];
for (int j = 0; j < 3; j++)
fP[j] = pLab[j];
}

Double_t mass = TDatabasePDG::Instance()->GetParticle(fPDGcode)->Mass();
Expand All @@ -120,21 +123,21 @@ void GeneratorSpectators::GenerateEvent() {
ExtractFermi(fPDGcode, ddp);
fP0 = TMath::Sqrt(fP[0] * fP[0] + fP[1] * fP[1] + fP[2] * fP[2] +
mass * mass);
for (int i = 0; i < 3; i++)
dddp[i] = ddp[i];
for (int j = 0; j < 3; j++)
dddp[j] = ddp[j];
dddp0 = TMath::Sqrt(dddp[0] * dddp[0] + dddp[1] * dddp[1] +
dddp[2] * dddp[2] + mass * mass);

TVector3 b(fP[0] / fP0, fP[1] / fP0, fP[2] / fP0);
TLorentzVector pFermi(dddp[0], dddp[1], dddp[2], dddp0);
pFermi.Boost(b);

for (int i = 0; i < 3; i++)
fP[i] = pFermi[i];
for (int j = 0; j < 3; j++)
fP[j] = pFermi[j];
}

if (fDebug)
printf(" ### Particle momentum = (%f, %f, %f)\n", fP[0], fP[1], fP[2]);
printf(" Particle momentum = (%f, %f, %f)\n", fP[0], fP[1], fP[2]);

Double_t energy = TMath::Sqrt(fP[0] * fP[0] + fP[1] * fP[1] +
fP[2] * fP[2] + mass * mass);
Expand Down Expand Up @@ -207,7 +210,7 @@ void GeneratorSpectators::FermiTwoGaussian(Float_t A) {
fProbintn[i] = fProbintp[i];
}
if (fDebug)
printf(" Initialization of Fermi momenta distribution \n");
printf(" Initialization of Fermi momenta distribution \n");
}

void GeneratorSpectators::ExtractFermi(Int_t id, Double_t *ddp) {
Expand Down Expand Up @@ -242,19 +245,32 @@ void GeneratorSpectators::ExtractFermi(Int_t id, Double_t *ddp) {

void GeneratorSpectators::BeamCrossing(Double_t *pLab)
{
// Applying beam crossing angle
pLab[1] = pLab[2]*TMath::Sin(fBeamCrossAngle)+pLab[1]*TMath::Cos(fBeamCrossAngle);
pLab[2] = pLab[2] * TMath::Cos(fBeamCrossAngle) -
pLab[1] * TMath::Sin(fBeamCrossAngle);
// Sample beam crossing angle if set (angle can be negative, so -999 marks "not set")
if (fBeamCrossAngleMin > -999. && fBeamCrossAngleMax > -999.) {
fBeamCrossAngle = gRandom->Uniform(fBeamCrossAngleMin, fBeamCrossAngleMax);
}

// Applying beam crossing angle: rotate the plane selected by fBeamCrossPlane
// (=1 -> horizontal, x-z plane; =2 -> vertical, y-z plane)
Int_t idx = (fBeamCrossPlane == 1) ? 0 : 1;
Double_t p0 = pLab[idx];
Double_t z0 = pLab[2];
pLab[idx] = z0 * TMath::Sin(fBeamCrossAngle) + p0 * TMath::Cos(fBeamCrossAngle);
pLab[2] = z0 * TMath::Cos(fBeamCrossAngle) - p0 * TMath::Sin(fBeamCrossAngle);

if (fDebug) {
printf(" Beam crossing angle = %f mrad -> ", fBeamCrossAngle * 1000.);
printf(" p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
printf("p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
}
}

void GeneratorSpectators::BeamDivergence(Double_t *pLab)
{
// Sample beam divergence if set
if (fBeamDivMin >= 0. && fBeamDivMax > 0.) {
fBeamDiv = gRandom->Uniform(fBeamDivMin, fBeamDivMax);
}

// Applying beam divergence and crossing angle
Double_t pmq = 0.;
for (int i = 0; i < 3; i++)
Expand Down Expand Up @@ -292,7 +308,7 @@ void GeneratorSpectators::BeamDivergence(Double_t *pLab)

if (fDebug) {
printf(" Beam divergence = %f mrad -> ", fBeamDiv * 1000.);
printf(" p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
printf("p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]);
}
}

Expand All @@ -316,10 +332,11 @@ void GeneratorSpectators::AddAngle(Double_t theta1, Double_t phi1,

Double_t rtetsum = TMath::ACos(cz);
Double_t tetsum = conv*rtetsum;
Double_t fisum = 0;
Double_t fisum = 0.;

if (tetsum == 0. || tetsum == 180.) {
fisum = 0.;
angleSum[0] = tetsum;
angleSum[1] = fisum;
return;
}

Expand All @@ -332,6 +349,7 @@ void GeneratorSpectators::AddAngle(Double_t theta1, Double_t phi1,
fisum = conv*TMath::ACos(temp);
if (cy < 0)
fisum = 360. - fisum;

angleSum[0] = tetsum;
angleSum[1] = fisum;
}
Expand All @@ -357,7 +375,7 @@ void GeneratorSpectators::InitParameterizations() {
Double_t GeneratorSpectators::ImpParFunc(Double_t *x, Double_t *par) {
// Function for parameterization of impact parameter distribution
// from fit to Pb-Pb at 5.02 TeV
if (x[0] < 5.e-2 or x[0] > 16.5)
if (x[0] < 5.e-2 || x[0] > 16.5)
return 0.0;

if (x[0] < 13.8)
Expand Down
18 changes: 15 additions & 3 deletions GeneratorSpectators/GeneratorSpectators.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,18 @@ class GeneratorSpectators : public TGenerator {
}
void SetFermi(Bool_t flag = kTRUE) { fFermiflag = flag; }
void SetDivergence(Float_t bmdiv = 0.) { fBeamDiv = bmdiv; }
void SetSampleDivergence(Float_t min = -999., Float_t max = -999.) {
fBeamDivMin = min;
fBeamDivMax = max;
}
void SetCrossing(Float_t xingangle = 0., Int_t xingplane = 2) {
fBeamCrossAngle = xingangle;
fBeamCrossPlane = xingplane;
}
void SetSampleCrossing(Float_t min = -999., Float_t max = -999.) {
fBeamCrossAngleMin = min;
fBeamCrossAngleMax = max;
}

// Getters
Double_t GetFermi2p(Int_t key) const { return fProbintp[key]; }
Expand All @@ -54,10 +62,10 @@ class GeneratorSpectators : public TGenerator {

protected:
Bool_t fDebug; // Debugging flag
Float_t fImpactParameter; // Impact parameter, if <=0 sample from realistic distribution
Float_t fImpactParameter; // Impact parameter, if <= 0 sample from realistic distribution
Bool_t fNpartFluctuation; // Enable gaussian fluctuation of sampled number of particles
Int_t fNpart; // Number of particles to be generated
// if>0 overwrite the impact-parameter information
// if > 0 overwrite the impact-parameter information
Int_t fPDGcode; // Particle to be generated - can be n (2112) or p (2212)
Float_t fPtot; // Nucleon momentum
Float_t fPseudoRapidity; // Pseudorapidity: =0->track director cosines, !=0->pc.eta
Expand All @@ -66,7 +74,11 @@ class GeneratorSpectators : public TGenerator {
Float_t fCosz; // Director cos of the track - z direction
Bool_t fFermiflag; // Fermi momentum flag (true -> Fermi smearing)
Float_t fBeamDiv; // Beam divergence (angle in rad)
Float_t fBeamCrossAngle; // Beam crossing angle (angle in rad)
Float_t fBeamDivMin; // Minimum beam divergence, if >= 0 sample from uniform distribution
Float_t fBeamDivMax; // Maximum beam divergence, if > 0 sample from uniform distribution
Float_t fBeamCrossAngle; // Beam crossing angle (angle in rad), can be negative
Float_t fBeamCrossAngleMin; // Minimum beam crossing angle, if > -999 sample from uniform distribution
Float_t fBeamCrossAngleMax; // Maximum beam crossing angle, if > -999 sample from uniform distribution
Int_t fBeamCrossPlane; // Beam crossing plane
// (=1 -> horizontal, =2 -> vertical plane)

Expand Down