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(_h["costh_lowQ"], 1, 1, 1);
 48      book(_h["costh_highQ"], 1, 1, 2);     
 49      book(_h["costh_lowQ_noEfwd"], 2, 1, 1);
 50      book(_h["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
 71     const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 72    
 73   
 74    //DIS kinematics
 75      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 76      if ( dk.failed() ) vetoEvent;
 77      double y   = dk.y();
 78      double w2  = dk.W2();
 79      double Q2  = dk.Q2();
 80      
 81      
 82      
 83      bool cut ;
 84      
 85      cut = Q2 > 12 && y < 0.6 && w2 > 3000 ;
 86      
 87      
 88      if ( !cut ) vetoEvent ; 
 89      
 90      const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
 91      if ( dl.failed() ) vetoEvent;
 92      /*
 93      cout << "  scattered lepton angle " << 180.- dl.out().momentum().angle(dl.in().momentum())/degree << endl;
 94      cout << " in  lepton " << dl.in().momentum() << endl;
 95      cout << " out lepton " << dl.out().momentum() << endl;
 96      */
 97      
 98      const FinalState& fs = apply<FinalState>(event, "FS");
 99      Particles particles; particles.reserve(fs.size());
100      ConstGenParticlePtr dislepGP = dl.out().genParticle();
101      for (const Particle& p : cfs.particles()) {
102        ConstGenParticlePtr loopGP = p.genParticle();
103        if (loopGP == dislepGP) continue;
104        particles.push_back(p);
105      }
106
107
108      double efwd = 0.;
109
110      for (const Particle& p : particles) {
111        const double th = 180. - p.momentum().angle(dl.in().momentum())/degree;
112
113        if (inRange(th, 4.4, 15.0)) { 
114           efwd += p.E();
115           //cout << " angle " << th << " pid " << p.pid() << " Efwd = " << efwd << endl;
116           } 
117      }
118
119      bool evcut[2];
120      evcut[0] =   efwd > 0.5;
121      
122       // fill the counter
123      _Nevt_after_cuts -> fill();
124      if (Q2 > 100 )  _Nevt_afterh_cuts -> fill();
125       
126      if ( evcut[0] && (Q2 < 80 ) ) _Nevt_afterfwd_cuts -> fill();
127      if ( evcut[0] && (Q2 > 100 ) ) _Nevt_afterhfwd_cuts -> fill();
128
129      double n_charg = 0; 
130
131      // Boost to Breit
132      const LorentzTransform breitboost = dk.boostBreit();
133       
134      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
135        const Particle& p = particles[ip1];
136        const FourMomentum BreMom = breitboost.transform(p.momentum());
137       // cout << BreMom.pz() << endl;
138        double x = cos(BreMom.theta());
139        
140        if (Q2 < 80 ) _h["costh_lowQ"] ->fill(x);
141        if (Q2 > 100 ) _h["costh_highQ"]  ->fill(x);
142        
143        if (Q2 < 80  && evcut[0]) _h["costh_lowQ_noEfwd"] ->fill(x);
144        if (Q2 > 100 && evcut[0]) _h["costh_highQ_noEfwd"]  ->fill(x);
145
146
147        if ( BreMom.pz() > 0. ) continue;
148        double pcal= sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2()) ;
149        double xp = 2*pcal/(sqrt(Q2));
150        double xi = log(1/xp);
151        
152        double charge = p.charge() ;
153       // cout << " charge " << charge << endl;
154       
155        if (charge > 0 ) {
156          if (Q2 < 80 ) _h["xp_posCharge_lowQ"] -> fill(xp);
157          if (Q2 > 100 ) _h["xp_posCharge_highQ"] -> fill(xp);
158        } else {
159         if (Q2 < 80 ) _h["xp_negCharge_lowQ"] -> fill(xp);
160         if (Q2 > 100 ) _h["xp_negCharge_highQ"] -> fill(xp);
161        }
162    
163        if (Q2 < 80 ) _h["ksi_lowQ"] -> fill(xi);
164        if (Q2 > 100 ) _h["ksi_highQ"] -> fill (xi);
165 
166        n_charg = n_charg + 1; 
167
168     }
169     _h["Mult_vrs_Q2_nchrg"] -> fill(Q2,n_charg) ;
170     _h["Mult_vrs_Q2_count"] -> fill(Q2) ;
171     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_nchrg"] -> fill(Q2,n_charg) ;
172     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_count"] -> fill(Q2) ;
173}
174
175    /// Normalise histograms etc., after the run
176    void finalize() {
177     
178   
179      normalize(_h["xp_posCharge_lowQ"]);
180      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_posCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
181
182      normalize(_h["xp_negCharge_lowQ"]);
183      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_negCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
184
185      scale(_h["costh_lowQ"], 1.0/ *_Nevt_after_cuts);
186      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["costh_highQ"], 1.0/ *_Nevt_afterh_cuts);
187      //cout << " after fwd cuts " << dbl(*_Nevt_afterfwd_cuts) << endl;
188      
189      if(dbl(*_Nevt_afterfwd_cuts)>0) scale(_h["costh_lowQ_noEfwd"], 1.0/ *_Nevt_afterfwd_cuts);
190      if(dbl(*_Nevt_afterhfwd_cuts)>0) scale(_h["costh_highQ_noEfwd"], 1.0/ *_Nevt_afterhfwd_cuts);
191      
192      scale(_h["ksi_lowQ"], 1.0/ *_Nevt_after_cuts);
193      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["ksi_highQ"], 1.0/ *_Nevt_afterh_cuts);
194      //cout << " Nevt " << dbl(*_Nevt_after_cuts) << endl;
195
196      divide(_h["Mult_vrs_Q2_nchrg"], _h["Mult_vrs_Q2_count"], _s["Mult_vrs_Q2"]);
197      divide(_h["Mult_vrs_Q2_noEfwd_nchrg"], _h["Mult_vrs_Q2_noEfwd_count"], _s["Mult_vrs_Q2_noEfwd"]);
198    }
199
200    ///@}
201
202
203    /// @name Histograms
204    ///@{
205    map<string, Histo1DPtr> _h;
206    map<string, Profile1DPtr> _p;
207    map<string, CounterPtr> _c;
208    map<string, Scatter2DPtr> _s;
209    CounterPtr _Nevt_after_cuts;
210    CounterPtr _Nevt_afterfwd_cuts;
211    CounterPtr _Nevt_afterh_cuts;
212    CounterPtr _Nevt_afterhfwd_cuts;
213    ///@}
214
215
216  };
217
218
219  RIVET_DECLARE_PLUGIN(H1_1995_I394793);
220
221}