C++ TGeoVolume::GetName方法代码示例

int main() {

  TGeoNode *node =  NULL;
  TGeoVolume *vol = NULL;
  
  LMCgeomN *g = new LMCgeomN("Telescope");
  TObjArray *oa = g->GetGeoManager()->GetListOfNodes();
  for (int i=0; i<oa->GetEntries(); i++) {
    node =  (TGeoNode*)oa->At(i);
    vol = node->GetVolume();
    cout << "= " << node->GetName() << "  " << vol->GetName()  <<endl;
    TObjArray *vnodes = vol->GetNodes();
    cout << vnodes->GetEntries() << endl;
    for (int j=0; j<vnodes->GetEntries(); j++) {
      node = (TGeoNode*)vnodes->At(j);
      cout << "== " << node->GetName() << endl;
      vol = node->GetVolume();
      TObjArray *vnodes1 = vol->GetNodes();
      for (int k=0; k<vnodes1->GetEntries(); k++) {
	node = (TGeoNode*)vnodes1->At(k);
	cout << "=== " << node->GetName() << endl;
	vol = node->GetVolume();
	TObjArray *vnodes2 = vol->GetNodes();
	if(!vnodes2) continue;
	for (int q=0; q<vnodes2->GetEntries(); q++) {
	   node = (TGeoNode*)vnodes2->At(q);
	   cout << "==== " << node->GetName() << endl;
	 }
	
      }
    }
  }

  return 0;
}

TGeoHMatrix &KVSpectroDetector::GetActiveVolToFocalMatrix(Int_t i ) const{ 
	// Returns the matrix which transforms coordinates form the reference
	// frame of the active volume 'i' to the reference frame of the focal
	// plan.
	
	static TGeoHMatrix mm;
	mm.Clear();

	TGeoVolume *vol = (TGeoVolume *)fActiveVolumes->At(i);

	if( vol ){

		Int_t prev_id = gGeoManager->GetCurrentNodeId();
		// To be sure that the matrices will be calculated again
		gGeoManager->CdTop();

		gGeoManager->CdNode( vol->GetUniqueID() );

		// matrix 'volume to target'
		TGeoMatrix *vol_to_tgt = gGeoManager->GetCurrentMatrix();

		if( fFocalToTarget ) mm = fFocalToTarget->Inverse()*(*vol_to_tgt);
		else                 mm = *vol_to_tgt;

		mm.SetName( Form("%s_to_focal",vol->GetName()) );

		// just to not create problems if this method is called during a tracking
		gGeoManager->CdNode( prev_id );
	}
	else mm.SetName("Identity_matrix");

		return mm;
}

//________________________________________________________________________________
Int_t StarMCHits::Init() {
  cout << "StarMCHits::Init() -I- Get Detectors" <<endl;
  if (! fDetectors ) delete fDetectors;
  fDetectors = 0;
  assert(StMaker::GetChain());
  fDetectors = StMaker::GetChain()->GetDataBase("VmcGeometry/Index");
  assert(fDetectors);
  // Make list of detector elements
  TDataSetIter next( fDetectors , 99);
  TDataSet *set = 0;
  if (fDetList) delete fDetList;
  fDetList = new THashList(100,0);
  Int_t N = 0;
  while ((set = next())) {
    StarVMCDetector *det = dynamic_cast<StarVMCDetector *>(set);
    if (! det ) continue;
    if (TString(det->GetName()) == "FPCT") continue; // ignore fpd
    if (TString(det->GetName()) == "BRSG") continue; // ignore tfr
    fDetList->Add(det);
    N++;
  }
  fDetList->Rehash(N);
  cout << "StarMCHits::Init() -I- Get Volume Info" << endl;
  //  TObjArray *UniqueVolumes = gGeoManager->GetListOfUVolumes();
  TObjArray *Volumes = gGeoManager->GetListOfUVolumes();
  Int_t Nids = Volumes->GetEntriesFast();
  if (! fVolUserInfo ) fVolUserInfo = new TObjArray(256);
  for (Int_t i = 0; i < Nids; i++) {
    TGeoVolume *vol  = (TGeoVolume *) Volumes->At(i);
    if (! vol) continue;
    Int_t uid = vol->GetNumber();
#if 0
    cout << "Volume:\t" << i << "\t" << vol->GetName() << "\t" << vol->GetTitle() << "\t" << uid;
#endif
    TString title(vol->GetName(),4);
    TObject *det = fDetList->FindObject(title.Data());
#if 0
    if (det) cout << "\tDetector: " << det->GetName();// << "\t" << det->GetTitle();
#endif
    fVolUserInfo->AddAtAndExpand(det,uid);
#if 0
    cout << endl;
#endif
  }
  return 0;
}

//_____________________________________________________________________________
void StarMCHits::MakeDetectorDescriptors() {
  if (! gGeoManager) {
    cout << "No gGeoManager" << endl;
    return;
  }
  TDataSet *Detectors = StMaker::GetChain()->GetDataBase("VmcGeometry/Index");
  if (! Detectors) {
    cout << "No Detectors found in VmcGeometry/Index" << endl;
  }
  // Make List of sensitive volumes
  TObjArray *vols = gGeoManager->GetListOfVolumes();
  Int_t nvol = vols->GetEntriesFast();
  Int_t nSensVol = 0;
  Int_t nsize = 100;
  TArrayI Indx(nsize); Int_t *indx = Indx.GetArray();
  for (Int_t i = 0; i < nvol; i++) {
    TGeoVolume *vol = (TGeoVolume *) vols->At(i);
    if (! vol) continue;
    TGeoMedium *med = vol->GetMedium();
    if (! med) continue;
    Int_t       isvol = (Int_t) med->GetParam(0);
    if (! isvol) continue;
    indx[nSensVol] = i;
    nSensVol++;
    if (nSensVol >= nsize) {
      nsize *= 2;
      Indx.Set(nsize); 
      indx = Indx.GetArray();
    }
    TString Path(MakeDetectorDescriptor(vol->GetName()));
    if (Detectors) {
      // Check consistency 
      StarVMCDetector *det = (StarVMCDetector *) Detectors->Find(vol->GetName());
      if (! det) {
	cout << "Detector description for " << vol->GetName() << "\t" << vol->GetTitle() << " is missing" << endl;
      } else {
	TString FMT = det->GetFMT();
	cout << "Found path:\t" << Path.Data() << endl;
	cout << "Set   path:\t" << FMT.Data();
	if (Path == FMT) cout << " are the same" << endl;
	else             cout << " are the different" << endl;
      }
    }
    
  }
}

//______________________________________________________________________________
void AddText(TPaveText *pave, TObject *pf, Int_t iaxis)
{
   char line[128];
   TGeoPatternFinder *finder = (TGeoPatternFinder*)pf;
   if (!pave || !pf) return;
   for (Int_t i=0; i<128; i++) line[i] = ' ';
   TGeoVolume *volume = finder->GetVolume();
   TGeoShape *sh = volume->GetShape();
   sprintf(line, "Division of %s on axis %d (%s)", volume->GetName(), iaxis,sh->GetAxisName(iaxis));
   TText *text = pave->AddText(line);
   text->SetTextColor(3);
   text->SetTextAlign(12);
   AddText(pave, "fNdiv",finder->GetNdiv(),"number of divisions");
   AddText(pave, "fStart",finder->GetStart(),"start divisioning position");
   AddText(pave, "fStep",finder->GetStep(),"division step");
}

void scan()
{
  // Load some root geometry
  TGeoVolume *top = gGeoManager->GetTopVolume();
  TGeoIterator iter(top);
  TGeoNode *current;
  // Inspect shape properties
  std::cout << "Top volume: " << top->GetName() << std::endl;
  top->InspectShape();
  TString path;
  while ((current = iter.Next())) {
    iter.GetPath(path);
    std::cout << "=IG====================================" std::cout << path << std::endl;
    std::cout << "=IG====================================" current->InspectNode();
  }
}

GeoHandler& GeoHandler::collect(DetElement element, GeometryInfo& info) {
  m_data->clear();
  i_collect(element.placement().ptr(), 0, Region(), LimitSet());
  for (Data::const_reverse_iterator i = m_data->rbegin(); i != m_data->rend(); ++i) {
    const Data::mapped_type& mapped = (*i).second;
    for (Data::mapped_type::const_iterator j = mapped.begin(); j != mapped.end(); ++j) {
      const TGeoNode* n = *j;
      TGeoVolume* v = n->GetVolume();
      if (v) {
        Material m(v->GetMedium());
        Volume vol = Ref_t(v);
        // Note : assemblies and the world do not have a real volume nor a material
        if (info.volumeSet.find(vol) == info.volumeSet.end()) {
          info.volumeSet.insert(vol);
          info.volumes.push_back(vol);
        }
        if (m.isValid())
          info.materials.insert(m);
        if (dynamic_cast<Volume::Object*>(v)) {
          VisAttr vis = vol.visAttributes();
          //Region      reg = vol.region();
          //LimitSet    lim = vol.limitSet();
          //SensitiveDetector det = vol.sensitiveDetector();

          if (vis.isValid())
            info.vis.insert(vis);
          //if ( lim.isValid() ) info.limits[lim.ptr()].insert(v);
          //if ( reg.isValid() ) info.regions[reg.ptr()].insert(v);
          //if ( det.isValid() ) info.sensitives[det.ptr()].insert(v);
        }
        collectSolid(info, v->GetName(), n->GetName(), v->GetShape(), n->GetMatrix());
      }
    }
  }
  return *this;
}

//________________________________________________________________________________
void StarMCHits::Step() {
  //  static Int_t Idevt0 = -1;
  static Double_t Gold = 0;
#if 0
  if (Debug() && gMC->IsA()->InheritsFrom("TGeant3TGeo")) {
    TGeant3TGeo *geant3 = (TGeant3TGeo *)gMC;
    geant3->Gdebug();
  }
#endif
  //  cout << "Call StarMCHits::Step" << endl;
  TGeoNode *nodeT = gGeoManager->GetCurrentNode();
  assert(nodeT);
  TGeoVolume *volT = nodeT->GetVolume();
  assert(volT);
  const TGeoMedium   *med = volT->GetMedium(); 
  /*   fParams[0] = isvol;
       fParams[1] = ifield;
       fParams[2] = fieldm;
       fParams[3] = tmaxfd;
       fParams[4] = stemax;
       fParams[5] = deemax;
       fParams[6] = epsil;
       fParams[7] = stmin; */
  Int_t Isvol = (Int_t) med->GetParam(0);
  fCurrentDetector = 0;
  if (Isvol <= 0) return;
  fCurrentDetector = (StarVMCDetector *) fVolUserInfo->At(volT->GetNumber());
  if (! fCurrentDetector) {
    volT = nodeT->GetMotherVolume();
    fCurrentDetector = (StarVMCDetector *) fVolUserInfo->At(volT->GetNumber());
    if (! fCurrentDetector) {
      TString path(gGeoManager->GetPath());
      TObjArray *obj = path.Tokenize("_/");
      Int_t N = obj->GetEntries();
      for (Int_t i = N-2; i >= 0; i -= 2) {
	TObjString *o = (TObjString  *) obj->At(i);
	const Char_t *name = o->GetName();
	volT = gGeoManager->GetVolume(name);
	assert (volT);
	fCurrentDetector = (StarVMCDetector *) fVolUserInfo->At(volT->GetNumber());
	if (fCurrentDetector) break;
      }
      delete obj;
    }
  }
  if (Isvol && ! fCurrentDetector && Debug()) {
    cout << "Active medium:" << med->GetName() << "\t for volume " << volT->GetName() 
	 << " has no detector description" << endl;
  }
  //  Int_t Idevt =  gMC->CurrentEvent();
  gMC->TrackPosition(fHit.Current.Global.xyzT);
  gMC->TrackMomentum(fHit.Current.Global.pxyzE);
  TGeoHMatrix  *matrixC = gGeoManager->GetCurrentMatrix();
  fHit.Current.Global2Local(matrixC);
  if (gMC->IsTrackEntering()) {
    fHit.Detector= fCurrentDetector;
    fHit.Entry = fHit.Current;
    fHit.Sleng = gMC->TrackLength();
    fHit.Charge = (Int_t) gMC->TrackCharge();
    fHit.Mass = gMC->TrackMass();
    fHit.AdEstep = fHit.AStep = 0;
    return;
  }
  Double_t GeKin = fHit.Current.Global.pxyzE.E() - fHit.Mass;
  fHit.Sleng = gMC->TrackLength();
  if (fHit.Sleng == 0.) Gold = GeKin;
  Double_t dEstep = gMC->Edep();
  Double_t Step = gMC->TrackStep();
  fHit.iPart = gMC->TrackPid();
  fHit.iTrack = StarVMCApplication::Instance()->GetStack()->GetCurrentTrackId(); // GetCurrentTrackNumber() + 1 to be consistent with g2t
  // - - - - - - - - - - - - - energy correction - - - - - - - - - -
  if (gMC->IsTrackStop() && TMath::Abs(fHit.iPart) == kElectron) {
    TArrayI proc;
    Int_t Nproc = gMC->StepProcesses(proc);
    Int_t Mec = 0;
    for (Int_t i = 0; i < Nproc; i++) if (proc[i] == kPAnnihilation || proc[i] == kPStop) Mec = proc[i];
    Int_t Ngkine = gMC->NSecondaries();
    if (fHit.iPart == kElectron && Ngkine == 0 && Mec == kPStop) dEstep = Gold;
    else {
      if (fHit.iPart == kPositron && Ngkine < 2 && Mec == kPAnnihilation) {
	dEstep = Gold + 2*fHit.Mass;
	if (Ngkine == 1) {
	  TLorentzVector x;
	  TLorentzVector p;
	  Int_t IpartSec;
	  gMC->GetSecondary(0,IpartSec,x,p);
	  dEstep -= p.E();
	}
      }
    }
  }
  // - - - - - - - - - - - - - - - - user - - - - - - - - - - - - - - -
  // user step
  // - - - - - - - - - - - - - - - sensitive - - - - - - - - - - - - -
  fHit.AdEstep += dEstep;  
  fHit.AStep   += Step;
  if (fHit.AdEstep == 0) return;
  if (! gMC->IsTrackExiting() && ! gMC->IsTrackStop()) return;
  fHit.Exit     = fHit.Current;
//.........这里部分代码省略.........

//////////////////////////////////
// main function
int main(int argc, char * argv[])
{
  int axis= 0;

  double axis1_start= 0.;
  double axis1_end= 0.;

  double axis2_start= 0.;
  double axis2_end= 0.;

  double pixel_width= 0;
  double pixel_axis= 1.;

  if( argc < 5 )
  {
    std::cerr<< std::endl;
    std::cerr<< "Need to give rootfile, volumename, axis and number of axis"<< std::endl;
    std::cerr<< "USAGE : ./XRayBenchmarkFromROOTFile [rootfile] [VolumeName] [ViewDirection(Axis)]"
             << "[PixelWidth(OutputImageSize)] [--usolids|--vecgeom(Default:usolids)] [--novoxel(Default:voxel)]"
             << std::endl;
    std::cerr<< "  ex) ./XRayBenchmarkFromROOTFile cms2015.root BSCTrap y 95"<< std::endl;
    std::cerr<< "      ./XRayBenchmarkFromROOTFile cms2015.root PLT z 500 --vecgeom --novoxel"<< std::endl<< std::endl;
    return 1;
  }

  TGeoManager::Import( argv[1] );
  std::string testvolume( argv[2] );

  if( strcmp(argv[3], "x")==0 )
    axis= 1;
  else if( strcmp(argv[3], "y")==0 )
    axis= 2;
  else if( strcmp(argv[3], "z")==0 )
    axis= 3;
  else
  {
    std::cerr<< "Incorrect axis"<< std::endl<< std::endl;
    return 1;
  }

  pixel_width= atof(argv[4]);

  for(auto i= 5; i< argc; i++)
  {
    if( ! strcmp(argv[i], "--usolids") )
      usolids= true;
    if( ! strcmp(argv[i], "--vecgeom") )
      usolids= false;
    if( ! strcmp(argv[i], "--novoxel") )
      voxelize = false;
  }

  int found = 0;
  TGeoVolume * foundvolume = NULL;
  // now try to find shape with logical volume name given on the command line
  TObjArray *vlist = gGeoManager->GetListOfVolumes( );
  for( auto i = 0; i < vlist->GetEntries(); ++i )
  {
    TGeoVolume * vol = reinterpret_cast<TGeoVolume*>(vlist->At( i ));
    std::string fullname(vol->GetName());
    
    std::size_t founds = fullname.compare(testvolume);
    if ( founds==0 ){
      found++;
      foundvolume = vol;

      std::cerr << "("<< i<< ")found matching volume " << foundvolume->GetName()
        << " of type " << foundvolume->GetShape()->ClassName() << "\n";
    }
  }

  std::cerr << "volume found " << found << " times \n\n";

  // if volume not found take world
  if( ! foundvolume ) {
      std::cerr << "specified volume not found; xraying complete detector\n";
      foundvolume = gGeoManager->GetTopVolume();
  }

  if( foundvolume ) {
    foundvolume->GetShape()->InspectShape();
    std::cerr << "volume capacity " 
          << foundvolume->GetShape()->Capacity() << "\n";

    // get bounding box to generate x-ray start positions
    double dx = ((TGeoBBox*)foundvolume->GetShape())->GetDX()*1.05;
    double dy = ((TGeoBBox*)foundvolume->GetShape())->GetDY()*1.05;
    double dz = ((TGeoBBox*)foundvolume->GetShape())->GetDZ()*1.05;
    double origin[3]= {0., };
    origin[0]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[0];
    origin[1]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[1];
    origin[2]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[2];
    
    TGeoMaterial * matVacuum = new TGeoMaterial("Vacuum",0,0,0);
    TGeoMedium * vac = new TGeoMedium("Vacuum",1,matVacuum);

    TGeoVolume* boundingbox= gGeoManager->MakeBox("BoundingBox", vac,
            std::abs(origin[0]) + dx,
//.........这里部分代码省略.........

//////////////////////////////////
// main function
int main(int argc, char * argv[])
{
  int axis= 0;

  double axis1_start= 0.;
  double axis1_end= 0.;

  double axis2_start= 0.;
  double axis2_end= 0.;

  double pixel_axis= 1.;

  if( argc < 5 )
  {
    std::cerr<< std::endl;
    std::cerr<< "Need to give rootfile, volumename, direction phi and direction theta (in degrees)"<< std::endl;
    return 1;
  }

  TGeoManager::Import( argv[1] );
  std::string testvolume( argv[2] );

  //double directionphi   = atof(argv[3])*vecgeom::kDegToRad;
  //double directiontheta = atof(argv[4])*vecgeom::kDegToRad;

  for(auto i= 5; i< argc; i++)
  {
    if( ! strcmp(argv[i], "--usolids") )
      usolids= true;
    if( ! strcmp(argv[i], "--vecgeom") )
      usolids= false;
    if( ! strcmp(argv[i], "--novoxel") )
      voxelize = false;
  }

  int found = 0;
  TGeoVolume * foundvolume = NULL;
  // now try to find shape with logical volume name given on the command line
  TObjArray *vlist = gGeoManager->GetListOfVolumes( );
  for( auto i = 0; i < vlist->GetEntries(); ++i )
  {
    TGeoVolume * vol = reinterpret_cast<TGeoVolume*>(vlist->At( i ));
    std::string fullname(vol->GetName());
    
    std::size_t founds = fullname.compare(testvolume);
    if ( founds==0 ){
      found++;
      foundvolume = vol;

      std::cerr << "("<< i<< ")found matching volume " << foundvolume->GetName()
        << " of type " << foundvolume->GetShape()->ClassName() << "\n";
    }
  }

  std::cerr << "volume found " << found << " times \n\n";

  // if volume not found take world
  if( ! foundvolume ) {
      std::cerr << "specified volume not found; xraying complete detector\n";
      foundvolume = gGeoManager->GetTopVolume();
  }

  if( foundvolume ) {
    foundvolume->GetShape()->InspectShape();
    std::cerr << "volume capacity " 
          << foundvolume->GetShape()->Capacity() << "\n";

    // get bounding box to generate x-ray start positions
    double dx = ((TGeoBBox*)foundvolume->GetShape())->GetDX()*1.5;
    double dy = ((TGeoBBox*)foundvolume->GetShape())->GetDY()*1.5;
    double dz = ((TGeoBBox*)foundvolume->GetShape())->GetDZ()*1.5;
    double origin[3]= {0., };
    origin[0]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[0];
    origin[1]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[1];
    origin[2]= ((TGeoBBox*)foundvolume->GetShape())->GetOrigin()[2];
    
    TGeoMaterial * matVacuum = new TGeoMaterial("Vacuum",0,0,0);
    TGeoMedium * vac = new TGeoMedium("Vacuum",1,matVacuum);

    TGeoVolume* boundingbox= gGeoManager->MakeBox("BoundingBox", vac,
            std::abs(origin[0]) + dx,
            std::abs(origin[1]) + dy,
            std::abs(origin[2]) + dz );

    // TGeoManager * geom = boundingbox->GetGeoManager();
    std::cout << gGeoManager->CountNodes() << "\n";

    if(! voxelize ) DeleteROOTVoxels();

//    TGeoManager * mg1 = gGeoManager;
    gGeoManager = 0;

    TGeoManager * mgr2 = new TGeoManager();

//    delete gGeoManager;
//    gGeoManager = new TGeoManager();
    boundingbox->AddNode( foundvolume, 1);
    mgr2->SetTopVolume( boundingbox );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  pipe->AddNode(pipe3, 0);
  TGeoVolume* pipevac2 = MakeVacuum(2, nSects02, z02, rin02, rout02, vacuum,     &infoFile); 
  pipevac2->SetLineColor(kCyan);
  pipe->AddNode(pipevac2, 0);
  //*/
  
  // define some rotation & translation for pipe4-pipe7 & pipevac3-pipevac5
  // Initial phi   (Euler rotation angle about Z)
  // Initial theta (Euler rotation angle about new X)
  // Initial psi   (Euler rotation angle about new Z)
  cout<<endl<<"Rotation angle (@[0;0;1800]mm)="<< rotangle <<"deg";
  TGeoRotation    *r1 = new TGeoRotation("r1", 90., rotangle, 0.); // deg - Euler angles
  r1->RegisterYourself();
  Double_t shift = rotangle/10.+0.1; // cm - volume boolean operation correction - TODO
  Double_t dx=0., dy=0., dz=1800.-10.*shift; // mm
  TGeoCombiTrans  *c1 = new TGeoCombiTrans("c1", dx/10., dy/10., dz/10., r1); // mm->cm
  c1->RegisterYourself();
  TGeoTranslation *t1 = new TGeoTranslation("t1", 0., 0., 180.-shift); // cm
  t1->RegisterYourself();
  Double_t phi, theta, psi;
  c1->GetRotation()->GetAngles(phi, theta, psi);
  
  infoFile << endl << "Beam pipe section: " << pipe4name << endl;
  infoFile << "Traslation(dx,dy,dz):    "<< dx <<" "<< dy <<" "<< dz <<" mm" << endl;
  infoFile << "Rotation(phi,theta,psi): " << phi <<" "<< theta <<" "<< psi <<" deg"<< endl;
  infoFile << setw(2) << "i" << setw(10) << "Z,mm" << setw(10) << "Rin,mm" << setw(10) << "Rout,mm" << setw(10) << "h,mm" << endl;
  // create shape components for pipe4
  TGeoPcon* p4 = MakeShape(nSects4, "p4", z4, rin4,  rout4, &infoFile);
  TGeoBBox* b  = new TGeoBBox("b", 7., 7., shift); // cm 0.13
  // create a composite for pipe4
  //*
  TGeoCompositeShape *cs1 = new TGeoCompositeShape("cs1", "p4:r1-b");
  TGeoVolume *pipe4 = new TGeoVolume("pipe4",cs1, pipeMedium);
  pipe4->SetLineColor(kGreen+2);
  pipe->AddNode(pipe4, 0, t1);
  //*/
  // create shape components for pipevac3
  TGeoPcon* p03 = MakeShape(nSects03, "p03", z03, rin03, rout03, &infoFileEmpty);
  // create a composite for pipevac3
  TGeoCompositeShape *cs2 = new TGeoCompositeShape("cs2", "p03:r1-b");
  TGeoVolume *pipevac3 = new TGeoVolume("pipevac3",cs2, vacuum);
  pipevac3->SetLineColor(kCyan);
  pipe->AddNode(pipevac3, 0, t1);
  //*
  infoFile << endl << "Beam pipe section: " << pipe5name << endl;
  infoFile << "Traslation(dx,dy,dz):    "<< dx <<" "<< dy <<" "<< dz <<" mm" << endl;
  infoFile << "Rotation(phi,theta,psi): " << phi <<" "<< theta <<" "<< psi <<" deg"<< endl;
  infoFile << setw(2) << "i" << setw(10) << "Z,mm" << setw(10) << "Rin,mm" << setw(10) << "Rout,mm" << setw(10) << "h,mm" << endl;
  TGeoVolume* pipe5    = MakePipe  (5, nSects5,  z5,  rin5,  rout5,  pipeMedium, &infoFile); 
  pipe5->SetLineColor(kGreen);
  pipe->AddNode(pipe5, 0, c1);
  TGeoVolume* pipevac4 = MakeVacuum(4, nSects04, z04, rin04, rout04, vacuum,     &infoFile); 
  pipevac4->SetLineColor(kCyan);
  pipe->AddNode(pipevac4, 0, c1);
  
  infoFile << endl << "Beam pipe section: " << pipe6name << endl;
  infoFile << "Traslation(dx,dy,dz):    "<< dx <<" "<< dy <<" "<< dz <<" mm" << endl;
  infoFile << "Rotation(phi,theta,psi): " << phi <<" "<< theta <<" "<< psi <<" deg"<< endl;
  infoFile << setw(2) << "i" << setw(10) << "Z,mm" << setw(10) << "Rin,mm" << setw(10) << "Rout,mm" << setw(10) << "h,mm" << endl;
  TGeoVolume* pipe6    = MakePipe  (6, nSects6,  z6,  rin6,  rout6,  pipeMedium, &infoFile); 
  pipe6->SetLineColor(kBlue);
  pipe->AddNode(pipe6, 0, c1);
  
  infoFile << endl << "Beam pipe section: " << pipe7name << endl;
  infoFile << "Traslation(dx,dy,dz):    "<< dx <<" "<< dy <<" "<< dz <<" mm" << endl;
  infoFile << "Rotation(phi,theta,psi): " << phi <<" "<< theta <<" "<< psi <<" deg"<< endl;
  infoFile << setw(2) << "i" << setw(10) << "Z,mm" << setw(10) << "Rin,mm" << setw(10) << "Rout,mm" << setw(10) << "h,mm" << endl;
  TGeoVolume* pipe7    = MakePipe  (7, nSects7,  z7,  rin7,  rout7,  pipeMedium, &infoFile); 
  pipe7->SetLineColor(kGreen);
  pipe->AddNode(pipe7, 0, c1);
  TGeoVolume* pipevac5 = MakeVacuum(5, nSects05, z05, rin05, rout05, vacuum,     &infoFile); 
  pipevac5->SetLineColor(kCyan);
  pipe->AddNode(pipevac5, 0, c1);
  //*/
  
  // -----   End   --------------------------------------------------

  // ---------------   Finish   -----------------------------------------------
  top->AddNode(pipe, 1);
  cout << endl << endl;
  gGeoMan->CloseGeometry();
  gGeoMan->CheckOverlaps(0.001);
  gGeoMan->PrintOverlaps();
  gGeoMan->Test();
  
  // visualize it with ray tracing, OGL/X3D viewer
  //top->Raytrace();
  top->Draw("ogl");
  //top->Draw("x3d");

  TFile* rootFile = new TFile(rootFileName, "RECREATE");
  top->Write();
  cout << endl;
  cout << "Geometry " << top->GetName() << " written to " 
       << rootFileName << endl;
  rootFile->Close();
  infoFile.close();


}

//////////////////////////////////
// main function
int main(int argc, char *argv[])
{
  if (argc < 3) {
    std::cerr << std::endl;
    std::cerr << "Need to give rootfile + volumename" << std::endl;
    return 1;
  }

  TGeoManager::Import(argv[1]);
  std::string testvolume(argv[2]);

  int found               = 0;
  TGeoVolume *foundvolume = NULL;
  // now try to find shape with logical volume name given on the command line
  TObjArray *vlist = gGeoManager->GetListOfVolumes();
  for (auto i = 0; i < vlist->GetEntries(); ++i) {
    TGeoVolume *vol = reinterpret_cast<TGeoVolume *>(vlist->At(i));
    std::string fullname(vol->GetName());
    std::size_t founds = fullname.compare(testvolume);
    if (founds == 0) {
      found++;
      foundvolume = vol;
      std::cerr << "(" << i << ")found matching volume " << foundvolume->GetName() << " of type "
                << foundvolume->GetShape()->ClassName() << "\n";
    }
  }
  std::cerr << "volume found " << found << " times \n\n";

  // if volume not found take world
  if (!foundvolume) {
    std::cerr << "specified volume not found; exiting\n";
    return 1;
  }

  if (foundvolume) {
    // convert current gGeoManager to a VecGeom geometry
    VPlacedVolume const *vecgeompvol = RootGeoManager::Instance().Convert(foundvolume)->Place();

    for (int i = 0; i < 20; ++i) {
      Vector3D<Precision> point = vecgeompvol->GetUnplacedVolume()->SamplePointOnSurface();

      if (vecgeompvol->Inside(point) != vecgeom::kSurface) {
        std::cerr << " WARNING : Inside does not report surface state \n";
      }
      Vector3D<Precision> dir = volumeUtilities::SampleDirection();
      bool contained          = vecgeompvol->Contains(point);
      bool exiting            = ExitingMethod2(vecgeompvol, point, dir);
      double DO               = vecgeompvol->DistanceToOut(point, dir);
      double DI               = vecgeompvol->DistanceToIn(point, dir);

      std::cerr << i << "  " << point << "  "
                << contained
                //                  << "   ExitingM1  " << ExitingMethod1(vecgeompvol,point,dir)
                << "   ExitingM2  " << exiting << " DI " << DI << " DO " << DO << " SO "
                << vecgeompvol->SafetyToOut(point) << " SI " << vecgeompvol->SafetyToIn(point) << "\n";

      if (exiting && DO > vecgeom::kTolerance) {
        std::cout << " WARNING FOR DO  : should be zero \n";
        std::cout << "./CompareDistances " << argv[1] << " " << argv[2] << std::setprecision(20) << " " << point[0]
                  << " " << point[1] << " " << point[2] << " " << dir[0] << " " << dir[1] << " " << dir[2] << "\n";
      }
      if (contained && !exiting && DO < vecgeom::kTolerance) {
        std::cout << " FATAL WARNING FOR DO : should be finite \n";
        std::cout << "./CompareDistances " << argv[1] << " " << argv[2] << std::setprecision(20) << " " << point[0]
                  << " " << point[1] << " " << point[2] << " " << dir[0] << " " << dir[1] << " " << dir[2] << "\n";
      }

      if (!exiting && DI > vecgeom::kTolerance) {
        std::cout << " WARNING FOR DI : should be zero \n";
        std::cout << "./CompareDistances " << argv[1] << " " << argv[2] << std::setprecision(20) << " " << point[0]
                  << " " << point[1] << " " << point[2] << " " << dir[0] << " " << dir[1] << " " << dir[2] << "\n";
      }
      if (!exiting && !contained && DI > 1E20) {
        std::cout << " FATAL WARNING FOR DI : should be zero \n";
        std::cout << "./CompareDistances " << argv[1] << " " << argv[2] << std::setprecision(20) << " " << point[0]
                  << " " << point[1] << " " << point[2] << " " << dir[0] << " " << dir[1] << " " << dir[2] << "\n";
      }

      std::cout << "\n\n";
    } // end for
    return 0;
  } // end if found volume
  return 1;
}

//.........这里部分代码省略.........
   Double_t larr_VolumeIsotopes[lcin_Z][lcin_A] = {0.}; //[Z][A], no map in pure ROOT

   for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
   {
      TGeoMaterial *lgeo_Mat = gGeoManager->GetMaterial(i);
      Int_t    idx     = gGeoManager->GetMaterial(i)->GetIndex();

      if( lgeo_Mat->IsMixture() )
      {
         TGeoMixture * lgeo_Mix = dynamic_cast <TGeoMixture*> ( lgeo_Mat );
         Int_t lint_Nelements = lgeo_Mix->GetNelements();

         for ( Int_t j=0; j<lint_Nelements; j++) 
         {
            Int_t lint_Z = TMath::Nint( (Double_t) lgeo_Mix->GetZmixt()[j] );
            Int_t lint_A = TMath::Nint( (Double_t) lgeo_Mix->GetAmixt()[j] );
            Double_t ldou_Fraction = lgeo_Mix->GetWmixt()[j];
            Double_t ldou_Density = lgeo_Mix->GetDensity() * density_unit_to_SI;

            larr_MassIsotopes[ lint_Z ][ lint_A ] += volume[idx] * ldou_Fraction * ldou_Density;
            larr_VolumeIsotopes[ lint_Z ][ lint_A ] += volume[idx] * ldou_Fraction;
         }
      }
   }

   //
   // print out volume/mass for each `material'
   //

   Double_t ldou_MinimumVolume = 1e-20;

   cout << endl
        << " Geometry: \"" <<  gFileName << "\"" << endl
        << " TopVolume: \"" << topvol->GetName() << "\"" 
        << endl;

   cout <<endl << "materials:" << endl;
   cout << setw(5) << "index"
        << setw(15) << "name"
        << setprecision(6) 
        << setw(14) << "volume (m^3)"
        << setw(14) << "mass (kg)"
        << setw(14) << "mass (%)"
        <<  endl;

   double total_mass_materials = 0;
   for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
   {
     Int_t    idx     = gGeoManager->GetMaterial(i)->GetIndex();
     Double_t density = gGeoManager->GetMaterial(i)->GetDensity() * density_unit_to_SI;
     Double_t mass_material = density * volume[idx];
     if ( volume[idx] > ldou_MinimumVolume ) {
       total_mass_materials += mass_material;
     }
   }


   for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
   {
      Int_t    idx     = gGeoManager->GetMaterial(i)->GetIndex();
      Double_t density = gGeoManager->GetMaterial(i)->GetDensity() * density_unit_to_SI;

      mass[idx] = density * volume[idx];

      if( volume[idx] > ldou_MinimumVolume ) {
        cout << setw(5) << i