1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/DISLepton.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief A Study of the Fragmentation of Quarks in ep Collisions at HERA (H1)
 12  class H1_1995_I394793 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1995_I394793);
 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      const FinalState fs(Cuts::abseta < 4.9);
 33      declare(fs, "FS");
 34      const ChargedFinalState cfs;
 35      declare(cfs, "CFS");
 36      
 37      // book a counter
 38      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
 39      book(_Nevt_afterfwd_cuts, "TMP/Nevt_afterfwd_cuts");
 40      book(_Nevt_afterh_cuts, "TMP/Nevt_afterh_cuts");
 41      book(_Nevt_afterhfwd_cuts, "TMP/Nevt_afterhfwd_cuts");
 42
 43
 44      // Book histograms
 45      // specify custom binning
 46      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 47      book(_d["costh_lowQ"]        , 1, 1, 1);
 48      book(_d["costh_highQ"]       , 1, 1, 2);
 49      book(_d["costh_lowQ_noEfwd"] , 2, 1, 1);
 50      book(_d["costh_highQ_noEfwd"], 2, 1, 2);
 51      book(_h["xp_posCharge_lowQ"], 3, 1, 1);
 52      book(_h["xp_negCharge_lowQ"], 3, 1, 2);
 53      book(_h["xp_posCharge_highQ"], 3, 1, 3);
 54      book(_h["xp_negCharge_highQ"], 3, 1, 4);
 55      book(_h["ksi_lowQ"], 4, 1, 1);
 56      book(_h["ksi_highQ"], 4, 1, 2);
 57      book(_s["Mult_vrs_Q2"], 5, 1, 1);
 58      book(_s["Mult_vrs_Q2_noEfwd"], 6, 1, 1);
 59      
 60      book(_h["Mult_vrs_Q2_nchrg"],"TMP/Mult_vrs_Q2_nchrg", refData(5,1,1));
 61      book(_h["Mult_vrs_Q2_noEfwd_nchrg"],"TMP/Mult_vrs_Q2_noEfwd_nchrg", refData(6,1,1));
 62      book(_h["Mult_vrs_Q2_count"],"TMP/Mult_vrs_Q2_count", refData(5,1,1));
 63      book(_h["Mult_vrs_Q2_noEfwd_count"],"TMP/Mult_vrs_Q2_noEfwd_count", refData(6,1,1));
 64
 65    }
 66 
 67
 68    /// Perform the per-event analysis
 69    void analyze(const Event& event) {
 70      if(_edges.empty()) _edges = _d["costh_lowQ"]->xEdges();
 71
 72      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 73    
 74   
 75      //DIS kinematics
 76      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 77      if ( dk.failed() ) vetoEvent;
 78      double y   = dk.y();
 79      double w2  = dk.W2();
 80      double Q2  = dk.Q2();
 81      
 82      
 83      
 84      bool cut = Q2 > 12 && y < 0.6 && w2 > 3000 ;
 85      
 86      
 87      if ( !cut ) vetoEvent ; 
 88      
 89      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 90      if ( dl.failed() ) vetoEvent;
 91      /*
 92      cout << "  scattered lepton angle " << 180.- dl.out().momentum().angle(dl.in().momentum())/degree << endl;
 93      cout << " in  lepton " << dl.in().momentum() << endl;
 94      cout << " out lepton " << dl.out().momentum() << endl;
 95      */
 96      
 97      const FinalState& fs = apply<FinalState>(event, "FS");
 98      Particles particles; particles.reserve(fs.size());
 99      ConstGenParticlePtr dislepGP = dl.out().genParticle();
100      for (const Particle& p : cfs.particles()) {
101        ConstGenParticlePtr loopGP = p.genParticle();
102        if (loopGP == dislepGP) continue;
103        particles.push_back(p);
104      }
105
106
107      double efwd = 0.;
108
109      for (const Particle& p : particles) {
110        const double th = 180. - p.momentum().angle(dl.in().momentum())/degree;
111
112        if (inRange(th, 4.4, 15.0)) { 
113           efwd += p.E();
114           //cout << " angle " << th << " pid " << p.pid() << " Efwd = " << efwd << endl;
115           } 
116      }
117
118      bool evcut[2];
119      evcut[0] =   efwd > 0.5;
120      
121       // fill the counter
122      _Nevt_after_cuts -> fill();
123      if (Q2 > 100 )  _Nevt_afterh_cuts -> fill();
124       
125      if ( evcut[0] && (Q2 < 80 ) ) _Nevt_afterfwd_cuts -> fill();
126      if ( evcut[0] && (Q2 > 100 ) ) _Nevt_afterhfwd_cuts -> fill();
127
128      double n_charg = 0; 
129
130      // Boost to Breit
131      const LorentzTransform breitboost = dk.boostBreit();
132       
133      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
134        const Particle& p = particles[ip1];
135        const FourMomentum BreMom = breitboost.transform(p.momentum());
136       // cout << BreMom.pz() << endl;
137        double x = cos(BreMom.theta());
138
139        size_t idx = _axis.index(x);
140        string edge = "OTHER";
141        if(idx && idx <= _edges.size()) edge=_edges[idx-1];
142        
143        if (Q2 < 80 )  _d["costh_lowQ" ] ->fill(edge);
144        if (Q2 > 100 ) _d["costh_highQ"]  ->fill(edge);
145        
146        if (Q2 < 80  && evcut[0]) _d["costh_lowQ_noEfwd" ] ->fill(edge);
147        if (Q2 > 100 && evcut[0]) _d["costh_highQ_noEfwd"]  ->fill(edge);
148
149
150        if ( BreMom.pz() > 0. ) continue;
151        double pcal= sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2()) ;
152        double xp = 2*pcal/(sqrt(Q2));
153        double xi = log(1/xp);
154        
155        double charge = p.charge() ;
156       // cout << " charge " << charge << endl;
157       
158        if (charge > 0 ) {
159          if (Q2 < 80 ) _h["xp_posCharge_lowQ"] -> fill(xp);
160          if (Q2 > 100 ) _h["xp_posCharge_highQ"] -> fill(xp);
161        } else {
162         if (Q2 < 80 ) _h["xp_negCharge_lowQ"] -> fill(xp);
163         if (Q2 > 100 ) _h["xp_negCharge_highQ"] -> fill(xp);
164        }
165    
166        if (Q2 < 80 ) _h["ksi_lowQ"] -> fill(xi);
167        if (Q2 > 100 ) _h["ksi_highQ"] -> fill (xi);
168 
169        n_charg = n_charg + 1; 
170
171     }
172     _h["Mult_vrs_Q2_nchrg"] -> fill(Q2,n_charg) ;
173     _h["Mult_vrs_Q2_count"] -> fill(Q2) ;
174     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_nchrg"] -> fill(Q2,n_charg) ;
175     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_count"] -> fill(Q2) ;
176}
177
178    /// Normalise histograms etc., after the run
179    void finalize() {
180     
181   
182      normalize(_h["xp_posCharge_lowQ"]);
183      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_posCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
184
185      normalize(_h["xp_negCharge_lowQ"]);
186      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_negCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
187
188      scale(_d["costh_lowQ"], 10.0/ *_Nevt_after_cuts);
189      if(dbl(*_Nevt_afterh_cuts)>0) scale(_d["costh_highQ"], 10.0/ *_Nevt_afterh_cuts);
190      
191      if(dbl(*_Nevt_afterfwd_cuts)>0) scale(_d["costh_lowQ_noEfwd"], 10.0/ *_Nevt_afterfwd_cuts);
192      if(dbl(*_Nevt_afterhfwd_cuts)>0) scale(_d["costh_highQ_noEfwd"], 10.0/ *_Nevt_afterhfwd_cuts);
193      
194      scale(_h["ksi_lowQ"], 1.0/ *_Nevt_after_cuts);
195      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["ksi_highQ"], 1.0/ *_Nevt_afterh_cuts);
196
197      divide(_h["Mult_vrs_Q2_nchrg"], _h["Mult_vrs_Q2_count"], _s["Mult_vrs_Q2"]);
198      divide(_h["Mult_vrs_Q2_noEfwd_nchrg"], _h["Mult_vrs_Q2_noEfwd_count"], _s["Mult_vrs_Q2_noEfwd"]);
199    }
200
201    ///@}
202
203
204    /// @name Histograms
205    ///@{
206    map<string, Histo1DPtr> _h;
207    map<string, BinnedHistoPtr<string> > _d;
208    map<string, Profile1DPtr> _p;
209    map<string, CounterPtr> _c;
210    map<string, Estimate1DPtr> _s;
211    CounterPtr _Nevt_after_cuts;
212    CounterPtr _Nevt_afterfwd_cuts;
213    CounterPtr _Nevt_afterh_cuts;
214    CounterPtr _Nevt_afterhfwd_cuts;
215    vector<string> _edges;
216    YODA::Axis<double> _axis = YODA::Axis<double>{-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,
217                                                  0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
218    ///@}
219
220
221  };
222
223
224  RIVET_DECLARE_PLUGIN(H1_1995_I394793);
225
226}