Skip to content

Commit 0fa433f

Browse files
sevdokimshahor02
sevdokim
authored and
shahor02
committed
CPV: fix for gain calibration algorithm
1 parent 08c14dd commit 0fa433f

3 files changed

Lines changed: 58 additions & 22 deletions

File tree

‎Detectors/CPV/calib/include/CPVCalibration/GainCalibrator.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@ class AmplitudeSpectrum
4242
double getRMS() const { return (mSumA2 / mNEntries) - ((mSumA * mSumA) / (mNEntries * mNEntries)); }; // return RMS of distribution
4343
uint32_t getNEntries() const { return mNEntries; }
4444
const uint32_t* getBinContent() { return mBinContent.data(); } // since C++17 std::array::data is constexpr
45-
void fillBinData(TH1F* h);
45+
void dumpToHisto(TH1F* h);
46+
int nEventsInRange(float lR, float rR);
4647

4748
private:
4849
uint32_t mNEntries;

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

Lines changed: 52 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -58,21 +58,43 @@ AmplitudeSpectrum& AmplitudeSpectrum::operator+=(const AmplitudeSpectrum& rhs)
5858
void AmplitudeSpectrum::fill(float amplitude)
5959
{
6060
if ((lRange <= amplitude) && (amplitude < rRange)) {
61-
int bin = (amplitude - lRange) / nBins;
61+
int bin = ((amplitude - lRange) / (rRange - lRange)) * nBins;
6262
mBinContent[bin]++;
6363
mNEntries++;
6464
mSumA += amplitude;
6565
mSumA2 += amplitude * amplitude;
6666
}
6767
}
6868
//_____________________________________________________________________________
69-
void AmplitudeSpectrum::fillBinData(TH1F* h)
69+
void AmplitudeSpectrum::dumpToHisto(TH1F* h)
7070
{
71-
for (uint16_t i = 0; i < nBins; i++) {
72-
h->SetBinContent(i + 1, float(mBinContent[i]));
73-
h->SetBinError(i + 1, sqrt(float(mBinContent[i])));
71+
int rebin = nBins / h->GetNbinsX();
72+
if (h != nullptr) {
73+
for (int iBin = 0; iBin < h->GetNbinsX(); iBin++) {
74+
float binC = 0.;
75+
for (uint16_t i = iBin * rebin; i < (iBin + 1) * rebin; i++) {
76+
binC += mBinContent[i];
77+
}
78+
h->SetBinContent(iBin + 1, binC);
79+
}
7480
}
7581
}
82+
int AmplitudeSpectrum::nEventsInRange(float lR, float rR)
83+
{
84+
if (lR < lRange) {
85+
lR = lRange;
86+
}
87+
if (rR > rRange) {
88+
rR = rRange;
89+
}
90+
int first = ((lR - lRange) / (rRange - lRange)) * nBins;
91+
int last = ((rR - lRange) / (rRange - lRange)) * nBins;
92+
int nEvents = 0;
93+
for (int i = first; i < last; i++) {
94+
nEvents += mBinContent[i];
95+
}
96+
return nEvents;
97+
}
7698
//_____________________________________________________________________________
7799
// GainCalibData
78100
//_____________________________________________________________________________
@@ -127,14 +149,16 @@ void GainCalibrator::configParameters()
127149
mMaxAllowedCoeff = cpvParams.gainMaxAllowedCoeff;
128150
mFitRangeL = cpvParams.gainFitRangeL;
129151
mFitRangeR = cpvParams.gainFitRangeR;
152+
mUpdateTFInterval = cpvParams.gainCheckForCalibrationInterval;
153+
130154
// adjust fit ranges to descrete binned values
131155
if (mFitRangeL < AmplitudeSpectrum::lRange) {
132156
mFitRangeL = AmplitudeSpectrum::lRange;
133157
}
134158
if (mFitRangeR > AmplitudeSpectrum::rRange) {
135159
mFitRangeR = AmplitudeSpectrum::rRange;
136160
}
137-
double binWidth = (AmplitudeSpectrum::rRange - AmplitudeSpectrum::lRange) / AmplitudeSpectrum::nBins;
161+
float binWidth = (AmplitudeSpectrum::rRange - AmplitudeSpectrum::lRange) / AmplitudeSpectrum::nBins;
138162
mFitRangeL = AmplitudeSpectrum::lRange + std::floor((mFitRangeL - AmplitudeSpectrum::lRange) / binWidth) * binWidth;
139163
mFitRangeR = AmplitudeSpectrum::lRange + std::floor((mFitRangeR - AmplitudeSpectrum::lRange) / binWidth) * binWidth;
140164

@@ -146,6 +170,7 @@ void GainCalibrator::configParameters()
146170
LOG(info) << "mMaxAllowedCoeff = " << mMaxAllowedCoeff;
147171
LOG(info) << "mFitRangeL = " << mFitRangeL;
148172
LOG(info) << "mFitRangeR = " << mFitRangeR;
173+
LOG(info) << "mUpdateTFInterval = " << mUpdateTFInterval;
149174
}
150175
//_____________________________________________________________________________
151176
void GainCalibrator::initOutput()
@@ -170,21 +195,27 @@ void GainCalibrator::finalizeSlot(GainTimeSlot& slot)
170195
TF1* fLandau = new TF1("fLandau", "landau", mFitRangeL, mFitRangeR);
171196
fLandau->SetParLimits(0, 0., 1.E6);
172197
fLandau->SetParLimits(1, mDesiredLandauMPV / mMaxAllowedCoeff, mDesiredLandauMPV / mMinAllowedCoeff);
173-
fLandau->SetParLimits(1, 0., 1.E3);
198+
fLandau->SetParLimits(2, 0., 1.E3);
174199
TH1F h("histoCPVMaxAmplSpectrum", "", AmplitudeSpectrum::nBins, AmplitudeSpectrum::lRange, AmplitudeSpectrum::rRange);
200+
h.Rebin(std::ceil(mMaxAllowedCoeff)); // rebin histogram in order to avoid descrete structures in it
175201
// double binWidth = (AmplitudeSpectrum::rRange - AmplitudeSpectrum::lRange) / AmplitudeSpectrum::nBins;
176202
// size_t nBinsToFit = (mFitRangeR - mFitRangeL) / binWidth;
177203
// uint32_t xMin = (mFitRangeL - AmplitudeSpectrum::lRange) / binWidth;
178204
// uint32_t xMax = (mFitRangeR - AmplitudeSpectrum::lRange) / binWidth;
179205
int nCalibratedChannels = 0;
180-
for (int i = 0; i < Geometry::kNCHANNELS; i++) {
206+
int badChi2Channels = 0;
207+
for (unsigned short i = 0; i < Geometry::kNCHANNELS; i++) {
208+
newGains->setGain(i, mPreviousGains.get()->getGain(i)); // copy previous gains first
181209
// print some info
182210
if ((i % 500) == 0) {
183211
LOG(info) << "GainCalibrator::finalizeSlot() : checking channel " << i;
184212
}
185-
if (slot.getContainer()->mAmplitudeSpectra[i].getNEntries() > mMinEvents) { // we are ready to fit
213+
if (slot.getContainer()->mAmplitudeSpectra[i].getNEntries() > mMinEvents) { // we are ready to fit
214+
if (slot.getContainer()->mAmplitudeSpectra[i].nEventsInRange(mFitRangeL, mFitRangeR) < mMinEvents) { // actually not enough events in fit range
215+
continue;
216+
}
186217
h.Reset();
187-
slot.getContainer()->mAmplitudeSpectra[i].fillBinData(&h);
218+
slot.getContainer()->mAmplitudeSpectra[i].dumpToHisto(&h);
188219
// set some starting values
189220
double mean = slot.getContainer()->mAmplitudeSpectra[i].getMean();
190221
double rms = slot.getContainer()->mAmplitudeSpectra[i].getRMS();
@@ -193,10 +224,17 @@ void GainCalibrator::finalizeSlot(GainTimeSlot& slot)
193224
auto fitResult = h.Fit(fLandau, "SQL0N", "", mFitRangeL, mFitRangeR);
194225
// auto fitResult = o2::math_utils::fit<uint32_t>(nBinsToFit, &(slot.getContainer()->mAmplitudeSpectra[i].getBinContent()[xMin]), xMin, xMax, *fLandau);
195226
if (fitResult->Chi2() / fitResult->Ndf() > mToleratedChi2PerNDF) {
196-
// in case of bad fit -> do something sofisticated. but what?
227+
badChi2Channels++;
228+
// in case of bad fit -> do something sofisticated. but what?
197229
// continue;
198-
LOG(info) << "GainCalibrator::finalizeSlot() : bad chi2/ndf in fit of spectrum in channel " << i;
199-
fitResult->Print("V");
230+
if (badChi2Channels < 20) {
231+
LOG(info) << "GainCalibrator::finalizeSlot() : bad chi2/ndf in fit of spectrum in channel " << i;
232+
fitResult->Print("V");
233+
} else if (badChi2Channels == 20) {
234+
LOG(info) << "GainCalibrator::finalizeSlot() : bad chi2/ndf in fit of spectrum in channel " << i;
235+
fitResult->Print("V");
236+
LOG(info) << "GainCalibrator::finalizeSlot() : bad chi2/ndf is reported 20 times. Muting ";
237+
}
200238
}
201239
// calib coeffs are defined as mDesiredLandauMPV/actualMPV*previousCoeff
202240
float coeff = mDesiredLandauMPV / fLandau->GetParameter(1) * mPreviousGains.get()->getGain(i);
@@ -213,7 +251,7 @@ void GainCalibrator::finalizeSlot(GainTimeSlot& slot)
213251
}
214252
delete fLandau;
215253
LOG(info) << "GainCalibrator::finalizeSlot() : succesfully calibrated " << nChannelsCalibrated << "channels";
216-
254+
LOG(info) << "GainCalibrator::finalizeSlot() : bad chi2/ndf is occured " << badChi2Channels << "times";
217255
// prepare new ccdb entries for sending
218256
mGainsVec.push_back(*newGains);
219257
// metadata for o2::cpv::CalibParams

‎Detectors/CPV/calib/testWorkflow/GainCalibratorSpec.h‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -38,21 +38,18 @@ class CPVGainCalibratorSpec : public o2::framework::Task
3838
void init(o2::framework::InitContext& ic) final
3939
{
4040
o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest);
41-
// auto slotL = ic.options().get<uint32_t>("tf-per-slot");
42-
// auto delay = ic.options().get<uint32_t>("max-delay");
4341
auto updateInterval = ic.options().get<uint32_t>("updateInterval"); // in TF
4442
bool updateAtTheEndOfRunOnly = ic.options().get<bool>("updateAtTheEndOfRunOnly");
4543
mCalibrator = std::make_unique<o2::cpv::GainCalibrator>();
46-
mCalibrator->setSlotLength(0);
47-
mCalibrator->setMaxSlotsDelay(1000);
44+
mCalibrator->setSlotLength(0); // infinite TF slot
45+
mCalibrator->setMaxSlotsDelay(10000);
4846
if (updateAtTheEndOfRunOnly) {
4947
mCalibrator->setUpdateAtTheEndOfRunOnly();
5048
}
5149
mCalibrator->setCheckIntervalInfiniteSlot(updateInterval);
52-
mCalibrator->setUpdateTFInterval(updateInterval);
5350
LOG(info) << "CPVGainCalibratorSpec initialized";
54-
LOG(info) << "tf-per-slot = 0 (this calibrator works only in single infinite slot mode)";
55-
LOG(info) << "max-delay = 1000 (inconfigurable for this calibrator)";
51+
LOG(info) << "tf-per-slot = 0 (inconfigurable, this calibrator works only in single infinite slot mode)";
52+
LOG(info) << "max-delay = 10000 (inconfigurable for this calibrator)";
5653
LOG(info) << "updateInterval = " << updateInterval;
5754
LOG(info) << "updateAtTheEndOfRunOnly = " << updateAtTheEndOfRunOnly;
5855
}

0 commit comments

Comments
 (0)