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

void GenerateNonStd(unsigned int Npart)
{
  frodo *fr = frodo::instance();

  static TRandom randy(0);

  fr->HBDParticles.clear();
  
  double s = randy.Rndm();
  //cout<<"s= "<<s<<endl;
  double s1 = 0;
  if(s < 0.5) s1 = 1;
  else s1 = -1;
  
  double a = randy.Rndm();
  //cout<<"a= "<<a <<endl;
  double Cx = 2*fr->Ring[0].C()*a - fr->Ring[0].C();
  double Cy = s1*sqrt(fr->Ring[0].C()*fr->Ring[0].C()-Cx*Cx);
  fr->CX.push_back(Coordinates(Cx));
  fr->CY.push_back(Coordinates(Cy));
  //cout<<"ring center is: ( "<< Cx <<" , " << Cy << " ) " <<endl;

  for (unsigned int i=0; i<Npart; i++)
    {
      //cout<<"Npart= " <<i<<endl;
      double S = randy.Rndm();
      //cout<<"S= "<<S<<endl;
      double s2 = 0;
      if(S < 0.5) s2 = 1;
      else s2 = -1;

      double b = randy.Rndm();
      //cout<<"b= "<<b<<endl;
      
      double x = Cx + 2*fr->Ring[0].R()*b - fr->Ring[0].R();      
      double y = Cy + s2*sqrt(fr->Ring[0].R()*fr->Ring[0].R() - (x-Cx)*(x-Cx));
      
      if(fabs(x) <= 1000 || fabs(y) <= 1125)
	{
	  i--;
	  //cout<<"reject!"<<endl;
	  continue;
	}
      double q = -1;
      while (q<=0) q=randy.Exp(80000);
      //cout<<"q= " <<q<<endl;
      
      fr->HBDParticles.push_back( AParticle(x,y,q) );
    }
  //cout<<"# of particles: " <<fr->HBDParticles.size()<<endl;
  
  
}

void drawing_pion_decay(){

	double Pi_lifetime = 0.26033e-6; //s
	double c = 3e8;//m/s
	TCanvas *c1 = new TCanvas("test", "test");
	TRandom *r = new TRandom();
	

	TView *view = TView::CreateView(1, 0, 0);
	view->ShowAxis();
	view->SetRange(-5, -5, 0, 5, 5, 10);
	
	

	for (int i = 0; i < 100; ++i)
	{
		double px = 0.3;
		double py = 0.3;
		double pz = -1; // GeV
		double energy = sqrt(px*px+py*py+pz*pz+M_pion*M_pion);
		
		//life time
		double t = r->Exp(Pi_lifetime);
		//decay length
		double gamma = energy/M_pion;
		double length = c*t*gamma/1e3;
		// decay position
		double vx = 0;
		double vy = 0;
		double vz = 10 - length;

		TLorentzVector pion(px, py, pz, energy);

		double decay_particle_mass[2] = {M_muon, M_neu};
		TGenPhaseSpace event;
		event.SetDecay(pion, 2, decay_particle_mass);
		event.Generate();

		TLorentzVector Muon = *(event.GetDecay(0));
		TLorentzVector Neu = *(event.GetDecay(1));
		
		plot_particle(Muon, vx, vy, vz, 2);
		plot_particle(Neu, vx, vy, vz, 4);


	}


}

void principal(Int_t n=10, Int_t m=10000) 
{
  // 
  // Principal Components Analysis (PCA) example
  // 
  // Example of using TPrincipal as a stand alone class. 
  // 
  // We create n-dimensional data points, where c = trunc(n / 5) + 1
  // are  correlated with the rest n - c randomly distributed variables. 
  //
  // Here's the plot of the eigenvalues Begin_Html
  // <IMG SRC="gif/principal_eigen.gif">
  // End_Html
  //Authors: Rene Brun, Christian Holm Christensen
    
  Int_t c = n / 5 + 1;

  cout << "*************************************************" << endl; 
  cout << "*         Principal Component Analysis          *" << endl;
  cout << "*                                               *" << endl;
  cout << "*  Number of variables:           " << setw(4) << n 
       << "          *" << endl;
  cout << "*  Number of data points:         " << setw(8) << m
       << "      *" << endl;
  cout << "*  Number of dependent variables: " << setw(4) << c
       << "          *" << endl;
  cout << "*                                               *" << endl;
  cout << "*************************************************" << endl; 
  
      
  // Initilase the TPrincipal object. Use the empty string for the
  // final argument, if you don't wan't the covariance
  // matrix. Normalising the covariance matrix is a good idea if your
  // variables have different orders of magnitude. 
  TPrincipal* principal = new TPrincipal(n,"ND");
  
  // Use a pseudo-random number generator
  TRandom* random = new TRandom;
  
  // Make the m data-points
  // Make a variable to hold our data
  // Allocate memory for the data point
  Double_t* data = new Double_t[n];
  for (Int_t i = 0; i < m; i++) {

    // First we create the un-correlated, random variables, according
    // to one of three distributions 
    for (Int_t j = 0; j < n - c; j++) {
       if (j % 3 == 0)
          data[j] = random->Gaus(5,1);
       else if (j % 3 == 1)
          data[j] = random->Poisson(8);
       else
          data[j] = random->Exp(2);
    }

    // Then we create the correlated variables
    for (Int_t j = 0 ; j < c; j++) {
       data[n - c + j] = 0;
       for (Int_t k = 0; k < n - c - j; k++)
          data[n - c + j] += data[k];
    }
    
    // Finally we're ready to add this datapoint to the PCA
    principal->AddRow(data);
  }
    
  // We delete the data after use, since TPrincipal got it by now. 
  delete [] data;
  
  // Do the actual analysis
  principal->MakePrincipals();
  
  // Print out the result on
  principal->Print();

  // Test the PCA 
  principal->Test();

  // Make some histograms of the orginal, principal, residue, etc data 
  principal->MakeHistograms();
  
  // Make two functions to map between feature and pattern space 
  principal->MakeCode();

  // Start a browser, so that we may browse the histograms generated
  // above 
  TBrowser* b = new TBrowser("principalBrowser", principal);
  
}