Skip to content

Commit

Permalink
Merge pull request #26335 from CTPPS/proton_reco_step2_algo
Browse files Browse the repository at this point in the history 
PPS: proton reconstruction, algorithm update
  • Loading branch information
@cmsbuild
cmsbuild committed Apr 12, 2019
2 parents c769947 + dfa0859 commit c4df47f
Show file tree
Hide file tree
Showing 12 changed files with 491 additions and 190 deletions.
14 changes: 10 additions & 4 deletions CalibPPS/ESProducers/plugins/CTPPSInterpolatedOpticalFunctionsESSource.cc
View file
Expand Up @@ -159,16 +159,22 @@ std::shared_ptr<LHCInterpolatedOpticalFunctionsSetCollection> CTPPSInterpolatedO
const auto rpId = rp_p.first;
const auto &rp_it2 = ofs2.find(rpId);
if (rp_it2 == ofs2.end())
throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "Mismatch between ofs1 and ofs2.";
throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "RP mismatch between ofs1 and ofs2.";

const auto &of1 = rp_p.second;
const auto &of2 = rp_it2->second;

LHCInterpolatedOpticalFunctionsSet iof;
iof.m_z = of1.getScoringPlaneZ();
const size_t num_xi_vals1 = of1.getXiValues().size();
const size_t num_xi_vals2 = of2.getXiValues().size();

if (num_xi_vals1 != num_xi_vals2)
throw cms::Exception("CTPPSInterpolatedOpticalFunctionsESSource") << "Size mismatch between ofs1 and ofs2.";

const size_t num_xi_vals = of1.getXiValues().size();
const size_t num_xi_vals = num_xi_vals1;

LHCInterpolatedOpticalFunctionsSet iof;
iof.m_z = of1.getScoringPlaneZ();
iof.m_fcn_values.resize(LHCInterpolatedOpticalFunctionsSet::nFunctions);
iof.m_xi_values.resize(num_xi_vals);

for (size_t fi = 0; fi < of1.getFcnValues().size(); ++fi)
Expand Down
4 changes: 2 additions & 2 deletions CalibPPS/ESProducers/plugins/CTPPSLHCInfoESSource.cc
View file
Expand Up @@ -54,12 +54,12 @@ void CTPPSLHCInfoESSource::setIntervalFor(const edm::eventsetup::EventSetupRecor
const edm::IOVSyncValue& iosv, edm::ValidityInterval& oValidity)
{
// TODO: the IOV specified in config is temporarily replaced with a hardcoded one
m_insideValidityRange = true;

edm::IOVSyncValue beg(edm::EventID(270293, 0, 0), edm::Timestamp(1461016800ULL << 32));
edm::IOVSyncValue end(edm::EventID(290872, 0, 0), edm::Timestamp(1483138800ULL << 32));
oValidity = edm::ValidityInterval(beg, end);

m_insideValidityRange = (beg < iosv && iosv < end);

/*
if (edm::contains(m_validityRange, iosv.eventID()))
{
Expand Down
130 changes: 77 additions & 53 deletions CalibPPS/ESProducers/plugins/CTPPSOpticalFunctionsESSource.cc
View file
Expand Up @@ -6,7 +6,6 @@
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESProducer.h"
#include "FWCore/Framework/interface/EventSetupRecordIntervalFinder.h"
#include "FWCore/Framework/interface/ESProducts.h"

#include "CondFormats/CTPPSReadoutObjects/interface/LHCOpticalFunctionsSetCollection.h"
#include "CondFormats/DataRecord/interface/CTPPSOpticsRcd.h"
Expand All @@ -22,51 +21,67 @@ class CTPPSOpticalFunctionsESSource: public edm::ESProducer, public edm::EventSe
CTPPSOpticalFunctionsESSource(const edm::ParameterSet &);
~CTPPSOpticalFunctionsESSource() override = default;

edm::ESProducts<std::unique_ptr<LHCOpticalFunctionsSetCollection>> produce(const CTPPSOpticsRcd &);
std::unique_ptr<LHCOpticalFunctionsSetCollection> produce(const CTPPSOpticsRcd &);
static void fillDescriptions(edm::ConfigurationDescriptions&);

private:
void setIntervalFor(const edm::eventsetup::EventSetupRecordKey&, const edm::IOVSyncValue&, edm::ValidityInterval&) override;

edm::EventRange m_validityRange;
bool m_insideValidityRange;

struct FileInfo
{
double xangle;
std::string fileName;
double m_xangle;
std::string m_fileName;
};
std::vector<FileInfo> m_fileInfo;

struct RPInfo
{
std::string dirName;
double scoringPlaneZ;
std::string m_dirName;
double m_scoringPlaneZ;
};

struct Entry
{
edm::EventRange m_validityRange;
std::vector<FileInfo> m_fileInfo;
std::unordered_map<unsigned int, RPInfo> m_rpInfo;
};
std::unordered_map<unsigned int, RPInfo> m_rpInfo;

std::vector<Entry> m_entries;

bool m_currentEntryValid;
unsigned int m_currentEntry;
};

//----------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------

CTPPSOpticalFunctionsESSource::CTPPSOpticalFunctionsESSource(const edm::ParameterSet& conf) :
m_validityRange(conf.getParameter<edm::EventRange>("validityRange")),
m_insideValidityRange(false)
m_currentEntryValid(false),
m_currentEntry(0)
{
for (const auto &pset : conf.getParameter<std::vector<edm::ParameterSet>>("opticalFunctions"))
for (const auto &entry_pset : conf.getParameter<std::vector<edm::ParameterSet>>("configuration"))
{
const double &xangle = pset.getParameter<double>("xangle");
const std::string &fileName = pset.getParameter<edm::FileInPath>("fileName").fullPath();
m_fileInfo.push_back({xangle, fileName});
}
edm::EventRange validityRange = entry_pset.getParameter<edm::EventRange>("validityRange");

for (const auto &pset : conf.getParameter<std::vector<edm::ParameterSet>>("scoringPlanes"))
{
const unsigned int rpId = pset.getParameter<unsigned int>("rpId");
const std::string dirName = pset.getParameter<std::string>("dirName");
const double z = pset.getParameter<double>("z");
const RPInfo entry = {dirName, z};
m_rpInfo.emplace(rpId, entry);
std::vector<FileInfo> fileInfo;
for (const auto &pset : entry_pset.getParameter<std::vector<edm::ParameterSet>>("opticalFunctions"))
{
const double &xangle = pset.getParameter<double>("xangle");
const std::string &fileName = pset.getParameter<edm::FileInPath>("fileName").fullPath();
fileInfo.push_back({xangle, fileName});
}

std::unordered_map<unsigned int, RPInfo> rpInfo;
for (const auto &pset : entry_pset.getParameter<std::vector<edm::ParameterSet>>("scoringPlanes"))
{
const unsigned int rpId = pset.getParameter<unsigned int>("rpId");
const std::string dirName = pset.getParameter<std::string>("dirName");
const double z = pset.getParameter<double>("z");
const RPInfo entry = {dirName, z};
rpInfo.emplace(rpId, entry);
}

m_entries.push_back({validityRange, fileInfo, rpInfo});
}

setWhatProduced(this);
Expand All @@ -78,56 +93,60 @@ CTPPSOpticalFunctionsESSource::CTPPSOpticalFunctionsESSource(const edm::Paramete
void CTPPSOpticalFunctionsESSource::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &key,
const edm::IOVSyncValue& iosv, edm::ValidityInterval& oValidity)
{
if (edm::contains(m_validityRange, iosv.eventID()))
for (unsigned int idx = 0; idx < m_entries.size(); ++idx)
{
m_insideValidityRange = true;
oValidity = edm::ValidityInterval(edm::IOVSyncValue(m_validityRange.startEventID()), edm::IOVSyncValue(m_validityRange.endEventID()));
} else {
m_insideValidityRange = false;
const auto &entry = m_entries[idx];

if (iosv.eventID() < m_validityRange.startEventID())
// is within an entry ?
if (edm::contains(entry.m_validityRange, iosv.eventID()))
{
edm::RunNumber_t run = m_validityRange.startEventID().run();
edm::LuminosityBlockNumber_t lb = m_validityRange.startEventID().luminosityBlock();
edm::EventID endEvent = (lb > 1) ? edm::EventID(run, lb-1, 0) : edm::EventID(run-1, edm::EventID::maxLuminosityBlockNumber(), 0);

oValidity = edm::ValidityInterval(edm::IOVSyncValue::beginOfTime(), edm::IOVSyncValue(endEvent));
} else {
edm::RunNumber_t run = m_validityRange.startEventID().run();
edm::LuminosityBlockNumber_t lb = m_validityRange.startEventID().luminosityBlock();
edm::EventID beginEvent = (lb < edm::EventID::maxLuminosityBlockNumber()-1) ? edm::EventID(run, lb+1, 0) : edm::EventID(run+1, 0, 0);

oValidity = edm::ValidityInterval(edm::IOVSyncValue(beginEvent), edm::IOVSyncValue::endOfTime());
m_currentEntryValid = true;
m_currentEntry = idx;
oValidity = edm::ValidityInterval(edm::IOVSyncValue(entry.m_validityRange.startEventID()), edm::IOVSyncValue(entry.m_validityRange.endEventID()));
return;
}
}

// not within any entry
m_currentEntryValid = false;
m_currentEntry = 0;

edm::LogInfo("") << "No configuration entry found for event " << iosv.eventID() << ", no optical functions will be available.";

const edm::EventID start(iosv.eventID().run(), iosv.eventID().luminosityBlock(), iosv.eventID().event());
const edm::EventID end(iosv.eventID().run(), iosv.eventID().luminosityBlock(), iosv.eventID().event());
oValidity = edm::ValidityInterval(edm::IOVSyncValue(start), edm::IOVSyncValue(end));
}

//----------------------------------------------------------------------------------------------------

edm::ESProducts<std::unique_ptr<LHCOpticalFunctionsSetCollection> >
std::unique_ptr<LHCOpticalFunctionsSetCollection>
CTPPSOpticalFunctionsESSource::produce(const CTPPSOpticsRcd &)
{
// fill the output
// prepare output, empty by default
auto output = std::make_unique<LHCOpticalFunctionsSetCollection>();

if (m_insideValidityRange)
// fill the output
if (m_currentEntryValid)
{
for (const auto &fi : m_fileInfo)
const auto &entry = m_entries[m_currentEntry];

for (const auto &fi : entry.m_fileInfo)
{
std::unordered_map<unsigned int, LHCOpticalFunctionsSet> xa_data;

for (const auto &rpi : m_rpInfo)
for (const auto &rpi : entry.m_rpInfo)
{
LHCOpticalFunctionsSet fcn(fi.fileName, rpi.second.dirName, rpi.second.scoringPlaneZ);
LHCOpticalFunctionsSet fcn(fi.m_fileName, rpi.second.m_dirName, rpi.second.m_scoringPlaneZ);
xa_data.emplace(rpi.first, std::move(fcn));
}

output->emplace(fi.xangle, xa_data);
output->emplace(fi.m_xangle, xa_data);
}
}

// commit the output
return edm::es::products(std::move(output));
return output;
}

//----------------------------------------------------------------------------------------------------
Expand All @@ -137,20 +156,25 @@ CTPPSOpticalFunctionsESSource::fillDescriptions(edm::ConfigurationDescriptions&
{
edm::ParameterSetDescription desc;

desc.add<edm::EventRange>("validityRange", edm::EventRange())->setComment("interval of validity");
edm::ParameterSetDescription config_desc;

config_desc.add<edm::EventRange>("validityRange", edm::EventRange())->setComment("interval of validity");

edm::ParameterSetDescription of_desc;
of_desc.add<double>("xangle")->setComment("half crossing angle value in urad");
of_desc.add<edm::FileInPath>("fileName")->setComment("ROOT file with optical functions");
std::vector<edm::ParameterSet> of;
desc.addVPSet("opticalFunctions", of_desc, of)->setComment("list of optical functions at different crossing angles");
config_desc.addVPSet("opticalFunctions", of_desc, of)->setComment("list of optical functions at different crossing angles");

edm::ParameterSetDescription sp_desc;
sp_desc.add<unsigned int>("rpId")->setComment("associated detector DetId");
sp_desc.add<std::string>("dirName")->setComment("associated path to the optical functions file");
sp_desc.add<double>("z")->setComment("longitudinal position at scoring plane/detector");
std::vector<edm::ParameterSet> sp;
desc.addVPSet("scoringPlanes", sp_desc, sp)->setComment("list of sensitive planes/detectors stations");
config_desc.addVPSet("scoringPlanes", sp_desc, sp)->setComment("list of sensitive planes/detectors stations");

std::vector<edm::ParameterSet> config;
desc.addVPSet("configuration", config_desc, sp)->setComment("list of configuration blocks");

descriptions.add("ctppsOpticalFunctionsESSource", desc);
}
Expand Down
18 changes: 10 additions & 8 deletions CalibPPS/ESProducers/python/ctppsOpticalFunctions_cff.py
View file
Expand Up @@ -2,25 +2,27 @@

from CalibPPS.ESProducers.ctppsLHCInfo_cff import *

from CalibPPS.ESProducers.ctppsOpticalFunctionsESSource_cfi import ctppsOpticalFunctionsESSource as _optics_tmp
from CalibPPS.ESProducers.ctppsOpticalFunctionsESSource_cfi import *

# optical functions for 2016 data
ctppsOpticalFunctionsESSource_2016 = _optics_tmp.clone(
validityRange = cms.EventRange("270293:min - 290872:max"),
# add 2016 pre-TS2 configuration
config_2016_preTS2 = cms.PSet(
validityRange = cms.EventRange("273725:min - 280385:max"),

opticalFunctions = cms.VPSet(
cms.PSet( xangle = cms.double(185), fileName = cms.FileInPath("CalibPPS/ESProducers/data/optical_functions_2016.root") )
),

scoringPlanes = cms.VPSet(
# z in cm
cms.PSet( rpId = cms.uint32(0x76100000), dirName = cms.string("XRPH_C6L5_B2"), z = cms.double(-20382.6) ), # RP 002
cms.PSet( rpId = cms.uint32(0x76180000), dirName = cms.string("XRPH_D6L5_B2"), z = cms.double(-21255.1) ), # RP 003
cms.PSet( rpId = cms.uint32(0x77100000), dirName = cms.string("XRPH_C6R5_B1"), z = cms.double(+20382.6) ), # RP 102
cms.PSet( rpId = cms.uint32(0x77180000), dirName = cms.string("XRPH_D6R5_B1"), z = cms.double(+21255.1) ), # RP 103
cms.PSet( rpId = cms.uint32(0x76100000), dirName = cms.string("XRPH_C6L5_B2"), z = cms.double(-20382.6) ), # RP 002, strip
cms.PSet( rpId = cms.uint32(0x76180000), dirName = cms.string("XRPH_D6L5_B2"), z = cms.double(-21255.1) ), # RP 003, strip
cms.PSet( rpId = cms.uint32(0x77100000), dirName = cms.string("XRPH_C6R5_B1"), z = cms.double(+20382.6) ), # RP 102, strip
cms.PSet( rpId = cms.uint32(0x77180000), dirName = cms.string("XRPH_D6R5_B1"), z = cms.double(+21255.1) ), # RP 103, strip
)
)

ctppsOpticalFunctionsESSource.configuration.append(config_2016_preTS2)

# optics interpolation between crossing angles
from CalibPPS.ESProducers.ctppsInterpolatedOpticalFunctionsESSource_cfi import *
ctppsInterpolatedOpticalFunctionsESSource.lhcInfoLabel = ctppsLHCInfoLabel
8 changes: 6 additions & 2 deletions DataFormats/ProtonReco/interface/ForwardProton.h
View file
Expand Up @@ -38,8 +38,9 @@ namespace reco
ForwardProton();
/// constructor from refit parameters, fitted vertex and momentum, and longitudinal fractional momentum loss
ForwardProton( double chi2, double ndof, const Point& vtx, const Vector& momentum, float xi,
const CovarianceMatrix& cov, ReconstructionMethod method,
const CTPPSLocalTrackLiteRefVector& local_tracks, bool valid );
const CovarianceMatrix& cov, ReconstructionMethod method,
const CTPPSLocalTrackLiteRefVector& local_tracks, bool valid,
const float time=0., const float time_err=0. );

/// fitted vertex position
const Point& vertex() const { return vertex_; }
Expand Down Expand Up @@ -107,8 +108,11 @@ namespace reco

/// time of proton arrival at forward stations
float time() const { return time_; }
void setTime(float time) { time_ = time; }

/// uncertainty on time of proton arrival at forward stations
float timeError() const { return time_err_; }
void setTimeError(float time_err) { time_err_ = time_err; }

/// set the flag for the fit validity
void setValidFit( bool valid = true ) { valid_fit_ = valid; }
Expand Down
5 changes: 3 additions & 2 deletions DataFormats/ProtonReco/src/ForwardProton.cc
View file
Expand Up @@ -17,9 +17,10 @@ ForwardProton::ForwardProton() :

ForwardProton::ForwardProton( double chi2, double ndof, const Point& vtx, const Vector& momentum, float xi,
const CovarianceMatrix& cov, ReconstructionMethod method,
const CTPPSLocalTrackLiteRefVector& local_tracks, bool valid ) :
const CTPPSLocalTrackLiteRefVector& local_tracks, bool valid,
const float time, const float time_err ) :
vertex_( vtx ), momentum_( momentum ),
time_( 0. ), time_err_( 0. ), xi_( xi ),
time_( time ), time_err_( time_err ), xi_( xi ),
covariance_( cov ), chi2_( chi2 ), ndof_( ndof ),
valid_fit_( valid ), method_( method ), contributing_local_tracks_( local_tracks )
{}
Expand Down
6 changes: 4 additions & 2 deletions PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
View file
Expand Up @@ -514,9 +514,11 @@ def miniAOD_customizeData(process):
runOnData( process, outputModules = [] )
process.load("RecoCTPPS.TotemRPLocal.ctppsLocalTrackLiteProducer_cff")
process.load("RecoCTPPS.ProtonReconstruction.ctppsProtons_cff")
process.load("Geometry.VeryForwardGeometry.geometryRPFromDB_cfi")
task = getPatAlgosToolsTask(process)
task.add(process.ctppsLocalTrackLiteProducer)
task.add(process.ctppsProtons)
from Configuration.Eras.Modifier_ctpps_2016_cff import ctpps_2016
ctpps_2016.toModify(task, func=lambda t: t.add(process.ctppsLocalTrackLiteProducer))
ctpps_2016.toModify(task, func=lambda t: t.add(process.ctppsProtons))

def miniAOD_customizeAllData(process):
miniAOD_customizeCommon(process)
Expand Down
4 changes: 4 additions & 0 deletions RecoCTPPS/ProtonReconstruction/BuildFile.xml
View file
Expand Up @@ -8,6 +8,10 @@

<use name="CondFormats/RunInfo"/>
<use name="CondFormats/CTPPSReadoutObjects"/>
<use name="CondFormats/DataRecord"/>

<use name="Geometry/Records"/>
<use name="Geometry/VeryForwardGeometryBuilder"/>

<export>
<lib name="1"/>
Expand Down
6 changes: 2 additions & 4 deletions RecoCTPPS/ProtonReconstruction/interface/ProtonReconstructionAlgorithm.h
View file
Expand Up @@ -33,13 +33,11 @@ class ProtonReconstructionAlgorithm
void release();

/// run proton reconstruction using single-RP strategy
void reconstructFromSingleRP(const CTPPSLocalTrackLiteRefVector& input,
reco::ForwardProtonCollection& output,
reco::ForwardProton reconstructFromSingleRP(const CTPPSLocalTrackLiteRef &track,
const LHCInfo& lhcInfo, std::ostream& os) const;

/// run proton reconstruction using multiple-RP strategy
void reconstructFromMultiRP(const CTPPSLocalTrackLiteRefVector& input,
reco::ForwardProtonCollection& output,
reco::ForwardProton reconstructFromMultiRP(const CTPPSLocalTrackLiteRefVector &tracks,
const LHCInfo& lhcInfo, std::ostream& os) const;

private:
Expand Down

0 comments on commit c4df47f

Please sign in to comment.