Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improved method for EB digi simulation #40855

Merged
merged 3 commits into from
Apr 1, 2023

Conversation

jhakala
Copy link
Contributor

@jhakala jhakala commented Feb 23, 2023

PR description:

This PR implements a new method for simulating EB digis which should be more accurate and thus will improve the quality of the MC for timing studies. The method is based on simulating the response of the APD+readout electronics in a way that takes into account the depth of PCaloHit energy deposits, so pulse shape fluctuations (sometimes referred to as the "Landau term") are modeled. In addition to improving the EB digis, a new collection is added in the EcalTimeDigis which stores a finely-sampled analog waveform, meant to emulate the output of the preamplifier prior to its sampling by the ADC.

In addition to implementing the method, I added a new CondFormat object, such that the "component shapes" used in the method can be accessed via the EventSetup from the CondDB.

The output will change in the following ways:

  • depending on a cfg parameter, a new EBUnsuppressedDigi collection will be added in the DIGI step. (default is that it is omitted.)
  • depending on a cfg parameter, a new array is added to the EcalTimeDigis which has the finely-sampled simulated waveform (it has 10 BX with 25 double values per BX [1 ns sampling] for each EB crystal.) (default is that this is omitted too)

This doesn't depend on any other PRs or externals. Here are some slides on the method: https://indico.cern.ch/event/1187720/#5-status-of-simulation-of-timi
Those slides are slightly out of date but I can provide more documentation if desired.

PR validation:

I've tested this with my own configuration files based on the ecal_devel workflows, in particular runTheMatrix workflow 24234.61, but with customizations added to enable the new functionality. I've looked at the new digis and waveforms from the method and they look as expected (slightly sharper waveforms with more pronounced fluctuations in pulse shapes.) I ran some of the recommended scram b runtests but it looked like they were failing due to unrelated reasons.

@cmsbuild
Copy link
Contributor

-code-checks

Logs: https://cmssdt.cern.ch/SDT/code-checks/cms-sw-PR-40855/34312

  • This PR adds an extra 100KB to repository

  • There are other open Pull requests which might conflict with changes you have proposed:

Code check has found code style and quality issues which could be resolved by applying following patch(s)

@cmsbuild
Copy link
Contributor

+code-checks

Logs: https://cmssdt.cern.ch/SDT/code-checks/cms-sw-PR-40855/34327

  • This PR adds an extra 100KB to repository

  • There are other open Pull requests which might conflict with changes you have proposed:

@cmsbuild
Copy link
Contributor

A new Pull Request was created by @jhakala (John Hakala) for master.

It involves the following packages:

  • CondCore/EcalPlugins (db)
  • CondFormats/DataRecord (db, alca)
  • CondFormats/EcalObjects (db, alca)
  • Configuration/Eras (operations)
  • Configuration/StandardSequences (operations)
  • DataFormats/EcalDigi (simulation)
  • SimCalorimetry/Configuration (simulation)
  • SimCalorimetry/EcalSimAlgos (simulation)
  • SimCalorimetry/EcalSimProducers (simulation)
  • SimGeneral/MixingModule (simulation)

@perrotta, @malbouis, @rappoccio, @yuanchao, @civanch, @mdhildreth, @tvami, @cmsbuild, @saumyaphor4252, @ggovi, @francescobrivio, @ChrisMisan, @fabiocos, @davidlange6 can you please review it and eventually sign? Thanks.
@VourMa, @felicepantaleo, @argiro, @Martin-Grunewald, @thomreis, @slomeo, @makortel, @JanFSchulte, @dgulhan, @missirol, @seemasharmafnal, @simonepigazzini, @GiacomoSguazzoni, @tocheng, @VinInn, @rovere, @mmusich, @mtosi, @fabiocos, @rchatter, @AnnikaStein, @sameasy this is something you requested to watch as well.
@perrotta, @dpiparo, @rappoccio you are the release manager for this.

cms-bot commands are listed here

@civanch
Copy link
Contributor

civanch commented Feb 23, 2023

please test

@cmsbuild
Copy link
Contributor

-1

Failed Tests: HeaderConsistency ClangBuild
Summary: https://cmssdt.cern.ch/SDT/jenkins-artifacts/pull-request-integration/PR-40d5f3/30850/summary.html
COMMIT: e0b676c
CMSSW: CMSSW_13_1_X_2023-02-23-1100/el8_amd64_gcc11
User test area: For local testing, you can use /cvmfs/cms-ci.cern.ch/week1/cms-sw/cmssw/40855/30850/install.sh to create a dev area with all the needed externals and cmssw changes.

Clang Build

I found compilation warning while trying to compile with clang. Command used:

USER_CUDA_FLAGS='--expt-relaxed-constexpr' USER_CXXFLAGS='-Wno-register -fsyntax-only' scram build -k -j 32 COMPILER='llvm compile'

See details on the summary page.

@cmsbuild
Copy link
Contributor

+code-checks

Logs: https://cmssdt.cern.ch/SDT/code-checks/cms-sw-PR-40855/34332

  • This PR adds an extra 36KB to repository

  • There are other open Pull requests which might conflict with changes you have proposed:

@perrotta
Copy link
Contributor

perrotta commented Apr 1, 2023

+1

@perrotta
Copy link
Contributor

perrotta commented Apr 1, 2023

merge

  • The last update was just an innocuous rebase: keep all previously given signatures for good

@cmsbuild cmsbuild merged commit dd08066 into cms-sw:master Apr 1, 2023
@srimanob
Copy link
Contributor

Hi @jhakala
I am wondering

  • if you will provide the follow up PR on the workflow?
  • What is the plan for 13_1 HLT production? Will this method becomes default, or it will be alternative (i.e. requests specific samples with it)

FYI @AdrianoDee @SohamBhattacharya

@srimanob
Copy link
Contributor

srimanob commented Jun 6, 2023

+Upgrade

void setSample(unsigned int i, const float sam) { data_[i] = sam; }
void setSampleOfInterest(int i) { sampleOfInterest_ = i; }

/// Gets the BX==0 sample. If =-1 then it means that only OOT hits are present
int sampleOfInterest() const { return sampleOfInterest_; }
std::vector<float> waveform() const { return waveform_; }

static const unsigned int WAVEFORMSAMPLES = 250;
Copy link
Contributor

@VinInn VinInn Jun 6, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why all capital?
constexpr is preferred.

@@ -2,11 +2,12 @@

namespace {
constexpr unsigned int MAXSAMPLES = 10;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto.


class CaloSubdetectorGeometry;

class EcalTimeMapDigitizer {
public:
struct time_average {
static const unsigned short time_average_capacity = 10;
static const unsigned short time_average_capacity = 10; // this corresponds to the number of BX
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

constexpr

else { // fill with dummy values, since this code is only reached before the actual digi simulation
m_thresh = 0.00013;
time_interval = 1.0;
aVec.reserve(500);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

aVec.resize(500,0.);

for (unsigned int i(0); i != shapeArray.size(); ++i) {
shapeArray[i] = shapeArray[i] / maxelt;
if (normalize) {
for (unsigned int i(0); i != shapeArray.size(); ++i) {
Copy link
Contributor

@VinInn VinInn Jun 6, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for (auto & s: shapeArray) s/=maxelt;

for (unsigned int bin(0); bin != result.waveform_capacity; ++bin) {
if (ComponentShapeCollection::toDepthBin((*it).depth()) <= ComponentShapeCollection::maxDepthBin()) {
result.waveform[bin] +=
(*(shapes()->at((*it).depth())))(binTime - jitter - 25 * (bunchCrossing - m_minBunch)) * (*it).energy();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use operator[] not at for production


private:
bool m_addToBarrel;
bool m_separateDigi;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const ?

void fillCollection(bool useDBShape);
void fillCollection(edm::ConsumesCollector iC);

bool m_useDBShape;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const?


private:
edm::ESGetToken<EcalSimComponentShape, EcalSimComponentShapeRcd> espsToken_;
int shapeIndex_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const ?

double time = (*payload).time_interval;
std::vector<std::vector<float> > EBshapes = (*payload).barrel_shapes;
char nameBuffer[50];
for (auto EBshape : EBshapes) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for (auto const &EBshape : EBshapes)

// --- get the shapes from the user provided txt files
if (!EBSimComponentShapeFiles_.empty()) {
int iShape(0);
for (auto filename : EBSimComponentShapeFiles_) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

auto const & filename

for (auto filename : EBSimComponentShapeFiles_) {
std::ifstream shapeEBFile;
EBshapes.emplace_back();
shapeEBFile.open(EBSimComponentShapeFiles_[iShape].c_str());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shapeEBFile.open(filename.c_str());

shapeEBFile.open(EBSimComponentShapeFiles_[iShape].c_str());
float ww;
while (shapeEBFile >> ww)
EBshapes[iShape].push_back(ww);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

EBshapes.back().push_back(ww);

or auto & EBshape = EBshapes.back() outside the while-loop and then
EBshape.push_back(ww); here

GlobalPoint layerPos =
(cellGeometry)->getPosition(double(layer) + 0.5); //depth in mm in the middle of the layer position
GlobalPoint layerPos = (cellGeometry)->getPosition();
//(cellGeometry)->getPosition(double(layer) + 0.5); //depth in mm in the middle of the layer position // JCH : I am not sure this is doing what it's supposed to, probably unimplemented since CaloCellGeometry returns the same value regardless of this double
return layerPos.mag() * cm / c_light;
Copy link
Contributor

@VinInn VinInn Jun 6, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this looks like a constant: maybe precompute and store it in cellGeometry?

@@ -13,25 +13,32 @@
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "CalibFormats/CaloObjects/interface/CaloTSamplesBase.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "SimCalorimetry/EcalSimAlgos/interface/ComponentShapeCollection.h"

class CaloSubdetectorGeometry;

class EcalTimeMapDigitizer {
public:
struct time_average {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TimeAverage according to CMS naming rules

@@ -66,14 +76,16 @@ class EcalTimeMapDigitizer {

typedef std::vector<unsigned int> VecInd;

explicit EcalTimeMapDigitizer(EcalSubdetector myDet);
explicit EcalTimeMapDigitizer(EcalSubdetector myDet, ComponentShapeCollection* componentShapes);

virtual ~EcalTimeMapDigitizer();
Copy link
Contributor

@VinInn VinInn Jun 6, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why virtual?

}
}

double ComponentShape::timeToRise() const {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

static constexpr
possibly inline

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.