Rivet analyses

title: STAR_2020_I1792394

Measurement of Central Exclusive Production of charged hadron pairs h+h- (h=pi,K,p) at sqrt(s)=200 GeV with forward proton tagging in Roman Pots

Experiment: STAR (RHIC)

Inspire ID: 1792394

Status: VALIDATED

Authors: - Rafal Sikora

References: - arXiv: 2004.11078

Beams: p+ p+

Beam energies: (100.0, 100.0)GeV

Run details: none listed

The differential fiducial cross sections are measured for the process of CEP of h+h- pairs (h=pi,K,p). For the details of analysis and definition of the fiducial region see the arXiv.

Source code:STAR_2020_I1792394.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/FinalState.hh"

namespace Rivet {

/// @brief CEP of h+h- (h=pi,K,p) at sqrt(s)=200 GeV with forward proton tagging class STAR_2020_I1792394 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(STAR_2020_I1792394);


/// @name Analysis methods
/// @{

enum CENTRAL_PARTICLES_PID { _PION, _KAON, _PROTON, _nAllowedPids };
enum PARTICLE_CHARGE { _PLUS, _MINUS, _nSigns };
enum PARTICLE_DIRECTION { _E, _W, _nBeamDirections }; // E = negative p_z, W = positive p_Z

const double minPt[_nAllowedPids] = {0.2 * GeV, 0.3 * GeV, 0.4 * GeV};
const double maxMinPt[_nAllowedPids] = {9e9 * GeV, 0.7 * GeV, 1.1 * GeV};

/// Book histograms and initialise projections before the run
void init() {

  // all final-state particles
  const FinalState fs(Cuts::NOCUT);
  declare(fs, "FS_all");

  // all final-state particles within STAR acceptance for this
  // measurement (reconstructed in the TPC and TOF)
  Cut centralCuts = Cuts::abscharge > 0 && Cuts::abseta < 0.7 && Cuts::pT > 0.2 * GeV
      && (Cuts::abspid == PID::PIPLUS || Cuts::abspid == PID::KPLUS || Cuts::abspid == PID::PROTON);
  const FinalState fs_central(centralCuts);
  declare(fs_central, "FS_central");

  // forward-scattered beam particles detectable in Roman Pots
  // Checking the ID is not needed
  Cut forwardCuts = Cuts::abscharge > 0
      && Cuts::abseta > 5.0; // inclusive cut to select forward particles
  const FinalState fs_forward(forwardCuts);
  declare(fs_forward, "FS_forward");


  // Book histograms with binning taken from HEPdata
  book(_h["m_pipi"], "d01-x01-y01");
  _scaleFactor["m_pipi"] = 1.0;
  book(_h["m_kk"], "d02-x01-y01");
  _scaleFactor["m_kk"] = 1.0;
  book(_h["m_ppbar"], "d03-x01-y01");
  _scaleFactor["m_ppbar"] = 1.0e3;
  book(_h["y_pipi"], "d04-x01-y01");
  _scaleFactor["y_pipi"] = 1.0;
  book(_h["y_kk"], "d05-x01-y01");
  _scaleFactor["y_kk"] = 1.0;
  book(_h["y_ppbar"], "d06-x01-y01");
  _scaleFactor["y_ppbar"] = 1.0e3;
  book(_h["deltaPhi_pipi"], "d07-x01-y01");
  _scaleFactor["deltaPhi_pipi"] = 1.0;
  book(_h["deltaPhi_kk"], "d08-x01-y01");
  _scaleFactor["deltaPhi_kk"] = 1.0e3;
  book(_h["deltaPhi_ppbar"], "d09-x01-y01");
  _scaleFactor["deltaPhi_ppbar"] = 1.0e3;
  book(_h["tSum_pipi"], "d10-x01-y01");
  _scaleFactor["tSum_pipi"] = 1.0;
  book(_h["tSum_kk"], "d11-x01-y01");
  _scaleFactor["tSum_kk"] = 1.0;
  book(_h["tSum_ppbar"], "d12-x01-y01");
  _scaleFactor["tSum_ppbar"] = 1.0e3;
}


/// Perform the per-event analysis
void analyze(const Event& event) {

  // Retrieve all final-state particles
  const FinalState& fs = apply<FinalState>(event, "FS_all");
  // Veto event if number of particles in the final state is different from 4
  if (fs.size() != 4) return;

  // Retrieve accepted centrally produced particles
  const FinalState& fs_central = apply<FinalState>(event, "FS_central");
  // Veto event if number of centrally produced particles is different from 2
  if (fs_central.size() != 2) return;

  // Retrieve forward-scattered particles
  const FinalState& fs_forward = apply<FinalState>(event, "FS_forward");
  // Veto event if number of forward particles is different from 2
  if (fs_forward.size() != 2) return;

  // Continue checking forward particles (intact beam particles)
  // Storing forward particles in an array with cell ID indicating the direction (p_z)
  bool forwardParticlesInFiducialRegion[_nBeamDirections] = {false};
  Particle forwardParticles[_nBeamDirections];
  Vector3 forwardParticle2Vec[_nBeamDirections];
  for (const Particle& p : fs_forward.particles()) {
    const int dir = p.pz() > 0 ? _W : _E;
    forwardParticle2Vec[dir] = Vector3(p.px(), p.py(), 0.0);
    forwardParticles[dir] = p;
    forwardParticlesInFiducialRegion[dir] = p.px() > -0.2 && fabs(p.py()) > 0.2 && fabs(p.py()) < 0.4
        && (pow(p.px() + 0.3, 2) + pow(p.py(), 2)) < 0.25;
  }
  if (!forwardParticlesInFiducialRegion[_E] || !forwardParticlesInFiducialRegion[_W]) return;

  // Storing central particles in an array with cell ID indicating the charge
  Particle csParticles[_nSigns];
  int totalCharge = 0;
  for (const Particle& p : fs_central.particles()) {
    csParticles[p.charge() > 0 ? _PLUS : _MINUS] = p;
    totalCharge += p.charge();
  }
  // Checking the charge conservation, just in case
  if (totalCharge != 0) return;

  // Determine PID of the central pair
  const int pid = (csParticles[_PLUS].pid() == PID::PIPLUS && csParticles[_MINUS].pid() == PID::PIMINUS)
      ? _PION
      : ((csParticles[_PLUS].pid() == PID::KPLUS && csParticles[_MINUS].pid() == PID::KMINUS)
             ? _KAON
             : ((csParticles[_PLUS].pid() == PID::PROTON && csParticles[_MINUS].pid() == PID::ANTIPROTON)
                    ? _PROTON
                    : _nAllowedPids));
  // skip event if particles in a pair are of different ID (should not happen)
  if (pid == _nAllowedPids) return;

  // Checking if central particles pass selection (in principle important for KK and ppbar)
  bool centralParticlesWithinFiducialRegion = true;
  for (int i = 0; i < _nSigns; ++i)
    if (csParticles[i].pT() < minPt[pid]
        || min(csParticles[i].pT(), csParticles[1 - i].pT()) > maxMinPt[pid]) {
      centralParticlesWithinFiducialRegion = false;
      break;
    }
  if (!centralParticlesWithinFiducialRegion) return;


  //-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
  // At this point event satisfies the definition of the fiducial region for events accepted
  // in the CEP measurement at STAR at 200 GeV

  const FourMomentum centralState4Mom = csParticles[_PLUS].momentum() + csParticles[_MINUS].momentum();
  const double invMass = centralState4Mom.mass() / GeV;
  const double rapidity = centralState4Mom.rapidity();
  const double deltaPhi = fabs(forwardParticle2Vec[_W].angle(forwardParticle2Vec[_E])) / degree;

  // We need beam particles to get momentum transfers
  /* // Fragment below did not work for Pythia, unfortunately (four momenta were [0,0,0,0]); using a workaround
  Particle beamParticles[_nBeamDirections];
  const ParticlePair & beams = Rivet::Beam().beams();
  beamParticles[_W] = beams.first;
  beamParticles[_E] = beams.second;
  */
  // workaround - at this point we know that process is exclusive (2 forward protons + 2 central state particles)
  // assume that both beams are of the same type and collision in symmetric (lab frame = c.m.s. frame)
  const double sqrt_s =
      (centralState4Mom + forwardParticles[_E].momentum() + forwardParticles[_W].momentum()).mass();
  const double beamParticleMass = forwardParticles[_W].momentum().mass();
  const double fabsPz = sqrt(sqrt_s * sqrt_s / 4. - beamParticleMass * beamParticleMass);
  FourMomentum beamParticles4Mom[_nBeamDirections];
  beamParticles4Mom[_W] = FourMomentum(sqrt(beamParticleMass * beamParticleMass + fabsPz * fabsPz), 0.,
                                       0., fabsPz);
  beamParticles4Mom[_E] = FourMomentum(sqrt(beamParticleMass * beamParticleMass + fabsPz * fabsPz), 0.,
                                       0., -fabsPz);
  // end of workaround

  double t[_nBeamDirections];
  for (int dir = 0; dir < _nBeamDirections; ++dir)
    t[dir] = (beamParticles4Mom[dir] - forwardParticles[dir].momentum()).mass2() / (GeV * GeV);
  const double tSum = fabs(t[_E] + t[_W]);

  if (pid == _PION) {
    _h["m_pipi"]->fill(invMass);
    _h["y_pipi"]->fill(rapidity);
    _h["deltaPhi_pipi"]->fill(deltaPhi);
    _h["tSum_pipi"]->fill(tSum);
  }
  else if (pid == _KAON) {
    _h["m_kk"]->fill(invMass);
    _h["y_kk"]->fill(rapidity);
    _h["deltaPhi_kk"]->fill(deltaPhi);
    _h["tSum_kk"]->fill(tSum);
  }
  else {
    _h["m_ppbar"]->fill(invMass);
    _h["y_ppbar"]->fill(rapidity);
    _h["deltaPhi_ppbar"]->fill(deltaPhi);
    _h["tSum_ppbar"]->fill(tSum);
  }
}


/// Normalise histograms etc., after the run
void finalize() {

  const double scalingFactor = crossSection() / nanobarn / sumOfWeights();
  // scale to cross section
  for (auto& hist : _h) scale(hist.second, scalingFactor * _scaleFactor[hist.first]);
}

/// @}


/// @name Histograms
/// @{
map<string, Histo1DPtr> _h;
map<string, double>
    _scaleFactor; // map with scale factors to ensure cross section units in agreement with HEPdata
/// @}

};

RIVET_DECLARE_PLUGIN(STAR_2020_I1792394);

} ```