macros/common/G4_CEmc_Spacal.C

0001 #ifndef MACRO_G4CEMCSPACAL_C
0002 #define MACRO_G4CEMCSPACAL_C
0003
0004 #include <GlobalVariables.C>
0005 #include <QA.C>
0006
0007 #include <phparameter/PHParameterUtils.h>
0008
0009 #include <g4detectors/PHG4CylinderCellReco.h>
0010 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h>
0011 #include <g4detectors/PHG4CylinderSubsystem.h>
0012 #include <g4detectors/PHG4FullProjSpacalCellReco.h>
0013 #include <g4detectors/PHG4SpacalSubsystem.h>
0014
0015 #include <g4calo/RawTowerBuilder.h>
0016 #include <g4calo/RawTowerDigitizer.h>
0017
0018 #include <g4eval/CaloEvaluator.h>
0019
0020 #include <g4main/PHG4Reco.h>
0021 #include <g4main/PHG4Utils.h>
0022
0023 #include <calowaveformsim/CaloWaveformSim.h>
0024
0025 #include <calobase/TowerInfoDefs.h>
0026
0027 #include <caloreco/RawClusterBuilderGraph.h>
0028 #include <caloreco/RawClusterBuilderTemplate.h>
0029 #include <caloreco/RawClusterPositionCorrection.h>
0030 #include <caloreco/RawTowerCalibration.h>
0031 #include <caloreco/CaloTowerBuilder.h>
0032 #include <caloreco/CaloTowerCalib.h>
0033 #include <caloreco/CaloWaveformProcessing.h>
0034 #include <caloreco/CaloTowerStatus.h>
0035
0036
0037 #include <simqa_modules/QAG4SimulationCalorimeter.h>
0038
0039 #include <fun4all/Fun4AllServer.h>
0040
0041 double
0042 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
0043
0044 R__LOAD_LIBRARY(libcalo_reco.so)
0045 R__LOAD_LIBRARY(libg4calo.so)
0046 R__LOAD_LIBRARY(libCaloWaveformSim.so)
0047 R__LOAD_LIBRARY(libg4detectors.so)
0048 R__LOAD_LIBRARY(libg4eval.so)
0049 R__LOAD_LIBRARY(libsimqa_modules.so)
0050
0051 namespace Enable
0052 {
0053   bool CEMC = false;
0054   bool CEMC_ABSORBER = false;
0055   bool CEMC_OVERLAPCHECK = false;
0056   bool CEMC_CELL = false;
0057   bool CEMC_TOWER = false;
0058   bool CEMC_CLUSTER = false;
0059   bool CEMC_EVAL = false;
0060   bool CEMC_QA = false;
0061   bool CEMC_G4Hit = true;
0062   bool CEMC_TOWERINFO = false;
0063   int CEMC_VERBOSITY = 0;
0064 }  // namespace Enable
0065
0066 namespace G4CEMC
0067 {
0068   int Min_cemc_layer = 1;
0069   int Max_cemc_layer = 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   //   2D azimuthal projective SPACAL (slow)
0081   int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal;
0082
0083   enum enu_Cemc_clusterizer
0084   {
0085     kCemcGraphClusterizer,
0086
0087     kCemcTemplateClusterizer
0088   };
0089
0090   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0091   enu_Cemc_clusterizer Cemc_clusterizer = kCemcTemplateClusterizer;
0092   //! graph clusterizer, RawClusterBuilderGraph
0093   // enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
0094   bool useTowerInfoV2 = true;
0095
0096 }  // namespace G4CEMC
0097
0098 // black hole parameters are set in CEmc function
0099 // needs a dummy argument to play with current G4Setup_sPHENIX.C
0100 void CEmcInit(const int i = 0)
0101 {
0102 }
0103
0104 //! EMCal main setup macro
0105 double
0106 CEmc(PHG4Reco *g4Reco, double radius, const int crossings)
0107 {
0108   return CEmc_2DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */ crossings);
0109 }
0110
0111 //! 2D full projective SPACAL
0112 double
0113 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
0114 {
0115   bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
0116   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
0117
0118   double emc_inner_radius = 92;  // emc inner radius from engineering drawing
0119   double cemcthickness = 22.50000 - no_overlapp;
0120
0121   // max radius is 116 cm;
0122   double emc_outer_radius = emc_inner_radius + cemcthickness;  // outer radius
0123   assert(emc_outer_radius < 116);
0124
0125   if (radius > emc_inner_radius)
0126   {
0127     cout << "inconsistency: preshower radius+thickness: " << radius
0128          << " larger than emc inner radius: " << emc_inner_radius << endl;
0129     gSystem->Exit(-1);
0130   }
0131
0132   // the radii are only to determined the thickness of the cemc
0133   radius = emc_inner_radius;
0134
0135   // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
0136   PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0137   cyl->set_double_param("radius", radius);
0138   cyl->set_string_param("material", "G4_TEFLON");
0139   cyl->set_double_param("thickness", 1.5 - no_overlapp);
0140   cyl->SuperDetector("CEMC_ELECTRONICS");
0141   cyl->OverlapCheck(OverlapCheck);
0142   if (AbsorberActive) cyl->SetActive();
0143   g4Reco->registerSubsystem(cyl);
0144
0145   radius += 1.5;
0146   cemcthickness -= 1.5 + no_overlapp;
0147
0148   // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
0149   cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0150   cyl->SuperDetector("CEMC_SPT");
0151   cyl->set_double_param("radius", radius + cemcthickness - 0.5);
0152   cyl->set_string_param("material", "SS310");  // SS310 Stainless Steel
0153   cyl->set_double_param("thickness", 0.5 - no_overlapp);
0154   cyl->OverlapCheck(OverlapCheck);
0155   if (AbsorberActive) cyl->SetActive();
0156   g4Reco->registerSubsystem(cyl);
0157
0158   // this is the z extend and outer radius of the support structure and therefore the z extend
0159   // and radius of the surrounding black holes
0160   double sptlen = PHG4Utils::GetLengthForRapidityCoverage(radius + cemcthickness);
0161   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, sptlen);
0162   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -sptlen);
0163   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius + cemcthickness);
0164
0165   cemcthickness -= 0.5 + no_overlapp;
0166
0167   int ilayer = 0;
0168   PHG4SpacalSubsystem *cemc;
0169
0170   cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0171
0172   cemc->set_int_param("virualize_fiber", 0);
0173   cemc->set_int_param("azimuthal_seg_visible", 1);
0174   cemc->set_int_param("construction_verbose", 0);
0175   cemc->Verbosity(0);
0176   cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0177   cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2023ProjTilted/"));
0178   cemc->set_double_param("radius", radius);            // overwrite minimal radius
0179   cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness
0180   if(G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal) cemc->set_int_param("saveg4hit", Enable::CEMC_G4Hit);
0181
0182   cemc->SetActive();
0183   cemc->SuperDetector("CEMC");
0184   if (AbsorberActive) cemc->SetAbsorberActive();
0185   cemc->OverlapCheck(OverlapCheck);
0186
0187   g4Reco->registerSubsystem(cemc);
0188
0189   if (ilayer > G4CEMC::Max_cemc_layer)
0190   {
0191     cout << "layer discrepancy, current layer " << ilayer
0192          << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
0193   }
0194
0195   radius += cemcthickness;
0196   radius += no_overlapp;
0197
0198   return radius;
0199 }
0200
0201 void CEMC_Cells()
0202 {
0203   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0204
0205   Fun4AllServer *se = Fun4AllServer::instance();
0206
0207   if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0208   {
0209     PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
0210     cemc_cells->Detector("CEMC");
0211     cemc_cells->Verbosity(verbosity);
0212     for (int i = G4CEMC::Min_cemc_layer; i <= G4CEMC::Max_cemc_layer; i++)
0213     {
0214       //          cemc_cells->etaphisize(i, 0.024, 0.024);
0215       const double radius = 95;
0216       cemc_cells->cellsize(i, 2 * M_PI / 256. * radius, 2 * M_PI / 256. * radius);
0217     }
0218     se->registerSubsystem(cemc_cells);
0219   }
0220   else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0221   {
0222     if (!Enable::CEMC_G4Hit) return;
0223     PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0224     cemc_cells->Detector("CEMC");
0225     cemc_cells->Verbosity(verbosity);
0226     cemc_cells->get_light_collection_model().load_data_file(
0227         string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
0228         "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
0229     se->registerSubsystem(cemc_cells);
0230   }
0231   else
0232   {
0233     cout << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0234          << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
0235     gSystem->Exit(-1);
0236     return;
0237   }
0238
0239   return;
0240 }
0241
0242 void CEMC_Towers()
0243 {
0244   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0245
0246   Fun4AllServer *se = Fun4AllServer::instance();
0247
0248   if (!Enable::CEMC_TOWERINFO)
0249   {
0250     if (Enable::CEMC_G4Hit)
0251     {
0252       RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0253       TowerBuilder->Detector("CEMC");
0254       TowerBuilder->set_sim_tower_node_prefix("SIM");
0255       TowerBuilder->Verbosity(verbosity);
0256       se->registerSubsystem(TowerBuilder);
0257     }
0258     double sampling_fraction = 1;
0259     //      sampling_fraction = 0.02244; //from production: /gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.3/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_e-_eta0_8GeV.root
0260     //    sampling_fraction = 2.36081e-02;  //from production: /gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.5/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_e-_eta0_8GeV.root
0261     //    sampling_fraction = 1.90951e-02; // 2017 Tilt porjective SPACAL, 8 GeV photon, eta = 0.3 - 0.4
0262     sampling_fraction = 2e-02;                 // 2017 Tilt porjective SPACAL, tower-by-tower calibration
0263     const double photoelectron_per_GeV = 500;  // 500 photon per total GeV deposition
0264
0265     RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
0266     TowerDigitizer->Detector("CEMC");
0267     TowerDigitizer->Verbosity(verbosity);
0268     TowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
0269     TowerDigitizer->set_variable_pedestal(true);  // read ped central and width from calibrations file comment next 2 lines if true
0270                                                   //   TowerDigitizer->set_pedstal_central_ADC(0);
0271                                                   //   TowerDigitizer->set_pedstal_width_ADC(8);  // eRD1 test beam setting
0272     TowerDigitizer->set_photonelec_ADC(1);        // not simulating ADC discretization error
0273     TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0274     TowerDigitizer->set_variable_zero_suppression(true);                                                          // read zs values from calibrations file comment next line if true
0275                                                                                                                   //   TowerDigitizer->set_zero_suppression_ADC(16);  // eRD1 test beam setting
0276     if (!Enable::CEMC_G4Hit) TowerDigitizer->set_towerinfo(RawTowerDigitizer::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0277     if (Enable::CDB)
0278     {
0279       PHParameterUtils::FillPHParametersFromCDB(TowerDigitizer->GetParameters(), "EMCTOWERCALIB");
0280     }
0281     else
0282     {
0283       TowerDigitizer->GetParameters().ReadFromFile("CEMC", "xml", 0, 0,
0284                                                    string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));  // calibration database
0285     }
0286     se->registerSubsystem(TowerDigitizer);
0287
0288     RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0289     TowerCalibration->Detector("CEMC");
0290     TowerCalibration->set_usetowerinfo_v2(G4CEMC::useTowerInfoV2);
0291     TowerCalibration->Verbosity(verbosity);
0292     if (!Enable::CEMC_G4Hit) TowerCalibration->set_towerinfo(RawTowerCalibration::ProcessTowerType::kTowerInfoOnly);  // just use towerinfo
0293     if (G4CEMC::TowerDigi == RawTowerDigitizer::kNo_digitization)
0294     {
0295       // just use sampling fraction set previously
0296       TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
0297     }
0298     else
0299     {
0300       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kTower_by_tower_calibration);
0301       if (Enable::CDB)
0302       {
0303         PHParameterUtils::FillPHParametersFromCDB(TowerCalibration->GetCalibrationParameters(), "EMCTOWERCALIB");
0304       }
0305       else
0306       {
0307         TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC", "xml", 0, 0,
0308                                                                   string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));  // calibration database
0309       }
0310       TowerCalibration->set_variable_GeV_ADC(true);  // read GeV per ADC from calibrations file comment next line if true
0311       //    TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715);                                                             // overall energy scale based on 4-GeV photon simulations
0312       TowerCalibration->set_variable_pedestal(true);  // read pedestals from calibrations file comment next line if true
0313       //    TowerCalibration->set_pedstal_ADC(0);
0314     }
0315     se->registerSubsystem(TowerCalibration);
0316   }
0317   else
0318   {
0319     CaloWaveformSim *caloWaveformSim = new CaloWaveformSim();
0320     caloWaveformSim->set_detector_type(CaloTowerDefs::CEMC);
0321     caloWaveformSim->set_detector("CEMC");
0322     caloWaveformSim->set_nsamples(12);
0323     caloWaveformSim->set_pedestalsamples(12);
0324     caloWaveformSim->set_timewidth(0.2);
0325     caloWaveformSim->set_peakpos(6);
0326     // pedestal scale down for different beam configurations according to Blair
0327     if(Input::BEAM_CONFIGURATION == Input::pp_ZEROANGLE)
0328     {
0329       caloWaveformSim->set_pedestal_scale(0.69);
0330     }
0331     if(Input::BEAM_CONFIGURATION == Input::pp_COLLISION)
0332     {
0333       caloWaveformSim->set_pedestal_scale(0.77);
0334     }
0335     // caloWaveformSim->Verbosity(2);
0336     // caloWaveformSim->set_noise_type(CaloWaveformSim::NOISE_NONE);
0337     se->registerSubsystem(caloWaveformSim);
0338
0339     CaloTowerBuilder *ca2 = new CaloTowerBuilder();
0340     ca2->set_detector_type(CaloTowerDefs::CEMC);
0341     ca2->set_nsamples(12);
0342     ca2->set_dataflag(false);
0343     ca2->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0344     ca2->set_builder_type(CaloTowerDefs::kWaveformTowerSimv1);
0345     // match our current ZS threshold ~60ADC for emcal
0346     ca2->set_softwarezerosuppression(true, 60);
0347     se->registerSubsystem(ca2);
0348
0349     CaloTowerStatus *statusEMC = new CaloTowerStatus("CEMCSTATUS");
0350     statusEMC->set_detector_type(CaloTowerDefs::CEMC);
0351     statusEMC->set_time_cut(1);
0352     se->registerSubsystem(statusEMC);
0353
0354     CaloTowerCalib *calibEMC = new CaloTowerCalib("CEMCCALIB");
0355     calibEMC->set_detector_type(CaloTowerDefs::CEMC);
0356     calibEMC->set_outputNodePrefix("TOWERINFO_CALIB_");
0357     se->registerSubsystem(calibEMC);
0358   }
0359   return;
0360 }
0361
0362 void CEMC_Clusters()
0363 {
0364   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0365
0366   Fun4AllServer *se = Fun4AllServer::instance();
0367
0368   if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcTemplateClusterizer)
0369   {
0370     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0371     ClusterBuilder->Detector("CEMC");
0372     ClusterBuilder->Verbosity(verbosity);
0373     ClusterBuilder->set_threshold_energy(0.030);  // This threshold should be the same as in CEMCprof_Thresh**.root file below
0374     std::string emc_prof = getenv("CALIBRATIONROOT");
0375     emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0376     ClusterBuilder->LoadProfile(emc_prof);
0377     if (!Enable::CEMC_G4Hit) ClusterBuilder->set_UseTowerInfo(1);  // just use towerinfo
0378     //    ClusterBuilder->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
0379     se->registerSubsystem(ClusterBuilder);
0380   }
0381   else if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcGraphClusterizer)
0382   {
0383     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
0384     ClusterBuilder->Detector("CEMC");
0385     ClusterBuilder->Verbosity(verbosity);
0386     se->registerSubsystem(ClusterBuilder);
0387   }
0388   else
0389   {
0390     cout << "CEMC_Clusters - unknown clusterizer setting!" << endl;
0391     exit(1);
0392   }
0393
0394   RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
0395   if (!Enable::CEMC_G4Hit || Enable::CEMC_TOWERINFO) clusterCorrection->set_UseTowerInfo(1);  // just use towerinfo
0396   //    clusterCorrection->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
0397
0398 //  if (Enable::CDB)
0399 //  {
0400 //    clusterCorrection->Get_eclus_CalibrationParameters().ReadFromCDB("CEMCRECALIB");
0401 //    clusterCorrection->Get_ecore_CalibrationParameters().ReadFromCDB("CEMC_ECORE_RECALIB");
0402 //  }
0403 //  else
0404 //  {
0405 //    clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
0406 //                                                                      // raw location
0407 //                                                                      string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0408 //    clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
0409 //                                                                      // raw location
0410 //                                                                      string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0411 //  }
0412
0413   clusterCorrection->Verbosity(verbosity);
0414   se->registerSubsystem(clusterCorrection);
0415
0416   return;
0417 }
0418 void CEMC_Eval(const std::string &outputfile)
0419 {
0420   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0421
0422   Fun4AllServer *se = Fun4AllServer::instance();
0423
0424   CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile);
0425   eval->Verbosity(verbosity);
0426   eval->set_use_towerinfo(Enable::CEMC_TOWERINFO);
0427   se->registerSubsystem(eval);
0428
0429   return;
0430 }
0431
0432 void CEMC_QA()
0433 {
0434   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::CEMC_VERBOSITY);
0435
0436   Fun4AllServer *se = Fun4AllServer::instance();
0437   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("CEMC");
0438   qa->Verbosity(verbosity);
0439   if(Enable::CEMC_TOWERINFO) qa->set_flags(QAG4SimulationCalorimeter::kProcessTowerinfo);
0440   se->registerSubsystem(qa);
0441
0442   return;
0443 }
0444
0445 #endif