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