24 #include "TGraphErrors.h"
37 nfitparameters_(iConfig.getUntrackedParameter<int>(
"numberOfFitParameters", 2)),
38 fitfunction_(iConfig.getUntrackedParameter<std::
string>(
"fitFunctionRootFormula",
"pol1")),
39 listofdetids_(conf_.getUntrackedParameter<std::
vector<uint32_t> >(
"listOfDetIDs")),
40 ignoreMode_(conf_.getUntrackedParameter<bool>(
"ignoreMode",
false)),
41 reject_plateaupoints_(iConfig.getUntrackedParameter<bool>(
"suppressPlateauInFit",
true)),
42 reject_single_entries_(iConfig.getUntrackedParameter<bool>(
"suppressPointsWithOneEntryOrLess",
true)),
43 plateau_max_slope_(iConfig.getUntrackedParameter<double>(
"plateauSlopeMax", 1.0)),
44 reject_first_point_(iConfig.getUntrackedParameter<bool>(
"rejectVCalZero",
true)),
45 reject_badpoints_frac_(iConfig.getUntrackedParameter<double>(
"suppressZeroAndPlateausInFitFrac", 0)),
46 bookBIGCalibPayload_(iConfig.getUntrackedParameter<bool>(
"saveFullPayloads",
false)),
47 savePixelHists_(iConfig.getUntrackedParameter<bool>(
"savePixelLevelHists",
false)),
48 chi2Threshold_(iConfig.getUntrackedParameter<double>(
"minChi2NDFforHistSave", 10)),
49 chi2ProbThreshold_(iConfig.getUntrackedParameter<double>(
"minChi2ProbforHistSave", 0.05)),
50 maxGainInHist_(iConfig.getUntrackedParameter<double>(
"maxGainInHist", 10)),
51 maxChi2InHist_(iConfig.getUntrackedParameter<double>(
"maxChi2InHist", 25)),
52 saveALLHistograms_(iConfig.getUntrackedParameter<bool>(
"saveAllHistograms",
false)),
54 filldb_(iConfig.getUntrackedParameter<bool>(
"writeDatabase",
false)),
55 writeSummary_(iConfig.getUntrackedParameter<bool>(
"writeSummary",
true)),
56 recordName_(conf_.getParameter<std::
string>(
"record")),
58 appendMode_(conf_.getUntrackedParameter<bool>(
"appendMode",
true)),
67 useVcalHigh_(conf_.getParameter<bool>(
"useVCALHIGH")),
68 scalarVcalHigh_VcalLow_(conf_.getParameter<double>(
"vcalHighToLowConversionFac")) {
73 ::putenv((
char *)
"CORAL_AUTH_USER=me");
74 ::putenv((
char *)
"CORAL_AUTH_PASSWORD=test");
78 graph_ =
new TGraphErrors();
80 summary_.open(
"SummaryPerDetID.txt");
82 for (
int ii = 0;
ii < 10;
ii++)
95 for (
int icol = 1; icol <= ncols; ++icol) {
98 for (
int irow = 1; irow <= nrows; ++irow) {
103 result.push_back(val);
127 for (std::map<uint32_t, std::map<std::string, MonitorElement *> >::const_iterator idet =
bookkeeper_.begin();
130 if (detid == idet->first)
136 std::ostringstream summarytext;
139 summarytext <<
"\t Following: values per column: column #, gain, pedestal, chi2\n";
140 for (uint32_t
i = 0;
i < gainvec.size();
i++)
141 summarytext <<
"\t " <<
i <<
" \t" << gainvec[
i] <<
" \t" << pedvec[
i] <<
" \t" << chi2vec[
i] <<
"\n";
142 summarytext <<
"\t list of pixels with high chi2 (chi2> " <<
chi2Threshold_ <<
"): \n";
147 summarytext <<
"\t " << ipix->first <<
"\n";
148 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") << summarytext.str() << std::endl;
154 summary_ <<
"Number of pixel tagged with status :" << endl;
155 for (
int ii = 0;
ii < 9;
ii++)
178 <<
"PLEASE do not fill the database directly from the gain calibration analyzer. This function is currently "
179 "disabled and no DB payloads will be produced!"
184 float lowmeanval = 255;
185 float highmeanval = 0;
189 makehistopersistent =
true;
193 double yerrvals[256];
194 double xvalsall[257];
195 double yvalsall[256];
196 double yerrvalsall[256];
199 bool use_point =
true;
203 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, 0);
209 for (uint32_t
ii = 0;
ii < ipix->getnpoints() &&
ii < 200;
ii++) {
217 yerrvalsall[
ii] = yvalsall[
ii] = 0;
220 yvalsall[
ii] = ipix->getsum(
ii) / (float)ipix->getnentries(
ii);
221 yerrvalsall[
ii] = ipix->getsumsquares(
ii) / (float)(ipix->getnentries(
ii));
222 yerrvalsall[
ii] -=
pow(yvalsall[
ii], 2);
223 yerrvalsall[
ii] =
sqrt(yerrvalsall[ii]) /
sqrt(ipix->getnentries(ii));
225 if (yvalsall[ii] < lowmeanval)
226 lowmeanval = yvalsall[
ii];
227 if (yvalsall[ii] > highmeanval)
228 highmeanval = yvalsall[
ii];
233 double plateauval = 0;
234 bool noPlateau =
false;
235 if (nallpoints >= 4) {
236 for (
int ii = nallpoints - 1;
ii > nallpoints - 5; --
ii)
237 plateauval += yvalsall[
ii];
239 for (
int ii = nallpoints - 1;
ii > nallpoints - 5; --
ii) {
240 if (fabs(yvalsall[
ii] - plateauval) > 5) {
247 int NbofPointsInPlateau = 0;
248 for (
int ii = 0;
ii < nallpoints; ++
ii)
249 if (fabs(yvalsall[
ii] - plateauval) < 10 || yvalsall[
ii] == 0)
250 NbofPointsInPlateau++;
252 if (NbofPointsInPlateau >= (nallpoints - 2)) {
254 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
256 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
266 for (
int ii = 0;
ii < nallpoints; ++
ii) {
275 if (yvalsall[
ii] > maxgoodvalinfit && !noPlateau)
277 if (
ii > 1 && fabs(yvalsall[
ii] - yvalsall[
ii - 1]) < 5. && yvalsall[
ii] > 0.8 * maxgoodvalinfit &&
291 float chi2,
slope, intercept, prob, slopeerror, intercepterror;
293 slope = intercept = slopeerror = intercepterror = 0;
299 for (
int ii = 0;
ii < nallpoints && npoints < 4 && yvalsall[
ii] < plateauval * 0.97; ++
ii) {
300 if (yvalsall[
ii] > 0) {
301 if (
ii > 0 && yvalsall[
ii] - yvalsall[
ii - 1] < 0.1)
312 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
314 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
317 std::ostringstream pixelinfo;
318 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
321 tempname += pixelinfo.str();
324 detid, pixelinfo.str(), tempname, 105 * nallpoints, xvalsall[0], xvalsall[nallpoints - 1] * 1.05);
325 for (
int ii = 0;
ii < nallpoints; ++
ii)
333 func_->SetParameter(0, 50.);
334 func_->SetParameter(1, 0.25);
335 for (
int ipointtemp = 0; ipointtemp <
npoints; ++ipointtemp) {
336 graph_->SetPoint(ipointtemp, xvals[ipointtemp], yvals[ipointtemp]);
337 graph_->SetPointError(ipointtemp, 0, yerrvals[ipointtemp]);
340 slope =
func_->GetParameter(1);
341 slopeerror =
func_->GetParError(1);
342 intercept =
func_->GetParameter(0);
343 intercepterror =
func_->GetParError(0);
344 chi2 =
func_->GetChisquare() / ((float)npoints -
func_->GetNpar());
345 prob = TMath::Prob(
func_->GetChisquare(), npoints -
func_->GetNpar());
347 while ((
isnan(slope) ||
isnan(intercept)) && ntimes < 10) {
349 makehistopersistent =
true;
351 edm::LogWarning(
"SiPixelGainCalibrationAnalysis") <<
"impossible to fit values, try " << ntimes <<
": ";
354 <<
"vcal " << xvals[
ii] <<
" response: " << yvals[
ii] <<
"+/-" << yerrvals[
ii] << std::endl;
357 slope =
func_->GetParameter(1);
358 slopeerror =
func_->GetParError(1);
359 intercept =
func_->GetParameter(0);
360 intercepterror =
func_->GetParError(0);
361 chi2 =
func_->GetChisquare() / ((float)npoints -
func_->GetNpar());
362 prob = TMath::Prob(
func_->GetChisquare(), npoints -
func_->GetNpar());
373 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
375 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
391 makehistopersistent =
true;
404 makehistopersistent =
true;
406 <<
"For DETID " << detid <<
"pixel row,col " << ipix->row() <<
"," << ipix->col() <<
" Gain was measured to be "
407 << slope <<
" which is outside the range of the summary plot (" <<
maxGainInHist_ <<
") !!!! " << std::endl;
409 bookkeeper_[detid][
"dynamicrange_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, highmeanval - lowmeanval);
410 bookkeeper_[detid][
"plateau_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, highmeanval);
411 bookkeeper_[detid][
"gain_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
slope);
412 bookkeeper_[detid][
"errorgain_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, slopeerror);
414 bookkeeper_[detid][
"ped_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, intercept);
415 bookkeeper_[detid][
"errorped_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, intercepterror);
417 bookkeeper_[detid][
"chi2_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
chi2);
419 bookkeeper_[detid][
"prob_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, prob);
421 bookkeeper_[detid][
"lowpoint_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, xvals[0]);
422 bookkeeper_[detid][
"highpoint_1d"]->Fill(xvals[npoints - 1]);
423 bookkeeper_[detid][
"highpoint_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, xvals[npoints - 1]);
424 bookkeeper_[detid][
"nfitpoints_1d"]->Fill(npoints);
425 bookkeeper_[detid][
"endpoint_1d"]->Fill((255 - intercept) * slope);
426 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
432 if (makehistopersistent) {
433 std::ostringstream pixelinfo;
434 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
437 tempname += pixelinfo.str();
444 detid, pixelinfo.str(), tempname, 105 * nallpoints, xvalsall[0], xvalsall[nallpoints - 1] * 1.05);
447 <<
"now saving histogram for pixel " << tempname <<
", gain = " << slope <<
", pedestal = " << intercept
448 <<
", chi2/NDF=" << chi2 <<
"(prob:" << prob <<
"), fit status " <<
status;
449 for (
int ii = 0;
ii < nallpoints; ++
ii) {
457 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col();
483 bookDQMHistogram1D(detid,
"GainChi2Prob1d",
"P(#chi^{2},NDOF) for " + tempname, 100, 0., 1.0);
489 detid,
"GainEndPoint1d",
"point where fit meets ADC=255 for " + tempname, 256, 0., 256. *
scalarVcalHigh_VcalLow_);
495 bookDQMHistogram1D(detid,
"GainNPoints1d",
"number of fit point for " + tempname, 20, 0., 20);
497 detid,
"GainDynamicRange2d",
"Difference lowest and highest points on gain curve for " + tempname);
std::map< uint32_t, std::map< std::string, MonitorElement * > > bookkeeper_pixels_
bool reject_single_entries_
#define DEFINE_FWK_MODULE(type)
static const double slope[3]
bool doFits(uint32_t detid, std::vector< SiPixelCalibDigi >::const_iterator ipix) override
std::string translateDetIdToString(uint32_t detid)
Log< level::Error, false > LogError
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
SiPixelGainCalibrationAnalysis(const edm::ParameterSet &iConfig)
double reject_badpoints_frac_
MonitorElement * bookDQMHistogram1D(uint32_t detid, std::string name, std::string title, int nchX, double lowX, double highX)
bool reject_plateaupoints_
double scalarVcalHigh_VcalLow_
double chi2ProbThreshold_
void calibrationSetup(const edm::EventSetup &iSetup) override
Abs< T >::type abs(const T &t)
static std::vector< short > vCalValues_
void calibrationEnd() override
std::vector< uint32_t > listofdetids_
~SiPixelGainCalibrationAnalysis() override
void newDetID(uint32_t detid) override
MonitorElement * bookDQMHistoPlaquetteSummary2D(uint32_t detid, std::string name, std::string title)
Log< level::Info, false > LogInfo
std::map< uint32_t, std::map< std::string, MonitorElement * > > bookkeeper_
std::vector< float > CalculateAveragePerColumn(uint32_t detid, std::string label)
std::string calibrationMode_
bool setDQMDirectory(std::string dirName)
Log< level::Warning, false > LogWarning
bool checkCorrectCalibrationType() override
Power< A, B >::type pow(const A &a, const B &b)