// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/RivetBoost.hh" #include "Rivet/Projections/DISKinematics.hh" #include "FinalStateBreit.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Math/Constants.hh" #include "Rivet/Tools/ParticleIdUtils.hh" #include "InvariantJetMass.hh" #include "Cparam.hh" #include "Rivet/Projections/Thrust.hh" #include "Broadening.hh" namespace Rivet { ///@brief ZEUS Event shapes in deep inelastic scattering at HERA ///@author Alexander Radovic edit by Jack Pegler ///based on ZEUS collaboration paper DESY-06-042 class ZEUS_2007_S6591329 : public Analysis { public: /// Constructor ZEUS_2007_S6591329() : Analysis("ZEUS_2007_S6591329"){ setBeams(POSITRON,PROTON); } public: void init() { /// Get Projections: const FinalState fs; addProjection(fs,"FS"); const DISKinematics dk; addProjection(dk, "DK"); const FinalStateBreit fsb= FinalStateBreit(dk,1,.25); addProjection(fsb, "FSB"); const InvariantJetMass invariantjetmassb(fsb); addProjection(invariantjetmassb, "IJM"); const Cparam cparamb(fsb); addProjection(cparamb, "CPM"); const Thrust thrustb(fsb); addProjection(thrustb, "Thrust"); const Broadening broadeningb(fsb); addProjection(broadeningb, "Broadening"); ///book histograms: ///DIS Variables _histQ = bookHistogram1D("Q2",100, 0, 20480); _histX = bookHistogram1D("X",100, 0,1); _histY = bookHistogram1D("Y",100, 0,1); ///Invariant Jet Mass _difM.addHistogram(80.0, 160.0, bookHistogram1D(7,1,1)); _difM.addHistogram(160.0, 320.0, bookHistogram1D(8,1,1)); _difM.addHistogram(320.0, 640.0, bookHistogram1D(9,1,1)); _difM.addHistogram(640.0, 1280.0, bookHistogram1D(10,1,1)); _difM.addHistogram(1280.0, 2560.0, bookHistogram1D(11,1,1)); _difM.addHistogram(2560.0, 5120.0, bookHistogram1D(12,1,1)); _difM.addHistogram(5120.0, 10240.0, bookHistogram1D(13,1,1)); _difM.addHistogram(10240.0, 20480.0, bookHistogram1D(14,1,1)); ///C-Paramater _difC.addHistogram(80.0, 160.0, bookHistogram1D(15,1,1)); _difC.addHistogram(160.0, 320.0, bookHistogram1D(16,1,1)); _difC.addHistogram(320.0, 640.0, bookHistogram1D(17,1,1)); _difC.addHistogram(640.0, 1280.0, bookHistogram1D(18,1,1)); _difC.addHistogram(1280.0, 2560.0, bookHistogram1D(19,1,1)); _difC.addHistogram(2560.0, 5120.0, bookHistogram1D(20,1,1)); _difC.addHistogram(5120.0, 10240.0, bookHistogram1D(21,1,1)); _difC.addHistogram(10240.0, 20480.0, bookHistogram1D(22,1,1)); ///Thrust _difT.addHistogram(80.0, 160.0, bookHistogram1D(23,1,1)); _difT.addHistogram(160.0, 320.0, bookHistogram1D(24,1,1)); _difT.addHistogram(320.0, 640.0, bookHistogram1D(25,1,1)); _difT.addHistogram(640.0, 1280.0, bookHistogram1D(26,1,1)); _difT.addHistogram(1280.0, 2560.0, bookHistogram1D(27,1,1)); _difT.addHistogram(2560.0, 5120.0, bookHistogram1D(28,1,1)); _difT.addHistogram(5120.0, 10240.0, bookHistogram1D(29,1,1)); _difT.addHistogram(10240.0, 20480.0, bookHistogram1D(30,1,1)); ///Broadening _difB.addHistogram(80.0, 160.0, bookHistogram1D(31,1,1)); _difB.addHistogram(160.0, 320.0, bookHistogram1D(32,1,1)); _difB.addHistogram(320.0, 640.0, bookHistogram1D(33,1,1)); _difB.addHistogram(640.0, 1280.0, bookHistogram1D(34,1,1)); _difB.addHistogram(1280.0, 2560.0, bookHistogram1D(35,1,1)); _difB.addHistogram(2560.0, 5120.0, bookHistogram1D(36,1,1)); _difB.addHistogram(5120.0, 10240.0, bookHistogram1D(37,1,1)); _difB.addHistogram(10240.0, 20480.0, bookHistogram1D(38,1,1)); ///and profiles: _profQ2vsX = bookProfile1D("Q2vsX",100,0,1); _profQ2vsY = bookProfile1D("Q2vsY",100,0,1); _profYvsX = bookProfile1D("YvsX",100,0,1); } void analyze(const Event& event) { const double weight = event.weight(); ///Apply the Projections needed more than once: const FinalState& fs = applyProjection(event, "FS"); const DISKinematics& dk = applyProjection(event, "DK"); const FinalStateBreit& fsb = applyProjection(event, "FSB"); ///Get final state information: ///Get DIS kinematics: ///Q^2: double q2 = dk.Q2(); ///bjorken x: double x = dk.x(); ///inlasticity: double y = dk.y(); ///W^2 double w2 = dk.W2(); ///Apply other Projections: ///Invariant Jet Mass const InvariantJetMass& invariantjetmassb = applyProjection(event, "IJM"); ///C-Parameter const Cparam& cparamb = applyProjection(event, "CPM"); ///Thrust in current hemisphere of Breit frame const Thrust& thrustb = applyProjection(event, "Thrust"); ///Broadening in current hemisphere of Breit Frame const Broadening& broadeningb = applyProjection(event, "Broadening"); ///before filling, double check we are within Q, x and y limits ///- Should have been done in MCGenerator if(((q2>80)&&(q2<20480))&&((x>.0024)&&(x<.6))&&((y>.04)&&(y<.9))){ _histQ -> fill(q2, weight); _histX -> fill(x, weight); _histY -> fill(y, weight); _difM.fill(160, invariantjetmassb.invariantjetmass(), weight); _difM.fill(320, invariantjetmassb.invariantjetmass(), weight); _difM.fill(640, invariantjetmassb.invariantjetmass(), weight); _difM.fill(1280, invariantjetmassb.invariantjetmass(), weight); _difM.fill(2560, invariantjetmassb.invariantjetmass(), weight); _difM.fill(5120, invariantjetmassb.invariantjetmass(), weight); _difM.fill(10240, invariantjetmassb.invariantjetmass(), weight); _difM.fill(20480, invariantjetmassb.invariantjetmass(), weight); _difC.fill(160, cparamb.cparam(), weight); _difC.fill(320, cparamb.cparam(), weight); _difC.fill(640, cparamb.cparam(), weight); _difC.fill(1280, cparamb.cparam(), weight); _difC.fill(2560, cparamb.cparam(), weight); _difC.fill(5120, cparamb.cparam(), weight); _difC.fill(10240, cparamb.cparam(), weight); _difC.fill(20480, cparamb.cparam(), weight); _difT.fill(160, thrustb.thrust(), weight); _difT.fill(320, thrustb.thrust(), weight); _difT.fill(640, thrustb.thrust(), weight); _difT.fill(1280, thrustb.thrust(), weight); _difT.fill(2560, thrustb.thrust(), weight); _difT.fill(5120, thrustb.thrust(), weight); _difT.fill(10240, thrustb.thrust(), weight); _difT.fill(20480, thrustb.thrust(), weight); _difB.fill(160, broadeningb.broadening(), weight); _difB.fill(320, broadeningb.broadening(), weight); _difB.fill(640, broadeningb.broadening(), weight); _difB.fill(1280, broadeningb.broadening(), weight); _difB.fill(2560, broadeningb.broadening(), weight); _difB.fill(5120, broadeningb.broadening(), weight); _difB.fill(10240, broadeningb.broadening(), weight); _difB.fill(20480, broadeningb.broadening(), weight); ///put data to all the plots that arent for a specific Q2 bin: _profQ2vsX -> fill(x, q2, weight); _profQ2vsY -> fill(y, q2, weight); _profYvsX -> fill(x, y, weight); } } void finalize() { foreach (AIDA::IHistogram1D* hist, _difM.getHistograms()) { normalize(hist); } foreach (AIDA::IHistogram1D* hist, _difC.getHistograms()) { normalize(hist); } foreach (AIDA::IHistogram1D* hist, _difT.getHistograms()) { normalize(hist); } foreach (AIDA::IHistogram1D* hist, _difB.getHistograms()) { normalize(hist); } } private: AIDA::IHistogram1D *_histQ; AIDA::IHistogram1D *_histX; AIDA::IHistogram1D *_histY; AIDA::IProfile1D *_profQ2vsX; AIDA::IProfile1D *_profQ2vsY; AIDA::IProfile1D *_profYvsX; BinnedHistogram _difM; BinnedHistogram _difC; BinnedHistogram _difT; BinnedHistogram _difB; }; // This global object acts as a hook for the plugin system AnalysisBuilder plugin_ZEUS_2007_S6591329; }