CMS 3D CMS Logo

SiPixelSCurveCalibrationAnalysis.cc
Go to the documentation of this file.
2 #include "TMath.h"
3 
4 #include <fstream>
5 #include <iostream>
6 
14 #include <sstream>
15 
16 // initialize static members
19 
23 }
24 
26  std::ofstream myfile;
27  myfile.open(thresholdfilename_.c_str());
28  for (detIDHistogramMap::iterator thisDetIdHistoGrams = histograms_.begin(); thisDetIdHistoGrams != histograms_.end();
29  ++thisDetIdHistoGrams) {
30  // loop over det id (det id = number (unsigned int) of pixel module
31  const MonitorElement *sigmahist = (*thisDetIdHistoGrams).second[kSigmas];
32  const MonitorElement *thresholdhist = (*thisDetIdHistoGrams).second[kThresholds];
33  uint32_t detid = (*thisDetIdHistoGrams).first;
34  std::string name = sigmahist->getTitle();
35  std::string rocname = name.substr(0, name.size() - 7);
36  rocname += "_ROC";
37  int total_rows = sigmahist->getNbinsY();
38  int total_columns = sigmahist->getNbinsX();
39  // loop over all rows on columns on all ROCs
40  for (int irow = 0; irow < total_rows; ++irow) {
41  for (int icol = 0; icol < total_columns; ++icol) {
42  float threshold_error = sigmahist->getBinContent(icol + 1, irow + 1); // +1 because root bins start at 1
43  if (writeZeroes_ || (!writeZeroes_ && threshold_error > 0)) {
44  // changing from offline to online numbers
45  int realfedID = -1;
46  for (int fedid = 0; fedid <= 40; ++fedid) {
48  if (converter.hasDetUnit(detid)) {
49  realfedID = fedid;
50  break;
51  }
52  }
53  if (realfedID == -1) {
54  std::cout << "error: could not obtain real fed ID" << std::endl;
55  }
56  sipixelobjects::DetectorIndex detector = {detid, irow, icol};
58  SiPixelFrameConverter formatter(theCablingMap_.product(), realfedID);
59  formatter.toCabling(cabling, detector);
60  // cabling should now contain cabling.roc and cabling.dcol and
61  // cabling.pxid however, the coordinates now need to be converted from
62  // dcl,pxid to the row,col coordinates used in the calibration info
64  loc.dcol = cabling.dcol;
65  loc.pxid = cabling.pxid;
66  // FIX to adhere to new cabling map. To be replaced with
67  // CalibTracker/SiPixelTools detid - > hardware id classes ASAP.
68  // const sipixelobjects::PixelFEDCabling *theFed=
69  // theCablingMap.product()->fed(realfedID); const
70  // sipixelobjects::PixelFEDLink * link =
71  // theFed->link(cabling.link); const sipixelobjects::PixelROC
72  // *theRoc = link->roc(cabling.roc);
73  sipixelobjects::LocalPixel locpixel(loc);
74  sipixelobjects::CablingPathToDetUnit path = {static_cast<unsigned int>(realfedID),
75  static_cast<unsigned int>(cabling.link),
76  static_cast<unsigned int>(cabling.roc)};
78  // END of FIX
79  int newrow = locpixel.rocRow();
80  int newcol = locpixel.rocCol();
81  myfile << rocname << theRoc->idInDetUnit() << " " << newcol << " " << newrow << " "
82  << thresholdhist->getBinContent(icol + 1, irow + 1) << " "
83  << threshold_error; // +1 because root bins start at 1
84  myfile << "\n";
85  }
86  }
87  }
88  }
89  myfile.close();
90 }
91 
92 // used for TMinuit fitting
93 void chi2toMinimize(int &npar, double *grad, double &fcnval, double *xval, int iflag) {
95  // setup function parameters
96  for (int i = 0; i < npar; i++)
97  theFormula->SetParameter(i, xval[i]);
98  fcnval = 0;
99  // compute Chi2 of all points
100  const std::vector<short> *theVCalValues = SiPixelSCurveCalibrationAnalysis::getVcalValues();
101  for (uint32_t i = 0; i < theVCalValues->size(); i++) {
102  float chi = (SiPixelSCurveCalibrationAnalysis::efficiencies_[i] - theFormula->Eval((*theVCalValues)[i]));
104  fcnval += chi * chi;
105  }
106 }
107 
109  edm::LogInfo("SiPixelSCurveCalibrationAnalysis") << "Setting up calibration paramters.";
110  std::vector<uint32_t> anEmptyDefaultVectorOfUInts;
111  std::vector<uint32_t> detIDsToSaveVector_;
112  useDetectorHierarchyFolders_ = iConfig.getUntrackedParameter<bool>("useDetectorHierarchyFolders", true);
113  saveCurvesThatFlaggedBad_ = iConfig.getUntrackedParameter<bool>("saveCurvesThatFlaggedBad", false);
114  detIDsToSaveVector_ =
115  iConfig.getUntrackedParameter<std::vector<uint32_t>>("detIDsToSave", anEmptyDefaultVectorOfUInts);
116  maxCurvesToSave_ = iConfig.getUntrackedParameter<uint32_t>("maxCurvesToSave", 1000);
117  write2dHistograms_ = iConfig.getUntrackedParameter<bool>("write2dHistograms", true);
118  write2dFitResult_ = iConfig.getUntrackedParameter<bool>("write2dFitResult", true);
119  printoutthresholds_ = iConfig.getUntrackedParameter<bool>("writeOutThresholdSummary", true);
120  thresholdfilename_ = iConfig.getUntrackedParameter<std::string>("thresholdOutputFileName", "thresholds.txt");
121  minimumChi2prob_ = iConfig.getUntrackedParameter<double>("minimumChi2prob", 0);
122  minimumThreshold_ = iConfig.getUntrackedParameter<double>("minimumThreshold", -10);
123  maximumThreshold_ = iConfig.getUntrackedParameter<double>("maximumThreshold", 300);
124  minimumSigma_ = iConfig.getUntrackedParameter<double>("minimumSigma", 0);
125  maximumSigma_ = iConfig.getUntrackedParameter<double>("maximumSigma", 100);
126  minimumEffAsymptote_ = iConfig.getUntrackedParameter<double>("minimumEffAsymptote", 0);
127  maximumEffAsymptote_ = iConfig.getUntrackedParameter<double>("maximumEffAsymptote", 1000);
128  maximumSigmaBin_ = iConfig.getUntrackedParameter<double>("maximumSigmaBin", 10);
129  maximumThresholdBin_ = iConfig.getUntrackedParameter<double>("maximumThresholdBin", 255);
130 
131  writeZeroes_ = iConfig.getUntrackedParameter<bool>("alsoWriteZeroThresholds", false);
132 
133  // convert the vector into a map for quicker lookups.
134  for (unsigned int i = 0; i < detIDsToSaveVector_.size(); i++)
135  detIDsToSave_.insert(std::make_pair(detIDsToSaveVector_[i], true));
136 }
137 
139  // do nothing
140 }
141 
143  const uint32_t &row,
144  const uint32_t &col,
145  sCurveErrorFlag errorFlag,
146  const std::vector<float> &efficiencies,
147  const std::vector<float> &errors) {
149  edm::LogWarning("SiPixelSCurveCalibrationAnalysis")
150  << "WARNING: Request to save curve for [detid](col/row): [" << detid << "](" << col << "/" << row
151  << ") denied. Maximum number of saved curves (defined in .cfi) "
152  "exceeded.";
153  return;
154  }
155  std::ostringstream rootName;
156  rootName << "SCurve_row_" << row << "_col_" << col;
157  std::ostringstream humanName;
158  humanName << translateDetIdToString(detid) << "_" << rootName.str() << "_ErrorFlag_" << (int)errorFlag;
159 
160  unsigned int numberOfVCalPoints = vCalPointsAsFloats_.size() - 1; // minus one is necessary since the lower
161  // edge of the last bin must be added
162  if (efficiencies.size() != numberOfVCalPoints || errors.size() != numberOfVCalPoints) {
163  edm::LogError("SiPixelSCurveCalibrationAnalysis")
164  << "Error saving single curve histogram! Number of Vcal values (" << numberOfVCalPoints
165  << ") does not match number of efficiency points or error points!";
166  return;
167  }
168  setDQMDirectory(detid);
169  float *vcalValuesToPassToCrappyRoot = &vCalPointsAsFloats_[0];
170  MonitorElement *aBadHisto =
171  bookDQMHistogram1D(detid,
172  rootName.str(),
173  humanName.str(),
174  numberOfVCalPoints,
175  vcalValuesToPassToCrappyRoot); // ROOT only takes an input as array. :(
176  // HOORAY FOR CINT!
178  for (unsigned int iBin = 0; iBin < numberOfVCalPoints; ++iBin) {
179  int rootBin = iBin + 1; // root bins start at 1
180  aBadHisto->setBinContent(rootBin, efficiencies[iBin]);
181  aBadHisto->setBinError(rootBin, errors[iBin]);
182  }
183 }
184 
186  edm::LogInfo("SiPixelSCurveCalibrationAnalysis")
187  << "Calibration Settings: VCalLow: " << vCalValues_[0] << " VCalHigh: " << vCalValues_[vCalValues_.size() - 1]
188  << " nVCal: " << vCalValues_.size() << " nTriggers: " << nTriggers_;
191  // build the vCal values as a vector of floats if we want to save single
192  // curves
193  const std::vector<short> *theVCalValues = this->getVcalValues();
194  unsigned int numberOfVCalPoints = theVCalValues->size();
195  edm::LogWarning("SiPixelSCurveCalibrationAnalysis")
196  << "WARNING: Option set to save indiviual S-Curves - max number: " << maxCurvesToSave_
197  << " This can lead to large memory consumption! (Got " << numberOfVCalPoints << " VCal Points";
198  for (unsigned int i = 0; i < numberOfVCalPoints; i++) {
199  vCalPointsAsFloats_.push_back(static_cast<float>((*theVCalValues)[i]));
200  edm::LogInfo("SiPixelSCurveCalibrationAnalysis") << "Adding calibration Vcal: " << (*theVCalValues)[i];
201  }
202  // must add lower edge of last bin to the vector
203  vCalPointsAsFloats_.push_back(vCalPointsAsFloats_[numberOfVCalPoints - 1] + 1);
204  }
205 
206  fitFunction_ = new TF1("sCurve",
207  "0.5*[2]*(1+TMath::Erf( (x-[0]) / ([1]*sqrt(2)) ) )",
208  vCalValues_[0],
209  vCalValues_[vCalValues_.size() - 1]);
210 }
211 
213  if (calibrationMode_ == "SCurve")
214  return true;
215  else if (calibrationMode_ == "unknown") {
216  edm::LogInfo("SiPixelSCurveCalibrationAnalysis")
217  << "calibration mode is: " << calibrationMode_ << ", continuing anyway...";
218  return true;
219  } else {
220  // edm::LogDebug("SiPixelSCurveCalibrationAnalysis") << "unknown
221  // calibration mode for SCurves, should be \"SCurve\" and is \"" <<
222  // calibrationMode_ << "\"";
223  }
224  return false;
225 }
226 
228  float &threshold,
229  float &sigma) {
231  bool allZeroSoFar = true;
232  int turnOnBin = -1;
233  int saturationBin = -1;
234  for (uint32_t iVcalPt = 0; iVcalPt < eff.size(); iVcalPt++) {
235  if (allZeroSoFar && eff[iVcalPt] != 0) {
236  turnOnBin = iVcalPt;
237  allZeroSoFar = false;
239  } else if (eff[iVcalPt] > 0.90) {
240  saturationBin = iVcalPt;
241  short turnOnVcal = vCalValues_[turnOnBin];
242  short saturationVcal = vCalValues_[saturationBin];
243  short delta = saturationVcal - turnOnVcal;
244  sigma = delta * 0.682;
245  if (sigma < 1) // check to make sure sigma guess is larger than our X resolution.
246  // Hopefully prevents Minuit from getting stuck at boundary
247  sigma = 1;
248  threshold = turnOnVcal + (0.5 * delta);
249  return errOK;
250  }
251  }
252  return output;
253 }
254 
256  float sigma,
257  float amplitude) {
258  // check if nonsensical
259  if (threshold > vCalValues_[vCalValues_.size() - 1] || threshold < vCalValues_[0] ||
260  sigma > vCalValues_[vCalValues_.size() - 1] - vCalValues_[0])
261  return errFitNonPhysical;
262 
263  if (threshold < minimumThreshold_ || threshold > maximumThreshold_ || sigma < minimumSigma_ ||
264  sigma > maximumSigma_ || amplitude < minimumEffAsymptote_ || amplitude > maximumEffAsymptote_)
265  return errFlaggedBadByUser;
266 
267  return errOK;
268 }
269 
270 void calculateEffAndError(int nADCResponse, int nTriggers, float &eff, float &error) {
271  eff = (float)nADCResponse / (float)nTriggers;
272  double effForErrorCalculation = eff;
273  if (eff <= 0 || eff >= 1)
274  effForErrorCalculation = 0.5 / (double)nTriggers;
275  error = TMath::Sqrt(effForErrorCalculation * (1 - effForErrorCalculation) / (double)nTriggers);
276 }
277 
278 // book histograms when new DetID is encountered in Event Record
280  edm::LogInfo("SiPixelSCurveCalibrationAnalysis")
281  << "Found a new DetID (" << detid << ")! Checking to make sure it has not been added.";
282  // ensure that this DetID has not been added yet
283  sCurveHistogramHolder tempMap;
284  std::pair<detIDHistogramMap::iterator, bool> insertResult;
285  insertResult = histograms_.insert(std::make_pair(detid, tempMap));
286  if (insertResult.second) // indicates successful insertion
287  {
288  edm::LogInfo("SiPixelSCurveCalibrationAnalysisHistogramReport")
289  << "Histogram Map.insert() returned true! Booking new histogrames for "
290  "detID: "
291  << detid;
292  // use detector hierarchy folders if desired
294  setDQMDirectory(detid);
295 
296  std::string detIdName = translateDetIdToString(detid);
297  if (write2dHistograms_) {
298  MonitorElement *D2sigma = bookDQMHistoPlaquetteSummary2D(detid, "ScurveSigmas", detIdName + " Sigmas");
299  MonitorElement *D2thresh = bookDQMHistoPlaquetteSummary2D(detid, "ScurveThresholds", detIdName + " Thresholds");
300  MonitorElement *D2chi2 = bookDQMHistoPlaquetteSummary2D(detid, "ScurveChi2Prob", detIdName + " Chi2Prob");
301  insertResult.first->second.insert(std::make_pair(kSigmas, D2sigma));
302  insertResult.first->second.insert(std::make_pair(kThresholds, D2thresh));
303  insertResult.first->second.insert(std::make_pair(kChi2s, D2chi2));
304  }
305  if (write2dFitResult_) {
306  MonitorElement *D2FitResult = bookDQMHistoPlaquetteSummary2D(detid, "ScurveFitResult", detIdName + " Fit Result");
307  insertResult.first->second.insert(std::make_pair(kFitResults, D2FitResult));
308  }
309  MonitorElement *D1sigma =
310  bookDQMHistogram1D(detid, "ScurveSigmasSummary", detIdName + " Sigmas Summary", 100, 0, maximumSigmaBin_);
312  detid, "ScurveThresholdSummary", detIdName + " Thresholds Summary", 255, 0, maximumThresholdBin_);
313  MonitorElement *D1chi2 =
314  bookDQMHistogram1D(detid, "ScurveChi2ProbSummary", detIdName + " Chi2Prob Summary", 101, 0, 1.01);
315  MonitorElement *D1FitResult =
316  bookDQMHistogram1D(detid, "ScurveFitResultSummary", detIdName + " Fit Result Summary", 10, -0.5, 9.5);
317  insertResult.first->second.insert(std::make_pair(kSigmaSummary, D1sigma));
318  insertResult.first->second.insert(std::make_pair(kThresholdSummary, D1thresh));
319  insertResult.first->second.insert(std::make_pair(kChi2Summary, D1chi2));
320  insertResult.first->second.insert(std::make_pair(kFitResultSummary, D1FitResult));
321  }
322 }
323 
324 bool SiPixelSCurveCalibrationAnalysis::doFits(uint32_t detid, std::vector<SiPixelCalibDigi>::const_iterator calibDigi) {
325  sCurveErrorFlag errorFlag = errOK;
326  uint32_t nVCalPts = calibDigi->getnpoints();
327  // reset and fill static datamembers with vector of points and errors
328  efficiencies_.resize(0);
329  effErrors_.resize(0);
330  for (uint32_t iVcalPt = 0; iVcalPt < nVCalPts; iVcalPt++) {
331  float eff;
332  float error;
333  calculateEffAndError(calibDigi->getnentries(iVcalPt), nTriggers_, eff, error);
334  edm::LogInfo("SiPixelSCurveCalibrationAnalysis")
335  << "Eff: " << eff << " Error: " << error << " nEntries: " << calibDigi->getnentries(iVcalPt)
336  << " nTriggers: " << nTriggers_ << " VCalPt " << vCalValues_[iVcalPt];
337  efficiencies_.push_back(eff);
338  effErrors_.push_back(error);
339  }
340 
341  // estimate the S-Curve parameters
342  float thresholdGuess = -1.0;
343  float sigmaGuess = -1.0;
344  errorFlag = estimateSCurveParameters(efficiencies_, thresholdGuess, sigmaGuess);
345 
346  // these -1.0 default values will only be filled if the curve is all zeroes,
347  // or doesn't turn on, WHICH INDICATES A SERIOUS PROBLEM
348  Double_t sigma = -1.0;
349  Double_t sigmaError = -1.0;
350  Double_t threshold = -1.0;
351  Double_t thresholdError = -1.0;
352  Double_t amplitude = -1.0;
353  Double_t amplitudeError = -1.0;
354  Double_t chi2 = -1.0;
355  // calculate NDF
356  Int_t nDOF = vCalValues_.size() - 3;
357  Double_t chi2probability = 0;
358 
359  if (errorFlag == errOK) // only do fit if curve is fittable
360  {
361  // set up minuit fit
362  TMinuit *gMinuit = new TMinuit(3);
363  gMinuit->SetPrintLevel(-1); // save ourselves from gigabytes of stdout
364  gMinuit->SetFCN(chi2toMinimize);
365 
366  // define threshold parameters - choose step size 1, max 300, min -50
367  gMinuit->DefineParameter(0, "Threshold", (Double_t)thresholdGuess, 1, -50, 300);
368  // sigma
369  gMinuit->DefineParameter(1, "Sigma", (Double_t)sigmaGuess, 0.1, 0, 255);
370  // amplitude
371  gMinuit->DefineParameter(2, "Amplitude", 1, 0.1, -0.001, 200);
372 
373  // Do Chi2 minimazation
374  gMinuit->Migrad();
375  gMinuit->GetParameter(0, threshold, thresholdError);
376  gMinuit->GetParameter(1, sigma, sigmaError);
377  gMinuit->GetParameter(2, amplitude, amplitudeError);
378 
379  // get Chi2
380  Double_t params[3] = {threshold, sigma, amplitude};
381  gMinuit->Eval(3, nullptr, chi2, params, 0);
382  // calculate Chi2 proability
383  if (nDOF <= 0)
384  chi2probability = 0;
385  else
386  chi2probability = TMath::Prob(chi2, nDOF);
387 
388  // check to make sure output makes sense (i.e. threshold > 0)
389  if (chi2probability > minimumChi2prob_)
390  errorFlag = fittedSCurveSanityCheck(threshold, sigma, amplitude);
391  else
392  errorFlag = errBadChi2Prob;
393 
394  edm::LogInfo("SiPixelSCurveCalibrationAnalysis")
395  << "Fit finished with errorFlag: " << errorFlag << " - threshold: " << threshold << " sigma: " << sigma
396  << " chi2: " << chi2 << " nDOF: " << nDOF << " chi2Prob: " << chi2probability
397  << " chi2MinUser: " << minimumChi2prob_;
398 
399  delete gMinuit;
400  }
401  // get row and column for this pixel
402  uint32_t row = calibDigi->row();
403  uint32_t col = calibDigi->col();
404 
405  // get iterator to histogram holder for this detid
406  detIDHistogramMap::iterator thisDetIdHistoGrams;
407  thisDetIdHistoGrams = histograms_.find(detid);
408  if (thisDetIdHistoGrams != histograms_.end()) {
409  edm::LogInfo("SiPixelSCurveCalibrationAnalysisHistogramReport")
410  << "Filling histograms for [detid](col/row): [" << detid << "](" << col << "/" << row
411  << ") ErrorFlag: " << errorFlag;
412  // always fill fit result
413  (*thisDetIdHistoGrams).second[kFitResultSummary]->Fill(errorFlag);
414  if (write2dFitResult_)
415  (*thisDetIdHistoGrams)
416  .second[kFitResults]
417  ->setBinContent(col + 1,
418  row + 1,
419  errorFlag); // +1 because root bins start at 1
420 
421  // fill sigma/threshold result
422  (*thisDetIdHistoGrams).second[kSigmaSummary]->Fill(sigma);
423  (*thisDetIdHistoGrams).second[kThresholdSummary]->Fill(threshold);
424  if (write2dHistograms_) {
425  (*thisDetIdHistoGrams)
426  .second[kSigmas]
427  ->setBinContent(col + 1,
428  row + 1,
429  sigma); // +1 because root bins start at 1
430  (*thisDetIdHistoGrams)
431  .second[kThresholds]
432  ->setBinContent(col + 1,
433  row + 1,
434  threshold); // +1 because root bins start at 1
435  }
436  // fill chi2
437  (*thisDetIdHistoGrams).second[kChi2Summary]->Fill(chi2probability);
438  if (write2dHistograms_)
439  (*thisDetIdHistoGrams).second[kChi2s]->Fill(col, row, chi2probability);
440  }
441  // save individual curves, if requested
443  bool thisDetIDinList = false;
444  if (detIDsToSave_.find(detid) != detIDsToSave_.end()) // see if we want to save this histogram
445  thisDetIDinList = true;
446 
447  if (errorFlag != errOK || thisDetIDinList) {
448  edm::LogError("SiPixelSCurveCalibrationAnalysis") << "Saving error histogram for [detid](col/row): [" << detid
449  << "](" << col << "/" << row << ") ErrorFlag: " << errorFlag;
450  buildACurveHistogram(detid, row, col, errorFlag, efficiencies_, effErrors_);
451  }
452  }
453 
454  return true;
455 }
std::map< sCurveHistogramType, MonitorElement * > sCurveHistogramHolder
void calculateEffAndError(int nADCResponse, int nTriggers, float &eff, float &error)
int toCabling(sipixelobjects::ElectronicIndex &cabling, const sipixelobjects::DetectorIndex &detector) const
void calibrationSetup(const edm::EventSetup &iSetup) override
std::string translateDetIdToString(uint32_t detid)
Log< level::Error, false > LogError
identify pixel inside single ROC
Definition: LocalPixel.h:7
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * bookDQMHistogram1D(uint32_t detid, std::string name, std::string title, int nchX, double lowX, double highX)
bool doFits(uint32_t detid, std::vector< SiPixelCalibDigi >::const_iterator ipix) override
T const * product() const
Definition: ESHandle.h:86
const sipixelobjects::PixelROC * findItem(const sipixelobjects::CablingPathToDetUnit &path) const final
sCurveErrorFlag estimateSCurveParameters(const std::vector< float > &eff, float &threshold, float &sigma)
virtual int getNbinsY() const
get # of bins in Y-axis
double collumn and pixel ID in double collumn representation
Definition: LocalPixel.h:19
virtual std::string getTitle() const
get MonitorElement title
MonitorElement * bookDQMHistoPlaquetteSummary2D(uint32_t detid, std::string name, std::string title)
Log< level::Info, false > LogInfo
virtual void setBinContent(int binx, double content)
set content of bin (1-D)
sCurveErrorFlag fittedSCurveSanityCheck(float threshold, float sigma, float amplitude)
void chi2toMinimize(int &npar, double *grad, double &fcnval, double *xval, int iflag)
col
Definition: cuy.py:1009
static const std::vector< short > * getVcalValues()
virtual int getNbinsX() const
get # of bins in X-axis
Definition: errors.py:1
Definition: output.py:1
virtual void setBinError(int binx, double error)
set uncertainty on content of bin (1-D)
Log< level::Warning, false > LogWarning
edm::ESHandle< SiPixelFedCablingMap > theCablingMap_
unsigned int idInDetUnit() const
id of this ROC in DetUnit etermined by token path
Definition: PixelROC.h:37
void buildACurveHistogram(const uint32_t &detid, const uint32_t &row, const uint32_t &col, sCurveErrorFlag errorFlag, const std::vector< float > &efficiencies, const std::vector< float > &errors)
virtual double getBinContent(int binx) const
get content of bin (1-D)