Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_TPC.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_TPC.C
1 #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 00, 0)
2 
3 #include <fun4all/Fun4AllDstOutputManager.h>
4 #include <fun4all/Fun4AllDummyInputManager.h>
5 #include <fun4all/Fun4AllInputManager.h>
6 #include <fun4all/Fun4AllOutputManager.h>
7 #include <fun4all/Fun4AllServer.h>
8 #include <fun4all/SubsysReco.h>
9 
10 #include <fun4all/Fun4AllServer.h>
11 
12 #include <g4eval/SvtxEvaluator.h>
13 #include <g4eval/TrkrEvaluator.h>
14 
15 #include <g4detectors/PHG4BlockSubsystem.h>
16 #include <g4main/PHG4Reco.h>
17 #include <phgeom/PHGeomFileImport.h>
18 
19 #include <g4tpc/PHG4TpcDigitizer.h>
20 #include <g4tpc/PHG4TpcElectronDrift.h>
21 #include <g4tpc/PHG4TpcPadPlane.h>
22 #include <g4tpc/PHG4TpcPadPlaneReadout.h>
23 #include <g4tpc/PHG4TpcSubsystem.h>
24 #include <trackreco/PHGenFitTrkFitter.h>
25 #include <trackreco/PHGenFitTrkProp.h>
26 #include <trackreco/PHHoughSeeding.h>
27 #include <trackreco/PHInitVertexing.h>
28 #include <trackreco/PHTrackSeeding.h>
29 #include <trackreco/PHTruthTrackSeeding.h>
30 #include <trackreco/PHTruthVertexing.h>
31 
32 #include <g4eval/PHG4DSTReader.h>
33 
34 #include <g4histos/G4HitNtuple.h>
35 
36 #include <g4main/PHG4ParticleGenerator.h>
37 #include <g4main/PHG4ParticleGun.h>
38 #include <g4main/PHG4Reco.h>
39 #include <g4main/PHG4SimpleEventGenerator.h>
40 #include <g4main/PHG4TruthSubsystem.h>
41 
42 #include <phool/recoConsts.h>
43 
44 #include <tpc2019/TpcPrototypeClusterizer.h>
45 #include <tpc2019/TpcPrototypeGenFitTrkFinder.h>
46 #include <tpc2019/TpcPrototypeGenFitTrkFitter.h>
47 
48 R__LOAD_LIBRARY(libcalo_reco.so)
49 R__LOAD_LIBRARY(libfun4all.so)
50 R__LOAD_LIBRARY(libg4tpc.so)
51 R__LOAD_LIBRARY(libg4detectors.so)
52 R__LOAD_LIBRARY(libg4eval.so)
53 R__LOAD_LIBRARY(libg4histos.so)
54 R__LOAD_LIBRARY(libg4testbench.so)
55 R__LOAD_LIBRARY(libg4tpc.so)
56 R__LOAD_LIBRARY(libg4intt.so)
57 R__LOAD_LIBRARY(libg4mvtx.so)
58 R__LOAD_LIBRARY(libg4hough.so)
59 R__LOAD_LIBRARY(libg4eval.so)
60 R__LOAD_LIBRARY(libintt.so)
61 R__LOAD_LIBRARY(libmvtx.so)
62 R__LOAD_LIBRARY(libtpc2019.so)
63 R__LOAD_LIBRARY(libtrack_reco.so)
64 #endif
65 
66 int n_tpc_layer_inner = 0;
67 int tpc_layer_rphi_count_inner = 1152;
68 int n_tpc_layer_mid = 16;
69 int n_tpc_layer_outer = 0;
70 int n_gas_layer = n_tpc_layer_inner + n_tpc_layer_mid + n_tpc_layer_outer;
71 
72 int Fun4All_G4_TPC(int nEvents = 100, bool eventDisp = false, int verbosity = 0)
73 {
74  gSystem->Load("libfun4all");
75  gSystem->Load("libg4detectors");
76  gSystem->Load("libg4testbench");
77  gSystem->Load("libg4histos");
78  gSystem->Load("libg4eval.so");
79  gSystem->Load("libqa_modules");
80  gSystem->Load("libg4tpc");
81  gSystem->Load("libtrack_io.so");
82  gSystem->Load("libfun4all.so");
83  gSystem->Load("libg4detectors.so");
84  gSystem->Load("libtpc2019.so");
85  gSystem->Load("libg4eval.so");
86  gSystem->Load("libfun4all.so");
87  gSystem->Load("libg4detectors.so");
88  gSystem->Load("libg4hough.so");
89  gSystem->Load("libtrack_reco.so");
90 
91  bool bh_on = false; // the surrounding boxes need some further thinking
92  bool dstreader = true;
93  bool dstoutput = false;
94 
95  const double TPCDriftLength = 40;
96 
98  // Make the Server
100  Fun4AllServer *se = Fun4AllServer::instance();
101  se->Verbosity(1);
102  recoConsts *rc = recoConsts::instance();
103  // only set this if you want a fixed random seed to make
104  // results reproducible for testing
105  rc->set_IntFlag("RANDOMSEED", 12345678);
106 
107  // simulated setup sits at eta=1, theta=40.395 degrees
108  double theta = 90 - 5;
109  double phi = 180 + 360 / 12 / 2;
110  // shift in x with respect to midrapidity setup
111  double add_place_z = -TPCDriftLength * .5;
112  // Test beam generator
113  PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
114  gen->add_particles("proton", 1); // mu-,e-,anti_proton,pi-
115  gen->set_vertex_distribution_mean(0.0, 0.0, add_place_z);
116  gen->set_vertex_distribution_width(0.0, .7, .7); // Rough beam profile size @ 16 GeV measured by Abhisek
117  gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
118  PHG4SimpleEventGenerator::Gaus,
119  PHG4SimpleEventGenerator::Gaus); // Gauss beam profile
120  double angle = theta * TMath::Pi() / 180.;
121  double eta = -1. * TMath::Log(TMath::Tan(angle / 2.));
122  gen->set_eta_range(eta - 0.001, eta + 0.001); // 1mrad angular divergence
123  gen->set_phi_range(TMath::Pi() * phi / 180 - 0.001, TMath::Pi() * phi / 180 + 0.001); // 1mrad angular divergence
124  const double momentum = 120;
125  gen->set_p_range(momentum, momentum, momentum * 2e-2); // 2% momentum smearing
126  se->registerSubsystem(gen);
127 
128  PHG4Reco *g4Reco = new PHG4Reco();
129  g4Reco->set_field(0);
130  // g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
131 
132  //----------------------------------------
133  // TPC G4
134  //----------------------------------------
135 
136  PHG4TpcSubsystem *tpc = new PHG4TpcSubsystem("TPC");
137  tpc->SetActive();
138  tpc->SuperDetector("TPC");
139  tpc->set_double_param("steplimits", 1);
140  tpc->set_double_param("tpc_length", TPCDriftLength * 2); // 2x 40 cm drift
141  // By default uses "sPHENIX_TPC_Gas", defined in PHG4Reco. That is 90:10 Ne:C4
142  tpc->SetAbsorberActive();
143  g4Reco->registerSubsystem(tpc);
144 
145  if (bh_on)
146  {
147  // BLACKHOLE, box surrounding the prototype to check for leakage
148  PHG4BlockSubsystem *bh[5];
149  // surrounding outer hcal
150  // top
151  bh[0] = new PHG4BlockSubsystem("bh1", 1);
152  bh[0]->set_double_param("size_x", 270.);
153  bh[0]->set_double_param("size_y", 0.01);
154  bh[0]->set_double_param("size_z", 165.);
155  bh[0]->set_double_param("place_x", 270. / 2.);
156  bh[0]->set_double_param("place_y", 125. / 2.);
157  // bottom
158  bh[1] = new PHG4BlockSubsystem("bh2", 2);
159  bh[1]->set_double_param("size_x", 270.);
160  bh[1]->set_double_param("size_y", 0.01);
161  bh[1]->set_double_param("size_z", 165.);
162  bh[1]->set_double_param("place_x", 270. / 2.);
163  bh[1]->set_double_param("place_y", -125. / 2.);
164  // right side
165  bh[2] = new PHG4BlockSubsystem("bh3", 3);
166  bh[2]->set_double_param("size_x", 200.);
167  bh[2]->set_double_param("size_y", 125.);
168  bh[2]->set_double_param("size_z", 0.01);
169  bh[2]->set_double_param("place_x", 200. / 2.);
170  bh[2]->set_double_param("place_z", 165. / 2.);
171  // left side
172  bh[3] = new PHG4BlockSubsystem("bh4", 4);
173  bh[3]->set_double_param("size_x", 270.);
174  bh[3]->set_double_param("size_y", 125.);
175  bh[3]->set_double_param("size_z", 0.01);
176  bh[3]->set_double_param("place_x", 270. / 2.);
177  bh[3]->set_double_param("place_z", -165. / 2.);
178  // back
179  bh[4] = new PHG4BlockSubsystem("bh5", 5);
180  bh[4]->set_double_param("size_x", 0.01);
181  bh[4]->set_double_param("size_y", 125.);
182  bh[4]->set_double_param("size_z", 165.);
183  bh[4]->set_double_param("place_x", 270.);
184  for (int i = 0; i < 5; i++)
185  {
186  bh[i]->BlackHole();
187  bh[i]->SetActive();
188  bh[i]->SuperDetector("BlackHole");
189  bh[i]->OverlapCheck(true);
190  g4Reco->registerSubsystem(bh[i]);
191  }
192  }
193  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
194  g4Reco->registerSubsystem(truth);
195 
196  g4Reco->save_DST_geometry(false);
197  se->registerSubsystem(g4Reco);
198 
199  //swap out with the test beam geometry for the analysis stage
200  PHGeomFileImport *import = new PHGeomFileImport("TpcPrototypeGeometry.gdml");
201  se->registerSubsystem(import);
202 
203  //=========================
204  // setup Tpc readout for filling cells
205  // g4tpc/PHG4TpcElectronDrift uses
206  // g4tpc/PHG4TpcPadPlaneReadout
207  //=========================
208 
209  PHG4TpcPadPlane *padplane = new PHG4TpcPadPlaneReadout();
210 
211  // The pad plane readout default is set in PHG4TpcPadPlaneReadout
212  //only build mid-layer and to size
213  padplane->set_int_param("tpc_minlayer_inner", 0);
214  padplane->set_int_param("ntpc_layers_inner", 0);
215  padplane->set_int_param("ntpc_layers_outer", 0);
216  padplane->set_double_param("maxdriftlength", TPCDriftLength);
217  padplane->set_int_param("ntpc_phibins_mid", 16 * 8 * 12);
218 
219  PHG4TpcElectronDrift *edrift = new PHG4TpcElectronDrift();
220  edrift->Detector("TPC");
221  // fudge factors to get drphi 150 microns (in mid and outer Tpc) and dz 500 microns cluster resolution
222  // They represent effects not due to ideal gas properties and ideal readout plane behavior
223  // defaults are 0.12 and 0.15, they can be changed here to get a different resolution
224  edrift->set_double_param("added_smear_trans", 0.12);
225  edrift->set_double_param("added_smear_long", 0.15);
226  edrift->registerPadPlane(padplane);
227  edrift->Verbosity(verbosity);
228  ;
229  se->registerSubsystem(edrift);
230 
231  // Tpc
232  //====
233  PHG4TpcDigitizer *digitpc = new PHG4TpcDigitizer();
234  digitpc->SetTpcMinLayer(0);
235  double ENC = 670.0; // standard
236  digitpc->SetENC(ENC);
237  double ADC_threshold = 4.0 * ENC;
238  digitpc->SetADCThreshold(ADC_threshold); // 4 * ENC seems OK
239 
240  cout << " Tpc digitizer: Setting ENC to " << ENC << " ADC threshold to " << ADC_threshold << endl;
241 
242  se->registerSubsystem(digitpc);
243 
244  // For the Tpc
245  //==========
246  TpcPrototypeClusterizer *tpcclusterizer = new TpcPrototypeClusterizer();
247  tpcclusterizer->Verbosity(verbosity);
248  ;
249  se->registerSubsystem(tpcclusterizer);
250 
251  //-------------
252  // Tracking
253  //------------
254 
255  // This should be true for everything except testing wirh truth track seeding!
256  const bool use_track_prop = true;
257  if (use_track_prop)
258  {
259  // Find all clusters associated with each seed track
260  TpcPrototypeGenFitTrkFinder *finder = new TpcPrototypeGenFitTrkFinder();
261  finder->Verbosity(verbosity);
262  finder->set_do_evt_display(eventDisp);
263  finder->set_do_eval(true);
264  se->registerSubsystem(finder);
265  }
266  else
267  {
268  //--------------------------------------------------
269  // Track finding using truth information
270  //--------------------------------------------------
271 
272  PHInitVertexing *init_vtx = new PHTruthVertexing("PHTruthVertexing");
273  init_vtx->Verbosity(verbosity);
274  ;
275  se->registerSubsystem(init_vtx);
276 
277  // For each truth particle, create a track and associate clusters with it using truth information
278  PHTruthTrackSeeding *pat_rec = new PHTruthTrackSeeding("PHTruthTrackSeeding");
279  pat_rec->Verbosity(verbosity);
280  ;
281  pat_rec->set_min_clusters_per_track(10);
282  se->registerSubsystem(pat_rec);
283  }
284  //
285  //------------------------------------------------
286  // Fitting of tracks using Kalman Filter
287  //------------------------------------------------
288 
289  TpcPrototypeGenFitTrkFitter *kalman = new TpcPrototypeGenFitTrkFitter();
290  kalman->Verbosity(verbosity);
291  kalman->set_do_evt_display(eventDisp);
292  kalman->set_do_eval(true);
293  se->registerSubsystem(kalman);
294 
295  //----------------
296  // Tracking evaluation
297  //----------------
298 
299  SvtxEvaluator *eval;
300  eval = new SvtxEvaluator("SVTXEVALUATOR", "G4TPC_eval.root", "SvtxTrackMap", 0, 0, n_gas_layer);
301  eval->do_cluster_eval(true);
302  eval->do_g4hit_eval(true);
303  eval->do_hit_eval(true); // enable to see the hits that includes the chamber physics...
304  eval->do_gpoint_eval(true);
305  eval->do_eval_light(false);
306  eval->scan_for_embedded(false); // take all tracks if false - take only embedded tracks if true
307  eval->Verbosity(verbosity);
308  ;
309  se->registerSubsystem(eval);
310 
311  //----------------------
312  // save a comprehensive evaluation file
313  //----------------------
314  if (dstreader)
315  {
316  PHG4DSTReader *ana = new PHG4DSTReader(string("G4TPC_DSTReader.root"));
317  ana->set_save_particle(true);
318  ana->set_load_all_particle(true);
319  ana->set_load_active_particle(true);
320  ana->set_save_vertex(true);
321 
322  ana->AddNode("ABSORBER_TPC");
323  ana->AddNode("TPC");
324  if (bh_on)
325  ana->AddNode("BlackHole"); // add a G4Hit node
326 
327  se->registerSubsystem(ana);
328  }
329 
330  if (dstoutput)
331  {
332  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", "G4TPC.root");
333  se->registerOutputManager(out);
334  }
335 
336  Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
337  se->registerInputManager(in);
338  if (nEvents <= 0)
339  {
340  return 0;
341  }
342 
343  gSystem->ListLibraries();
344  se->run(nEvents);
345 
346  se->End();
347 
348  // std::cout << "All done" << std::endl;
349  delete se;
350  // return 0;
351  gSystem->Exit(0);
352  return 0;
353 }
354 
355 // for using QuickTest to check if macro loads
356 void RunLoadTest() {}