calo/tower_slope_macro/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 <g4detectors/PHG4CylinderSubsystem.h>
0012 #include <g4detectors/PHG4HcalCellReco.h>
0013 #include <g4detectors/PHG4InnerHcalSubsystem.h>
0014
0015 #include <g4eval/CaloEvaluator.h>
0016
0017 #include <g4main/PHG4Reco.h>
0018
0019 #include <caloreco/RawClusterBuilderGraph.h>
0020 #include <caloreco/RawClusterBuilderTemplate.h>
0021 #include <caloreco/RawTowerCalibration.h>
0022 #include <qa_modules/QAG4SimulationCalorimeter.h>
0023
0024 #include <fun4all/Fun4AllServer.h>
0025
0026 R__LOAD_LIBRARY(libcalo_reco.so)
0027 R__LOAD_LIBRARY(libg4calo.so)
0028 R__LOAD_LIBRARY(libg4detectors.so)
0029 R__LOAD_LIBRARY(libg4eval.so)
0030 R__LOAD_LIBRARY(libqa_modules.so)
0031
0032 void HCalInner_SupportRing(PHG4Reco *g4Reco);
0033
0034 namespace Enable
0035 {
0036   bool HCALIN = false;
0037   bool HCALIN_ABSORBER = false;
0038   bool HCALIN_OVERLAPCHECK = false;
0039   bool HCALIN_CELL = false;
0040   bool HCALIN_TOWER = false;
0041   bool HCALIN_CLUSTER = false;
0042   bool HCALIN_EVAL = false;
0043   bool HCALIN_QA = false;
0044   bool HCALIN_SUPPORT = false;
0045   int HCALIN_VERBOSITY = 0;
0046 }  // namespace Enable
0047
0048 namespace G4HCALIN
0049 {
0050   double support_ring_outer_radius = 178.0 - 0.001;
0051   double support_ring_z_ring2 = (2150 + 2175) / 2. / 10.;
0052   double dz = 25. / 10.;
0053
0054   //Inner HCal absorber material selector:
0055   //false - old version, absorber material is SS310
0056   //true - default Choose if you want Aluminum
0057   bool inner_hcal_material_Al = true;
0058
0059   int inner_hcal_eic = 0;
0060
0061   // Digitization (default photon digi):
0062   RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0063   // directly pass the energy of sim tower to digitized tower
0064   // kNo_digitization
0065   // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
0066   // kSimple_photon_digitization
0067   // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
0068   // kSiPM_photon_digitization
0069
0070   enum enu_HCalIn_clusterizer
0071   {
0072     kHCalInGraphClusterizer,
0073
0074     kHCalInTemplateClusterizer
0075   };
0076
0077   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0078   enu_HCalIn_clusterizer HCalIn_clusterizer = kHCalInTemplateClusterizer;
0079   //! graph clusterizer, RawClusterBuilderGraph
0080   //enu_HCalIn_clusterizer HCalIn_clusterizer = kHCalInGraphClusterizer;
0081 }  // namespace G4HCALIN
0082
0083 // Init is called by G4Setup.C
0084 void HCalInnerInit(const int iflag = 0)
0085 {
0086   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4HCALIN::support_ring_outer_radius);
0087   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4HCALIN::support_ring_z_ring2 + G4HCALIN::dz / 2.);
0088   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -G4HCALIN::support_ring_z_ring2 - G4HCALIN::dz / 2.);
0089   if (iflag == 1)
0090   {
0091     G4HCALIN::inner_hcal_eic = 1;
0092   }
0093 }
0094
0095 double HCalInner(PHG4Reco *g4Reco,
0096                  double radius,
0097                  const int crossings)
0098 {
0099   bool AbsorberActive = Enable::ABSORBER || Enable::HCALIN_ABSORBER;
0100   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::HCALIN_OVERLAPCHECK;
0101   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0102
0103   // all sizes are in cm!
0104   PHG4InnerHcalSubsystem *hcal = new PHG4InnerHcalSubsystem("HCALIN");
0105   // these are the parameters you can change with their default settings
0106   // hcal->set_string_param("material","SS310");
0107   if (G4HCALIN::inner_hcal_material_Al)
0108   {
0109     if (verbosity > 0)
0110     {
0111       cout << "HCalInner - construct inner HCal absorber with G4_Al" << endl;
0112     }
0113     hcal->set_string_param("material", "G4_Al");
0114   }
0115   else
0116   {
0117     if (verbosity > 0)
0118     {
0119       cout << "HCalInner - construct inner HCal absorber with SS310" << endl;
0120     }
0121     hcal->set_string_param("material", "SS310");
0122   }
0123   // hcal->set_double_param("inner_radius", 117.27);
0124   //-----------------------------------------
0125   // the light correction can be set in a single call
0126   // hcal->set_double_param("light_balance_inner_corr", NAN);
0127   // hcal->set_double_param("light_balance_inner_radius", NAN);
0128   // hcal->set_double_param("light_balance_outer_corr", NAN);
0129   // hcal->set_double_param("light_balance_outer_radius", NAN);
0130   // hcal->SetLightCorrection(NAN,NAN,NAN,NAN);
0131   //-----------------------------------------
0132   // hcal->set_double_param("outer_radius", 134.42);
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_neighbor", 0.1);
0141   // hcal->set_double_param("scinti_inner_gap", 0.85);
0142   // hcal->set_double_param("scinti_outer_gap", 1.22 * (5.0 / 4.0));
0143   // hcal->set_double_param("scinti_outer_radius", 133.3);
0144   // hcal->set_double_param("scinti_tile_thickness", 0.7);
0145   // hcal->set_double_param("size_z", 175.94 * 2);
0146   // hcal->set_double_param("steplimits", NAN);
0147   // hcal->set_double_param("tilt_angle", 36.15);
0148
0149   // hcal->set_int_param("light_scint_model", 1);
0150   // hcal->set_int_param("ncross", 0);
0151   // hcal->set_int_param("n_towers", 64);
0152   // hcal->set_int_param("n_scinti_plates_per_tower", 4);
0153   // hcal->set_int_param("n_scinti_tiles", 12);
0154
0155   // hcal->set_string_param("material", "SS310");
0156
0157   hcal->SetActive();
0158   hcal->SuperDetector("HCALIN");
0159   if (AbsorberActive)
0160   {
0161     hcal->SetAbsorberActive();
0162   }
0163   hcal->OverlapCheck(OverlapCheck);
0164
0165   g4Reco->registerSubsystem(hcal);
0166
0167   radius = hcal->get_double_param("outer_radius");
0168
0169   HCalInner_SupportRing(g4Reco);
0170
0171   radius += no_overlapp;
0172   return radius;
0173 }
0174
0175 //! A rough version of the inner HCal support ring, from Richie's CAD drawing. - Jin
0176 void HCalInner_SupportRing(PHG4Reco *g4Reco)
0177 {
0178   bool AbsorberActive = Enable::SUPPORT || Enable::HCALIN_SUPPORT;
0179
0180   const double z_ring1 = (2025 + 2050) / 2. / 10.;
0181   const double innerradius_sphenix = 116.;
0182   const double innerradius_ephenix_hadronside = 138.;
0183   const double z_rings[] =
0184       {-G4HCALIN::support_ring_z_ring2, -z_ring1, z_ring1, G4HCALIN::support_ring_z_ring2};
0185
0186   PHG4CylinderSubsystem *cyl;
0187
0188   for (int i = 0; i < 4; i++)
0189   {
0190     double innerradius = innerradius_sphenix;
0191     if (z_rings[i] > 0 && G4HCALIN::inner_hcal_eic == 1)
0192     {
0193       innerradius = innerradius_ephenix_hadronside;
0194     }
0195     cyl = new PHG4CylinderSubsystem("HCALIN_SPT_N1", i);
0196     cyl->set_double_param("place_z", z_rings[i]);
0197     cyl->SuperDetector("HCALIN_SPT");
0198     cyl->set_double_param("radius", innerradius);
0199     cyl->set_int_param("lengthviarapidity", 0);
0200     cyl->set_double_param("length", G4HCALIN::dz);
0201     cyl->set_string_param("material", "SS310");
0202     cyl->set_double_param("thickness", G4HCALIN::support_ring_outer_radius - innerradius);
0203     if (AbsorberActive)
0204     {
0205       cyl->SetActive();
0206     }
0207     g4Reco->registerSubsystem(cyl);
0208   }
0209
0210   return;
0211 }
0212
0213 void HCALInner_Cells()
0214 {
0215   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0216
0217   Fun4AllServer *se = Fun4AllServer::instance();
0218
0219   PHG4HcalCellReco *hc = new PHG4HcalCellReco("HCALIN_CELLRECO");
0220   hc->Detector("HCALIN");
0221   //  hc->Verbosity(2);
0222   // check for energy conservation - needs modified "infinite" timing cuts
0223   // 0-999999999
0224   //  hc->checkenergy();
0225   // timing cuts with their default settings
0226   // hc->set_double_param("tmin",0.);
0227   // hc->set_double_param("tmax",60.0);
0228   // or all at once:
0229   // hc->set_timing_window(0.0,60.0);
0230   // this sets all cells to a fixed energy for debugging
0231   // hc->set_fixed_energy(1.);
0232   se->registerSubsystem(hc);
0233
0234   return;
0235 }
0236
0237 void HCALInner_Towers()
0238 {
0239   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0240   Fun4AllServer *se = Fun4AllServer::instance();
0241
0242   HcalRawTowerBuilder *TowerBuilder = new HcalRawTowerBuilder("HcalInRawTowerBuilder");
0243   TowerBuilder->Detector("HCALIN");
0244   TowerBuilder->set_sim_tower_node_prefix("SIM");
0245   // this sets specific decalibration factors
0246   // for a given cell
0247   // TowerBuilder->set_cell_decal_factor(1,10,0.1);
0248   // for a whole tower
0249   // TowerBuilder->set_tower_decal_factor(0,10,0.2);
0250   TowerBuilder->Verbosity(verbosity);
0251   se->registerSubsystem(TowerBuilder);
0252
0253   // From 2016 Test beam sim
0254   RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("HcalInRawTowerDigitizer");
0255   TowerDigitizer->Detector("HCALIN");
0256   //  TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0257   TowerDigitizer->set_digi_algorithm(G4HCALIN::TowerDigi);
0258   TowerDigitizer->set_pedstal_central_ADC(0);
0259   TowerDigitizer->set_pedstal_width_ADC(1);  // From Jin's guess. No EMCal High Gain data yet! TODO: update
0260   TowerDigitizer->set_photonelec_ADC(32. / 5.);
0261   TowerDigitizer->set_photonelec_yield_visible_GeV(32. / 5 / (0.4e-3));
0262   TowerDigitizer->set_zero_suppression_ADC(-0);  // no-zero suppression
0263   se->registerSubsystem(TowerDigitizer);
0264
0265   //Default sampling fraction for SS310
0266   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
0267
0268   if (G4HCALIN::inner_hcal_material_Al) visible_sample_fraction_HCALIN = 0.162166;  //for "G4_Al", Abhisek Sen <sen.abhisek@gmail.com>
0269
0270   RawTowerCalibration *TowerCalibration = new RawTowerCalibration("HcalInRawTowerCalibration");
0271   TowerCalibration->Detector("HCALIN");
0272   //  TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0273   //  TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0274
0275   //  TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0276
0277   TowerCalibration->set_calib_algorithm(RawTowerCalibration::kDbfile_tbt_gain_corr);
0278   TowerCalibration->set_UseConditionsDB(false);
0279   // see documentation in github/jtaou/coresoft
0280   // ...[calib/towerslope] /macros/G4_CEmc_Spacal.C
0281   // for more info about TowerCalibration dbfile mode
0282
0283   //  TowerCalibration->set_CalibrationFileName("HCALIN_GainsCalib1.22.txt");
0284   TowerCalibration->set_CalibrationFileName("HCALIN_GainsCalib1.12_PosEtaOnly_NegEta1.0.txt");
0285
0286
0287   if (G4HCALIN::TowerDigi == RawTowerDigitizer::kNo_digitization)
0288   {
0289     // 0.176 extracted from electron sims (edep(scintillator)/edep(total))
0290     TowerCalibration->set_calib_const_GeV_ADC(1. / 0.176);
0291   }
0292   else
0293   {
0294     TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / visible_sample_fraction_HCALIN);
0295   }
0296   TowerCalibration->set_pedstal_ADC(0);
0297   se->registerSubsystem(TowerCalibration);
0298
0299   return;
0300 }
0301
0302 void HCALInner_Clusters()
0303 {
0304   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0305
0306   Fun4AllServer *se = Fun4AllServer::instance();
0307
0308   if (G4HCALIN::HCalIn_clusterizer == G4HCALIN::kHCalInTemplateClusterizer)
0309   {
0310     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("HcalInRawClusterBuilderTemplate");
0311     ClusterBuilder->Detector("HCALIN");
0312     ClusterBuilder->SetCylindricalGeometry();  // has to be called after Detector()
0313     ClusterBuilder->Verbosity(verbosity);
0314     se->registerSubsystem(ClusterBuilder);
0315   }
0316   else if (G4HCALIN::HCalIn_clusterizer == G4HCALIN::kHCalInGraphClusterizer)
0317   {
0318     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("HcalInRawClusterBuilderGraph");
0319     ClusterBuilder->Detector("HCALIN");
0320     ClusterBuilder->Verbosity(verbosity);
0321     se->registerSubsystem(ClusterBuilder);
0322   }
0323   else
0324   {
0325     cout << "HCalIn_Clusters - unknown clusterizer setting!" << endl;
0326     exit(1);
0327   }
0328   return;
0329 }
0330
0331 void HCALInner_Eval(const std::string &outputfile)
0332 {
0333   int verbosity = std::max(Enable::VERBOSITY, Enable::HCALIN_VERBOSITY);
0334   Fun4AllServer *se = Fun4AllServer::instance();
0335
0336   CaloEvaluator *eval = new CaloEvaluator("HCALINEVALUATOR", "HCALIN", outputfile);
0337   eval->Verbosity(verbosity);
0338   se->registerSubsystem(eval);
0339
0340   return;
0341 }
0342
0343 void HCALInner_QA()
0344 {
0345   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::HCALIN_VERBOSITY);
0346
0347   Fun4AllServer *se = Fun4AllServer::instance();
0348   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("HCALIN");
0349   qa->Verbosity(verbosity);
0350   se->registerSubsystem(qa);
0351
0352   return;
0353 }
0354
0355 #endif