C++ PHG4CylinderSubsystem::set_double_param方法代码示例

//! A rough version of the inner HCal support ring, from Richie's CAD drawing. - Jin
void HCalInner_SupportRing(PHG4Reco* g4Reco,
			   const int absorberactive = 0) {
  
  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  const double z_ring1 = (2025 + 2050) / 2. / 10.;
  const double z_ring2 = (2150 + 2175) / 2. / 10.;
  const double dz = 25. / 10.;
  const double innerradius_sphenix = 116.;
  const double innerradius_ephenix_hadronside = 138.;
  const double maxradius = 178.0 - 0.001; // avoid touching the outer HCal envelop volumne
  const double z_rings[] =
    { -z_ring2, -z_ring1, z_ring1, z_ring2 };

  PHG4CylinderSubsystem *cyl;

  for (int i = 0; i < 4; i++)
    {
      double innerradius = innerradius_sphenix;
      if ( z_rings[i] > 0 && inner_hcal_eic == 1)
      {
        innerradius = innerradius_ephenix_hadronside;
      }
      cyl = new PHG4CylinderSubsystem("HCALIN_SPT_N1", i);
      cyl->set_double_param("place_z",z_rings[i]);
      cyl->SuperDetector("HCALIN_SPT");
      cyl->set_double_param("radius",innerradius);
      cyl->set_int_param("lengthviarapidity",0);
      cyl->set_double_param("length",dz);
      cyl->set_string_param("material","SS310");
      cyl->set_double_param("thickness",maxradius - innerradius);
      if (absorberactive)
	{
	  cyl->SetActive();
	}
      g4Reco->registerSubsystem(cyl);
    }
  
  return;
}

//! EMCal setup macro - 1D azimuthal projective SPACAL
double
CEmc_1DProjectiveSpacal(PHG4Reco* g4Reco, double radius, const int crossings, const int absorberactive = 0)
{
  double emc_inner_radius = 95.; // emc inner radius from engineering drawing
  double cemcthickness = 12.7;
  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius

  
  if (radius > emc_inner_radius) {
    cout << "inconsistency: preshower radius+thickness: " << radius 
	 << " larger than emc inner radius: " <<  emc_inner_radius
	 << endl;
    gSystem->Exit(-1);
  }
  
  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

//  boundary check
  if (radius> emc_inner_radius - 1.5 - no_overlapp)
    {
      cout <<"G4_CEmc_Spacal.C::CEmc() - expect radius < "<<emc_inner_radius - 1.5 - no_overlapp<<" to install SPACAL"<<endl;
      exit(1);
    }
  radius = emc_inner_radius - 1.5 - no_overlapp;

  // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
  cyl->SuperDetector("CEMC_ELECTRONICS");
  cyl->set_double_param("radius",radius);
  cyl->set_string_param("material","G4_TEFLON");
  cyl->set_double_param("thickness",1.5);
  if (absorberactive)  cyl->SetActive();
  g4Reco->registerSubsystem( cyl );

  radius += 1.5;
  radius += no_overlapp;


  int ilayer = Min_cemc_layer;
  PHG4SpacalSubsystem *cemc;
  cemc = new PHG4SpacalSubsystem("CEMC", ilayer);

  cemc ->get_geom().set_radius(emc_inner_radius);
  cemc ->get_geom().set_thickness(cemcthickness);
  cemc ->get_geom().set_construction_verbose(1);

  cemc->SetActive();
  cemc->SuperDetector("CEMC");
  if (absorberactive)  cemc->SetAbsorberActive();
  cemc->OverlapCheck(overlapcheck);
  
  g4Reco->registerSubsystem( cemc );

  if (ilayer > Max_cemc_layer)
    {
      cout << "layer discrepancy, current layer " << ilayer
           << " max cemc layer: " << Max_cemc_layer << endl;
    }

  radius += cemcthickness;
  radius += no_overlapp;

  // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
  cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
  cyl->SuperDetector("CEMC_SPT");
  cyl->set_double_param("radius",radius);
cyl->set_string_param("material","SS310"); // SS310 Stainless Steel
 cyl->set_double_param("thickness",0.5);
  if (absorberactive)
    cyl->SetActive();
  g4Reco->registerSubsystem(cyl);

  radius += 0.5;
  radius += no_overlapp;

  return radius;
}

//! test for 2D projective SPACAL using 1D projective modules. For evaluation only
double
CEmc_Proj(PHG4Reco* g4Reco, double radius, const int crossings, const int absorberactive = 0)
{
  double emc_inner_radius = 95.; // emc inner radius from engineering drawing
  double cemcthickness = 12.9+1.5;
  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius
  
  if (radius > emc_inner_radius) {
    cout << "inconsistency: preshower radius+thickness: " << radius 
	 << " larger than emc inner radius: " <<  emc_inner_radius
	 << endl;
    gSystem->Exit(-1);
  }
  
  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  cout <<"Get_Spacal_Tilt() = "<<Get_Spacal_Tilt()<<endl;

  // the radii are only to determined the thickness of the cemc
  //double emc_inner_radius = radius; // emc inner radius from engineering spreadsheet
  // double cemcthickness = 12.9+1.5;
  //  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius
  radius = emc_outer_radius;
  
  int ilayer = Min_cemc_layer;
  PHG4SpacalSubsystem *cemc;
  cemc = new PHG4SpacalSubsystem("CEMC", ilayer);

  cemc ->get_geom().set_radius(emc_inner_radius);
  cemc ->get_geom().set_thickness(cemcthickness);
  cemc ->get_geom().set_construction_verbose(2);

  cemc ->get_geom().set_config(PHG4CylinderGeom_Spacalv1::kProjective_PolarTaper);
  cemc ->get_geom().set_azimuthal_tilt(Get_Spacal_Tilt());
//  cemc ->get_geom().set_azimuthal_tilt(-2.3315/2./95*1.9);
//  cemc ->get_geom().set_azimuthal_tilt(6.28/256*2);
//  cemc ->get_geom().set_azimuthal_tilt(-6.28/256*4);
//  cemc ->get_geom().set_azimuthal_seg_visible(true);
//  cemc ->get_geom().set_virualize_fiber(false);
//  cemc ->get_geom().set_assembly_spacing(0.001);
  cemc ->get_geom().set_polar_taper_ratio(1.128);
//  cemc ->get_geom().set_polar_taper_ratio(1.123);
//  cemc ->get_geom().set_polar_taper_ratio(1.117);
//  cemc ->get_geom().set_absorber_mat("G4_AIR");
//  cemc ->get_geom().set_azimuthal_n_sec(256/2);

  cemc->SetActive();
  cemc->SuperDetector("CEMC");
  if (absorberactive)  cemc->SetAbsorberActive();
  cemc->OverlapCheck(overlapcheck);

  g4Reco->registerSubsystem( cemc );

  if (ilayer > Max_cemc_layer)
    {
      cout << "layer discrepancy, current layer " << ilayer
           << " max cemc layer: " << Max_cemc_layer << endl;
    }

  radius += cemcthickness;
  radius += no_overlapp;

  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("EMCELECTRONICS", 0);
  cyl->set_double_param("radius",radius);
  cyl->set_string_param("material","G4_TEFLON");
  cyl->set_double_param("thickness",0.5);
  if (absorberactive)  cyl->SetActive();
  g4Reco->registerSubsystem( cyl );
  radius += 0.5;
  radius += no_overlapp;

  return radius;
}

double Svtx(PHG4Reco* g4Reco, double radius,
            const int absorberactive = 0,
            int verbosity = 0)
{
  if (n_maps_layer > 0)
    {
      bool maps_overlapcheck = false;  // set to true if you want to check for overlaps
      
      // MAPS inner barrel layers
      //======================================================
      
      double maps_layer_radius[3] = {24.61, 32.59, 39.88}; // mm - numbers from Walt 6 Aug 2018
      
      // D. McGlinchey 6Aug2018 - type no longer is used, included here because I was too lazy to remove it from the code
      int stave_type[3] = {0, 0, 0};
      int staves_in_layer[3] = {12, 16, 20};       // Number of staves per layer in sPHENIX MVTX
      double phi_tilt[3] = {0.300, 0.305, 0.300}; // radians - numbers from Walt 6 Aug 2018
      
      for (int ilayer = 0; ilayer < n_maps_layer; ilayer++)
	{
	  if (verbosity)
	    cout << "Create Maps layer " << ilayer << " with radius " << maps_layer_radius[ilayer] << " mm, stave type " << stave_type[ilayer]
		 << " pixel size 30 x 30 microns "
		 << " active pixel thickness 0.0018 microns" << endl;
	  
	  PHG4MapsSubsystem* lyr = new PHG4MapsSubsystem("MAPS", ilayer, stave_type[ilayer]);
	  lyr->Verbosity(verbosity);
	  
	  lyr->set_double_param("layer_nominal_radius", maps_layer_radius[ilayer]);  // thickness in cm
	  lyr->set_int_param("N_staves", staves_in_layer[ilayer]);       // uses fixed number of staves regardless of radius, if set. Otherwise, calculates optimum number of staves
	  
	  // The cell size is used only during pixilization of sensor hits, but it is convemient to set it now because the geometry object needs it
	  lyr->set_double_param("pixel_x", 0.0030);          // pitch in cm
	  lyr->set_double_param("pixel_z", 0.0030);          // length in cm
	  lyr->set_double_param("pixel_thickness", 0.0018);  // thickness in cm
	  lyr->set_double_param("phitilt", phi_tilt[ilayer]);
	  
	  lyr->set_int_param("active", 1);
	  lyr->OverlapCheck(maps_overlapcheck);
	  
	  //lyr->set_string_param("stave_geometry_file", "/phenix/hhj3/dcm07e/sPHENIX/macros/macros/g4simulations/mvtx_stave_v01.gdml");
	  lyr->set_string_param("stave_geometry_file", string(getenv("CALIBRATIONROOT")) + string("/Tracking/geometry/mvtx_stave_v01.gdml"));

	  g4Reco->registerSubsystem(lyr);
	  
	  radius = maps_layer_radius[ilayer];
	}
    }
  
  if (n_intt_layer > 0)
    {
      //-------------------
      // INTT ladders
      //-------------------
      
      bool intt_overlapcheck = false;  // set to true if you want to check for overlaps
      
      // instantiate the Silicon tracker subsystem and register it
      // We make one instance of PHG4TrackerSubsystem for all four layers of tracker
      // dimensions are in mm, angles are in radians
      
      // PHG4SiliconTrackerSubsystem creates the detetor layer using PHG4SiliconTrackerDetector
      // and instantiates the appropriate PHG4SteppingAction
      const double intt_radius_max = 140.;  // including stagger radius (mm)
      
      // The length of vpair is used to determine the number of layers
      std::vector<std::pair<int, int>> vpair;  // (sphxlayer, inttlayer)
      for (int i = 0; i < n_intt_layer; i++)
	{
	  // We want the sPHENIX layer numbers for the INTT to be from n_maps_layer to n_maps_layer+n_intt_layer - 1
	  vpair.push_back(std::make_pair(n_maps_layer + i, i));  // sphxlayer=n_maps_layer+i corresponding to inttlayer=i
	  if (verbosity) cout << "Create strip tracker layer " << vpair[i].second << " as  sphenix layer  " << vpair[i].first << endl;
	}
      
      PHG4SiliconTrackerSubsystem* sitrack = new PHG4SiliconTrackerSubsystem("SILICON_TRACKER", vpair);
      sitrack->Verbosity(verbosity);
      sitrack->SetActive(1);
      sitrack->OverlapCheck(intt_overlapcheck);
      g4Reco->registerSubsystem(sitrack);
      
      // Update the laddertype and ladder spacing configuration
      for(int i=0;i<n_intt_layer;i++)
	{
	  sitrack->set_int_param(i, "laddertype", laddertype[i]);
	  sitrack->set_int_param(i, "nladder", nladder[i]);
	  sitrack->set_double_param(i,"sensor_radius", sensor_radius[i]);  // expecting cm
	  sitrack->set_double_param(i,"offsetphi",offsetphi[i]);  // expecting degrees
	}
      
      // outer radius marker (translation back to cm)
      radius = intt_radius_max * 0.1;
    }
  
  //  int verbosity = 1;
  
  // time projection chamber layers --------------------------------------------
  
  // switch ONLY for backward compatibility with 40 layer hits files!
  if (tpc_layers_40)
    {
//.........这里部分代码省略.........

double Svtx(PHG4Reco* g4Reco,
	    double radius,
	    const int absorberactive = 0,
	    int verbosity = 0)
{

  float svtx_inner_radius = 2.3;
  
  if (radius > svtx_inner_radius) {
    cout << "inconsistency: radius: " << radius 
	 << " larger than SVTX inner radius: " << svtx_inner_radius << endl;
    gSystem->Exit(-1);
  }
  
  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  PHG4CylinderSubsystem *cyl;

  //======================================================================================================
  // The thicknesses from Yasuyuki on June 12, 2014 are as follows:
  // For Si 1mm = 1.07% X_0
  // For Cu 1mm = 6.96% X_0
  // The thickness of the tracking layers is:
  // Pixels:         1.3% X_0  (0.21% sensor +  1.07% support)  sensor = 200 mc Si, support = 154 mc Cu
  // Stripixel:      5%   X_0  (0.67% sensor + 4.3% support)    sensor = 624 mc Si, support = 618 mc Cu
  // Outer strips:   2%   X_0 (conservative)  (0.34% sensor + 1.66% support) sensor = 320 mc Si, support = 238 mc Cu
  //=======================================================================================================

  //double svxrad[7] = {2.3, 3.2, 3.9, 19.6, 24.5, 34.4, 39.3}; // ALICE ITS upgrade layer radii in cm
  double svxrad[7] = {svtx_inner_radius, 3.2, 3.9, 19.6, 24.5, 34.4, 64.0}; // ALICE ITS upgrade with outer layer pushed out
  double si_thickness[7] = {0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005};  // ALICE ITS upgrade Si thickness is 50 microns
  double length[7] = {27.0, 27.0, 27.0, 88.0, 88.0, 148.0, 148.0}; // ALICE ladder lengths (stave - 2cm)

  // ALICE ITS total thickness (% of X_0 of 0.3, 0.3, 0.3, 0.8, 0.8, 0.8, 0.8
  // Pixel chip thickness of 50 um (% of X_0 of 0.05 x 1.07 = 0.053%) in all layers
  // so inner 3 layers support thickness = 0.3 - 0.053 = 0.25%, outer 4 layers = 0.8 - 0.053 = 0.75%
  // Support thickness equivalent for Cu in inner layer = 0.25%/6.96% x 1 mm = 0.036 mm = 0.0036 cm
  // Support thickness equivalent for Cu in outer layer = 0.75%/6.96% x 1 mm = 0.108 mm = 0.0108 cm
  double support_thickness[7] = {0.0036, 0.0036, 0.0036, 0.0108, 0.0108, 0.0108, 0.0108};

  // here is our silicon:
  double inner_radius = radius;
  for (int ilayer = Min_si_layer; ilayer <= Max_si_layer; ilayer++)
    {
      cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
      radius = svxrad[ilayer];
      // protect against installing layer with radius < inner radius from argument
      if (radius < inner_radius)
	{
	  cout << "current radius " << radius << " smaller than inner radius "
	       << inner_radius << endl;
	  gSystem->Exit(-1);
	}
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Si");
      cyl->set_double_param("thickness",si_thickness[ilayer]);
      cyl->SetActive();
      cyl->SuperDetector("SVTX");
      g4Reco->registerSubsystem( cyl );

      radius += si_thickness[ilayer] + no_overlapp;
      cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Cu");
      cyl->set_double_param("thickness",support_thickness[ilayer]);
      if (absorberactive)  cyl->SetActive();
      cyl->SuperDetector("SVTXSUPPORT");
      g4Reco->registerSubsystem( cyl );
    }
  if (ilayer != (Max_si_layer+1)) // coming out of the loop, layer is layer+1
    {
      cout << "layer number mismatch for Max_si_layer, Max_si_layer "
           << Max_si_layer << " should be " << ilayer << endl;
      gSystem->Exit(-1);
    }
  radius += support_thickness[Max_si_layer] + no_overlapp;

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

double Magnet(PHG4Reco* g4Reco,
              double radius,
              const int crossings = 0,
              const int absorberactive = 0,
              int verbosity = 0) {

  double magnet_inner_cryostat_wall_radius = 142;
  double magnet_inner_cryostat_wall_thickness = 1;
  double magnet_outer_cryostat_wall_radius = 174.5;
  double magnet_outer_cryostat_wall_thickness = 2.5;
  double magnet_coil_radius = 150.8;
  double magnet_coil_thickness = 9.38;
  double magnet_length = 379.;
  double coil_length = 361.5;
  if (radius > magnet_inner_cryostat_wall_radius) {
    cout << "inconsistency: radius: " << radius
         << " larger than Magnet inner radius: " << magnet_inner_cryostat_wall_radius << endl;
    gSystem->Exit(-1);
  }

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  radius = magnet_inner_cryostat_wall_radius;
  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("MAGNET", 0);
  cyl->set_double_param("radius",magnet_inner_cryostat_wall_radius);
  cyl->set_int_param("lengthviarapidity",0);
  cyl->set_double_param("length",magnet_length);
  cyl->set_double_param("thickness",magnet_inner_cryostat_wall_thickness);
  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
  cyl->SuperDetector("MAGNET");
  if (absorberactive)  cyl->SetActive();
  g4Reco->registerSubsystem( cyl );

  cyl = new PHG4CylinderSubsystem("MAGNET", 1);
  cyl->set_double_param("radius",magnet_coil_radius);
  cyl->set_int_param("lengthviarapidity",0);
  cyl->set_double_param("length",coil_length);
  cyl->set_double_param("thickness",magnet_coil_thickness);
  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
  cyl->SuperDetector("MAGNET");
  if (absorberactive)  cyl->SetActive();
  g4Reco->registerSubsystem( cyl );

  cyl = new PHG4CylinderSubsystem("MAGNET", 2);
  cyl->set_double_param("radius",magnet_outer_cryostat_wall_radius);
  cyl->set_int_param("lengthviarapidity",0);
  cyl->set_double_param("length",magnet_length);
  cyl->set_double_param("thickness",magnet_outer_cryostat_wall_thickness);
  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
  cyl->SuperDetector("MAGNET");
  if (absorberactive)  cyl->SetActive();
  g4Reco->registerSubsystem( cyl );

  radius = magnet_outer_cryostat_wall_radius + magnet_outer_cryostat_wall_thickness; // outside of magnet

  if (verbosity > 0) {
    cout << "========================= G4_Magnet.C::Magnet() ===========================" << endl;
    cout << " MAGNET Material Description:" << endl;
    cout << "  inner radius = " << magnet_inner_cryostat_wall_radius << " cm" << endl;
    cout << "  outer radius = " << magnet_outer_cryostat_wall_radius + magnet_outer_cryostat_wall_thickness << " cm" << endl;
    cout << "  length = " << magnet_length << " cm" << endl;
    cout << "===========================================================================" << endl;
  }

  radius += no_overlapp;

  return radius;
}

//.........这里部分代码省略.........
  stringline >> fieldstrength;
  if (stringline.fail()) { // conversion to double fails -> we have a string

    if (field.find("sPHENIX.root") != string::npos) {
      g4Reco->set_field_map(field, PHFieldConfig::Field3DCartesian);
    } else {
      g4Reco->set_field_map(field, PHFieldConfig::kField2D);
    }
  } else {
    g4Reco->set_field(fieldstrength); // use const soleniodal field
  }
  g4Reco->set_field_rescale(magfield_rescale);
  
  double radius = 0.;

  //----------------------------------------
  // PIPE
  if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
  
  //----------------------------------------
  // TRACKING
  if (do_tracking) radius = Tracking(g4Reco, radius, absorberactive);

  //----------------------------------------
  // PSTOF
  
  if (do_pstof) radius = PSTOF(g4Reco, radius, absorberactive);

  //----------------------------------------
  // CEMC
//
  if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
//  if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
  
  //----------------------------------------
  // HCALIN
  
  if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);

  //----------------------------------------
  // MAGNET
  
  if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);

  //----------------------------------------
  // HCALOUT
  
  if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);

  //----------------------------------------
  // sPHENIX forward flux return door
  if (do_plugdoor) PlugDoor(g4Reco, absorberactive);

  //----------------------------------------
  // BLACKHOLE
  
  // swallow all particles coming out of the backend of sPHENIX
  PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
blackhole->set_double_param("radius",radius + 10); // add 10 cm

  blackhole->set_int_param("lengthviarapidity",0);
  blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
  blackhole->BlackHole();
  blackhole->set_double_param("thickness",0.1); // it needs some thickness
  blackhole->SetActive(); // always see what leaks out
  blackhole->OverlapCheck(overlapcheck);
  g4Reco->registerSubsystem(blackhole);

  //----------------------------------------
  // FORWARD BLACKHOLEs
  // +Z
  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
  blackhole->SuperDetector("BH_FORWARD_PLUS");
  blackhole->set_double_param("radius",0); // add 10 cm
  blackhole->set_int_param("lengthviarapidity",0);
  blackhole->set_double_param("length",0.1); // make it cover the world in length
  blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1  - no_overlapp);
  blackhole->BlackHole();
  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
  blackhole->SetActive(); // always see what leaks out
  blackhole->OverlapCheck(overlapcheck);
  g4Reco->registerSubsystem(blackhole);

  blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
  blackhole->SuperDetector("BH_FORWARD_NEG");
  blackhole->set_double_param("radius",0); // add 10 cm
  blackhole->set_int_param("lengthviarapidity",0);
  blackhole->set_double_param("length",0.1); // make it cover the world in length
  blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1  + no_overlapp);
  blackhole->BlackHole();
  blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
  blackhole->SetActive(); // always see what leaks out
  blackhole->OverlapCheck(overlapcheck);
  g4Reco->registerSubsystem(blackhole);

  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
  g4Reco->registerSubsystem(truth);
  se->registerSubsystem( g4Reco );
  return 0;
}

//! Babar DIRC (Without most of support structure)
//! Ref: I. Adam et al. The DIRC particle identification system for the BaBar experiment.
//! Nucl. Instrum. Meth., A538:281-357, 2005. doi:10.1016/j.nima.2004.08.129.
double
DIRCSetup(PHG4Reco* g4Reco)
{
  const double radiator_R = 83.65;
  const double length = 470;
  const double z_shift = -115;
  const double z_start = z_shift + length / 2.;
  const double z_end = z_shift - length / 2.;

  PHG4SectorSubsystem *dirc;
  dirc = new PHG4SectorSubsystem("DIRC");
  dirc->get_geometry().set_normal_polar_angle(3.14159265358979323846/2);
  dirc->get_geometry().set_normal_start(
                                        83.65 * PHG4Sector::Sector_Geometry::Unit_cm());
  dirc->get_geometry().set_min_polar_angle(atan2(radiator_R, z_start));
  dirc->get_geometry().set_max_polar_angle(atan2(radiator_R, z_end));
  dirc->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
  dirc->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
  dirc->get_geometry().set_material("Quartz");
  dirc->get_geometry().set_N_Sector(12);
  dirc->OverlapCheck(overlapcheck);
  dirc->get_geometry().AddLayer("Radiator", "Quartz",
                                1.7 * PHG4Sector::Sector_Geometry::Unit_cm(), true);
  g4Reco->registerSubsystem(dirc);

  PHG4CylinderSubsystem *cyl;

  //  The cylinder skins provide most of the strength
  //  and stiffness of the CST. The thickness of the inner
  //  and outer skins is 1.27 and 0.76 mm, respectively

  // Inner skin:
  cyl = new PHG4CylinderSubsystem("DIRC_CST_Inner_Skin", 10);
  cyl->set_double_param("radius",81.71);
  cyl->set_int_param("lengthviarapidity",0);
  cyl->set_double_param("length",length);
  cyl->set_string_param("material","G4_Al");
  cyl->set_double_param("thickness",0.127);
  cyl->set_double_param("place_x",0.);
  cyl->set_double_param("place_y",0.);
  cyl->set_double_param("place_z",z_shift);
  cyl->SetActive(0);
  cyl->SuperDetector("DIRC");
  cyl->OverlapCheck(overlapcheck);

  g4Reco->registerSubsystem(cyl);

  // Outer skin:
  cyl = new PHG4CylinderSubsystem("DIRC_CST_Outer_Skin", 11);
  cyl->set_double_param("radius",89.25 - 0.076);
  cyl->set_int_param("lengthviarapidity",0);
  cyl->set_double_param("length",length);
  cyl->set_string_param("material","G4_Al");
  cyl->set_double_param("thickness",0.076);
  cyl->set_double_param("place_x",0.);
  cyl->set_double_param("place_y",0.);
  cyl->set_double_param("place_z",z_shift);
  cyl->SetActive(0);
  cyl->SuperDetector("DIRC");
  cyl->OverlapCheck(overlapcheck);

  g4Reco->registerSubsystem(cyl);

  // Done
  return 89.25;

}

//! 2D full projective SPACAL
double
CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings,
                        const int absorberactive = 0)
{
  double emc_inner_radius = 92;  // emc inner radius from engineering drawing
  double cemcthickness = 24.00000 - no_overlapp;

  //max radius is 116 cm;
  double emc_outer_radius = emc_inner_radius + cemcthickness;  // outer radius
  assert(emc_outer_radius < 116);

  if (radius > emc_inner_radius)
  {
    cout << "inconsistency: preshower radius+thickness: " << radius
         << " larger than emc inner radius: " << emc_inner_radius << endl;
    gSystem->Exit(-1);
  }

  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");

  // the radii are only to determined the thickness of the cemc
  radius = emc_inner_radius;

  //---------------
  // Load libraries
  //---------------

  // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
  cyl->set_double_param("radius", radius);
  cyl->set_string_param("material", "G4_TEFLON");
  cyl->set_double_param("thickness", 1.5 - no_overlapp);
  cyl->SuperDetector("CEMC_ELECTRONICS");
  cyl->OverlapCheck(overlapcheck);
  if (absorberactive) cyl->SetActive();
  g4Reco->registerSubsystem(cyl);

  radius += 1.5;
  cemcthickness -= 1.5 + no_overlapp;

  // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
  cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
  cyl->SuperDetector("CEMC_SPT");
  cyl->set_double_param("radius", radius + cemcthickness - 0.5);
  cyl->set_string_param("material", "SS310");  // SS310 Stainless Steel
  cyl->set_double_param("thickness", 0.5 - no_overlapp);
  cyl->OverlapCheck(overlapcheck);
  if (absorberactive)
    cyl->SetActive();
  g4Reco->registerSubsystem(cyl);

  cemcthickness -= 0.5 + no_overlapp;

  //---------------
  // Load libraries
  //---------------

  int ilayer = 0;
  PHG4SpacalSubsystem *cemc;

  const bool use_2015_design = false;
  if (use_2015_design)
  {
    cemc = new PHG4SpacalSubsystem("CEMC", ilayer);

    cemc->set_int_param("config", PHG4CylinderGeom_Spacalv1::kFullProjective_2DTaper_SameLengthFiberPerTower);
    cemc->set_double_param("radius", radius);            // overwrite minimal radius
    cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness
    cemc->set_int_param("azimuthal_n_sec", 32);
    //    cemc->set_int_param("construction_verbose", 2);

    cemc->SetActive();
    cemc->SuperDetector("CEMC");
    if (absorberactive)
      cemc->SetAbsorberActive();
    cemc->OverlapCheck(overlapcheck);
  }

  else
  {
    cemc = new PHG4SpacalSubsystem("CEMC", ilayer);

    cemc->set_int_param("virualize_fiber", 0);
    cemc->set_int_param("azimuthal_seg_visible", 1);
    cemc->set_int_param("construction_verbose", 0);
    cemc->Verbosity(0);

    cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
    cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2017ProjTilted/"));
    cemc->set_double_param("radius", radius);            // overwrite minimal radius
    cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness

    cemc->SetActive();
    cemc->SuperDetector("CEMC");
    if (absorberactive)
//.........这里部分代码省略.........

double Svtx(PHG4Reco* g4Reco, double radius, 
	    const int absorberactive = 0,
	    int verbosity = 0) {
  
  float svtx_inner_radius = 2.71;
  
  if (radius > svtx_inner_radius) {
    cout << "inconsistency: radius: " << radius 
	 << " larger than SVTX inner radius: " << svtx_inner_radius << endl;
    gSystem->Exit(-1);
  }
 
  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  PHG4CylinderSubsystem *cyl;
 
  //======================================================================================================
  // The thicknesses from Yasuyuki on June 12, 2014 are as follows:
  // For Si 1mm = 1.07% X_0
  // For Cu 1mm = 6.96% X_0
  // The thickness of the tracking layers is:
  // Pixels:         1.3% X_0  (0.21% sensor +  1.07% support)  sensor = 200 mc Si, support = 154 mc Cu
  // Stripixel:      5%   X_0  (0.67% sensor + 4.3% support)    sensor = 624 mc Si, support = 618 mc Cu
  // Outer strips:   2%   X_0 (conservative)  (0.34% sensor + 1.66% support) sensor = 320 mc Si, support = 238 mc Cu
  //=======================================================================================================

  double si_thickness[7] = {0.02, 0.02, 0.032, 0.032, 0.032, 0.032, 0.032};
  double svxrad[7] = {svtx_inner_radius, 4.63, 9.5, 10.5, 44.5, 45.5, 80.0}; // provides 98 MeV Upsilon resolution
  // Thicknesses (in % X_0) of 1.3,1.3,2.7/2,2.7/2,2.0/2,2.0/2,2.0 - YA's most conservative case
  double support_thickness[7] = {0.0154, 0.0154, 0.0338/2.0, 0.0338/2.0, 0.0238/2.0, 0.0238/2.0, 0.0238};
  double length[7] = {20., 20., -1, -1., - 1., - 1., -1}; // -1 use eta coverage to determine length

  // here is our silicon:
  double inner_radius = radius;
  for (int ilayer = Min_si_layer; ilayer <= Max_si_layer; ilayer++)
    {
      cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
      radius = svxrad[ilayer];
      // protect against installing layer with radius < inner radius from argument
      if (radius < inner_radius)
	{
	  cout << "current radius " << radius << " smaller than inner radius "
	       << inner_radius << endl;
	  gSystem->Exit(-1);
	}
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Si");
      cyl->set_double_param("thickness",si_thickness[ilayer]);
      cyl->SetActive();
      cyl->SuperDetector("SVTX");
      g4Reco->registerSubsystem( cyl );

      radius += si_thickness[ilayer] + no_overlapp;
      cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Cu");
      cyl->set_double_param("thickness",support_thickness[ilayer]);
      if (absorberactive)  cyl->SetActive();
      cyl->SuperDetector("SVTXSUPPORT");
      g4Reco->registerSubsystem( cyl );
    }
  if (ilayer != (Max_si_layer+1)) // coming out of the loop, layer is layer+1
    {
      cout << "layer number mismatch for Max_si_layer, Max_si_layer "
           << Max_si_layer << " should be " << ilayer << endl;
      gSystem->Exit(-1);
    }
  radius += support_thickness[Max_si_layer] + no_overlapp;

  if (verbosity > 0) {
    cout << "============================ G4_Svtx.C::Svtx() ============================" << endl;
    cout << " SVTX Material Description:" << endl;
    for (int ilayer = Min_si_layer; ilayer <= Max_si_layer; ilayer++) {
      cout << "  layer " << ilayer
	   << "  radius " << svxrad[ilayer]
	   << "  zlength " << length[ilayer]
	   << "  thickness (Si) " << si_thickness[ilayer]
//.........这里部分代码省略.........

double Svtx(PHG4Reco* g4Reco, double radius, 
	    const int absorberactive = 0,
	    int verbosity = 0)
{

  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libg4detectors.so");
  gSystem->Load("libg4testbench.so");

  //---------------------------------
  // Inner Cylinder layers for pixels
  //--------------------------------- 

  PHG4CylinderSubsystem *cyl;
 
  //======================================================================================================
  // The thicknesses from Yasuyuki on June 12, 2014 are as follows:
  // For Si 1mm = 1.07% X_0
  // For Cu 1mm = 6.96% X_0
  // The thickness of the tracking layers is:
  // Pixels:         1.3% X_0  (0.21% sensor +  1.07% support)  sensor = 200 mc Si, support = 154 mc Cu
  //=======================================================================================================

  double si_thickness[2] = {0.02, 0.02};
  double svxrad[2] = {2.71, 4.63};
  double support_thickness[2] = {0.0154, 0.0154};
  double length[2] = {20., 20.};

  // here is our silicon:
  double inner_radius = radius;
  for (int ilayer = Min_si_layer; ilayer < 2; ilayer++)
    {
      cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
      cyl->Verbosity(verbosity);
      radius = svxrad[ilayer];
      // protect against installing layer with radius < inner radius from argument
      if (radius < inner_radius)
	{
	  cout << "current radius " << radius << " smaller than inner radius "
	       << inner_radius << endl;
	  gSystem->Exit(-1);
	}
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Si");
      cyl->set_double_param("thickness",si_thickness[ilayer]);
      cyl->SetActive();
      cyl->SuperDetector("SVTX");
      g4Reco->registerSubsystem( cyl );

      radius += si_thickness[ilayer] + no_overlapp;
      cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
      cyl->Verbosity(verbosity);
      cyl->set_double_param("radius",radius);
      if (length[ilayer] > 0)
        {
          cyl->set_int_param("lengthviarapidity",0);
          cyl->set_double_param("length",length[ilayer]);
        }
      else
	{
          cyl->set_int_param("lengthviarapidity",1);
	}
      cyl->set_string_param("material","G4_Cu");
      cyl->set_double_param("thickness",support_thickness[ilayer]);
      if (absorberactive)  cyl->SetActive();
      cyl->SuperDetector("SVTXSUPPORT");
      g4Reco->registerSubsystem( cyl );
    }

  //--------------------------------
  // Outer Silicon tracking subsytem
  //--------------------------------
  
  bool overlapcheck = false; // set to true if you want to check for overlaps
  
  // instantiate the Silicon tracker subsystem and register it
  // We make one instance of PHG4TrackerSubsystem per layer of tracker
  // dimensions are in mm, angles are in radians
 
  bool option_double_layer[5] = {false, true, false, true, false};
  double layer_radius[5] = {85.0, 85.0, 400.0, 400.0, 800.0};
  int N_strips_sensor_phi[5] = {256, 256, 512, 512, 1536}; 
  double radius_stagger[5] = {10.0, 10.0, 10.0, 10.0, 10.0};
  int N_staggers[5] = {4, 4, 4, 4, 2};
  bool add_lower_roc[5] = {false, false, true, true, true};
  double strip_tilt[5] = {0.0156, -0.0156, 0.0156, -0.0156, 0.0};  // radians, usually 1:64 tilt
  //double strip_tilt[num_si_layers] = {0.0, 0.0, 0.0, 0.0, 0.0};  // radians, usually 1:64 tilt

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