Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_G4_Prototype2.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_G4_Prototype2.C
1 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
2 #include <caloreco/RawTowerCalibration.h>
3 
4 #include <fun4all/SubsysReco.h>
5 #include <fun4all/Fun4AllDstOutputManager.h>
6 #include <fun4all/Fun4AllDummyInputManager.h>
7 #include <fun4all/Fun4AllInputManager.h>
8 #include <fun4all/Fun4AllOutputManager.h>
9 #include <fun4all/Fun4AllServer.h>
10 
11 #include <g4calo/RawTowerBuilder.h>
12 #include <g4calo/RawTowerDigitizer.h>
13 
14 #include <g4caloprototype/Prototype2RawTowerBuilder.h>
15 #include <g4detectors/PHG4BlockSubsystem.h>
16 #include <g4detectors/PHG4FullProjSpacalCellReco.h>
17 #include <g4caloprototype/PHG4Prototype2HcalCellReco.h>
18 #include <g4caloprototype/PHG4Prototype2InnerHcalSubsystem.h>
19 #include <g4caloprototype/PHG4Prototype2OuterHcalSubsystem.h>
20 #include <g4caloprototype/PHG4SpacalPrototypeSubsystem.h>
21 
22 #include <g4eval/PHG4DSTReader.h>
23 
24 #include <g4histos/G4HitNtuple.h>
25 
26 #include <g4main/PHG4ParticleGun.h>
27 #include <g4main/PHG4Reco.h>
28 #include <g4main/PHG4SimpleEventGenerator.h>
29 #include <g4main/PHG4TruthSubsystem.h>
30 
31 #include <phool/recoConsts.h>
32 
33 #include <qa_modules/QAHistManagerDef.h>
34 #include <qa_modules/QAG4SimulationCalorimeter.h>
35 
36 R__LOAD_LIBRARY(libcalo_reco.so)
37 R__LOAD_LIBRARY(libfun4all.so)
38 R__LOAD_LIBRARY(libg4calo.so)
39 R__LOAD_LIBRARY(libg4caloprototype.so)
40 R__LOAD_LIBRARY(libg4eval.so)
41 R__LOAD_LIBRARY(libg4histos.so)
42 R__LOAD_LIBRARY(libg4testbench.so)
43 R__LOAD_LIBRARY(libqa_modules.so)
44 
45 #endif
46 
47 int Fun4All_G4_Prototype2(int nEvents = 1)
48 {
49 
50  gSystem->Load("libfun4all");
51  gSystem->Load("libg4caloprototype");
52  gSystem->Load("libg4testbench");
53  gSystem->Load("libg4histos");
54  gSystem->Load("libg4eval");
55  gSystem->Load("libqa_modules");
56  gSystem->Load("libg4calo");
57  gSystem->Load("libcalo_reco");
58 
59  bool cemc_on = true;
60  bool cemc_cell = cemc_on && true;
61  bool cemc_twr = cemc_cell && true;
62  bool cemc_digi = cemc_twr && true;
63  bool cemc_twrcal = cemc_digi && true;
64 
65  bool ihcal_on = true;
66  bool ihcal_cell = ihcal_on && true;
67  bool ihcal_twr = ihcal_cell && true;
68  bool ihcal_digi = ihcal_twr && true;
69  bool ihcal_twrcal = ihcal_digi && true;
70 
71  bool ohcal_on = true;
72  bool ohcal_cell = ohcal_on && true;
73  bool ohcal_twr = ohcal_cell && true;
74  bool ohcal_digi = ohcal_twr && true;
75  bool ohcal_twrcal = ohcal_digi && true;
76 
77  bool cryo_on = true;
78  bool bh_on = true;
79  bool hit_ntuple = false;
80  bool dstreader = true;
81  bool dstoutput = false;
82 
84  // Make the Server
86  Fun4AllServer *se = Fun4AllServer::instance();
87  // se->Verbosity(1);
88  recoConsts *rc = recoConsts::instance();
89  // only set this if you want a fixed random seed to make
90  // results reproducible for testing
91  // rc->set_IntFlag("RANDOMSEED",12345);
92 
93  // Test beam generator
94  PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
95  gen->add_particles("pi-", 1); // mu-,e-,anti_proton,pi-
96  gen->set_vertex_distribution_mean(0.0, 0.0, 0);
97  gen->set_vertex_distribution_width(0.0, .7, .7); // Rough beam profile size @ 16 GeV measured by Abhisek
98  gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
99  PHG4SimpleEventGenerator::Gaus, PHG4SimpleEventGenerator::Gaus); // Gauss beam profile
100  gen->set_eta_range(-.001, .001); // 1mrad angular divergence
101  gen->set_phi_range(-.001, .001); // 1mrad angular divergence
102  const double momentum = 32;
103  gen->set_p_range(momentum,momentum, momentum*2e-2); // 2% momentum smearing
104  se->registerSubsystem(gen);
105 
106  // Simple single particle generator
107  PHG4ParticleGun *gun = new PHG4ParticleGun();
108  // gun->set_name("anti_proton");
109  // gun->set_name("geantino");
110  gun->set_name("proton");
111  gun->set_vtx(0, 0, 0);
112  gun->set_mom(120, 0, 0);
113  // gun->AddParticle("geantino",1.7776,-0.4335,0.);
114  // gun->AddParticle("geantino",1.7709,-0.4598,0.);
115  // gun->AddParticle("geantino",2.5621,0.60964,0.);
116  // gun->AddParticle("geantino",1.8121,0.253,0.);
117  // se->registerSubsystem(gun);
118 
119  PHG4Reco* g4Reco = new PHG4Reco();
120  g4Reco->set_field(0);
121  // g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
122 
123  //----------------------------------------
124  // EMCal G4
125  //----------------------------------------
126  if (cemc_on)
127  {
128  PHG4SpacalPrototypeSubsystem *cemc = new PHG4SpacalPrototypeSubsystem("CEMC");
129  cemc->SetActive();
130  cemc->SuperDetector("CEMC");
131  cemc->SetAbsorberActive();
132  cemc->OverlapCheck(true);
133  cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
134  cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/Prototype2/Geometry/") );
135  // cemc->set_double_param("z_rotation_degree", 15); // rotation around CG
136  cemc->set_double_param("xpos", (116.77 + 137.0)*.5 - 26.5 - 10.2); // location in cm of EMCal CG. Updated with final positioning of EMCal
137  cemc->set_double_param("ypos", 4); // put it some where in UIUC blocks
138  cemc->set_double_param("zpos", 4); // put it some where in UIUC blocks
139  g4Reco->registerSubsystem(cemc);
140  }
141  //----------------------------------------
142  // HCal G4
143  //----------------------------------------
144  if (ihcal_on)
145  {
146  PHG4Prototype2InnerHcalSubsystem *innerhcal = new PHG4Prototype2InnerHcalSubsystem("HCalIn");
147  innerhcal->SetActive();
148  innerhcal->SetAbsorberActive();
149  innerhcal->SetAbsorberTruth(1);
150  innerhcal->OverlapCheck(true);
151  innerhcal->SuperDetector("HCALIN");
152  g4Reco->registerSubsystem(innerhcal);
153  }
154  if (ohcal_on)
155  {
156  PHG4Prototype2OuterHcalSubsystem *outerhcal = new PHG4Prototype2OuterHcalSubsystem("HCalOut");
157  outerhcal->SetActive();
158  outerhcal->SetAbsorberActive();
159  outerhcal->SetAbsorberTruth(1);
160  outerhcal->OverlapCheck(true);
161  outerhcal->SuperDetector("HCALOUT");
162  g4Reco->registerSubsystem(outerhcal);
163  }
164  // Cryostat from engineering drawing
165  if (cryo_on)
166  {
167  PHG4BlockSubsystem *cryo1 = new PHG4BlockSubsystem("cryo1",1);
168  cryo1->set_double_param("size_x",0.95);
169  cryo1->set_double_param("size_y",60.96);
170  cryo1->set_double_param("size_z",60.96);
171  cryo1->set_double_param("place_x",141.96+0.95/2.);
172  cryo1->set_string_param("material","G4_Al");
173  cryo1->SetActive(); // it is an active volume - save G4Hits
174  cryo1->SuperDetector("CRYO");
175  g4Reco->registerSubsystem(cryo1);
176 
177  PHG4BlockSubsystem *cryo2 = new PHG4BlockSubsystem("cryo2",2);
178  cryo2->set_double_param("size_x",8.89);
179  cryo2->set_double_param("size_y",60.96);
180  cryo2->set_double_param("size_z",60.96);
181  cryo2->set_double_param("place_x",150.72+8.89/2.);
182  cryo2->set_string_param("material","G4_Al");
183  cryo2->SetActive(); // it is an active volume - save G4Hits
184  cryo2->SuperDetector("CRYO");
185  g4Reco->registerSubsystem(cryo2);
186 
187  PHG4BlockSubsystem *cryo3 = new PHG4BlockSubsystem("cryo3",3);
188  cryo3->set_double_param("size_x",2.54);
189  cryo3->set_double_param("size_y",60.96);
190  cryo3->set_double_param("size_z",60.96);
191  cryo3->set_double_param("place_x",173.93+2.54/2.);
192  cryo3->set_string_param("material","G4_Al");
193  cryo3->SetActive(); // it is an active volume - save G4Hits
194  cryo3->SuperDetector("CRYO");
195  g4Reco->registerSubsystem(cryo3);
196  }
197  // BLACKHOLE, box surrounding the prototype to check for leakage
198  if (bh_on)
199  {
200  PHG4BlockSubsystem *bh[5];
201  // surrounding outer hcal
202  // top
203  bh[0] = new PHG4BlockSubsystem("bh1",1);
204  bh[0]->set_double_param("size_x",270.);
205  bh[0]->set_double_param("size_y",0.01);
206  bh[0]->set_double_param("size_z",165.);
207  bh[0]->set_double_param("place_x",270./2.);
208  bh[0]->set_double_param("place_y",125./2.);
209  // bottom
210  bh[1] = new PHG4BlockSubsystem("bh2",2);
211  bh[1]->set_double_param("size_x",270.);
212  bh[1]->set_double_param("size_y",0.01);
213  bh[1]->set_double_param("size_z",165.);
214  bh[1]->set_double_param("place_x",270./2.);
215  bh[1]->set_double_param("place_y",-125./2.);
216  // right side
217  bh[2] = new PHG4BlockSubsystem("bh3",3);
218  bh[2]->set_double_param("size_x",270.);
219  bh[2]->set_double_param("size_y",125.);
220  bh[2]->set_double_param("size_z",0.01);
221  bh[2]->set_double_param("place_x",270./2.);
222  bh[2]->set_double_param("place_z",165./2.);
223  // left side
224  bh[3] = new PHG4BlockSubsystem("bh4",4);
225  bh[3]->set_double_param("size_x",270.);
226  bh[3]->set_double_param("size_y",125.);
227  bh[3]->set_double_param("size_z",0.01);
228  bh[3]->set_double_param("place_x",270./2.);
229  bh[3]->set_double_param("place_z",-165./2.);
230  // back
231  bh[4] = new PHG4BlockSubsystem("bh5",5);
232  bh[4]->set_double_param("size_x",0.01);
233  bh[4]->set_double_param("size_y",125.);
234  bh[4]->set_double_param("size_z",165.);
235  bh[4]->set_double_param("place_x",270.);
236  for (int i=0; i<5; i++)
237  {
238  bh[i]->BlackHole();
239  bh[i]->SetActive();
240  bh[i]->SuperDetector("BlackHole");
241  bh[i]->OverlapCheck(true);
242  g4Reco->registerSubsystem(bh[i]);
243  }
244  }
245  PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
246  g4Reco->registerSubsystem(truth);
247 
248  se->registerSubsystem( g4Reco );
249  //----------------------------------------
250  // EMCal digitization
251  //----------------------------------------
252 
253  if (cemc_cell)
254  {
255  PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
256  cemc_cells->Detector("CEMC");
257  cemc_cells->set_timing_window(0.,60.);
258  cemc_cells->get_light_collection_model().load_data_file(
259  string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype2Module.xml"),
260  "data_grid_light_guide_efficiency","data_grid_fiber_trans");
261 
262  se->registerSubsystem(cemc_cells);
263  }
264  if (cemc_twr)
265  {
266  RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
267  TowerBuilder->Detector("CEMC");
268  TowerBuilder->set_sim_tower_node_prefix("SIM");
269  se->registerSubsystem(TowerBuilder);
270  }
271 
272 // these parameters are used in multiple modules
273  const double ADC_per_photoelectron_LG = 0.24; // From Sean Stoll, Mar 29
274  const double ADC_per_photoelectron_HG = 3.8; // From Sean Stoll, Mar 29
275  const double photoelectron_per_GeV = 500; //500 photon per total GeV deposition
276 
277  if (cemc_digi)
278  {
279  const double sampling_fraction = 0.0233369; // +/- 8.22211e-05 from 15 Degree indenting 8 GeV electron showers
280 
281  // low gains
282  RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerLG");
283  TowerDigitizer->Detector("CEMC");
284  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
285  TowerDigitizer->set_digi_algorithm(
286  RawTowerDigitizer::kSimple_photon_digitalization);
287  TowerDigitizer->set_pedstal_central_ADC(0);
288  TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
289  TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_LG);
290  TowerDigitizer->set_photonelec_yield_visible_GeV(
291  photoelectron_per_GeV / sampling_fraction);
292  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
293  se->registerSubsystem(TowerDigitizer);
294  // high gains
295  TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerHG");
296  TowerDigitizer->Detector("CEMC");
297  TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
298  TowerDigitizer->set_digi_algorithm(
299  RawTowerDigitizer::kSimple_photon_digitalization);
300  TowerDigitizer->set_pedstal_central_ADC(0);
301  TowerDigitizer->set_pedstal_width_ADC(15); // From John Haggerty, Mar 29
302  TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_HG);
303  TowerDigitizer->set_photonelec_yield_visible_GeV(
304  photoelectron_per_GeV / sampling_fraction);
305  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
306  se->registerSubsystem(TowerDigitizer);
307  }
308  if (cemc_twrcal)
309  {
310  RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationLG");
311  TowerCalibration->Detector("CEMC");
312  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
313  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
314  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
315  TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_LG / photoelectron_per_GeV);
316  TowerCalibration->set_pedstal_ADC(0);
317  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
318  se->registerSubsystem(TowerCalibration);
319 
320  TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationHG");
321  TowerCalibration->Detector("CEMC");
322  TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
323  TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
324  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
325  TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_HG / photoelectron_per_GeV);
326  TowerCalibration->set_pedstal_ADC(0);
327  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
328  se->registerSubsystem(TowerCalibration);
329  }
330 
331  //----------------------------------------
332  // HCal towering
333  //----------------------------------------
334  if (ihcal_cell)
335  {
336  PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALinCellReco");
337  hccell->Detector("HCALIN");
338  se->registerSubsystem(hccell);
339  }
340 
341  if (ohcal_cell)
342  {
343  PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALoutCellReco");
344  hccell->Detector("HCALOUT");
345  se->registerSubsystem(hccell);
346  }
347  if (ihcal_twr)
348  {
349  Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALinRawTowerBuilder");
350  hcaltwr->Detector("HCALIN");
351  hcaltwr->set_sim_tower_node_prefix("SIM");
352  se->registerSubsystem(hcaltwr);
353  }
354  if (ohcal_twr)
355  {
356  Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALoutRawTowerBuilder");
357  hcaltwr->Detector("HCALOUT");
358  hcaltwr->set_sim_tower_node_prefix("SIM");
359  se->registerSubsystem(hcaltwr);
360  }
361 
362  //----------------------------------------
363  // HCal digitization
364  //----------------------------------------
365  // From: Abhisek Sen [mailto:sen.abhisek@gmail.com]
366  // Sent: Tuesday, April 19, 2016 10:55 PM
367  // To: Huang, Jin <jhuang@bnl.gov>; Haggerty, John <haggerty@bnl.gov>
368 
369  // HCALIN:
370  // 1/5 pixel / HG ADC channel
371  // 32/5 pixel / LG ADC channel
372  // 0.4 MeV/ LG ADC
373  // 0.4/32 MeV/ HG ADC
374 
375  // HCALOUT:
376  // 1/5 pixel / HG ADC channel
377  // 16/5 pixel / LG ADC channel
378  // 0.2 MeV/ LG ADC
379  // 0.2/16 MeV/ HG ADC
380 
381  if (ihcal_digi)
382  {
383  RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerLG");
384  TowerDigitizer->Detector("HCALIN");
385  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
386  TowerDigitizer->set_digi_algorithm(
387  RawTowerDigitizer::kSimple_photon_digitalization);
388  TowerDigitizer->set_pedstal_central_ADC(0);
389  TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
390  TowerDigitizer->set_photonelec_ADC(32. / 5.);
391  TowerDigitizer->set_photonelec_yield_visible_GeV(32. / 5 / (0.4e-3));
392  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
393  se->registerSubsystem(TowerDigitizer);
394 
395  TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerHG");
396  TowerDigitizer->Detector("HCALIN");
397  TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
398  TowerDigitizer->set_digi_algorithm(
399  RawTowerDigitizer::kSimple_photon_digitalization);
400  TowerDigitizer->set_pedstal_central_ADC(0);
401  TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
402  TowerDigitizer->set_photonelec_ADC(1. / 5.);
403  TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.4e-3 / 32));
404  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
405  se->registerSubsystem(TowerDigitizer);
406  }
407  if (ohcal_digi)
408  {
409  RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerLG");
410  TowerDigitizer->Detector("HCALOUT");
411  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
412  TowerDigitizer->set_digi_algorithm(
413  RawTowerDigitizer::kSimple_photon_digitalization);
414  TowerDigitizer->set_pedstal_central_ADC(0);
415  TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
416  TowerDigitizer->set_photonelec_ADC(16. / 5.);
417  TowerDigitizer->set_photonelec_yield_visible_GeV(16. / 5 / (0.2e-3));
418  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
419  se->registerSubsystem(TowerDigitizer);
420 
421  TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerHG");
422  TowerDigitizer->Detector("HCALOUT");
423  TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
424  TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
425  TowerDigitizer->set_pedstal_central_ADC(0);
426  TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
427  TowerDigitizer->set_photonelec_ADC(1. / 5.);
428  TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.2e-3 / 16));
429  TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
430  se->registerSubsystem(TowerDigitizer);
431  }
432  //----------------------------------------
433  // HCal calibration
434  //----------------------------------------
435  // 32 GeV Pi+ scan
436  const double visible_sample_fraction_HCALIN = 7.19505e-02 ; // 1.34152e-02
437  const double visible_sample_fraction_HCALOUT = 0.0313466 ; // +/- 0.0067744
438 
439  if (ihcal_twrcal)
440  {
441  RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationLG");
442  TowerCalibration->Detector("HCALIN");
443  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
444  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
445  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
446  TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / visible_sample_fraction_HCALIN);
447  TowerCalibration->set_pedstal_ADC(0);
448  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
449  se->registerSubsystem(TowerCalibration);
450 
451  TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationHG");
452  TowerCalibration->Detector("HCALIN");
453  TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
454  TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
455  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
456  TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / 32 / visible_sample_fraction_HCALIN);
457  TowerCalibration->set_pedstal_ADC(0);
458  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
459  se->registerSubsystem(TowerCalibration);
460  }
461  if (ohcal_twrcal)
462  {
463  RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationLG");
464  TowerCalibration->Detector("HCALOUT");
465  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
466  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
467  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
468  TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / visible_sample_fraction_HCALOUT);
469  TowerCalibration->set_pedstal_ADC(0);
470  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
471  se->registerSubsystem(TowerCalibration);
472 
473  TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationHG");
474  TowerCalibration->Detector("HCALOUT");
475  TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
476  TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
477  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
478  TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / 16 / visible_sample_fraction_HCALOUT);
479  TowerCalibration->set_pedstal_ADC(0);
480  TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
481  se->registerSubsystem(TowerCalibration);
482  }
483  //----------------------
484  // QA Histograms
485  //----------------------
486  if (cemc_on)
487  {
488  se->registerSubsystem(new QAG4SimulationCalorimeter("CEMC",QAG4SimulationCalorimeter::kProcessG4Hit));
489  }
490  if (ihcal_on)
491  {
492  se->registerSubsystem(new QAG4SimulationCalorimeter("HCALIN",QAG4SimulationCalorimeter::kProcessG4Hit));
493  }
494  if (ohcal_on)
495  {
496  se->registerSubsystem(new QAG4SimulationCalorimeter("HCALOUT",QAG4SimulationCalorimeter::kProcessG4Hit));
497  }
498  //----------------------
499  // G4HitNtuple
500  //----------------------
501  if (hit_ntuple)
502  {
503  G4HitNtuple *hit = new G4HitNtuple("G4HitNtuple","g4hitntuple.root");
504  hit->AddNode("HCALIN", 0);
505  hit->AddNode("HCALOUT", 1);
506  hit->AddNode("CRYO", 2);
507  hit->AddNode("BlackHole", 3);
508  hit->AddNode("ABSORBER_HCALIN", 10);
509  hit->AddNode("ABSORBER_HCALOUT", 11);
510  se->registerSubsystem(hit);
511  }
512  //----------------------
513  // save a comprehensive evaluation file
514  //----------------------
515  if (dstreader)
516  {
517  PHG4DSTReader* ana = new PHG4DSTReader(
518  string("DSTReader.root"));
519  ana->set_save_particle(true);
520  ana->set_load_all_particle(false);
521  ana->set_load_active_particle(false);
522  ana->set_save_vertex(true);
523  ana->set_tower_zero_sup(-1000); // no zero suppression
524 
525  // ana->AddNode("CEMC");
526  // if (absorberactive)
527  // {
528  // ana->AddNode("ABSORBER_CEMC");
529  // }
530  ana->AddTower("SIM_CEMC");
531  ana->AddTower("RAW_LG_CEMC");
532  ana->AddTower("CALIB_LG_CEMC");// Low gain CEMC
533  ana->AddTower("RAW_HG_CEMC");
534  ana->AddTower("CALIB_HG_CEMC");// High gain CEMC
535 
536  ana->AddTower("SIM_HCALOUT");
537  ana->AddTower("SIM_HCALIN");
538 
539  ana->AddTower("RAW_LG_HCALIN");
540  ana->AddTower("RAW_HG_HCALIN");
541  ana->AddTower("RAW_LG_HCALOUT");
542  ana->AddTower("RAW_HG_HCALOUT");
543 
544  ana->AddTower("CALIB_LG_HCALIN");
545  ana->AddTower("CALIB_HG_HCALIN");
546  ana->AddTower("CALIB_LG_HCALOUT");
547  ana->AddTower("CALIB_HG_HCALOUT");
548 
549  ana->AddNode("BlackHole");// add a G4Hit node
550 
551  se->registerSubsystem(ana);
552  }
553  if (dstoutput)
554  {
555  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT","G4Prototype2New.root");
556  se->registerOutputManager(out);
557  }
558 
559  Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
560  se->registerInputManager( in );
561  if (nEvents <= 0)
562  {
563  return 0;
564  }
565  se->run(nEvents);
566 
567  se->End();
568 
569  QAHistManagerDef::saveQARootFile("G4Prototype2_qa.root");
570 
571 
572  // std::cout << "All done" << std::endl;
573  delete se;
574  gSystem->Exit(0);
575  return 0;
576 
577 }
578