macros/common/G4_HcalIn_ref.C

0001 // Inner HCal reconstruction macro
0002 #ifndef MACRO_G4HCALINREF_C
0003 #define MACRO_G4HCALINREF_C
0004
0005 #include <GlobalVariables.C>
0006 #include <QA.C>
0007
0008 #include <g4calo/HcalRawTowerBuilder.h>
0009 #include <g4calo/RawTowerDigitizer.h>
0010
0011 #include <g4ihcal/PHG4IHCalSubsystem.h>
0012
0013 #include <g4detectors/PHG4CylinderSubsystem.h>
0014 #include <g4detectors/PHG4HcalCellReco.h>
0015 #include <g4detectors/PHG4InnerHcalSubsystem.h>
0016
0017 #include <g4eval/CaloEvaluator.h>
0018
0019 #include <g4main/PHG4Reco.h>
0020
0021 #include <calowaveformsim/CaloWaveformSim.h>
0022
0023 #include <calobase/TowerInfoDefs.h>
0024
0025 #include <caloreco/CaloTowerBuilder.h>
0026 #include <caloreco/CaloTowerCalib.h>
0027 #include <caloreco/CaloTowerStatus.h>
0028 #include <caloreco/CaloWaveformProcessing.h>
0029 #include <caloreco/RawClusterBuilderGraph.h>
0030 #include <caloreco/RawClusterBuilderTemplate.h>
0031 #include <caloreco/RawTowerCalibration.h>
0032
0033 #include <simqa_modules/QAG4SimulationCalorimeter.h>
0034
0035 #include <fun4all/Fun4AllServer.h>
0036
0037 R__LOAD_LIBRARY(libcalo_reco.so)
0038 R__LOAD_LIBRARY(libg4calo.so)
0039 R__LOAD_LIBRARY(libCaloWaveformSim.so)
0040 R__LOAD_LIBRARY(libg4detectors.so)
0041 R__LOAD_LIBRARY(libg4eval.so)
0042 R__LOAD_LIBRARY(libg4ihcal.so)
0043 R__LOAD_LIBRARY(libsimqa_modules.so)
0044
0045 namespace Enable
0046 {
0047   bool HCALIN = false;
0048   bool HCALIN_ABSORBER = false;
0049   bool HCALIN_OVERLAPCHECK = false;
0050   bool HCALIN_CELL = false;
0051   bool HCALIN_TOWER = false;
0052   bool HCALIN_CLUSTER = false;
0053   bool HCALIN_EVAL = false;
0054   bool HCALIN_QA = false;
0055   bool HCALIN_SUPPORT = false;
0056   bool HCALIN_OLD = false;
0057   bool HCALIN_G4Hit = true;
0058   bool HCALIN_TOWERINFO = false;
0059   int HCALIN_VERBOSITY = 0;
0060 }  // namespace Enable
0061
0062 namespace G4HCALIN
0063 {
0064   double inch = 2.54;
0065   double support_ring_outer_radius = 177.323;
0066   double support_ring_z = 175.375 * inch / 2;
0067   double dz = 4. * inch;
0068
0069   double phistart = NAN;
0070   double tower_emin = NAN;
0071   int light_scint_model = -1;
0072   int tower_energy_source = -1;
0073
0074   // Inner HCal absorber material selector:
0075   // false - old version, absorber material is SS310
0076   // true - default Choose if you want Aluminum
0077   bool inner_hcal_material_Al = true;
0078
0079   int inner_hcal_eic = 0;
0080
0081   // Digitization (default photon digi):
0082   RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0083   // directly pass the energy of sim tower to digitized tower
0084   // kNo_digitization
0085   // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
0086   // kSimple_photon_digitization
0087   // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
0088   // kSiPM_photon_digitization
0089
0090   enum enu_HCalIn_clusterizer
0091   {
0092     kHCalInGraphClusterizer,
0093
0094     kHCalInTemplateClusterizer
0095   };
0096
0097   bool useTowerInfoV2 = true;
0098   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0099   enu_HCalIn_clusterizer HCalIn_clusterizer = kHCalInTemplateClusterizer;
0100   //! graph clusterizer, RawClusterBuilderGraph
0101   // enu_HCalIn_clusterizer HCalIn_clusterizer = kHCalInGraphClusterizer;
0102 }  // namespace G4HCALIN
0103
0104 // Init is called by G4Setup.C
0105 void HCalInnerInit(const int iflag = 0)
0106 {
0107   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4HCALIN::support_ring_outer_radius);
0108   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4HCALIN::support_ring_z + G4HCALIN::dz / 2.);
0109   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4HCALIN::support_ring_z - G4HCALIN::dz / 2.);
0110   if (iflag == 1)
0111   {
0112     G4HCALIN::inner_hcal_eic = 1;
0113   }
0114 }
0115
0116 double HCalInner(PHG4Reco *g4Reco,
0117                  double radius,
0118                  const int crossings)
0119 {
0120   bool AbsorberActive = Enable::ABSORBER || Enable::HCALIN_ABSORBER;
0121   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::HCALIN_OVERLAPCHECK;
0122   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0123
0124   PHG4DetectorSubsystem *hcal = nullptr;
0125   //  Mephi Maps
0126   //  Maps are different for old/new but how to set is identical
0127   //  here are the ones for the gdml based inner hcal
0128   //  use hcal->set_string_param("MapFileName",""); to disable map
0129   //  hcal->set_string_param("MapFileName",std::string(getenv("CALIBRATIONROOT")) + "/HCALIN/tilemap/ihcalgdmlmap09212022.root");
0130   //  hcal->set_string_param("MapHistoName","ihcalcombinedgdmlnormtbyt");
0131   //  use hcal->set_string_param("MapFileName",""); to disable map
0132   //  hcal->set_string_param("MapFileName",std::string(getenv("CALIBRATIONROOT")) + "/HCALIN/tilemap/ihcalgdmlmap09212022.root");
0133   //  hcal->set_string_param("MapHistoName","ihcalcombinedgdmlnormtbyt");
0134
0135   // all sizes are in cm!
0136   if (Enable::HCALIN_OLD)
0137   {
0138     hcal = new PHG4InnerHcalSubsystem("HCALIN");
0139     // these are the parameters you can change with their default settings
0140     // hcal->set_string_param("material","SS310");
0141     if (G4HCALIN::inner_hcal_material_Al)
0142     {
0143       if (verbosity > 0)
0144       {
0145         std::cout << "HCalInner - construct inner HCal absorber with G4_Al" << std::endl;
0146       }
0147       hcal->set_string_param("material", "G4_Al");
0148     }
0149     else
0150     {
0151       if (verbosity > 0)
0152       {
0153         std::cout << "HCalInner - construct inner HCal absorber with SS310" << std::endl;
0154       }
0155       hcal->set_string_param("material", "SS310");
0156     }
0157     // hcal->set_double_param("inner_radius", 117.27);
0158     //-----------------------------------------
0159     // the light correction can be set in a single call
0160     // hcal->set_double_param("light_balance_inner_corr", NAN);
0161     // hcal->set_double_param("light_balance_inner_radius", NAN);
0162     // hcal->set_double_param("light_balance_outer_corr", NAN);
0163     // hcal->set_double_param("light_balance_outer_radius", NAN);
0164     // hcal->SetLightCorrection(NAN,NAN,NAN,NAN);
0165     //-----------------------------------------
0166     // hcal->set_double_param("outer_radius", 134.42);
0167     // hcal->set_double_param("place_x", 0.);
0168     // hcal->set_double_param("place_y", 0.);
0169     // hcal->set_double_param("place_z", 0.);
0170     // hcal->set_double_param("rot_x", 0.);
0171     // hcal->set_double_param("rot_y", 0.);
0172     // hcal->set_double_param("rot_z", 0.);
0173     // hcal->set_double_param("scinti_eta_coverage", 1.1);
0174     // hcal->set_double_param("scinti_gap_neighbor", 0.1);
0175     // hcal->set_double_param("scinti_inner_gap", 0.85);
0176     // hcal->set_double_param("scinti_outer_gap", 1.22 * (5.0 / 4.0));
0177     // hcal->set_double_param("scinti_outer_radius", 133.3);
0178     // hcal->set_double_param("scinti_tile_thickness", 0.7);
0179     // hcal->set_double_param("size_z", 175.94 * 2);
0180     // hcal->set_double_param("steplimits", NAN);
0181     // hcal->set_double_param("tilt_angle", 36.15);
0182
0183     // hcal->set_int_param("light_scint_model", 1);
0184     // hcal->set_int_param("ncross", 0);
0185     // hcal->set_int_param("n_towers", 64);
0186     // hcal->set_int_param("n_scinti_plates_per_tower", 4);
0187     // hcal->set_int_param("n_scinti_tiles", 12);
0188
0189     // hcal->set_string_param("material", "SS310");
0190   }
0191   else
0192   {
0193     hcal = new PHG4IHCalSubsystem("HCALIN");
0194     // hcal->set_string_param("GDMPath", "mytestgdml.gdml"); // try other gdml file
0195   }
0196   if (G4HCALIN::light_scint_model >= 0)
0197   {
0198     hcal->set_int_param("light_scint_model", G4HCALIN::light_scint_model);
0199   }
0200   hcal->SetActive();
0201   hcal->SuperDetector("HCALIN");
0202   if (AbsorberActive)
0203   {
0204     hcal->SetAbsorberActive();
0205   }
0206   if (!std::isfinite(G4HCALIN::phistart))
0207   {
0208     if (Enable::HCALIN_OLD)
0209     {
0210       G4HCALIN::phistart = 0.0328877688;  // offet in phi (from zero) extracted from geantinos
0211     }
0212     else
0213     {
0214       G4HCALIN::phistart = 0.0445549893;  // offet in phi (from zero) extracted from geantinos
0215     }
0216   }
0217   hcal->set_int_param("saveg4hit", Enable::HCALIN_G4Hit);
0218   hcal->set_double_param("phistart", G4HCALIN::phistart);
0219   hcal->OverlapCheck(OverlapCheck);
0220
0221   g4Reco->registerSubsystem(hcal);
0222
0223   radius = hcal->get_double_param("outer_radius");
0224
0225   // HCalInner_SupportRing(g4Reco);
0226
0227   radius += no_overlapp;
0228   return radius;
0229 }
0230
0231 void HCALInner_Cells()
0232 {
0233   if (!Enable::HCALIN_G4Hit) return;
0234   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0235
0236   Fun4AllServer *se = Fun4AllServer::instance();
0237
0238   PHG4HcalCellReco *hc = new PHG4HcalCellReco("HCALIN_CELLRECO");
0239   hc->Detector("HCALIN");
0240   hc->Verbosity(verbosity);
0241   // check for energy conservation - needs modified "infinite" timing cuts
0242   // 0-999999999
0243   //  hc->checkenergy();
0244   // timing cuts with their default settings
0245   // hc->set_double_param("tmin",0.);
0246   // hc->set_double_param("tmax",60.0);
0247   // or all at once:
0248   // hc->set_timing_window(0.0,60.0);
0249   // this sets all cells to a fixed energy for debugging
0250   // hc->set_fixed_energy(1.);
0251   se->registerSubsystem(hc);
0252
0253   return;
0254 }
0255
0256 void HCALInner_Towers()
0257 {
0258   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0259   Fun4AllServer *se = Fun4AllServer::instance();
0260
0261   if (!Enable::HCALIN_TOWERINFO)
0262   {
0263     if (Enable::HCALIN_G4Hit)
0264     {
0265       HcalRawTowerBuilder *TowerBuilder = new HcalRawTowerBuilder("HcalInRawTowerBuilder");
0266       TowerBuilder->Detector("HCALIN");
0267       TowerBuilder->set_sim_tower_node_prefix("SIM");
0268       if (!std::isfinite(G4HCALIN::phistart))
0269       {
0270         if (Enable::HCALIN_OLD)
0271         {
0272           G4HCALIN::phistart = 0.0328877688;  // offet in phi (from zero) extracted from geantinos
0273         }
0274         else
0275         {
0276           G4HCALIN::phistart = 0.0445549893;  // offet in phi (from zero) extracted from geantinos
0277         }
0278       }
0279       TowerBuilder->set_double_param("phistart", G4HCALIN::phistart);
0280       if (std::isfinite(G4HCALIN::tower_emin))
0281       {
0282         TowerBuilder->set_double_param("emin", G4HCALIN::tower_emin);
0283       }
0284       if (G4HCALIN::tower_energy_source >= 0)
0285       {
0286         TowerBuilder->set_int_param("tower_energy_source", G4HCALIN::tower_energy_source);
0287       }
0288       // this sets specific decalibration factors
0289       // for a given cell
0290       // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0291       // for a whole tower
0292       // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0293       TowerBuilder->Verbosity(verbosity);
0294       se->registerSubsystem(TowerBuilder);
0295     }
0296     // From 2016 Test beam sim
0297     RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HcalInRawTowerDigitizer");
0298     TowerDigitizer->Detector("HCALIN");
0299     //  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0300     TowerDigitizer->set_digi_algorithm(G4HCALIN::TowerDigi);
0301     TowerDigitizer->set_pedstal_central_ADC(0);
0302     TowerDigitizer->set_pedstal_width_ADC(1);  // From Jin's guess. No EMCal High Gain data yet! TODO: update
0303     TowerDigitizer->set_photonelec_ADC(32. / 5.);
0304     TowerDigitizer->set_photonelec_yield_visible_GeV(32. / 5 / (0.4e-3));
0305     TowerDigitizer->set_zero_suppression_ADC(-0);  // no-zero suppression
0306     TowerDigitizer->Verbosity(verbosity);
0307     if (!Enable::HCALIN_G4Hit) TowerDigitizer->set_towerinfo(RawTowerDigitizer::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0308     se->registerSubsystem(TowerDigitizer);
0309
0310     // Default sampling fraction for SS310
0311     double visible_sample_fraction_HCALIN = 0.0631283;  //, /gpfs/mnt/gpfs04/sphenix/user/jinhuang/prod_analysis/hadron_shower_res_nightly/./G4Hits_sPHENIX_pi-_eta0_16GeV-0000.root_qa.rootQA_Draw_HCALIN_G4Hit.pdf
0312
0313     if (G4HCALIN::inner_hcal_material_Al) visible_sample_fraction_HCALIN = 0.162166;  // for "G4_Al", Abhisek Sen <sen.abhisek@gmail.com>
0314
0315     RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HcalInRawTowerCalibration");
0316     TowerCalibration->Detector("HCALIN");
0317     TowerCalibration->set_usetowerinfo_v2(G4HCALIN::useTowerInfoV2);
0318     //  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0319     //  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0320     TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0321     if (G4HCALIN::TowerDigi == RawTowerDigitizer::kNo_digitization)
0322     {
0323       // 0.176 extracted from electron sims (edep(scintillator)/edep(total))
0324       TowerCalibration->set_calib_const_GeV_ADC(1. / 0.176);
0325     }
0326     else
0327     {
0328       TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / visible_sample_fraction_HCALIN);
0329     }
0330     TowerCalibration->set_pedstal_ADC(0);
0331     TowerCalibration->Verbosity(verbosity);
0332     if (!Enable::HCALIN_G4Hit) TowerCalibration->set_towerinfo(RawTowerCalibration::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0333     se->registerSubsystem(TowerCalibration);
0334   }
0335   else
0336   {
0337     CaloWaveformSim *caloWaveformSim = new CaloWaveformSim();
0338     caloWaveformSim->set_detector_type(CaloTowerDefs::HCALIN);
0339     caloWaveformSim->set_detector("HCALIN");
0340     caloWaveformSim->set_nsamples(12);
0341     caloWaveformSim->set_pedestalsamples(12);
0342     caloWaveformSim->set_timewidth(0.2);
0343     caloWaveformSim->set_peakpos(6);
0344     // caloWaveformSim->Verbosity(2);
0345     // caloWaveformSim->set_noise_type(CaloWaveformSim::NOISE_NONE);
0346     se->registerSubsystem(caloWaveformSim);
0347
0348     CaloTowerBuilder *ca2 = new CaloTowerBuilder();
0349     ca2->set_detector_type(CaloTowerDefs::HCALIN);
0350     ca2->set_nsamples(12);
0351     ca2->set_dataflag(false);
0352     ca2->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0353     ca2->set_builder_type(CaloTowerDefs::kWaveformTowerSimv1);
0354     ca2->set_softwarezerosuppression(true, 30);
0355     se->registerSubsystem(ca2);
0356
0357     CaloTowerStatus *statusHCALIN = new CaloTowerStatus("HCALINSTATUS");
0358     statusHCALIN->set_detector_type(CaloTowerDefs::HCALIN);
0359     statusHCALIN->set_time_cut(2);
0360     se->registerSubsystem(statusHCALIN);
0361
0362     CaloTowerCalib *calibIHCal = new CaloTowerCalib("HCALINCALIB");
0363     calibIHCal->set_detector_type(CaloTowerDefs::HCALIN);
0364     calibIHCal->set_outputNodePrefix("TOWERINFO_CALIB_");
0365     se->registerSubsystem(calibIHCal);
0366   }
0367   return;
0368 }
0369
0370 void HCALInner_Clusters()
0371 {
0372   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0373
0374   Fun4AllServer *se = Fun4AllServer::instance();
0375
0376   if (G4HCALIN::HCalIn_clusterizer == G4HCALIN::kHCalInTemplateClusterizer)
0377   {
0378     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("HcalInRawClusterBuilderTemplate");
0379     ClusterBuilder->Detector("HCALIN");
0380     ClusterBuilder->SetCylindricalGeometry();  // has to be called after Detector()
0381     ClusterBuilder->Verbosity(verbosity);
0382     if (!Enable::HCALIN_G4Hit || Enable::HCALIN_TOWERINFO) ClusterBuilder->set_UseTowerInfo(1);  // just use towerinfo
0383     se->registerSubsystem(ClusterBuilder);
0384   }
0385   else if (G4HCALIN::HCalIn_clusterizer == G4HCALIN::kHCalInGraphClusterizer)
0386   {
0387     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("HcalInRawClusterBuilderGraph");
0388     ClusterBuilder->Detector("HCALIN");
0389     ClusterBuilder->Verbosity(verbosity);
0390     // if (!Enable::HCALIN_G4Hit) ClusterBuilder->set_UseTowerInfo(1);  // just use towerinfo
0391     se->registerSubsystem(ClusterBuilder);
0392   }
0393   else
0394   {
0395     std::cout << "HCalIn_Clusters - unknown clusterizer setting!" << std::endl;
0396     exit(1);
0397   }
0398   return;
0399 }
0400
0401 void HCALInner_Eval(const std::string &outputfile, int start_event = 0)
0402 {
0403   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0404   Fun4AllServer *se = Fun4AllServer::instance();
0405
0406   CaloEvaluator *eval = new CaloEvaluator("HCALINEVALUATOR", "HCALIN", outputfile);
0407   eval->set_event(start_event);
0408   eval->Verbosity(verbosity);
0409   eval->set_use_towerinfo(Enable::HCALIN_TOWERINFO);
0410   se->registerSubsystem(eval);
0411
0412   return;
0413 }
0414
0415 void HCALInner_QA()
0416 {
0417   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::HCALIN_VERBOSITY);
0418
0419   Fun4AllServer *se = Fun4AllServer::instance();
0420   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("HCALIN");
0421   if (Enable::HCALIN_TOWERINFO) qa->set_flags(QAG4SimulationCalorimeter::kProcessTowerinfo);
0422   qa->Verbosity(verbosity);
0423   se->registerSubsystem(qa);
0424
0425   return;
0426 }
0427
0428 #endif