Skip to content

CPV fix simulation workflow #5887

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.

Sign up for GitHub

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

Merged
davidrohr merged 3 commits into from
Apr 11, 2021
Merged

CPV fix simulation workflow #5887

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8e4543a
Fixed CPV simulation workflow
sevdokim Apr 8, 2021
5acd63c
applied clang-format
sevdokim Apr 8, 2021
3e937b0
fix for ~statement should be inside braces~ mumbling
sevdokim Apr 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 21 additions & 19 deletions Detectors/CPV/base/include/CPVBase/CPVSimParams.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,34 +30,36 @@ struct CPVSimParams : public o2::conf::ConfigurableParamHelper<CPVSimParams> {
float mPadSizeX = 1.13; ///< overall size of CPV active size
float mPadSizeZ = 2.1093; ///< in phi and z directions
float mDetR = 0.1; ///< Relative energy fluctuation in track for 100 e-
float mdEdx = 4.0; ///< Average energy loss in CPV;
float mdEdx = 400.0; ///< Average energy loss in CPV (arbitrary units);
int mNgamz = 5; ///< Ionization size in Z
int mNgamx = 9; ///< Ionization size in Phi
float mCPVGasThickness = 1.3; ///< width of ArC02 gas gap
float mA = 1.0; ///< Parameter to model CPV response
float mB = 0.7; ///< Parameter to model CPV response

//Parameters used in electronic noise calculation and thresholds (Digitizer)
float mReadoutTime = 5.; ///< Read-out time in ns for default simulaionts
float mDeadTime = 20.; ///< PHOS dead time (includes Read-out time i.e. mDeadTime>=mReadoutTime)
float mReadoutTimePU = 2000.; ///< Read-out time in ns if pileup simulation on in DigitizerSpec
float mDeadTimePU = 30000.; ///< PHOS dead time if pileup simulation on in DigitizerSpec
bool mApplyDigitization = true; ///< if energy digitization should be applied
float mZSthreshold = 0.01; ///< Zero Suppression threshold
float mADCWidth = 0.005; ///< Widht of ADC channel used for energy digitization
float mNoise = 0.01; ///< charge noise in one pad
float mCoeffToNanoSecond = 1.e+9; ///< Conversion for time units
float mSortingDelta = 0.1; ///< used in sorting clusters inverse sorting band in cm
float mReadoutTime = 5.; ///< Read-out time in ns for default simulaionts
float mDeadTime = 20.; ///< PHOS dead time (includes Read-out time i.e. mDeadTime>=mReadoutTime)
float mReadoutTimePU = 2000.; ///< Read-out time in ns if pileup simulation on in DigitizerSpec
float mDeadTimePU = 30000.; ///< PHOS dead time if pileup simulation on in DigitizerSpec
bool mApplyDigitization = true; ///< if energy digitization should be applied
//float mZSthreshold = 0.01; ///< Zero Suppression threshold
float mZSnSigmas = 3.; ///< Zero Suppression threshold
//float mADCWidth = 0.005; ///< Widht of ADC channel used for energy digitization
//float mNoise = 0.01; ///< charge noise in one pad
//float mCoeffToNanoSecond = 1.e+9; ///< Conversion for time units
float mSortingDelta = 0.1; ///< used in sorting clusters inverse sorting band in cm

//Parameters used in clusterization
float mDigitMinEnergy = 0.01; ///< Minimal amplitude of a digit to be used in cluster
float mClusteringThreshold = 0.050; ///< Seed digit minimal amplitude
float mUnfogingEAccuracy = 1.e-3; ///< Accuracy of energy calculation in unfoding prosedure (GeV)
float mUnfogingXZAccuracy = 1.e-1; ///< Accuracy of position calculation in unfolding procedure (cm)
float mLocalMaximumCut = 0.030; ///< Threshold to separate local maxima
float mLogWeight = 4.5; ///< weight in cluster center of gravity calculation
int mNMaxIterations = 10; ///< Maximal number of iterations in unfolding procedure
bool mUnfoldClusters = false; ///< Perform cluster unfolding?
//float mDigitMinEnergy = 0.01; ///< Minimal amplitude of a digit to be used in cluster
float mDigitMinEnergy = 5.; ///< Minimal amplitude of a digit to be used in cluster
float mClusteringThreshold = 10.; ///< Seed digit minimal amplitude
float mUnfogingEAccuracy = 1.e-3; ///< Accuracy of energy calculation in unfoding prosedure (GeV)
float mUnfogingXZAccuracy = 1.e-1; ///< Accuracy of position calculation in unfolding procedure (cm)
float mLocalMaximumCut = 0.030; ///< Threshold to separate local maxima
float mLogWeight = 4.5; ///< weight in cluster center of gravity calculation
int mNMaxIterations = 10; ///< Maximal number of iterations in unfolding procedure
bool mUnfoldClusters = false; ///< Perform cluster unfolding?

inline float CellWr() const { return 0.5 * mPadSizeX; } ///< Distance between wires (2 wires above 1 pad)

Expand Down
13 changes: 9 additions & 4 deletions Detectors/CPV/calib/include/CPVCalib/Pedestals.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "TObject.h"
Copy link
Collaborator

@shahor02 shahor02 Apr 9, 2021

Choose a reason for hiding this comment

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

If this is a calibration data format which should be stored in the ccdb, please consider moving it to DataFormats/Detectorts/CPV (can be other PR, I guess there are other such classes).

Copy link
Contributor Author

@sevdokim sevdokim Apr 9, 2021

Choose a reason for hiding this comment

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

@shahor02 actually there is plenty of work: I need to adjust calibration procedure as well. I think configuring ccdb path should naturally be included in that PR (which is not ready yet). I hope it's OK for you.

Copy link
Collaborator

@shahor02 shahor02 Apr 9, 2021

Choose a reason for hiding this comment

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

@sevdokim Sure, but in the DPL workflows keep it simple: eventually the DPL framework will take care of CCDB related stuff.

Copy link
Contributor Author

@sevdokim sevdokim Apr 10, 2021

Choose a reason for hiding this comment

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

@shahor02 to be honest, I didn't get you last sentence. Can you please detalize it?
Particularly, what is the reason to move calibration objects headers to DataFormats/Detectorts/CPV? Will it help CCDB framework somehow to take care of them? What is CCDB related stuff? In our code we connect to CCDB and exctract
calibration objects from there. Maybe there is a global instance of CCDB which can be used and we don't have to create a new one?
PS: sorry for this ignorance. I want to improve my understanding in order to produce effective code and safe your time for reviews in the future.

Copy link
Collaborator

@shahor02 shahor02 Apr 10, 2021

Choose a reason for hiding this comment

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

@sevdokim We try to have all data (including the calibration) format classes (i.e. those classes which ends up in the DPL messages and/or stored in the files) in DataFormats directories, with no dependence on Detectors sim/rec/calib processing classes. Mostly in order to minimize the risk of circular dependencies.

Concerning the "CCDB stuff": eventually, the detector code should not deal with CCDB server directly, instead it will declare in the DPL device specifications what objects it will need (in the same way as it currently subscribes to data) and the framework will take care of shipping the objects matching to the timestamp of the TF being processed.
This is not available at the moment so we have to configure/query the ccdb server. Note that this will work for DPL workflows only. For simple executables, like digi2raw, we will always deal with CCDB directly.


class TH1;
class TH1F;

namespace o2
{
Expand All @@ -45,22 +46,26 @@ class Pedestals
/// \param cellID Absolute ID of cell
/// \return pedestal for the cell
short getPedestal(short cellID) const { return short(mPedestals.at(cellID)); }
float getPedSigma(short cellID) const { return mPedSigmas.at(cellID); }

/// \brief Set pedestal
/// \param cellID Absolute ID of cell
/// \param c is the pedestal (expected to be in range <254)
void setPedestal(short cellID, short c) { mPedestals[cellID] = static_cast<unsigned char>(c); }
void setPedestal(short cellID, short c) { mPedestals[cellID] = (c > 0 && c < 511) ? c : mPedestals[cellID]; }
void setPedSigma(short cellID, float c) { mPedSigmas[cellID] = (c > 0) ? c : mPedSigmas[cellID]; }

/// \brief Set pedestals from 1D histogram with cell absId in x axis
/// \param 1D(NCHANNELS) histogram with calibration coefficients
/// \return Is successful
bool setPedestals(TH1* h);
bool setPedSigmas(TH1F* h);

private:
static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0
std::array<unsigned char, NCHANNELS> mPedestals; ///< Container for pedestals
static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0
std::array<short, NCHANNELS> mPedestals; ///< Container for pedestals
std::array<float, NCHANNELS> mPedSigmas; ///< Container for pedestal sigmas

ClassDefNV(Pedestals, 1);
ClassDefNV(Pedestals, 2);
};

} // namespace cpv
Expand Down
2 changes: 1 addition & 1 deletion Detectors/CPV/calib/src/CalibParams.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ using namespace o2::cpv;
CalibParams::CalibParams(short /*dummy*/)
{
//produce reasonable objest for test purposes
mGainCalib.fill(0.01);
mGainCalib.fill(1.);
}

bool CalibParams::setGain(TH2* h, short module)
Expand Down
38 changes: 33 additions & 5 deletions Detectors/CPV/calib/src/Pedestals.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,15 +11,17 @@
#include "CPVCalib/Pedestals.h"
#include "FairLogger.h"
#include <TH1.h>
#include <TH1F.h>

using namespace o2::cpv;

Pedestals::Pedestals(int /*dummy*/)
{
//produce reasonable objest for test purposes
mPedestals.fill(40);
mPedestals.fill(200); //typical pedestal value
mPedSigmas.fill(1.5); //typical pedestal sigma
}

//______________________________________________________________________________
bool Pedestals::setPedestals(TH1* h)
{
if (!h) {
Expand All @@ -33,11 +35,37 @@ bool Pedestals::setPedestals(TH1* h)
}

for (short i = 1; i <= NCHANNELS; i++) {
if (h->GetBinContent(i) > 255) {
LOG(ERROR) << "pedestal value too large:" << h->GetBinContent(i) << "can not be stored in char";
if (h->GetBinContent(i) > 511) {
LOG(ERROR) << "setPedestals : pedestal value = " << h->GetBinContent(i)
<< " in channel " << i
<< " exceeds max possible value 511 (limited by CPV electronics)";
continue;
}
mPedestals[i] = short(h->GetBinContent(i));
}
return true;
}
//_______________________________________________________________________________
bool Pedestals::setPedSigmas(TH1F* h)
{
if (!h) {
LOG(ERROR) << "no input histogam";
return false;
}

if (h->GetNbinsX() != NCHANNELS) {
LOG(ERROR) << "Wrong dimentions of input histogram:" << h->GetNbinsX() << " instead of " << NCHANNELS;
return false;
}

for (short i = 1; i <= NCHANNELS; i++) {
if (h->GetBinContent(i) < 0) {
LOG(ERROR) << "pedestal sigma = " << h->GetBinContent(i)
<< " in channel " << i
<< " cannot be less than 0";
continue;
}
mPedestals[i] = char(h->GetBinContent(i));
mPedSigmas[i] = float(h->GetBinContent(i));
}
return true;
}
13 changes: 8 additions & 5 deletions Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
#include "DataFormatsCPV/Digit.h"
#include "CPVBase/Geometry.h"
#include "CPVCalib/CalibParams.h"
#include "CPVCalib/Pedestals.h"
#include "CPVCalib/BadChannelMap.h"
#include "CPVBase/Hit.h"
#include "SimulationDataFormat/MCCompLabel.h"
#include "SimulationDataFormat/MCTruthContainer.h"
Expand All @@ -39,15 +41,16 @@ class Digitizer : public TObject
int source, int entry, double dt);

protected:
float uncalibrate(float e, int absId);
float simulateNoise();
float simulatePedestalNoise(int absId);

private:
static constexpr short NCHANNELS = 23040; //128*60*3: toatl number of CPV channels
std::unique_ptr<CalibParams> mCalibParams; /// Calibration coefficients
std::array<Digit, NCHANNELS> mArrayD;

ClassDefOverride(Digitizer, 2);
std::unique_ptr<Pedestals> mPedestals; /// Pedestals
std::unique_ptr<BadChannelMap> mBadMap; /// Bad channel map
std::array<Digit, NCHANNELS> mArrayD; ///array of digits (for inner use)
std::array<float, NCHANNELS> mDigitThresholds;
ClassDefOverride(Digitizer, 3);
};
} // namespace cpv
} // namespace o2
Expand Down
2 changes: 2 additions & 0 deletions Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "DataFormatsCPV/Digit.h"
#include "DataFormatsCPV/TriggerRecord.h"
#include "CPVCalib/CalibParams.h"
#include "CPVCalib/Pedestals.h"

namespace o2
{
Expand Down Expand Up @@ -78,6 +79,7 @@ class RawWriter
FileFor_t mFileFor = FileFor_t::kFullDet; ///< Granularity of the output files
std::string mOutputLocation = "./"; ///< Rawfile name
std::unique_ptr<CalibParams> mCalibParams; ///< CPV calibration
std::unique_ptr<Pedestals> mPedestals; ///< CPV pedestals
std::vector<char> mPayload; ///< Payload to be written
gsl::span<o2::cpv::Digit> mDigits; ///< Digits input vector - must be in digitized format including the time response
std::unique_ptr<o2::raw::RawFileWriter> mRawWriter; ///< Raw writer
Expand Down
1 change: 0 additions & 1 deletion Detectors/CPV/simulation/src/Detector.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,6 @@ Bool_t Detector::ProcessHits(FairVolume* v)

// Now calculate pad response
double qpad = padResponseFunction(qhit, zg, xg);
qpad += cpvparam.mNoise * rnor2;
if (qpad < 0) {
continue;
}
Expand Down
80 changes: 58 additions & 22 deletions Detectors/CPV/simulation/src/Digitizer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ void Digitizer::init()
mCalibParams.reset(new CalibParams(1)); // test default calibration
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
} else {
LOG(INFO) << "[CPVDigitizer] can not gey calibration object from ccdb yet";
LOG(INFO) << "[CPVDigitizer] can not get calibration object from ccdb yet. Using default";
mCalibParams.reset(new CalibParams(1)); // test default calibration
// o2::ccdb::CcdbApi ccdb;
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
Expand All @@ -41,6 +42,45 @@ void Digitizer::init()
// }
}
}
if (!mPedestals) {
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
mPedestals.reset(new Pedestals(1)); // test default calibration
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
} else {
LOG(INFO) << "[CPVDigitizer] can not get pedestal object from ccdb yet. Using default";
mPedestals.reset(new Pedestals(1)); // test default calibration
// o2::ccdb::CcdbApi ccdb;
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
// mPedestals = ccdb.retrieveFromTFileAny<o2::cpv::Pedestals>("CPV/Calib", metadata, mEventTime);
// if (!mPedestals) {
// LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb";
// }
}
}
if (!mBadMap) {
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
mBadMap.reset(new BadChannelMap(1)); // test default calibration
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
} else {
LOG(INFO) << "[CPVDigitizer] can not get bad channel map object from ccdb yet. Using default";
mBadMap.reset(new BadChannelMap(1)); // test default calibration
// o2::ccdb::CcdbApi ccdb;
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
// mBadMap = ccdb.retrieveFromTFileAny<o2::cpv::BadChannelMap>("CPV/Calib", metadata, mEventTime);
// if (!mBadMap) {
// LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb";
// }
}
}

//signal thresolds for digits
//note that digits are calibrated objects
for (int i = 0; i < NCHANNELS; i++) {
mDigitThresholds[i] = o2::cpv::CPVSimParams::Instance().mZSnSigmas *
mPedestals->getPedSigma(i) * mCalibParams->getGain(i);
}
}

//_______________________________________________________________________
Expand All @@ -60,14 +100,11 @@ void Digitizer::processHits(const std::vector<Hit>* hits, const std::vector<Digi
mArrayD[i].reset();
}

if (digitsBg.size() == 0) { // no digits provided: try simulate noise
if (digitsBg.size() == 0) { // no digits provided: try simulate pedestal noise (do it only once)
for (int i = NCHANNELS; i--;) {
float energy = simulateNoise();
energy = uncalibrate(energy, i);
if (energy > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
mArrayD[i].setAmplitude(energy);
mArrayD[i].setAbsId(i);
}
float amplitude = simulatePedestalNoise(i);
mArrayD[i].setAmplitude(amplitude);
mArrayD[i].setAbsId(i);
}
} else { //if digits exist, no noise should be added
for (auto& dBg : digitsBg) { //digits are sorted and unique
Expand All @@ -79,12 +116,12 @@ void Digitizer::processHits(const std::vector<Hit>* hits, const std::vector<Digi
for (auto& h : *hits) {
int i = h.GetDetectorID();
if (mArrayD[i].getAmplitude() > 0) {
mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + uncalibrate(h.GetEnergyLoss(), i));
mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + h.GetEnergyLoss());
} else {
mArrayD[i].setAmplitude(uncalibrate(h.GetEnergyLoss(), i));
mArrayD[i].setAmplitude(h.GetEnergyLoss());
mArrayD[i].setAbsId(i);
}
if (mArrayD[i].getAmplitude() > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) {
int labelIndex = mArrayD[i].getLabel();
if (labelIndex == -1) { //no digit or noisy
labelIndex = labels.getIndexedSize();
Expand All @@ -109,23 +146,22 @@ void Digitizer::processHits(const std::vector<Hit>* hits, const std::vector<Digi
}
}

//finalize output digits
for (int i = 0; i < NCHANNELS; i++) {
if (mArrayD[i].getAmplitude() > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
if (!mBadMap->isChannelGood(i)) {
continue; //bad channel -> skip this digit
}
if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) {
digitsOut.push_back(mArrayD[i]);
}
}
}

float Digitizer::simulateNoise() { return gRandom->Gaus(0., o2::cpv::CPVSimParams::Instance().mNoise); }

//_______________________________________________________________________
float Digitizer::uncalibrate(const float e, const int absId)
float Digitizer::simulatePedestalNoise(int absId)
{
// Decalibrate CPV digit, i.e. transform from amplitude to ADC counts a factor read from CDB
float calib = mCalibParams->getGain(absId);
if (calib > 0) {
return e / calib;
} else {
return 0; // TODO apply de-calibration from OCDB
//this function is to simulate pedestal and its noise (ADC counts)
if (absId < 0 || absId >= NCHANNELS) {
return 0.;
}
return gRandom->Gaus(0, mPedestals->getPedSigma(absId) * mCalibParams->getGain(absId));
}
24 changes: 23 additions & 1 deletion Detectors/CPV/simulation/src/RawWriter.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,28 @@ void RawWriter::digitsToRaw(gsl::span<o2::cpv::Digit> digitsbranch, gsl::span<o2
}
}

if (!mPedestals) {
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
mPedestals = std::make_unique<o2::cpv::Pedestals>(1); // test default calibration
LOG(INFO) << "[RawWriter] No reading calibration from ccdb requested, set default";
} else {
LOG(INFO) << "[RawWriter] getting calibration object from ccdb";
o2::ccdb::CcdbApi ccdb;
std::map<std::string, std::string> metadata;
ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
Copy link
Collaborator

@shahor02 shahor02 Apr 9, 2021

Choose a reason for hiding this comment

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

could you please make the ccdb patch a settable parameter?

auto tr = triggerbranch.begin();
double eventTime = -1;
// if(tr!=triggerbranch.end()){
// eventTime = (*tr).getBCData().getTimeNS() ;
// }
//add copy constructor if necessary
// mPedestals = std::make_unique<o2::cpv::Pedestals>(ccdb.retrieveFromTFileAny<o2::cpv::Pedestals>("CPV/Calib", metadata, eventTime));
if (!mPedestals) {
LOG(FATAL) << "[RawWriter] can not get calibration object from ccdb";
}
}
}

for (auto trg : triggerbranch) {
processTrigger(digitsbranch, trg);
}
Expand All @@ -84,7 +106,7 @@ bool RawWriter::processTrigger(const gsl::span<o2::cpv::Digit> digitsbranch, con
o2::cpv::Geometry::absIdToHWaddress(absId, ccId, dil, gas, pad);

//Convert Amp to ADC counts
short charge = dig.getAmplitude() / mCalibParams->getGain(absId);
short charge = dig.getAmplitude() / mCalibParams->getGain(absId) + mPedestals->getPedestal(absId);
if (charge > 4095) {
charge = 4095;
}
Expand Down
5 changes: 5 additions & 0 deletions Detectors/CPV/testsimulation/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,3 +21,8 @@ o2_add_test_root_macro(plot_hit_cpv.C
o2_add_test_root_macro(plot_dig_cpv.C
PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation
LABELS cpv)

o2_add_test_root_macro(plot_clu_cpv.C
PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation
LABELS cpv)

Loading