1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief CEP of h+h- (h=pi,K,p) at sqrt(s)=200 GeV with forward proton tagging
 10  class STAR_2020_I1792394 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(STAR_2020_I1792394);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    enum CENTRAL_PARTICLES_PID {_PION, _KAON, _PROTON, _nAllowedPids};
 21    enum PARTICLE_CHARGE { _PLUS, _MINUS, _nSigns };
 22    enum PARTICLE_DIRECTION { _E, _W, _nBeamDirections }; // E = negative p_z, W = positive p_Z
 23
 24    const double minPt[_nAllowedPids] = { 0.2*GeV, 0.3*GeV, 0.4*GeV };
 25    const double maxMinPt[_nAllowedPids] = { 9e9*GeV, 0.7*GeV, 1.1*GeV };
 26
 27    /// Book histograms and initialise projections before the run
 28    void init() {
 29
 30      // all final-state particles
 31      const FinalState fs( Cuts::NOCUT );
 32      declare(fs, "FS_all");
 33
 34      // all final-state particles within STAR acceptance for this
 35      // measurement (reconstructed in the TPC and TOF)
 36      Cut centralCuts = Cuts::abscharge > 0
 37                     && Cuts::abseta < 0.7
 38                     && Cuts::pT > 0.2*GeV
 39                     && (  Cuts::abspid == PID::PIPLUS
 40                        || Cuts::abspid == PID::KPLUS
 41                        || Cuts::abspid == PID::PROTON );
 42      const FinalState fs_central( centralCuts );
 43      declare(fs_central, "FS_central");
 44
 45      // forward-scattered beam particles detectable in Roman Pots
 46      // Checking the ID is not needed
 47      Cut forwardCuts = Cuts::abscharge > 0
 48                     && Cuts::abseta > 5.0; // inclusive cut to select forward particles
 49      const FinalState fs_forward(forwardCuts);
 50      declare(fs_forward, "FS_forward");
 51
 52
 53      // Book histograms with binning taken from HEPdata
 54      book(_h["m_pipi"], "d01-x01-y01");            _scaleFactor["m_pipi"] = 1.0;
 55      book(_h["m_kk"], "d02-x01-y01");              _scaleFactor["m_kk"] = 1.0;
 56      book(_h["m_ppbar"], "d03-x01-y01");           _scaleFactor["m_ppbar"] = 1.0e3;
 57      book(_h["y_pipi"], "d04-x01-y01");            _scaleFactor["y_pipi"] = 1.0;
 58      book(_h["y_kk"], "d05-x01-y01");              _scaleFactor["y_kk"] = 1.0;
 59      book(_h["y_ppbar"], "d06-x01-y01");           _scaleFactor["y_ppbar"] = 1.0e3;
 60      book(_h["deltaPhi_pipi"], "d07-x01-y01");     _scaleFactor["deltaPhi_pipi"] = 1.0;
 61      book(_h["deltaPhi_kk"], "d08-x01-y01");       _scaleFactor["deltaPhi_kk"] = 1.0e3;
 62      book(_h["deltaPhi_ppbar"], "d09-x01-y01");    _scaleFactor["deltaPhi_ppbar"] = 1.0e3;
 63      book(_h["tSum_pipi"], "d10-x01-y01");         _scaleFactor["tSum_pipi"] = 1.0;
 64      book(_h["tSum_kk"], "d11-x01-y01");           _scaleFactor["tSum_kk"] = 1.0;
 65      book(_h["tSum_ppbar"], "d12-x01-y01");        _scaleFactor["tSum_ppbar"] = 1.0e3;
 66
 67    }
 68
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72
 73      // Retrieve all final-state particles
 74      const FinalState & fs = apply<FinalState>(event, "FS_all");
 75      // Veto event if number of particles in the final state is different from 4
 76      if(fs.size() != 4)
 77          return;
 78
 79      // Retrieve accepted centrally produced particles
 80      const FinalState & fs_central = apply<FinalState>(event, "FS_central");
 81      // Veto event if number of centrally produced particles is different from 2
 82      if(fs_central.size() != 2)
 83          return;
 84
 85      // Retrieve forward-scattered particles
 86      const FinalState & fs_forward = apply<FinalState>(event, "FS_forward");
 87      // Veto event if number of forward particles is different from 2
 88      if(fs_forward.size() != 2)
 89          return;
 90
 91      // Continue checking forward particles (intact beam particles)
 92      // Storing forward particles in an array with cell ID indicating the direction (p_z)
 93      bool forwardParticlesInFiducialRegion[_nBeamDirections] = {false};
 94      Particle forwardParticles[_nBeamDirections];
 95      Vector3 forwardParticle2Vec[_nBeamDirections];
 96      for(const Particle & p : fs_forward.particles()){
 97          const int dir = p.pz() > 0 ? _W : _E;
 98          forwardParticle2Vec[dir] = Vector3(p.px(), p.py(), 0.0);
 99          forwardParticles[dir] = p;
100          forwardParticlesInFiducialRegion[dir] = p.px() > -0.2
101                                                && fabs(p.py()) > 0.2
102                                                && fabs(p.py()) < 0.4
103                                                && (pow( p.px() + 0.3, 2) + pow( p.py(), 2 )) < 0.25;
104      }
105      if( !forwardParticlesInFiducialRegion[_E] || !forwardParticlesInFiducialRegion[_W] )
106          return;
107
108      // Storing central particles in an array with cell ID indicating the charge
109      Particle csParticles[_nSigns];
110      int totalCharge = 0;
111      for(const Particle & p : fs_central.particles()){
112          csParticles[ p.charge()>0 ? _PLUS : _MINUS] = p;
113          totalCharge += p.charge();
114      }
115      // Checking the charge conservation, just in case
116      if( totalCharge != 0 )
117          return;
118
119      // Determine PID of the central pair
120      const int pid = ( csParticles[_PLUS].pid()==PID::PIPLUS && csParticles[_MINUS].pid()==PID::PIMINUS ) ? _PION :
121                     (( csParticles[_PLUS].pid()==PID::KPLUS  && csParticles[_MINUS].pid()==PID::KMINUS) ? _KAON :
122                     (( csParticles[_PLUS].pid()==PID::PROTON && csParticles[_MINUS].pid()==PID::ANTIPROTON ) ? _PROTON : _nAllowedPids) );
123      // skip event if particles in a pair are of different ID (should not happen)
124      if(pid == _nAllowedPids)
125          return;
126
127      // Checking if central particles pass selection (in principle important for KK and ppbar)
128      bool centralParticlesWithinFiducialRegion = true;
129      for(int i=0; i<_nSigns; ++i)
130          if( csParticles[i].pT() < minPt[pid] || min(csParticles[i].pT(), csParticles[1-i].pT()) > maxMinPt[pid] ){
131              centralParticlesWithinFiducialRegion = false;
132              break;
133          }
134      if( !centralParticlesWithinFiducialRegion )
135          return;
136
137
138      //-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
139      // At this point event satisfies the definition of the fiducial region for events accepted
140      // in the CEP measurement at STAR at 200 GeV
141
142      const FourMomentum centralState4Mom = csParticles[_PLUS].momentum() + csParticles[_MINUS].momentum();
143      const double invMass = centralState4Mom.mass()/GeV;
144      const double rapidity = centralState4Mom.rapidity();
145      const double deltaPhi = fabs( forwardParticle2Vec[_W].angle( forwardParticle2Vec[_E] ) )/degree;
146
147      // We need beam particles to get momentum transfers
148      /* // Fragment below did not work for Pythia, unfortunately (four momenta were [0,0,0,0]); using a workaround
149      Particle beamParticles[_nBeamDirections];
150      const ParticlePair & beams = Rivet::Beam().beams();
151      beamParticles[_W] = beams.first;
152      beamParticles[_E] = beams.second;
153      */
154      // workaround - at this point we know that process is exclusive (2 forward protons + 2 central state particles)
155      // assume that both beams are of the same type and collision in symmetric (lab frame = c.m.s. frame)
156      const double sqrt_s = (centralState4Mom + forwardParticles[_E].momentum() + forwardParticles[_W].momentum()).mass();
157      const double beamParticleMass = forwardParticles[_W].momentum().mass();
158      const double fabsPz = sqrt( sqrt_s*sqrt_s/4. - beamParticleMass*beamParticleMass );
159      FourMomentum beamParticles4Mom[_nBeamDirections];
160      beamParticles4Mom[_W] = FourMomentum( sqrt(beamParticleMass*beamParticleMass + fabsPz*fabsPz), 0., 0., fabsPz );
161      beamParticles4Mom[_E] = FourMomentum( sqrt(beamParticleMass*beamParticleMass + fabsPz*fabsPz), 0., 0., -fabsPz );
162      // end of workaround
163
164      double t[_nBeamDirections];
165      for(int dir=0; dir<_nBeamDirections; ++dir)
166          t[dir] = (beamParticles4Mom[dir] - forwardParticles[dir].momentum()).mass2()/(GeV*GeV);
167      const double tSum = fabs(t[_E] + t[_W]);
168
169      if( pid==_PION ){
170        _h["m_pipi"]->fill( invMass );
171        _h["y_pipi"]->fill( rapidity );
172        _h["deltaPhi_pipi"]->fill( deltaPhi );
173        _h["tSum_pipi"]->fill( tSum );
174      } else if( pid==_KAON ){
175        _h["m_kk"]->fill( invMass );
176        _h["y_kk"]->fill( rapidity );
177        _h["deltaPhi_kk"]->fill( deltaPhi );
178        _h["tSum_kk"]->fill( tSum );
179      } else{
180        _h["m_ppbar"]->fill( invMass );
181        _h["y_ppbar"]->fill( rapidity );
182        _h["deltaPhi_ppbar"]->fill( deltaPhi );
183        _h["tSum_ppbar"]->fill( tSum );
184      }
185
186
187    }
188
189
190    /// Normalise histograms etc., after the run
191    void finalize() {
192
193      const double scalingFactor = crossSection()/nanobarn/sumOfWeights();
194      // scale to cross section
195      for(auto &hist : _h)
196        scale(hist.second, scalingFactor*_scaleFactor[hist.first]);
197    }
198
199    /// @}
200
201
202    /// @name Histograms
203    /// @{
204    map<string, Histo1DPtr> _h;
205    map<string, double> _scaleFactor; // map with scale factors to ensure cross section units in agreement with HEPdata
206    /// @}
207
208
209  };
210
211
212  RIVET_DECLARE_PLUGIN(STAR_2020_I1792394);
213
214}