macros/common/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 <g4ohcal/PHG4OHCalSubsystem.h>
0011
0012 #include <g4detectors/PHG4HcalCellReco.h>
0013 #include <g4detectors/PHG4OuterHcalSubsystem.h>
0014
0015 #include <g4eval/CaloEvaluator.h>
0016
0017 #include <g4main/PHG4Reco.h>
0018
0019 #include <calowaveformsim/CaloWaveformSim.h>
0020
0021 #include <calobase/TowerInfoDefs.h>
0022
0023 #include <caloreco/CaloTowerBuilder.h>
0024 #include <caloreco/CaloTowerCalib.h>
0025 #include <caloreco/CaloTowerStatus.h>
0026 #include <caloreco/CaloWaveformProcessing.h>
0027 #include <caloreco/RawClusterBuilderGraph.h>
0028 #include <caloreco/RawClusterBuilderTemplate.h>
0029 #include <caloreco/RawTowerCalibration.h>
0030
0031 #include <simqa_modules/QAG4SimulationCalorimeter.h>
0032
0033 #include <fun4all/Fun4AllServer.h>
0034
0035 R__LOAD_LIBRARY(libcalo_reco.so)
0036 R__LOAD_LIBRARY(libg4calo.so)
0037 R__LOAD_LIBRARY(libCaloWaveformSim.so)
0038 R__LOAD_LIBRARY(libg4detectors.so)
0039 R__LOAD_LIBRARY(libg4eval.so)
0040 R__LOAD_LIBRARY(libg4ohcal.so)
0041 R__LOAD_LIBRARY(libsimqa_modules.so)
0042
0043 namespace Enable
0044 {
0045   bool HCALOUT = false;
0046   bool HCALOUT_ABSORBER = false;
0047   bool HCALOUT_OVERLAPCHECK = false;
0048   bool HCALOUT_CELL = false;
0049   bool HCALOUT_TOWER = false;
0050   bool HCALOUT_CLUSTER = false;
0051   bool HCALOUT_EVAL = false;
0052   bool HCALOUT_QA = false;
0053   bool HCALOUT_OLD = false;
0054   bool HCALOUT_RING = false;
0055   bool HCALOUT_G4Hit = true;
0056   bool HCALOUT_TOWERINFO = false;
0057   int HCALOUT_VERBOSITY = 0;
0058 }  // namespace Enable
0059
0060 namespace G4HCALOUT
0061 {
0062   double outer_radius = 269.317 + 5;
0063   double size_z = 639.240 + 10;
0064   double phistart = NAN;
0065   double tower_emin = NAN;
0066   int light_scint_model = -1;
0067   int tower_energy_source = -1;
0068
0069   // Digitization (default photon digi):
0070   RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0071   // directly pass the energy of sim tower to digitized tower
0072   // kNo_digitization
0073   // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
0074   // kSimple_photon_digitization
0075   // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
0076   // kSiPM_photon_digitization
0077
0078   enum enu_HCalOut_clusterizer
0079   {
0080     kHCalOutGraphClusterizer,
0081     kHCalOutTemplateClusterizer
0082   };
0083
0084   bool useTowerInfoV2 = true;
0085
0086   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0087   enu_HCalOut_clusterizer HCalOut_clusterizer = kHCalOutTemplateClusterizer;
0088   //! graph clusterizer, RawClusterBuilderGraph
0089   // enu_HCalOut_clusterizer HCalOut_clusterizer = kHCalOutGraphClusterizer;
0090 }  // namespace G4HCALOUT
0091
0092 // Init is called by G4Setup.C
0093 void HCalOuterInit()
0094 {
0095   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4HCALOUT::outer_radius);
0096   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4HCALOUT::size_z / 2.);
0097   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4HCALOUT::size_z / 2.);
0098 }
0099
0100 double HCalOuter(PHG4Reco *g4Reco,
0101                  double radius,
0102                  const int crossings)
0103 {
0104   bool AbsorberActive = Enable::ABSORBER || Enable::HCALOUT_ABSORBER;
0105   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::HCALOUT_OVERLAPCHECK;
0106   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0107
0108   PHG4DetectorSubsystem *hcal = nullptr;
0109   //  Mephi Maps
0110   //  Maps are different for old/new but how to set is identical
0111   //  here are the ones for the old outer hcal since the new maps do not exist yet
0112   //  use hcal->set_string_param("MapFileName",""); to disable map
0113   //  hcal->set_string_param("MapFileName",std::string(getenv("CALIBRATIONROOT")) + "/HCALOUT/tilemap/oHCALMaps092021.root");
0114   //  hcal->set_string_param("MapHistoName","hCombinedMap");
0115
0116   if (Enable::HCALOUT_OLD)
0117   {
0118     hcal = new PHG4OuterHcalSubsystem("HCALOUT");
0119     // hcal->set_double_param("inner_radius", 183.3);
0120     //-----------------------------------------
0121     // the light correction can be set in a single call
0122     // hcal->set_double_param("light_balance_inner_corr", NAN);
0123     // hcal->set_double_param("light_balance_inner_radius", NAN);
0124     // hcal->set_double_param("light_balance_outer_corr", NAN);
0125     // hcal->set_double_param("light_balance_outer_radius", NAN);
0126     // hcal->set_double_param("magnet_cutout_radius", 195.31);
0127     // hcal->set_double_param("magnet_cutout_scinti_radius", 195.96);
0128     // hcal->SetLightCorrection(NAN,NAN,NAN,NAN);
0129     //-----------------------------------------
0130     // hcal->set_double_param("outer_radius", G4HCALOUT::outer_radius);
0131     // hcal->set_double_param("place_x", 0.);
0132     // hcal->set_double_param("place_y", 0.);
0133     // hcal->set_double_param("place_z", 0.);
0134     // hcal->set_double_param("rot_x", 0.);
0135     // hcal->set_double_param("rot_y", 0.);
0136     // hcal->set_double_param("rot_z", 0.);
0137     // hcal->set_double_param("scinti_eta_coverage", 1.1);
0138     // hcal->set_double_param("scinti_gap", 0.85);
0139     // hcal->set_double_param("scinti_gap_neighbor", 0.1);
0140     // hcal->set_double_param("scinti_inner_radius",183.89);
0141     // hcal->set_double_param("scinti_outer_radius",263.27);
0142     // hcal->set_double_param("scinti_tile_thickness", 0.7);
0143     // hcal->set_double_param("size_z", G4HCALOUT::size_z);
0144     // hcal->set_double_param("steplimits", NAN);
0145     // hcal->set_double_param("tilt_angle", -11.23);
0146
0147     // hcal->set_int_param("light_scint_model", 1);
0148     // hcal->set_int_param("magnet_cutout_first_scinti", 8);
0149     // hcal->set_int_param("ncross", 0);
0150     // hcal->set_int_param("n_towers", 64);
0151     // hcal->set_int_param("n_scinti_plates_per_tower", 5);
0152     // hcal->set_int_param("n_scinti_tiles", 12);
0153
0154     // hcal->set_string_param("material", "Steel_1006");
0155   }
0156   else
0157   {
0158     hcal = new PHG4OHCalSubsystem("HCALOUT");
0159     if (Enable::HCALOUT_RING)
0160     {
0161       std::string gdmlfile_no_ring = std::string(getenv("CALIBRATIONROOT")) + "/HcalGeo/OuterHCalAbsorberTiles_merged.gdml";
0162       hcal->set_string_param("GDMPath", gdmlfile_no_ring);
0163     }
0164     // hcal->set_string_param("GDMPath", "mytestgdml.gdml"); // try other gdml file
0165     // common setting with tracking, we likely want to move to the cdb with this
0166     hcal->set_string_param("IronFieldMapPath", G4MAGNET::magfield_OHCAL_steel);
0167     hcal->set_double_param("IronFieldMapScale", G4MAGNET::magfield_rescale);
0168   }
0169
0170   if (G4HCALOUT::light_scint_model >= 0)
0171   {
0172     hcal->set_int_param("light_scint_model", G4HCALOUT::light_scint_model);
0173   }
0174   // hcal->set_int_param("field_check", 1); // for validating the field in HCal
0175   hcal->SetActive();
0176   hcal->SuperDetector("HCALOUT");
0177   if (AbsorberActive)
0178   {
0179     hcal->SetAbsorberActive();
0180   }
0181   hcal->OverlapCheck(OverlapCheck);
0182   if (!std::isfinite(G4HCALOUT::phistart))
0183   {
0184     if (Enable::HCALOUT_OLD)
0185     {
0186       G4HCALOUT::phistart = 0.026598397;  // offet in phi (from zero) extracted from geantinos
0187     }
0188     else
0189     {
0190       G4HCALOUT::phistart = 0.0240615415;  // offet in phi (from zero) extracted from geantinos
0191     }
0192   }
0193   hcal->set_int_param("saveg4hit", Enable::HCALOUT_G4Hit);
0194   hcal->set_double_param("phistart", G4HCALOUT::phistart);
0195   g4Reco->registerSubsystem(hcal);
0196
0197   if (!Enable::HCALOUT_OLD)
0198   {
0199     // HCal support rings, approximated as solid rings
0200     // note there is only one ring on either side, but to allow part of the ring inside the HCal envelope two rings are used
0201     const double inch = 2.54;
0202     const double support_ring_outer_radius = 74.061 * inch;
0203     const double innerradius = 56.188 * inch;
0204     const double hcal_envelope_radius = 182.423 - 5.;
0205     const double support_ring_z = 175.375 * inch / 2.;
0206     const double support_ring_dz = 4. * inch;
0207     const double z_rings[] =
0208         {-support_ring_z, support_ring_z};
0209     PHG4CylinderSubsystem *cyl;
0210     PHG4CylinderSubsystem *cylout;
0211
0212     for (int i = 0; i < 2; i++)
0213     {
0214       // rings outside of HCal envelope
0215       cyl = new PHG4CylinderSubsystem("HCAL_SPT_N1", i);
0216       cyl->set_double_param("place_z", z_rings[i]);
0217       cyl->SuperDetector("HCALIN_SPT");
0218       cyl->set_double_param("radius", innerradius);
0219       cyl->set_int_param("lengthviarapidity", 0);
0220       cyl->set_double_param("length", support_ring_dz);
0221       cyl->set_string_param("material", "G4_Al");
0222       cyl->set_double_param("thickness", hcal_envelope_radius - 0.1 - innerradius);
0223       cyl->set_double_param("start_phi_rad", 1.867);
0224       cyl->set_double_param("delta_phi_rad", 5.692);
0225       cyl->OverlapCheck(Enable::OVERLAPCHECK);
0226       if (AbsorberActive)
0227       {
0228         cyl->SetActive();
0229       }
0230       g4Reco->registerSubsystem(cyl);
0231
0232       // rings inside outer HCal envelope
0233       // only use if we want to use the old version of the ring instead of the gdml implementation
0234       if (Enable::HCALOUT_RING)
0235       {
0236         cylout = new PHG4CylinderSubsystem("HCAL_SPT_N1", i + 2);
0237         cylout->set_double_param("place_z", z_rings[i]);
0238         cylout->SuperDetector("HCALIN_SPT");
0239         cylout->set_double_param("radius", hcal_envelope_radius + 0.1);  // add a mm to avoid overlaps
0240         cylout->set_int_param("lengthviarapidity", 0);
0241         cylout->set_double_param("length", support_ring_dz);
0242         cylout->set_string_param("material", "G4_Al");
0243         cylout->set_double_param("thickness", support_ring_outer_radius - (hcal_envelope_radius + 0.1));
0244         cylout->set_double_param("start_phi_rad", 1.867);
0245         cylout->set_double_param("delta_phi_rad", 5.692);
0246         if (AbsorberActive)
0247         {
0248           cylout->SetActive();
0249         }
0250         cylout->SetMotherSubsystem(hcal);
0251         cylout->OverlapCheck(OverlapCheck);
0252         g4Reco->registerSubsystem(cylout);
0253       }
0254     }
0255   }
0256
0257   radius = hcal->get_double_param("outer_radius");
0258
0259   radius += no_overlapp;
0260
0261   return radius;
0262 }
0263
0264 void HCALOuter_Cells()
0265 {
0266   if (!Enable::HCALOUT_G4Hit) return;
0267   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0268
0269   Fun4AllServer *se = Fun4AllServer::instance();
0270
0271   PHG4HcalCellReco *hc = new PHG4HcalCellReco("HCALOUT_CELLRECO");
0272   hc->Detector("HCALOUT");
0273   hc->Verbosity(verbosity);
0274   // check for energy conservation - needs modified "infinite" timing cuts
0275   // 0-999999999
0276   //  hc->checkenergy();
0277   // timing cuts with their default settings
0278   // hc->set_double_param("tmin",0.);
0279   // hc->set_double_param("tmax",60.0);
0280   // or all at once:
0281   // hc->set_timing_window(0.0,60.0);
0282   // this sets all cells to a fixed energy for debugging
0283   // hc->set_fixed_energy(1.);
0284   se->registerSubsystem(hc);
0285
0286   return;
0287 }
0288
0289 void HCALOuter_Towers()
0290 {
0291   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0292   Fun4AllServer *se = Fun4AllServer::instance();
0293
0294   if (!Enable::HCALOUT_TOWERINFO)
0295   {
0296     // build the raw tower anyways for the geom nodes
0297   if (Enable::HCALOUT_G4Hit)
0298     {
0299       HcalRawTowerBuilder *TowerBuilder = new HcalRawTowerBuilder("HcalOutRawTowerBuilder");
0300       TowerBuilder->Detector("HCALOUT");
0301       TowerBuilder->set_sim_tower_node_prefix("SIM");
0302       if (!isfinite(G4HCALOUT::phistart))
0303       {
0304         if (Enable::HCALOUT_OLD)
0305         {
0306           G4HCALOUT::phistart = 0.026598397;  // offet in phi (from zero) extracted from geantinos
0307         }
0308         else
0309         {
0310           G4HCALOUT::phistart = 0.0240615415;  // offet in phi (from zero) extracted from geantinos
0311         }
0312       }
0313       TowerBuilder->set_double_param("phistart", G4HCALOUT::phistart);
0314       if (isfinite(G4HCALOUT::tower_emin))
0315       {
0316         TowerBuilder->set_double_param("emin", G4HCALOUT::tower_emin);
0317       }
0318       if (G4HCALOUT::tower_energy_source >= 0)
0319       {
0320         TowerBuilder->set_int_param("tower_energy_source", G4HCALOUT::tower_energy_source);
0321       }
0322       // this sets specific decalibration factors
0323       // for a given cell
0324       // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0325       // for a whole tower
0326       // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0327       // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0328       // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0329       TowerBuilder->Verbosity(verbosity);
0330       se->registerSubsystem(TowerBuilder);
0331     }
0332     // From 2016 Test beam sim
0333     RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HcalOutRawTowerDigitizer");
0334     TowerDigitizer->Detector("HCALOUT");
0335     //  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0336     TowerDigitizer->set_digi_algorithm(G4HCALOUT::TowerDigi);
0337     TowerDigitizer->set_pedstal_central_ADC(0);
0338     TowerDigitizer->set_pedstal_width_ADC(1);  // From Jin's guess. No EMCal High Gain data yet! TODO: update
0339     TowerDigitizer->set_photonelec_ADC(16. / 5.);
0340     TowerDigitizer->set_photonelec_yield_visible_GeV(16. / 5 / (0.2e-3));
0341     TowerDigitizer->set_zero_suppression_ADC(-0);  // no-zero suppression
0342     TowerDigitizer->Verbosity(verbosity);
0343     if (!Enable::HCALOUT_G4Hit) TowerDigitizer->set_towerinfo(RawTowerDigitizer::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0344     se->registerSubsystem(TowerDigitizer);
0345
0346     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
0347
0348     RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HcalOutRawTowerCalibration");
0349     TowerCalibration->Detector("HCALOUT");
0350     TowerCalibration->set_usetowerinfo_v2(G4HCALOUT::useTowerInfoV2);
0351
0352     //  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0353     //  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0354     TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0355     if (G4HCALOUT::TowerDigi == RawTowerDigitizer::kNo_digitization)
0356     {
0357       // 0.033 extracted from electron sims (edep(scintillator)/edep(total))
0358       TowerCalibration->set_calib_const_GeV_ADC(1. / 0.033);
0359     }
0360     else
0361     {
0362       TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / visible_sample_fraction_HCALOUT);
0363     }
0364     TowerCalibration->set_pedstal_ADC(0);
0365     TowerCalibration->Verbosity(verbosity);
0366     if (!Enable::HCALOUT_G4Hit) TowerCalibration->set_towerinfo(RawTowerCalibration::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0367     se->registerSubsystem(TowerCalibration);
0368   }
0369   // where I use waveformsim
0370   else
0371   {
0372     CaloWaveformSim *caloWaveformSim = new CaloWaveformSim();
0373     caloWaveformSim->set_detector_type(CaloTowerDefs::HCALOUT);
0374     caloWaveformSim->set_detector("HCALOUT");
0375     caloWaveformSim->set_nsamples(12);
0376     caloWaveformSim->set_pedestalsamples(12);
0377     caloWaveformSim->set_timewidth(0.2);
0378     caloWaveformSim->set_peakpos(6);
0379     // caloWaveformSim->Verbosity(2);
0380     // caloWaveformSim->set_noise_type(CaloWaveformSim::NOISE_NONE);
0381     se->registerSubsystem(caloWaveformSim);
0382
0383     CaloTowerBuilder *ca2 = new CaloTowerBuilder();
0384     ca2->set_detector_type(CaloTowerDefs::HCALOUT);
0385     ca2->set_nsamples(12);
0386     ca2->set_dataflag(false);
0387     ca2->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0388     ca2->set_builder_type(CaloTowerDefs::kWaveformTowerSimv1);
0389     ca2->set_softwarezerosuppression(true, 30);
0390     se->registerSubsystem(ca2);
0391
0392     CaloTowerStatus *statusHCALOUT = new CaloTowerStatus("HCALOUTSTATUS");
0393     statusHCALOUT->set_detector_type(CaloTowerDefs::HCALOUT);
0394     statusHCALOUT->set_time_cut(2);
0395     se->registerSubsystem(statusHCALOUT);
0396
0397     CaloTowerCalib *calibOHCal = new CaloTowerCalib("HCALOUTCALIB");
0398     calibOHCal->set_detector_type(CaloTowerDefs::HCALOUT);
0399     calibOHCal->set_outputNodePrefix("TOWERINFO_CALIB_");
0400     se->registerSubsystem(calibOHCal);
0401   }
0402
0403   return;
0404 }
0405
0406 void HCALOuter_Clusters()
0407 {
0408   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0409
0410   Fun4AllServer *se = Fun4AllServer::instance();
0411
0412   if (G4HCALOUT::HCalOut_clusterizer == G4HCALOUT::kHCalOutTemplateClusterizer)
0413   {
0414     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("HcalOutRawClusterBuilderTemplate");
0415     ClusterBuilder->Detector("HCALOUT");
0416     ClusterBuilder->SetCylindricalGeometry();  // has to be called after Detector()
0417     ClusterBuilder->Verbosity(verbosity);
0418     if (!Enable::HCALOUT_G4Hit || Enable::HCALOUT_TOWERINFO) ClusterBuilder->set_UseTowerInfo(1);  // just use towerinfo
0419     se->registerSubsystem(ClusterBuilder);
0420   }
0421   else if (G4HCALOUT::HCalOut_clusterizer == G4HCALOUT::kHCalOutGraphClusterizer)
0422   {
0423     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("HcalOutRawClusterBuilderGraph");
0424     ClusterBuilder->Detector("HCALOUT");
0425     ClusterBuilder->Verbosity(verbosity);
0426     // if (!Enable::HCALOUT_G4Hit) ClusterBuilder->set_UseTowerInfo(1);  // just use towerinfo
0427     se->registerSubsystem(ClusterBuilder);
0428   }
0429   else
0430   {
0431     std::cout << "HCALOuter_Clusters - unknown clusterizer setting!" << std::endl;
0432     exit(1);
0433   }
0434
0435   return;
0436 }
0437
0438 void HCALOuter_Eval(const std::string &outputfile, int start_event = 0)
0439 {
0440   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALOUT_VERBOSITY);
0441
0442   Fun4AllServer *se = Fun4AllServer::instance();
0443
0444   CaloEvaluator *eval = new CaloEvaluator("HCALOUTEVALUATOR", "HCALOUT", outputfile);
0445   eval->set_event(start_event);
0446   eval->Verbosity(verbosity);
0447   eval->set_use_towerinfo(Enable::HCALOUT_TOWERINFO);
0448   se->registerSubsystem(eval);
0449
0450   return;
0451 }
0452
0453 void HCALOuter_QA()
0454 {
0455   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::HCALOUT_VERBOSITY);
0456
0457   Fun4AllServer *se = Fun4AllServer::instance();
0458   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("HCALOUT");
0459   if (Enable::HCALOUT_TOWERINFO) qa->set_flags(QAG4SimulationCalorimeter::kProcessTowerinfo);
0460   qa->Verbosity(verbosity);
0461   se->registerSubsystem(qa);
0462
0463   return;
0464 }
0465
0466 #endif