1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/DISFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11  namespace H1_2007_I746380_PROJECTIONS {
 12
 13
 14    /// Projection to find the largest gaps and the masses of the two
 15    /// systems separated by the gap. Based on the HZTools gap-finding
 16    /// method (hzhadgap.F). Note that gaps are found in the HCM frame.
 17    ///
 18    /// @author Christine O. Rasmussen.
 19    class RapidityGap : public Projection {
 20    public:
 21
 22      /// Type of DIS boost to apply
 23      enum Frame { HCM, LAB, XCM };
 24
 25      RapidityGap() {
 26        setName("RapidityGap");
 27        declare(DISKinematics(), "DISKIN");
 28        declare(DISFinalState(DISFrame::HCM), "DISFS");
 29      }
 30
 31      RIVET_DEFAULT_PROJ_CLONE(RapidityGap);
 32
 33      /// Import to avoid warnings about overload-hiding
 34      using Projection::operator =;
 35
 36      double M2X() const { return _M2X; }
 37      double M2Y() const { return _M2Y; }
 38      double t() const { return _t; }
 39      double gap() const { return _gap; }
 40      double gapUpp() const { return _gapUpp; }
 41      double gapLow() const { return _gapLow; }
 42      double EpPzX(Frame f) const {
 43        if (f == LAB) return _ePpzX_LAB;
 44        else if (f == XCM) return _ePpzX_XCM;
 45        else return _ePpzX_HCM;
 46      }
 47      double EmPzX(Frame f) const {
 48        if (f == LAB) return _eMpzX_LAB;
 49        else if (f == XCM) return _eMpzX_XCM;
 50        else return _eMpzX_HCM;
 51      }
 52      FourMomentum pX(Frame f) const {
 53        if (f == LAB) return _momX_LAB;
 54        else if (f == XCM) return _momX_XCM;
 55        else return _momX_HCM;
 56      }
 57      FourMomentum pY(Frame f) const {
 58        if (f == LAB) return _momY_LAB;
 59        else if (f == XCM) return _momY_XCM;
 60        else return _momY_HCM;
 61      }
 62      const Particles& systemX(Frame f) const {
 63        if (f == LAB) return _pX_LAB;
 64        else if (f == XCM) return _pX_XCM;
 65        else return _pX_HCM;
 66      }
 67      const Particles& systemY(Frame f) const {
 68        if (f == LAB) return _pY_LAB;
 69        else if (f == XCM) return _pY_XCM;
 70        else return _pY_HCM;
 71      }
 72
 73    protected:
 74
 75      virtual CmpState compare(const Projection& p) const {
 76        const RapidityGap& other = pcast<RapidityGap>(p);
 77        return mkNamedPCmp(other, "DISKIN") || mkNamedPCmp(other, "DISFS");
 78      }
 79
 80      virtual void project(const Event& e){
 81        const DISKinematics& dk = apply<DISKinematics>(e, "DISKIN");
 82        const Particles& p      = apply<DISFinalState>(e, "DISFS").particles(cmpMomByEta);
 83        findgap(p, dk);
 84      }
 85
 86      void clearAll(){
 87        _M2X = _M2Y = _t = _gap = 0.;
 88        _gapUpp = _gapLow = -8.;
 89        _ePpzX_HCM = _eMpzX_HCM =_ePpzX_LAB = _eMpzX_LAB = _ePpzX_XCM = _eMpzX_XCM = 0.;
 90        _momX_HCM.setPE(0., 0., 0., 0.);
 91        _momY_HCM.setPE(0., 0., 0., 0.);
 92        _momX_XCM.setPE(0., 0., 0., 0.);
 93        _momY_XCM.setPE(0., 0., 0., 0.);
 94        _momX_LAB.setPE(0., 0., 0., 0.);
 95        _momY_LAB.setPE(0., 0., 0., 0.);
 96        _pX_HCM.clear();
 97        _pY_HCM.clear();
 98        _pX_XCM.clear();
 99        _pY_XCM.clear();
100        _pX_LAB.clear();
101        _pY_LAB.clear();
102      }
103
104      void findgap(const Particles& particles, const DISKinematics& diskin){
105
106        clearAll();
107
108        // Begin by finding largest gap and gapedges between all final
109        // state particles in HCM frame.
110        int nP  = particles.size();
111        int dir = diskin.orientation();
112        for (int i = 0; i < nP-1; ++i){
113          double tmpGap = abs(particles[i+1].eta() - particles[i].eta());
114          if (tmpGap > _gap) {
115            _gap    = tmpGap;
116            _gapLow = (dir > 0) ? particles[i].eta() : dir * particles[i+1].eta();
117            _gapUpp = (dir > 0) ? particles[i+1].eta() : dir * particles[i].eta();
118          }
119        }
120
121        // Define the two systems X and Y.
122        Particles tmp_pX, tmp_pY;
123        for (const Particle& ip : particles) {
124          if (dir * ip.eta() > _gapLow) tmp_pX.push_back(ip);
125          else tmp_pY.push_back(ip);
126        }
127
128        Particles pX, pY;
129        pX = (dir < 0) ? tmp_pY : tmp_pX;
130        pY = (dir < 0) ? tmp_pX : tmp_pY;
131
132        // Find variables related to HCM frame.
133        // Note that HCM has photon along +z, as opposed to
134        // H1 where proton is along +z. This results in a sign change
135        // as compared to H1 papers!
136
137        // X - side
138        FourMomentum momX;
139        for (const Particle& jp : pX) {
140          momX  += jp.momentum();
141          _ePpzX_HCM += jp.E() - jp.pz(); // Sign + => -
142          _eMpzX_HCM += jp.E() + jp.pz(); // Sign - => +
143        }
144        _momX_HCM = momX;
145        _pX_HCM   = pX;
146        _M2X      = _momX_HCM.mass2();
147
148        // Y - side
149        FourMomentum momY;
150        for (const Particle& kp : pY) momY += kp.momentum();
151        _momY_HCM = momY;
152        _pY_HCM   = pY;
153        _M2Y      = _momY_HCM.mass2();
154
155        // Find variables related to LAB frame
156        const LorentzTransform hcmboost   = diskin.boostHCM();
157        const LorentzTransform hcminverse = hcmboost.inverse();
158        _momX_LAB = hcminverse.transform(_momX_HCM);
159        _momY_LAB = hcminverse.transform(_momY_HCM);
160
161        // Find momenta in XCM frame. Note that it is HCM frame that is
162        // boosted, resulting in a sign change later!
163        const bool doXCM = (momX.betaVec().mod2() < 1.);
164        if (doXCM) {
165          const LorentzTransform xcmboost =
166            LorentzTransform::mkFrameTransformFromBeta(momX.betaVec());
167          _momX_XCM = xcmboost.transform(momX);
168          _momY_XCM = xcmboost.transform(momY);
169        }
170
171        for (const Particle& jp : pX) {
172          // Boost from HCM to LAB.
173          FourMomentum lab = hcminverse.transform(jp.momentum());
174          _ePpzX_LAB += lab.E() + dir * lab.pz();
175          _eMpzX_LAB += lab.E() - dir * lab.pz();
176          Particle plab = jp;
177          plab.setMomentum(lab);
178          _pX_LAB.push_back(plab);
179          // Set XCM. Note that since HCM frame is boosted to XCM frame,
180          // we have a sign change
181          if (doXCM) {
182            const LorentzTransform xcmboost =
183              LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
184            FourMomentum xcm = xcmboost.transform(jp.momentum());
185            _ePpzX_XCM += xcm.E() - xcm.pz(); // Sign + => -
186            _eMpzX_XCM += xcm.E() + xcm.pz(); // Sign - => +
187            Particle pxcm = jp;
188            pxcm.setMomentum(xcm);
189            _pX_XCM.push_back(pxcm);
190          }
191        }
192
193        for (const Particle& jp : pY) {
194          // Boost from HCM to LAB
195          FourMomentum lab = hcminverse.transform(jp.momentum());
196          Particle plab = jp;
197          plab.setMomentum(lab);
198          _pY_LAB.push_back(plab);
199          // Boost from HCM to XCM
200          if (doXCM) {
201            const LorentzTransform xcmboost =
202              LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
203            FourMomentum xcm = xcmboost.transform(jp.momentum());
204            Particle pxcm = jp;
205            pxcm.setMomentum(xcm);
206            _pY_XCM.push_back(pxcm);
207          }
208        }
209
210        // Boost p from LAB to HCM frame to find t.
211        ///
212        /// @note Currently can only handle gap on proton side.
213        ///
214        /// @todo Expand to also handle gap on photon side
215        const FourMomentum proton = hcmboost.transform(diskin.beamHadron().momentum());
216        FourMomentum pPom         = proton - _momY_HCM;
217        _t                        = pPom * pPom;
218
219      }
220
221    private:
222
223      double _M2X, _M2Y, _t;
224      double _gap, _gapUpp, _gapLow;
225      double _ePpzX_LAB, _eMpzX_LAB, _ePpzX_HCM, _eMpzX_HCM, _ePpzX_XCM, _eMpzX_XCM;
226      FourMomentum _momX_HCM, _momY_HCM,_momX_LAB, _momY_LAB, _momX_XCM, _momY_XCM;
227      Particles _pX_HCM, _pY_HCM, _pX_LAB, _pY_LAB, _pX_XCM, _pY_XCM;
228
229    };
230
231    /// Projection to boost system X (photon+Pomeron) particles into its rest frame.
232    ///
233    /// @author Ilkka Helenius
234    class BoostedXSystem : public FinalState {
235    public:
236
237      BoostedXSystem(const FinalState& fs) {
238        setName("BoostedXSystem");
239        declare(fs,"FS");
240        declare(RapidityGap(), "RAPGAP");
241      }
242
243      // Return the boost to XCM frame.
244      const LorentzTransform& boost() const { return _boost; }
245
246      RIVET_DEFAULT_PROJ_CLONE(BoostedXSystem);
247
248      /// Import to avoid warnings about overload-hiding
249      using Projection::operator =;
250
251    protected:
252
253      // Apply the projection on the supplied event.
254      void project(const Event& e){
255
256        const RapidityGap& rg = apply<RapidityGap>(e, "RAPGAP");
257
258        // Total momentum of the system X.
259        const FourMomentum pX = rg.pX(RapidityGap::HCM);
260
261        // Reset the boost. Is there a separate method for this?
262        _boost = combine(_boost, _boost.inverse());
263
264        // Define boost only when numerically safe, otherwise negligible.
265        if (pX.betaVec().mod2() < 1.)
266          _boost = LorentzTransform::mkFrameTransformFromBeta(pX.betaVec());
267
268        // Boost the particles from system X.
269        _theParticles.clear();
270        _theParticles.reserve(rg.systemX(RapidityGap::HCM).size());
271        for (const Particle& p : rg.systemX(RapidityGap::HCM)) {
272          Particle temp = p;
273          temp.setMomentum(_boost.transform(temp.momentum()));
274          _theParticles.push_back(temp);
275        }
276
277      }
278
279      // Compare projections.
280      CmpState compare(const Projection& p) const {
281        const BoostedXSystem& other = pcast<BoostedXSystem>(p);
282        return mkNamedPCmp(other, "RAPGAP") || mkNamedPCmp(other, "FS");
283      }
284
285    private:
286
287      LorentzTransform _boost;
288
289    };
290
291  }
292
293
294
295  /// @brief H1 diffractive dijets
296  ///
297  /// Diffractive dijets H1 with 920 GeV p and 27.5 GeV e
298  /// Note tagged protons!
299  ///
300  /// @author Christine O. Rasmussen
301  class H1_2007_I746380 : public Analysis {
302  public:
303
304    typedef H1_2007_I746380_PROJECTIONS::RapidityGap RapidityGap;
305    typedef H1_2007_I746380_PROJECTIONS::BoostedXSystem BoostedXSystem;
306
307    /// Constructor
308    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I746380);
309
310    /// @name Analysis methods
311    /// @{
312
313    // Book projections and histograms
314    void init() {
315
316      declare(DISKinematics(), "Kinematics");
317      const DISFinalState& disfs = declare(DISFinalState(DISFrame::HCM), "DISFS");
318      const BoostedXSystem& disfsXcm = declare( BoostedXSystem(disfs), "BoostedXFS");
319      declare(FastJets(disfsXcm, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
320                       JetMuons::ALL, JetInvisibles::NONE, nullptr), "DISFSJets");
321      declare(RapidityGap(), "RapidityGap");
322
323      // Book histograms from REF data
324      book(_h_DIS_dsigdzPom, 1, 1, 1);
325      book(_h_DIS_dsigdlogXpom, 2, 1, 1);
326      book(_h_DIS_dsigdW, 3, 1, 1);
327      book(_h_DIS_dsigdQ2, 4, 1, 1);
328      book(_h_DIS_dsigdEtJet1, 5, 1, 1);
329      book(_h_DIS_dsigdAvgEta, 6, 1, 1);
330      book(_h_DIS_dsigdDeltaEta, 7, 1, 1);
331
332      book(_h_PHO_dsigdzPom, 8, 1, 1);
333      book(_h_PHO_dsigdxGam, 9, 1, 1);
334      book(_h_PHO_dsigdlogXpom, 10, 1, 1);
335      book(_h_PHO_dsigdW, 11, 1, 1);
336      book(_h_PHO_dsigdEtJet1, 12, 1, 1);
337      book(_h_PHO_dsigdAvgEta, 13, 1, 1);
338      book(_h_PHO_dsigdDeltaEta, 14, 1, 1);
339      book(_h_PHO_dsigdMjets, 15, 1, 1);
340
341      isDIS  = false;
342      nVeto0 = 0;
343      nVeto1 = 0;
344      nVeto2 = 0;
345      nVeto3 = 0;
346      nVeto4 = 0;
347      nVeto5 = 0;
348      nPHO   = 0;
349      nDIS   = 0;
350    }
351
352    // Do the analysis
353    void analyze(const Event& event) {
354
355      // Event weight
356      isDIS  = false;
357
358      // Projections - special handling of events where no proton found:
359      const RapidityGap&    rg = apply<RapidityGap>(event, "RapidityGap");
360      const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
361      const BoostedXSystem& disfsXcm = apply<BoostedXSystem>( event, "BoostedXFS");
362
363      // Determine kinematics: H1 has +z = proton direction
364      int dir   = kin.orientation();
365      double W2 = kin.W2();
366      double W  = sqrt(W2);
367      double y  = kin.y();
368      double Q2 = kin.Q2();
369
370      // Separate into DIS and PHO regimes else veto
371      if (!inRange(W, 165.*GeV, 242.*GeV)) vetoEvent;
372      if (Q2 < 0.01*GeV2) {
373        isDIS = false;
374        ++nPHO;
375      } else if (inRange(Q2, 4.0*GeV2, 80.*GeV2)) {
376        isDIS = true;
377        ++nDIS;
378      } else {
379        vetoEvent;
380      }
381      ++nVeto0;
382
383      // Find diffractive variables as defined in paper.
384      const double M2Y  = rg.M2Y();
385      const double M2X  = rg.M2X();
386      const double abst = abs(rg.t());
387      const double xPom = (isDIS) ? (Q2 + M2X) / (Q2 + W2) :
388        rg.EpPzX(RapidityGap::LAB) / (2. * kin.beamHadron().E());
389
390      // Veto if outside allowed region
391      if (sqrt(M2Y) > 1.6*GeV)    vetoEvent;
392      ++nVeto1;
393      if (abst > 1.0*GeV2) vetoEvent;
394      ++nVeto2;
395      if (xPom > 0.03)     vetoEvent;
396      ++nVeto3;
397
398      // Jet selection. Note jets are found in photon-proton (XCM)
399      // frame, but eta cut is applied in lab frame!
400      Cut jetcuts = Cuts::Et > 4.* GeV;
401      Jets jets   = apply<FastJets>(event, "DISFSJets").jets(jetcuts, cmpMomByEt);
402      // Veto if not dijets and if Et_j1 < 5.0
403      if (jets.size() < 2)       vetoEvent;
404      if (jets[0].Et() < 5.*GeV) vetoEvent;
405      ++nVeto4;
406      // Find Et_jet1 and deltaEta* in XCM frame
407      double EtJet1       = jets[0].Et() * GeV;
408      double etaXCMJet1   = jets[0].eta();
409      double etaXCMJet2   = jets[1].eta();
410      double deltaEtaJets = abs(etaXCMJet1 - etaXCMJet2);
411
412      // Transform from XCM to HCM
413      const LorentzTransform xcmboost = disfsXcm.boost();
414      for (int i = 0; i < 2; ++i) jets[i].transformBy(xcmboost.inverse());
415      // Find mass of jets and EpPz, EmPz of jets
416      FourMomentum momJets = jets[0].momentum() + jets[1].momentum();
417      double M2jets   = momJets.mass2();
418      double EpPzJets = 0.;
419      double EmPzJets = 0.;
420      // DIS variables are found in XCM frame, so boost back again
421      if (isDIS){
422        for (int i = 0; i < 2; ++i) jets[i].transformBy(xcmboost);
423      }
424      // Note sign change wrt. H1 because photon is in +z direction
425      // Jets in HCM so no need to consider orientation.
426      for (int i = 0; i < 2; ++i){
427        EpPzJets += jets[i].E() - jets[i].pz(); // Sign: + => -
428        EmPzJets += jets[i].E() + jets[i].pz(); // Sign: - => +
429      }
430
431      // Transform the jets from HCM to LAB frame where eta cut is
432      // applied for photoproduction.
433      const LorentzTransform hcmboost = kin.boostHCM();
434      for (int i = 0; i < 2; ++i) jets[i].transformBy(hcmboost.inverse());
435      double etaLabJet1 = dir * jets[0].eta();
436      double etaLabJet2 = dir * jets[1].eta();
437      double etaMin     = (isDIS) ? -3. : -1.;
438      double etaMax     = (isDIS) ? 0. : 2.;
439      double eta1       = (isDIS) ? etaXCMJet1 : etaLabJet1;
440      double eta2       = (isDIS) ? etaXCMJet2 : etaLabJet2;
441      if (!inRange(eta1, etaMin, etaMax)) vetoEvent;
442      if (!inRange(eta2, etaMin, etaMax)) vetoEvent;
443      ++nVeto5;
444
445      // Pseudorapidity distributions are examined in lab frame:
446      double avgEtaJets   = 0.5 * (etaLabJet1 + etaLabJet2);
447
448      // Derive xPom and xGam values from the jet kinematics.
449      double zPomJets, xGamJets;
450      if (isDIS) {
451        zPomJets = (Q2 + M2jets) / (Q2 + M2X);
452        xGamJets = EmPzJets / rg.EmPzX(RapidityGap::XCM);
453      } else {
454        // Boost E_p, E_e to HCM frame
455        FourMomentum lep = hcmboost.transform(kin.beamLepton().momentum());
456        FourMomentum had = hcmboost.transform(kin.beamHadron().momentum());
457        zPomJets = EpPzJets / (2. * xPom * had.E());
458        xGamJets = EmPzJets / (2. * y * lep.E());
459      }
460
461      // Now fill histograms
462      if (isDIS){
463        _h_DIS_dsigdzPom     ->fill(zPomJets);
464        _h_DIS_dsigdlogXpom  ->fill(log10(xPom));
465        _h_DIS_dsigdW        ->fill(W);
466        _h_DIS_dsigdQ2       ->fill(Q2);
467        _h_DIS_dsigdEtJet1   ->fill(EtJet1);
468        _h_DIS_dsigdAvgEta   ->fill(avgEtaJets);
469        _h_DIS_dsigdDeltaEta ->fill(deltaEtaJets);
470      } else {
471        _h_PHO_dsigdzPom     ->fill(zPomJets);
472        _h_PHO_dsigdxGam     ->fill(xGamJets);
473        _h_PHO_dsigdlogXpom  ->fill(log10(xPom));
474        _h_PHO_dsigdW        ->fill(W);
475        _h_PHO_dsigdEtJet1   ->fill(EtJet1);
476        _h_PHO_dsigdAvgEta   ->fill(avgEtaJets);
477        _h_PHO_dsigdDeltaEta ->fill(deltaEtaJets);
478        _h_PHO_dsigdMjets    ->fill(sqrt(M2jets));
479      }
480
481    }
482
483    // Finalize
484    void finalize() {
485      // Normalise to cross section
486      const double norm = crossSection()/picobarn/sumOfWeights();
487
488      scale( _h_DIS_dsigdzPom    , norm);
489      scale( _h_DIS_dsigdlogXpom , norm);
490      scale( _h_DIS_dsigdW       , norm);
491      scale( _h_DIS_dsigdQ2      , norm);
492      scale( _h_DIS_dsigdEtJet1  , norm);
493      scale( _h_DIS_dsigdAvgEta  , norm);
494      scale( _h_DIS_dsigdDeltaEta, norm);
495
496      scale( _h_PHO_dsigdzPom    , norm);
497      scale( _h_PHO_dsigdxGam    , norm);
498      scale( _h_PHO_dsigdlogXpom , norm);
499      scale( _h_PHO_dsigdW       , norm);
500      scale( _h_PHO_dsigdEtJet1  , norm);
501      scale( _h_PHO_dsigdAvgEta  , norm);
502      scale( _h_PHO_dsigdDeltaEta, norm);
503      scale( _h_PHO_dsigdMjets   , norm);
504
505      const double dPHO = nPHO;
506      MSG_INFO("H1_2007_I746380");
507      MSG_INFO("Cross section = " << crossSection()/picobarn << " pb");
508      MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
509      MSG_INFO("Number of PHO = " << nPHO << ", number of DIS = " << nDIS);
510      MSG_INFO("Events passing electron veto   = " << nVeto0 << " (" << nVeto0/dPHO * 100. << "%)" );
511      MSG_INFO("Events passing MY              = " << nVeto1 << " (" << nVeto1/dPHO * 100. << "%)" );
512      MSG_INFO("Events passing t veto          = " << nVeto2 << " (" << nVeto2/dPHO * 100. << "%)" );
513      MSG_INFO("Events passing xPom            = " << nVeto3 << " (" << nVeto3/dPHO * 100. << "%)" );
514      MSG_INFO("Events passing jet Et veto     = " << nVeto4 << " (" << nVeto4/dPHO * 100. << "%)" );
515      MSG_INFO("Events passing jet eta veto    = " << nVeto5 << " (" << nVeto5/dPHO * 100. << "%)" );
516
517    }
518
519    /// @}
520
521
522  private:
523
524    /// @name Histograms
525    /// @{
526    // Book histograms from REF data
527    Histo1DPtr _h_DIS_dsigdzPom    ;
528    Histo1DPtr _h_DIS_dsigdlogXpom ;
529    Histo1DPtr _h_DIS_dsigdW       ;
530    Histo1DPtr _h_DIS_dsigdQ2      ;
531    Histo1DPtr _h_DIS_dsigdEtJet1  ;
532    Histo1DPtr _h_DIS_dsigdAvgEta  ;
533    Histo1DPtr _h_DIS_dsigdDeltaEta;
534
535    Histo1DPtr _h_PHO_dsigdzPom    ;
536    Histo1DPtr _h_PHO_dsigdxGam    ;
537    Histo1DPtr _h_PHO_dsigdlogXpom ;
538    Histo1DPtr _h_PHO_dsigdW       ;
539    Histo1DPtr _h_PHO_dsigdEtJet1  ;
540    Histo1DPtr _h_PHO_dsigdAvgEta  ;
541    Histo1DPtr _h_PHO_dsigdDeltaEta;
542    Histo1DPtr _h_PHO_dsigdMjets   ;
543    /// @}
544
545    bool isDIS;
546    int  nVeto0, nVeto1, nVeto2, nVeto3, nVeto4, nVeto5;
547    int nPHO, nDIS;
548  };
549
550  RIVET_DECLARE_PLUGIN(H1_2007_I746380);
551
552}