EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_FastMom.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_FastMom.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 
34 #include <g4lblvtx/SimpleNtuple.h>
36 #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 = 5; // 1 = particle generator, 2 = particle gun, 3 = simple event generator, 4 = pythia8 e+p collision, 5 = particle generator flat in pT
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 generator flat in pT
94  gen_pT->set_name(std::string(genpar)); // geantino, pi-, pi+, mu-, mu+, e-., e+, proton, ... (currently passed as an input)
95  gen_pT->set_vtx(0,0,0); // Vertex generation range
96  gen_pT->set_pT_range(.00001,30.); // Momentum generation range in GeV/c
97  gen_pT->set_z_range(0.,0.);
98  gen_pT->set_eta_range(-4,4); // Detector coverage corresponds to |η|< 4
99  gen_pT->set_phi_range(0.,2.*TMath::Pi());
100  // --------------------------------------------------------------------------------------
101  // Particle Gun Setup
102  PHG4ParticleGun *gun = new PHG4ParticleGun();
103  gun->set_name(std::string(genpar)); // geantino, pi-, pi+, mu-, mu+, e-., e+, proton, ...
104  gun->set_vtx(0,0,0);
105  gun->set_mom(0,1,0);
106  // --------------------------------------------------------------------------------------
107  // Simple event generator
109  segen->add_particles(std::string(genpar),100); // 100 pion option
111  segen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
112  segen->set_vertex_distribution_width(0.0, 0.0, 5.0);
114  segen->set_vertex_size_parameters(0.0, 0.0);
115  segen->set_eta_range(-4.,4.);
116  segen->set_phi_range(-TMath::Pi(),TMath::Pi());
117  segen->set_p_range(.00001,30.);
118  //segen->Embed(2);
119  segen->Verbosity(0);
120  // --------------------------------------------------------------------------------------
121  // Pythia 8 events
122  bool do_pythia8_jets = false;
123  if (particle_gen==1){se->registerSubsystem( gen); cout << "Using particle generator" << endl;}
124  else if(particle_gen==2){se->registerSubsystem( gun); cout << "Using particle gun" << endl;}
125  else if(particle_gen==3){se->registerSubsystem(segen); cout << "Using simple event generator" << endl;}
126  else if(particle_gen==4){
127  do_pythia8_jets = true;
128  gSystem->Load("libPHPythia8.so");
129 
130  PHPythia8 *pythia8 = new PHPythia8(); // see coresoftware/generators/PHPythia8 for example config
131  pythia8->set_config_file("phpythia8.cfg"); // example configure files : https://github.com/sPHENIX-Collaboration/coresoftware/tree/master/generators/PHPythia8
132  se->registerSubsystem(pythia8);
133 
134  // 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
135  HepMCNodeReader *hr = new HepMCNodeReader();
136  se->registerSubsystem(hr);
137  }
138  else if (particle_gen==5){se->registerSubsystem(gen_pT); cout << "Using particle generator flat in pT" << endl;}
139  // ======================================================================================================
140  PHG4Reco *g4Reco = new PHG4Reco();
141  //g4Reco->SetWorldMaterial("G4_Galactic");
142  // ======================================================================================================
143  // Magnetic field setting
144  TString B_label;
145  if(magnetic_field==1){ // uniform 1.5T
146  B_label = "_B_1.5T";
147  g4Reco->set_field(1.5);
148  }
149  else if(magnetic_field==2){ // uniform 3.0T
150  B_label = "_B_3.0T";
151  g4Reco->set_field(3.0);
152  }
153  else if(magnetic_field==3){ // sPHENIX 1.4T map
154  B_label = "_sPHENIX";
155  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"), PHFieldConfig::kField2D);
156  g4Reco->set_field_rescale(-1.4/1.5);
157  }
158  else if(magnetic_field==4){ // Beast 3.0T map
159  B_label = "_Beast";
160  g4Reco->set_field_map(string(getenv("CALIBRATIONROOT")) + string("/Field/Map/mfield.4col.dat"), PHFieldConfig::kFieldBeast);
161  }
162  else{ // The user did not provide a valid B field setting
163  cout << "User did not provide a valid magnetic field setting. Set 'magnetic_field'. Bailing out!" << endl;
164  }
165  // ======================================================================================================
166  // Physics list (default list is "QGSP_BERT")
167  //g4Reco->SetPhysicsList("FTFP_BERT_HP"); // This list is slower and only useful for hadronic showers.
168  // ======================================================================================================
169  load_AllSi_geom(g4Reco, det_ver); // Loading All-Si Tracker and beampipe geometries
170 
171  // ======================================================================================================
172  if(do_projections){
174  cyl = new PHG4CylinderSubsystem(projname1,0);
175  cyl->set_double_param("length", length1);
176  cyl->set_double_param("radius", projradius1); // dirc radius
177  cyl->set_double_param("thickness", 0.1); // needs some thickness
178  cyl->set_string_param("material", "G4_AIR");
179  cyl->SetActive(1);
180  cyl->SuperDetector(projname1);
181  cyl->BlackHole();
182  cyl->set_color(1,0,0,0.7); //reddish
183  g4Reco->registerSubsystem(cyl);
184 
185  cyl = new PHG4CylinderSubsystem(projname2,0);
186  cyl->set_double_param("length", thinness);
187  cyl->set_double_param("radius", 2); // beampipe needs to fit here
188  cyl->set_double_param("thickness", projradius2); //
189  cyl->set_string_param("material", "G4_AIR");
190  cyl->set_double_param("place_z", projzpos2);
191  cyl->SetActive(1);
192  cyl->SuperDetector(projname2);
193  cyl->BlackHole();
194  cyl->set_color(0,1,1,0.3); //reddish
195  g4Reco->registerSubsystem(cyl);
196 
197  cyl = new PHG4CylinderSubsystem(projname3,0);
198  cyl->set_double_param("length", thinness);
199  cyl->set_double_param("radius", 2); // beampipe needs to fit here
200  cyl->set_double_param("thickness", projradius3); //
201  cyl->set_string_param("material", "G4_AIR");
202  cyl->set_double_param("place_z", projzpos3);
203  cyl->SetActive(1);
204  cyl->SuperDetector(projname3);
205  cyl->BlackHole();
206  cyl->set_color(0,1,1,0.3); //reddish
207  g4Reco->registerSubsystem(cyl);
208  }
209  // ======================================================================================================
210 
211  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
212  g4Reco->registerSubsystem(truth);
213 
214  se->registerSubsystem(g4Reco);
215 
216  // ======================================================================================================
217  if(do_pythia8_jets){
218  Bbc_Reco();
219  Global_Reco();
220  }
221 
222  // ======================================================================================================
223  // Fast pattern recognition and full Kalman filter
224  // output evaluation file for truth track and reco tracks are PHG4TruthInfoContainer
225  double um_to_cm = 1E-04; // Conversion factor from um to cm
226  char nodename[100];
227  PHG4TrackFastSim *kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
228  kalman->set_use_vertex_in_fitting(false);
229  kalman->set_sub_top_node_name("SVTX");
230  kalman->set_trackmap_out_name("SvtxTrackMap");
231 
232  add_AllSi_to_kalman( kalman , pixel_size , det_ver ); // Add All-Silicon tracker to Kalman filter
233 
234  // Projections
235  if(do_projections){
236  kalman->add_cylinder_state(projname1, projradius1); // projection on cylinder (DIRC)
237  kalman->add_zplane_state (projname2, projzpos2 ); // projection on vertical planes
238  kalman->add_zplane_state (projname3, projzpos3 ); // projection on vertical planes
239  }
240 
241  if(DISPLACED_VERTEX){
242  // use very loose vertex constraint (1cm in sigma) to allow reco of displaced vertex
243  kalman->set_use_vertex_in_fitting(true);
244  kalman->set_vertex_xy_resolution(1);
245  kalman->set_vertex_z_resolution(1);
246  kalman->enable_vertexing(true);
247  kalman->set_vertex_min_ndf(2); // Only tracks with number of degrees of freedom greater than this value are used to fit the vertex
248  }
249  /*
250  else{
251  kalman->set_use_vertex_in_fitting(false);
252  kalman->set_vertex_xy_resolution(0);
253  kalman->set_vertex_z_resolution(0);
254  kalman->enable_vertexing(false); // this is false by default
255  kalman->set_vertex_min_ndf(2);
256  }
257  */
258  //kalman -> Verbosity(10);
259  se->registerSubsystem(kalman);
260  // -----------------------------------------------------
261  // INFO: The resolution numbers above correspond to:
262  // 20e-4/sqrt(12) cm = 5.8e-4 cm, to simulate 20x20 um
263 
264  // ======================================================================================================
265  TrackFastSimEval *fast_sim_eval = new TrackFastSimEval("FastTrackingEval");
266  fast_sim_eval->set_filename(TString(outputFile)+B_label+"_FastTrackingEval.root");
267  if(do_projections){
268  fast_sim_eval->AddProjection(projname1);
269  fast_sim_eval->AddProjection(projname2);
270  fast_sim_eval->AddProjection(projname3);
271  }
272  se->registerSubsystem(fast_sim_eval);
273 
274  // ======================================================================================================
275  // resonstruct jets after the tracking
276  if(do_pythia8_jets) Jet_Reco();
277 
278  SimpleNtuple *hits = new SimpleNtuple("Hits");
279  add_AllSi_hits(hits,det_ver); // Add All-Silicon tracker hits
280  se->registerSubsystem(hits);
281 
282  // ======================================================================================================
283  if(do_pythia8_jets)
284  Jet_Eval(string(outputFile) + "_g4jet_eval.root");
285  // ======================================================================================================
286 
288  // IOManagers...
290  const std::string dst_name = std::string(outputFile)+"_G4LBLVtx.root";
291  //Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",TString(outputFile)+"_G4LBLVtx.root");
292  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",dst_name);
293  out->Verbosity(0);
294  se->registerOutputManager(out);
295 
297  se->registerInputManager(in);
298 
299  if(do_pythia8_jets){
300  gSystem->Load("libmyjetanalysis");
301  std::string jetoutputFile = std::string(outputFile) + std::string("_electrons+jets.root");
302  MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r10","AntiKt_Truth_r10",jetoutputFile.data());
303  //MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r04","AntiKt_Truth_r04",jetoutputFile.data());
304  //MyJetAnalysis_AllSi *myJetAnalysis = new MyJetAnalysis_AllSi("AntiKt_Track_r02","AntiKt_Truth_r02",jetoutputFile.data());
305  se->registerSubsystem(myJetAnalysis);
306  }
307 
308  if (nEvents <= 0) return;
309 
310  se->run(nEvents);
311  se->End();
312  delete se;
313 
314  gSystem->Exit(0);
315 }