File indexing completed on 2025-08-05 08:19:42
0001 #ifndef MACRO_G4CEMCSPACAL_C
0002 #define MACRO_G4CEMCSPACAL_C
0003
0004 #include <GlobalVariables.C>
0005 #include <QA.C>
0006
0007 #include <g4detectors/PHG4CylinderCellReco.h>
0008 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h>
0009 #include <g4detectors/PHG4CylinderSubsystem.h>
0010 #include <g4detectors/PHG4FullProjSpacalCellReco.h>
0011 #include <g4detectors/PHG4SpacalSubsystem.h>
0012
0013 #include <g4calo/RawTowerBuilder.h>
0014 #include <g4calo/RawTowerDigitizer.h>
0015
0016 #include <g4eval/CaloEvaluator.h>
0017
0018 #include <g4main/PHG4Reco.h>
0019 #include <g4main/PHG4Utils.h>
0020
0021 #include <caloreco/RawClusterBuilderGraph.h>
0022 #include <caloreco/RawClusterBuilderTemplate.h>
0023 #include <caloreco/RawClusterPositionCorrection.h>
0024 #include <caloreco/RawTowerCalibration.h>
0025 #include <qa_modules/QAG4SimulationCalorimeter.h>
0026
0027 #include <fun4all/Fun4AllServer.h>
0028
0029 double
0030 CEmc_1DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
0031
0032 double
0033 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
0034
0035 R__LOAD_LIBRARY(libcalo_reco.so)
0036 R__LOAD_LIBRARY(libg4calo.so)
0037 R__LOAD_LIBRARY(libg4detectors.so)
0038 R__LOAD_LIBRARY(libg4eval.so)
0039 R__LOAD_LIBRARY(libqa_modules.so)
0040
0041 namespace Enable
0042 {
0043 bool CEMC = false;
0044 bool CEMC_ABSORBER = false;
0045 bool CEMC_OVERLAPCHECK = false;
0046 bool CEMC_CELL = false;
0047 bool CEMC_TOWER = false;
0048 bool CEMC_CLUSTER = false;
0049 bool CEMC_EVAL = false;
0050 bool CEMC_QA = false;
0051 int CEMC_VERBOSITY = 0;
0052 }
0053
0054 namespace G4CEMC
0055 {
0056 int Min_cemc_layer = 1;
0057 int Max_cemc_layer = 1;
0058
0059
0060 RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073 int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal;
0074
0075 enum enu_Cemc_clusterizer
0076 {
0077 kCemcGraphClusterizer,
0078
0079 kCemcTemplateClusterizer
0080 };
0081
0082
0083 enu_Cemc_clusterizer Cemc_clusterizer = kCemcTemplateClusterizer;
0084
0085
0086
0087 }
0088
0089
0090
0091 void CEmcInit(const int i = 0)
0092 {
0093 }
0094
0095
0096 double
0097 CEmc(PHG4Reco *g4Reco, double radius, const int crossings)
0098 {
0099 if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0100 {
0101 return CEmc_1DProjectiveSpacal( g4Reco, radius, crossings);
0102 }
0103 else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0104 {
0105 return CEmc_2DProjectiveSpacal( g4Reco, radius, crossings);
0106 }
0107 else
0108 {
0109 std::cout
0110 << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0111 << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0112 exit(-1);
0113 return 0;
0114 }
0115 }
0116
0117
0118 double
0119 CEmc_1DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
0120 {
0121 bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
0122 bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
0123
0124 double emc_inner_radius = 95.;
0125 double cemcthickness = 12.7;
0126 double emc_outer_radius = emc_inner_radius + cemcthickness;
0127
0128 if (radius > emc_inner_radius)
0129 {
0130 cout << "inconsistency: pstof outer radius: " << radius
0131 << " larger than emc inner radius: " << emc_inner_radius
0132 << endl;
0133 gSystem->Exit(-1);
0134 }
0135
0136
0137 if (radius > emc_inner_radius - 1.5 - no_overlapp)
0138 {
0139 cout << "G4_CEmc_Spacal.C::CEmc() - expect radius < " << emc_inner_radius - 1.5 - no_overlapp << " to install SPACAL" << endl;
0140 exit(1);
0141 }
0142 radius = emc_inner_radius - 1.5 - no_overlapp;
0143
0144
0145 PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0146 cyl->SuperDetector("CEMC_ELECTRONICS");
0147 cyl->set_double_param("radius", radius);
0148 cyl->set_string_param("material", "G4_TEFLON");
0149 cyl->set_double_param("thickness", 1.5);
0150 if (AbsorberActive) cyl->SetActive();
0151 g4Reco->registerSubsystem(cyl);
0152
0153 radius += 1.5;
0154 radius += no_overlapp;
0155
0156 int ilayer = G4CEMC::Min_cemc_layer;
0157 PHG4SpacalSubsystem *cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0158 cemc->set_double_param("radius", emc_inner_radius);
0159 cemc->set_double_param("thickness", cemcthickness);
0160
0161 cemc->SetActive();
0162 cemc->SuperDetector("CEMC");
0163 if (AbsorberActive) cemc->SetAbsorberActive();
0164 cemc->OverlapCheck(OverlapCheck);
0165
0166 g4Reco->registerSubsystem(cemc);
0167
0168 if (ilayer > G4CEMC::Max_cemc_layer)
0169 {
0170 cout << "layer discrepancy, current layer " << ilayer
0171 << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
0172 }
0173
0174 radius += cemcthickness;
0175 radius += no_overlapp;
0176
0177
0178 cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0179 cyl->SuperDetector("CEMC_SPT");
0180 cyl->set_double_param("radius", radius);
0181 cyl->set_string_param("material", "SS310");
0182 cyl->set_double_param("thickness", 0.5);
0183 if (AbsorberActive) cyl->SetActive();
0184 g4Reco->registerSubsystem(cyl);
0185 radius += 0.5;
0186
0187
0188 BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 149.47);
0189 BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -149.47);
0190 BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius);
0191 radius += no_overlapp;
0192
0193 return radius;
0194 }
0195
0196
0197 double
0198 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
0199 {
0200 bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
0201 bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
0202
0203 double emc_inner_radius = 92;
0204 double cemcthickness = 24.00000 - no_overlapp;
0205
0206
0207 double emc_outer_radius = emc_inner_radius + cemcthickness;
0208 assert(emc_outer_radius < 116);
0209
0210 if (radius > emc_inner_radius)
0211 {
0212 cout << "inconsistency: preshower radius+thickness: " << radius
0213 << " larger than emc inner radius: " << emc_inner_radius << endl;
0214 gSystem->Exit(-1);
0215 }
0216
0217
0218 radius = emc_inner_radius;
0219
0220
0221 PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0222 cyl->set_double_param("radius", radius);
0223 cyl->set_string_param("material", "G4_TEFLON");
0224 cyl->set_double_param("thickness", 1.5 - no_overlapp);
0225 cyl->SuperDetector("CEMC_ELECTRONICS");
0226 cyl->OverlapCheck(OverlapCheck);
0227 if (AbsorberActive) cyl->SetActive();
0228 g4Reco->registerSubsystem(cyl);
0229
0230 radius += 1.5;
0231 cemcthickness -= 1.5 + no_overlapp;
0232
0233
0234 cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0235 cyl->SuperDetector("CEMC_SPT");
0236 cyl->set_double_param("radius", radius + cemcthickness - 0.5);
0237 cyl->set_string_param("material", "SS310");
0238 cyl->set_double_param("thickness", 0.5 - no_overlapp);
0239 cyl->OverlapCheck(OverlapCheck);
0240 if (AbsorberActive) cyl->SetActive();
0241 g4Reco->registerSubsystem(cyl);
0242
0243
0244
0245 double sptlen = PHG4Utils::GetLengthForRapidityCoverage(radius + cemcthickness);
0246 BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, sptlen);
0247 BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -sptlen);
0248 BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius + cemcthickness);
0249
0250 cemcthickness -= 0.5 + no_overlapp;
0251
0252 int ilayer = 0;
0253 PHG4SpacalSubsystem *cemc;
0254
0255 cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0256
0257 cemc->set_int_param("virualize_fiber", 0);
0258 cemc->set_int_param("azimuthal_seg_visible", 1);
0259 cemc->set_int_param("construction_verbose", 0);
0260 cemc->Verbosity(0);
0261
0262 cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0263 cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2018ProjTilted/"));
0264 cemc->set_double_param("radius", radius);
0265 cemc->set_double_param("thickness", cemcthickness);
0266
0267 cemc->SetActive();
0268 cemc->SuperDetector("CEMC");
0269 if (AbsorberActive) cemc->SetAbsorberActive();
0270 cemc->OverlapCheck(OverlapCheck);
0271
0272 g4Reco->registerSubsystem(cemc);
0273
0274 if (ilayer > G4CEMC::Max_cemc_layer)
0275 {
0276 cout << "layer discrepancy, current layer " << ilayer
0277 << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
0278 }
0279
0280 radius += cemcthickness;
0281 radius += no_overlapp;
0282
0283 return radius;
0284 }
0285
0286 void CEMC_Cells()
0287 {
0288 int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0289
0290 Fun4AllServer *se = Fun4AllServer::instance();
0291
0292 if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0293 {
0294 PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
0295 cemc_cells->Detector("CEMC");
0296 cemc_cells->Verbosity(verbosity);
0297 for (int i = G4CEMC::Min_cemc_layer; i <= G4CEMC::Max_cemc_layer; i++)
0298 {
0299
0300 const double radius = 95;
0301 cemc_cells->cellsize(i, 2 * M_PI / 256. * radius, 2 * M_PI / 256. * radius);
0302 }
0303 se->registerSubsystem(cemc_cells);
0304 }
0305 else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0306 {
0307 PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0308 cemc_cells->Detector("CEMC");
0309 cemc_cells->Verbosity(verbosity);
0310 cemc_cells->get_light_collection_model().load_data_file(
0311 string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
0312 "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
0313 se->registerSubsystem(cemc_cells);
0314 }
0315 else
0316 {
0317 cout << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0318 << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
0319 gSystem->Exit(-1);
0320 return;
0321 }
0322
0323 return;
0324 }
0325
0326 void CEMC_Towers()
0327 {
0328 int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0329
0330 Fun4AllServer *se = Fun4AllServer::instance();
0331
0332 RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0333 TowerBuilder->Detector("CEMC");
0334 TowerBuilder->set_sim_tower_node_prefix("SIM");
0335 TowerBuilder->Verbosity(verbosity);
0336 se->registerSubsystem(TowerBuilder);
0337
0338 double sampling_fraction = 1;
0339 if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0340 {
0341 sampling_fraction = 0.0234335;
0342 }
0343 else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0344 {
0345
0346
0347
0348 sampling_fraction = 2e-02;
0349 }
0350 else
0351 {
0352 std::cout
0353 << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
0354 << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0355 exit(-1);
0356 return;
0357 }
0358
0359 const double photoelectron_per_GeV = 500;
0360
0361 bool doSimple = true;
0362
0363 RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
0364 TowerDigitizer->Detector("CEMC");
0365 TowerDigitizer->Verbosity(verbosity);
0366 TowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
0367 TowerDigitizer->set_variable_pedestal(true);
0368
0369
0370 TowerDigitizer->set_photonelec_ADC(1);
0371 TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0372 TowerDigitizer->set_variable_zero_suppression(true);
0373
0374 TowerDigitizer->GetParameters().ReadFromFile("CEMC", "xml", 0, 0,
0375 string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));
0376
0377
0378
0379
0380 TowerDigitizer->set_UseConditionsDB(false);
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395 TowerDigitizer->set_DoTowerDecal(true,"emcal_newPatternCinco.root",false);
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421 se->registerSubsystem(TowerDigitizer);
0422
0423
0424 RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0425 TowerCalibration->Detector("CEMC");
0426 TowerCalibration->Verbosity(verbosity);
0427
0428
0429
0430
0431 if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0432 {
0433 if (G4CEMC::TowerDigi == RawTowerDigitizer::kNo_digitization)
0434 {
0435
0436 TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
0437 }
0438 else
0439 {
0440 TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0441 TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV);
0442 TowerCalibration->set_pedstal_ADC(0);
0443 }
0444 }
0445 else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0446 {
0447 if (G4CEMC::TowerDigi == RawTowerDigitizer::kNo_digitization)
0448 {
0449
0450 TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
0451 }
0452 else
0453 {
0454
0455 if (!doSimple)
0456 {
0457
0458
0459
0460 TowerCalibration->set_calib_algorithm(RawTowerCalibration::kTower_by_tower_calibration);
0461 TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC", "xml", 0, 0,
0462 string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));
0463 TowerCalibration->set_variable_GeV_ADC(true);
0464
0465 TowerCalibration->set_variable_pedestal(true);
0466
0467
0468
0469 }
0470 else
0471 {
0472
0473 TowerCalibration->set_calib_algorithm(RawTowerCalibration::kDbfile_tbt_gain_corr);
0474 TowerCalibration->set_UseConditionsDB(false);
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485 TowerCalibration->set_CalibrationFileName("emcal_corr1_00.root");
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495 TowerCalibration->set_calib_const_GeV_ADC(0.87 * 1./ photoelectron_per_GeV / 0.9715);
0496 TowerCalibration->set_pedstal_ADC(0);
0497
0498
0499
0500
0501
0502 }
0503
0504 }
0505 }
0506 else
0507 {
0508 cout << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
0509 << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
0510 gSystem->Exit(-1);
0511 return;
0512 }
0513 se->registerSubsystem(TowerCalibration);
0514
0515 return;
0516 }
0517
0518 void CEMC_Clusters()
0519 {
0520 int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0521
0522 Fun4AllServer *se = Fun4AllServer::instance();
0523
0524 if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcTemplateClusterizer)
0525 {
0526 RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0527 ClusterBuilder->Detector("CEMC");
0528 ClusterBuilder->Verbosity(verbosity);
0529 ClusterBuilder->set_threshold_energy(0.030);
0530 std::string emc_prof = getenv("CALIBRATIONROOT");
0531 emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0532 ClusterBuilder->LoadProfile(emc_prof);
0533 se->registerSubsystem(ClusterBuilder);
0534 }
0535 else if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcGraphClusterizer)
0536 {
0537 RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
0538 ClusterBuilder->Detector("CEMC");
0539 ClusterBuilder->Verbosity(verbosity);
0540 se->registerSubsystem(ClusterBuilder);
0541 }
0542 else
0543 {
0544 cout << "CEMC_Clusters - unknown clusterizer setting!" << endl;
0545 exit(1);
0546 }
0547
0548 RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
0549
0550 clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
0551
0552 string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0553
0554 clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
0555
0556 string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0557
0558 clusterCorrection->Verbosity(verbosity);
0559 se->registerSubsystem(clusterCorrection);
0560
0561 return;
0562 }
0563 void CEMC_Eval(const std::string &outputfile)
0564 {
0565 int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0566
0567 Fun4AllServer *se = Fun4AllServer::instance();
0568
0569 CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile);
0570 eval->Verbosity(verbosity);
0571 se->registerSubsystem(eval);
0572
0573 return;
0574 }
0575
0576 void CEMC_QA()
0577 {
0578 int verbosity = std::max(Enable::QA_VERBOSITY, Enable::CEMC_VERBOSITY);
0579
0580 Fun4AllServer *se = Fun4AllServer::instance();
0581 QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("CEMC");
0582 qa->Verbosity(verbosity);
0583 se->registerSubsystem(qa);
0584
0585 return;
0586 }
0587
0588 #endif