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(DISFinalState::BoostFrame::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        _h_binned["Q2xbj"].add( 2., 4.22, book(_h["1011"], 10, 1, 1));        
 47        _h_binned["Q2xbj"].add( 4.22, 10., book(_h["1111"], 11, 1, 1));
 48        _h_binned["Q2xbj"].add( 10., 17.8, book(_h["1211"], 12, 1, 1));
 49        _h_binned["Q2xbj"].add( 17.8, 31.6, book(_h["1311"], 13, 1, 1));
 50        _h_binned["Q2xbj"].add( 31.6, 100., book(_h["1411"], 14, 1, 1));
 51        book(_h["1511"], 15, 1, 1);
 52        book(_h["1611"], 16, 1, 1);
 53        book(_h["1711"], 17, 1, 1);
 54        book(_h["1811"], 18, 1, 1);
 55        book(_h["1911"], 19, 1, 1);
 56        book(_h["2011"], 20, 1, 1);
 57        book(_h["2111"], 21, 1, 1);
 58        book(_h["2211"], 22, 1, 1);
 59
 60        _h_binned["Q2phi"].add( 2., 10., book(tmp, 23, 1, 1));
 61        _h_binned["Q2phi"].add( 10., 100., book(tmp, 24, 1, 1));
 62        book(_h["2511"], 25, 1, 1);
 63        book(_h["2611"], 26, 1, 1);
 64        book(_h["2711"], 27, 1, 1);
 65        book(_h["2811"], 28, 1, 1);
 66        _h_binned["Q2xgam"].add( 2., 5., book(_h["2911"], 29, 1, 1));
 67        _h_binned["Q2xgam"].add( 5., 10., book(_h["2912"], 29, 1, 2));
 68        _h_binned["Q2xgam"].add( 10., 100., book(_h["2913"], 29, 1, 3));
 69        book(_h["3011"], 30, 1, 1);
 70        _h_binned["Q2xglue"].add( 2., 5., book(_h["3111"], 31, 1, 1));
 71        _h_binned["Q2xglue"].add( 5., 10., book(_h["3112"], 31, 1, 2));
 72        _h_binned["Q2xglue"].add( 10., 100., book(_h["3113"], 31, 1, 3));
 73    }
 74
 75
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78        
 79        const FinalState& fs = apply<FinalState>(event, "FS");
 80        const size_t numParticles = fs.particles().size();
 81        
 82        Jets jets_fs = apply<JetAlg>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
 83        double jet_radius = 1.0;
 84        
 85        const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
 86        
 87        if (numParticles < 2){
 88            MSG_DEBUG("Failed leptonic event cut");
 89            vetoEvent;
 90        }
 91        
 92        Particles Dstar;
 93        for(const Particle& p : filter_select(ufs.particles(), Cuts::pT > 1.5*GeV and Cuts::pT < 15*GeV and Cuts::abseta < 1.5 and Cuts::abspid==413)) {
 94            Dstar.push_back(p);
 95        }
 96        if(Dstar.size() == 0){ // Cut on Dstar
 97            MSG_DEBUG("Failed Dstar cut");
 98            vetoEvent;
 99        }
100            
101        const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
102        const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
103        
104       
105        double Q2 = dk.Q2();
106        double y  = dk.y();
107
108        if(y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100){ // Cut on event kinematics
109            MSG_DEBUG("Failed kinematics cut");
110            vetoEvent; 
111        } 
112        
113        // Extract the particles other than the lepton
114        Particles particles;
115        particles.reserve(fs.particles().size());
116        ConstGenParticlePtr dislepGP = dl.out().genParticle();
117        for (const Particle& p : fs.particles()) {
118           ConstGenParticlePtr loopGP = p.genParticle();
119           if (loopGP == dislepGP) continue;
120           particles.push_back(p);
121        }
122
123        
124        
125        const LorentzTransform hcmboost = dk.boostHCM(); // Hadron cm system 
126        const LorentzTransform breitboost = dk.boostBreit(); //Breit system
127        const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit
128        
129        double xbj  = dk.x();
130        double W2 = dk.W2();
131        double pT_cm(0); 
132        
133        _h["411"] -> fill(Q2);
134        _h["511"] -> fill(xbj);
135        _h["611"] -> fill(std::sqrt(W2)); 
136        
137        _h_binned["Q2xbj"].fill(Q2, xbj);      
138        
139        for(const Particle& p : Dstar){
140            _h["711"] -> fill(p.pT());
141            _h["811"] -> fill(p.eta());
142            _h["911"] -> fill((p.E() - p.pz())/(2*y*dk.beamLepton().E()));
143            
144            const FourMomentum hcmMom = hcmboost.transform(p.momentum()); 
145            if(pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();
146            
147            if(hcmMom.pT() > 2){              //Cut for 1711 and 1811 histograms
148                _h["1711"] -> fill(p.pT());
149                _h["1811"] -> fill(p.eta());
150            }            
151        }
152
153        if(pT_cm > 2){
154            _h["1511"] -> fill(Q2);
155            _h["1611"] -> fill(xbj);
156        }     
157        
158/*
159        FourMomentum gammaZ;
160        const FourMomentum leptonOUT = dl.out();
161        const FourMomentum leptonIN = dl.in();
162        gammaZ = leptonIN - leptonOUT ;
163        cout << " LAB: gammaZ " <<  gammaZ <<  endl; 
164        cout << " HCM: gammaZ " <<   hcmboost.transform(gammaZ) <<  endl; 
165        cout << " LAB: p-beam " <<  dk.beamHadron().momentum() <<  endl; 
166        cout << " HCM: p-beam " <<   hcmboost.transform(dk.beamHadron().momentum()) <<  endl; 
167        cout << " Breit: gammaZ " <<   breitboost.transform(gammaZ) <<  endl; 
168        cout << " Breit: p-beam " <<   breitboost.transform(dk.beamHadron().momentum()) <<  endl; 
169        for(ConstGenParticlePtr p: HepMCUtils::particles(event.genEvent())) {
170          const PdgId pid = p->pdg_id();
171          if (abs(pid) == 23) {
172          cout<< " HEPMC: photon/Z " << p->momentum() << endl;
173          }
174        }        
175*/
176        bool two_jets_with_cut(false); 
177        
178        Jets jet_cut; 
179        
180        for(const Jet& j : jets_fs){
181            double etalab = labboost.transform(j.momentum()).eta() ;
182            if(j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5){  // Cut on jet energy in Breit sys.
183                jet_cut.push_back(j);                            }     
184        }  
185
186        if(jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4 ) two_jets_with_cut = true;
187        
188        double delta_phi;
189        FourMomentum jet1;
190        FourMomentum jet2;
191        
192        Jets jet_DJ, jet_OJ;
193        bool found_DJ(false), found_OJ(false);
194        //cout << " check D and other jets " << endl;        
195        if(two_jets_with_cut){
196            jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys. 
197            jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys. 
198
199            delta_phi = deltaPhi(jet1,jet2)/degree ;
200            _h["1911"] -> fill(Q2);     
201            _h["2011"] -> fill(xbj); 
202            _h["2111"] -> fill(jet_cut[0].momentum().Et()/GeV); 
203            _h["2211"] -> fill(FourMomentum(jet1+jet2).mass()/GeV);  // Jets invariant mass in Breit sys.
204            //cout << " dphi " << delta_phi <<" " << deltaPhi(jet_cut[0].momentum(),jet_cut[1].momentum())/degree <<  endl;
205            _h_binned["Q2phi"].fill(Q2, delta_phi);
206            for (const Jet& jet : jet_cut) {
207                for(const Particle & p : Dstar) {
208                    if(deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius  ) {
209                        jet_DJ.push_back(jet);
210                        found_DJ = true;
211                    }
212                    else {
213                    jet_OJ.push_back(jet);
214                    found_OJ = true;
215                    }
216                }                
217
218               if( found_DJ &&  found_OJ ) break ;  
219            }
220            //cout << " DJ jet size " << jet_DJ.size() << "found Dj " << found_DJ <<  " OJ jet size " << jet_DJ.size()  << " found OJ " << found_OJ << endl;
221            if(jet_DJ.size()>0 && jet_OJ.size()>0 )
222            {            
223                double eta_DJ_breit = jet_DJ[0].momentum().eta();
224                double eta_OJ_breit = jet_OJ[0].momentum().eta(); 
225                
226                _h["2511"] -> fill(eta_DJ_breit);
227                _h["2611"] -> fill(eta_OJ_breit);
228                _h["2711"] -> fill(abs(eta_DJ_breit - eta_OJ_breit));    
229                
230                double x_gamma, x_gluon;
231                double E_star_p_z_star(0); // E() - pz() sum for for x_gamma in \gamma p cm
232                double E_star_p_z_had(0);
233
234                // for jets in had cms: boost first back to lab and then to hcm
235               Jet jet_DJ_hcm, jet_OJ_hcm;
236               jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
237               jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));
238
239// Need to change sign: by default hcmboost has gamma* in +z dir, and p in -z dir, but here we need: gamma* -in -z and proton in +z.
240// so - -> +: watch out for E-pz -> E+pz  and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
241// this was already noted in  H1_2007_I746380.cc
242
243
244   
245               double Et_DJ_hcm = jet_DJ_hcm.Et();
246               double eta_DJ_hcm = jet_DJ_hcm.eta();
247               double Et_OJ_hcm = jet_OJ_hcm.Et();
248               double eta_OJ_hcm = jet_OJ_hcm.eta();
249
250                //observed fraction of the photon momentum carried by the parton involved in the hard subprocess 
251                
252                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();
253
254                for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
255                  const Particle& p = particles[ip1];
256                  E_star_p_z_had += (hcmboost.transform(p.momentum())).E() + (hcmboost.transform(p.momentum())).pz();
257                }
258                x_gamma = E_star_p_z_star/E_star_p_z_had;
259                x_gluon = (Et_OJ_hcm*exp(-eta_OJ_hcm) + Et_DJ_hcm*exp(-eta_DJ_hcm))/(2.*hcmboost.transform(dk.beamHadron().momentum()).E());
260                //observed fraction of the proton momentum carried by the gluon 
261                
262                //cout << " xgamm " << x_gamma << " x_glu " << x_gluon << endl;
263                _h["2811"] -> fill(x_gamma); 
264                _h_binned["Q2xgam"].fill(Q2, x_gamma);  
265                
266                _h["3011"] -> fill(log10(x_gluon));
267                _h_binned["Q2xglue"].fill(Q2, log10(x_gluon));                  
268            }
269        }        
270    }
271
272
273    /// Normalise histograms etc., after the run
274    void finalize() {
275        double norm = crossSection()/nanobarn/sumW(); 
276        //double norm = crossSection()/nanobarn/sumOfWeights(); 
277        //cout << " SumOfWeights " << sumW() << " "<< sumOfWeights() << endl;
278       
279        scale(_h["411"], norm);
280        scale(_h["511"], norm);
281        scale(_h["611"], norm);
282        scale(_h["711"], norm);
283        scale(_h["811"], norm);
284        scale(_h["911"], norm);
285        _h_binned["Q2xbj"].scale(norm, this);
286        scale(_h["1511"], norm);
287        scale(_h["1611"], norm);
288        scale(_h["1711"], norm);
289        scale(_h["1811"], norm);
290        scale(_h["1911"], norm);
291        scale(_h["2011"], norm);
292        scale(_h["2111"], norm);
293        scale(_h["2211"], norm);
294        _h_binned["Q2phi"].scale(norm*180./M_PI, this);
295        scale(_h["2511"], norm);
296        scale(_h["2611"], norm);
297        scale(_h["2711"], norm);
298        scale(_h["2811"], norm);
299        _h_binned["Q2xgam"].scale(norm, this);
300        
301        scale(_h["3011"], norm);
302        // new scaling needed, since x bins are in log10(x)
303        vector<YODA::HistoBin1D>& bins = _h["3011"] -> bins(); 
304        for (auto & b : bins) {          
305           double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
306           // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x 
307           double factor = pow(10,b.xMid())/2.3026;
308           b.scaleW(scale_new/factor) ;
309        }
310        
311        _h_binned["Q2xglue"].scale(norm, this);
312        // new scaling needed, since x bins are in log10(x)
313        for (Histo1DPtr histo : _h_binned["Q2xglue"].histos()) {
314           vector<YODA::HistoBin1D>& bins = histo -> bins(); 
315           for (auto & b : bins) {             
316              double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
317              // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x 
318              double factor = pow(10,b.xMid())/2.3026;
319              b.scaleW(scale_new/factor) ;
320           }
321
322        }
323
324    }
325
326    ///@}
327
328
329    /// @name Histograms
330    ///@{
331    map<string, Histo1DPtr> _h;
332    map<string, Profile1DPtr> _p;
333    map<string, CounterPtr> _c;
334    map<string, BinnedHistogram> _h_binned;
335    ///@}
336
337
338  };
339
340
341  RIVET_DECLARE_PLUGIN(H1_2007_I736052);
342
343}