26 #include "TGraphErrors.h" 39 nfitparameters_(iConfig.getUntrackedParameter<
int>(
"numberOfFitParameters", 2)),
40 fitfunction_(iConfig.getUntrackedParameter<
std::
string>(
"fitFunctionRootFormula",
"pol1")),
41 listofdetids_(conf_.getUntrackedParameter<
std::vector<uint32_t> >(
"listOfDetIDs")),
42 ignoreMode_(conf_.getUntrackedParameter<
bool>(
"ignoreMode",
false)),
43 reject_plateaupoints_(iConfig.getUntrackedParameter<
bool>(
"suppressPlateauInFit",
true)),
44 reject_single_entries_(iConfig.getUntrackedParameter<
bool>(
"suppressPointsWithOneEntryOrLess",
true)),
45 plateau_max_slope_(iConfig.getUntrackedParameter<double>(
"plateauSlopeMax", 1.0)),
46 reject_first_point_(iConfig.getUntrackedParameter<
bool>(
"rejectVCalZero",
true)),
47 reject_badpoints_frac_(iConfig.getUntrackedParameter<double>(
"suppressZeroAndPlateausInFitFrac", 0)),
48 bookBIGCalibPayload_(iConfig.getUntrackedParameter<
bool>(
"saveFullPayloads",
false)),
49 savePixelHists_(iConfig.getUntrackedParameter<
bool>(
"savePixelLevelHists",
false)),
50 chi2Threshold_(iConfig.getUntrackedParameter<double>(
"minChi2NDFforHistSave", 10)),
51 chi2ProbThreshold_(iConfig.getUntrackedParameter<double>(
"minChi2ProbforHistSave", 0.05)),
52 maxGainInHist_(iConfig.getUntrackedParameter<double>(
"maxGainInHist", 10)),
53 maxChi2InHist_(iConfig.getUntrackedParameter<double>(
"maxChi2InHist", 25)),
54 saveALLHistograms_(iConfig.getUntrackedParameter<
bool>(
"saveAllHistograms",
false)),
56 filldb_(iConfig.getUntrackedParameter<
bool>(
"writeDatabase",
false)),
57 writeSummary_(iConfig.getUntrackedParameter<
bool>(
"writeSummary",
true)),
58 recordName_(conf_.getParameter<
std::
string>(
"record")),
60 appendMode_(conf_.getUntrackedParameter<
bool>(
"appendMode",
true)),
69 useVcalHigh_(conf_.getParameter<
bool>(
"useVCALHIGH")),
70 scalarVcalHigh_VcalLow_(conf_.getParameter<double>(
"vcalHighToLowConversionFac")) {
75 ::putenv((
char *)
"CORAL_AUTH_USER=me");
76 ::putenv((
char *)
"CORAL_AUTH_PASSWORD=test");
80 graph_ =
new TGraphErrors();
82 summary_.open(
"SummaryPerDetID.txt");
84 for (
int ii = 0;
ii < 10;
ii++)
97 for (
int icol = 1; icol <=
ncols; ++icol) {
100 for (
int irow = 1; irow <= nrows; ++irow) {
105 result.push_back(val);
129 for (
std::map<uint32_t, std::map<std::string, MonitorElement *> >::const_iterator idet =
bookkeeper_.begin();
132 if (detid == idet->first)
138 std::ostringstream summarytext;
141 summarytext <<
"\t Following: values per column: column #, gain, pedestal, chi2\n";
142 for (uint32_t
i = 0;
i < gainvec.size();
i++)
143 summarytext <<
"\t " <<
i <<
" \t" << gainvec[
i] <<
" \t" << pedvec[
i] <<
" \t" << chi2vec[
i] <<
"\n";
144 summarytext <<
"\t list of pixels with high chi2 (chi2> " <<
chi2Threshold_ <<
"): \n";
149 summarytext <<
"\t " << ipix->first <<
"\n";
150 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") << summarytext.str() << std::endl;
156 summary_ <<
"Number of pixel tagged with status :" << endl;
157 for (
int ii = 0;
ii < 9;
ii++)
180 <<
"PLEASE do not fill the database directly from the gain calibration analyzer. This function is currently " 181 "disabled and no DB payloads will be produced!" 186 float lowmeanval = 255;
187 float highmeanval = 0;
191 makehistopersistent =
true;
195 double yerrvals[256];
196 double xvalsall[257];
197 double yvalsall[256];
198 double yerrvalsall[256];
201 bool use_point =
true;
205 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, 0);
211 for (uint32_t
ii = 0;
ii < ipix->getnpoints() &&
ii < 200;
ii++) {
219 yerrvalsall[
ii] = yvalsall[
ii] = 0;
222 yvalsall[
ii] = ipix->getsum(
ii) / (
float)ipix->getnentries(
ii);
223 yerrvalsall[
ii] = ipix->getsumsquares(
ii) / (
float)(ipix->getnentries(
ii));
224 yerrvalsall[
ii] -=
pow(yvalsall[
ii], 2);
225 yerrvalsall[
ii] =
sqrt(yerrvalsall[ii]) /
sqrt(ipix->getnentries(ii));
227 if (yvalsall[ii] < lowmeanval)
228 lowmeanval = yvalsall[
ii];
229 if (yvalsall[ii] > highmeanval)
230 highmeanval = yvalsall[
ii];
235 double plateauval = 0;
236 bool noPlateau =
false;
237 if (nallpoints >= 4) {
238 for (
int ii = nallpoints - 1;
ii > nallpoints - 5; --
ii)
239 plateauval += yvalsall[
ii];
241 for (
int ii = nallpoints - 1;
ii > nallpoints - 5; --
ii) {
242 if (fabs(yvalsall[
ii] - plateauval) > 5) {
249 int NbofPointsInPlateau = 0;
250 for (
int ii = 0;
ii < nallpoints; ++
ii)
251 if (fabs(yvalsall[
ii] - plateauval) < 10 || yvalsall[
ii] == 0)
252 NbofPointsInPlateau++;
254 if (NbofPointsInPlateau >= (nallpoints - 2)) {
256 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
258 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
268 for (
int ii = 0;
ii < nallpoints; ++
ii) {
277 if (yvalsall[
ii] > maxgoodvalinfit && !noPlateau)
279 if (
ii > 1 && fabs(yvalsall[
ii] - yvalsall[
ii - 1]) < 5. && yvalsall[
ii] > 0.8 * maxgoodvalinfit &&
293 float chi2,
slope, intercept,
prob, slopeerror, intercepterror;
295 slope = intercept = slopeerror = intercepterror = 0;
301 for (
int ii = 0;
ii < nallpoints && npoints < 4 && yvalsall[
ii] < plateauval * 0.97; ++
ii) {
302 if (yvalsall[
ii] > 0) {
303 if (
ii > 0 && yvalsall[
ii] - yvalsall[
ii - 1] < 0.1)
314 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
316 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
319 std::ostringstream pixelinfo;
320 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
323 tempname += pixelinfo.str();
326 detid, pixelinfo.str(), tempname, 105 * nallpoints, xvalsall[0], xvalsall[nallpoints - 1] * 1.05);
327 for (
int ii = 0;
ii < nallpoints; ++
ii)
335 func_->SetParameter(0, 50.);
336 func_->SetParameter(1, 0.25);
337 for (
int ipointtemp = 0; ipointtemp <
npoints; ++ipointtemp) {
338 graph_->SetPoint(ipointtemp, xvals[ipointtemp], yvals[ipointtemp]);
339 graph_->SetPointError(ipointtemp, 0, yerrvals[ipointtemp]);
342 slope =
func_->GetParameter(1);
343 slopeerror =
func_->GetParError(1);
344 intercept =
func_->GetParameter(0);
345 intercepterror =
func_->GetParError(0);
347 prob = TMath::Prob(
func_->GetChisquare(), npoints -
func_->GetNpar());
349 while ((
isnan(slope) ||
isnan(intercept)) && ntimes < 10) {
351 makehistopersistent =
true;
353 edm::LogWarning(
"SiPixelGainCalibrationAnalysis") <<
"impossible to fit values, try " << ntimes <<
": ";
356 <<
"vcal " << xvals[
ii] <<
" response: " << yvals[
ii] <<
"+/-" << yerrvals[
ii] << std::endl;
359 slope =
func_->GetParameter(1);
360 slopeerror =
func_->GetParError(1);
361 intercept =
func_->GetParameter(0);
362 intercepterror =
func_->GetParError(0);
364 prob = TMath::Prob(
func_->GetChisquare(), npoints -
func_->GetNpar());
375 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
377 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col() <<
" status_" << status << endl;
393 makehistopersistent =
true;
406 makehistopersistent =
true;
408 <<
"For DETID " << detid <<
"pixel row,col " << ipix->row() <<
"," << ipix->col() <<
" Gain was measured to be " 409 << slope <<
" which is outside the range of the summary plot (" <<
maxGainInHist_ <<
") !!!! " << std::endl;
411 bookkeeper_[detid][
"dynamicrange_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, highmeanval - lowmeanval);
412 bookkeeper_[detid][
"plateau_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, highmeanval);
413 bookkeeper_[detid][
"gain_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
slope);
414 bookkeeper_[detid][
"errorgain_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, slopeerror);
416 bookkeeper_[detid][
"ped_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, intercept);
417 bookkeeper_[detid][
"errorped_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, intercepterror);
419 bookkeeper_[detid][
"chi2_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
chi2);
421 bookkeeper_[detid][
"prob_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
prob);
423 bookkeeper_[detid][
"lowpoint_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, xvals[0]);
424 bookkeeper_[detid][
"highpoint_1d"]->Fill(xvals[npoints - 1]);
425 bookkeeper_[detid][
"highpoint_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1, xvals[npoints - 1]);
426 bookkeeper_[detid][
"nfitpoints_1d"]->Fill(npoints);
427 bookkeeper_[detid][
"endpoint_1d"]->Fill((255 - intercept) * slope);
428 bookkeeper_[detid][
"status_2d"]->setBinContent(ipix->col() + 1, ipix->row() + 1,
status);
434 if (makehistopersistent) {
435 std::ostringstream pixelinfo;
436 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
439 tempname += pixelinfo.str();
446 detid, pixelinfo.str(), tempname, 105 * nallpoints, xvalsall[0], xvalsall[nallpoints - 1] * 1.05);
449 <<
"now saving histogram for pixel " << tempname <<
", gain = " << slope <<
", pedestal = " << intercept
450 <<
", chi2/NDF=" << chi2 <<
"(prob:" << prob <<
"), fit status " <<
status;
451 for (
int ii = 0;
ii < nallpoints; ++
ii) {
459 summary_ <<
"row_" << ipix->row() <<
" col_" << ipix->col();
485 bookDQMHistogram1D(detid,
"GainChi2Prob1d",
"P(#chi^{2},NDOF) for " + tempname, 100, 0., 1.0);
491 detid,
"GainEndPoint1d",
"point where fit meets ADC=255 for " + tempname, 256, 0., 256. *
scalarVcalHigh_VcalLow_);
497 bookDQMHistogram1D(detid,
"GainNPoints1d",
"number of fit point for " + tempname, 20, 0., 20);
499 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_
static const double slope[3]
bool doFits(uint32_t detid, std::vector< SiPixelCalibDigi >::const_iterator ipix) override
std::string translateDetIdToString(uint32_t detid)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
SiPixelGainCalibrationAnalysis(const edm::ParameterSet &iConfig)
MonitorElement * bookDQMHistoPlaquetteSummary2D(uint32_t detid, std::string name, std::string title)
double reject_badpoints_frac_
#define DEFINE_FWK_MODULE(type)
bool reject_plateaupoints_
double scalarVcalHigh_VcalLow_
double chi2ProbThreshold_
void calibrationSetup(const edm::EventSetup &iSetup) override
MonitorElement * bookDQMHistogram1D(uint32_t detid, std::string name, std::string title, int nchX, double lowX, double highX)
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
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)
bool checkCorrectCalibrationType() override
Power< A, B >::type pow(const A &a, const B &b)