CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/RecoJets/src/CompMethods.cc

Go to the documentation of this file.
00001 #include <memory>
00002 #include <string>
00003 #include <fstream>
00004 #include <iostream>
00005 #include <cmath>
00006 
00007 #include "Validation/RecoJets/interface/CompMethods.h"
00008 
00009 using std::cerr;
00010 using std::cout;
00011 using std::endl;
00012 
00013 double Gauss(double *x, double *par){
00014   double arg=0;
00015   if(par[2]!=0){
00016     arg = (x[0]-par[1])/par[2];
00017   }
00018   return par[0]*TMath::Exp(-0.5*arg*arg);
00019 }
00020 
00021 // ------------------------------------------------------------------------------------------------------------------
00022 
00023 StabilizedGauss::StabilizedGauss(const char* funcName, int funcType, double lowerBound, double upperBound):
00024   funcType_  (funcType  ),
00025   funcName_  (funcName  ),
00026   lowerBound_(lowerBound),
00027   upperBound_(upperBound)
00028 {
00029   if(funcType_==0){
00030     func_ = new TF1(funcName_, Gauss, lowerBound_, upperBound_, 3);  
00031     func_->SetParNames( "Const", "Mean", "Sigma" );
00032   }
00033   else{
00034     std::cout << "Sorry: not yet implemented" << std::endl;
00035   }
00036 }
00037 
00038 void
00039 StabilizedGauss::fit(TH1F& hist)
00040 {
00041   //set start values for first iteration
00042   if(funcType_==0){
00043     double maxValue=hist.GetBinCenter(hist.GetMaximumBin());
00044     func_->SetParameter(1, maxValue);
00045     func_->SetParameter(2, hist.GetRMS());
00046   }
00047 
00048   //set parameter limits
00049   if(funcType_==0){
00050     func_->SetParLimits(1, lowerBound_, upperBound_);
00051     func_->SetParLimits(2, 0., 5.*hist.GetRMS());
00052   }
00053 
00054   //do the fit
00055   mean_ = func_->GetParameter(1);
00056   sigma_= func_->GetParameter(2);
00057 
00058   hist.Fit( "func", "RE0", "", (mean_-2.*sigma_), (mean_+2.*sigma_) );
00059   if(hist.GetFunction("func")){
00060     //get mean and sigma 
00061     //from first iteration
00062     mean_ = hist.GetFunction("func")->GetParameter(1);
00063     sigma_= hist.GetFunction("func")->GetParameter(2);
00064 
00065     //set start values for 
00066     //second iteration
00067     func_->SetParameter(1, mean_ );
00068     func_->SetParameter(2, sigma_);
00069     hist.Fit( func_, "MEL", "", (mean_-1.5*sigma_), (mean_+1.5*sigma_) );
00070   }
00071   else{
00072     std::cout << "sorry... no fit function found..." << std::endl;
00073   }
00074 }
00075 
00076 // ------------------------------------------------------------------------------------------------------------------
00077 
00078 std::pair<int,int>
00079 MaximalValue::contour(TH1F& hist, double& frac)
00080 {
00081   int idx=hist.GetMaximumBin(), jdx=hist.GetMaximumBin();
00082   if(0<=frac && frac<=1){
00083     while( hist.GetBinContent(idx)/hist.GetMaximum()>frac) --idx;
00084     while( hist.GetBinContent(jdx)/hist.GetMaximum()>frac) ++jdx;
00085   }
00086   return std::pair<int, int>(idx, jdx);
00087 }
00088 
00089 // ------------------------------------------------------------------------------------------------------------------
00090 
00091 double 
00092 Quantile::spreadError(TH1F& hist)
00093 {
00094   quantiles(hist, 0.25+err_); double outer=distance(hist);
00095   quantiles(hist, 0.25-err_); double inner=distance(hist);
00096   return std::fabs(outer-inner)/2;
00097 }