use strict;
use warnings;
use SOOT ':all';
use Math::Trig;
$gSystem->Load("libGeom");
$gSystem->Load("libGeomBuilder");
$gSystem->Load("libGeomPainter");
SOOT->UpdateClasses();
# use TGeo classes to draw a model of a nucleus
#
# Author: Otto Schaile
my $nProtons = shift || 40;
my $nNeutrons = shift || 60;
my $NeutronRadius = 60;
my $ProtonRadius = 60;
my $NucleusRadius;
my $distance = 60;
my $vol = $nProtons + $nNeutrons;
$vol = 3 * $vol / (4 * pi);
$NucleusRadius = $distance * $vol**(1./3.);
my $geom = TGeoManager->new("nucleus", "Model of a nucleus");
$geom->SetNsegments(40);
my $matEmptySpace = TGeoMaterial->new("EmptySpace", 0, 0, 0);
my $matProton = TGeoMaterial->new("Proton" , .938, 1., 10000.);
my $matNeutron = TGeoMaterial->new("Neutron" , .935, 0., 10000.);
my $EmptySpace = TGeoMedium->new("Empty", 1, $matEmptySpace);
my $Proton = TGeoMedium->new("Proton", 1, $matProton);
my $Neutron = TGeoMedium->new("Neutron",1, $matNeutron);
# the space where the nucleus lives (top container volume)
my $worldx = 200.;
my $worldy = 200.;
my $worldz = 200.;
my $top = $geom->MakeBox("WORLD", $EmptySpace, $worldx, $worldy, $worldz);
$geom->SetTopVolume($top);
my $proton = $geom->MakeSphere("proton", $Proton, 0., $ProtonRadius);
my $neutron = $geom->MakeSphere("neutron", $Neutron, 0., $NeutronRadius);
$proton->SetLineColor(kRed);
$neutron->SetLineColor(kBlue);
my ($x, $y, $z);
my $i = 0;
while ($i < $nProtons) {
$x = $gRandom->Gaus(0.0, 1.0);
$y = $gRandom->Gaus(0.0, 1.0);
$z = $gRandom->Gaus(0.0, 1.0);
printf "%f %f %f\n", $x, $y, $z;
if (sqrt($x**2 + $y**2 + $z**2) < 1) {
$x = (2 * $x - 1) * $NucleusRadius;
$y = (2 * $y - 1) * $NucleusRadius;
$z = (2 * $z - 1) * $NucleusRadius;
my $trans = TGeoTranslation->new($x*1.0, $y*1.0, $z*1.0)->keep;
$top->AddNode($proton, $i, $trans);
$i++;
}
}
$i = 0;
while ($i < $nNeutrons) {
$x = $gRandom->Gaus(0.0, 1.0);
$y = $gRandom->Gaus(0.0, 1.0);
$z = $gRandom->Gaus(0.0, 1.0);
if (sqrt($x**2 + $y**2 + $z**2) < 1) {
$x = (2 * $x - 1) * $NucleusRadius;
$y = (2 * $y - 1) * $NucleusRadius;
$z = (2 * $z - 1) * $NucleusRadius;
my $trans = TGeoTranslation->new($x*1.0, $y*1.0, $z*1.0)->keep;
$top->AddNode($neutron, $i + $nProtons, $trans);
$i++;
}
}
$geom->CloseGeometry();
$geom->SetVisLevel(4);
$top->Draw("ogl");
$gApplication->Run;