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