EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_FGEM_fsPHENIX.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_FGEM_fsPHENIX.C
1 #ifndef MACRO_G4FGEMFSPHENIX_C
2 #define MACRO_G4FGEMFSPHENIX_C
3 
4 #include <GlobalVariables.C>
5 
6 #include <G4_FEMC.C>
7 #include <G4_FHCAL.C>
8 
9 #include <g4detectors/PHG4SectorSubsystem.h>
10 
11 #include <g4main/PHG4Reco.h>
12 
13 #include <g4trackfastsim/PHG4TrackFastSim.h>
14 #include <g4trackfastsim/PHG4TrackFastSimEval.h>
15 
16 #include <fun4all/Fun4AllServer.h>
17 
18 R__LOAD_LIBRARY(libg4detectors.so)
19 R__LOAD_LIBRARY(libg4eval.so)
20 R__LOAD_LIBRARY(libg4trackfastsim.so)
21 
22 int make_GEM_station(string name, PHG4Reco *g4Reco, double zpos, double etamin,
23  double etamax, const int N_Sector = 8);
24 void AddLayers_MiniTPCDrift(PHG4SectorSubsystem *gem);
25 
26 namespace Enable
27 {
28  bool FGEM = false;
29  bool FGEM_OVERLAPCHECK = false;
30  bool FGEM_TRACK = false;
31  bool FGEM_EVAL = false;
32  int FGEM_VERBOSITY = 0;
33 } // namespace Enable
34 
35 void FGEM_Init()
36 {
37  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 130.);
38  BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 280.);
39  TRACKING::TrackNodeName = "SvtxTrackMap"; // node name for tracks
40 }
41 
42 void FGEMSetup(PHG4Reco *g4Reco, const int N_Sector = 8, //
43  const double min_eta = 1.45 //
44 )
45 {
46  const double tilt = .1;
47 
48  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::FGEM_OVERLAPCHECK;
49 
50  string name;
51  double etamax;
52  double etamin;
53  double zpos;
54  PHG4SectorSubsystem *gem;
55 
56  make_GEM_station("FGEM_0", g4Reco, 17, 1.01, 2.7, N_Sector);
57  make_GEM_station("FGEM_1", g4Reco, 62, 2.15, 4.0, N_Sector);
58 
60 
61  name = "FGEM_2";
62  etamax = 4;
63  etamin = min_eta;
64  zpos = 1.2e2;
65 
66  gem = new PHG4SectorSubsystem(name);
67 
68  gem->get_geometry().set_normal_polar_angle(tilt);
69  gem->get_geometry().set_normal_start(zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
70  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
71  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
72  gem->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
73  gem->get_geometry().set_material("G4_METHANE");
74  gem->get_geometry().set_N_Sector(N_Sector);
75  gem->OverlapCheck(OverlapCheck);
76  AddLayers_MiniTPCDrift(gem);
77  gem->get_geometry().AddLayers_HBD_GEM();
78  g4Reco->registerSubsystem(gem);
79 
81 
82  name = "FGEM_3";
83  etamax = 4;
84  etamin = min_eta;
85  zpos = 1.6e2;
86  gem = new PHG4SectorSubsystem(name);
87 
88  gem->SuperDetector(name);
89  gem->get_geometry().set_normal_polar_angle(tilt);
90  gem->get_geometry().set_normal_start(zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
91  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
92  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
93  gem->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
94  gem->get_geometry().set_material("G4_METHANE");
95  gem->get_geometry().set_N_Sector(N_Sector);
96  gem->OverlapCheck(OverlapCheck);
97  AddLayers_MiniTPCDrift(gem);
98  gem->get_geometry().AddLayers_HBD_GEM();
99  g4Reco->registerSubsystem(gem);
100 
101  gem = new PHG4SectorSubsystem(name + "_LowerEta");
102  gem->SuperDetector(name);
103 
104  zpos = zpos - (zpos * sin(tilt) + zpos * cos(tilt) * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt)) * sin(tilt);
105 
106  gem->get_geometry().set_normal_polar_angle((PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta) + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
107  gem->get_geometry().set_normal_start(
108  zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
109  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
110  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
111  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
112  gem->get_geometry().set_material("G4_METHANE");
113  gem->get_geometry().set_N_Sector(N_Sector);
114  gem->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
115 
116  AddLayers_MiniTPCDrift(gem);
117  gem->get_geometry().AddLayers_HBD_GEM();
118  gem->OverlapCheck(OverlapCheck);
119  g4Reco->registerSubsystem(gem);
120 
122 
123  name = "FGEM_4";
124  etamax = 4;
125  etamin = min_eta;
126  zpos = 2.75e2;
127  gem = new PHG4SectorSubsystem(name);
128 
129  gem->SuperDetector(name);
130  gem->get_geometry().set_normal_polar_angle(tilt);
131  gem->get_geometry().set_normal_start(zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
132  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
133  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
134  gem->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
135  gem->get_geometry().set_material("G4_METHANE");
136  gem->get_geometry().set_N_Sector(N_Sector);
137  gem->OverlapCheck(OverlapCheck);
138  AddLayers_MiniTPCDrift(gem);
139  gem->get_geometry().AddLayers_HBD_GEM();
140  g4Reco->registerSubsystem(gem);
141 
142  zpos = zpos - (zpos * sin(tilt) + zpos * cos(tilt) * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt)) * sin(tilt);
143 
144  gem = new PHG4SectorSubsystem(name + "_LowerEta");
145  gem->SuperDetector(name);
146 
147  gem->get_geometry().set_normal_polar_angle((PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta) + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
148  gem->get_geometry().set_normal_start(
149  zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
150  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
151  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
152  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
153  gem->get_geometry().set_material("G4_METHANE");
154  gem->get_geometry().set_N_Sector(N_Sector);
155  gem->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
156 
157  AddLayers_MiniTPCDrift(gem);
158  gem->get_geometry().AddLayers_HBD_GEM();
159  gem->OverlapCheck(OverlapCheck);
160  g4Reco->registerSubsystem(gem);
161 
163 }
164 
166 void AddLayers_MiniTPCDrift(PHG4SectorSubsystem *gem)
167 {
168  assert(gem);
169 
170  const double cm = PHG4Sector::Sector_Geometry::Unit_cm();
171  const double mm = .1 * cm;
172  const double um = 1e-3 * mm;
173 
174  // const int N_Layers = 70; // used for mini-drift TPC timing digitalization
175  const int N_Layers = 1; // simplified setup
176  const double thickness = 2 * cm;
177 
178  gem->get_geometry().AddLayer("EntranceWindow", "G4_MYLAR", 25 * um, false, 100);
179  gem->get_geometry().AddLayer("Cathode", "G4_GRAPHITE", 10 * um, false, 100);
180 
181  for (int d = 1; d <= N_Layers; d++)
182  {
183  stringstream s;
184  s << "DriftLayer_";
185  s << d;
186 
187  gem->get_geometry().AddLayer(s.str(), "G4_METHANE", thickness / N_Layers, true);
188  }
189 }
190 
191 int make_GEM_station(string name, PHG4Reco *g4Reco, double zpos, double etamin,
192  double etamax, const int N_Sector = 8)
193 {
194  // cout
195  // << "make_GEM_station - GEM construction with PHG4SectorSubsystem - make_GEM_station_EdgeReadout of "
196  // << name << endl;
197 
198  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::FGEM_OVERLAPCHECK;
199 
200  double polar_angle = 0;
201 
202  if (zpos < 0)
203  {
204  zpos = -zpos;
205  polar_angle = TMath::Pi();
206  }
207  if (etamax < etamin)
208  {
209  double t = etamax;
210  etamax = etamin;
211  etamin = t;
212  }
213 
214  PHG4SectorSubsystem *gem;
215  gem = new PHG4SectorSubsystem(name);
216 
217  gem->SuperDetector(name);
218 
219  gem->get_geometry().set_normal_polar_angle(polar_angle);
220  gem->get_geometry().set_normal_start(zpos * PHG4Sector::Sector_Geometry::Unit_cm());
221  gem->get_geometry().set_min_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
222  gem->get_geometry().set_max_polar_angle(PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
223  gem->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
224  gem->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge());
225  gem->get_geometry().set_N_Sector(N_Sector);
226  gem->get_geometry().set_material("G4_METHANE");
227  gem->OverlapCheck(OverlapCheck);
228 
229  AddLayers_MiniTPCDrift(gem);
230  gem->get_geometry().AddLayers_HBD_GEM();
231  g4Reco->registerSubsystem(gem);
232  return 0;
233 }
234 
235 void FGEM_FastSim_Reco()
236 {
237  int verbosity = std::max(Enable::VERBOSITY, Enable::FGEM_VERBOSITY);
238 
239  //---------------
240  // Fun4All server
241  //---------------
242 
243  Fun4AllServer *se = Fun4AllServer::instance();
244 
245  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
246  kalman->Verbosity(verbosity);
247 
248  kalman->set_use_vertex_in_fitting(true);
249  kalman->set_vertex_xy_resolution(50E-4);
250  kalman->set_vertex_z_resolution(50E-4);
251 
252  kalman->set_sub_top_node_name("SVTX");
253  kalman->set_trackmap_out_name(TRACKING::TrackNodeName);
254 
255  // MAPS in MVTX detector
256  kalman->add_phg4hits(
257  "G4HIT_MVTX", // const std::string& phg4hitsNames,
258  PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype,
259  5e-4, // const float radres,
260  5e-4, // const float phires,
261  5e-4, // const float lonres,
262  1, // const float eff,
263  0 // const float noise
264  );
265 
266  // GEM0, 70um azimuthal resolution, 1cm radial strips
267  kalman->add_phg4hits(
268  "G4HIT_FGEM_0", // const std::string& phg4hitsNames,
269  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
270  1. / sqrt(12.), // const float radres,
271  70e-4, // const float phires,
272  100e-4, // const float lonres,
273  1, // const float eff,
274  0 // const float noise
275  );
276  // GEM1, 70um azimuthal resolution, 1cm radial strips
277  kalman->add_phg4hits(
278  "G4HIT_FGEM_1", // const std::string& phg4hitsNames,
279  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
280  1. / sqrt(12.), // const float radres,
281  70e-4, // const float phires,
282  100e-4, // const float lonres,
283  1, // const float eff,
284  0 // const float noise
285  );
286 
287  // TPC
288  kalman->add_phg4hits(
289  "G4HIT_TPC", // const std::string& phg4hitsNames,
290  PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype,
291  1, // const float radres,
292  200e-4, // const float phires,
293  500e-4, // const float lonres,
294  1, // const float eff,
295  0 // const float noise
296  );
297 
298  // GEM2, 70um azimuthal resolution, 1cm radial strips
299  kalman->add_phg4hits(
300  "G4HIT_FGEM_2", // const std::string& phg4hitsNames,
301  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
302  1. / sqrt(12.), // const float radres,
303  70e-4, // const float phires,
304  100e-4, // const float lonres,
305  1, // const float eff,
306  0 // const float noise
307  );
308  // GEM3, 70um azimuthal resolution, 1cm radial strips
309  kalman->add_phg4hits(
310  "G4HIT_FGEM_3", // const std::string& phg4hitsNames,
311  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
312  1. / sqrt(12.), // const float radres,
313  70e-4, // const float phires,
314  100e-4, // const float lonres,
315  1, // const float eff,
316  0 // const float noise
317  );
318  // GEM4, 70um azimuthal resolution, 1cm radial strips
319  kalman->add_phg4hits(
320  "G4HIT_FGEM_4", // const std::string& phg4hitsNames,
321  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
322  1. / sqrt(12.), // const float radres,
323  70e-4, // const float phires,
324  100e-4, // const float lonres,
325  1, // const float eff,
326  0 // const float noise
327  );
328 
329  // Saved track states (projections)
330 
331  if (Enable::FEMC)
332  {
333  kalman->add_state_name("FEMC");
334  }
335  if (Enable::FHCAL)
336  {
337  kalman->add_state_name("FHCAL");
338  }
339  se->registerSubsystem(kalman);
340 }
341 
342 void FGEM_FastSim_Eval(const std::string &outputfile)
343 {
344  int verbosity = std::max(Enable::VERBOSITY, Enable::FGEM_VERBOSITY);
345  Fun4AllServer *se = Fun4AllServer::instance();
346 
347  PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
348  fast_sim_eval->set_trackmapname(TRACKING::TrackNodeName);
349  fast_sim_eval->set_filename(outputfile);
350  se->registerSubsystem(fast_sim_eval);
351 
352  return;
353 }
354 #endif