//____________________________________________________________________________
void get_target_mass ( 
  string geometry_file, string topvolname = "", 
  double length_unit=units::cm, double density_unit=units::g_cm3, 
  Bool_t checkoverlaps = kFALSE)
{
  PDGLibrary* pdglib = PDGLibrary::Instance(); // get message out of the way

   gFileName = geometry_file;
   // import geometry
   TGeoManager *tgeo = new TGeoManager("TGeo","TGeo");
   tgeo->Import( geometry_file.c_str() );

   // set top volume 
   if( ! set_top_volume(topvolname) ) return;

   // draw
   gGeoManager->GetTopVolume()->Draw();

   // method fails if there are volume overlaps, check it if the user wants to... (takes some time)
   if(checkoverlaps) {
      Double_t overlap_acc = 0.05;
      gGeoManager->CheckOverlaps( overlap_acc );
   }

   get_mass(length_unit, density_unit);
}

void OpNoviceGeom()
{
// Create root geometry corresponding to OpNovice example
  double fExpHall_x, fExpHall_y, fExpHall_z;
  double fTank_x, fTank_y, fTank_z;
  double fBubble_x, fBubble_y, fBubble_z;
  double a, z, density;
  int nelements;

  fExpHall_x = fExpHall_y = fExpHall_z = 10000.0;
  fTank_x    = fTank_y    = fTank_z    =  5000.0;
  fBubble_x  = fBubble_y  = fBubble_z  =  500.0;
  TGeoManager *geom = new TGeoManager("OpNoviceGeom", "GEANT4 OpNovice example");
// Air
//
  TGeoElement *N = new TGeoElement("Nitrogen", "N", z=7 , a=14.01);
  TGeoElement *O = new TGeoElement("Oxygen"  , "O", z=8 , a=16.00);
  TGeoMixture *air = new TGeoMixture("Air", 2, density=1.29E-3);
  air->AddElement(N, 0.7);
  air->AddElement(O, 0.3);
  TGeoMedium *medair = new TGeoMedium("Air", 1, air);
// Water
//
  TGeoElement *H = new TGeoElement("Hydrogen", "H", z=1 , a=1.01);
  TGeoMixture *water = new TGeoMixture("Water", 2, density=1.0); 
  water->AddElement(H,2);
  water->AddElement(O,1);
  TGeoMedium *medwater = new TGeoMedium("Water", 2, water);
// The experimental Hall
//
  TGeoBBox *expHall_box = new TGeoBBox("World",fExpHall_x,fExpHall_y,fExpHall_z);
  TGeoVolume *expHall_log = new TGeoVolume("World", expHall_box, medair);
  expHall_log->SetLineColor(1);
  expHall_log->SetVisContainers();
  geom->SetTopVolume(expHall_log);
// The Water Tank
//
  TGeoBBox *waterTank_box = new TGeoBBox("Tank",fTank_x,fTank_y,fTank_z);
  TGeoVolume *waterTank_log = new TGeoVolume("Tank", waterTank_box, medwater);
  waterTank_log->SetLineColor(kBlue);
  waterTank_log->SetTransparency(70);
  waterTank_log->SetVisContainers();
  expHall_log->AddNode(waterTank_log, 0);
// The Air Bubble
//
  TGeoBBox *bubbleAir_box = new TGeoBBox("Bubble",fBubble_x,fBubble_y,fBubble_z);
  TGeoVolume *bubbleAir_log = new TGeoVolume("Bubble", bubbleAir_box, medair);
  bubbleAir_log->SetLineColor(kCyan);
  bubbleAir_log->SetTransparency(70);
  waterTank_log->AddNode(bubbleAir_log, 0, new TGeoTranslation(0,2500,0));
  
  geom->CloseGeometry();
  geom->SetTopVisible(true);
  geom->Export("OpNoviceGeom.root");
}

void LCDDImp::init() {
  if (!m_world.isValid()) {
    TGeoManager* mgr = m_manager;
    Box worldSolid("world_x", "world_y", "world_z");
    printout(INFO,"LCDD"," *********** created World volume with size : %7.0f %7.0f %7.0f",
             worldSolid->GetDX(), worldSolid->GetDY(), worldSolid->GetDZ());

    m_materialAir = material("Air");
    m_materialVacuum = material("Vacuum");

    Volume world_vol("world_volume", worldSolid, m_materialAir);
    m_world = DetElement(new WorldObject(*this,"world"));
    m_worldVol = world_vol;
    // Set the world volume to invisible.
    VisAttr worldVis("WorldVis");
    worldVis.setAlpha(1.0);
    worldVis.setVisible(false);
    worldVis.setShowDaughters(true);
    worldVis.setColor(1.0, 1.0, 1.0);
    worldVis.setLineStyle(VisAttr::SOLID);
    worldVis.setDrawingStyle(VisAttr::WIREFRAME);
    //m_worldVol.setVisAttributes(worldVis);
    m_worldVol->SetVisibility(kFALSE);
    m_worldVol->SetVisDaughters(kTRUE);
    m_worldVol->SetVisContainers(kTRUE);
    add(worldVis);

#if 0
    Tube trackingSolid(0.,"tracking_region_radius","2*tracking_region_zmax",2*M_PI);
    Volume tracking("tracking_volume",trackingSolid, m_materialAir);
    m_trackers = TopDetElement("tracking");
    m_trackingVol = tracking;
    PlacedVolume pv = m_worldVol.placeVolume(tracking);
    m_trackers.setPlacement(pv);
    m_world.add(m_trackers);
#endif
    m_detectors.append(m_world);
    m_manager->SetTopVolume(m_worldVol.ptr());
    m_world.setPlacement(PlacedVolume(mgr->GetTopNode()));
  }
}

void LCDDImp::endDocument() {
  TGeoManager* mgr = m_manager;
  if (!mgr->IsClosed()) {
#if 0
    Region trackingRegion("TrackingRegion");
    trackingRegion.setThreshold(1);
    trackingRegion.setStoreSecondaries(true);
    add(trackingRegion);
    m_trackingVol.setRegion(trackingRegion);
    // Set the tracking volume to invisible.
    VisAttr trackingVis("TrackingVis");
    trackingVis.setVisible(false);
    m_trackingVol.setVisAttributes(trackingVis);
    add(trackingVis);
#endif
    /// Since we allow now for anonymous shapes,
    /// we will rename them to use the name of the volume they are assigned to
    mgr->CloseGeometry();
    ShapePatcher patcher(m_volManager, m_world);
    patcher.patchShapes();
    mapDetectorTypes();
  }
}

void AddLadder()
{
  SvxTGeo *geo = new SvxTGeo;
  geo->ReadParFile("parfiles/svxPISA.par");

  // Make an empty 100x100x100 cm^3 space in the experiment hall.
  geo->MakeTopVolume(200/2, 200/2, 200/2);

  // Place VTX sensors in the volume.
  geo->AddSensors();

  // Get handles for further manipulation
  TGeoManager *mgr = geo->GeoManager();
  TGeoVolume *top = mgr->GetTopVolume();

  // Add a ladder with new index 20 and the same geometry as B1L9.
  double xyz[3] = {0};
  geo->GetSensorXYZ(1, 9, 0, xyz);
  geo->AddLadder(1, 20, xyz[0], xyz[1], 0., geo->GetLadderPhiTilt(1,9));
  
  // Now rotate B1L20 by angular difference between B1L9 and B1L8.
  double dphi = geo->SensorPhiRad(1,9,0) - geo->SensorPhiRad(1,8,0);
  geo->RotateLadder(1, 20, 0., 0., dphi);

  // Done building model.
  // Close geometry to check for problems (overlapping boundaries)
  mgr->CloseGeometry();

  // Press j,k to zoom; u,i to look up/down; h,l to look left, right.
  TCanvas *c = new TCanvas("c", "svx model", 1400, 1000);
  c->SetFillColor(kBlack);
  top->Draw();

  geo->WriteParFile("parfiles/svxPISA.newladders.par");
  
  return;
}

void event3D()
{
  gSystem->Load("libGeom");
   TGeoManager *geom = new TGeoManager("simple1", "Simple geometry");

   //--- 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("Root Material",2, matAl);

     //--- make the top container volume
   Double_t worldx = 110.;
   Double_t worldy = 50.;
   Double_t worldz = 5.;

   TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 100., 100., 100.);
   geom->SetTopVolume(top);

   gGeoManager->GetVolume("TOP")->InvisibleAll();

   //Measurements (in cm)
   
   Double_t target_rad = 0.3;
   Double_t target_len = 48.58;

   Double_t RTPC_len = 41.;
   
   Double_t foil1_R1 = 2.;
   Double_t foil1_R2 = 2.0000018;

   Double_t foil2_R1 = 3.;
   Double_t foil2_R2 = 3.0000018;

   Double_t GEM1_R1 = 7.;
   Double_t GEM1_R2 = 7.0005;

   Double_t GEM2_R1 = 7.3;
   Double_t GEM2_R2 = 7.3005;

   Double_t GEM3_R1 = 7.6;
   Double_t GEM3_R2 = 7.6005;

   Double_t Z_offset = 0;

   Int_t trans_lvl = 70;
   
   TGeoTube *target = new TGeoTube("target_tube",0, target_rad, target_len/2);
   
   TGeoVolume *target_vol = new TGeoVolume("target_vol",target, Al); //(*)

   top->AddNode(target_vol,1,new TGeoTranslation(0,0, Z_offset));
   gGeoManager->GetVolume("target_vol")->SetTransparency(0);   

   TGeoTube *foil1 = new TGeoTube("foil1_tube", foil1_R1, foil1_R2, RTPC_len/2);
   
   TGeoVolume *foil1_vol = new TGeoVolume("foil1_vol", foil1, Al); //(*)

   top->AddNode(foil1_vol,1,new TGeoTranslation(0,0, Z_offset));


   TGeoTube *foil2 = new TGeoTube("foil2_tube", foil2_R1, foil2_R2, RTPC_len/2);
   
   TGeoVolume *foil2_vol = new TGeoVolume("foil2_vol", foil2, Al); //(*)

   top->AddNode(foil2_vol,1,new TGeoTranslation(0,0, Z_offset));


   TGeoTube *gem1 = new TGeoTube("gem1_tube", GEM1_R1, GEM1_R2, RTPC_len/2);
   TGeoVolume *gem1_vol = new TGeoVolume("gem1_vol", gem1, Al); //(*)
   gem1_vol->SetLineColor(kOrange);
   top->AddNode(gem1_vol,1,new TGeoTranslation(0,0, Z_offset));

   TGeoTube *gem2 = new TGeoTube("gem2_tube", GEM2_R1, GEM2_R2, RTPC_len/2);
   TGeoVolume *gem2_vol = new TGeoVolume("gem2_vol", gem2, Al); //(*)
   gem2_vol->SetLineColor(kGreen);
   top->AddNode(gem2_vol,1,new TGeoTranslation(0,0, Z_offset));

   TGeoTube *gem3 = new TGeoTube("gem3_tube", GEM3_R1, GEM3_R2, RTPC_len/2);
   TGeoVolume *gem3_vol = new TGeoVolume("gem3_vol", gem3, Al); //(*)
   top->AddNode(gem3_vol,1,new TGeoTranslation(0,0, Z_offset));
   gem3_vol->SetLineColor(kRed);
   gGeoManager->GetVolume("gem3_vol")->SetTransparency(trans_lvl);   

   
   //--- draw the ROOT box.
   // by default the picture will appear in the standard ROOT TPad.
   //if you have activated the following line in system.rootrc,
   //it will appear in the GL viewer
   //#Viewer3D.DefaultDrawOption:   ogl
   
   geom->CloseGeometry();
   
   geom->SetVisLevel(4);
  
   //   top->SetVisibility(kFALSE);
   top->Draw("ogl");
   gWorld = top;
}

void lego()
{
  // Drawing a figure, made of lego block, using ROOT geometry class.
  // Name: lego.C
  // Author: Soon Gi Kwon([email protected]), Dept. of Physics, Univ. of Seoul
  // Reviewed by Sunman Kim ([email protected])
  // Supervisor: Prof. Inkyu Park ([email protected])
  //
  // How to run: .x lego.C in ROOT terminal, then use OpenGL
  //
  // This macro was created for the evaluation of Computational Physics course in 2006.
  // We thank to Prof. Inkyu Park for his special lecture on ROOT and to all of ROOT team
  //

   TGeoManager *geom = new TGeoManager("geom","My first 3D geometry");


   TGeoMaterial *vacuum=new TGeoMaterial("vacuum",0,0,0);
   TGeoMaterial *Fe=new TGeoMaterial("Fe",55.845,26,7.87);



   TGeoMedium *Air=new TGeoMedium("Vacuum",0,vacuum);
   TGeoMedium *Iron=new TGeoMedium("Iron",1,Fe);


 // create volume
   TGeoVolume *top=geom->MakeBox("top",Air,100,100,100);
   geom->SetTopVolume(top);
   geom->SetTopVisible(0);
   // If you want to see the boundary, please input the number, 1 instead of 0.
   // Like this, geom->SetTopVisible(1);


//----------------------------------------------------------------------

TGeoVolume *ha1=geom->MakeSphere("ha1",Iron,0,10,80,90,0,360);
 ha1->SetLineColor(41);
top->AddNodeOverlap(ha1,1,new TGeoCombiTrans(0,0,4,new TGeoRotation("ha1",0,0,0)));

TGeoVolume *ha2=geom->MakeSphere("ha2",Iron,0,7,90,180,0,360);
 ha2->SetLineColor(41);
top->AddNodeOverlap(ha2,1,new TGeoCombiTrans(0,0,4,new TGeoRotation("ha2",0,180,0)));

TGeoVolume *ha3=geom->MakeSphere("ha3",Iron,0,7.3,80,90,0,360);
 ha3->SetLineColor(2);
top->AddNodeOverlap(ha3,1,new TGeoCombiTrans(0,0,4.8,new TGeoRotation("ha3",0,0,0)));


TGeoVolume *h1=geom->MakeTubs("h1",Iron,0,6,4.5,0,0);
 h1->SetLineColor(5);
top->AddNodeOverlap(h1,1,new TGeoCombiTrans(0,0,0,new TGeoRotation("h1",0,0,0)));

TGeoVolume *h2=geom->MakeSphere("h2",Iron,0,7.5,0,52.5,0,360);
 h2->SetLineColor(5);
top->AddNodeOverlap(h2,1,new TGeoCombiTrans(0,0,0,new TGeoRotation("h2",0,0,0)));

TGeoVolume *h3=geom->MakeSphere("h3",Iron,0,7.5,0,52.5,0,360);
 h3->SetLineColor(5);
top->AddNodeOverlap(h3,1,new TGeoCombiTrans(0,0,0,new TGeoRotation("h3",180,180,0)));

TGeoVolume *h4=geom->MakeTubs("h4",Iron,2.5,3.5,1.5,0,0);
 h4->SetLineColor(5);
top->AddNodeOverlap(h4,1,new TGeoCombiTrans(0,0,7.5,new TGeoRotation("h4",0,0,0)));



TGeoVolume *t1_1=geom->MakeTubs("t1_1",Iron,0,0.8,1,0,360);
 t1_1->SetLineColor(12);
top->AddNodeOverlap(t1_1,1,new TGeoCombiTrans(-5,2,1.5,new TGeoRotation("t1_1",-90,90,0)));

TGeoVolume *t2_1=geom->MakeTubs("t2_1",Iron,0,0.8,1,0,360);
 t2_1->SetLineColor(12);
top->AddNodeOverlap(t2_1,1,new TGeoCombiTrans(-5,-2,1.5,new TGeoRotation("t2_1",-90,90,0)));

TGeoVolume *fb1=geom->MakeTubs("fb1",Iron,2,2.3,1,100,260);
 fb1->SetLineColor(12);
top->AddNodeOverlap(fb1,1,new TGeoCombiTrans(-5,0,-1,new TGeoRotation("fb1",90,90,90)));



TGeoVolume *m1=geom->MakeBox("m1",Iron,7,8,4);
 m1->SetLineColor(2);
top->AddNodeOverlap(m1,1,new TGeoCombiTrans(0,0,-17,new TGeoRotation("m1",90,90,0)));

TGeoVolume *m2=geom->MakeTubs("m2",Iron,0,1,7,90,180);
 m2->SetLineColor(2);
top->AddNodeOverlap(m2,1,new TGeoCombiTrans(-3,0,-9,new TGeoRotation("m2",0,90,0)));

TGeoVolume *m3=geom->MakeTubs("m3",Iron,0,1,7,0,90);
 m3->SetLineColor(2);
top->AddNodeOverlap(m3,1,new TGeoCombiTrans(3,0,-9,new TGeoRotation("m3",0,90,0)));

TGeoVolume *m4=geom->MakeBox("m4",Iron,3,7,0.5);
 m4->SetLineColor(2);
top->AddNodeOverlap(m4,1,new TGeoCombiTrans(0,0,-8.5,new TGeoRotation("m4",90,0,90)));

TGeoVolume *m5=geom->MakeTubs("m5",Iron,0,1.5,1.2,0,0);
 m5->SetLineColor(5);
top->AddNodeOverlap(m5,1,new TGeoCombiTrans(0,0,-7.8,new TGeoRotation("m5",0,0,0)));
//.........这里部分代码省略.........

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

  double fFeedback = 0.;
  if (sumIonsTotal > 0.) fFeedback = sumIonsDrift / sumIonsTotal;
  std::cout << "Fraction of ions drifting back: " << fFeedback << "\n"; 

  const double neMean = hElectrons->GetMean();
  std::cout << "Mean number of electrons: " << neMean << "\n";
  const double niMean = hIons->GetMean();
  std::cout << "Mean number of ions: " << niMean << "\n";

  std::cout << "Mean number of electrons on plastic: "
            << sumElectronsPlastic / nEvents << "\n";
  std::cout << "Mean number of ions on plastic: "
            << sumIonsPlastic / nEvents << "\n";
 
  std::cout << "Electron endpoints:\n";
  const double fUpperMetal = sumElectronsUpperMetal / sumElectronsTotal;
  const double fPlastic = sumElectronsPlastic / sumElectronsTotal;
  const double fLowerMetal = sumElectronsLowerMetal / sumElectronsTotal;
  const double fTransfer = sumElectronsTransfer / sumElectronsTotal;
  const double fOther = sumElectronsOther / sumElectronsTotal;
  std::cout << "    upper metal: " << fUpperMetal * 100. << "%\n";
  std::cout << "    plastic:     " << fPlastic * 100. << "%\n";
  std::cout << "    lower metal: " << fLowerMetal * 100. << "%\n";
  std::cout << "    transfer:    " << fTransfer * 100. << "%\n";
  std::cout << "    other:       " << fOther * 100. << "%\n";

  TCanvas* cD = new TCanvas();
  const bool plotGeo = true;
  if (plotGeo && plotDrift) {
    // Build the geometry in Root.
    TGeoManager* geoman = new TGeoManager("world", "geometry");
    TGeoMaterial* matVacuum = new TGeoMaterial("Vacuum", 0, 0, 0);
    TGeoMedium* medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
    TGeoMaterial* matKapton = new TGeoMaterial("Kapton", 12, 6, 1.42);
    TGeoMedium* medKapton = new TGeoMedium("Kapton", 2, matKapton);
    TGeoMaterial* matCopper = new TGeoMaterial("Copper", 63, 29, 8.94);
    TGeoMedium* medCopper = new TGeoMedium("Copper", 3, matCopper);
    TGeoVolume* volTop = geoman->MakeBox("TOP", 
                                         medVacuum, pitch, pitch, 0.02);
    volTop->SetVisibility(0);
    TGeoBBox* shpKapton = new TGeoBBox("K", pitch / 2., 
                                            pitch / 2., 
                                            kapton / 2.);
    TGeoPcon* shpHole = new TGeoPcon("H", 0., 360., 3);
    shpHole->DefineSection(0, -kapton / 2., 0., outdia / 2.);
    shpHole->DefineSection(1,           0., 0., middia / 2.);
    shpHole->DefineSection(2,  kapton / 2., 0., outdia / 2.);

    TGeoCompositeShape* shpGem = new TGeoCompositeShape("G", "K - H");
    TGeoVolume* volKapton = new TGeoVolume("Kapton", shpGem, medKapton);
    volKapton->SetLineColor(kGreen);
    volKapton->SetTransparency(50);

    TGeoBBox* shpMetal = new TGeoBBox("M", pitch / 2., 
                                           pitch / 2., 
                                           metal / 2.);
    TGeoTube* shpTube = new TGeoTube("T", 0., outdia / 2., metal / 2.);
    TGeoCompositeShape* shpElectrode = new TGeoCompositeShape("E", "M - T");
    TGeoVolume* volElectrode = new TGeoVolume("Electrode", 
                                              shpElectrode, medCopper);
    volElectrode->SetLineColor(kBlue);
    volElectrode->SetTransparency(50);

int visualizationTracker(float minZ, float maxZ, float minX, float maxX, float theta, float phi){
    gSystem->Load("libGeom");
//++++++++++++++++++++ Set up stuff ++++++++++++++++++++//
    TGeoManager *geom = new TGeoManager("simple1", "Simple geometry");
//--- define some materials and media
    TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
//--- make the top container volume
    TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 500., 500., 500.);
//TGeoVolume *toptop = geom->MakeBox("TOPTOP", Vacuum, 1000., 1000., 500.);
    geom->SetTopVolume(top);

    int count = 0;
    for (int i = 0; i < _nEntries; ++i){
        _inTree->GetEntry(i);
        if (isRightSubDet()&&(_zVal >= minZ && _zVal < maxZ)&&(_xVal >= minX && _xVal < maxX)/*&&(_rVal <= 12)&&(_rVal >=8)*/){
            char modName[192];
            sprintf(modName, "testModule%i", i);
            TGeoVolume* testMod = geom->MakeBox( modName, Vacuum, 90., 90., 40. );
            getModule( geom, top, testMod );
            count++;
        }
    }

    if(count == 0) return -1;

    getBeamVisuals(geom, top, minZ, maxZ);

//--- close the geometry
    geom->CloseGeometry();
// -- draw
    geom->SetVisLevel(4);

    TCanvas * c = new TCanvas();
    c->SetTheta(theta);
    c->SetPhi(phi);
    top->Draw();

//--- putting words on canvas...
    bool with0T = true;

    //can play with these numbers
    double widthofeach = 0.07;
    double textsize = 0.05;

    double xmax = 2*widthofeach;
    if (with0T) xmax = widthofeach;

    TPaveText* pt = new TPaveText(0,0,xmax,1,"brNDC");
    pt->SetBorderSize(0);
    pt->SetFillStyle(0);
    pt->SetTextAlign(22);
    pt->SetTextFont(42);
    pt->SetTextSize(0.1);
    TText *text = pt->AddText(0,0,TString("#font[42]{"+_line1+"}"));
    text->SetTextSize(textsize);
    text->SetTextAngle(90);
    pt->Draw();

    TPaveText *pt2 = new TPaveText(widthofeach, 0, 2*widthofeach, 1, "brNDC");
    pt2->SetBorderSize(0);
    pt2->SetFillStyle(0);
    pt2->SetTextAlign(22);
    pt2->SetTextFont(42);
    pt2->SetTextSize(0.1);
    TText *text2 = pt2->AddText(0,0,TString("#font[42]{"+_line2+"}"));
    text2->SetTextSize(textsize);
    text2->SetTextAngle(90);
    pt2->Draw();

    TPaveText *pt3 = new TPaveText(2*widthofeach, 0, 3*widthofeach, 1, "brNDC");
    pt3->SetBorderSize(0);
    pt3->SetFillStyle(0);
    pt3->SetTextAlign(22);
    pt3->SetTextFont(42);
    pt3->SetTextSize(0.1);
    TText *text3 = pt3->AddText(0,0,TString("#font[42]{"+_line3+"}"));
    text3->SetTextSize(textsize);
    text3->SetTextAngle(90);
    pt3->Draw();

    string str = string("i") + to_string(_i) + string(".gif");
    c->SaveAs(TString(str));
    gSystem->Exec(TString("mv "+str+" images/"+str));
    delete c;
    cout << "Created image " << str << endl;
    return 0;
}

void rootgeom()
{
//Definition of a simple geometry (the 4 ROOT characters)
//Author: Andrei Gheata
   
   // gStyle->SetCanvasPreferGL(true);
   gSystem->Load("libGeom");
   TGeoManager *geom = new TGeoManager("simple1", "Simple geometry");
   
   //--- 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("Root Material",2, matAl);

   //--- define the transformations
   TGeoTranslation *tr1 = new TGeoTranslation(20., 0, 0.);
   TGeoTranslation *tr2 = new TGeoTranslation(10., 0., 0.);
   TGeoTranslation *tr3 = new TGeoTranslation(10., 20., 0.);
   TGeoTranslation *tr4 = new TGeoTranslation(5., 10., 0.);
   TGeoTranslation *tr5 = new TGeoTranslation(20., 0., 0.);
   TGeoTranslation *tr6 = new TGeoTranslation(-5., 0., 0.);
   TGeoTranslation *tr7 = new TGeoTranslation(7.5, 7.5, 0.);
   TGeoRotation   *rot1 = new TGeoRotation("rot1", 90., 0., 90., 270., 0., 0.);
   TGeoCombiTrans *combi1 = new TGeoCombiTrans(7.5, -7.5, 0., rot1);
   TGeoTranslation *tr8 = new TGeoTranslation(7.5, -5., 0.);
   TGeoTranslation *tr9 = new TGeoTranslation(7.5, 20., 0.);
   TGeoTranslation *tr10 = new TGeoTranslation(85., 0., 0.);
   TGeoTranslation *tr11 = new TGeoTranslation(35., 0., 0.);
   TGeoTranslation *tr12 = new TGeoTranslation(-15., 0., 0.);
   TGeoTranslation *tr13 = new TGeoTranslation(-65., 0., 0.);
   
   TGeoTranslation  *tr14 = new TGeoTranslation(0,0,-100);
   TGeoCombiTrans *combi2 = new TGeoCombiTrans(0,0,100, 
                                   new TGeoRotation("rot2",90,180,90,90,180,0));
   TGeoCombiTrans *combi3 = new TGeoCombiTrans(100,0,0, 
                                   new TGeoRotation("rot3",90,270,0,0,90,180));
   TGeoCombiTrans *combi4 = new TGeoCombiTrans(-100,0,0, 
                                   new TGeoRotation("rot4",90,90,0,0,90,0));
   TGeoCombiTrans *combi5 = new TGeoCombiTrans(0,100,0, 
                                   new TGeoRotation("rot5",0,0,90,180,90,270));
   TGeoCombiTrans *combi6 = new TGeoCombiTrans(0,-100,0, 
                                   new TGeoRotation("rot6",180,0,90,180,90,90));
   
   //--- make the top container volume
   Double_t worldx = 110.;
   Double_t worldy = 50.;
   Double_t worldz = 5.;
   TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 270., 270., 120.);
   geom->SetTopVolume(top);
   TGeoVolume *replica = geom->MakeBox("REPLICA", Vacuum,120,120,120);
   replica->SetVisibility(kFALSE);
   TGeoVolume *rootbox = geom->MakeBox("ROOT", Vacuum, 110., 50., 5.);
   rootbox->SetVisibility(kFALSE);
   
   //--- make letter 'R'
   TGeoVolume *R = geom->MakeBox("R", Vacuum, 25., 25., 5.);
   R->SetVisibility(kFALSE);
   TGeoVolume *bar1 = geom->MakeBox("bar1", Al, 5., 25, 5.);
   bar1->SetLineColor(kRed);
   R->AddNode(bar1, 1, tr1);
   TGeoVolume *bar2 = geom->MakeBox("bar2", Al, 5., 5., 5.);
   bar2->SetLineColor(kRed);
   R->AddNode(bar2, 1, tr2);
   R->AddNode(bar2, 2, tr3);
   TGeoVolume *tub1 = geom->MakeTubs("tub1", Al, 5., 15., 5., 90., 270.);
   tub1->SetLineColor(kRed);
   R->AddNode(tub1, 1, tr4);
   TGeoVolume *bar3 = geom->MakeArb8("bar3", Al, 5.);
   bar3->SetLineColor(kRed);
   TGeoArb8 *arb = (TGeoArb8*)bar3->GetShape();
   arb->SetVertex(0, 15., -5.);
   arb->SetVertex(1, 0., -25.);
   arb->SetVertex(2, -10., -25.);
   arb->SetVertex(3, 5., -5.);
   arb->SetVertex(4, 15., -5.);
   arb->SetVertex(5, 0., -25.);
   arb->SetVertex(6, -10., -25.);
   arb->SetVertex(7, 5., -5.);
   R->AddNode(bar3, 1, gGeoIdentity);
   
   //--- make letter 'O'
   TGeoVolume *O = geom->MakeBox("O", Vacuum, 25., 25., 5.);
   O->SetVisibility(kFALSE);
   TGeoVolume *bar4 = geom->MakeBox("bar4", Al, 5., 7.5, 5.);
   bar4->SetLineColor(kYellow);
   O->AddNode(bar4, 1, tr5);
   O->AddNode(bar4, 2, tr6);
   TGeoVolume *tub2 = geom->MakeTubs("tub1", Al, 7.5, 17.5, 5., 0., 180.);
   tub2->SetLineColor(kYellow);
   O->AddNode(tub2, 1, tr7);
   O->AddNode(tub2, 2, combi1);
   
   //--- make letter 'T'
   TGeoVolume *T = geom->MakeBox("T", Vacuum, 25., 25., 5.);
   T->SetVisibility(kFALSE);
   TGeoVolume *bar5 = geom->MakeBox("bar5", Al, 5., 20., 5.);
   bar5->SetLineColor(kBlue);
   T->AddNode(bar5, 1, tr8);
//.........这里部分代码省略.........

void
show
(
 char* inpDir  = "",         // MuDST directory
 char* inpFile = "show.lis", // MuDST file(s);                      
 char* outFile = "show.root",// output tree file
 Int_t nFiles  = 50,         // # of MuDST file(s)
 Int_t nEvents = 100         // # of events
 )
  // remeber to adjust dbase timestamp below !!!! 
  // what a ... design
{ 
  //gErrorIgnoreLevel=1999;

  // load root/root4star libraries
  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();

  // load more libraries :)
  gSystem->Load("libmysqlclient");
  gSystem->Load("libGeom");
  gSystem->Load("StDbLib");
  gSystem->Load("StDbBroker");
  gSystem->Load("St_db_Maker");

  // load even more libraries (EEMC stuff) 
  gSystem->Load("StEEmcUtil");
  gSystem->Load("StEEmcDbMaker");
  gSystem->Load("StEEmcPoolTTM");

  // create the chain    
  chain = new StChain("StChain"); 
  //

  now = new TDatime;
  // for display
  TCanvas    *c1   = new TCanvas("eemc","eemc",10,10,1000,1000);
  TPaveLabel *tlab = new TPaveLabel(-0.99,+0.99,+0.99,+0.90,
				    "EEMC TOWERS & TPC TRACKS     Piotr A Zolnierczuk (IU)");

  eventInfo = new TPaveText (-0.99,-0.99,+0.0 ,-0.75);
  dateInfo  = new TPaveLabel(+0.60,-0.99,+0.99,-0.95,now->AsString());

  TGeoManager  *gm    = new TGeoManager("eemc", "eemc tower display");
  TGeoVolume   *top   = gm->MakeBox("star",0, 200., 200., 350.);
  TGeoVolume   *smbox = gm->MakeBox("smbox1",0, 2., 2., 2.);
  smbox->SetLineColor(kRed);
  // eemc 
  eemc  = new EEmcTTDisplay();
  eemc->SetMagneticField(0.5); // in Tesla
  eemc->SetShowExtrapolatedTracks(true);
 
  TGeoTranslation *etra = new TGeoTranslation(0.0,0.0,0.5*(eemc->getZ1()+eemc->getZ2()));
  top->AddNode(smbox, 1,NULL);
  top->AddNode(eemc(),1,etra);
  gm->SetTopVolume(top);
  gm->CloseGeometry();
  gm->SetVisLevel(4);
  gm->SetVisOption(0);

  c1->SetTheta(90);
  c1->SetPhi(0);

  top->Draw();
  tlab->Draw();

  gPad->Update();
  
  // now we add Makers to the chain...  some of that is black magic :) 
  muDstMk  = new StMuDstMaker(0,0,inpDir,inpFile,"",nFiles);       // muDST main chain
  StMuDbReader  *db       = StMuDbReader::instance();              // need the database
  St_db_Maker   *dbMk        = new St_db_Maker("StarDb", "MySQL:StarDb");   // need another db(?) 
  new StEEmcDbMaker("eemcDb");                     // need EEMC database  

  // now comment in/out/change the below if you want it your way
  dbMk->setTimeStampDay(20040331);   // format: yyyymmdd

  // finally after so many lines we arrive at the good stuff

  ttm = new  EEmcTTMMaker ("TTM",muDstMk,eemcDbMk);
  ttm->Summary(cout);    // 

  StMuDebug::setLevel(0);

  chain->Init();

  StEEmcDb *eemcDb = (StEEmcDb*)chain->GetDataSet("StEEmcDb");
  eemcDb->setSectors(1,12);            // request EEMC DB for sectors you need (dafault:1-12)
  eemcDb->setPreferedFlavor("onlped","eemcPMTped"); // request alternative flavor

  chain->ls(3);

  //---------------------------------------------------
  next();
}

int main(int argc, char * argv[])
{

   // generate benchmark cases
   TransformationMatrix const * identity = new TransformationMatrix(0,0,0,0,0,0);

   double L = 10.;
   double Lz = 10.;
   const double Sqrt2 = sqrt(2.);

   TGeoManager * geom = new TGeoManager("simple1", "ToyDetector");
   TGeoMaterial * matVacuum = new TGeoMaterial("Vacuum",0,0,0);

   TGeoMedium * vac = new TGeoMedium("Vacuum",1,matVacuum);
   TGeoVolume * world = geom->MakeBox("world", vac, L, L, Lz  );

   geom->SetTopVolume( world );

   TGeoVolume * boxlevel2 = geom->MakeBox( "boxlevel2", vac, Sqrt2*L/2./2., Sqrt2*L/2./2.,Lz);
   TGeoVolume * boxlevel3 = geom->MakeBox( "boxlevel3", vac, L/2./2.,L/2./2.,Lz);
   TGeoVolume * boxlevel1 = geom->MakeBox( "boxlevel1", vac, L/2.,L/2., Lz);

   boxlevel2->AddNode( boxlevel3, 1, new TGeoRotation("mat1",0,0,-45));
   boxlevel1->AddNode( boxlevel2, 1, new TGeoRotation("mat2",0,0,45));
   world->AddNode( boxlevel1, 1, new TGeoTranslation("trans1",-L/2., 0,0 ) );
   world->AddNode( boxlevel1, 2, new TGeoTranslation("trans2",L/2., 0,0 ) );
   geom->CloseGeometry();
   delete world->GetVoxels();
   world->SetVoxelFinder(0);

    // perform basic tests
   TGeoNavigator * nav = geom->GetCurrentNavigator(); // new TGeoNavigator(geom);
   StopWatch timer;

   timer.Start();
   int stepsdone=0;
   for(int n=0;n<1000;n++)
     {
   for(int i=0;i<100000;i++)
   // testing the NavigationAndStepInterface
   {
      int localstepsdone=0;
      double distancetravelled=0.;
      Vector3D p;
      PhysicalVolume::samplePoint( p, L, L, Lz, 1. );

      //std::cerr << p << std::endl;
      //setup point in world
      Vector3D d(1,0,0);

      TGeoNode const * vol;
      nav->SetCurrentPoint( p.x, p.y, p.z );
      nav->SetCurrentDirection( d.x, d.y, d.z);
      vol=nav->FindNode(p.x,p.y,p.z);

      while( vol!=NULL )
      {
         localstepsdone++;

         // do one step ( this will internally adjust the current point and so on )
         vol = nav->FindNextBoundaryAndStep(Utils::kInfinity);

         distancetravelled+=nav->GetStep();
         //double const * p = nav->GetCurrentPoint();
         //double const * pl = nav->GetLastPoint();
         //double const * cd = nav->GetCurrentDirection();

         //std::cerr << " proposed step: " << nav->GetStep();
         //std::cerr << " current point " << p[0] << " " << p[1] << " " << p[2] << std::endl;
         //std::cerr << " last point " << pl[0] << " " << pl[1] << " " << pl[2] << std::endl;
         //std::cerr << " current dir " << cd[0] << " " << cd[1] << " " << cd[2] << std::endl;
      }
   //   std::cerr << localstepsdone << " " << distancetravelled << std::endl;
      stepsdone+=localstepsdone;
   }
     }
   timer.Stop();
   std::cout << " time for 100000 particles " << timer.getDeltaSecs( ) << std::endl;
   std::cout << " average steps done " << stepsdone / 100000. << std::endl;
   std::cout << " time per step " << timer.getDeltaSecs()/stepsdone << std::endl;
}

void geom_Al()
{
// Set the range
//---------------------------------------------------------------------------//

  // 0.314 Mev Ranges (g/cm2):
  /* 0.0025 0.0094 0.0181 0.0255 0.0336 0.0403 0.0477 0.0566 0.0654 0.0721
   * 0.0810 0.0993 */

  // 0.521 Mev Ranges (g/cm2):
  /* 0.0025 0.0094 0.0180 0.0255 0.0335 0.0405 0.0475 0.0566 0.0653 0.0721
   * 0.0807 0.0992 0.1111 0.1259 0.1439 0.1596 0.1825 0.2125 */

  // 1.033 Mev Ranges (g/cm2):
  /* 0.0025 0.0094 0.0180 0.0255 0.0336 0.0402 0.0476 0.0562 0.0654 0.0723
   * 0.0808 0.0990 0.1110 0.1257 0.1440 0.1593 0.1821 0.2122 0.2225 0.2452
   * 0.2521 0.2908 0.3141 0.3533 0.4188 0.4814 */

  // Range of dose depth calculation in cm
  double range = 0.0993;
//---------------------------------------------------------------------------//


  // Set up manager of the geometry world
  gSystem->Load( "libGeom" );
  TGeoManager* geom = new TGeoManager(
    "geom_Al",
    "Geometry for electron energy deposition in Al test problems.");

//---------------------------------------------------------------------------//
// Material and Mixture Definitions
//---------------------------------------------------------------------------//

  // Al ( density given as 2.7 not 2.69890 g/cm^3 )
  density = 2.7;
  TGeoMaterial* mat_1 = new TGeoMaterial( "mat_1", 26.9815385, 13, -density );

  // Void material
  TGeoMaterial* void_mat = new TGeoMaterial( "void", 0, 0, 0 );

  // Graveyard (terminal)
  TGeoMaterial* graveyard_mat = new TGeoMaterial( "graveyard", 0, 0, 0 );

//---------------------------------------------------------------------------//
// Medium Definitions
//---------------------------------------------------------------------------//

  // Al
  TGeoMedium* med_1 = new TGeoMedium( "med_1", 1, mat_1 );

  // Void
  TGeoMedium* void_med = new TGeoMedium( "void_med", 2, void_mat );

  // Graveyard
  TGeoMedium* graveyard_med = new TGeoMedium( "graveyard", 3, graveyard_mat );

//---------------------------------------------------------------------------//
// Volume Definitions
//---------------------------------------------------------------------------//

  //----------------------------
  // Create the calorimeter foil
  //----------------------------

  // Thickness of the calorimeter foil (given in g/cm^2)
  double calorimeter_thickness = 5.050E-03/density;
  // Half the height
  double half_cal_thickness = calorimeter_thickness/2.0;

  // Height (given - largest dimension)
  double height = 8.255;
  // Half side length
  double half_height = height/2.0;

  TGeoVolume* calorimeter_foil =
    geom->MakeBox( "Calorimeter_Foil", med_1, half_height, half_height, half_cal_thickness );

  calorimeter_foil->SetUniqueID( 1 );
  calorimeter_foil->SetLineColor( 2 );

  // //----------------------
  // // Create the front foil
  // //----------------------

  // Calculate the range = 0.4814;
  double range_cm = range/density;

  // Calculate the thickness of the front foil
  double front_thickness = range_cm - half_cal_thickness;

  TGeoVolume* front_foil;
  if (front_thickness > 0.0 )
  {
    // Height (given - largest dimension)
    height = 8.89;
    // Half side length
    half_height = height/2.0;

    front_foil =
      geom->MakeBox( "Front_Foil", med_1, half_height, half_height, front_thickness/2.0 );
//.........这里部分代码省略.........

void LCDDImp::dump() const {
  TGeoManager* mgr = m_manager;
  mgr->SetVisLevel(4);
  mgr->SetVisOption(1);
  m_worldVol->Draw("ogl");
}

//.........这里部分代码省略.........
//                for(int reps=0;reps<NREPS;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<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();

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 *too