C++ TH1::SetAxisRange方法代码示例

/*============================================================================*/
void gaus1peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[5],epar[5],x[4],y[4];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  setcanvas(1);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x2;
  x[2] = x3;
  x[3] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x2);
  y[1] = hist->GetBinContent(bin2);
  Int_t bin3 = hist->FindBin(x3);
  y[2] = hist->GetBinContent(bin3);
  Int_t bin4 = hist->FindBin(x4);
  y[3] = hist->GetBinContent(bin4);
  TGraph *g = new TGraph(4,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[3]=fpol1->GetParameter(0);
  par[4]=fpol1->GetParameter(1);

  //--**-- Gaussian Peak estimation without background --**--//
  TF1 *fgaus = new TF1("GAUS","gaus",x2,x3);
  hist->Fit(fgaus,"RQN");
  fgaus->GetParameters(&par[0]);

  //--**-- Final Peak Fit with Background --**--//
  TF1 *func = new TF1("FGAUS","gaus(0)+pol1(3)",x1,x4);
  func->SetParameters(par);
  hist->Fit(func,"R+QN");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  Double_t fwhm = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N0 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r0 = epar[0]/par[0];
  Double_t r2 = epar[2]/par[2];
  Double_t eN0= N0*TMath::Sqrt(r0*r0+r2*r2);
  printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
          par[1],epar[1],fwhm,efwhm,N0,eN0);
  //printf("%11.4f %11.4f %11.0f %11.0f\n",
  //        par[1],epar[1],N0,eN0);
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
}

/*============================================================================*/
void gaus2peakfit(Char_t *s, Float_t x1, Float_t x2, Float_t x3, Float_t x4)
{
  Double_t par[8],epar[8],x[2],y[2];
  TH1 *hist;
  hist = (TH1 *) gROOT->FindObject(s);
  TCanvas *c1=(TCanvas*) gROOT->FindObject("c1");
  if(c1==NULL)setcanvas(1);
  c1->Clear();
  hist->SetAxisRange(x1-30,x4+30);
  hist->Draw();

  //--**-- Linear background estimation --**--//
  x[0] = x1;
  x[1] = x4;
  Int_t bin1 = hist->FindBin(x1);
  y[0] = hist->GetBinContent(bin1);
  Int_t bin2 = hist->FindBin(x4);
  y[1] = hist->GetBinContent(bin2);
  TGraph *g = new TGraph(2,x,y);
  TF1 *fpol1 = new TF1("POL1","pol1",x1,x4);
  g->Fit(fpol1,"RQN");
  par[6]=fpol1->GetParameter(0);
  par[7]=fpol1->GetParameter(1);

  //--**-- Two Gaussian Peak estimation without background --**--//
  fgaus1 = new TF1("m1","gaus",x1,x2);
  fgaus2 = new TF1("m2","gaus",x3,x4);
  hist->Fit(fgaus1,"R+QN");
  hist->Fit(fgaus2,"R+QN");
  fgaus1->GetParameters(&par[0]);
  fgaus2->GetParameters(&par[3]);

  //--**-- Final Peak Fit with Background --**--//
  func = new TF1("m","gaus(0)+gaus(3)+pol1(6)",x1,x4);
  func->SetParameters(par);

  
  hist->Fit(func,"R+QN");
  func->SetLineWidth(0.5);
  func->SetLineStyle(1);
  func->SetLineColor(4);
  func->SetFillColor(4);
  func->Draw("same");
  func->GetParameters(par);
  epar[0]=func->GetParError(0);
  epar[1]=func->GetParError(1);
  epar[2]=func->GetParError(2);
  epar[3]=func->GetParError(3);
  epar[4]=func->GetParError(4);
  epar[5]=func->GetParError(5);
 
  Double_t fwhm1 = par[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm1 = epar[2]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N10 = par[0]*(TMath::Sqrt(TMath::TwoPi())*par[2]);
  Double_t r10 = epar[0]/par[0];
  Double_t r12 = epar[2]/par[2];
  Double_t eN10= N10*TMath::Sqrt(r10*r10+r12*r12);

  Double_t fwhm2 = par[5]*TMath::Sqrt(8*TMath::Log(2));
  Double_t efwhm2 = epar[5]*TMath::Sqrt(8*TMath::Log(2));
  Double_t N20 = par[3]*(TMath::Sqrt(TMath::TwoPi())*par[5]);
  Double_t r20 = epar[3]/par[3];
  Double_t r22 = epar[5]/par[5];
  Double_t eN20= N20*TMath::Sqrt(r20*r20+r22*r22);

  //printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
  //        par[1],epar[1],fwhm1,efwhm1,N10,eN10);
  //printf("Peak = %f +- %f; FFHM = %f +- %f; Area = %f +- %f\n",
  //        par[4],epar[4],fwhm2,efwhm2,N20,eN20);
  printf("%11.4f %11.4f %11.0f %11.0f\n",
          par[1],epar[1],N10,eN10);
  printf("%11.4f %11.4f %11.0f %11.0f\n",
          par[4],epar[4],N20,eN20);
}

void Template::makeResidualsControlPlot(const string& tag, unsigned int rebin)
/*****************************************************************/
{
    if(numberOfDimensions()>0 && m_template->GetNbinsX()%rebin!=0) return;
    if(numberOfDimensions()>1 && m_template->GetNbinsY()%rebin!=0) return;
    if(numberOfDimensions()>2 && m_template->GetNbinsZ()%rebin!=0) return;

    stringstream cpName, cpRawName, resMapName, resDistName, relErrDistName;
    cpName << m_name << "_cp";
    cpRawName << m_name << "_rawcp";
    resMapName << m_name << "_resmap_" << tag << "_rebin" << rebin;
    resDistName << m_name << "_resdist_" << tag << "_rebin" << rebin;
    relErrDistName << m_name << "_relerrdist_" << tag << "_rebin" << rebin;
    TH1* cpTmp = NULL;
    TH1* cpRawTmp = NULL;
    if(rebin==1)
    {
        cpTmp = dynamic_cast<TH1*>(m_template->Clone(cpName.str().c_str()));
        cpRawTmp = dynamic_cast<TH1*>(m_rawTemplate->Clone(cpRawName.str().c_str()));
    }
    else
    {
        if(numberOfDimensions()==2)
        {
            cpTmp = dynamic_cast<TH2F*>(m_template)->Rebin2D(rebin, rebin, cpName.str().c_str());
            cpRawTmp = dynamic_cast<TH2F*>(m_rawTemplate)->Rebin2D(rebin, rebin, cpName.str().c_str());
        }
        else if(numberOfDimensions()==3)
        {
            cpTmp = dynamic_cast<TH3F*>(m_template)->Rebin3D(rebin, rebin, rebin, cpName.str().c_str());
            cpRawTmp = dynamic_cast<TH3F*>(m_rawTemplate)->Rebin3D(rebin, rebin, rebin, cpName.str().c_str());
        }
    }
    TH1* resMap = dynamic_cast<TH1*>(cpTmp->Clone(resMapName.str().c_str()));
    TH1D* resDist = new TH1D(resDistName.str().c_str(), resDistName.str().c_str(), 30, -3, 3);
    resDist->StatOverflows();
    TH1D* relErrDist = new TH1D(relErrDistName.str().c_str(), relErrDistName.str().c_str(), 200, -1, 1);
    relErrDist->StatOverflows();
    unsigned int nbins1 = cpTmp->GetNbinsX();
    unsigned int nbins2 = cpTmp->GetNbinsY();
    unsigned int nbins3 = cpTmp->GetNbinsZ();
    for(unsigned int b1=1;b1<=nbins1;b1++)
    {
        for(unsigned int b2=1;b2<=nbins2;b2++)
        {
            for(unsigned int b3=1;b3<=nbins3;b3++)
            {
                double tmpValue = cpTmp->GetBinContent(b1,b2,b3);
                double tmpRawValue = cpRawTmp->GetBinContent(b1,b2,b3);
                double tmpRawError = cpRawTmp->GetBinError(b1,b2,b3);
                if(tmpValue>0. && tmpRawValue>0. && tmpRawError>0.)
                {
                    double res = (tmpRawValue-tmpValue)/tmpRawError;
                    double relErr = (tmpRawValue-tmpValue)/tmpRawValue;
                    if(numberOfDimensions()==2)
                    {
                        resMap->SetBinContent(b1,b2,b3, res);
                        resMap->SetBinError(b1,b2,b3, 0.);
                    }
                    resDist->Fill(res);
                    relErrDist->Fill(relErr);
                }
            }
        }
    }
    stringstream plotMapName, plotDistName, plotErrDistName;
    plotMapName << "control_" << getName() << "_resMap" << "_" << tag << "_rebin" << rebin;
    plotDistName << "control_" << getName() << "_resDist" << "_" << tag << "_rebin" << rebin;
    plotErrDistName << "control_" << getName() << "_relErrDist" << "_" << tag << "_rebin" << rebin;
    if(numberOfDimensions()==2)
    {
        TCanvas* c = new TCanvas(plotMapName.str().c_str(),plotMapName.str().c_str(), 700,700);
        resMap->SetContour(99);
        resMap->SetAxisRange(-3., 3., "z");
        resMap->Draw("color z");
        resMap->SetXTitle(getVariable(0).c_str());
        resMap->SetYTitle(getVariable(1).c_str());
        addControlPlot(c);
    }
    TCanvas* c2 = new TCanvas(plotDistName.str().c_str(),plotDistName.str().c_str(), 700,700);
    resDist->SetLineColor(kBlack);
    resDist->SetLineWidth(2);
    resDist->SetMarkerColor(kBlack);
    resDist->SetMarkerStyle(20);
    resDist->SetXTitle("(raw-template)/error_{raw}");
    resDist->Draw();
    addControlPlot(c2);

    TCanvas* c3 = new TCanvas(plotErrDistName.str().c_str(),plotErrDistName.str().c_str(), 700,700);
    relErrDist->SetLineColor(kBlack);
    relErrDist->SetLineWidth(2);
    relErrDist->SetMarkerColor(kBlack);
    relErrDist->SetMarkerStyle(20);
    relErrDist->SetXTitle("(raw-template)/raw");
    relErrDist->Draw();
    addControlPlot(c3);

    if(cpTmp) delete cpTmp;
    if(cpRawTmp) delete cpRawTmp;
}

void makeMuEff(const string filename="") 
{
	if (filename.length()<1)
	{
		cout << "Need input root file!" << endl;
		return;
	}

	vector< pair<double,double> >  vJetBins, vHtBins;
	vJetBins.push_back(make_pair(2,2));
	vJetBins.push_back(make_pair(3,5));
	vJetBins.push_back(make_pair(6,7));
	vJetBins.push_back(make_pair(8,1000));
	vHtBins.push_back(make_pair(500,8000));

	TFile f(filename.c_str());
	if (f.IsZombie())
	{
		cout << "File with name " << filename << " cannot be opened or does not exist!" << endl;
		return;
	}

	const string numname("muEff_num");
	const string denname("muEff_den");

	TCanvas *c = new TCanvas();
	size_t dotpos = filename.find(".");
	string substr = filename.substr(0,dotpos);
	stringstream epsname_o, epsname, epsname_c;
	epsname_o << substr << ".eps[";
	epsname << substr << ".eps";
	epsname_c << substr << ".eps]";
	c->Print(epsname_o.str().c_str());
	gStyle->SetOptStat(0);
	gStyle->SetPaintTextFormat("1.1f");
	for (unsigned jetbin=0; jetbin < vJetBins.size(); ++jetbin)
	{
		double minjet= vJetBins.at(jetbin).first;
		double maxjet= vJetBins.at(jetbin).second;
		for (unsigned htbin=0; htbin < vHtBins.size(); ++htbin)
		{
			double minht= vHtBins.at(htbin).first;
			double maxht= vHtBins.at(htbin).second;
			stringstream folder, numhistname, denhistname;
			folder << "muoneff/Hist/Njet" << minjet << "to" << maxjet << "Ht" << minht << "to" << maxht << "/";
			numhistname << folder.str() << numname;
			denhistname << folder.str() << denname;
			TH1* numhist = dynamic_cast<TH1*>( f.Get(numhistname.str().c_str()));
			TH1* denhist = dynamic_cast<TH1*> (f.Get(denhistname.str().c_str()));

			if (numhist == NULL) 
			{
				cout << "Numerator hist " << numhistname.str() << " is not found in file " << f.GetName() << endl;
				return;
			}

			if (denhist == NULL) 
			{
				cout << "Denominator hist " << denhistname.str() << " is not found in file " << f.GetName() << endl;
				return;
			}
			numhist->Draw("colz");
			c->Print(epsname.str().c_str());
			denhist->Draw("colz");
			c->Print(epsname.str().c_str());
			numhist->Sumw2();
			numhist->Divide(denhist);
			stringstream title;
			title << "Jets " << minjet << "-" << maxjet << ", HT>" << minht <<  ": Muon Reco+ID efficiency from " << substr << endl;
			numhist->SetTitle(title.str().c_str());
			//numhist->GetXaxis()->SetRangeUser(0,500);
			numhist->SetAxisRange(0,500,"X");
			numhist->SetAxisRange(0,3.5,"Y");
			//numhist->GetYaxis()->SetRangeUser(0,3.5);
			numhist->Draw("colzTEXT90E");
			c->Print(epsname.str().c_str());
		}
	}
	c->Print(epsname_c.str().c_str());

}