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

void peaks(Int_t np=10) {
   npeaks = TMath::Abs(np);
   TH1F *h = new TH1F("h","test",500,0,1000);
   //generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
   Int_t p;
   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 1;
      par[3*p+3] = 10+gRandom->Rndm()*980;
      par[3*p+4] = 3+2*gRandom->Rndm();
   }
   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
   f->SetNpx(1000);
   f->SetParameters(par);
   TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
   c1->Divide(1,2);
   c1->cd(1);
   h->FillRandom("f",200000);
   h->Draw();
   TH1F *h2 = (TH1F*)h->Clone("h2");
   //Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,2,"",0.10);
   printf("Found %d candidate peaks to fit\n",nfound);
   //Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,20,"same");
   if (hb) c1->Update();
   if (np <0) return;

   //estimate linear background using a fitting method
   c1->cd(2);
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h->Fit("fline","qn");
   //Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Double_t *xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
      par[3*npeaks+4] = 3;
      npeaks++;
   }
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   //we may have more than the default 25 parameters
   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h2->Fit("fit");
}

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 draw()
{
  TF1* f = new TF1("f",_gaussCrystalball,0,5,8);
  f->SetParameters(39,1.42,0.22,172,2.506,0.14,0.64,1.7);
  f->SetNpx(500);
  f->Draw();
}

void example(double E0 = 50, int nevents = 100000)
{
   TStopwatch timer;

   // 12C* -> 3(4He)
   // compound nucleus = carbon
   KVNucleus CN(6, 12);
   CN.SetExcitEnergy(E0);

   // decay products
   KVEvent decay;
   KVNucleus* n = decay.AddParticle();
   n->SetZandA(2, 4);
   n = decay.AddParticle();
   n->SetZandA(2, 4);
   n = decay.AddParticle();
   n->SetZandA(2, 4);

   MicroStat::mdweight gps;
   Double_t etot = E0 + decay.GetChannelQValue();

   if (etot <= 0) {
      printf("Break-up channel is not allowed\n");
      return;
   }
   gps.SetWeight(&decay, etot);
   gps.initGenerateEvent(&decay);

   TH1F* h1 = new TH1F("h1", "Kinetic energy of alpha particle 3", 200, 0, etot * 2. / 3.);
   h1->Sumw2();

   KVEvent event;

   while (nevents--) {
      gps.GenerateEvent(&decay, &event);
      h1->Fill(event.GetParticle(3)->GetKE());
      gps.resetGenerateEvent();
   }

   h1->Draw();
   TF1* EDis = new TF1("EDis", edist, 0., etot, 3);
   EDis->SetNpx(500);
   EDis->SetParLimits(0, 0, 1.e+08);
   EDis->SetParLimits(1, 0, 2 * etot);
   EDis->FixParameter(2, 3);
   gStyle->SetOptFit(1);
   h1->Fit(EDis, "EM");

   timer.Print();
}

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 FindPeak(TH1 *hm, int * i, char * namefile){
  int np =5, p, max = 0;
  int npeaks = TMath::Abs(np);
  double par[3000];
  par[0] = 0.8;
  par[1] = -0.6/1000;

  for (p=0;p<npeaks;p++) {
    par[3*p+2] = 1;
    par[3*p+3] = 10+gRandom->Rndm()*980;
    par[3*p+4] = 3+2*gRandom->Rndm();
  }
 
  TSpectrum *s = new TSpectrum(2*npeaks,1);
  int nfound = s->Search(hm,2,"",0.10);
  printf("Found %d candidate peaks to fit\n",nfound);

  TH1 *hb = s->Background(hm,20,"same");
  if (np <0) return;

  // loope over peaks
  TF1 *fline = new TF1("fline","pol1",0,1000);
  hm->Fit("fline","qn");
  par[0] = fline->GetParameter(0);
  par[1] = fline->GetParameter(1);
  npeaks = 0;
  float *xpeaks = s->GetPositionX();
  for (p=0;p<nfound;p++) {
    float xp = xpeaks[p];
    int bin = hm->GetXaxis()->FindBin(xp);
    float yp = hm->GetBinContent(bin);
    if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
    par[3*npeaks+2] = yp;
    par[3*npeaks+3] = xp;
    par[3*npeaks+4] = 3;
    npeaks++;
  }
  printf("Found %d useful peaks to fit\n",npeaks); 
  printf("Now fitting: Be patient\n");
  if (max < npeaks) max = npeaks;
  TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
  TVirtualFitter::Fitter(hm,10+3*npeaks);
  fit->SetParameters(par);
  fit->SetNpx(1000);
  hm->Fit("fit"); 
}

void drawShapes (TGraph ** tg_sample_par, int Npar)
  {
    TGraph * log_tg_sample_par0 = makeLog (tg_sample_par[0]) ;
    //PG loop on candidate masses
    for (double mass = 300 ; mass < 1100 ; mass += 25)
      {
        TF1 * func = new TF1 ("func",crystalBallLowHigh, 200, 2000, 7) ;
        func->SetParameter (0, TMath::Exp (log_tg_sample_par0->Eval (mass))) ;
        for (int iPar = 1 ; iPar < Npar ; ++iPar)
          func->SetParameter (iPar, tg_sample_par[iPar]->Eval (mass)) ;
        func->SetLineWidth (1) ;
        func->SetLineColor (kBlue + 1) ;
        func->SetNpx (10000) ;
        func->Draw ("same") ;
      } //PG loop on candidate masses
    return ;
  }

void drawOriginals (TGraph ** tg_sample_par, int color, int Npar)
  {
    //PG loop on candidate masses
    double masses[5] = {350, 500, 650, 800, 1000} ;
    for (int i = 0 ; i < 5 ; ++i)
      {
        mass = masses[i] ;
        TF1 * func = new TF1 ("func",crystalBallLowHigh, 200, 2000, 7) ;
        for (int iPar = 0 ; iPar < Npar ; ++iPar)
          func->SetParameter (iPar, tg_sample_par[iPar]->Eval (mass)) ;
        func->SetLineWidth (1) ;
        func->SetLineColor (color) ;
        func->SetNpx (10000) ;
        func->Draw ("same") ;
      } //PG loop on candidate masses
    return ;
  }

void GetRandomTest(){
	double k0=1.39;
	double k1 = 0.425;
	double theta0 = 3.41;
	double theta1 = 1.30;
	int iNpart=2;
	double k_=k0+k1*(iNpart-2);
        double theta_=theta0+theta1*TMath::Log(iNpart-1);
	TF1 *f = new TF1("f","TMath::GammaDist(x,[0],0,[1])",0,200);
	f->SetParameters(k1,theta_);
	cout<<"Value at 0 = "<<f->Eval(0)<<endl;
	cout<<"Value at 1e-11 = "<<f->Eval(1e-11)<<endl;
	cout<<"Integral 1= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	f->SetRange(1e-12,200);
	cout<<"Integral 2= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	cout<<"fXmin = "<<f->GetXmin()<<"\tfXmax = "<<f->GetXmax()<<"\tfNpx = "<<f->GetNpx()<<endl;
	f->SetNpx(1e5);
	TCanvas *c1 = new TCanvas();
	c1->SetLogy();
	cout<<"f mean = "<<f->Mean(0,200)<<endl;
	cout<<"math mean = "<<f->GetParameter(0)*f->GetParameter(1)<<endl;
	TH1D* h = new TH1D("h","h",1000,0,200);
	for(int i=0;i<1e6;i++){
		double para = f->GetRandom();
		h->Fill(para);
	}
	h->Scale(1.0/h->Integral()*1000/200);
	h->GetYaxis()->SetRangeUser(1e-10,1);
	h->SetMarkerStyle(24);
	h->SetMarkerColor(4);
	h->SetMarkerSize(1.1);
	TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
        leg->SetFillColor(0);
        leg->SetFillStyle(0);
        leg->SetBorderSize(0);
        leg->SetTextFont(42);
        leg->SetTextSize(0.03);
        leg->AddEntry(f,"function","lp");
        leg->AddEntry(h,"filled histogram","lp");
	h->Draw("P");
	f->Draw("same");
        leg->Draw("same");
	cout<<"h mean = "<<h->GetMean(1)<<endl;
	c1->Print("TestGetRandom.png");
	}

Double_t fitfulllang( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){
    cout << hs << " does not exist\n";
    return 0;
  }

  double aa = h->GetEntries();//normalization

  // find peak:
  int ipk = h->GetMaximumBin();
  double xpk = h->GetBinCenter(ipk);
  double sm = xpk / 9; // sigma
  double ns = sm; // noise

  // fit range:
  int ib0 = ipk/2;
  int ib9 = h->GetNbinsX() - 1;

  double x0 = h->GetBinLowEdge(ib0);
  double x9 = h->GetBinLowEdge(ib9) + h->GetBinWidth(ib9);

  // create a TF1 with the range from x0 to x9 and 4 parameters
  TF1 *fitFcn = new TF1( "fitFcn", fitLandauGauss, x0, x9, 4 );

  fitFcn->SetParName( 0, "peak" );
  fitFcn->SetParName( 1, "sigma" );
  fitFcn->SetParName( 2, "area" );
  fitFcn->SetParName( 3, "smear" );

  fitFcn->SetNpx(500);
  fitFcn->SetLineWidth(4);
  fitFcn->SetLineColor(kMagenta);

  // set start values:
  fitFcn->SetParameter( 0, xpk ); // peak position, defined above
  fitFcn->SetParameter( 1, sm ); // width
  fitFcn->SetParameter( 2, aa ); // area
  fitFcn->SetParameter( 3, ns ); // noise

  h->Fit("fitFcn", "NQR", "ep" );// R = range from fitFcn
  return fitFcn->GetParameter(0);
}

void v()
{

    TF1 *signal = new TF1("DGauss", DGauss, 0, 800, 7);
    signal->SetParameter(0, 50);
    signal->SetParameter(1, 1);
    signal->SetParameter(2, 2);
    signal->SetParameter(3, 100);
    signal->SetParameter(4, -20);

    signal->SetParameter(5, 22.374);
    signal->SetParameter(6, 5.55);
    signal->SetNpx(10000);


    signal->Draw();

}

// This is the main program
void slits() {
  float r,ns;

  // request user input
  cout << "slit width / g ? ";
  scanf("%f",&r);
  cout << "# of slits? ";
  scanf("%f",&ns);
  cout <<"interference pattern for "<< ns
       <<" slits, width/distance: "<<r<<endl;

  // define function and set options
  TF1 *Fnslit  = new TF1("Fnslit",nslit,-5.001,5.,2);
  Fnslit->SetNpx(500);
    
  // set parameters, as read in above
  Fnslit->SetParameter(0,r);
  Fnslit->SetParameter(1,ns);
    
  // draw the interference pattern for a grid with n slits
  Fnslit->Draw();
}

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

//.........这里部分代码省略.........
      hDeltaPhiBack12->Fill((Phi1-Phi2)*TMath::RadToDeg());
      hThetaRelBack12->Fill(thetarel12);
    }

    if(TMath::Abs(DeltaT13-corr13) < 500){
      hDeltaTheta13->Fill((Theta1-Theta3)*TMath::RadToDeg());
      hDeltaPhi13->Fill((Phi1-Phi3)*TMath::RadToDeg());
      hThetaRel13->Fill(thetarel13);
    }
    else if(TMath::Abs(DeltaT13-corr13) > 1000 && TMath::Abs(DeltaT13-corr13) < 6000){
      hDeltaThetaBack13->Fill((Theta1-Theta3)*TMath::RadToDeg());
      hDeltaPhiBack13->Fill((Phi1-Phi3)*TMath::RadToDeg());
      hThetaRelBack13->Fill(thetarel13);
    }
  }
  
  h->SetStats(0);

  hDeltaThetaBack12->Sumw2();
  hDeltaPhiBack12->Sumw2();
  hThetaRelBack12->Sumw2();
  hDeltaThetaBack12->Scale(0.1);
  hDeltaPhiBack12->Scale(0.1);
  hThetaRelBack12->Scale(0.1);
  hDeltaThetaBack13->Sumw2();
  hDeltaPhiBack13->Sumw2();
  hThetaRelBack13->Sumw2();
  hDeltaThetaBack13->Scale(0.1);
  hDeltaPhiBack13->Scale(0.1);
  hThetaRelBack13->Scale(0.1);

  Float_t val,eval;
  TCanvas *c1=new TCanvas();
  TF1 *ff = new TF1("ff","[0]*[4]/[2]/sqrt(2*TMath::Pi())*TMath::Exp(-(x-[1])*(x-[1])*0.5/[2]/[2]) + [3]*[4]/6/[2]");
  ff->SetParName(0,"signal");
  ff->SetParName(1,"mean");
  ff->SetParName(2,"sigma");
  ff->SetParName(3,"background");
  ff->SetParName(4,"bin width");
  ff->SetParameter(0,42369);
  ff->SetParameter(1,0);
  ff->SetParLimits(2,10,1000);
  ff->SetParameter(2,150); // fix witdh if needed
  ff->SetParameter(3,319);
  ff->FixParameter(4,20000./nbint); // bin width
  
  ff->SetNpx(1000);
  
  h->Fit(ff);
  
  val = ff->GetParameter(2);
  eval = ff->GetParError(2);
  
  printf("significance = %f\n",ff->GetParameter(0)/sqrt(ff->GetParameter(0) + ff->GetParameter(3)));

  h->Draw();
  
  TF1 *func1 = (TF1 *)  h->GetListOfFunctions()->At(0);
  
  func1->SetLineColor(2);
  h->SetLineColor(4);
  
  TPaveText *text = new TPaveText(1500,(h->GetMinimum()+(h->GetMaximum()-h->GetMinimum())*0.6),9500,h->GetMaximum());
  text->SetFillColor(0);
  sprintf(title,"width = %5.1f #pm %5.1f",func1->GetParameter(2),func1->GetParError(2));
  text->AddText(title);
  sprintf(title,"signal (S) = %5.1f #pm %5.1f",func1->GetParameter(0),func1->GetParError(0));
  text->AddText(title);
  sprintf(title,"background (B) (3#sigma) = %5.1f #pm %5.1f",func1->GetParameter(3),func1->GetParError(3));
  text->AddText(title);
  sprintf(title,"significance (S/#sqrt{S+B}) = %5.1f",func1->GetParameter(0)/sqrt(func1->GetParameter(0)+func1->GetParameter(3)));
  text->AddText(title);
  
  text->SetFillStyle(0);
  text->SetBorderSize(0);
  
  text->Draw("SAME");
  
  printf("n_day = %f\nn_dayGR = %f\n",nsec*1./86400,nsecGR*1./86400);

  text->AddText(Form("rate = %f #pm %f per day",func1->GetParameter(0)*86400/nsecGR,func1->GetParError(0)*86400/nsecGR));

  TFile *fo = new TFile("output-SAVO-010203.root","RECREATE");
  h->Write();
  h2->Write();
  hDeltaTheta12->Write();
  hDeltaPhi12->Write();
  hThetaRel12->Write();
  hDeltaThetaBack12->Write();
  hDeltaPhiBack12->Write();
  hThetaRelBack12->Write();
  hDeltaTheta13->Write();
  hDeltaPhi13->Write();
  hThetaRel13->Write();
  hDeltaThetaBack13->Write();
  hDeltaPhiBack13->Write();
  hThetaRelBack13->Write();
  fo->Close();
  
}

void proc_pileup(const AliVdMMetaData& vdmMetaData,
                 AliVdMScanData& allData,
                 const char* classAC,
                 const char* classAnotC,
                 const char* classCnotA,
                 const std::vector<Double_t>& par0)
{
  typedef std::map<Short_t, TGraphErrors> map_t; // BCID -> TGraphErrors
  map_t gAnotC, gCnotA;  // one-arm/two-arm ratios

  // (1) fill one-arm/two-arm ratio graphs for all BCIDs
  for (Int_t iScan=0; iScan<4; ++iScan) {
    AliVdMTree& vtAND   = allData.GetMap(iScan)[classAC];
    AliVdMTree& vtAnotC = allData.GetMap(iScan)[classAnotC];
    AliVdMTree& vtCnotA = allData.GetMap(iScan)[classCnotA];
    vtAND.Zip3([&gAnotC,&gCnotA](const AliVdMTree::DefaultBranchData& d,
                                 AliVdMTree::branchMap_t& mapAC,
                                 AliVdMTree::branchMap_t& mapAnotC,
                                 AliVdMTree::branchMap_t& mapCnotA)
               {
                 AliVdMTree::ValErr v1 = mapAnotC["rate"];
                 v1 /= mapAC["rate"];
                 // v1 /= mapAC["relBkgd"];
                 v1 *= mapAnotC["relBkgd"];
                 if (!v1.isInf() && v1.val()) {
                   const Int_t m1 = gAnotC[d.BCID()].GetN();
                   gAnotC[d.BCID()].SetPoint     (m1, mapAC["mu"].val(), v1.val());
                   gAnotC[d.BCID()].SetPointError(m1, mapAC["mu"].err(), v1.err());
                 }
                 AliVdMTree::ValErr v2 = mapCnotA["rate"];
                 v2 /= mapAC["rate"];
                 // v2 /= mapAC["relBkgd"];
                 v2 *= mapCnotA["relBkgd"];
                 if (!v2.isInf() && v2.val()) {
                   const Int_t m2 = gCnotA[d.BCID()].GetN();
                   gCnotA[d.BCID()].SetPoint     (m2, mapAC["mu"].val(), v2.val());
                   gCnotA[d.BCID()].SetPointError(m2, mapAC["mu"].err(), v2.err());
                 }
               },
               vtAnotC, vtCnotA);
  }

  // (2) fit model
  AliVdMPileup pileupModel;

  TString pn = TString::Format("pileup_%s.pdf", classAC);
  TCanvas *c1 = new TCanvas;
  c1->SaveAs(pn+"[");

  const AliTriggerBCMask& bcMask = vdmMetaData.GetTriggerBCMask();
  const Int_t             nBCs   = bcMask.GetNUnmaskedBCs();

  TH1 *hPar[5] = {
    SetAttr(new TH1D("hrA",      ";BCID;r_{A}",           nBCs,0,nBCs), kRed),
    SetAttr(new TH1D("hrC",      ";BCID;r_{C}",           nBCs,0,nBCs), kBlue),
    SetAttr(new TH1D("hbkgdA",   ";BCID;bkgd_{A}",        nBCs,0,nBCs), kRed),
    SetAttr(new TH1D("hbkgdC",   ";BCID;bkgd_{C}",        nBCs,0,nBCs), kBlue),
    SetAttr(new TH1D("hChi2NDF", ";BCID;#chi^{2}/n.d.f.", nBCs,0,nBCs), kBlue)
  };

  for (Int_t bc=0, counter=0; bc<3564; ++bc) {
    if (bcMask.GetMask(bc))
      continue;

    const TString binLabel = TString::Format("%d", bc);
    for (Int_t i=0; i<5; ++i)
      hPar[i]->GetXaxis()->SetBinLabel(1+counter, binLabel);

    c1->Clear();
    c1->SetLogx();
    c1->SetLogy();
    TH1 *hf = c1->DrawFrame(1e-6, 0.01, 0.5, 20);
    hf->SetTitle(TString::Format("BCID=%d %s;two-arm #mu;one-arm/two-arm", bc, classAC));
    pileupModel.DoFit(&gAnotC[bc], &gCnotA[bc], &par0[0]);

    SetAttr(&gAnotC[bc], kBlue);
    SetAttr(&gCnotA[bc], kRed);
    gAnotC[bc].Draw("PE");
    gCnotA[bc].Draw("PE");

    TF1 *fAnotC = SetAttr(new TF1("fAnotC", &pileupModel, &AliVdMPileup::fcnAnotC, 1e-6, 0.5, 5), kBlue);
    fAnotC->SetParameters(pileupModel.GetPar());
    fAnotC->SetNpx(1000);
    fAnotC->Draw("same");

    TF1 *fCnotA = SetAttr(new TF1("fCnotA", &pileupModel, &AliVdMPileup::fcnCnotA, 1e-6, 0.5, 5), kRed);
    fCnotA->SetParameters(pileupModel.GetPar());
    fCnotA->SetNpx(1000);
    fCnotA->Draw("same");

    TLegend *leg = new TLegend(0.6, 0.75, 0.9, 0.9);
    leg->AddEntry(&gCnotA[bc], "AnotC/AandC", "PEL");
    leg->AddEntry(&gAnotC[bc], "CnotA/AandC", "PEL");
    leg->Draw();

    TPaveText *pt = new TPaveText(0.6, 0.4, 0.9, 0.7, "NDC NB");
    pt->SetFillStyle(0);
    pt->AddText(TString::Format("#chi^{2}/n.d.f = %.0f/%.0f = %.2f", pileupModel.GetChi2(), pileupModel.GetNDF(), pileupModel.GetChi2()/pileupModel.GetNDF()));
    {
      double curval,err, lowlim, uplim;
//.........这里部分代码省略.........