1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3// #include "Rivet/Projections/ChargedFinalState.hh"
 4
 5namespace Rivet {
 6
 7  /// Generic analysis looking at various distributions of final state particles
 8  class MC_PDFS : public Analysis {
 9  public:
10
11    /// Constructor
12    MC_PDFS()
13      : Analysis("MC_PDFS")
14    {    }
15
16
17  public:
18
19    /// @name Analysis methods
20    /// @{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24      // Projections
25      // declare(ChargedFinalState((Cuts::etaIn(-5.0, 5.0) && Cuts::pT >=  500*MeV)), "CFS");
26
27      // Histograms
28      book(_histPdfX ,"PdfX", logspace(50, 0.000001, 1.0));
29      book(_histPdfXmin ,"PdfXmin", logspace(50, 0.000001, 1.0));
30      book(_histPdfXmax ,"PdfXmax", logspace(50, 0.000001, 1.0));
31      book(_histPdfQ ,"PdfQ", 50, 0.0, 30.0);
32      book(_histPdfXQ,"PdfXQ", logspace(50, 0.000001, 1.0), linspace(50, 0.0, 30.0));
33      //book( _histPdfTrackptVsX ,"PdfTrackptVsX", logspace(50, 0.000001, 1.0));
34      //book( _histPdfTrackptVsQ ,"PdfTrackptVsQ", 50, 0.0, 30.0);
35    }
36
37
38    /// Perform the per-event analysis
39    void analyze(const Event& event) {
40
41      // This analysis needs a valid HepMC PDF info object to do anything
42      if (event.genEvent()->pdf_info() == 0) vetoEvent;
43      PdfInfo pdfi = *(event.genEvent()->pdf_info());
44
45      MSG_DEBUG("PDF Q = " << pdfi.scale<< " for (id, x) = "
46                << "(" << pdfi.parton_id[0] << ", " << pdfi.x[0] << ") "
47                << "(" << pdfi.parton_id[1] << ", " << pdfi.x[1] << ")");
48      _histPdfX->fill(pdfi.x[0]);
49      _histPdfX->fill(pdfi.x[1]);
50      _histPdfXmin->fill(std::min(pdfi.x[0], pdfi.x[1]));
51      _histPdfXmax->fill(std::max(pdfi.x[0], pdfi.x[1]));
52      _histPdfQ->fill(pdfi.scale); // always in GeV?
53      _histPdfXQ->fill(pdfi.x[0], pdfi.scale); // always in GeV?
54      _histPdfXQ->fill(pdfi.x[1], pdfi.scale); // always in GeV?
55
56      // const FinalState& cfs = apply<FinalState>(event, "CFS");
57      // for (const Particle& p : cfs.particles()) {
58      //   if (fabs(eta) < 2.5 && p.pT() > 10*GeV) {
59      //     _histPdfTrackptVsX->fill(pdfi.x1(), p.pT()/GeV);
60      //     _histPdfTrackptVsX->fill(pdfi.x2(), p.pT()/GeV);
61      //     _histPdfTrackptVsQ->fill(pdfi.scalePDF(), p.pT()/GeV);
62      //   }
63      // }
64
65    }
66
67
68
69    /// Finalize
70    void finalize() {
71      scale(_histPdfX, 1/sumOfWeights());
72      scale(_histPdfXmin, 1/sumOfWeights());
73      scale(_histPdfXmax, 1/sumOfWeights());
74      scale(_histPdfQ, 1/sumOfWeights());
75    }
76
77    /// @}
78
79
80  private:
81
82    /// @name Histograms
83    /// @{
84    Histo1DPtr _histPdfX, _histPdfXmin, _histPdfXmax, _histPdfQ;
85    Histo2DPtr _histPdfXQ;
86    // Profile1DPtr   _histPdfTrackptVsX, _histPdfTrackptVsQ;
87    /// @}
88
89  };
90
91
92  RIVET_DECLARE_PLUGIN(MC_PDFS);
93
94}