ROOT_528-00b_version: tutorials/gl/nucleus.C Source File

00001 // Model of a nucleus built from TGeo classes.
00002 // 
00003 // Author: Otto Schaile
00004 
00005 void nucleus(Int_t nProtons  = 40,Int_t  nNeutrons = 60) 
00006 {
00007    Double_t NeutronRadius = 60,
00008             ProtonRadius = 60,
00009             NucleusRadius,
00010             distance = 60;
00011    Double_t vol = nProtons + nNeutrons;
00012    vol = 3 * vol / (4 * TMath::Pi());
00013 
00014    NucleusRadius = distance * TMath::Power(vol, 1./3.);
00015 //   cout << "NucleusRadius: " << NucleusRadius << endl;
00016 
00017    TGeoManager * geom = new TGeoManager("nucleus", "Model of a nucleus");
00018    geom->SetNsegments(40);
00019    TGeoMaterial *matEmptySpace = new TGeoMaterial("EmptySpace", 0, 0, 0);
00020    TGeoMaterial *matProton     = new TGeoMaterial("Proton"    , .938, 1., 10000.);
00021    TGeoMaterial *matNeutron    = new TGeoMaterial("Neutron"   , .935, 0., 10000.);
00022 
00023    TGeoMedium *EmptySpace = new TGeoMedium("Empty", 1, matEmptySpace);
00024    TGeoMedium *Proton     = new TGeoMedium("Proton", 1, matProton);
00025    TGeoMedium *Neutron    = new TGeoMedium("Neutron",1, matNeutron);
00026 
00027 //  the space where the nucleus lives (top container volume)
00028 
00029    Double_t worldx = 200.;
00030    Double_t worldy = 200.;
00031    Double_t worldz = 200.;
00032  
00033    TGeoVolume *top = geom->MakeBox("WORLD", EmptySpace, worldx, worldy, worldz); 
00034    geom->SetTopVolume(top);
00035 
00036    TGeoVolume * proton  = geom->MakeSphere("proton",  Proton,  0., ProtonRadius); 
00037    TGeoVolume * neutron = geom->MakeSphere("neutron", Neutron, 0., NeutronRadius); 
00038    proton->SetLineColor(kRed);
00039    neutron->SetLineColor(kBlue);
00040 
00041    Double_t x, y, z, dummy;
00042    Int_t i = 0; 
00043    while ( i<  nProtons) {
00044       gRandom->Rannor(x, y);
00045       gRandom->Rannor(z,dummy);
00046       if ( TMath::Sqrt(x*x + y*y + z*z) < 1) {
00047          x = (2 * x - 1) * NucleusRadius;
00048          y = (2 * y - 1) * NucleusRadius;
00049          z = (2 * z - 1) * NucleusRadius;
00050          top->AddNode(proton, i, new TGeoTranslation(x, y, z));
00051          i++;
00052       }
00053    }
00054    i = 0; 
00055    while ( i <  nNeutrons) {
00056       gRandom->Rannor(x, y);
00057       gRandom->Rannor(z,dummy);
00058       if ( TMath::Sqrt(x*x + y*y + z*z) < 1) {
00059          x = (2 * x - 1) * NucleusRadius;
00060          y = (2 * y - 1) * NucleusRadius;
00061          z = (2 * z - 1) * NucleusRadius;
00062          top->AddNode(neutron, i + nProtons, new TGeoTranslation(x, y, z));
00063          i++;
00064       }
00065    }
00066    geom->CloseGeometry();
00067    geom->SetVisLevel(4);
00068    top->Draw("ogl");
00069 }