Skip to content

Commit

Permalink
Merge pull request #41027 from fabiocos/fc-BTLamplitude
Browse files Browse the repository at this point in the history 
MTD digitization: add position dependence to BTL amplitudes
  • Loading branch information
@cmsbuild
cmsbuild committed Mar 20, 2023
2 parents 276a056 + e580431 commit ce98e5b
Show file tree
Hide file tree
Showing 7 changed files with 121 additions and 49 deletions.
1 change: 1 addition & 0 deletions SimFastTiming/FastTimingCommon/BuildFile.xml
View file
Expand Up @@ -9,6 +9,7 @@
<use name="CLHEP"/>
<use name="CommonTools/Utils"/>
<use name="vdt_headers"/>
<use name="rootmath"/>
<export>
<lib name="1"/>
</export>
1 change: 1 addition & 0 deletions SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h
View file
Expand Up @@ -46,6 +46,7 @@ class BTLDeviceSim {
const float LightCollSlopeR_;
const float LightCollSlopeL_;
const float PDE_;
const float LCEpositionSlope_;
};

#endif
1 change: 1 addition & 0 deletions SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h
View file
Expand Up @@ -63,6 +63,7 @@ class BTLElectronicsSim {
const bool smearTimeForOOTtails_;
const float Npe_to_pC_;
const float Npe_to_V_;
const std::vector<double> sigmaRelTOFHIRenergy_;

// adc/tdc bitwidths
const uint32_t adcNbits_, tdcNbits_;
Expand Down
2 changes: 2 additions & 0 deletions SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py
View file
Expand Up @@ -15,6 +15,7 @@
LightCollectionSlopeR = cms.double(0.075), # [ns/cm]
LightCollectionSlopeL = cms.double(0.075), # [ns/cm]
PhotonDetectionEff = cms.double(0.20),
LCEpositionSlope = cms.double(0.071), # [1/cm] LCE variation vs longitudinal position shift
),
ElectronicsSimulation = cms.PSet(
bxTime = cms.double(25), # [ns]
Expand All @@ -36,6 +37,7 @@
SmearTimeForOOTtails = cms.bool(True),
Npe_to_pC = cms.double(0.016), # [pC]
Npe_to_V = cms.double(0.0064),# [V]
SigmaRelTOFHIRenergy = cms.vdouble(0.139,-4.35e-05,3.315e-09,-1.20e-13,1.67e-18), # [%] coefficients of 4th degree Chebyshev polynomial parameterization

# n bits for the ADC
adcNbits = cms.uint32(10),
Expand Down
7 changes: 4 additions & 3 deletions SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc
View file
Expand Up @@ -22,7 +22,8 @@ BTLDeviceSim::BTLDeviceSim(const edm::ParameterSet& pset, edm::ConsumesCollector
LightCollEff_(pset.getParameter<double>("LightCollectionEff")),
LightCollSlopeR_(pset.getParameter<double>("LightCollectionSlopeR")),
LightCollSlopeL_(pset.getParameter<double>("LightCollectionSlopeL")),
PDE_(pset.getParameter<double>("PhotonDetectionEff")) {}
PDE_(pset.getParameter<double>("PhotonDetectionEff")),
LCEpositionSlope_(pset.getParameter<double>("LCEpositionSlope")) {}

void BTLDeviceSim::getEventSetup(const edm::EventSetup& evs) {
geom_ = &evs.getData(geomToken_);
Expand Down Expand Up @@ -97,7 +98,7 @@ void BTLDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t, f
// --- Right side
if (iBXR > 0 && iBXR < mtd_digitizer::kNumberOfBX) {
// Accumulate the energy of simHits in the same crystal
(simHitIt->second).hit_info[0][iBXR] += Npe;
(simHitIt->second).hit_info[0][iBXR] += Npe * (1. + LCEpositionSlope_ * convertMmToCm(hit.localPosition().x()));

// Store the time of the first SimHit in the i-th BX
if ((simHitIt->second).hit_info[1][iBXR] == 0 || tR < (simHitIt->second).hit_info[1][iBXR])
Expand All @@ -107,7 +108,7 @@ void BTLDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t, f
// --- Left side
if (iBXL > 0 && iBXL < mtd_digitizer::kNumberOfBX) {
// Accumulate the energy of simHits in the same crystal
(simHitIt->second).hit_info[2][iBXL] += Npe;
(simHitIt->second).hit_info[2][iBXL] += Npe * (1. - LCEpositionSlope_ * convertMmToCm(hit.localPosition().x()));

// Store the time of the first SimHit in the i-th BX
if ((simHitIt->second).hit_info[3][iBXL] == 0 || tL < (simHitIt->second).hit_info[3][iBXL])
Expand Down
15 changes: 14 additions & 1 deletion SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc
View file
Expand Up @@ -6,6 +6,8 @@
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandGaussQ.h"

#include "Math/ChebyshevPol.h"

using namespace mtd;

BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
Expand All @@ -27,6 +29,7 @@ BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset, edm::Consume
smearTimeForOOTtails_(pset.getParameter<bool>("SmearTimeForOOTtails")),
Npe_to_pC_(pset.getParameter<double>("Npe_to_pC")),
Npe_to_V_(pset.getParameter<double>("Npe_to_V")),
sigmaRelTOFHIRenergy_(pset.getParameter<std::vector<double>>("SigmaRelTOFHIRenergy")),
adcNbits_(pset.getParameter<uint32_t>("adcNbits")),
tdcNbits_(pset.getParameter<uint32_t>("tdcNbits")),
adcSaturation_MIP_(pset.getParameter<double>("adcSaturation_MIP")),
Expand Down Expand Up @@ -141,7 +144,17 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
finalToA1 += cosPhi_ * smearing_thr1_uncorr + sinPhi_ * smearing_thr2_uncorr;
finalToA2 += sinPhi_ * smearing_thr1_uncorr + cosPhi_ * smearing_thr2_uncorr;

chargeColl[iside] = Npe * Npe_to_pC_; // the p.e. number is here converted to pC
//Smear the energy according to TOFHIR energy branch measured resolution
float tofhir_ampnoise_relsigma = ROOT::Math::Chebyshev4(Npe,
sigmaRelTOFHIRenergy_[0],
sigmaRelTOFHIRenergy_[1],
sigmaRelTOFHIRenergy_[2],
sigmaRelTOFHIRenergy_[3],
sigmaRelTOFHIRenergy_[4]);
float smearing_tofhir = CLHEP::RandGaussQ::shoot(hre, 0., tofhir_ampnoise_relsigma);
// the amplitude resolution already includes the photostatistics fluctuation, use the original average deposit
chargeColl[iside] = (it->second).hit_info[2 * iside][iBX] * Npe_to_pC_ *
(1. + smearing_tofhir); // the p.e. number is here converted to pC

toa1[iside] = finalToA1;
toa2[iside] = finalToA2;
Expand Down
143 changes: 98 additions & 45 deletions Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc
View file
Expand Up @@ -110,10 +110,11 @@ class BtlLocalRecoValidation : public DQMEDAnalyzer {
MonitorElement* meHitTvsEta_;
MonitorElement* meHitTvsZ_;
MonitorElement* meHitLongPos_;
MonitorElement* meHitLongPosErr_;

MonitorElement* meTimeRes_;
MonitorElement* meTimeResVsE_;
MonitorElement* meEnergyRes_;
MonitorElement* meEnergyRelResVsE_;

MonitorElement* meLongPosPull_;
MonitorElement* meLongPosPullvsE_;
Expand Down Expand Up @@ -179,6 +180,11 @@ class BtlLocalRecoValidation : public DQMEDAnalyzer {

// --- UncalibratedRecHits histograms

MonitorElement* meUncEneLVsX_;
MonitorElement* meUncEneRVsX_;
MonitorElement* meUncTimeLVsX_;
MonitorElement* meUncTimeRVsX_;

static constexpr int nBinsQ_ = 20;
static constexpr float binWidthQ_ = 30.;
static constexpr int nBinsQEta_ = 3;
Expand Down Expand Up @@ -212,9 +218,8 @@ BtlLocalRecoValidation::BtlLocalRecoValidation(const edm::ParameterSet& iConfig)
mtdtopoToken_(esConsumes<MTDTopology, MTDTopologyRcd>()),
cpeToken_(esConsumes<MTDClusterParameterEstimator, MTDCPERecord>(edm::ESInputTag("", "MTDCPEBase"))) {
btlRecHitsToken_ = consumes<FTLRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsTag"));
if (uncalibRecHitsPlots_)
btlUncalibRecHitsToken_ =
consumes<FTLUncalibratedRecHitCollection>(iConfig.getParameter<edm::InputTag>("uncalibRecHitsTag"));
btlUncalibRecHitsToken_ =
consumes<FTLUncalibratedRecHitCollection>(iConfig.getParameter<edm::InputTag>("uncalibRecHitsTag"));
btlSimHitsToken_ = consumes<CrossingFrame<PSimHit> >(iConfig.getParameter<edm::InputTag>("simHitsTag"));
btlRecCluToken_ = consumes<FTLClusterCollection>(iConfig.getParameter<edm::InputTag>("recCluTag"));
mtdTrackingHitToken_ = consumes<MTDTrackingDetSetVector>(iConfig.getParameter<edm::InputTag>("trkHitTag"));
Expand Down Expand Up @@ -303,7 +308,6 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
meHitTime_->Fill(recHit.time());
meHitTimeError_->Fill(recHit.timeError());
meHitLongPos_->Fill(recHit.position());
meHitLongPosErr_->Fill(recHit.positionError());

meOccupancy_->Fill(global_point.z(), global_point.phi());

Expand Down Expand Up @@ -339,7 +343,9 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
const auto& global_point_sim = thedet->toGlobal(local_point_sim);

meTimeRes_->Fill(time_res);
meTimeResVsE_->Fill(recHit.energy(), time_res);
meEnergyRes_->Fill(energy_res);
meEnergyRelResVsE_->Fill(recHit.energy(), energy_res / recHit.energy());

meLongPosPull_->Fill(longpos_res / recHit.positionError());
meLongPosPullvsEta_->Fill(std::abs(global_point_sim.eta()), longpos_res / recHit.positionError());
Expand Down Expand Up @@ -565,16 +571,49 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
meNevents_->Fill(0.5);

// --- Loop over the BTL Uncalibrated RECO hits
if (uncalibRecHitsPlots_) {
auto btlUncalibRecHitsHandle = makeValid(iEvent.getHandle(btlUncalibRecHitsToken_));
auto btlUncalibRecHitsHandle = makeValid(iEvent.getHandle(btlUncalibRecHitsToken_));

for (const auto& uRecHit : *btlUncalibRecHitsHandle) {
BTLDetId detId = uRecHit.id();
for (const auto& uRecHit : *btlUncalibRecHitsHandle) {
BTLDetId detId = uRecHit.id();

// --- Skip UncalibratedRecHits not matched to SimHits
if (m_btlSimHits.count(detId.rawId()) != 1)
continue;
// --- Skip UncalibratedRecHits not matched to SimHits
if (m_btlSimHits.count(detId.rawId()) != 1)
continue;

// --- Combine the information from the left and right BTL cell sides

float nHits = 0.;
float hit_amplitude = 0.;
float hit_time = 0.;

// left side:
if (uRecHit.amplitude().first > 0.) {
hit_amplitude += uRecHit.amplitude().first;
hit_time += uRecHit.time().first;
nHits += 1.;
}
// right side:
if (uRecHit.amplitude().second > 0.) {
hit_amplitude += uRecHit.amplitude().second;
hit_time += uRecHit.time().second;
nHits += 1.;
}

hit_amplitude /= nHits;
hit_time /= nHits;

// --- Fill the histograms

if (hit_amplitude < hitMinAmplitude_)
continue;

meUncEneRVsX_->Fill(uRecHit.position(), uRecHit.amplitude().first - hit_amplitude);
meUncEneLVsX_->Fill(uRecHit.position(), uRecHit.amplitude().second - hit_amplitude);

meUncTimeRVsX_->Fill(uRecHit.position(), uRecHit.time().first - hit_time);
meUncTimeLVsX_->Fill(uRecHit.position(), uRecHit.time().second - hit_time);

if (uncalibRecHitsPlots_) {
DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
const MTDGeomDet* thedet = geom->idToDet(geoId);
if (thedet == nullptr)
Expand All @@ -587,33 +626,6 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
const auto& global_point = thedet->toGlobal(local_point);

// --- Combine the information from the left and right BTL cell sides

float nHits = 0.;
float hit_amplitude = 0.;
float hit_time = 0.;

// left side:
if (uRecHit.amplitude().first > 0.) {
hit_amplitude += uRecHit.amplitude().first;
hit_time += uRecHit.time().first;
nHits += 1.;
}
// right side:
if (uRecHit.amplitude().second > 0.) {
hit_amplitude += uRecHit.amplitude().second;
hit_time += uRecHit.time().second;
nHits += 1.;
}

hit_amplitude /= nHits;
hit_time /= nHits;

// --- Fill the histograms

if (hit_amplitude < hitMinAmplitude_)
continue;

float time_res = hit_time - m_btlSimHits[detId.rawId()].time;

// amplitude histograms
Expand Down Expand Up @@ -645,9 +657,8 @@ void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventS
qBin = ibin;

meTimeResEtavsQ_[etaBin][qBin]->Fill(time_res);

} // uRecHit loop
}
}
} // uRecHit loop}
}

// ------------ method for histogram booking ------------
Expand Down Expand Up @@ -692,11 +703,19 @@ void BtlLocalRecoValidation::bookHistograms(DQMStore::IBooker& ibook,
ibook.bookProfile("BtlHitTvsEta", "BTL RECO ToA vs #eta;#eta_{RECO};ToA_{RECO} [ns]", 50, -1.6, 1.6, 0., 100.);
meHitTvsZ_ =
ibook.bookProfile("BtlHitTvsZ", "BTL RECO ToA vs Z;Z_{RECO} [cm];ToA_{RECO} [ns]", 50, -260., 260., 0., 100.);
meHitLongPos_ = ibook.book1D("BtlLongPos", "BTL RECO hits longitudinal position;long. pos._{RECO}", 100, -10, 10);
meHitLongPosErr_ =
ibook.book1D("BtlLongPosErr", "BTL RECO hits longitudinal position error; long. pos. error_{RECO}", 100, -1, 1);
meHitLongPos_ = ibook.book1D("BtlLongPos", "BTL RECO hits longitudinal position;long. pos._{RECO}", 50, -5, 5);
meTimeRes_ = ibook.book1D("BtlTimeRes", "BTL time resolution;T_{RECO}-T_{SIM}", 100, -0.5, 0.5);
meTimeResVsE_ = ibook.bookProfile(
"BtlTimeResvsE", "BTL time resolution vs hit energy;E_{RECO} [MeV];T_{RECO}-T_{SIM}", 50, 0., 20., -0.5, 0.5, "S");
meEnergyRes_ = ibook.book1D("BtlEnergyRes", "BTL energy resolution;E_{RECO}-E_{SIM}", 100, -0.5, 0.5);
meEnergyRelResVsE_ = ibook.bookProfile("BtlEnergyRelResvsE",
"BTL relative energy resolution vs hit energy;E_{RECO} [MeV];E_{RECO}-E_{SIM}",
50,
0.,
20.,
-0.15,
0.15,
"S");
meLongPosPull_ = ibook.book1D("BtlLongPosPull",
"BTL longitudinal position pull;X^{loc}_{RECO}-X^{loc}_{SIM}/#sigma_{xloc_{RECO}}",
100,
Expand Down Expand Up @@ -901,6 +920,40 @@ void BtlLocalRecoValidation::bookHistograms(DQMStore::IBooker& ibook,

// --- UncalibratedRecHits histograms

meUncEneLVsX_ = ibook.bookProfile("BTLUncEneLVsX",
"BTL uncalibrated hit amplitude left - average vs X;X [cm];Delta(Q_{left}) [pC]",
20,
-5.,
5.,
-640.,
640.,
"S");
meUncEneRVsX_ = ibook.bookProfile("BTLUncEneRVsX",
"BTL uncalibrated hit amplitude right - average vs X;X [cm];Delta(Q_{right}) [pC]",
20,
-5.,
5.,
-640.,
640.,
"S");

meUncTimeLVsX_ = ibook.bookProfile("BTLUncTimeLVsX",
"BTL uncalibrated hit time left - average vs X;X [cm];Delta(T_{left}) [MeV]",
20,
-5.,
5.,
-25.,
25.,
"S");
meUncTimeRVsX_ = ibook.bookProfile("BTLUncTimeRVsX",
"BTL uncalibrated hit time right - average vs X;X [cm];Delta(T_{right}) [MeV]",
20,
-5.,
5.,
-25.,
25.,
"S");

if (uncalibRecHitsPlots_) {
for (unsigned int ihistoQ = 0; ihistoQ < nBinsQ_; ++ihistoQ) {
std::string hname = Form("TimeResQ_%d", ihistoQ);
Expand Down

0 comments on commit ce98e5b

Please sign in to comment.