float pi = 3.1415926;
float pAlpha = 18;  //MeV/c
float mp = 931.5; //MeV/c2
float mn = 934.5; //MeV/c2
  
float mBe    = 4*mp + 4*mn + 4.942;
float mAlpha = 2*mp + 2*mn + 2.425;

TVector3 sorteiaAlpha1()
{
  // sorteia direcao de decaimento de uma das alphas
  // o decaimento é isotrópico
  float thetaAlpha1 = gRandom->Uniform(-pi,pi);
  float phiAlpha1   = gRandom->Uniform(0,2*pi);

  // cria os vetores de decaimento das alphas no 
  // centro de massa do berílio
  TVector3 pAlpha1CM;
  pAlpha1CM.SetMagThetaPhi(pAlpha,thetaAlpha1,phiAlpha1);

  return pAlpha1CM;
}

void simula(float VMAX, int N)
{
  gStyle->SetPalette(1,0);
  
  // A distr. de probabilidades eh dada por uma
  // distr de Maxwell-Boltzmann
  TF1 *prob = new TF1("prob","(x^2)*exp(-[0]*x^2)",0,0.2);
  prob->SetParameter(0,1.0/(VMAX*VMAX));
  
  TH1F *h1 = new TH1F("Ealpha","Energia das particulas alpha (MeV)", 1000, 0, 20);
  TH1F *h2 = new TH1F("Aalpha","Angulo das particulas alpha (MeV)", 1000, 0, 2*pi);
  TH1F *h3 = new TH1F("EBe","Energia dos berilios (MeV)", 1000, 0, 20);

  TH2F *h4 = new TH2F("EA","Energia em funcao do angulo", 1000, 0, 2*pi, 1000, 0, 20);

  for(int i = 0; i<N; i++)
  {    
    // cria o vetor de momento do berilio no espaco  
    TVector3 pBe;
    float v = prob->GetRandom(); // em unidades de C
    float p = mBe*v;             // em MeV/c
    pBe.SetMagThetaPhi(p,0,0);
    float EBe8 = pBe.Mag2()/(2.0*mBe); // em MeV

    //sorteia particula alpha1 e calcula alpha2
    TVector3 pAlpha1CM = sorteiaAlpha1();
    TVector3 pAlpha2CM = -pAlpha1CM;
    
    // converte para o referencial do laboratorio
    TVector3 pAlpha1 = pAlpha1CM + pBe;
    TVector3 pAlpha2 = pAlpha2CM + pBe;
    
    // calcula energia cinetica das alphas
    // E = p²/2m  
    float EAlpha1 = pAlpha1.Mag2()/(2.0*mAlpha);
    float EAlpha2 = pAlpha2.Mag2()/(2.0*mAlpha);

    // calcula os angulos relativos entre as
    // alphas e o berílio (angle1 e angle2)
    float angle1 = pAlpha1.Angle(pBe);
    float angle2 = pAlpha2.Angle(pBe);
    
    h1->Fill(EAlpha1);
    h1->Fill(EAlpha2);
    
    h2->Fill(angle1);
    h2->Fill(angle2);    

    h3->Fill(EBe8);    

    h4->Fill(angle1,EAlpha1);
    h4->Fill(angle2,EAlpha2);
  }
  
  TCanvas *c1 = new TCanvas();
  h3->Draw();
  TCanvas *c2 = new TCanvas();
  h1->Draw();
  TCanvas *c3 = new TCanvas();
  h2->Draw();
  TCanvas *c4 = new TCanvas();
  h4->Draw("col");
  
}