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

//---------------------
// FitSection
//---------------------
TF1* FitSection(int order, TTree *Bfield, const char *Bi, Double_t z_min_inches, Double_t z_max_inches, Double_t r, const char *fname)
{
	cout<<"Fitting "<<fname<<" at r="<<r<<endl;

	// The following should be defined based on the grid spacing/location
	// of points in the Bfield tree. 
	Double_t delta_z = 2.54; // in cm
	Double_t delta_r = 2.54; // in cm

	// Determine proper histogram limits 
	Double_t z_min = z_min_inches*delta_z - delta_z/2.0;
	Double_t z_max = z_max_inches*delta_z + delta_z/2.0;
	Double_t epsilon = delta_z/1000.0; // used to prevent round-off problems
	UInt_t NbinsZ = floor((z_max - z_min + epsilon)/delta_z);

	// Note that here we explicitly set the limits to be -1 to +1 so that the
	// histogram conforms to the requirements of the chebyshev_FindBestFunction
	// function.
	TH1D *Bi_vs_znorm = new TH1D("Bi_vs_znorm", "Bi as function of modified z coordinate", NbinsZ, z_min, z_max);

	// Make varexp that properly converts from the lab z-coordinates,
	// to the "Chebyshev compatible" coordinates where z_min=-1 and
	// z_max=+1.
	stringstream cut;
	cut<<Bi<<"*(abs(r-"<<r<<")<"<<delta_r/10.0<<")"; // select a single r bin
	Bfield->Project("Bi_vs_znorm", "z", cut.str().c_str());

	TF1 *func = chebyshev_Fit(Bi_vs_znorm, order);
	TF1 *myfunc = (TF1*)func->Clone(fname);

	delete Bi_vs_znorm;
	
	return myfunc;
}

void fitmmp(TH1 *hmmp, bool verbose = false) {
  TString options;
  if (verbose) options = "";
  else options = "Q";
  TF1 *fbg = f_pol4("fbg",0.4,2,true);
  hmmp->Fit(fbg,options.Data(),"",0.42,2);
  //printf("%s\n",gMinuit->fCstatu.Data());
  if (true) { //gMinuit->fCstatu.Contains("CONVER")) {
    TF1 *fbg2 = f_pol4("fbg",0.4,2,false);
    fbg2->SetParameters(fbg->GetParameters());
    //fbg2->Draw("same");
    fbg2->SetParameter(5,0);
    fbg2->SetParameter(6,0);
    TH1 *htmp = (TH1*)hmmp->Clone("hmmp_bgsub");
    htmp->Add(fbg2,-1);
    //htmp->Draw();
    TF1 *fgaus = new TF1("fgaus","gaus",0.4,2);
    htmp->Fit(fgaus,options.Data(),"",0.74,0.85);
    TF1 *f = f_pol4gaus("f",0.4,2,fbg2,fgaus);
    f->SetNpx(500);
    hmmp->Fit(f,options.Data(),"",0.42,2);
    fgaus->SetRange(0.4,2);
    for (int i = 0; i < 3; i++) fgaus->SetParameter(i,f->GetParameter(i+5));
    for (int i = 0; i < 5; i++) fbg2->SetParameter(i,f->GetParameter(i));
    fbg2->SetLineStyle(2);
    fbg2->SetLineColor(kRed+1);
    fgaus->SetLineStyle(2);
    fgaus->SetLineColor(kGreen+1);
    hmmp->GetListOfFunctions()->Add(fbg2->Clone());
    hmmp->GetListOfFunctions()->Add(fgaus->Clone());
    delete fbg2;
    delete htmp;
    delete fgaus;
    delete f;
  } else hmmp->GetListOfFunctions()->Delete();
  delete fbg;
}

void fitDstar(char *infname = "/data/wangj/MC2015/Dntuple/PbPb/ntD_Pythia8_5020GeV_DstarD0kpipipi_755patch3_GEN_SIM_PU_20151120_Dstar5p_tkPt2_20151126_Evt_All.root") {

    TFile *inf = new TFile(infname);
    //TTree *ntmix=(TTree*)inf->Get("ntDD0kpipipipi");
    TTree *ntmix=(TTree*)inf->Get("ntDD0kpipi");
    TH1D *h = new TH1D("h","",100,0.139,0.159);
    TCut cutTrk = "";//"trk1PixelHit>=2&&trk1StripHit>=10&&trk1Chi2ndf<5&&trk2PixelHit>=2&&trk2StripHit>=10&&trk2Chi2ndf<5";

    TCanvas *c = new TCanvas("c","",750,750);
    ntmix->Draw("Dmass-DtktkResmass>>h","abs(DtktkResmass-1.86486)<0.015&&Dpt>10&&(DsvpvDistance/DsvpvDisErr)>0.&&Dchi2cl>0.05&&Dalpha<1."&&cutTrk,"",10000000);
    //ntmix->Draw("Dmass-DtktkResmass>>h","abs(DtktkResmass-1.86486)<0.015&&Dpt>10&&Dtrk1Pt>0.5&&DRestrk1Pt>0.5&&DRestrk2Pt>0.5&&DRestrk3Pt>0.5&&DRestrk4Pt>0.5&&Dchi2cl>0.1&&Dalpha<0.2&&(DsvpvDistance/DsvpvDisErr)>.0"&&cutTrk,"",10000000);
    h->Sumw2();
    //   TF1 *f = new TF1("f","(1-exp(-(x-[8])/[0]))*(((x)/[8])**[1]+[2]*(((x)/[8])-1))*[3]+[4]*(TMath::Voigt(x-[5],[6],[7]))");
    TF1* f = new TF1("f","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*((1-[8])*TMath::Gaus(x,[6],[7])/(sqrt(2*3.14159)*[7])+[8]*TMath::Gaus(x,[6],[9])/(sqrt(2*3.14159)*[9]))",minmass3prong,maxmass3prong);
//   TF1 *f = new TF1("f","(1-exp(-(x-0.13957018)/[0]))*(((x)/0.13957018)**[1]+[2]*(((x)/0.13957018)-1))*[3]+[4]*(TMath::Gaus(x,[5],[6]))");
    //TF1 *f = new TF1("f","(1-exp(-(x-1.86486-0.13957018)/[0]))*(((x-1.86486)/0.13957018)**[1]+[2]*(((x-1.86486)/0.13957018)-1))*[3]+[4]*TMath::Gaus(x,[5],[6])");
    f->SetLineColor(4);
//   f->SetParameters(-2.3825e6,-7.99713e-1,-1.42957,-5.50069e10,5.33573,1.45491e-1,2.78677e-6,1.43145e-3,0.13957018);
//   f->SetParameters(-7.3e5,-2.2e1,5.24e-1,-7.18e9,2e2,1.45491e-1,9e-4,0.1,8e-4);
    f->SetParameters(0,0,0,0,0,2e2,1.45491e-1,9e-4,0.1,8e-4);
    f->FixParameter(9,15e-4);
    f->FixParameter(6,0.145491);
    f->FixParameter(7,8e-4);
    f->SetParLimits(8,0,1);
    h->Fit("f","LL");
    h->Fit("f","LL");
    h->Fit("f","LL","",0.142,0.147);
    f->ReleaseParameter(6);
    f->ReleaseParameter(7);
    f->ReleaseParameter(9);
    f->SetParLimits(7,1e-4,9e-4);
    f->SetParLimits(9,1e-4,9e-4);
    h->Fit("f","LL","",0.142,0.148);
    h->Fit("f","LL","",0.142,0.16);
    h->Fit("f","LL","",0.142,0.16);
    h->Fit("f","LL","",0.141,0.16);
    h->Fit("f","LL","",0.141,0.16);
    h->Fit("f","LL","",0.141,0.16);

    h->SetXTitle("M_{K#pi#pi#pi#pi}-M_{K#pi#pi#pi} (GeV/c^{2})");
    h->SetYTitle("Entries");
    h->SetStats(0);
    h->SetAxisRange(1,h->GetMaximum()*1.3,"Y");
    TF1 *f2 = (TF1*)f->Clone("f2");
    f2->SetParameter(5,0);
    f2->SetRange(0.141,0.16);
    TF1 *f3 = (TF1*)f->Clone("f3");
    f3->SetParameter(0,0);
    f3->SetParameter(1,0);
    f3->SetParameter(2,0);
    f3->SetParameter(3,0);
    f3->SetParameter(4,0);

    f->SetLineColor(4);
    f2->SetLineColor(4);
    f2->SetLineStyle(2);
    f3->SetLineStyle(2);
    f2->Draw("same");
    f3->SetLineColor(2);
    f3->SetFillStyle(3004);
    f3->SetFillColor(2);
    f3->Draw("same");

    c->SaveAs("canvasDstar.pdf");
}

void mk_sigaccanplots(string flavor = "112", bool tdrstyle = false)
{
	setTDRStyle(tdrstyle);
	// gStyle->SetOptFit(1100); // chi2 and prob, not parameters
	gStyle->SetOptFit(0); // chi2 and prob, not parameters
	// gStyle->SetOptStat("irme"); // integral, RMS, mean, # of entries
	gStyle->SetStatFontSize(0.005);
	gStyle->SetStatY(0.4);
	float unused = 0.9, simtxt_y = 0.97;
	if (tdrstyle == false) {
		simtxt_y = 0.93;
	}
	TLatex *tex = new TLatex(unused, simtxt_y, "CMS simulation at #sqrt{s} = 8 TeV");
	tex->SetNDC();
	tex->SetTextAlign(12); // Left-adjusted
	tex->SetTextFont(42);
	float textsiz = 0.06;
	tex->SetTextSize(0.05);
	tex->SetLineWidth(2);
  // for(int k=0; k<3; k++){	
  // for(int k=0; k<2; k++){	
  for(int k=2; k<3; k++){	
  //Which plots to make
  //define input variables
  vector <TFile* > filelist_kin;

  // string uncert="btagup";
  // if(k==1) uncert="btagdown";
  string uncert="up";
  if(k==1) uncert="down";
  if(k==2) uncert="";
  
  // filelist_kin.push_back(TFile::Open(("Acc_Gluino_NoLumi_pt110_8TeVxsec_BtagMap_2500GeV"+uncert+".root").c_str()));
  // filelist_kin.push_back(TFile::Open("Acc_Gluino_NoLumi_pt60_pt110_8TeVxsec_BtagMap.root"));
  // filelist_kin.push_back(TFile::Open(("Acc_RPV112" + uncert + ".root").c_str()));
  filelist_kin.push_back(TFile::Open(("Acc_RPV" + flavor + uncert + "_Sph4.root").c_str()));
//   filelist_kin.push_back(TFile::Open(("Acc_Gluino_NoLumi_pt60_pt110_8TeVxsec_BtagMap" + uncert + ".root").c_str()));
  
  TFile f1(("RPV_" + flavor + "_accandwdithtest"+uncert+".root").c_str(), "recreate");
  f1.cd();
  
  
  
  string nameIN;
  string cuts;
  string prefix;
  string postfix;
  string folder;
  string ptcut;
  folder="plots_paper/";
  postfix=".pdf";
  /////////////Plots for each Mass separatly//////////////
  /////////////----------------------------------------/////////////
  
  for (int p = 0; p < 2; p++) {
    // ptcut="112";
    // if(p==1) ptcut="" + flavor;
     ptcut="60";
    if (p==1) {
		ptcut="110";
		textsiz = 0.045;
    }
    // string histname = string("GausWidth_vs_Mass_112") + ptcut;
      string histname = string("GausWidth_vs_Mass_" + flavor) + "_" + ptcut;
      cout<< histname << endl;
      TGraphErrors *h_GluinoHist_Fit = (TGraphErrors*) filelist_kin[0]->Get(histname.c_str())->Clone();
      histname = string("GausAcceptance_vs_Mass_" + flavor) + "_" + ptcut;
      // histname = string("GausAcceptance_vs_Mass_112") + ptcut;
      cout<< histname << endl;
      TGraphErrors *h_GluinoHist_MCcomb = (TGraphErrors*) filelist_kin[0]->Get(histname.c_str())->Clone();
      
      // histname = string("FullAcceptance_vs_Mass_" + flavor) + "_" + ptcut;
      histname = string("GausMean_vs_Mass_" + flavor) + "_" + ptcut;
      // histname = string("GausAcceptance_vs_Mass_112") + ptcut;
      cout<< histname << endl;
      // TGraphErrors *h_FullAccept = (TGraphErrors*) filelist_kin[0]->Get(histname.c_str())->Clone();
      TGraph *h_FullAccept = (TGraph*) filelist_kin[0]->Get(histname.c_str())->Clone();
     
      histname = string("GausMeanOffset_vs_Mass_" + flavor) + "_" + ptcut;
      cout<< histname << endl;
      TGraph *h_MeanOffset = (TGraph*) filelist_kin[0]->Get(histname.c_str())->Clone();
     
      TCanvas * cGluinoFitsOpti = new TCanvas(("RPV_"+flavor +ptcut).c_str(), ("RPV_" + ptcut+"_"+cuts).c_str(), 800, 600);
      //h_GluinoHist_Fit->SetFillColor(kOrange-2);
      // h_GluinoHist_Fit->SetLineColor(kBlack);
		string title;
		string systematic = "pile-up";
		// title="Gaussian Width vs. Mass for Light-ptcut RPV";
		string tag = "qqb", sphericity = " Sphericity #geq 0.4";
		if (flavor.compare("112") == 0)
			tag = "qqq";
		else if (ptcut.compare("60") == 0)
			sphericity = "";
		string titlepart = "Hadronic RPV #tilde{g} #rightarrow " + tag;
		string titleln2 = "\\Delta = 110 GeV";
		string titleln3 = "6^{th}-jet p_{T} #geq " + ptcut + " GeV";
		title = titlepart;
		// title = "Width for " + titlepart;
		/*
		if(k==0){
//.........这里部分代码省略.........

void fitmmp(TH1 *hmmp, float step_x0 = 2, float step_x1 = 2.1, float wval = 0, bool verbose = false) {
  float fnrangemax = 1.1;
  float mmplow = 0.4;
  //float bgp2 = 2.1e8;
  float gsig = -0.0065+0.013*wval;
  float pol2cutoff = 4; //wval under which to use pol2
  float fitcutoff = 1.1;
  TString options;
  if (verbose) options = "";
  else options = "Q";
  TF1 *fbg = f_pol4("fbg",0.4,fnrangemax,true);
  fbg->FixParameter(7,step_x0);
  fbg->FixParameter(8,step_x1);
  if (wval<pol2cutoff) {
    fitcutoff = 1.1;
    fbg->FixParameter(3,0);
    fbg->FixParameter(4,0);
  } else fitcutoff = 2;
  hmmp->Fit(fbg,options.Data(),"",mmplow,fitcutoff);
  if (true) { //gMinuit->fCstatu.Contains("CONVER")) {
    TF1 *fbg2 = f_pol4("fbg",0.4,fnrangemax,false);
    fbg2->SetParameters(fbg->GetParameters());
    fbg2->FixParameter(7,step_x0);
    fbg2->FixParameter(8,step_x1);
    if (wval<pol2cutoff) {
      fbg2->FixParameter(3,0);
      fbg2->FixParameter(4,0);
    }
    fbg2->SetParameter(5,0);
    fbg2->SetParameter(6,0);
    TH1 *htmp = (TH1*)hmmp->Clone("hmmp_bgsub");
    htmp->Add(fbg2,-1);
    TF1 *fgaus = new TF1("fgaus","gaus",0.4,fnrangemax);
    fgaus->FixParameter(2,gsig);
    fgaus->FixParameter(1,0.783);
    float gmax = 0.85;
    if (gmax > step_x0) gmax = step_x0;
    htmp->Fit(fgaus,options.Data(),"",0.74,gmax);
    TF1 *f = f_pol4gaus("f",0.4,fnrangemax,fbg2,fgaus);
    f->FixParameter(8,step_x0);
    f->FixParameter(9,step_x1);
    if (wval<pol2cutoff) {
      f->FixParameter(3,0);
      f->FixParameter(4,0);
    }
    f->SetNpx(500);
    f->FixParameter(6,0.783);
    f->FixParameter(7,gsig);
    hmmp->Fit(f,options.Data(),"",mmplow,fitcutoff);

    //f->FixParameter(2,0.9*f->GetParameter(2));
    //hmmp->Fit(f,options.Data(),"",mmplow,fitcutoff);

    fgaus->SetRange(0.4,fnrangemax);
    for (int i = 0; i < 3; i++) fgaus->SetParameter(i,f->GetParameter(i+5));
    for (int i = 0; i < 5; i++) fbg2->SetParameter(i,f->GetParameter(i));
    fbg2->SetLineStyle(2);
    fbg2->SetLineColor(kRed+1);
    fgaus->SetLineStyle(2);
    fgaus->SetLineColor(kGreen+1);
    hmmp->GetListOfFunctions()->Add(fbg2->Clone());
    hmmp->GetListOfFunctions()->Add(fgaus->Clone());
    float ymax = 1.1*hmmp->GetMaximum();
    hmmp->GetYaxis()->SetRangeUser(0,ymax);
    TLine *l1 = new TLine(step_x0,0,step_x0,ymax);
    TLine *l2 = new TLine(step_x1,0,step_x1,ymax);
    l1->SetLineColor(kBlue+1);
    l2->SetLineColor(kBlue+1);
    hmmp->GetListOfFunctions()->Add(l1);
    hmmp->GetListOfFunctions()->Add(l2);
    delete fbg2;
    delete htmp;
    delete fgaus;
    delete f;
  } else hmmp->GetListOfFunctions()->Delete();
  delete fbg;
}

void fit(float bgpar2smudge=1.0) {
  TF1 *fsig = new TF1("fsig",&d_sig,0.4,2,5);
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  //int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\tstepx0\tstepx1\txsectFn\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      delete [] wq;
      int wbin = h3->GetYaxis()->FindBin(wval);
      float wlow = h3->GetYaxis()->GetBinLowEdge(wbin);
      float step_x0 = sqrt(wlow*wlow+MASS_P*MASS_P-2*wlow*MASS_P);
      float whigh = h3->GetYaxis()->GetBinLowEdge(wbin+1);
      float step_x1 = sqrt(whigh*whigh+MASS_P*MASS_P-2*whigh*MASS_P);
      fitmmp(h,step_x0,step_x1,wval);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      //int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);

      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      fsig->SetParameter(0,ftmp->GetParameter(5));
      fsig->SetParameter(1,ftmp->GetParameter(6));
      fsig->SetParameter(2,ftmp->GetParameter(7));
      fsig->SetParameter(3,step_x0);
      fsig->SetParameter(4,step_x1);
      fsig->SetLineWidth(2);
      fsig->SetLineColor(kBlue+1);
      h->GetListOfFunctions()->Add((TF1*)fsig->Clone("fsig"));
      //fsig->Print();
      double Nfn = 0;
      for (int b = 1; b < h->GetNbinsX(); b++) {
        double x = h->GetXaxis()->GetBinCenter(b);
        Nfn += fsig->Eval(x);
      }
      //printf("**** %.3e\t\%.3e\n",Nfn,N);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double xsectFn = Nfn/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile,"\t%.3e",xsectFn/(1e6));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// root -l       computeDrellYanPtllWeight.C+
// root -l -b -q computeDrellYanPtllWeight.C+
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void computeDrellYanPtllWeight(TString fname = "h_pt2l_mm")
{
  gInterpreter->ExecuteMacro("PaperStyle.C");

  TFile* file = TFile::Open("figures/Control/01_DY/" + fname + ".root");

  TH1D* ratio = (TH1D*)file->Get("ratio");


  // Draw the ratio
  //----------------------------------------------------------------------------
  TCanvas* c1 = new TCanvas("c1", "c1");

  ratio->SetMinimum(0.85);
  ratio->SetMaximum(1.20);

  ratio->Draw("ep");

  TString ytitle = Form("data / MC ratio / %.0f GeV", ratio->GetBinWidth(0));

  SetAxis(ratio, ratio->GetXaxis()->GetTitle(), ytitle, 1.7, 1.8);


  // Draw the old function
  //----------------------------------------------------------------------------
  TF1* fOld = new TF1("fOld", "[3]*(0.95 - [0]*TMath::Erf((x-[1])/[2]))", 0, 90);

  fOld->SetLineColor  (kGreen+2);
  fOld->SetMarkerColor(kGreen+2);

  // https://github.com/latinos/PlotsConfigurations/blob/master/Configurations/ggH/nuisances_iteos.py#L969-L977
  fOld->SetParameter(0,     0.1);
  fOld->SetParameter(1,    14.0);
  fOld->SetParameter(2,     8.8);
  fOld->SetParameter(3, 1.08683);

  fOld->Draw("same");


  // Draw the old function error band
  //------------------------------------------------------------------------------
  if (errorband)
    {
      TF1* fOld_down = (TF1*)fOld->Clone("fOld_down");
      TF1* fOld_up   = (TF1*)fOld->Clone("fOld_up");

      fOld_down->SetParameter(0, 13.6);
      fOld_down->SetParameter(1,  8.6);

      fOld_up->SetParameter(0, 14.4);
      fOld_up->SetParameter(1,  9.0);

      fOld_down->Draw("same");
      fOld_up  ->Draw("same");
    }


  // Update the fit parameters of the old function
  //
  //   1  p0           6.85257e-02   2.99341e-01   1.42129e-05   4.16760e-04
  //   2  p1           1.24518e+01   4.76906e+01   3.55312e-03   6.41625e-07
  //   3  p2           5.40627e+00   7.83907e+01   6.49546e-03  -1.96376e-06
  //   4  p3           1.05396e+00   3.48880e-01   1.58779e-05  -3.97115e-05
  //
  //----------------------------------------------------------------------------
  TF1* fNew = new TF1("fNew", "[3]*(0.95 - [0]*TMath::Erf((x-[1])/[2]))", 0, 90);

  fNew->SetLineColor  (kRed+1);
  fNew->SetMarkerColor(kRed+1);

  fNew->SetParameter(0, 0.1);
  fNew->SetParameter(1,  10);
  fNew->SetParameter(2,   1);
  fNew->SetParameter(3,   1);

  ratio->Fit(fNew, "mlr0");

  fNew->Draw("same");


  // Draw the updated error band
  //------------------------------------------------------------------------------
  if (errorband)
    {
      TF1* fNew_down = (TF1*)fNew->Clone("fNew_down");
      TF1* fNew_up   = (TF1*)fNew->Clone("fNew_up");

      fNew_down->SetParameter(0, 0.97 * fNew->GetParameter(0));
      fNew_down->SetParameter(1, 0.97 * fNew->GetParameter(1));

      fNew_up->SetParameter(0, 1.03 * fNew->GetParameter(0));
      fNew_up->SetParameter(1, 1.03 * fNew->GetParameter(1));

      fNew_down->Draw("same");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  leg->SetFillStyle(0);
  leg->AddEntry(h,"Data","pl");
  leg->AddEntry(f,"Fit","l");
  leg->AddEntry(mass,"D^{0}+#bar{D^{#lower[0.2]{0}}} Signal","f");
  leg->AddEntry(massSwap,"K-#pi swapped","f");
  leg->AddEntry(background,"Combinatorial","l");
  leg->Draw("same");

  TLatex Tl;
  Tl.SetNDC();
  Tl.SetTextAlign(12);
  Tl.SetTextSize(0.04);
  Tl.SetTextFont(42);
  Tl.DrawLatex(0.18,0.93, "#scale[1.25]{CMS} Performance");
  Tl.DrawLatex(0.65,0.93, Form("%s #sqrt{s_{NN}} = 5.02 TeV",collisionsystem.Data()));

  TLatex* tex;

  tex = new TLatex(0.22,0.78,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax));
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  tex = new TLatex(0.22,0.83,"|y| < 1.0");
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  h->GetFunction(Form("f%d",count))->Delete();
  TH1F* histo_copy_nofitfun = ( TH1F * ) h->Clone("histo_copy_nofitfun");
  histo_copy_nofitfun->Draw("esame");

 TH1D* hTest = new TH1D("hTest","",nbinsmasshisto,minhisto,maxhisto);

 for (int m=0;m<yieldtotal;m++){
 double r = f->GetRandom();
 hTest->Fill(r);
 
 }

 TF1* ffaketotal=(TF1*)f->Clone("ffake");
 TF1* ffakemass=(TF1*)mass->Clone("ffakemass");
 TF1* ffakebackground=(TF1*)background->Clone("ffakebackground");
 TF1* ffakemassSwap=(TF1*)massSwap->Clone("ffakemassSwap");
 
 Double_t yieldtotal_original = ffaketotal->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldmass_original = ffakemass->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldbackground_original = ffakebackground->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldswapped_original = ffakemassSwap->Integral(minhisto,maxhisto)/binwidthmass;

 TH1D* hTestFake = new TH1D("hTestFake","",nbinsmasshisto,minhisto,maxhisto);
  ffakemass->SetParameter(2,ffaketotal->GetParameter(2)*0.8);
  ffakemass->SetParameter(10,ffaketotal->GetParameter(10)*0.8);
 
Double_t yieldmass_modified= ffakemass->Integral(minhisto,maxhisto)/binwidthmass;

cout<<"mass original="<<yieldmass_original<<endl;
cout<<"mass modified="<<yieldmass_modified<<endl;

 for (int m=0;m<yieldmass_original*scalefactor;m++){
   double r = ffakemass->GetRandom();
   hTestFake->Fill(r);

//.........这里部分代码省略.........
	    if (!ga) {
	      ga = new TGraphErrors(0);
	      ga->SetMarkerStyle(g->GetMarkerStyle());
	      ga->SetMarkerColor(g->GetMarkerColor());
	      ga->SetLineColor(g->GetLineColor());
	      gamap[ct][cm][cs][ipt] = ga;
	    }
	    int n = ga->GetN();
	    ga->SetPoint(n, 0.01*a, g->GetY()[i]);
	    ga->SetPointError(n, 0, g->GetEY()[i]);
	  } // for i 

	} // for ialpha

      } // for isample
    } // for imethod
  } // for itype

  cout << "Drawing plots vs pT for each alpha" << endl << flush;

  // 2x6 plots
  for (int itype = 0; itype != ntypes; ++itype) {
    for (int  imethod = 0; imethod != nmethods; ++imethod) {

      const char *ct = types[itype];
      const char *cm = methods[imethod];

      int ipad = ntypes*imethod + itype + 1; assert(ipad<=6);
      c1->cd(ipad);
      gPad->SetLogx();
      h1->SetMaximum(itype<2 ? 1.15 : 1.08);
      h1->SetMinimum(itype<2 ? 0.85 : 0.93);
      h1->SetYTitle(Form("Response (%s)",ct));
      h1->DrawClone("AXIS");
      tex->DrawLatex(0.20,0.85,texlabel[cm]);
      tex->DrawLatex(0.20,0.80,"|#eta| < 1.3, #alpha=0.1--0.3");
      TLegend *leg = tdrLeg(0.60,0.75,0.90,0.90);

      for (int  isample = 0; isample != nsamples; ++isample) {
	for (int  ialpha = 0; ialpha != nalphas; ++ialpha) {

	  const char *cs = samples[isample];
	  const int a = alphas[ialpha];
	  TGraphErrors *g = gemap[ct][cm][cs][a]; assert(g);

	  // Clean out points with very large uncertainty for plot readability
	  for (int i = g->GetN()-1; i != -1; --i) {
	    if (g->GetEY()[i]>0.02) g->RemovePoint(i);
	  }

	  g->Draw("SAME Pz");

	  if (ialpha==0) leg->AddEntry(g,texlabel[cs],"P");
	}
      } // for isample

      // Individual plots for JEC paper
      if ( true ) { // paper

	TH1D *h = new TH1D(Form("h_5%s_%s",ct,cm),
			   Form(";p_{T} (GeV);Response (%s)",ct),
			   1270,30,1300);
	h->GetXaxis()->SetMoreLogLabels();
	h->GetXaxis()->SetNoExponent();
	h->SetMinimum(0.88);
	h->SetMaximum(1.13);

void MakePi0Analysis(){

  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadRightMargin(0.01);
  gStyle->SetPadTopMargin(0.09);
  gStyle->SetPadBottomMargin(0.11);
  //gStyle->SetOptStat(0);
  //gStyle->SetOptTitle(1);
  //gStyle->SetPadTickX(1);
  //gStyle->SetPadTickY(1);
  TGaxis::SetMaxDigits(3);
  
  Double_t bins[] = {1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0, 3.2, 3.4, 3.6,3.8,4.0,4.5, 5.0, 5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,10.0,11,12}; //Nicolas LHC12d PHI
  const Int_t binnum = sizeof(bins)/sizeof(Double_t) - 1;
  
  TString cwd = gSystem->WorkingDirectory();

  TH1F* fHistPeakMean  = new TH1F("fHistPeakMean","",binnum,bins);
  TH1F* fHistPeakWidth = new TH1F("fHistPeakWidth","",binnum,bins);
  TH1F* fHistRawYield  = new TH1F("fHistRawYield","",binnum,bins);
  TH1F* fHistRawYieldPerEvent  = new TH1F("fHistRawYieldPerEvent","",binnum,bins);
  
  TCanvas *c1[binnum];
  TCanvas *c2[binnum];
  TCanvas *c3;

  TFile* ef     = TFile::Open("../AnalysisResults.root");
  TList* list   = (TList*)ef->Get("list_kINT7_Pi0");
  TList *listEv = (TList*)ef->Get("list_kINT7_Event");
  
  TH1F* fHistEvents = (TH1F*)listEv->FindObject("fHistAnalyzedEvents");
  Double_t nEvt = fHistEvents->GetEntries();
  
  TString tofName[] = {"","_TOFcut1","_TOFcut2"};
  TString fitName[] = {"CrystalBall","AsymmGauss"};
  TString modName[] = {"","_M1","_M3"};
  
  Int_t ntof = 3;
  
  Double_t signal_range_min=0.11;
  Double_t signal_range_max=0.15;
  Double_t bg_range_min=0.050;
  Double_t bg_range_max=0.250;
  
  for(Int_t itof=0; itof<ntof; ++itof){
    for(Int_t iMod=0; iMod<1; iMod++){
      for(Int_t iFit=0; iFit<2; iFit++){
	for(Int_t iPol=0; iPol<3; iPol++){
	  
	  Int_t pol = iPol;

	  TString fNameSame = "fHistMassTwoGammas"+tofName[itof]+modName[iMod];
	  TString fNameMix  = "fHistMixMassTwoGammas"+tofName[itof]+modName[iMod];
	  
	  THnSparse *fTHnSparseRecPi0;
	  THnSparse *fTHnSparseRecMixPi0;
	  THnSparse *fTHnSparseGeneratePi0;
	  
	  TH2* hPi0A08    = 0x0;
	  TH2* hMixPi0A08 = 0x0;
	  hPi0A08    = (TH2*)list->FindObject(fNameSame)->Clone();
	  hMixPi0A08 = (TH2*)list->FindObject(fNameMix)->Clone();

	  Int_t sbin = 0;
	  Int_t lbin = 0;
      
	  TH1F* fHistFinalRatio[binnum];
	  TH1F* fHistFinalSame[binnum];
	  TH1F* fHistFinalBG[binnum];
	  TH1F* fHistFinalSignal[binnum];
	  TH1F* fHistOnlyFitSignal[binnum];
	  
	  TF1* fFitFinalRatio[binnum];
	  TF1* fFitFinalSignal[binnum];
	  TF1* fFitOnlyFitSignal[binnum];

	  for(Int_t i=0; i<binnum; i++){    
	    
	    Double_t bin_width = (bins[i+1]-bins[i]);
	    Double_t bin_mean  = (bins[i+1]+bins[i])/2;
	    
	    Int_t rebin = 1;
	    
	    if(bin_mean>5.0){
	      rebin=10;
	    }
	    else{
	      rebin=5;
	    }
	    
	    Int_t sproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i]);
	    Int_t lproj_bin = hPi0A08->GetYaxis()->FindBin(bins[i+1])-1;
	    
	    TH1* fHistBasicSame = (TH1*)hPi0A08   ->ProjectionX(Form("fHistBasicSame_No%d",i+1),sproj_bin,lproj_bin,""); 
	    TH1* fHistBasicMix  = (TH1*)hMixPi0A08->ProjectionX(Form("fHistBasicMix_No%d",i+1),sproj_bin,lproj_bin,""); 
	    fHistBasicSame->Rebin(rebin);
	    fHistBasicMix->Rebin(rebin);
	    fHistBasicSame->Sumw2();
	    fHistBasicMix->Sumw2();
	    fHistBasicSame->GetYaxis()->SetTitle(Form("dN/dM per %.0f MeV/c^{2}",fHistBasicSame->GetBinWidth(1)*1000));
//.........这里部分代码省略.........