emc_pi0_tbt/towInfoClustering/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 <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_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
0031
0032 R__LOAD_LIBRARY(libcalo_reco.so)
0033 R__LOAD_LIBRARY(libg4calo.so)
0034 R__LOAD_LIBRARY(libg4detectors.so)
0035 R__LOAD_LIBRARY(libg4eval.so)
0036 R__LOAD_LIBRARY(libqa_modules.so)
0037
0038 namespace Enable
0039 {
0040   bool CEMC = false;
0041   bool CEMC_ABSORBER = false;
0042   bool CEMC_OVERLAPCHECK = false;
0043   bool CEMC_CELL = false;
0044   bool CEMC_TOWER = false;
0045   bool CEMC_CLUSTER = false;
0046   bool CEMC_EVAL = false;
0047   bool CEMC_QA = false;
0048   int CEMC_VERBOSITY = 0;
0049 }  // namespace Enable
0050
0051 namespace G4CEMC
0052 {
0053   int Min_cemc_layer = 1;
0054   int Max_cemc_layer = 1;
0055
0056   // Digitization (default photon digi):
0057   RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization;
0058   // directly pass the energy of sim tower to digitized tower
0059   // kNo_digitization
0060   // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
0061   // kSimple_photon_digitization
0062   // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
0063   // kSiPM_photon_digitization
0064
0065   //   2D azimuthal projective SPACAL (slow)
0066   int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal;
0067
0068   enum enu_Cemc_clusterizer
0069   {
0070     kCemcGraphClusterizer,
0071
0072     kCemcTemplateClusterizer
0073   };
0074
0075   //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0076   enu_Cemc_clusterizer Cemc_clusterizer = kCemcTemplateClusterizer;
0077   //! graph clusterizer, RawClusterBuilderGraph
0078   //enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
0079
0080 }  // namespace G4CEMC
0081
0082 // black hole parameters are set in CEmc function
0083 // needs a dummy argument to play with current G4Setup_sPHENIX.C
0084 void CEmcInit(const int i = 0)
0085 {
0086 }
0087
0088 //! EMCal main setup macro
0089 double
0090 CEmc(PHG4Reco *g4Reco, double radius, const int crossings)
0091 {
0092   return CEmc_2DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */ crossings);
0093 }
0094
0095 //! 2D full projective SPACAL
0096 double
0097 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
0098 {
0099   bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
0100   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
0101
0102   double emc_inner_radius = 92;  // emc inner radius from engineering drawing
0103   double cemcthickness = 22.50000 - no_overlapp;
0104
0105   //max radius is 116 cm;
0106   double emc_outer_radius = emc_inner_radius + cemcthickness;  // outer radius
0107   assert(emc_outer_radius < 116);
0108
0109   if (radius > emc_inner_radius)
0110   {
0111     cout << "inconsistency: preshower radius+thickness: " << radius
0112          << " larger than emc inner radius: " << emc_inner_radius << endl;
0113     gSystem->Exit(-1);
0114   }
0115
0116   // the radii are only to determined the thickness of the cemc
0117   radius = emc_inner_radius;
0118
0119   // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
0120   PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0121   cyl->set_double_param("radius", radius);
0122   cyl->set_string_param("material", "G4_TEFLON");
0123   cyl->set_double_param("thickness", 1.5 - no_overlapp);
0124   cyl->SuperDetector("CEMC_ELECTRONICS");
0125   cyl->OverlapCheck(OverlapCheck);
0126   if (AbsorberActive) cyl->SetActive();
0127   g4Reco->registerSubsystem(cyl);
0128
0129   radius += 1.5;
0130   cemcthickness -= 1.5 + no_overlapp;
0131
0132   // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
0133   cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0134   cyl->SuperDetector("CEMC_SPT");
0135   cyl->set_double_param("radius", radius + cemcthickness - 0.5);
0136   cyl->set_string_param("material", "SS310");  // SS310 Stainless Steel
0137   cyl->set_double_param("thickness", 0.5 - no_overlapp);
0138   cyl->OverlapCheck(OverlapCheck);
0139   if (AbsorberActive) cyl->SetActive();
0140   g4Reco->registerSubsystem(cyl);
0141
0142   // this is the z extend and outer radius of the support structure and therefore the z extend
0143   // and radius of the surrounding black holes
0144   double sptlen = PHG4Utils::GetLengthForRapidityCoverage(radius + cemcthickness);
0145   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, sptlen);
0146   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -sptlen);
0147   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius + cemcthickness);
0148
0149   cemcthickness -= 0.5 + no_overlapp;
0150
0151   int ilayer = 0;
0152   PHG4SpacalSubsystem *cemc;
0153
0154   cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0155
0156   cemc->set_int_param("virualize_fiber", 0);
0157   cemc->set_int_param("azimuthal_seg_visible", 1);
0158   cemc->set_int_param("construction_verbose", 0);
0159   cemc->Verbosity(0);
0160   if (Enable::XPLOAD)
0161   {
0162     cemc->UseCDB("CEMC_GEOMETRY");
0163   }
0164   else
0165   {
0166     cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0167     cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2023ProjTilted/"));
0168   }
0169   cemc->set_double_param("radius", radius);            // overwrite minimal radius
0170   cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness
0171
0172   cemc->SetActive();
0173   cemc->SuperDetector("CEMC");
0174   if (AbsorberActive) cemc->SetAbsorberActive();
0175   cemc->OverlapCheck(OverlapCheck);
0176
0177   g4Reco->registerSubsystem(cemc);
0178
0179   if (ilayer > G4CEMC::Max_cemc_layer)
0180   {
0181     cout << "layer discrepancy, current layer " << ilayer
0182          << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
0183   }
0184
0185   radius += cemcthickness;
0186   radius += no_overlapp;
0187
0188   return radius;
0189 }
0190
0191 void CEMC_Cells()
0192 {
0193   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0194
0195   Fun4AllServer *se = Fun4AllServer::instance();
0196
0197   if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0198   {
0199     PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
0200     cemc_cells->Detector("CEMC");
0201     cemc_cells->Verbosity(verbosity);
0202     for (int i = G4CEMC::Min_cemc_layer; i <= G4CEMC::Max_cemc_layer; i++)
0203     {
0204       //          cemc_cells->etaphisize(i, 0.024, 0.024);
0205       const double radius = 95;
0206       cemc_cells->cellsize(i, 2 * M_PI / 256. * radius, 2 * M_PI / 256. * radius);
0207     }
0208     se->registerSubsystem(cemc_cells);
0209   }
0210   else if (G4CEMC::Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0211   {
0212     PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0213     cemc_cells->Detector("CEMC");
0214     cemc_cells->Verbosity(verbosity);
0215     cemc_cells->get_light_collection_model().load_data_file(
0216         string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
0217         "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
0218     se->registerSubsystem(cemc_cells);
0219   }
0220   else
0221   {
0222     cout << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0223          << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
0224     gSystem->Exit(-1);
0225     return;
0226   }
0227
0228   return;
0229 }
0230
0231 void CEMC_Towers()
0232 {
0233   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0234
0235   Fun4AllServer *se = Fun4AllServer::instance();
0236
0237   RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0238   TowerBuilder->Detector("CEMC");
0239   TowerBuilder->set_sim_tower_node_prefix("SIM");
0240   TowerBuilder->Verbosity(verbosity);
0241   se->registerSubsystem(TowerBuilder);
0242
0243   double sampling_fraction = 1;
0244   //      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
0245   //    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
0246   //    sampling_fraction = 1.90951e-02; // 2017 Tilt porjective SPACAL, 8 GeV photon, eta = 0.3 - 0.4
0247   sampling_fraction = 2e-02;                 // 2017 Tilt porjective SPACAL, tower-by-tower calibration
0248   const double photoelectron_per_GeV = 500;  //500 photon per total GeV deposition
0249
0250   RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
0251   TowerDigitizer->Detector("CEMC");
0252   TowerDigitizer->Verbosity(verbosity);
0253   TowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
0254   TowerDigitizer->set_variable_pedestal(true);  //read ped central and width from calibrations file comment next 2 lines if true
0255                                                 //  TowerDigitizer->set_pedstal_central_ADC(0);
0256                                                 //  TowerDigitizer->set_pedstal_width_ADC(8);  // eRD1 test beam setting
0257   TowerDigitizer->set_photonelec_ADC(1);        //not simulating ADC discretization error
0258   TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0259   TowerDigitizer->set_variable_zero_suppression(true);  //read zs values from calibrations file comment next line if true
0260                                                         //  TowerDigitizer->set_zero_suppression_ADC(16);  // eRD1 test beam setting
0261   if (Enable::XPLOAD)
0262   {
0263     TowerDigitizer->GetParameters().ReadFromCDB("EMCTOWERCALIB");
0264   }
0265   else
0266   {
0267     TowerDigitizer->GetParameters().ReadFromFile("CEMC", "xml", 0, 0,
0268                                                  string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));  // calibration database
0269   }
0270   se->registerSubsystem(TowerDigitizer);
0271
0272   RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0273   TowerCalibration->Detector("CEMC");
0274   TowerCalibration->Verbosity(verbosity);
0275
0276   if (G4CEMC::TowerDigi == RawTowerDigitizer::kNo_digitization)
0277   {
0278     // just use sampling fraction set previously
0279     TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
0280   }
0281   else
0282   {
0283     TowerCalibration->set_calib_algorithm(RawTowerCalibration::kTower_by_tower_calibration);
0284     if (Enable::XPLOAD)
0285     {
0286       TowerCalibration->GetCalibrationParameters().ReadFromCDB("EMCTOWERCALIB");
0287     }
0288     else
0289     {
0290       TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC", "xml", 0, 0,
0291                                                                 string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/"));  // calibration database
0292     }
0293     TowerCalibration->set_variable_GeV_ADC(true);  //read GeV per ADC from calibrations file comment next line if true
0294     //    TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715);                                                             // overall energy scale based on 4-GeV photon simulations
0295     TowerCalibration->set_variable_pedestal(true);  //read pedestals from calibrations file comment next line if true
0296     //    TowerCalibration->set_pedstal_ADC(0);
0297   }
0298   se->registerSubsystem(TowerCalibration);
0299
0300   return;
0301 }
0302
0303 void CEMC_Clusters()
0304 {
0305   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0306
0307   Fun4AllServer *se = Fun4AllServer::instance();
0308
0309   if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcTemplateClusterizer)
0310   {
0311     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0312     ClusterBuilder->Detector("CEMC");
0313     ClusterBuilder->Verbosity(verbosity);
0314     ClusterBuilder->set_threshold_energy(0.030);  // This threshold should be the same as in CEMCprof_Thresh**.root file below
0315     std::string emc_prof = getenv("CALIBRATIONROOT");
0316     emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0317     ClusterBuilder->LoadProfile(emc_prof);
0318     //    ClusterBuilder->set_UseTowerInfo(0);
0319     se->registerSubsystem(ClusterBuilder);
0320
0321     RawClusterBuilderTemplate *ClusterBuilder2 = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate2");
0322     ClusterBuilder2->Detector("CEMC");
0323     ClusterBuilder2->Verbosity(verbosity);
0324     ClusterBuilder2->set_threshold_energy(0.030);  // This threshold should be the same as in CEMCprof_Thresh**.root file below
0325     std::string emc_prof2 = getenv("CALIBRATIONROOT");
0326     emc_prof2 += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0327     ClusterBuilder2->LoadProfile(emc_prof2);
0328     ClusterBuilder2->set_UseTowerInfo(1);
0329     se->registerSubsystem(ClusterBuilder2);
0330
0331   }
0332   else if (G4CEMC::Cemc_clusterizer == G4CEMC::kCemcGraphClusterizer)
0333   {
0334     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
0335     ClusterBuilder->Detector("CEMC");
0336     ClusterBuilder->Verbosity(verbosity);
0337     se->registerSubsystem(ClusterBuilder);
0338   }
0339   else
0340   {
0341     cout << "CEMC_Clusters - unknown clusterizer setting!" << endl;
0342     exit(1);
0343   }
0344
0345   RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
0346   if (Enable::XPLOAD)
0347   {
0348     clusterCorrection->Get_eclus_CalibrationParameters().ReadFromCDB("CEMCRECALIB");
0349     clusterCorrection->Get_ecore_CalibrationParameters().ReadFromCDB("CEMC_ECORE_RECALIB");
0350   }
0351   else
0352   {
0353     clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
0354                                                                       //raw location
0355                                                                       string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0356     clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
0357                                                                       //raw location
0358                                                                       string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0359   }
0360
0361   clusterCorrection->Verbosity(verbosity);
0362   se->registerSubsystem(clusterCorrection);
0363
0364
0365   RawClusterPositionCorrection *clusterCorrection2 = new RawClusterPositionCorrection("CEMC");
0366
0367   clusterCorrection2->set_UseTowerInfo(1);
0368   if (Enable::XPLOAD)
0369   {
0370     clusterCorrection2->Get_eclus_CalibrationParameters().ReadFromCDB("CEMCRECALIB");
0371     clusterCorrection2->Get_ecore_CalibrationParameters().ReadFromCDB("CEMC_ECORE_RECALIB");
0372   }
0373   else
0374   {
0375     clusterCorrection2->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
0376                                                                       //raw location
0377                                                                       string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0378     clusterCorrection2->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
0379                                                                       //raw location
0380                                                                       string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
0381   }
0382
0383   clusterCorrection2->Verbosity(verbosity);
0384   se->registerSubsystem(clusterCorrection2);
0385
0386
0387   return;
0388 }
0389 void CEMC_Eval(const std::string &outputfile)
0390 {
0391   int verbosity = std::max(Enable::VERBOSITY, Enable::CEMC_VERBOSITY);
0392
0393   Fun4AllServer *se = Fun4AllServer::instance();
0394
0395   CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile);
0396   eval->Verbosity(verbosity);
0397   se->registerSubsystem(eval);
0398
0399   return;
0400 }
0401
0402 void CEMC_QA()
0403 {
0404   int verbosity = std::max(Enable::QA_VERBOSITY, Enable::CEMC_VERBOSITY);
0405
0406   Fun4AllServer *se = Fun4AllServer::instance();
0407   QAG4SimulationCalorimeter *qa = new QAG4SimulationCalorimeter("CEMC");
0408   qa->Verbosity(verbosity);
0409   se->registerSubsystem(qa);
0410
0411   return;
0412 }
0413
0414 #endif