From a21c3d57b7b0b3839986d636224e1d55efb4f508 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Fri, 5 Dec 2025 11:07:13 +0100
Subject: [PATCH 1/8] Add functions to generate protons

---
 GeneratorSpectators/GeneratorSpectators.C | 26 +++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index 28f2421..7eb541d 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -46,6 +46,17 @@ FairGenerator *GeneratorSingleNeutron() {
   return wrap;
 }
 
+FairGenerator *GeneratorSingleProton() {
+  auto wrap = new GeneratorSpectatorsO2();
+  auto spec = wrap->getGenerator();
+  spec->SetNpart(1);
+  spec->SetParticle(2212);
+  spec->SetMomentum(5360. / 2.);
+  spec->SetDirection(0, 0., 0., 1.);
+  spec->SetDivergence(0.000032);  // beam divergence in murad
+  return wrap;
+}
+
 FairGenerator *GeneratorNeutrons(int nNeutrons = -1, float b = -1., bool useFluctuation = false) {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -58,3 +69,18 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1, float b = -1., bool useFluc
   return wrap;
 }
 
+FairGenerator *GeneratorProtons(int nNeutrons = -1, float b = -1., bool useFluctuation = false) {
+  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(nNeutrons);
+  spec->SetMomentum(5360. / 2.);
+  spec->SetDirection(0, 0., 0., 1.);
+  spec->SetDivergence(0.000032);  // beam divergence in murad
+  return wrap;
+}
+

From 0d7ccc7f88b4fb85578d5f6e41fcf55f1fdd2a5b Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Wed, 6 May 2026 15:06:30 +0200
Subject: [PATCH 2/8] Fully expose configuration in function

---
 GeneratorSpectators/GeneratorSpectators.C | 50 ++++++++++++++++-------
 1 file changed, 35 insertions(+), 15 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index 7eb541d..2beeb66 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -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) {
@@ -36,40 +39,57 @@ private:
   GeneratorSpectators *spec = nullptr;
 };
 
-FairGenerator *GeneratorSingleNeutron() {
+
+FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
+                                      float beamDiv = 0.000032)
+{
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
   spec->SetNpart(1);
-  spec->SetMomentum(5360. / 2.);
+  spec->SetParticle(2112);
+  spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
-  spec->SetDivergence(0.000032);  // beam divergence in murad
+  spec->SetDivergence(beamDiv);
   return wrap;
 }
 
-FairGenerator *GeneratorSingleProton() {
+FairGenerator *GeneratorSingleProton(float mom = 2680.,
+                                     float beamDiv = 0.000032)
+{
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
   spec->SetNpart(1);
   spec->SetParticle(2212);
-  spec->SetMomentum(5360. / 2.);
+  spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
-  spec->SetDivergence(0.000032);  // beam divergence in murad
+  spec->SetDivergence(beamDiv);
   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)
+{
   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(0.000032);  // beam divergence in murad
+  spec->SetDivergence(beamDiv);
   return wrap;
 }
 
-FairGenerator *GeneratorProtons(int nNeutrons = -1, float b = -1., bool useFluctuation = false) {
+FairGenerator *GeneratorProtons(int nNeutrons = -1,
+                                float b = -1.,
+                                bool useFluctuation = false,
+                                float mom = 2680.,
+                                float beamDiv = 0.000032)
+{
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
   spec->SetParticle(2212);
@@ -78,9 +98,9 @@ FairGenerator *GeneratorProtons(int nNeutrons = -1, float b = -1., bool useFluct
   // At the moment, those for the neutrons are used
   spec->SetNpartFluctuation(useFluctuation);
   spec->SetNpart(nNeutrons);
-  spec->SetMomentum(5360. / 2.);
+  spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
-  spec->SetDivergence(0.000032);  // beam divergence in murad
+  spec->SetDivergence(beamDiv);
   return wrap;
 }
 

From f4f7d79356d026962f6c772a132577b2e2936c33 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Wed, 6 May 2026 15:36:21 +0200
Subject: [PATCH 3/8] Add possibility to sample beam divergence and beam
 crossing angle from an uniform distribution

---
 GeneratorSpectators/GeneratorSpectators.C   | 40 +++++++++++++++++--
 GeneratorSpectators/GeneratorSpectators.cxx | 43 ++++++++++++++-------
 GeneratorSpectators/GeneratorSpectators.h   | 16 +++++++-
 3 files changed, 78 insertions(+), 21 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index 2beeb66..b5321a8 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -41,7 +41,12 @@ private:
 
 
 FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
-                                      float beamDiv = 0.000032)
+                                      float beamDiv = 0.000032,
+                                      float beamDivMin = -1.,
+                                      float beamDivMax = -1.,
+                                      float beamCrossAngle = 0.,
+                                      float beamCrossAngleMin = -1.,
+                                      float beamCrossAngleMax = -1.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -50,11 +55,19 @@ FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
   spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
+  spec->SetSampleDivergence(beamDivMin, beamDivMax);
+  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
 
 FairGenerator *GeneratorSingleProton(float mom = 2680.,
-                                     float beamDiv = 0.000032)
+                                     float beamDiv = 0.000032,
+                                     float beamDivMin = -1.,
+                                     float beamDivMax = -1.,
+                                     float beamCrossAngle = 0.,
+                                     float beamCrossAngleMin = -1.,
+                                     float beamCrossAngleMax = -1.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -63,6 +76,9 @@ FairGenerator *GeneratorSingleProton(float mom = 2680.,
   spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
+  spec->SetSampleDivergence(beamDivMin, beamDivMax);
+  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
 
@@ -70,7 +86,12 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
                                  float b = -1.,
                                  bool useFluctuation = false,
                                  float mom = 2680.,
-                                 float beamDiv = 0.000032)
+                                 float beamDiv = 0.000032,
+                                 float beamDivMin = -1.,
+                                 float beamDivMax = -1.,
+                                 float beamCrossAngle = 0.,
+                                 float beamCrossAngleMin = -1.,
+                                 float beamCrossAngleMax = -1.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -81,6 +102,9 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
   spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
+  spec->SetSampleDivergence(beamDivMin, beamDivMax);
+  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
 
@@ -88,7 +112,12 @@ FairGenerator *GeneratorProtons(int nNeutrons = -1,
                                 float b = -1.,
                                 bool useFluctuation = false,
                                 float mom = 2680.,
-                                float beamDiv = 0.000032)
+                                float beamDiv = 0.000032,
+                                float beamDivMin = -1.,
+                                float beamDivMax = -1.,
+                                float beamCrossAngle = 0.,
+                                float beamCrossAngleMin = -1.,
+                                float beamCrossAngleMax = -1.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -101,6 +130,9 @@ FairGenerator *GeneratorProtons(int nNeutrons = -1,
   spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
+  spec->SetSampleDivergence(beamDivMin, beamDivMax);
+  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
 
diff --git a/GeneratorSpectators/GeneratorSpectators.cxx b/GeneratorSpectators/GeneratorSpectators.cxx
index 5737196..f595705 100644
--- a/GeneratorSpectators/GeneratorSpectators.cxx
+++ b/GeneratorSpectators/GeneratorSpectators.cxx
@@ -36,7 +36,9 @@ GeneratorSpectators::GeneratorSpectators()
   SetDirection();
   SetFermi();
   SetDivergence();
+  SetSampleDivergence();
   SetCrossing();
+  SetSampleCrossing();
 
   for(Int_t i=0; i<201; i++){
      fProbintp[i] = 0;
@@ -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 MeV/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() {
@@ -93,16 +96,16 @@ void GeneratorSpectators::GenerateEvent() {
 
     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 >= 0. && fBeamCrossAngleMax > 0.)) {
       BeamCrossing(pLab);
       for (int i = 0; i < 3; i++)
         fP[i] = pLab[i];
     }
-    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];
@@ -134,7 +137,7 @@ void GeneratorSpectators::GenerateEvent() {
     }
 
     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);
@@ -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) {
@@ -242,6 +245,11 @@ void GeneratorSpectators::ExtractFermi(Int_t id, Double_t *ddp) {
 
 void GeneratorSpectators::BeamCrossing(Double_t *pLab)
 {
+  // Sample beam crossing angle if set
+  if (fBeamCrossAngleMin >= 0. && fBeamCrossAngleMax > 0.) {
+    fBeamCrossAngle = gRandom->Uniform(fBeamCrossAngleMin, fBeamCrossAngleMax);
+  }
+
   // Applying beam crossing angle
   pLab[1] = pLab[2]*TMath::Sin(fBeamCrossAngle)+pLab[1]*TMath::Cos(fBeamCrossAngle);
   pLab[2] = pLab[2] * TMath::Cos(fBeamCrossAngle) -
@@ -249,12 +257,17 @@ void GeneratorSpectators::BeamCrossing(Double_t *pLab)
 
   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++)
@@ -292,7 +305,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]);
   }
 }
 
diff --git a/GeneratorSpectators/GeneratorSpectators.h b/GeneratorSpectators/GeneratorSpectators.h
index 561e47d..8a49433 100644
--- a/GeneratorSpectators/GeneratorSpectators.h
+++ b/GeneratorSpectators/GeneratorSpectators.h
@@ -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 = -1., Float_t max = -1.) {
+    fBeamDivMin = min;
+    fBeamDivMax = max;
+  }
   void SetCrossing(Float_t xingangle = 0., Int_t xingplane = 2) {
     fBeamCrossAngle = xingangle;
     fBeamCrossPlane = xingplane;
   }
+  void SetSampleCrossing(Float_t min = -1., Float_t max = -1.) {
+    fBeamCrossAngleMin = min;
+    fBeamCrossAngleMax = max;
+  }
 
   // Getters
   Double_t GetFermi2p(Int_t key) const { return fProbintp[key]; }
@@ -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
@@ -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  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)
+  Float_t  fBeamCrossAngleMin;  // Minimum beam crossing angle, if >= 0 sample from uniform distribution
+  Float_t  fBeamCrossAngleMax;  // Maximum beam crossing angle, if > 0 sample from uniform distribution
   Int_t    fBeamCrossPlane;     // Beam crossing plane
                                 // (=1 -> horizontal, =2 -> vertical plane)
 

From 5d270fc92312eb0d91ab78eef925ec390f39d2c1 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Thu, 2 Jul 2026 15:38:51 +0200
Subject: [PATCH 4/8] Fix sampling of beam-crossing angle

---
 GeneratorSpectators/GeneratorSpectators.C   | 32 ++++++++++-----------
 GeneratorSpectators/GeneratorSpectators.cxx |  6 ++--
 GeneratorSpectators/GeneratorSpectators.h   | 10 +++----
 3 files changed, 24 insertions(+), 24 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index b5321a8..6a23869 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -42,11 +42,11 @@ private:
 
 FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
                                       float beamDiv = 0.000032,
-                                      float beamDivMin = -1.,
-                                      float beamDivMax = -1.,
+                                      float beamDivMin = -999.,
+                                      float beamDivMax = -999.,
                                       float beamCrossAngle = 0.,
-                                      float beamCrossAngleMin = -1.,
-                                      float beamCrossAngleMax = -1.)
+                                      float beamCrossAngleMin = -999.,
+                                      float beamCrossAngleMax = -999.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -63,11 +63,11 @@ FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
 
 FairGenerator *GeneratorSingleProton(float mom = 2680.,
                                      float beamDiv = 0.000032,
-                                     float beamDivMin = -1.,
-                                     float beamDivMax = -1.,
+                                     float beamDivMin = -999.,
+                                     float beamDivMax = -999.,
                                      float beamCrossAngle = 0.,
-                                     float beamCrossAngleMin = -1.,
-                                     float beamCrossAngleMax = -1.)
+                                     float beamCrossAngleMin = -999.,
+                                     float beamCrossAngleMax = -999.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -87,11 +87,11 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
                                  bool useFluctuation = false,
                                  float mom = 2680.,
                                  float beamDiv = 0.000032,
-                                 float beamDivMin = -1.,
-                                 float beamDivMax = -1.,
+                                 float beamDivMin = -999.,
+                                 float beamDivMax = -999.,
                                  float beamCrossAngle = 0.,
-                                 float beamCrossAngleMin = -1.,
-                                 float beamCrossAngleMax = -1.)
+                                 float beamCrossAngleMin = -999.,
+                                 float beamCrossAngleMax = -999.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -113,11 +113,11 @@ FairGenerator *GeneratorProtons(int nNeutrons = -1,
                                 bool useFluctuation = false,
                                 float mom = 2680.,
                                 float beamDiv = 0.000032,
-                                float beamDivMin = -1.,
-                                float beamDivMax = -1.,
+                                float beamDivMin = -999.,
+                                float beamDivMax = -999.,
                                 float beamCrossAngle = 0.,
-                                float beamCrossAngleMin = -1.,
-                                float beamCrossAngleMax = -1.)
+                                float beamCrossAngleMin = -999.,
+                                float beamCrossAngleMax = -999.)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
diff --git a/GeneratorSpectators/GeneratorSpectators.cxx b/GeneratorSpectators/GeneratorSpectators.cxx
index f595705..fff6d59 100644
--- a/GeneratorSpectators/GeneratorSpectators.cxx
+++ b/GeneratorSpectators/GeneratorSpectators.cxx
@@ -100,7 +100,7 @@ void GeneratorSpectators::GenerateEvent() {
     }
 
     // Beam divergence and crossing angle
-    if (TMath::Abs(fBeamCrossAngle) > 0. || (fBeamCrossAngleMin >= 0. && fBeamCrossAngleMax > 0.)) {
+    if (TMath::Abs(fBeamCrossAngle) > 0. || (fBeamCrossAngleMin > -999. && fBeamCrossAngleMax > -999.)) {
       BeamCrossing(pLab);
       for (int i = 0; i < 3; i++)
         fP[i] = pLab[i];
@@ -245,8 +245,8 @@ void GeneratorSpectators::ExtractFermi(Int_t id, Double_t *ddp) {
 
 void GeneratorSpectators::BeamCrossing(Double_t *pLab)
 {
-  // Sample beam crossing angle if set
-  if (fBeamCrossAngleMin >= 0. && fBeamCrossAngleMax > 0.) {
+  // 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);
   }
 
diff --git a/GeneratorSpectators/GeneratorSpectators.h b/GeneratorSpectators/GeneratorSpectators.h
index 8a49433..0792cfe 100644
--- a/GeneratorSpectators/GeneratorSpectators.h
+++ b/GeneratorSpectators/GeneratorSpectators.h
@@ -42,7 +42,7 @@ 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 = -1., Float_t max = -1.) {
+  void SetSampleDivergence(Float_t min = -999., Float_t max = -999.) {
     fBeamDivMin = min;
     fBeamDivMax = max;
   }
@@ -50,7 +50,7 @@ class GeneratorSpectators : public TGenerator {
     fBeamCrossAngle = xingangle;
     fBeamCrossPlane = xingplane;
   }
-  void SetSampleCrossing(Float_t min = -1., Float_t max = -1.) {
+  void SetSampleCrossing(Float_t min = -999., Float_t max = -999.) {
     fBeamCrossAngleMin = min;
     fBeamCrossAngleMax = max;
   }
@@ -76,9 +76,9 @@ class GeneratorSpectators : public TGenerator {
   Float_t  fBeamDiv;            // Beam divergence (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)
-  Float_t  fBeamCrossAngleMin;  // Minimum beam crossing angle, if >= 0 sample from uniform distribution
-  Float_t  fBeamCrossAngleMax;  // Maximum beam crossing angle, 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)
 

From 3d8ce0136323188b544a71ac6bcff53143cd5187 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Thu, 2 Jul 2026 16:21:55 +0200
Subject: [PATCH 5/8] Fix crossing-angle calculation and use plane information

---
 GeneratorSpectators/GeneratorSpectators.cxx | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.cxx b/GeneratorSpectators/GeneratorSpectators.cxx
index fff6d59..4090b2e 100644
--- a/GeneratorSpectators/GeneratorSpectators.cxx
+++ b/GeneratorSpectators/GeneratorSpectators.cxx
@@ -250,10 +250,13 @@ void GeneratorSpectators::BeamCrossing(Double_t *pLab)
     fBeamCrossAngle = gRandom->Uniform(fBeamCrossAngleMin, fBeamCrossAngleMax);
   }
 
-  // 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);
+  // 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.);

From a44de8bed07f357b25c18dd84c3ac05824976026 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Thu, 2 Jul 2026 16:36:58 +0200
Subject: [PATCH 6/8] Fix minor issues

---
 GeneratorSpectators/GeneratorSpectators.C   |  4 +--
 GeneratorSpectators/GeneratorSpectators.cxx | 28 +++++++++++----------
 2 files changed, 17 insertions(+), 15 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index 6a23869..f200e22 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -108,7 +108,7 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
   return wrap;
 }
 
-FairGenerator *GeneratorProtons(int nNeutrons = -1,
+FairGenerator *GeneratorProtons(int nProtons = -1,
                                 float b = -1.,
                                 bool useFluctuation = false,
                                 float mom = 2680.,
@@ -126,7 +126,7 @@ FairGenerator *GeneratorProtons(int nNeutrons = -1,
   // TODO: implement parameterizations for the number of protons produced
   // At the moment, those for the neutrons are used
   spec->SetNpartFluctuation(useFluctuation);
-  spec->SetNpart(nNeutrons);
+  spec->SetNpart(nProtons);
   spec->SetMomentum(mom);
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
diff --git a/GeneratorSpectators/GeneratorSpectators.cxx b/GeneratorSpectators/GeneratorSpectators.cxx
index 4090b2e..ffbd27d 100644
--- a/GeneratorSpectators/GeneratorSpectators.cxx
+++ b/GeneratorSpectators/GeneratorSpectators.cxx
@@ -91,8 +91,8 @@ 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: ");
@@ -102,13 +102,13 @@ void GeneratorSpectators::GenerateEvent() {
     // Beam divergence and crossing angle
     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. || (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();
@@ -123,8 +123,8 @@ 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);
 
@@ -132,8 +132,8 @@ void GeneratorSpectators::GenerateEvent() {
       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)
@@ -332,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;
   }
 
@@ -348,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;
 }
@@ -373,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)

From 2dd4c82efed2a59ed580b0fbb36bd51fdc794a87 Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Thu, 2 Jul 2026 16:39:44 +0200
Subject: [PATCH 7/8] Make crossing plane configurable

---
 GeneratorSpectators/GeneratorSpectators.C | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.C b/GeneratorSpectators/GeneratorSpectators.C
index f200e22..e162012 100644
--- a/GeneratorSpectators/GeneratorSpectators.C
+++ b/GeneratorSpectators/GeneratorSpectators.C
@@ -46,7 +46,8 @@ FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
                                       float beamDivMax = -999.,
                                       float beamCrossAngle = 0.,
                                       float beamCrossAngleMin = -999.,
-                                      float beamCrossAngleMax = -999.)
+                                      float beamCrossAngleMax = -999.,
+                                      int beamCrossPlane = 2)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -56,7 +57,7 @@ FairGenerator *GeneratorSingleNeutron(float mom = 2680.,
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
   spec->SetSampleDivergence(beamDivMin, beamDivMax);
-  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetCrossing(beamCrossAngle, beamCrossPlane);
   spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
@@ -67,7 +68,8 @@ FairGenerator *GeneratorSingleProton(float mom = 2680.,
                                      float beamDivMax = -999.,
                                      float beamCrossAngle = 0.,
                                      float beamCrossAngleMin = -999.,
-                                     float beamCrossAngleMax = -999.)
+                                     float beamCrossAngleMax = -999.,
+                                     int beamCrossPlane = 2)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -77,7 +79,7 @@ FairGenerator *GeneratorSingleProton(float mom = 2680.,
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
   spec->SetSampleDivergence(beamDivMin, beamDivMax);
-  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetCrossing(beamCrossAngle, beamCrossPlane);
   spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
@@ -91,7 +93,8 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
                                  float beamDivMax = -999.,
                                  float beamCrossAngle = 0.,
                                  float beamCrossAngleMin = -999.,
-                                 float beamCrossAngleMax = -999.)
+                                 float beamCrossAngleMax = -999.,
+                                 int beamCrossPlane = 2)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -103,7 +106,7 @@ FairGenerator *GeneratorNeutrons(int nNeutrons = -1,
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
   spec->SetSampleDivergence(beamDivMin, beamDivMax);
-  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetCrossing(beamCrossAngle, beamCrossPlane);
   spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }
@@ -117,7 +120,8 @@ FairGenerator *GeneratorProtons(int nProtons = -1,
                                 float beamDivMax = -999.,
                                 float beamCrossAngle = 0.,
                                 float beamCrossAngleMin = -999.,
-                                float beamCrossAngleMax = -999.)
+                                float beamCrossAngleMax = -999.,
+                                int beamCrossPlane = 2)
 {
   auto wrap = new GeneratorSpectatorsO2();
   auto spec = wrap->getGenerator();
@@ -131,7 +135,7 @@ FairGenerator *GeneratorProtons(int nProtons = -1,
   spec->SetDirection(0, 0., 0., 1.);
   spec->SetDivergence(beamDiv);
   spec->SetSampleDivergence(beamDivMin, beamDivMax);
-  spec->SetCrossing(beamCrossAngle, 2);
+  spec->SetCrossing(beamCrossAngle, beamCrossPlane);
   spec->SetSampleCrossing(beamCrossAngleMin, beamCrossAngleMax);
   return wrap;
 }

From 398f029a6c7636075c23923562dcf1b0926fe83e Mon Sep 17 00:00:00 2001
From: Fabio Catalano <fabio.catalano@cern.ch>
Date: Fri, 3 Jul 2026 10:46:26 +0200
Subject: [PATCH 8/8] Minor fix

---
 GeneratorSpectators/GeneratorSpectators.cxx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/GeneratorSpectators/GeneratorSpectators.cxx b/GeneratorSpectators/GeneratorSpectators.cxx
index ffbd27d..f61e5a8 100644
--- a/GeneratorSpectators/GeneratorSpectators.cxx
+++ b/GeneratorSpectators/GeneratorSpectators.cxx
@@ -54,7 +54,7 @@ void GeneratorSpectators::Init() {
   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, Fermi flag: %d\n", fPtot, fFermiflag);
+  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);