1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DISFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/DISKinematics.hh"
  7#include "Rivet/Projections/DISLepton.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)
 14  class H1_2007_I736052 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I736052);
 19
 20
 21    /// @name Analysis methods
 22    ///@{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27        declare(DISLepton(), "Lepton");
 28        declare(DISKinematics(), "Kinematics");
 29
 30        const FinalState fs;
 31        declare(fs, "FS");
 32        const UnstableParticles ufs;
 33        declare(ufs, "UFS");
 34
 35        double jet_radius = 1.0;
 36        const DISFinalState DISfs(DISFrame::BREIT);
 37        declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme, jet_radius), "jets_fs");
 38
 39        Histo1DPtr tmp;
 40        book(_h["411"], 4, 1, 1);
 41        book(_h["511"], 5, 1, 1);
 42        book(_h["611"], 6, 1, 1);
 43        book(_h["711"], 7, 1, 1);
 44        book(_h["811"], 8, 1, 1);
 45        book(_h["911"], 9, 1, 1);
 46        book(_h_binned["Q2xbj"], {2., 4.22, 10., 17.8, 31.6, 100.},
 47                                 {"d10-x01-y01", "d11-x01-y01", "d12-x01-y01",
 48                                  "d13-x01-y01", "d14-x01-y01"});
 49        book(_h["1511"], 15, 1, 1);
 50        book(_h["1611"], 16, 1, 1);
 51        book(_h["1711"], 17, 1, 1);
 52        book(_h["1811"], 18, 1, 1);
 53        book(_h["1911"], 19, 1, 1);
 54        book(_h["2011"], 20, 1, 1);
 55        book(_h["2111"], 21, 1, 1);
 56        book(_h["2211"], 22, 1, 1);
 57
 58        book(_h_binned["Q2phi"], {2., 10., 100.}, {"d23-x01-y01", "d24-x01-y01"});
 59        book(_h["2511"], 25, 1, 1);
 60        book(_h["2611"], 26, 1, 1);
 61        book(_h["2711"], 27, 1, 1);
 62        book(_h["2811"], 28, 1, 1);
 63        book(_h_binned["Q2xgam"], {2., 5., 10., 100.}, {"d29-x01-y01", "d29-x01-y02", "d29-x01-y03"});
 64        book(_h["3011"], 30, 1, 1);
 65        book(_h_binned["Q2xglue"], {2., 5., 10., 100.}, {"d31-x01-y01", "d31-x01-y02", "d31-x01-y03"});
 66    }
 67
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event& event) {
 71
 72        const FinalState& fs = apply<FinalState>(event, "FS");
 73        const size_t numParticles = fs.particles().size();
 74
 75        Jets jets_fs = apply<JetFinder>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
 76        double jet_radius = 1.0;
 77
 78        const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
 79
 80        if (numParticles < 2){
 81            MSG_DEBUG("Failed leptonic event cut");
 82            vetoEvent;
 83        }
 84
 85        Particles Dstar;
 86        for(const Particle& p : select(ufs.particles(), Cuts::pT > 1.5*GeV and Cuts::pT < 15*GeV and Cuts::abseta < 1.5 and Cuts::abspid==413)) {
 87            Dstar.push_back(p);
 88        }
 89        if(Dstar.size() == 0){ // Cut on Dstar
 90            MSG_DEBUG("Failed Dstar cut");
 91            vetoEvent;
 92        }
 93
 94        const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 95        const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 96
 97
 98        double Q2 = dk.Q2();
 99        double y  = dk.y();
100
101        if(y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100){ // Cut on event kinematics
102            MSG_DEBUG("Failed kinematics cut");
103            vetoEvent;
104        }
105
106        // Extract the particles other than the lepton
107        Particles particles;
108        particles.reserve(fs.particles().size());
109        ConstGenParticlePtr dislepGP = dl.out().genParticle();
110        for (const Particle& p : fs.particles()) {
111           ConstGenParticlePtr loopGP = p.genParticle();
112           if (loopGP == dislepGP) continue;
113           particles.push_back(p);
114        }
115
116
117
118        const LorentzTransform hcmboost = dk.boostHCM(); // Hadron cm system
119        const LorentzTransform breitboost = dk.boostBreit(); //Breit system
120        const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit
121
122        double xbj  = dk.x();
123        double W2 = dk.W2();
124        double pT_cm(0);
125
126        _h["411"] -> fill(Q2);
127        _h["511"] -> fill(xbj);
128        _h["611"] -> fill(std::sqrt(W2));
129
130        _h_binned["Q2xbj"]->fill(Q2, xbj);
131
132        for(const Particle& p : Dstar){
133            _h["711"] -> fill(p.pT());
134            _h["811"] -> fill(p.eta());
135            _h["911"] -> fill((p.E() - p.pz())/(2*y*dk.beamLepton().E()));
136
137            const FourMomentum hcmMom = hcmboost.transform(p.momentum());
138            if(pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();
139
140            if(hcmMom.pT() > 2){              //Cut for 1711 and 1811 histograms
141                _h["1711"] -> fill(p.pT());
142                _h["1811"] -> fill(p.eta());
143            }
144        }
145
146        if(pT_cm > 2){
147            _h["1511"] -> fill(Q2);
148            _h["1611"] -> fill(xbj);
149        }
150
151        bool two_jets_with_cut(false);
152
153        Jets jet_cut;
154
155        for(const Jet& j : jets_fs){
156            double etalab = labboost.transform(j.momentum()).eta() ;
157            if(j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5){  // Cut on jet energy in Breit sys.
158                jet_cut.push_back(j);                            }
159        }
160
161        if(jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4 ) two_jets_with_cut = true;
162
163        double delta_phi;
164        FourMomentum jet1;
165        FourMomentum jet2;
166
167        Jets jet_DJ, jet_OJ;
168        bool found_DJ(false), found_OJ(false);
169        //cout << " check D and other jets " << endl;
170        if (two_jets_with_cut){
171            jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys.
172            jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys.
173
174            delta_phi = deltaPhi(jet1,jet2);
175            _h["1911"] -> fill(Q2);
176            _h["2011"] -> fill(xbj);
177            _h["2111"] -> fill(jet_cut[0].momentum().Et()/GeV);
178            _h["2211"] -> fill(FourMomentum(jet1+jet2).mass()/GeV);  // Jets invariant mass in Breit sys.
179            _h_binned["Q2phi"]->fill(Q2, delta_phi);
180            for (const Jet& jet : jet_cut) {
181                for(const Particle & p : Dstar) {
182                    if(deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius  ) {
183                        jet_DJ.push_back(jet);
184                        found_DJ = true;
185                    }
186                    else {
187                    jet_OJ.push_back(jet);
188                    found_OJ = true;
189                    }
190                }
191
192               if( found_DJ &&  found_OJ ) break ;
193            }
194            if(jet_DJ.size()>0 && jet_OJ.size()>0 ) {
195                double eta_DJ_breit = jet_DJ[0].momentum().eta();
196                double eta_OJ_breit = jet_OJ[0].momentum().eta();
197
198                _h["2511"] -> fill(eta_DJ_breit);
199                _h["2611"] -> fill(eta_OJ_breit);
200                _h["2711"] -> fill(abs(eta_DJ_breit - eta_OJ_breit));
201
202                double x_gamma, x_gluon;
203                double E_star_p_z_star(0); // E() - pz() sum for for x_gamma in \gamma p cm
204                double E_star_p_z_had(0);
205
206                // for jets in had cms: boost first back to lab and then to hcm
207               Jet jet_DJ_hcm, jet_OJ_hcm;
208               jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
209               jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));
210
211               // Need to change sign: by default hcmboost has gamma* in +z dir,
212               // and p in -z dir, but here we need: gamma* -in -z and proton in +z.
213               // so - -> +: watch out for E-pz -> E+pz  and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
214               // this was already noted in  H1_2007_I746380
215
216               double Et_DJ_hcm = jet_DJ_hcm.Et();
217               double eta_DJ_hcm = jet_DJ_hcm.eta();
218               double Et_OJ_hcm = jet_OJ_hcm.Et();
219               double eta_OJ_hcm = jet_OJ_hcm.eta();
220
221                //observed fraction of the photon momentum carried by the parton involved in the hard subprocess
222
223                E_star_p_z_star = jet_DJ_hcm.momentum().E() + jet_DJ_hcm.momentum().pz() + jet_OJ_hcm.momentum().E() + jet_OJ_hcm.momentum().pz();
224
225                for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
226                  const Particle& p = particles[ip1];
227                  E_star_p_z_had += (hcmboost.transform(p.momentum())).E() + (hcmboost.transform(p.momentum())).pz();
228                }
229                x_gamma = E_star_p_z_star/E_star_p_z_had;
230                x_gluon = (Et_OJ_hcm*exp(-eta_OJ_hcm) + Et_DJ_hcm*exp(-eta_DJ_hcm))/(2.*hcmboost.transform(dk.beamHadron().momentum()).E());
231                //observed fraction of the proton momentum carried by the gluon
232
233                _h["2811"]->fill(x_gamma);
234                _h_binned["Q2xgam"]->fill(Q2, x_gamma);
235
236                _h["3011"]->fill(log10(x_gluon));
237                _h_binned["Q2xglue"]->fill(Q2, log10(x_gluon));
238            }
239        }
240    }
241
242
243    /// Normalise histograms etc., after the run
244    void finalize() {
245
246        const double norm = crossSection()/nanobarn/sumW();
247        scale(_h, norm);
248        scale(_h_binned, norm);
249        divByGroupWidth(_h_binned);
250
251        // new scaling needed, since x bins are in log10(x)
252        for (auto& b : _h["3011"]->bins()) {
253           const double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
254           // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
255           const double factor = pow(10,b.xMid())/2.3026;
256           b.scaleW(scale_new/factor) ;
257        }
258
259        // new scaling needed, since x bins are in log10(x)
260        for (auto& histo : _h_binned["Q2xglue"]->bins()) {
261           for (auto& b : histo->bins()) {
262              const double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
263              // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
264              const double factor = pow(10,b.xMid())/2.3026;
265              b.scaleW(scale_new/factor) ;
266           }
267        }
268
269    }
270
271    ///@}
272
273
274    /// @name Histograms
275    ///@{
276    map<string, Histo1DPtr> _h;
277    map<string, Histo1DGroupPtr> _h_binned;
278    ///@}
279
280
281  };
282
283
284  RIVET_DECLARE_PLUGIN(H1_2007_I736052);
285
286}