27 #include "TGraphErrors.h"
40 nfitparameters_(iConfig.getUntrackedParameter<int>(
"numberOfFitParameters",2)),
41 fitfunction_(iConfig.getUntrackedParameter<std::
string>(
"fitFunctionRootFormula",
"pol1")),
42 listofdetids_(conf_.getUntrackedParameter<std::vector<uint32_t> >(
"listOfDetIDs")),
43 ignoreMode_(conf_.getUntrackedParameter<bool>(
"ignoreMode",
false)),
44 reject_plateaupoints_(iConfig.getUntrackedParameter<bool>(
"suppressPlateauInFit",
true)),
45 reject_single_entries_(iConfig.getUntrackedParameter<bool>(
"suppressPointsWithOneEntryOrLess",
true)),
46 plateau_max_slope_(iConfig.getUntrackedParameter<double>(
"plateauSlopeMax",1.0)),
47 reject_first_point_(iConfig.getUntrackedParameter<bool>(
"rejectVCalZero",
true)),
48 reject_badpoints_frac_(iConfig.getUntrackedParameter<double>(
"suppressZeroAndPlateausInFitFrac",0)),
49 bookBIGCalibPayload_(iConfig.getUntrackedParameter<bool>(
"saveFullPayloads",
false)),
50 savePixelHists_(iConfig.getUntrackedParameter<bool>(
"savePixelLevelHists",
false)),
51 chi2Threshold_(iConfig.getUntrackedParameter<double>(
"minChi2NDFforHistSave",10)),
52 chi2ProbThreshold_(iConfig.getUntrackedParameter<double>(
"minChi2ProbforHistSave",0.05)),
53 maxGainInHist_(iConfig.getUntrackedParameter<double>(
"maxGainInHist",10)),
54 maxChi2InHist_(iConfig.getUntrackedParameter<double>(
"maxChi2InHist",25)),
55 saveALLHistograms_(iConfig.getUntrackedParameter<bool>(
"saveAllHistograms",
false)),
58 filldb_(iConfig.getUntrackedParameter<bool>(
"writeDatabase",
false)),
59 writeSummary_(iConfig.getUntrackedParameter<bool>(
"writeSummary",
true)),
60 recordName_(conf_.getParameter<std::
string>(
"record")),
62 appendMode_(conf_.getUntrackedParameter<bool>(
"appendMode",
true)),
67 gainlow_(10.),gainhi_(0.),pedlow_(255.),pedhi_(0.),
68 useVcalHigh_(conf_.getParameter<bool>(
"useVCALHIGH")),
69 scalarVcalHigh_VcalLow_(conf_.getParameter<double>(
"vcalHighToLowConversionFac"))
75 ::putenv((
char*)
"CORAL_AUTH_USER=me");
76 ::putenv((
char*)
"CORAL_AUTH_PASSWORD=test");
79 graph_ =
new TGraphErrors();
81 summary_.open(
"SummaryPerDetID.txt");
98 for(
int icol=1; icol<=ncols; ++icol){
101 for(
int irow=1; irow<=nrows; ++irow){
106 result.push_back(val);
137 if(detid==idet->first)
143 std::ostringstream summarytext;
146 summarytext <<
"\t Following: values per column: column #, gain, pedestal, chi2\n";
147 for(uint32_t
i=0;
i<gainvec.size();
i++)
148 summarytext <<
"\t " <<
i <<
" \t" << gainvec[
i] <<
" \t" << pedvec[
i] <<
" \t" << chi2vec[
i] <<
"\n";
149 summarytext <<
"\t list of pixels with high chi2 (chi2> " <<
chi2Threshold_ <<
"): \n";
153 summarytext <<
"\t " << ipix->first <<
"\n";
154 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") << summarytext.str() << std::endl;
161 summary_<<
"Number of pixel tagged with status :"<<endl;
186 edm::LogError(
"SiPixelGainCalibration::fillDatabase()") <<
"PLEASE do not fill the database directly from the gain calibration analyzer. This function is currently disabled and no DB payloads will be produced!" << std::endl;
193 float lowmeanval=255;
198 makehistopersistent=
true;
202 double yerrvals[256];
203 double xvalsall[257];
204 double yvalsall[256];
205 double yerrvalsall[256];
212 bookkeeper_[
detid][
"status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,0);
218 for(uint32_t
ii=0;
ii< ipix->getnpoints() &&
ii<200;
ii++){
227 yerrvalsall[
ii]=yvalsall[
ii]=0;
230 yvalsall[
ii]=ipix->getsum(
ii)/(float)ipix->getnentries(
ii);
231 yerrvalsall[
ii]=ipix->getsumsquares(
ii)/(float)(ipix->getnentries(
ii));
232 yerrvalsall[
ii]-=
pow(yvalsall[
ii],2);
233 yerrvalsall[
ii]=
sqrt(yerrvalsall[ii])/
sqrt(ipix->getnentries(ii));
235 if(yvalsall[ii]<lowmeanval)
236 lowmeanval=yvalsall[
ii];
237 if(yvalsall[ii]>highmeanval)
238 highmeanval=yvalsall[
ii];
246 for(
int ii=nallpoints-1;
ii>nallpoints-5; --
ii) plateauval+=yvalsall[
ii];
248 for(
int ii=nallpoints-1;
ii>nallpoints-5; --
ii){
249 if(fabs(yvalsall[
ii]-plateauval)>5){
256 int NbofPointsInPlateau=0;
257 for(
int ii=0;
ii<nallpoints; ++
ii)
258 if(fabs(yvalsall[
ii]-plateauval)<10 || yvalsall[
ii]==0) NbofPointsInPlateau++;
260 if(NbofPointsInPlateau>=(nallpoints-2)){
264 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
274 for(
int ii=0;
ii<nallpoints; ++
ii){
284 if(yvalsall[
ii]>maxgoodvalinfit && !noPlateau)
292 xvals[npoints]=xvalsall[
ii];
293 yvals[npoints]=yvalsall[
ii];
294 yerrvals[npoints]=yerrvalsall[
ii];
299 float chi2,
slope,intercept,
prob,slopeerror,intercepterror;
301 slope=intercept=slopeerror=intercepterror=0;
308 for(
int ii=0;
ii<nallpoints && npoints<4 && yvalsall[
ii]<plateauval*0.97; ++
ii){
310 if(
ii>0 && yvalsall[
ii]-yvalsall[
ii-1]<0.1)
312 xvals[npoints]=xvalsall[
ii];
313 yvals[npoints]=yvalsall[
ii];
314 yerrvals[npoints]=yerrvalsall[
ii];
323 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
326 std::ostringstream pixelinfo;
327 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
330 tempname+=pixelinfo.str();
333 for(
int ii=0;
ii<nallpoints; ++
ii)
341 func_->SetParameter(0,50.);
342 func_->SetParameter(1,0.25);
343 for(
int ipointtemp=0; ipointtemp<npoints; ++ipointtemp){
344 graph_->SetPoint(ipointtemp,xvals[ipointtemp],yvals[ipointtemp]);
345 graph_->SetPointError(ipointtemp,0,yerrvals[ipointtemp]);
347 Int_t tempresult =
graph_->Fit(
"func",
"FQ0N");
348 slope=
func_->GetParameter(1);
349 slopeerror=
func_->GetParError(1);
350 intercept=
func_->GetParameter(0);
351 intercepterror=
func_->GetParError(0);
352 chi2=
func_->GetChisquare()/((float)npoints-
func_->GetNpar());
353 prob= TMath::Prob(
func_->GetChisquare(),npoints-
func_->GetNpar());
355 while((
isnan(slope) ||
isnan(intercept) )&& ntimes<10){
357 makehistopersistent=
true;
359 edm::LogWarning(
"SiPixelGainCalibrationAnalysis") <<
"impossible to fit values, try " << ntimes <<
": " ;
360 for(
int ii=0;
ii<npoints; ++
ii){
361 edm::LogWarning(
"SiPixelGainCalibrationAnalysis")<<
"vcal " << xvals[
ii] <<
" response: " << yvals[
ii] <<
"+/-" << yerrvals[
ii] << std::endl;
363 tempresult =
graph_->Fit(
"func",
"FQ0NW");
364 slope=
func_->GetParameter(1);
365 slopeerror=
func_->GetParError(1);
366 intercept =
func_->GetParameter(0);
367 intercepterror=
func_->GetParError(0);
368 chi2=
func_->GetChisquare()/((float)npoints-
func_->GetNpar());
369 prob= TMath::Prob(
func_->GetChisquare(),npoints-
func_->GetNpar());
382 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
398 makehistopersistent=
true;
412 makehistopersistent=
true;
413 edm::LogWarning(
"SiPixelGainCalibration") <<
"For DETID " << detid <<
"pixel row,col " << ipix->row() <<
"," << ipix->col() <<
" Gain was measured to be " << slope <<
" which is outside the range of the summary plot (" <<
maxGainInHist_ <<
") !!!! " << std::endl;
415 bookkeeper_[
detid][
"dynamicrange_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval-lowmeanval);
416 bookkeeper_[
detid][
"plateau_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval);
418 bookkeeper_[
detid][
"errorgain_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,slopeerror);
420 bookkeeper_[
detid][
"ped_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,intercept);
421 bookkeeper_[
detid][
"errorped_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,intercepterror);
423 bookkeeper_[
detid][
"chi2_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,chi2);
427 bookkeeper_[
detid][
"lowpoint_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,xvals[0]);
429 bookkeeper_[
detid][
"highpoint_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,xvals[npoints-1]);
438 if(makehistopersistent){
439 std::ostringstream pixelinfo;
440 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
443 tempname+=pixelinfo.str();
451 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") <<
"now saving histogram for pixel " << tempname <<
", gain = " << slope <<
", pedestal = " << intercept <<
", chi2/NDF=" << chi2 <<
"(prob:" << prob <<
"), fit status " <<
status;
452 for(
int ii=0;
ii<nallpoints; ++
ii){
461 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col();
~SiPixelGainCalibrationAnalysis()
bool reject_single_entries_
virtual bool checkCorrectCalibrationType()
std::map< uint32_t, std::map< std::string, MonitorElement * > > bookkeeper_pixels_
virtual void calibrationEnd()
#define DEFINE_FWK_MODULE(type)
static const double slope[3]
std::string translateDetIdToString(uint32_t detid)
virtual void calibrationSetup(const edm::EventSetup &iSetup)
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_
virtual bool doFits(uint32_t detid, std::vector< SiPixelCalibDigi >::const_iterator ipix)
bool reject_plateaupoints_
double scalarVcalHigh_VcalLow_
double chi2ProbThreshold_
MonitorElement * bookDQMHistogram1D(uint32_t detid, std::string name, std::string title, int nchX, double lowX, double highX)
virtual void newDetID(uint32_t detid)
static std::vector< short > vCalValues_
std::vector< uint32_t > listofdetids_
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)
Power< A, B >::type pow(const A &a, const B &b)