ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_hFarFwdBeamLine_EIC.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_hFarFwdBeamLine_EIC.C
1 #ifndef MACRO_G4HFARFWDBEAMLINE_EIC_C
2 #define MACRO_G4HFARFWDBEAMLINE_EIC_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <g4detectors/BeamLineMagnetSubsystem.h>
7 #include <g4detectors/PHG4BlockSubsystem.h>
8 #include <g4detectors/PHG4ConeSubsystem.h>
9 #include <g4detectors/PHG4CylinderSubsystem.h>
10 
11 #include <eicg4zdc/EICG4ZDCHitTree.h>
12 #include <eicg4zdc/EICG4ZDCNtuple.h>
13 #include <eicg4zdc/EICG4ZDCSubsystem.h>
14 
15 #include <eicg4b0/EICG4B0Subsystem.h>
16 #include <eicg4b0ecal/EICG4B0ECALSubsystem.h>
17 #include <eicg4rp/EICG4RPSubsystem.h>
18 
19 #include <eiceval/FarForwardEvaluator.h>
20 
21 #include <g4main/PHG4Reco.h>
22 
23 #include <TSystem.h>
24 
25 #include <fun4all/Fun4AllServer.h>
26 
27 R__LOAD_LIBRARY(libg4detectors.so)
28 
29 float PosFlip(float pos);
30 float AngleFlip(float angle);
31 float MagFieldFlip(float Bfield);
32 
33 
34 // This creates the Enable Flag to be used in the main steering macro
35 namespace Enable
36 {
37  bool HFARFWD_MAGNETS = false;
38  bool HFARFWD_VIRTUAL_DETECTORS = false;
39 
40  bool HFARFWD_PIPE = false;
41  bool HFARFWD_OVERLAPCHECK = false;
42  int HFARFWD_VERBOSITY = 0;
43 
44 // Detector configuration options
45  bool ZDC_DISABLE_BLACKHOLE = false;
46  bool B0_DISABLE_HITPLANE = false;
47  bool B0_FULLHITPLANE = false;
48  bool B0_VAR_PIPE_HOLE = false;
49  bool B0_CIRCLE_PIPE_HOLE = false;
50  bool RP_DISABLE_HITPLANE = false;
51  bool B0ECALTOWERS = true; //Set to 'false' for nice PackMan views. Set 'true' for physics studies.
52 
53  //enabled automatically in hFarFwdBeamLineInit(), unless overridden by user
54  bool HFARFWD_MAGNETS_IP6 = false;
55  bool HFARFWD_MAGNETS_IP8 = false;
56 
57  //enabled automatically in hFarFwdBeamLineInit(), unless overridden by user
58  bool HFARFWD_VIRTUAL_DETECTORS_IP6 = false;
59  bool HFARFWD_VIRTUAL_DETECTORS_IP8 = false;
60 
61  bool FFR_EVAL = false;
62 
63 } // namespace Enable
64 
65 namespace hFarFwdBeamLine
66 {
67  double starting_z = 500; //cm as center-forward interface
68  double enclosure_z_max = NAN;
69  double enclosure_r_max = NAN;
70  double enclosure_center = NAN;
71 
72  PHG4CylinderSubsystem *hFarFwdBeamLineEnclosure(nullptr);
73 
74  BeamLineMagnetSubsystem *B0Magnet = (nullptr);
75  double B0Magnet_x = NAN;
76  double B0Magnet_y = NAN;
77  double B0Magnet_z = NAN;
78 } // namespace hFarFwdBeamLine
79 
80 void hFarFwdBeamLineInit()
81 {
82  Enable::HFARFWD_MAGNETS_IP6 |= Enable::HFARFWD_MAGNETS and Enable::IP6;
83  Enable::HFARFWD_MAGNETS_IP8 |= Enable::HFARFWD_MAGNETS and Enable::IP8;
84 
85  Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 |= Enable::HFARFWD_VIRTUAL_DETECTORS and Enable::IP6;
86  Enable::HFARFWD_VIRTUAL_DETECTORS_IP8 |= Enable::HFARFWD_VIRTUAL_DETECTORS and Enable::IP8;
87 
88  if (Enable::HFARFWD_MAGNETS_IP6 && Enable::HFARFWD_MAGNETS_IP8)
89  {
90  cout << "You cannot have magnets for both IP6 and IP8 ON at the same time" << endl;
91  gSystem->Exit(1);
92  }
93 
94  if (Enable::HFARFWD_MAGNETS_IP6)
95  {
96  hFarFwdBeamLine::enclosure_z_max = 4500.;
97  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, hFarFwdBeamLine::starting_z);
98  hFarFwdBeamLine::enclosure_r_max = 200.;
99  }
100 
101  if (Enable::HFARFWD_MAGNETS_IP8)
102  {
103  hFarFwdBeamLine::enclosure_z_max = 4500.;
104  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, hFarFwdBeamLine::starting_z);
105  hFarFwdBeamLine::enclosure_r_max = 200.;
106  }
107 
108  hFarFwdBeamLine::enclosure_center = 0.5 * (hFarFwdBeamLine::starting_z + hFarFwdBeamLine::enclosure_z_max);
109 
110  BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, hFarFwdBeamLine::enclosure_z_max);
111  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, hFarFwdBeamLine::enclosure_r_max);
112 }
113 
114 void hFarFwdDefineMagnets(PHG4Reco *g4Reco)
115 {
116  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
117  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
118 
119  hFarFwdBeamLine::hFarFwdBeamLineEnclosure = new PHG4CylinderSubsystem("hFarFwdBeamLineEnclosure");
120  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("place_z", hFarFwdBeamLine::enclosure_center);
121  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("radius", 0);
122  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("thickness", hFarFwdBeamLine::enclosure_r_max); // This is intentionally made large 25cm radius
123  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_double_param("length", hFarFwdBeamLine::enclosure_z_max - hFarFwdBeamLine::starting_z);
124  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_string_param("material", "G4_Galactic");
125  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_color(.5, .5, .5, 0.2);
126  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->OverlapCheck(overlapCheck);
127  hFarFwdBeamLine::hFarFwdBeamLineEnclosure->SetActive();
128  if (verbosity) hFarFwdBeamLine::hFarFwdBeamLineEnclosure->Verbosity(verbosity);
129  g4Reco->registerSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
130 
131  if (verbosity > 0)
132  {
133  std::cout << "hFarFwdBeamLine::hFarFwdBeamLineEnclosure CanBeMotherSubsystem = " << hFarFwdBeamLine::hFarFwdBeamLineEnclosure->CanBeMotherSubsystem() << std::endl;
134  }
135 
136  string magFile;
137  if (Enable::HFARFWD_MAGNETS_IP6)
138  magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farFwdBeamLineMagnets_v2.0.dat";
139  else if (Enable::HFARFWD_MAGNETS_IP8)
140  magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip8_35mrad_h_farFwdBeamLineMagnets.dat";
141  else
142  {
143  cout << " You have to enable either the IP6 or IP8 Magnet configuration to define magnets! " << endl;
144  gSystem->Exit(1);
145  }
146 
147  // make magnet active volume if you want to study the hits
148  bool magnet_active = true;
149  int absorberactive = 0;
150 
151  // if you insert numbers it only displays those magnets, do not comment out the set declaration
152  set<int> magnetlist;
153  //magnetlist.insert(7);
154 
155  BeamLineMagnetSubsystem *bl = nullptr;
156  std::ifstream infile(magFile);
157  if (infile.is_open())
158  {
159  double biggest_z = 0.;
160  int imagnet = 0;
161  std::string line;
162  while (std::getline(infile, line))
163  {
164  if (!line.compare(0, 1, "B") ||
165  !line.compare(0, 1, "Q") ||
166  !line.compare(0, 1, "S"))
167  {
168  std::istringstream iss(line);
169  string magname;
170  double x;
171  double y;
172  double z;
173  double inner_radius_zin;
174  double inner_radius_zout;
175  double outer_magnet_diameter;
176  double length;
177  double angle;
178  double dipole_field_x;
179  double fieldgradient;
180  if (!(iss >> magname >> x >> y >> z >> inner_radius_zin >> inner_radius_zout >> outer_magnet_diameter >> length >> angle >> dipole_field_x >> fieldgradient))
181  {
182  cout << "could not decode " << line << endl;
183  gSystem->Exit(1);
184  }
185  else
186  {
187  //------------------------
188  //Select only the magnet component in the far forward region
189  if (z < 0.0)
190  continue;
191 
192  string magtype;
193  if (inner_radius_zin != inner_radius_zout)
194  {
195  cout << "inner radius at front of magnet " << inner_radius_zin
196  << " not equal radius at back of magnet " << inner_radius_zout
197  << " needs change in code (replace tube by cone for beamline)" << endl;
198  gSystem->Exit(1);
199  }
200  if (verbosity > 0)
201  {
202  cout << endl
203  << endl
204  << "\tID number " << imagnet << endl;
205  cout << "magname: " << magname << endl;
206  cout << "x: " << x << endl;
207  cout << "y: " << y << endl;
208  cout << "z: " << z << endl;
209  cout << "inner_radius_zin: " << inner_radius_zin << endl;
210  cout << "inner_radius_zout: " << inner_radius_zout << endl;
211  cout << "outer_magnet_diameter: " << outer_magnet_diameter << endl;
212  cout << "length: " << length << endl;
213  cout << "angle: " << angle << endl;
214  cout << "dipole_field_x: " << dipole_field_x << endl;
215  cout << "fieldgradient: " << fieldgradient << endl;
216  }
217  if (!magname.compare(0, 1, "B"))
218  {
219  magtype = "DIPOLE";
220  }
221  else if (!magname.compare(0, 1, "Q"))
222  {
223  magtype = "QUADRUPOLE";
224  }
225  else if (!magname.compare(0, 1, "S"))
226  {
227  magtype = "SEXTUPOLE";
228  }
229  else
230  {
231  cout << "cannot decode magnet name " << magname << endl;
232  gSystem->Exit(1);
233  }
234  // convert to our units (cm, deg)
235  x *= 100.;
236  y *= 100.;
237  z *= 100.;
238  length *= 100.;
239  inner_radius_zin *= 100.;
240  outer_magnet_diameter *= 100.;
241  angle = (angle / TMath::Pi() * 180.) / 1000.; // given in mrad
242 
243  //------------------------
244  // Linearly scaling down the magnetic field for lower energy proton
245  if( Enable::HFARFWD_ION_ENERGY != 275 ) {
246  float scaleFactor = Enable::HFARFWD_ION_ENERGY / 275. ;
247  dipole_field_x = dipole_field_x*scaleFactor;
248  fieldgradient = fieldgradient * scaleFactor;
249  }
250 
251  if (magnetlist.empty() || magnetlist.find(imagnet) != magnetlist.end())
252  {
253  bl = new BeamLineMagnetSubsystem("BEAMLINEMAGNET", imagnet);
254  bl->set_double_param("field_y", MagFieldFlip(dipole_field_x));
255  bl->set_double_param("fieldgradient", MagFieldFlip(fieldgradient));
256  bl->set_string_param("magtype", magtype);
257  bl->set_double_param("length", length);
258  bl->set_double_param("place_x", PosFlip(x)); // relative position to mother vol.
259  bl->set_double_param("place_y", y); // relative position to mother vol.
260  bl->set_double_param("place_z", z - hFarFwdBeamLine::enclosure_center); // relative position to mother vol.
261  bl->set_double_param("field_global_position_x", PosFlip(x)); // abs. position to world for field manager
262  bl->set_double_param("field_global_position_y", y); // abs. position to world for field manager
263  bl->set_double_param("field_global_position_z", z); // abs. position to world for field manager
264  bl->set_double_param("rot_y", AngleFlip(angle));
265  bl->set_double_param("field_global_rot_y", AngleFlip(angle)); // abs. rotation to world for field manager
266  bl->set_double_param("inner_radius", inner_radius_zin);
267  bl->set_double_param("outer_radius", outer_magnet_diameter / 2.);
268  bl->SetActive(magnet_active);
269  bl->SetAbsorberActive();
270  bl->BlackHole();
271  bl->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
272  if (absorberactive)
273  {
274  bl->SetAbsorberActive();
275  }
276  bl->OverlapCheck(overlapCheck);
277  bl->SuperDetector("BEAMLINEMAGNET");
278  if (verbosity) bl->Verbosity(verbosity);
279  g4Reco->registerSubsystem(bl);
280 
281  // rag the B0 magnet
282  if (imagnet == 0)
283  { //To tell the B0 Calorimeter the global coordinates of the B0 Magnet
284  hFarFwdBeamLine::B0Magnet = bl;
285  hFarFwdBeamLine::B0Magnet_x = PosFlip(x);
286  hFarFwdBeamLine::B0Magnet_y = y;
287  hFarFwdBeamLine::B0Magnet_z = z;
288  }
289  }
290  imagnet++;
291  if (fabs(z) + length > biggest_z)
292  {
293  biggest_z = fabs(z) + length;
294  }
295  }
296  }
297  }
298  infile.close();
299  }
300 }
301 
302 void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco)
303 {
304  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
305  if (Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARFWD_VIRTUAL_DETECTORS_IP8)
306  {
307  cout << "You cannot have detectors enabled for both IP6 and IP8 ON at the same time" << endl;
308  gSystem->Exit(1);
309  }
310 
311  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
312 
313  auto *detZDCsurrogate = new PHG4BlockSubsystem("zdcTruth");
314  const double detZDCsurrogate_size_z = 0.1;
315  detZDCsurrogate->SuperDetector("ZDCsurrogate");
316  detZDCsurrogate->set_double_param("place_x", PosFlip(-96.24));
317  detZDCsurrogate->set_double_param("place_y", 0);
318  detZDCsurrogate->set_double_param("place_z", 3700 - hFarFwdBeamLine::enclosure_center);
319  detZDCsurrogate->set_double_param("rot_y", AngleFlip(0.025 * TMath::RadToDeg()));
320  detZDCsurrogate->set_double_param("size_x", 60);
321  detZDCsurrogate->set_double_param("size_y", 60);
322  detZDCsurrogate->set_double_param("size_z", detZDCsurrogate_size_z);
323  detZDCsurrogate->set_string_param("material", "G4_Si");
324  detZDCsurrogate->SetActive();
325  detZDCsurrogate->set_color(1, 0, 0, 0.5);
326  detZDCsurrogate->OverlapCheck(overlapCheck);
327  if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole();
328  if (verbosity) detZDCsurrogate->Verbosity(verbosity);
329  detZDCsurrogate->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
330  g4Reco->registerSubsystem(detZDCsurrogate);
331 
332  if (Enable::ZDC_DISABLE_BLACKHOLE)
333  {
334  EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC");
335  detZDC->SetActive();
336  detZDC->set_double_param("place_z", 3700. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center);
337  detZDC->set_double_param("place_x", PosFlip(-96.24));
338  detZDC->set_double_param("rot_y", AngleFlip(0.025));
339  detZDC->OverlapCheck(overlapCheck);
340  detZDC->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
341  g4Reco->registerSubsystem(detZDC);
342  }
343 
344  const int offMomDetNr = 2;
345  const double om_zCent[offMomDetNr] = {3450, 3650};
346  const double om_xCent[offMomDetNr] = {-162, -171};
347  for (int i = 0; i < offMomDetNr; i++)
348  {
349  auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i), i);
350  detOM->SuperDetector("offMomTruth");
351  detOM->set_double_param("place_x", PosFlip(om_xCent[i]));
352  detOM->set_double_param("place_y", 0);
353  detOM->set_double_param("place_z", om_zCent[i] - hFarFwdBeamLine::enclosure_center);
354  detOM->set_double_param("rot_y", AngleFlip(0.045 * TMath::RadToDeg()));
355  detOM->set_double_param("size_x", 50);
356  detOM->set_double_param("size_y", 35);
357  detOM->set_double_param("size_z", 0.03);
358  detOM->set_string_param("material", "G4_Si");
359  detOM->SetActive();
360  if (verbosity) detOM->Verbosity(verbosity);
361  detOM->OverlapCheck(overlapCheck);
362  detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
363  g4Reco->registerSubsystem(detOM);
364  }
365 
366 
367 
368  //----------------------
369  // Roman Pots
370  //----------------------
371 
372  if( ! Enable::RP_DISABLE_HITPLANE )
373  {
374  string paramFile = string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP6.dat";
375  int Nlayers = GetParameterFromFile <int> (paramFile, "Number_layers");
376 
377  for( int layer = 0; layer < Nlayers; layer++ ) {
378  auto *detRP = new EICG4RPSubsystem( Form("rpTruth_%d", layer), layer );
379  detRP->SuperDetector( "rpTruth" );
380  detRP->set_double_param( "FFenclosure_center", hFarFwdBeamLine::enclosure_center );
381  detRP->set_int_param( "layerNumber", layer + 1 );
382  detRP->SetBeamConfig( (Enable::BEAM_COLLISION_SETTING).Data() );
383  detRP->SetIonBeamEnergy( Enable::HFARFWD_ION_ENERGY );
384  detRP->SetElectronBeamEnergy( Enable::HFARBWD_E_ENERGY );
385  detRP->SetParametersFromFile( paramFile );
386  detRP->OverlapCheck( overlapCheck );
387  detRP->SetMotherSubsystem( hFarFwdBeamLine::hFarFwdBeamLineEnclosure );
388  detRP->SetActive( true );
389  if( verbosity ) { detRP->Verbosity( verbosity ); }
390  g4Reco->registerSubsystem( detRP );
391  }
392  }
393 
394 
395  //---------------------------------
396  // B0 implementation
397  // Three choices: 1. Realistic detector; 2. Circulat plane; 3. hit plane with realistic detector goemetry
398 
399  double b0tr_z = 0; //Subsystem position relative to B0 magnet (for iterator)
400  const int b0DetNr = 4;
401  const double b0Mag_zCent = 640;
402  const double b0Mag_zLen = 120;
403  const double b0tr[4]={10,40,70,100};
404 // const double b0tr[4]={10,45,80,115}; //Tracker layers when no ECAL
405  const double b0Cu_zLen = .2; //B0 dead material length
406  const double b0Si_zLen = .1; //B0 Si length
407  const double b0Ecal_zLen = 10; //B0 Ecal length
408  double pipe_hole_r = 3.5; //detector cut off for beam pipe
409  double pipe_hole = 2.5;
410  const double cable_hole = 2.0;
411  const double cable_x = -17.0;
412  double pipe_x = -1.; //pipe hole position
413  const double d_radius = 7.0; //detector cut off Packman
414  const double b0_radius = 19.0; //outer radius of B0-detector
415  const double b0_magradius = 20.0; //inner radius of B0-magnet
416  const double spanning_angle = 240; //spanning angle Packman
417  const double b0Ecal_z = 48;//B0 ECal position (relative to the B0-magnet)
418  double start_angle = 60; //start angle Packman
419  const double cross_angle = 0.025;
420 
421  if (Enable::B0_DISABLE_HITPLANE) {
422 
423  // Choice 1 realistic detector
424 // const double b0tr[4]={10,45,80,115};
425  //const double b0tr[4]={0,30,60,90};
426  //const double b0tr[5]={0,25,50,75,100};
427  cout << "Realistic B0"<<endl;
428  for (int i = 0; i < b0DetNr; i++)
429  {
430  if (Enable::B0_VAR_PIPE_HOLE){
431  pipe_hole = b0tr[i]*cross_angle;
432  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x;
433  }
434  else if (Enable::B0_CIRCLE_PIPE_HOLE){
435  pipe_hole = 0.1;
436  pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2;
437  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
438  }
439  else {
440  pipe_hole = b0tr[b0DetNr-1]*cross_angle;
441  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
442  }
443  cout <<"Starting B0 Tracker layer "<<i+1<<endl;
444  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
445  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
446  auto *detB0 = new EICG4B0Subsystem(Form("b0Truth_%d", i), i);
447  detB0->SuperDetector(Form("b0Truth_%d", i));
448  detB0->set_double_param("place_x", 0);
449  detB0->set_double_param("place_y", 0);
450  // detB0->set_int_param("ispipe", 0); //for future pipe implementation
451  detB0->set_double_param("pipe_hole", pipe_hole);
452  detB0->set_double_param("cable_hole", cable_hole);
453  detB0->set_double_param("outer_radius", b0_radius);
454  detB0->set_double_param("d_radius", d_radius);
455  detB0->set_double_param("length", b0Si_zLen);
456  detB0->set_string_param("material", "G4_Si");
457  detB0->set_double_param("startAngle",start_angle);
458  detB0->set_double_param("spanningAngle",spanning_angle);
459  detB0->set_double_param("detid",i);
460  detB0->set_double_param("pipe_x", pipe_x);
461  detB0->set_double_param("pipe_y", 0);
462  detB0->set_double_param("pipe_z", 0);
463  detB0->set_double_param("pipe_hole_r", pipe_hole_r);
464  detB0->set_double_param("cable_x", cable_x);
465  detB0->set_double_param("cable_y", 0);
466  detB0->set_double_param("cable_z", 0);
467  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
468  detB0->SetActive(true);
469  if (verbosity)
470  detB0->Verbosity(verbosity);
471  detB0->OverlapCheck(overlapCheck);
472  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
473  g4Reco->registerSubsystem(detB0);
474 // For B0 Tracking Implementation
475  if (Enable::B0TRACKING){
476  if (B0TRACKING::FastKalmanFilter)
477  {
478  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
479  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
480  G4B0TRACKING::PositionResolution, // const float radres,
481  G4B0TRACKING::PositionResolution, // const float phires,
482  0, // const float lonres, *ignored in plane detector*
483  1, // const float eff,
484  0); // const float noise
485  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
486  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
487  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
488  G4B0TRACKING::PositionResolution, // const float radres,
489  G4B0TRACKING::PositionResolution, // const float phires,
490  0, // const float lonres, *ignored in plane detector*
491  1, // const float eff,
492  0); // const float noise
493  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
494  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
495  }
496  }
497 
498  auto *detB0e = new EICG4B0Subsystem(Form("b0Dead_%d", i), i);
499  detB0e->SuperDetector("b0Dead");
500  // detB0e->set_int_param("ispipe", 0); //for future pipe implementation
501  detB0e->set_double_param("pipe_hole", pipe_hole);
502  detB0e->set_double_param("place_x", 0);
503  detB0e->set_double_param("place_y", 0);
504  detB0e->set_double_param("d_radius", d_radius);
505  detB0e->set_double_param("pipe_x", pipe_x);
506  detB0e->set_double_param("pipe_y", 0);
507  detB0e->set_double_param("pipe_z", 0);
508  detB0e->set_double_param("pipe_hole_r", pipe_hole_r);
509  detB0e->set_double_param("cable_x", cable_x);
510  detB0e->set_double_param("cable_y", 0);
511  detB0e->set_double_param("cable_z", 0);
512  detB0e->set_double_param("outer_radius", b0_radius);
513  detB0e->set_double_param("length", b0Cu_zLen);
514  detB0e->set_string_param("material", "G4_Cu");
515  detB0e->set_double_param("detid",i);
516  detB0e->set_double_param("startAngle",start_angle);
517  detB0e->set_double_param("spanningAngle",spanning_angle);
518  detB0e->set_double_param("place_z", b0tr_z +(b0Cu_zLen+b0Si_zLen)/2) ; // relative to B0 magnet
519  detB0e->SetActive(false);
520  if (verbosity)
521  detB0e->Verbosity(verbosity);
522  detB0e->OverlapCheck(overlapCheck);
523  detB0e->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
524  g4Reco->registerSubsystem(detB0e);
525  }
526 
527  if (Enable::B0ECAL) {
528  pipe_hole = b0Mag_zLen*cross_angle;
529  pipe_x = - cross_angle*b0Mag_zCent - hFarFwdBeamLine::B0Magnet_x;
530  if (Enable::B0_CIRCLE_PIPE_HOLE){
531  pipe_hole = 0.1;
532  pipe_hole_r = pipe_hole_r + b0Mag_zLen*cross_angle/2;
533  }
534  cout <<"Starting B0 ECAL "<<endl;
535  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
536  if (Enable::B0ECALTOWERS){ //Use this option to do physics studies
537 // pipe_x=-1.25;
538 // pipe_hole=3.0;
539  cout << hFarFwdBeamLine::B0Magnet_x<<endl;
540  ostringstream mapping_b0ecal;
541  mapping_b0ecal << getenv("CALIBRATIONROOT") << "/B0Ecal/mapping/B0ECAL_mapping_v2.txt"; // Specify the mapping file for B0 ECal Towers here
542  // mapping_b0ecal << "B0ECAL_mapping_v2.txt"; // Specify the mapping file for B0 ECal Towers here
543  //cout <<"Will use B0 mapping file "<< mapping_b0ecal.str()<<endl;
544  auto *B0Ecal = new EICG4B0ECALSubsystem("B0ECAL");
545  B0Ecal->SetTowerMappingFile(mapping_b0ecal.str());
546  B0Ecal->SuperDetector("B0ECAL");
547  B0Ecal->set_double_param("pipe_hole", pipe_hole);
548  B0Ecal->set_double_param("place_x", 0);
549  B0Ecal->set_double_param("place_y", 0);
550  B0Ecal->set_double_param("place_z", b0Ecal_z);
551  B0Ecal->set_double_param("pipe_x", pipe_x);
552  B0Ecal->set_double_param("pipe_y", 0);
553  B0Ecal->set_double_param("pipe_z", 0);
554  B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r);
555  B0Ecal->set_double_param("cable_x", cable_x);
556  B0Ecal->set_double_param("cable_y", 0);
557  B0Ecal->set_double_param("cable_z", 0);
558  B0Ecal->set_double_param("length", b0Ecal_zLen);
559  B0Ecal->set_double_param("outer_radius", b0_radius);
560  B0Ecal->set_double_param("d_radius", d_radius);
561  B0Ecal->set_string_param("material", "G4_PbWO4");
562  B0Ecal->set_double_param("startAngle",start_angle);
563  B0Ecal->set_double_param("spanningAngle",spanning_angle);
564  B0Ecal->set_double_param("detid",0);
565  B0Ecal->set_double_param("global_x",hFarFwdBeamLine::B0Magnet_x);
566  B0Ecal->set_double_param("global_y",hFarFwdBeamLine::B0Magnet_y);
567  B0Ecal->set_double_param("global_z",hFarFwdBeamLine::B0Magnet_z);
568  B0Ecal->set_int_param("lightyield",1); //Note additional parameter for storing Light Yield in B0 Ecal
569  B0Ecal->SetActive(true);
570  if (verbosity)
571  B0Ecal->Verbosity(verbosity);
572  B0Ecal->OverlapCheck(overlapCheck);
573  B0Ecal->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
574  g4Reco->registerSubsystem(B0Ecal);
575  }
576  else { //Use this option to have a circular packman-shape of the B0 ECal for plots.
577  auto *B0Ecal = new EICG4B0Subsystem(Form("b0Truth_%d", 2*b0DetNr), 2*b0DetNr);
578  B0Ecal->SuperDetector("b0Truth");
579  B0Ecal->set_double_param("pipe_hole", pipe_hole);
580  B0Ecal->set_double_param("place_x", 0);
581  B0Ecal->set_double_param("place_y", 0);
582  B0Ecal->set_double_param("place_z", b0Ecal_z);
583  B0Ecal->set_double_param("pipe_x", pipe_x);
584  B0Ecal->set_double_param("pipe_y", 0);
585  B0Ecal->set_double_param("pipe_z", 0);
586  B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r);
587  B0Ecal->set_double_param("cable_x", cable_x);
588  B0Ecal->set_double_param("cable_y", 0);
589  B0Ecal->set_double_param("cable_z", 0);
590  B0Ecal->set_double_param("length", b0Ecal_zLen);
591  B0Ecal->set_double_param("outer_radius", b0_radius);
592  B0Ecal->set_double_param("d_radius", d_radius);
593  B0Ecal->set_string_param("material", "G4_PbWO4");
594  B0Ecal->set_double_param("startAngle",start_angle);
595  B0Ecal->set_double_param("spanningAngle",spanning_angle);
596  B0Ecal->set_double_param("detid",2*b0DetNr);
597  B0Ecal->SetActive(true);
598  if (verbosity)
599  B0Ecal->Verbosity(verbosity);
600  B0Ecal->OverlapCheck(overlapCheck);
601  B0Ecal->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
602  g4Reco->registerSubsystem(B0Ecal);
603  }
604 
605  auto *B0Ecale = new EICG4B0Subsystem(Form("b0Dead_%d", b0DetNr), b0DetNr); //B0 ECal dead layer is the same subsystem as other four dead layers
606  B0Ecale->SuperDetector("b0Dead");
607  // B0Ecale->set_int_param("ispipe", 0); //for future pipe implementation
608  B0Ecale->set_double_param("pipe_hole", pipe_hole);
609  B0Ecale->set_double_param("place_x", 0);
610  B0Ecale->set_double_param("place_y", 0);
611  B0Ecale->set_double_param("place_z", b0Ecal_z + (b0Ecal_zLen + b0Cu_zLen)/2);
612  B0Ecale->set_double_param("pipe_x", pipe_x);
613  B0Ecale->set_double_param("pipe_y", 0);
614  B0Ecale->set_double_param("pipe_z", 0);
615  B0Ecale->set_double_param("pipe_hole_r", pipe_hole_r);
616  B0Ecale->set_double_param("cable_x", cable_x);
617  B0Ecale->set_double_param("cable_y", 0);
618  B0Ecale->set_double_param("cable_z", 0);
619  B0Ecale->set_double_param("length", b0Cu_zLen);
620  B0Ecale->set_double_param("d_radius", d_radius);
621  B0Ecale->set_double_param("outer_radius", b0_radius);
622  B0Ecale->set_string_param("material", "G4_Cu");
623  B0Ecale->set_double_param("startAngle",start_angle);
624  B0Ecale->set_double_param("spanningAngle",spanning_angle);
625  B0Ecale->set_double_param("detid",b0DetNr+1);
626  //B0Ecale->SetActive(true);
627  B0Ecale->SetActive(false);
628  if (verbosity)
629  B0Ecale->Verbosity(verbosity);
630  B0Ecale->OverlapCheck(overlapCheck);
631  B0Ecale->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
632  g4Reco->registerSubsystem(B0Ecale);
633  }
634  } else {
635 
636  if (Enable::B0_FULLHITPLANE) {
637 
638  // Choice 2 circular hit planes
639  cout << "Circular hit planes"<<endl;
640 
641  for (int i = 0; i < b0DetNr; i++)
642  {
643  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
644  auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), i);
645  detB0->SuperDetector("b0Truth");
646  //detB0->SuperDetector(Form("b0Truth_%d", i));
647  detB0->set_double_param("radius", 0);
648  detB0->set_double_param("thickness", 20);
649  detB0->set_double_param("length", 0.1);
650  detB0->set_string_param("material", "G4_Si");
651  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
652  detB0->SetActive(true);
653  if (verbosity) detB0->Verbosity(verbosity);
654  detB0->OverlapCheck(overlapCheck);
655 
656  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
657 
658  g4Reco->registerSubsystem(detB0);
659  if (Enable::B0TRACKING){
660  if (B0TRACKING::FastKalmanFilter)
661  {
662  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
663  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
664  G4B0TRACKING::PositionResolution, // const float radres,
665  G4B0TRACKING::PositionResolution, // const float phires,
666  0, // const float lonres, *ignored in plane detector*
667  1, // const float eff,
668  0); // const float noise
669  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
670  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
671  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
672  G4B0TRACKING::PositionResolution, // const float radres,
673  G4B0TRACKING::PositionResolution, // const float phires,
674  0, // const float lonres, *ignored in plane detector*
675  1, // const float eff,
676  0); // const float noise
677  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
678  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
679  }
680  }
681 
682  }
683 
684  } else {
685 
688  cout << "Realistic hit planes"<<endl;
689 
690  for (int i = 0; i < b0DetNr; i++) {
691  if (Enable::B0_VAR_PIPE_HOLE){
692  pipe_hole = b0tr[i]*cross_angle;
693  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x;
694  }
695  else if (Enable::B0_CIRCLE_PIPE_HOLE){
696  pipe_hole = 0.1;
697  pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2;
698  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
699  }
700  else {
701  pipe_hole = b0tr[b0DetNr-1]*cross_angle;
702  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
703  }
704  cout <<"Starting B0 Tracker layer "<<i+1<<endl;
705  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
706  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
707  auto *detB0 = new EICG4B0Subsystem(Form("b0Truth_%d", i), i);
708  detB0->SuperDetector(Form("b0Truth_%d", i));
709  detB0->set_double_param("place_x", 0);
710  detB0->set_double_param("place_y", 0);
711  // detB0->set_int_param("ispipe", 0); //for future pipe implementation
712  detB0->set_double_param("pipe_hole", pipe_hole);
713  detB0->set_double_param("cable_hole", cable_hole);
714  detB0->set_double_param("outer_radius", b0_radius);
715  detB0->set_double_param("d_radius", d_radius);
716  detB0->set_double_param("length", b0Si_zLen);
717  detB0->set_string_param("material", "G4_Si");
718  detB0->set_double_param("startAngle",start_angle);
719  detB0->set_double_param("spanningAngle",spanning_angle);
720  detB0->set_double_param("detid",i);
721  detB0->set_double_param("pipe_x", pipe_x);
722  detB0->set_double_param("pipe_y", 0);
723  detB0->set_double_param("pipe_z", 0);
724  detB0->set_double_param("pipe_hole_r", pipe_hole_r);
725  detB0->set_double_param("cable_x", cable_x);
726  detB0->set_double_param("cable_y", 0);
727  detB0->set_double_param("cable_z", 0);
728  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
729  detB0->SetActive(true);
730  if (verbosity)
731  detB0->Verbosity(verbosity);
732  detB0->OverlapCheck(overlapCheck);
733  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
734  g4Reco->registerSubsystem(detB0);
735  if (Enable::B0TRACKING){
736  if (B0TRACKING::FastKalmanFilter)
737  {
738  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
739  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
740  G4B0TRACKING::PositionResolution, // const float radres,
741  G4B0TRACKING::PositionResolution, // const float phires,
742  0, // const float lonres, *ignored in plane detector*
743  1, // const float eff,
744  0); // const float noise
745  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
746  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
747  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
748  G4B0TRACKING::PositionResolution, // const float radres,
749  G4B0TRACKING::PositionResolution, // const float phires,
750  0, // const float lonres, *ignored in plane detector*
751  1, // const float eff,
752  0); // const float noise
753  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
754  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
755  }
756  }
757  }
758 
759  }
760  }
761 
762 }
763 
764 void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco)
765 {
766 //--------------------------------------------------------
767 // The IP8 detector position is implemented by Wenliang Li (billlee@jlab.org)
768 // on July 07, 2021
769 // Reference of this implementation: https://indico.bnl.gov/event/10974/contributions/51160/
770 
771  bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARFWD_OVERLAPCHECK;
772  if (Enable::HFARFWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARFWD_VIRTUAL_DETECTORS_IP8)
773  {
774  cout << "You cannot have detectors enabled for both IP6 and IP8 ON at the same time" << endl;
775  gSystem->Exit(1);
776  }
777 
778  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
779 
780  const int offMomDetNr = 2;
781  const double om_xCent[offMomDetNr] = {46, 49};
782  const double om_zCent[offMomDetNr] = {3250, 3450};
783 
784  for (int i = 0; i < offMomDetNr; i++)
785  {
786  auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i), i);
787  detOM->SuperDetector("offMomTruth");
788  detOM->set_double_param("place_x", PosFlip(om_xCent[i]));
789  detOM->set_double_param("place_y", 0);
790  detOM->set_double_param("place_z", PosFlip(om_zCent[i] - hFarFwdBeamLine::enclosure_center));
791  detOM->set_double_param("rot_y", AngleFlip(-0.045 * TMath::RadToDeg()));
792  detOM->set_double_param("size_x", 40); // Original design specification
793  detOM->set_double_param("size_y", 35); // Original design specification
794  detOM->set_double_param("size_z", 0.03);
795  detOM->set_string_param("material", "G4_Si");
796  detOM->OverlapCheck(overlapCheck);
797  detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
798  detOM->SetActive();
799  detOM->set_color(0, 0, 1, 0.5);
800  if (verbosity) detOM->Verbosity(verbosity);
801  g4Reco->registerSubsystem(detOM);
802  }
803 
804  auto *detZDCsurrogate = new PHG4BlockSubsystem("zdcTruth");
805  const double detZDCsurrogate_size_z = 0.1;
806  detZDCsurrogate->SuperDetector("ZDCsurrogate");
807  detZDCsurrogate->set_double_param("place_x", PosFlip(120));
808  detZDCsurrogate->set_double_param("place_y", 0);
809  detZDCsurrogate->set_double_param("place_z", 3350 - hFarFwdBeamLine::enclosure_center);
810  detZDCsurrogate->set_double_param("rot_y", AngleFlip(-0.035 * TMath::RadToDeg()));
811  detZDCsurrogate->set_double_param("size_x", 60);
812  detZDCsurrogate->set_double_param("size_y", 60);
813  detZDCsurrogate->set_double_param("size_z", detZDCsurrogate_size_z);
814  detZDCsurrogate->set_string_param("material", "G4_Si");
815  detZDCsurrogate->SetActive();
816  detZDCsurrogate->OverlapCheck(overlapCheck);
817  detZDCsurrogate->set_color(1, 0, 0, 0.5);
818  if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole();
819  if (verbosity) detZDCsurrogate->Verbosity(verbosity);
820  detZDCsurrogate->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
821  g4Reco->registerSubsystem(detZDCsurrogate);
822 
823  if (Enable::ZDC_DISABLE_BLACKHOLE)
824  {
825 
826  EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC");
827  detZDC->SetActive();
828  detZDC->set_double_param("place_z", 3350. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center);
829  detZDC->set_double_param("place_x", PosFlip(120.0));
830  detZDC->set_double_param("rot_y", AngleFlip(-0.035));
831  detZDC->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure);
832  detZDC->OverlapCheck(overlapCheck);
833  g4Reco->registerSubsystem(detZDC);
834 
835  }
836 
837 
838  //----------------------
839  // Roman Pots: Both sets before and near the secondary focus
840  //----------------------
841 
842  if( ! Enable::RP_DISABLE_HITPLANE )
843  {
844  string paramFile = string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP8.dat";
845  int Nlayers = GetParameterFromFile <int> (paramFile, "Number_layers");
846 
847  for( int layer = 0; layer < Nlayers; layer++ ) {
848  auto *detRP = new EICG4RPSubsystem( Form("rpTruth_%d", layer), layer );
849  detRP->SuperDetector( "rpTruth" );
850  detRP->set_double_param( "FFenclosure_center", hFarFwdBeamLine::enclosure_center );
851  detRP->set_int_param( "layerNumber", layer + 1 );
852  detRP->SetBeamConfig( (Enable::BEAM_COLLISION_SETTING).Data() );
853  detRP->SetIonBeamEnergy( Enable::HFARFWD_ION_ENERGY );
854  detRP->SetElectronBeamEnergy( Enable::HFARBWD_E_ENERGY );
855  detRP->SetParametersFromFile( paramFile );
856  detRP->OverlapCheck( overlapCheck );
857  detRP->SetMotherSubsystem( hFarFwdBeamLine::hFarFwdBeamLineEnclosure );
858  detRP->SetActive( true );
859  if( verbosity ) { detRP->Verbosity( verbosity ); }
860  g4Reco->registerSubsystem( detRP );
861  }
862  }
863 
864  if (verbosity > 0)
865  {
866  std::cout << "B0Magnet can be mother = " << hFarFwdBeamLine::B0Magnet->CanBeMotherSubsystem() << std::endl;
867  }
868 
869 /* const int b0DetNr = 4;
870  const double b0Mag_zCent = 610;
871  const double b0Mag_zLen = 120;
872  for (int i = 0; i < b0DetNr; i++)
873  {
874  auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), i);
875  detB0->SuperDetector("b0Truth");
876  detB0->set_double_param("radius", 0);
877  detB0->set_double_param("thickness", 20);
878  detB0->set_double_param("length", 0.1);
879  detB0->set_string_param("material", "G4_Si");
880  detB0->set_double_param("place_y", 0);
881  detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2));
882  detB0->OverlapCheck(overlapCheck);
883  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
884  detB0->SetActive(true);
885  if (verbosity)
886  detB0->Verbosity(verbosity);
887  g4Reco->registerSubsystem(detB0);
888  }*/
889  //---------------------------------
890  // B0 implementation
891  // Three choices: 1. Realistic detector; 2. Circulat plane; 3. hit plane with realistic detector goemetry
892  double b0tr_z = 0; //Subsystem position relative to B0 magnet (for iterator)
893  const int b0DetNr = 4;
894  const double b0Mag_zCent = 610;
895  const double b0Mag_zLen = 120;
896  const double b0tr[4]={10,40,70,100};
897  const double b0Cu_zLen = .2; //B0 dead material length
898  const double b0Si_zLen = .1; //B0 Si length
899  const double b0Ecal_zLen = 10; //B0 Ecal length
900  double pipe_hole_r = 3.5; //detector cut off for beam pipe
901  double pipe_hole = 2.5;
902  const double cable_hole = 2.0;
903  const double cable_x = 21.5;
904  double pipe_x = -1.; //pipe hole position
905  const double d_radius = 7.0; //detector cut off Packman
906  const double b0_radius = 23.5; //outer radius of B0-detector
907  const double b0_magradius = 24.5; //inner radius of B0-magnet
908  const double spanning_angle = 240; //spanning angle Packman
909  const double b0Ecal_z = 48;//B0 ECal position (relative to the B0-magnet)
910  double start_angle = -120; //start angle Packman
911  const double cross_angle = 0.035;
912 
913  if (Enable::B0_DISABLE_HITPLANE) {
914 
915  // Choice 1 realistic detector
916 // const double b0tr[4]={10,45,80,115};
917  //const double b0tr[4]={0,30,60,90};
918  //const double b0tr[5]={0,25,50,75,100};
919  cout << "Realistic B0"<<endl;
920  for (int i = 0; i < b0DetNr; i++)
921  {
922  if (Enable::B0_VAR_PIPE_HOLE){
923  pipe_hole = b0tr[i]*cross_angle;
924  pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x;
925  }
926  else if (Enable::B0_CIRCLE_PIPE_HOLE){
927  pipe_hole = 0.1;
928  pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2;
929  pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
930  }
931  else {
932  pipe_hole = b0tr[b0DetNr-1]*cross_angle;
933  pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
934  }
935  cout <<"Starting B0 Tracker layer "<<i+1<<endl;
936  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
937  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
938  auto *detB0 = new EICG4B0Subsystem(Form("b0Truth_%d", i), i);
939  detB0->SuperDetector(Form("b0Truth_%d", i));
940  detB0->set_double_param("place_x", 0);
941  detB0->set_double_param("place_y", 0);
942  // detB0->set_int_param("ispipe", 0); //for future pipe implementation
943  detB0->set_double_param("pipe_hole", pipe_hole);
944  detB0->set_double_param("cable_hole", cable_hole);
945  detB0->set_double_param("outer_radius", b0_radius);
946  detB0->set_double_param("d_radius", d_radius);
947  detB0->set_double_param("length", b0Si_zLen);
948  detB0->set_string_param("material", "G4_Si");
949  detB0->set_double_param("startAngle",start_angle);
950  detB0->set_double_param("spanningAngle",spanning_angle);
951  detB0->set_double_param("detid",i);
952  detB0->set_double_param("pipe_x", pipe_x);
953  detB0->set_double_param("pipe_y", 0);
954  detB0->set_double_param("pipe_z", 0);
955  detB0->set_double_param("pipe_hole_r", pipe_hole_r);
956  detB0->set_double_param("cable_x", cable_x);
957  detB0->set_double_param("cable_y", 0);
958  detB0->set_double_param("cable_z", 0);
959  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
960  detB0->SetActive(true);
961  if (verbosity)
962  detB0->Verbosity(verbosity);
963  detB0->OverlapCheck(overlapCheck);
964  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
965  g4Reco->registerSubsystem(detB0);
966 // For B0 Tracking Implementation
967  if (Enable::B0TRACKING){
968  if (B0TRACKING::FastKalmanFilter)
969  {
970  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
971  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
972  G4B0TRACKING::PositionResolution, // const float radres,
973  G4B0TRACKING::PositionResolution, // const float phires,
974  0, // const float lonres, *ignored in plane detector*
975  1, // const float eff,
976  0); // const float noise
977  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
978  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
979  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
980  G4B0TRACKING::PositionResolution, // const float radres,
981  G4B0TRACKING::PositionResolution, // const float phires,
982  0, // const float lonres, *ignored in plane detector*
983  1, // const float eff,
984  0); // const float noise
985  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
986  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
987  }
988  }
989 
990  auto *detB0e = new EICG4B0Subsystem(Form("b0Dead_%d", i), i);
991  detB0e->SuperDetector("b0Dead");
992  // detB0e->set_int_param("ispipe", 0); //for future pipe implementation
993  detB0e->set_double_param("pipe_hole", pipe_hole);
994  detB0e->set_double_param("place_x", 0);
995  detB0e->set_double_param("place_y", 0);
996  detB0e->set_double_param("d_radius", d_radius);
997  detB0e->set_double_param("pipe_x", pipe_x);
998  detB0e->set_double_param("pipe_y", 0);
999  detB0e->set_double_param("pipe_z", 0);
1000  detB0e->set_double_param("pipe_hole_r", pipe_hole_r);
1001  detB0e->set_double_param("cable_x", cable_x);
1002  detB0e->set_double_param("cable_y", 0);
1003  detB0e->set_double_param("cable_z", 0);
1004  detB0e->set_double_param("outer_radius", b0_radius);
1005  detB0e->set_double_param("length", b0Cu_zLen);
1006  detB0e->set_string_param("material", "G4_Cu");
1007  detB0e->set_double_param("detid",i);
1008  detB0e->set_double_param("startAngle",start_angle);
1009  detB0e->set_double_param("spanningAngle",spanning_angle);
1010  detB0e->set_double_param("place_z", b0tr_z +(b0Cu_zLen+b0Si_zLen)/2) ; // relative to B0 magnet
1011  detB0e->SetActive(false);
1012  if (verbosity)
1013  detB0e->Verbosity(verbosity);
1014  detB0e->OverlapCheck(overlapCheck);
1015  detB0e->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1016  g4Reco->registerSubsystem(detB0e);
1017  }
1018 
1019  if (Enable::B0ECAL) {
1020  pipe_hole = b0Mag_zLen*cross_angle;
1021  pipe_x = cross_angle*b0Mag_zCent - hFarFwdBeamLine::B0Magnet_x;
1022  if (Enable::B0_CIRCLE_PIPE_HOLE){
1023  pipe_hole = 0.1;
1024  pipe_hole_r = pipe_hole_r + b0Mag_zLen*cross_angle/2;
1025  }
1026  cout <<"Starting B0 ECAL "<<endl;
1027  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
1028  if (Enable::B0ECALTOWERS){ //Use this option to do physics studies
1029 // pipe_x=-1.25;
1030 // pipe_hole=3.0;
1031  cout << hFarFwdBeamLine::B0Magnet_x<<endl;
1032  ostringstream mapping_b0ecal;
1033  mapping_b0ecal << getenv("CALIBRATIONROOT") << "/B0Ecal/mapping/B0ECAL_mapping_ip8_v1.txt"; // Specify the mapping file for B0 ECal Towers here
1034 // mapping_b0ecal << "B0ECAL_mapping_ip8_v1.txt"; // Specify the mapping file for B0 ECal Towers here
1035  //cout <<"Will use B0 mapping file "<< mapping_b0ecal.str()<<endl;
1036  auto *B0Ecal = new EICG4B0ECALSubsystem("B0ECAL");
1037  B0Ecal->SetTowerMappingFile(mapping_b0ecal.str());
1038  B0Ecal->SuperDetector("B0ECAL");
1039  B0Ecal->set_double_param("pipe_hole", pipe_hole);
1040  B0Ecal->set_double_param("place_x", 0);
1041  B0Ecal->set_double_param("place_y", 0);
1042  B0Ecal->set_double_param("place_z", b0Ecal_z);
1043  B0Ecal->set_double_param("pipe_x", pipe_x);
1044  B0Ecal->set_double_param("pipe_y", 0);
1045  B0Ecal->set_double_param("pipe_z", 0);
1046  B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r);
1047  B0Ecal->set_double_param("cable_x", cable_x);
1048  B0Ecal->set_double_param("cable_y", 0);
1049  B0Ecal->set_double_param("cable_z", 0);
1050  B0Ecal->set_double_param("length", b0Ecal_zLen);
1051  B0Ecal->set_double_param("outer_radius", b0_radius);
1052  B0Ecal->set_double_param("d_radius", d_radius);
1053  B0Ecal->set_string_param("material", "G4_PbWO4");
1054  B0Ecal->set_double_param("startAngle",start_angle);
1055  B0Ecal->set_double_param("spanningAngle",spanning_angle);
1056  B0Ecal->set_double_param("detid",0);
1057  B0Ecal->set_double_param("global_x",hFarFwdBeamLine::B0Magnet_x);
1058  B0Ecal->set_double_param("global_y",hFarFwdBeamLine::B0Magnet_y);
1059  B0Ecal->set_double_param("global_z",hFarFwdBeamLine::B0Magnet_z);
1060  B0Ecal->set_int_param("lightyield",1); //Note additional parameter for storing Light Yield in B0 Ecal
1061  B0Ecal->SetActive(true);
1062  if (verbosity)
1063  B0Ecal->Verbosity(verbosity);
1064  B0Ecal->OverlapCheck(overlapCheck);
1065  B0Ecal->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1066  g4Reco->registerSubsystem(B0Ecal);
1067  }
1068  else { //Use this option to have a circular packman-shape of the B0 ECal for plots.
1069  auto *B0Ecal = new EICG4B0Subsystem(Form("b0Truth_%d", 2*b0DetNr), 2*b0DetNr);
1070  B0Ecal->SuperDetector("b0Truth");
1071  B0Ecal->set_double_param("pipe_hole", pipe_hole);
1072  B0Ecal->set_double_param("place_x", 0);
1073  B0Ecal->set_double_param("place_y", 0);
1074  B0Ecal->set_double_param("place_z", b0Ecal_z);
1075  B0Ecal->set_double_param("pipe_x", pipe_x);
1076  B0Ecal->set_double_param("pipe_y", 0);
1077  B0Ecal->set_double_param("pipe_z", 0);
1078  B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r);
1079  B0Ecal->set_double_param("cable_x", cable_x);
1080  B0Ecal->set_double_param("cable_y", 0);
1081  B0Ecal->set_double_param("cable_z", 0);
1082  B0Ecal->set_double_param("length", b0Ecal_zLen);
1083  B0Ecal->set_double_param("outer_radius", b0_radius);
1084  B0Ecal->set_double_param("d_radius", d_radius);
1085  B0Ecal->set_string_param("material", "G4_PbWO4");
1086  B0Ecal->set_double_param("startAngle",start_angle);
1087  B0Ecal->set_double_param("spanningAngle",spanning_angle);
1088  B0Ecal->set_double_param("detid",2*b0DetNr);
1089  B0Ecal->SetActive(true);
1090  if (verbosity)
1091  B0Ecal->Verbosity(verbosity);
1092  B0Ecal->OverlapCheck(overlapCheck);
1093  B0Ecal->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1094  g4Reco->registerSubsystem(B0Ecal);
1095  }
1096 
1097  auto *B0Ecale = new EICG4B0Subsystem(Form("b0Dead_%d", b0DetNr), b0DetNr); //B0 ECal dead layer is the same subsystem as other four dead layers
1098  B0Ecale->SuperDetector("b0Dead");
1099  // B0Ecale->set_int_param("ispipe", 0); //for future pipe implementation
1100  B0Ecale->set_double_param("pipe_hole", pipe_hole);
1101  B0Ecale->set_double_param("place_x", 0);
1102  B0Ecale->set_double_param("place_y", 0);
1103  B0Ecale->set_double_param("place_z", b0Ecal_z + (b0Ecal_zLen + b0Cu_zLen)/2);
1104  B0Ecale->set_double_param("pipe_x", pipe_x);
1105  B0Ecale->set_double_param("pipe_y", 0);
1106  B0Ecale->set_double_param("pipe_z", 0);
1107  B0Ecale->set_double_param("pipe_hole_r", pipe_hole_r);
1108  B0Ecale->set_double_param("cable_x", cable_x);
1109  B0Ecale->set_double_param("cable_y", 0);
1110  B0Ecale->set_double_param("cable_z", 0);
1111  B0Ecale->set_double_param("length", b0Cu_zLen);
1112  B0Ecale->set_double_param("d_radius", d_radius);
1113  B0Ecale->set_double_param("outer_radius", b0_radius);
1114  B0Ecale->set_string_param("material", "G4_Cu");
1115  B0Ecale->set_double_param("startAngle",start_angle);
1116  B0Ecale->set_double_param("spanningAngle",spanning_angle);
1117  B0Ecale->set_double_param("detid",b0DetNr+1);
1118  //B0Ecale->SetActive(true);
1119  B0Ecale->SetActive(false);
1120  if (verbosity)
1121  B0Ecale->Verbosity(verbosity);
1122  B0Ecale->OverlapCheck(overlapCheck);
1123  B0Ecale->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1124  g4Reco->registerSubsystem(B0Ecale);
1125  }
1126  } else {
1127 
1128  if (Enable::B0_FULLHITPLANE) {
1129 
1130  // Choice 2 circular hit planes
1131  cout << "Circular hit planes"<<endl;
1132 
1133  for (int i = 0; i < b0DetNr; i++)
1134  {
1135  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
1136  auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), i);
1137  detB0->SuperDetector("b0Truth");
1138  //detB0->SuperDetector(Form("b0Truth_%d", i));
1139  detB0->set_double_param("radius", 0);
1140  detB0->set_double_param("thickness", 20);
1141  detB0->set_double_param("length", 0.1);
1142  detB0->set_string_param("material", "G4_Si");
1143  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
1144  detB0->SetActive(true);
1145  if (verbosity) detB0->Verbosity(verbosity);
1146  detB0->OverlapCheck(overlapCheck);
1147 
1148  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1149 
1150  g4Reco->registerSubsystem(detB0);
1151  if (Enable::B0TRACKING){
1152  if (B0TRACKING::FastKalmanFilter)
1153  {
1154  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
1155  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
1156  G4B0TRACKING::PositionResolution, // const float radres,
1157  G4B0TRACKING::PositionResolution, // const float phires,
1158  0, // const float lonres, *ignored in plane detector*
1159  1, // const float eff,
1160  0); // const float noise
1161  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
1162  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
1163  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
1164  G4B0TRACKING::PositionResolution, // const float radres,
1165  G4B0TRACKING::PositionResolution, // const float phires,
1166  0, // const float lonres, *ignored in plane detector*
1167  1, // const float eff,
1168  0); // const float noise
1169  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
1170  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
1171  }
1172  }
1173 
1174  }
1175 
1176  } else {
1177 
1180  cout << "Realistic hit planes"<<endl;
1181 
1182  for (int i = 0; i < b0DetNr; i++) {
1183  if (Enable::B0_VAR_PIPE_HOLE){
1184  pipe_hole = b0tr[i]*cross_angle;
1185  pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x;
1186  }
1187  else if (Enable::B0_CIRCLE_PIPE_HOLE){
1188  pipe_hole = 0.1;
1189  pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2;
1190  pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
1191  }
1192  else {
1193  pipe_hole = b0tr[b0DetNr-1]*cross_angle;
1194  pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x;
1195  }
1196  cout <<"Starting B0 Tracker layer "<<i+1<<endl;
1197  cout <<"Pipe Hole: "<< pipe_hole<<"\t"<<pipe_x<<endl;
1198  b0tr_z = b0tr[i] - b0Mag_zLen / 2;
1199  auto *detB0 = new EICG4B0Subsystem(Form("b0Truth_%d", i), i);
1200  detB0->SuperDetector(Form("b0Truth_%d", i));
1201  detB0->set_double_param("place_x", 0);
1202  detB0->set_double_param("place_y", 0);
1203  // detB0->set_int_param("ispipe", 0); //for future pipe implementation
1204  detB0->set_double_param("pipe_hole", pipe_hole);
1205  detB0->set_double_param("cable_hole", cable_hole);
1206  detB0->set_double_param("outer_radius", b0_radius);
1207  detB0->set_double_param("d_radius", d_radius);
1208  detB0->set_double_param("length", b0Si_zLen);
1209  detB0->set_string_param("material", "G4_Si");
1210  detB0->set_double_param("startAngle",start_angle);
1211  detB0->set_double_param("spanningAngle",spanning_angle);
1212  detB0->set_double_param("detid",i);
1213  detB0->set_double_param("pipe_x", pipe_x);
1214  detB0->set_double_param("pipe_y", 0);
1215  detB0->set_double_param("pipe_z", 0);
1216  detB0->set_double_param("pipe_hole_r", pipe_hole_r);
1217  detB0->set_double_param("cable_x", cable_x);
1218  detB0->set_double_param("cable_y", 0);
1219  detB0->set_double_param("cable_z", 0);
1220  detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet
1221  detB0->SetActive(true);
1222  if (verbosity)
1223  detB0->Verbosity(verbosity);
1224  detB0->OverlapCheck(overlapCheck);
1225  detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet);
1226  g4Reco->registerSubsystem(detB0);
1227  if (Enable::B0TRACKING){
1228  if (B0TRACKING::FastKalmanFilter)
1229  {
1230  B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
1231  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
1232  G4B0TRACKING::PositionResolution, // const float radres,
1233  G4B0TRACKING::PositionResolution, // const float phires,
1234  0, // const float lonres, *ignored in plane detector*
1235  1, // const float eff,
1236  0); // const float noise
1237  B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
1238  B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames,
1239  B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
1240  G4B0TRACKING::PositionResolution, // const float radres,
1241  G4B0TRACKING::PositionResolution, // const float phires,
1242  0, // const float lonres, *ignored in plane detector*
1243  1, // const float eff,
1244  0); // const float noise
1245  B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z);
1246  B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i));
1247  }
1248  }
1249  }
1250 
1251  }
1252  }
1253 }
1254 
1255 
1256 void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco)
1257 {
1258  int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY);
1259  //exit window
1260  PHG4CylinderSubsystem *exitWin = new PHG4CylinderSubsystem("exitWin", 0);
1261  exitWin->set_double_param("radius", 3.2);
1262  exitWin->set_double_param("thickness", 11.8);
1263  exitWin->set_double_param("length", 0.15);
1264  exitWin->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg()));
1265  exitWin->set_string_param("material", "G4_STAINLESS-STEEL");
1266  exitWin->set_double_param("place_x", PosFlip(12.5));
1267  exitWin->set_double_param("place_y", 0);
1268  exitWin->set_double_param("place_z", 500);
1269  exitWin->SetActive(false);
1270  g4Reco->registerSubsystem(exitWin);
1271 
1272  //B0 magnet pipe
1273  PHG4CylinderSubsystem *pipeB0 = new PHG4CylinderSubsystem("beamPipeB0", 0);
1274  pipeB0->set_double_param("radius", 2.8);
1275  pipeB0->set_double_param("thickness", 0.25);
1276  pipeB0->set_double_param("length", 195);
1277  pipeB0->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg()));
1278  pipeB0->set_string_param("material", "G4_Al");
1279  pipeB0->set_double_param("place_x", PosFlip(14.748));
1280  pipeB0->set_double_param("place_y", 0);
1281  pipeB0->set_double_param("place_z", 590);
1282  pipeB0->SetActive(false);
1283  g4Reco->registerSubsystem(pipeB0);
1284 
1285  //Quad pipes
1286  const int nSecQ = 5; //B0apf, Q1apf, Q1bpf, Q2pf, B1pf
1287  const string nm[nSecQ] = {"B0apf", "Q1apf", "Q1bpf", "Q2pf", "B1pf"};
1288  const double qlen[nSecQ] = {160, 150, 220, 440, 330};
1289  const double qir[nSecQ] = {4, 5.1, 7, 12, 12.2};
1290  const double qor[nSecQ] = {4.1, 5.2, 7.2, 12.2, 12.4};
1291  const double qrot[nSecQ] = {25, 19.5, 15, 15, 34}; //mrad
1292  const double qxC[nSecQ] = {19.8, 24.47, 30.05, 39.5, 48};
1293  const double qyC[nSecQ] = {0, 0, 0, 0, 0};
1294  const double qzC[nSecQ] = {770, 922.8, 1106.3, 1416.7, 1806.7};
1295  for (int i = 0; i < nSecQ; i++)
1296  {
1297  PHG4CylinderSubsystem *pipe = new PHG4CylinderSubsystem(Form("beamPipe%s", nm[i].c_str()), 0);
1298  pipe->set_double_param("radius", qir[i]);
1299  pipe->set_double_param("thickness", qor[i] - qir[i]);
1300  pipe->set_double_param("length", qlen[i]);
1301  pipe->set_double_param("rot_y", AngleFlip(-qrot[i] / 1000 * TMath::RadToDeg()));
1302  pipe->set_string_param("material", "G4_Al");
1303  pipe->set_double_param("place_x", PosFlip(qxC[i]));
1304  pipe->set_double_param("place_y", qyC[i]);
1305  pipe->set_double_param("place_z", qzC[i]);
1306  pipe->SetActive(false);
1307  g4Reco->registerSubsystem(pipe);
1308  }
1309 
1310  //Electron pipe
1311  PHG4CylinderSubsystem *pipeElectron = new PHG4CylinderSubsystem("beamPipeElectron", 0);
1312  pipeElectron->set_double_param("radius", 1);
1313  pipeElectron->set_double_param("thickness", 1);
1314  pipeElectron->set_double_param("length", 3000);
1315  pipeElectron->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg()));
1316  pipeElectron->set_string_param("material", "G4_Al");
1317  pipeElectron->set_double_param("place_x", PosFlip(0));
1318  pipeElectron->set_double_param("place_y", 0);
1319  pipeElectron->set_double_param("place_z", 2000);
1320  pipeElectron->SetActive(false);
1321  //g4Reco->registerSubsystem(pipeElectron);
1322 
1323  //ZDC pipe
1324  PHG4CylinderSubsystem *pipeZDC = new PHG4CylinderSubsystem("beamPipeZDC", 0);
1325  pipeZDC->set_double_param("radius", 16.5);
1326  pipeZDC->set_double_param("thickness", 0.1);
1327  pipeZDC->set_double_param("length", 170);
1328  pipeZDC->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg()));
1329  pipeZDC->set_string_param("material", "G4_Al");
1330  pipeZDC->set_double_param("place_x", PosFlip(59));
1331  pipeZDC->set_double_param("place_y", 0);
1332  pipeZDC->set_double_param("place_z", 2041.59);
1333  pipeZDC->SetActive(false);
1334  g4Reco->registerSubsystem(pipeZDC);
1335 
1336  //Roman Pot pipe
1337  const int nSec = 2;
1338  const double len[nSec] = {850, 1150};
1339  const double ir1[nSec] = {17, 17};
1340  const double or1[nSec] = {17.1, 17.1};
1341  const double ir2[nSec] = {17, 7};
1342  const double or2[nSec] = {17.1, 7.1};
1343  const double xC[nSec] = {83, 130};
1344  const double yC[nSec] = {0, 0};
1345  const double zC[nSec] = {2550, 3550};
1346  for (int i = 0; i < nSec; i++)
1347  {
1348  PHG4ConeSubsystem *pipe = new PHG4ConeSubsystem(Form("beamPipeRP%d", i), 0);
1349  pipe->set_string_param("material", "G4_STAINLESS-STEEL");
1350  pipe->set_double_param("place_x", PosFlip(xC[i]));
1351  pipe->set_double_param("place_y", yC[i]);
1352  pipe->set_double_param("place_z", zC[i]);
1353  pipe->set_double_param("length", len[i] / 2);
1354  pipe->set_double_param("rmin1", ir1[i]);
1355  pipe->set_double_param("rmin2", ir2[i]);
1356  pipe->set_double_param("rmax1", or1[i]);
1357  pipe->set_double_param("rmax2", or2[i]);
1358  pipe->set_double_param("rot_y", AngleFlip(-0.047 * TMath::RadToDeg()));
1359  g4Reco->registerSubsystem(pipe);
1360  }
1361 }
1362 
1363 float PosFlip(float pos) {
1364  if(Enable::HFARFWD_MAGNETS_IP6) {
1365  return pos;
1366  } else {
1367  return pos;
1368  }
1369 }
1370 
1371 float AngleFlip(float angle){
1372  if(Enable::HFARFWD_MAGNETS_IP6) {
1373  return angle;
1374  } else {
1375  return angle;
1376  }
1377 }
1378 
1379 float MagFieldFlip(float Bfield){
1380  if(Enable::HFARFWD_MAGNETS_IP6) {
1381  return Bfield;
1382  } else {
1383  return Bfield;
1384  }
1385 }
1386 
1387 
1388 //------------------------------------------
1389 
1390 void FFR_Eval(const std::string &outputfile)
1391 {
1392 
1393  string ip_str;
1394 
1395 
1396  if(Enable::IP6) {
1397  ip_str = "IP6";
1398  } else {
1399  ip_str = "IP8";
1400  }
1401 
1402  int verbosity = std::max(Enable::VERBOSITY, Enable::EEMC_VERBOSITY);
1403 
1404  Fun4AllServer *se = Fun4AllServer::instance();
1405 
1406  FarForwardEvaluator *eval = new FarForwardEvaluator("FARFORWARDEVALUATOR", "FFR", outputfile.c_str(), ip_str);
1407 
1408  eval->Verbosity(verbosity);
1409  se->registerSubsystem(eval);
1410 
1411  return;
1412 }
1413 
1414 
1415 #endif