macros/common/G4_HcalOut_ref.C

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