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

void plot_diff_xsec()
{
  TF1* fXsec = new TF1("fXsec", "diff_xsec(x)", 1, 0);
  fXsec->SetLineWidth(3);
  fXsec->SetTitle("#nu_{#mu}-e differential cross section (E_{#nu} = 2 GeV)");
  fXsec->Draw();
  fXsec->GetXaxis()->SetTitle("cos#theta");
  fXsec->GetXaxis()->CenterTitle();
  fXsec->GetYaxis()->SetTitle("d#sigma/dcos#theta (cm^{-2})");
  fXsec->GetYaxis()->CenterTitle();
  fXsec->GetYaxis()->SetTitleOffset(1.4);
  gPad->SetLogy(1);
  
  double totXsec = fXsec->Integral(0,1);
  for(int i = 0; i < 100; i++)
  {
    cout << fXsec->Integral(0.01*i,1)/totXsec << endl;
  }
  cout << totXsec << " " << fXsec->Integral(.99,1) << endl;
  
  ROOT::Math::GSLIntegrator ig(1.E-6,1.E-6,1000); 
  ROOT::Math::WrappedTF1 wf(*fXsec);
  ig.SetFunction(wf);
  cout << ig.Integral(.9,1) << endl;
  
  //~ TCanvas* c2 = new TCanvas("c2");
  //~ TF1* fCXsec = new TF1("fCXsec", "cumul_xsec(x)", 0, 90);
  //~ fCXsec->Draw();
  //~ cout << cumul_xsec(45) << endl;
}

void setGraphOptions(TF1 &f, 
		     const char *xname, 
		     const char *yname)
{
  TAxis::TAxis *axis1;
  TAxis::TAxis *axis2;
  axis1 = f.GetXaxis();
  axis2 = f.GetYaxis();
  axis1->SetTitle(xname);
  axis2->SetTitle(yname);
  setAxisOptions(axis1);
  setAxisOptions(axis2);
}

void wilenbrock_pseudoscalar()
{

  setTDRStyle();

  TCanvas* c = new TCanvas();
  TF1* func = new TF1("Signal + Interferences",Somme,350,900,1);
  func->SetTitle("Scalar : Signal + interferences");
  func->SetMinimum(-9e-9);
  //func->GetXaxis()->SetTitle("m_{t#bar{t}}");
  double m=400;
  func->SetParameter(0,m);
  func->SetLineColor(TColor::GetColor("#542437"));
  func->SetLineWidth(2);
  func->SetNpx(500);
  func->Draw();
  func->GetXaxis()->SetTitle("m_{t#bar{t}}");
  func->GetYaxis()->SetTitle("#sigma(#hat{s}) - #sigma_{QCD}(#hat{s})");
  TF1* func2 = new TF1("Signal + Interferences",Somme,350,900,1);
  m=500;
  func2->SetParameter(0,m);
  func2->SetLineColor(TColor::GetColor("#C02942"));
  func2->SetLineWidth(2);
  func2->SetNpx(500);
  func2->Draw("SAME");
  TF1* func3 = new TF1("Signal + Interferences",Somme,350,900,1);
  m=600;
  func3->SetParameter(0,m);
  func3->SetLineColor(TColor::GetColor("#53777A"));
  func3->SetLineWidth(2);
  func3->SetNpx(500);
  func3->Draw("SAME");
  TF1* func4 = new TF1("Signal + Interferences",Somme,350,900,1);
  m=700;
  func4->SetParameter(0,m);
  func4->SetLineColor(TColor::GetColor("#D95B43"));
  func4->SetLineWidth(2);
  func4->SetNpx(500);
  func4->Draw("SAME");
  TF1* func5 = new TF1("Signal + Interferences",Somme,350,900,1);
  m=800;
  func5->SetParameter(0,m);
  func5->SetLineColor(TColor::GetColor("#ECD078"));
  func5->SetLineWidth(2);
  func5->SetNpx(500);
  func5->Draw("SAME");

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

void FFT()
{

//This tutorial illustrates the Fast Fourier Transforms interface in ROOT.
//FFT transform types provided in ROOT:
// - "C2CFORWARD" - a complex input/output discrete Fourier transform (DFT) 
//                  in one or more dimensions, -1 in the exponent
// - "C2CBACKWARD"- a complex input/output discrete Fourier transform (DFT) 
//                  in one or more dimensions, +1 in the exponent
// - "R2C"        - a real-input/complex-output discrete Fourier transform (DFT)
//                  in one or more dimensions,
// - "C2R"        - inverse transforms to "R2C", taking complex input 
//                  (storing the non-redundant half of a logically Hermitian array) 
//                  to real output
// - "R2HC"       - a real-input DFT with output in ¡Èhalfcomplex¡É format, 
//                  i.e. real and imaginary parts for a transform of size n stored as
//                  r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
// - "HC2R"       - computes the reverse of FFTW_R2HC, above
// - "DHT"        - computes a discrete Hartley transform
// Sine/cosine transforms:
//  DCT-I  (REDFT00 in FFTW3 notation)
//  DCT-II (REDFT10 in FFTW3 notation)
//  DCT-III(REDFT01 in FFTW3 notation)
//  DCT-IV (REDFT11 in FFTW3 notation)
//  DST-I  (RODFT00 in FFTW3 notation)
//  DST-II (RODFT10 in FFTW3 notation)
//  DST-III(RODFT01 in FFTW3 notation)
//  DST-IV (RODFT11 in FFTW3 notation)
//First part of the tutorial shows how to transform the histograms
//Second part shows how to transform the data arrays directly
//Authors: Anna Kreshuk and Jens Hoffmann


//********* Histograms ********//


   //prepare the canvas for drawing
   TCanvas *myc = new TCanvas("myc", "Fast Fourier Transform", 800, 600);
   myc->SetFillColor(45);
   TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.67,0.49,0.99);
   TPad *c1_2 = new TPad("c1_2", "c1_2",0.51,0.67,0.99,0.99);
   TPad *c1_3 = new TPad("c1_3", "c1_3",0.01,0.34,0.49,0.65);
   TPad *c1_4 = new TPad("c1_4", "c1_4",0.51,0.34,0.99,0.65);
   TPad *c1_5 = new TPad("c1_5", "c1_5",0.01,0.01,0.49,0.32);
   TPad *c1_6 = new TPad("c1_6", "c1_6",0.51,0.01,0.99,0.32);
   c1_1->Draw();
   c1_2->Draw();
   c1_3->Draw();
   c1_4->Draw();
   c1_5->Draw();
   c1_6->Draw();
   c1_1->SetFillColor(30);
   c1_1->SetFrameFillColor(42);
   c1_2->SetFillColor(30);
   c1_2->SetFrameFillColor(42);
   c1_3->SetFillColor(30);
   c1_3->SetFrameFillColor(42);
   c1_4->SetFillColor(30);
   c1_4->SetFrameFillColor(42);
   c1_5->SetFillColor(30);
   c1_5->SetFrameFillColor(42);
   c1_6->SetFillColor(30);
   c1_6->SetFrameFillColor(42);
   
   c1_1->cd();
   TH1::AddDirectory(kFALSE);
     
   //A function to sample
   TF1 *fsin = new TF1("fsin", "sin(x)*sin(x)/(x*x)", 0, 4*TMath::Pi());
   fsin->Draw();
   
   Int_t n=25;
   TH1D *hsin = new TH1D("hsin", "hsin", n+1, 0, 4*TMath::Pi());
   Double_t x;
   
   //Fill the histogram with function values
   for (Int_t i=0; i<=n; i++){
      x = (Double_t(i)/n)*(4*TMath::Pi());
      hsin->SetBinContent(i+1, fsin->Eval(x));
   }
   hsin->Draw("same");
   fsin->GetXaxis()->SetLabelSize(0.05);
   fsin->GetYaxis()->SetLabelSize(0.05);
   
   c1_2->cd();
   //Compute the transform and look at the magnitude of the output
   TH1 *hm =0;
   TVirtualFFT::SetTransform(0);
   hm = hsin->FFT(hm, "MAG");
   hm->SetTitle("Magnitude of the 1st transform");
   hm->Draw();
   //NOTE: for "real" frequencies you have to divide the x-axes range with the range of your function 
   //(in this case 4*Pi); y-axes has to be rescaled by a factor of 1/SQRT(n) to be right: this is not done automatically!
   
   hm->SetStats(kFALSE);
   hm->GetXaxis()->SetLabelSize(0.05);
   hm->GetYaxis()->SetLabelSize(0.05);
   c1_3->cd();   
   //Look at the phase of the output   
   TH1 *hp = 0;
//.........这里部分代码省略.........

void betaMS(Double_t z = 1000, const TString &opt="pdf") {

  gROOT->Reset();
  PGlobals::Initialize();
  
  // gStyle->SetPadTopMargin(0.10);
  // gStyle->SetPadRightMargin(0.05);
  // gStyle->SetPadBottomMargin(0.15);
  // gStyle->SetPadLeftMargin(0.15);

  Double_t E0  = 1000;  // MeV
  Double_t gamma = E0/0.511; 
  Double_t X0  = 35.0E4; // um (Berylium)
  Double_t sx0 = 10.0;   // um
  Double_t emitxn  = 1.00; // um
  Double_t emitx   = emitxn/gamma; // um
  Double_t sd0 = 1000. * (emitx/sx0); // mrad
  Double_t kp  = 0.0595; // um^-1  for n0 = 10^17 cm^-3
  Double_t kbeta  = kp/TMath::Sqrt(2 * gamma);
  Double_t betamatch = (1/kbeta) / 1000.; // mm

  TF1 *ftheta = new TF1("ftheta",theta,-0.01*z,z,4);
  ftheta->SetParameters(E0,X0,sx0,sd0);
  ftheta->SetParNames("Energy","X0","sx0","sd0");
  ftheta->SetTitle("");
  ftheta->GetXaxis()->SetTitle("#Deltaz [#mum]");
  ftheta->GetYaxis()->SetTitle("#Theta_{0} [mrad]");

  TF1 *fsigmax = new TF1("fsigmax",Sigmax,-0.01*z,z,4);
  fsigmax->SetParameters(E0,X0,sx0,sd0);
  fsigmax->SetParNames("Energy","X0","sx0","sd0");
  fsigmax->SetTitle("");
  fsigmax->GetXaxis()->SetTitle("#Deltaz [#mum]");
  fsigmax->GetYaxis()->SetTitle("#sigma_{x} [#mum]");

  TF1 *fsigmaxp = new TF1("fsigmaxp",Sigmaxp,-0.01*z,z,4);
  fsigmaxp->SetParameters(E0,X0,sx0,sd0);
  fsigmaxp->SetParNames("Energy","X0","sx0","sd0");
  fsigmaxp->SetTitle("");
  fsigmaxp->GetXaxis()->SetTitle("#Deltaz [#mum]");
  fsigmaxp->GetYaxis()->SetTitle("#sigma_{x'} [mrad]");

  TF1 *femit = new TF1("femit",emittance,-0.01*z,z,4);
  femit->SetParameters(E0,X0,sx0,sd0);
  femit->SetParNames("Energy","X0","sx0","sd0");
  femit->SetTitle("");
  femit->GetXaxis()->SetTitle("#Deltaz [#mum]");
  femit->GetYaxis()->SetTitle("#epsilon_{x,n} [#mum]");

  TF1 *fbeta = new TF1("fbeta",beta,-0.01*z,z,4);
  fbeta->SetParameters(E0,X0,sx0,sd0);
  fbeta->SetParNames("Energy","X0","sx0","sd0");
  fbeta->SetTitle("");
  fbeta->GetXaxis()->SetTitle("#Deltaz [#mum]");
  fbeta->GetYaxis()->SetTitle("#beta [mm]");
  
  TF1 *femitfin = new TF1("femitfin",finalemittance,-0.01*z,z,5);
  femitfin->SetParameters(E0,X0,sx0,sd0,betamatch);
  femitfin->SetParNames("Energy","X0","sx0","sd0","betam");
  femitfin->SetTitle("");
  femitfin->GetXaxis()->SetTitle("#Deltaz [#mum]");
  femitfin->GetYaxis()->SetTitle("#epsilon_{x,n,F} [#mum]");

  const UInt_t NPad = 6;
  TCanvas *C = new TCanvas("C","",640,1024);
  // Setup Pad layout:
  TPad **pad = new TPad*[NPad];
  Float_t bMargin = 0.10;
  Float_t tMargin = 0.04;
  Float_t lMargin = 0.20;
  Float_t rMargin = 0.08;
  Float_t vSpacing = 0.017;
  PGlobals::CanvasPartition(C,NPad,lMargin,rMargin,bMargin,tMargin,vSpacing);
 
  // Define the frames for plotting
  Int_t fonttype = 43;
  Int_t fontsize = 22;
  Int_t tfontsize = 26;
  Float_t txoffset = 5.0;
  Float_t lxoffset = 0.02;
  Float_t tyoffset = 3.0;
  Float_t lyoffset = 0.01;
  Float_t tylength = 0.02;
  Float_t txlength = 0.04;

  TF1 *func[NPad];
  func[5] = ftheta;
  func[4] = fsigmax;
  func[3] = fsigmaxp;
  func[2] = femit;
  func[1] = fbeta;
  func[0] = femitfin;

  for(Int_t i=0;i<NPad;i++) {
    char name[16];
    sprintf(name,"pad_%i",i);
    pad[i] = (TPad*) gROOT->FindObject(name);
    pad[i]->SetFrameLineWidth(2);  
    pad[i]->SetTickx(1);
    pad[i]->SetTicky(1);
//.........这里部分代码省略.........

void FFT()
{
   // Histograms
   // =========
   //prepare the canvas for drawing
   TCanvas *myc = new TCanvas("myc", "Fast Fourier Transform", 800, 600);
   myc->SetFillColor(45);
   TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.67,0.49,0.99);
   TPad *c1_2 = new TPad("c1_2", "c1_2",0.51,0.67,0.99,0.99);
   TPad *c1_3 = new TPad("c1_3", "c1_3",0.01,0.34,0.49,0.65);
   TPad *c1_4 = new TPad("c1_4", "c1_4",0.51,0.34,0.99,0.65);
   TPad *c1_5 = new TPad("c1_5", "c1_5",0.01,0.01,0.49,0.32);
   TPad *c1_6 = new TPad("c1_6", "c1_6",0.51,0.01,0.99,0.32);
   c1_1->Draw();
   c1_2->Draw();
   c1_3->Draw();
   c1_4->Draw();
   c1_5->Draw();
   c1_6->Draw();
   c1_1->SetFillColor(30);
   c1_1->SetFrameFillColor(42);
   c1_2->SetFillColor(30);
   c1_2->SetFrameFillColor(42);
   c1_3->SetFillColor(30);
   c1_3->SetFrameFillColor(42);
   c1_4->SetFillColor(30);
   c1_4->SetFrameFillColor(42);
   c1_5->SetFillColor(30);
   c1_5->SetFrameFillColor(42);
   c1_6->SetFillColor(30);
   c1_6->SetFrameFillColor(42);

   c1_1->cd();
   TH1::AddDirectory(kFALSE);

   //A function to sample
   TF1 *fsin = new TF1("fsin", "sin(x)+sin(2*x)+sin(0.5*x)+1", 0, 4*TMath::Pi());
   fsin->Draw();

   Int_t n=25;
   TH1D *hsin = new TH1D("hsin", "hsin", n+1, 0, 4*TMath::Pi());
   Double_t x;

   //Fill the histogram with function values
   for (Int_t i=0; i<=n; i++){
      x = (Double_t(i)/n)*(4*TMath::Pi());
      hsin->SetBinContent(i+1, fsin->Eval(x));
   }
   hsin->Draw("same");
   fsin->GetXaxis()->SetLabelSize(0.05);
   fsin->GetYaxis()->SetLabelSize(0.05);

   c1_2->cd();
   //Compute the transform and look at the magnitude of the output
   TH1 *hm =0;
   TVirtualFFT::SetTransform(0);
   hm = hsin->FFT(hm, "MAG");
   hm->SetTitle("Magnitude of the 1st transform");
   hm->Draw();
   //NOTE: for "real" frequencies you have to divide the x-axes range with the range of your function
   //(in this case 4*Pi); y-axes has to be rescaled by a factor of 1/SQRT(n) to be right: this is not done automatically!

   hm->SetStats(kFALSE);
   hm->GetXaxis()->SetLabelSize(0.05);
   hm->GetYaxis()->SetLabelSize(0.05);
   c1_3->cd();
   //Look at the phase of the output
   TH1 *hp = 0;
   hp = hsin->FFT(hp, "PH");
   hp->SetTitle("Phase of the 1st transform");
   hp->Draw();
   hp->SetStats(kFALSE);
   hp->GetXaxis()->SetLabelSize(0.05);
   hp->GetYaxis()->SetLabelSize(0.05);

   //Look at the DC component and the Nyquist harmonic:
   Double_t re, im;
   //That's the way to get the current transform object:
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
   c1_4->cd();
   //Use the following method to get just one point of the output
   fft->GetPointComplex(0, re, im);
   printf("1st transform: DC component: %f\n", re);
   fft->GetPointComplex(n/2+1, re, im);
   printf("1st transform: Nyquist harmonic: %f\n", re);

   //Use the following method to get the full output:
   Double_t *re_full = new Double_t[n];
   Double_t *im_full = new Double_t[n];
   fft->GetPointsComplex(re_full,im_full);

   //Now let's make a backward transform:
   TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &n, "C2R M K");
   fft_back->SetPointsComplex(re_full,im_full);
   fft_back->Transform();
   TH1 *hb = 0;
   //Let's look at the output
   hb = TH1::TransformHisto(fft_back,hb,"Re");
   hb->SetTitle("The backward transform result");
   hb->Draw();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      gPhase->SetPoint(i,zp,(phase-zg)*kp);
    else
      gPhase->SetPoint(i,zp,(phase-zg));

    v = PFunc::PhaseVelocity2(zp,phase2,gamma,nb,n0);
    phase2 += (v - 1) * Dz;
    // cout << " z = " << zp << "  phase = " << phase << endl;
    if(opt.Contains("units"))
      gPhase2->SetPoint(i,zp,(phase2-zg)*kp);
    else
      gPhase2->SetPoint(i,zp,(phase2-zg));
  }
  
  TCanvas *C = new TCanvas("C","Wake phase velocity",1000,750);
  C->Divide(2,2);

  TLegend *Leg = new TLegend(0.6,0.20,0.85,0.35);
  PlasmaGlob::SetPaveStyle(Leg);
  Leg->SetTextAlign(22);
  Leg->SetTextColor(kGray+3);
  Leg->SetLineColor(1);
  Leg->SetBorderSize(1);
  Leg->SetFillColor(0);
  Leg->SetFillStyle(1001);
  //Leg-> SetNColumns(2);
  Leg->AddEntry(fPhaseVsZ,"PRL 107,145002","L");
  Leg->AddEntry(fPhaseVsZ2,"PRL 107,145003","L");
  Leg->SetTextColor(kGray+3);
  

  C->cd(1);
  fPhaseVsZ->GetYaxis()->SetTitle("(v_{p} - c)/c");
  if(opt.Contains("units")) 
    fPhaseVsZ->GetXaxis()->SetTitle("z [m]"); 
  else
    fPhaseVsZ->GetXaxis()->SetTitle("z [c/#omega_{p}]"); 
  fPhaseVsZ->GetYaxis()->SetRangeUser(-3e-4,0.);
  fPhaseVsZ->GetXaxis()->CenterTitle();
  fPhaseVsZ->GetYaxis()->CenterTitle();
  fPhaseVsZ->SetLineWidth(2);
  fPhaseVsZ->Draw("C");
  fPhaseVsZ2->SetLineWidth(2);
  fPhaseVsZ2->SetLineStyle(2);
  fPhaseVsZ2->Draw("C same");
  textZetag->Draw();
  Leg->Draw();

  C->cd(2);
  fPhaseVsZg->GetYaxis()->SetTitle("(v_{p} - c)/c");
  if(opt.Contains("units"))
    fPhaseVsZg->GetXaxis()->SetTitle("#zeta [m]"); 
  else
    fPhaseVsZg->GetXaxis()->SetTitle("#zeta [c/#omega_{p}]"); 
  fPhaseVsZg->GetYaxis()->SetRangeUser(-3e-4,0.);
  fPhaseVsZg->GetYaxis()->SetNdivisions(505);
  fPhaseVsZg->GetXaxis()->CenterTitle();
  fPhaseVsZg->SetLineWidth(2);
  fPhaseVsZg->GetYaxis()->CenterTitle();
  fPhaseVsZg->Draw("C");
  fPhaseVsZg2->SetLineWidth(2);
  fPhaseVsZg2->SetLineStyle(2);
  fPhaseVsZg2->Draw("C same");
  textZeta->Draw();
  
  C->cd(3);
  gPhase->GetYaxis()->SetTitle("#Delta#zeta [c/#omega_{p}]");

void fitBvar(TString collsyst="PbPb", TString inputfile ="", TString npfile="ROOTfiles/NPFitPbPb.root", float centMin=0, float centMax=100, TString outputfile="outHisto")
{
  collisionsystem = collsyst;
  infname = outputfile;
  centmin = centMin;
  centmax = centMax;
  if(collsyst != "PbPb") isPbPb = false;

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

  TF1* fit (float ptmin, float ptmax, int s, int b, int widVar); // widVar=0,1

  TCanvas* c = new TCanvas("c","",600,600);

  TF1* bmass = new TF1("bmass","[0]",7,50);
  bmass->SetTitle(";B^{+} p_{T} (GeV/c);fraction of default yield");
  bmass->SetMinimum(0.9);
  bmass->SetMaximum(1.1);
  if(isPbPb)
    {
    bmass->SetMinimum(0.9);
    bmass->SetMaximum(1.1);
    }
  bmass->SetParameter(0,1);
  bmass->SetLineWidth(1);
  bmass->SetLineColor(kRed);
  bmass->SetLineStyle(2);
  bmass->GetXaxis()->SetTitleOffset(1.3);
  bmass->GetYaxis()->SetTitleOffset(1.8);
  bmass->GetXaxis()->SetLabelOffset(0.007);
  bmass->GetYaxis()->SetLabelOffset(0.007);
  bmass->GetXaxis()->SetTitleSize(0.045);
  bmass->GetYaxis()->SetTitleSize(0.045);
  bmass->GetXaxis()->SetTitleFont(42);
  bmass->GetYaxis()->SetTitleFont(42);
  bmass->GetXaxis()->SetLabelFont(42);
  bmass->GetYaxis()->SetLabelFont(42);
  bmass->GetXaxis()->SetLabelSize(0.04);
  bmass->GetYaxis()->SetLabelSize(0.04);
  bmass->Draw();

  TLatex* Title = new TLatex(0.1,0.94, Form("Fit Variation for %s",collisionsystem.Data()));
  Title->SetNDC();
  Title->SetTextAlign(12);
  Title->SetTextSize(0.04);
  Title->SetTextFont(42);
  Title->Draw("Same");

  TLegend* leg = new TLegend(0.2,0.67,0.4,0.87,NULL,"brNDC");
  leg->SetBorderSize(0);
  leg->SetTextSize(0.04);
  leg->SetTextFont(42);
  leg->SetFillStyle(0);

    getNPFnPar(npfile, NPpar);
    std::cout<<"NP parameter 0: "<<NPpar[0]<<std::endl;
    std::cout<<"NP parameter 1: "<<NPpar[1]<<std::endl;

    TString inputf;
    if(nBins == 1) inputf = Form("%s_integrated.root",inputfile.Data());
    else inputf = Form("%s.root",inputfile.Data());

    TFile* data = new TFile(inputf.Data());  
    TH1D* h_def = (TH1D*)data->Get("hPt"); 
    TH1D* hwidvar;
    int cnt =0;

    for(int s=0;s<nsig;s++)
    {
	bkgmax[0]=sigmax[s]+2;

	for(int b=0;b<nbkg;b++)
	{

	  if(s==0 || b==0)
	  {
	      TH1D* hvar = new TH1D(Form("h_%s",bkgname[b].Data()),"",nBins,ptBins);
	      if(s==0 && b==0) hwidvar = new TH1D("h_widvar","",nBins,ptBins);

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

		  double def_y = h_def->GetBinContent(i+1);
		  double def_err = h_def->GetBinError(i+1);
		 
		  TF1* f = fit(ptBins[i],ptBins[i+1],s,b,0);
		  cout<<"YIELD / YIELDERR:  "<< yield <<" "<< yieldErr << endl;
		  cout << "DEF YIELD:  "<< def_y*(ptBins[i+1]-ptBins[i]) <<endl;
		  double y = yield/(ptBins[i+1]-ptBins[i]);
		  double err = yieldErr/(ptBins[i+1]-ptBins[i]);

		  double y_fr = y/def_y;
		  cout << "YIELD FRACTION: " << y_fr << endl;
		  double err_fr = sqrt(pow(err/def_y,2)+pow(def_err*y/(def_y*def_y),2));
//.........这里部分代码省略.........

void rsLimit7TeV(){

  setTDRStyle();

//=========Macro generated from canvas: cLimit/Limit
//=========  (Mon Feb 22 22:44:48 2010) by ROOT version5.18/00a
//   TCanvas *cLimit = new TCanvas("cLimit", "Limit",450,40,800,550);
//  TCanvas *cLimit = new TCanvas("cLimit", "Limit",100,122,600,600);
  cLimit = new TCanvas("cLimit","cLimit",800,600);

   gStyle->SetOptStat(0);
   //   cLimit->Range(595.5973,-0.03483694,1345.283,0.2539571);
   cLimit->SetFillColor(0);
   cLimit->SetBorderMode(0);
   cLimit->SetBorderSize(2);
   cLimit->SetLeftMargin(0.139262);
   cLimit->SetRightMargin(0.0604027);
   cLimit->SetTopMargin(0.0804196);
   cLimit->SetBottomMargin(0.120629);
   cLimit->SetFrameBorderMode(0);
   cLimit->SetFrameBorderMode(0);
   
   TGraph *graph = new TGraph(11);
   graph->SetName("Graph");
   graph->SetTitle("");
   graph->SetFillColor(1);


   Int_t ci;   // for color index setting
   ci = TColor::GetColor("#ff0000");
   graph->SetLineColor(ci);
   graph->SetLineWidth(3);

   ci = TColor::GetColor("#ff0000");
   graph->SetMarkerColor(ci);
   graph->SetMarkerStyle(20);
   graph->SetMarkerSize(1.0);
   graph->SetPoint(0,  863,0.01);
   graph->SetPoint(1, 1132,0.02);
   graph->SetPoint(2, 1274,0.03);
   graph->SetPoint(3, 1395,0.04);
   graph->SetPoint(4, 1503,0.05);
   graph->SetPoint(5, 1597,0.06);
   graph->SetPoint(6, 1676,0.07);
   graph->SetPoint(7, 1742,0.08);
   graph->SetPoint(8, 1801,0.09);
   graph->SetPoint(9, 1844,0.1);
   graph->SetPoint(10,1881,0.11);
   
   cout << " Don't forget to change this..." << endl;
   TH1 *Graph6 = new TH1F("Graph6","",100,863,1881);
   Graph6->SetMinimum(0);
   Graph6->SetMaximum(0.12);
   Graph6->SetDirectory(0);
   Graph6->SetStats(0);
   Graph6->GetXaxis()->SetTitle("M_{1} [GeV]");
   Graph6->GetYaxis()->SetTitle("Coupling k/#bar{M}_{Pl}");
   Graph6->GetXaxis()->SetLabelFont(42);
   Graph6->GetYaxis()->SetLabelFont(42);
   Graph6->GetYaxis()->SetTitleOffset(1.8);
   graph->SetHistogram(Graph6);
   
   graph->GetYaxis()->SetTitleOffset(1.19);
   graph->Draw("al");


   TGraph* graph2 = new TGraph(11);
   graph2->SetMarkerColor(ci);
   graph2->SetMarkerStyle(20);
   graph2->SetMarkerSize(0.0);
   graph2->SetPoint(0,  845,0.01);
   graph2->SetPoint(1, 1133,0.02);
   graph2->SetPoint(2, 1275,0.03);
   graph2->SetPoint(3, 1396,0.04);
   graph2->SetPoint(4, 1504,0.05);
   graph2->SetPoint(5, 1598,0.06);
   graph2->SetPoint(6, 1677,0.07);
   graph2->SetPoint(7, 1743,0.08);
   graph2->SetPoint(8, 1801,0.09);
   graph2->SetPoint(9, 1844,0.1);
   graph2->SetPoint(10,1881,0.11);
   graph2->SetLineStyle(kDashed);
   graph2->SetLineColor(ci);
   graph2->SetLineWidth(3);
   graph2->GetXaxis()->SetLabelFont(42);
   graph2->GetYaxis()->SetLabelFont(42);
   graph2->Draw("plsame");


//    graph = new TGraph(3);
//    graph->SetName("Graph");
//    graph->SetTitle("");
//    graph->SetFillColor(1);

//    ci = TColor::GetColor("#0000ff");
//    graph->SetLineColor(ci);
//    graph->SetLineWidth(3);

//    ci = TColor::GetColor("#0000ff");
//    graph->SetMarkerColor(ci);
//.........这里部分代码省略.........

Float_t FitOmega(TH1D *hist, char *plotFilePrefix)
{
	Float_t fPeak, fWidth, SumLo, SumHi;
	Float_t xLo = 0.5;  // lower value of x-axis for drawing histogram
	Float_t xHi = 1.0;   // upper value of x-axis for drawing histogram
	Float_t PeakLo = 0.7; // lower limit on the peak range
	Float_t PeakHi = 0.9; // upper limit on the peak range
	Float_t Nsigma = 3.0;  // number of sigma of the peak
	
	Int_t ibg = 3; // parameter index in the par[] array where the background parameters start

	Int_t i, j, k;
	Float_t x, pval;
	Float_t yield;
	
    cout << "Fitting histogram " << hist->GetName() << endl;
    
	TH1D *hBgFit; // histogram of background
		
	// create canvas
	char title[100];
	char xtitle[100];
	char ytitle[100];
	sprintf(title,"Fitting Analysis"); // canvas title
	TCanvas *can1 = new TCanvas("can1",title,0,0,600,600); // create the canvas
	
//	gStyle->SetOptStat(1111);
	gStyle->SetOptFit(1111);
	can1->SetBorderMode(0);  //Bordermode (-1=down, 0 = no border, 1=up)
	can1->SetBorderSize(5);  
	can1->SetFillStyle(4000);
    can1->Divide(2,2);
	
	sprintf(xtitle,"Invariant Mass (GeV)"); // set the x-axis title
	sprintf(ytitle,"Counts"); // set the y-axis title
		
	// histogram of background
	hBgFit = (TH1D*)hist->Clone("hBgFit"); // clone original hist. into background temp. hist.
	hBgFit->SetName("hBgFit");
	hBgFit->SetTitle("Background");

	// fit the histogram
	Double_t par[7]={1.0,0.782,1.0,1.0,1.0,1.0,1.0};
    TF1 *g1 = new TF1("g1",gaussFit,PeakLo,PeakHi,3); // declare fit fcn
    TF1 *pol = new TF1("pol",polFit,xLo,xHi,4);
    TF1 *t1 = new TF1("t1",totFit,xLo,xHi,7);
		
    can1->cd(1);
	g1->SetParameters(&par[0]);    // set parameters for initial peak fit
	hist->Fit("g1","WR");           // fit the peak
	g1->GetParameters(&par[0]);    // get parameters from initial peak fit
    g1->Draw("same");

    can1->cd(2);
	hist->Fit("pol","WR+");         // fit the background
	pol->GetParameters(&par[ibg]); // get parameters fromt background fit
    pol->Draw("same");

	t1->SetParameters(par);       // set the parameters from initial fits for total fit
	t1->SetLineWidth(5);          // make fit line thicker
	t1->SetLineColor(1);          // set the fit line color
    can1->cd(3);
	hist->Fit("t1","R");          // fit spectrum with total fit function
	t1->GetParameters(&par[0]);   //get final parameters
	
    can1->cd();
	fPeak = t1->GetParameter(1);  // get peak centroid
	fWidth = t1->GetParameter(2); // get peak width
	SumLo = fPeak - Nsigma*fWidth; // calc lower edge at Nsigma away from centroid
	SumHi = fPeak + Nsigma*fWidth; // calc upper edge at Nsigma away from centroid
    cout << "fPeak:  " << fPeak << endl;
			
	pol->SetParameters(&par[ibg]); // set the pfinal parameters for the background function
	hBgFit->Add(pol,-1.0); // subtract background function
	for(k=1; k<hBgFit->GetNbinsX(); k++){
		x = hBgFit->GetBinCenter(k); // read bin center value on the x-axis
		if(x>=SumLo && x<=SumHi){ // check that x is in the peak summation region
				pval = hBgFit->GetBinContent(k); // get the number of counts in the bin
			if(pval<0.0) pval = 0.0; // if neg. counts, set to zero
		}else{
			pval =0.0;
		}
		hBgFit->SetBinContent(k,pval); // refill the histogram
	}
	yield = hBgFit->Integral(hBgFit->FindBin(SumLo),hBgFit->FindBin(SumHi)); // sum total counts in peak

    //g1, pol, t1 TODO
    Double_t binWidth = hist->GetXaxis()->GetBinWidth(0);
    cout << "Bin width:  " << binWidth << endl;
    cout << "SumLo:  " << SumLo << "\nSumHi:  " << SumHi << endl;
    TCanvas *can2 = new TCanvas("can2","Fit Functions",0,0,600,600); // create the canvas
    can2->Divide(2);
    can2->cd(1);
    g1->GetXaxis()->SetRangeUser(SumLo,SumHi);
    g1->Draw();
    can2->cd(2);
    pol->GetXaxis()->SetRangeUser(SumLo,SumHi);
    pol->Draw();
    Double_t peakSum = g1->Integral(SumLo,SumHi);
    Double_t bgSum = pol->Integral(SumLo,SumHi);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    gNonPrJpsiTrad_pt6530y1216->SetMarkerStyle(25);
    gNonPrJpsiTrad_pt6530y1624->SetMarkerStyle(27);
    gNonPrJpsiTrad_pt365y1624->SetMarkerStyle(28);

    // pr
    gPrJpsiTrad->SetMarkerSize(1.2);

    gPrJpsiTrad_pt6530y012->SetMarkerSize(1.2);
    gPrJpsiTrad_pt6530y1216->SetMarkerSize(1.2);
    gPrJpsiTrad_pt6530y1624->SetMarkerSize(2.0);
    gPrJpsiTrad_pt365y1624->SetMarkerSize(1.7);

    // nonPr
    gNonPrJpsiTrad->SetMarkerSize(2.0);

    gNonPrJpsiTrad_pt6530y012->SetMarkerSize(1.2);
    gNonPrJpsiTrad_pt6530y1216->SetMarkerSize(1.2);
    gNonPrJpsiTrad_pt6530y1624->SetMarkerSize(2.0);
    gNonPrJpsiTrad_pt365y1624->SetMarkerSize(1.7);

    // same for the raa denominator and nominator yields histograms
    for(int ih=0; ih<nInHist; ih++) // for each kinematic range
    {
        ahRatio_pr_pp[ih]->SetMarkerStyle(20);
        ahRatio_npr_pp[ih]->SetMarkerStyle(20);
        ahRatio_pr_aa[ih]->SetMarkerStyle(20);
        ahRatio_npr_aa[ih]->SetMarkerStyle(20);
    }


    //-------------------------------------------
    TF1 *f4 = new TF1("f4","1",0,400);
    f4->SetLineWidth(1);
    f4->GetXaxis()->SetTitle("N_{part}");
    f4->GetYaxis()->SetTitle("R_{AA}");
    f4->GetYaxis()->SetRangeUser(0.0,1.5);
    f4->GetXaxis()->CenterTitle(kTRUE);

    TLegend *leg11a = new TLegend(0.6,0.5,0.85,0.65);
    leg11a->SetFillStyle(0);
    leg11a->SetFillColor(0);
    leg11a->SetBorderSize(0);
    leg11a->SetMargin(0.2);
    leg11a->SetTextSize(0.04);
    leg11a->AddEntry(gPrJpsi_pt6530y012,"|y|<1.2","P");
    leg11a->AddEntry(gPrJpsi_pt6530y1216,"1.2<|y|<1.6","P");

    TLegend *leg11b = new TLegend(0.5,0.52,0.8,0.65);
    leg11b->SetFillStyle(0);
    leg11b->SetFillColor(0);
    leg11b->SetBorderSize(0);
    leg11b->SetMargin(0.2);
    leg11b->SetTextSize(0.04);
    leg11b->AddEntry(gPrJpsi_pt6530y1624,"6.5<p_{T}<30 GeV/c","P");
    leg11b->AddEntry(gPrJpsi_pt365y1624,"3<p_{T}<6.5 GeV/c","P");

    TLegend *leg21a = new TLegend(0.6,0.5,0.85,0.65);
    leg21a->SetFillStyle(0);
    leg21a->SetFillColor(0);
    leg21a->SetBorderSize(0);
    leg21a->SetMargin(0.2);
    leg21a->SetTextSize(0.04);
    leg21a->AddEntry(gNonPrJpsi_pt6530y012,"|y|<1.2","P");
    leg21a->AddEntry(gNonPrJpsi_pt6530y1216,"1.2<|y|<1.6","P");

    TLegend *leg21b = new TLegend(0.55,0.5,0.8,0.65);

//.........这里部分代码省略.........
   
   Int_t n=1024;
   TH1D *hsin = new TH1D("hsin", "hsin", n+1, 0, 1023);
   Double_t x;
   //hsin->Fit( model,"MLR" );
   //Fill the histogram with function values
   for (Int_t i=0; i<=n; i++){
     /*
     if( i >= n/2 )
       {
	 x = (Double_t(i-(n/2+1))/n)*(160*TMath::Pi());
       }
     else
       {
	 x = -80*TMath::Pi()+(Double_t(i)/n)*(160*TMath::Pi());
       }
     */
     x = (Double_t(i)/n)*(1024);
     //std::cout << "n: " << i << " x: " << x << std::endl;
     hsin->SetBinContent(i+1, fsin->Eval(x));
   }
   
   hsin->Fit( model2,"MLR" );
   //TFile* fn = new TFile("/Users/cmorgoth/Software/git/TimingAna_New/CIT_Laser_022015_69_ana.root", "READ");
   TFile* fn = new TFile("/Users/cmorgoth/Work/data/LaserDataAtCaltech/02282015/CIT_Laser_022015_69_ana.root", "READ");
   TH1F* pulse = (TH1F*)fn->Get("CH2pulse");
   //hsin->Draw("same");
   hsin->SetLineColor(kGreen-4);
   hsin->Draw();
   model->Draw("same");
   //pulse->SetAxisRange(650, 780, "X");
   pulse->Scale(22.0);
   pulse->Draw("same");
   fsin->GetXaxis()->SetLabelSize(0.05);
   fsin->GetYaxis()->SetLabelSize(0.05);

   c1_2->cd();
   //Compute the transform and look at the magnitude of the output
   TH1 *hm =0;
   TVirtualFFT::SetTransform(0);
   //hm = hsin->FFT(hm, "MAG");
   hm = pulse->FFT(hm, "MAG");
   hm->SetTitle("Magnitude of the 1st transform");
   //hm->Draw();
   double sf = 5e3;//to go from sample to picosecons and also from Hz to MHz
   double range = sf*(double)n/(1023.);
   int n_bin_fft = hm->GetNbinsX();
   TH1F* hmr = new TH1F( "hmr" ,"Magnitude of the 1st transform Rescaled", n_bin_fft, 0, range);
   for( int i = 1; i <= n_bin_fft; i++)
     {
       double bc = hm->GetBinContent( i )/sqrt( n );
       hmr->SetBinContent( i, bc );
     }
   hmr->SetXTitle("f (MHz)");
   hmr->Draw();
   //Transfor to the theoretical function
   TH1 *hm2 =0;
   TVirtualFFT::SetTransform(0);
   hm2 = hsin->FFT(hm2, "MAG");
   hm2->SetLineColor(2);
   //hm2->Draw("same");
   TH1F* hmr2 = new TH1F( "hmr2" ,"Magnitude of the 1st transform Rescaled", n_bin_fft, 0, range);
   for( int i = 1; i <= n_bin_fft; i++)
     {
       double bc = hm2->GetBinContent( i )/sqrt( n );
       hmr2->SetBinContent( i, bc );

void normalDist() {

   TF1 *pdfunc  = new TF1("pdf","ROOT::Math::normal_pdf(x, [0],[1])",-5,5);
   TF1 *cdfunc  = new TF1("cdf","ROOT::Math::normal_cdf(x, [0],[1])",-5,5);
   TF1 *ccdfunc = new TF1("cdf_c","ROOT::Math::normal_cdf_c(x, [0])",-5,5);
   TF1 *qfunc   = new TF1("quantile","ROOT::Math::normal_quantile(x, [0])",0,1);
   TF1 *cqfunc  = new TF1("quantile_c","ROOT::Math::normal_quantile_c(x, [0])",0,1);

   pdfunc->SetParameters(1.0,0.0);  // set sigma to 1 and mean to zero
   pdfunc->SetTitle("");
   pdfunc->SetLineColor(kBlue);

   pdfunc->GetXaxis()->SetLabelSize(0.06);
   pdfunc->GetXaxis()->SetTitle("x");
   pdfunc->GetXaxis()->SetTitleSize(0.07);
   pdfunc->GetXaxis()->SetTitleOffset(0.55);
   pdfunc->GetYaxis()->SetLabelSize(0.06);

   cdfunc->SetParameters(1.0,0.0);  // set sigma to 1 and mean to zero
   cdfunc->SetTitle("");
   cdfunc->SetLineColor(kRed);

   cdfunc->GetXaxis()->SetLabelSize(0.06);
   cdfunc->GetXaxis()->SetTitle("x");
   cdfunc->GetXaxis()->SetTitleSize(0.07);
   cdfunc->GetXaxis()->SetTitleOffset(0.55);

   cdfunc->GetYaxis()->SetLabelSize(0.06);
   cdfunc->GetYaxis()->SetTitle("p");
   cdfunc->GetYaxis()->SetTitleSize(0.07);
   cdfunc->GetYaxis()->SetTitleOffset(0.55);

   ccdfunc->SetParameters(1.0,0.0);  // set sigma to 1 and mean to zero
   ccdfunc->SetTitle("");
   ccdfunc->SetLineColor(kGreen);

   qfunc->SetParameter(0, 1.0);  // set sigma to 1
   qfunc->SetTitle("");
   qfunc->SetLineColor(kRed);
   qfunc->SetNpx(1000); // to get more precision for p close to 0 or 1

   qfunc->GetXaxis()->SetLabelSize(0.06);
   qfunc->GetXaxis()->SetTitle("p");
   qfunc->GetYaxis()->SetLabelSize(0.06);
   qfunc->GetXaxis()->SetTitleSize(0.07);
   qfunc->GetXaxis()->SetTitleOffset(0.55);
   qfunc->GetYaxis()->SetTitle("x");
   qfunc->GetYaxis()->SetTitleSize(0.07);
   qfunc->GetYaxis()->SetTitleOffset(0.55);

   cqfunc->SetParameter(0, 1.0);  // set sigma to 1
   cqfunc->SetTitle("");
   cqfunc->SetLineColor(kGreen);
   cqfunc->SetNpx(1000);

   TCanvas * c1 = new TCanvas("c1","Normal Distributions",100,10,600,800);

   c1->Divide(1,3);
   c1->cd(1);

   pdfunc->Draw();
   TLegend *legend1 = new TLegend(0.583893,0.601973,0.885221,0.854151);
   legend1->AddEntry(pdfunc,"normal_pdf","l");
   legend1->Draw();

   c1->cd(2);
   cdfunc->Draw();
   ccdfunc->Draw("same");
   TLegend *legend2 = new TLegend(0.585605,0.462794,0.886933,0.710837);
   legend2->AddEntry(cdfunc,"normal_cdf","l");
   legend2->AddEntry(ccdfunc,"normal_cdf_c","l");
   legend2->Draw();

   c1->cd(3);
   qfunc->Draw();
   cqfunc->Draw("same");
   TLegend *legend3 = new TLegend(0.315094,0.633668,0.695179,0.881711);
   legend3->AddEntry(qfunc,"normal_quantile","l");
   legend3->AddEntry(cqfunc,"normal_quantile_c","l");
   legend3->Draw();
}

//.........这里部分代码省略.........
  gPrJpsiTrad_mb->SetMarkerStyle(27);
  gNonPrJpsiTrad_mb->SetMarkerStyle(27);
  
  // marker size
  gPrJpsiTrad->SetMarkerSize(1.2);
  gPrJpsiTrad_y_y->SetMarkerSize(1.7);
  gNonPrJpsiTrad->SetMarkerSize(1.7);
  gNonPrJpsiTrad_y_y->SetMarkerSize(1.7);

  //mb
  gPrJpsiTrad_mb->SetMarkerSize(1.5);
  gNonPrJpsiTrad_mb->SetMarkerSize(1.5);

 // same for the raa denominator and nominator yields histograms
  for(int ih=0; ih<nInHist;ih++) // for each kinematic range
  {
    ahRatio_pr_pp[ih]->SetMarkerStyle(20);
    ahRatio_npr_pp[ih]->SetMarkerStyle(20);
    ahRatio_pr_aa[ih]->SetMarkerStyle(20);
    ahRatio_npr_aa[ih]->SetMarkerStyle(20);
    if(ih==2)
      {
        ahRatio_pr_pp[ih]->SetMarkerColor(kCyan+2);
        ahRatio_npr_pp[ih]->SetMarkerColor(kCyan+2);
        ahRatio_pr_aa[ih]->SetMarkerColor(kCyan+2);
        ahRatio_npr_aa[ih]->SetMarkerColor(kCyan+2);
      }
  }


  //-------------------------------------------
  TF1 *f4 = new TF1("f4","1",0,2.4);
  f4->SetLineWidth(1);
  f4->GetXaxis()->SetTitle("|y|");
  f4->GetXaxis()->SetNdivisions(-6);
  f4->GetYaxis()->SetTitle("R_{AA}");
  f4->GetYaxis()->SetRangeUser(0.0,1.5);
  f4->GetXaxis()->CenterTitle(kTRUE);

  //---------------- general stuff
  TLatex *lcent = new TLatex(1.1,1.03,"Cent. 0-100%");
  lcent->SetTextFont(42);
  lcent->SetTextSize(0.05);

  TLatex *lpt = new TLatex(1.1,0.9,"6.5 < p_{T} < 30 GeV/c");
  lpt->SetTextFont(42);
  lpt->SetTextSize(0.05);

  TLatex *lPr = new TLatex(0.2,1.35,"Prompt J/#psi");
  lPr->SetTextFont(42);
  lPr->SetTextSize(0.05);

  TLatex *lNpr = new TLatex(0.2,1.35,"Non-prompt J/#psi");
  lNpr->SetTextFont(42);
  lNpr->SetTextSize(0.05);

 //  -----------------------for comparison purposes
  // axis for the yields
  TF1 *fBin = new TF1("fBin","1",0,8);
  fBin->SetLineWidth(1);
  fBin->GetXaxis()->SetTitle("Bin number");
  fBin->GetYaxis()->SetTitle("Yield ratio");
  fBin->GetYaxis()->SetRangeUser(0.5,2);
  fBin->GetXaxis()->CenterTitle(kTRUE);

  TLatex *lRatio = new TLatex(0.5,1.7,"Yield: Ev-by-Ev/Trad. correction");