Skip to content

Commit e059886

Browse files
sevdokimsevdokim
authored
CPV fix simulation workflow (#5887)
* Fixed CPV simulation workflow * applied clang-format * fix for ~statement should be inside braces~ mumbling Co-authored-by: sevdokim <[email protected]>
1 parent bbc7397 commit e059886

File tree

15 files changed

+359
-84
lines changed

15 files changed

+359
-84
lines changed

‎Detectors/CPV/base/include/CPVBase/CPVSimParams.h

+21-19
Original file line numberDiff line numberDiff line change
@@ -30,34 +30,36 @@ struct CPVSimParams : public o2::conf::ConfigurableParamHelper<CPVSimParams> {
3030
float mPadSizeX = 1.13; ///< overall size of CPV active size
3131
float mPadSizeZ = 2.1093; ///< in phi and z directions
3232
float mDetR = 0.1; ///< Relative energy fluctuation in track for 100 e-
33-
float mdEdx = 4.0; ///< Average energy loss in CPV;
33+
float mdEdx = 400.0; ///< Average energy loss in CPV (arbitrary units);
3434
int mNgamz = 5; ///< Ionization size in Z
3535
int mNgamx = 9; ///< Ionization size in Phi
3636
float mCPVGasThickness = 1.3; ///< width of ArC02 gas gap
3737
float mA = 1.0; ///< Parameter to model CPV response
3838
float mB = 0.7; ///< Parameter to model CPV response
3939

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

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

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

‎Detectors/CPV/calib/include/CPVCalib/Pedestals.h

+9-4
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#include "TObject.h"
2323

2424
class TH1;
25+
class TH1F;
2526

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

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

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

5963
private:
60-
static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0
61-
std::array<unsigned char, NCHANNELS> mPedestals; ///< Container for pedestals
64+
static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0
65+
std::array<short, NCHANNELS> mPedestals; ///< Container for pedestals
66+
std::array<float, NCHANNELS> mPedSigmas; ///< Container for pedestal sigmas
6267

63-
ClassDefNV(Pedestals, 1);
68+
ClassDefNV(Pedestals, 2);
6469
};
6570

6671
} // namespace cpv

‎Detectors/CPV/calib/src/CalibParams.cxx

+1-1
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ using namespace o2::cpv;
2222
CalibParams::CalibParams(short /*dummy*/)
2323
{
2424
//produce reasonable objest for test purposes
25-
mGainCalib.fill(0.01);
25+
mGainCalib.fill(1.);
2626
}
2727

2828
bool CalibParams::setGain(TH2* h, short module)

‎Detectors/CPV/calib/src/Pedestals.cxx

+33-5
Original file line numberDiff line numberDiff line change
@@ -11,15 +11,17 @@
1111
#include "CPVCalib/Pedestals.h"
1212
#include "FairLogger.h"
1313
#include <TH1.h>
14+
#include <TH1F.h>
1415

1516
using namespace o2::cpv;
1617

1718
Pedestals::Pedestals(int /*dummy*/)
1819
{
1920
//produce reasonable objest for test purposes
20-
mPedestals.fill(40);
21+
mPedestals.fill(200); //typical pedestal value
22+
mPedSigmas.fill(1.5); //typical pedestal sigma
2123
}
22-
24+
//______________________________________________________________________________
2325
bool Pedestals::setPedestals(TH1* h)
2426
{
2527
if (!h) {
@@ -33,11 +35,37 @@ bool Pedestals::setPedestals(TH1* h)
3335
}
3436

3537
for (short i = 1; i <= NCHANNELS; i++) {
36-
if (h->GetBinContent(i) > 255) {
37-
LOG(ERROR) << "pedestal value too large:" << h->GetBinContent(i) << "can not be stored in char";
38+
if (h->GetBinContent(i) > 511) {
39+
LOG(ERROR) << "setPedestals : pedestal value = " << h->GetBinContent(i)
40+
<< " in channel " << i
41+
<< " exceeds max possible value 511 (limited by CPV electronics)";
42+
continue;
43+
}
44+
mPedestals[i] = short(h->GetBinContent(i));
45+
}
46+
return true;
47+
}
48+
//_______________________________________________________________________________
49+
bool Pedestals::setPedSigmas(TH1F* h)
50+
{
51+
if (!h) {
52+
LOG(ERROR) << "no input histogam";
53+
return false;
54+
}
55+
56+
if (h->GetNbinsX() != NCHANNELS) {
57+
LOG(ERROR) << "Wrong dimentions of input histogram:" << h->GetNbinsX() << " instead of " << NCHANNELS;
58+
return false;
59+
}
60+
61+
for (short i = 1; i <= NCHANNELS; i++) {
62+
if (h->GetBinContent(i) < 0) {
63+
LOG(ERROR) << "pedestal sigma = " << h->GetBinContent(i)
64+
<< " in channel " << i
65+
<< " cannot be less than 0";
3866
continue;
3967
}
40-
mPedestals[i] = char(h->GetBinContent(i));
68+
mPedSigmas[i] = float(h->GetBinContent(i));
4169
}
4270
return true;
4371
}

‎Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h

+8-5
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,8 @@
1414
#include "DataFormatsCPV/Digit.h"
1515
#include "CPVBase/Geometry.h"
1616
#include "CPVCalib/CalibParams.h"
17+
#include "CPVCalib/Pedestals.h"
18+
#include "CPVCalib/BadChannelMap.h"
1719
#include "CPVBase/Hit.h"
1820
#include "SimulationDataFormat/MCCompLabel.h"
1921
#include "SimulationDataFormat/MCTruthContainer.h"
@@ -39,15 +41,16 @@ class Digitizer : public TObject
3941
int source, int entry, double dt);
4042

4143
protected:
42-
float uncalibrate(float e, int absId);
43-
float simulateNoise();
44+
float simulatePedestalNoise(int absId);
4445

4546
private:
4647
static constexpr short NCHANNELS = 23040; //128*60*3: toatl number of CPV channels
4748
std::unique_ptr<CalibParams> mCalibParams; /// Calibration coefficients
48-
std::array<Digit, NCHANNELS> mArrayD;
49-
50-
ClassDefOverride(Digitizer, 2);
49+
std::unique_ptr<Pedestals> mPedestals; /// Pedestals
50+
std::unique_ptr<BadChannelMap> mBadMap; /// Bad channel map
51+
std::array<Digit, NCHANNELS> mArrayD; ///array of digits (for inner use)
52+
std::array<float, NCHANNELS> mDigitThresholds;
53+
ClassDefOverride(Digitizer, 3);
5154
};
5255
} // namespace cpv
5356
} // namespace o2

‎Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h

+2
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include "DataFormatsCPV/Digit.h"
2727
#include "DataFormatsCPV/TriggerRecord.h"
2828
#include "CPVCalib/CalibParams.h"
29+
#include "CPVCalib/Pedestals.h"
2930

3031
namespace o2
3132
{
@@ -78,6 +79,7 @@ class RawWriter
7879
FileFor_t mFileFor = FileFor_t::kFullDet; ///< Granularity of the output files
7980
std::string mOutputLocation = "./"; ///< Rawfile name
8081
std::unique_ptr<CalibParams> mCalibParams; ///< CPV calibration
82+
std::unique_ptr<Pedestals> mPedestals; ///< CPV pedestals
8183
std::vector<char> mPayload; ///< Payload to be written
8284
gsl::span<o2::cpv::Digit> mDigits; ///< Digits input vector - must be in digitized format including the time response
8385
std::unique_ptr<o2::raw::RawFileWriter> mRawWriter; ///< Raw writer

‎Detectors/CPV/simulation/src/Detector.cxx

-1
Original file line numberDiff line numberDiff line change
@@ -261,7 +261,6 @@ Bool_t Detector::ProcessHits(FairVolume* v)
261261

262262
// Now calculate pad response
263263
double qpad = padResponseFunction(qhit, zg, xg);
264-
qpad += cpvparam.mNoise * rnor2;
265264
if (qpad < 0) {
266265
continue;
267266
}

‎Detectors/CPV/simulation/src/Digitizer.cxx

+58-22
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,8 @@ void Digitizer::init()
3131
mCalibParams.reset(new CalibParams(1)); // test default calibration
3232
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
3333
} else {
34-
LOG(INFO) << "[CPVDigitizer] can not gey calibration object from ccdb yet";
34+
LOG(INFO) << "[CPVDigitizer] can not get calibration object from ccdb yet. Using default";
35+
mCalibParams.reset(new CalibParams(1)); // test default calibration
3536
// o2::ccdb::CcdbApi ccdb;
3637
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
3738
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
@@ -41,6 +42,45 @@ void Digitizer::init()
4142
// }
4243
}
4344
}
45+
if (!mPedestals) {
46+
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
47+
mPedestals.reset(new Pedestals(1)); // test default calibration
48+
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
49+
} else {
50+
LOG(INFO) << "[CPVDigitizer] can not get pedestal object from ccdb yet. Using default";
51+
mPedestals.reset(new Pedestals(1)); // test default calibration
52+
// o2::ccdb::CcdbApi ccdb;
53+
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
54+
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
55+
// mPedestals = ccdb.retrieveFromTFileAny<o2::cpv::Pedestals>("CPV/Calib", metadata, mEventTime);
56+
// if (!mPedestals) {
57+
// LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb";
58+
// }
59+
}
60+
}
61+
if (!mBadMap) {
62+
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
63+
mBadMap.reset(new BadChannelMap(1)); // test default calibration
64+
LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default";
65+
} else {
66+
LOG(INFO) << "[CPVDigitizer] can not get bad channel map object from ccdb yet. Using default";
67+
mBadMap.reset(new BadChannelMap(1)); // test default calibration
68+
// o2::ccdb::CcdbApi ccdb;
69+
// std::map<std::string, std::string> metadata; // do we want to store any meta data?
70+
// ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
71+
// mBadMap = ccdb.retrieveFromTFileAny<o2::cpv::BadChannelMap>("CPV/Calib", metadata, mEventTime);
72+
// if (!mBadMap) {
73+
// LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb";
74+
// }
75+
}
76+
}
77+
78+
//signal thresolds for digits
79+
//note that digits are calibrated objects
80+
for (int i = 0; i < NCHANNELS; i++) {
81+
mDigitThresholds[i] = o2::cpv::CPVSimParams::Instance().mZSnSigmas *
82+
mPedestals->getPedSigma(i) * mCalibParams->getGain(i);
83+
}
4484
}
4585

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

63-
if (digitsBg.size() == 0) { // no digits provided: try simulate noise
103+
if (digitsBg.size() == 0) { // no digits provided: try simulate pedestal noise (do it only once)
64104
for (int i = NCHANNELS; i--;) {
65-
float energy = simulateNoise();
66-
energy = uncalibrate(energy, i);
67-
if (energy > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
68-
mArrayD[i].setAmplitude(energy);
69-
mArrayD[i].setAbsId(i);
70-
}
105+
float amplitude = simulatePedestalNoise(i);
106+
mArrayD[i].setAmplitude(amplitude);
107+
mArrayD[i].setAbsId(i);
71108
}
72109
} else { //if digits exist, no noise should be added
73110
for (auto& dBg : digitsBg) { //digits are sorted and unique
@@ -79,12 +116,12 @@ void Digitizer::processHits(const std::vector<Hit>* hits, const std::vector<Digi
79116
for (auto& h : *hits) {
80117
int i = h.GetDetectorID();
81118
if (mArrayD[i].getAmplitude() > 0) {
82-
mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + uncalibrate(h.GetEnergyLoss(), i));
119+
mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + h.GetEnergyLoss());
83120
} else {
84-
mArrayD[i].setAmplitude(uncalibrate(h.GetEnergyLoss(), i));
121+
mArrayD[i].setAmplitude(h.GetEnergyLoss());
85122
mArrayD[i].setAbsId(i);
86123
}
87-
if (mArrayD[i].getAmplitude() > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
124+
if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) {
88125
int labelIndex = mArrayD[i].getLabel();
89126
if (labelIndex == -1) { //no digit or noisy
90127
labelIndex = labels.getIndexedSize();
@@ -109,23 +146,22 @@ void Digitizer::processHits(const std::vector<Hit>* hits, const std::vector<Digi
109146
}
110147
}
111148

149+
//finalize output digits
112150
for (int i = 0; i < NCHANNELS; i++) {
113-
if (mArrayD[i].getAmplitude() > o2::cpv::CPVSimParams::Instance().mZSthreshold) {
151+
if (!mBadMap->isChannelGood(i)) {
152+
continue; //bad channel -> skip this digit
153+
}
154+
if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) {
114155
digitsOut.push_back(mArrayD[i]);
115156
}
116157
}
117158
}
118159

119-
float Digitizer::simulateNoise() { return gRandom->Gaus(0., o2::cpv::CPVSimParams::Instance().mNoise); }
120-
121-
//_______________________________________________________________________
122-
float Digitizer::uncalibrate(const float e, const int absId)
160+
float Digitizer::simulatePedestalNoise(int absId)
123161
{
124-
// Decalibrate CPV digit, i.e. transform from amplitude to ADC counts a factor read from CDB
125-
float calib = mCalibParams->getGain(absId);
126-
if (calib > 0) {
127-
return e / calib;
128-
} else {
129-
return 0; // TODO apply de-calibration from OCDB
162+
//this function is to simulate pedestal and its noise (ADC counts)
163+
if (absId < 0 || absId >= NCHANNELS) {
164+
return 0.;
130165
}
166+
return gRandom->Gaus(0, mPedestals->getPedSigma(absId) * mCalibParams->getGain(absId));
131167
}

‎Detectors/CPV/simulation/src/RawWriter.cxx

+23-1
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,28 @@ void RawWriter::digitsToRaw(gsl::span<o2::cpv::Digit> digitsbranch, gsl::span<o2
6060
}
6161
}
6262

63+
if (!mPedestals) {
64+
if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) {
65+
mPedestals = std::make_unique<o2::cpv::Pedestals>(1); // test default calibration
66+
LOG(INFO) << "[RawWriter] No reading calibration from ccdb requested, set default";
67+
} else {
68+
LOG(INFO) << "[RawWriter] getting calibration object from ccdb";
69+
o2::ccdb::CcdbApi ccdb;
70+
std::map<std::string, std::string> metadata;
71+
ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation
72+
auto tr = triggerbranch.begin();
73+
double eventTime = -1;
74+
// if(tr!=triggerbranch.end()){
75+
// eventTime = (*tr).getBCData().getTimeNS() ;
76+
// }
77+
//add copy constructor if necessary
78+
// mPedestals = std::make_unique<o2::cpv::Pedestals>(ccdb.retrieveFromTFileAny<o2::cpv::Pedestals>("CPV/Calib", metadata, eventTime));
79+
if (!mPedestals) {
80+
LOG(FATAL) << "[RawWriter] can not get calibration object from ccdb";
81+
}
82+
}
83+
}
84+
6385
for (auto trg : triggerbranch) {
6486
processTrigger(digitsbranch, trg);
6587
}
@@ -84,7 +106,7 @@ bool RawWriter::processTrigger(const gsl::span<o2::cpv::Digit> digitsbranch, con
84106
o2::cpv::Geometry::absIdToHWaddress(absId, ccId, dil, gas, pad);
85107

86108
//Convert Amp to ADC counts
87-
short charge = dig.getAmplitude() / mCalibParams->getGain(absId);
109+
short charge = dig.getAmplitude() / mCalibParams->getGain(absId) + mPedestals->getPedestal(absId);
88110
if (charge > 4095) {
89111
charge = 4095;
90112
}

‎Detectors/CPV/testsimulation/CMakeLists.txt

+5
Original file line numberDiff line numberDiff line change
@@ -21,3 +21,8 @@ o2_add_test_root_macro(plot_hit_cpv.C
2121
o2_add_test_root_macro(plot_dig_cpv.C
2222
PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation
2323
LABELS cpv)
24+
25+
o2_add_test_root_macro(plot_clu_cpv.C
26+
PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation
27+
LABELS cpv)
28+

0 commit comments

Comments
 (0)