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    DEFAULT_RIVET_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    
 55    //const FinalState& fsall = apply<FinalState>(event, "FS");
 56    const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 57
 58    const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
 59    const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
 60
 61    double x = dk.x();
 62    double y = dk.y();
 63    const double Q2 = dk.Q2();
 64    
 65     double PT[] = { 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0};
 66    // Extract the particles other than the lepton
 67
 68      if(x<0.01||x>0.1) vetoEvent;
 69      if(y<0.2||y>0.8) vetoEvent;
 70      if(Q2 > 7220 || Q2 < 180) vetoEvent;
 71
 72      _Nevt_after_cuts -> fill();
 73        
 74      Particles particles;
 75      particles.reserve(cfs.particles().size());
 76      
 77      ConstGenParticlePtr dislepGP = dl.out().genParticle();
 78      for (const Particle& p : cfs.particles()) {
 79          ConstGenParticlePtr loopGP = p.genParticle();
 80
 81          if (loopGP == dislepGP) continue;
 82          particles.push_back(p);
 83      }
 84    
 85    const LorentzTransform hcmboost = dk.boostHCM();
 86    for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
 87        const Particle& p = particles[ip1];
 88       
 89       // calculate zh        
 90        double zh = 2.*x/Q2* (dk.beamHadron().E()*p.momentum().E() - dk.beamHadron().pz()*p.momentum().pz()) ;
 91       // cout << " zh " << zh << endl;
 92       // Boost to hcm
 93       
 94       if (zh < 0.2 ) continue ;
 95
 96         const FourMomentum hcmMom = hcmboost.transform(p.momentum());      
 97                  
 98         const double phi =mapAngleMPiToPi(hcmMom.phi())/degree ;
 99         
100         
101         
102// Filling histograms with values of cos(phi) and cos(2phi) wrt the corresponding momentum cuts
103
104         for (size_t i = 0; i < 8; ++i) {
105            if(hcmMom.pT() > PT[i] ) { 
106             _p["cosphi"]->fill(i+1,cos(hcmMom.phi()));
107             _p["cos2phi"]->fill(i+1,cos(2.*hcmMom.phi()));
108             }
109         }
110         
111         
112                    
113         if(hcmMom.pT() > PT[1] ) { _h["A1"] -> fill(phi); }          
114                     
115         if(hcmMom.pT() > PT[3] ) { _h["A2"] -> fill(phi); }
116         
117         if(hcmMom.pT() > PT[5] ) { _h["A3"] -> fill(phi); }
118                   
119         if(hcmMom.pT() > PT[7] ) { _h["A4"] -> fill(phi); }
120        
121     }
122        
123  
124  }
125
126
127
128    /// Normalise histograms etc., after the run
129    void finalize() {
130
131      // correct binwidth in degree to correct for binning from degree to rad by: binwidth/(2PI/10.)
132      double norm = dbl(*_Nevt_after_cuts) ;
133   //   cout << " Nev " << norm << " bin_width= " <<_h["A2"]->bin(0).xWidth() << endl;
134      double degTOrad_width = _h["A1"]->bin(0).xWidth()*10./2./M_PI ;
135      if (norm > 1 ) {
136         scale(_h["A1"], degTOrad_width/norm); 
137         scale(_h["A2"], degTOrad_width/norm); 
138         scale(_h["A3"], degTOrad_width/norm);
139         scale(_h["A4"], degTOrad_width/norm); 
140      }     
141      
142    }
143
144    ///@}
145
146
147    /// @name Histograms
148    ///@{
149    map<string, Histo1DPtr> _h;
150    map<string, Profile1DPtr> _p;
151    map<string, CounterPtr> _c;
152    CounterPtr _Nevt_after_cuts;
153    ///@}
154    
155
156
157  };
158
159
160  DECLARE_RIVET_PLUGIN(ZEUS_2000_I524911);
161
162}