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

/*============================================================================*/
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 GFHistManager::ColourFuncs(GFHistArray *hists) const
{
  // adjust colour of funcs to match hist, but only if exactly one function per hist
  if (!hists) return;
  for (Int_t iH = 0; iH < hists->GetEntriesFast(); ++iH) {
    TH1 *h = hists->At(iH);
    if (!h) continue;

    // look for _the_ TF1 (not > 1!)
    TF1 *func = NULL;
    TIter nextprim(h->GetListOfFunctions());
    while (TObject* next = nextprim()) {
      if (next->InheritsFrom(TF1::Class())) {
	if (func) { // there is already a TF1, so...
	  func = NULL; // remove it again...
	  break;       // ...and stop searching for  more!
	} else {
	  func = static_cast<TF1*>(next);
	}
      }
    }
// if exactly 1 found, adjust line style/colour
    if (func) {
      func->SetLineColor(h->GetLineColor());
      func->SetLineStyle(h->GetLineStyle());
    }
  } 
}

void Draw_vn( TH1D * h_vn_pt, TString MBorDtrig = "MBtrig", TString EPorSP = "EP", TString vnname = "v2", TString Ytitle = "v_{2}", int cent_low = 0, int cent_high = 100, double ptlow = 1.0, double pthigh = 35.0, TString fitoption = "poly3bkg")
{
	TCanvas * cfg_vn = new TCanvas(Form("cfg_%s_%s", MBorDtrig.Data(), vnname.Data()));

	h_vn_pt->GetXaxis()->SetRangeUser(ptlow+0.05, pthigh-0.05);
	h_vn_pt->GetYaxis()->SetRangeUser(-0.2, 0.4);
	h_vn_pt->GetXaxis()->SetTitle("D^{0} p_{T} (GeV/c)");
	h_vn_pt->GetYaxis()->SetTitle(Ytitle);

	h_vn_pt->GetYaxis()->SetTitleSize(0.05);
	h_vn_pt->GetXaxis()->SetTitleSize(0.05);

	if( EPorSP == "EP" )
	{
		h_vn_pt->SetMarkerColor(6.0);
		h_vn_pt->SetLineColor(6.0);
		h_vn_pt->SetMarkerStyle(22);
		h_vn_pt->SetMarkerSize(1.0);
	}
	else if( EPorSP == "SP" )
	{
		h_vn_pt->SetMarkerColor(4.0);
		h_vn_pt->SetLineColor(4.0);
		h_vn_pt->SetMarkerStyle(21);
		h_vn_pt->SetMarkerSize(1.3);
	}

	h_vn_pt->Draw();

	TLatex Tl; 
	Tl.SetNDC();
	Tl.SetTextAlign(12);
	Tl.SetTextSize(0.05);
	Tl.SetTextFont(42);
	Tl.DrawLatex(0.125,0.965, "#font[61]{CMS} #scale[0.8]{Preliminary}");
	Tl.DrawLatex(0.57,0.965, "#scale[0.8]{PbPb #sqrt{s_{NN}} = 5.02 TeV}");

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

	tex = new TLatex(0.20,0.78,Form("Cent. %d-%d%%", cent_low, cent_high));
	tex->SetNDC();
	tex->SetTextFont(42);
	tex->SetTextSize(0.04);
	tex->SetLineWidth(2);
	tex->Draw();

	TF1 * fun = new TF1("fun", "0.0", 0, 100);
	fun->SetLineColor(1.0);
	fun->SetLineStyle(3);
	fun->SetLineWidth(1);
	fun->Draw("same");

	cfg_vn->SaveAs(Form("Plots_vn/cfg_vnvsmass_%s_%s_cent%dto%d_%s.pdf", MBorDtrig.Data(), vnname.Data(), cent_low, cent_high, fitoption.Data()));
}

	void CreateDrawAndSaveHistogramBreitWignerFit(TH1* &histo, TString outputdir, TString outputname, bool saveoutput, bool close, double range=1){
		setPandaStyle();

		TString name = TString(histo->GetName());
		TString title = TString(histo->GetTitle());

		TCanvas * canvas = new TCanvas("c_"+name, title, 0,0,1500,1000);
		histo->Draw();

//		canvas->Update();

		TF1 * bw;

		bw = GetVoigtFit(histo, range);

		bw->SetLineColor(kRed);
		bw->SetLineStyle(7);
		bw->SetLineWidth(3);
		bw->Draw("SAME");

		PandaSmartLabel("L");

		if (saveoutput){
			canvas->Print(outputdir + "root-files/" + outputname + ".root");
			canvas->Print(outputdir + "png-files/" + outputname + ".png");
			canvas->Print(outputdir + "pdf-files/" + outputname + ".pdf");
		}

		if (close) canvas->Close();

	}

void solveNumerically()
{
  
  Double_t n0=1.35;
  Double_t m=-1*(1.78-1.35)/140;
  Double_t z0=-100;
  Double_t dMax=n0+m*z0;


  TCanvas *can = new TCanvas("can","can",600,400);
  // TF1 *lefty = new TF1("lefty",strangeLeft,0,dMax,5);
  // lefty->SetParameters(n0,m,z0,40,-30);
  // lefty->SetLineColor(8);
  // lefty->SetLineStyle(1);
  // lefty->SetNpx(1000);
  // lefty->Draw("");

  // TF1 *righty = new TF1("righty",strangeRight,0,dMax,5);
  // righty->SetParameters(n0,m,z0,40,-30);
  // righty->SetLineColor(kViolet);
  // righty->SetLineStyle(1);
  // righty->SetNpx(1000);
  // righty->Draw("same");


  TH1F*framey = can->DrawFrame(0,-50,dMax*1.1,0);
  TF1 *endy = new TF1("endy",strangeZ,0.02,dMax,4);
  endy->SetParameters(n0,m,z0,200);
  endy->SetLineColor(kViolet);
  endy->SetLineStyle(1);
  endy->SetNpx(1000);
  endy->Draw("same");

}

void problem_1() {
    //Declare variable for energy
    Float_t energy;

    //declare a histogram for the energy
    TH1D *energy_hist = new TH1D("Energy","Energy",50,100,160);
    //declare the canvas to display hist
    TCanvas *c1 = new TCanvas("c1","c1",10,10,900,900);
    TF1 *fitFcn = new TF1("fitFcn",fitFunction,90,180,6);

    //get the input file
    TFile *inputFile = new TFile("data.root");
    //get the correct tree from the input file
    //depending on the tree name change "data"
    TTree *data_tree = (TTree*)inputFile->Get("data");

    //Get the branch named E1 and put it into the varialbe energy
    data_tree->SetBranchAddress("E", &energy);

    //Get the number of events for the for loop
    int count = data_tree->GetEntries();

    for (int i = 0; i < count; i++) {
        //get the current event i
        data_tree->GetEntry(i);
        //Fill the histogram
        energy_hist->Fill(energy);
    }
    //label the axis and draw the histogram after it's filled
    energy_hist->GetXaxis()->SetTitle("E (GeV)");
    energy_hist->Draw();
    // first try without starting values for the parameters
    // this defaults to 1 for each param.
    energy_hist->Fit("fitFcn");
    // this results in an ok fit for the polynomial function however
    // the non-linear part (Lorentzian) does not respond well
    // second try: set start values for some parameters
    fitFcn->SetParameter(4,1.66); // width
    fitFcn->SetParameter(5,126.5); // peak
    energy_hist->Fit("fitFcn");

    // improve the picture:
    TF1 *backFcn = new TF1("backFcn",background,100,160,3);
    backFcn->SetLineColor(1);
    backFcn->SetLineStyle(3);
    TF1 *signalFcn = new TF1("signalFcn",gaussian,100,160,3);
    signalFcn->SetLineColor(4);
    Double_t par[6];

    // writes the fit results into the par array
    fitFcn->GetParameters(par);
    backFcn->SetParameters(par);
    backFcn->Draw("same");
    signalFcn->SetParameters(&par[4]);
    signalFcn->Draw("same");

    gStyle->SetOptFit(1111);

}

	  TF1 * getDoubleFit(TH1 * hist, bool autorange, double inner, double outer, int which=1){
		  /**
		   * @brief performes a double gaussian Fit
		   */

		  Double_t parameters[6] = {0,0,0,0,0,0};
		  Double_t parameters_new[6] = {0,0,0,0,0,0};


		  double center = hist->GetMean();

		  if (autorange){
			  double inner = fabs(center - hist->GetXaxis()->GetXmax()) /10;
			  double outer = fabs(center - hist->GetXaxis()->GetXmax()) /10*8;
		  }

		  TF1 * fit1 = new TF1("fit1", "gaus", center-inner, center+inner );

		  hist->Fit(fit1, "Q0R");
		  fit1->GetParameters(&parameters[0]);

		  TF1 * fit2 = new TF1("fit2", "gaus", center-outer, center+outer );

		  hist->Fit(fit2, "Q0R");
		  fit2->GetParameters(&parameters[3]);

		  TF1 * fitproper = new TF1("fitproper", "gaus(0)+gaus(3)", center-outer, center+outer);
		  fitproper->SetParameters(parameters);
		  fitproper->SetParName(0, "Const.(inner)");
		  fitproper->SetParName(1, "Mean (inner)");
		  fitproper->SetParName(2, "Sigma (inner)");
		  fitproper->SetParName(3, "Const.(outer)");
		  fitproper->SetParName(4, "Mean (outer)");
		  fitproper->SetParName(5, "Sigma (outer)");

		  hist->Fit(fitproper, "Q0R");
		  fitproper->GetParameters(&parameters_new[0]);


		  fitproper->SetLineColor(hist->GetLineColor());
		  fitproper->SetLineWidth(hist->GetLineWidth());
		  fitproper->SetLineStyle(2);

		  if (which==1){
		  	  TF1 * Gausinner = new TF1("Gausinner", "gaus(0)", center-inner, center+inner);
			  Gausinner->SetParameters(parameters_new);
		  	  return Gausinner;
		  }
		  else{
			  TF1 * Gausouter = new TF1("Gausouter", "gaus(3)", center-outer, center+outer);
			  Gausouter->SetParameters(parameters_new);
			  return Gausouter;
		  }

	  }

void Draw_vn(TH1D * vn_cent0to10, TH1D * vn_cent10to30, TH1D * vn_cent30to50, TH1D * vn_cent50to70, TString vnname = "v2_SP", TString Ytitle = "v_{2}{SP}")
{
	TCanvas * cfg_vn_centrality = new TCanvas(Form("cfg_%s_centrality", vnname.Data()), Form("cfg_%s_centrality", vnname.Data()));

	vn_cent10to30->SetLineColor(4.0);
	vn_cent10to30->SetMarkerColor(4.0);

	vn_cent30to50->SetLineColor(2.0);
	vn_cent30to50->SetMarkerColor(2.0);

	vn_cent50to70->SetLineColor(8.0);
	vn_cent50to70->SetMarkerColor(8.0);

	vn_cent0to10->GetYaxis()->SetTitle(Ytitle);

	vn_cent0to10->Draw();
	vn_cent10to30->Draw("same");
	vn_cent30to50->Draw("same");
	vn_cent50to70->Draw("same");

    TF1 * fun = new TF1("fun", "0.0", 0, 100);
    fun->SetLineColor(1.0);
    fun->SetLineStyle(3);
    fun->SetLineWidth(1);
    fun->Draw("same");

	TLegend * leg = new TLegend(0.45, 0.67, 0.65, 0.90);
	leg->SetTextSize(0.04);
	leg->SetTextFont(42);
	leg->SetBorderSize(0);
	leg->SetFillStyle(0);
	leg->AddEntry( vn_cent0to10, "Cent. 0-10%");
	leg->AddEntry( vn_cent10to30, "Cent. 10-30%");
	leg->AddEntry( vn_cent30to50, "Cent. 30-50%");
	leg->AddEntry( vn_cent50to70, "Cent. 50-70%");
	leg->Draw();

    TLatex Tl;
    Tl.SetNDC();
    Tl.SetTextAlign(12);
    Tl.SetTextSize(0.05);
    Tl.SetTextFont(42);
    Tl.DrawLatex(0.125,0.965, "#font[61]{CMS} #scale[0.8]{Preliminary}");
    Tl.DrawLatex(0.57,0.965, "#scale[0.8]{PbPb #sqrt{s_{NN}} = 5.02 TeV}");

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

	cfg_vn_centrality->SaveAs(Form("Plots_vn/cfg_%s_centrality.pdf", vnname.Data()));
}

TF1* plotmax(TString tagger, TString TaggerStrength, float Etamin, float Etamax, TString opt = "" , int col = 1, float lineWidth = 1, int lineStyle = 1)
{
  TF1* f = GetMistagmax(tagger,TaggerStrength,Etamin,Etamax);
  if( f != NULL )
  {
    f->SetLineColor(col);
    f->SetLineWidth(lineWidth);
    f->SetLineStyle(lineStyle);
    f->Draw(opt);
  }
  //else cout << "NULL pointer returned... Function seems not to exist" << endl;
  return f;
}

void Draw_vn_ratio( TH1D * vn_default, TH1D * vn_poly3bkg_floatwidth, TString vnname, TString method, TString MBorDtrig, int cent_low, int cent_high, float ptlow, float pthigh)
{
	TCanvas * cfg_vn_floatwidth = new TCanvas(Form("cfg_vn_floatwidth_%s_cent%dto%d", (vnname+"_"+method+"_"+MBorDtrig).Data(), cent_low, cent_high), Form("cfg_vn_floatwidth_%s_cent%dto%d", (vnname+"_"+method+"_"+MBorDtrig).Data(), cent_low, cent_high), 400, 600);
	
	TPad *pad1 = new TPad("pad1","top pad",0.0,0.4,1.0,1.0);
	pad1->SetBottomMargin(0.01);
	pad1->Draw();
	
	TPad *pad2 = new TPad("pad2","bottom pad",0.0,0.0,1.0,0.4);
	pad2->SetTopMargin(0.0);
	pad2->Draw();

	pad1->cd();
	vn_default->Draw();
	vn_poly3bkg_floatwidth->SetMarkerColor(4);
	vn_poly3bkg_floatwidth->SetLineColor(4);
	vn_poly3bkg_floatwidth->Draw("same");

    TF1 * fun = new TF1("fun", "0.0", 0, 100);
    fun->SetLineColor(1.0);
    fun->SetLineStyle(3);
    fun->SetLineWidth(1);
    fun->Draw("same");

	TLegend *  leg = new TLegend(0.45, 0.70, 0.65, 0.90);
	leg->SetTextSize(0.04);
	leg->SetTextFont(42);
	leg->SetBorderSize(0);
	leg->SetFillStyle(0);
	leg->AddEntry( (TObject*)0, Form("Cent. %d-%d%%", cent_low, cent_high) , "");
	leg->AddEntry( vn_default, "Default");
	leg->AddEntry( vn_poly3bkg_floatwidth, "Float width");
	leg->Draw();

	pad2->cd();
	TH1D * h_vn_ratio_poly3bkg_floatwidth_default = (TH1D *) vn_poly3bkg_floatwidth->Clone("h_vn_ratio_poly3bkg_floatwidth_default");

	h_vn_ratio_poly3bkg_floatwidth_default->Divide( vn_poly3bkg_floatwidth, vn_default, 1.0, 1.0, "B");

	h_vn_ratio_poly3bkg_floatwidth_default->GetYaxis()->SetRangeUser(0.7,1.3);
	h_vn_ratio_poly3bkg_floatwidth_default->GetYaxis()->SetTitle("*/Default");
	h_vn_ratio_poly3bkg_floatwidth_default->Draw();

    TF1 * fun1 = new TF1("fun1", "1.0", 0, 100);
    fun1->SetLineColor(1.0);
    fun1->SetLineStyle(3);
    fun1->SetLineWidth(1);
    fun1->Draw("same");

	cfg_vn_floatwidth->SaveAs(Form("Plots_vn/plots_PDFvariation/cfg_vn_floatwidth_%s_cent%dto%d.pdf", (vnname+"_"+method+"_"+MBorDtrig).Data(), cent_low, cent_high));
}

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 drawSupport(){

#ifdef __CINT__
    gROOT->LoadMacro("AtlasUtils.C");
#endif
    SetAtlasStyle();


    char name[200];

    //sprintf(name,"Mydecorr_0_Calib1.root");
    sprintf(name,"Mydecorr_0_Calib1_1.root"); 
    //sprintf(name,"support.root");
    fin = TFile::Open(name);

    sprintf(name, "hcentrality");
    TH1* hcent = (TH1D*)fin->Get(name);

    sprintf(name, "hntrk_Fcal_%d",0);
    hntrk_Fcal[0] = (TH2*)fin->Get(name);
    hntrk_Fcal[0]->Rebin2D(5,5);

    sprintf(name, "hntrk_Fcal_%d",1);
    hntrk_Fcal[1] = (TH2*)fin->Get(name);
    sprintf(name, "hFCalAC");
    TH2* hFCalAC = (TH2*)fin->Get(name);

    TF1 *fa = new TF1("fa","[0] + [1]*x",0,7);
    double par2[] = {-1300,1125};
    fa->SetLineColor(kRed);
    fa->SetParameters(par2);
    fa->SetLineWidth(2);
    fa->SetLineStyle(2);

    TF1 *fa2 = new TF1("fa2","[0] + [1]*x",0,7);
    double par[] = {4.1,-1};
    fa2->SetLineColor(kRed);
    fa2->SetParameters(par);
    fa2->SetLineWidth(2);
    fa2->SetLineStyle(2);

    for (int ic=0; ic<NCENT; ic++) {
        sprintf(name,"hntrkloose_ic%d", ic); hntrkloose[ic] =(TH1D*)fin->Get(name);
        cout<<name<< " "<<hntrkloose[ic]->GetBinContent(3325)<< endl;
        sprintf(name,"hntrktight_ic%d", ic); hntrktight[ic] =(TH1D*)fin->Get(name);
    }

    for (int ic=0; ic<20; ic++) {
        sprintf(name,"hntrkcut_ic%d", ic); hntrkcut[ic] =(TH1D*)fin->Get(name);
    }



       sprintf(name,"sfcalAC");
       can = SetCan2(name, 2,2);
       hFCalAC->Draw("colz");
       hFCalAC->GetXaxis()->SetTitle("#Sigma E_{T}^{FCal A side} [TeV]");
       hFCalAC->GetYaxis()->SetTitle("#Sigma E_{T}^{FCal C side} [TeV]");
       hFCalAC->Rebin2D(10,10);
       hFCalAC->GetXaxis()->SetRangeUser(0,3.5);
       hFCalAC->GetYaxis()->SetRangeUser(0,3.5);
       gPad->SetLogz();
       fa2->Draw("same");
       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);


       sprintf(name,"sfcal");
       can = SetCan2(name, 2,2);
       hntrk_Fcal[0]->Draw("colz");
       hntrk_Fcal[0]->GetXaxis()->SetTitle("#Sigma E_{T}^{FCal} [TeV]");
       hntrk_Fcal[0]->GetYaxis()->SetTitle("N_{track}^{Loose}  [p_{T}>0.5GeV]");

       gPad->SetLogz();
       fa->Draw("same");
       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);



       sprintf(name,"sCentrality");
       can = SetCan2(name, 2,2);
       hcent->Draw();

       hcent->GetXaxis()->SetTitle("Centrality %");
       hcent->GetYaxis()->SetTitle("Nevents");

       sprintf(name,"%s.pdf", can->GetName()); can->SaveAs(name);



       int Nx=5, Ny=4;
       sprintf(name,"sntrkLoose");
       can = SetCan2D2(name,Nx, Ny);
       for (int ix=0; ix<Nx; ix++) {
       for (int iy=0; iy<Ny; iy++) {
       int ic = ix + iy*Nx;
       can->cd(ic+1);
       if(ic>=NCENT) continue;
       hntrkloose[ic] ->Rebin(10);
       hntrkloose[ic] ->DrawClone();
       gPad->SetLogy();
//.........这里部分代码省略.........

void fitMKVoight(TH1 *h33 , Double_t low, Double_t high, Double_t p0, Double_t p1, Double_t p2, Double_t p3, Double_t initialPar, Double_t width, Double_t Gamma, Double_t factor, Int_t draw_opt){
  
  //double nEnt = h33->GetEntries();
  double nEnt = h33->GetMaximum();
  double in_par[8] = {100, initialPar, width , Gamma, p0, p1, p2}; //{A,mean,sigma,Gamma}
  TF1 *fitter = new TF1("fitter",myFuncBWG,low,high,8);
  fitter->SetParameters(&in_par[0]);
  fitter->SetParLimits(0,10,nEnt); fitter->SetParLimits(1,initialPar-width,initialPar+width);
  fitter->SetLineColor(kBlue);
  h33->Fit("fitter","REM+","same");
  //h33->Draw("E");
  TF1 *backFcn = new TF1("backFcn", "pol2",low,high);
  TF1 *signalFcn = new TF1("signalFcn", myVoight,low,high,4);
  signalFcn->SetLineColor(kBlue);
  signalFcn->SetLineWidth(2);
  Double_t par[8];
  fitter->GetParameters(par);
  backFcn->SetParameters(&par[4]);
  backFcn->SetLineStyle(2);
  backFcn->SetLineColor(6);
  backFcn->SetLineWidth(1);
  
  
  signalFcn->SetParameters(par);
  signalFcn->SetLineStyle(2);
  signalFcn->SetLineColor(4);
  signalFcn->SetLineWidth(1);
  // backFcn->Draw("same");
  //signalFcn->Draw("same");
  //Double_t Intg = abs(signalFcn->Integral(par[1]-factor*par[2],par[1]+factor*par[2]));
  ROOT::Math::GSLIntegrator ig(1.E-6,1.E-6,1000);
  ROOT::Math::WrappedTF1 wf(*fitter);
  ig.SetFunction(wf);
  double Intg = ig.Integral(par[1]-factor*par[2],par[1]+factor*par[2]);
  Double_t Intb = abs(backFcn->Integral(par[1]-factor*par[2],par[1]+factor*par[2]));
  
  
  Double_t binw = h33->GetBinWidth(1);
  Int_t yield = Intg/binw;
  Int_t bckgd = Intb/binw;
  //Double_t ratio = double(yield)/TMath::Sqrt(double(bckgd));
  Double_t ratio = double(yield)/double(yield+bckgd);
  
  cout << yield << "\t" << ratio << endl;
  TAxis *x=h33->GetXaxis();
  TAxis *y=h33->GetYaxis();
  
  Double_t startx=x->GetXmin()+0.45*(x->GetXmax()-x->GetXmin());
  /*
   Double_t starty0=0.5*h33->GetMaximum();
   Double_t starty1=0.56*h33->GetMaximum();
   Double_t starty2=0.62*h33->GetMaximum();
   Double_t starty3=0.68*h33->GetMaximum();
   Double_t starty4=0.74*h33->GetMaximum();
   Double_t starty5=0.8*h33->GetMaximum();
   */
  Double_t starty0=0.35*h33->GetMaximum();
  Double_t starty1=0.43*h33->GetMaximum();
  Double_t starty2=0.49*h33->GetMaximum();
  Double_t starty3=0.56*h33->GetMaximum();
  Double_t starty4=0.63*h33->GetMaximum();
  Double_t starty5=0.7*h33->GetMaximum();
  
  
  
  double meanError = fitter->GetParError(1);
  double sigmaError = fitter->GetParError(2);
  if (draw_opt ==1) {
    
    TLatex *sum = new TLatex(startx*0.93, starty5,Form("Yield: %i",yield));
    TLatex *sum12 = new TLatex(startx*0.93, starty4,Form("Background: %i",bckgd));
    TLatex *sum0=new TLatex(startx*0.93, starty3,Form("Range: #pm %2.1f #sigma",factor));
    TLatex *sum2=new TLatex(startx*0.93, starty2,Form("Mean: %4.4f #pm %.4f GeV",par[1], meanError));
    TLatex *sum3=new TLatex(startx*0.93, starty1,Form("#sigma: %5.4f #pm %.4f GeV",par[2], sigmaError));
    //TLatex *ra = new TLatex(startx*0.93, starty0,Form("#frac{S}{#sqrt{B}}= %.1f", ratio));
    TLatex *ra = new TLatex(startx*0.93, starty0,Form("#frac{S}{S+B} = %1.4f", ratio));
    
    
    sum->SetTextSize(0.04);
    sum->SetTextColor(kBlue);//2
    sum->Draw("same");
    sum12->SetTextSize(0.04);
    sum12->SetTextColor(kBlue);//6
    sum12->Draw("same");
    
    ra->SetTextSize(0.04);
    ra->SetTextColor(kBlue);//was black before;
    ra->Draw("same");
    
    
    sum0->SetTextSize(0.04);
    sum0->SetTextColor(kBlue);//2
    sum0->Draw("same");
    sum2->SetTextSize(0.04);
    sum2->SetTextColor(kBlue);//4
    sum2->Draw("same");
    sum3->SetTextSize(0.04);
    sum3->SetTextColor(kBlue);//4
    sum3->Draw("same");
    
//.........这里部分代码省略.........

// Function for the computation of channeling efficiency at various incoming angle
Int_t AnalyseChannelingEfficiency(TTree *fTree,Float_t fChannelingMinimum = 35., Float_t fChannelingMaximum = 70.){
    //**//Channeling Gaussian Fit Function
    TF1 *vChanneling = new TF1("vChanneling","gaus",fChannelingMinimum,fChannelingMaximum);
    vChanneling->SetParNames("Const","Mean","Sigma");
    vChanneling->SetLineColor(4);
    vChanneling->SetLineStyle(2);

    TH2D *hChannelingPlot = new TH2D("hChannelingPlot","Deflection Angle vs. Incoming Angle;Horizontal Incoming Angle [#murad];Horizontal Deflection Angle [#murad]",21,-10.5,10.5,256,-127.5,128.5);
    
    TH1F *hChannelingEfficiency = new TH1F("hChannelingEfficiency","G4Channeling;Horizontal Incoming Angle [#murad];Efficiency [%]",21,-10.5,10.5);

    fTree->Draw("-(angXout-angXin):-angXin>>hChannelingPlot");
    
    Double_t vNormalizationToAmorphous = 0.965; // Normalization for channeling efficiency, see PRSTAB 11, 063501 (2008)
    
    for(int i=2;i<=21;i++){
        TH1D* h1 = hChannelingPlot->ProjectionY("h1",i,i);
        h1->Fit(vChanneling,"QR");
        Double_t *vChannelingParameters;
        vChannelingParameters = vChanneling->GetParameters();
        hChannelingEfficiency->SetBinContent(i,ComputeEfficiency(h1,vChannelingParameters)/vNormalizationToAmorphous);
        h1->Delete();
    }
    hChannelingEfficiency->SetLineColor(3);
    hChannelingEfficiency->SetLineStyle(4);
    hChannelingEfficiency->SetMarkerColor(3);
    hChannelingEfficiency->SetFillStyle(0);
    hChannelingEfficiency->SetMarkerStyle(20);
    hChannelingEfficiency->Draw("PL");

    TGraph* gRoughExperimentalData = new TGraph(11);
    gRoughExperimentalData->SetPoint(	0	,	-10	,	20	);
    gRoughExperimentalData->SetPoint(	1	,	-8	,	38	);
    gRoughExperimentalData->SetPoint(	2	,	-6	,	56	);
    gRoughExperimentalData->SetPoint(	3	,	-4	,	72	);
    gRoughExperimentalData->SetPoint(	4	,	-2	,	80	);
    gRoughExperimentalData->SetPoint(	5	,	0	,	84	);
    gRoughExperimentalData->SetPoint(	6	,	2	,	82	);
    gRoughExperimentalData->SetPoint(	7	,	4	,	78	);
    gRoughExperimentalData->SetPoint(	8	,	6	,	66	);
    gRoughExperimentalData->SetPoint(	9	,	8	,	52	);
    gRoughExperimentalData->SetPoint(	10	,	10	,	37	);

    gRoughExperimentalData->SetLineColor(4);
    gRoughExperimentalData->SetLineStyle(3);
    gRoughExperimentalData->SetFillStyle(0);
    gRoughExperimentalData->SetFillColor(0);
    gRoughExperimentalData->SetMarkerColor(4);
    gRoughExperimentalData->SetMarkerStyle(21);
    gRoughExperimentalData->SetTitle("Phys. Lett. B 680, 129");

    gRoughExperimentalData->Draw("sameCP");

    TLegend *aLegend = new TLegend(0.30,0.15,0.55,0.3);
    aLegend->AddEntry(hChannelingEfficiency);
    aLegend->AddEntry(gRoughExperimentalData);
    aLegend->SetFillStyle(0);
    aLegend->SetLineColor(0);
    aLegend->Draw();

    return 0;
}

TF1* fitMass(TH1D* hData, TH1D* hMCSignal, TH1D* hMCSwapped)
{
	Double_t setparam0=100.;
	Double_t setparam1=1.865;
	Double_t setparam2=0.03;
	Double_t setparam10=0.005;
	Double_t setparam8=0.1;
	Double_t setparam9=0.1;
	Double_t fixparam1=1.865;
	Double_t minhisto=1.7;
	Double_t maxhisto=2.0;

	TF1* f = new TF1("fMass","[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.1415927)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.1415927)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.1415927)*[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.005,0.05);
	f->SetParLimits(2,0.01,0.1);
	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);

	hMCSignal->Fit("fMass","q","",minhisto,maxhisto);
	hMCSignal->Fit("fMass","q","",minhisto,maxhisto);
	f->ReleaseParameter(1);
	hMCSignal->Fit("fMass","L q","",minhisto,maxhisto);
	hMCSignal->Fit("fMass","L q","",minhisto,maxhisto);
	hMCSignal->Fit("fMass","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("fMass","L q","",minhisto,maxhisto);
	hMCSwapped->Fit("fMass","L q","",minhisto,maxhisto);
	hMCSwapped->Fit("fMass","L q","",minhisto,maxhisto);
	hMCSwapped->Fit("fMass","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);

	hData->Fit("fMass","q","",minhisto,maxhisto);
	hData->Fit("fMass","q","",minhisto,maxhisto);
	f->ReleaseParameter(1);
	f->SetParLimits(1,1.86,1.87);
	f->ReleaseParameter(11);
	f->SetParLimits(11,-0.2,0.2);
	hData->Fit("fMass","L q","",minhisto,maxhisto);
	hData->Fit("fMass","L q","",minhisto,maxhisto);
	hData->Fit("fMass","L q","",minhisto,maxhisto);
	hData->Fit("fMass","L m","",minhisto,maxhisto);

	TF1* background = new TF1("fBackground","[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("fSignal","[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.1415927)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.1415927)*[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("fBackground","[0]*(1-[2])*Gaus(x,[1],[3]*(1+[4]))/(sqrt(2*3.1415927)*[3]*(1+[4]))");
	massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8),f->GetParameter(11));
//.........这里部分代码省略.........