EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_Magnet.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_Magnet.C
1 #pragma once
2 #include <fun4all/SubsysReco.h>
3 #include <fun4all/Fun4AllServer.h>
4 #include <fun4all/Fun4AllInputManager.h>
5 #include <fun4all/Fun4AllDummyInputManager.h>
6 #include <g4detectors/PHG4CylinderSubsystem.h>
7 #include <g4histos/G4HitNtuple.h>
8 #include <g4main/PHG4ParticleGun.h>
9 #include <g4main/PHG4Reco.h>
10 #include <phool/recoConsts.h>
11 #include <eicdetectors/BeastMagnetSubsystem.h>
12 
13 R__LOAD_LIBRARY(libfun4all.so)
14 R__LOAD_LIBRARY(libg4detectors.so)
15 R__LOAD_LIBRARY(libg4histos.so)
16 R__LOAD_LIBRARY(libeicdetectors.so)
17 
18 void Fun4All_G4_Magnet(int nEvents = -1)
19 {
20 
22  // Make the Server
24  Fun4AllServer *se = Fun4AllServer::instance();
25  recoConsts *rc = recoConsts::instance();
26 
27 // set magnet you want to run to true
28  bool use_babar = true;
29  bool use_beast = false;
30  bool use_cleo = false;
31 
32 // replace fieldmap with constant solenoidal field from x/y/z = 0
33  bool use_solenoid_field = false;
34 
35 // make magnet active volume if you want to study the hits
36  bool magnet_active=false;
37 
38 // if you want to fix the random seed to reproduce results
39 // set this flag
40 // nail this down so I know what the first event looks like...
41 // rc->set_IntFlag("RANDOMSEED",12345);
42 
43  double pz = 1.; // some forward momentum to prevent geantino loopers
44 
45 // ParticleGun shoots right into the original MyDetector volume
46  PHG4ParticleGun *gun = new PHG4ParticleGun();
47 // gun->set_name("pi-");
48  gun->set_name("chargedgeantino");
49  gun->set_vtx(0, 0, 0);
50  gun->set_mom(0, 0.5, pz);
51  gun->AddParticle("chargedgeantino",0,1,pz);
52  gun->AddParticle("chargedgeantino",0,2,pz);
53  gun->AddParticle("chargedgeantino",0,3,pz);
54  se->registerSubsystem(gun);
55 
56 //
57 // Geant4 setup
58 //
59  PHG4Reco* g4Reco = new PHG4Reco();
60 // setup of G4:
61 // no saving of geometry: it takes time and we do not do tracking
62 // so we do not need the geometry
63  g4Reco->save_DST_geometry(false);
64 
65  if (use_beast)
66  {
67  if (use_solenoid_field)
68  {
69  g4Reco->set_field(3);
70  }
71  else
72  {
73  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/mfield.4col.dat"), PHFieldConfig::kFieldBeast);
74 // g4Reco->set_field_rescale(0.5);
75  }
76  BeastMagnetSubsystem *beast = new BeastMagnetSubsystem();
77  beast->set_string_param("GDMPath",(string(getenv("CALIBRATIONROOT")) + string("/Magnet/BeastSolenoid.gdml")));
78  beast->set_string_param("TopVolName","SOLENOID");
79  beast->SetActive(magnet_active);
80  beast->SuperDetector("MAGNET");
81  g4Reco->registerSubsystem(beast);
82  }
83  else if (use_cleo)
84  {
85  if (use_solenoid_field)
86  {
87  g4Reco->set_field(2.02);
88  }
89  else
90  {
91  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/SolenoidMag3D.TABLE"), PHFieldConfig::kFieldCleo);
92 // g4Reco->set_field_rescale(0.5);
93  }
94  double magnet_inner_radius = 137.;
95  double magnet_thickness = 20.;
96  double magnet_length = 350.;
97  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("MAGNET",0);
98  cyl->set_double_param("radius",magnet_inner_radius);
99  cyl->set_double_param("length",magnet_length);
100  cyl->set_double_param("thickness",magnet_thickness);
101  cyl->set_string_param("material","G4_Al");
102  cyl->set_color(0,1,0,1);
103  cyl->SuperDetector("MAGNET");
104  cyl->SetActive(magnet_active);
105  g4Reco->registerSubsystem( cyl );
106  }
107  else if (use_babar)
108  {
109  if (use_solenoid_field)
110  {
111  g4Reco->set_field(1.4);
112  }
113  else
114  {
115  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"), PHFieldConfig::kField2D);
116  g4Reco->set_field_rescale(-1.4/1.5);
117  }
118  double magnet_inner_cryostat_wall_radius = 142;
119  double magnet_inner_cryostat_wall_thickness = 1;
120  double magnet_outer_cryostat_wall_radius = 174.5;
121  double magnet_outer_cryostat_wall_thickness = 2.5;
122  double magnet_coil_radius = 150.8;
123  double magnet_coil_thickness = 9.38;
124  double magnet_length = 379.;
125  double coil_length = 361.5;
126  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("MAGNET",0);
127  cyl->set_double_param("radius",magnet_inner_cryostat_wall_radius);
128  cyl->set_int_param("lengthviarapidity",0);
129  cyl->set_double_param("length",magnet_length);
130  cyl->set_double_param("thickness",magnet_inner_cryostat_wall_thickness);
131  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
132  cyl->SuperDetector("MAGNET");
133  cyl->SetActive(magnet_active);
134  cyl->set_color(0,0,1,1);
135  g4Reco->registerSubsystem( cyl );
136 
137  cyl = new PHG4CylinderSubsystem("MAGNET", 1);
138  cyl->set_double_param("radius",magnet_coil_radius);
139  cyl->set_int_param("lengthviarapidity",0);
140  cyl->set_double_param("length",coil_length);
141  cyl->set_double_param("thickness",magnet_coil_thickness);
142  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
143  cyl->SuperDetector("MAGNET");
144  cyl->SetActive(magnet_active);
145  cyl->set_color(0,0,1,1);
146  g4Reco->registerSubsystem( cyl );
147 
148  cyl = new PHG4CylinderSubsystem("MAGNET", 2);
149  cyl->set_double_param("radius",magnet_outer_cryostat_wall_radius);
150  cyl->set_int_param("lengthviarapidity",0);
151  cyl->set_double_param("length",magnet_length);
152  cyl->set_double_param("thickness",magnet_outer_cryostat_wall_thickness);
153  cyl->set_string_param("material","Al5083"); // use 1 radiation length Al for magnet thickness
154  cyl->SuperDetector("MAGNET");
155  cyl->SetActive(magnet_active);
156  cyl->set_color(0,0,1,1);
157  g4Reco->registerSubsystem(cyl);
158  }
159 
160  se->registerSubsystem( g4Reco );
161 
163  // Fun4All modules
165 
166  G4HitNtuple *hits = new G4HitNtuple("Hits");
167  hits->AddNode("MAGNET",0);
168  se->registerSubsystem(hits);
169 
170 // this (dummy) input manager just drives the event loop
171  Fun4AllInputManager *in = new Fun4AllDummyInputManager( "Dummy");
172  se->registerInputManager( in );
173 // events = 0 => run forever, default is -1, do not run event loop
174  if (nEvents <= 0)
175  {
176  return 0;
177  }
178  se->run(nEvents);
179  se->End();
180  delete se;
181  gSystem->Exit(0);
182 }