Rivet analyses

// -*- 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);

}