Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/DISLepton.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {


  /// @brief Measurement of azimuthal asymmetries in deep inelastic scattering (ZEUS)
  class ZEUS_2000_I524911 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2000_I524911);


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

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

      // Initialise and register projections
      declare(DISLepton(), "Lepton");
      declare(DISKinematics(), "Kinematics");

      // The basic final-state projection:
      // all final-state particles within
      // the given eta acceptance

      const ChargedFinalState cfs;
      declare(cfs, "CFS");

      book(_h["A1"], 1, 1, 1);
      book(_h["A2"], 1, 1, 2);
      book(_h["A3"], 1, 1, 3);
      book(_h["A4"], 1, 1, 4);

      book(_p["cosphi"],2, 1, 1) ;
      book(_p["cos2phi"],2, 1, 2) ;

      // counter pointer to store the no. of events
      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");



    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if(_edges.empty()) _edges = _p["cosphi"]->xEdges();
    //const FinalState& fsall = apply<FinalState>(event, "FS");
    const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");

    const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
    const DISLepton& dl = apply<DISLepton>(event,"Lepton");

    double x = dk.x();
    double y = dk.y();
    const double Q2 = dk.Q2();

    // Extract the particles other than the lepton

      if(x<0.01||x>0.1) vetoEvent;
      if(y<0.2||y>0.8) vetoEvent;
      if(Q2 > 7220 || Q2 < 180) vetoEvent;

      _Nevt_after_cuts -> fill();

      Particles particles;
      particles.reserve(cfs.particles().size());

      ConstGenParticlePtr dislepGP = dl.out().genParticle();
      for (const Particle& p : cfs.particles()) {
          ConstGenParticlePtr loopGP = p.genParticle();

          if (loopGP == dislepGP) continue;
          particles.push_back(p);
      }

    const LorentzTransform hcmboost = dk.boostHCM();
    for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
        const Particle& p = particles[ip1];

       // calculate zh
        double zh = 2.*x/Q2* (dk.beamHadron().E()*p.momentum().E() - dk.beamHadron().pz()*p.momentum().pz()) ;
       // cout << " zh " << zh << endl;
       // Boost to hcm

       if (zh < 0.2 ) continue ;

         const FourMomentum hcmMom = hcmboost.transform(p.momentum());

         const double phi_rad = mapAngleMPiToPi(hcmMom.phi());
         const double phi_deg = phi_rad/degree;

        // Filling histograms with values of cos(phi) and cos(2phi)
        // wrt the corresponding momentum cuts
         for (size_t i = 0; i < 8; ++i) {
            if(hcmMom.pT() > _ptCut[i] ) {
              _p["cosphi"] ->fill(_edges[i], cos(phi_rad));
              _p["cos2phi"]->fill(_edges[i], cos(2.*phi_rad));
            }
         }

         if(hcmMom.pT() > _ptCut[1] ) { _h["A1"] -> fill(phi_deg); }

         if(hcmMom.pT() > _ptCut[3] ) { _h["A2"] -> fill(phi_deg); }

         if(hcmMom.pT() > _ptCut[5] ) { _h["A3"] -> fill(phi_deg); }

         if(hcmMom.pT() > _ptCut[7] ) { _h["A4"] -> fill(phi_deg); }

     }


  }



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

      // correct binwidth in degree to correct for binning from degree to rad by: binwidth/(2PI/10.)
      double norm = dbl(*_Nevt_after_cuts) ;
      double degTOrad_width = _h["A1"]->bin(1).xWidth()*10./2./M_PI ;
      if (norm > 1 ) {
         scale(_h["A1"], degTOrad_width/norm);
         scale(_h["A2"], degTOrad_width/norm);
         scale(_h["A3"], degTOrad_width/norm);
         scale(_h["A4"], degTOrad_width/norm);
      }

    }

    ///@}


    /// @name Histograms
    ///@{
    map<string, Histo1DPtr> _h;
    map<string, BinnedProfilePtr<string>> _p;
    map<string, CounterPtr> _c;
    CounterPtr _Nevt_after_cuts;
    vector<string> _edges;
    vector<double> _ptCut = { 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0};
    ///@}



  };


  RIVET_DECLARE_PLUGIN(ZEUS_2000_I524911);

}