Merge pull request #31453 from CTPPS/pps_systematics_rebased

PPS: direct simulation with xangle distributions

CalibPPS/ESProducers/plugins/CTPPSLHCInfoESSource.cc

@@ -20,6 +20,7 @@ class CTPPSLHCInfoESSource : public edm::ESProducer, public edm::EventSetupRecor
 public:
   CTPPSLHCInfoESSource(const edm::ParameterSet &);
   edm::ESProducts<std::unique_ptr<LHCInfo>> produce(const LHCInfoRcd &);
+  static void fillDescriptions(edm::ConfigurationDescriptions &);
 
 private:
   void setIntervalFor(const edm::eventsetup::EventSetupRecordKey &,
@@ -30,6 +31,7 @@ class CTPPSLHCInfoESSource : public edm::ESProducer, public edm::EventSetupRecor
 
   edm::EventRange m_validityRange;
   double m_beamEnergy;
+  double m_betaStar;
   double m_xangle;
 
   bool m_insideValidityRange;
@@ -42,6 +44,7 @@ CTPPSLHCInfoESSource::CTPPSLHCInfoESSource(const edm::ParameterSet &conf)
     : m_label(conf.getParameter<std::string>("label")),
       m_validityRange(conf.getParameter<edm::EventRange>("validityRange")),
       m_beamEnergy(conf.getParameter<double>("beamEnergy")),
+      m_betaStar(conf.getParameter<double>("betaStar")),
       m_xangle(conf.getParameter<double>("xangle")),
       m_insideValidityRange(false) {
   setWhatProduced(this, m_label);
@@ -50,6 +53,22 @@ CTPPSLHCInfoESSource::CTPPSLHCInfoESSource(const edm::ParameterSet &conf)
 
 //----------------------------------------------------------------------------------------------------
 
+void CTPPSLHCInfoESSource::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
+  edm::ParameterSetDescription desc;
+
+  desc.add<std::string>("label", "")->setComment("label of the LHCInfo record");
+
+  desc.add<edm::EventRange>("validityRange", edm::EventRange())->setComment("interval of validity");
+
+  desc.add<double>("beamEnergy", 0.)->setComment("beam energy");
+  desc.add<double>("betaStar", 0.)->setComment("beta*");
+  desc.add<double>("xangle", 0.)->setComment("crossing angle");
+
+  descriptions.add("ctppsLHCInfoESSource", desc);
+}
+
+//----------------------------------------------------------------------------------------------------
+
 void CTPPSLHCInfoESSource::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &key,
                                           const edm::IOVSyncValue &iosv,
                                           edm::ValidityInterval &oValidity) {
@@ -85,9 +104,11 @@ edm::ESProducts<std::unique_ptr<LHCInfo>> CTPPSLHCInfoESSource::produce(const LH
 
   if (m_insideValidityRange) {
     output->setEnergy(m_beamEnergy);
+    output->setBetaStar(m_betaStar);
     output->setCrossingAngle(m_xangle);
   } else {
     output->setEnergy(0.);
+    output->setBetaStar(0.);
     output->setCrossingAngle(0.);
   }
 

CalibPPS/ESProducers/plugins/CTPPSLHCInfoRandomXangleESSource.cc

@@ -66,7 +66,12 @@ CTPPSLHCInfoRandomXangleESSource::CTPPSLHCInfoRandomXangleESSource(const edm::Pa
   const auto &xangleHistogramObject = conf.getParameter<std::string>("xangleHistogramObject");
 
   TFile *f_in = TFile::Open(xangleHistogramFile.c_str());
+  if (!f_in)
+    throw cms::Exception("PPS") << "Cannot open input file '" << xangleHistogramFile << "'.";
+
   TH1D *h_xangle = (TH1D *)f_in->Get(xangleHistogramObject.c_str());
+  if (!h_xangle)
+    throw cms::Exception("PPS") << "Cannot load input object '" << xangleHistogramObject << "'.";
 
   double s = 0.;
   for (int bi = 1; bi <= h_xangle->GetNbinsX(); ++bi)

Validation/CTPPS/plugins/CTPPSProtonReconstructionSimulationValidator.cc

@@ -78,20 +78,20 @@ class CTPPSProtonReconstructionSimulationValidator : public edm::one::EDAnalyzer
 
     PlotGroup()
         : h_de_xi(new TH1D("", ";#xi_{reco} - #xi_{simu}", 100, 0., 0.)),
-          p_de_xi_vs_xi_simu(new TProfile("", ";#xi_{simu};#xi_{reco} - #xi_{simu}", 19, 0.015, 0.205)),
+          p_de_xi_vs_xi_simu(new TProfile("", ";#xi_{simu};#xi_{reco} - #xi_{simu}", 50, 0., 0.25)),
           h_xi_reco_vs_xi_simu(new TH2D("", ";#xi_{simu};#xi_{reco}", 100, 0., 0.30, 100, 0., 0.30)),
 
           h_de_th_x(new TH1D("", ";#theta_{x,reco} - #theta_{x,simu}", 100, 0., 0.)),
-          p_de_th_x_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{x,reco} - #theta_{x,simu}", 19, 0.015, 0.205)),
+          p_de_th_x_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{x,reco} - #theta_{x,simu}", 50, 0., 0.25)),
 
           h_de_th_y(new TH1D("", ";#theta_{y,reco} - #theta_{y,simu}", 100, 0., 0.)),
-          p_de_th_y_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{y,reco} - #theta_{y,simu}", 19, 0.015, 0.205)),
+          p_de_th_y_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{y,reco} - #theta_{y,simu}", 50, 0., 0.25)),
 
           h_de_vtx_y(new TH1D("", ";vtx_{y,reco} - vtx_{y,simu}   (mm)", 100, 0., 0.)),
-          p_de_vtx_y_vs_xi_simu(new TProfile("", ";#xi_{simu};vtx_{y,reco} - vtx_{y,simu} (mm)", 19, 0.015, 0.205)),
+          p_de_vtx_y_vs_xi_simu(new TProfile("", ";#xi_{simu};vtx_{y,reco} - vtx_{y,simu} (mm)", 50, 0., 0.25)),
 
           h_de_t(new TH1D("", ";t_{reco} - t_{simu}", 100, -1., +1.)),
-          p_de_t_vs_xi_simu(new TProfile("", ";xi_{simu};t_{reco} - t_{simu}", 19, 0.015, 0.205)),
+          p_de_t_vs_xi_simu(new TProfile("", ";xi_{simu};t_{reco} - t_{simu}", 50, 0., 0.25)),
           p_de_t_vs_t_simu(new TProfile("", ";t_{simu};t_{reco} - t_{simu}", 20, 0., 5.)) {}
 
     static TGraphErrors profileToRMSGraph(TProfile *p, const char *name = "") {

Validation/CTPPS/python/simu_config/base_cff.py

@@ -17,7 +17,8 @@
   label = cms.string(""),
   validityRange = cms.EventRange("0:min - 999999:max"),
   beamEnergy = cms.double(6500),  # GeV
-  xangle = cms.double(-1)  # murad
+  xangle = cms.double(-1),  # murad
+  betaStar = cms.double(-1)
 )
 
 # beam parameters as determined by PPS
@@ -76,7 +77,8 @@
 
 # default source
 source = cms.Source("EmptySource",
-  firstRun = cms.untracked.uint32(1)
+  firstRun = cms.untracked.uint32(1),
+  numberEventsInLuminosityBlock = cms.untracked.uint32(10)
 )
 
 # particle generator
@@ -117,23 +119,34 @@
 
 #----------------------------------------------------------------------------------------------------
 
-def SetLevel1(process):
-  process.ctppsBeamParametersESSource.vtxStddevX = 0E-4
-  process.ctppsBeamParametersESSource.vtxStddevZ = 0
+def SetSmearingLevel1(obj):
+  obj.vtxStddevX = 0E-4
+  obj.vtxStddevZ = 0
+
+  obj.beamDivX45 = 0E-6
+  obj.beamDivX56 = 0E-6
+  obj.beamDivY45 = 0E-6
+  obj.beamDivY56 = 0E-6
 
-  process.ctppsBeamParametersESSource.beamDivX45 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivX56 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivY45 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivY56 = 0E-6
+def SetLevel1(process):
+  if hasattr(process, "ctppsBeamParametersESSource"):
+    SetSmearingLevel1(process.ctppsBeamParametersESSource)
+  else:
+    SetSmearingLevel1(process.ctppsBeamParametersFromLHCInfoESSource)
 
   process.ctppsDirectProtonSimulation.roundToPitch = False
 
+def SetSmearingLevel2(obj):
+  obj.beamDivX45 = 0E-6
+  obj.beamDivX56 = 0E-6
+  obj.beamDivY45 = 0E-6
+  obj.beamDivY56 = 0E-6
 
 def SetLevel2(process):
-  process.ctppsBeamParametersESSource.beamDivX45 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivX56 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivY45 = 0E-6
-  process.ctppsBeamParametersESSource.beamDivY56 = 0E-6
+  if hasattr(process, "ctppsBeamParametersESSource"):
+    SetSmearingLevel2(process.ctppsBeamParametersESSource)
+  else:
+    SetSmearingLevel2(process.ctppsBeamParametersFromLHCInfoESSource)
 
   process.ctppsDirectProtonSimulation.roundToPitch = False
 
@@ -159,3 +172,30 @@ def UseCrossingAngle(xangle, process):
   process.ctppsLHCInfoESSource.xangle = xangle
   process.ctppsBeamParametersESSource.halfXangleX45 = xangle * 1E-6
   process.ctppsBeamParametersESSource.halfXangleX56 = xangle * 1E-6
+
+def UseCrossingAngleHistgoram(process, f, obj):
+  process.load("CalibPPS.ESProducers.ctppsLHCInfoRandomXangleESSource_cfi")
+  process.ctppsLHCInfoRandomXangleESSource.generateEveryNEvents = 1
+  process.ctppsLHCInfoRandomXangleESSource.xangleHistogramFile = f
+  process.ctppsLHCInfoRandomXangleESSource.xangleHistogramObject = obj
+  process.ctppsLHCInfoRandomXangleESSource.beamEnergy = ctppsLHCInfoESSource.beamEnergy
+  process.ctppsLHCInfoRandomXangleESSource.betaStar = ctppsLHCInfoESSource.betaStar
+
+  del process.ctppsLHCInfoESSource
+
+  process.load("CalibPPS.ESProducers.ctppsBeamParametersFromLHCInfoESSource_cfi")
+  process.ctppsBeamParametersFromLHCInfoESSource.beamDivX45 = process.ctppsBeamParametersESSource.beamDivX45
+  process.ctppsBeamParametersFromLHCInfoESSource.beamDivX56 = process.ctppsBeamParametersESSource.beamDivX56
+  process.ctppsBeamParametersFromLHCInfoESSource.beamDivY45 = process.ctppsBeamParametersESSource.beamDivY45
+  process.ctppsBeamParametersFromLHCInfoESSource.beamDivY56 = process.ctppsBeamParametersESSource.beamDivY56
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetX45 = process.ctppsBeamParametersESSource.vtxOffsetX45
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetX56 = process.ctppsBeamParametersESSource.vtxOffsetX56
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetY45 = process.ctppsBeamParametersESSource.vtxOffsetY45
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetY56 = process.ctppsBeamParametersESSource.vtxOffsetY56
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetZ45 = process.ctppsBeamParametersESSource.vtxOffsetZ45
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxOffsetZ56 = process.ctppsBeamParametersESSource.vtxOffsetZ56
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxStddevX = process.ctppsBeamParametersESSource.vtxStddevX
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxStddevY = process.ctppsBeamParametersESSource.vtxStddevY
+  process.ctppsBeamParametersFromLHCInfoESSource.vtxStddevZ = process.ctppsBeamParametersESSource.vtxStddevZ
+
+  del process.ctppsBeamParametersESSource

Validation/CTPPS/python/simu_config/year_2016_postTS2_cff.py

@@ -42,3 +42,7 @@
 # defaults
 def SetDefaults(process):
   UseCrossingAngle(140, process)
+
+# xangle distribution
+def UseCrossingAngleDistribution(process, f):
+  UseCrossingAngleHistgoram(process, f, "h_xangle_2016_postTS2")

Validation/CTPPS/python/simu_config/year_2016_preTS2_cff.py

@@ -42,3 +42,7 @@
 # defaults
 def SetDefaults(process):
   UseCrossingAngle(185, process)
+
+# xangle distribution
+def UseCrossingAngleDistribution(process, f):
+  UseCrossingAngleHistgoram(process, f, "h_xangle_2016_preTS2")

Validation/CTPPS/python/simu_config/year_2017_postTS2_cff.py

@@ -18,3 +18,6 @@
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "2 * 0.130"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "2 * 0.130"
 
+# xangle distribution
+def UseCrossingAngleDistribution(process, f):
+  UseCrossingAngleHistgoram(process, f, "h_xangle_2017_postTS2")

Validation/CTPPS/python/simu_config/year_2017_preTS2_cff.py

@@ -18,3 +18,6 @@
 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")

Validation/CTPPS/python/simu_config/year_2018_TS1_TS2_cff.py

@@ -11,3 +11,7 @@
 # timing resolution
 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")

Validation/CTPPS/python/simu_config/year_2018_postTS2_cff.py

@@ -11,3 +11,7 @@
 # timing resolution
 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")

Validation/CTPPS/python/simu_config/year_2018_preTS1_cff.py

@@ -12,3 +12,7 @@
 ctppsLocalTrackLiteProducer.includeDiamonds = False
 ctppsDirectProtonSimulation.timeResolutionDiamonds45 = "999"
 ctppsDirectProtonSimulation.timeResolutionDiamonds56 = "999"
+
+# xangle distribution
+def UseCrossingAngleDistribution(process, f):
+  UseCrossingAngleHistgoram(process, f, "h_xangle_2018_preTS1")