Skip to content

Commit

Permalink
CPV fix simulation workflow (#5887)
Browse files Browse the repository at this point in the history 
* Fixed CPV simulation workflow

* applied clang-format

* fix for ~statement should be inside braces~ mumbling

Co-authored-by: sevdokim <[email protected]>
  • Loading branch information
@sevdokim
sevdokim and sevdokim authored Apr 11, 2021
1 parent bbc7397 commit e059886
Show file tree
Hide file tree
Showing 15 changed files with 359 additions and 84 deletions.
40 changes: 21 additions & 19 deletions Detectors/CPV/base/include/CPVBase/CPVSimParams.h
View file
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
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "TObject.h"

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
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
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
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
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
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
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
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
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
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

0 comments on commit e059886

Please sign in to comment.