Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_FGEM_FastSim.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_FGEM_FastSim.C
1 // $Id: G4_FGEM_fsPHENIX.C,v 1.2 2014/01/22 01:44:13 jinhuang Exp $
2 
11 using namespace std;
12 
13 void
14 FGEM_Init()
15 {
16 
17 }
18 
19 void
20 FGEMSetup(PHG4Reco* g4Reco, const int N_Sector = 8, //
21  const double min_eta = 1.45 //
22  )
23 {
24 
25  const double tilt = .1;
26 
27  string name;
28  double etamax;
29  double etamin;
30  double zpos;
31  PHG4SectorSubsystem *gem;
32 
33  make_GEM_station("FGEM_0", g4Reco, 17, 1.01, 2.7, N_Sector);
34  make_GEM_station("FGEM_1", g4Reco, 62, 2.15, 4.0, N_Sector);
35 
37 
38  name = "FGEM_2";
39  etamax = 4;
40  etamin = min_eta;
41  zpos = 1.2e2;
42 
43  gem = new PHG4SectorSubsystem(name.c_str());
44 
45  gem->get_geometry().set_normal_polar_angle(tilt);
46  gem->get_geometry().set_normal_start(
47  zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
48  gem->get_geometry().set_min_polar_angle(
49  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
50  gem->get_geometry().set_max_polar_angle(
51  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
52  gem->get_geometry().set_max_polar_edge(
53  PHG4Sector::Sector_Geometry::FlatEdge());
54  gem->get_geometry().set_material("G4_METHANE");
55  gem->get_geometry().set_N_Sector(N_Sector);
56  gem->OverlapCheck(overlapcheck);
57  AddLayers_MiniTPCDrift(gem);
58  gem->get_geometry().AddLayers_HBD_GEM();
59  g4Reco->registerSubsystem(gem);
60 
62 
63  name = "FGEM_3";
64  etamax = 4;
65  etamin = min_eta;
66  zpos = 1.6e2;
67  gem = new PHG4SectorSubsystem(name.c_str());
68 
69  gem->SuperDetector(name);
70  gem->get_geometry().set_normal_polar_angle(tilt);
71  gem->get_geometry().set_normal_start(
72  zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
73  gem->get_geometry().set_min_polar_angle(
74  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
75  gem->get_geometry().set_max_polar_angle(
76  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
77  gem->get_geometry().set_max_polar_edge(
78  PHG4Sector::Sector_Geometry::FlatEdge());
79  gem->get_geometry().set_material("G4_METHANE");
80  gem->get_geometry().set_N_Sector(N_Sector);
81  gem->OverlapCheck(overlapcheck);
82  AddLayers_MiniTPCDrift(gem);
83  gem->get_geometry().AddLayers_HBD_GEM();
84  g4Reco->registerSubsystem(gem);
85 
86  gem = new PHG4SectorSubsystem(name + "_LowerEta");
87  gem->SuperDetector(name);
88 
89  zpos = zpos
90  - (zpos * sin(tilt)
91  + zpos * cos(tilt)
92  * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt))
93  * sin(tilt);
94 
95  gem->get_geometry().set_normal_polar_angle(
96  (PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta)
97  + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
98  gem->get_geometry().set_normal_start(
99  zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
100  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
101  gem->get_geometry().set_min_polar_angle(
102  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
103  gem->get_geometry().set_max_polar_angle(
104  PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
105  gem->get_geometry().set_material("G4_METHANE");
106  gem->get_geometry().set_N_Sector(N_Sector);
107  gem->get_geometry().set_min_polar_edge(
108  PHG4Sector::Sector_Geometry::FlatEdge());
109 
110  AddLayers_MiniTPCDrift(gem);
111  gem->get_geometry().AddLayers_HBD_GEM();
112  gem->OverlapCheck(overlapcheck);
113  g4Reco->registerSubsystem(gem);
114 
116 
117  name = "FGEM_4";
118  etamax = 4;
119  etamin = min_eta;
120  zpos = 2.75e2;
121  gem = new PHG4SectorSubsystem(name.c_str());
122 
123  gem->SuperDetector(name);
124  gem->get_geometry().set_normal_polar_angle(tilt);
125  gem->get_geometry().set_normal_start(
126  zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
127  gem->get_geometry().set_min_polar_angle(
128  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
129  gem->get_geometry().set_max_polar_angle(
130  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
131  gem->get_geometry().set_max_polar_edge(
132  PHG4Sector::Sector_Geometry::FlatEdge());
133  gem->get_geometry().set_material("G4_METHANE");
134  gem->get_geometry().set_N_Sector(N_Sector);
135  gem->OverlapCheck(overlapcheck);
136  AddLayers_MiniTPCDrift(gem);
137  gem->get_geometry().AddLayers_HBD_GEM();
138  g4Reco->registerSubsystem(gem);
139 
140  zpos = zpos
141  - (zpos * sin(tilt)
142  + zpos * cos(tilt)
143  * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt))
144  * sin(tilt);
145 
146  gem = new PHG4SectorSubsystem(name + "_LowerEta");
147  gem->SuperDetector(name);
148 
149  gem->get_geometry().set_normal_polar_angle(
150  (PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta)
151  + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
152  gem->get_geometry().set_normal_start(
153  zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
154  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
155  gem->get_geometry().set_min_polar_angle(
156  PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
157  gem->get_geometry().set_max_polar_angle(
158  PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
159  gem->get_geometry().set_material("G4_METHANE");
160  gem->get_geometry().set_N_Sector(N_Sector);
161  gem->get_geometry().set_min_polar_edge(
162  PHG4Sector::Sector_Geometry::FlatEdge());
163 
164  AddLayers_MiniTPCDrift(gem);
165  gem->get_geometry().AddLayers_HBD_GEM();
166  gem->OverlapCheck(overlapcheck);
167  g4Reco->registerSubsystem(gem);
168 
170 
171 }
172 
174 void
175 AddLayers_MiniTPCDrift(PHG4SectorSubsystem *gem)
176 {
177  assert(gem);
178 
179  const double cm = PHG4Sector::Sector_Geometry::Unit_cm();
180  const double mm = .1 * cm;
181  const double um = 1e-3 * mm;
182 
183 // const int N_Layers = 70; // used for mini-drift TPC timing digitalization
184  const int N_Layers = 1; // simplified setup
185  const double thickness = 2 * cm;
186 
187  gem->get_geometry().AddLayer("EntranceWindow", "G4_MYLAR", 25 * um, false,
188  100);
189  gem->get_geometry().AddLayer("Cathode", "G4_GRAPHITE", 10 * um, false, 100);
190 
191  for (int d = 1; d <= N_Layers; d++)
192  {
193  stringstream s;
194  s << "DriftLayer_";
195  s << d;
196 
197  gem->get_geometry().AddLayer(s.str(), "G4_METHANE", thickness / N_Layers,
198  true);
199 
200  }
201 }
202 
203 int
204 make_GEM_station(string name, PHG4Reco* g4Reco, double zpos, double etamin,
205  double etamax, const int N_Sector = 8)
206 {
207 
208 // cout
209 // << "make_GEM_station - GEM construction with PHG4SectorSubsystem - make_GEM_station_EdgeReadout of "
210 // << name << endl;
211 
212  double polar_angle = 0;
213 
214  if (zpos < 0)
215  {
216  zpos = -zpos;
217  polar_angle = TMath::Pi();
218 
219  }
220  if (etamax < etamin)
221  {
222  double t = etamax;
223  etamax = etamin;
224  etamin = t;
225  }
226 
227  PHG4SectorSubsystem *gem;
228  gem = new PHG4SectorSubsystem(name.c_str());
229 
230  gem->SuperDetector(name);
231 
232  gem->get_geometry().set_normal_polar_angle(polar_angle);
233  gem->get_geometry().set_normal_start(
234  zpos * PHG4Sector::Sector_Geometry::Unit_cm());
235  gem->get_geometry().set_min_polar_angle(
236  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
237  gem->get_geometry().set_max_polar_angle(
238  PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
239  gem->get_geometry().set_max_polar_edge(
240  PHG4Sector::Sector_Geometry::FlatEdge());
241  gem->get_geometry().set_min_polar_edge(
242  PHG4Sector::Sector_Geometry::FlatEdge());
243  gem->get_geometry().set_N_Sector(N_Sector);
244  gem->get_geometry().set_material("G4_METHANE");
245  gem->OverlapCheck(overlapcheck);
246 
247  AddLayers_MiniTPCDrift(gem);
248  gem->get_geometry().AddLayers_HBD_GEM();
249  g4Reco->registerSubsystem(gem);
250 
251 }
252 
253 
254 void FGEM_FastSim_Reco(int verbosity = 0) {
255 
256  //---------------
257  // Load libraries
258  //---------------
259 
260  gSystem->Load("libfun4all.so");
261  gSystem->Load("libg4hough.so");
262 
263  //---------------
264  // Fun4All server
265  //---------------
266 
267  Fun4AllServer *se = Fun4AllServer::instance();
268 
269  PHG4TrackFastSim* kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
270  kalman->Verbosity(0);
271 
272  kalman->set_use_vertex_in_fitting(true);
273  kalman->set_vertex_xy_resolution(50E-4);
274  kalman->set_vertex_z_resolution(50E-4);
275 
276  kalman->set_detector_type(PHG4TrackFastSim::Vertical_Plane); // Vertical_Plane, Cylinder
277  kalman->set_phi_resolution(50E-4);
278  kalman->set_r_resolution(1.);
279 
280  // now automatic loading magnetic field. See https://github.com/sPHENIX-Collaboration/coresoftware/pull/349
281 // kalman->set_mag_field_file_name("fieldmap.root");
282 // kalman->set_mag_field_re_scaling_factor(1.);
283 
284  kalman->set_pat_rec_hit_finding_eff(1.);
285  kalman->set_pat_rec_noise_prob(0.);
286 
287  std::string phg4hits_names[] = {"G4HIT_FGEM_0","G4HIT_FGEM_1","G4HIT_FGEM_2","G4HIT_FGEM_3","G4HIT_FGEM_4"};
288  kalman->set_phg4hits_names(phg4hits_names, 5);
289  kalman->set_sub_top_node_name("SVTX");
290  kalman->set_trackmap_out_name("SvtxTrackMap");
291 
292  kalman->set_fit_alg_name("KalmanFitterRefTrack");//
293  kalman->set_primary_assumption_pid(13);
294  kalman->set_do_evt_display(false);
295 
296  se->registerSubsystem(kalman);
297 }
298 
299 void FGEM_FastSim_Eval(const char* outname) {
300 
301  //---------------
302  // Load libraries
303  //---------------
304 
305  //gSystem->Load("libfun4all.so");
306  //gSystem->Load("libg4hough.so");
307  gSystem->Load("/phenix/hhj/yuhw//GitHub//HaiwangYu/install//lib/libFastTrackingEval.so");
308 
309 
310  //---------------
311  // Fun4All server
312  //---------------
313 
314  Fun4AllServer *se = Fun4AllServer::instance();
315 
316 
317  FastTrackingEval *fast_sim_eval = new FastTrackingEval("FastTrackingEval");
318  fast_sim_eval->set_filename( Form("%s",outname) );
319  se->registerSubsystem( fast_sim_eval );
320 
321 }