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

//.........这里部分代码省略.........
  //GFI2 position
  //dx = -145.270039; //Justyna
  //dy = 0.000000;    //Justyna
  //dz = 730.318956;  //Justyna
  //dx = -145.212570; //Justyna new
  //dy = -0.010000;    //Justyna new
  //dz = 730.336197;  //Justyna new
  //dx = -147.486; //dE tracker
  //dy = 0.070;
  //dz = 729.798;
  
  //LABPOS(GFI2,-147.135037,0.069976,729.680342)
  dx = -147.135037; //dE tracker, correction due to wrong angle
  dy = 0.069976;
  dz = 729.680342;
  
  /*
   dx = 141.800000;
   dy = 0.000000;
   dz = 727.300000;
   */
  // Rotation:
  thx = -106.700000;    phx = 0.000000;
  //   thx = 106.700000;    phx = 0.000000;
  thy = 90.000000;    phy = 90.000000;
  thz = -16.700000;    phz = 0.000000;
  TGeoRotation *pMatrix5 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
  TGeoCombiTrans*
  pMatrix4 = new TGeoCombiTrans("", dx,dy,dz,pMatrix5);
  
  
  // World definition
  TGeoVolume* pWorld = gGeoManager->GetTopVolume();
  pWorld->SetVisLeaves(kTRUE);
  
  
  
  
  // SHAPES, VOLUMES AND GEOMETRICAL HIERARCHY
  // Shape: GFIBoxWorld type: TGeoBBox
  dx = 29.00000;
  dy = 29.00000;
  dz = 0.050000;
  TGeoShape *pGFIBoxWorld = new TGeoBBox("GFIBoxWorld", dx,dy,dz);
  // Volume: GFILogWorld
  TGeoVolume*
  pGFILogWorld = new TGeoVolume("GFILogWorld",pGFIBoxWorld, pMed2);
  pGFILogWorld->SetVisLeaves(kTRUE);
  
  // Global positioning
  
  TGeoCombiTrans *pGlobal1 = GetGlobalPosition(pMatrix2);
  TGeoCombiTrans *pGlobal2 = GetGlobalPosition(pMatrix4);
    
  TGeoVolumeAssembly *pw1 = new TGeoVolumeAssembly("GFI");
  pWorld->AddNode(pw1, 0);
    
  pw1->AddNode( pGFILogWorld, 0, pGlobal1 );
  pw1->AddNode( pGFILogWorld, 1, pGlobal2 );
  
  
  
  // Shape: GFIBox type: TGeoBBox
  dx = 25.000000;
  dy = 25.000000;
  dz = 0.050000;

//.........这里部分代码省略.........
    Int_t nel, numed;



    // TRANSFORMATION MATRICES
    // Combi transformation:
    dx = -417.359574; //Justyna
    dy = 2.400000;    //Justyna
    dz = 960.777114;  //Justyna

    //    dx = -421.33683; //Christoph
    //    dy = 2.12;    //Christoph
    //    dz = 958.387337;  //Christoph

    /*    dx = -171.1;  //position directrly (15cm) after DCH2
     dy = 2.400000;
     dz = 548.95;*/
    // dz = 0.;
    /*   // Rotation:
     thx = -121.000000;    phx = 0.000000;
     thy = 90.000000;    phy = 90.000000;
     thz = -31.000000;    phz = 0.000000;*/  //this

    //   TGeoRotation *pMatrix3 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);  //this
    TGeoRotation *gRot = new TGeoRotation();
    gRot->RotateX(0.);
    gRot->RotateY(-31.000000);
    gRot->RotateZ(0.);

    TGeoCombiTrans*
    //   pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,pMatrix3);  //this
    pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,gRot);  //this







    /* PREVIOUS!!!
     // TRANSFORMATION MATRICES
     // Combi transformation:
     dx = 419.700000;
     dy = 0.000000;
     dz = 952.400000;
     // dz = 0.;
     // Rotation:
     thx = 121.000000;    phx = 0.000000;
     thy = 90.000000;    phy = 90.000000;
     thz = 31.000000;    phz = 0.000000;
     TGeoRotation *pMatrix3 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
     TGeoCombiTrans*
     pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,pMatrix3);
     */

    //Top Volume
    TGeoVolume* pWorld = gGeoManager->GetTopVolume();
    pWorld->SetVisLeaves(kTRUE);


    TGeoVolumeAssembly *ptof = new TGeoVolumeAssembly("TOF");


    // SHAPES, VOLUMES AND GEOMETRICAL HIERARCHY
    // Shape: TOFBox type: TGeoBBox
    dx = 94.450000;
    dy = 73.450000;
    dz = 0.500000;
    TGeoShape *pTOFBox = new TGeoBBox("TOFBox", dx,dy,dz);
    // Volume: TOFLog
    TGeoVolume*
    pTOFLog = new TGeoVolume("TOFLog",pTOFBox, pMed34);
    pTOFLog->SetVisLeaves(kTRUE);

    TGeoCombiTrans *t0 = new TGeoCombiTrans("t0");
    ptof->AddNode(pTOFLog, 0, t0);

    TGeoCombiTrans *pGlobal = GetGlobalPosition(pMatrix2);
    if (pGlobal) {
        pWorld->AddNode(ptof, 0, pGlobal);
    } else {
        pWorld->AddNode(ptof, 0, pMatrix2);
    }

//  AddSensitiveVolume(pTOFLog);
//  fNbOfSensitiveVol+=1;



    // ---------------   Finish   -----------------------------------------------
    gGeoMan->CloseGeometry();
    gGeoMan->CheckOverlaps(0.001);
    gGeoMan->PrintOverlaps();
    gGeoMan->Test();

    TFile* geoFile = new TFile(geoFileName, "RECREATE");
    top->Write();
    geoFile->Close();
    // --------------------------------------------------------------------------
}

//.........这里部分代码省略.........
  dz = 761.870346;
  
  // Rotation:
  thx = -106.700000;    phx = 0.000000;
  thy = 90.000000;    phy = 90.000000;
  thz = -16.700000;    phz = 0.000000;
  
  /*    dx = -171.1;
   dy = 2.400000;
   dz = 548.95;
   // dz = 0.;
   // Rotation:
   thx = -121.000000;    phx = 0.000000;
   thy = 90.000000;    phy = 90.000000;
   thz = -31.000000;    phz = 0.000000;*/
  
  TGeoRotation *pMatrix3 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
  TGeoCombiTrans*
  pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,pMatrix3);
  
  /*
   // TRANSFORMATION MATRICES
   // Combi transformation:
   dx = 151.000000;
   dy = 0.000000;
   dz = 758.000000;
   // Rotation:
   thx = 106.700000;    phx = 0.000000;
   thy = 90.000000;    phy = 90.000000;
   thz = 16.700000;    phz = 0.000000;
   TGeoRotation *pMatrix3 = new TGeoRotation("",thx,phx,thy,phy,thz,phz);
   TGeoCombiTrans*
   pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,pMatrix3);
   */
  
  // Shape: World type: TGeoBBox
  TGeoVolume* pWorld = gGeoManager->GetTopVolume();
  pWorld->SetVisLeaves(kTRUE);
  
  // Create a global Mother Volume
  /*
   dx = 200.000000;
   dy = 200.000000;
   dz = 200.000000;
   TGeoShape *pBoxWorld = new TGeoBBox("mTofBoxWorld", dx,dy,dz);
   TGeoVolume*
   pWorld  = new TGeoVolume("mTofBoxLogWorld",pBoxWorld, pMed2);
   pWorld->SetVisLeaves(kTRUE);
   TGeoCombiTrans *pGlobalc = GetGlobalPosition();
   
   // add the sphere as Mother Volume
   pAWorld->AddNode(pWorld, 0, pGlobalc);
   */
  
  
  // SHAPES, VOLUMES AND GEOMETRICAL HIERARCHY
  // Shape: mTOFBox type: TGeoBBox
  dx = 24.000000;
  dy = 24.000000;
  //dz = 0.250000;	//wrong: should be 0.5->1cm total
  dz = 0.500000;
  
  /*   dx = 94.450000;  //TFW size
   dy = 73.450000;
   dz = 0.500000;*/
  
  TGeoShape *pmTOFBox = new TGeoBBox("mTOFBox", dx,dy,dz);
  // Volume: mTOFLog
  TGeoVolume *
  pmTOFLog = new TGeoVolume("mTOFLog",pmTOFBox, pMed34);
  pmTOFLog->SetVisLeaves(kTRUE);
  
  TGeoVolumeAssembly *pmTof = new TGeoVolumeAssembly("mTOF");
  TGeoCombiTrans *t0 = new TGeoCombiTrans("t0");
  pmTof->AddNode(pmTOFLog, 0, t0);
  
  TGeoCombiTrans *pGlobal = GetGlobalPosition(pMatrix2);
  if (pGlobal){
    pWorld->AddNode(pmTof, 0, pGlobal);
  }else{
    pWorld->AddNode(pmTof, 0, pMatrix2);
  }
  
  
//  AddSensitiveVolume(pmTOFLog);
//  fNbOfSensitiveVol+=1;
  
  
  
  // ---------------   Finish   -----------------------------------------------
  gGeoMan->CloseGeometry();
  gGeoMan->CheckOverlaps(0.001);
  gGeoMan->PrintOverlaps();
  gGeoMan->Test();

  TFile* geoFile = new TFile(geoFileName, "RECREATE");
  top->Write();
  geoFile->Close();
  // --------------------------------------------------------------------------
}

void create_dch_geo(const char* geoTag)
{

  fGlobalTrans->SetTranslation(0.0,0.0,0.0);

  // -------   Load media from media file   -----------------------------------
  FairGeoLoader*    geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  TString geoPath = gSystem->Getenv("VMCWORKDIR");
  TString medFile = geoPath + "/geometry/media_r3b.geo";
  geoFace->setMediaFile(medFile);
  geoFace->readMedia();
  gGeoMan = gGeoManager;
  // --------------------------------------------------------------------------



  // -------   Geometry file name (output)   ----------------------------------
  TString geoFileName = geoPath + "/geometry/dch_";
  geoFileName = geoFileName + geoTag + ".geo.root";
  // --------------------------------------------------------------------------



  // -----------------   Get and create the required media    -----------------
  FairGeoMedia*   geoMedia = geoFace->getMedia();
  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();

  FairGeoMedium* mAir      = geoMedia->getMedium("Air");
  if ( ! mAir ) Fatal("Main", "FairMedium Air not found");
  geoBuild->createMedium(mAir);
  TGeoMedium* pMed2 = gGeoMan->GetMedium("Air");
  if ( ! pMed2 ) Fatal("Main", "Medium Air not found");

  FairGeoMedium* mAl      = geoMedia->getMedium("aluminium");
  if ( ! mAl ) Fatal("Main", "FairMedium aluminium not found");
  geoBuild->createMedium(mAl);
  TGeoMedium* pMed21 = gGeoMan->GetMedium("aluminium");
  if ( ! pMed21 ) Fatal("Main", "Medium aluminium not found");

  FairGeoMedium* mMylar      = geoMedia->getMedium("mylar");
  if ( ! mMylar ) Fatal("Main", "FairMedium mylar not found");
  geoBuild->createMedium(mMylar);
  TGeoMedium* pMed15 = gGeoMan->GetMedium("mylar");
  if ( ! pMed15 ) Fatal("Main", "Medium mylar not found");

  FairGeoMedium* mHe      = geoMedia->getMedium("helium");
  if ( ! mHe ) Fatal("Main", "FairMedium helium not found");
  geoBuild->createMedium(mHe);
  TGeoMedium* pMed4 = gGeoMan->GetMedium("helium");
  if ( ! pMed4 ) Fatal("Main", "Medium helium not found");

  FairGeoMedium* mDCH      = geoMedia->getMedium("DCHgas");
  if ( ! mDCH ) Fatal("Main", "FairMedium DCHgas not found");
  geoBuild->createMedium(mDCH);
  TGeoMedium* pMed33 = gGeoMan->GetMedium("DCHgas");
  if ( ! pMed33 ) Fatal("Main", "Medium DCHgas not found");
  // --------------------------------------------------------------------------



  // --------------   Create geometry and top volume  -------------------------
  gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
  gGeoMan->SetName("DCHgeom");
  TGeoVolume* top = new TGeoVolumeAssembly("TOP");
  gGeoMan->SetTopVolume(top);
  // --------------------------------------------------------------------------



  Double_t dx,dy,dz;
  Double_t theta, phi;
  Double_t thx, phx, thy, phy, thz, phz;
  Double_t alpha;



  // Define DCH Geometry
  TGeoVolume* topDCH = new TGeoVolumeAssembly("DCH");


  //Active volume	(envelope box enclosing the sense wires)	
  Double_t actGasDx = 99.07/2.; // [cm]	//99.07326
  Double_t actGasDy = 77.60/2.;  // [cm]	//77.59584
  Double_t actGasDz = 5.20/2.;     // [cm]
  //distance of x and y planes centres is 4 cm. Distance between the first x halfplane and the second y halfplane is 4cm + 2*0.69*cos30 cm.

  TGeoShape* pActGasBox = new TGeoBBox("ActGasBox",
				       actGasDx,
				       actGasDy,
				       actGasDz);
  TGeoVolume*
    pActGasDchLog = new TGeoVolume("ActGASBoxLog",pActGasBox, pMed33);
  pActGasDchLog->SetVisLeaves(kTRUE);
  pActGasDchLog->SetVisContainers(kTRUE);


  // Gas box
  Double_t gasDx = 106.4/2.; // [cm]
  Double_t gasDy = 83.4/2.;  // [cm]
//.........这里部分代码省略.........

/**
 * Initialise ROOT geometry objects from GEAR objects
 * 
 * @param geomName name of ROOT geometry object
 * @param dumpRoot dump automatically generated ROOT geometry file for further inspection
 */
void EUTelGeometryTelescopeGeoDescription::initializeTGeoDescription( std::string const & geomName, bool dumpRoot = false ) {
	if( _isGeoInitialized ) {
		streamlog_out( WARNING3 ) << "EUTelGeometryTelescopeGeoDescription: Geometry already initialized, using old initialization" << std::endl;
		return;
	} else {
    		_geoManager = std::make_unique<TGeoManager>("Telescope", "v0.1");
			_geoManager->SetBit(kCanDelete);
	}

	if( !_geoManager ) {
		streamlog_out( ERROR3 ) << "Can't instantiate ROOT TGeoManager " << std::endl;
		return;
	}
   
    
    // Create top world volume containing telescope geometry
    
    // Create air mixture
    // see http://pdg.lbl.gov/2013/AtomicNuclearProperties/HTML_PAGES/104.html
    double air_density = 1.2e-3;         // g/cm^3
    double air_radlen  = 36.62;          // g/cm^2 //Must be -ve to use this value rather than internal one calculated.

    TGeoMixture* pMatAir = new TGeoMixture("AIR",3,air_density);
    pMatAir->DefineElement(0, 14.007, 7.,  0.755267 );     //Nitrogen
    pMatAir->DefineElement(1, 15.999, 8.,  0.231781 );     //Oxygen
    pMatAir->DefineElement(2, 39.948, 18., 0.012827 );     //Argon
    pMatAir->DefineElement(3, 12.011, 6.,  0.000124 );     //Carbon
    pMatAir->SetRadLen( air_radlen );
    // Medium: medium_World_AIR
    TGeoMedium* pMedAir = new TGeoMedium("medium_World_AIR", 3, pMatAir );

    // The World is the 10 x 10m x 10m box filled with air mixture
    Double_t dx,dy,dz;
    dx = 5000.000000; // [mm]
    dy = 5000.000000; // [mm]
    dz = 5000.000000; // [mm]
    TGeoShape *pBoxWorld = new TGeoBBox("Box_World", dx,dy,dz);
    // Volume: volume_World
    TGeoVolume* pvolumeWorld = new TGeoVolume("volume_World",pBoxWorld, pMedAir);
    pvolumeWorld->SetLineColor(4);
    pvolumeWorld->SetLineWidth(3);
    pvolumeWorld->SetVisLeaves(kTRUE);

   // Set top volume of geometry
   gGeoManager->SetTopVolume( pvolumeWorld );
   
   IntVec::const_iterator itrPlaneId;
   for ( itrPlaneId = _sensorIDVec.begin(); itrPlaneId != _sensorIDVec.end(); ++itrPlaneId ) {
       translateSiPlane2TGeo(pvolumeWorld, *itrPlaneId );
   }
 
    _geoManager->CloseGeometry();
    _isGeoInitialized = true;
    // Dump ROOT TGeo object into file
    if ( dumpRoot ) _geoManager->Export( geomName.c_str() );

   for(auto& mapEntry: _planePath) {
		auto pathName = mapEntry.second;
		auto sensorID = mapEntry.first;
    	_geoManager->cd( pathName.c_str() );
		_TGeoMatrixMap[sensorID] = _geoManager->GetCurrentNode()->GetMatrix();
	} 

    return;
}

//.........这里部分代码省略.........
	stTranslationName.append( strId.str() );
	TGeoTranslation* pMatrixTrans = new TGeoTranslation( stTranslationName.c_str(), xc, yc, zc );
	//ALL clsses deriving from TGeoMatrix are not owned by the ROOT geometry manager, invoking RegisterYourself() transfers
	//ownership and thus ROOT will clean up
	pMatrixTrans->RegisterYourself();      

	//Create TGeoRotation object. 
	//Translations are of course just positional changes in the global frame.
	//Note that each subsequent rotation is using the new coordinate system of the last transformation all the way back to the global frame.
	//The way to think about this is that each rotation is the multiplication of the last rotation matrix by a new one.
	//The order is:
	//Integer Z rotation and reflections.
	//Z rotations specified by in degrees.
	//X rotations 
	//Y rotations
	TGeoRotation* pMatrixRotRefCombined = new TGeoRotation();
	//We have to ensure that we retain a right handed coordinate system, i.e. if we only flip the x or y axis, we have to also flip the z-axis. If we flip both we have to flip twice.	
	double integerRotationsAndReflections[9]={rotRef1,rotRef2,0,rotRef3,rotRef4,0,0,0, determinant};
	pMatrixRotRefCombined->SetMatrix(integerRotationsAndReflections);
	std::cout << "Rotating plane " << SensorId << " to gamma: " << gamma << std::endl;
	pMatrixRotRefCombined->RotateZ(gamma);//Z Rotation (degrees)//This will again rotate a vector around z axis usign the right hand rule.  
	pMatrixRotRefCombined->RotateX(alpha);//X Rotations (degrees)//This will rotate a vector usign the right hand rule round the x-axis
	pMatrixRotRefCombined->RotateY(beta);//Y Rotations (degrees)//Same again for Y axis
	pMatrixRotRefCombined->RegisterYourself();//We must allow the matrix to be used by the TGeo manager.
	// Combined translation and orientation
	TGeoCombiTrans* combi = new TGeoCombiTrans( *pMatrixTrans, *pMatrixRotRefCombined );
	//This is to print to screen the rotation and translation matrices used to transform from local to global frame.
	streamlog_out(MESSAGE9) << "THESE MATRICES ARE USED TO TAKE A POINT IN THE LOCAL FRAME AND MOVE IT TO THE GLOBAL FRAME."  << std::endl;   
	streamlog_out(MESSAGE9) << "SensorID: " << SensorId << " Rotation/Reflection matrix for this object."  << std::endl;   
	const double* rotationMatrix =  combi->GetRotationMatrix();	
	streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[0]<<"  "<<rotationMatrix[1]<<"   "<<rotationMatrix[2]<< std::endl;
	streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[3]<<"  "<<rotationMatrix[4]<<"   "<<rotationMatrix[5]<< std::endl;
	streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[6]<<"  "<<rotationMatrix[7]<<"   "<<rotationMatrix[8]<< std::endl;

	//streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[0] << std::setw(10) <<rotationMatrix[1]<< std::setw(10) <<rotationMatrix[2]<< std::setw(10)<< std::endl<< std::endl; 
	//streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[3] << std::setw(10) <<rotationMatrix[4]<< std::setw(10) <<rotationMatrix[5]<< std::setw(10)<< std::endl<< std::endl; 
	//streamlog_out(MESSAGE9) << std::setw(10) <<rotationMatrix[6] << std::setw(10) <<rotationMatrix[7]<< std::setw(10) <<rotationMatrix[8]<< std::setw(10)<< std::endl<< std::endl; 
	const double* translationMatrix =  combi->GetTranslation();	
	streamlog_out(MESSAGE9) << "SensorID: " << SensorId << " Translation vector for this object."  << std::endl;   
	streamlog_out(MESSAGE9) << std::setw(10) <<translationMatrix[0] << std::setw(10) <<translationMatrix[1]<< std::setw(10) <<translationMatrix[2]<< std::setw(10)<< std::endl; 

	combi->RegisterYourself();   
	
	// Construct object medium. Required for radiation length determination
	// assume SILICON, though all information except of radiation length is ignored
	double a       = 28.085500;     
	double z       = 14.000000;
	double density = 2.330000;
	double radl    = siPlaneRadLength( SensorId );
	double absl    = 45.753206;
	std::string stMatName = "materialSensor";
	stMatName.append( strId.str() );
	TGeoMaterial* pMat = new TGeoMaterial( stMatName.c_str(), a, z, density, -radl, absl );
	pMat->SetIndex( 1 );
	// Medium: medium_Sensor_SILICON
	int numed   = 0;  // medium number
	double par[8];
	par[0]  = 0.000000; // isvol
	par[1]  = 0.000000; // ifield
	par[2]  = 0.000000; // fieldm
	par[3]  = 0.000000; // tmaxfd
	par[4]  = 0.000000; // stemax
	par[5]  = 0.000000; // deemax
	par[6]  = 0.000000; // epsil
	par[7]  = 0.000000; // stmin
	std::string stMedName = "mediumSensor";
	stMedName.append( strId.str() );
	TGeoMedium* pMed = new TGeoMedium( stMedName.c_str(), numed, pMat, par );

	// Construct object shape
	// Shape: Box type: TGeoBBox
	// TGeo requires half-width of box side
	Double_t dx = siPlaneXSize( SensorId ) / 2.;
	Double_t dy = siPlaneYSize( SensorId ) / 2.;
	Double_t dz = siPlaneZSize( SensorId ) / 2.;
	TGeoShape *pBoxSensor = new TGeoBBox( "BoxSensor", dx, dy, dz );

	std::cout << "Box for sensor: " << SensorId << " is: " << dx << "|" << dy  << "|" << dz << '\n';

	// Geometry navigation package requires following names for objects that have an ID  name:ID
	std::string stVolName = "volume_SensorID:";
	stVolName.append( strId.str() );

	_planePath.insert( std::make_pair(SensorId, "/volume_World_1/"+stVolName+"_1") );

	TGeoVolume* pvolumeSensor = new TGeoVolume( stVolName.c_str(), pBoxSensor, pMed );
	pvolumeSensor->SetVisLeaves( kTRUE );
	pvolumeWorld->AddNode(pvolumeSensor, 1/*(SensorId)*/, combi);

	//this line tells the pixel geometry manager to load the pixel geometry into the plane			
	streamlog_out(DEBUG1) << " sensorID: " << SensorId << " " << stVolName << std::endl;   
	std::string name = geoLibName(SensorId);

	if( name == "CAST" ) {
		_pixGeoMgr->addCastedPlane( SensorId, siPlaneXNpixels(SensorId), siPlaneYNpixels(SensorId), siPlaneXSize(SensorId), siPlaneYSize(SensorId), siPlaneZSize(SensorId), siPlaneRadLength(SensorId), stVolName);
	} else {
		_pixGeoMgr->addPlane( SensorId, name, stVolName);
		updatePlaneInfo(SensorId);
	}
}

/**
 * Initialise ROOT geometry objects from GEAR objects
 * 
 * @param geomName name of ROOT geometry object
 * @param dumpRoot dump automatically generated ROOT geometry file for further inspection
 */
void EUTelGeometryTelescopeGeoDescription::initializeTGeoDescription( std::string& geomName, bool dumpRoot = false ) {
//    #ifdef USE_TGEO
    // get access to ROOT's geometry manager
    
	if( _isGeoInitialized )
	{
		streamlog_out( WARNING3 ) << "EUTelGeometryTelescopeGeoDescription: Geometry already initialized, using old initialization" << std::endl;
		return;
	}
	else
	{
    		_geoManager = new TGeoManager("Telescope", "v0.1");
	}

	if( !_geoManager )
	{
		streamlog_out( ERROR3 ) << "Can't instantiate ROOT TGeoManager " << std::endl;
		return;
	}
   
    
    // Create top world volume containing telescope/DUT geometry
    
    
    // Create air mixture
    // see http://pdg.lbl.gov/2013/AtomicNuclearProperties/HTML_PAGES/104.html
    double air_density = 1.2e-3;         // g/cm^3
    double air_radlen  = 36.62;          // g/cm^2
    TGeoMixture* pMatAir = new TGeoMixture("AIR",3,air_density);
    pMatAir->DefineElement(0, 14.007, 7.,  0.755267 );     //Nitrogen
    pMatAir->DefineElement(1, 15.999, 8.,  0.231781 );     //Oxygen
    pMatAir->DefineElement(2, 39.948, 18., 0.012827 );     //Argon
    pMatAir->DefineElement(3, 12.011, 6.,  0.000124 );     //Carbon
    pMatAir->SetRadLen( air_radlen );
    // Medium: medium_World_AIR
    TGeoMedium* pMedAir = new TGeoMedium("medium_World_AIR", 3, pMatAir );

    // The World is the 10 x 10m x 10m box filled with air mixture
    Double_t dx,dy,dz;
    dx = 5000.000000; // [mm]
    dy = 5000.000000; // [mm]
    dz = 5000.000000; // [mm]
    TGeoShape *pBoxWorld = new TGeoBBox("Box_World", dx,dy,dz);
    // Volume: volume_World
    TGeoVolume* pvolumeWorld = new TGeoVolume("volume_World",pBoxWorld, pMedAir);
    pvolumeWorld->SetLineColor(4);
    pvolumeWorld->SetLineWidth(3);
    pvolumeWorld->SetVisLeaves(kTRUE);

   // Set top volume of geometry
   gGeoManager->SetTopVolume( pvolumeWorld );
   
 
   
   // Iterate over registered GEAR objects and construct their TGeo representation
   
   const Double_t PI = 3.141592653589793;
   const Double_t DEG = 180./PI;
   
   double xc, yc, zc;   // volume center position 
   double alpha, beta, gamma;
   
   IntVec::const_iterator itrPlaneId;
   for ( itrPlaneId = _sensorIDVec.begin(); itrPlaneId != _sensorIDVec.end(); ++itrPlaneId ) {
       
       std::stringstream strId;
       strId << *itrPlaneId;
       
       // Get sensor center position
       xc = siPlaneXPosition( *itrPlaneId );
       yc = siPlaneYPosition( *itrPlaneId );
       zc = siPlaneZPosition( *itrPlaneId );
       
       // Get sensor orientation
       alpha = siPlaneXRotation( *itrPlaneId ); // [rad]
       beta  = siPlaneYRotation( *itrPlaneId ); // [rad]
       gamma = siPlaneZRotation( *itrPlaneId ); // [rad]
       
       // Spatial translations of the sensor center
       string stTranslationName = "matrixTranslationSensor";
       stTranslationName.append( strId.str() );
       TGeoTranslation* pMatrixTrans = new TGeoTranslation( stTranslationName.c_str(), xc, yc, zc );
       //ALL clsses deriving from TGeoMatrix are not owned by the ROOT geometry manager, invoking RegisterYourself() transfers
       //ownership and thus ROOT will clean up
       pMatrixTrans->RegisterYourself();      
       
       // Spatial rotation around sensor center
       // TGeoRotation requires Euler angles in degrees
       string stRotationName = "matrixRotationSensorX";
       stRotationName.append( strId.str() );
       TGeoRotation* pMatrixRotX = new TGeoRotation( stRotationName.c_str(), 0.,  alpha*DEG, 0.);                // around X axis
       stRotationName = "matrixRotationSensorY";
       stRotationName.append( strId.str() );
       TGeoRotation* pMatrixRotY = new TGeoRotation( stRotationName.c_str(), 90., beta*DEG,  0.);                // around Y axis (combination of rotation around Z axis and new X axis)
//.........这里部分代码省略.........

void create_sfi_geo(const char* geoTag)
{

  //fGlobalTrans->SetTranslation(0.0,0.0,0.0);

  // -------   Load media from media file   -----------------------------------
  FairGeoLoader*    geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");
  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  TString geoPath = gSystem->Getenv("VMCWORKDIR");
  TString medFile = geoPath + "/geometry/media_r3b.geo";
  geoFace->setMediaFile(medFile);
  geoFace->readMedia();
  gGeoMan = gGeoManager;
  // --------------------------------------------------------------------------



  // -------   Geometry file name (output)   ----------------------------------
  TString geoFileName = geoPath + "/geometry/sfi_";
  geoFileName = geoFileName + geoTag + ".geo.root";
  // --------------------------------------------------------------------------



  // -----------------   Get and create the required media    -----------------
  FairGeoMedia*   geoMedia = geoFace->getMedia();
  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();

  FairGeoMedium* mAir      = geoMedia->getMedium("Air");
  if ( ! mAir ) Fatal("Main", "FairMedium Air not found");
  geoBuild->createMedium(mAir);
  TGeoMedium* pMed2 = gGeoMan->GetMedium("Air");
  if ( ! pMed2 ) Fatal("Main", "Medium Air not found");

  FairGeoMedium* mVac      = geoMedia->getMedium("vacuum");
  if ( ! mVac ) Fatal("Main", "FairMedium vacuum not found");
  geoBuild->createMedium(mVac);
  TGeoMedium* pMed1 = gGeoMan->GetMedium("vacuum");
  if ( ! pMed1 ) Fatal("Main", "Medium vacuum not found");
  
  FairGeoMedium* mGfi      = geoMedia->getMedium("plasticForGFI");
  if ( ! mGfi ) Fatal("Main", "FairMedium plasticForGFI not found");
  geoBuild->createMedium(mGfi);
  TGeoMedium* pMed35 = gGeoMan->GetMedium("plasticForGFI");
  if ( ! pMed35 ) Fatal("Main", "Medium plasticForGFI not found");

  FairGeoMedium* mAl      = geoMedia->getMedium("aluminium");
  if ( ! mAl ) Fatal("Main", "FairMedium aluminium not found");
  geoBuild->createMedium(mAl);
  TGeoMedium* pMed21 = gGeoMan->GetMedium("aluminium");
  if ( ! pMed21 ) Fatal("Main", "Medium aluminium not found");
  // --------------------------------------------------------------------------



  // --------------   Create geometry and top volume  -------------------------
  gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");
  gGeoMan->SetName("GFIgeom");
  TGeoVolume* top = new TGeoVolumeAssembly("TOP");
  gGeoMan->SetTopVolume(top);
  // --------------------------------------------------------------------------
  
  //LABPOS(GFI1,-73.274339,0.069976,513.649524)
  Float_t dx = -73.274339; //dE tracker, correction due to wrong angle
  Float_t dy = 0.069976;
  Float_t dz = 513.649524;
  
  TGeoRotation *pMatrix3 = new TGeoRotation();
  //pMatrix3->RotateY(-16.7);
  TGeoCombiTrans*
  pMatrix2 = new TGeoCombiTrans("", dx,dy,dz,pMatrix3);

  //LABPOS(GFI2,-147.135037,0.069976,729.680342)
  dx = -147.135037; //dE tracker, correction due to wrong angle
  dy = 0.069976;
  dz = 729.680342;

  TGeoRotation *pMatrix5 = new TGeoRotation();
  //pMatrix5->RotateY(-16.7);
  TGeoCombiTrans*
  pMatrix4 = new TGeoCombiTrans("", dx,dy,dz,pMatrix5);
  
  
  // World definition
  TGeoVolume* pWorld = gGeoManager->GetTopVolume();
  pWorld->SetVisLeaves(kTRUE);
  
  // SHAPES, VOLUMES AND GEOMETRICAL HIERARCHY 
  
  // Volume: GFILogWorld
  
  TGeoVolume*   pGFILogWorld = new TGeoVolumeAssembly("GFILogWorld");
  pGFILogWorld->SetVisLeaves(kTRUE);
  
  // Global positioning
  pWorld->AddNode( pGFILogWorld, 0, pMatrix2 );
  
  Float_t detector_size = 5.120000;
  Float_t fiber_thickness = 0.020000; 
  
//.........这里部分代码省略.........