CMS 3D CMS Logo

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