1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8	// @brief ZEUS dijet photoproduction in different xgamma regions
  9  
 10	class ZEUS_1997_I450085 : public Analysis {
 11	public:
 12
 13		// Constructor
 14		RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_1997_I450085);
 15		
 16		// Initialization, books projections and histograms
 17		void init() {
 18			
 19			// Projections
 20			// checking recombination scheme and radius checked with original code from M.Wing
 21			FinalState fs;
 22			declare(FastJets(fs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme, 1.0), "Jets");
 23			declare(DISKinematics(), "Kinematics");
 24
 25			//Histograms
 26			// Table 1
 27			book(_h_etabar_all[0], 1, 1, 1);
 28			book(_h_etabar_all[1], 2, 1, 1);
 29			book(_h_etabar_all[2], 3, 1, 1);
 30			book(_h_etabar_all[3], 4, 1, 1);
 31			// Table 6
 32			book(_h_etabar[1][0], 21, 1, 1);
 33			book(_h_etabar[1][1], 22, 1, 1);
 34			book(_h_etabar[1][2], 23, 1, 1);
 35			book(_h_etabar[1][3], 24, 1, 1);
 36			// Table 7
 37			book(_h_etabar[0][0], 25, 1, 1);
 38			book(_h_etabar[0][1], 26, 1, 1);
 39			book(_h_etabar[0][2], 27, 1, 1);
 40			book(_h_etabar[0][3], 28, 1, 1);
 41
 42		}
 43		// Analysis
 44		void analyze(const Event & event) {
 45
 46			// Determine kinematics, including event orientation since ZEUS coord system is for +z = proton direction
 47			const DISKinematics & kin = apply<DISKinematics>(event, "Kinematics");
 48		    if (kin.failed()) vetoEvent;
 49			const int orientation = kin.orientation();
 50
 51			// Q2 and inelasticity cuts
 52			if (kin.Q2() > 4 * GeV2) vetoEvent;
 53			if (!inRange(kin.y(), 0.2, 0.8)) vetoEvent;
 54
 55			// Jet calculation
 56			const Jets jets = apply<FastJets>(event, "Jets") \
 57				.jets(Cuts::Et > 6 * GeV && Cuts::etaIn(-1.375 * orientation, 1.875 * orientation), cmpMomByEt);
 58			MSG_DEBUG("Jet multiplicity = " << jets.size());
 59			//Dijet event selection
 60			if (jets.size() < 2) vetoEvent;
 61			const Jet & j1 = jets[0];
 62		    const Jet & j2 = jets[1];
 63	
 64			//Jet eta, average eta and eta difference calculation
 65			const double eta1 = orientation * j1.eta(), eta2 = orientation * j2.eta();
 66			const double etabar = (eta1 + eta2) / 2;
 67			const double etadiff = eta1 - eta2;
 68
 69			//Cut in pseudorapidity difference
 70			if (abs(etadiff) > 0.5) vetoEvent;
 71
 72			// Calculation of x_gamma^obs
 73			/// note Assuming Ee is the lab frame positron momentum, not in proton rest frame cf.
 74		    const double xyobs = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2 * kin.y() * kin.beamLepton().E());
 75			const size_t i_xyobs = (xyobs < 0.75) ? 0 : 1;
 76
 77			//Classify events according to minimum jet energies
 78			size_t iE_min = 0;
 79			if (j1.Et() > 8 * GeV && j2.Et() > 8 * GeV)
 80				iE_min = 1;
 81			if (j1.Et() > 11 * GeV && j2.Et() > 11 * GeV)
 82				iE_min = 2;
 83			if (j1.Et() > 15 * GeV && j2.Et() > 15 * GeV)
 84				iE_min = 3;
 85
 86			//Fill histograms
 87			for (size_t isel = 0; isel <= iE_min; ++isel) {
 88
 89				//T1
 90				_h_etabar_all[isel]->fill(etabar);
 91				//T6, T7
 92				if (xyobs < 0.3) vetoEvent;
 93				_h_etabar[i_xyobs][isel]->fill(etabar);
 94			}
 95			
 96
 97		}
 98		
 99		// Finalize
100		void finalize() {
101			const double sf = crossSection() / nanobarn / sumOfWeights();
102			for (size_t ix = 0; ix < 2; ++ix) {
103				for (auto& h : _h_etabar[ix]) scale(h, sf);
104
105			}
106			for (auto& h : _h_etabar_all) scale(h, sf);
107
108		}
109
110
111	private:
112		//name Histograms
113		Histo1DPtr _h_etabar_all[4], _h_etabar[2][4];
114
115	};
116	
117	RIVET_DECLARE_PLUGIN(ZEUS_1997_I450085);
118
119}