EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_simple_vertex.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_simple_vertex.C
1 #pragma once
2 #include <phgenfit/Track.h>
3 #include <fun4all/Fun4AllDstInputManager.h>
4 #include <fun4all/Fun4AllDstOutputManager.h>
5 #include <fun4all/Fun4AllDummyInputManager.h>
6 #include <fun4all/Fun4AllInputManager.h>
7 #include <fun4all/Fun4AllNoSyncDstInputManager.h>
8 #include <fun4all/Fun4AllOutputManager.h>
9 #include <fun4all/Fun4AllServer.h>
10 #include <fun4all/SubsysReco.h>
11 #include <g4detectors/PHG4DetectorSubsystem.h>
12 #include <g4detectors/PHG4CylinderSubsystem.h>
13 #include <g4histos/G4HitNtuple.h>
14 #include <g4main/PHG4ParticleGenerator.h>
15 #include <g4main/PHG4ParticleGeneratorBase.h>
16 #include <g4main/PHG4Reco.h>
17 #include <g4main/PHG4TruthSubsystem.h>
18 #include <g4trackfastsim/PHG4TrackFastSim.h>
19 #include <g4trackfastsim/PHG4TrackFastSimEval.h>
20 #include <phool/recoConsts.h>
21 
22 R__LOAD_LIBRARY(libfun4all.so)
23 R__LOAD_LIBRARY(libg4detectors.so)
24 R__LOAD_LIBRARY(libg4lblvtx.so)
25 R__LOAD_LIBRARY(libg4trackfastsim.so)
26 
27  void Fun4All_G4_simple_vertex(
28  int nEvents = -1, // number of events
29  bool vtx_lyr_1 = true,
30  bool vtx_lyr_2 = true,
31  bool vtx_lyr_3 = true,
32  double vtx_matBud = 0.05, //% X/X0
33  double pix_size = 10.,
34  TString out_name = "out_vtx_study") // output filename
35 {
36  std::string outputFile = std::string(out_name)+"_FastSimEval.root";
37  // ======================================================================================================
38  // Make the Server
39  Fun4AllServer *se = Fun4AllServer::instance();
40  // ======================================================================================================
41  // Particle Generator Setup
42  PHG4ParticleGenerator *gen = new PHG4ParticleGenerator();
43  gen->set_name(std::string("pi-")); // geantino, pi-, pi+, mu-, mu+, e-., e+, proton, ... (currently passed as an input)
44  gen->set_vtx(0,0,0); // Vertex generation range
45  gen->set_mom_range(0.,8.); // Momentum generation range in GeV/c
46  gen->set_z_range(0.,0.);
47  gen->set_eta_range(0,1);
48  gen->set_phi_range(0,2*TMath::Pi());
49  se->registerSubsystem(gen);
50  // ======================================================================================================
51  PHG4Reco *g4Reco = new PHG4Reco();
52  g4Reco->SetWorldMaterial("G4_Galactic");
53  g4Reco->set_field(3.0);
54  // ======================================================================================================
55  // Detector setup
56  double si_r_pos[] = {3.64,4.45,5.26,21.,22.68,39.3,43.23};
57  const int nTrckLayers = sizeof(si_r_pos)/sizeof(*si_r_pos);
58  double si_z_length[] = {14.,14.,14.,18.,20.,35.,38.};
59  for(int i = 0 ; i < nTrckLayers ; i++) si_z_length[i] *= 3.;
60  double si_thick_vtx = vtx_matBud/100.*9.37;
61  double si_thick_bar = 0.55/100.*9.37;
62 
63  PHG4CylinderSubsystem *cyl;
64  for (int ilayer = 0; ilayer < nTrckLayers ; ilayer++)
65  {
66  if(
67  (ilayer==0&&vtx_lyr_1)||
68  (ilayer==1&&vtx_lyr_2)||
69  (ilayer==2&&vtx_lyr_3)||
70  (ilayer>2)
71  ){
72  cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
73  cyl->set_string_param("material" , "G4_Si" );
74  cyl->set_double_param("radius" , si_r_pos[ilayer]);
75 
76  if(ilayer<2)
77  cyl->set_double_param("thickness", si_thick_vtx);
78  else
79  cyl->set_double_param("thickness", si_thick_bar);
80 
81  cyl->set_double_param("place_z" , 0 );
82  cyl->set_double_param("length" , si_z_length[ilayer] );
83  cyl->SetActive();
84  cyl->SuperDetector("SVTX");
85  g4Reco->registerSubsystem(cyl);
86  }
87  }
88 
89  //---------------------------
90  // mid-rapidity beryllium pipe
91  double be_pipe_radius = 3.1000;
92  double be_pipe_thickness = 3.1762 - be_pipe_radius; // 760 um for sPHENIX
93  double be_pipe_length_plus = 66.8; // +z beam pipe extend.
94  double be_pipe_length_neg = -79.8; // -z beam pipe extend.
95  double be_pipe_length = be_pipe_length_plus - be_pipe_length_neg;
96  double be_pipe_center = 0.5 * (be_pipe_length_plus + be_pipe_length_neg);
97 
98  cyl = new PHG4CylinderSubsystem("BE_PIPE", 1);
99  cyl->set_double_param("radius", be_pipe_radius);
100  cyl->set_int_param("lengthviarapidity", 0);
101  cyl->set_double_param("length", be_pipe_length);
102  cyl->set_double_param("place_z", be_pipe_center);
103  cyl->set_string_param("material", "G4_Be");
104  cyl->set_double_param("thickness", be_pipe_thickness);
105  cyl->SuperDetector("PIPE");
106  g4Reco->registerSubsystem(cyl);
107  //---------------------------
108 
109  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
110  g4Reco->registerSubsystem(truth);
111 
112  se->registerSubsystem(g4Reco);
113 
114  //---------------------------
115  // fast pattern recognition and full Kalman filter
116  // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
117  //---------------------------
118  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
119  kalman->set_use_vertex_in_fitting(false);
120  kalman->set_sub_top_node_name("SVTX");
121  kalman->set_trackmap_out_name("SvtxTrackMap");
122 
123  // add Si Trtacker
124  kalman->add_phg4hits(
125  "G4HIT_SVTX", // const std::string& phg4hitsNames,
126  PHG4TrackFastSim::Cylinder,
127  999., // radial-resolution [cm]
128  pix_size/10000./sqrt(12.), // azimuthal-resolution [cm]
129  pix_size/10000./sqrt(12.), // z-resolution [cm]
130  1, // efficiency,
131  0 // noise hits
132  );
133  //kalman->Verbosity(10);
134  kalman->set_use_vertex_in_fitting(false);
135  kalman->set_vertex_xy_resolution(0);
136  kalman->set_vertex_z_resolution(0);
137  kalman->enable_vertexing(false); // this is false by default
138  kalman->set_vertex_min_ndf(2);
139 
140  se->registerSubsystem(kalman);
141 
142  PHG4TrackFastSimEval *fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval");
143  fast_sim_eval->set_filename(outputFile);
144  se->registerSubsystem(fast_sim_eval);
145 
146  // ======================================================================================================
147  // IOManagers...
148  const std::string dst_name = std::string(out_name)+"_G4LBLVtx.root";
149  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",dst_name);
150  out->Verbosity(0);
151  se->registerOutputManager(out);
152 
153  Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
154  se->registerInputManager(in);
155 
156  if (nEvents <= 0) return;
157 
158  se->run(nEvents);
159  se->End();
160  delete se;
161 
162  gSystem->Exit(0);
163 }