Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:20:34

0001 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0002 #include <caloreco/RawTowerCalibration.h>
0003 
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllDstOutputManager.h>
0006 #include <fun4all/Fun4AllDummyInputManager.h>
0007 #include <fun4all/Fun4AllInputManager.h>
0008 #include <fun4all/Fun4AllOutputManager.h>
0009 #include <fun4all/Fun4AllServer.h>
0010 
0011 #include <g4calo/RawTowerBuilder.h>
0012 #include <g4calo/RawTowerDigitizer.h>
0013 
0014 #include <g4caloprototype/Prototype2RawTowerBuilder.h>
0015 #include <g4detectors/PHG4BlockSubsystem.h>
0016 #include <g4detectors/PHG4FullProjSpacalCellReco.h>
0017 #include <g4caloprototype/PHG4Prototype2HcalCellReco.h>
0018 #include <g4caloprototype/PHG4Prototype2InnerHcalSubsystem.h>
0019 #include <g4caloprototype/PHG4Prototype2OuterHcalSubsystem.h>
0020 #include <g4caloprototype/PHG4SpacalPrototypeSubsystem.h>
0021 
0022 #include <g4eval/PHG4DSTReader.h>
0023 
0024 #include <g4histos/G4HitNtuple.h>
0025 
0026 #include <g4main/PHG4ParticleGun.h>
0027 #include <g4main/PHG4Reco.h>
0028 #include <g4main/PHG4SimpleEventGenerator.h>
0029 #include <g4main/PHG4TruthSubsystem.h>
0030 
0031 #include <phool/recoConsts.h>
0032 
0033 #include <qa_modules/QAHistManagerDef.h>
0034 #include <qa_modules/QAG4SimulationCalorimeter.h>
0035 
0036 R__LOAD_LIBRARY(libcalo_reco.so)
0037 R__LOAD_LIBRARY(libfun4all.so)
0038 R__LOAD_LIBRARY(libg4calo.so)
0039 R__LOAD_LIBRARY(libg4caloprototype.so)
0040 R__LOAD_LIBRARY(libg4eval.so)
0041 R__LOAD_LIBRARY(libg4histos.so)
0042 R__LOAD_LIBRARY(libg4testbench.so)
0043 R__LOAD_LIBRARY(libqa_modules.so)
0044 
0045 #endif
0046 
0047 int Fun4All_G4_Prototype2(int nEvents = 1)
0048 {
0049 
0050   gSystem->Load("libfun4all");
0051   gSystem->Load("libg4caloprototype");
0052   gSystem->Load("libg4testbench");
0053   gSystem->Load("libg4histos");
0054   gSystem->Load("libg4eval");
0055   gSystem->Load("libqa_modules");
0056   gSystem->Load("libg4calo");
0057   gSystem->Load("libcalo_reco");
0058 
0059   bool cemc_on = true;
0060   bool cemc_cell = cemc_on && true;
0061   bool cemc_twr = cemc_cell && true;
0062   bool cemc_digi = cemc_twr && true;
0063   bool cemc_twrcal = cemc_digi && true;
0064 
0065   bool ihcal_on = true;
0066   bool ihcal_cell = ihcal_on && true;
0067   bool ihcal_twr = ihcal_cell && true;
0068   bool ihcal_digi = ihcal_twr && true;
0069   bool ihcal_twrcal = ihcal_digi && true;
0070 
0071   bool ohcal_on = true;
0072   bool ohcal_cell = ohcal_on && true;
0073   bool ohcal_twr = ohcal_cell && true;
0074   bool ohcal_digi = ohcal_twr && true;
0075   bool ohcal_twrcal =  ohcal_digi && true;
0076 
0077   bool cryo_on = true;
0078   bool bh_on = true;
0079   bool hit_ntuple = false;
0080   bool dstreader = true;
0081   bool dstoutput = false;
0082 
0083   ///////////////////////////////////////////
0084   // Make the Server
0085   //////////////////////////////////////////
0086   Fun4AllServer *se = Fun4AllServer::instance();
0087   //  se->Verbosity(1);
0088   recoConsts *rc = recoConsts::instance();
0089   // only set this if you want a fixed random seed to make
0090   // results reproducible for testing
0091   //   rc->set_IntFlag("RANDOMSEED",12345);
0092 
0093   // Test beam generator
0094   PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0095   gen->add_particles("pi-", 1); // mu-,e-,anti_proton,pi-
0096   gen->set_vertex_distribution_mean(0.0, 0.0, 0);
0097   gen->set_vertex_distribution_width(0.0, .7, .7); // Rough beam profile size @ 16 GeV measured by Abhisek
0098   gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0099                     PHG4SimpleEventGenerator::Gaus, PHG4SimpleEventGenerator::Gaus); // Gauss beam profile
0100   gen->set_eta_range(-.001, .001); // 1mrad angular divergence
0101   gen->set_phi_range(-.001, .001); // 1mrad angular divergence
0102   const double momentum = 32;
0103   gen->set_p_range(momentum,momentum, momentum*2e-2); // 2% momentum smearing
0104   se->registerSubsystem(gen);
0105 
0106   // Simple single particle generator
0107   PHG4ParticleGun *gun = new PHG4ParticleGun();
0108   //  gun->set_name("anti_proton");
0109   //  gun->set_name("geantino");
0110   gun->set_name("proton");
0111   gun->set_vtx(0, 0, 0);
0112   gun->set_mom(120, 0, 0);
0113   // gun->AddParticle("geantino",1.7776,-0.4335,0.);
0114   // gun->AddParticle("geantino",1.7709,-0.4598,0.);
0115   // gun->AddParticle("geantino",2.5621,0.60964,0.);
0116   // gun->AddParticle("geantino",1.8121,0.253,0.);
0117   //  se->registerSubsystem(gun);
0118 
0119   PHG4Reco* g4Reco = new PHG4Reco();
0120   g4Reco->set_field(0);
0121   //  g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
0122 
0123   //----------------------------------------
0124   // EMCal G4
0125   //----------------------------------------
0126   if (cemc_on)
0127     {
0128       PHG4SpacalPrototypeSubsystem *cemc = new PHG4SpacalPrototypeSubsystem("CEMC");
0129       cemc->SetActive();
0130       cemc->SuperDetector("CEMC");
0131       cemc->SetAbsorberActive();
0132       cemc->OverlapCheck(true);
0133       cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0134       cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/Prototype2/Geometry/") ); 
0135       //  cemc->set_double_param("z_rotation_degree", 15); // rotation around CG
0136       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
0137       cemc->set_double_param("ypos", 4); // put it some where in UIUC blocks
0138       cemc->set_double_param("zpos", 4); // put it some where in UIUC blocks
0139       g4Reco->registerSubsystem(cemc);
0140     }
0141   //----------------------------------------
0142   // HCal G4
0143   //----------------------------------------
0144   if (ihcal_on)
0145     {
0146       PHG4Prototype2InnerHcalSubsystem *innerhcal = new PHG4Prototype2InnerHcalSubsystem("HCalIn");
0147       innerhcal->SetActive();
0148       innerhcal->SetAbsorberActive();
0149       innerhcal->SetAbsorberTruth(1);
0150       innerhcal->OverlapCheck(true);
0151       innerhcal->SuperDetector("HCALIN");
0152       g4Reco->registerSubsystem(innerhcal);
0153     }
0154   if (ohcal_on)
0155     {
0156       PHG4Prototype2OuterHcalSubsystem *outerhcal = new PHG4Prototype2OuterHcalSubsystem("HCalOut");
0157       outerhcal->SetActive();
0158       outerhcal->SetAbsorberActive();
0159       outerhcal->SetAbsorberTruth(1);
0160       outerhcal->OverlapCheck(true);
0161       outerhcal->SuperDetector("HCALOUT");
0162       g4Reco->registerSubsystem(outerhcal);
0163     }
0164   // Cryostat from engineering drawing
0165   if (cryo_on)
0166     {
0167       PHG4BlockSubsystem *cryo1 = new PHG4BlockSubsystem("cryo1",1);
0168       cryo1->set_double_param("size_x",0.95);
0169       cryo1->set_double_param("size_y",60.96);
0170       cryo1->set_double_param("size_z",60.96);
0171       cryo1->set_double_param("place_x",141.96+0.95/2.);
0172       cryo1->set_string_param("material","G4_Al");
0173       cryo1->SetActive(); // it is an active volume - save G4Hits
0174       cryo1->SuperDetector("CRYO");
0175       g4Reco->registerSubsystem(cryo1);
0176 
0177       PHG4BlockSubsystem *cryo2 = new PHG4BlockSubsystem("cryo2",2);
0178       cryo2->set_double_param("size_x",8.89);
0179       cryo2->set_double_param("size_y",60.96);
0180       cryo2->set_double_param("size_z",60.96);
0181       cryo2->set_double_param("place_x",150.72+8.89/2.);
0182       cryo2->set_string_param("material","G4_Al");
0183       cryo2->SetActive(); // it is an active volume - save G4Hits
0184       cryo2->SuperDetector("CRYO");
0185       g4Reco->registerSubsystem(cryo2);
0186 
0187       PHG4BlockSubsystem *cryo3 = new PHG4BlockSubsystem("cryo3",3);
0188       cryo3->set_double_param("size_x",2.54);
0189       cryo3->set_double_param("size_y",60.96);
0190       cryo3->set_double_param("size_z",60.96);
0191       cryo3->set_double_param("place_x",173.93+2.54/2.);
0192       cryo3->set_string_param("material","G4_Al");
0193       cryo3->SetActive(); // it is an active volume - save G4Hits
0194       cryo3->SuperDetector("CRYO");
0195       g4Reco->registerSubsystem(cryo3);
0196     }
0197   // BLACKHOLE, box surrounding the prototype to check for leakage
0198   if (bh_on)
0199     {
0200       PHG4BlockSubsystem *bh[5];
0201       // surrounding outer hcal
0202       // top
0203       bh[0] = new PHG4BlockSubsystem("bh1",1);
0204       bh[0]->set_double_param("size_x",270.);
0205       bh[0]->set_double_param("size_y",0.01);
0206       bh[0]->set_double_param("size_z",165.);
0207       bh[0]->set_double_param("place_x",270./2.);
0208       bh[0]->set_double_param("place_y",125./2.);
0209       // bottom
0210       bh[1] = new PHG4BlockSubsystem("bh2",2);
0211       bh[1]->set_double_param("size_x",270.);
0212       bh[1]->set_double_param("size_y",0.01);
0213       bh[1]->set_double_param("size_z",165.);
0214       bh[1]->set_double_param("place_x",270./2.);
0215       bh[1]->set_double_param("place_y",-125./2.);
0216       // right side
0217       bh[2] = new PHG4BlockSubsystem("bh3",3);
0218       bh[2]->set_double_param("size_x",270.);
0219       bh[2]->set_double_param("size_y",125.);
0220       bh[2]->set_double_param("size_z",0.01);
0221       bh[2]->set_double_param("place_x",270./2.);
0222       bh[2]->set_double_param("place_z",165./2.);
0223       // left side
0224       bh[3] = new PHG4BlockSubsystem("bh4",4);
0225       bh[3]->set_double_param("size_x",270.);
0226       bh[3]->set_double_param("size_y",125.);
0227       bh[3]->set_double_param("size_z",0.01);
0228       bh[3]->set_double_param("place_x",270./2.);
0229       bh[3]->set_double_param("place_z",-165./2.);
0230       // back
0231       bh[4] = new PHG4BlockSubsystem("bh5",5);
0232       bh[4]->set_double_param("size_x",0.01);
0233       bh[4]->set_double_param("size_y",125.);
0234       bh[4]->set_double_param("size_z",165.);
0235       bh[4]->set_double_param("place_x",270.);
0236       for (int i=0; i<5; i++)
0237     {
0238       bh[i]->BlackHole();
0239       bh[i]->SetActive();
0240       bh[i]->SuperDetector("BlackHole");
0241       bh[i]->OverlapCheck(true);
0242       g4Reco->registerSubsystem(bh[i]);
0243     }
0244     }
0245   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0246   g4Reco->registerSubsystem(truth);
0247 
0248   se->registerSubsystem( g4Reco );
0249   //----------------------------------------
0250   // EMCal digitization
0251   //----------------------------------------
0252 
0253   if (cemc_cell)
0254     {
0255       PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0256       cemc_cells->Detector("CEMC");
0257       cemc_cells->set_timing_window(0.,60.);
0258       cemc_cells->get_light_collection_model().load_data_file(
0259                                   string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype2Module.xml"),
0260                                   "data_grid_light_guide_efficiency","data_grid_fiber_trans");
0261 
0262       se->registerSubsystem(cemc_cells);
0263     }
0264   if (cemc_twr)
0265     {
0266       RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0267       TowerBuilder->Detector("CEMC");
0268       TowerBuilder->set_sim_tower_node_prefix("SIM");
0269       se->registerSubsystem(TowerBuilder);
0270     }
0271 
0272 // these parameters are used in multiple modules
0273   const double ADC_per_photoelectron_LG = 0.24; // From Sean Stoll, Mar 29
0274   const double ADC_per_photoelectron_HG = 3.8; // From Sean Stoll, Mar 29
0275   const double photoelectron_per_GeV = 500; //500 photon per total GeV deposition
0276 
0277   if (cemc_digi)
0278     {
0279       const double sampling_fraction = 0.0233369; //  +/-   8.22211e-05  from 15 Degree indenting 8 GeV electron showers
0280 
0281       // low gains
0282       RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerLG");
0283       TowerDigitizer->Detector("CEMC");
0284       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0285       TowerDigitizer->set_digi_algorithm(
0286                      RawTowerDigitizer::kSimple_photon_digitalization);
0287       TowerDigitizer->set_pedstal_central_ADC(0);
0288       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0289       TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_LG);
0290       TowerDigitizer->set_photonelec_yield_visible_GeV(
0291                                photoelectron_per_GeV / sampling_fraction);
0292       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0293       se->registerSubsystem(TowerDigitizer);
0294       // high gains
0295       TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerHG");
0296       TowerDigitizer->Detector("CEMC");
0297       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0298       TowerDigitizer->set_digi_algorithm(
0299                      RawTowerDigitizer::kSimple_photon_digitalization);
0300       TowerDigitizer->set_pedstal_central_ADC(0);
0301       TowerDigitizer->set_pedstal_width_ADC(15); // From John Haggerty, Mar 29
0302       TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_HG);
0303       TowerDigitizer->set_photonelec_yield_visible_GeV(
0304                                photoelectron_per_GeV / sampling_fraction);
0305       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0306       se->registerSubsystem(TowerDigitizer);
0307     }
0308   if (cemc_twrcal)
0309     {
0310       RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationLG");
0311       TowerCalibration->Detector("CEMC");
0312       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0313       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0314       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0315       TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_LG / photoelectron_per_GeV);
0316       TowerCalibration->set_pedstal_ADC(0);
0317       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0318       se->registerSubsystem(TowerCalibration);
0319 
0320       TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationHG");
0321       TowerCalibration->Detector("CEMC");
0322       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0323       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0324       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0325       TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_HG / photoelectron_per_GeV);
0326       TowerCalibration->set_pedstal_ADC(0);
0327       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0328       se->registerSubsystem(TowerCalibration);
0329     }
0330 
0331   //----------------------------------------
0332   // HCal towering
0333   //----------------------------------------
0334   if (ihcal_cell)
0335     {
0336       PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALinCellReco");
0337       hccell->Detector("HCALIN");
0338       se->registerSubsystem(hccell);
0339     }
0340 
0341   if (ohcal_cell)
0342     {
0343       PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALoutCellReco");
0344       hccell->Detector("HCALOUT");
0345       se->registerSubsystem(hccell);
0346     }
0347   if (ihcal_twr)
0348     {
0349       Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALinRawTowerBuilder");
0350       hcaltwr->Detector("HCALIN");
0351       hcaltwr->set_sim_tower_node_prefix("SIM");
0352       se->registerSubsystem(hcaltwr);
0353     }
0354   if (ohcal_twr)
0355     {
0356       Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALoutRawTowerBuilder");
0357       hcaltwr->Detector("HCALOUT");
0358       hcaltwr->set_sim_tower_node_prefix("SIM");
0359       se->registerSubsystem(hcaltwr);
0360     }
0361 
0362   //----------------------------------------
0363   // HCal digitization
0364   //----------------------------------------
0365   //  From: Abhisek Sen [mailto:sen.abhisek@gmail.com]
0366   //  Sent: Tuesday, April 19, 2016 10:55 PM
0367   //  To: Huang, Jin <jhuang@bnl.gov>; Haggerty, John <haggerty@bnl.gov>
0368 
0369   //  HCALIN:
0370   //     1/5 pixel / HG ADC channel
0371   //     32/5 pixel / LG ADC channel
0372   //     0.4 MeV/ LG ADC
0373   //     0.4/32 MeV/ HG ADC
0374 
0375   //  HCALOUT:
0376   //     1/5 pixel / HG ADC channel
0377   //     16/5 pixel / LG ADC channel
0378   //     0.2 MeV/ LG ADC
0379   //     0.2/16 MeV/ HG ADC
0380 
0381   if (ihcal_digi)
0382     {
0383       RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerLG");
0384       TowerDigitizer->Detector("HCALIN");
0385       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0386       TowerDigitizer->set_digi_algorithm(
0387                      RawTowerDigitizer::kSimple_photon_digitalization);
0388       TowerDigitizer->set_pedstal_central_ADC(0);
0389       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0390       TowerDigitizer->set_photonelec_ADC(32. / 5.);
0391       TowerDigitizer->set_photonelec_yield_visible_GeV(32. / 5 / (0.4e-3));
0392       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0393       se->registerSubsystem(TowerDigitizer);
0394 
0395       TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerHG");
0396       TowerDigitizer->Detector("HCALIN");
0397       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0398       TowerDigitizer->set_digi_algorithm(
0399                      RawTowerDigitizer::kSimple_photon_digitalization);
0400       TowerDigitizer->set_pedstal_central_ADC(0);
0401       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0402       TowerDigitizer->set_photonelec_ADC(1. / 5.);
0403       TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.4e-3 / 32));
0404       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0405       se->registerSubsystem(TowerDigitizer);
0406     }
0407   if (ohcal_digi)
0408     {
0409       RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerLG");
0410       TowerDigitizer->Detector("HCALOUT");
0411       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0412       TowerDigitizer->set_digi_algorithm(
0413                      RawTowerDigitizer::kSimple_photon_digitalization);
0414       TowerDigitizer->set_pedstal_central_ADC(0);
0415       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0416       TowerDigitizer->set_photonelec_ADC(16. / 5.);
0417       TowerDigitizer->set_photonelec_yield_visible_GeV(16. / 5 / (0.2e-3));
0418       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0419       se->registerSubsystem(TowerDigitizer);
0420 
0421       TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerHG");
0422       TowerDigitizer->Detector("HCALOUT");
0423       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0424       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0425       TowerDigitizer->set_pedstal_central_ADC(0);
0426       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0427       TowerDigitizer->set_photonelec_ADC(1. / 5.);
0428       TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.2e-3 / 16));
0429       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0430       se->registerSubsystem(TowerDigitizer);
0431     }
0432   //----------------------------------------
0433   // HCal calibration
0434   //----------------------------------------
0435   // 32 GeV Pi+ scan
0436   const double visible_sample_fraction_HCALIN = 7.19505e-02 ; // 1.34152e-02
0437   const double visible_sample_fraction_HCALOUT = 0.0313466 ; //  +/-   0.0067744
0438 
0439   if (ihcal_twrcal)
0440     {
0441       RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationLG");
0442       TowerCalibration->Detector("HCALIN");
0443       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0444       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0445       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0446       TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / visible_sample_fraction_HCALIN);
0447       TowerCalibration->set_pedstal_ADC(0);
0448       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0449       se->registerSubsystem(TowerCalibration);
0450 
0451       TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationHG");
0452       TowerCalibration->Detector("HCALIN");
0453       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0454       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0455       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0456       TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / 32 / visible_sample_fraction_HCALIN);
0457       TowerCalibration->set_pedstal_ADC(0);
0458       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0459       se->registerSubsystem(TowerCalibration);
0460     }
0461   if (ohcal_twrcal)
0462     {
0463       RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationLG");
0464       TowerCalibration->Detector("HCALOUT");
0465       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0466       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0467       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0468       TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / visible_sample_fraction_HCALOUT);
0469       TowerCalibration->set_pedstal_ADC(0);
0470       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0471       se->registerSubsystem(TowerCalibration);
0472 
0473       TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationHG");
0474       TowerCalibration->Detector("HCALOUT");
0475       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0476       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0477       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0478       TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / 16 / visible_sample_fraction_HCALOUT);
0479       TowerCalibration->set_pedstal_ADC(0);
0480       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0481       se->registerSubsystem(TowerCalibration);
0482     }
0483   //----------------------
0484   // QA Histograms
0485   //----------------------
0486   if (cemc_on)
0487     {
0488       se->registerSubsystem(new QAG4SimulationCalorimeter("CEMC",QAG4SimulationCalorimeter::kProcessG4Hit));
0489     }
0490   if (ihcal_on)
0491     {
0492       se->registerSubsystem(new QAG4SimulationCalorimeter("HCALIN",QAG4SimulationCalorimeter::kProcessG4Hit));
0493     }
0494   if (ohcal_on)
0495     {
0496       se->registerSubsystem(new QAG4SimulationCalorimeter("HCALOUT",QAG4SimulationCalorimeter::kProcessG4Hit));
0497     }
0498   //----------------------
0499   // G4HitNtuple
0500   //----------------------
0501   if (hit_ntuple)
0502     {
0503       G4HitNtuple *hit = new G4HitNtuple("G4HitNtuple","g4hitntuple.root");
0504       hit->AddNode("HCALIN", 0);
0505       hit->AddNode("HCALOUT", 1);
0506       hit->AddNode("CRYO", 2);
0507       hit->AddNode("BlackHole", 3);
0508       hit->AddNode("ABSORBER_HCALIN", 10);
0509       hit->AddNode("ABSORBER_HCALOUT", 11);
0510       se->registerSubsystem(hit);
0511     }
0512   //----------------------
0513   // save a comprehensive  evaluation file
0514   //----------------------
0515   if (dstreader)
0516     {
0517       PHG4DSTReader* ana = new PHG4DSTReader(
0518                          string("DSTReader.root"));
0519       ana->set_save_particle(true);
0520       ana->set_load_all_particle(false);
0521       ana->set_load_active_particle(false);
0522       ana->set_save_vertex(true);
0523       ana->set_tower_zero_sup(-1000); // no zero suppression
0524 
0525       //  ana->AddNode("CEMC");
0526       //  if (absorberactive)
0527       //    {
0528       //      ana->AddNode("ABSORBER_CEMC");
0529       //    }
0530       ana->AddTower("SIM_CEMC");
0531       ana->AddTower("RAW_LG_CEMC");
0532       ana->AddTower("CALIB_LG_CEMC");// Low gain CEMC
0533       ana->AddTower("RAW_HG_CEMC");
0534       ana->AddTower("CALIB_HG_CEMC");// High gain CEMC
0535 
0536       ana->AddTower("SIM_HCALOUT");
0537       ana->AddTower("SIM_HCALIN");
0538 
0539       ana->AddTower("RAW_LG_HCALIN");
0540       ana->AddTower("RAW_HG_HCALIN");
0541       ana->AddTower("RAW_LG_HCALOUT");
0542       ana->AddTower("RAW_HG_HCALOUT");
0543 
0544       ana->AddTower("CALIB_LG_HCALIN");
0545       ana->AddTower("CALIB_HG_HCALIN");
0546       ana->AddTower("CALIB_LG_HCALOUT");
0547       ana->AddTower("CALIB_HG_HCALOUT");
0548 
0549       ana->AddNode("BlackHole");// add a G4Hit node
0550 
0551       se->registerSubsystem(ana);
0552     }
0553   if (dstoutput)
0554     {
0555       Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT","G4Prototype2New.root");
0556       se->registerOutputManager(out);
0557     }
0558 
0559   Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0560   se->registerInputManager( in );
0561   if (nEvents <= 0)
0562     {
0563       return 0;
0564     }
0565   se->run(nEvents);
0566 
0567   se->End();
0568 
0569   QAHistManagerDef::saveQARootFile("G4Prototype2_qa.root");
0570 
0571 
0572   //   std::cout << "All done" << std::endl;
0573   delete se;
0574   gSystem->Exit(0);
0575   return 0;
0576 
0577 }
0578