stressShapes.cxx
// Author : Mihaela Gheata 12-01-03
#ifndef __CINT__
#include <TRandom3.h>
#include <TROOT.h>
#include <TH1.h>
#include <TMath.h>
#include <TGeoManager.h>
#include <TGeoVolume.h>
#include <TGeoPcon.h>
#include <TGeoMatrix.h>
#include <TBenchmark.h>
#include <TApplication.h>
void stressShapes();
int main(int argc, char **argv)
{
TApplication theApp("App", &argc, argv);
stressShapes();
return 0;
}
#endif
//--- This macro creates a simple geometry based on all shapes known
//--- by TGeo. The first test generates 1 million random points inside
//--- the bounding box of each shape and computes the volume of the
//--- shape as Vbbox*Ninside/Ntotal.
//--- The second test tracks 100K random rays in the geometry, histogramming
//--- the length of all segments passing through each different shape.
//--- It computes mean, RMS and sum of lengths of all segments inside a
//--- given shape and compares with reference values.
//
// This test program is automatically created by $ROOTSYS/test/Makefile.
// To run it in batch, execute stressGeom.
// To run this test with interactive CINT, do
// root > .x stressShapes.cxx++
// or
// root > .x stressShapes.cxx
void sample_volume(Int_t ivol)
{
const Double_t vshape[16] = {40000.0, 36028.3, 39978.70, 48001.3, 28481.2,
8726.2, 42345.4, 9808.2, 12566.8, 64655.6, 37730.4, 23579.7,
25559.5, 18418.3, 49960.2, 47771.1};
gRandom = new TRandom3();
TGeoVolume *vol = (TGeoVolume*)gGeoManager->GetListOfVolumes()->At(ivol);
TGeoShape *shape = vol->GetShape();
Double_t dx = ((TGeoBBox*)shape)->GetDX();
Double_t dy = ((TGeoBBox*)shape)->GetDY();
Double_t dz = ((TGeoBBox*)shape)->GetDZ();
Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
Double_t ratio;
Double_t point[3];
Double_t ngen=1000000;
Double_t iin=0;
Double_t i;
for (i=0; i<ngen; i++) {
point[0] = ox-dx+2*dx*gRandom->Rndm();
point[1] = oy-dy+2*dy*gRandom->Rndm();
point[2] = oz-dz+2*dz*gRandom->Rndm();
if (shape->Contains(point)) iin++;
}
ratio = Double_t(iin)/Double_t(ngen);
Double_t vbox = 8*dx*dy*dz;
Double_t vv = vbox*ratio;
Double_t dvv = TMath::Abs(vv-vshape[ivol-1]);
Double_t sigma = vv/TMath::Sqrt(iin+1);
char result[16];
sprintf(result, "FAILED");
if (dvv<2*sigma) sprintf(result, "OK");
printf("---> testing %-4s ............... %s\n", vol->GetName(), result);
}
void length()
{
const Double_t rms[16] = {6.284, 10.79, 9.545, 14.15, 11.45,
5.871, 7.673, 5.935, 7.61, 5.334, 6.581, 4.954,
7.718, 3.238, 19.09, 14.77};
const Double_t mean[16] = {19.34, 22.53, 18.87, 21.95, 23.29,
16.73, 15.09, 9.516, 12.68, 8.852, 9.518, 7.432,
8.881, 6.489, 28.29, 26.05};
TObjArray *vlist = gGeoManager->GetListOfVolumes();
TGeoVolume *volume;
Int_t nvolumes = vlist->GetEntriesFast();
Double_t len[17];
TList *hlist = new TList();
TH1F *hist;
Int_t i=0;
memset(len, 0, nvolumes*sizeof(Double_t));
for (i=0; i<nvolumes; i++) {
volume = (TGeoVolume*)(vlist->At(i));
hist = new TH1F(volume->GetName(), "lengths inside", 100, 0, 100);
hist->SetBit(TH1::kCanRebin);
hlist->Add(hist);
}
Int_t nrays = 100000;
Double_t dir[3];
TGeoNode *startnode, *endnode;
Int_t istep=0, icrt;
Int_t itot=0;
Int_t n10=nrays/10;
Double_t theta,phi, step;
while (itot<nrays) {
itot++;
if (n10) {
if ((itot%n10) == 0) printf(" %i percent\n", Int_t(100*itot/nrays));
}
phi = 2*TMath::Pi()*gRandom->Rndm();
theta= TMath::ACos(1.-2.*gRandom->Rndm());
dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
dir[2]=TMath::Cos(theta);
gGeoManager->InitTrack(0,0,0, dir[0], dir[1], dir[2]);
startnode = gGeoManager->GetCurrentNode();
if (gGeoManager->IsOutside()) startnode=0;
icrt = 0;
if (startnode) icrt =vlist->IndexOf(startnode->GetVolume());
// find where we end-up
gGeoManager->FindNextBoundary();
step = gGeoManager->GetStep();
endnode = gGeoManager->Step();
while (step<1E10) {
while (!gGeoManager->IsEntering()) {
istep++;
if (istep>10000) break;
gGeoManager->SetStep(1E-3);
endnode = gGeoManager->Step();
step += 1E-3;
}
if (istep>10000) break;
len[icrt] += step;
hist = (TH1F*)(hlist->At(icrt));
hist->Fill(step);
// now see if we can make an other step
if (endnode==0 && step>1E10) break;
istep = 0;
// generate an extra step to cross boundary
startnode = endnode;
icrt = 0;
if (startnode) icrt =vlist->IndexOf(startnode->GetVolume());
gGeoManager->FindNextBoundary();
step = gGeoManager->GetStep();
endnode = gGeoManager->Step();
}
}
// draw all segments
Double_t drms, dmean;
for (i=1; i<nvolumes; i++) {
volume = (TGeoVolume*)(vlist->At(i));
hist = (TH1F*)(hlist->At(i));
char result[16];
drms = TMath::Abs(rms[i-1]-hist->GetRMS());
dmean = TMath::Abs(mean[i-1]-hist->GetMean());
sprintf(result, "FAILED");
if (dmean<0.01) {
if (drms<0.01) sprintf(result,"OK");
}
printf(" %-4s : mean_len=%7.4g RMS=%7.4g total_len=%11.4g ... %s\n",
volume->GetName(), hist->GetMean(), hist->GetRMS(), len[i],result);
}
hlist->Delete();
delete hlist;
}
void stressShapes()
{
// New geometry test suite. Creates a geometry containing all shape
// types. Loop over all volumes and compute the following :
// - generate 1 million random points and count how many are inside
// each shape -> compute volume of each shape
// - generate 10000 random directions and propagate from the center
// of each volume -> compute total step length to exit current shape
#ifdef __CINT__
gSystem->Load("libGeom");
#endif
gBenchmark = new TBenchmark();
gBenchmark->Start("stressShapes");
TGeoManager *geom = new TGeoManager("stressShapes", "arbitrary shapes");
TGeoMaterial *mat;
TGeoMixture *mix;
//---> create some materials
mat = new TGeoMaterial("Vacuum",0,0,0);
mat->SetUniqueID(0);
mat = new TGeoMaterial("Be", 9.01,4,1.848);
mat->SetUniqueID(1);
mat = new TGeoMaterial("Al", 26.98,13,2.7);
mat->SetUniqueID(2);
mat = new TGeoMaterial("Fe", 55.85,26,7.87);
mat->SetUniqueID(3);
mat = new TGeoMaterial("Cu", 63.55,29,8.96);
mat->SetUniqueID(4);
mat = new TGeoMaterial("C",12.01,6,2.265);
mat->SetUniqueID(5);
mat = new TGeoMaterial("Pb",207.19,82,11.35);
mat->SetUniqueID(6);
mat = new TGeoMaterial("Si",28.09,14,2.33);
mat->SetUniqueID(7);
mix = new TGeoMixture("scint",2, 1.03200 );
mix->DefineElement(0,1.008,1,0.7749078E-01);
mix->DefineElement(1,12,6,0.9225092);
mix->SetUniqueID(8);
mat = new TGeoMaterial("Li",6.94,3,0.534);
mat->SetUniqueID(9);
mat = new TGeoMaterial("N",14.01,7,0.808);
mat->SetUniqueID(10);
mat = new TGeoMaterial("Ne",20.18,10,1.207);
mat->SetUniqueID(11);
mat = new TGeoMaterial("Ts",183.85,74,19.3);
mat->SetUniqueID(12);
mat = new TGeoMaterial("CFRP",12.01,6,2.3);
mat->SetUniqueID(13);
mix = new TGeoMixture("H2O",2, 1.00000 );
mix->DefineElement(0,1.01,1,0.1120977);
mix->DefineElement(1,16,8,0.8879023);
mix->SetUniqueID(14);
mix = new TGeoMixture("Ethane_Gas",2, 0.135600E-02);
mix->DefineElement(0,12.01,6,0.7985373);
mix->DefineElement(1,1.01,1,0.2014628);
mix->SetUniqueID(15);
mat = new TGeoMaterial("RHONEYCM",26.98,13,0.601);
mat->SetUniqueID(16);
//---> create mediums
TGeoMedium *med0 = new TGeoMedium("Vacuum",0,0,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med1 = new TGeoMedium("Be",1,1,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med2 = new TGeoMedium("Al",2,2,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med3 = new TGeoMedium("Fe",3,3,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med4 = new TGeoMedium("Cu",4,4,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med5 = new TGeoMedium("C",5,5,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med6 = new TGeoMedium("Pb",6,6,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med7 = new TGeoMedium("Si",7,7,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med8 = new TGeoMedium("scint",8,8,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med9 = new TGeoMedium("Li",9,9,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med10 = new TGeoMedium("N",10,10,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med11 = new TGeoMedium("Ne",11,11,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med12 = new TGeoMedium("Ts",12,12,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med13 = new TGeoMedium("CFRP",13,13,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med14 = new TGeoMedium("H2O",14,14,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med15 = new TGeoMedium("Ethane gas",15,15,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
TGeoMedium *med16 = new TGeoMedium("RHONEYCM",16,16,0,0,0,20,0.1000000E+11,0.212,0.1000000E-02,1.150551);
//---> create volumes
TGeoVolume *TOP = geom->MakeBox("TOP", med0, 200, 200, 200);
TGeoVolume *BOX = geom->MakeBox("BOX", med1, 10,20,25);
TGeoVolume *TRD1 = geom->MakeTrd1("TRD1", med2, 20,10,20,15);
TGeoVolume *TRD2 = geom->MakeTrd2("TRD2", med3, 20,5,10,25,25);
TGeoVolume *PARA = geom->MakePara("PARA", med4, 10, 20, 30, 15, 15, 120);
Double_t v[16];
v[0]=-22; v[1]=-18; v[2]=-18; v[3]=22; v[4]=22; v[5]=18; v[6]=18; v[7]=-22;
v[8]=-12; v[9]=-8; v[10]=-8; v[11]=12; v[12]=12; v[13]=8; v[14]=8; v[15]=-12;
TGeoVolume *ARB8 = geom->MakeArb8("ARB8", med5, 15, v);
TGeoVolume *SPHE = geom->MakeSphere("SPHE", med6, 5, 15, 45, 180, 0, 270);
TGeoVolume *TUBE = geom->MakeTube("TUBE", med7, 20, 25, 30);
TGeoVolume *TUBS = geom->MakeTubs("TUBS", med8, 10,15,20, 45, 270);
TGeoVolume *ELTU = geom->MakeEltu("ELTU", med9, 10,20,10);
TGeoVolume *CTUB = geom->MakeCtub("CTUB", med10, 25, 30, 50, 0,270, 0, 0.5, -0.5*TMath::Sqrt(3.),0.5,0,0.5*TMath::Sqrt(3.));
TGeoVolume *CONE = geom->MakeCone("CONE", med11, 30, 25, 30, 10, 15);
TGeoVolume *CONS = geom->MakeCons("CONS", med12, 30, 25, 30, 10, 15, -45, 180);
TGeoVolume *PCON = geom->MakePcon("PCON", med13, 0, 360, 3);
TGeoPcon *pcon = (TGeoPcon*)(PCON->GetShape());
pcon->DefineSection(0,-25, 10, 15);
pcon->DefineSection(1,0, 10, 15);
pcon->DefineSection(2,25, 25, 30);
TGeoVolume *PGON = geom->MakePgon("PGON", med14, 0, 270, 4, 3);
pcon = (TGeoPcon*)(PGON->GetShape());
pcon->DefineSection(0,-25, 10, 15);
pcon->DefineSection(1,0, 10, 15);
pcon->DefineSection(2,25, 15, 20);
TGeoVolume *TRAP = geom->MakeTrap("TRAP", med15, 25, 15, 30, 20, 15,10,15, 20,15,10,15);
TGeoVolume *GTRA = geom->MakeGtra("GTRA", med16, 25, 15, 30, 30, 20, 15,10,15, 20,15,10,15);
//---> create nodes
geom->SetTopVolume(TOP);
TOP->AddNode(BOX, 1);
TOP->AddNode(TRD1, 2, new TGeoTranslation(100, 0, 0));
TOP->AddNode(TRD2, 3, new TGeoTranslation(-100, 0, 0));
TOP->AddNode(PARA, 4, new TGeoTranslation(0, 100, 0));
TOP->AddNode(ARB8, 5, new TGeoTranslation(0, -100, 0));
TOP->AddNode(SPHE, 6, new TGeoTranslation(0, 0, 100));
TOP->AddNode(TUBE, 7, new TGeoTranslation(0, 0, -100));
TOP->AddNode(TUBS, 8, new TGeoTranslation(100, 0, 100));
TOP->AddNode(ELTU, 9, new TGeoTranslation(100, 0, -100));
TOP->AddNode(CTUB, 10, new TGeoTranslation(-100, 0, 100));
TOP->AddNode(CONE, 11, new TGeoTranslation(-100, 0, -100));
TOP->AddNode(CONS, 12, new TGeoTranslation(0, 100, 100));
TOP->AddNode(PCON, 13, new TGeoTranslation(0, -100, 100));
TOP->AddNode(PGON, 14, new TGeoTranslation(100, 100, 100));
TOP->AddNode(TRAP, 15, new TGeoTranslation(-100, -100, -100));
TOP->AddNode(GTRA, 16, new TGeoTranslation(100, 100, 0));
//---> close geometry
geom->CloseGeometry("d");
geom->DefaultColors();
TIter next(gGeoManager->GetListOfVolumes());
TGeoVolume *vol = (TGeoVolume*)next();
Int_t ivol=1;
printf("=== testing shapes ...\n");
while ((vol=(TGeoVolume*)next())) {
sample_volume(ivol);
ivol++;
}
printf("=== testing global tracking ...\n");
length();
// print ROOTMARKs
printf("\n");
gBenchmark->Stop("stressShapes");
gBenchmark->Print("stressShapes");
Float_t ct = gBenchmark->GetCpuTime("stressShapes");
Float_t cp_brun = 5.55;
Float_t rootmarks = 860*cp_brun/ct;
printf("*******************************************************************\n");
printf("* ROOTMARKS =%6.1f * Root%-8s %d/%d CP = %7.2fs\n",rootmarks,gROOT->GetVersion(),gROOT->GetVersionDate(),gROOT->GetVersionTime(),ct);
printf("*******************************************************************\n");
}
