ECCE @ EIC 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  // if no material is given, world material is chosen (typically G4_AIR)
73  cyl = new PHG4CylinderSubsystem("MyCylinder1");
74  cyl->set_double_param("radius", 2); // 80 cm
75  cyl->set_double_param("thickness", 0.01); // does not matter (but > 0)
76  cyl->SuperDetector("MyCylinder1");
77  cyl->set_double_param("length", 90.);
78  cyl->SetActive();
79  cyl->SaveAllHits(); // save all hits, also zero energy hits (which are normally discarde)
80  g4Reco->registerSubsystem(cyl);
81 
82  cyl = new PHG4CylinderSubsystem("MyCylinder2");
83  cyl->set_double_param("radius", 70); // 80 cm
84  cyl->set_double_param("thickness", 0.01); // does not matter (but > 0)
85  cyl->SuperDetector("MyCylinder2");
86  cyl->set_double_param("length", 90.);
87  cyl->SetActive();
88  cyl->SaveAllHits();
89  g4Reco->registerSubsystem(cyl);
90 
91  cyl = new PHG4CylinderSubsystem("MyPlane1");
92  cyl->set_double_param("radius", 0); // 80 cm
93  cyl->set_double_param("thickness", 75); // for a cylindrical plane thickness is diameter
94  cyl->set_double_param("length", 0.01); // for a cylindrical plane length is depth
95  cyl->set_double_param("place_z", 100. + 0.005); // position in z, 1/2 depth needs to be added to put front at exactly 100 cm
96  cyl->SuperDetector("MyPlane1");
97  cyl->SetActive();
98  cyl->SaveAllHits();
99  g4Reco->registerSubsystem(cyl);
100 
101  // Black hole swallows everything - prevent loopers from returning
102  // to inner detectors
103  cyl = new PHG4CylinderSubsystem("BlackHole", 0);
104  cyl->set_double_param("radius", 80); // 80 cm
105  cyl->set_double_param("thickness", 0.1); // does not matter (but > 0)
106  cyl->SetActive();
107  cyl->BlackHole(); // eats everything
108  g4Reco->registerSubsystem(cyl);
109 
110  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
111  g4Reco->registerSubsystem(truth);
112 
113  se->registerSubsystem(g4Reco);
114 
115  //---------------------------
116  // fast pattern recognition and full Kalman filter
117  // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
118  //---------------------------
119  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
120  kalman->set_use_vertex_in_fitting(false);
121  kalman->set_sub_top_node_name("SVTX");
122  kalman->set_trackmap_out_name("SvtxTrackMap");
123 
124  // add Si Trtacker
125  kalman->add_phg4hits(
126  "G4HIT_SVTX", // const std::string& phg4hitsNames,
127  PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype,
128  300e-4, // radial-resolution [cm]
129  30e-4, // azimuthal-resolution [cm]
130  1, // z-resolution [cm]
131  1, // efficiency,
132  0 // noise hits
133  );
134 
135  kalman->add_cylinder_state("MyCylinder1", 2.);
136  kalman->add_cylinder_state("MyCylinder2", 100.);
137  kalman->add_zplane_state("MyPlane1", 100.);
138 
139  se->registerSubsystem(kalman);
140 
141  PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
142  fast_sim_eval->set_filename(evalfile);
143  fast_sim_eval->AddProjection("MyCylinder1");
144  fast_sim_eval->AddProjection("MyCylinder2");
145  fast_sim_eval->AddProjection("MyPlane1");
146 
147  se->registerSubsystem(fast_sim_eval);
148  //---------------------------
149 
150  //---------------------------
151  // output DST file for further offlien analysis
152  //---------------------------
153  if (!outfile.empty())
154  {
155  Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outfile);
156  se->registerOutputManager(out);
157  }
158  Fun4AllInputManager *in = new Fun4AllDummyInputManager("ZEN");
159  se->registerInputManager(in);
160 
161  if (nEvents > 0)
162  {
163  se->run(nEvents);
164  // finish job - close and save output files
165  se->End();
166  std::cout << "All done" << std::endl;
167 
168  // cleanup - delete the server and exit
169  delete se;
170  gSystem->Exit(0);
171  }
172  return 0;
173 }
174 
175 PHG4ParticleGenerator *get_gen(const char *name = "PGENERATOR")
176 {
177  Fun4AllServer *se = Fun4AllServer::instance();
178  PHG4ParticleGenerator *pgun = (PHG4ParticleGenerator *) se->getSubsysReco(name);
179  return pgun;
180 }
181 
182 #endif