Merge pull request #31639 from CTPPS/pps_xangle_beta_distributions

PPS: xangle-beta* distributions

CalibPPS/ESProducers/plugins/CTPPSLHCInfoRandomXangleESSource.cc

@@ -16,7 +16,7 @@
 #include "CLHEP/Random/JamesRandom.h"
 
 #include "TFile.h"
-#include "TH1D.h"
+#include "TH2D.h"
 
 //----------------------------------------------------------------------------------------------------
 
@@ -39,15 +39,15 @@ class CTPPSLHCInfoRandomXangleESSource : public edm::ESProducer, public edm::Eve
   unsigned int m_generateEveryNEvents;
 
   double m_beamEnergy;
-  double m_betaStar;
 
   std::unique_ptr<CLHEP::HepRandomEngine> m_engine;
 
   struct BinData {
-    double min, max, xangle;
+    double min, max;
+    double xangle, betaStar;
   };
 
-  std::vector<BinData> binData;
+  std::vector<BinData> xangleBetaStarBins;
 };
 
 //----------------------------------------------------------------------------------------------------
@@ -59,35 +59,42 @@ CTPPSLHCInfoRandomXangleESSource::CTPPSLHCInfoRandomXangleESSource(const edm::Pa
       m_generateEveryNEvents(conf.getParameter<unsigned int>("generateEveryNEvents")),
 
       m_beamEnergy(conf.getParameter<double>("beamEnergy")),
-      m_betaStar(conf.getParameter<double>("betaStar")),
 
       m_engine(new CLHEP::HepJamesRandom(conf.getParameter<unsigned int>("seed"))) {
-  const auto &xangleHistogramFile = conf.getParameter<std::string>("xangleHistogramFile");
-  const auto &xangleHistogramObject = conf.getParameter<std::string>("xangleHistogramObject");
+  // get input beta* vs. xangle histogram
+  const auto &xangleBetaStarHistogramFile = conf.getParameter<std::string>("xangleBetaStarHistogramFile");
+  const auto &xangleBetaStarHistogramObject = conf.getParameter<std::string>("xangleBetaStarHistogramObject");
 
-  TFile *f_in = TFile::Open(xangleHistogramFile.c_str());
+  edm::FileInPath fip(xangleBetaStarHistogramFile);
+  std::unique_ptr<TFile> f_in(TFile::Open(fip.fullPath().c_str()));
   if (!f_in)
-    throw cms::Exception("PPS") << "Cannot open input file '" << xangleHistogramFile << "'.";
+    throw cms::Exception("PPS") << "Cannot open input file '" << xangleBetaStarHistogramFile << "'.";
 
-  TH1D *h_xangle = (TH1D *)f_in->Get(xangleHistogramObject.c_str());
-  if (!h_xangle)
-    throw cms::Exception("PPS") << "Cannot load input object '" << xangleHistogramObject << "'.";
+  TH2D *h_xangle_beta_star = (TH2D *)f_in->Get(xangleBetaStarHistogramObject.c_str());
+  if (!h_xangle_beta_star)
+    throw cms::Exception("PPS") << "Cannot load input object '" << xangleBetaStarHistogramObject << "'.";
 
-  double s = 0.;
-  for (int bi = 1; bi <= h_xangle->GetNbinsX(); ++bi)
-    s += h_xangle->GetBinContent(bi);
-
-  double cw = 0.;
-  for (int bi = 1; bi <= h_xangle->GetNbinsX(); ++bi) {
-    double xangle = h_xangle->GetBinCenter(bi);
-    double w = h_xangle->GetBinContent(bi) / s;
-
-    binData.push_back({cw, cw + w, xangle});
-
-    cw += w;
+  // parse histogram
+  double sum = 0.;
+  for (int x = 1; x <= h_xangle_beta_star->GetNbinsX(); ++x) {
+    for (int y = 1; y <= h_xangle_beta_star->GetNbinsY(); ++y)
+      sum += h_xangle_beta_star->GetBinContent(x, y);
   }
 
-  delete f_in;
+  double cw = 0;
+  for (int x = 1; x <= h_xangle_beta_star->GetNbinsX(); ++x) {
+    for (int y = 1; y <= h_xangle_beta_star->GetNbinsY(); ++y) {
+      const double c = h_xangle_beta_star->GetBinContent(x, y);
+      const double xangle = h_xangle_beta_star->GetXaxis()->GetBinCenter(x);
+      const double betaStar = h_xangle_beta_star->GetYaxis()->GetBinCenter(y);
+
+      if (c > 0.) {
+        const double rc = c / sum;
+        xangleBetaStarBins.push_back({cw, cw + rc, xangle, betaStar});
+        cw += rc;
+      }
+    }
+  }
 
   setWhatProduced(this, m_label);
   findingRecord<LHCInfoRcd>();
@@ -104,11 +111,10 @@ void CTPPSLHCInfoRandomXangleESSource::fillDescriptions(edm::ConfigurationDescri
 
   desc.add<unsigned int>("generateEveryNEvents", 1)->setComment("how often to generate new xangle");
 
-  desc.add<std::string>("xangleHistogramFile", "")->setComment("ROOT file with xangle distribution");
-  desc.add<std::string>("xangleHistogramObject", "")->setComment("xangle distribution object in the ROOT file");
+  desc.add<std::string>("xangleBetaStarHistogramFile", "")->setComment("ROOT file with xangle distribution");
+  desc.add<std::string>("xangleBetaStarHistogramObject", "")->setComment("xangle distribution object in the ROOT file");
 
   desc.add<double>("beamEnergy", 0.)->setComment("beam energy");
-  desc.add<double>("betaStar", 0.)->setComment("beta*");
 
   descriptions.add("ctppsLHCInfoRandomXangleESSource", desc);
 }
@@ -128,23 +134,23 @@ void CTPPSLHCInfoRandomXangleESSource::setIntervalFor(const edm::eventsetup::Eve
 edm::ESProducts<std::unique_ptr<LHCInfo>> CTPPSLHCInfoRandomXangleESSource::produce(const LHCInfoRcd &) {
   auto output = std::make_unique<LHCInfo>();
 
+  double xangle = 0., betaStar = 0.;
   const double u = CLHEP::RandFlat::shoot(m_engine.get(), 0., 1.);
-
-  double xangle = 0.;
-  for (const auto &d : binData) {
+  for (const auto &d : xangleBetaStarBins) {
     if (d.min <= u && u <= d.max) {
       xangle = d.xangle;
+      betaStar = d.betaStar;
       break;
     }
   }
 
   output->setEnergy(m_beamEnergy);
   output->setCrossingAngle(xangle);
-  output->setBetaStar(m_betaStar);
+  output->setBetaStar(betaStar);
 
   return edm::es::products(std::move(output));
 }
 
 //----------------------------------------------------------------------------------------------------
 
-DEFINE_FWK_EVENTSETUP_SOURCE(CTPPSLHCInfoRandomXangleESSource);
+DEFINE_FWK_EVENTSETUP_SOURCE(CTPPSLHCInfoRandomXangleESSource);

Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc

@@ -16,13 +16,16 @@
 
 #include "TFile.h"
 #include "TH1D.h"
+#include "TH2D.h"
 
 //----------------------------------------------------------------------------------------------------
 
 class CTPPSLHCInfoPlotter : public edm::one::EDAnalyzer<> {
 public:
   explicit CTPPSLHCInfoPlotter(const edm::ParameterSet &);
 
+  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
+
 private:
   void analyze(const edm::Event &, const edm::EventSetup &) override;
   void endJob() override;
@@ -34,6 +37,7 @@ class CTPPSLHCInfoPlotter : public edm::one::EDAnalyzer<> {
   TH1D *h_beamEnergy_;
   TH1D *h_xangle_;
   TH1D *h_betaStar_;
+  TH2D *h2_betaStar_vs_xangle_;
 
   TH1D *h_fill_;
   TH1D *h_run_;
@@ -53,19 +57,39 @@ CTPPSLHCInfoPlotter::CTPPSLHCInfoPlotter(const edm::ParameterSet &iConfig)
       h_beamEnergy_(new TH1D("h_beamEnergy", ";beam energy   (GeV)", 81, -50., 8050.)),
       h_xangle_(new TH1D("h_xangle", ";(half) crossing angle   (#murad)", 201, -0.5, 200.5)),
       h_betaStar_(new TH1D("h_betaStar", ";#beta^{*}   (m)", 101, -0.005, 1.005)),
+      h2_betaStar_vs_xangle_(new TH2D("h2_betaStar_vs_xangle",
+                                      ";(half) crossing angle   (#murad);#beta^{*}   (m)",
+                                      201,
+                                      -0.5,
+                                      200.5,
+                                      101,
+                                      -0.005,
+                                      1.005)),
 
       h_fill_(new TH1D("h_fill", ";fill", 4001, 3999.5, 8000.5)),
       h_run_(new TH1D("h_run", ";run", 6000, 270E3, 330E3)) {}
 
 //----------------------------------------------------------------------------------------------------
 
+void CTPPSLHCInfoPlotter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
+  edm::ParameterSetDescription desc;
+
+  desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
+  desc.add<std::string>("outputFile", "")->setComment("output file");
+
+  descriptions.add("ctppsLHCInfoPlotter", desc);
+}
+
+//----------------------------------------------------------------------------------------------------
+
 void CTPPSLHCInfoPlotter::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
   edm::ESHandle<LHCInfo> hLHCInfo;
   iSetup.get<LHCInfoRcd>().get(lhcInfoLabel_, hLHCInfo);
 
   h_beamEnergy_->Fill(hLHCInfo->energy());
   h_xangle_->Fill(hLHCInfo->crossingAngle());
   h_betaStar_->Fill(hLHCInfo->betaStar());
+  h2_betaStar_vs_xangle_->Fill(hLHCInfo->crossingAngle(), hLHCInfo->betaStar());
 
   h_fill_->Fill(hLHCInfo->fillNumber());
   h_run_->Fill(iEvent.id().run());
@@ -79,6 +103,7 @@ void CTPPSLHCInfoPlotter::endJob() {
   h_beamEnergy_->Write();
   h_xangle_->Write();
   h_betaStar_->Write();
+  h2_betaStar_vs_xangle_->Write();
 
   h_fill_->Write();
   h_run_->Write();

Validation/CTPPS/python/simu_config/base_cff.py

@@ -173,13 +173,14 @@ def UseCrossingAngle(xangle, process):
   process.ctppsBeamParametersESSource.halfXangleX45 = xangle * 1E-6
   process.ctppsBeamParametersESSource.halfXangleX56 = xangle * 1E-6
 
-def UseCrossingAngleHistgoram(process, f, obj):
+default_xangle_beta_star_file = "CalibPPS/ESProducers/data/xangle_beta_distributions/version1.root"
+
+def UseXangleBetaStarHistogram(process, f, obj):
   process.load("CalibPPS.ESProducers.ctppsLHCInfoRandomXangleESSource_cfi")
-  process.ctppsLHCInfoRandomXangleESSource.generateEveryNEvents = 1
-  process.ctppsLHCInfoRandomXangleESSource.xangleHistogramFile = f
-  process.ctppsLHCInfoRandomXangleESSource.xangleHistogramObject = obj
+  process.ctppsLHCInfoRandomXangleESSource.generateEveryNEvents = 10 # this is to be synchronised with source.numberEventsInLuminosityBlock
+  process.ctppsLHCInfoRandomXangleESSource.xangleBetaStarHistogramFile = f
+  process.ctppsLHCInfoRandomXangleESSource.xangleBetaStarHistogramObject = obj
   process.ctppsLHCInfoRandomXangleESSource.beamEnergy = ctppsLHCInfoESSource.beamEnergy
-  process.ctppsLHCInfoRandomXangleESSource.betaStar = ctppsLHCInfoESSource.betaStar
 
   del process.ctppsLHCInfoESSource
 

Validation/CTPPS/python/simu_config/year_2016_postTS2_cff.py

@@ -38,11 +38,13 @@
 ctppsDirectProtonSimulation.empiricalAperture45="6.10374E-05+(([xi]<0.113491)*0.00795942+([xi]>=0.113491)*0.01935)*([xi]-0.113491)"
 ctppsDirectProtonSimulation.empiricalAperture56="([xi]-0.110)/130.0"
 
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2016_postTS2/h2_betaStar_vs_xangle")
 
 # defaults
 def SetDefaults(process):
-  UseCrossingAngle(140, process)
-
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2016_postTS2")
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2016_preTS2_cff.py

@@ -38,11 +38,13 @@
 ctppsDirectProtonSimulation.empiricalAperture45="3.76296E-05+(([xi]<0.117122)*0.00712775+([xi]>=0.117122)*0.0148651)*([xi]-0.117122)"
 ctppsDirectProtonSimulation.empiricalAperture56="1.85954E-05+(([xi]<0.14324)*0.00475349+([xi]>=0.14324)*0.00629514)*([xi]-0.14324)"
 
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(185, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2016_preTS2/h2_betaStar_vs_xangle")
 
 # defaults
 def SetDefaults(process):
-  UseCrossingAngle(185, process)
-
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2016_preTS2")
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2017_cff.py

@@ -54,7 +54,3 @@
   rp_56_N = cms.uint32(103),
   rp_56_F = cms.uint32(123)
 )
-
-# defaults
-def SetDefaults(process):
-  UseCrossingAngle(140, process)

Validation/CTPPS/python/simu_config/year_2017_postTS2_cff.py

@@ -18,6 +18,13 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "2 * 0.130"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "2 * 0.130"
 
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2017_postTS2")
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2017_postTS2/h2_betaStar_vs_xangle")
+
+# defaults
+def SetDefaults(process):
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2017_preTS2_cff.py

@@ -18,6 +18,13 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "2 * (0.0025*(x-3) + 0.080)"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "2 * (0.0050*(x-3) + 0.060)"
 
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2017_preTS2")
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2017_preTS2/h2_betaStar_vs_xangle")
+
+# defaults
+def SetDefaults(process):
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2018_TS1_TS2_cff.py

@@ -12,6 +12,13 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "2 * ( (x<10)*(-0.0086*(x-10) + 0.100) + (x>=10)*(0.100) )"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "2 * ( (x<8) *(-0.0100*(x-8)  + 0.100) + (x>=8) *(-0.0027*(x-8) + 0.100) )"
 
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2018_TS1_TS2")
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2018_TS1_TS2/h2_betaStar_vs_xangle")
+
+# defaults
+def SetDefaults(process):
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2018_cff.py

@@ -66,7 +66,3 @@
   rp_56_N = cms.uint32(103),
   rp_56_F = cms.uint32(123)
 )
-
-# defaults
-def SetDefaults(process):
-  UseCrossingAngle(140, process)

Validation/CTPPS/python/simu_config/year_2018_postTS2_cff.py

@@ -12,6 +12,13 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "2 * (-0.0031 * (x - 3) + 0.16)"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "2 * ( (x<10)*(-0.0057*(x-10) + 0.110) + (x>=10)*(-0.0022*(x-10) + 0.110) )"
 
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2018_postTS2")
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2018_postTS2/h2_betaStar_vs_xangle")
+
+# defaults
+def SetDefaults(process):
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/python/simu_config/year_2018_preTS1_cff.py

@@ -13,6 +13,13 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "999"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "999"
 
-# xangle distribution
-def UseCrossingAngleDistribution(process, f):
-  UseCrossingAngleHistgoram(process, f, "h_xangle_2018_preTS1")
+# xangle/beta* options
+def UseDefaultXangleBetaStar(process):
+  UseCrossingAngle(140, process)
+
+def UseDefaultXangleBetaStarDistribution(process):
+  UseXangleBetaStarHistogram(process, default_xangle_beta_star_file, "2018_preTS1/h2_betaStar_vs_xangle")
+
+# defaults
+def SetDefaults(process):
+  UseDefaultXangleBetaStarDistribution(process)

Validation/CTPPS/test/simu/template_cfg.py

@@ -6,7 +6,7 @@
 # load config
 import Validation.CTPPS.simu_config.year_$CONFIG_cff as config
 process.load("Validation.CTPPS.simu_config.year_$CONFIG_cff")
-config.SetDefaults(process)
+config.UseXangleBetaStarDistribution(process,"../../../../CalibPPS/ESProducers/data/xangle_beta_distributions/version1.root")
 
 # minimal logger settings
 process.MessageLogger = cms.Service("MessageLogger",