C++ TRandom3::Gaus方法代码示例

void write_ntuple_to_file(){

    TFile ofile("conductivity_experiment.root","RECREATE");

    // Initialise the TNtuple
    TNtuple cond_data("cond_data",
                      "Example N-Tuple",
                      "Potential:Current:Temperature:Pressure");

    // Fill it randomly to fake the acquired data
    TRandom3 rndm;
    float pot,cur,temp,pres;
    for (int i=0;i<10000;++i){
        pot=rndm.Uniform(0.,10.);      // get voltage
        temp=rndm.Uniform(250.,350.);  // get temperature
        pres=rndm.Uniform(0.5,1.5);    // get pressure
        cur=pot/(10.+0.05*(temp-300.)-0.2*(pres-1.)); // current
// add some random smearing (measurement errors)
        pot*=rndm.Gaus(1.,0.01); // 1% error on voltage
        temp+=rndm.Gaus(0.,0.3); // 0.3 abs. error on temp.
        pres*=rndm.Gaus(1.,0.02);// 1% error on pressure
        cur*=rndm.Gaus(1.,0.01); // 1% error on current
// write to ntuple
        cond_data.Fill(pot,cur,temp,pres);
        }

    // Save the ntuple and close the file
    cond_data.Write();
    ofile.Close();
}

void macro6(){

    TH1F* sig_h=new TH1F("sig_h","Signal Histo",50,0,10);
    TH1F* gaus_h1=new TH1F("gaus_h1","Gauss Histo 1",30,0,10);
    TH1F* gaus_h2=new TH1F("gaus_h2","Gauss Histo 2",30,0,10);
    TH1F* bkg_h=new TH1F("exp_h","Exponential Histo",50,0,10);

    // simulate the measurements
    TRandom3 rndgen;
    for (int imeas=0;imeas<4000;imeas++){
        bkg_h->Fill(rndgen.Exp(4));
        if (imeas%4==0) gaus_h1->Fill(rndgen.Gaus(5,2));
        if (imeas%4==0) gaus_h2->Fill(rndgen.Gaus(5,2));
        if (imeas%10==0)sig_h->Fill(rndgen.Gaus(5,.5));}

    // Format Histograms
    TH1F* histos[4]={sig_h,bkg_h,gaus_h1,gaus_h2};
    for (int i=0;i<4;++i){
        histos[i]->Sumw2(); // *Very* Important
        format_h(histos[i],i+1);
        }

    // Sum
    TH1F* sum_h= new TH1F(*bkg_h);
    sum_h->Add(sig_h,1.);
    sum_h->SetTitle("Exponential + Gaussian");
    format_h(sum_h,kBlue);

    TCanvas* c_sum= new TCanvas();
    sum_h->Draw("hist");
    bkg_h->Draw("SameHist");
    sig_h->Draw("SameHist");

    // Divide
    TH1F* dividend=new TH1F(*gaus_h1);
    dividend->Divide(gaus_h2);

    // Graphical Maquillage
    dividend->SetTitle(";X axis;Gaus Histo 1 / Gaus Histo 2");
    format_h(dividend,kOrange);
    gaus_h1->SetTitle(";;Gaus Histo 1 and Gaus Histo 2");
    gStyle->SetOptStat(0);

    TCanvas* c_divide= new TCanvas();
    c_divide->Divide(1,2,0,0);
    c_divide->cd(1);
    c_divide->GetPad(1)->SetRightMargin(.01);
    gaus_h1->DrawNormalized("Hist");
    gaus_h2->DrawNormalized("HistSame");

    c_divide->cd(2);
    dividend->GetYaxis()->SetRangeUser(0,2.49);
    c_divide->GetPad(2)->SetGridy();
    c_divide->GetPad(2)->SetRightMargin(.01);
    dividend->Draw();
}

int main(int argc, char *argv[])
{
    std::auto_ptr<TRint> application(new TRint("histInMemory", 0, 0));

    TH1F *hist1 = new TH1F("hist1", "hist1", 50, 0, 10);
    TH1F *hist2 = new TH1F("hist2", "hist2", 50, 0, 10);
    TH1F *hist3 = new TH1F("hist3", "hist3", 50, 0, 10);
    TH1F *hist4 = new TH1F("hist4", "hist4", 50, 0, 10);
    TH1F *hist5 = new TH1F("hist5", "hist5", 50, 0, 10);

    {
        TRandom3 *random = new TRandom3(1);
        for(int i = 0; 10000 > i; ++i)
        {
            hist1->Fill(random->Gaus(5, 1));
            hist2->Fill(random->Poisson(5));
            hist3->Fill(random->Uniform(0, 10));
            hist4->Fill(random->Gaus(4, 1));
            hist5->Fill(random->Poisson(3));
        }
        delete random;
    }

    TFile *output = new TFile("output.root", "RECREATE");
    TDirectory *group1 = output->mkdir("group1");
    group1->cd();

    hist1->Write();
    hist2->Write();

    TDirectory *group2 = output->mkdir("group2");
    group2->cd();

    hist3->Write();
    hist4->Write();
    hist5->Write();

    delete output;

    application->Run(true);

    delete hist1;
    delete hist2;
    delete hist3;
    delete hist4;
    delete hist5;

    return 0;
}

//____________________________________________________________________
void fitCircle(Int_t n=10000) {
   //generates n points around a circle and fit them
   TCanvas *c1 = new TCanvas("c1","c1",600,600);
   c1->SetGrid();
   gr = new TGraph(n);
   if (n> 999) gr->SetMarkerStyle(1);
   else        gr->SetMarkerStyle(3);
   TRandom3 r;
   Double_t x,y;
   for (Int_t i=0;i<n;i++) {
      r.Circle(x,y,r.Gaus(4,0.3));
      gr->SetPoint(i,x,y);
   }
   c1->DrawFrame(-5,-5,5,5);
   gr->Draw("p");

   //Fit a circle to the graph points
   TVirtualFitter::SetDefaultFitter("Minuit");  //default is Minuit
   TVirtualFitter *fitter = TVirtualFitter::Fitter(0, 3);
   fitter->SetFCN(myfcn);

   fitter->SetParameter(0, "x0",   0, 0.1, 0,0);
   fitter->SetParameter(1, "y0",   0, 0.1, 0,0);
   fitter->SetParameter(2, "R",    1, 0.1, 0,0);

   Double_t arglist[1] = {0};
   fitter->ExecuteCommand("MIGRAD", arglist, 0);

   //Draw the circle on top of the points
   TArc *arc = new TArc(fitter->GetParameter(0),
      fitter->GetParameter(1),fitter->GetParameter(2));
   arc->SetLineColor(kRed);
   arc->SetLineWidth(4);
   arc->Draw();
}

void smearEmEnergy( const EnergyScaleCorrection_class& egcor, TLorentzVector& p, const ZGTree& myTree, float r9 ) {

    float smearSigma =  egcor.getSmearingSigma(myTree.run, fabs(p.Eta())<1.479, r9, p.Eta(), p.Pt(), 0., 0.);
    float cor = myRandom_.Gaus( 1., smearSigma );
    p.SetPtEtaPhiM( p.Pt()*cor, p.Eta(), p.Phi(), p.M() );

}

void plot_distortion(){
    
    TRandom3* rand = new TRandom3();

//    TF1 *func = new TF1("func","TMath::Abs([0]+[1]*x)",0,12);
//
//    rand->SetSeed(0);
//    Double_t slope=0.2*rand->Gaus(0,1);
//    Double_t anchor_point=3.5;
//    //want offset to be set by requiring func at anchor point to be 1
//    Double_t offset=(1-slope*anchor_point);
//    func->SetParameter(0,offset);
//    func->SetParameter(1,slope);
    
    //FN distortion
    TF1 *func = new TF1("func","[0]/(0.2*pow(x,0.1))+[1]",0,12);
    rand->SetSeed(0);
    Double_t scaling =0.2*rand->Gaus(0,1);
    func->SetParameter(0,scaling);
    func->SetParameter(1,0);
    //set offset so that func(FinalEnergy)=1;
    func->SetParameter(1,1-1*func->Eval(12));
    
    func->Draw();
    
    
   }

void gtime2(Int_t nsteps = 200, Int_t np=5000) {
   if (np > 5000) np = 5000;
   Int_t color[5000];
   Double_t cosphi[5000], sinphi[5000], speed[5000];
   TRandom3 r;
   Double_t xmin = 0, xmax = 10, ymin = -10, ymax = 10;
   TGraphTime *g = new TGraphTime(nsteps,xmin,ymin,xmax,ymax);
   g->SetTitle("TGraphTime demo 2;X;Y");
   Int_t i,s;
   Double_t phi,fact = xmax/Double_t(nsteps);
   for (i=0;i<np;i++) { //calculate some object parameters
      speed[i]  = r.Uniform(0.5,1);
      phi       = r.Gaus(0,TMath::Pi()/6.);
      cosphi[i] = fact*speed[i]*TMath::Cos(phi);
      sinphi[i] = fact*speed[i]*TMath::Sin(phi);
      Double_t rc = r.Rndm();
      color[i] = kRed;
      if (rc > 0.3) color[i] = kBlue;
      if (rc > 0.7) color[i] = kYellow;
   }
   for (s=0;s<nsteps;s++) { //fill the TGraphTime step by step
      for (i=0;i<np;i++) {
         Double_t xx = s*cosphi[i];
         if (xx < xmin) continue;
         Double_t yy = s*sinphi[i];
         TMarker *m = new TMarker(xx,yy,25);
         m->SetMarkerColor(color[i]);
         m->SetMarkerSize(1.5 -s/(speed[i]*nsteps));
         g->Add(m,s);
      }
      g->Add(new TPaveLabel(.70,.92,.98,.99,Form("shower at %5.3f nsec",3.*s/nsteps),"brNDC"),s);
   }
   g->Draw();
}

//____________________________________________________________________
void fitCircle(Int_t n=10000) {
   //generates n points around a circle and fit them
   TCanvas *c1 = new TCanvas("c1","c1",600,600);
   c1->SetGrid();
   gr = new TGraph(n);
   if (n> 999) gr->SetMarkerStyle(1);
   else        gr->SetMarkerStyle(3);
   TRandom3 r;
   Double_t x,y;
   for (Int_t i=0;i<n;i++) {
      r.Circle(x,y,r.Gaus(4,0.3));
      gr->SetPoint(i,x,y);
   }
   c1->DrawFrame(-5,-5,5,5);
   gr->Draw("p");


   auto chi2Function = [&](const Double_t *par) {
      //minimisation function computing the sum of squares of residuals
      // looping at the graph points
      Int_t np = gr->GetN();
      Double_t f = 0;
      Double_t *x = gr->GetX();
      Double_t *y = gr->GetY();
      for (Int_t i=0;i<np;i++) {
         Double_t u = x[i] - par[0];
         Double_t v = y[i] - par[1];
         Double_t dr = par[2] - std::sqrt(u*u+v*v);
         f += dr*dr;
      }
      return f;
   };

   // wrap chi2 funciton in a function object for the fit
   // 3 is the number of fit parameters (size of array par)
   ROOT::Math::Functor fcn(chi2Function,3);
   ROOT::Fit::Fitter  fitter;


   double pStart[3] = {0,0,1};
   fitter.SetFCN(fcn, pStart);
   fitter.Config().ParSettings(0).SetName("x0");
   fitter.Config().ParSettings(1).SetName("y0");
   fitter.Config().ParSettings(2).SetName("R");

   // do the fit 
   bool ok = fitter.FitFCN();
   if (!ok) {
      Error("line3Dfit","Line3D Fit failed");
   }   

   const ROOT::Fit::FitResult & result = fitter.Result();
   result.Print(std::cout);

   //Draw the circle on top of the points
   TArc *arc = new TArc(result.Parameter(0),result.Parameter(1),result.Parameter(2));
   arc->SetLineColor(kRed);
   arc->SetLineWidth(4);
   arc->Draw();
}

// Code from http://www.hongliangjie.com/2012/12/19/how-to-generate-gamma-random-variables/
// Parameter b could be theta...
double gsl_ran_gamma(const double a, const double b, TRandom3 &rand)
{
  
  if (a < 1)
  {
    double u = rand.Uniform(1);
    return gsl_ran_gamma(1.0 + a, b, rand) * pow (u, 1.0 / a);
  }

  double x, v, u;
  double d = a - 1.0 / 3.0;
  double c = (1.0 / 3.0) / sqrt (d);
  
  while (1) 
  {
    do 
    {
      x = rand.Gaus(0, 1.0);
      v = 1.0 + c * x;
    }
    while (v <= 0);
      
    v = v * v * v;
    u = rand.Uniform(1);

    if (u < 1 - 0.0331 * x * x * x * x) 
      break;

    if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
      break;
  }
    
  return b * d * v;
}

int main()
{
	// Manual input:
	if (true)
	{
		vector<KLV> gen(4); // pt eta phi m
		gen[1].p4.SetCoordinates(30,  1.0, 0.5, 0);
		gen[3].p4.SetCoordinates(25, -1.0, 2.5, 0);
		gen[2].p4.SetCoordinates(15,  2.1, 1.1, 0);
		gen[0].p4.SetCoordinates(10,  2.3, 1.0, 0);
		vector<KLV> rec(3);
		rec[2].p4.SetCoordinates(20,  1.3, 0.4, 0);
		rec[0].p4.SetCoordinates(18, -0.9, 2.2, 0);
		rec[1].p4.SetCoordinates(13,  2.2, 1.1, 0);
		test(gen, rec);
	}

	// Automatic input
	TRandom3 rnd;
	for (size_t i = 0; i < 4; ++i)
	{
		vector<KLV> gen, rec;
		// GEN stuff
		for (size_t j = 0; j < (size_t)rnd.Poisson(10); ++j)
		{
			gen.push_back(KLV());
			gen.back().p4.SetCoordinates(rnd.Uniform(5, 100), rnd.Gaus(0, 5), rnd.Uniform(-M_PI, +M_PI), rnd.Gaus(0, 10));
			// RECO efficiency
			if (rnd.Uniform(0, 1) > 0.05)
			{
				rec.push_back(KLV());
				rec.back().p4.SetCoordinates(
					gen.back().p4.pt() * fabs(rnd.Gaus(0.8, 0.1)),
					gen.back().p4.eta() * fabs(rnd.Gaus(1, 0.1)),
					gen.back().p4.phi() + rnd.Uniform(0, 0.1),
					gen.back().p4.mass() * rnd.Gaus(1, 0.1));
			}
		}
		// RECO noise
		for (size_t j = 0; j < (size_t)rnd.Poisson(5); ++j)
		{
			rec.push_back(KLV());
			rec.back().p4.SetCoordinates(rnd.Uniform(1, 10), rnd.Gaus(0, 5), rnd.Uniform(-M_PI, +M_PI), rnd.Gaus(0, 10));
		}
		test(gen, rec);
	}
}

int main() {

     //Teil a)
    int array[100];
    for(int i = 0; i < 100; i++) {
        array[i] = i+1;
        //cout << array[i] << endl;
    }


    //Teil b)
    double array2d[100][8];
    
    TRandom3 randomgen;
    for(int i = 0; i < 100; i++) {
        //0.Spalte
        array2d[i][0] = i+1;

        //1. Spalte
        array2d[i][1] = randomgen.Rndm();

        //2. Spalte
        array2d[i][2] = randomgen.Rndm() * 10;

        //3. Spalte
        array2d[i][3] = randomgen.Rndm() * 10 + 20;

        //4. Spalte
        array2d[i][4]= randomgen.Gaus();

        //5. Spalte
        array2d[i][5] = randomgen.Gaus(5, 2);

        //6. Spalte
        array2d[i][6] = array2d[i][0] * array2d[i][0];

        //7. Spalte
        array2d[i][7] = TMath::Cos(array2d[i][0]); 
    }


    return 0;
}

//ks-test for Gaussian Random samples
void test(Int_t num)
{
  Double_t *d=new Double_t[num];
  Double_t p,t;
  TRandom3 r;
  for (int i=0;i<num;i++)
    d[i]=r.Gaus();
  ROOT::Math::GoFTest *tt=new ROOT::Math::GoFTest(num,d,ROOT::Math::GoFTest::kGaussian);
  tt->KolmogorovSmirnovTest(p,t);
  cout<<"p:"<<p<<",t:"<<t<<endl;
}

int Ngenerated(double muIn)
{

	double sigma = sqrt(muIn + background);
	double mean = muIn + background;
	double nn = rnd.Gaus(mean, sigma);
	int ng = int(round(nn));

	return ng;

}

void generateGaussRandom(){

  TRandom3 *gr = new TRandom3(12345);
  cout<<"generating gaus rand for use "<< Nrand<<endl; 
  for(int j=0;j<Nrand;j++){
    if(j%100000000==0) cout<<" j " << j <<endl; 
    float tmp = gr->Gaus();
    vrandgaus.push_back(tmp);
  }
  
}

Double_t Reconstruct( Double_t xt, TRandom3& R )
{
   // apply some Gaussian smearing + bias and efficiency corrections to fake reconstruction
   const Double_t cutdummy = -99999.0;
   Double_t xeff = 0.3 + (1.0 - 0.3)/20.0*(xt + 10.0);  // efficiency
   Double_t x    = R.Rndm();
   if (x > xeff) return cutdummy;
   else {
     Double_t xsmear= R.Gaus(-2.5,0.2); // bias and smear
     return xt+xsmear;
   }
}