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

TF1* fit(TH1D* h, TH1D* hMCSignal, TH1D* hMCSwapped, Float_t ptmin, Float_t ptmax, Int_t j)
{
  TCanvas* c = new TCanvas(Form("c_%.0f_%.0f_%d",ptmin,ptmax,j),"",600,600);
  TF1* f = new TF1(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.14159)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.14159)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.14159)*[8]*(1+[11])))+[3]+[4]*x+[5]*x*x+[6]*x*x*x", 1.7, 2.0);
  
  f->SetParLimits(4,-1000,1000);
  f->SetParLimits(10,0.001,0.05);
  f->SetParLimits(2,0.01,0.5);
  f->SetParLimits(8,0.02,0.2);
  f->SetParLimits(7,0,1);
  f->SetParLimits(9,0,1);
  
  f->SetParameter(0,setparam0);
  f->SetParameter(1,setparam1);
  f->SetParameter(2,setparam2);
  f->SetParameter(10,setparam10);
  f->SetParameter(9,setparam9);

  f->FixParameter(8,setparam8);
  f->FixParameter(7,1);
  f->FixParameter(1,fixparam1);
  f->FixParameter(3,0);
  f->FixParameter(4,0);
  f->FixParameter(5,0);
  f->FixParameter(6,0);
  f->FixParameter(11,0);
  h->GetEntries();
  
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L m","",minhisto,maxhisto);
  
  f->FixParameter(1,f->GetParameter(1));
  f->FixParameter(2,f->GetParameter(2));
  f->FixParameter(10,f->GetParameter(10));
  f->FixParameter(9,f->GetParameter(9));
  f->FixParameter(7,0);
  f->ReleaseParameter(8);
  f->SetParameter(8,setparam8);
  
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L m","",minhisto,maxhisto);
  
  f->FixParameter(7,hMCSignal->Integral(0,1000)/(hMCSwapped->Integral(0,1000)+hMCSignal->Integral(0,1000)));
  f->FixParameter(8,f->GetParameter(8));
  f->ReleaseParameter(3);
  f->ReleaseParameter(4);
  f->ReleaseParameter(5);
  f->ReleaseParameter(6);

  f->SetLineColor(kRed);
  
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  f->SetParLimits(1,1.86,1.87);
  /*
  f->ReleaseParameter(11);
  f->SetParLimits(11,-1.,1.);
  */
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L m","",minhisto,maxhisto);
  
  TF1* background = new TF1(Form("background_%.0f_%.0f_%d",ptmin,ptmax,j),"[0]+[1]*x+[2]*x*x+[3]*x*x*x");
  background->SetParameter(0,f->GetParameter(3));
  background->SetParameter(1,f->GetParameter(4));
  background->SetParameter(2,f->GetParameter(5));
  background->SetParameter(3,f->GetParameter(6));
  background->SetLineColor(4);
  background->SetRange(minhisto,maxhisto);
  background->SetLineStyle(2);
  
  TF1* mass = new TF1(Form("mass_%.0f_%.0f_%d",ptmin,ptmax,j),"[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.14159)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.14159)*[5]*(1+[6]))))");
  mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10),f->GetParameter(11));
  mass->SetParError(0,f->GetParError(0));
  mass->SetParError(1,f->GetParError(1));
  mass->SetParError(2,f->GetParError(2));
  mass->SetParError(3,f->GetParError(7));
  mass->SetParError(4,f->GetParError(9));
  mass->SetParError(5,f->GetParError(10));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);
  
  TF1* massSwap = new TF1(Form("massSwap_%.0f_%.0f_%d",ptmin,ptmax,j),"[0]*(1-[2])*Gaus(x,[1],[3]*(1+[4]))/(sqrt(2*3.14159)*[3]*(1+[4]))");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8),f->GetParameter(11));
  massSwap->SetParError(0,f->GetParError(0));
  massSwap->SetParError(1,f->GetParError(1));
  massSwap->SetParError(2,f->GetParError(7));
  massSwap->SetParError(3,f->GetParError(8));
  massSwap->SetFillColor(kGreen+4);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          //hlumi0->SetBinContent(i, lum0);
          //}

          outdir->cd();
          if (h) hpt_notimedep = (TH1D*)h->Clone("hpt_notimedep");
          if (hpre && h) htimedep = (TH1D*)hpre->Clone("htimedep");
          if (hpre && h) htimedep->Divide(hpre,h);//,1,1,"B");

          // Figure out trigger luminosities
          double lumi = 0;
          if (hlumi) lumi = hlumi->GetBinContent(1);
          double lumi0 = 0;
          if (hlumi0) lumi0 = hlumi0->GetBinContent(1);
          if (htimedep && lumi && lumi0) {
            htimedep->Scale(lumi / lumi0);
          }

          // Find proper pT range and fit
          double minpt = 0.;
          double maxpt = 6500.;
          if (hsel) {
            for (int i = 1; i != hsel->GetNbinsX()+1; ++i) {
              if (hsel->GetBinContent(i)!=0 &&
                  hsel->GetBinLowEdge(i)>=_jp_xmin57) {
                if (minpt<20) minpt = hsel->GetBinLowEdge(i);
                maxpt = hsel->GetBinLowEdge(i+1);
              }
            }
          }
          TF1 *ftmp = new TF1("ftmp","[0]",minpt,maxpt);
          ftmp->SetParameter(0,1);
          if (htimedep && htimedep->Integral()>0) htimedep->Fit(ftmp,"QRN");

          if (htimedep && ftmp->GetParameter(0)>0)
            ktime = 1./ftmp->GetParameter(0);

          if (htimedep) {
            outdir->cd();
            htimefit = (TH1D*)hsel->Clone("htimefit");
            hpt_withtimedep = (TH1D*)h->Clone("hpt_withtimedep");

            for (int i = 1; i != htimefit->GetNbinsX()+1; ++i) {

              if (hsel->GetBinContent(i)!=0) {
                htimefit->SetBinContent(i, ftmp->GetParameter(0));
                htimefit->SetBinError(i, ftmp->GetParError(0));
              }

              // Calculate with time dependence here to add ktime fit error
              hpt_withtimedep->SetBinContent(i, hpt_notimedep->GetBinContent(i)
                                             * htimefit->GetBinContent(i));
              double err1 = hpt_notimedep->GetBinError(i)
                / hpt_notimedep->GetBinContent(i);
              double err2 = htimefit->GetBinError(i)
                / htimefit->GetBinContent(i);
              hpt_withtimedep->SetBinError(i, hpt_notimedep->GetBinContent(i)
                                           * sqrt(pow(err1,2) + pow(err2,2)));
            }
          }
        } // dotimedep


        if (!(hpt->GetNbinsX()==peff->GetNbinsX() || isoth || isgen) ||
            !(hpt->GetNbinsX()==hlumi->GetNbinsX() || isoth || isgen)) {
          cerr << "Hist " << hpt->GetName() << " " << dir->GetName()
               << " Nbins=" << hpt->GetNbinsX() << endl << flush;

//
// Analyze_ECvsP - fit the EC/P vs P distributions
//
//                  fAna = output from eg2a DMS
//                  target = target name
//
void Analyze_ECvsP(char *fAna="Ana.root", char *target)
{
    Int_t i, j;
    
    const Int_t NPARAM = 3;
    TH1D *h1D[NSECTORS][NPARAM];
    TH2D *h2D[NSECTORS];
    
    char HistName2D[50];
    char strname[50];
    char *yname[NPARAM] = {"Amplitude","Mean","Sigma"};
    
    // open text file for the yields
    char OutFile[100];
    sprintf(OutFile,"ECvsP_Fit_%s.dat",target);
    ofstream fout(OutFile);
    
    TCanvas *can[NSECTORS]; // Canvas to plot histogram
    
    Double_t par[3]={1.0,1.0,1.0};
    TF1 *pol;
    TF1 *sig;
    
    // data files contain the trees
    printf("Analyzing file %s\n",fAna);
    TFile *fm = new TFile(fAna,"READ");
    TDirectory *tmp = fm->GetDirectory("ElectronID");
    
    for(j=0; j<NSECTORS; j++){
        sprintf(cname,"can%i",j+1);
        sprintf(cname,"Canvas, Sector %i",j+1);
        can[j] = new TCanvas(cname,ctitle,50*j,0,600,600);
        can[j]->SetBorderMode(1);  //Bordermode (-1=down, 0 = no border, 1=up)
        can[j]->SetBorderSize(5);
        gStyle->SetOptStat(1111);
        can[j]->SetFillStyle(4000);
        can[j]->Divide(2,2);
        
        sprintf(HistName2D,"ECtotP_VS_P_Sector%i",j+1);
        h2D[j] = (TH2D*)tmp->Get(HistName2D);

        can[j]->cd(1);
        gPad->SetLeftMargin(Lmar);
        gPad->SetRightMargin(Rmar);
        gPad->SetFillColor(0);

        h2D[j]->GetXaxis()->CenterTitle();
        h2D[j]->GetYaxis()->CenterTitle();
        h2D[j]->GetYaxis()->SetTitleOffset(yoff);
    
        h2D[j]->SetAxisRange(0.0,3.0,"X");
        h2D[j]->Draw("colz");
    
        h2D[j]->FitSlicesY(0,0,-1,20);
    
        fout<<j+1<<"\t";
        
        for(i=0; i<NPARAM; i++){
            can[j]->cd(i+2);
            gPad->SetLeftMargin(Lmar);
            gPad->SetRightMargin(Rmar);
            gPad->SetFillColor(0);

            sprintf(strname,"%s_%i",HistName2D,i);
            h1D[j][i] = (TH1D*)gDirectory->Get(strname);
            h1D[j][i]->GetXaxis()->CenterTitle();
            h1D[j][i]->GetYaxis()->CenterTitle();
            h1D[j][i]->GetYaxis()->SetTitle(yname[i]);
            h1D[j][i]->GetYaxis()->SetTitleOffset(yoff);
            h1D[j][i]->SetAxisRange(0.0,2.5,"X");
            if (i==0){
                h1D[j][i]->SetAxisRange(0.0,40.0,"Y");
            }
            if (i==1) {
                h1D[j][i]->SetAxisRange(0.15,0.35,"Y");
                pol = new TF1("pol","pol2",0.5,3.0);
                h1D[j][i]->Fit("pol","R");
                fout<<pol->GetParameter(0)<<"\t"<<pol->GetParameter(1)<<"\t"<<pol->GetParameter(2)<<"\t";
            }
            if (i==2) {
                par[0] = 10.0;
                par[1] = 0.1;
                par[2] = 0.1;
                h1D[j][i]->SetAxisRange(0.0,0.3,"Y");
                sig = new TF1("sig",SigmaFit,0.5,3.0,2);
                sig->SetParameters(par);
                h1D[j][i]->Fit("sig","R");
                fout<<sig->GetParameter(0)<<"\t"<<sig->GetParameter(1)<<endl;
            }
            h1D[j][i]->Draw();
        }
       	sprintf(OutCan,"Analyze_ECvsP_S%i_%s.gif",j+1,target);
        can[j]->Print(OutCan);
        sprintf(OutCan,"Analyze_ECvsP_S%i_%s.eps",j+1,target);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  float sideBandYield = 0;
  float mixedSideBandYield = 0;
  for (int i = 0; i < massNBins; ++i)
  {
    const float binCenter = sig->GetXaxis()->GetBinCenter(i+1);
    if(binCenter < sideBandTo || binCenter > sideBandFrom)
    {
      sideBandYield += sig->GetBinContent(i+1);
      mixedSideBandYield += bkgMEUS->GetBinContent(i+1);
    }

  }
  TH1D* sideBandMarker = new TH1D("sideBandMarker","marker for the sideband in the form of a histogram", massNBins, massFrom, massTo);
  const float binWidth = (massTo - massFrom) / static_cast<float>(massNBins);
  for(int i = 0; i < massNBins; ++i)
    sideBandMarker->SetBinContent(i+1, 1000);
  for(float bin = sideBandTo; bin < sideBandFrom; bin += binWidth)
  {
    // cout << bin << endl;
    sideBandMarker->Fill(bin + 0.5*binWidth, -1e6);
  }
  sideBandMarker->SetFillColor(kRed);
  sideBandMarker->SetFillStyle(3004);
  sideBandMarker->SetLineColorAlpha(kRed,0);

  bkgSELS->Scale(1./3.);
  cout << "Scaling mixed event bg by " << sideBandYield / mixedSideBandYield << endl;
  bkgMEUS->Scale(sideBandYield / (mixedSideBandYield) );

  TF1 *line = new TF1("line", "[0] + [1]*x");
  TF1 *gaussPlusLine = new TF1("gaussPlusLine", "[0] + [1]*x + [2]*TMath::Gaus(x,[3],[4],1)");

  bkgMEUS->Fit(line, "", "", 2.1, 2.5);
  double offset = line->GetParameter(0);
  double slope = line->GetParameter(1);
  const double LcMass = 2.28646;

  gaussPlusLine->FixParameter(0,offset);
  gaussPlusLine->FixParameter(1,slope);
  gaussPlusLine->SetParameter(3,LcMass);

  gaussPlusLine->SetParameter(2,1);
  gaussPlusLine->SetParameter(4,sigma);

  sig->Fit(gaussPlusLine, "", "", 2.1, 2.5);

  // -- Drawing
  sig->GetXaxis()->SetTitle("#font[12]{m}_{pK#pi} (GeV/#font[12]{c}^{2})");
  sig->GetXaxis()->CenterTitle();
  sig->GetXaxis()->SetLabelSize(0.04);
  sig->GetXaxis()->SetTitleSize(0.055);
  sig->GetXaxis()->SetTitleOffset(0.8);
  // sig->GetXaxis()->SetTickLength(0.001);
  sig->GetYaxis()->SetTitle("#font[12]{N}/(10 MeV/#font[12]{c}^{2})");
  sig->GetYaxis()->CenterTitle();
  sig->GetYaxis()->SetLabelSize(0.04);
  sig->GetYaxis()->SetTitleSize(0.055);
  sig->GetYaxis()->SetTitleOffset(0.8);
  // sig->GetYaxis()->SetRangeUser(17, 95);

  sig->SetMarkerStyle(20);
  sig->SetLineColor(kBlack);

  bkgMEUS->SetMarkerStyle(21);
  bkgMEUS->SetMarkerColor(kOrange+5);
  bkgMEUS->SetLineColor(kOrange+5);

//.........这里部分代码省略.........
		//}
		//else{
			gfit->SetParameters((Double_t)N,m,s);
			gfit->SetRange(m-2*s,m+1*s); //fit within 2 std devs
			//if(m > p) gfit->SetRange(p-2*s,p+1*s); //high tail
			//else gfit->SetRange(p-1*s,p+2*s); //low tail
		//}
		
		//formatting
		gfit->SetLineColor(kRed);
		gfit->SetMarkerColor(kRed);
		gfit->SetLineWidth(2);
		//fit
		h_res->Fit(gfit,"LNQR");
	}
	
	//store parameters
	theRes->setStats(stats,stat_err);
	if(do_fit) theRes->setFit(gfit);
	//store histo
	h_res->SetDirectory(0);
	theRes->setHist(h_res);
	
	std::stringstream muname, signame, musigname, aLname, nLname, aRname, nRname, Nname, chiname;
	muname.precision(2);
	signame.precision(2);
	musigname.precision(3);
	aLname.precision(2);
	nLname.precision(2);
	aRname.precision(2);
	nRname.precision(2);
	chiname.precision(5);
	if (do_fit) {
		muname << fixed << "#mu = " << gfit->GetParameter(1) << " #pm " << gfit->GetParError(1);
		signame << fixed << "#sigma = " << gfit->GetParameter(2) << " #pm " << gfit->GetParError(2);
		musigname << fixed << "#sigma/#mu = " << gfit->GetParameter(2)/gfit->GetParameter(1) << " #pm " << 
		gfit->GetParameter(2)/gfit->GetParameter(1) * sqrt( Power(gfit->GetParError(1),2)/Power(gfit->GetParameter(1),2) +  Power(gfit->GetParError(2),2)/Power(gfit->GetParameter(2),2) );
		//aLname << fixed << "a_{L} = " << gfit->GetParameter(3) << " #pm " << gfit->GetParError(3);
		//nLname << fixed << "n_{L} = " << gfit->GetParameter(4) << " #pm " << gfit->GetParError(4);
		//aRname << fixed << "a_{R} = " << gfit->GetParameter(5) << " #pm " << gfit->GetParError(5);
		//nRname << fixed << "n_{R} = " << gfit->GetParameter(6) << " #pm " << gfit->GetParError(6);
		chiname << fixed << "#chi^{2}/ndf = " << gfit->GetChisquare()/gfit->GetNDF();
	}
	else {
		muname << fixed << "Mean = " << m << " #pm " << me;
		signame << fixed << "RMS = " << s << " #pm " << se;
		musigname << fixed << "RMS/Mean = " << s/m << " #pm " << s/m*sqrt((me*me)/(m*m)+(se*se)/(s*s));
	}
	Nname << "N = " << N; 

	//plotting
	if (do_show){
		can = new TCanvas((oname.str()).c_str(),(oname.str()).c_str(),700,500);
		can->cd();
		pad = new TPad("graph","",0,0,1,1);
		pad->SetMargin(0.12,0.05,0.15,0.05);
		pad->Draw();
		pad->cd();
		
		//formatting
		h_res->SetStats(kTRUE);
		gStyle->SetOptStat("mr");
		h_res->GetYaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		h_res->GetYaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));
		h_res->GetXaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		h_res->GetXaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));

//.........这里部分代码省略.........
  } //for entries

  TCanvas* c_prova = new TCanvas("c_prova", "c_prova", 800, 600);
  c_prova->cd();
  


  Double_t Lower[nBins];

  Lower[0] = 80.;  
  Lower[1] = 100.;  
  Lower[2] = 125.;  
  Lower[3] = 160.;  
  Lower[4] = 190.;  
  Lower[5] = 215.;  


  TH1F* h1_resolution_vs_pt = new TH1F("resolution_vs_pt", "Jet Energy Resolution vs. Average Jet p_{T}", nBins, Lower);
  h1_resolution_vs_pt->SetXTitle("Average p_{T}^{RECO} [GeV/c]");
  h1_resolution_vs_pt->SetYTitle("#sigma(p_{T})/p_{T}");


  for(int iBin=0; iBin<nBins; ++iBin) {

    int nPoints = 6;
    Float_t x[nPoints], y[nPoints];
    Float_t err_x[nPoints], err_y[nPoints];

    x[0] = 10.;
    x[1] = 12.;
    x[2] = 15.;
    x[3] = 17.;
    x[4] = 20.;
    x[5] = 30.;

    y[0] = sqrt(2)*h1_asymmetry_less10[iBin]->GetRMS();
    y[1] = sqrt(2)*h1_asymmetry_less12[iBin]->GetRMS();
    y[2] = sqrt(2)*h1_asymmetry_less15[iBin]->GetRMS();
    y[3] = sqrt(2)*h1_asymmetry_less17[iBin]->GetRMS();
    y[4] = sqrt(2)*h1_asymmetry_less20[iBin]->GetRMS();
    y[5] = sqrt(2)*h1_asymmetry_less30[iBin]->GetRMS();

    err_x[0] = 0.; 
    err_x[1] = 0.; 
    err_x[2] = 0.; 
    err_x[3] = 0.; 
    err_x[4] = 0.; 
    err_x[5] = 0.; 

    err_y[0] = y[0]/sqrt(h1_asymmetry_less10[iBin]->GetEntries()); 
    err_y[1] = y[1]/sqrt(h1_asymmetry_less12[iBin]->GetEntries()); 
    err_y[2] = y[2]/sqrt(h1_asymmetry_less15[iBin]->GetEntries()); 
    err_y[3] = y[3]/sqrt(h1_asymmetry_less17[iBin]->GetEntries()); 
    err_y[4] = y[4]/sqrt(h1_asymmetry_less20[iBin]->GetEntries()); 
    err_y[5] = y[5]/sqrt(h1_asymmetry_less30[iBin]->GetEntries()); 


    TGraphErrors* resolution_vs_pt3 = new TGraphErrors(nPoints, x, y, err_x, err_y);
    resolution_vs_pt3->SetMarkerStyle(20);

    TF1* line = new TF1("line", "[0] + [1]*x");
    line->SetParameter(0, y[0]);
    line->SetParameter(1, (y[1]-y[0])/(x[1]-x[0]));
    line->SetRange(0., 32.);

    resolution_vs_pt3->Fit(line, "R");

    char plotName[70];
    sprintf(plotName, "asymmetry/asymmetry_less10_%d.eps", iBin);
    h1_asymmetry_less10[iBin]->Draw();
    c_prova->SaveAs(plotName);

    sprintf(plotName, "asymmetry/asymmetry_less20_%d.eps", iBin);
    h1_asymmetry_less20[iBin]->Draw();
    c_prova->SaveAs(plotName);

    sprintf(plotName, "asymmetry/asymmetry_less30_%d.eps", iBin);
    h1_asymmetry_less30[iBin]->Draw();
    c_prova->SaveAs(plotName);


    sprintf(plotName, "asymmetry/resolution_vs_pt3_%d.eps", iBin);
    resolution_vs_pt3->Draw("AP");
    c_prova->SaveAs(plotName);


    h1_resolution_vs_pt->SetBinContent(iBin+1, line->GetParameter(0));
    h1_resolution_vs_pt->SetBinError(iBin+1, line->GetParError(0));

  } //for pt bins


  TFile* outFile = new TFile("asymmetryOutputFile.root", "RECREATE");

  outFile->cd();
  h1_resolution_vs_pt->Write();
  outFile->Close();
  outFile->Write();

} //asymmetry

//.........这里部分代码省略.........
	  vector<Double_t> myt0_err;
	  vector<Double_t> x_t0;
	  vector<Double_t> x_err_t0;

	  int count_valid_points(0);

	  for(Int_t myttrig = 0; myttrig < 5; myttrig++){

	    if(Ares_mean_err[mywheel+2][mysector-1][mystation-1][mySL-1][myttrig] > 0.2) continue;

	    mymean.push_back(Ares_mean[mywheel+2][mysector-1][mystation-1][mySL-1][myttrig]);
	    mymean_err.push_back(Ares_mean_err[mywheel+2][mysector-1][mystation-1][mySL-1][myttrig]);
	    x_err.push_back(0.);
	    x_mean.push_back(deltaTtrig[myttrig]);

	    myt0.push_back(t0res_mean[mywheel+2][mysector-1][mystation-1][mySL-1][myttrig]);
	    myt0_err.push_back(0.5);
	    x_err_t0.push_back(0.);
	    x_t0.push_back(deltaTtrig[myttrig]);

	    count_valid_points++;
	  }

	  //cout << " res mean " << t0res_mean[mywheel+2][mysector-1][mystation-1][mySL-1][2] << endl;

	  if(count_valid_points < 2){
	    fout << mywheel << " " << mysector << " " << mystation << " " << mySL << endl;
	    continue;
	  }

	  TGraphErrors mygr(mymean.size(),&x_mean[0],&mymean[0],&x_err[0],&mymean_err[0]);
	  mygr.Fit("pol1");
	  TF1 * mygrFit = mygr.GetFunction("pol1");
	  double par0res = mygrFit->GetParameter("p0");
	  double par1res = mygrFit->GetParameter("p1");
	  double ttrigCorrRes  = 99999;
	  if(par1res != 0) ttrigCorrRes = -par0res/par1res;

	  ttrigCorr[mywheel+2][mysector-1][mystation-1][mySL-1] = par1res;

          string stringName1;
          stringstream theStream1;
          theStream1 <<  "Wh" << mywheel << "_Sect" << mysector << "_St"<< mystation << "_SL" << mySL;
          theStream1 >> stringName1;

	  mygr.SetTitle(stringName1.c_str());

	  c1.cd();
	  mygr.Draw("AP");

	  if(count == 0) c1.SaveAs("resplots.ps(");
	  else c1.SaveAs("resplots.ps");

	  TGraphErrors mygrt0(myt0.size(),&x_t0[0],&myt0[0],&x_err_t0[0],&myt0_err[0]);
	  mygrt0.Fit("pol1");

	  TF1 * mygrt0Fit = mygrt0.GetFunction("pol1");
	  double par0t0 = mygrt0Fit->GetParameter("p0");
	  double par1t0 = mygrt0Fit->GetParameter("p1");
	  double ttrigCorrT0    = 99999;
	  if(par1t0 != 0) ttrigCorrT0 = -par0t0/par1t0;

	  if(mySL != 2) {
	    hCorr->Fill(ttrigCorrRes, ttrigCorrT0);
	    hPar1_phi->Fill(-par1res);
	    hPar1VsSect_phi->Fill(mysector,-par1res);

void SimpleFitter()
{
   
   TFile*output=new TFile("histoX.root","read");    
   TH1D*h=(TH1D*)output->Get("h");

   h->Sumw2();
   TF1 *f = new TF1("f","[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2]) +[3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5])+[6]+[7]*x+[8]*x*x+[9]*x*x*x");
   f->SetLineColor(4);
   f->SetParameters(21,3.68,0.0057,1.34714e+00,3.86,0.006,-2.21334e+04,1.01661e+04,-1.05364e+03,1.95474e+04);
   f->FixParameter(1,3.686);
   f->FixParameter(2,0.00357);
   f->FixParameter(4,3.8725);
   f->FixParameter(5,0.005);

   h->Fit("f","LL");
   h->Fit("f","");
   h->Fit("f","LL");
   
   h->Fit("f","LL","",3.65,3.94);
   h->Fit("f","LL","",3.65,3.94);

   f->ReleaseParameter(1);
   f->ReleaseParameter(2);
   f->ReleaseParameter(4);
   f->ReleaseParameter(5);
   h->Fit("f","LL","",3.65,3.94);
   h->SetXTitle("m(J/#psi#pi^{+}#pi^{-}) [GeV]");
   h->SetYTitle("Entries");
   h->SetStats(0);

   h->SetAxisRange(0,h->GetMaximum()*1.3	,"Y");
   TF1 *f2 = new TF1("f2","[6]+[7]*x+[8]*x*x+[9]*x*x*x");
   f2->SetParameter(6,f->GetParameter(6));
   f2->SetParameter(7,f->GetParameter(7));
   f2->SetParameter(8,f->GetParameter(8));
   f2->SetParameter(9,f->GetParameter(9));
   TF1 *f3 = new TF1("f3","[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])");
   f3->SetParameter(0,f->GetParameter(0));
   f3->SetParameter(1,f->GetParameter(1));
   f3->SetParameter(2,f->GetParameter(2));

   TF1 *f4 = new TF1("f4","[3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5])");
   f4->SetParameter(3,f->GetParameter(3));
   f4->SetParameter(4,f->GetParameter(4));
   f4->SetParameter(5,f->GetParameter(5));
      
   f->SetLineColor(4);
   f2->SetLineColor(4);
   f3->SetRange(3.65,3.94);
   f4->SetRange(3.65,3.94);
   f2->SetRange(3.65,3.94);
   f2->SetLineStyle(2);
   f3->SetLineStyle(2);
   f4->SetLineStyle(2);
   f2->Draw("same");
   f3->SetLineColor(2);
   f3->SetFillStyle(3004);
   f3->SetFillColor(2);
   f3->Draw("same");
   f4->SetLineColor(2);
   f4->SetFillStyle(3004);
   f4->SetFillColor(2);
   f4->Draw("same");

   TLatex *l = new TLatex(3.7,70./80*h->GetMaximum(),"#psi(2S)");
   l->Draw();
   TLatex *l2 = new TLatex(3.875,50./80*h->GetMaximum(),"X(3872)");
   l2->Draw();
   TLatex *l3 = new TLatex(3.812,70./80*h->GetMaximum(),"CMS Preliminary");
   l3->Draw();
   TLatex *l4 = new TLatex(3.78,60./80*h->GetMaximum(),"pp #sqrt{s_{NN}}=5.02 TeV");
   l4->Draw();
   
   cout<<"********************* SIGNIFICANCE *********************"<<endl;
   
   double meanX=f4->GetParameter(4);
   double sigmaX=f4->GetParameter(5);
   double signalX=f4->Integral(meanX-3*sigmaX,meanX+3*sigmaX)/h->GetBinWidth(0);
   double bkgX=f2->Integral(meanX-3*sigmaX,meanX+3*sigmaX)/h->GetBinWidth(0);
   double signalpsi2=f4->Integral(meanX-3*sigmaX,meanX+3*sigmaX)/h->GetBinWidth(0);
   double bkgpsi2=f2->Integral(meanX-3*sigmaX,meanX+3*sigmaX)/h->GetBinWidth(0);
   double significanceX=signalX/TMath::Sqrt(signalX+bkgX);   
   
   double meanpsi2=f3->GetParameter(1);
   double sigmapsi2=f3->GetParameter(2);
   double signalpsi2=f3->Integral(meanpsi2-3*sigmapsi2,meanpsi2+3*sigmapsi2)/h->GetBinWidth(0);
   double bkgpsi2=f2->Integral(meanpsi2-3*sigmapsi2,meanpsi2+3*sigmapsi2)/h->GetBinWidth(0);
   double signalpsi2=f3->Integral(meanpsi2-3*sigmapsi2,meanpsi2+3*sigmapsi2)/h->GetBinWidth(0);
   double bkgpsi2=f2->Integral(meanpsi2-3*sigmapsi2,meanpsi2+3*sigmapsi2)/h->GetBinWidth(0);
   double significancepsi2=signalpsi2/TMath::Sqrt(signalpsi2+bkgpsi2);


   cout<<"significance X ="<<significanceX<<endl;
   cout<<"significance psi(2S)="<<significancepsi2<<endl;
   
   cout<<"yield X="<<signalX<<endl;
   cout<<"yield 2S="<<signalpsi2<<endl;
   cout<<"yield X/yield 2S="<<signalX/signalpsi2<<endl;

//.........这里部分代码省略.........

int piRandom() {

  TRandom                  r0;
  TRandom1                  r1;
  TRandom2                 r2;
  TRandom3                 r3;
  //Random<GSLRngRand>       r1;
//   Random<GSLRngTaus>       r2;
//   Random<GSLRngRanLux>     r3;

  double n = NEVT;
  int nloop = NLOOP;
  double dy = 15/std::sqrt(n);

  TH1D * h0 = new TH1D("h0","TRandom delta",100,-dy,dy);
  TH1D * h1 = new TH1D("h1","TRandom1 delta",100,-dy,dy);
  TH1D * h2 = new TH1D("h2","TRandom2 delta",100,-dy,dy);
  TH1D * h3 = new TH1D("h3","TRandom3 delta",100,-dy,dy);

  double sigma = std::sqrt( PI * (4 - PI)/n );
  std::cout << "**********************************************************" << std::endl;
  std::cout << " Generate " << int(n) << " for " << nloop      << " times " << std::endl;
  std::cout << "**********************************************************" << std::endl;
  std::cout << "\tExpected Sigma = " << sigma << std::endl;

#define INTERACTIVE
#ifdef INTERACTIVE

  double del, err;
  TCanvas * c1 = new TCanvas("c1_piRandom","PI Residuals");
  gStyle->SetOptFit(1111);
  gStyle->SetOptLogy();
  c1->Divide(2,2);
  c1->cd(1);
  generate(r0,h0);
  h0->Fit("gaus");
  h0->Draw();
  TF1 * fg = (TF1*) h0->FindObject("gaus");
  if (fg) {
    del = (fg->GetParameter(2)-sigma);
    err = fg->GetParError(2);
  }
  else { del = -999; err = 1; }

  char text[20];
  sprintf(text,"Spread %8.4f",del/err);
  TPaveLabel * pl0 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
  pl0->Draw();




  c1->cd(2);
  generate(r1,h1);
  h1->Fit("gaus");
  h1->Draw();
  fg = (TF1*) h1->FindObject("gaus");
  if (fg) {
    del = (fg->GetParameter(2)-sigma);
    err = fg->GetParError(2);
  }
  else { del = -999; err = 1; }

  sprintf(text,"Spread %8.4f",del/err);
  TPaveLabel * pl1 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
  pl1->Draw();


  c1->cd(3);
  generate(r2,h2);
  h2->Fit("gaus");
  h2->Draw();
  fg = (TF1*) h2->FindObject("gaus");
  if (fg) {
    del = (fg->GetParameter(2)-sigma);
    err = fg->GetParError(2);
  }
  else { del = -999; err = 1; }

  sprintf(text,"Spread %8.4f",del/err);
  TPaveLabel * pl2 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
  pl2->Draw();


  c1->cd(4);
  generate(r3,h3);
  h3->Fit("gaus");
  h3->Draw();
  fg = (TF1*) h3->FindObject("gaus");
  if (fg) {
    del = (fg->GetParameter(2)-sigma);
    err = fg->GetParError(2);
  }
  else { del = -999; err = 1; }

  sprintf(text,"Spread %8.4f",del/err);
  TPaveLabel * pl3 = new TPaveLabel(0.6,0.3,0.9,0.4,text,"brNDC");
  pl3->Draw();

#else
//.........这里部分代码省略.........

void ABCD2vari_p3(const char* ISO) { //Only iso is external, all the regions go in a loop!!! to be fed in the macro at once

  
  TString PREFIX = "./";

  //Get the fractions
  //gStyle->SetOptStat(0);
  TFile* fB = new TFile(PREFIX+"factor_qcd_"+TString(ISO)+"_B.root", "open");
  fB->cd();
  TH1D* hRatioB = (TH1D*)gDirectory->Get("hfactor");
  hRatioB->Fit("pol0");
  //2D needs, preliminary
  TH1D* hRatioB2030 = (TH1D*)gDirectory->Get("hfactor2030");
  hRatioB2030->Fit("pol0");
  TH1D* hRatioB3045 = (TH1D*)gDirectory->Get("hfactor3045");
  hRatioB3045->Fit("pol0");
  TH1D* hRatioB4560 = (TH1D*)gDirectory->Get("hfactor4560");
  hRatioB4560->Fit("pol0");
  TH1D* hRatioB60120 = (TH1D*)gDirectory->Get("hfactor60120");
  hRatioB60120->Fit("pol0");
  TH1D* hRatioB120200 = (TH1D*)gDirectory->Get("hfactor120200");
  hRatioB120200->Fit("pol0");
  TH1D* hRatioB2001500 = (TH1D*)gDirectory->Get("hfactor2001500");
  hRatioB2001500->Fit("pol0");

  TF1 *myfitB = (TF1*)hRatioB->GetFunction("pol0");
  Double_t factor_B = myfitB->GetParameter(0);
  Double_t factor_B_err = myfitB->GetParError(0);
  //2D needs, preliminary
  TF1 *myfitB2030 = (TF1*)hRatioB2030->GetFunction("pol0");
  Double_t factor_B2030 = myfitB2030->GetParameter(0);
  Double_t factor_B_err2030 = myfitB2030->GetParError(0);
  TF1 *myfitB3045 = (TF1*)hRatioB3045->GetFunction("pol0");
  Double_t factor_B3045 = myfitB3045->GetParameter(0);
  Double_t factor_B_err3045 = myfitB3045->GetParError(0);
  TF1 *myfitB4560 = (TF1*)hRatioB4560->GetFunction("pol0");
  Double_t factor_B4560 = myfitB4560->GetParameter(0);
  Double_t factor_B_err4560 = myfitB4560->GetParError(0);
  TF1 *myfitB60120 = (TF1*)hRatioB60120->GetFunction("pol0");
  Double_t factor_B60120 = myfitB60120->GetParameter(0);
  Double_t factor_B_err60120 = myfitB60120->GetParError(0);
  TF1 *myfitB120200 = (TF1*)hRatioB120200->GetFunction("pol0");
  Double_t factor_B120200 = myfitB120200->GetParameter(0);
  Double_t factor_B_err120200 = myfitB120200->GetParError(0);
  TF1 *myfitB2001500 = (TF1*)hRatioB2001500->GetFunction("pol0");
  Double_t factor_B2001500 = myfitB2001500->GetParameter(0);
  Double_t factor_B_err2001500 = myfitB2001500->GetParError(0);


  //repeat for region D
  TFile* fD = new TFile(PREFIX+"factor_qcd_"+TString(ISO)+"_D.root", "open");
  fD->cd();
  TH1D* hRatioD = (TH1D*)gDirectory->Get("hfactor");
  hRatioD->Fit("pol0");

  //2D needs, preliminary
  TH1D* hRatioD2030 = (TH1D*)gDirectory->Get("hfactor2030");
  hRatioD2030->Fit("pol0");
  TH1D* hRatioD3045 = (TH1D*)gDirectory->Get("hfactor3045");
  hRatioD3045->Fit("pol0");
  TH1D* hRatioD4560 = (TH1D*)gDirectory->Get("hfactor4560");
  hRatioD4560->Fit("pol0");
  TH1D* hRatioD60120 = (TH1D*)gDirectory->Get("hfactor60120");
  hRatioD60120->Fit("pol0");
  TH1D* hRatioD120200 = (TH1D*)gDirectory->Get("hfactor120200");
  hRatioD120200->Fit("pol0");
  TH1D* hRatioD2001500 = (TH1D*)gDirectory->Get("hfactor2001500");
  hRatioD2001500->Fit("pol0");

  TF1 *myfitD = (TF1*)hRatioD->GetFunction("pol0");
  Double_t factor_D = myfitD->GetParameter(0);
  Double_t factor_D_err = myfitD->GetParError(0);
  //2D needs, preliminary
  TF1 *myfitD2030 = (TF1*)hRatioD2030->GetFunction("pol0");
  Double_t factor_D2030 = myfitD2030->GetParameter(0);
  Double_t factor_D_err2030 = myfitD2030->GetParError(0);
  TF1 *myfitD3045 = (TF1*)hRatioD3045->GetFunction("pol0");
  Double_t factor_D3045 = myfitD3045->GetParameter(0);
  Double_t factor_D_err3045 = myfitB3045->GetParError(0);
  TF1 *myfitD4560 = (TF1*)hRatioD4560->GetFunction("pol0");
  Double_t factor_D4560 = myfitD4560->GetParameter(0);
  Double_t factor_D_err4560 = myfitD4560->GetParError(0);
  TF1 *myfitD60120 = (TF1*)hRatioD60120->GetFunction("pol0");
  Double_t factor_D60120 = myfitD60120->GetParameter(0);
  Double_t factor_D_err60120 = myfitD60120->GetParError(0);
  TF1 *myfitD120200 = (TF1*)hRatioD120200->GetFunction("pol0");
  Double_t factor_D120200 = myfitD120200->GetParameter(0);
  Double_t factor_D_err120200 = myfitD120200->GetParError(0);
  TF1 *myfitD2001500 = (TF1*)hRatioD2001500->GetFunction("pol0");
  Double_t factor_D2001500 = myfitD2001500->GetParameter(0);
  Double_t factor_D_err2001500 = myfitD2001500->GetParError(0);

  // directory with data
  TString protocol = "file://";
  TString dirname = "/scratch/lustreC/a/asvyatko/DY2013/rootfiles/";

  // EWK
  TFileCollection* c6 = new TFileCollection("WJets","WJets");
  c6->Add(protocol+dirname+"WJets_PU"+"/*.root");

//.........这里部分代码省略.........

//.........这里部分代码省略.........
//      for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}
      TH1D * hpm2   = (TH1D*)hpm->Clone(Form("Bg1_%d",iPID)) ;
      TH1D * hpcopy = (TH1D*)hp ->Clone(Form("hpcopy_%d",iPID)) ;
      TH1D * hp2    = (TH1D*)hp ->Clone(Form("hp2_%d",iPID)) ;
      hpcopy->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
      hp2   ->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
      hpcopy->Divide(hpm) ;
      hpcopy->SetTitle(Form("%3.1f<p_{T}<%3.1f GeV/c",xa[i-1],xa[i])) ;
      hpcopy->SetMarkerStyle(20) ;
      hpcopy->SetMarkerSize(0.7) ;
      hpcopy->GetXaxis()->SetRangeUser(0.05,0.25) ;
      
//      fit1->SetParameters(0.0002+0.0001*i*i,0.136,0.011,0.0002,-0.002,0.0) ;
      fit1->SetParameters(0.0002+0.0001*i*i,0.136,0.005,0.0002,-0.002,0.0) ;
//      fit1->SetParameters(0.0002,0.136,0.011,0.0002,-0.002,0.0) ;
      if(cen==0)
        fit1->SetParameters(0.001,0.146,0.005,0.12,-0.002,0.0) ;
      fit1->SetParLimits(0,0.000,2.*hpcopy->GetMaximum()) ;
      fit1->SetParLimits(1,0.125,0.145) ;
      fit1->SetParLimits(2,0.0045,0.010) ;

      Double_t rangeMin=0.07 ;//TMath::Max(0.06,0.11-0.01*i) ;
      Double_t rangeMax=0.22; //TMath::Min(0.25,0.18+0.01*i) ;
//       if(pt>10){
//         rangeMax=0.18; 
//         rangeMin=0.10; 
//       }
//      Double_t rangeMin=TMath::Max(0.06,0.12-0.01*i) ;
//      Double_t rangeMax=TMath::Min(0.25,0.16+0.01*i) ;

      hpcopy->Fit(fit1,"Q","",rangeMin,rangeMax) ;
      hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;

      fit2->SetParameters(fit1->GetParameters()) ;
      fit2->SetParameter(3,3.2) ;
      fit2->SetParameter(4,1.56) ;
      fit2->SetParameter(5,fit1->GetParameter(3)) ;
      fit2->SetParameter(6,fit1->GetParameter(4)) ;
      fit2->SetParameter(7,fit1->GetParameter(5)) ;      
      fit2->SetParLimits(0,0.000,2.*hpcopy->GetMaximum()) ;
      fit2->SetParLimits(1,0.125,0.145) ;
      fit2->SetParLimits(2,0.0045,0.010) ;
      
      hpcopy->Fit(fit2,"+QN","",rangeMin,rangeMax) ;
      hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;
      c1[iPID]->cd(i) ;
      hpcopy->Draw() ;
      if(i<17)
        c1[iPID]->Update() ;
      else
        c1b[iPID]->Update() ;

      if(i<17)
        c2[iPID]->cd(i) ;
      else
        c2b[iPID]->cd(i-16) ;
      
      fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4),fit1->GetParameter(5),fit1->GetParameter(6)); 
      fbg2->SetParameters(fit2->GetParameter(5),fit2->GetParameter(6),fit2->GetParameter(7),fit2->GetParameter(8)); 
      
      Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ;
      Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ;
      Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;
      Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;
      Double_t errStat     = hpm->Integral(intBinMin,intBinMax); 

//.........这里部分代码省略.........
	     -0.545069,
	     -0.281830,
	     0.026143,
	     2.549494,
	     0.,
	     0.
	   };

	   for(int ipar=0; ipar <NRETURN; ipar++)
	     myFitPar[ipar] = tempPar[ipar];
	 } // end if EE, Et=40 GeV


      Double_t sigFitPar[NPAR]={0};
      for(int ipar=0; ipar<NPAR; ipar++)
	sigFitPar[ipar] = myFitPar[ipar];

      Double_t bkgFitPar[NPAR]={0};
      for(int ipar=0; ipar<NPAR; ipar++)
	bkgFitPar[ipar] = myFitPar[ipar+NPAR];


      Double_t sumFitPar[NPAR]={0};
      for(int ipar=0; ipar<NPAR; ipar++)
	sumFitPar[ipar] = myFitPar[ipar+NPAR*2];


      TF1* fsig = new TF1("fsig", exp_conv, fit_lo, fit_hi, 11);
      fsig->SetParameters(sigFitPar);       
      fsig->SetParameter(0,1.0);


      // changing the signal tail
      Double_t current = fsig->GetParameter(1);
      cout << "Current parameter = " << current  << endl;
      current = fCentralTail[ieta][ipt];
      fsig->SetParameter(1, current);
      cout << "Now Current parameter = " << fsig->GetParameter(1)  << endl;

      mySigPDFnorm = fsig->Integral(fit_lo,fit_hi);
      cout << "mySigPDFnorm = " << mySigPDFnorm << endl;
      
      TF1* fbkg = new TF1("fbkg", expinv_power, fit_lo, fit_hi, 11);
      fbkg->SetParameters(bkgFitPar);
      fbkg->SetParameter(4,1.0);

      myBkgPDFnorm = fbkg->Integral(fit_lo, fit_hi);
      cout << "myBkgPDFnorm = " << myBkgPDFnorm << endl;
 
      TF1* fsum = new TF1("fsum",mysum_norm, fit_lo, fit_hi,11);
      fsum->SetParameters(sumFitPar);

      cout << "Using nsig_input = " << nsig_input << endl;
      cout << "Using nbkg_input = " << nbkg_input << endl;
      cout << "Using nsig_input5GeV = " << nsig_input5GeV << endl;

      fsum->SetParameter(0, nsig_input*hdata_data[ieta][ipt]->GetBinWidth(2));
      fsum->SetParameter(4, nbkg_input*hdata_data[ieta][ipt]->GetBinWidth(2));
      fsum->SetLineColor(2);
      fsum->SetNpx(2500);

      //        cout << "fsum integral = " << fsum->Integral(fit_lo,fit_hi) << endl;
      //        for(int ipar=0; ipar<NPAR; ipar++)cout << "fsum par " << ipar << " = " << fsum->GetParameter(ipar) << endl;

//       FILE *infile =  fopen("/afs/cern.ch/user/s/syu/scratch0/LxplusArea/EGdata_comb3Iso_et.dat","r");  
//       Double_t xdata, xdata1, xdata2; // combined isolation, pt, eta

TF1 *fit(TTree *nt,TTree *ntMC,double ptmin,double ptmax){   
   //cout<<cut.Data()<<endl;
   static int count=0;
   count++;
   TCanvas *c= new TCanvas(Form("c%d",count),"",600,600);
   TH1D *h = new TH1D(Form("h%d",count),"",50,5,6);
   TH1D *hMC = new TH1D(Form("hMC%d",count),"",50,5,6);

   TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
   TF1 *f = new TF1(Form("f%d",count),"[0]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*("+iNP+")");
   nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   clean0(h);
   h->Draw();
   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.05);
   f->SetParLimits(8,0.01,0.05);
   f->SetParLimits(7,0,1);
   f->SetParLimits(5,0,1000);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(8,setparam3);
   f->FixParameter(1,fixparam1);
   h->GetEntries();

   hMC->Fit(Form("f%d",count),"q","",5,6);
   hMC->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L m","",5,6);

   cout<<"======= MC ======="<<endl;
   cout<<f->GetParameter(2)<<"  "<<f->GetParameter(8)<<endl;
   cout<<"===== MC end ====="<<endl;
   cout<<endl;

   f->FixParameter(1,f->GetParameter(1));
   //f->FixParameter(2,f->GetParameter(2));
   f->FixParameter(7,f->GetParameter(7));
   //f->FixParameter(8,f->GetParameter(8));
   
   h->Fit(Form("f%d",count),"q","",5,6);
   h->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L m","",5,6);
   h->SetMarkerSize(0.8);
   h->SetMarkerStyle(20);
   cout <<h->GetEntries()<<endl;

   cout<<"======= data ======="<<endl;
   cout<<f->GetParameter(2)<<"  "<<f->GetParameter(8)<<endl;
   cout<<"===== data end ====="<<endl;
   cout<<endl;
   cout<<"======= chi2 ======="<<endl;
   cout<<f->GetChisquare()<<endl;
   cout<<"===== chi2 end ====="<<endl;

   // function for background shape plotting. take the fit result from f
   TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
   background->SetParameter(0,f->GetParameter(3));
   background->SetParameter(1,f->GetParameter(4));
   background->SetParameter(2,f->GetParameter(5));
   background->SetParameter(3,f->GetParameter(6));
   background->SetLineColor(4);
   background->SetRange(5,6);
   background->SetLineStyle(2);
   
   // function for signal shape plotting. take the fit result from f
   TF1 *Bkpi = new TF1(Form("fBkpi",count),"[0]*("+iNP+")");
   Bkpi->SetParameter(0,f->GetParameter(5));
   Bkpi->SetLineColor(kGreen+1);
   Bkpi->SetFillColor(kGreen+1);
//   Bkpi->SetRange(5.00,5.28);
   Bkpi->SetRange(5.00,6.00);
   Bkpi->SetLineStyle(1);
   Bkpi->SetFillStyle(3004);

   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1(Form("fmass",count),"[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(8));
   mass->SetParError(0,f->GetParError(0));
   mass->SetParError(1,f->GetParError(1));
   mass->SetParError(2,f->GetParError(2));
   mass->SetParError(7,f->GetParError(7));
   mass->SetParError(8,f->GetParError(8));
   mass->SetLineColor(2);
   mass->SetLineStyle(2);

//   cout <<mass->Integral(0,1.2)<<" "<<mass->IntegralError(0,1.2)<<endl;
   h->SetMarkerStyle(24);
   h->SetStats(0);
   h->Draw("e");
   h->SetXTitle("M_{B} (GeV/c^{2})");
//.........这里部分代码省略.........

void plot()
{
  TGraph * gr1 = new TGraph();
  gr1->SetMarkerStyle(21);
  int nPoints1 = 0;

  TGraph * gr2 = new TGraph();
  gr2->SetMarkerStyle(21);
  int nPoints2 = 0;

  TGraph * gr3 = new TGraph();
  gr3->SetMarkerStyle(21);
  int nPoints3 = 0;

  TGraph * gr4 = new TGraph();
  gr4->SetMarkerStyle(21);
  int nPoints4 = 0;

  TGraph * gr5 = new TGraph();
  gr5->SetMarkerStyle(21);
  int nPoints5 = 0;

  while (M1[nPoints1].bath!=-999.0) {
    gr1->SetPoint(nPoints1, M1[nPoints1].ambient, M1[nPoints1].slope);
    nPoints1++;

    gr2->SetPoint(nPoints2, M1[nPoints2].bath, M1[nPoints2].ambient);
    nPoints2++;

    gr3->SetPoint(nPoints3, M1[nPoints3].bottom, M1[nPoints3].ambient);
    nPoints3++;

    gr4->SetPoint(nPoints4, M1[nPoints4].top, M1[nPoints4].ambient);
    nPoints4++;

    gr5->SetPoint(nPoints5, M1[nPoints5].top, M1[nPoints5].slope);
    nPoints5++;
  }

  TF1 * fit = new TF1("fit", myfunction, 5, 35, 2);
  fit->SetLineColor(2);

  gr3->Draw("AP");
  gr3->Fit(fit, "NR");
  fit->Draw("same");

  double Tbottom = fit->GetParameter(0) / (1.0 - fit->GetParameter(1));
  cout << Tbottom << endl;

  gr4->Draw("AP");
  gr4->Fit(fit, "NR");
  fit->Draw("same");

  double Ttop = fit->GetParameter(0) / (1.0 - fit->GetParameter(1));
  cout << Ttop << endl;

  gr5->Draw("AP");

  gr5->Fit(fit, "NR");
  fit->Draw("same");

  cout << fit->Eval(Tbottom) << endl;
  cout << fit->Eval(Ttop) << endl;
  cout << fit->Eval(Ttop + (Ttop-Tbottom)) << endl;
}

void fitD(TString inputdata="/data/dmeson2015/DataDntuple/nt_20160112_DfinderData_pp_20160111_dPt0tkPt1_D0Dstar3p5p_DCSJSON_v2.root", TString inputmc="/afs/cern.ch/work/w/wangj/public/Dmeson/ntD_20151110_DfinderMC_20151110_EvtMatching_Pythia_D0pt15p0_Pthat15_TuneZ2_5020GeV_GENSIM_75x_1015_20151110_ppGlobaTrackingPPmenuHFlowpuv11_MBseed_twang-Pythia_1107.root", TString trgselection="((HLT_DmesonPPTrackingGlobal_Dpt15_v1&&Dpt>25&&Dpt<40)||(HLT_DmesonPPTrackingGlobal_Dpt30_v1&&Dpt>40&&Dpt<60)||(HLT_DmesonPPTrackingGlobal_Dpt50_v1&&Dpt>60))",  TString cut="Dy>-1.&&Dy<1.&&(Dtrk1highPurity&&Dtrk2highPurity)&&(DsvpvDistance/DsvpvDisErr)>3.5&&Dchi2cl>0.05&&Dalpha<0.12&&Dtrk1Pt>1.5&&Dtrk2Pt>1.5", TString selmcgen="((GisSignal==1||GisSignal==2)&&(Gy>-1&&Gy<1))", int isMC=0, Double_t luminosity=26., int doweight=0, TString collsyst="PbPb", TString outputfile="mytest.root")
{
  collisionsystem=collsyst;
  seldata = Form("%s&&%s",trgselection.Data(),cut.Data());
  selmc = Form("%s",cut.Data());

  gStyle->SetTextSize(0.05);
  gStyle->SetTextFont(42);
  gStyle->SetPadRightMargin(0.043);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadBottomMargin(0.145);
  gStyle->SetTitleX(.0f);

  void clean0 (TH1D* h);
  TF1* fit (TTree* nt, TTree* ntMC, double ptmin, double ptmax, int isMC);

  if(!doweight) weight="1";
  TFile* inf = new TFile(inputdata.Data());
  TFile* infMC = new TFile(inputmc.Data());

  TTree* nt = (TTree*) inf->Get("ntDkpi");
  TTree* HltTree= (TTree*) inf->Get("ntHlt");
  HltTree->AddFriend(nt);
  nt->AddFriend(HltTree);
  TTree* ntHid = (TTree*) inf->Get("ntHi");
  nt->AddFriend(ntHid);
  
  TTree* ntMC = (TTree*)infMC->Get("ntDkpi");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  TTree* ntHi = (TTree*)infMC->Get("ntHi");
  
  ntGen->AddFriend(ntMC);
  ntGen->AddFriend(ntHi);
  ntMC->AddFriend(ntGen);
  ntMC->AddFriend(ntHi);
  ntHi->AddFriend(ntMC);
  
  TH1D* hPt = new TH1D("hPt","",nBins,ptBins);
  TH1D* hPtRecoTruth = new TH1D("hPtRecoTruth","",nBins,ptBins);
  TH1D* hPtMC = new TH1D("hPtMC","",nBins,ptBins);
  TH1D* hPtGen = new TH1D("hPtGen","",nBins,ptBins);

  TH1D* hMean = new TH1D("hMean","",nBins,ptBins);                       
  TH1D* hSigmaGaus1 = new TH1D("hSigmaGaus1","",nBins,ptBins); 
  TH1D* hSigmaGaus2 = new TH1D("hSigmaGaus2","",nBins,ptBins); 
  TH1D* hRelMagnGaus1Gaus2 = new TH1D("hRelMagnGaus1Gaus2","",nBins,ptBins); 
  
  for(int i=0;i<nBins;i++)
    {
      TF1* f = fit(nt,ntMC,ptBins[i],ptBins[i+1],isMC);
      double yield = f->Integral(minhisto,maxhisto)/binwidthmass;
      double yieldErr = f->Integral(minhisto,maxhisto)/binwidthmass*f->GetParError(0)/f->GetParameter(0);
      hPt->SetBinContent(i+1,yield/(ptBins[i+1]-ptBins[i]));
      hPt->SetBinError(i+1,yieldErr/(ptBins[i+1]-ptBins[i]));
      hMean->SetBinContent(i+1,f->GetParameter(1));
      hMean->SetBinError(i+1,f->GetParError(1));
      hSigmaGaus1->SetBinContent(i+1,f->GetParameter(2));
      hSigmaGaus1->SetBinError(i+1,f->GetParError(2));
      hSigmaGaus2->SetBinContent(i+1,f->GetParameter(5));
      hSigmaGaus2->SetBinError(i+1,f->GetParError(5));
      hRelMagnGaus1Gaus2->SetBinContent(i+1,f->GetParameter(4));
      hRelMagnGaus1Gaus2->SetBinError(i+1,f->GetParError(4));
    }  

  ntMC->Project("hPtMC","Dpt",TCut(weight)*(TCut(selmc.Data())&&"(Dgen==23333)"));
  divideBinWidth(hPtMC);
  ntMC->Project("hPtRecoTruth","Dpt",TCut(selmc.Data())&&"(Dgen==23333)");
  divideBinWidth(hPtRecoTruth);
  ntGen->Project("hPtGen","Gpt",TCut(weight)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);

  TCanvas* cPt =  new TCanvas("cPt","",600,600);
  cPt->SetLogy();
  hPt->SetXTitle("D^{0} p_{T} (GeV/c)");
  hPt->SetYTitle("Uncorrected dN(D^{0})/dp_{T}");
  hPt->Sumw2();
  hPt->Draw();
  if(isMC)
    {
      hPtMC->Draw("same hist");
      TLegend* legPt = myLegend(0.55,0.80,0.90,0.94);
      legPt->AddEntry(hPt,"Signal extraction","pl");
      legPt->AddEntry(hPtMC,"Matched reco","lf");
      legPt->Draw("same");  
    }
  hPtMC->Sumw2();
  TH1D* hEff = (TH1D*)hPtMC->Clone("hEff");
  hEff->SetTitle(";D^{0} p_{T} (GeV/c);Efficiency");
  hEff->Sumw2();
  hEff->Divide(hPtGen);
  TCanvas* cEff = new TCanvas("cEff","",600,600);
  hEff->Draw();
  
  TH1D* hPtCor = (TH1D*)hPt->Clone("hPtCor");
  hPtCor->SetTitle(";D^{0} p_{T} (GeV/c);Corrected dN(D^{0})/dp_{T}");
  hPtCor->Divide(hEff);
  TCanvas* cPtCor=  new TCanvas("cCorResult","",600,600);
  cPtCor->SetLogy();
  hPtCor->Draw();
//.........这里部分代码省略.........