1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/DISLepton.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Measurement of azimuthal asymmetries in deep inelastic scattering (ZEUS)
 12  class ZEUS_2000_I524911 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2000_I524911);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26      declare(DISLepton(), "Lepton");
 27      declare(DISKinematics(), "Kinematics");
 28
 29      // The basic final-state projection:
 30      // all final-state particles within
 31      // the given eta acceptance
 32
 33      const ChargedFinalState cfs;
 34      declare(cfs, "CFS");
 35
 36      book(_h["A1"], 1, 1, 1);
 37      book(_h["A2"], 1, 1, 2);
 38      book(_h["A3"], 1, 1, 3);
 39      book(_h["A4"], 1, 1, 4);
 40
 41      book(_p["cosphi"],2, 1, 1) ;
 42      book(_p["cos2phi"],2, 1, 2) ;
 43
 44      // counter pointer to store the no. of events
 45      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
 46
 47
 48
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      if(_edges.empty()) _edges = _p["cosphi"]->xEdges();
 55    //const FinalState& fsall = apply<FinalState>(event, "FS");
 56    const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 57
 58    const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 59    const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 60
 61    double x = dk.x();
 62    double y = dk.y();
 63    const double Q2 = dk.Q2();
 64
 65    // Extract the particles other than the lepton
 66
 67      if(x<0.01||x>0.1) vetoEvent;
 68      if(y<0.2||y>0.8) vetoEvent;
 69      if(Q2 > 7220 || Q2 < 180) vetoEvent;
 70
 71      _Nevt_after_cuts -> fill();
 72
 73      Particles particles;
 74      particles.reserve(cfs.particles().size());
 75
 76      ConstGenParticlePtr dislepGP = dl.out().genParticle();
 77      for (const Particle& p : cfs.particles()) {
 78          ConstGenParticlePtr loopGP = p.genParticle();
 79
 80          if (loopGP == dislepGP) continue;
 81          particles.push_back(p);
 82      }
 83
 84    const LorentzTransform hcmboost = dk.boostHCM();
 85    for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
 86        const Particle& p = particles[ip1];
 87
 88       // calculate zh
 89        double zh = 2.*x/Q2* (dk.beamHadron().E()*p.momentum().E() - dk.beamHadron().pz()*p.momentum().pz()) ;
 90       // cout << " zh " << zh << endl;
 91       // Boost to hcm
 92
 93       if (zh < 0.2 ) continue ;
 94
 95         const FourMomentum hcmMom = hcmboost.transform(p.momentum());
 96
 97         const double phi_rad = mapAngleMPiToPi(hcmMom.phi());
 98         const double phi_deg = phi_rad/degree;
 99
100        // Filling histograms with values of cos(phi) and cos(2phi)
101        // wrt the corresponding momentum cuts
102         for (size_t i = 0; i < 8; ++i) {
103            if(hcmMom.pT() > _ptCut[i] ) {
104              _p["cosphi"] ->fill(_edges[i], cos(phi_rad));
105              _p["cos2phi"]->fill(_edges[i], cos(2.*phi_rad));
106            }
107         }
108
109         if(hcmMom.pT() > _ptCut[1] ) { _h["A1"] -> fill(phi_deg); }
110
111         if(hcmMom.pT() > _ptCut[3] ) { _h["A2"] -> fill(phi_deg); }
112
113         if(hcmMom.pT() > _ptCut[5] ) { _h["A3"] -> fill(phi_deg); }
114
115         if(hcmMom.pT() > _ptCut[7] ) { _h["A4"] -> fill(phi_deg); }
116
117     }
118
119
120  }
121
122
123
124    /// Normalise histograms etc., after the run
125    void finalize() {
126
127      // correct binwidth in degree to correct for binning from degree to rad by: binwidth/(2PI/10.)
128      double norm = dbl(*_Nevt_after_cuts) ;
129      double degTOrad_width = _h["A1"]->bin(1).xWidth()*10./2./M_PI ;
130      if (norm > 1 ) {
131         scale(_h["A1"], degTOrad_width/norm);
132         scale(_h["A2"], degTOrad_width/norm);
133         scale(_h["A3"], degTOrad_width/norm);
134         scale(_h["A4"], degTOrad_width/norm);
135      }
136
137    }
138
139    ///@}
140
141
142    /// @name Histograms
143    ///@{
144    map<string, Histo1DPtr> _h;
145    map<string, BinnedProfilePtr<string>> _p;
146    map<string, CounterPtr> _c;
147    CounterPtr _Nevt_after_cuts;
148    vector<string> _edges;
149    vector<double> _ptCut = { 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0};
150    ///@}
151
152
153
154  };
155
156
157  RIVET_DECLARE_PLUGIN(ZEUS_2000_I524911);
158
159}