C++ TGeoRotation类代码示例

// function to convert this transformation to a TGeo transformation
// mainly used for the benchmark comparisons with ROOT
TGeoMatrix * Transformation3D::ConvertToTGeoMatrix() const
{
  if( fIdentity ){
      return new TGeoIdentity();
  }
  if( fHasTranslation && ! fHasRotation ) {
      return new TGeoTranslation(fTranslation[0], fTranslation[1], fTranslation[2]);
  }
  if( fHasRotation && ! fHasTranslation ) {
      TGeoRotation * tmp = new TGeoRotation();
      tmp->SetMatrix( Rotation() );
      return tmp;
  }
  if( fHasTranslation && fHasRotation )
  {
      TGeoRotation * tmp = new TGeoRotation();
      tmp->SetMatrix( Rotation() );
      return  new TGeoCombiTrans(fTranslation[0], fTranslation[1],
                     fTranslation[2], tmp);
  }
  return 0;
}

void Draw3D(AliITSOnlineCalibrationSPDhandler *h){

  TGeoHMatrix m2t[240];
  for(Int_t imod=0; imod<240; imod++){
    int vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(imod);
    AliITSAlignMille2Module::SensVolMatrix(vid,&m2t[imod]);
  }

  delete gGeoManager;

  new TGeoManager("SPD","active");

  TGeoMaterial *vacuum = new TGeoMaterial("vacuum",0,0,0);
  TGeoMedium *none = new TGeoMedium("Vacuum",0,vacuum);
  TGeoVolume *top = gGeoManager->MakeBox("TOP",none,500,500,500);
  gGeoManager->SetTopVolume(top);

  TGeoVolume *ladder = gGeoManager->MakeBox("ladder",none,0.6375,0.001/2,3.48);

  Int_t nActive[2]={0,0};
  for(Int_t imod=0; imod<240; imod++){
    TGeoRotation *rot  = new TGeoRotation();
    rot->SetMatrix(m2t[imod].GetRotationMatrix());
    TGeoCombiTrans *matrix = new TGeoCombiTrans(m2t[imod].GetTranslation()[0],m2t[imod].GetTranslation()[1],m2t[imod].GetTranslation()[2],rot);
    if((40960-h->GetNrBad(imod))>0) {
      top->AddNode(ladder,imod,matrix);
      if(imod<80) nActive[0]++;
      else nActive[1]++;
    }
  }

  printf("  \n\n   Number of Active SPD modules (->Total)  : inner %i (80) outer %i (160) \n\n\n",nActive[0],nActive[1]);
  gGeoManager->CloseGeometry();
  top->Draw("ogl");
  gPad->GetView()->ShowAxis();

}

TGeoCombiTrans* GetGlobalPosition(TGeoCombiTrans *fRef)
{
  if (fLocalTrans == kTRUE ) {

    if ( ( fThetaX == 0 )  && ( fThetaY==0 )  && ( fThetaZ == 0 )
	 &&
	 ( fX == 0 ) && ( fY == 0 ) && ( fZ == 0 )
	 )  return fRef;


    // X axis
    Double_t xAxis[3] = { 1. , 0. , 0. };
    Double_t yAxis[3] = { 0. , 1. , 0. };
    Double_t zAxis[3] = { 0. , 0. , 1. };
    // Reference Rotation
    fRefRot = fRef->GetRotation();

    if (fRefRot) {
      Double_t mX[3] = {0.,0.,0.};
      Double_t mY[3] = {0.,0.,0.};
      Double_t mZ[3] = {0.,0.,0.};

      fRefRot->LocalToMasterVect(xAxis,mX);
      fRefRot->LocalToMasterVect(yAxis,mY);
      fRefRot->LocalToMasterVect(zAxis,mZ);

      Double_t a[4]={ mX[0],mX[1],mX[2], fThetaX };
      Double_t b[4]={ mY[0],mY[1],mY[2], fThetaY };
      Double_t c[4]={ mZ[0],mZ[1],mZ[2], fThetaZ };

      ROOT::Math::AxisAngle aX(a,a+4);
      ROOT::Math::AxisAngle aY(b,b+4);
      ROOT::Math::AxisAngle aZ(c,c+4);

      ROOT::Math::Rotation3D fMatX( aX );
      ROOT::Math::Rotation3D fMatY( aY );
      ROOT::Math::Rotation3D fMatZ( aZ );

      ROOT::Math::Rotation3D  fRotXYZ = (fMatZ * (fMatY * fMatX));

      //cout << fRotXYZ << endl;

      Double_t fRotable[9]={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
      fRotXYZ.GetComponents(
			    fRotable[0],fRotable[3],fRotable[6],
			    fRotable[1],fRotable[4],fRotable[7],
			    fRotable[2],fRotable[5],fRotable[8]
			    );
      TGeoRotation *pRot = new TGeoRotation();
      pRot->SetMatrix(fRotable);
      TGeoCombiTrans *pTmp = new TGeoCombiTrans(*fGlobalTrans,*pRot);

      // ne peut pas etre applique ici
      // il faut differencier trans et rot dans la multi.
      TGeoRotation rot_id;
      rot_id.SetAngles(0.0,0.0,0.0);

      TGeoCombiTrans c1;
      c1.SetRotation(rot_id);
      const Double_t *t = pTmp->GetTranslation();
      c1.SetTranslation(t[0],t[1],t[2]);

      TGeoCombiTrans c2;
      c2.SetRotation(rot_id);
      const Double_t *tt = fRefRot->GetTranslation();
      c2.SetTranslation(tt[0],tt[1],tt[2]);

      TGeoCombiTrans cc = c1 * c2 ;

      TGeoCombiTrans c3;
      c3.SetRotation(pTmp->GetRotation());
      TGeoCombiTrans c4;
      c4.SetRotation(fRefRot);

      TGeoCombiTrans ccc = c3 * c4;

      TGeoCombiTrans pGlobal;
      pGlobal.SetRotation(ccc.GetRotation());
      const Double_t *allt = cc.GetTranslation();
      pGlobal.SetTranslation(allt[0],allt[1],allt[2]);

      return  ( new TGeoCombiTrans( pGlobal ) );

    }else{

      cout << "-E- R3BDetector::GetGlobalPosition() \
	      No. Ref. Transformation defined ! " << endl;
      cout << "-E- R3BDetector::GetGlobalPosition() \
	      cannot create Local Transformation " << endl;
      return NULL;
    } //! fRefRot

  } else {
    // Lab Transf.
    if ( ( fPhi == 0 )  && ( fTheta==0 )  && ( fPsi == 0 )

//.........这里部分代码省略.........

  //------------------ wrapping Alu------------------------------------------
  TGeoShape* padle_h_box1 = new TGeoBBox("padle_h_box1",
					 neuLAND_paddle_dimx, 
					 neuLAND_paddle_dimy + neuLAND_wrapping1_dim, 
					 neuLAND_paddle_dimz + neuLAND_wrapping1_dim);
  TGeoShape* padle_h_box2 = new TGeoBBox("padle_h_box2",
					 neuLAND_paddle_dimx, 
					 neuLAND_paddle_dimy, 
					 neuLAND_paddle_dimz);
  
  // Create a composite shape
  TGeoCompositeShape *wrapping1 = new TGeoCompositeShape("diffbox", "padle_h_box1 - padle_h_box2");
  TGeoVolume *bvol1 = new TGeoVolume("wrapping1", wrapping1, pMedAl);

  //------------------ wrapping Tape------------------------------------------
  TGeoShape* padle_h_box3 = new TGeoBBox("padle_h_box3",
					 neuLAND_paddle_dimx,
					 neuLAND_paddle_dimy + neuLAND_wrapping1_dim + neuLAND_wrapping2_dim, 
					 neuLAND_paddle_dimz + neuLAND_wrapping1_dim + neuLAND_wrapping2_dim);
  TGeoShape* padle_h_box4 = new TGeoBBox("padle_h_box4",
					 neuLAND_paddle_dimx, 
					 neuLAND_paddle_dimy + neuLAND_wrapping1_dim, 
					 neuLAND_paddle_dimz + neuLAND_wrapping1_dim);

  // Create a composite shape
  TGeoCompositeShape *wrapping2 = new TGeoCompositeShape("diffbox", "padle_h_box3 - padle_h_box4");
  TGeoVolume *bvol2 = new TGeoVolume("wrapping2", wrapping2, pMed38);



  // Make the elementary assembly of the whole structure
  TGeoVolume *aLand = new TGeoVolumeAssembly("ALAND");

  Double_t total_dimx = neuLAND_paddle_dimx;
  Double_t total_dimy = neuLAND_paddle_dimy +
    neuLAND_wrapping1_dim + 
    neuLAND_wrapping2_dim + 
    neuLAND_gap_dim;
  Double_t total_dimz = neuLAND_paddle_dimz +
    neuLAND_wrapping1_dim +
    neuLAND_wrapping2_dim +
    neuLAND_gap_dim;
  
  //paddles
  TGeoRotation *zeroRotation = new TGeoRotation();
  zeroRotation->RotateX(0.);
  zeroRotation->RotateY(0.);
  zeroRotation->RotateZ(0.);

  TGeoRotation *rot1 = new TGeoRotation();
  rot1->RotateX(0.);
  rot1->RotateY(0.);
  rot1->RotateZ(90.);

  Double_t xx = 0.;
  Double_t yy = 0.;
  Double_t zz = 0.;
  
  aLand->AddNode(padle_h_box5, 1, new TGeoCombiTrans(xx, yy, zz, zeroRotation));
  aLand->AddNode(bvol1, 1, new TGeoCombiTrans(xx, yy, zz, zeroRotation));
  aLand->AddNode(bvol2, 1, new TGeoCombiTrans(xx, yy, zz, zeroRotation));

  TGeoVolume *cell = new TGeoVolumeAssembly("CELL");

  Int_t nindex = 0, i = 0;
  tx = 0.;
  tz = -neuLAND_depth_dim + total_dimz;
  for(tz = -neuLAND_depth_dim + total_dimz; tz < neuLAND_depth_dim; tz += total_dimz*2) {
    i += 1;
    for (ty = -total_dimx + total_dimy; ty < total_dimx; ty += total_dimy*2) {
      nindex++;
      if (i % 2 == 1) {
        cell->AddNode(aLand, nindex, new TGeoCombiTrans(tx, ty, tz, zeroRotation));
      } else {
        cell->AddNode(aLand, nindex, new TGeoCombiTrans(ty, tx, tz, rot1));
      }
    }
  }

  tx = 0.0;
  ty = 0.0;
  tz = fZ;
  TGeoCombiTrans *t0 = new TGeoCombiTrans(tx, ty, tz, zeroRotation);

  top->AddNode(cell, 1, /*GetGlobalPosition(t0)*/t0);



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

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

void create_tof_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/tof_";
    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* mTof      = geoMedia->getMedium("plasticForTOF");
    if ( ! mTof ) Fatal("Main", "FairMedium plasticForTOF not found");
    geoBuild->createMedium(mTof);
    TGeoMedium* pMed34 = gGeoMan->GetMedium("plasticForTOF");
    if ( ! pMed34 ) Fatal("Main", "Medium plasticForTOF not found");
    // --------------------------------------------------------------------------



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



    // out-of-file geometry definition
    Double_t dx,dy,dz;
    Double_t a;
    //   Double_t thx, phx, thy, phy, thz, phz;
    Double_t z, density, w;
    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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  Double_t alDy = 103.4/2.; // [cm]
  //Double_t alDz = 4.06;     // [cm]
  Double_t alDz = 7.50;     // [cm]	//modified

  TGeoShape* pAlBox = new TGeoBBox("AlBox",
				   alDx,
				   alDy,
				   alDz);

  TGeoVolume*
    pAlDchLog = new TGeoVolume("ALBoxLog",pAlBox, pMed21);
  pAlDchLog->SetVisLeaves(kTRUE);



  // Mylar Entrance exit windows
  Double_t mylDx= gasDx; //[cm]
  Double_t mylDy= gasDy; //[cm]
  Double_t mylDz= 0.0006; //[cm]
  TGeoShape* pMylarBox = new TGeoBBox("MylarBox",
				      mylDx,
				      mylDy,
				      mylDz);

  TGeoVolume*
    pMylDchLog = new TGeoVolume("MYLBoxLog",pMylarBox, pMed15);
  pMylDchLog->SetVisLeaves(kTRUE);



  // First assembly
  TGeoVolume *dch1 = new TGeoVolumeAssembly("DCH1");

  TGeoRotation *rot = new TGeoRotation();
  rot->RotateX(0.);
  rot->RotateY(0.);
  rot->RotateZ(0.);
  //Double_t tx = -3.5;	//correct values? active-area vs Al frame offsets?
  //Double_t ty = -5.;
  Double_t tx = +3.5;	//corrected values, active-area vs Al frame offsets. F.Wamers.
  Double_t ty = -3.5;
  Double_t tz = 0.;

  TGeoCombiTrans*
    pTransfo1 = new TGeoCombiTrans("", 0.,0.,0.,rot);
  TGeoCombiTrans*
    pTransfo2 = new TGeoCombiTrans("", tx,ty,tz,rot);
  TGeoCombiTrans*
    pTransfo3 = new TGeoCombiTrans("", -tx,-ty,-tz,rot); //new, in order to compensate Al vs Gas

  dch1->AddNode(pAlDchLog,0,pTransfo3);	//1->3
  pAlDchLog->AddNode(pGasDchLog,0,pTransfo2);	//place gas in aluminum, shifted
  pGasDchLog->AddNode(pActGasDchLog, 0,pTransfo1);	//place active in gas, centrally
  // Mylar Windows front+back
  //dch1->AddNode(pMylDchLog,0,new TGeoCombiTrans("", tx,ty,-alDz-mylDz,rot));
  //dch1->AddNode(pMylDchLog,1,new TGeoCombiTrans("", tx,ty, alDz+mylDz,rot));
  dch1->AddNode(pMylDchLog,0,new TGeoCombiTrans("", 0.,0.,-alDz-mylDz,rot));
  dch1->AddNode(pMylDchLog,1,new TGeoCombiTrans("", 0.,0., alDz+mylDz,rot));



  // Global Positioning
  //in agreement with the s318 tracker, those are supposed to be the centres of the active volumes!!!
  //Double_t pDch1x = -123.22 ; //Justyna 
  //Double_t pDch1y = 3.6 ;     //Justyna
  //Double_t pDch1z = 444.13 ;  //Justyna

void assembly()
{
//--- Definition of a simple geometry
    gSystem->Load("libGeom");
    TGeoManager *geom = new TGeoManager("Assemblies",
                                        "Geometry using assemblies");
    Int_t i;
    //--- define some materials
    TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
//   //--- define some media
    TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
    TGeoMedium *Al = new TGeoMedium("Aluminium",2, matAl);

    //--- make the top container volume
    TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 1000., 1000., 100.);
    geom->SetTopVolume(top);

    // Make the elementary assembly of the whole structure
    TGeoVolume *tplate = new TGeoVolumeAssembly("TOOTHPLATE");

    Int_t ntooth = 5;
    Double_t xplate = 25;
    Double_t yplate = 50;
    Double_t xtooth = 10;
    Double_t ytooth = 0.5*yplate/ntooth;
    Double_t dshift = 2.*xplate + xtooth;
    Double_t xt,yt;

    TGeoVolume *plate = geom->MakeBox("PLATE", Al, xplate,yplate,1);
    plate->SetLineColor(kBlue);
    TGeoVolume *tooth = geom->MakeBox("TOOTH", Al, xtooth,ytooth,1);
    tooth->SetLineColor(kBlue);
    tplate->AddNode(plate,1);
    for (i=0; i<ntooth; i++) {
        xt = xplate+xtooth;
        yt = -yplate + (4*i+1)*ytooth;
        tplate->AddNode(tooth, i+1, new TGeoTranslation(xt,yt,0));
        xt = -xplate-xtooth;
        yt = -yplate + (4*i+3)*ytooth;
        tplate->AddNode(tooth, ntooth+i+1, new TGeoTranslation(xt,yt,0));
    }

    TGeoRotation *rot1 = new TGeoRotation();
    rot1->RotateX(90);
    TGeoRotation *rot;
    // Make a hexagone cell out of 6 tooth plates. These can zip together
    // without generating overlaps (they are self-contained)
    TGeoVolume *cell = new TGeoVolumeAssembly("CELL");
    for (i=0; i<6; i++) {
        Double_t phi =  60.*i;
        Double_t phirad = phi*TMath::DegToRad();
        Double_t xp = dshift*TMath::Sin(phirad);
        Double_t yp = -dshift*TMath::Cos(phirad);
        rot = new TGeoRotation(*rot1);
        rot->RotateZ(phi);
        cell->AddNode(tplate,i+1,new TGeoCombiTrans(xp,yp,0,rot));
    }

    // Make a row as an assembly of cells, then combine rows in a honeycomb
    // structure. This again works without any need to define rows as
    // "overlapping"
    TGeoVolume *row = new TGeoVolumeAssembly("ROW");
    Int_t ncells = 5;
    for (i=0; i<ncells; i++) {
        Double_t ycell = (2*i+1)*(dshift+10);
        row->AddNode(cell, ncells+i+1, new TGeoTranslation(0,ycell,0));
        row->AddNode(cell,ncells-i,new TGeoTranslation(0,-ycell,0));
    }

    Double_t dxrow = 3.*(dshift+10.)*TMath::Tan(30.*TMath::DegToRad());
    Double_t dyrow = dshift+10.;
    Int_t nrows = 5;
    for (i=0; i<nrows; i++) {
        Double_t xrow = 0.5*(2*i+1)*dxrow;
        Double_t yrow = 0.5*dyrow;
        if ((i%2)==0) yrow = -yrow;
        top->AddNode(row, nrows+i+1, new TGeoTranslation(xrow,yrow,0));
        top->AddNode(row, nrows-i, new TGeoTranslation(-xrow,-yrow,0));
    }

    //--- close the geometry
    geom->CloseGeometry();

    geom->SetVisLevel(4);
    geom->SetVisOption(0);
    top->Draw();
}

void AddMirrors(AOpticalComponent* opt)
{
  // dummy hexagonal prism to cut a spherical mirror
  TGeoPgon* mirCut = new TGeoPgon("mirCut", 0., 360., 6, 2);
  mirCut->DefineSection(0, -100*mm, 0, kMirrorD/2.);
  mirCut->DefineSection(1,  100*mm, 0, kMirrorD/2.);

  double theta = TMath::ASin(kMirrorD/TMath::Sqrt(3)/kMirrorR)*TMath::RadToDeg();
  TGeoSphere* mirSphere = new TGeoSphere("mirSphere", kMirrorR, kMirrorR + kMirrorT, 180. - theta, 180.);
  TGeoTranslation* transZ = new TGeoTranslation("transZ", 0, 0, kMirrorR);
  transZ->RegisterYourself();
  TGeoCompositeShape* mirComposite = new TGeoCompositeShape("mirComposite", "mirSphere:transZ*mirCut");
  AMirror* mirror = new AMirror("mirror", mirComposite);

  const int kNMirror = 88;
  double dx = kMirrorD/TMath::Sqrt(3);
  double dy = kMirrorD/2.;
  double x[kNMirror] = {0, 0, 0, 0, 0, 0, 0, 0,
                        1.5*dx, 1.5*dx, 1.5*dx, 1.5*dx, 1.5*dx,
                        1.5*dx, 1.5*dx, 1.5*dx, 1.5*dx, 1.5*dx,
                        -1.5*dx, -1.5*dx, -1.5*dx, -1.5*dx, -1.5*dx,
                        -1.5*dx, -1.5*dx, -1.5*dx, -1.5*dx, -1.5*dx,
                        3*dx, 3*dx, 3*dx, 3*dx, 3*dx,
                        3*dx, 3*dx, 3*dx, 3*dx,
                        -3*dx, -3*dx, -3*dx, -3*dx, -3*dx,
                        -3*dx, -3*dx, -3*dx, -3*dx,
                        4.5*dx, 4.5*dx, 4.5*dx, 4.5*dx,
                        4.5*dx, 4.5*dx, 4.5*dx, 4.5*dx,
                        -4.5*dx, -4.5*dx, -4.5*dx, -4.5*dx,
                        -4.5*dx, -4.5*dx, -4.5*dx, -4.5*dx,
                        6*dx, 6*dx, 6*dx, 6*dx,
                        6*dx, 6*dx, 6*dx,
                        -6*dx, -6*dx, -6*dx, -6*dx,
                        -6*dx, -6*dx, -6*dx,
                        7.5*dx, 7.5*dx, 7.5*dx,
                        7.5*dx, 7.5*dx, 7.5*dx,
                        -7.5*dx, -7.5*dx, -7.5*dx,
                        -7.5*dx, -7.5*dx, -7.5*dx};

  double y[kNMirror] = {2*dy, 4*dy, 6*dy, 8*dy, -2*dy, -4*dy, -6*dy, -8*dy,
                        1*dy, 3*dy, 5*dy, 7*dy, 9*dy,
                        -1*dy, -3*dy, -5*dy, -7*dy, -9*dy,
                        1*dy, 3*dy, 5*dy, 7*dy, 9*dy,
                        -1*dy, -3*dy, -5*dy, -7*dy, -9*dy,
                        0*dy, 2*dy, 4*dy, 6*dy, 8*dy,
                        -2*dy, -4*dy, -6*dy, -8*dy,
                        0*dy, 2*dy, 4*dy, 6*dy, 8*dy,
                        -2*dy, -4*dy, -6*dy, -8*dy,
                        1*dy, 3*dy, 5*dy, 7*dy,
                        -1*dy, -3*dy, -5*dy, -7*dy,
                        1*dy, 3*dy, 5*dy, 7*dy,
                        -1*dy, -3*dy, -5*dy, -7*dy,
                        0*dy, 2*dy, 4*dy, 6*dy,
                        -2*dy, -4*dy, -6*dy,
                        0*dy, 2*dy, 4*dy, 6*dy,
                        -2*dy, -4*dy, -6*dy,
                        1*dy, 3*dy, 5*dy,
                        -1*dy, -3*dy, -5*dy,
                        1*dy, 3*dy, 5*dy,
                        -1*dy, -3*dy, -5*dy};

  for(int i = 0; i < kNMirror; i++){
    double r2d = TMath::RadToDeg();
    double r2 = TMath::Power(x[i], 2) + TMath::Power(y[i], 2);
    double z = kF - TMath::Sqrt(TMath::Power(kF, 2) - r2);

    // each mirror center is relocated from the origin (0, 0, 0) to (x, y, z)
    TGeoTranslation* trans = new TGeoTranslation(Form("mirTrans%d", i), x[i], y[i], z);

    // and is rotated to compose a DC optics
    double phi = TMath::ATan2(y[i], x[i])*r2d;
    TGeoRotation* rot = new TGeoRotation(Form("mirRot%d", i), - phi + 90., 0, 0);
    theta = TMath::ATan2(TMath::Sqrt(r2), 2*kF - z)*r2d;
    TGeoRotation* rot2 = new TGeoRotation("", phi - 90., theta, 0);
    rot->MultiplyBy(rot2, 0);
    
    // make a matrix from translation and rotation matrices
    TGeoCombiTrans* combi = new TGeoCombiTrans(*trans, *rot);

    // finally add this mirror to the world
    opt->AddNode(mirror, i + 1, combi);
  } // i
}

void EUTelGeometryTelescopeGeoDescription::translateSiPlane2TGeo(TGeoVolume* pvolumeWorld, int SensorId ) {
	double xc, yc, zc;   // volume center position 
	double alpha, beta, gamma;
	double rotRef1, rotRef2, rotRef3, rotRef4;

	std::stringstream strId;
	strId << SensorId;

	// Get sensor center position
	xc = siPlaneXPosition( SensorId );
	yc = siPlaneYPosition( SensorId );
	zc = siPlaneZPosition( SensorId );

	// Get sensor orientation
	alpha = siPlaneXRotation( SensorId ); //  in degrees !
	beta  = siPlaneYRotation( SensorId ); // 
	gamma = siPlaneZRotation( SensorId ); // 

	rotRef1 = siPlaneRotation1( SensorId );
	rotRef2 = siPlaneRotation2( SensorId );
	rotRef3 = siPlaneRotation3( SensorId );
	rotRef4 = siPlaneRotation4( SensorId );

	//We must check that the input is correct. Since this is a combination of initial rotations and reflections the determinate must be 1 or -1
	float determinant = rotRef1*rotRef4 - rotRef2*rotRef3  ;
	if(determinant==1 or determinant==-1) { 
		streamlog_out(DEBUG5) << "SensorID: " << SensorId << ". Determinant =  " <<determinant <<"  This is the correct determinate for this transformation." << std::endl;   
	} else {
		streamlog_out(ERROR5) << "SensorID: " << SensorId << ". Determinant =  " <<determinant << std::endl;   
		throw(lcio::Exception("The initial rotation and reflection matrix does not have determinant of 1 or -1. Gear file input must be wrong.")); 	
	}
	//Create spatial TGeoTranslation object.
	std::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();      

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

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