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

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

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 write(int n) {

   TRandom R;
   TStopwatch timer;

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

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

   XYZTVectorF *v1 = new XYZTVectorF();
   t1.Branch("LV branch","ROOT::Math::XYZTVectorF",&v1);

   timer.Start();
   for (int i = 0; i < n; ++i) {
      double Px = R.Gaus(0,10);
      double Py = R.Gaus(0,10);
      double Pz = R.Gaus(0,10);
      double E  = R.Gaus(100,10);
      v1->SetCoordinates(Px,Py,Pz,E);
      t1.Fill();
   }

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

// Throw toy experiments to predict observed yields:
// - Throw nPredictions mean values (prediction) from the background estimates
//   considering their uncertainties
//   - Per prediction, throw nExperiments observed yields from a Poisson
//     distribution with mean value prediction
void ToyExperiments::run(unsigned int nPredictions, unsigned int nExperiments) const {
  std::cout << "Predicting event yields in " << Parameters::nBins() << " bins from " << nPredictions*nExperiments << " toy experiments ...  " << std::flush;

  // Minimal value (in standard deviations) allowed for
  //correlated fluctuation
  const double minCorr = findMinValidRandomNumberForCorrelatedUncertainties();

  // Throw mean values
  for(unsigned int p = 0; p < nPredictions; ++p) {
    // Throw one (normalized) random number for correlated
    // uncertainties that is valid in all bins
    double rCorr = rand_->Gaus(0.,1.);
    while( rCorr <= minCorr ) {
      rCorr = rand_->Gaus(0.,1.);
    }

    // Loop over all bins and get individual predictions
    for(unsigned int bin = 0; bin < Parameters::nBins(); ++bin) {
      double prediction = -1.;
      bool negativePrediction = true;
      while( negativePrediction ) {
	// Throw one (normalized) random number for uncorrelated
	// uncertainties that is valid in this bin only
	double rUncorr = rand_->Gaus(0.,1.);
	
	// Scale the normalized random numbers by the uncertainties' size
	// to obtain variation of yield
	double uncorrVar = rUncorr * uncorrelatedUncerts_.at(bin);
	double corrVar   = rCorr   * correlatedUncerts_.at(bin);
	
	// Add variations to yield
	prediction = meanPredictions_.at(bin) + uncorrVar + corrVar;
	
	// Check if prediction is positive
	if( prediction >= 0. ) {
	  negativePrediction = false;
	}
      }

      // Throw predicted yields from Poisson with
      // this mean
      for(unsigned int e = 0; e < nExperiments; ++e) {
	predictedYields_.at(bin)->Fill(rand_->Poisson(prediction));
      }
    }
  }
  std::cout << "ok" << std::endl;

  for(unsigned int bin = 0; bin < Parameters::nBins(); ++bin) {
    if( predictedYields_.at(bin)->GetBinContent(Parameters::maxYield()+1) > 0 ) {
      std::cerr << "\n\nWARNING: Overflows in yield histograms!" << std::endl;
      std::cerr << "This is probably safe, but better increase Parameters::maxYield.\n\n" << std::endl;
    }
  }
}

void th2polyHoneycomb(){
   gStyle->SetCanvasPreferGL(true);
   TH2Poly *hc = new TH2Poly();
   hc->Honeycomb(0,0,.1,25,25);

   TRandom ran;
   for (int i = 0; i<30000; i++) {
      hc->Fill(ran.Gaus(2.,1), ran.Gaus(2.,1));
   }

   hc->Draw("gllego");
}

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

// PValue for finding a local p-value as observed in 'bin' or worse 
void ToyExperiments::printGlobalPValueOfLocalFluctuation(unsigned int bin, unsigned int nExperiments) const {
  std::cout << "Determining global p-value for observed fluctuation in bin " << bin << " from " << nExperiments << " toy experiments ...  " << std::flush;

  // Find the predicted yields that correspond
  // to the local p-value 'localPValue'
  std::vector<unsigned int> limitYields = yields(localPValue(bin,observedYields_.at(bin)));

  TH1* hIsAbovePValue = new TH1D("hIsAbovePValue","",2,0,2);

  const double minCorr = findMinValidRandomNumberForCorrelatedUncertainties();

  for(unsigned int p = 0; p < nExperiments; ++p) {
    bool isAbovePValue = false;
    double rCorr = rand_->Gaus(0.,1.);
    while( rCorr <= minCorr ) {
      rCorr = rand_->Gaus(0.,1.);
    }
    for(unsigned int b = 0; b < Parameters::nBins(); ++b) {
      double prediction = -1.;
      bool negativePrediction = true;
      while( negativePrediction ) {
	double rUncorr = rand_->Gaus(0.,1.);
	double uncorrVar = rUncorr * uncorrelatedUncerts_.at(b);
	double corrVar   = rCorr   * correlatedUncerts_.at(b);
	prediction = meanPredictions_.at(b) + uncorrVar + corrVar;
	if( prediction >= 0. ) {
	  negativePrediction = false;
	}
      }
      double predictedYield = rand_->Poisson(prediction);
      if( predictedYield >= limitYields.at(b) ) {
	isAbovePValue = true;
	break;
      }      
    }
    if( isAbovePValue ) {
      hIsAbovePValue->Fill(1);
    } else {
      hIsAbovePValue->Fill(0);
    }
  }
  std::cout << "ok" << std::endl;

  double lpv      = localPValue(bin,observedYields_.at(bin));
  double gpUncorr = 1. - pow(1.-lpv,Parameters::nBins());
  double gpCorr   = hIsAbovePValue->Integral(2,2)/hIsAbovePValue->Integral(1,2);

  std::cout << "\n\n----- Global p-value for observed fluctuation in bin " << bin << " -----" << std::endl;
  std::cout << "  local p-value                           : " << lpv << " (" << TMath::NormQuantile(1.-lpv) << "sig)" << std::endl;
  std::cout << "  global p-value (without correlations)   : " << gpUncorr  << " (" << TMath::NormQuantile(1.-gpUncorr) << "sig)" << std::endl;
  std::cout << "  global p-value (including correlations) : " << gpCorr  << " (" << TMath::NormQuantile(1.-gpCorr) << "sig)" << std::endl;
}

void JEC_fit_Uncertainty(int N)
{
  TF1 *func[1000];
  TFile *inf = new TFile("L3Graphs_test_Icone5.root");
  TGraphErrors *g = (TGraphErrors*)inf->Get("Correction_vs_CaloPt");
  TGraphErrors *vg[1000];
  int i,k;
  double x[20],y[20],ex[20],ey[20];
  double vx[20],vy[20],vex[20],vey[20];
  for(i=0;i<g->GetN();i++)
    {
      g->GetPoint(i,x[i],y[i]);
      ex[i]=g->GetErrorX(i);
      ey[i]=g->GetErrorY(i); 
    }  
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  for(k=0;k<N;k++)
    {
      for(i=0;i<g->GetN();i++)
        {	
          vx[i] = rnd->Gaus(x[i],ex[i]);
          //vx[i] = x[i];
          vy[i] = rnd->Gaus(y[i],ey[i]);
          vex[i] = ex[i];
          vey[i] = ey[i];
        }
      vg[k] = new TGraphErrors(g->GetN(),vx,vy,vex,vey);
      func[k] = new TF1("func","[0]+[1]/(pow(log10(x),[2])+[3])",1,2000);
      func[k]->SetParameters(1,3,6,5);
      vg[k]->Fit(func[k],"RQ");     	
    }
  
  TCanvas *c = new TCanvas("c","c");
  gPad->SetLogx();
  g->SetMarkerStyle(20);
  g->SetMaximum(3.5);
  g->Draw("AP");
  for(k=0;k<N;k++)
    {
      func[k]->SetLineColor(5);
      func[k]->SetLineWidth(1);
      cout<<func[k]->GetChisquare()<<endl;
      vg[k]->SetMarkerColor(2);
      vg[k]->SetLineColor(2);
      vg[k]->SetMarkerStyle(21);
      //if (func[k]->GetChisquare()<0.1)
        //vg[k]->Draw("sameP");
      func[k]->Draw("same");  
    }  	 
}  

void candlehisto()
{
   TCanvas *c1 = new TCanvas("c1","Candle Presets",800,600);
   c1->Divide(3,2);

   TRandom *rand = new TRandom();
   TH2I *h1 = new TH2I("h1","Sin",18,0,360,100,-1.5,1.5);
   h1->GetXaxis()->SetTitle("Deg");

   float myRand;
   for (int i = 0; i < 360; i+=10) {
      for (int j = 0; j < 100; j++) {
         myRand = rand->Gaus(sin(i*3.14/180),0.2);
         h1->Fill(i,myRand);
      }
   }

   for (int i = 1; i < 7; i++) {
      c1->cd(i);
      char str[16];
      sprintf(str,"CANDLEX%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }

   TCanvas *c2 = new TCanvas("c2","Violin Presets",800,300);
   c2->Divide(2,1);

   for (int i = 1; i < 3; i++) {
      c2->cd(i);
      char str[16];
      sprintf(str,"VIOLINX%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetFillColor(kGray+2);
   }

   TCanvas *c3 = new TCanvas("c3","Playing with candle and violin-options",800,600);
   c3->Divide(3,2);
   char myopt[6][16] = {"1000000","2000000","3000000","1112111","112111","112111"};
   for (int i = 0; i < 6; i++) {
      c3->cd(i+1);
      char str[16];
      sprintf(str, "candlex(%s)",myopt[i]);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetFillColor(kYellow);
      if (i == 4) {
         TH2I * myhist2 = (TH2I*)h1->DrawCopy("candlex(1000000) same");
         myhist2->SetFillColor(kRed);
      }
      if (i == 5) {
         myhist->SetBarWidth(0.2);
         myhist->SetBarOffset(0.25);
         TH2I * myhist2 = (TH2I*)h1->DrawCopy("candlex(2000000) same");
         myhist2->SetFillColor(kRed);
         myhist2->SetBarWidth(0.6);
         myhist2->SetBarOffset(-0.5);
      }
      myhist->SetTitle(str);
   }
}

void tv3Write() {
   //creates the Tree
   TVector3 *v = new TVector3();
   TVector3::Class()->IgnoreTObjectStreamer();
   TFile *f = new TFile("v3.root","recreate");
   TTree *T = new TTree("T","v3 Tree");
   T->Branch("v3","TVector3",&v,32000,1);
   TRandom r;
   for (Int_t i=0;i<10000;i++) {
      v->SetXYZ(r.Gaus(0,1),r.Landau(0,1),r.Gaus(100,10));
      T->Fill();
   }
   T->Write();
   T->Print();
   delete f;
}

void candleplotoption()
{
   TCanvas *c1 = new TCanvas("c1","Candle Presets",1200,800);
   c1->Divide(3,2);

   TRandom *randnum = new TRandom();
   TH2I *h1 = new TH2I("h1","Sin",18,0,360,300,-1.5,1.5);
   h1->GetXaxis()->SetTitle("Deg");
   float myRand;
   for (int i = 0; i < 360; i+=10) {
      for (int j = 0; j < 100; j++) {
         myRand = randnum->Gaus(sin(i*3.14/180),0.2);
         h1->Fill(i,myRand);
      }
   }
   for (int i = 1; i < 7; i++) {
      c1->cd(i);
      char str[16];
      sprintf(str,"candlex%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }

   TCanvas *c2 = new TCanvas("c2","Candle Individual",1200,800);
   c2->Divide(4,4);
   char myopt[16][8] = {"0","1","11","21","31","30","111","311","301","1111","2321","12111","112111","212111","312111"};
   for (int i = 0; i < 15; i++) {
      c2->cd(i+1);
      char str[16];
      sprintf(str, "candlex(%s)",myopt[i]);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }
}

int ConvertEvent::get_electron_numbers(double energy)
{
    if (energy <= 0)
        return 0;
    int n = int(round(random.Gaus(energy / work_function, TMath::Sqrt(fano * energy / work_function))));
    return n > 0 ? n : 0;
}

void Interpolation()
{
   ROOT::R::TRInterface &r=ROOT::R::TRInterface::Instance();
   //Creating points
   TRandom rg;
   std::vector<Double_t> x(10),y(10);
   for(int i=0;i<10;i++)
   {
      x[i]=i;
      y[i]=rg.Gaus();
   }

   r["x"]=x;
   r["y"]=y;


   // do plotting only in non-batch mode
   if (!gROOT->IsBatch() )  {

      r<<"dev.new()";//Required to activate new window for plot
      //Plot parameter. Plotting using two rows and one column
      r<<"par(mfrow = c(2,1))";

      //plotting the points
      r<<"plot(x, y, main = 'approx(.) and approxfun(.)')";

      //The function "approx" returns a list with components x and y
      //containing n coordinates which interpolate the given data points according to the method (and rule) desired.
      r<<"points(approx(x, y), col = 2, pch = '*')";
      r<<"points(approx(x, y, method = 'constant'), col = 4, pch = '*')";
   }
   else {
      r << "print('Interpolated points')";
      r << "print(approx(x,y,n=20))";
   }

   //The function "approxfun" returns a function performing (linear or constant)
   //interpolation of the given data.
   //For a given set of x values, this function will return the corresponding interpolated values.
   r<<"f <- approxfun(x, y)";
   //using approxfun with const method
   r<<"fc <- approxfun(x, y, method = 'const')";

   if (!gROOT->IsBatch() ) {
      r<<"curve(f(x), 0, 11, col = 'green2')";
      r<<"points(x, y)";

      r<<"curve(fc(x), 0, 10, col = 'darkblue', add = TRUE)";
      // different interpolation on left and right side :
      r<<"plot(approxfun(x, y, rule = 2:1), 0, 11,col = 'tomato', add = TRUE, lty = 3, lwd = 2)";
      r<<"dev.off()";//Required to close new window for plot
   }
   else {
      r << "x2=x+0.5";
      r << "print('Result of approxfun with default method')";
      r << "print(paste('x = ',x,'  f(x) = ',f(x2)))";
      r << "print('Result of approxfun with const method')";
      r << "print(paste('x = ',x,'  f(x) = ',fc(x2)))";
   }
}

void MixingResult::randomise () {//randomizer
  static TRandom donram(42); 
  std::cout << "Randomising " << getName() << "; old value "
	    << measurement;
  double delta = donram.Gaus()*error;
  measurement += delta;
  std::cout << " delta " << delta << " new value " << measurement << std::endl; 
}

// Example of a canvas showing two histograms with different scales.
// The second histogram is drawn in a transparent pad
void transpad() {
   TCanvas *c1 = new TCanvas("c1","transparent pad",200,10,700,500);
   TPad *pad1 = new TPad("pad1","",0,0,1,1);
   TPad *pad2 = new TPad("pad2","",0,0,1,1);
   pad2->SetFillStyle(4000); //will be transparent
   pad1->Draw();
   pad1->cd();

   TH1F *h1 = new TH1F("h1","h1",100,-3,3);
   TH1F *h2 = new TH1F("h2","h2",100,-3,3);
   TRandom r;
   for (Int_t i=0;i<100000;i++) {
      Double_t x1 = r.Gaus(-1,0.5);
      Double_t x2 = r.Gaus(1,1.5);
      if (i <1000) h1->Fill(x1);
      h2->Fill(x2);
   }
   h1->Draw();
   pad1->Update(); //this will force the generation of the "stats" box
   TPaveStats *ps1 = (TPaveStats*)h1->GetListOfFunctions()->FindObject("stats");
   ps1->SetX1NDC(0.4); ps1->SetX2NDC(0.6);
   pad1->Modified();
   c1->cd();
   
   //compute the pad range with suitable margins
   Double_t ymin = 0;
   Double_t ymax = 2000;
   Double_t dy = (ymax-ymin)/0.8; //10 per cent margins top and bottom
   Double_t xmin = -3;
   Double_t xmax = 3;
   Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
   pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
   pad2->Draw();
   pad2->cd();
   h2->SetLineColor(kRed);
   h2->Draw("sames");
   pad2->Update();
   TPaveStats *ps2 = (TPaveStats*)h2->GetListOfFunctions()->FindObject("stats");
   ps2->SetX1NDC(0.65); ps2->SetX2NDC(0.85);
   ps2->SetTextColor(kRed);
   
   // draw axis on the right side of the pad
   TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
   axis->SetLabelColor(kRed);
   axis->Draw();
}