Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConstructGeometry.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ConstructGeometry.C
1 // $Id: $
2 
11 #include <fun4all/Fun4AllDstOutputManager.h>
12 #include <fun4all/Fun4AllDummyInputManager.h>
13 #include <fun4all/Fun4AllInputManager.h>
14 #include <fun4all/Fun4AllOutputManager.h>
15 #include <fun4all/Fun4AllServer.h>
16 #include <fun4all/SubsysReco.h>
17 
18 #include <fun4all/Fun4AllServer.h>
19 
20 #include <g4eval/SvtxEvaluator.h>
21 #include <g4eval/TrkrEvaluator.h>
22 
23 #include <g4detectors/PHG4BlockSubsystem.h>
24 #include <g4detectors/PHG4CylinderSubsystem.h>
25 #include <g4main/PHG4Reco.h>
26 
27 #include <g4eval/PHG4DSTReader.h>
28 
29 #include <g4histos/G4HitNtuple.h>
30 
31 #include <g4main/PHG4ParticleGenerator.h>
32 #include <g4main/PHG4ParticleGun.h>
33 #include <g4main/PHG4Reco.h>
34 #include <g4main/PHG4SimpleEventGenerator.h>
35 #include <g4main/PHG4TruthSubsystem.h>
36 
37 #include <phgeom/PHGeomUtility.h>
38 #include <phool/recoConsts.h>
39 
40 #include <TMath.h>
41 #include <TSystem.h>
42 
43 R__LOAD_LIBRARY(libcalo_reco.so)
44 R__LOAD_LIBRARY(libfun4all.so)
45 R__LOAD_LIBRARY(libg4tpc.so)
46 R__LOAD_LIBRARY(libg4detectors.so)
47 R__LOAD_LIBRARY(libg4eval.so)
48 R__LOAD_LIBRARY(libg4histos.so)
49 R__LOAD_LIBRARY(libg4testbench.so)
50 R__LOAD_LIBRARY(libg4tpc.so)
51 R__LOAD_LIBRARY(libg4intt.so)
52 R__LOAD_LIBRARY(libg4mvtx.so)
53 R__LOAD_LIBRARY(libg4hough.so)
54 R__LOAD_LIBRARY(libg4eval.so)
55 R__LOAD_LIBRARY(libintt.so)
56 R__LOAD_LIBRARY(libmvtx.so)
57 R__LOAD_LIBRARY(libtpc2019.so)
58 R__LOAD_LIBRARY(libtrack_reco.so)
59 
60 void ConstructGeometry()
61 {
62  const double placementR = 0.5 * (40 + 60);
63  const double rotaitonZ = TMath::Pi() + TMath::Pi() * 2 / 12. / 2.;
64  // const double rotaitonZ = 0;
65  const double driftLength = 40;
66  const double cageRadius = 20;
67 
68  gSystem->Load("libfun4all");
69  gSystem->Load("libg4detectors");
70  gSystem->Load("libg4testbench");
71  gSystem->Load("libg4histos");
72  gSystem->Load("libg4eval.so");
73  gSystem->Load("libqa_modules");
74  gSystem->Load("libg4tpc");
75  gSystem->Load("libtrack_io.so");
76  gSystem->Load("libfun4all.so");
77  gSystem->Load("libg4detectors.so");
78  gSystem->Load("libtpc2019.so");
79  gSystem->Load("libg4eval.so");
80  gSystem->Load("libfun4all.so");
81  gSystem->Load("libg4detectors.so");
82  gSystem->Load("libg4hough.so");
83  gSystem->Load("libtrack_reco.so");
84 
86  // Make the Server
88  Fun4AllServer *se = Fun4AllServer::instance();
89  se->Verbosity(1);
90  recoConsts *rc = recoConsts::instance();
91  // only set this if you want a fixed random seed to make
92  // results reproducible for testing
93  // rc->set_IntFlag("RANDOMSEED",12345678);
94 
95  // simulated setup sits at eta=1, theta=40.395 degrees
96  double theta = 90;
97  double phi = 180 + 360 / 12 / 2;
98  // shift in x with respect to midrapidity setup
99  double add_place_z = -driftLength * .5;
100  // Test beam generator
101  PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
102  gen->add_particles("proton", 1); // mu-,e-,anti_proton,pi-
103  gen->set_vertex_distribution_mean(0.0, 0.0, add_place_z);
104  gen->set_vertex_distribution_width(0.0, .0, .0); // Rough beam profile size @ 16 GeV measured by Abhisek
105  gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
106  PHG4SimpleEventGenerator::Gaus,
107  PHG4SimpleEventGenerator::Gaus); // Gauss beam profile
108  double angle = theta * TMath::Pi() / 180.;
109  double eta = -1. * TMath::Log(TMath::Tan(angle / 2.));
110  gen->set_eta_range(eta - 0.001, eta + 0.001); // 1mrad angular divergence
111  gen->set_phi_range(TMath::Pi() * phi / 180 - 0.001, TMath::Pi() * phi / 180 + 0.001); // 1mrad angular divergence
112  const double momentum = 120;
113  gen->set_p_range(momentum, momentum, momentum * 2e-2); // 2% momentum smearing
114  se->registerSubsystem(gen);
115 
116  PHG4Reco *g4Reco = new PHG4Reco();
117  g4Reco->set_field(0);
118  // g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
119 
120  {
121  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_GasVol", 0);
122  cyl->set_double_param("length", driftLength);
123  cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
124  cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
125  cyl->set_double_param("place_z", -driftLength / 2);
126  cyl->set_double_param("radius", 0.0);
127  cyl->set_int_param("lengthviarapidity", 0);
128  cyl->set_string_param("material", "sPHENIX_TPC_Gas");
129  cyl->set_double_param("thickness", cageRadius);
130  cyl->SuperDetector("TPC");
131  cyl->SetActive();
132  cyl->OverlapCheck(1);
133  g4Reco->registerSubsystem(cyl);
134  }
135  {
136  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_FieldCage", 1);
137  cyl->set_double_param("length", driftLength);
138  cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
139  cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
140  cyl->set_double_param("place_z", -driftLength / 2);
141  cyl->set_double_param("radius", cageRadius + 1e-4);
142  cyl->set_int_param("lengthviarapidity", 0);
143  cyl->set_string_param("material", "G4_Cu");
144  cyl->set_double_param("thickness", 0.00347);
145  cyl->SuperDetector("TPC_Support");
146  cyl->SetActive();
147  cyl->OverlapCheck(1);
148  g4Reco->registerSubsystem(cyl);
149  }
150 
151  for (int sign = -1; sign <= 1; sign += 2)
152  {
153  const double endcap_thickness = 0.5;
154 
155  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("TPC_EndCap", 4 + sign);
156  cyl->set_double_param("length", endcap_thickness);
157  cyl->set_double_param("place_x", placementR * TMath::Cos(rotaitonZ));
158  cyl->set_double_param("place_y", placementR * TMath::Sin(rotaitonZ));
159  cyl->set_double_param("place_z", -driftLength / 2 + sign * ((driftLength / 2) + endcap_thickness / 2 + 1e-4));
160  cyl->set_double_param("radius", 0.0);
161  cyl->set_int_param("lengthviarapidity", 0);
162  cyl->set_string_param("material", "G4_Al");
163  cyl->set_double_param("thickness", cageRadius);
164  cyl->SuperDetector("TPC_Support");
165  cyl->SetActive();
166  cyl->OverlapCheck(1);
167  g4Reco->registerSubsystem(cyl);
168  }
169 
170  for (int i = 0; i < 3; ++i)
171  {
172  const double GEMLocation = 0;
173  const double GEMSpacing = -20;
174 
175  PHG4BlockSubsystem *cyl = new PHG4BlockSubsystem("GEM", i);
176  cyl->set_double_param("size_x", 0.5);
177  cyl->set_double_param("size_y", 40);
178  cyl->set_double_param("size_z", 40);
179  cyl->set_double_param("place_x", (GEMLocation + GEMSpacing * i) * TMath::Cos(rotaitonZ));
180  cyl->set_double_param("place_y", (GEMLocation + GEMSpacing * i) * TMath::Sin(rotaitonZ));
181  cyl->set_double_param("place_z", -driftLength / 2);
182  cyl->set_double_param("rot_z", -rotaitonZ * 180 / TMath::Pi());
183  cyl->set_int_param("lengthviarapidity", 0);
184  cyl->set_string_param("material", "ePHEINX_TPC_Gas"); // ePHENIX TPC is Ar + CO2 gas
185  cyl->SuperDetector("GEM");
186  cyl->SetActive();
187  cyl->OverlapCheck(1);
188  g4Reco->registerSubsystem(cyl);
189  }
190 
191  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
192  g4Reco->registerSubsystem(truth);
193 
194  se->registerSubsystem(g4Reco);
195 
196  Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
197  se->registerInputManager(in);
198 
199  char cmd[100];
200  g4Reco->InitRun(se->topNode());
201  g4Reco->ApplyDisplayAction();
202  sprintf(cmd, "/control/execute vis.mac");
203  g4Reco->ApplyCommand(cmd);
204 
205  se->run(1);
206  g4Reco->ApplyCommand("/vis/scene/add/axes 0 0 0 50 cm");
207  g4Reco->ApplyCommand("/vis/viewer/zoom 1");
208 
209  PHGeomUtility::ExportGeomtry(se->topNode(), "TpcPrototypeGeometry.root");
210  PHGeomUtility::ExportGeomtry(se->topNode(), "TpcPrototypeGeometry.gdml");
211 }