EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_FastMom_GEM.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_FastMom_GEM.C
1 #pragma once
2 #include <phgenfit/Track.h>
10 #include <fun4all/SubsysReco.h>
13 #include <phpythia8/PHPythia8.h>
14 #include <g4histos/G4HitNtuple.h>
16 #include <g4main/HepMCNodeReader.h>
19 #include <g4main/PHG4ParticleGun.h>
20 #include <g4main/PHG4Reco.h>
25 #include <phool/recoConsts.h>
26 
27 #include "G4_Jets.C"
28 #include "G4_Bbc.C"
29 #include "G4_Global.C"
30 #include "G4_Pipe_EIC.C"
31 #include "G4_AllSi.C"
32 #include "G4_GEM.C"
33 
35 #include <g4lblvtx/SimpleNtuple.h>
37 #include <myjetanalysis/MyJetAnalysis_AllSi.h>
38 R__LOAD_LIBRARY(libmyjetanalysis.so)
39 R__LOAD_LIBRARY(libfun4all.so)
40 R__LOAD_LIBRARY(libg4detectors.so)
41 R__LOAD_LIBRARY(libg4lblvtx.so)
42 R__LOAD_LIBRARY(libg4trackfastsim.so)
43 R__LOAD_LIBRARY(libPHPythia8.so)
44 
46  int nEvents = -1, // number of events
47  const char *outputFile = "out_allSi", // output filename
48  const char *genpar = "pi-", // particle species to simulate with the simple generators
49  const int det_ver = 2, // version of detector geometry
50  const double pixel_size = 10.) // pixel length (um)
51 {
52  // ======================================================================================================
53  // Input from the user
54  const int particle_gen = 1; // 1 = particle generator, 2 = particle gun, 3 = simple particle generator, 4 = pythia8 e+p collision
55  const int magnetic_field = 4; // 1 = uniform 1.5T, 2 = uniform 3.0T, 3 = sPHENIX 1.4T map, 4 = Beast 3.0T map
56  bool DISPLACED_VERTEX = false; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes
57  bool do_projections = false; // Project momentum vectors to surfaces defined below
58  // ======================================================================================================
59  // Parameters for projections
60  string projname1 = "DIRC"; // Cylindrical surface object name
61  double projradius1 = 50.; // [cm]
62  double length1 = 400.; // [cm]
63  // ---
64  double thinness = 0.1; // black hole thickness, needs to be taken into account for the z positions
65  // ---
66  string projname2 = "FOR"; // Forward plane object name
67  double projzpos2 = 130+thinness/2.; // [cm]
68  double projradius2 = 50.; // [cm]
69  // ---
70  string projname3 = "BACK"; // Backward plane object name
71  double projzpos3 = -(130+thinness/2.);// [cm]
72  double projradius3 = 50.; // [cm]
73  // ======================================================================================================
74  // Make the Server
76  // If you want to fix the random seed for reproducibility
77  // recoConsts *rc = recoConsts::instance();
78  // rc->set_IntFlag("RANDOMSEED", 12345);
79  // ======================================================================================================
80  // Particle Generation
81  if(particle_gen<4) cout << "Particle that will be generated: " << std::string(genpar) << endl;
82  // --------------------------------------------------------------------------------------
83  // Particle Generator Setup
85  gen->set_name(std::string(genpar)); // geantino, pi-, pi+, mu-, mu+, e-., e+, proton, ... (currently passed as an input)
86  gen->set_vtx(0,0,0); // Vertex generation range
87  gen->set_mom_range(.00001,30.); // Momentum generation range in GeV/c
88  gen->set_z_range(0.,0.);
89  gen->set_eta_range(-4,4); // Detector coverage corresponds to |η|< 4
90  gen->set_phi_range(0.,2.*TMath::Pi());
91  // --------------------------------------------------------------------------------------
92  // Particle Gun Setup
93  PHG4ParticleGun *gun = new PHG4ParticleGun();
94  gun->set_name(std::string(genpar)); // geantino, pi-, pi+, mu-, mu+, e-., e+, proton, ...
95  gun->set_vtx(0,0,0);
96  gun->set_mom(0,1,0);
97  // --------------------------------------------------------------------------------------
98  // Simple event generator
100  segen->add_particles(std::string(genpar),100); // 100 pion option
102  segen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
103  segen->set_vertex_distribution_width(0.0, 0.0, 5.0);
105  segen->set_vertex_size_parameters(0.0, 0.0);
106  segen->set_eta_range(-4.,4.);
107  segen->set_phi_range(0.,2.*TMath::Pi());
108  segen->set_p_range(.00001,30.);
109  //segen->Embed(2);
110  segen->Verbosity(0);
111 
112  // --------------------------------------------------------------------------------------
113  // Pythia 8 events
114  bool do_pythia8_jets = false;
115  if (particle_gen==1){se->registerSubsystem( gen); cout << "Using particle generator" << endl;}
116  else if(particle_gen==2){se->registerSubsystem( gun); cout << "Using particle gun" << endl;}
117  else if(particle_gen==3){se->registerSubsystem(segen); cout << "Using simple event generator" << endl;}
118  else if(particle_gen==4){
119  do_pythia8_jets = true;
120  gSystem->Load("libPHPythia8.so");
121 
122  PHPythia8 *pythia8 = new PHPythia8(); // see coresoftware/generators/PHPythia8 for example config
123  pythia8->set_config_file("phpythia8.cfg"); // example configure files : https://github.com/sPHENIX-Collaboration/coresoftware/tree/master/generators/PHPythia8
124  se->registerSubsystem(pythia8);
125 
126  // read-in HepMC events to Geant4 if there are any pythia8 produces hepmc records, so this is needed to read the above generated pythia8 events
127  HepMCNodeReader *hr = new HepMCNodeReader();
128  se->registerSubsystem(hr);
129  }
130 
131  // ======================================================================================================
132  PHG4Reco *g4Reco = new PHG4Reco();
133  //g4Reco->SetWorldMaterial("G4_Galactic");
134  // ======================================================================================================
135  // Magnetic field setting
136  TString B_label;
137  if(magnetic_field==1){ // uniform 1.5T
138  B_label = "_B_1.5T";
139  g4Reco->set_field(1.5);
140  }
141  else if(magnetic_field==2){ // uniform 3.0T
142  B_label = "_B_3.0T";
143  g4Reco->set_field(3.0);
144  }
145  else if(magnetic_field==3){ // sPHENIX 1.4T map
146  B_label = "_sPHENIX";
147  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"), PHFieldConfig::kField2D);
148  g4Reco->set_field_rescale(-1.4/1.5);
149  }
150  else if(magnetic_field==4){ // Beast 3.0T map
151  B_label = "_Beast";
152  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/mfield.4col.dat"), PHFieldConfig::kFieldBeast);
153  }
154  else{ // The user did not provide a valid B field setting
155  cout << "User did not provide a valid magnetic field setting. Set 'magnetic_field'. Bailing out!" << endl;
156  }
157  // ======================================================================================================
158  // Physics list (default list is "QGSP_BERT")
159  //g4Reco->SetPhysicsList("FTFP_BERT_HP"); // This list is slower and only useful for hadronic showers.
160  // ======================================================================================================
161  load_AllSi_geom(g4Reco, det_ver); // Loading All-Si Tracker and beampipe geometries
162  EGEM_Init(); // Loading backward GEM geometry
163  FGEM_Init(); // Loading forward GEM geometry
164  EGEMSetup(g4Reco);
165  FGEMSetup(g4Reco);
166  //addPassiveMaterial(g4Reco);
167 
168  // ======================================================================================================
169  if(do_projections){
171  cyl = new PHG4CylinderSubsystem(projname1,0);
172  cyl->set_double_param("length", length1);
173  cyl->set_double_param("radius", projradius1); // dirc radius
174  cyl->set_double_param("thickness", 0.1); // needs some thickness
175  cyl->set_string_param("material", "G4_AIR");
176  cyl->SetActive(1);
177  cyl->SuperDetector(projname1);
178  cyl->BlackHole();
179  cyl->set_color(1,0,0,0.7); //reddish
180  g4Reco->registerSubsystem(cyl);
181 
182  cyl = new PHG4CylinderSubsystem(projname2,0);
183  cyl->set_double_param("length", thinness);
184  cyl->set_double_param("radius", 2); // beampipe needs to fit here
185  cyl->set_double_param("thickness", projradius2); //
186  cyl->set_string_param("material", "G4_AIR");
187  cyl->set_double_param("place_z", projzpos2);
188  cyl->SetActive(1);
189  cyl->SuperDetector(projname2);
190  cyl->BlackHole();
191  cyl->set_color(0,1,1,0.3); //reddish
192  g4Reco->registerSubsystem(cyl);
193 
194  cyl = new PHG4CylinderSubsystem(projname3,0);
195  cyl->set_double_param("length", thinness);
196  cyl->set_double_param("radius", 2); // beampipe needs to fit here
197  cyl->set_double_param("thickness", projradius3); //
198  cyl->set_string_param("material", "G4_AIR");
199  cyl->set_double_param("place_z", projzpos3);
200  cyl->SetActive(1);
201  cyl->SuperDetector(projname3);
202  cyl->BlackHole();
203  cyl->set_color(0,1,1,0.3); //reddish
204  g4Reco->registerSubsystem(cyl);
205  }
206  // ======================================================================================================
207 
208  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
209  g4Reco->registerSubsystem(truth);
210 
211  se->registerSubsystem(g4Reco);
212 
213  // ======================================================================================================
214  if(do_pythia8_jets){
215  Bbc_Reco();
216  Global_Reco();
217  }
218 
219  // ======================================================================================================
220  // Fast pattern recognition and full Kalman filter
221  // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
222  double um_to_cm = 1E-04; // Conversion factor from um to cm
223  char nodename[100];
224  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
225  kalman->set_use_vertex_in_fitting(false);
226  kalman->set_sub_top_node_name("SVTX");
227  kalman->set_trackmap_out_name("SvtxTrackMap");
228 
229  add_AllSi_to_kalman( kalman , pixel_size , det_ver ); // Add All-Silicon tracker to Kalman filter
230 
231  //-------------------------
232  // Adding GEMs to the Kalman filter
233  // BACKWARD GEM, 70um azimuthal resolution, 1cm radial strips
234  kalman->add_phg4hits("G4HIT_EGEM", // const std::string& phg4hitsNames,
235  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
236  50e-4,//1. / sqrt(12.), // const float radres,
237  50e-4, // const float phires,
238  999., // longitudinal (z) resolution [cm] (this number is not used in vertical plane geometry)
239  1, // const float eff,
240  0 // const float noise
241  );
242  // FORWARD GEM2, 70um azimuthal resolution, 1cm radial strips
243  kalman->add_phg4hits("G4HIT_FGEM", // const std::string& phg4hitsNames,
244  PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype,
245  50e-4,//1. / sqrt(12.), // const float radres,
246  50e-4, // const float phires,
247  999., // longitudinal (z) resolution [cm] (this number is not used in vertical plane geometry)
248  1, // const float eff,
249  0 // const float noise
250  );
251 
252  //-------------------------
253  // Projections
254  if(do_projections){
255  kalman->add_cylinder_state(projname1, projradius1); // projection on cylinder (DIRC)
256  kalman->add_zplane_state (projname2, projzpos2 ); // projection on vertical planes
257  kalman->add_zplane_state (projname3, projzpos3 ); // projection on vertical planes
258  }
259 
260  if(DISPLACED_VERTEX){
261  // use very loose vertex constraint (1cm in sigma) to allow reco of displaced vertex
262  kalman->set_use_vertex_in_fitting(true);
263  kalman->set_vertex_xy_resolution(1);
264  kalman->set_vertex_z_resolution(1);
265  kalman->enable_vertexing(true);
266  kalman->set_vertex_min_ndf(2); // Only tracks with number of degrees of freedom greater than this value are used to fit the vertex
267  }
268  /*
269  else{
270  kalman->set_use_vertex_in_fitting(false);
271  kalman->set_vertex_xy_resolution(0);
272  kalman->set_vertex_z_resolution(0);
273  kalman->enable_vertexing(false); // this is false by default
274  kalman->set_vertex_min_ndf(2);
275  }
276  */
277  //kalman->Verbosity(10);
278  se->registerSubsystem(kalman);
279  // -----------------------------------------------------
280  // INFO: The resolution numbers above correspond to:
281  // 20e-4/sqrt(12) cm = 5.8e-4 cm, to simulate 20x20 um
282 
283  // ======================================================================================================
284  TrackFastSimEval *fast_sim_eval = new TrackFastSimEval("FastTrackingEval");
285  fast_sim_eval->set_filename(TString(outputFile)+B_label+"_FastTrackingEval.root");
286  if(do_projections){
287  fast_sim_eval->AddProjection(projname1);
288  fast_sim_eval->AddProjection(projname2);
289  fast_sim_eval->AddProjection(projname3);
290  }
291  se->registerSubsystem(fast_sim_eval);
292 
293  // ======================================================================================================
294  // resonstruct jets after the tracking
295  if(do_pythia8_jets) Jet_Reco();
296 
297  SimpleNtuple *hits = new SimpleNtuple("Hits");
298  add_AllSi_hits(hits,det_ver); // Add All-Silicon tracker hits
299  hits->AddNode("EGEM");
300  hits->AddNode("FGEM");
301  se->registerSubsystem(hits);
302 
303  // ======================================================================================================
304  if(do_pythia8_jets)
305  Jet_Eval(string(outputFile) + "_g4jet_eval.root");
306  // ======================================================================================================
307 
309  // IOManagers...
311  const std::string dst_name = std::string(outputFile)+"_G4LBLVtx.root";
312  //Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",TString(outputFile)+"_G4LBLVtx.root");
313  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",dst_name);
314  out->Verbosity(0);
315  se->registerOutputManager(out);
316 
318  se->registerInputManager(in);
319 
320  if(do_pythia8_jets){
321  gSystem->Load("libmyjetanalysis");
322  std::string jetoutputFile = std::string(outputFile) + std::string("_electrons+jets.root");
323  MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r10","AntiKt_Truth_r10",jetoutputFile.data());
324  //MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r04","AntiKt_Truth_r04",jetoutputFile.data());
325  //MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r02","AntiKt_Truth_r02",jetoutputFile.data());
326  se->registerSubsystem(myJetAnalysis);
327  }
328 
329  if (nEvents <= 0) return;
330 
331  se->run(nEvents);
332  se->End();
333  delete se;
334 
335  gSystem->Exit(0);
336 }