Reference Manual

Go to the documentation of this file.
00001 #include "DQMServices/Diagnostic/interface/HDQMfitUtilities.h"
00002 
00003 namespace HDQMUtil{
00004   //-----------------------------------------------------------------------------------------------
00005   double langaufun(double *x, double *par){
00006     //Fit parameters:
00007     //par[0]=Width (scale) parameter of Landau density
00008     //par[1]=Most Probable (MP, location) parameter of Landau density
00009     //par[2]=Total area (integral -inf to inf, normalization constant)
00010     //par[3]=Width (sigma) of convoluted Gaussian function
00011     //
00012     //In the Landau distribution (represented by the CERNLIB approximation), 
00013     //the maximum is located at x=-0.22278298 with the location parameter=0.
00014     //This shift is corrected within this function, so that the actual
00015     //maximum is identical to the MP parameter.
00016 
00017     // Numeric constants
00018     double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
00019     double mpshift  = -0.22278298;       // Landau maximum location
00020 
00021     // Control constants
00022     double np = 100.0;      // number of convolution steps
00023     double sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
00024 
00025     // Variables
00026     double xx;
00027     double mpc;
00028     double fland;
00029     double sum = 0.0;
00030     double xlow,xupp;
00031     double step;
00032     double i;
00033 
00034 
00035     // MP shift correction
00036     mpc = par[1] - mpshift * par[0]; 
00037 
00038     // Range of convolution integral
00039     xlow = x[0] - sc * par[3];
00040     xupp = x[0] + sc * par[3];
00041 
00042     step = (xupp-xlow) / np;
00043 
00044     // Landau Distribution Production
00045     for(i=1.0; i<=np/2; i++) {
00046       xx = xlow + (i-.5) * step;
00047       fland = TMath::Landau(xx,mpc,par[0]) / par[0];
00048       sum += fland * TMath::Gaus(x[0],xx,par[3]);
00049 
00050       xx = xupp - (i-.5) * step;
00051       fland = TMath::Landau(xx,mpc,par[0]) / par[0];
00052       sum += fland * TMath::Gaus(x[0],xx,par[3]);
00053     }
00054 
00055     return (par[2] * step * sum * invsq2pi / par[3]);
00056   }
00057 
00058   //-----------------------------------------------------------------------------------------------
00059   int32_t langaupro(double *params, double &maxx, double &FWHM) {
00060     edm::LogInfo("fitUtility") << "inside langaupro " << std::endl;
00061     // Seaches for the location (x value) at the maximum of the 
00062     // Landau and its full width at half-maximum.
00063     //
00064     // The search is probably not very efficient, but it's a first try.
00065 
00066     double p,x,fy,fxr,fxl;
00067     double step;
00068     double l,lold,dl;
00069     int32_t i = 0;
00070     const int32_t MAXCALLS = 10000;
00071     const double dlStop = 1e-3; // relative change < .001
00072 
00073     // Search for maximum
00074     p = params[1] - 0.1 * params[0];
00075     step = 0.05 * params[0];
00076     lold = -2.0;
00077     l    = -1.0;
00078 
00079     dl = (l-lold)/lold;    // FIXME catch divide by zero
00080     while ( (TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
00081       i++;
00082 
00083       lold = l;
00084       x = p + step;
00085       l = langaufun(&x,params);
00086       dl = (l-lold)/lold; // FIXME catch divide by zero
00087         
00088       if (l < lold)
00089         step = -step/10;
00090  
00091       p += step;
00092     }
00093 
00094     if (i == MAXCALLS)
00095       return (-1);
00096 
00097     maxx = x;
00098 
00099     fy = l/2;
00100 
00101 
00102     // Search for right x location of fy
00103     p = maxx + params[0];
00104     step = params[0];
00105     lold = -2.0;
00106     l    = -1e300;
00107     i    = 0;
00108 
00109     dl = (l-lold)/lold;   // FIXME catch divide by zero
00110     while ( ( TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
00111       i++;
00112 
00113       lold = l;
00114       x = p + step;
00115       l = TMath::Abs(langaufun(&x,params) - fy);
00116       dl = (l-lold)/lold; // FIXME catch divide by zero
00117  
00118       if (l > lold)
00119         step = -step/10;
00120  
00121       p += step;
00122     }
00123 
00124     if (i == MAXCALLS)
00125       return (-2);
00126 
00127     fxr = x;
00128 
00129 
00130     // Search for left x location of fy
00131     p = maxx - 0.5 * params[0];
00132     step = -params[0];
00133     lold = -2.0;
00134     l    = -1e300;
00135     i    = 0;
00136 
00137     dl = (l-lold)/lold;    // FIXME catch divide by zero
00138     while ( ( TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
00139       i++;
00140 
00141       lold = l;
00142       x = p + step;
00143       l = TMath::Abs(langaufun(&x,params) - fy);
00144       dl = (l-lold)/lold; // FIXME catch divide by zero
00145  
00146       if (l > lold)
00147         step = -step/10;
00148  
00149       p += step;
00150     }
00151 
00152     if (i == MAXCALLS)
00153       return (-3);
00154 
00155 
00156     fxl = x;
00157 
00158     FWHM = fxr - fxl;
00159     return (0);
00160   }
00161 
00162 
00163   //-----------------------------------------------------------------------------------------------
00164   double Gauss(double *x, double *par)
00165   {
00166     // The noise function: a gaussian
00167 
00168     double arg = 0;
00169     if (par[2]) arg = (x[0] - par[1])/par[2];
00170 
00171     double noise = par[0]*TMath::Exp(-0.5*arg*arg);
00172     return noise;
00173   }
00174 
00175 }
00176 
00177 //-----------------------------------------------------------------------------------------------
00178 HDQMfitUtilities::HDQMfitUtilities():langausFit(0),gausFit(0){
00179   init();
00180 }
00181 
00182 //-----------------------------------------------------------------------------------------------
00183 void HDQMfitUtilities::init(){
00184   pLanGausS[0]=0 ;   pLanGausS[1]=0;   pLanGausS[2]=0;   pLanGausS[3]=0;
00185   epLanGausS[0]=0;  epLanGausS[1]=0;  epLanGausS[2]=0;  epLanGausS[3]=0;
00186   pGausS[0]=0 ;   pGausS[1]=0;   pGausS[2]=0;
00187   epGausS[0]=0;  epGausS[1]=0;  epGausS[2]=0;
00188   pLanConv[0]=0;  pLanConv[1]=0;
00189 
00190   //if ( langausFit!=0 ) delete langausFit;
00191   //if ( gausFit!=0 ) delete gausFit;
00192 
00193   chi2GausS = -9999.;
00194   nDofGausS = -9999;
00195 }
00196 
00197 //-----------------------------------------------------------------------------------------------
00198 HDQMfitUtilities::~HDQMfitUtilities(){
00199   if ( langausFit!=0 ) delete langausFit;
00200   if ( gausFit!=0 ) delete gausFit;
00201 }
00202 
00203   
00204  
00205 //-----------------------------------------------------------------------------------------------
00206 double HDQMfitUtilities::doLanGaussFit(TH1F* htoFit){
00207   init();
00208  
00209   // if (htoFit->GetEntries()!=0) {
00210   // Check for the entries excluding over/underflows
00211   if (htoFit->Integral()!=0) {
00212      edm::LogInfo("fitUtility")<<"Fitting "<< htoFit->GetTitle() <<std::endl;
00213     // Setting fit range and start values
00214     double fr[2];
00215     double sv[4], pllo[4], plhi[4];
00216     fr[0]=0.5*htoFit->GetMean();
00217     fr[1]=3.0*htoFit->GetMean();
00218               
00219     // (EM) parameters setting good for signal only 
00220     int32_t imax=htoFit->GetMaximumBin();
00221     double xmax=htoFit->GetBinCenter(imax);
00222     double ymax=htoFit->GetBinContent(imax);
00223     int32_t i[2];
00224     int32_t iArea[2];
00225               
00226     i[0]=htoFit->GetXaxis()->FindBin(fr[0]);
00227     i[1]=htoFit->GetXaxis()->FindBin(fr[1]);
00228                   
00229     iArea[0]=htoFit->GetXaxis()->FindBin(fr[0]);
00230     iArea[1]=htoFit->GetXaxis()->FindBin(fr[1]);
00231     double AreaFWHM=htoFit->Integral(iArea[0],iArea[1],"width");
00232                   
00233     sv[1]=xmax;
00234     sv[2]=htoFit->Integral(i[0],i[1],"width");
00235     sv[3]=AreaFWHM/(4*ymax);
00236     sv[0]=sv[3];
00237                   
00238     plhi[0]=25.0; plhi[1]=200.0; plhi[2]=1000000.0; plhi[3]=50.0;
00239     pllo[0]=1.5 ; pllo[1]=10.0 ; pllo[2]=1.0      ; pllo[3]= 1.0;
00240                   
00241     Char_t FunName[100];
00242     sprintf(FunName,"FitfcnLG_%s",htoFit->GetName());  
00243                   
00244     langausFit = new TF1(FunName,HDQMUtil::langaufun,fr[0],fr[1],4);
00245     langausFit->SetParameters(sv);
00246     langausFit->SetParNames("Width","MP","Area","GSigma");
00247                   
00248     for (int32_t i=0; i<4; i++) {
00249       langausFit->SetParLimits(i,pllo[i],plhi[i]);
00250     }  
00251                   
00252     try{
00253       htoFit->Fit(langausFit,"R0");  // "R" fit in a range,"0" quiet fit
00254       
00255       langausFit->SetRange(fr[0],fr[1]);
00256       langausFit->GetParameters(pLanGausS);
00257       std::memcpy((void*) epLanGausS, (void*) langausFit->GetParErrors(), 4*sizeof(double));
00258       
00259       chi2GausS =langausFit->GetChisquare();  // obtain chi^2
00260       nDofGausS = langausFit->GetNDF();           // obtain ndf
00261       
00262       double sPeak, sFWHM;
00263       HDQMUtil::langaupro(pLanGausS,sPeak,sFWHM);
00264       pLanConv[0]=sPeak;
00265       pLanConv[1]=sFWHM;
00266       edm::LogInfo("fitUtility") << "langaupro:  max  " << sPeak << std::endl;
00267       edm::LogInfo("fitUtility") << "langaupro:  FWHM " << sFWHM << std::endl;
00268     }
00269     catch(cms::Exception& iException){
00270       edm::LogError("fitUtility") << "problem in fitting " << htoFit->GetTitle() << " \n\tDefault values of the parameters will be used";
00271       pLanGausS[0]=-9999; pLanGausS[1]=-9999; pLanGausS[2]=-9999; pLanGausS[3]=-9999;
00272       epLanGausS[0]=-9999; epLanGausS[1]=-9999; epLanGausS[2]=-9999; epLanGausS[3]=-9999;
00273       pLanConv[0]=-9999;   pLanConv[1]=-9999;   
00274       chi2GausS=-9999;  nDofGausS=-9999;    
00275     }
00276   }
00277   else {  
00278     pLanGausS[0]=-9999; pLanGausS[1]=-9999; pLanGausS[2]=-9999; pLanGausS[3]=-9999;
00279     epLanGausS[0]=-9999; epLanGausS[1]=-9999; epLanGausS[2]=-9999; epLanGausS[3]=-9999;
00280     pLanConv[0]=-9999;   pLanConv[1]=-9999;   
00281     chi2GausS=-9999;  nDofGausS=-9999;    
00282   }
00283 
00284   return htoFit->GetEntries();
00285 
00286 }
00287 
00288   
00289 //-----------------------------------------------------------------------------------------------
00290 double HDQMfitUtilities::doGaussFit(TH1F* htoFit){
00291   init();
00292   // if (htoFit->GetEntries()!=0) {
00293   if (htoFit->Integral()!=0) {
00294     
00295     // Setting fit range and start values
00296     double fr[2];
00297     double sv[3], pllo[3], plhi[3];
00298     fr[0]=htoFit->GetMean()-5*htoFit->GetRMS();
00299     fr[1]=htoFit->GetMean()+5*htoFit->GetRMS();
00300                   
00301     int32_t imax=htoFit->GetMaximumBin();
00302     double xmax=htoFit->GetBinCenter(imax);
00303     double ymax=htoFit->GetBinContent(imax);
00304     int32_t i[2];
00305     int32_t iArea[2];
00306                   
00307     i[0]=htoFit->GetXaxis()->FindBin(fr[0]);
00308     i[1]=htoFit->GetXaxis()->FindBin(fr[1]);
00309                   
00310     iArea[0]=htoFit->GetXaxis()->FindBin(fr[0]);
00311     iArea[1]=htoFit->GetXaxis()->FindBin(fr[1]);
00312     double AreaFWHM=htoFit->Integral(iArea[0],iArea[1],"width");
00313                   
00314     sv[2]=AreaFWHM/(4*ymax);
00315     sv[1]=xmax;
00316     sv[0]=htoFit->Integral(i[0],i[1],"width");
00317                   
00318     plhi[0]=1000000.0; plhi[1]=10.0; plhi[2]=10.;
00319     pllo[0]=1.5 ; pllo[1]=0.1; pllo[2]=0.3;
00320     Char_t FunName[100];
00321     sprintf(FunName,"FitfcnLG_%s",htoFit->GetName());
00322     gausFit = new TF1(FunName,HDQMUtil::Gauss,fr[0],fr[1],3);
00323     gausFit->SetParameters(sv);
00324     gausFit->SetParNames("Constant","GaussPeak","Sigma");
00325                   
00326     for (int32_t i=0; i<3; i++) {
00327       gausFit->SetParLimits(i,pllo[i],plhi[i]);
00328     }
00329  
00330     try{
00331       htoFit->Fit(gausFit,"R0");
00332       
00333       gausFit->SetRange(fr[0],fr[1]);
00334       gausFit->GetParameters(pGausS);
00335       std::memcpy((void*) epGausS, (void*) gausFit->GetParErrors(), 3*sizeof(double));
00336       
00337       chi2GausS =langausFit->GetChisquare(); // obtain chi^2
00338       nDofGausS = langausFit->GetNDF();// obtain ndf
00339     }
00340     catch(cms::Exception& iException){
00341       edm::LogError("fitUtility") << "problem in fitting " << htoFit->GetTitle() << " \n\tDefault values of the parameters will be used";
00342       pGausS[0]=-9999; pGausS[1]=-9999; pGausS[2]=-9999;
00343       epGausS[0]=-9999; epGausS[1]=-9999; epGausS[2]=-9999;
00344       chi2GausS=-9999;  nDofGausS=-9999;    
00345     }
00346                  
00347   }
00348   else {
00349     pGausS[0]=-9999; pGausS[1]=-9999; pGausS[2]=-9999;
00350     epGausS[0]=-9999; epGausS[1]=-9999; epGausS[2]=-9999;
00351     chi2GausS=-9999;  nDofGausS=-9999;    
00352   }
00353   
00354   return htoFit->GetEntries();
00355 }
00356 
00357 //-----------------------------------------------------------------------------------------------
00358 double HDQMfitUtilities::getLanGaussPar(std::string s){
00359   if(s=="landau_width")
00360     return pLanGausS[0];
00361   else if(s=="mpv")
00362     return pLanGausS[1];
00363   else if(s=="area")
00364     return pLanGausS[2];
00365   else if(s=="gauss_sigma")
00366     return pLanGausS[3];
00367   else
00368     return -99999;
00369 }
00370 
00371 double HDQMfitUtilities::getLanGaussParErr(std::string s){
00372   if(s=="landau_width")
00373     return epLanGausS[0];
00374   else if(s=="mpv")
00375     return epLanGausS[1];
00376   else if(s=="area")
00377     return epLanGausS[2];
00378   else if(s=="gauss_sigma")
00379     return epLanGausS[3];
00380   else
00381     return -99999;
00382 }
00383 
00384 double HDQMfitUtilities::getLanGaussConv(std::string s) {
00385   if(s=="mpv")
00386     return pLanConv[0];
00387   else if(s=="fwhm")
00388     return pLanConv[1];
00389   else
00390     return -99999;
00391 }
00392 
00393 double HDQMfitUtilities::getGaussPar(std::string s) {
00394   if(s=="area")
00395     return pGausS[0];
00396   else if(s=="mean")
00397     return pGausS[1];
00398   else if(s=="sigma")
00399     return pGausS[2];
00400   else
00401     return -99999;
00402 }
00403 
00404 double HDQMfitUtilities::getGaussParErr(std::string s) {
00405   if(s=="area")
00406     return epGausS[0];
00407   else if(s=="mean")
00408     return epGausS[1];
00409   else if(s=="sigma")
00410     return epGausS[2];
00411   else
00412     return -99999;
00413 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQMServices/Diagnostic/src/HDQMfitUtilities.cc
