diff --git a/r3bgen/R3BCALIFATestGenerator.cxx b/r3bgen/R3BCALIFATestGenerator.cxx index c24f0121d..9fbdc0750 100644 --- a/r3bgen/R3BCALIFATestGenerator.cxx +++ b/r3bgen/R3BCALIFATestGenerator.cxx @@ -20,6 +20,7 @@ #include #include #include + #include #include #include @@ -33,31 +34,13 @@ R3BCALIFATestGenerator::R3BCALIFATestGenerator() R3BCALIFATestGenerator::R3BCALIFATestGenerator(Int_t pdgid, Int_t mult) : fPDGType(pdgid) , fMult(mult) - , fPDGMass(0) - , fPtMin(0) - , fPtMax(0) - , fEtaMin(0) - , fEtaMax(0) - , fYMin(0) - , fYMax(0) - , fPMin(0) - , fPMax(0) - , fX(0) - , fY(0) - , fZ(0) - , fX1(0) - , fY1(0) - , fZ1(0) - , fX2(0) - , fY2(0) - , fZ2(0) { } Bool_t R3BCALIFATestGenerator::Init() { // Initialize generator - R3BLOG_IF(fatal, fPhiMax - fPhiMin > 360, "phi range is too wide: " << fPhiMin << " 360., "phi range is too wide: " << fPhiMin << " 1, "beta of fragment larger than 1!"); + R3BLOG_IF(fatal, fBetaOfEmittingFragment > 1., "beta of fragment larger than 1!"); + R3BLOG_IF(fatal, + fBeamAngularDivIsSet && !fLorentzBoostIsSet, + "SetBeamAngularDivergence has no effect without SetLorentzBoost"); - for (Int_t i = 0; i < fGammaEnergies.size(); i++) + sumBranchingRatios = 0.; + for (size_t i = 0; i < fGammaEnergies.size(); ++i) { - if (fGammaBranchingRatios[i] > 1) + if (fGammaBranchingRatios[i] > 1.) { R3BLOG(fatal, "gamma branching ratio in position " << i << " larger than 1!"); } sumBranchingRatios += fGammaBranchingRatios[i]; } - R3BLOG_IF(warn, sumBranchingRatios > 1, "gamma branching ratio sum larger than 1!"); + R3BLOG_IF(warn, sumBranchingRatios > 1., "gamma branching ratio sum larger than 1!"); - // Check for particle type TDatabasePDG* pdgBase = TDatabasePDG::Instance(); TParticlePDG* particle = pdgBase->GetParticle(fPDGType); R3BLOG_IF(fatal, !particle, "PDG code " << fPDGType << " not defined."); @@ -88,6 +80,85 @@ Bool_t R3BCALIFATestGenerator::Init() return kTRUE; } +void R3BCALIFATestGenerator::SampleVertex(Double_t& vx, Double_t& vy, Double_t& vz) const +{ + vx = 0.; + vy = 0.; + vz = 0.; + + if (fBoxVtxIsSet) + { + vx = gRandom->Uniform(fX1, fX2); + vy = gRandom->Uniform(fY1, fY2); + vz = gRandom->Uniform(fZ1, fZ2); + return; + } + + if (fPointVtxIsSet) + { + if (fDispersionVtxIsSet) + { + vx = gRandom->Gaus(fX, fDX); + vy = gRandom->Gaus(fY, fDY); + vz = gRandom->Uniform(fZ - fDZ, fZ + fDZ); + } + else + { + vx = fX; + vy = fY; + vz = fZ; + } + } +} + +void R3BCALIFATestGenerator::SampleBeamDirection(Double_t& nx, Double_t& ny, Double_t& nz) const +{ + nx = 0.; + ny = 0.; + nz = 1.; + + if (!fBeamAngularDivIsSet) + return; + + // Sample projected angles in the x-z and y-z planes from a Gaussian + // centered on zero (nominal beam axis) + const Double_t thetaX = gRandom->Gaus(0., fBeamDivX); + const Double_t thetaY = gRandom->Gaus(0., fBeamDivY); + nx = TMath::Sin(thetaX); + ny = TMath::Sin(thetaY); + // Ensure unit vector; clamp argument to avoid sqrt of negative from + // floating-point rounding when divergence angles are very large. + const Double_t nz2 = 1. - nx * nx - ny * ny; + nz = (nz2 > 0.) ? TMath::Sqrt(nz2) : 0.; +} + +void R3BCALIFATestGenerator::ApplyLorentzBoost(Double_t& px, + Double_t& py, + Double_t& pz, + Double_t nx, + Double_t ny, + Double_t nz) const +{ + // General Lorentz boost along unit vector n = (nx, ny, nz): + // E' = gamma * (E + beta * (n . p)) + // p' = p + n * [(gamma - 1) * (n . p) + gamma * beta * E] + // + // fGammaFactor stores sqrt(1 - beta^2) = 1/gamma, matching the + // convention used in SetFragmentVelocity. + + const Double_t E = (fPDGType == 22) ? TMath::Sqrt(px * px + py * py + pz * pz) + : TMath::Sqrt(px * px + py * py + pz * pz + fPDGMass * fPDGMass); + + const Double_t gamma = 1. / fGammaFactor; + const Double_t beta = fBetaOfEmittingFragment; + const Double_t pDotN = px * nx + py * ny + pz * nz; + const Double_t coeff = (gamma - 1.) * pDotN + gamma * beta * E; + + px += nx * coeff; + py += ny * coeff; + pz += nz * coeff; +} + Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) { // Generate one event: produce primary particles emitted from one vertex. @@ -96,11 +167,21 @@ Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) // if SetCosTheta() function is used, the distribution will be uniform in // cos(theta) - Double32_t pabs = 0, phi, pt = 0, theta = 0, eta, y, mt, px, py, pz = 0; + Double32_t pabs = 0., phi = 0., pt = 0., theta = 0., eta = 0., y = 0., mt = 0.; + Double32_t px = 0., py = 0., pz = 0.; + + // Sample the boost direction once for the whole event: all de-excitation + // photons from a single fragment share the same beam direction. + Double_t beamNx = 0., beamNy = 0., beamNz = 1.; + if (fLorentzBoostIsSet) + SampleBeamDirection(beamNx, beamNy, beamNz); // Generate particles - for (Int_t k = 0; k < fMult; k++) + for (Int_t k = 0; k < fMult; ++k) { + Double_t vx = 0., vy = 0., vz = 0.; + SampleVertex(vx, vy, vz); + phi = gRandom->Uniform(fPhiMin, fPhiMax) * TMath::DegToRad(); if (fPRangeIsSet) @@ -118,7 +199,7 @@ Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) else if (fEtaRangeIsSet) { eta = gRandom->Uniform(fEtaMin, fEtaMax); - theta = 2 * TMath::ATan(TMath::Exp(-eta)); + theta = 2. * TMath::ATan(TMath::Exp(-eta)); } else if (fYRangeIsSet) { @@ -135,7 +216,9 @@ Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) pt = pabs * TMath::Sin(theta); } else if (fPtRangeIsSet) + { pz = pt / TMath::Tan(theta); + } } px = pt * TMath::Cos(phi); @@ -151,56 +234,42 @@ Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) if (fNuclearDecayChainIsSet) { LOG_IF(fatal, fPDGType != 22) << "PDG code " << fPDGType << " is not a gamma!"; - for (auto i = 0; i < fGammaEnergies.size(); i++) + for (size_t i = 0; i < fGammaEnergies.size(); ++i) { - auto br = gRandom->Uniform(); + const auto br = gRandom->Uniform(); if (br < fGammaBranchingRatios[i]) { pabs = fGammaEnergies[i]; + phi = gRandom->Uniform(fPhiMin, fPhiMax) * TMath::DegToRad(); theta = acos(gRandom->Uniform(cos(fThetaMin * TMath::DegToRad()), cos(fThetaMax * TMath::DegToRad()))); pz = pabs * TMath::Cos(theta); pt = pabs * TMath::Sin(theta); - px = pt * TMath::Cos(phi); py = pt * TMath::Sin(phi); - auto gammaMomentum = TMath::Sqrt(px * px + py * py + pz * pz); + const auto gammaMomentum = TMath::Sqrt(px * px + py * py + pz * pz); px = px * fGammaEnergies[i] / gammaMomentum; py = py * fGammaEnergies[i] / gammaMomentum; pz = pz * fGammaEnergies[i] / gammaMomentum; - if (sumBranchingRatios <= 1) - break; if (fLorentzBoostIsSet) { - // Lorentz transformation Pz(lab) = gamma * Pz(cm) + gamma * beta * E - // As each Lorentz transformation can be performed sequentially, - // we can separate the gamma factor corresponding to each direction - auto gammaMom = TMath::Sqrt(px * px + py * py + pz * pz); - pz = (pz + fBetaOfEmittingFragment * gammaMom) / fGammaFactor; + ApplyLorentzBoost(px, py, pz, beamNx, beamNy, beamNz); } - primGen->AddTrack(fPDGType, px, py, pz, fX, fY, fZ); + + primGen->AddTrack(fPDGType, px, py, pz, vx, vy, vz); + + if (sumBranchingRatios <= 1.) + break; } } - return kTRUE; + continue; } - if (fPDGType == 22 && fLorentzBoostIsSet) - { // for gamma-rays - // Lorentz transformation Pz(lab) = gamma * Pz(cm) + gamma * beta * E - // As each Lorentz transformation can be performed sequentially, - // we can separate the gamma factor corresponding to each direction - Double32_t gammaMomentum = TMath::Sqrt(px * px + py * py + pz * pz); - pz = (pz + fBetaOfEmittingFragment * gammaMomentum) / fGammaFactor; - } - else if (fLorentzBoostIsSet) - { // for any massive particle - // Lorentz transformation Pz(lab) = gamma * Pz(cm) + gamma * beta * E - // As each Lorentz transformation can be performed sequentially, - // we can separate the gamma factor corresponding to each direction - Double32_t particleEnergy = TMath::Sqrt(px * px + py * py + pz * pz + fPDGMass * fPDGMass); - pz = (pz + fBetaOfEmittingFragment * particleEnergy) / fGammaFactor; + if (fLorentzBoostIsSet) + { + ApplyLorentzBoost(px, py, pz, beamNx, beamNy, beamNz); } if (fDebug) @@ -208,19 +277,27 @@ Bool_t R3BCALIFATestGenerator::ReadEvent(FairPrimaryGenerator* primGen) std::ostringstream oss; oss << "CALIFATestGen: kf=" << fPDGType << ", p=(" << std::fixed << std::setprecision(2) << px << ", " << py << ", " << pz << ") GeV" - << ", x=(" << std::setprecision(1) << fX << ", " << fY << ", " << fZ << ") cm"; + << ", x=(" << std::setprecision(2) << vx << ", " << vy << ", " << vz << ") cm"; + if (fLorentzBoostIsSet && fBeamAngularDivIsSet) + { + oss << ", beam_n=(" << std::setprecision(4) << beamNx << ", " << beamNy << ", " << beamNz << ")"; + } + R3BLOG(info, oss.str()); } - primGen->AddTrack(fPDGType, px, py, pz, fX, fY, fZ); + + primGen->AddTrack(fPDGType, px, py, pz, vx, vy, vz); } + return kTRUE; } void R3BCALIFATestGenerator::SetFragmentVelocity(double beta, double dispersion) { - // Sets the velocity and gamma factor of the fragment emitting the gammas + // Keep original convention: sample once when the setter is called. + // This matches the current class behavior and avoids unexpected API changes. R3BLOG(info, "Set a beta of " << beta << " with a sigma dispersion of " << dispersion); fBetaOfEmittingFragment = gRandom->Gaus(beta, dispersion); - fGammaFactor = TMath::Sqrt(1 - fBetaOfEmittingFragment * fBetaOfEmittingFragment); + fGammaFactor = TMath::Sqrt(1. - fBetaOfEmittingFragment * fBetaOfEmittingFragment); } void R3BCALIFATestGenerator::SetDecayChainPoint(double gammaEnergy, double branchingRatio) @@ -229,4 +306,5 @@ void R3BCALIFATestGenerator::SetDecayChainPoint(double gammaEnergy, double branc fGammaEnergies.push_back(gammaEnergy / 1000.); // GeV fGammaBranchingRatios.push_back(branchingRatio); } + ClassImp(R3BCALIFATestGenerator) diff --git a/r3bgen/R3BCALIFATestGenerator.h b/r3bgen/R3BCALIFATestGenerator.h index 567a9653c..97ec33707 100644 --- a/r3bgen/R3BCALIFATestGenerator.h +++ b/r3bgen/R3BCALIFATestGenerator.h @@ -15,9 +15,24 @@ *@author H. Alvarez Pol * * The R3BCALIFATestGenerator generates gammas with different options - * for testing CALIFA. Copies from FairBoxGenerator (I tried to derive from it, - * but requires a modification of their data members to protected). - * Derived from FairGenerator. + * for testing CALIFA. Originally copied from FairBoxGenerator. + * + * Updated to include vertex distributions: + * - SetXYZ(x,y,z) : central vertex position [cm] + * - SetDxDyDz(sx,sy,sz): Gaussian sigma for x/y, uniform half-width for z + * i.e. z is sampled in [z-sz, z+sz] + * + * For backward compatibility, SetBoxXYZ(...) is retained, which samples a + * uniform box distribution. + * + * Beam angular divergence: + * - SetBeamAngularDivergence(sx_mrad, sy_mrad): a Gaussian + * distribution the fragment beam around the z-axis, characterised by + * 1-sigma projected angles (in mrad) in the x-z and y-z planes. + * The Lorentz boost is applied along the sampled beam direction from this Gaussian + * rather than a fixed z-axis, giving correct lab-frame angles of + * the gammas seen by CALIFA. + * Requires SetLorentzBoost() to be set; has no effect otherwise. **/ #pragma once @@ -39,60 +54,61 @@ class R3BCALIFATestGenerator : public FairGenerator **@param pdgid Particle type (PDG encoding) **@param mult Multiplicity (default is 1) **/ + R3BCALIFATestGenerator(Int_t pdgid, Int_t mult = 1); // Destructor - ~R3BCALIFATestGenerator() = default; + ~R3BCALIFATestGenerator() override = default; // Modifiers - void SetPDGType(Int_t pdg) { fPDGType = pdg; }; - - void SetMultiplicity(Int_t mult) { fMult = mult; }; + void SetPDGType(Int_t pdg) { fPDGType = pdg; } + void SetMultiplicity(Int_t mult) { fMult = mult; } - void SetPRange(double pmin = 0, double pmax = 10) + void SetPRange(double pmin = 0., double pmax = 10.) { fPMin = pmin; fPMax = pmax; fPRangeIsSet = true; } - void SetPtRange(double ptmin = 0, double ptmax = 10) + void SetPtRange(double ptmin = 0., double ptmax = 10.) { fPtMin = ptmin; fPtMax = ptmax; fPtRangeIsSet = true; - }; + } - void SetPhiRange(double phimin = 0, double phimax = 360) + void SetPhiRange(double phimin = 0., double phimax = 360.) { fPhiMin = phimin; fPhiMax = phimax; - }; + } - void SetEtaRange(double etamin = -5, double etamax = 7) + void SetEtaRange(double etamin = -5., double etamax = 7.) { fEtaMin = etamin; fEtaMax = etamax; fEtaRangeIsSet = true; - }; + } - void SetYRange(double ymin = -5, double ymax = 7) + void SetYRange(double ymin = -5., double ymax = 7.) { fYMin = ymin; fYMax = ymax; fYRangeIsSet = true; - }; + } - void SetThetaRange(double thetamin = 0, double thetamax = 180) + void SetThetaRange(double thetamin = 0., double thetamax = 180.) { fThetaMin = thetamin; fThetaMax = thetamax; fThetaRangeIsSet = true; - }; + } - void SetCosTheta() { fCosThetaIsSet = true; }; + void SetCosTheta() { fCosThetaIsSet = true; } - void SetXYZ(Double32_t x = 0, Double32_t y = 0, Double32_t z = 0) + // Vertex configuration: same conventions as R3BINCLRootGenerator + void SetXYZ(Double32_t x = 0., Double32_t y = 0., Double32_t z = 0.) { fX = x; fY = y; @@ -100,12 +116,30 @@ class R3BCALIFATestGenerator : public FairGenerator fPointVtxIsSet = true; } - void SetBoxXYZ(Double32_t x1 = 0, - Double32_t y1 = 0, - Double32_t z1 = 0, - Double32_t x2 = 0, - Double32_t y2 = 0, - Double32_t z2 = 0) + /** + * Set dispersion around the central vertex point. + * - x sampled as Gaussian(x, sx) + * - y sampled as Gaussian(y, sy) + * - z sampled uniformly in [z-sz, z+sz] + */ + void SetDxDyDz(Double32_t sx = 0., Double32_t sy = 0., Double32_t sz = 0.) + { + fDX = sx; + fDY = sy; + fDZ = sz; + fDispersionVtxIsSet = true; + } + + /** + * Backward-compatible box vertex. + * Samples x,y,z uniformly in the specified box. + */ + void SetBoxXYZ(Double32_t x1 = 0., + Double32_t y1 = 0., + Double32_t z1 = 0., + Double32_t x2 = 0., + Double32_t y2 = 0., + Double32_t z2 = 0.) { fX1 = x1; fY1 = y1; @@ -113,7 +147,7 @@ class R3BCALIFATestGenerator : public FairGenerator fX2 = x2; fY2 = y2; fZ2 = z2; - fBoxVtxIsSet = kTRUE; + fBoxVtxIsSet = true; } void SetDebug(bool debug = false) { fDebug = debug; } @@ -122,15 +156,30 @@ class R3BCALIFATestGenerator : public FairGenerator { SetFragmentVelocity(beta); fLorentzBoostIsSet = true; - }; + } void SetFragmentVelocity(double beta = 0., double dispersion = 0.); - void SetNuclearDecayChain() { fNuclearDecayChainIsSet = true; }; + /** + * Set 1-sigma projected angular divergence of the beam around the z-axis. + * @param sigmaX_mrad Gaussian sigma of the x-z projection angle [mrad] + * @param sigmaY_mrad Gaussian sigma of the y-z projection angle [mrad] + * + * Per event a random beam direction is sampled: + * n = (sin(thetaX), sin(thetaY), sqrt(1 - sin²(thetaX) - sin²(thetaY))) + * and the Lorentz boost is applied along n instead of the fixed z-axis. + * Requires SetLorentzBoost() to be set; ignored otherwise. + */ + void SetBeamAngularDivergence(double sigmaX_mrad = 0., double sigmaY_mrad = 0.) + { + fBeamDivX = sigmaX_mrad * 1e-3; // convert mrad -> rad + fBeamDivY = sigmaY_mrad * 1e-3; + fBeamAngularDivIsSet = true; + } - void SetDecayChainPoint(double gammaEnergy = 0, double branchingRatio = 0); + void SetNuclearDecayChain() { fNuclearDecayChainIsSet = true; } + void SetDecayChainPoint(double gammaEnergy = 0., double branchingRatio = 0.); - /** Initializer **/ Bool_t Init() override; /** Creates an event with given type and multiplicity. @@ -139,41 +188,65 @@ class R3BCALIFATestGenerator : public FairGenerator Bool_t ReadEvent(FairPrimaryGenerator* primGen) override; private: - Int_t fPDGType; // Particle type (PDG encoding) - Int_t fMult; // Multiplicity - - double fPDGMass = 0.; // Particle mass [GeV] - double fPtMin = 0., fPtMax = 0.; // Transverse momentum range [GeV] - double fPhiMin = 0., fPhiMax = 360.; // Azimuth angle range [degree] - double fEtaMin = 0., fEtaMax = 0.; // Pseudorapidity range in lab system - double fYMin = 0., fYMax = 0.; // Rapidity range in lab system - double fPMin = 0., fPMax = 0.; // Momentum range in lab system - double fThetaMin = 0., fThetaMax = 180.; // Polar angle range in lab system [degree] - double fX = 0., fY = 0., fZ = 0.; // Point vertex coordinates [cm] - double fX1 = 0., fY1 = 0., fZ1 = 0., fX2 = 0., fY2 = 0., fZ2 = 0.; // Box vertex coords (x1,y1,z1)->(x2,y2,z2) - double sumBranchingRatios = 0.; - - bool fEtaRangeIsSet = false; // True if eta range is set - bool fYRangeIsSet = false; // True if rapidity range is set - bool fThetaRangeIsSet = true; // True if theta range is set - bool fCosThetaIsSet = false; // True if uniform distribution in - // cos(theta) is set (default -> not set) - bool fPtRangeIsSet = false; // True if transverse momentum range is set - bool fPRangeIsSet = false; // True if abs.momentum range is set - bool fPointVtxIsSet = false; // True if point vertex is set - bool fBoxVtxIsSet = false; // True if box vertex is set - bool fDebug = false; // Debug switch + void SampleVertex(Double_t& vx, Double_t& vy, Double_t& vz) const; + + /** Sample a beam direction unit vector for the current event. + * When fBeamAngularDivIsSet, spreads the beam by Gaussian projected angles. + * Otherwise returns the nominal z-axis (0, 0, 1). */ + void SampleBeamDirection(Double_t& nx, Double_t& ny, Double_t& nz) const; + + /** Apply a Lorentz boost along the unit vector (nx, ny, nz) to the + * momentum (px, py, pz), using fBetaOfEmittingFragment / fGammaFactor. + * Handles both massive particles and photons (fPDGType == 22). */ + void ApplyLorentzBoost(Double_t& px, Double_t& py, Double_t& pz, Double_t nx, Double_t ny, Double_t nz) const; + + Int_t fPDGType = 22; + Int_t fMult = 1; + + double fPDGMass = 0.; + double fPtMin = 0., fPtMax = 0.; + double fPhiMin = 0., fPhiMax = 360.; + double fEtaMin = 0., fEtaMax = 0.; + double fYMin = 0., fYMax = 0.; + double fPMin = 0., fPMax = 0.; + double fThetaMin = 0., fThetaMax = 180.; + + // Central vertex coordinates [cm] + double fX = 0., fY = 0., fZ = 0.; + // INCL-like dispersions: x/y Gaussian sigma, z uniform half-width [cm] + double fDX = 0., fDY = 0., fDZ = 0.; + // Explicit box vertex coordinates [cm] + double fX1 = 0., fY1 = 0., fZ1 = 0., fX2 = 0., fY2 = 0., fZ2 = 0.; - // SPECIFIC OF CALIFA - std::vector fGammaEnergies; // Gamma energies for the nuclear decay chain in GeV UNITS - std::vector fGammaBranchingRatios; // Gamma branching ratios for the nuclear decay chain + double sumBranchingRatios = 0.; - double fBetaOfEmittingFragment = 0.; // Velocity of the fragment emitting the gammas - double fGammaFactor = 1.; // Velocity of the fragment emitting the gammas + bool fEtaRangeIsSet = false; // True if eta range is set + bool fYRangeIsSet = false; // True if rapidity range is set + bool fThetaRangeIsSet = true; // True if theta range is set + bool fCosThetaIsSet = false; // True if uniform distribution in + // cos(theta) is set (default -> not set) + bool fPtRangeIsSet = false; // True if transverse momentum range is set + bool fPRangeIsSet = false; // True if abs.momentum range is set + bool fPointVtxIsSet = true; // True if point vertex is set + bool fDispersionVtxIsSet = true; // True if x,y,z dispersion around vertex is set + bool fBoxVtxIsSet = false; // True if box vertex is set + bool fDebug = false; // Debug switch + + // CALIFA-specific + std::vector fGammaEnergies; // gamma energies in GeV + std::vector fGammaBranchingRatios; // branching ratios + + double fBetaOfEmittingFragment = 0.; // Velocity of the fragment emitting the photons + double fGammaFactor = 1.; // Lorentz gamma of the fragment emitting the photons bool fLorentzBoostIsSet = false; // True if Lorentz Boost is set bool fNuclearDecayChainIsSet = false; // True if a nuclear decay chain is set + // Beam angular divergence (projected 1-sigma angles in x-z and y-z planes) + double fBeamDivX = 0.; // [rad] + double fBeamDivY = 0.; // [rad] + bool fBeamAngularDivIsSet = false; + public: - ClassDefOverride(R3BCALIFATestGenerator, 2); // NOLINT + ClassDefOverride(R3BCALIFATestGenerator, 4); // NOLINT };