Use the thread safe form of TH*::Fit

The version of TH*::Fit which uses a string to designate the function
to use is not thread safe. The string is used to lookup a global
instance of the TF1 and that instance can be used simultaneously
but multiple Fits to store their results. However, the version of
TH*::Fit which takes an explicit pointer to a TF1 is safe.

Alignment/LaserAlignment/plugins/TkLasBeamFitter.cc

@@ -676,19 +676,19 @@ void TkLasBeamFitter::fitter(TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatri
     TF1 tecPlus("tecPlus", tecPlusFunction, zMin, zMax, nFitParams );
     tecPlus.SetParameter( 1, 0 ); // slope
     tecPlus.SetParameter( nFitParams - 1, 0 ); // BS 
-    lasData->Fit("tecPlus", "R"); // "R", "RV" or "RQ"
+    lasData->Fit(&tecPlus, "R"); // "R", "RV" or "RQ"
   }
   else if(beam.isTecInternal(-1)){
     TF1 tecMinus("tecMinus", tecMinusFunction, zMin, zMax, nFitParams );
     tecMinus.SetParameter( 1, 0 ); // slope
     tecMinus.SetParameter( nFitParams - 1, 0 ); // BS 
-    lasData->Fit("tecMinus", "R");
+    lasData->Fit(&tecMinus, "R");
   }
   else{
     TF1 at("at", atFunction, zMin, zMax, nFitParams );
     at.SetParameter( 1, 0 ); // slope
     at.SetParameter( nFitParams - 1, 0 ); // BS 
-    lasData->Fit("at","R");
+    lasData->Fit(&at,"R");
   }
 
   // get values and errors for offset and slope

Alignment/MuonAlignmentAlgorithms/src/MuonResidualsFitter.cc

@@ -503,7 +503,7 @@ void MuonResidualsFitter::histogramChi2GaussianFit(int which, double &fit_mean,
   f1->SetParameter(0, hist->GetEntries());
   f1->SetParameter(1, 0);
   f1->SetParameter(2, hist->GetRMS());
-  hist->Fit("f1","RQ");
+  hist->Fit(f1,"RQ");
 
   fit_mean  = f1->GetParameter(1);
   fit_sigma = f1->GetParameter(2);

CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzerYousi.cc

@@ -451,7 +451,7 @@ EcalLaserAnalyzerYousi::endJob() {
       min_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) - K*RMS;
       max_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) + K*RMS;
       gs1 = new TF1("gs1", "gaus", min_r, max_r);
-      C_APD[j]->Fit("gs1", "RQI");
+      C_APD[j]->Fit(gs1, "RQI");
       Sigma = gs1->GetParameter(2);
       K = K*1.5;
       iCount++;
@@ -471,7 +471,7 @@ EcalLaserAnalyzerYousi::endJob() {
       min_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) - K*RMS;
       max_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) + K*RMS;
       gs2 = new TF1("gs2", "gaus", min_r, max_r);
-      C_APDPN[j]->Fit("gs2", "RQI");
+      C_APDPN[j]->Fit(gs2, "RQI");
       Sigma = gs2->GetParameter(2);
       K = K*1.5;
       iCount++;

CalibCalorimetry/EcalPedestalOffsets/src/TPedValues.cc

@@ -240,7 +240,7 @@ int TPedValues::makePlots (TFile * rootFile, const std::string & dirName,
           sprintf (name,"XTL%04d_GAIN%02d",endcapCrystalNumbers[xtl],gainHuman) ;
           graph.GetXaxis()->SetTitle("DAC value");
           graph.GetYaxis()->SetTitle("Average pedestal ADC");
-          graph.Fit("fitFunction","RWQ");
+          graph.Fit(&fitFunction,"RWQ");
           graph.Write (name);
 
           double slope1 = fitFunction.GetParameter(0);

CalibMuon/DTCalibration/src/DTTimeBoxFitter.cc

@@ -84,7 +84,7 @@ pair<double, double> DTTimeBoxFitter::fitTimeBox(TH1F *hTimeBox) {
   if(theVerbosityLevel >= 2)
     option = "";
 
-  hTimeBox->Fit("IntGauss", option.c_str(), "",xFitMin, xFitMax);
+  hTimeBox->Fit(fIntGaus, option.c_str(), "",xFitMin, xFitMax);
 
   // Get fitted parameters
   double mean =  fIntGaus->GetParameter("Mean");

CalibTracker/SiPixelGainCalibration/plugins/SiPixelGainCalibrationAnalysis.cc

@@ -343,7 +343,7 @@ SiPixelGainCalibrationAnalysis::doFits(uint32_t detid, std::vector<SiPixelCalibD
     graph_->SetPoint(ipointtemp,xvals[ipointtemp],yvals[ipointtemp]);
     graph_->SetPointError(ipointtemp,0,yerrvals[ipointtemp]);
   }
-  Int_t tempresult = graph_->Fit("func","FQ0N");
+  Int_t tempresult = graph_->Fit(func_,"FQ0N");
   slope=func_->GetParameter(1);
   slopeerror=func_->GetParError(1);
   intercept=func_->GetParameter(0);
@@ -359,7 +359,7 @@ SiPixelGainCalibrationAnalysis::doFits(uint32_t detid, std::vector<SiPixelCalibD
     for(int ii=0; ii<npoints; ++ii){
       edm::LogWarning("SiPixelGainCalibrationAnalysis")<< "vcal " << xvals[ii] << " response: " << yvals[ii] << "+/-" << yerrvals[ii] << std::endl; 
     } 
-    tempresult = graph_->Fit("func","FQ0NW");
+    tempresult = graph_->Fit(func_,"FQ0NW");
     slope=func_->GetParameter(1);
     slopeerror=func_->GetParError(1);
     intercept = func_->GetParameter(0);

CalibTracker/SiStripChannelGain/plugins/SiStripGainCosmicCalculator.cc

@@ -222,24 +222,24 @@ std::pair<double,double> SiStripGainCosmicCalculator::getPeakOfLandau( TH1F * in
 //  delete landaufit;
 //
   // perform fit with standard landau
-  inputHisto->Fit("landau","0Q");
-  TF1 * fitfunction = (TF1*) inputHisto->GetListOfFunctions()->First();
+  // make our own copy to avoid problems with threads
+  std::unique_ptr<TF1> fitfunction( new TF1("landaufit","landau") );
+  inputHisto->Fit(fitfunction.get(),"0Q");
   adcs = fitfunction->GetParameter("MPV");
   error = fitfunction->GetParError(1); // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
   double chi2 = fitfunction->GetChisquare();
   double ndf = fitfunction->GetNDF();
   double chi2overndf = chi2 / ndf;
   // in case things went wrong, try to refit in smaller range
   if(adcs< 2. || (error/adcs)>1.8 ){
-     inputHisto->Fit("landau","0Q",0,0.,400.);
-     TF1 * fitfunction2 = (TF1*) inputHisto->GetListOfFunctions()->First();
+     inputHisto->Fit(fitfunction.get(),"0Q",0,0.,400.);
      std::cout<<"refitting landau for histogram "<<inputHisto->GetTitle()<<std::endl;
      std::cout<<"initial error/adcs ="<<error<<" / "<<adcs<<std::endl;
-     std::cout<<"new     error/adcs ="<<fitfunction2->GetParError(1)<<" / "<<fitfunction2->GetParameter("MPV")<<std::endl;
-     adcs = fitfunction2->GetParameter("MPV");
-     error = fitfunction2->GetParError(1); // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
-     chi2 = fitfunction2->GetChisquare();
-     ndf = fitfunction2->GetNDF();
+     std::cout<<"new     error/adcs ="<<fitfunction->GetParError(1)<<" / "<<fitfunction->GetParameter("MPV")<<std::endl;
+     adcs = fitfunction->GetParameter("MPV");
+     error = fitfunction->GetParError(1); // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
+     chi2 = fitfunction->GetChisquare();
+     ndf = fitfunction->GetNDF();
      chi2overndf = chi2 / ndf;
    }
    // if still wrong, give up

CalibTracker/SiStripChannelGain/plugins/SiStripGainFromCalibTree.cc

@@ -376,7 +376,7 @@ void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResul
    // perform fit with standard landau
    TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
    MyLandau->SetParameter(1,300);
-   InputHisto->Fit("MyLandau","0QR WW");
+   InputHisto->Fit(MyLandau,"0QR WW");
 
    // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
    FitResults[0]         = MyLandau->GetParameter(1);  //MPV

CalibTracker/SiStripChannelGain/plugins/SiStripGainFromData.cc

@@ -1255,15 +1255,14 @@ void SiStripGainFromData::getPeakOfLandau(TH1* InputHisto, double* FitResults, d
    TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
    MyLandau->SetParameter("MPV",300);
 
-   InputHisto->Fit("MyLandau","QR WW");
-   TF1 * fitfunction = (TF1*) InputHisto->GetListOfFunctions()->First();
+   InputHisto->Fit(MyLandau,"QR WW");
 
    // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
-    adcs        = fitfunction->GetParameter("MPV");
-    adcs_err    = fitfunction->GetParError(1);
-    width       = fitfunction->GetParameter(2);
-    width_err   = fitfunction->GetParError(2);
-    chi2overndf = fitfunction->GetChisquare() / fitfunction->GetNDF();
+    adcs        = MyLandau->GetParameter("MPV");
+    adcs_err    = MyLandau->GetParError(1);
+    width       = MyLandau->GetParameter(2);
+    width_err   = MyLandau->GetParError(2);
+    chi2overndf = MyLandau->GetChisquare() / MyLandau->GetNDF();
 
     // if still wrong, give up
     if(adcs<2. || chi2overndf>MaxChi2OverNDF){

CalibTracker/SiStripLorentzAngle/src/LA_Results.cc

@@ -95,9 +95,11 @@ summarize_ensembles(Book& book) const {
       book.fill( p.measured.second, name+"_merr",         500, 0, 0.01);
       book.fill( (p.measured.first-p.reco.first)/p.measured.second, name+"_pull", 500, -10,10);
     }
-    book[name+"_measure"]->Fit("gaus","LLQ");
-    book[name+"_merr"]->Fit("gaus","LLQ");
-    book[name+"_pull"]->Fit("gaus","LLQ");
+    //Need our own copy for thread safety
+    TF1 gaus("mygaus","gaus");
+    book[name+"_measure"]->Fit(&gaus,"LLQ");
+    book[name+"_merr"]->Fit(&gaus,"LLQ");
+    book[name+"_pull"]->Fit(&gaus,"LLQ");
   }
 }
 
@@ -137,7 +139,9 @@ offset_slope(const std::vector<LA_Filler_Fitter::EnsembleSummary>& ensembles) {
       yerr.push_back(ensemble.meanMeasured.second);
     }
     TGraphErrors graph(x.size(),&(x[0]),&(y[0]),&(xerr[0]),&(yerr[0]));
-    graph.Fit("pol1");
+    //Need our own copy for thread safety
+    TF1 pol1("mypol1","pol1");
+    graph.Fit(&pol1);
     const TF1*const fit = graph.GetFunction("pol1");
 
     return std::make_pair( std::make_pair(fit->GetParameter(0), fit->GetParError(0)),

CalibTracker/SiStripLorentzAngle/src/SymmetryFit.cc

@@ -1,6 +1,7 @@
 #include "CalibTracker/SiStripLorentzAngle/interface/SymmetryFit.h"
 #include <cmath>
 #include <cassert>
+#include <memory>
 #include "boost/foreach.hpp"
 
 TH1* SymmetryFit::symmetryChi2(std::string basename, const std::vector<TH1*>& candidates, const std::pair<unsigned,unsigned> range)
@@ -128,12 +129,12 @@ int SymmetryFit::fit() {
       p[0] > chi2_->GetBinCenter(chi2_->GetNbinsX()))
     return 7;
 
-  TF1* f = fitfunction();
+  std::unique_ptr<TF1> f( fitfunction() );
   f->SetParameter(0, p[0]);  f->SetParLimits(0, p[0], p[0]);
   f->SetParameter(1, p[1]);  f->SetParLimits(1, p[1], p[1]);
   f->SetParameter(2, p[2]);  f->SetParLimits(2, p[2], p[2]);
   f->SetParameter(3, ndf_);  f->SetParLimits(3, ndf_,ndf_); //Fixed
-  chi2_->Fit(f,"WQ");
+  chi2_->Fit(f.get(),"WQ");
   return 0;
 }
 
@@ -151,12 +152,14 @@ TF1* SymmetryFit::fitfunction()
 std::vector<double> SymmetryFit::pol2_from_pol2(TH1* hist) {
   std::vector<double> v;
 
-  int status = hist->Fit("pol2","WQ");
+  //Need our own copy for thread safety
+  TF1 func("mypol2","pol2"); 
+  int status = hist->Fit(&func,"WQ");
   if(!status) {
     std::vector<double> p;
-    p.push_back(hist->GetFunction("pol2")->GetParameter(0));
-    p.push_back(hist->GetFunction("pol2")->GetParameter(1));
-    p.push_back(hist->GetFunction("pol2")->GetParameter(2));
+    p.push_back(func.GetParameter(0));
+    p.push_back(func.GetParameter(1));
+    p.push_back(func.GetParameter(2));
     if(p[2]>0) {
       v.push_back( -0.5*p[1]/p[2] );
       v.push_back( 1./sqrt(p[2]) );
@@ -169,13 +172,14 @@ std::vector<double> SymmetryFit::pol2_from_pol2(TH1* hist) {
 std::vector<double> SymmetryFit::pol2_from_pol3(TH1* hist) {
   std::vector<double> v;
 
-  int status = hist->Fit("pol3","WQ");
+  std::unique_ptr<TF1> func( new TF1("mypol3","pol3")); 
+  int status = hist->Fit(func.get(),"WQ");
   if(!status) {
     std::vector<double> p;
-    p.push_back(hist->GetFunction("pol3")->GetParameter(0));
-    p.push_back(hist->GetFunction("pol3")->GetParameter(1));
-    p.push_back(hist->GetFunction("pol3")->GetParameter(2));
-    p.push_back(hist->GetFunction("pol3")->GetParameter(3));
+    p.push_back(func->GetParameter(0));
+    p.push_back(func->GetParameter(1));
+    p.push_back(func->GetParameter(2));
+    p.push_back(func->GetParameter(3));
     double radical = p[2]*p[2] - 3*p[1]*p[3] ;
     if(radical>0) {
       double x0 = ( -p[2] + sqrt(radical) ) / ( 3*p[3] ) ;

Calibration/EcalCalibAlgos/src/PhiSymmetryCalibration.cc

@@ -464,15 +464,16 @@ void PhiSymmetryCalibration::getKfactors()
   std::vector<TGraph*>  k_barl_graph(kBarlRings);
   std::vector<TCanvas*> k_barl_plot(kBarlRings);
 
+  //Create our own TF1 to avoid threading problems
+  TF1 mypol1("mypol1","pol1");
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     for (int imiscal=0; imiscal<kNMiscalBinsEB; imiscal++) {
       int middlebin =  int (kNMiscalBinsEB/2);
       epsilon_T_eb[imiscal] = etsum_barl_miscal_[imiscal][ieta]/etsum_barl_miscal_[middlebin][ieta] - 1.;
       epsilon_M_eb[imiscal] = miscalEB_[imiscal] - 1.;
     }
     k_barl_graph[ieta] = new TGraph (kNMiscalBinsEB,epsilon_M_eb,epsilon_T_eb);
-    k_barl_graph[ieta]->Fit("pol1");
-
+    k_barl_graph[ieta]->Fit(&mypol1);
 
     ostringstream t;
     t<< "k_barl_" << ieta+1; 
@@ -504,7 +505,7 @@ void PhiSymmetryCalibration::getKfactors()
       epsilon_M_ee[imiscal] = miscalEE_[imiscal] - 1.;
     }
     k_endc_graph[ring] = new TGraph (kNMiscalBinsEE,epsilon_M_ee,epsilon_T_ee);
-    k_endc_graph[ring]->Fit("pol1");
+    k_endc_graph[ring]->Fit(&mypol1);
 
     ostringstream t;
     t<< "k_endc_"<< ring+1;

Calibration/Tools/src/ZIterativeAlgorithmWithFit.cc

@@ -365,11 +365,11 @@ ZIterativeAlgorithmWithFit::~ZIterativeAlgorithmWithFit()
 }
 
 void ZIterativeAlgorithmWithFit::gausfit(TH1F * histoou,double* par,double* errpar,float nsigmalow, float nsigmaup, double* myChi2, int* iterations) {
-  TF1 *gausa = new TF1("gausa","gaus",histoou->GetMean()-3*histoou->GetRMS(),histoou->GetMean()+3*histoou->GetRMS());
+  std::unique_ptr<TF1> gausa{ new TF1("gausa","gaus",histoou->GetMean()-3*histoou->GetRMS(),histoou->GetMean()+3*histoou->GetRMS()) };
 
   gausa->SetParameters(histoou->GetMaximum(),histoou->GetMean(),histoou->GetRMS());
 
-  histoou->Fit("gausa","qR0N");
+  histoou->Fit(gausa.get(),"qR0N");
 
   double p1    = gausa->GetParameter(1);
   double sigma = gausa->GetParameter(2);
@@ -383,7 +383,6 @@ void ZIterativeAlgorithmWithFit::gausfit(TH1F * histoou,double* par,double* errp
   double xmax_fit=p1+nsigmaup*sigma;
 
   int iter=0;
-  TF1* fitFunc;
 
   while ((chi2>1. && iter<5) || iter<2 )
     {
@@ -394,12 +393,12 @@ void ZIterativeAlgorithmWithFit::gausfit(TH1F * histoou,double* par,double* errp
 
       char suffix[20];
       sprintf (suffix,"_iter_%d",iter); 
-      fitFunc = new TF1("FitFunc"+TString(suffix),"gaus",xmin_fit,xmax_fit);
+      std::unique_ptr<TF1> fitFunc{ new TF1("FitFunc"+TString(suffix),"gaus",xmin_fit,xmax_fit) };
       fitFunc->SetParameters(nor,p1,sigma);
       fitFunc->SetLineColor((int)(iter+1));
       fitFunc->SetLineWidth(1);
       //histoou->Fit("FitFunc","lR+","");
-      histoou->Fit("FitFunc"+TString(suffix),"qR0+","");
+      histoou->Fit(fitFunc.get(),"qR0+","");
 
       histoou->GetXaxis()->SetRangeUser(xmi,xma);
       histoou->GetXaxis()->SetLabelSize(0.055);

DQM/BeamMonitor/plugins/BeamMonitor.cc

@@ -33,6 +33,7 @@ V00-03-25
 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
 #include <numeric>
 #include <math.h>
+#include <memory>
 #include <TMath.h>
 #include <iostream>
 #include <TStyle.h>
@@ -854,10 +855,10 @@ void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &n
       sprintf(tmpTitle,"%s %i %s %i","Fitted Primary Vertex (cm) of LS: ",beginLumiOfPVFit_," to ",endLumiOfPVFit_);
       pvResults->setAxisTitle(tmpTitle,1);
 
-      TF1 *fgaus = new TF1("fgaus","gaus");
+      std::unique_ptr<TF1> fgaus{ new TF1("fgaus","gaus") };
       double mean,width,meanErr,widthErr;
       fgaus->SetLineColor(4);
-      h_PVx[0]->getTH1()->Fit("fgaus","QLM0");
+      h_PVx[0]->getTH1()->Fit(fgaus.get(),"QLM0");
       mean = fgaus->GetParameter(1);
       width = fgaus->GetParameter(2);
       meanErr = fgaus->GetParError(1);
@@ -892,10 +893,10 @@ void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &n
       tmphisto = static_cast<TH1D *>((h_PVx[0]->getTH1())->Clone("tmphisto"));
       h_PVx[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
       h_PVx[1] = dbe_->book1D(tmpfile,h_PVx[0]->getTH1F());
-      h_PVx[1]->getTH1()->Fit("fgaus","QLM");
+      h_PVx[1]->getTH1()->Fit(fgaus.get(),"QLM");
 
 
-      h_PVy[0]->getTH1()->Fit("fgaus","QLM0");
+      h_PVy[0]->getTH1()->Fit(fgaus.get(),"QLM0");
       mean = fgaus->GetParameter(1);
       width = fgaus->GetParameter(2);
       meanErr = fgaus->GetParError(1);
@@ -921,10 +922,10 @@ void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &n
       h_PVy[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
       h_PVy[1]->update();
       h_PVy[1] = dbe_->book1D(tmpfile,h_PVy[0]->getTH1F());
-      h_PVy[1]->getTH1()->Fit("fgaus","QLM");
+      h_PVy[1]->getTH1()->Fit(fgaus.get(),"QLM");
 
 
-      h_PVz[0]->getTH1()->Fit("fgaus","QLM0");
+      h_PVz[0]->getTH1()->Fit(fgaus.get(),"QLM0");
       mean = fgaus->GetParameter(1);
       width = fgaus->GetParameter(2);
       meanErr = fgaus->GetParError(1);
@@ -950,7 +951,7 @@ void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &n
       h_PVz[1]->getTH1()->SetBins(tmphisto->GetNbinsX(),tmphisto->GetXaxis()->GetXmin(),tmphisto->GetXaxis()->GetXmax());
       h_PVz[1]->update();
       h_PVz[1] = dbe_->book1D(tmpfile,h_PVz[0]->getTH1F());
-      h_PVz[1]->getTH1()->Fit("fgaus","QLM");
+      h_PVz[1]->getTH1()->Fit(fgaus.get(),"QLM");
 
     }//check if found min Vertices
   }//do PVfit
@@ -1165,10 +1166,10 @@ void BeamMonitor::FitAndFill(const LuminosityBlock& lumiSeg,int &lastlumi,int &n
 
 	double mean = bs.z0();
 	double width = bs.sigmaZ();
-	TF1 *fgaus = new TF1("fgaus","gaus");
+	std::unique_ptr<TF1> fgaus{ new TF1("fgaus","gaus") };
 	fgaus->SetParameters(mean,width);
 	fgaus->SetLineColor(4);
-	h_trk_z0->getTH1()->Fit("fgaus","QLRM","",mean-3*width,mean+3*width);
+	h_trk_z0->getTH1()->Fit(fgaus.get(),"QLRM","",mean-3*width,mean+3*width);
       }
 
       fitResults->Reset();