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

UChar_t KVSpectroDetector::GetPosition( Double_t *XYZf, Int_t idx ){
	// Returns in the array 'XYZf', the coordinates (in cm) of a point randomly drawn in the 
	// active volume with index 'idx'. We assume that the shape of this
	// volume is a TGeoBBox. These coordinates are given in the focal-plane
	// of reference.
	// The function returns a binary code where:
	//   - bit 0 = 1 if X is OK  (001);
	//   - bit 1 = 1 if Y is OK  (010);
	//   - bit 2 = 1 if Z is OK  (100);
	//
	//  If no coordinates are OK, the returned value is null (000) and if 
	// X, Y and Z are OK then the returned value is equal 7 (111).

	TGeoVolume *vol = GetActiveVolume(idx);
	if( !vol ||  !PositionIsOK() ) return 0;

   	const TGeoBBox *box = (TGeoBBox *)vol->GetShape();
   	Double_t dx = box->GetDX();
   	Double_t dy = box->GetDY();
   	Double_t dz = box->GetDZ();
   	Double_t ox = (box->GetOrigin())[0];
   	Double_t oy = (box->GetOrigin())[1];
   	Double_t oz = (box->GetOrigin())[2];
   	Double_t xyz[3];
  	xyz[0] = ox+( TestBit(kRdmPos) ? dx*(2*gRandom->Rndm()-1) : 0.);
   	xyz[1] = oy+( TestBit(kRdmPos) ? dy*(2*gRandom->Rndm()-1) : 0.);
   	xyz[2] = oz+( TestBit(kRdmPos) ? dz*(2*gRandom->Rndm()-1) : 0.);
	if( ActiveVolumeToFocal( xyz , XYZf, idx ) ) return 7;
	return 0;
}

 long dump(TGeoNode* ideal, TGeoNode* aligned,int level, VIDs volids) const {
   char fmt[128];
   string opt_info;
   PlacedVolume pv(ideal);
   bool sensitive = false;
   if ( m_printPositions || m_printVolIDs )  {
     stringstream log;
     if ( m_printPositions )  {
       const double* trans = ideal->GetMatrix()->GetTranslation();
       ::snprintf(fmt, sizeof(fmt), "Pos: (%f,%f,%f) ",trans[0],trans[1],trans[2]);
       log << fmt;
     }
     // Top level volume! have no volume ids
     if ( m_printVolIDs && ideal && ideal->GetMotherVolume() )  {
       VIDs vid = pv.volIDs();
       if ( !vid.empty() )  {
         sensitive = true;
         log << " VolID: ";
         volids.std::vector<VID>::insert(volids.end(),vid.begin(),vid.end());
         for(VIDs::const_iterator i=volids.begin(); i!=volids.end(); ++i)  {
           ::snprintf(fmt, sizeof(fmt), "%s:%2d ",(*i).first.c_str(), (*i).second);
           log << fmt;
         }
       }
     }
     opt_info = log.str();
   }
   TGeoVolume* volume = ideal->GetVolume();
   if ( !m_printSensitivesOnly || (m_printSensitivesOnly && sensitive) )  {
     char sens = pv.volume().isSensitive() ? 'S' : ' ';
     if ( ideal == aligned )  {
       ::snprintf(fmt,sizeof(fmt),"%03d %%-%ds %%s (%%s: %%s) \t[%p] %c %%s",
                  level+1,2*level+1,(void*)ideal, sens);
     }
     else  {
       ::snprintf(fmt,sizeof(fmt),"%03d %%-%ds %%s (%%s: %%s) Ideal:%p Aligned:%p %c %%s",
                  level+1,2*level+1,(void*)ideal,(void*)aligned, sens);
     }
     printout(INFO,"+++",fmt,"",
              aligned->GetName(),
              volume->GetTitle(),
              volume->GetShape()->IsA()->GetName(),
              opt_info.c_str());
   }
   for (Int_t idau = 0, ndau = aligned->GetNdaughters(); idau < ndau; ++idau)  {
     TGeoNode*   ideal_daughter   = ideal->GetDaughter(idau);
     const char* daughter_name    = ideal_daughter->GetName();
     TGeoNode*   aligned_daughter = volume->GetNode(daughter_name);
     dump(ideal_daughter, aligned_daughter, level+1, volids);
   }
   return 1;
 }

void KVSpectroDetector::GetDeltaXYZf(Double_t *DXYZf, Int_t idx){
	// Returns in the DXYZf array the errors of each coordinate of the position returned by
	// GetPosition(...) in the focal-plane frame of reference.
	//
	// In this mother class, the surface of the detector is assumed to be perpendicular to the
	// Z-axis. To be modified in the child classes.

	DXYZf[0] =  DXYZf[1] =  DXYZf[2] = -1;
	TGeoVolume *vol = GetActiveVolume(idx);
	if( !vol ||  !PositionIsOK() ) return;

   	const TGeoBBox *box = (TGeoBBox *)vol->GetShape();
   	DXYZf[0] = box->GetDX();
   	DXYZf[1] = box->GetDY();
   	DXYZf[2] = box->GetDZ();
}

//______________________________________________________________________________
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");
}

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;
}

//.........这里部分代码省略.........
    						distances[j]=daughter->DistanceToIn( conventionalpoints[j], conventionaldirs[j], steps[j]);
    					}
    				}
    				timer.Stop();
    				double t3 = timer.getDeltaSecs();

    				// now unplaced version
    				timer.Start();
    				for(int reps=0;reps<1000;reps++)
    				{
    					for(auto j=0;j<np;++j)
    		    			{
    		    				Vector3D localp, localdir;
    		    				tm->MasterToLocal(conventionalpoints[j], localp);
    		    				tm->MasterToLocalVec(conventionaldirs[j], localdir);
    		    				distances2[j]=unplaceddaughter->DistanceToIn( localp, localdir, steps[j]);
    		    			}
    				}
    				timer.Stop();
    				double t4 = timer.getDeltaSecs();

    				// now unplaced version
    				timer.Start();
    				for(int reps=0;reps<1000;reps++)
    				{
    					for(auto j=0;j<np;++j)
    		      			{
    		       				Vector3D localp, localdir;
    		       				sm->MasterToLocal(conventionalpoints[j], localp);
    		       				sm->MasterToLocalVec(conventionaldirs[j], localdir);
    		       				distances2[j]=unplaceddaughter->DistanceToIn( localp, localdir, steps[j]);
    		      			}
    				}
    				timer.Stop();
    				double t5 = timer.getDeltaSecs();

    				// now unplaced version but inlined matrices
    				timer.Start();
    				for(int reps=0;reps<1000;reps++)
    				{
    					for(auto j=0;j<np;++j)
    					{
    		        	     Vector3D localp, localdir;
    		        	     // this inlines I think
    		           		 tm->MasterToLocal<-1,-1>(conventionalpoints[j], localp);
    		           		 tm->MasterToLocalVec<-1>(conventionaldirs[j], localdir);
    		           		 distances2[j]=unplaceddaughter->DistanceToIn( localp, localdir, 1E30);
    		            }
    				}
    				timer.Stop();
    				double t6 = timer.getDeltaSecs();

    				std::cerr << "SCALAR " << tm->isTranslation() << " " << tm->isRotation() << "("<<tm->getNumberOfZeroEntries()<<")" << " " << t3 <<  " " << t4 << " " << t5 << " " << t6 << std::endl;

    				TGeoMatrix * rootmatrix= new TGeoCombiTrans(TransCases[t][0], TransCases[t][1], TransCases[t][2],
						   new TGeoRotation("rot1",EulerAngles[r][0], EulerAngles[r][1], EulerAngles[r][2]));
    				TGeoManager *geom = new TGeoManager("","");
    		     TGeoVolume * vol = geom->MakeBox("abox",0,10,15,20);
    		     TGeoShape *  rootbox=vol->GetShape();

    		     // now the scalar version from ROOTGeantV
    		     timer.Start();
    		     for(int reps=0;reps<1000;reps++)
    		     {
    		    	 for(auto j=0;j<np;++j)
    		    	 {
    		    		 Vector3D localp, localdir;
    		    		 // this inlines I think
    		    		 rootmatrix->MasterToLocal( &conventionalpoints[j].x, &localp.x );
    		        	 rootmatrix->MasterToLocalVect( &conventionaldirs[j].x, &localdir.x );
    		             distances[j]=rootbox->DistFromOutside( &localp.x, &localdir.x, 3,1e30, 0);
    		         }
    		     }
    		     timer.Stop();
    		     double t7 = timer.getDeltaSecs();

    		     // now the VECTOR version from ROOT
    		     // now the scalar version from ROOTGeantV
    		     timer.Start();
    		     for(int reps=0;reps<1000;reps++)
    		     {
    		    	 rootmatrix->MasterToLocalCombined_v( reinterpret_cast<StructOfCoord const &>(points), reinterpret_cast<StructOfCoord &>(intermediatepoints),
    		    			     		    			 reinterpret_cast<StructOfCoord const &>(dirs), reinterpret_cast<StructOfCoord &>(intermediatedirs), np );
    		         rootbox->DistFromOutsideSOA_v( reinterpret_cast<StructOfCoord const &>(intermediatepoints),
    		        		 	 reinterpret_cast<StructOfCoord const &>(intermediatedirs), 3, steps, 0, distances2, np);
    		     }
    		     timer.Stop();
    		     double t8 = timer.getDeltaSecs();
    		     std::cerr << "RSCAL " << tm->isTranslation() << " " << tm->isRotation() << "("<<tm->getNumberOfZeroEntries()<<")" << " " << t7 << std::endl;
    		     std::cerr << "RVEC " << tm->isTranslation() << " " << tm->isRotation() << "("<<tm->getNumberOfZeroEntries()<<")" << " " << t8 << std::endl;

    		     cmpresults( distances, distances2, np );

    		    delete tm;
    			delete sm;
    		  }

    _mm_free(distances);
    return 1;
}

//////////////////////////////////
// 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 );
//.........这里部分代码省略.........

//////////////////////////////////
// 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;
}

//.........这里部分代码省略.........
//                   for(auto j=0;j<np;++j)
//                         {
//                             Vector3D localp, localdir;
//                             sm->MasterToLocal(conventionalpoints[j], localp);
//                             sm->MasterToLocalVec(conventionaldirs[j], localdir);
//                             distances2[j]=unplaceddaughter->DistanceToIn( localp, localdir, steps[j]);
//                         }
//                }
//                timer.Stop();
//                double t5 = timer.getDeltaSecs();
//
//                // now unplaced version but inlined matrices
//                timer.Start();
//                for(int reps=0;reps<NREPS;reps++)
//                {
//                   for(auto j=0;j<np;++j)
//                   {
//                          Vector3D localp, localdir;
//                          // this inlines I think
//                            tm->MasterToLocal<-1,-1>(conventionalpoints[j], localp);
//                            tm->MasterToLocalVec<-1>(conventionaldirs[j], localdir);
//                            distances2[j]=unplaceddaughter->DistanceToIn( localp, localdir, 1E30);
//                      }
//                }
//                timer.Stop();
//                double t6 = timer.getDeltaSecs();
//
//                std::cerr << "SCALAR " << tm->isTranslation() << " " << tm->isRotation() << "("<<tm->getNumberOfZeroEntries()<<")" << " " << t3 <<  " " << t4 << " " << t5 << " " << t6 << std::endl;
//
              TGeoMatrix * rootmatrix= new TGeoCombiTrans(TransCases[t][0], TransCases[t][1], TransCases[t][2],
                     new TGeoRotation("rot1",EulerAngles[r][0], EulerAngles[r][1], EulerAngles[r][2]));
              TGeoManager *geom = new TGeoManager("","");
               TGeoVolume * vol = geom->MakeTubs("atube",0, tp->GetRmin(), tp->GetRmax(), tp->GetDZ(),tp->GetSPhi() *360/(2.*M_PI), tp->GetSPhi()+360*tp->GetDPhi()/(2.*M_PI));
               TGeoShape *  roottube=vol->GetShape();

               // now the scalar version from ROOTGeantV
           //    timer.Start();
           //    for(int reps=0;reps<NREPS;reps++)
           //    {
           //       for(auto j=0;j<np;++j)
           //       {
           //          Vector3D localp, localdir;
           //          rootmatrix->MasterToLocal( &conventionalpoints[j].x, &localp.x );
           //           rootmatrix->MasterToLocalVect( &conventionaldirs[j].x, &localdir.x );
           //            distancesROOTSCALAR[j]=roottube->DistFromOutside( &localp.x, &localdir.x, 3, Utils::kInfinity, 0);
           //        }
           //    }
           //    timer.Stop();
           //    double t7 = timer.getDeltaSecs();

               // now the scalar version from ROOTGeantV
               timer.Start();
               for(int reps=0;reps<NREPS;reps++)
               {
                  for(auto j=0;j<np;++j)
                  {
                     Vector3D localp, localdir;
                     rootmatrix->MasterToLocal( &conventionalpoints[j].x, &localp.x );
                      rootmatrix->MasterToLocalVect( &conventionaldirs[j].x, &localdir.x );
                      distancesROOTSCALAR[j]=roottube->DistFromInside( &conventionalpoints[j].x, &conventionaldirs[j].x, 3, Utils::kInfinity, 0);
                  }
               }
               timer.Stop();
               double t7 = timer.getDeltaSecs();

               // now the VECTOR version from ROOT