C++ TH1D::FillRandom方法代码示例

void  makeGaussianSignals(SigData_t& m_sigdata)
{
  //std::vector<TH1D *> vgsh(NUMCHAN);
  std::vector<TH1D *> vcdfh(NUMCHAN);

  if( m_sigdata.find("gs") == m_sigdata.end() ) {
    cerr << "Gaussian signal data not found, not making CDF signal!" << endl;
    return;
  }

  for (int ichan=0; ichan<NUMCHAN; ichan++) {
    TH1D *cdfh;

    TString channame(channames[ichan]);
    TString name;

    TH1D * gsh = m_sigdata["gs"].at(ichan);

    assert(gsh) ;

#if 0
    name = "Signalgs_"+channame;
    gsh = (TH1D *)tch->Clone(name.Data());

    assert(gsh);

    gsh->SetTitle("Gaussian signal");
    
    gsh->Reset();

    TF1 *g = (TF1 *)gROOT->GetFunction("gaus");
    g->SetParameters(1,gaussian_mean_mass_gev,gaussian_mass_sigma_gev);
    gsh->FillRandom("gaus",100000);

    // norm to 1pb signal with 1/fb integrated luminosity
    double norm = 1000 * gseffxacc[ichan]/gsh->Integral(0,gsh->GetNbinsX()+1,"width");

    //gsh->Scale(norm/eff_fudgefactor); // kludge: pre-undo the fudge in the next module
    gsh->Scale(norm);

    vgsh[ichan] = gsh;
#endif

    // New CDF bump, same as Gauss but set to CDF (obs/theor)*(LHC theor) = 3.43pb
    cdfh = (TH1D *)gsh->Clone("CDFbump");

    cdfh->Scale(3.43);

    vcdfh[ichan] = cdfh;

    cdfh->Draw();

    gsh->Draw("same");

  } // channel loop

  //m_sigdata["gs"]  = vgsh;
  m_sigdata["cdf"] = vcdfh;
}                                                           // makeGaussianSignals

//=== FUNCTION DEFINITIONS ======================================================================================
TTree* toyHist(int iNToys, bool iDoMu,std::string iName,TH1D *iData,TH1D *iW,TH1D *iEWK,TH1D *iAntiData,TH1D *iAntiW,TH1D *iAntiEWK,double *iOriginalResults) {
    TCanvas *c = 0;
    TTree *tree = new TTree(("Toys"+iName).c_str(),("Toys"+iName).c_str());
    float *vals = new float[16];
    for(int i0 = 0; i0 < 16; i0++) vals[i0] = 0.;
    tree->Branch("nsig"      ,&vals[0] ,"Val0/F");
    tree->Branch("newk"      ,&vals[1] ,"Val1/F");
    tree->Branch("nqcd"      ,&vals[2] ,"Val2/F");
    tree->Branch("nantisig"  ,&vals[3] ,"Val3/F");
    tree->Branch("nantiewk"  ,&vals[4] ,"Val4/F");
    tree->Branch("nantiqcd"  ,&vals[5] ,"Val5/F");
    tree->Branch("sigma"     ,&vals[6] ,"Val6/F");
    tree->Branch("a1"        ,&vals[7] ,"Val7/F");
    tree->Branch("nsig_e"    ,&vals[8] ,"Val8/F");
    tree->Branch("newk_e"    ,&vals[9] ,"Val9/F");
    tree->Branch("nqcd_e"    ,&vals[10],"Val10/F");
    tree->Branch("nantisig_e",&vals[11],"Val11/F");
    tree->Branch("nantiewk_e",&vals[12],"Val12/F");
    tree->Branch("nantiqcd_e",&vals[13],"Val13/F");
    tree->Branch("sigma_e"   ,&vals[14],"Val14/F");
    tree->Branch("a1_e"      ,&vals[15],"Val15/F");
    //Set The first entry of the tree to be the default values
    for(int i1 = 0; i1 < 16; i1++) vals[i1] = float(iOriginalResults[i1]);
    tree->Fill();

    TH1D *pData       = 0;
    TH1D *pAntiData   = 0;
    for(int i0 = 0; i0 < iNToys; i0++) {
        std::stringstream pSS;
        pSS << "Hist" << i0;
        pData       = (TH1D*) iData    ->Clone((pSS.str()+"A").c_str());
        pAntiData   = (TH1D*) iAntiData->Clone((pSS.str()+"B").c_str());
        for(int i1 = 0; i1 < pData    ->GetNbinsX()+1; i1++) pData    ->SetBinContent(i1,0);
        for(int i1 = 0; i1 < pAntiData->GetNbinsX()+1; i1++) pAntiData->SetBinContent(i1,0);
        pData    ->FillRandom(iData,    iData    ->Integral());
        pAntiData->FillRandom(iAntiData,iAntiData->Integral());
        cout <<"==> " << pData->Integral() << " -- " << iData->Integral() << endl;
        double* pVals = fitHist(c,iDoMu,0,pSS.str(),pData,iW,iEWK,pAntiData,iAntiW,iAntiEWK,100,100,100,100);
        for(int i1 = 0; i1 < 16; i1++) vals[i1] = float(pVals[i1]);
        tree->Fill();
        delete pData;
        delete pAntiData;
    }
    return tree;
}

void combinedFit() {

  TH1D * hB = new TH1D("hB","histo B",100,0,100);
  TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);

  TF1 * fB = new TF1("fB","expo",0,100);
  fB->SetParameters(1,-0.05);
  hB->FillRandom("fB");

  TF1 * fS = new TF1("fS","gaus",0,100);
  fS->SetParameters(1,30,5);

  hSB->FillRandom("fB",2000);
  hSB->FillRandom("fS",1000);

  // perform now global fit

  TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);

  ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
  ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);

  ROOT::Fit::DataOptions opt;
  ROOT::Fit::DataRange rangeB;
  // set the data range
  rangeB.SetRange(10,90);
  ROOT::Fit::BinData dataB(opt,rangeB);
  ROOT::Fit::FillData(dataB, hB);

  ROOT::Fit::DataRange rangeSB;
  rangeSB.SetRange(10,50);
  ROOT::Fit::BinData dataSB(opt,rangeSB);
  ROOT::Fit::FillData(dataSB, hSB);

  ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
  ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);

  GlobalChi2 globalChi2(chi2_B, chi2_SB);

  ROOT::Fit::Fitter fitter;

  const int Npar = 6;
  double par0[Npar] = { 5,5,-0.1,100, 30,10};

  // create before the parameter settings in order to fix or set range on them
  fitter.Config().SetParamsSettings(6,par0);
  // fix 5-th parameter
  fitter.Config().ParSettings(4).Fix();
  // set limits on the third and 4-th parameter
  fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4);
  fitter.Config().ParSettings(3).SetLimits(0,10000);
  fitter.Config().ParSettings(3).SetStepSize(5);

  fitter.Config().MinimizerOptions().SetPrintLevel(0);
  fitter.Config().SetMinimizer("Minuit2","Migrad");

  // fit FCN function directly
  // (specify optionally data size and flag to indicate that is a chi2 fit)
  fitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms",
                             10,10,700,700);
  c1->Divide(1,2);
  c1->cd(1);
  gStyle->SetOptFit(1111);

  fB->SetFitResult( result, iparB);
  fB->SetRange(rangeB().first, rangeB().second);
  fB->SetLineColor(kBlue);
  hB->GetListOfFunctions()->Add(fB);
  hB->Draw();

  c1->cd(2);
  fSB->SetFitResult( result, iparSB);
  fSB->SetRange(rangeSB().first, rangeSB().second);
  fSB->SetLineColor(kRed);
  hSB->GetListOfFunctions()->Add(fSB);
  hSB->Draw();


}

void first_script() {
  TH1D* hist = new TH1D("myhist", "Gaussian Histogram (#sigma = 1)", 50, -5, 5);
  hist->FillRandom("gaus", 10000);
  hist->Draw();
}

void test()
{
//Illustrates TVirtualFitter::GetConfidenceIntervals
//This method computes confidence intervals for the fitted function
//Author: Anna Kreshuk

   TCanvas *myc = new TCanvas("myc",
      "Confidence intervals on the fitted function",1200, 500);
   myc->Divide(3,1);

/////1. A graph
   //Create and fill a graph
   Int_t ngr = 100;
   TGraph *gr = new TGraph(ngr);
   gr->SetName("GraphNoError");
   Double_t x, y;
   Int_t i;
   for (i=0; i<ngr; i++){
      x = gRandom->Uniform(-1, 1);
      y = -1 + 2*x + gRandom->Gaus(0, 1);
      gr->SetPoint(i, x, y);
   }
   //Create the fitting function
   TF1 *fpol = new TF1("fpol", "pol1", -1, 1);
   fpol->SetLineWidth(2);
   gr->Fit(fpol, "Q");

   //Create a TGraphErrors to hold the confidence intervals
   TGraphErrors *grint = new TGraphErrors(ngr);
   grint->SetTitle("Fitted line with .95 conf. band");
   for (i=0; i<ngr; i++)
      grint->SetPoint(i, gr->GetX()[i], 0);
   //Compute the confidence intervals at the x points of the created graph
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint);
   //Now the "grint" graph contains function values as its y-coordinates
   //and confidence intervals as the errors on these coordinates
   //Draw the graph, the function and the confidence intervals
   myc->cd(1);
   grint->SetLineColor(kRed);
   grint->Draw("ap");
   gr->SetMarkerStyle(5);
   gr->SetMarkerSize(0.7);
   gr->Draw("psame");

/////2. A histogram
   myc->cd(2);
   //Create, fill and fit a histogram
   Int_t nh=5000;
   TH1D *h = new TH1D("h",
      "Fitted gaussian with .95 conf.band", 100, -3, 3);
   h->FillRandom("gaus", nh);
   TF1 *f = new TF1("fgaus", "gaus", -3, 3);
   f->SetLineWidth(2);
   h->Fit(f, "Q");
   h->Draw();

   //Create a histogram to hold the confidence intervals
   TH1D *hint = new TH1D("hint",
      "Fitted gaussian with .95 conf.band", 100, -3, 3);
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(hint);
   //Now the "hint" histogram has the fitted function values as the
   //bin contents and the confidence intervals as bin errors
   hint->SetStats(kFALSE);
   hint->SetFillColor(2);
   hint->Draw("e3 same");

/////3. A 2d graph
   //Create and fill the graph
   Int_t ngr2 = 100;
   Double_t z, rnd, e=0.3;
   TGraph2D *gr2 = new TGraph2D(ngr2);
   gr2->SetName("Graph2DNoError");
   TF2  *f2 = new TF2("f2",
      "1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+250",-6,6,-6,6);
   f2->SetParameters(1,1);
   for (i=0; i<ngr2; i++){
      f2->GetRandom2(x,y);
      // Generate a random number in [-e,e]
      rnd = 2*gRandom->Rndm()*e-e;
      z = f2->Eval(x,y)*(1+rnd);
      gr2->SetPoint(i,x,y,z);
   }
   //Create a graph with errors to store the intervals
   TGraph2DErrors *grint2 = new TGraph2DErrors(ngr2);
   for (i=0; i<ngr2; i++)
      grint2->SetPoint(i, gr2->GetX()[i], gr2->GetY()[i], 0);

   //Fit the graph
   f2->SetParameters(0.5,1.5);
   gr2->Fit(f2, "Q");
   //Compute the confidence intervals
   (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint2);
   //Now the "grint2" graph contains function values as z-coordinates
   //and confidence intervals as their errors
   //draw
   myc->cd(3);
   f2->SetNpx(30);
   f2->SetNpy(30);
   f2->SetFillColor(kBlue);
   f2->Draw("surf4");
//.........这里部分代码省略.........

//___________________________________________________________________________
Double_t* Ifit( TFile* fin, std::string name, std::string output=".", std::string xTitle="", std::string yTitle="Events", int rebin=1, int channel=0, int fit_data=1 )
{
    Double_t* fitted = new Double_t[8];

    TCanvas *c1 = new TCanvas("HF1", "Histos1", 258,92,748,702);
    c1->Range(-104.4905,-2560.33,537.9965,11563.46);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetBorderSize(2);
    c1->SetLeftMargin(0.1626344);
    c1->SetRightMargin(0.05913978);
    c1->SetTopMargin(0.05349183);
    c1->SetBottomMargin(0.1812779);
    c1->SetFrameBorderMode(0);
    c1->SetFrameBorderMode(0);

    double count=0;
    dataColl.clear();
    sigColl.clear();
    bkgColl.clear();

    totalColl.clear();
    ctauColl.clear();

    //Get data from looping tree

    //TFile *fin = new TFile("results/15Dec_LepJet_MCDATA/Template_EvtChi2_Top_Hadronic_Mass.root");
    //TFile *fin = new TFile("results/15Dec_LepJet_MCDATA/Template_EvtChi2_Top_Leptonic_Mbl.root");


    //   TH1D *hsig = new TH1D();
    //   TH1D *hbkg = new TH1D();
    TH1D *hsig_toymc = new TH1D();
    TH1D *hbkg_toymc = new TH1D();

    char hname[30];
    std::string ch;
    if( channel == 1 )
        ch="_El";
    else if( channel == 2)
        ch="_Mu";
    else
        ch="";

    TH1D * hsig = (TH1D*)((TH1D*)fin->Get(("SigMC"+ch).c_str()))->Clone();
    TH1D * hbkg = (TH1D*)((TH1D*)fin->Get(("BkgMC"+ch).c_str()))->Clone();
    TH1D *hEGdata;
    hEGdata = (TH1D*)((TH1D*)fin->Get(("DATA"+ch).c_str()))->Clone();

    // hsig->Sumw2();
    //hbkg->Sumw2();

    // hsig->Rebin(2);
    // hbkg->Rebin(2);
    // hEGdata->Rebin(2);
    // hbkg_toymc->Rebin(2);
    // hsig_toymc->Rebin(2);

    hsig->Rebin(rebin);
    hbkg->Rebin(rebin);
    hEGdata->Rebin(rebin);

    // normalize template
    hsig->Scale(1./hsig->Integral());
    hbkg->Scale(1./hbkg->Integral());  
    if(fit_data==0){
        hsig_toymc->Scale(1./hsig_toymc->Integral());
        hbkg_toymc->Scale(1./hbkg_toymc->Integral());
    }

  TH1D *hsum = new TH1D();
  int ntemplate = 1000.;
  float sigfrac = 0.5;
  TH1D *hsum_norm = new TH1D();
  TH1D *hdata = new TH1D();

  int ndata=0;
  if ( fit_data>0 ) {
    hdata = (TH1D*)hEGdata->Clone();
    ndata = hdata->GetEntries();
  }else { //generate toymc
    hsum = (TH1D*)hsig_toymc->Clone();
    hsum->Scale(toymc_sig);
    hsum->Add(hbkg_toymc,toymc_bkg);
    
    hsum_norm = (TH1D*)hsum->Clone();  
    hsum_norm->Scale(1./hsum->Integral());
    hdata = (TH1D*)hsum_norm->Clone();
    //ndata = (int) gRandom->Poisson(hsum->Integral());
    ndata=toymc_sig+toymc_bkg;
    hdata->FillRandom(hsum_norm, ndata);
  }
  if(ndata==0) {
    printf(" ---  no events in the fit \n");
    fitted[0] = 0.;
    fitted[1] = 0.;
    fitted[2] = 0.;
    fitted[3] = 0.;
    fitted[4] = 0.;
//.........这里部分代码省略.........

void first_script3(int nbins, int nevents) {
  TH1D* hist =
      new TH1D("myhist", "Gaussian Histogram (#sigma = 1)", nbins, -5, 5);
  hist->FillRandom("gaus", nevents);
  hist->Draw();
}