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