C++ TRandom::Uniform方法代码示例

void makePoints(Int_t n, Double_t *x, Double_t *y, Double_t *e, Int_t p)
{
  Int_t i;
  TRandom r;

  if (p==2) {
    for (i=0; i<n; i++) {
      x[i] = r.Uniform(-2, 2);
      y[i]=TMath::Sin(x[i]) + TMath::Sin(2*x[i]) + r.Gaus()*0.1;
      e[i] = 0.1;
    }
  }
  if (p==3) {
    for (i=0; i<n; i++) {
      x[i] = r.Uniform(-2, 2);
      y[i] = -7 + 2*x[i]*x[i] + x[i]*x[i]*x[i]+ r.Gaus()*0.1;
      e[i] = 0.1;
    }
  }
  if (p==4) {
    for (i=0; i<n; i++) {
      x[i] = r.Uniform(-2, 2);
      y[i]=-2 + TMath::Exp(-x[i]) + r.Gaus()*0.1;
      e[i] = 0.1;
    }
  }
}

double write(int n) {



  TRandom R;
  TStopwatch timer;


  TFile f1("mathcoreLV.root","RECREATE");

  // create tree
  TTree t1("t1","Tree with new LorentzVector");

  std::vector<ROOT::Math::XYZTVector>  tracks;
  std::vector<ROOT::Math::XYZTVector> * pTracks = &tracks;
  t1.Branch("tracks","std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > >",&pTracks);

  double M = 0.13957;  // set pi+ mass

  timer.Start();
  double sum = 0;
  for (int i = 0; i < n; ++i) {
    int nPart = R.Poisson(5);
    pTracks->clear();
    pTracks->reserve(nPart);
    for (int j = 0; j < nPart; ++j) {
      double px = R.Gaus(0,10);
      double py = R.Gaus(0,10);
      double pt = sqrt(px*px +py*py);
      double eta = R.Uniform(-3,3);
      double phi = R.Uniform(0.0 , 2*TMath::Pi() );
      RhoEtaPhiVector vcyl( pt, eta, phi);
      // set energy
      double E = sqrt( vcyl.R()*vcyl.R() + M*M);
      XYZTVector q( vcyl.X(), vcyl.Y(), vcyl.Z(), E);
      // fill track vector
      pTracks->push_back(q);
      // evaluate sum of components to check
      sum += q.x()+q.y()+q.z()+q.t();
    }
    t1.Fill();
  }

  f1.Write();
  timer.Stop();
  std::cout << " Time for new Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;

  t1.Print();
  return sum;
}

 void Build() {
    TRandom rand;
    for (int i = 0; i < nelems; ++i) {
       fV0[i] = rand.Uniform();
       fV1[i] = rand.Uniform();
       fV2[i] = rand.Uniform();
       fV3[i] = rand.Uniform();
       fV4[i] = rand.Uniform();
       fV5[i] = rand.Uniform();
       fV6[i] = rand.Uniform();
       fV7[i] = rand.Uniform();
       fV8[i] = rand.Uniform();
       fV9[i] = rand.Uniform();
    }
 }

void fmpmt(void){

  TTree *t = new TTree("tree","fine mesh");
  double adc;
  t->Branch("adc",&adc,"adc/D");

  double Mult(double seed){
    do{
      double _fac=2.79360;
      TRandom _rand;
      _rand.SetSeed(0);
      double _judge = _rand.Uniform(0,1);
      double _x = _rand.Uniform(0.001,10);
    } while (TMath::Gaus(_x, _fac, _fac/20) < _judge);
    return seed*_x;
  }


  double prob = 0.7;
  double npe = 0.1;
  double judge, tmp;
  TRandom rand;
  rand.SetSeed(0);

  for(int i=0; i<4000; ++i){
    if(i%100==0){
      cout<<"=====  "<<i<<"  ====="<<endl;
    }
    do{
      tmp = rand.Uniform(0,4);
      judge = rand.Uniform(0,1);
    } while (TMath::Poisson(tmp, npe) < judge);
    for(int j=0; j<19; ++j){
      if(j==0){
        judge = rand.Uniform(0,1);
        if(judge>prob){
          continue;
        }
      }
      tmp=Mult(tmp);
    }
    adc=tmp;
    t->Fill();
  }
  t->Draw("adc");


}

void graph2derrorsfit()
{
   TCanvas *c1 = new TCanvas("c1");

   Double_t rnd, x, y, z, ex, ey, ez;
   Double_t e = 0.3;
   Int_t nd = 500;

   TRandom r;
   TF2  *f2 = new TF2("f2","1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200",-6,6,-6,6);
   f2->SetParameters(1,1);
   TGraph2DErrors *dte = new TGraph2DErrors(nd);

   // Fill the 2D graph
   for (Int_t i=0; i<nd; i++) {
      f2->GetRandom2(x,y);      
      rnd = r.Uniform(-e,e); // Generate a random number in [-e,e]
      z = f2->Eval(x,y)*(1+rnd);
      dte->SetPoint(i,x,y,z);
      ex = 0.05*r.Rndm();
      ey = 0.05*r.Rndm();
      ez = TMath::Abs(z*rnd);
      dte->SetPointError(i,ex,ey,ez);
   }

   f2->SetParameters(0.5,1.5);
   dte->Fit(f2);
   TF2 *fit2 = (TF2*)dte->FindObject("f2");
   fit2->SetTitle("Minuit fit result on the Graph2DErrors points");
   fit2->Draw("surf1");
   dte->Draw("axis p0");
}

void multicolor(Int_t isStack=0) {
   TCanvas *c1 = new TCanvas;
   Int_t nbins = 20;
   TH2F *h1 = new TH2F("h1","h1",nbins,-4,4,nbins,-4,4);
   h1->SetFillColor(kBlue);
   TH2F *h2 = new TH2F("h2","h2",nbins,-4,4,nbins,-4,4);
   h2->SetFillColor(kRed);
   TH2F *h3 = new TH2F("h3","h3",nbins,-4,4,nbins,-4,4);
   h3->SetFillColor(kYellow);
   THStack *hs = new THStack("hs","three plots");
   hs->Add(h1);
   hs->Add(h2);
   hs->Add(h3);
   TRandom r;
   Int_t i;
   for (i=0;i<20000;i++) h1->Fill(r.Gaus(),r.Gaus());
   for (i=0;i<200;i++) {
      Int_t ix = (Int_t)r.Uniform(0,nbins);
      Int_t iy = (Int_t)r.Uniform(0,nbins);
      Int_t bin = h1->GetBin(ix,iy);
      Double_t val = h1->GetBinContent(bin);
      if (val <= 0) continue;
      if (!isStack) h1->SetBinContent(bin,0);
      if (r.Rndm() > 0.5) {
         if (!isStack) h2->SetBinContent(bin,0);
         h3->SetBinContent(bin,val);
      } else {
         if (!isStack) h3->SetBinContent(bin,0);
         h2->SetBinContent(bin,val);
      }
   }
   hs->Draw("lego1");
}

void spheres(Int_t nspheres=20, Int_t npoints=100000) {
   // generate random points uniformly distributed over the surface
   // of spheres generated at random positions.

   TCanvas *c1 = new TCanvas("c1","spheres",600,600);
   c1->SetFillColor(kBlack);
   TView *view = new TView(1);
   Double_t k=0.8;
   view->SetRange( -k, -k, -k, k, k, k);
   
   //generate sphere coordinates and radius
   TRandom r;
   if (nspheres <=0) nspheres = 10;
   Double_t *xs = new Double_t[nspheres];
   Double_t *ys = new Double_t[nspheres];
   Double_t *zs = new Double_t[nspheres];
   Double_t *rs = new Double_t[nspheres];
   Int_t i;
   for (i=0;i<nspheres;i++) {
      xs[i] = r.Uniform(-1,1);
      ys[i] = r.Uniform(-1,1);
      zs[i] = r.Uniform(-1,1);
      rs[i] = r.Uniform(0.05,0.25);
   }

   //generate points
   TPolyMarker3D *pm = new TPolyMarker3D(npoints);
   pm->SetMarkerColor(kRed);
   Double_t x,y,z;
   for (Int_t j=0;j<npoints;j++) {
       i = (Int_t)r.Uniform(0,nspheres);  //choose sphere
       r.Sphere(x,y,z,rs[i]);
       pm->SetPoint(j,xs[i]+x, ys[i]+y,zs[i]+z);
   }
   
   delete [] xs;
   delete [] ys;
   delete [] zs;
   delete [] rs;
    
   pm->Draw();
}

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

void triangles(Int_t ntriangles=50) {
  TCanvas *c1 = new TCanvas("c1","triangles",10,10,700,700);
  TRandom r;
  Double_t dx = 0.2; Double_t dy = 0.2;
  Int_t ncolors = gStyle->GetNumberOfColors();
  Double_t x[4],y[4];
  TColor *c;
  Int_t ci;
  for (Int_t i=0;i<ntriangles;i++) {
     x[0] = r.Uniform(.05,.95); y[0] = r.Uniform(.05,.95);
     x[1] = x[0] + dx*r.Rndm(); y[1] = y[0] + dy*r.Rndm();
     x[2] = x[1] - dx*r.Rndm(); y[2] = y[1] - dy*r.Rndm();
     x[3] = x[0];               y[3] = y[0];
     TPolyLine *pl = new TPolyLine(4,x,y);
     pl->SetUniqueID(i);
     ci = ncolors*r.Rndm();
     c  = gROOT->GetColor(TColor::GetColorPalette(ci));
     c->SetAlpha(r.Rndm());
     pl->SetFillColor(ci);
     pl->Draw("f");
  }
  c1->AddExec("ex","TriangleClicked()");
}

void generate_imt_tree()
{
  // Create the file and the tree
  TFile hfile("ttree_read_imt.root", "RECREATE", "File for IMT test");
  TTree tree("TreeIMT", "TTree for IMT test");

  int nvbranches = 50, nabranches = 50;
  int nentries = 1000, nvelems = 100;
  std::vector<std::vector<Double_t>> vbranches(nvbranches);
  std::vector<A> abranches(nabranches);

  // Create the tree branches
  for (int i = 0; i < nvbranches; ++i) {
    vbranches[i] = std::vector<Double_t>(nvelems);

    std::string branchname("Vbranch");
    branchname += std::to_string(i);
    branchname += std::string("."); // make the top-level branch name be included in the sub-branch names

    tree.Branch(branchname.c_str(), &vbranches[i]);
  }
  for (int i = 0; i < nabranches; ++i) {
    std::string branchname("Abranch");
    branchname += std::to_string(i);
    branchname += std::string("."); // make the top-level branch name be included in the sub-branch names

    tree.Branch(branchname.c_str(), &abranches[i]);
  }

  // Fill the tree
  TRandom rand;
  for (int i = 0; i < nentries; i++) {
    for (int i = 0; i < nvbranches; ++i) {
      for (int j = 0; j < nvelems; ++j) {
        vbranches[i][j] = rand.Uniform();
      }
    }
    for (int i = 0; i < nabranches; ++i) {
      abranches[i].Build();
    }
    Int_t nb = tree.Fill();
  }

  // Write the file 
  hfile.Write();
  hfile.Close();
}

void apply_fake_rate(TRandom & rng, hxx_tree & data, double rate) {
    data.weight *= rate;

    //if (data.l1_pt <= 0.0) return;
    //if (data.l2_pt > 0.0) return;

    int n = data.jet_pt->size();

    if (n <= 0) return;

    int i = (int) (rng.Uniform() * n);
    //cout << i << " of " << n << "\n";

    double pt = data.jet_pt->at(i);

    if (pt > data.l1_pt) {
        data.l2_pt  = data.l1_pt;
        data.l2_eta = data.l1_eta;
        data.l2_phi = data.l1_phi;
        data.l1_pt  = pt;
        data.l1_eta = data.jet_eta->at(i);
        data.l1_phi = data.jet_phi->at(i);
    } else {
        data.l2_pt  = pt;
        data.l2_eta = data.jet_eta->at(i);
        data.l2_phi = data.jet_phi->at(i);
    }

    data.erase_jet(i);

    if (data.l2_pt <= 0.0) {
        data.mll = 0.0;
        return;
    }

    TLorentzVector v, v1, v2;
    v1.SetPtEtaPhiM(data.l1_pt, data.l1_eta, data.l1_phi, 0.0);
    v2.SetPtEtaPhiM(data.l2_pt, data.l2_eta, data.l2_phi, 0.0);
    v = v1 + v2;
    data.mll = v.M();

    //cout << data.mll << "\n";
    //cout << "   " << data.l1_pt << "\n";
    //cout << "   " << data.l2_pt << "\n";

}

double sampleXgenFromFunction1d(TRandom& rnd, TF1* function, double xMin, double xMax, double pMax)
{
  //std::cout << "<sampleXgenFromFunction1d>:" << std::endl;
  //std::cout << " function = " << function->GetTitle() << std::endl;
  //int numParameter = function->GetNpar();
  //for ( int iParameter = 0; iParameter < numParameter; ++iParameter ) {
  //  std::cout << " parameter #" << iParameter << " = " << function->GetParameter(iParameter) << std::endl;
  //}
  //assert(0);
  double x;
  bool isDone = false;
  while ( !isDone ) {
    x = rnd.Uniform(xMin, xMax);
    double u = rnd.Rndm();
    double f_x = function->Eval(x);
    //std::cout << "x = " << x << ": f(x) = " << f_x << " (pMax = " << pMax << ")" << std::endl;
    isDone = ((u*pMax) < f_x);
  }
  return x;
}

void test_mc_rejection_method_lin()
{
  const int    N     = 500000;
  const int    nbins = 300;
  const double ymax  = 1.1;
  const double xmin  = 0;
  const double xmax  = 10;
  const double dx    = xmax-xmin;

  TRandom rg;

  TF1 *  func = new TF1  ("func","1/(x+1)",0,10);
  TH1D * hgen = new TH1D ("hgen","generated",nbins,xmin,xmax);

  for(int i=0; i<N; i++) {
     cout << "..................." << i << endl;
     bool selected=false;
     while(1) {
       double xg = xmin + rg.Uniform() * dx;
       double yg = ymax * rg.Uniform();
       double yc = func->Eval(xg);
       selected = (yg<yc);
       if(selected) {
           hgen->Fill(xg);
           break;
       }
     }
  }

  double IF = func->Integral(xmin,xmax);
  double IH = hgen->Integral("width");
  double sc = IF/IH;
  hgen->Scale(sc);

  hgen->Draw();
  func->Draw("same");
}

TLorentzVector* applyResolution(TLorentzVector *l_gen){

  if(randomForSmearing.Uniform()<.9){
    double l_Perp =
      l_gen->Perp()+(randomForSmearing.Gaus(0,0.012*1.15*
					    l_gen->Perp()+0.00000*1.15* l_gen->Perp()* l_gen->Perp() ));
    double l_Theta = l_gen->Theta()+(randomForSmearing.Gaus(0,0.001));
    double l_Phi = l_gen->Phi()+(randomForSmearing.Gaus(0,0.001));
  }else{
    double l_Perp =
      l_gen->Perp()+randomForSmearing.Gaus(-l_gen->Perp()*.04,l_gen->Perp()*.08);
    //double l_Perp = l_gen->Perp();
    double l_Theta = l_gen->Theta()+(randomForSmearing.Gaus(0,0.001));
    double l_Phi = l_gen->Phi()+(randomForSmearing.Gaus(0,0.001));
  }
  double l_Px = l_Perp*cos(l_Phi);
  double l_Py = l_Perp*sin(l_Phi);
  double l_Pz = l_Perp/tan(l_Theta);
  double l_E = sqrt(l_Px*l_Px+l_Py*l_Py+l_Pz*l_Pz);

  return (new TLorentzVector( l_Px, l_Py, l_Pz, l_E ));


}

void FitterUtilsSimultaneousExpOfPolyTimesX::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma,
      RooRealVar& nKemu, RooRealVar& nSignal, RooRealVar& nPartReco,
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne,
      RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma,
      RooRealVar& l1Kee, RooRealVar& l2Kee, RooRealVar& l3Kee, RooRealVar& l4Kee, RooRealVar& l5Kee,
      RooRealVar& l1Kemu, RooRealVar& l2Kemu, RooRealVar& l3Kemu, RooRealVar& l4Kemu, RooRealVar& l5Kemu,
      RooRealVar const& l1KeeGen, RooRealVar const& l2KeeGen, RooRealVar const& l3KeeGen, RooRealVar const& l4KeeGen, RooRealVar const& l5KeeGen 

      )
{
  FitterUtilsExpOfPolyTimesX::initiateParams(nGenSignalZeroGamma, nGenSignalOneGamma, nGenSignalTwoGamma,
           nSignal, nPartReco, nComb, fracZero, fracOne, nJpsiLeak, constPartReco, fracPartRecoSigma,
           l1Kee, l2Kee, l3Kee, l4Kee, l5Kee,
           l1KeeGen, l2KeeGen, l3KeeGen, l4KeeGen, l5KeeGen ); 



  TRandom rand;
  rand.SetSeed();

  nKemu.setVal(rand.Uniform(nGenKemu-5*sqrt(nGenKemu), nGenKemu+5*sqrt(nGenKemu)));
  nKemu.setRange(nGenKemu-10*sqrt(nGenKemu), nGenKemu+10*sqrt(nGenKemu));

  l1Kemu.setVal(rand.Uniform( l1KeeGen.getVal() - 5*l1KeeGen.getError(), l1KeeGen.getVal() + 5*l1KeeGen.getError() ) );
  l1Kemu.setRange( l1KeeGen.getVal() - 10*l1KeeGen.getError(), l1KeeGen.getVal() + 10*l1KeeGen.getError() );

  l2Kemu.setVal(rand.Uniform( l2KeeGen.getVal() - 5*l2KeeGen.getError(), l2KeeGen.getVal() + 5*l2KeeGen.getError() ) );
  l2Kemu.setRange( l2KeeGen.getVal() - 10*l2KeeGen.getError(), l2KeeGen.getVal() + 10*l2KeeGen.getError() );

  l3Kemu.setVal(rand.Uniform( l3KeeGen.getVal() - 5*l3KeeGen.getError(), l3KeeGen.getVal() + 5*l3KeeGen.getError() ) );
  l3Kemu.setRange( l3KeeGen.getVal() - 10*l3KeeGen.getError(), l3KeeGen.getVal() + 10*l3KeeGen.getError() );

  l4Kemu.setVal(rand.Uniform( l4KeeGen.getVal() - 5*l4KeeGen.getError(), l4KeeGen.getVal() + 5*l4KeeGen.getError() ) );
  l4Kemu.setRange( l4KeeGen.getVal() - 10*l4KeeGen.getError(), l4KeeGen.getVal() + 10*l4KeeGen.getError() );

  l5Kemu.setVal(rand.Uniform( l5KeeGen.getVal() - 5*l5KeeGen.getError(), l5KeeGen.getVal() + 5*l5KeeGen.getError() ) );
  l5Kemu.setRange( l5KeeGen.getVal() - 10*l5KeeGen.getError(), l5KeeGen.getVal() + 10*l5KeeGen.getError() );

}