CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    SiPixelGainCalibrationAnalysis
00004 // Class:      SiPixelGainCalibrationAnalysis
00005 // 
00013 //
00014 // Original Author:  Freya Blekman
00015 //         Created:  Wed Nov 14 15:02:06 CET 2007
00016 // $Id: SiPixelGainCalibrationAnalysis.cc,v 1.34 2008/09/15 15:32:39 fblekman Exp $
00017 //
00018 //
00019 
00020 // user include files
00021 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00022 
00023 #include "SiPixelGainCalibrationAnalysis.h"
00024 #include <sstream>
00025 #include <vector>
00026 #include <math.h>
00027 #include "TGraphErrors.h"
00028 #include "TMath.h"
00029 //
00030 // constructors and destructor
00031 //
00032 SiPixelGainCalibrationAnalysis::SiPixelGainCalibrationAnalysis(const edm::ParameterSet& iConfig):
00033   SiPixelOfflineCalibAnalysisBase(iConfig),
00034   conf_(iConfig),
00035   bookkeeper_(),
00036   bookkeeper_pixels_(),
00037   nfitparameters_(iConfig.getUntrackedParameter<int>("numberOfFitParameters",2)),
00038   fitfunction_(iConfig.getUntrackedParameter<std::string>("fitFunctionRootFormula","pol1")),
00039   reject_plateaupoints_(iConfig.getUntrackedParameter<bool>("suppressPlateauInFit",true)),
00040   reject_single_entries_(iConfig.getUntrackedParameter<bool>("suppressPointsWithOneEntryOrLess",true)),
00041   plateau_max_slope_(iConfig.getUntrackedParameter<double>("plateauSlopeMax",1.0)),
00042   reject_first_point_(iConfig.getUntrackedParameter<bool>("rejectVCalZero",true)),
00043   reject_badpoints_frac_(iConfig.getUntrackedParameter<double>("suppressZeroAndPlateausInFitFrac",0)),
00044   bookBIGCalibPayload_(iConfig.getUntrackedParameter<bool>("saveFullPayloads",false)),
00045   savePixelHists_(iConfig.getUntrackedParameter<bool>("savePixelLevelHists",false)),
00046   chi2Threshold_(iConfig.getUntrackedParameter<double>("minChi2NDFforHistSave",10)),
00047   chi2ProbThreshold_(iConfig.getUntrackedParameter<double>("minChi2ProbforHistSave",0.05)),
00048   maxGainInHist_(iConfig.getUntrackedParameter<double>("maxGainInHist",10)),
00049   maxChi2InHist_(iConfig.getUntrackedParameter<double>("maxChi2InHist",25)),
00050   saveALLHistograms_(iConfig.getUntrackedParameter<bool>("saveAllHistograms",false)),
00051   filldb_(iConfig.getUntrackedParameter<bool>("writeDatabase",false)),
00052   recordName_(conf_.getParameter<std::string>("record")),
00053   appendMode_(conf_.getUntrackedParameter<bool>("appendMode",true)),
00054   listofdetids_(conf_.getUntrackedParameter<std::vector<uint32_t> >("listOfDetIDs")),
00055   theGainCalibrationDbInput_(0),
00056   theGainCalibrationDbInputOffline_(0),
00057   theGainCalibrationDbInputHLT_(0),
00058   theGainCalibrationDbInputService_(iConfig),
00059   gainlow_(10.),gainhi_(0.),pedlow_(255.),pedhi_(0.),
00060   useVcalHigh_(conf_.getParameter<bool>("useVCALHIGH")),
00061   scalarVcalHigh_VcalLow_(conf_.getParameter<double>("vcalHighToLowConversionFac"))
00062 {
00063   if(reject_single_entries_)
00064     min_nentries_=1;
00065   else
00066     min_nentries_=0;
00067   ::putenv("CORAL_AUTH_USER=me");
00068   ::putenv("CORAL_AUTH_PASSWORD=test");   
00069   edm::LogInfo("SiPixelGainCalibrationAnalysis") << "now using fit function " << fitfunction_ << ", which has " << nfitparameters_ << " free parameters. " << std::endl;
00070   func_= new TF1("func",fitfunction_.c_str(),0,256*scalarVcalHigh_VcalLow_);
00071   graph_ = new TGraphErrors();
00072 }
00073 
00074 SiPixelGainCalibrationAnalysis::~SiPixelGainCalibrationAnalysis()
00075 {
00076 }
00077 // member functions
00078 //
00079 // ------------ method called once each job just before starting event loop  ------------
00080 
00081 std::vector<float> SiPixelGainCalibrationAnalysis::CalculateAveragePerColumn(uint32_t detid, std::string label){
00082   std::vector<float> result;
00083   int ncols= bookkeeper_[detid][label]->getNbinsX();
00084   int nrows= bookkeeper_[detid][label]->getNbinsY();
00085   for(int icol=1; icol<=ncols; ++icol){
00086     float val=0;
00087     float ntimes =0;
00088     for(int irow=1; irow<=nrows; ++irow){
00089       val+= bookkeeper_[detid][label]->getBinContent(icol,irow);
00090       ntimes++;
00091     }
00092     val/= ntimes;
00093     result.push_back(val);
00094   }
00095   return result;
00096 }
00097 
00098 bool
00099 SiPixelGainCalibrationAnalysis::checkCorrectCalibrationType()
00100 {
00101   if(calibrationMode_=="GainCalibration")
00102     return true;
00103   else if(calibrationMode_=="unknown"){
00104     edm::LogInfo("SiPixelGainCalibrationAnalysis") <<  "calibration mode is: " << calibrationMode_ << ", continuing anyway..." ;
00105     return true;
00106   }
00107   else{
00108     //    edm::LogError("SiPixelGainCalibrationAnalysis") << "unknown calibration mode for Gain calibration, should be \"Gain\" and is \"" << calibrationMode_ << "\"";
00109   }
00110   return false;
00111 }
00112 
00113 void SiPixelGainCalibrationAnalysis::calibrationSetup(const edm::EventSetup&)
00114 {
00115 }
00116 //------- summary printing method. Very verbose.
00117 void
00118 SiPixelGainCalibrationAnalysis::printSummary(){
00119 
00120   uint32_t detid=0;
00121   for(std::map<uint32_t,std::map<std::string,MonitorElement *> >::const_iterator idet = bookkeeper_.begin(); idet != bookkeeper_.end(); ++idet){
00122     if(detid==idet->first)
00123       continue;// only do things once per detid
00124     detid=idet->first;
00125     std::vector<float> gainvec=CalculateAveragePerColumn(detid,"gain_2d");
00126     std::vector<float> pedvec =CalculateAveragePerColumn(detid,"ped_2d");
00127     std::vector<float> chi2vec = CalculateAveragePerColumn(detid,"chi2_2d");
00128     std::ostringstream summarytext;
00129 
00130     summarytext << "Summary for det ID " << detid << "(" << translateDetIdToString(detid) << ")\n";
00131     summarytext << "\t Following: values per column: column #, gain, pedestal, chi2\n";
00132     for(uint32_t i=0; i<gainvec.size(); i++)
00133       summarytext << "\t " << i << " \t" << gainvec[i] << " \t" << pedvec[i] << " \t" << chi2vec[i] << "\n";
00134     summarytext << "\t list of pixels with high chi2 (chi2> " << chi2Threshold_ << "): \n";
00135 
00136     
00137     for(std::map<std::string, MonitorElement *>::const_iterator ipix = bookkeeper_pixels_[detid].begin(); ipix!=bookkeeper_pixels_[detid].end(); ++ipix)
00138       summarytext << "\t " << ipix->first << "\n";
00139     edm::LogInfo("SiPixelGainCalibrationAnalysis") << summarytext.str() << std::endl;
00140 
00141   }
00142 
00143 }
00144 
00145 // ------------ method called once each job just after ending the event loop  ------------
00146 
00147 void 
00148 SiPixelGainCalibrationAnalysis::calibrationEnd() {
00149 
00150   //  printSummary();
00151   
00152   // this is where we loop over all histograms and save the database objects
00153   if(filldb_)
00154     fillDatabase();
00155 }
00156 //-----------method to fill the database
00157 void SiPixelGainCalibrationAnalysis::fillDatabase(){
00158  // only create when necessary.
00159   // process the minimum and maximum gain & ped values...
00160   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;
00161 
00162 }
00163 // ------------ method called to do fits to all objects available  ------------
00164 bool
00165 SiPixelGainCalibrationAnalysis::doFits(uint32_t detid, std::vector<SiPixelCalibDigi>::const_iterator ipix)
00166 {
00167   float lowmeanval=255;
00168   float highmeanval=0;
00169   bool makehistopersistent = saveALLHistograms_;
00170   std::vector<uint32_t>::const_iterator detidfinder=find(listofdetids_.begin(),listofdetids_.end(),detid);
00171   if(detidfinder!=listofdetids_.end())
00172     makehistopersistent=true;
00173   // first, fill the input arrays to the TLinearFitter.
00174   double xvals[257];
00175   double yvals[256];
00176   double yerrvals[256];
00177   double xvalsall[257];
00178   double yvalsall[256];
00179   double yerrvalsall[256];
00180   int npoints=0;
00181   int nallpoints=0;
00182   bool use_point=true;
00183   bookkeeper_[detid]["status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,0);
00184   for(uint32_t ii=0; ii< ipix->getnpoints() && ii<200; ii++){
00185     //    std::cout << ipix->getsum(ii) << " " << ipix->getnentries(ii) << " " << ipix->getsumsquares(ii) << std::endl;
00186     nallpoints++;
00187     use_point=true;
00188     if(useVcalHigh_){
00189       xvalsall[ii]=vCalValues_[ii]*scalarVcalHigh_VcalLow_;
00190     }
00191     else
00192       xvalsall[ii]=vCalValues_[ii];
00193     yerrvalsall[ii]=yvalsall[ii]=0; 
00194   
00195     if(ipix->getnentries(ii)>min_nentries_){
00196       yvalsall[ii]=ipix->getsum(ii)/(float)ipix->getnentries(ii);
00197       yerrvalsall[ii]=ipix->getsumsquares(ii)/(float)(ipix->getnentries(ii));
00198       yerrvalsall[ii]-=pow(yvalsall[ii],2);
00199       yerrvalsall[ii]=sqrt(yerrvalsall[ii])*sqrt(ipix->getnentries(ii));
00200       
00201       if(yvalsall[ii]<lowmeanval)
00202         lowmeanval=yvalsall[ii];
00203       if(yvalsall[ii]>highmeanval)
00204         highmeanval=yvalsall[ii];
00205     }
00206   }
00207   
00208   
00209   // calculate plateau value from last 3 full entries
00210   double plateauval=0;
00211   for(int ii=nallpoints-1; ii>=0 && npoints<10; --ii){
00212     if(yvalsall[ii]>0 && yerrvalsall[ii]>0){
00213       plateauval+=yvalsall[ii];
00214       npoints++;
00215     }
00216   }
00217   plateauval/=npoints;
00218   double maxgoodvalinfit=plateauval*(1.-reject_badpoints_frac_);
00219   if(maxgoodvalinfit<1)
00220     maxgoodvalinfit=255*(1.-reject_badpoints_frac_);
00221   npoints=0;
00222   for(int ii=0; ii<nallpoints; ++ii){
00223    
00224     // now selecting the appropriate points for the fit.
00225     use_point=true;
00226     if(reject_first_point_ && xvalsall[ii]<0.1)
00227       use_point=false;
00228     if(ipix->getnentries(ii)<=min_nentries_ && reject_single_entries_)
00229       use_point=false;
00230     if(ipix->getnentries(ii)==0 && reject_badpoints_)
00231       use_point=false;
00232     if(yvalsall[ii]>maxgoodvalinfit)
00233       use_point=false;
00234     if(ii>1 && fabs(yvalsall[ii]-yvalsall[ii-1])<1. && yvalsall[ii]>0.8*maxgoodvalinfit && reject_plateaupoints_)
00235       use_point=false;
00236     
00237     if(use_point){
00238       xvals[npoints]=xvalsall[ii];
00239       yvals[npoints]=yvalsall[ii];
00240       yerrvals[npoints]=yerrvalsall[ii];
00241       npoints++;
00242     }
00243   }
00244   int status=0;
00245   float chi2,slope,intercept,prob;
00246   prob=chi2=-1;
00247   slope=intercept=0;
00248   
00249 
00250   // now check on number of points. If bad just start taking the first few:
00251 
00252   if(npoints<4){
00253     npoints=0;
00254     for(int ii=0; ii<nallpoints; ++ii){
00255       if(yvalsall[ii]>0){
00256         if(ii>0 && yvalsall[ii]-yvalsall[ii-1]<0.1)
00257           continue;
00258         xvals[npoints]=xvalsall[ii];
00259         yvals[npoints]=yvalsall[ii];
00260         yerrvals[npoints]=yerrvalsall[ii];
00261         npoints++;
00262       }
00263     }
00264   }
00265   if(npoints<2){
00266     bookkeeper_[detid]["status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,-2);
00267     return false;
00268   }
00269     
00270   //  std::cout << "starting fit!" << std::endl;
00271   graph_->Set(npoints);
00272   
00273   func_->SetParameter(0,50.);
00274   func_->SetParameter(1,0.25);
00275   for(int ipointtemp=0; ipointtemp<npoints; ++ipointtemp){
00276     graph_->SetPoint(ipointtemp,xvals[ipointtemp],yvals[ipointtemp]);
00277     graph_->SetPointError(ipointtemp,0,yerrvals[ipointtemp]);
00278   }
00279   Int_t tempresult = graph_->Fit("func","FQ0N");
00280   slope=func_->GetParameter(1);
00281   intercept = func_->GetParameter(0);
00282   chi2=func_->GetChisquare()/((float)npoints-func_->GetNpar());
00283   prob= TMath::Prob(func_->GetChisquare(),npoints-func_->GetNpar());
00284   size_t ntimes=0;
00285   while((isnan(slope) || isnan(intercept) )&& ntimes<10){
00286     ntimes++;
00287     makehistopersistent=true;
00288     //    std::cout << slope << " " << intercept << " " << prob << std::endl;
00289     edm::LogWarning("SiPixelGainCalibrationAnalysis") << "impossible to fit values, try " << ntimes << ": " ;
00290     for(int ii=0; ii<npoints; ++ii){
00291       edm::LogWarning("SiPixelGainCalibrationAnalysis")<< "vcal " << xvals[ii] << " response: " << yvals[ii] << "+/-" << yerrvals[ii] << std::endl; 
00292     } 
00293     tempresult = graph_->Fit("func","FQ0NW");
00294     slope=func_->GetParameter(1);
00295     intercept = func_->GetParameter(0);
00296     chi2=func_->GetChisquare()/((float)npoints-func_->GetNpar());
00297     prob= TMath::Prob(func_->GetChisquare(),npoints-func_->GetNpar());
00298 
00299   }
00300   if(tempresult==0)
00301     status=1;
00302   else
00303     status=-2;
00304   if(slope!=0)
00305     slope = 1./slope;
00306   if(isnan(slope) || isnan(intercept)){
00307     bookkeeper_[detid]["status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,-2);
00308     return 0;
00309 
00310   }
00311   if(chi2>chi2Threshold_ && chi2Threshold_>=0)
00312     makehistopersistent=true;
00313   if(prob<chi2ProbThreshold_)
00314     makehistopersistent=true;
00315   if(status<=0)
00316     makehistopersistent=true;
00317 
00318 
00319   if(slope<gainlow_)
00320     gainlow_=slope;
00321   if(slope>gainhi_)
00322     gainhi_=slope;
00323   if(intercept>pedhi_)
00324     pedhi_=intercept;
00325   if(intercept<pedlow_)
00326     pedlow_=intercept;
00327   bookkeeper_[detid]["gain_1d"]->Fill(slope);
00328   if(slope>maxGainInHist_){
00329     makehistopersistent=true;
00330     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;
00331   }
00332   bookkeeper_[detid]["dynamicrange_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval-lowmeanval);
00333   bookkeeper_[detid]["plateau_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,highmeanval);
00334   bookkeeper_[detid]["gain_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,slope);
00335   bookkeeper_[detid]["ped_1d"]->Fill(intercept);
00336   bookkeeper_[detid]["ped_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,intercept);
00337   bookkeeper_[detid]["chi2_1d"]->Fill(chi2);
00338   bookkeeper_[detid]["chi2_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,chi2);
00339   bookkeeper_[detid]["prob_1d"]->Fill(prob);
00340   bookkeeper_[detid]["prob_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,prob);
00341   bookkeeper_[detid]["lowpoint_1d"]->Fill(xvals[0]);
00342   bookkeeper_[detid]["highpoint_1d"]->Fill(xvals[npoints-1]);
00343   bookkeeper_[detid]["nfitpoints_1d"]->Fill(npoints);
00344   bookkeeper_[detid]["endpoint_1d"]->Fill((255 - intercept)*slope);
00345   bookkeeper_[detid]["status_2d"]->setBinContent(ipix->col()+1,ipix->row()+1,status);
00346   
00347   if(!savePixelHists_)
00348     return true;
00349   if(makehistopersistent){
00350     std::ostringstream pixelinfo;
00351     pixelinfo << "GainCurve_row_" << ipix->row() << "_col_" << ipix->col();
00352     std::string tempname=translateDetIdToString(detid);
00353     tempname+="_";
00354     tempname+=pixelinfo.str();
00355     // and book the histo
00356     // fill the last value of the vcal array...   
00357 
00358     setDQMDirectory(detid);
00359     bookkeeper_pixels_[detid][pixelinfo.str()] =  bookDQMHistogram1D(detid,pixelinfo.str(),tempname,105*nallpoints,xvalsall[0],xvalsall[nallpoints-1]*1.05);
00360     
00361     edm::LogInfo("SiPixelGainCalibrationAnalysis") << "now saving histogram for pixel " << tempname << ", gain = " << slope << ", pedestal = " << intercept << ", chi2/NDF=" << chi2 << "(prob:" << prob << "), fit status " << status;
00362     for(int ii=0; ii<nallpoints; ++ii){
00363       //      std::cout << xvalsall[ii]<<","<<yvalsall[ii]<< " " << tempfloats[ii] << std::endl;
00364       bookkeeper_pixels_[detid][pixelinfo.str()]->Fill(xvalsall[ii],yvalsall[ii]);
00365     }
00366     
00367     //    addTF1ToDQMMonitoringElement(bookkeeper_pixels_[detid][pixelinfo.str()],func_);
00368   } 
00369   return true;
00370 }
00371 // ------------ method called to do fill new detids  ------------
00372 void 
00373 SiPixelGainCalibrationAnalysis::newDetID(uint32_t detid)
00374 {
00375   setDQMDirectory(detid);
00376   std::string tempname=translateDetIdToString(detid);
00377   bookkeeper_[detid]["gain_1d"] = bookDQMHistogram1D(detid,"Gain1d","gain for "+tempname,100,0.,maxGainInHist_);
00378   bookkeeper_[detid]["gain_2d"] = bookDQMHistoPlaquetteSummary2D(detid, "Gain2d","gain for "+tempname);
00379   bookkeeper_[detid]["ped_1d"] = bookDQMHistogram1D(detid,"Pedestal1d","pedestal for "+tempname,256,0.,256.0);
00380   bookkeeper_[detid]["ped_2d"] = bookDQMHistoPlaquetteSummary2D(detid,"Pedestal2d","pedestal for "+tempname);
00381   bookkeeper_[detid]["chi2_1d"] = bookDQMHistogram1D(detid,"GainChi2NDF1d","#chi^{2}/NDOF for "+tempname,100,0.,maxChi2InHist_);
00382   bookkeeper_[detid]["chi2_2d"] = bookDQMHistoPlaquetteSummary2D(detid,"GainChi2NDF2d","#chi^{2}/NDOF for "+tempname);
00383   bookkeeper_[detid]["prob_1d"] = bookDQMHistogram1D(detid,"GainChi2Prob1d","P(#chi^{2},NDOF) for "+tempname,100,0.,1.0);
00384   bookkeeper_[detid]["prob_2d"] = bookDQMHistoPlaquetteSummary2D(detid,"GainChi2Prob2d","P(#chi^{2},NDOF) for "+tempname);
00385   bookkeeper_[detid]["status_2d"] = bookDQMHistoPlaquetteSummary2D(detid,"GainFitResult2d","Fit result for "+tempname);
00386   bookkeeper_[detid]["endpoint_1d"]= bookDQMHistogram1D(detid,"GainEndPoint1d","point where fit meets ADC=255 for "+tempname,256,0.,256.*scalarVcalHigh_VcalLow_);
00387   bookkeeper_[detid]["lowpoint_1d"]= bookDQMHistogram1D(detid,"GainLowPoint1d","lowest fit point for "+tempname,256,0.,256.*scalarVcalHigh_VcalLow_);
00388   bookkeeper_[detid]["highpoint_1d"]= bookDQMHistogram1D(detid,"GainHighPoint1d","highest fit point for "+tempname,256,0.,256.*scalarVcalHigh_VcalLow_);
00389   bookkeeper_[detid]["nfitpoints_1d"]= bookDQMHistogram1D(detid,"GainNPoints1d","number of fit point for "+tempname,20,0.,20);
00390   bookkeeper_[detid]["dynamicrange_2d"]=bookDQMHistoPlaquetteSummary2D(detid,"GainDynamicRange2d","Difference lowest and highest points on gain curve for "+tempname);
00391   bookkeeper_[detid]["plateau_2d"]=bookDQMHistoPlaquetteSummary2D(detid,"GainSaturate2d","Highest points on gain curve for "+tempname);
00392 
00393 }
00394 //define this as a plug-in
00395 DEFINE_FWK_MODULE(SiPixelGainCalibrationAnalysis);