Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_Momentum_Projection_Detectors.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_Momentum_Projection_Detectors.C
1 #ifndef FUN4ALL_G4_MOMENTUM_PROJECTION_DETECTORS_C
2 #define FUN4ALL_G4_MOMENTUM_PROJECTION_DETECTORS_C
3 
4 #include <g4detectors/PHG4CylinderSubsystem.h>
5 
6 #include <g4trackfastsim/PHG4TrackFastSim.h>
7 #include <g4trackfastsim/PHG4TrackFastSimEval.h>
8 
9 #include <g4main/PHG4ParticleGenerator.h>
10 #include <g4main/PHG4ParticleGun.h>
11 #include <g4main/PHG4Reco.h>
12 #include <g4main/PHG4TruthSubsystem.h>
13 
14 #include <fun4all/Fun4AllDstOutputManager.h>
15 #include <fun4all/Fun4AllDummyInputManager.h>
16 #include <fun4all/Fun4AllInputManager.h>
17 #include <fun4all/Fun4AllOutputManager.h>
18 #include <fun4all/Fun4AllServer.h>
19 #include <fun4all/SubsysReco.h>
20 
21 #include <phool/recoConsts.h>
22 
23 R__LOAD_LIBRARY(libfun4all.so)
24 R__LOAD_LIBRARY(libg4testbench.so)
25 R__LOAD_LIBRARY(libg4detectors.so)
26 R__LOAD_LIBRARY(libg4trackfastsim.so)
27 
28 int Fun4All_G4_Momentum_Projection_Detectors(const int nEvents = 1000, const string &evalfile = "FastTrackingEval.root", const string &outfile = "")
29 {
31  // Make the Server
33  Fun4AllServer *se = Fun4AllServer::instance();
34  se->Verbosity(0);
35 
36  recoConsts *rc = recoConsts::instance();
37  // if you want to use a fixed seed for reproducible results
38  // rc->set_IntFlag("RANDOMSEED", 12345); // if you want to use a fixed seed
39  // PHG4ParticleGenerator generates particle
40  // distributions in eta/phi/mom range
41  PHG4ParticleGenerator *gen = new PHG4ParticleGenerator("PGENERATOR");
42  gen->set_name("pi-");
43  gen->set_vtx(0, 0, 0);
44  gen->set_eta_range(0.5, 1.5);
45  gen->set_mom_range(2, 2); // GeV/c
46  gen->set_phi_range(0., 90. / 180. * TMath::Pi()); // 0-90 deg
47  se->registerSubsystem(gen);
48 
49  PHG4Reco *g4Reco = new PHG4Reco();
50  g4Reco->set_field(1.5); // 1.5 T solenoidal field
51 
52  double si_thickness[6] = {0.02, 0.02, 0.0625, 0.032, 0.032, 0.032};
53  double svxrad[6] = {2.71, 4.63, 11.765, 25.46, 41.38, 63.66};
54  double length[6] = {20., 20., 36., -1., -1., -1.}; // -1 use eta coverage to determine length
55  PHG4CylinderSubsystem *cyl;
56  // here is our silicon:
57  for (int ilayer = 0; ilayer < 6; ilayer++)
58  {
59  cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
60  cyl->set_double_param("radius", svxrad[ilayer]);
61  cyl->set_string_param("material", "G4_Si");
62  cyl->set_double_param("thickness", si_thickness[ilayer]);
63  cyl->SetActive();
64  cyl->SuperDetector("SVTX");
65  if (length[ilayer] > 0)
66  {
67  cyl->set_double_param("length", length[ilayer]);
68  }
69  g4Reco->registerSubsystem(cyl);
70  }
71 
72  // projection volumes (thin cylinders)
73  // if no material is given, world material is chosen (typically G4_AIR)
74  double mycylinder1_radius = 2.; // in cm
75  cyl = new PHG4CylinderSubsystem("MyCylinder1");
76  cyl->set_double_param("radius", mycylinder1_radius);
77  cyl->set_double_param("thickness", 0.01); // does not matter (but > 0)
78  cyl->SuperDetector("MyCylinder1");
79  cyl->set_double_param("length", 90.);
80  cyl->SetActive();
81  cyl->SaveAllHits(); // save all hits, also zero energy hits (which are normally discarded)
82  g4Reco->registerSubsystem(cyl);
83 
84  double mycylinder2_radius = 70.; // in cm
85  cyl = new PHG4CylinderSubsystem("MyCylinder2");
86  cyl->set_double_param("radius", mycylinder2_radius);
87  cyl->set_double_param("thickness", 0.01); // does not matter (but > 0)
88  cyl->SuperDetector("MyCylinder2");
89  cyl->set_double_param("length", 90.);
90  cyl->SetActive();
91  cyl->SaveAllHits();
92  g4Reco->registerSubsystem(cyl);
93 
94  double myplane1_z = 100.;
95  cyl = new PHG4CylinderSubsystem("MyPlane1");
96  cyl->set_double_param("radius", 0); // 80 cm
97  cyl->set_double_param("thickness", 75); // for a cylindrical plane thickness is diameter
98  cyl->set_double_param("length", 0.01); // for a cylindrical plane length is depth
99  cyl->set_double_param("place_z", myplane1_z + 0.005); // position in z, 1/2 depth needs to be added to put front at exactly 100 cm
100  cyl->SuperDetector("MyPlane1");
101  cyl->SetActive();
102  cyl->SaveAllHits();
103  g4Reco->registerSubsystem(cyl);
104 
105  // Black hole swallows everything - prevent loopers from returning
106  // to inner detectors
107  cyl = new PHG4CylinderSubsystem("BlackHole", 0);
108  cyl->set_double_param("radius", 80); // 80 cm
109  cyl->set_double_param("thickness", 0.1); // does not matter (but > 0)
110  cyl->SetActive();
111  cyl->BlackHole(); // eats everything
112  g4Reco->registerSubsystem(cyl);
113 
114  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
115  g4Reco->registerSubsystem(truth);
116 
117  se->registerSubsystem(g4Reco);
118 
119  //---------------------------
120  // fast pattern recognition and full Kalman filter
121  // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
122  //---------------------------
123  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
124  kalman->set_use_vertex_in_fitting(false);
125  kalman->set_sub_top_node_name("SVTX");
126  kalman->set_trackmap_out_name("SvtxTrackMap");
127 
128  // add Si Trtacker
129  kalman->add_phg4hits(
130  "G4HIT_SVTX", // const std::string& phg4hitsNames,
131  PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype,
132  300e-4, // radial-resolution [cm]
133  30e-4, // azimuthal-resolution [cm]
134  1, // z-resolution [cm]
135  1, // efficiency,
136  0 // noise hits
137  );
138 
139  kalman->add_cylinder_state("MyCylinder1", mycylinder1_radius);
140  kalman->add_cylinder_state("MyCylinder2", mycylinder2_radius);
141  kalman->add_zplane_state("MyPlane1", myplane1_z);
142 
143  se->registerSubsystem(kalman);
144 
145  PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
146  fast_sim_eval->set_filename(evalfile);
147  fast_sim_eval->AddProjection("MyCylinder1");
148  fast_sim_eval->AddProjection("MyCylinder2");
149  fast_sim_eval->AddProjection("MyPlane1");
150 
151  se->registerSubsystem(fast_sim_eval);
152  //---------------------------
153 
154  //---------------------------
155  // output DST file for further offlien analysis
156  //---------------------------
157  if (!outfile.empty())
158  {
159  Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outfile);
160  se->registerOutputManager(out);
161  }
162  Fun4AllInputManager *in = new Fun4AllDummyInputManager("ZEN");
163  se->registerInputManager(in);
164 
165  if (nEvents > 0)
166  {
167  se->run(nEvents);
168  // finish job - close and save output files
169  se->End();
170  std::cout << "All done" << std::endl;
171 
172  // cleanup - delete the server and exit
173  delete se;
174  gSystem->Exit(0);
175  }
176  return 0;
177 }
178 
179 PHG4ParticleGenerator *get_gen(const char *name = "PGENERATOR")
180 {
181  Fun4AllServer *se = Fun4AllServer::instance();
182  PHG4ParticleGenerator *pgun = (PHG4ParticleGenerator *) se->getSubsysReco(name);
183  return pgun;
184 }
185 
186 #endif