1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/ParticleIdUtils.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7#include "Rivet/Projections/DISKinematics.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Transverse momentum spectra of charged particles in DIS
 13  class H1_1996_I424463 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I424463);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Book projections
 27      declare(DISLepton(), "Lepton");
 28      declare(DISKinematics(), "Kinematics");
 29      declare(ChargedFinalState(), "CFS");
 30      declare(FinalState(), "FS");
 31
 32      // Book histograms
 33      book(_NevAll, "TMP/Nev_all");
 34      int ixx = 0 ;
 35      for (size_t ix = 0; ix < 10; ++ix) {
 36        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
 37        for(unsigned int ih=0;ih<2;++ih) {
 38          book(_h_dndpt_eta_bin[ih][ix], ih*10 + ix+1 , 1, 1);
 39          book(_hdndeta_bin    [ih][ix], ih*10 + ix+29, 1, 1);
 40        }
 41        if (ix != 6 && ix != 9) {
 42          book(_hdndptmax_low_eta_bin[ixx], ixx+21, 1, 1);
 43          ixx=ixx+1;
 44        }
 45      }
 46    }
 47
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      if (_edgespT.empty()) {
 52        _edgespT    = _h_dndpt_eta_bin   [0][0]->xEdges();
 53        _edgesEta   = _hdndeta_bin       [0][0]->xEdges();
 54        _edgespTMax = _hdndptmax_low_eta_bin[0]->xEdges();
 55      }
 56      const FinalState& fs = apply<FinalState>(event, "FS");
 57      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 58      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 59
 60      // Get the DIS kinematics
 61      double x  = dk.x();
 62      double y = dk.y();
 63      double Q2 = dk.Q2()/GeV;
 64      double W2 = dk.W2()/GeV;
 65
 66      // Momentum of the scattered lepton
 67      FourMomentum leptonMom = dl.out().momentum();
 68      double enel = leptonMom.E()/GeV;
 69      double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
 70
 71      _NevAll -> fill() ;
 72
 73      bool cut = y > 0.05 && Q2 > 5. && Q2 < 100.&& enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
 74      if (!cut) vetoEvent;
 75
 76      int ibin[10] ;
 77      for (int i = 0; i < 10; i++) ibin[i] = 0;
 78      if ( 5. < Q2 && Q2 < 50. && x > 0.0001 && x < 0.0010) ibin[0] = 1;
 79      if ( 5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.0002) ibin[1] = 1;
 80      if ( 6. < Q2 && Q2 < 10. && x > 0.0002 && x < 0.0005) ibin[2] = 1;
 81      if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.0005) ibin[3] = 1;
 82      if (10. < Q2 && Q2 < 20. && x > 0.0005 && x < 0.0008) ibin[4] = 1;
 83      if (10. < Q2 && Q2 < 20. && x > 0.0008 && x < 0.0015) ibin[5] = 1;
 84      if (10. < Q2 && Q2 < 20. && x > 0.0015 && x < 0.0040) ibin[6] = 1;
 85      if (20. < Q2 && Q2 < 50. && x > 0.0005 && x < 0.0014) ibin[7] = 1;
 86      if (20. < Q2 && Q2 < 50. && x > 0.0014 && x < 0.0030) ibin[8] = 1;
 87      if (20. < Q2 && Q2 < 50. && x > 0.0030 && x < 0.0100) ibin[9] = 1;
 88      for (int i = 0; i < 10; i++) {
 89        if (ibin[i]==1) _Nevt_after_cuts[i]->fill();
 90      }
 91
 92      // Extract the particles other than the lepton
 93      /// @todo Improve to avoid HepMC digging
 94      Particles particles;
 95      particles.reserve(fs.particles().size());
 96      ConstGenParticlePtr dislepGP = dl.out().genParticle();
 97      for (const Particle& p : fs.particles()) {
 98        ConstGenParticlePtr loopGP = p.genParticle();
 99        if (loopGP == dislepGP) continue;
100        particles.push_back(p);
101      }
102
103      // Boost to hadronic CM
104      const LorentzTransform hcmboost = dk.boostHCM();
105
106      // Loop over the particles
107      int mult = 0 ;
108      double ptmax_high[10], ptmax_low[10];
109      for (int i = 0; i < 10; i++) {
110        ptmax_high[i] = 0.; ptmax_low[i] = 0.;
111      }
112      double EtSum = 0;
113      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
114        const Particle& p = particles[ip1];
115
116        double eta = p.momentum().pseudorapidity();
117        // Boost to hcm
118        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
119
120        // Apply safety cuts
121        if (eta > -5  && eta < 10.) {
122          mult = mult + 1;
123          double pThcm =hcmMom.pT() ;
124          double etahcm = hcmMom.pseudorapidity();
125          if (etahcm > 0. && etahcm < 2.0) {
126            EtSum = EtSum + hcmMom.Et();
127          }
128          if (PID::charge(p.pid()) != 0) {
129            if (etahcm > 0.5 && etahcm < 1.5) {
130              for ( int i=0; i< 10; i++) {
131                if (ibin[i]==1) {
132                  fillpT(1,i,pThcm);
133                  if (pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
134                }
135              }
136
137            }
138            if (etahcm > 1.5 && etahcm < 2.5){
139              for ( int i=0; i< 10; i++) {
140                if (ibin[i]==1) {
141                  fillpT(0,i,pThcm);
142                  if (pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
143                }
144              }
145            }
146            for ( int i=0; i< 10; i++) {
147              if (ibin[i]==1)               fillEta(1,i,etahcm);
148              if (ibin[i]==1 && pThcm > 1.) fillEta(0,i,etahcm);
149            }
150          }
151        }  // end of loop over the particles
152      }
153      int ii=0;
154      for ( int i=0; i< 10; i++) {
155        if (i != 6 && i != 9 ) {
156          if ( ibin[i]==1 && EtSum > 6. ) {
157            fillpTMax(ii,ptmax_low[i]);
158          }
159          ii=ii+1;
160        }
161      }
162    }
163
164
165    /// Normalise histograms etc., after the run
166    void finalize() {
167      MSG_DEBUG("All events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[0]->val());
168      MSG_DEBUG("Cut1 events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[1]->val());
169      int ii = 0;
170      for ( int i=0; i< 10; i++) {
171        if (_Nevt_after_cuts[i]->val()  != 0) {
172          for(unsigned int ih=0;ih<2;++ih) {
173            scale(_h_dndpt_eta_bin[ih][i], 1./ *_Nevt_after_cuts[i]);
174            size_t ioff=0;
175            if( (ih==0 && (i==6 || i==9)) ||
176                (ih==1 && (i<=4 || i==7)))    ioff=1;
177            else if (ih==1 && (i==5 || i==8)) ioff=2;
178            for(auto & b : _h_dndpt_eta_bin[ih][i]->bins()) {
179              const size_t idx = b.index()+ioff;
180              b.scaleW(1./_axispT.width(idx));
181            }
182            scale(_hdndeta_bin    [ih][i], 1./ *_Nevt_after_cuts[i]);
183            ioff=0;
184            if( (ih==0 && (i==5 || i==6 || i==8 || i==9)) ||
185                (ih==1 && (i==0 || i==4 || i==5 || i==8))) ioff=1;
186            else if (ih==1 && (i==6 || i==9))       ioff=2;
187            for(auto & b : _hdndeta_bin[ih][i]->bins()) {
188              const size_t idx = b.index()+ioff;
189              b.scaleW(1./_axisEta.width(idx));
190            }
191          }
192        }
193        if ( i !=6 && i!=9) {
194          if (_Nevt_after_cuts[i]->val()  != 0) {
195            normalize(_hdndptmax_low_eta_bin[ii]);
196            for(auto & b : _hdndptmax_low_eta_bin[ii]->bins()) {
197              const size_t idx = b.index();
198              b.scaleW(1./_axispTMax.width(idx));
199            }
200          }
201          ii=ii+1;
202        }
203      }
204      
205    }
206    /// @}
207
208    void fillEta(const unsigned int ix, const unsigned int iy, const double value) {
209      string edge = "OTHER";
210      const size_t idx = _axisEta.index(value);
211      if (idx && idx <= _edgesEta.size()) {
212        if ( ( (ix==0 && iy>=1&& iy<=3) || (ix==0 && (iy==1 || iy==3)) ) && idx == _edgesEta.size())
213          ;
214        else if (idx==1 && ( (ix==0 && (iy==5 || iy==6 || iy==8 || iy==9)) ||
215                             (ix==1 && (iy==0 || iy==4 || iy==5 || iy==8))))
216          ;
217        else if (idx<=2 && (ix==1 && (iy==6 || iy==9)) )
218          ;
219        else {
220          edge = _edgesEta[idx-1];
221        }
222      }
223      _hdndeta_bin[ix][iy]->fill(edge);
224    }
225    
226    void fillpT(const unsigned int ix, const unsigned int iy, const double value) {
227      string edge = "OTHER";
228      const size_t idx = _axispT.index(value);
229      if (idx && idx <= _edgespT.size()) {
230        if(((ix==0 && iy==8) || (ix==1 && (iy==2 || iy==4 || iy==6))) && idx == _edgespT.size())
231          ;
232        else if (idx==1 && ( (ix==0 && (iy==6 || iy==9)) ||
233                             (ix==1 && (iy<=4 || iy==7))))
234          ;
235        else if (idx<=2 && (ix==1 && (iy==5 || iy==8)))
236          ;
237        else
238          edge = _edgespT[idx-1];
239      }
240      _h_dndpt_eta_bin[ix][iy]->fill(edge);
241    }
242
243    void fillpTMax(const unsigned int ix, const double value) {
244      string edge = "OTHER";
245      const size_t idx = _axispTMax.index(value);
246      if (idx && idx <= _edgespTMax.size()) {
247        if(idx==_edgespTMax.size() && (ix==1 || ix==2 || ix==4))
248          ;
249        else if(idx>=_edgespTMax.size()-1 && ix==5)
250          ;
251        else
252          edge = _edgespTMax[idx-1];
253      }
254      _hdndptmax_low_eta_bin[ix]->fill(edge);
255    }
256
257    /// @name Histograms
258    /// @{
259    array<CounterPtr,10> _Nevt_after_cuts;
260    BinnedHistoPtr<string> _h_dndpt_eta_bin[2][10], _hdndeta_bin[2][10], _hdndptmax_low_eta_bin[8];
261    CounterPtr _NevAll ;
262
263    vector<string> _edgespT,_edgesEta,_edgespTMax;
264    YODA::Axis<double> _axispT    = YODA::Axis<double>{ 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.825, 2.125, 2.525, 3.125, 4.0, 5.0};
265    YODA::Axis<double> _axisEta   = YODA::Axis<double>{ 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.00};
266    YODA::Axis<double> _axispTMax = YODA::Axis<double>{0.005, 0.255, 0.505, 0.755, 1.005, 1.255, 1.505, 1.755, 2.065, 2.5, 3.125, 4.0, 5.0};
267    /// @}
268
269  };
270
271
272  RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
273
274}