EIC-IR/macros_fun4all/G4_IR_EIC.C

0001
0002 void IRInit() {
0003
0004 }
0005
0006 double IRSetup(PHG4Reco* g4Reco,
0007                const int absorberactive = 1,
0008                int verbosity = 0) {
0009
0010   /* Increase world size to fit extended IR magnets */
0011   g4Reco->SetWorldSizeZ(25000);
0012
0013   /* Define outer radius for all extneded IR magnets */
0014   double ir_magnet_outer_radius = 50;
0015
0016   /* Read IR configuration file- this needs to go somewhere else using parameters and a .root file to store them*/
0017   // string irfile = "../input_ir/magnet_quad_test.dat";
0018   // string irfile = "../input_ir/magnet_dipole_test.dat";
0019   // string irfile = "../input_ir/linac-ring-proton-magnets-Version3.01-21apr2016.dat";
0020   string irfile = "../input_ir/proton-magnets-250GeV-opt2.dat";
0021   // string irfile = "../input_ir/updated-magnets-2017.dat";
0022
0023   ifstream irstream(irfile.c_str());
0024
0025   while(!irstream.eof()){
0026     string str;
0027     getline(irstream, str);
0028     if(str[0] == '#') continue; //for comments
0029     stringstream ss(str);
0030
0031     string name;
0032     double center_z, center_x, center_y, aperture_radius, length, angle, B, gradient;
0033
0034     ss >> name >> center_z >> center_x >> center_y >> aperture_radius >> length >> angle >> B >> gradient;
0035
0036     if ( name == "" ) continue; //for empty lines
0037
0038     /* convert units from m to cm */
0039     center_x *= 100;
0040     center_y *= 100;
0041     center_z *= 100;
0042     aperture_radius *= 100;
0043     length *= 100;
0044
0045     /* convert angle from millirad to degrees */
0046     angle = (angle / 1000.) * (180./TMath::Pi());
0047
0048     /* Place IR component */
0049     cout << "New IR component: " << name << " at x = " << center_x << ", y = " << center_y << ", z = " << center_z << ", angle = " << angle << endl;
0050
0051     string volname = "IRMAGNET_";
0052     volname.append(name);
0053
0054     PHG4BeamlineMagnetSubsystem *ir_magnet_i = new PHG4BeamlineMagnetSubsystem(volname, 0);
0055     ir_magnet_i->set_int_param("lengthviarapidity",0);
0056     ir_magnet_i->set_double_param("length",length);
0057     ir_magnet_i->set_double_param("radius",aperture_radius);
0058     ir_magnet_i->set_double_param("thickness", ir_magnet_outer_radius - aperture_radius);
0059     ir_magnet_i->set_double_param("place_x",center_x);
0060     ir_magnet_i->set_double_param("place_y",center_y);
0061     ir_magnet_i->set_double_param("place_z",center_z);
0062     ir_magnet_i->set_double_param("rot_y",angle);
0063     ir_magnet_i->set_string_param("material","G4_Galactic"); // magnet yoke material is 'vacuum'
0064     // ir_magnet_i->set_string_param("material","G4_Fe"); // magnet yoke material is 'iron'
0065
0066     if ( B != 0 && gradient == 0.0 )
0067       {
0068         ir_magnet_i->set_string_param("magtype","dipole");
0069         ir_magnet_i->set_double_param("field_y",B);
0070       }
0071     else if ( B == 0 && gradient != 0.0 )
0072       {
0073         ir_magnet_i->set_string_param("magtype","quadrupole");
0074         ir_magnet_i->set_double_param("fieldgradient",gradient);
0075       }
0076     else if ( B == 0 && gradient == 0.0 )
0077       {
0078         ir_magnet_i->set_string_param("magtype","dipole");
0079         ir_magnet_i->set_double_param("fieldgradient",0);
0080         ir_magnet_i->set_double_param("field_y",0);
0081       }
0082     else
0083       {
0084         cout << "Error in G4_IR_EIC: Could not identify magnet type. Abort." << endl;
0085         return 1;
0086       }
0087
0088     ir_magnet_i->SetActive(false);
0089     ir_magnet_i->OverlapCheck(overlapcheck);
0090     g4Reco->registerSubsystem(ir_magnet_i);
0091   }
0092
0093   return 0;
0094 }