Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #33119 from nminafra/ecalTiming_11_3_X_kucc
Cross-correlation algorithm for the ECAL timing reconstruction
- Loading branch information
Showing
6 changed files
with
331 additions
and
3 deletions.
There are no files selected for viewing
43 changes: 43 additions & 0 deletions
43
RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,43 @@ | ||
#ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimingCCAlgo_HH | ||
#define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimingCCAlgo_HH | ||
|
||
/** \class EcalUncalibRecHitTimingCCAlgo | ||
* CrossCorrelation algorithm for timing reconstruction | ||
* | ||
* \author N. Minafra, J. King, C. Rogan | ||
*/ | ||
|
||
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h" | ||
#include "FWCore/MessageLogger/interface/MessageLogger.h" | ||
|
||
#include "CondFormats/EcalObjects/interface/EcalPedestals.h" | ||
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h" | ||
#include "DataFormats/EcalDigi/interface/EcalConstants.h" | ||
#include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h" | ||
|
||
class EcalUncalibRecHitTimingCCAlgo { | ||
public: | ||
EcalUncalibRecHitTimingCCAlgo(const float startTime, const float stopTime, const float targetTimePrecision); | ||
double computeTimeCC(const EcalDataFrame& dataFrame, | ||
const std::vector<double>& amplitudes, | ||
const EcalPedestals::Item* aped, | ||
const EcalMGPAGainRatio* aGain, | ||
const FullSampleVector& fullpulse, | ||
EcalUncalibratedRecHit& uncalibRecHit, | ||
float& errOnTime) const; | ||
|
||
private: | ||
const float startTime_; | ||
const float stopTime_; | ||
const float targetTimePrecision_; | ||
|
||
static constexpr int TIME_WHEN_NOT_CONVERGING = 100; | ||
static constexpr int MAX_NUM_OF_ITERATIONS = 30; | ||
static constexpr int MIN_NUM_OF_ITERATIONS = 2; | ||
static constexpr float GLOBAL_TIME_SHIFT = 100; | ||
|
||
FullSampleVector interpolatePulse(const FullSampleVector& fullpulse, const float t = 0) const; | ||
float computeCC(const std::vector<float>& samples, const FullSampleVector& sigmalTemplate, const float t) const; | ||
}; | ||
|
||
#endif |
158 changes: 158 additions & 0 deletions
158
RecoLocalCalo/EcalRecAlgos/src/EcalUncalibRecHitTimingCCAlgo.cc
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,158 @@ | ||
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h" | ||
|
||
EcalUncalibRecHitTimingCCAlgo::EcalUncalibRecHitTimingCCAlgo(const float startTime, | ||
const float stopTime, | ||
const float targetTimePrecision) | ||
: startTime_(startTime), stopTime_(stopTime), targetTimePrecision_(targetTimePrecision) {} | ||
|
||
double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFrame, | ||
const std::vector<double>& amplitudes, | ||
const EcalPedestals::Item* aped, | ||
const EcalMGPAGainRatio* aGain, | ||
const FullSampleVector& fullpulse, | ||
EcalUncalibratedRecHit& uncalibRecHit, | ||
float& errOnTime) const { | ||
constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES; | ||
|
||
double maxamplitude = -std::numeric_limits<double>::max(); | ||
float pulsenorm = 0.; | ||
|
||
std::vector<float> pedSubSamples(nsample); | ||
for (unsigned int iSample = 0; iSample < nsample; iSample++) { | ||
const EcalMGPASample& sample = dataFrame.sample(iSample); | ||
|
||
float amplitude = 0.; | ||
int gainId = sample.gainId(); | ||
|
||
double pedestal = 0.; | ||
double gainratio = 1.; | ||
|
||
if (gainId == 0 || gainId == 3) { | ||
pedestal = aped->mean_x1; | ||
gainratio = aGain->gain6Over1() * aGain->gain12Over6(); | ||
} else if (gainId == 1) { | ||
pedestal = aped->mean_x12; | ||
gainratio = 1.; | ||
} else if (gainId == 2) { | ||
pedestal = aped->mean_x6; | ||
gainratio = aGain->gain12Over6(); | ||
} | ||
|
||
amplitude = (static_cast<float>(sample.adc()) - pedestal) * gainratio; | ||
|
||
if (gainId == 0) { | ||
//saturation | ||
amplitude = (4095. - pedestal) * gainratio; | ||
} | ||
|
||
pedSubSamples[iSample] = amplitude; | ||
|
||
if (amplitude > maxamplitude) { | ||
maxamplitude = amplitude; | ||
} | ||
pulsenorm += fullpulse(iSample); | ||
} | ||
|
||
int ipulse = -1; | ||
for (auto const& amplit : amplitudes) { | ||
ipulse++; | ||
int bxp3 = ipulse - 2; | ||
int firstsamplet = std::max(0, bxp3); | ||
int offset = 7 - bxp3; | ||
|
||
for (unsigned int isample = firstsamplet; isample < nsample; ++isample) { | ||
auto const pulse = fullpulse(isample + offset); | ||
pedSubSamples[isample] = std::max(0., pedSubSamples[isample] - amplit * pulse / pulsenorm); | ||
} | ||
} | ||
|
||
// Start of time computation | ||
float tStart = startTime_ + GLOBAL_TIME_SHIFT; | ||
float tStop = stopTime_ + GLOBAL_TIME_SHIFT; | ||
float tM = (tStart + tStop) / 2; | ||
|
||
float distStart, distStop; | ||
int counter = 0; | ||
|
||
do { | ||
++counter; | ||
distStart = computeCC(pedSubSamples, fullpulse, tStart); | ||
distStop = computeCC(pedSubSamples, fullpulse, tStop); | ||
|
||
if (distStart > distStop) { | ||
tStop = tM; | ||
} else { | ||
tStart = tM; | ||
} | ||
tM = (tStart + tStop) / 2; | ||
|
||
} while (tStop - tStart > targetTimePrecision_ && counter < MAX_NUM_OF_ITERATIONS); | ||
|
||
tM -= GLOBAL_TIME_SHIFT; | ||
errOnTime = targetTimePrecision_; | ||
|
||
if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) { | ||
tM = TIME_WHEN_NOT_CONVERGING * ecalPh1::Samp_Period; | ||
//Negative error means that there was a problem with the CC | ||
errOnTime = -targetTimePrecision_ / ecalPh1::Samp_Period; | ||
} | ||
|
||
return -tM / ecalPh1::Samp_Period; | ||
} | ||
|
||
FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse, | ||
const float time) const { | ||
// t is in ns | ||
int shift = time / ecalPh1::Samp_Period; | ||
if (time < 0) | ||
shift -= 1; | ||
float tt = time / ecalPh1::Samp_Period - shift; | ||
|
||
FullSampleVector interpPulse; | ||
// 2nd poly with avg | ||
unsigned int numberOfSamples = fullpulse.size(); | ||
auto facM1orP2 = 0.25 * tt * (tt - 1); | ||
auto fac = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1); | ||
auto facP1 = (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt; | ||
for (unsigned int i = 1; i < numberOfSamples - 2; ++i) { | ||
float a = | ||
facM1orP2 * fullpulse[i - 1] + fac * fullpulse[i] + facP1 * fullpulse[i + 1] + facM1orP2 * fullpulse[i + 2]; | ||
if (a > 0) | ||
interpPulse[i] = a; | ||
else | ||
interpPulse[i] = 0; | ||
} | ||
interpPulse[0] = facM1orP2 * fullpulse[0] + facP1 * fullpulse[1] + facM1orP2 * fullpulse[2]; | ||
interpPulse[numberOfSamples - 2] = facM1orP2 * fullpulse[numberOfSamples - 3] + fac * fullpulse[numberOfSamples - 2] + | ||
facP1 * fullpulse[numberOfSamples - 1]; | ||
interpPulse[numberOfSamples - 1] = 2 * facM1orP2 * fullpulse[numberOfSamples - 2] - | ||
4 * facM1orP2 * fullpulse[numberOfSamples - 1] + | ||
facP1 * fullpulse[numberOfSamples - 1]; | ||
|
||
FullSampleVector interpPulseShifted; | ||
for (int i = 0; i < interpPulseShifted.size(); ++i) { | ||
if (i + shift >= 0 && i + shift < interpPulse.size()) | ||
interpPulseShifted[i] = interpPulse[i + shift]; | ||
else | ||
interpPulseShifted[i] = 0; | ||
} | ||
return interpPulseShifted; | ||
} | ||
|
||
float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples, | ||
const FullSampleVector& signalTemplate, | ||
const float time) const { | ||
constexpr int exclude = 1; | ||
float powerSamples = 0.; | ||
float powerTemplate = 0.; | ||
float cc = 0.; | ||
auto interpolated = interpolatePulse(signalTemplate, time); | ||
for (int i = exclude; i < int(samples.size() - exclude); ++i) { | ||
powerSamples += std::pow(samples[i], 2); | ||
powerTemplate += std::pow(interpolated[i], 2); | ||
cc += interpolated[i] * samples[i]; | ||
} | ||
|
||
float denominator = std::sqrt(powerTemplate * powerSamples); | ||
return cc / denominator; | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
97 changes: 97 additions & 0 deletions
97
RecoLocalCalo/EcalRecProducers/test/testEcalUncalibRechitProducerWithCC_cfg.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,97 @@ | ||
import FWCore.ParameterSet.Config as cms | ||
|
||
from Configuration.StandardSequences.Eras import eras | ||
|
||
process = cms.Process('RECO', eras.Run2_2018) | ||
|
||
# import of standard configurations | ||
process.load('Configuration.StandardSequences.Services_cff') | ||
process.load('FWCore.MessageService.MessageLogger_cfi') | ||
process.load('Configuration.StandardSequences.GeometryRecoDB_cff') | ||
process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') | ||
process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') | ||
process.load('Configuration.StandardSequences.Reconstruction_cff') | ||
|
||
|
||
|
||
process.source = cms.Source('PoolSource', | ||
fileNames = cms.untracked.vstring( | ||
'/store/data/Run2018D/EGamma/RAW/v1/000/323/414/00000/042D6023-E0A2-8649-8D86-445F752A8F6B.root', | ||
), | ||
) | ||
|
||
|
||
# Other statements | ||
from Configuration.AlCa.GlobalTag import GlobalTag | ||
process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_data', '') | ||
|
||
|
||
process.maxEvents = cms.untracked.PSet( | ||
input = cms.untracked.int32(1000) | ||
) | ||
|
||
#----------------------------------------- | ||
# CMSSW/Hcal non-DQM Related Module import | ||
#----------------------------------------- | ||
process.load('Configuration.StandardSequences.GeometryRecoDB_cff') | ||
process.load("RecoLocalCalo.Configuration.hcalLocalReco_cff") | ||
process.load("RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff") | ||
process.load("EventFilter.HcalRawToDigi.HcalRawToDigi_cfi") | ||
process.load("EventFilter.EcalRawToDigi.EcalUnpackerData_cfi") | ||
process.load("RecoLuminosity.LumiProducer.bunchSpacingProducer_cfi") | ||
|
||
# load both cpu plugins | ||
process.load("RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi") | ||
|
||
## | ||
## force HLT configuration for ecalMultiFitUncalibRecHit | ||
## | ||
|
||
process.ecalMultiFitUncalibRecHit.algoPSet = cms.PSet( | ||
# for crossCorrelationMethod | ||
timealgo = cms.string( "crossCorrelationMethod" ), | ||
crossCorrelationStartTime = cms.double(-25), | ||
crossCorrelationStopTime = cms.double(25), | ||
crossCorrelationTargetTimePrecision = cms.double(0.01), | ||
) | ||
|
||
## | ||
|
||
process.ecalDigis = process.ecalEBunpacker.clone() | ||
process.ecalDigis.InputLabel = cms.InputTag('rawDataCollector') | ||
|
||
process.out = cms.OutputModule( | ||
"PoolOutputModule", | ||
fileName = cms.untracked.string("test_uncalib.root") | ||
) | ||
|
||
process.finalize = cms.EndPath(process.out) | ||
|
||
process.bunchSpacing = cms.Path( | ||
process.bunchSpacingProducer | ||
) | ||
|
||
process.digiPath = cms.Path( | ||
process.ecalDigis | ||
) | ||
|
||
process.recoPath = cms.Path( | ||
process.ecalMultiFitUncalibRecHit | ||
*process.ecalRecHit | ||
) | ||
|
||
process.schedule = cms.Schedule( | ||
process.bunchSpacing, | ||
process.digiPath, | ||
process.recoPath, | ||
process.finalize | ||
) | ||
|
||
process.options = cms.untracked.PSet( | ||
numberOfThreads = cms.untracked.uint32(8), | ||
numberOfStreams = cms.untracked.uint32(8), | ||
SkipEvent = cms.untracked.vstring('ProductNotFound'), | ||
wantSummary = cms.untracked.bool(True) | ||
) | ||
|
||
|