C++ TF1::GetName方法代码示例

void Draw_KL_Test(){

  TChain* ch = new TChain("Tree");
  TChain* ch1 = new TChain("Tree");
  
  TH1D* his  = new TH1D("Klong6g","Klong6g",20,450,550);
  TH1D* his1 = new TH1D("Klong4g","Klong4g",20,450,550);
  TH1D* his2 = new TH1D("Klong4gAll","Klong4gAll",60,250,550);

  for( int i = 0; i< 68; i++){
    ch->Add(Form("klongRootFile/kl%d.root" ,4162+i));
    ch1->Add(Form("klongRootFile/ks%d.root",4162+i));
  }
  ch->Project(his->GetName()  ,"KlongMass[0]","CutCondition==0");
  ch1->Project(his1->GetName(),"KlongMass[0]","CutCondition==0");
  ch1->Project(his2->GetName(),"KlongMass[0]","CutCondition==0");
  
  TF1* func = new TF1("func","gaus(0)+expo(3)",0,550);
  func->SetParameter(1,498);
  func->SetParameter(2,5);
  TF1* func2 = new TF1("func2","gaus(0)",0,550);
  func2->SetParameter(1,498);
  func2->SetParameter(2,5);



  TCanvas* can = new TCanvas("can","",1200,600);
  can->Divide(2,1);
  can->cd(1);
  his2->Fit(func->GetName(),"","",450,550);
  his2->Draw();
  TF1* func1 = new TF1("Test","gaus",450,550);
  func1->SetParameter(0,func->GetParameter(0));
  func1->SetParameter(1,func->GetParameter(1));
  func1->SetParameter(2,func->GetParameter(2));

  can->cd(2);
  his1->SetLineColor(2);  
  his->Draw();
  his->Fit(func2->GetName(),"","",450,550);
  func->Draw("same");
  his1->Draw("same");

  std::cout<< func2->GetParameter(0) << " " 
	   << func->GetParameter(0)  << " " 
	   << func->GetParameter(0)/func2->GetParameter(0)<< std::endl;
  std::cout<< func2->Integral(450,550) << " " 
	   << func1->Integral(450,550)  << " " 
	   << func1->Integral(450,550)/func2->Integral(450,550)
	   << std::endl;
  //ch->Draw("KlongPt[0]:KlongMass[0]>>(400,200,600,50,0,20)","(CutCondition&(1|2|4|8))==0","colz");
  gPad->SetLogz();
  TText* text = new TText(0.5,0.5,"");
  TText* text1 = new TText(0.5,0.5,"");
  text->DrawTextNDC(0.5,0.5,Form("Integral:%2.3lf",func1->Integral(450,550)));
  text1->DrawTextNDC(0.5,0.6,Form("Integral:%2.3lf",func2->Integral(450,550)));
  
}

void FitWithConstant(TH1D* h1, TFile *out)
{
  Double_t fitRangeLow = 0.;
  Double_t fitRangeHigh = .4;
  TF1 *f = new TF1("f","[0]", fitRangeLow, fitRangeHigh);

  Int_t binLow = h1->FindBin(fitRangeLow);
  Int_t binHigh = h1->FindBin(fitRangeHigh);
  cout<<"\t\t********************Fitting with a constant value********************"<<endl<<endl;;
  h1->Fit(f, "R0");

  TString fitName = h1->GetName();
  fitName += "Fit";
  f->SetName(fitName);
  f->SetTitle(fitName);

  // TFile out("Compare.root","update");
  TDirectory *dir = out->GetDirectory("Fit");
  if(!dir) dir = out->mkdir("Fit");
  dir->cd();
  // f->SetDirectory(0); 
  f->Write(f->GetName(), TObject::kOverwrite);
  // out->Close();
}

//.........这里部分代码省略.........
            pad = (TPad*)canvas[i]->cd(nfile+1);
            histos[i]->SetLineColor(1);
            histos[i]->SetLineWidth(1);
            histos[i]->GetXaxis()->SetRangeUser(0.7, 1.2);
            histos[i]->GetYaxis()->SetRangeUser(0, 1.5*histos[i]->GetBinContent(histos[i]->GetMaximumBin()));
            //histos[i]->SetStats(kFALSE);
            
            //histos[i]->Draw("HIST");

            printf("mean PT: %f\n", meanPT);

            TF1 *fit = new TF1(Form("fitPTbin%d00particle%d", nfile*2+1, i), FitFunRelBW, 0.68, 1.05, 5);
            //TF1 *fit = new TF1(Form("fitPTbin%d00particle%d", nfile*2+1, i), "gaus(0)", 0.86, 0.92);


            fit->SetParNames("BW Area", "Mass", "Width", "PT", "Temp");
            fit->SetParameters(TMath::Power(10.0, (float)(nfile)/1.7), 0.89, 0.1, 0.5, 0.130);
            //fit->SetParNames("BW Area", "Mass", "Width");
            //fit->SetParameters(100, 0.89, 0.0474);
            //fit->SetParLimits(0, -10, 1.5e9);
            Float_t max = histos[i]->GetXaxis()->GetBinCenter(histos[i]->GetMaximumBin());
            //if(max < 0.91 && max > 0.892){
            //    fit->SetParLimits(1, max-0.001, max+0.001);
            //}else{
                fit->SetParLimits(1, 0.82, 0.98);
            //}
            //fit->SetParLimits(2, 0.005, 0.15);
            fit->FixParameter(2, 0.042);
            fit->FixParameter(3, meanPT);
            //fit->SetParLimits(4, 0.05, 0.2);
            fit->FixParameter(4, 0.100001);
            fit->SetLineColor(2);

            printf("%s\n", fit->GetName());

            histos[i]->Fit(Form("fitPTbin%d00particle%d", nfile*2+1, i), "BRIM", "SAME");
            TVirtualFitter *fitter = TVirtualFitter::GetFitter();
           
            histos[i]->SetStats(1);
            histos[i]->Draw();
            gPad->Update();
            pad->Update();
            st = (TPaveStats*)histos[i]->FindObject("stats");
            st->SetX1NDC(0.524);
            st->SetY1NDC(0.680);
            st->SetX2NDC(0.884);
            st->SetY2NDC(0.876);
            //fit->Draw("SAME");
            //histos[i]->Draw();
            gPad->Update();
            pad->Update();
            printf("\n");
    
            diffHistos[i] = (TH1D*)histos[i]->Clone(Form("diffPTbin%d00particl%d", nfile*2+1, i));
            diffHistos[i]->Add(fit, -1);
            diffCanvas[i]->cd(nfile+1);
            diffHistos[i]->Draw("HIST E");
            diffHistos[i]->Write();

            //counting bins
            Float_t integral = histos[i]->Integral(1, 500)*500.0;
            //integrals << integral <<" \n";
            histos[i]->Write();
            fit->Write();
            //Do mass and width vs. pT plots just for K*0
            if(i==PARTICLE_NUM){

bool leptonic_fitter_algebraic::fit( const TLorentzVector& B, const TH1& BITF, const TF1& Beff, 
				     const TLorentzVector& lep, 
				     double MEX, double MEY, const TF1& dnuPDF )
{
  if( _dbg > 19 ) cout<<"DBG20 Entered leptonic_fitter_algebraic::fit with B mass: "<<B.M()<<", l_m:"<<lep.M()<<", MET: "<<MEX<<" "<<MEY<<endl;
  if( B.M() <= 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a b-jet with an illegal (non-positive) mass!"); 
  if( lep.M() < 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a lepton with an illegal (negative) mass!"); 
  _converged = _swapped = false;
  _obsB = B;
  _obsL = lep;

  _BITF = &BITF;
  _Beff = &Beff;
  _dnuPDF = dnuPDF;

  _b_m2 = B.M2();

  double lep_b_angle = lep.Angle( B.Vect() );
  double cos_lep_b = TMath::Cos( lep_b_angle );
  double sin_lep_b = TMath::Sin( lep_b_angle );
  double b_p = B.P();
  double b_e = B.E();
  _denom = b_e - cos_lep_b * b_p;
  
  _lep_p = lep.P();
  _x0 = - _W_m2 / ( 2 * _lep_p );
  _y1 = - sin_lep_b * _x0 * b_p / _denom;
  _x1_0 = _x0 * b_e / _denom  -  _y1*_y1 / _x0;
  _Z2_0 = _x0*_x0 - _W_m2 - _y1*_y1;
  if( _dbg > 219 ) cout<<"DBG220 lfa updated lepton with: "<<lv2str( lep )<<" -> x0:"<<_x0<<", y1: "<<_y1<<", x1_0: "<<_x1_0<<", Z2_0: "<<_Z2_0<<endl;

  static double bnums[3];
  bnums[0] = B.X();
  bnums[1] = B.Y();
  bnums[2] = B.Z();
  TMatrixD bXYZ( 3, 1, bnums );
  _R_T = rotation( 2, lep.Phi() ); // R_z^T
  _R_T *= rotation( 1, lep.Theta() - 0.5*TMath::Pi() ); // R_z^T R_y^T
  TMatrixD rotation_vect( _R_T, TMatrixD::kTransposeMult, bXYZ ); // R_y R_z
  double* rotation_array = rotation_vect.GetMatrixArray();
  double phi_x = - TMath::ATan2( rotation_array[2], rotation_array[1] );
  if( _dbg > 99 ) cout<<"DBG100 lfa x rotation vector is:"<<rotation_array[0]<<" "<<rotation_array[1]<<" "<<rotation_array[2]<<" -> phi_x:"<<phi_x<<endl;
  _R_T *= rotation( 0, - phi_x ); // R_z^T R_y^T R_x^T

  // set up _Nu's non-zero elements so that \vec{nu} = Nu \vec{t} for any \vec{t} (since only t's 3nd component is used, and its always 1).
  _Nu[0][2] = MEX;
  _Nu[1][2] = MEY;

  double iVarMET = TMath::Power( TMath::Max( 1., dnuPDF.GetHistogram()->GetRMS() ), -2 );
  _invFlatVar[0][0] = _invFlatVar[1][1] = iVarMET; // set up the chi^2 distance with the right order of magnitude (generalizes to rotated covariance matrix)
  if( _dbg > 209 ) cout<<"DBG210 lfa "<<dnuPDF.GetName()<<" --> iVarMET:"<<iVarMET<<endl;

  // (re)define fit parameter, so all fits start off on an equal footing
  _mini->SetPrintLevel( _minimizer_print_level );
  _mini->Clear();
  _mini->SetFunction( _functor );
  leptonic_fitter_algebraic_object = this; // set the function in the functor pointing back to this object. Doubtfull that all this redirection is needed...
  _mini->SetTolerance( _tolerance );
  bool OK = _mini->SetLimitedVariable( 0, "sB", 1.0, 0.4, 0.1, 6.0 );
  //bool OK = _mini->SetVariable( 0, "sB", 1.0, 0.4 );
  if( ! OK ) {cerr<<"minimizer (@lfa) failed to SetVariable."<<endl; return false;}

  // define 1 sigma in terms of the function
  _mini->SetErrorDef( 0.5 ); // since this is a likelihood fit

  // do the minimization
  OK = _mini->Minimize(); 
  if( _dbg > 19 && ( ! OK || _dbg > 59 ) ) cout<<"DBG INFO: initial fit @lfa returned OK: "<<OK<<", has status: "<<_mini->Status()<<endl;

  _converged = OK; // use status somehow? depends on fitter?

  // read parameters
  const double *xs = _mini->X();
  for( int ip = 0; ip < 1; ++ip ) _params[ ip ] = xs[ ip ];

  // return all intermediate results to the minimum, in particular, the discriminant
  calc_MLL( _params, true );
  TMatrixD nu_vec( _Emat, TMatrixD::kMult, _tvec );
  update_nu_and_decay_chain( nu_vec );
  if( _dbg > 203 ) cout<<"DBG204 lfa finalized _genN: "<<lv2str(_genN)<<", _W: "<<lv2str(_W)<<", & _t: "<<lv2str(_T)<<endl;

  _MLL = _mini->MinValue();
  return true;
} 

//.........这里部分代码省略.........
      if (factor / background < 10) continue;

      // Get the total number of events in the gaussian 1 sigma
      Int_t sigRangeLow = hist->FindBin(mean - sigma);
      Int_t sigRangeHigh = hist->FindBin(mean + sigma);

      double sigRangeTotal = 0;
      for (Int_t bin = sigRangeLow; bin <= sigRangeHigh; bin++)
        sigRangeTotal += hist->GetBinContent(bin);

      // 2 * 1 sigma integral shoudl give ~ area under gaussian
      sigma /= sqrt(2 * sigRangeTotal);

      ptsX.push_back(hist->GetYaxis()->GetBinCenter(row));
      ptsY.push_back(mean);
      ptsErr.push_back(sigma);
    }

    if (ptsX.size() < 3) continue;

    std::vector<double> yvals = ptsY;
    std::sort(yvals.begin(), yvals.end());
    const double median = yvals[yvals.size()/2];
    double avgDeviation = 0;
    for (unsigned int i = 0; i < yvals.size(); i++)
      avgDeviation += fabs(yvals[i] - median);
    avgDeviation /= (double)yvals.size();

    std::vector<double> ptsXGood;
    std::vector<double> ptsYGood;
    std::vector<double> ptsErrGood;

    for (unsigned int i = 0; i < ptsX.size(); i++)
    {
      if (fabs(ptsY[i] - median) > 1.5*avgDeviation) continue;
      ptsXGood.push_back(ptsX[i]);
      ptsYGood.push_back(ptsY[i]);
      ptsErrGood.push_back(ptsErr[i]);
    }

    if (ptsXGood.size() < 3) continue;

    TGraphErrors* graph = new TGraphErrors(ptsXGood.size(),
                                           &(ptsXGood.at(0)),
                                           &(ptsYGood.at(0)), 0,
                                           &(ptsErrGood.at(0)));

    TF1* fitFunc = new TF1("f1", "1 ++ x");
    TF1* result = 0;

    graph->Fit(fitFunc, "Q0E").Get();
    result = graph->GetFunction(fitFunc->GetName());

    // Weight the angle by the slope uncertainty and the inverse of the chi2 normalized
    double weight = result->GetParError(1);
    const double chi2 = result->GetChisquare() / (double)result->GetNDF();
    fitChi2 += chi2;
    weight *= chi2;
    if (weight > 10 * DBL_MIN) weight = 1.0 / weight;
    else weight = 1.0;

    if (axis)
    {
      rotation -= weight * atan(result->GetParameter(1));
      offsetX = result->GetParameter(0);
    }
    else
    {
      rotation += weight * atan(result->GetParameter(1));
      offsetY = result->GetParameter(0);
    }

    angleWeights += weight;

    if (display)
    {
      TCanvas* can = new TCanvas("ResidualAlignment", "Residual Alignment", 900, 600);
      can->Divide(2);
      can->cd(1);
      hist->Draw("COLZ");
      can->cd(2);
      result->SetLineColor(46);
      result->SetLineWidth(2);
      graph->Draw("ap");
      result->Draw("SAME");
      can->Update();
      can->WaitPrimitive();
    }

    delete fitFunc;
    delete graph;
  }

  if (angleWeights > 10 * DBL_MIN)
    rotation /= angleWeights;
  std::cout << "relaxation: " << relaxation << std::endl;
  rotation *= relaxation;
  offsetX *= relaxation;
  offsetY *= relaxation;
}

//.........这里部分代码省略.........
	  tex->DrawLatex(0.55,0.27-0.045*isample,
			 Form("#chi^{2}/NDF = %1.1f / %d",
			      fk->GetChisquare(), fk->GetNDF()));
	  tex->SetTextColor(kBlack);

	  // Error band
	  const int n = fk->GetNpar();
	  TMatrixD emat(n,n);
	  gMinuit->mnemat(emat.GetMatrixArray(), n);
	  TF1 *fke = new TF1(Form("fk_%s_%s_%s",ct,cm,cs),
			     sr_fitError, 30, 1300, 1);
	  _sr_fitError_func = fk;
	  _sr_fitError_emat = &emat;

	  fke->SetLineStyle(kSolid);
	  fke->SetLineColor(fk->GetLineColor()-10);
	  fke->SetParameter(0,-1);
	  fke->DrawClone("SAME");
	  fke->SetParameter(0,+1);
	  fke->DrawClone("SAME");

	  fk->DrawClone("SAME");
	  gk->DrawClone("SAME Pz");

	  // Store soft radiation corrections in fsr subdirectory
	  assert(fin->cd(ct));
	  assert(gDirectory->cd(bin));
	  if (!gDirectory->FindObject("fsr")) gDirectory->mkdir("fsr");
	  assert(gDirectory->cd("fsr"));

	  TH1D *hk = (TH1D*)(isample==0 ? hpt2->Clone() : hpt1->Clone());
	  hk->SetName(Form("hkfsr_%s_%s",cm,cs));
	  TProfile *ppt = (isample==0 ? ppt2 : ppt1);
	  if (isample==3) { ppt = ppt4; } // pas-v6
	  for (int i = 1; i != hk->GetNbinsX()+1; ++i) {
	    double pt = ppt->GetBinContent(i);
	    if (pt>30 && pt<1300) {
	      hk->SetBinContent(i, fk->Eval(pt));
	      hk->SetBinError(i, fabs(fke->Eval(pt)-fk->Eval(pt)));
	    }
	    else {
	      hk->SetBinContent(i, 0);
	      hk->SetBinError(i, 0);
	    }
	  }
	  
	  hk->Write(hk->GetName(), TObject::kOverwrite);

	  // Factorize error matrix into eigenvectors
	  // Remember: A = Q*Lambda*Q^-1, where
	  // A is emat, Q is eigmat, and Lambda is a diagonal matrix with
	  // eigenvalues from eigvec on the diagonal. For eigenmatrix
	  // Q^-1 = Q^T, i.e. inverse matrix is the original transposed
	  TVectorD eigvec(n);
	  TMatrixD eigmat = emat.EigenVectors(eigvec);

	  // Eigenvectors are the columns and sum of eigenvectors squared
	  // equals original uncertainty. Calculate histograms from the
	  // eigenvectors and store them
	  TF1 *fkeig = (TF1*)fk->Clone(Form("%s_eig",fk->GetName()));
	  fkeig->SetLineStyle(kDotted);
	  for (int ieig = 0; ieig != n; ++ieig) {

	    // Eigenvector functions
	    for (int i = 0; i != n; ++i) {
	      fkeig->SetParameter(i, fk->GetParameter(i)
				  + eigmat[i][ieig] * sqrt(eigvec[ieig]));
	    }
	    fkeig->DrawClone("SAMEL");

	    // Eigenvector histograms evaluated at bin mean pT
	    TH1D *hke = (TH1D*)hk->Clone(Form("%s_eig%d",hk->GetName(),ieig));
	    hke->Reset();

	    for (int i = 0; i != gk->GetN(); ++i) {

	      double pt = gk->GetX()[i];
	      int ipt = hke->FindBin(pt);
	      // Need to store central value as well, because
	      // uncertainty sources are signed
	      hke->SetBinContent(ipt, fkeig->Eval(pt)-fk->Eval(pt));
	      hke->SetBinError(ipt, fabs(fkeig->Eval(pt)-fk->Eval(pt)));
	    }
	    hke->Write(hke->GetName(), TObject::kOverwrite);
	  }

	  cout << "." << flush;
	} // if tree
      } // for isample
    } // for imethod
  } // for itype
  
  c3->cd(0);
  //cmsPrel(_lumi, true);
  CMS_lumi(c3, 2, 33);
  c3->SaveAs("pdf/softrad_2x6_kfsr.pdf");

  fin->Close();
  curdir->cd();
} // softrad