macros/checkgeom.C

   1 #include "TGeoManager.h"
   2 #include "TGeoVolume.h"
   3 #include "TGeoNode.h"
   4 #include "TGeoBBox.h"
   5 #include "TRandom3.h"
   6 #include "TMath.h"
   7 #include "TH1.h"
   8 #include "TH2.h"
   9 #include "TCanvas.h"
  10 #include "TStopwatch.h"
  11 #include "TFile.h"
  12
  13 TStopwatch *timer = 0;
  14
  15 Int_t propagate_in_geom(Double_t *, Double_t *);
  16 void score(TGeoVolume *, Int_t, Double_t);
  17 Double_t timing_per_volume(TGeoVolume *);
  18 Double_t *val1 = 0;
  19 Double_t *val2 = 0;
  20 Bool_t   *flags = 0;
  21
  22 void checkgeom()
  23 {
  24    TGeoManager::Import("geometry.root");
  25    Int_t nvol = gGeoManager->GetListOfVolumes()->GetEntries();
  26    Int_t nuid = gGeoManager->GetListOfUVolumes()->GetEntries();
  27    timer = new TStopwatch();
  28    Int_t i;
  29    Double_t value;
  30    flags = new Bool_t[nuid];
  31    val1 = new Double_t[nuid];
  32    val2 = new Double_t[nuid];
  33    memset(flags, 0, nuid*sizeof(Bool_t));
  34    memset(val1, 0, nuid*sizeof(Double_t));
  35    memset(val2, 0, nuid*sizeof(Double_t));
  36    TGeoVolume *vol;
  37
  38 // STAGE 1: Overlap checking by sampling per volume
  39
  40    printf("====================================================================\n");
  41    printf("STAGE 1: Overlap checking by sampling per volume\n");
  42    printf("====================================================================\n");
  43    for (i=0; i<nvol; i++) {
  44       vol = (TGeoVolume*)gGeoManager->GetListOfVolumes()->At(i);
  45       Int_t uid = vol->GetNumber();
  46       if (flags[uid]) continue;
  47       flags[uid] = kTRUE;
  48       if (!vol->GetNdaughters()) continue;
  49       vol->CheckOverlaps(0.01, "s");
  50    }
  51
  52 // STAGE 2: Global overlap checking
  53    printf("====================================================================\n");
  54    printf("STAGE 2: Global overlap checking\n");
  55    printf("====================================================================\n");
  56    gGeoManager->CheckOverlaps(0.01);
  57
  58 // STAGE 3: How many crossings per volume in a realistic case ?
  59    // Ignore volumes with no daughters
  60
  61    // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
  62    Int_t ntracks = 1000000;
  63    printf("====================================================================\n");
  64    printf("STAGE 3: Propagating %i tracks starting from vertex\n and conting number of boundary crossings...\n", ntracks);
  65    printf("====================================================================\n");
  66    Int_t nbound = 0;
  67    Double_t theta, phi;
  68    Double_t point[3], dir[3];
  69    new TRandom3();
  70
  71    timer->Start();
  72    for (i=0; i<ntracks; i++) {
  73       phi = 2.*TMath::Pi()*gRandom->Rndm();
  74       theta= TMath::ACos(1.-2.*gRandom->Rndm());
  75       point[0] = point[1] = point[2] = 0;
  76       dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
  77       dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
  78       dir[2]=TMath::Cos(theta);
  79       if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
  80       nbound += propagate_in_geom(point,dir);
  81    }
  82    Double_t time1 = timer->CpuTime() *1.E6;
  83    Double_t time2 = time1/ntracks;
  84    Double_t time3 = time1/nbound;
  85    printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
  86    printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);
  87
  88 // STAGE 4: How much time per volume:
  89
  90    printf("====================================================================\n");
  91    printf("STAGE 4: How much navigation time per volume per next+safety call\n");
  92    printf("====================================================================\n");
  93    TGeoIterator next(gGeoManager->GetTopVolume());
  94    TGeoNode*current;
  95    TString path;
  96    vol = gGeoManager->GetTopVolume();
  97    memset(flags, 0, nuid*sizeof(Bool_t));
  98    score(vol, 1, timing_per_volume(vol));
  99    while ((current=next())) {
 100       vol = current->GetVolume();
 101       Int_t uid = vol->GetNumber();
 102       if (flags[uid]) continue;
 103       flags[uid] = kTRUE;
 104       next.GetPath(path);
 105       gGeoManager->cd(path.Data());
 106       score(vol,1,timing_per_volume(vol));
 107    }
 108
 109    // Draw some histos
 110    Double_t time_tot_pertrack = 0.;
 111    TCanvas *c1 = new TCanvas("c2","ncrossings",10,10,900,500);
 112    c1->SetGrid();
 113    c1->SetTopMargin(0.15);
 114    TFile *f = new TFile("histos.root", "RECREATE");
 115    TH1F *h = new TH1F("h","number of boundary crossings per volume",3,0,3);
 116    h->SetStats(0);
 117    h->SetFillColor(38);
 118    h->SetBit(TH1::kCanRebin);
 119
 120    memset(flags, 0, nuid*sizeof(Bool_t));
 121    for (i=0; i<nuid; i++) {
 122       vol = gGeoManager->GetVolume(i);
 123       if (!vol->GetNdaughters()) continue;
 124       time_tot_pertrack += val1[i]*val2[i];
 125       h->Fill(vol->GetName(), (Int_t)val1[i]);
 126    }
 127    time_tot_pertrack /= ntracks;
 128    h->LabelsDeflate();
 129    h->LabelsOption(">","X");
 130    h->Draw();
 131
 132
 133    TCanvas *c2 = new TCanvas("c3","time spent per volume in navigation",10,10,900,500);
 134    c2->SetGrid();
 135    c2->SetTopMargin(0.15);
 136    TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0,100,100,0,15);
 137    h2->SetStats(0);
 138    h2->SetMarkerStyle(2);
 139    TH1F *h1 = new TH1F("h1","percent of time spent per volume",3,0,3);
 140    h1->SetStats(0);
 141    h1->SetFillColor(38);
 142    h1->SetBit(TH1::kCanRebin);
 143    for (i=0; i<nuid; i++) {
 144       vol = gGeoManager->GetVolume(i);
 145       if (!vol->GetNdaughters()) continue;
 146       value = val1[i]*val2[i]/ntracks/time_tot_pertrack;
 147       h1->Fill(vol->GetName(), value);
 148       h2->Fill(vol->GetNdaughters(), val2[i]);
 149    }
 150    h1->LabelsDeflate();
 151    h1->LabelsOption(">","X");
 152    h1->Draw();
 153    TCanvas *c3 = new TCanvas("c4","timing vs. ndaughters",10,10,900,500);
 154    c3->SetGrid();
 155    c3->SetTopMargin(0.15);
 156    h2->Draw();
 157    f->Write();
 158    delete [] flags;
 159    delete [] val1;
 160    delete [] val2;
 161 }
 162
 163 Int_t propagate_in_geom(Double_t *start, Double_t *dir)
 164 {
 165 // Propagate from START along DIR from boundary to boundary until exiting
 166 // geometry. Fill array of hits.
 167    gGeoManager->InitTrack(start, dir);
 168    TGeoNode *current = 0;
 169    Int_t nzero = 0;
 170    Int_t nhits = 0;
 171    while (!gGeoManager->IsOutside()) {
 172       current = gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
 173       if (!current || gGeoManager->IsOutside()) return nhits;
 174       Double_t step = gGeoManager->GetStep();
 175       if (step<2.*TGeoShape::Tolerance()) {
 176          nzero++;
 177          continue;
 178       }
 179       else nzero = 0;
 180       if (nzero>3) {
 181       // Problems in trying to cross a boundary
 182          printf("Error in trying to cross boundary of %s\n", current->GetName());
 183          return nhits;
 184       }
 185       // Generate the hit
 186       nhits++;
 187       TGeoVolume *vol = current->GetVolume();
 188       score(vol,0,1.);
 189       Int_t iup = 1;
 190       TGeoNode *mother = gGeoManager->GetMother(iup++);
 191       while (mother && mother->GetVolume()->IsAssembly()) {
 192          score(mother->GetVolume(), 0, 1.);
 193          mother = gGeoManager->GetMother(iup++);
 194       }
 195    }
 196    return nhits;
 197 }
 198
 199 void score(TGeoVolume *vol, Int_t ifield, Double_t value)
 200 {
 201 // Score something for VOL
 202    Int_t uid = vol->GetNumber();
 203    switch (ifield) {
 204       case 0:
 205          val1[uid] += value;
 206          break;
 207       case 1:
 208          val2[uid] += value;
 209    }
 210 }
 211
 212 Double_t timing_per_volume(TGeoVolume *vol)
 213 {
 214 // Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
 215 // in the current path.
 216    timer->Reset();
 217    const TGeoShape *shape = vol->GetShape();
 218    TGeoBBox *box = (TGeoBBox *)shape;
 219    Double_t dx = box->GetDX();
 220    Double_t dy = box->GetDY();
 221    Double_t dz = box->GetDZ();
 222    Double_t ox = (box->GetOrigin())[0];
 223    Double_t oy = (box->GetOrigin())[1];
 224    Double_t oz = (box->GetOrigin())[2];
 225    Double_t point[3], dir[3], lpt[3], ldir[3];
 226    Double_t pstep = 0.;
 227    pstep = TMath::Max(pstep,dz);
 228    Double_t theta, phi;
 229    Int_t idaughter;
 230    timer->Start();
 231    Double_t dist;
 232    Bool_t inside;
 233    for (Int_t i=0; i<1000000; i++) {
 234       lpt[0] = ox-dx+2*dx*gRandom->Rndm();
 235       lpt[1] = oy-dy+2*dy*gRandom->Rndm();
 236       lpt[2] = oz-dz+2*dz*gRandom->Rndm();
 237       gGeoManager->GetCurrentMatrix()->LocalToMaster(lpt,point);
 238       gGeoManager->SetCurrentPoint(point[0],point[1],point[2]);
 239       phi = 2*TMath::Pi()*gRandom->Rndm();
 240       theta= TMath::ACos(1.-2.*gRandom->Rndm());
 241       ldir[0]=TMath::Sin(theta)*TMath::Cos(phi);
 242       ldir[1]=TMath::Sin(theta)*TMath::Sin(phi);
 243       ldir[2]=TMath::Cos(theta);
 244       gGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir,dir);
 245       gGeoManager->SetCurrentDirection(dir);
 246       gGeoManager->SetStep(pstep);
 247       gGeoManager->ResetState();
 248       inside = kTRUE;
 249       dist = TGeoShape::Big();
 250       if (!vol->IsAssembly()) {
 251          inside = vol->Contains(lpt);
 252          if (!inside) {
 253             dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
 254 //            if (dist>=pstep) continue;
 255          } else {
 256             vol->GetShape()->DistFromInside(lpt,ldir,3,pstep);
 257          }
 258
 259          if (!vol->GetNdaughters()) vol->GetShape()->Safety(lpt, inside);
 260       }
 261       if (vol->GetNdaughters()) {
 262          gGeoManager->Safety();
 263          gGeoManager->FindNextDaughterBoundary(point,dir,idaughter,kFALSE);
 264       }
 265    }
 266    timer->Stop();
 267    Double_t time_per_track = timer->CpuTime();
 268    Int_t uid = vol->GetNumber();
 269    Int_t ncrossings = (Int_t)val1[uid];
 270    if (!vol->GetNdaughters())
 271       printf("Time for volume %s (shape=%s): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->GetShape()->GetName(), time_per_track, vol->GetNdaughters(), ncrossings);
 272    else
 273       printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(), time_per_track, vol->GetNdaughters(), ncrossings);
 274    return time_per_track;
 275 }
 276