calo/tower_slope_macro/G4_HcalOut_ref.C

0001 #ifndef MACRO_G4HCALOUTREF_C
0002 #define MACRO_G4HCALOUTREF_C
0003
0004 #include <GlobalVariables.C>
0005 #include <QA.C>
0006
0007 #include <g4calo/HcalRawTowerBuilder.h>
0008 #include <g4calo/RawTowerDigitizer.h>
0009
0010 #include <g4detectors/PHG4HcalCellReco.h>
0011 #include <g4detectors/PHG4OuterHcalSubsystem.h>
0012
0013 #include <g4eval/CaloEvaluator.h>
0014
0015 #include <g4main/PHG4Reco.h>
0016
0017 #include <caloreco/RawClusterBuilderGraph.h>
0018 #include <caloreco/RawClusterBuilderTemplate.h>
0019 #include <caloreco/RawTowerCalibration.h>
0020 #include <qa_modules/QAG4SimulationCalorimeter.h>
0021
0022 #include <fun4all/Fun4AllServer.h>
0023
0024 R__LOAD_LIBRARY(libcalo_reco.so)
0025 R__LOAD_LIBRARY(libg4calo.so)
0026 R__LOAD_LIBRARY(libg4detectors.so)
0027 R__LOAD_LIBRARY(libg4eval.so)
0028 R__LOAD_LIBRARY(libqa_modules.so)
0029
0030 namespace Enable
0031 {
0032   bool HCALOUT = false;
0033   bool HCALOUT_ABSORBER = false;
0034   bool HCALOUT_OVERLAPCHECK = false;
0035   bool HCALOUT_CELL = false;
0036   bool HCALOUT_TOWER = false;
0037   bool HCALOUT_CLUSTER = false;
0038   bool HCALOUT_EVAL = false;
0039   bool HCALOUT_QA = false;
0040   int HCALOUT_VERBOSITY = 0;
0041 }  // namespace Enable
0042
0043 namespace G4HCALOUT
0044 {
0045   double outer_radius = 264.71;
0046   double size_z = 304.91 * 2;
0047
0048   // Digitization (default photon digi):
0049   RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0050   // directly pass the energy of sim tower to digitized tower
0051   // kNo_digitization
0052   // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
0053   // kSimple_photon_digitization
0054   // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
0055   // kSiPM_photon_digitization
0056
0057   enum enu_HCalOut_clusterizer
0058   {
0059     kHCalOutGraphClusterizer,
0060     kHCalOutTemplateClusterizer
0061   };
0062
0063   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0064   enu_HCalOut_clusterizer HCalOut_clusterizer = kHCalOutTemplateClusterizer;
0065   //! graph clusterizer, RawClusterBuilderGraph
0066   //enu_HCalOut_clusterizer HCalOut_clusterizer = kHCalOutGraphClusterizer;
0067 }  // namespace G4HCALOUT
0068
0069 // Init is called by G4Setup.C
0070 void HCalOuterInit()
0071 {
0072   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4HCALOUT::outer_radius);
0073   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4HCALOUT::size_z / 2.);
0074   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4HCALOUT::size_z / 2.);
0075 }
0076
0077 double HCalOuter(PHG4Reco *g4Reco,
0078                  double radius,
0079                  const int crossings)
0080 {
0081   bool AbsorberActive = Enable::ABSORBER || Enable::HCALOUT_ABSORBER;
0082   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::HCALOUT_OVERLAPCHECK;
0083   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0084
0085   PHG4OuterHcalSubsystem *hcal = new PHG4OuterHcalSubsystem("HCALOUT");
0086   // hcal->set_double_param("inner_radius", 183.3);
0087   //-----------------------------------------
0088   // the light correction can be set in a single call
0089   // hcal->set_double_param("light_balance_inner_corr", NAN);
0090   // hcal->set_double_param("light_balance_inner_radius", NAN);
0091   // hcal->set_double_param("light_balance_outer_corr", NAN);
0092   // hcal->set_double_param("light_balance_outer_radius", NAN);
0093   // hcal->set_double_param("magnet_cutout_radius", 195.31);
0094   // hcal->set_double_param("magnet_cutout_scinti_radius", 195.96);
0095   // hcal->SetLightCorrection(NAN,NAN,NAN,NAN);
0096   //-----------------------------------------
0097   // hcal->set_double_param("outer_radius", G4HCALOUT::outer_radius);
0098   // hcal->set_double_param("place_x", 0.);
0099   // hcal->set_double_param("place_y", 0.);
0100   // hcal->set_double_param("place_z", 0.);
0101   // hcal->set_double_param("rot_x", 0.);
0102   // hcal->set_double_param("rot_y", 0.);
0103   // hcal->set_double_param("rot_z", 0.);
0104   // hcal->set_double_param("scinti_eta_coverage", 1.1);
0105   // hcal->set_double_param("scinti_gap", 0.85);
0106   // hcal->set_double_param("scinti_gap_neighbor", 0.1);
0107   // hcal->set_double_param("scinti_inner_radius",183.89);
0108   // hcal->set_double_param("scinti_outer_radius",263.27);
0109   // hcal->set_double_param("scinti_tile_thickness", 0.7);
0110   // hcal->set_double_param("size_z", G4HCALOUT::size_z);
0111   // hcal->set_double_param("steplimits", NAN);
0112   // hcal->set_double_param("tilt_angle", -11.23);
0113
0114   // hcal->set_int_param("light_scint_model", 1);
0115   // hcal->set_int_param("magnet_cutout_first_scinti", 8);
0116   // hcal->set_int_param("ncross", 0);
0117   // hcal->set_int_param("n_towers", 64);
0118   // hcal->set_int_param("n_scinti_plates_per_tower", 5);
0119   // hcal->set_int_param("n_scinti_tiles", 12);
0120
0121   // hcal->set_string_param("material", "Steel_1006");
0122
0123   hcal->SetActive();
0124   hcal->SuperDetector("HCALOUT");
0125   if (AbsorberActive)
0126   {
0127     hcal->SetAbsorberActive();
0128   }
0129   hcal->OverlapCheck(OverlapCheck);
0130   g4Reco->registerSubsystem(hcal);
0131
0132   radius = hcal->get_double_param("outer_radius");
0133
0134   radius += no_overlapp;
0135
0136   return radius;
0137 }
0138
0139 void HCALOuter_Cells()
0140 {
0141   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0142
0143   Fun4AllServer *se = Fun4AllServer::instance();
0144
0145   PHG4HcalCellReco *hc = new PHG4HcalCellReco("HCALOUT_CELLRECO");
0146   hc->Detector("HCALOUT");
0147   //  hc->Verbosity(2);
0148   // check for energy conservation - needs modified "infinite" timing cuts
0149   // 0-999999999
0150   //  hc->checkenergy();
0151   // timing cuts with their default settings
0152   // hc->set_double_param("tmin",0.);
0153   // hc->set_double_param("tmax",60.0);
0154   // or all at once:
0155   // hc->set_timing_window(0.0,60.0);
0156   // this sets all cells to a fixed energy for debugging
0157   // hc->set_fixed_energy(1.);
0158   se->registerSubsystem(hc);
0159
0160   return;
0161 }
0162
0163 void HCALOuter_Towers()
0164 {
0165   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0166
0167   Fun4AllServer *se = Fun4AllServer::instance();
0168
0169   HcalRawTowerBuilder *TowerBuilder = new HcalRawTowerBuilder("HcalOutRawTowerBuilder");
0170   TowerBuilder->Detector("HCALOUT");
0171   TowerBuilder->set_sim_tower_node_prefix("SIM");
0172   // this sets specific decalibration factors
0173   // for a given cell
0174   // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0175   // for a whole tower
0176   // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0177   // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0178   // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0179   TowerBuilder->Verbosity(verbosity);
0180   se->registerSubsystem(TowerBuilder);
0181
0182   // From 2016 Test beam sim
0183   RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HcalOutRawTowerDigitizer");
0184   TowerDigitizer->Detector("HCALOUT");
0185   //  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0186   TowerDigitizer->set_digi_algorithm(G4HCALOUT::TowerDigi);
0187   TowerDigitizer->set_pedstal_central_ADC(0);
0188   TowerDigitizer->set_pedstal_width_ADC(1);  // From Jin's guess. No EMCal High Gain data yet! TODO: update
0189   TowerDigitizer->set_photonelec_ADC(16. / 5.);
0190   TowerDigitizer->set_photonelec_yield_visible_GeV(16. / 5 / (0.2e-3));
0191   TowerDigitizer->set_zero_suppression_ADC(-0);  // no-zero suppression
0192   se->registerSubsystem(TowerDigitizer);
0193
0194   const double visible_sample_fraction_HCALOUT = 3.38021e-02;  // /gpfs/mnt/gpfs04/sphenix/user/jinhuang/prod_analysis/hadron_shower_res_nightly/./G4Hits_sPHENIX_pi-_eta0_16GeV.root_qa.rootQA_Draw_HCALOUT_G4Hit.pdf
0195
0196   RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HcalOutRawTowerCalibration");
0197   TowerCalibration->Detector("HCALOUT");
0198   //  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0199   //  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0200
0201   //TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0202
0203
0204   // see documentation in github/jtaou/coresoft
0205   // ...[calib/towerslope] /macros/G4_CEmc_Spacal.C
0206   // for more info about TowerCalibration dbfile mode
0207   TowerCalibration->set_calib_algorithm(RawTowerCalibration::kDbfile_tbt_gain_corr);
0208   TowerCalibration->set_UseConditionsDB(false);
0209   //
0210   //TowerCalibration->set_CalibrationFileName("HCALOUT_GainsCalib1.22.txt");
0211   TowerCalibration->set_CalibrationFileName("HCALOUT_GainsCalib1.12_PosEtaOnly_NegEta1.0.txt");
0212
0213
0214   if (G4HCALOUT::TowerDigi == RawTowerDigitizer::kNo_digitization)
0215   {
0216     // 0.033 extracted from electron sims (edep(scintillator)/edep(total))
0217     TowerCalibration->set_calib_const_GeV_ADC(1. / 0.033);
0218   }
0219   else
0220   {
0221     TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / visible_sample_fraction_HCALOUT);
0222   }
0223   TowerCalibration->set_pedstal_ADC(0);
0224   se->registerSubsystem(TowerCalibration);
0225
0226   return;
0227 }
0228
0229 void HCALOuter_Clusters()
0230 {
0231   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0232
0233   Fun4AllServer *se = Fun4AllServer::instance();
0234
0235   if (G4HCALOUT::HCalOut_clusterizer == G4HCALOUT::kHCalOutTemplateClusterizer)
0236   {
0237     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("HcalOutRawClusterBuilderTemplate");
0238     ClusterBuilder->Detector("HCALOUT");
0239     ClusterBuilder->SetCylindricalGeometry();  // has to be called after Detector()
0240     ClusterBuilder->Verbosity(verbosity);
0241     se->registerSubsystem(ClusterBuilder);
0242   }
0243   else if (G4HCALOUT::HCalOut_clusterizer == G4HCALOUT::kHCalOutGraphClusterizer)
0244   {
0245     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("HcalOutRawClusterBuilderGraph");
0246     ClusterBuilder->Detector("HCALOUT");
0247     ClusterBuilder->Verbosity(verbosity);
0248     se->registerSubsystem(ClusterBuilder);
0249   }
0250   else
0251   {
0252     cout << "HCALOuter_Clusters - unknown clusterizer setting!" << endl;
0253     exit(1);
0254   }
0255
0256   return;
0257 }
0258
0259 void HCALOuter_Eval(const std::string &outputfile)
0260 {
0261   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0262
0263   Fun4AllServer *se = Fun4AllServer::instance();
0264
0265   CaloEvaluator *eval = new CaloEvaluator("HCALOUTEVALUATOR", "HCALOUT", outputfile);
0266   eval->Verbosity(verbosity);
0267   se->registerSubsystem(eval);
0268
0269   return;
0270 }
0271
0272 void HCALOuter_QA()
0273 {
0274   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::HCALOUT_VERBOSITY);
0275
0276   Fun4AllServer *se = Fun4AllServer::instance();
0277   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("HCALOUT");
0278   qa->Verbosity(verbosity);
0279   se->registerSubsystem(qa);
0280
0281   return;
0282 }
0283
0284 #endif