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)),
57 filldb_(iConfig.getUntrackedParameter<bool>(
"writeDatabase",
false)),
58 writeSummary_(iConfig.getUntrackedParameter<bool>(
"writeSummary",
true)),
59 recordName_(conf_.getParameter<std::
string>(
"record")),
61 appendMode_(conf_.getUntrackedParameter<bool>(
"appendMode",
true)),
66 gainlow_(10.),gainhi_(0.),pedlow_(255.),pedhi_(0.),
67 useVcalHigh_(conf_.getParameter<bool>(
"useVCALHIGH")),
68 scalarVcalHigh_VcalLow_(conf_.getParameter<double>(
"vcalHighToLowConversionFac"))
74 ::putenv((
char*)
"CORAL_AUTH_USER=me");
75 ::putenv((
char*)
"CORAL_AUTH_PASSWORD=test");
78 graph_ =
new TGraphErrors();
80 summary_.open(
"SummaryPerDetID.txt");
97 for(
int icol=1; icol<=ncols; ++icol){
100 for(
int irow=1; irow<=nrows; ++irow){
105 result.push_back(val);
136 if(detid==idet->first)
142 std::ostringstream summarytext;
145 summarytext <<
"\t Following: values per column: column #, gain, pedestal, chi2\n";
146 for(uint32_t
i=0;
i<gainvec.size();
i++)
147 summarytext <<
"\t " <<
i <<
" \t" << gainvec[
i] <<
" \t" << pedvec[
i] <<
" \t" << chi2vec[
i] <<
"\n";
148 summarytext <<
"\t list of pixels with high chi2 (chi2> " <<
chi2Threshold_ <<
"): \n";
152 summarytext <<
"\t " << ipix->first <<
"\n";
153 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") << summarytext.str() << std::endl;
160 summary_<<
"Number of pixel tagged with status :"<<endl;
185 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;
192 float lowmeanval=255;
197 makehistopersistent=
true;
201 double yerrvals[256];
202 double xvalsall[257];
203 double yvalsall[256];
204 double yerrvalsall[256];
211 bookkeeper_[
detid][
"status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,0);
217 for(uint32_t
ii=0;
ii< ipix->getnpoints() &&
ii<200;
ii++){
226 yerrvalsall[
ii]=yvalsall[
ii]=0;
229 yvalsall[
ii]=ipix->getsum(
ii)/(float)ipix->getnentries(
ii);
230 yerrvalsall[
ii]=ipix->getsumsquares(
ii)/(float)(ipix->getnentries(
ii));
231 yerrvalsall[
ii]-=
pow(yvalsall[
ii],2);
232 yerrvalsall[
ii]=
sqrt(yerrvalsall[ii])/
sqrt(ipix->getnentries(ii));
234 if(yvalsall[ii]<lowmeanval)
235 lowmeanval=yvalsall[
ii];
236 if(yvalsall[ii]>highmeanval)
237 highmeanval=yvalsall[
ii];
245 for(
int ii=nallpoints-1;
ii>nallpoints-5; --
ii) plateauval+=yvalsall[
ii];
247 for(
int ii=nallpoints-1;
ii>nallpoints-5; --
ii){
248 if(fabs(yvalsall[
ii]-plateauval)>5){
255 int NbofPointsInPlateau=0;
256 for(
int ii=0;
ii<nallpoints; ++
ii)
257 if(fabs(yvalsall[
ii]-plateauval)<10 || yvalsall[
ii]==0) NbofPointsInPlateau++;
259 if(NbofPointsInPlateau>=(nallpoints-2)){
263 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
273 for(
int ii=0;
ii<nallpoints; ++
ii){
283 if(yvalsall[
ii]>maxgoodvalinfit && !noPlateau)
298 float chi2,
slope,intercept,prob,slopeerror,intercepterror;
300 slope=intercept=slopeerror=intercepterror=0;
307 for(
int ii=0;
ii<nallpoints && npoints<4 && yvalsall[
ii]<plateauval*0.97; ++
ii){
309 if(
ii>0 && yvalsall[
ii]-yvalsall[
ii-1]<0.1)
322 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
325 std::ostringstream pixelinfo;
326 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
329 tempname+=pixelinfo.str();
332 for(
int ii=0;
ii<nallpoints; ++
ii)
340 func_->SetParameter(0,50.);
341 func_->SetParameter(1,0.25);
342 for(
int ipointtemp=0; ipointtemp<
npoints; ++ipointtemp){
343 graph_->SetPoint(ipointtemp,xvals[ipointtemp],yvals[ipointtemp]);
344 graph_->SetPointError(ipointtemp,0,yerrvals[ipointtemp]);
347 slope=
func_->GetParameter(1);
348 slopeerror=
func_->GetParError(1);
349 intercept=
func_->GetParameter(0);
350 intercepterror=
func_->GetParError(0);
351 chi2=
func_->GetChisquare()/((float)npoints-
func_->GetNpar());
352 prob= TMath::Prob(
func_->GetChisquare(),npoints-
func_->GetNpar());
354 while((
isnan(slope) ||
isnan(intercept) )&& ntimes<10){
356 makehistopersistent=
true;
358 edm::LogWarning(
"SiPixelGainCalibrationAnalysis") <<
"impossible to fit values, try " << ntimes <<
": " ;
360 edm::LogWarning(
"SiPixelGainCalibrationAnalysis")<<
"vcal " << xvals[
ii] <<
" response: " << yvals[
ii] <<
"+/-" << yerrvals[
ii] << std::endl;
363 slope=
func_->GetParameter(1);
364 slopeerror=
func_->GetParError(1);
365 intercept =
func_->GetParameter(0);
366 intercepterror=
func_->GetParError(0);
367 chi2=
func_->GetChisquare()/((float)npoints-
func_->GetNpar());
368 prob= TMath::Prob(
func_->GetChisquare(),npoints-
func_->GetNpar());
381 summary_<<
"row_"<<ipix->row()<<
" col_"<<ipix->col()<<
" status_"<<status<<endl;
397 makehistopersistent=
true;
411 makehistopersistent=
true;
412 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;
414 bookkeeper_[
detid][
"dynamicrange_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval-lowmeanval);
415 bookkeeper_[
detid][
"plateau_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval);
417 bookkeeper_[
detid][
"errorgain_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,slopeerror);
419 bookkeeper_[
detid][
"ped_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,intercept);
420 bookkeeper_[
detid][
"errorped_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,intercepterror);
422 bookkeeper_[
detid][
"chi2_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,chi2);
424 bookkeeper_[
detid][
"prob_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,prob);
426 bookkeeper_[
detid][
"lowpoint_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,xvals[0]);
428 bookkeeper_[
detid][
"highpoint_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,xvals[npoints-1]);
437 if(makehistopersistent){
438 std::ostringstream pixelinfo;
439 pixelinfo <<
"GainCurve_row_" << ipix->row() <<
"_col_" << ipix->col();
442 tempname+=pixelinfo.str();
450 edm::LogInfo(
"SiPixelGainCalibrationAnalysis") <<
"now saving histogram for pixel " << tempname <<
", gain = " << slope <<
", pedestal = " << intercept <<
", chi2/NDF=" << chi2 <<
"(prob:" << prob <<
"), fit status " <<
status;
451 for(
int ii=0;
ii<nallpoints; ++
ii){
460 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)
Abs< T >::type abs(const T &t)
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)
volatile std::atomic< bool > shutdown_flag false
Power< A, B >::type pow(const A &a, const B &b)