[Rivet-svn] r3320 - in trunk: data/plotinfo src/Analyses

[blackhole at projects.hepforge.org]

Author: buckley
Date: Tue Aug 23 12:30:31 2011
New Revision: 3320

Log:
Adding plotting of some PDF quantities if PdfInfo objects are available in the events

Modified:
   trunk/data/plotinfo/MC_GENERIC.plot
   trunk/src/Analyses/MC_GENERIC.cc

Modified: trunk/data/plotinfo/MC_GENERIC.plot
==============================================================================
--- trunk/data/plotinfo/MC_GENERIC.plot	Tue Aug 23 12:20:37 2011	(r3319)
+++ trunk/data/plotinfo/MC_GENERIC.plot	Tue Aug 23 12:30:31 2011	(r3320)
@@ -168,3 +168,44 @@
 LogY=0
 LegendYPos=0.5
 # END PLOT
+
+# BEGIN PLOT /MC_GENERIC/PdfX
+Title=PDF sampling $x$ value distribution
+XLabel=$x$
+YLabel=$1/N_\text{evt} \  \mathrm{d}N_\text{evt}/\mathrm{d}x$
+LogX=1
+LogY=1
+FullRange=1
+LegendXPos=0.1
+LegendYPos=0.5
+# END PLOT
+
+# BEGIN PLOT /MC_GENERIC/PdfXmin
+Title=PDF sampling $x_\text{min}$ value distribution
+XLabel=$x_\text{min}$
+# END PLOT
+
+# BEGIN PLOT /MC_GENERIC/PdfXmax
+Title=PDF sampling $x_\text{max}$ value distribution
+XLabel=$x_\text{max}$
+LegendXPos=0.3
+XMin=0.00005
+# END PLOT
+
+# BEGIN PLOT /MC_GENERIC/PdfQ
+Title=PDF sampling $Q$ value distribution
+XLabel=$Q$
+YLabel=$1/N_\text{evt} \  \mathrm{d}N_\text{evt}/\mathrm{d}Q$
+LogX=0
+LogY=1
+FullRange=1
+LegendYPos=0.9
+# END PLOT
+
+# # BEGIN PLOT /MC_GENERIC/PdfXQ
+# Title=PDF sampling $x, Q$ value distribution
+# XLabel=$x$
+# YLabel=$Q$
+# LogX=0
+# LogY=0
+# # END PLOT

Modified: trunk/src/Analyses/MC_GENERIC.cc
==============================================================================
--- trunk/src/Analyses/MC_GENERIC.cc	Tue Aug 23 12:20:37 2011	(r3319)
+++ trunk/src/Analyses/MC_GENERIC.cc	Tue Aug 23 12:30:31 2011	(r3320)
@@ -72,6 +72,14 @@
 
       _histPhi    = bookHistogram1D("Phi", 50, 0, TWOPI);
       _histPhiCh  = bookHistogram1D("PhiCh", 50, 0, TWOPI);
+
+      _histPdfX = bookHistogram1D("PdfX", logspace(0.000001, 1.0, 50));
+      _histPdfXmin = bookHistogram1D("PdfXmin", logspace(0.000001, 1.0, 50));
+      _histPdfXmax = bookHistogram1D("PdfXmax", logspace(0.000001, 1.0, 50));
+      _histPdfQ = bookHistogram1D("PdfQ", 50, 0.0, 30.0);
+      // _histPdfXQ = bookHistogram2D("PdfXQ", logspace(0.000001, 1.0, 50), linspace(0.0, 30.0, 50));
+      // _histPdfTrackptVsX = bookProfile1D("PdfTrackptVsX", logspace(0.000001, 1.0, 50));
+      // _histPdfTrackptVsQ = bookProfile1D("PdfTrackptVsQ", 50, 0.0, 30.0);
     }
 
 
@@ -85,6 +93,19 @@
         _histAllPIDs->fill(abs(gp->pdg_id()), weight);
       }
 
+      // Print and plot PDF info
+      if (event.genEvent().pdf_info() != 0) {
+        HepMC::PdfInfo pdfi = *event.genEvent().pdf_info();
+        MSG_DEBUG("PDF Q = " << pdfi.scalePDF() << " for (id, x) = "
+                  << "(" << pdfi.id1() << ", " << pdfi.x1() << ") "
+                  << "(" << pdfi.id2() << ", " << pdfi.x2() << ")");
+        _histPdfX->fill(pdfi.x1(), weight);
+        _histPdfX->fill(pdfi.x2(), weight);
+        _histPdfXmin->fill(std::min(pdfi.x1(), pdfi.x2()), weight);
+        _histPdfXmax->fill(std::max(pdfi.x1(), pdfi.x2()), weight);
+        _histPdfQ->fill(pdfi.scalePDF(), weight); // always in GeV?
+      }
+
       // Charged + neutral final state
       const FinalState& cnfs = applyProjection<FinalState>(event, "FS");
       MSG_DEBUG("Total multiplicity = " << cnfs.size());
@@ -132,6 +153,15 @@
         _histPtCh->fill(p.momentum().pT()/GeV, weight);
         _histECh->fill(p.momentum().E()/GeV, weight);
         _histPhiCh->fill(p.momentum().phi(), weight);
+
+        // if (event.genEvent().pdf_info() != 0) {
+        //   if (fabs(eta) < 2.5 && p.momentum().pT() > 10*GeV) {
+        //     HepMC::PdfInfo pdfi = *event.genEvent().pdf_info();
+        //     _histPdfTrackptVsX->fill(pdfi.x1(), p.momentum().pT()/GeV, weight);
+        //     _histPdfTrackptVsX->fill(pdfi.x2(), p.momentum().pT()/GeV, weight);
+        //     _histPdfTrackptVsQ->fill(pdfi.scalePDF(), p.momentum().pT()/GeV, weight);
+        //   }
+        // }
       }
 
 
@@ -180,6 +210,11 @@
       scale(_histPhi, 1/sumOfWeights());
       scale(_histPhiCh, 1/sumOfWeights());
 
+      scale(_histPdfX, 1/sumOfWeights());
+      scale(_histPdfXmin, 1/sumOfWeights());
+      scale(_histPdfXmax, 1/sumOfWeights());
+      scale(_histPdfQ, 1/sumOfWeights());
+
       histogramFactory().divide(histoPath("EtaPMRatio"), *_tmphistEtaPlus, *_tmphistEtaMinus);
       histogramFactory().divide(histoPath("EtaChPMRatio"), *_tmphistEtaChPlus, *_tmphistEtaChMinus);
       histogramFactory().divide(histoPath("RapidityPMRatio"), *_tmphistRapPlus, *_tmphistRapMinus);
@@ -197,8 +232,8 @@
     shared_ptr<LWH::Histogram1D> _tmphistRapPlus, _tmphistRapMinus;
     shared_ptr<LWH::Histogram1D> _tmphistRapChPlus, _tmphistRapChMinus;
 
+    /// @name Histograms
     //@{
-    /// Histograms
     AIDA::IHistogram1D *_histMult, *_histMultCh;
     AIDA::IHistogram1D *_histStablePIDs, *_histDecayedPIDs, *_histAllPIDs;
     AIDA::IHistogram1D *_histEtaPi, *_histEtaK, *_histEtaLambda;
@@ -208,6 +243,8 @@
     AIDA::IHistogram1D *_histPt, *_histPtCh;
     AIDA::IHistogram1D *_histE, *_histECh;
     AIDA::IHistogram1D *_histPhi, *_histPhiCh;
+    AIDA::IHistogram1D *_histPdfX, *_histPdfXmin, *_histPdfXmax, *_histPdfQ;
+    // AIDA::IProfile1D   *_histPdfTrackptVsX, *_histPdfTrackptVsQ;
     //@}
 
   };