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