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

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 plotsi()
{
  TCanvas *c1 = new TCanvas();
  TF1 *d = new TF1("d",si_dose,0.001,10.0,0);
  d->SetTitle("Si-1 LBSD diode;#DeltaV [V];#Phi [10^{12} cm^{-2}]");
  d->Draw();
  c1->SetLogx();
  c1->SetLogy();
  c1->SetGridx();
  c1->SetGridy();
  c1->Print("si_dose.pdf","pdf portrait");
}

void plotradfet()
{
  TCanvas *c1 = new TCanvas();
  TF1 *d = new TF1("d",radfet_dose,0.001,100.0,0);
  d->SetTitle("REM-501C Type K 250 nm RadFET;#DeltaV [V];Dose [Gy]");
  d->Draw();
  c1->SetLogx();
  c1->SetLogy();
  c1->SetGridx();
  c1->SetGridy();
  c1->Print("radfet_dose.pdf","pdf portrait");
}

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 FitWithConstant(TH1D* h1, TFile *out)
{
  Double_t fitRangeLow = 0.;
  Double_t fitRangeHigh = .4;
  TF1 *f = new TF1("f","[0]", fitRangeLow, fitRangeHigh);

  Int_t binLow = h1->FindBin(fitRangeLow);
  Int_t binHigh = h1->FindBin(fitRangeHigh);
  cout<<"\t\t********************Fitting with a constant value********************"<<endl<<endl;;
  h1->Fit(f, "R0");

  TString fitName = h1->GetName();
  fitName += "Fit";
  f->SetName(fitName);
  f->SetTitle(fitName);

  // TFile out("Compare.root","update");
  TDirectory *dir = out->GetDirectory("Fit");
  if(!dir) dir = out->mkdir("Fit");
  dir->cd();
  // f->SetDirectory(0); 
  f->Write(f->GetName(), TObject::kOverwrite);
  // out->Close();
}

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

   //gSystem->Load("libMathMore");
   // this is the way to force load of MathMore in Cling
   ROOT::Math::MathMoreLibrary::Load();

   int n=100;
   double a=-5.;
   double b=5.;
   //double r  = 3;
   TF1* pdf = new TF1("pdf", "ROOT::Math::tdistribution_pdf(x,3.0)", a,b);
   TF1* cum = new TF1("cum", "ROOT::Math::tdistribution_cdf(x,3.0)", a,b);

   TH1D* quant = new TH1D("quant", "", 9, 0, 0.9);

   for(int i=1; i < 10; i++)
      quant->Fill((i-0.5)/10.0, ROOT::Math::tdistribution_quantile((1.0*i)/10, 3.0 ) );

   double xx[10];
   xx[0] = -1.5;
   for(int i=1; i<9; i++)
      xx[i]= quant->GetBinContent(i);
   xx[9] = 1.5;
   TH1D* pdfq[10];
   //int nbin = n/10.0;
   for(int i=0; i < 9; i++) {
      int nbin = n * (xx[i+1]-xx[i])/3.0 + 1.0;
      TString name = "pdf";
      name += i;
      pdfq[i]= new TH1D(name, "", nbin,xx[i],xx[i+1] );
      for(int j=1; j<nbin; j++) {
         double x= j*(xx[i+1]-xx[i])/nbin + xx[i];
         pdfq[i]->SetBinContent(j, ROOT::Math::tdistribution_pdf(x,3));
      }
   }

   TCanvas *Canvas = new TCanvas("DistCanvas", "Student Distribution graphs", 10, 10, 800, 700);
   pdf->SetTitle("Student t distribution function");
   cum->SetTitle("Cumulative for Student t");
   quant->SetTitle("10-quantiles  for Student t");
   Canvas->Divide(2, 2);
   Canvas->cd(1);
   pdf->SetLineWidth(2);
   pdf->DrawCopy();
   Canvas->cd(2);
   cum->SetLineWidth(2);
   cum->SetLineColor(kRed);
   cum->Draw();
   Canvas->cd(3);
   quant->Draw();
   quant->SetLineWidth(2);
   quant->SetLineColor(kBlue);
   quant->SetStats(0);
   Canvas->cd(4);
   pdfq[0]->SetTitle("Student t & its quantiles");
   pdf->SetTitle("");
   pdf->Draw();
   //pdfq[0]->SetAxisRange(-1.5, 0, 1.5,1.0);
   pdfq[0]->SetTitle("Student t & its quantiles");
   for(int i=0; i < 9; i++) {
      pdfq[i]->SetStats(0);
      pdfq[i]->SetFillColor(i+1);
      pdfq[i]->Draw("same");
   }
   Canvas->Modified();
   Canvas->cd();
}

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 makeGenLevel_pTRatios(int erg_tev){
  char TREE_NAME[]="candTree";
  TString OUTPUT_NAME = Form("compareProductionModes_pT_%iTeV.root", erg_tev);
  char csuffix[] = "_pToverMZZ";
  char cxtitle[] = "p_{T} / m_{4l}";

  TFile* foutput = new TFile(OUTPUT_NAME, "recreate");
  TH1F* h_pt[nProdModes];
  TH1F* hratio_pt[nProdModes-1];
  TProfile* hpr_pt[nProdModes-1];
  TGraphAsymmErrors* tgratio_pt[nProdModes-1];
  TF1* fratio_pt[nProdModes];
  double xmin = 0;
  double xmax = 6;
  const int nbins = 23;
  double xbins[nbins+1];
  for (int bin=0; bin<=nbins; bin++){
    if (bin==0) xbins[bin] = xmin;
    else if (xbins[bin-1]<1.2) xbins[bin] = xbins[bin-1] + 0.1; // 12
    else if (xbins[bin-1]<2.4) xbins[bin] = xbins[bin-1] + 0.2; // 6
    else if (xbins[bin-1]<4) xbins[bin] = xbins[bin-1] + 0.4; // 4
    else if (xbins[bin-1]<xmax) xbins[bin] = xbins[bin-1] + 2.0; // 1
  }
  TString cytitle = "Rate / bin";
  for (int p = 0; p < nProdModes; p++){
    TString hname = Form("%s%s", prodName[p], csuffix);
//    TString htitle = Form("%s %i TeV", prodName[p], erg_tev);
    h_pt[p] = new TH1F(hname, "", nbins, xbins);
    h_pt[p]->Sumw2();
    h_pt[p]->GetXaxis()->SetTitle(cxtitle);
    h_pt[p]->GetYaxis()->SetTitle(cytitle);

    if (p>0){
      hname = Form("%s%s_profile", prodName[p], csuffix);
      hpr_pt[p-1] = new TProfile(hname, "", nbins, xbins);
      hpr_pt[p-1]->Sumw2();
    }
  }
  int genFinalState = -1;
  float GenHMass = 0;
  float GenHPt = 0;
  TString cinput_common = location_primaryTrees;
  if (erg_tev==7) cinput_common.Append("141217/");
  else if (erg_tev==8) cinput_common.Append("141021/");
  for (int p = 0; p < nProdModes; p++){
    TString cinput = cinput_common;
//    cinput.Append(Form("%s/ZZ4lAnalysis_%s%s", sampleName[p], sampleName[p], ".root"));
    cinput.Append(Form("ZZ4lAnalysis_%s%s", sampleName[p], ".root"));
    TFile* finput = new TFile(cinput, "read");
    cout << "Production " << prodName[p] << endl;
    TTree* tin[3] ={
      (TTree*)finput->Get(Form("ZZ4muTree/%s", TREE_NAME)),
      (TTree*)finput->Get(Form("ZZ4eTree/%s", TREE_NAME)),
      (TTree*)finput->Get(Form("ZZ2e2muTree/%s", TREE_NAME))
    };
    for (int f = 0; f < 3; f++){
      cout << "Channel " << f << endl;
      cout << "Nentries: " << tin[f]->GetEntries() << endl;
      tin[f]->SetBranchAddress("genFinalState", &genFinalState);
      tin[f]->SetBranchAddress("GenHMass", &GenHMass);
      tin[f]->SetBranchAddress("GenHPt", &GenHPt);
      for (int ev = 0; ev < tin[f]->GetEntries(); ev++){
        //			for (int ev = 0; ev < 1000; ev++){
        genFinalState = -1;
        GenHMass = 0;
        GenHPt = 0;
        tin[f]->GetEntry(ev);
        if (f != genFinalState) continue;
        double kd = GenHPt / GenHMass;
        h_pt[p]->Fill(kd);
        if (p==0){
          for (int pp=0; pp<nProdModes-1; pp++) hpr_pt[pp]->Fill(kd, kd);
        }
        else hpr_pt[p-1]->Fill(kd, kd);
      }
    }
    finput->Close();
    cout << "File closed" << endl;
  }
  foutput->cd();
/*  double zhxsec, whxsec;
  if (erg_tev == 7){ whxsec = 0.5688; zhxsec = 0.3299; }
  else if (erg_tev == 8){ whxsec = 0.6931; zhxsec = 0.4091; }
  h_pt[nProdModes]->Add(h_pt[nProdModes-1], whxsec);
  h_pt[nProdModes]->Add(h_pt[nProdModes-2], zhxsec);
*/  

  double max_plot = 0;
  for (int p = 0; p < nProdModes; p++){
    h_pt[p]->Scale(1. / h_pt[p]->Integral(0, h_pt[p]->GetNbinsX() + 1));

    h_pt[p]->GetYaxis()->SetRangeUser(0, 1);
    h_pt[p]->SetLineWidth(2);
    h_pt[p]->SetLineStyle(1);
    if (p==0) h_pt[p]->SetLineColor(kBlack);
    else if (p==3) h_pt[p]->SetLineColor(kRed);
    else if (p==4) h_pt[p]->SetLineColor(kGreen+2);
    else if (p==2) h_pt[p]->SetLineColor(kBlue);
    else if (p==1) h_pt[p]->SetLineColor(kViolet);
    if (p==0) h_pt[p]->SetMarkerColor(kBlack);
//.........这里部分代码省略.........

int compareL2L3Correction(const std::string& kalibri, const std::string& jetMETL2, const std::string& jetMETL3) {
//int compareL2L3Correction(const std::string& jetMETL2, const std::string& jetMETL3) {
    // -- Read calibration constants ---------------
    std::cout << "Reading calibration constants...\n";
    parmap corrKalibri  = readParameters(kalibri);
    parmap corrJetMETL2 = readParameters(jetMETL2);
    parmap corrJetMETL3 = readParameters(jetMETL3);
    if( (1 + corrKalibri.size() ) != ( corrJetMETL2.size() + corrJetMETL3.size() ) ) {
        std::cerr << "Linker Kompaktierer defekt. Kontaktieren Sie das Personal (68/111)." << std::endl;
        return 1;
    }


    // -- Global variables -------------------------
    std::vector<TF1*> fCorrKalibri(1+corrKalibri.size());
    std::vector<TF1*> fCorrJetMET(1+corrKalibri.size());
    TF1 * f = 0;
    char name[50];
    std::string outFileName = "comparisonL2L3_Summer08_IC5Calo.ps";
    bool printPar = false;


    // -- Create correction plots for JetMET -------
    std::cout << "Creating JetMET plots...\n";
    // L3 correction
    if( printPar ) std::cout << "  JetMETL3: " << std::flush;
    std::vector<double>& par = corrJetMETL3[0];
    f = new TF1("fJetMETL3",corrL3,par.at(2),par.at(3),4);
    for(int i = 0; i < 4; i++) {
        f->SetParameter(i,par.at(4+i));
        if( printPar ) std::cout << par.at(4+i) << " ";
    }
    if( printPar ) std::cout << "\n";
    fCorrJetMET.at(0) = f;

    // L2 correction
    for(unsigned int bin = 0; bin < corrJetMETL2.size(); bin++) {
        if( printPar ) std::cout << "  JetMETL2 " << bin << ": " << std::flush;
        par = corrJetMETL2[bin];
        sprintf(name,"fJetMETL2_%i",bin);
        f = new TF1(name,corrL2,par.at(2),par.at(3),3);
        sprintf(name,"Bin %i: %.3f <  #eta < %.3f",bin,par.at(0),par.at(1));
        f->SetTitle(name);
        for(int i = 0; i < 3; i++) {
            f->SetParameter(i,par.at(4+i));
            if( printPar ) std::cout << par.at(4+i) << " ";
        }
        if( printPar ) std::cout << "\n";
        fCorrJetMET.at(1+bin) = f;
    }

    // Nice plots
    for(unsigned i = 0; i < fCorrJetMET.size(); i++) {
        fCorrJetMET.at(i)->SetLineWidth(2);
        fCorrJetMET.at(i)->SetLineColor(4);
    }


    // -- Create correction plots for Kalibri ------
    std::cout << "Creating Kalibri plots...\n";
    int nTowerPar     = 0;
    int nJetPar       = 3;
    int nTrackPar     = 0;
    int nGlobalJetPar = 4;
    // Global correction
    if( printPar ) std::cout << "  Kalibri: " << std::flush;
    par = corrKalibri[0];
    f = new TF1("fKalibri_global",corrL3,par.at(2),par.at(3),4);
    for(int i = 0; i < nGlobalJetPar; i++) {
        f->SetParameter(i,par.at(4+nTowerPar+nJetPar+nTrackPar+i));
        if( printPar ) std::cout << par.at(4+nTowerPar+nJetPar+nTrackPar+i) << " ";
    }
//   f = new TF1("fKalibri_global",corrL3,4,2000,4);
//   f->SetParameter(0,0.998293);
//   f->SetParameter(1,5.43056);
//   f->SetParameter(2,3.3444);
//   f->SetParameter(3,2.39809);
    if( printPar ) std::cout << "\n";
    fCorrKalibri.at(0) = f;

    // Local correction
    double scale[3] = { 1., 0.1, 0.01 };
    for(unsigned int bin = 0; bin < corrKalibri.size(); bin++) {
        if( printPar ) std::cout << "  Kalibri " << bin << ": " << std::flush;
        par = corrKalibri[bin];
        sprintf(name,"fKalibri_%i",bin);
        f = new TF1(name,corrL2,par.at(2),par.at(3),3);
        sprintf(name,"Bin %i: %.3f <  #eta < %.3f",bin,par.at(0),par.at(1));
        f->SetTitle(name);
        for(int i = 0; i < nJetPar; i++) {
            f->SetParameter(i,scale[i]*par.at(4+nTowerPar+i));
            if( printPar ) std::cout << scale[i]*par.at(4+nTowerPar+i) << " ";
        }
        if( printPar ) std::cout << "\nDone";
        fCorrKalibri.at(1+bin) = f;
    }

    // Nice plots
    for(unsigned i = 0; i < fCorrKalibri.size(); i++) {
        fCorrKalibri.at(i)->SetLineWidth(2);
//.........这里部分代码省略.........

void checkRsnPIDqa(TString filename, TString foldername, Bool_t savePng,
		   TString plotTPCpi, TString plotTPCka, TString plotTPCpro,
		   TString plotTTOFpi, TString plotTOFka, TString plotTOFpro)
{
  //Open input file 
  TFile * fin = TFile::Open(filename.Data());
  if (!fin) return 0x0;
  
  //Access output of specific wagon
  TList * list = (TList*) fin->Get(foldername.Data());
  if (!list) return 0x0;

  //Set range for fit
  Float_t RangeFitMomMin = 0.1; //range in momentum where to check the mean and pull
  Float_t RangeFitMomMax = 2.0;
  Int_t xbinFitMin = 0;
  Int_t xbinFitMax = -1;
  Float_t RangeFitNsigmaPIDmin = -2.0; //range in Nsigma where the fit is to be performed
  Float_t RangeFitNsigmaPIDmax = 2.0;

  //Set range for visualisation
  Float_t RangeShowTPC[2] = {0.1, 2.0}; 
  Float_t RangeShowTOF[2] = {0.25, 2.0};
    
  //--------------------------
  // TPC PID Nsigma
  // fit with simple gaussian
  //--------------------------
  //Gaussian function
  TF1 *fGaus = new TF1("f","gaus", -7.0, 7.0);

  //--- pions
  TH2F * hTPCsigmaPi = (TH2F*)list->FindObject(plotTPCpi.Data());
  hTPCsigmaPi->RebinX(2);
  hTPCsigmaPi->SetTitle("TPC Pions");
  MakeUpHisto(hTPCsigmaPi,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  hTPCsigmaPi->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaPi->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  xbinFitMin = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMin);
  xbinFitMax = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMax);
  hTPCsigmaPi->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  TH1D * hTPCsigmaPi_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_mean");
  TH1D * hTPCsigmaPi_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_pull");
  MakeUpHisto(hTPCsigmaPi_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaPi_pull, "", "", 1, kRed+2, 2);

  //--- kaons
  TH2F * hTPCsigmaKa = (TH2F*)list->FindObject(plotTPCka.Data());
  hTPCsigmaKa->RebinX(2);
  hTPCsigmaKa->SetTitle("TPC Kaons");
  hTPCsigmaKa->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaKa->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  hTPCsigmaKa->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  MakeUpHisto(hTPCsigmaKa,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  TH1D * hTPCsigmaKa_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCka.Data())))->Clone("hNsigmaTPCka_mean");
  TH1D * hTPCsigmaKa_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCka.Data())))->Clone("hNsigmaTPCka_pull");
  MakeUpHisto(hTPCsigmaKa_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaKa_pull, "", "", 1, kRed+2, 2);

  //--- protons
  TH2F * hTPCsigmaPro = (TH2F*)list->FindObject(plotTPCpro.Data());
  hTPCsigmaPro->RebinX(2);
  hTPCsigmaPro->SetTitle("TPC Protons");
  MakeUpHisto(hTPCsigmaPro,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  hTPCsigmaPro->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaPro->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  hTPCsigmaPro->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  TH1D * hTPCsigmaPro_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_mean");
  TH1D * hTPCsigmaPro_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_pull");
  MakeUpHisto(hTPCsigmaPro_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaPro_pull, "", "", 1, kRed+2, 2);

   //--- plot TPC
  TLine *l11=new TLine(RangeShowTPC[0],0.,RangeShowTPC[1],0.); l11->SetLineWidth(1); l11->SetLineStyle(7);
  TLine *l12=new TLine(RangeShowTPC[0],1.,RangeShowTPC[1],1.); l12->SetLineWidth(1); l12->SetLineStyle(7);

  gStyle->SetOptStat(0);
  TCanvas *cPidPerformance4 = new TCanvas("cPIDperformance4","TPC PID",1200,500);
  cPidPerformance4->Divide(3,1);
  cPidPerformance4->cd(1);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaPi->DrawCopy("colz");
  hTPCsigmaPi_mean->DrawCopy("same");
  hTPCsigmaPi_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

  cPidPerformance4->cd(2);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaKa->DrawCopy("colz");
  hTPCsigmaKa_mean->DrawCopy("same");
  hTPCsigmaKa_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

  cPidPerformance4->cd(3);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaPro->DrawCopy("colz");
  hTPCsigmaPro_mean->DrawCopy("same");
  hTPCsigmaPro_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

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

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();
}