diff --git a/common/G4_AllSilicon.C b/common/G4_AllSilicon.C new file mode 100644 index 00000000..b4b01e95 --- /dev/null +++ b/common/G4_AllSilicon.C @@ -0,0 +1,46 @@ +#ifndef MACRO_G4ALLSILICON_C +#define MACRO_G4ALLSILICON_C + +#include + +#include + +R__LOAD_LIBRARY(libg4lblvtx.so) + +namespace Enable +{ + bool ALLSILICON = false; + bool ALLSILICON_ABSORBER = false; + bool ALLSILICON_OVERLAPCHECK = false; +} // namespace Enable + +namespace G4ALLSILICON +{ + namespace SETTING + { + int geomVersion = 2; + } // namespace SETTING +} // namespace G4FHCAL + +void AllSiliconInit() {} + +void AllSiliconSetup(PHG4Reco *g4Reco) +{ + bool AbsorberActive = Enable::ABSORBER || Enable::ALLSILICON_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::ALLSILICON_OVERLAPCHECK; + AllSiliconTrackerSubsystem *allsili = new AllSiliconTrackerSubsystem(); + + + allsili->set_string_param("GDMPath", string(getenv("CALIBRATIONROOT")) + Form("/AllSiliconTracker/genfitGeom_AllSi_v%d.gdml",G4ALLSILICON::SETTING::geomVersion)); + + allsili->AddAssemblyVolume("VST"); // Barrel + allsili->AddAssemblyVolume("FST"); // Forward disks + allsili->AddAssemblyVolume("BST"); // Backward disks + //allsili->AddAssemblyVolume("BEAMPIPE"); // Beampipe + allsili->SuperDetector("LBLVTX"); + allsili->OverlapCheck(OverlapCheck); + allsili->SetActive(); // this saves hits in the MimosaCore volumes + if (AbsorberActive) allsili->SetAbsorberActive(); // this saves hits in all volumes (in the absorber node) + g4Reco->registerSubsystem(allsili); +} +#endif diff --git a/common/G4_BECAL.C b/common/G4_BECAL.C new file mode 100644 index 00000000..3fc67477 --- /dev/null +++ b/common/G4_BECAL.C @@ -0,0 +1,162 @@ +#ifndef MACRO_G4BECAL_C +#define MACRO_G4BECAL_C + +#include + +#include + +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +R__LOAD_LIBRARY(libcalo_reco.so) +R__LOAD_LIBRARY(libg4calo.so) +R__LOAD_LIBRARY(libg4eiccalos.so) +R__LOAD_LIBRARY(libg4eval.so) + +namespace Enable +{ + bool BECAL = false; + bool BECAL_ABSORBER = false; + bool BECAL_CELL = false; + bool BECAL_TOWER = false; + bool BECAL_CLUSTER = false; + bool BECAL_EVAL = false; + bool BECAL_OVERLAPCHECK = false; + int BECAL_VERBOSITY = 0; +} // namespace Enable + + + +namespace G4BECAL +{ + + //double minz = -273.6*cm; + //double maxz = 142.4*cm; + double minz = -453; + double maxz = 371; + double topradius = 138; + double radius = 85; + + // this is default set to -1.5set_string_param("mapping_file", mapping_becal.str()); + becal->OverlapCheck(OverlapCheck); + becal->SetActive(); + becal->SuperDetector("BECAL"); + if (AbsorberActive) becal->SetAbsorberActive(); + + g4Reco->registerSubsystem(becal); + + return G4BECAL::topradius; + +} + +void BECAL_Cells(int verbosity = 0) +{ + return; +} + +void BECAL_Towers() +{ + + int verbosity = std::max(Enable::VERBOSITY, Enable::BECAL_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + ostringstream mapping_BECAL; + mapping_BECAL << getenv("CALIBRATIONROOT") << "/BarrelEcal/mapping/towerMap_BEMC_v001.txt"; + + const double photoelectron_per_GeV = 5000; + + RawTowerBuilderByHitIndexBECAL *tower_BECAL = new RawTowerBuilderByHitIndexBECAL("TowerBuilder_BECAL"); + tower_BECAL->Detector("BECAL"); + tower_BECAL->set_sim_tower_node_prefix("SIM"); + tower_BECAL->EminCut(1e-7); + tower_BECAL->GeometryTableFile(mapping_BECAL.str()); + tower_BECAL->Verbosity(verbosity); + se->registerSubsystem(tower_BECAL); + + RawTowerDigitizer *TowerDigitizer_BECAL = new RawTowerDigitizer("BECALRawTowerDigitizer"); + TowerDigitizer_BECAL->Detector("BECAL"); + TowerDigitizer_BECAL->Verbosity(verbosity); +// TowerDigitizer_BECAL->Verbosity(verbosity); + TowerDigitizer_BECAL->set_digi_algorithm(G4BECAL::TowerDigi); + TowerDigitizer_BECAL->set_raw_tower_node_prefix("RAW"); + TowerDigitizer_BECAL->set_pedstal_central_ADC(0); + TowerDigitizer_BECAL->set_pedstal_width_ADC(0); // eRD1 test beam setting + TowerDigitizer_BECAL->set_photonelec_ADC(1); // not simulating ADC discretization error + TowerDigitizer_BECAL->set_photonelec_yield_visible_GeV(photoelectron_per_GeV); + TowerDigitizer_BECAL->set_zero_suppression_ADC(0); // eRD1 test beam setting + se->registerSubsystem(TowerDigitizer_BECAL); + + RawTowerCalibration *TowerCalibration_BECAL = new RawTowerCalibration("BECALRawTowerCalibration"); + TowerCalibration_BECAL->Detector("BECAL"); + TowerCalibration_BECAL->Verbosity(verbosity); + TowerCalibration_BECAL->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + TowerCalibration_BECAL->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV); + TowerCalibration_BECAL->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration_BECAL); + +} + +void BECAL_Clusters() +{ + + return; +} + +void BECAL_Eval(const std::string &outputfile) +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::BECAL_VERBOSITY); + Fun4AllServer *se = Fun4AllServer::instance(); + + CaloEvaluator *eval = new CaloEvaluator("BECALEVALUATOR", "BECAL", outputfile.c_str()); + eval->set_do_cluster_eval(false); + eval->Verbosity(1); + se->registerSubsystem(eval); + + return; +} +#endif diff --git a/common/G4_Barrel_EIC.C b/common/G4_Barrel_EIC.C index 20b4bc72..cb9132ea 100644 --- a/common/G4_Barrel_EIC.C +++ b/common/G4_Barrel_EIC.C @@ -21,69 +21,82 @@ namespace Enable void BarrelInit() { } + +void BarrelFastKalmanFilterConfig(PHG4TrackFastSim * kalman_filter) +{ + + // import Kalman filter config (lines 226 to 246 here: https://github.com/eic/g4lblvtx/blob/master/macros/auxiliary_studies/simplified_geometry/Fun4All_G4_simplified_v2.C): + + // add Vertexing Layers + kalman_filter->add_phg4hits( + "G4HIT_SVTX", // const std::string& phg4hitsNames, + PHG4TrackFastSim::Cylinder, + 999., // radial-resolution [cm] + 10. / 10000. / sqrt(12.), // azimuthal-resolution [cm] + 10. / 10000. / sqrt(12.), // z-resolution [cm] + 1, // efficiency, + 0 // noise hits + ); + + // add Barrel Layers + kalman_filter->add_phg4hits( + "G4HIT_BARR", // const std::string& phg4hitsNames, + PHG4TrackFastSim::Cylinder, + 999., // radial-resolution [cm] + 10. / 10000. / sqrt(12.), // azimuthal-resolution [cm] + 10. / 10000. / sqrt(12.), // z-resolution [cm] + 1, // efficiency, + 0 // noise hits + ); + +} + + //-----------------------------------------------------------------------------------// void Barrel(PHG4Reco *g4Reco, int det_ver = 3) { - // import Geometry (lines 111 to 148 in https://github.com/eic/g4lblvtx/blob/master/macros/auxiliary_studies/simplified_geometry/Fun4All_G4_simplified_v2.C): - PHG4CylinderSubsystem * cyl(nullptr); + PHG4CylinderSubsystem *cyl(nullptr); //--------------------------- // Vertexing - double si_vtx_r_pos[] = {3.30,5.70}; - const int nVtxLayers = sizeof(si_vtx_r_pos)/sizeof(*si_vtx_r_pos); - for (int ilayer = 0; ilayer < nVtxLayers ; ilayer++){ - cyl = new PHG4CylinderSubsystem("SVTX", ilayer); - cyl->set_string_param("material" , "G4_Si" ); - cyl->set_double_param("radius" , si_vtx_r_pos[ilayer] ); - cyl->set_double_param("thickness", 0.05/100.*9.37 ); - cyl->set_double_param("place_z" , 0 ); - cyl->set_double_param("length" , 30.); - cyl->SetActive(); - cyl->SuperDetector("SVTX"); - g4Reco->registerSubsystem(cyl); + double si_vtx_r_pos[] = {3.30, 5.70, 8.10}; + const int nVtxLayers = sizeof(si_vtx_r_pos) / sizeof(*si_vtx_r_pos); + for (int ilayer = 0; ilayer < nVtxLayers; ilayer++) + { + cyl = new PHG4CylinderSubsystem("SVTX", ilayer); + cyl->set_string_param("material", "G4_Si"); + cyl->set_double_param("radius", si_vtx_r_pos[ilayer]); + cyl->set_double_param("thickness", 0.05 / 100. * 9.37); + cyl->set_double_param("place_z", 0); + cyl->set_double_param("length", 30.); + cyl->SetActive(); + cyl->SuperDetector("SVTX"); + g4Reco->registerSubsystem(cyl); } //--------------------------- // Barrel - double si_r_pos[] = {21.,22.68,39.3,43.23}; - const int nTrckLayers = sizeof(si_r_pos)/sizeof(*si_r_pos); - double si_z_length[] = {54.,60.,105.,114.}; - for (int ilayer = 0; ilayer < nTrckLayers ; ilayer++){ - cyl = new PHG4CylinderSubsystem("BARR", ilayer); - cyl->set_string_param("material" , "G4_Si" ); - cyl->set_double_param("radius" , si_r_pos[ilayer] ); - cyl->set_double_param("thickness", 0.55/100.*9.37 ); - cyl->set_double_param("place_z" , 0 ); - cyl->set_double_param("length" , si_z_length[ilayer]); - cyl->SetActive(); - cyl->SuperDetector("BARR"); - g4Reco->registerSubsystem(cyl); + double si_r_pos[] = {21., 22.68}; + const int nTrckLayers = sizeof(si_r_pos) / sizeof(*si_r_pos); + double si_z_length[] = {60., 60.}; + for (int ilayer = 0; ilayer < nTrckLayers; ilayer++) + { + cyl = new PHG4CylinderSubsystem("BARR", ilayer); + cyl->set_string_param("material", "G4_Si"); + cyl->set_double_param("radius", si_r_pos[ilayer]); + cyl->set_double_param("thickness", 0.05 / 100. * 9.37); + cyl->set_double_param("place_z", 0); + cyl->set_double_param("length", si_z_length[ilayer]); + cyl->SetActive(); + cyl->SuperDetector("BARR"); + g4Reco->registerSubsystem(cyl); } - // import Kalman filter config (lines 226 to 246 here: https://github.com/eic/g4lblvtx/blob/master/macros/auxiliary_studies/simplified_geometry/Fun4All_G4_simplified_v2.C): - - // add Vertexing Layers - TRACKING::FastKalmanFilter->add_phg4hits( - "G4HIT_SVTX", // const std::string& phg4hitsNames, - PHG4TrackFastSim::Cylinder, - 999., // radial-resolution [cm] - 10./10000./sqrt(12.), // azimuthal-resolution [cm] - 10./10000./sqrt(12.), // z-resolution [cm] - 1, // efficiency, - 0 // noise hits - ); + BarrelFastKalmanFilterConfig(TRACKING::FastKalmanFilter); - // add Barrel Layers - TRACKING::FastKalmanFilter->add_phg4hits( - "G4HIT_BARR", // const std::string& phg4hitsNames, - PHG4TrackFastSim::Cylinder, - 999., // radial-resolution [cm] - 10./10000./sqrt(12.), // azimuthal-resolution [cm] - 10./10000./sqrt(12.), // z-resolution [cm] - 1, // efficiency, - 0 // noise hits - ); + BarrelFastKalmanFilterConfig(TRACKING::FastKalmanFilterInnerTrack); + BarrelFastKalmanFilterConfig(TRACKING::FastKalmanFilterSiliconTrack); } #endif diff --git a/common/G4_DIRC.C b/common/G4_DIRC.C index ff236bb5..d0eb503e 100644 --- a/common/G4_DIRC.C +++ b/common/G4_DIRC.C @@ -1,139 +1,77 @@ +/*! + * \file G4_DIRC.C + * \brief Setup the gas DIRC detector full sim + * \author Cameron Dean + * \version $Revision: 0.1 $ + * \date $Date: 2021/08/17 $ + */ #ifndef MACRO_G4DIRC_C #define MACRO_G4DIRC_C #include -#include -#include -#include - +#include #include #include -#include - R__LOAD_LIBRARY(libg4detectors.so) - -/*! - * \file G4_DIRC.C - * \brief Macro setting up the barrel DIRC - * \author Jin Huang - * \version $Revision: 1.3 $ - * \date $Date: 2013/10/09 01:08:17 $ - */ +R__LOAD_LIBRARY(libg4eicdirc.so) namespace Enable { bool DIRC = false; bool DIRC_OVERLAPCHECK = false; + int DIRC_VERBOSITY = 0; + double DIRC_SCALE = 10; //DIRC class is in mm, ECCE is in cm } // namespace Enable -namespace G4DIRC -{ - // Grzegorz Kalicy - // -position in z around IP -168 cm to 287 cm - // -69cm radius for the barrel (70cm inner radius for the bars) - // -12 bar boxes, 10 long bars side-by-side in a bar box - // -bar length 425cm - // -Solid fused silica prism: 24 x 36 x 30 cm3 (H x W x L) - // -Additional longitudinal space for MCP-PMTs, readout cards, cables: ~13cm - // -radial thickness 7-8 cm, including mechanical support - // -~16-18% of a radiation length at normal incidence - - double radiator_R = 70; - double z_prism = 30; - double z_end = +168; - double z_start = -287 + z_prism; - double length = z_end - z_start; - double z_shift = 0.5 * (z_end + z_start); - double outer_skin_radius = 78; -} // namespace G4DIRC - void DIRCInit() { - BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4DIRC::outer_skin_radius); - BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4DIRC::z_end); - BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, G4DIRC::z_start - G4DIRC::z_prism); + BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 210.); + BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 280.); } -//! Babar DIRC (Without most of support structure) -//! Ref: I. Adam et al. The DIRC particle identification system for the BaBar experiment. -//! Nucl. Instrum. Meth., A538:281-357, 2005. doi:10.1016/j.nima.2004.08.129. -double DIRCSetup(PHG4Reco *g4Reco) +void DIRCSetup(PHG4Reco* g4Reco) { bool OverlapCheck = Enable::OVERLAPCHECK || Enable::DIRC_OVERLAPCHECK; + int verbosity = std::max(Enable::VERBOSITY, Enable::DIRC_VERBOSITY); + + G4EicDircSubsystem *dircSubsys = new G4EicDircSubsystem("hpDIRC"); + dircSubsys->SuperDetector("hpDIRC"); + dircSubsys->set_double_param("rMin", 74.1 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("rMin_inner", 60.0 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("length", (287 + 168) * Enable::DIRC_SCALE); + dircSubsys->set_double_param("NBars", 11); + dircSubsys->set_double_param("Radius", 75.0 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Prizm_width", 38.65 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Prizm_length", 30.0 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Prizm_height_at_lens", 3.7 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Bar_thickness", 1.7 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Bar_width", 3.5 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("BarL_length", 122.5 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("BarS_length", 56.0 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("Mirror_height", 2.0 * Enable::DIRC_SCALE); + dircSubsys->set_double_param("z_shift", -45 * Enable::DIRC_SCALE); + dircSubsys->set_int_param("Geom_type", 0); // 0-whole DIRC, 1-one bar box + dircSubsys->set_int_param("Lens_id", 3); // 3- 3-layer spherical lens + dircSubsys->set_int_param("MCP_rows", 6); + dircSubsys->set_int_param("MCP_columns", 4); + dircSubsys->set_int_param("NBoxes",12); + dircSubsys->set_int_param("Bar_pieces", 4); + + dircSubsys->OverlapCheck(OverlapCheck); + dircSubsys->Verbosity(verbosity); + dircSubsys->SetActive(); + + g4Reco->registerSubsystem(dircSubsys); - PHG4SectorSubsystem *dirc; - dirc = new PHG4SectorSubsystem("DIRC"); - dirc->get_geometry().set_normal_polar_angle(M_PI / 2); - dirc->get_geometry().set_normal_start(G4DIRC::radiator_R * PHG4Sector::Sector_Geometry::Unit_cm()); - dirc->get_geometry().set_min_polar_angle(atan2(G4DIRC::radiator_R, G4DIRC::z_end)); - dirc->get_geometry().set_max_polar_angle(atan2(G4DIRC::radiator_R, G4DIRC::z_start)); - dirc->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge()); - dirc->get_geometry().set_max_polar_edge(PHG4Sector::Sector_Geometry::FlatEdge()); - dirc->get_geometry().set_material("Quartz"); - dirc->get_geometry().set_N_Sector(12); - dirc->OverlapCheck(OverlapCheck); - dirc->get_geometry().AddLayer("Radiator", "Quartz", 1.7 * PHG4Sector::Sector_Geometry::Unit_cm(), true); - g4Reco->registerSubsystem(dirc); - - PHG4CylinderSubsystem *cyl; - - // The cylinder skins provide most of the strength - // and stiffness of the CST. The thickness of the inner - // and outer skins is 1.27 and 0.76 mm, respectively - - // Inner skin: - cyl = new PHG4CylinderSubsystem("DIRC_CST_Inner_Skin", 10); - cyl->set_double_param("radius", 69); - cyl->set_double_param("length", G4DIRC::length + G4DIRC::z_prism); - cyl->set_string_param("material", "G4_Al"); - cyl->set_double_param("thickness", 0.127); - cyl->set_double_param("place_x", 0.); - cyl->set_double_param("place_y", 0.); - cyl->set_double_param("place_z", G4DIRC::z_shift - G4DIRC::z_prism * 0.5); - cyl->SetActive(0); - cyl->SuperDetector("DIRC"); - cyl->OverlapCheck(OverlapCheck); - - g4Reco->registerSubsystem(cyl); - - // Outer skin: - cyl = new PHG4CylinderSubsystem("DIRC_CST_Outer_Skin", 11); - cyl->set_double_param("radius", G4DIRC::outer_skin_radius - 0.076); - cyl->set_double_param("length", G4DIRC::length); - cyl->set_string_param("material", "G4_Al"); - cyl->set_double_param("thickness", 0.076); - cyl->set_double_param("place_x", 0.); - cyl->set_double_param("place_y", 0.); - cyl->set_double_param("place_z", G4DIRC::z_shift); - cyl->SetActive(0); - cyl->SuperDetector("DIRC"); - cyl->OverlapCheck(OverlapCheck); - - g4Reco->registerSubsystem(cyl); - - // simple approximation for DIRC prism - PHG4ConeSubsystem *cone = new PHG4ConeSubsystem("DIRC_Prism"); - cone->set_color(0, 1, 0); - cone->SetR1(G4DIRC::radiator_R, G4DIRC::radiator_R + 20); - cone->SetR2(G4DIRC::radiator_R, G4DIRC::radiator_R + 2); - cone->SetZlength(0.5 * G4DIRC::z_prism); - cone->SetPlaceZ(G4DIRC::z_start - 0.5 * G4DIRC::z_prism); - cone->SetMaterial("Quartz"); - cone->SetActive(0); - cone->SuperDetector("DIRC"); - cone->OverlapCheck(OverlapCheck); - g4Reco->registerSubsystem(cone); - - // track projection to DIRC reference radiator R if (TRACKING::FastKalmanFilter) { - TRACKING::FastKalmanFilter->add_cylinder_state("DIRC", G4DIRC::radiator_R); + // project to an reference plane at z=170 cm + TRACKING::FastKalmanFilter-> add_zplane_state("DIRC", 185); TRACKING::ProjectionNames.insert("DIRC"); } - // Done - return G4DIRC::outer_skin_radius; + } #endif diff --git a/common/G4_DRCALO.C b/common/G4_DRCALO.C new file mode 100644 index 00000000..0c2e391d --- /dev/null +++ b/common/G4_DRCALO.C @@ -0,0 +1,234 @@ +#ifndef MACRO_G4DR_C +#define MACRO_G4DR_C + +#include + +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +R__LOAD_LIBRARY(libcalo_reco.so) +R__LOAD_LIBRARY(libg4calo.so) +R__LOAD_LIBRARY(libg4detectors.so) +R__LOAD_LIBRARY(libg4eval.so) + +namespace Enable +{ + bool DRCALO = false; + bool DRCALO_ABSORBER = false; + bool DRCALO_CELL = false; + bool DRCALO_TOWER = false; + bool DRCALO_CLUSTER = false; + bool DRCALO_EVAL = false; + bool DRCALO_OVERLAPCHECK = false; + int DRCALO_VERBOSITY = 0; +} // namespace Enable + +namespace G4DRCALO +{ + // from ForwardDualReadout/mapping/towerMap_DRCALO_v005.txt + double Gz0 = 400.; + double Gdz = 100.; + double outer_radius = 262.; + enum enu_FHcal_clusterizer + { + kFHcalGraphClusterizer, + kFHcalTemplateClusterizer + }; + //template clusterizer, as developed by Sasha Bazilevsky + enu_FHcal_clusterizer FHcal_clusterizer = kFHcalTemplateClusterizer; + // graph clusterizer + //enu_FHcal_clusterizer FHcal_clusterizer = kFHcalGraphClusterizer; + namespace SETTING + { + bool Tungsten = false; + bool Quartz = false; + bool PMMA = false; + bool FwdConfig = false; + bool FwdSquare = false; + } // namespace SETTING +} // namespace G4DRCALO + +void DRCALOInit() +{ + // simple way to check if only 1 of the settings is true + if (( (G4DRCALO::SETTING::Tungsten ? 1 : 0) + (G4DRCALO::SETTING::Quartz ? 1 : 0) + (G4DRCALO::SETTING::PMMA ? 1 : 0)) > 1) + { + cout << "use only G4DRCALO::SETTING::Tungsten=true or G4DRCALO::SETTING::Quartz=true or G4DRCALO::SETTING::PMMA=true" << endl; + gSystem->Exit(1); + } + if (( (G4DRCALO::SETTING::FwdSquare ? 1 : 0) + (G4DRCALO::SETTING::FwdConfig ? 1 : 0)) > 1) + { + cout << "use only G4DRCALO::SETTING::FwdSquare=true or G4DRCALO::SETTING::FwdConfig=true" << endl; + gSystem->Exit(1); + } + + BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4DRCALO::outer_radius); + BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4DRCALO::Gz0 + G4DRCALO::Gdz / 2.); +} + +void DRCALOSetup(PHG4Reco *g4Reco) +{ + const bool AbsorberActive = Enable::ABSORBER || Enable::DRCALO_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::DRCALO_OVERLAPCHECK; + Fun4AllServer *se = Fun4AllServer::instance(); + + /** Use dedicated DRCALO module */ + PHG4ForwardDualReadoutSubsystem *hhcal = new PHG4ForwardDualReadoutSubsystem("DRCALO"); + hhcal->OverlapCheck(OverlapCheck); + hhcal->SuperDetector("DRCALO"); + hhcal->SetActive(); + + ostringstream mapping_drcalo; + + // Switch to desired calo setup + if(G4DRCALO::SETTING::Tungsten) { // design with PMMA cerenkov fibers and tungsten absorber material + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_tungsten_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_tungsten.txt"; + } else if(G4DRCALO::SETTING::Quartz) { // design with Quartz cerenkov fibers + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_Quartz_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_Quartz.txt"; + } else if(G4DRCALO::SETTING::PMMA) { // design with PMMA cerenkov fibers + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_PMMA_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_PMMA.txt"; + } else { // NOTCHED design with PMMA cerenkov fiber + if (G4DRCALO::SETTING::FwdConfig) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default_FwdConfig.txt"; + else if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default.txt"; + } + + hhcal->SetTowerMappingFile(mapping_drcalo.str()); + + + if (AbsorberActive) hhcal->SetAbsorberActive(); + + g4Reco->registerSubsystem(hhcal); +} + +void DRCALO_Cells(int verbosity = 0) +{ + return; +} + +void DRCALO_Towers() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::DRCALO_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + ostringstream mapping_drcalo; + + // Switch to desired calo setup + if(G4DRCALO::SETTING::Tungsten) { // design with PMMA cerenkov fibers and tungsten absorber material + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_tungsten_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_tungsten.txt"; + } else if(G4DRCALO::SETTING::Quartz) { // design with Quartz cerenkov fibers + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_Quartz_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_Quartz.txt"; + } else if(G4DRCALO::SETTING::PMMA) { // design with PMMA cerenkov fibers + if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_PMMA_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_PMMA.txt"; + } else { // NOTCHED design with PMMA cerenkov fiber + if (G4DRCALO::SETTING::FwdConfig) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default_FwdConfig.txt"; + else if (G4DRCALO::SETTING::FwdSquare) + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default_FwdSquare.txt"; + else + mapping_drcalo << getenv("CALIBRATIONROOT") << "/DRCALO/mapping/towerMap_DRCALO_default.txt"; + } + cout << "using " << mapping_drcalo.str() << endl; + + + RawTowerBuilderDRCALO *tower_DRCALO = new RawTowerBuilderDRCALO("TowerBuilder_DRCALO"); + tower_DRCALO->Detector("DRCALO"); + tower_DRCALO->set_sim_tower_node_prefix("SIM"); + tower_DRCALO->GeometryTableFile(mapping_drcalo.str()); + + se->registerSubsystem(tower_DRCALO); + + cout << "def: using default for DRCALO towers" << endl; + RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("DRCALORawTowerDigitizer"); + TowerDigitizer->Detector("DRCALO"); + TowerDigitizer->Verbosity(1); + TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kNo_digitization); + se->registerSubsystem(TowerDigitizer); + + RawTowerCalibration *TowerCalibration = new RawTowerCalibration("DRCALORawTowerCalibration"); + TowerCalibration->Detector("DRCALO"); + TowerCalibration->Verbosity(verbosity); + TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + TowerCalibration->set_calib_const_GeV_ADC(1. / 0.03898); // calibrated with muons + TowerCalibration->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration); +} + +void DRCALO_Clusters() +{ + // int verbosity = std::max(Enable::VERBOSITY, Enable::DRCALO_VERBOSITY); + // Fun4AllServer *se = Fun4AllServer::instance(); + + // if (G4DRCALO::FHcal_clusterizer == G4DRCALO::kFHcalTemplateClusterizer) + // { + // RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("DRCALORawClusterBuilderTemplate"); + // ClusterBuilder->Detector("DRCALO"); + // ClusterBuilder->SetPlanarGeometry(); // has to be called after Detector() + // ClusterBuilder->Verbosity(verbosity); + // ClusterBuilder->set_threshold_energy(0.100); + // se->registerSubsystem(ClusterBuilder); + // } + // else if (G4DRCALO::FHcal_clusterizer == G4DRCALO::kFHcalTemplateClusterizer) + // { + // RawClusterBuilderFwd *ClusterBuilder = new RawClusterBuilderFwd("DRCALORawClusterBuilderFwd"); + // ClusterBuilder->Detector("DRCALO"); + // ClusterBuilder->Verbosity(verbosity); + // ClusterBuilder->set_threshold_energy(0.100); + // se->registerSubsystem(ClusterBuilder); + // } + // else + // { + // cout << "DRCALO_Clusters - unknown clusterizer setting " << G4DRCALO::FHcal_clusterizer << endl; + // gSystem->Exit(1); + // } + + return; +} + +void DRCALO_Eval(const std::string &outputfile) +{ + // int verbosity = std::max(Enable::VERBOSITY, Enable::DRCALO_VERBOSITY); + // Fun4AllServer *se = Fun4AllServer::instance(); + + // CaloEvaluator *eval = new CaloEvaluator("DRCALOEVALUATOR", "DRCALO", outputfile.c_str()); + // eval->Verbosity(verbosity); + // se->registerSubsystem(eval); + + return; +} +#endif diff --git a/common/G4_DSTReader_EICDetector.C b/common/G4_DSTReader_EICDetector.C index 19f0bed0..528549c5 100644 --- a/common/G4_DSTReader_EICDetector.C +++ b/common/G4_DSTReader_EICDetector.C @@ -67,81 +67,111 @@ void G4DSTreader_EICDetector(const string &outputFile = "G4sPHENIXCells.root") if (Enable::EGEM) { ana->AddNode("EGEM_0"); - ana->AddNode("EGEM_1"); } if (Enable::FGEM) { ana->AddNode("FGEM_0"); - ana->AddNode("FGEM_1"); } - if (Enable::FST) - { - ana->AddNode("FST_0"); - ana->AddNode("FST_1"); - ana->AddNode("FST_2"); - ana->AddNode("FST_3"); - ana->AddNode("FST_4"); - ana->AddNode("FST_5"); - } - - if (Enable::CEMC) + if (Enable::FTTL) { - ana->AddNode("CEMC"); - if (Enable::ABSORBER || Enable::CEMC_ABSORBER) - { - ana->AddNode("ABSORBER_CEMC"); - ana->AddNode("CEMC_ELECTRONICS_0"); - } + ana->AddNode("FTTL_0"); + ana->AddNode("FTTL_1"); } - - if (Enable::HCALIN) + if (Enable::ETTL) { - ana->AddNode("HCALIN"); - if (Enable::ABSORBER || Enable::HCALIN_ABSORBER) - { - ana->AddNode("ABSORBER_HCALIN"); - ana->AddNode("HCALIN_SPT"); - } + ana->AddNode("ETTL_0"); + ana->AddNode("ETTL_1"); } - - if (Enable::MAGNET) + if (Enable::CTTL) { - if (Enable::ABSORBER || Enable::MAGNET_ABSORBER) - ana->AddNode("MAGNET"); + ana->AddNode("CTTL_0"); } - - if (Enable::HCALOUT) + if (Enable::RWELL) { - ana->AddNode("HCALOUT"); - if (Enable::ABSORBER || Enable::HCALOUT_ABSORBER) - ana->AddNode("ABSORBER_HCALOUT"); + ana->AddNode("RWELL_0"); + ana->AddNode("RWELL_1"); + ana->AddNode("RWELL_2"); } - - if (Enable::FHCAL) - { - ana->AddNode("FHCAL"); - if (Enable::ABSORBER || Enable::FHCAL_ABSORBER) - ana->AddNode("ABSORBER_FHCAL"); - } - - if (Enable::FEMC) + if (Enable::BARREL) { - ana->AddNode("FEMC"); - if (Enable::ABSORBER || Enable::FEMC_ABSORBER) - ana->AddNode("ABSORBER_FEMC"); + ana->AddNode("SVTX"); + ana->AddNode("BARR"); } - - if (Enable::EEMC) + if (Enable::FST) { - ana->AddNode("EEMC"); - } + ana->AddNode("FST_0"); + ana->AddNode("FST_1"); + ana->AddNode("FST_2"); + ana->AddNode("FST_3"); + ana->AddNode("FST_4"); + ana->AddNode("EFST_0"); + ana->AddNode("EFST_1"); + ana->AddNode("EFST_2"); + ana->AddNode("EFST_3"); + } +// +// if (Enable::CEMC) +// { +// ana->AddNode("CEMC"); +// if (Enable::ABSORBER || Enable::CEMC_ABSORBER) +// { +// ana->AddNode("ABSORBER_CEMC"); +// ana->AddNode("CEMC_ELECTRONICS_0"); +// } +// } +// +// if (Enable::HCALIN) +// { +// ana->AddNode("HCALIN"); +// if (Enable::ABSORBER || Enable::HCALIN_ABSORBER) +// { +// ana->AddNode("ABSORBER_HCALIN"); +// ana->AddNode("HCALIN_SPT"); +// } +// } +// +// if (Enable::MAGNET) +// { +// if (Enable::ABSORBER || Enable::MAGNET_ABSORBER) +// ana->AddNode("MAGNET"); +// } +// +// if (Enable::HCALOUT) +// { +// ana->AddNode("HCALOUT"); +// if (Enable::ABSORBER || Enable::HCALOUT_ABSORBER) +// ana->AddNode("ABSORBER_HCALOUT"); +// } +// +// if (Enable::FHCAL) +// { +// ana->AddNode("FHCAL"); +// if (Enable::ABSORBER || Enable::FHCAL_ABSORBER) +// ana->AddNode("ABSORBER_FHCAL"); +// } +// +// if (Enable::FEMC) +// { +// ana->AddNode("FEMC"); +// if (Enable::ABSORBER || Enable::FEMC_ABSORBER) +// ana->AddNode("ABSORBER_FEMC"); +// } +// +// if (Enable::EEMC) +// { +// ana->AddNode("EEMC"); +// } if (Enable::DIRC) { ana->AddNode("DIRC"); } if (Enable::RICH) { - ana->AddNode("RICH"); + ana->AddNode("dRICh_0"); + } + if (Enable::mRICH) + { + ana->AddNode("mRICH"); } if (Enable::BLACKHOLE) diff --git a/common/G4_EEMC_hybrid.C b/common/G4_EEMC_hybrid.C new file mode 100644 index 00000000..ba1b4793 --- /dev/null +++ b/common/G4_EEMC_hybrid.C @@ -0,0 +1,319 @@ +#ifndef MACRO_G4EEMCHYBRID_C +#define MACRO_G4EEMCHYBRID_C + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +R__LOAD_LIBRARY(libcalo_reco.so) +R__LOAD_LIBRARY(libg4eiccalos.so) +R__LOAD_LIBRARY(libg4detectors.so) +R__LOAD_LIBRARY(libg4eval.so) + +namespace Enable +{ + bool EEMCH = false; + bool EEMCH_ABSORBER = false; + bool EEMCH_CELL = false; + bool EEMCH_TOWER = false; + bool EEMCH_CLUSTER = false; + bool EEMCH_EVAL = false; + bool EEMCH_OVERLAPCHECK = false; + int EEMCH_VERBOSITY = 0; +} // namespace Enable + + +namespace G4EEMCH +{ + int use_projective_geometry = 0; + + // double Gdz = 18. + 0.0001; // These 2 paras are only served as the dimension of the black hole + // double Gz0 = -170.; + double Gdz = 20. + 0.1; + double Gz0 = -180.; + + namespace SETTING + { + bool USEHYBRID = true; + bool USECEMCGeo = false; + } // namespace SETTING + + // Digitization (default photon digi): + RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kSimple_photon_digitization; + // directly pass the energy of sim tower to digitized tower + // kNo_digitization + // simple digitization with photon statistics, single amplitude ADC conversion and pedestal + // kSimple_photon_digitization + // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal + // kSiPM_photon_digitization + + enum enu_Eemc_clusterizer + { + kEemcGraphClusterizer, + kEemcTemplateClusterizer + }; + //default template clusterizer, as developed by Sasha Bazilevsky + enu_Eemc_clusterizer Eemc_clusterizer = kEemcTemplateClusterizer; + // graph clusterizer + //enu_Eemc_clusterizer Eemc_clusterizer = kEemcGraphClusterizer; + +} // namespace G4EEMC + + +void EEMCHInit() +{ + BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 77.); + // from towerMap_EEMC_v006.txt + BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, G4EEMCH::Gz0 - G4EEMCH::Gdz / 2.); +} + + +void EEMCHSetup(PHG4Reco *g4Reco) +{ + bool AbsorberActive = Enable::ABSORBER || Enable::EEMCH_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::EEMCH_OVERLAPCHECK; + int verbosity = std::max(Enable::VERBOSITY, Enable::EEMCH_VERBOSITY); + + /** Use dedicated EEMCH module */ + ostringstream mapping_eemc_1, mapping_eemc_2; + + if (Enable::EEMCH_VERBOSITY > 0) cout << "hybrid: " << G4EEMCH::SETTING::USEHYBRID << "\t CEMC:" << G4EEMCH::SETTING::USECEMCGeo << endl; + + PHG4HybridHomogeneousCalorimeterSubsystem *eemc_crystal = new PHG4HybridHomogeneousCalorimeterSubsystem("EEMC"); + eemc_crystal->SuperDetector("EEMC"); + eemc_crystal->SetActive(); + if (AbsorberActive) + eemc_crystal->SetAbsorberActive(); + + if (G4EEMCH::SETTING::USEHYBRID && !G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_crystal_200cm_SciGlassBarrel.txt"; + else if (G4EEMCH::SETTING::USEHYBRID && G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_crystal_200cm_CEMCBarrel.txt"; + else if (!G4EEMCH::SETTING::USEHYBRID && !G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_purecrystal_185cm.txt"; + else if (!G4EEMCH::SETTING::USEHYBRID && G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_purecrystal_200cm_CEMCBarrel.txt"; + else { + cout << "*******************************************************************************" << endl; + cout << "****** ATTENTION no EEMC set as your settings aren't correct ******" << endl; + cout << "*******************************************************************************" << endl; + return; + } + if (Enable::EEMCH_VERBOSITY > 0) cout << "setting EEMC crystal mapping: " << mapping_eemc_1.str() << endl; + eemc_crystal->set_string_param("mappingtower", mapping_eemc_1.str()); + eemc_crystal->OverlapCheck(OverlapCheck); + + g4Reco->registerSubsystem(eemc_crystal); + + if (G4EEMCH::SETTING::USEHYBRID){ + if (!G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_2 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_glass_200cm_SciGlassBarrel.txt"; + else if ( G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_2 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_glass_200cm_CEMCBarrel.txt"; + else { + cout << "*******************************************************************************" << endl; + cout << "****** requested hybrid option but no glass mapping set ******" << endl; + cout << "*******************************************************************************" << endl; + return; + } + + if (Enable::EEMCH_VERBOSITY > 0) cout << "setting EEMC glass mapping: " << mapping_eemc_2.str() << endl; + PHG4HybridHomogeneousCalorimeterSubsystem *eemc_glass = new PHG4HybridHomogeneousCalorimeterSubsystem("EEMC_glass"); + eemc_glass->SuperDetector("EEMC_glass"); + eemc_glass->SetActive(); + if (AbsorberActive) + eemc_glass->SetAbsorberActive(); + + eemc_glass->set_string_param("mappingtower", mapping_eemc_2.str()); + + eemc_glass->OverlapCheck(OverlapCheck); + g4Reco->registerSubsystem(eemc_glass); + + } + +} + +void EEMCH_Cells() +{} + +void EEMCH_Towers() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::EEMCH_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + ostringstream mapping_eemc_1, mapping_eemc_2; + if (G4EEMCH::SETTING::USEHYBRID && !G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_crystal_200cm_SciGlassBarrel.txt"; + else if (G4EEMCH::SETTING::USEHYBRID && G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_crystal_200cm_CEMCBarrel.txt"; + else if (!G4EEMCH::SETTING::USEHYBRID && !G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_purecrystal_185cm.txt"; + else if (!G4EEMCH::SETTING::USEHYBRID && G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_1 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_purecrystal_200cm_CEMCBarrel.txt"; + + // CMS lead tungstate barrel ECAL at 18 degree centrigrade: 4.5 photoelectrons per MeV + // lead tungsten test in Orsay is 15~20 p.e. per MeV, sci-glass is 5 p.e. per MeV + const double EEMC_photoelectron_per_GeV_crystal = 15000; + const double EEMC_photoelectron_per_GeV_glass = 5000; + + //the original values are [8, 16], no noise case[0, 0], really high case[80, 160] + const double crystal_pedestal_ADC = 0, crystal_zero_suppression_ADC = 0; + const double glass_pedestal_ADC = 0, glass_zero_suppression_ADC = 0; + + RawTowerBuilderByHitIndex *tower_EEMC_crystal = new RawTowerBuilderByHitIndex("TowerBuilder_EEMC_crystal"); + tower_EEMC_crystal->Detector("EEMC"); + tower_EEMC_crystal->set_sim_tower_node_prefix("SIM"); + tower_EEMC_crystal->GeometryTableFile(mapping_eemc_1.str()); + se->registerSubsystem(tower_EEMC_crystal); + + // Calorimeter digitization + RawTowerDigitizer *TowerDigitizer_EEMC_crystal = new RawTowerDigitizer("EEMCRawTowerDigitizer_crystal"); + TowerDigitizer_EEMC_crystal->Detector("EEMC"); + TowerDigitizer_EEMC_crystal->Verbosity(verbosity); + TowerDigitizer_EEMC_crystal->set_raw_tower_node_prefix("RAW"); + TowerDigitizer_EEMC_crystal->set_digi_algorithm(G4EEMCH::TowerDigi); + TowerDigitizer_EEMC_crystal->set_pedstal_central_ADC(0); + TowerDigitizer_EEMC_crystal->set_pedstal_width_ADC(crystal_pedestal_ADC); // eRD1 test beam setting + TowerDigitizer_EEMC_crystal->set_photonelec_ADC(1); //not simulating ADC discretization error + TowerDigitizer_EEMC_crystal->set_photonelec_yield_visible_GeV(EEMC_photoelectron_per_GeV_crystal); + TowerDigitizer_EEMC_crystal->set_zero_suppression_ADC(crystal_zero_suppression_ADC); // eRD1 test beam setting + se->registerSubsystem(TowerDigitizer_EEMC_crystal); + + // Calorimeter calibration + RawTowerCalibration *TowerCalibration_EEMC_crystal = new RawTowerCalibration("EEMCRawTowerCalibration_crystal"); + TowerCalibration_EEMC_crystal->Detector("EEMC"); + TowerCalibration_EEMC_crystal->Verbosity(verbosity); + TowerCalibration_EEMC_crystal->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + if (G4EEMCH::TowerDigi == RawTowerDigitizer::kNo_digitization) + TowerCalibration_EEMC_crystal->set_calib_const_GeV_ADC(1.); + else + TowerCalibration_EEMC_crystal->set_calib_const_GeV_ADC(1. / EEMC_photoelectron_per_GeV_crystal); + TowerCalibration_EEMC_crystal->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration_EEMC_crystal); + + if (G4EEMCH::SETTING::USEHYBRID){ + if (!G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_2 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_glass_200cm_SciGlassBarrel.txt"; + else if ( G4EEMCH::SETTING::USECEMCGeo) + mapping_eemc_2 << getenv("CALIBRATIONROOT") << "/CrystalCalorimeter/mapping/crystal_mapping/tower_map_glass_200cm_CEMCBarrel.txt"; + + RawTowerBuilderByHitIndex *tower_EEMC_glass = new RawTowerBuilderByHitIndex("TowerBuilder_EEMC_glass"); + tower_EEMC_glass->Detector("EEMC_glass"); + tower_EEMC_glass->set_sim_tower_node_prefix("SIM"); + tower_EEMC_glass->GeometryTableFile(mapping_eemc_2.str()); + se->registerSubsystem(tower_EEMC_glass); + + RawTowerDigitizer *TowerDigitizer_EEMC_glass = new RawTowerDigitizer("EEMCRawTowerDigitizer_glass"); + TowerDigitizer_EEMC_glass->Detector("EEMC_glass"); + TowerDigitizer_EEMC_glass->Verbosity(verbosity); + TowerDigitizer_EEMC_glass->set_raw_tower_node_prefix("RAW"); + TowerDigitizer_EEMC_glass->set_digi_algorithm(G4EEMCH::TowerDigi); + TowerDigitizer_EEMC_glass->set_pedstal_central_ADC(0); + TowerDigitizer_EEMC_glass->set_pedstal_width_ADC(glass_pedestal_ADC); // eRD1 test beam setting + TowerDigitizer_EEMC_glass->set_photonelec_ADC(1); //not simulating ADC discretization error + TowerDigitizer_EEMC_glass->set_photonelec_yield_visible_GeV(EEMC_photoelectron_per_GeV_glass); + TowerDigitizer_EEMC_glass->set_zero_suppression_ADC(glass_zero_suppression_ADC); // eRD1 test beam setting + se->registerSubsystem(TowerDigitizer_EEMC_glass); + + RawTowerCalibration *TowerCalibration_EEMC_glass = new RawTowerCalibration("EEMCRawTowerCalibration_glass"); + TowerCalibration_EEMC_glass->Detector("EEMC_glass"); + TowerCalibration_EEMC_glass->Verbosity(verbosity); + TowerCalibration_EEMC_glass->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + if (G4EEMCH::TowerDigi == RawTowerDigitizer::kNo_digitization) + TowerCalibration_EEMC_glass->set_calib_const_GeV_ADC(1.); + else + TowerCalibration_EEMC_glass->set_calib_const_GeV_ADC(1. / EEMC_photoelectron_per_GeV_glass); + TowerCalibration_EEMC_glass->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration_EEMC_glass); + } +} + + +void EEMCH_Clusters() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::EEMCH_VERBOSITY); + Fun4AllServer *se = Fun4AllServer::instance(); + + if (G4EEMCH::Eemc_clusterizer == G4EEMCH::kEemcTemplateClusterizer) + { + + RawClusterBuilderTemplate *ClusterBuilder_crystal = new RawClusterBuilderTemplate("EEMCRawClusterBuilderTemplate_crystal"); + ClusterBuilder_crystal->Detector("EEMC"); + ClusterBuilder_crystal->Verbosity(2); + se->registerSubsystem(ClusterBuilder_crystal); + + if (G4EEMCH::SETTING::USEHYBRID){ + RawClusterBuilderTemplate *ClusterBuilder_glass = new RawClusterBuilderTemplate("EEMCRawClusterBuilderTemplate_glass"); + ClusterBuilder_glass->Detector("EEMC_glass"); + ClusterBuilder_glass->Verbosity(verbosity); + se->registerSubsystem(ClusterBuilder_glass); + } + } + else if (G4EEMCH::Eemc_clusterizer == G4EEMCH::kEemcGraphClusterizer) + { + + RawClusterBuilderFwd *ClusterBuilder_crystal = new RawClusterBuilderFwd("EEMCRawClusterBuilderFwd_crystal"); + ClusterBuilder_crystal->Detector("EEMC"); + ClusterBuilder_crystal->Verbosity(verbosity); + ClusterBuilder_crystal->Verbosity(2); + se->registerSubsystem(ClusterBuilder_crystal); + + if (G4EEMCH::SETTING::USEHYBRID){ + RawClusterBuilderFwd *ClusterBuilder_glass = new RawClusterBuilderFwd("EEMCRawClusterBuilderFwd_glass"); + ClusterBuilder_glass->Detector("EEMC_glass"); + ClusterBuilder_glass->Verbosity(verbosity); + se->registerSubsystem(ClusterBuilder_glass); + } + } + else + { + cout << "EEMC_Clusters - unknown clusterizer setting " << G4EEMCH::Eemc_clusterizer << endl; + gSystem->Exit(1); + } + return; +} + + +void EEMCH_Eval(const std::string &outputfile) +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::EEMCH_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + string outputroot = outputfile; + string remove_this = ".root"; + size_t pos = outputroot.find(remove_this); + if (pos != string::npos){ + outputroot.erase(pos, remove_this.length()); + } + string outputrootc = outputroot+"_crystal.root"; + string outputrootg = outputroot+"_glass.root"; + + CaloEvaluator *eval_crystal = new CaloEvaluator("EEMCEVALUATOR", "EEMC", outputrootc.c_str()); + eval_crystal->Verbosity(verbosity); + eval_crystal->set_do_cluster_eval(true); + se->registerSubsystem(eval_crystal); + + if (G4EEMCH::SETTING::USEHYBRID){ + CaloEvaluator *eval_glass = new CaloEvaluator("EEMCGLASSEVALUATOR", "EEMC_glass", outputrootg.c_str()); + eval_glass->Verbosity(verbosity); + eval_glass->set_do_cluster_eval(true); + se->registerSubsystem(eval_glass); + } + return; +} +#endif diff --git a/common/G4_EventEvaluator.C b/common/G4_EventEvaluator.C index d864dc67..28da71b7 100644 --- a/common/G4_EventEvaluator.C +++ b/common/G4_EventEvaluator.C @@ -2,10 +2,10 @@ #define MACRO_EventEvaluator_C #include -#include +#include R__LOAD_LIBRARY(libfun4all.so) -R__LOAD_LIBRARY(libg4eval.so) +R__LOAD_LIBRARY(libeiceval.so) namespace Enable { @@ -23,8 +23,8 @@ void Event_Eval(const std::string &filename) { Fun4AllServer *se = Fun4AllServer::instance(); - EventEvaluator *eval = new EventEvaluator("EVENTEVALUATOR", filename); - eval->set_reco_tracing_energy_threshold(EVENT_EVALUATOR::EnergyThreshold); + EventEvaluatorEIC *eval = new EventEvaluatorEIC("EVENTEVALUATOR", filename); + eval->set_reco_tracing_energy_thresholdMC(EVENT_EVALUATOR::EnergyThreshold); eval->Verbosity(EVENT_EVALUATOR::Verbosity); if (Enable::TRACKING) @@ -41,8 +41,10 @@ void Event_Eval(const std::string &filename) if (Enable::HCALIN_CLUSTER) eval->set_do_HCALIN(true); if (Enable::HCALOUT_CLUSTER) eval->set_do_HCALOUT(true); if (Enable::FHCAL_CLUSTER) eval->set_do_FHCAL(true); - if (Enable::FHCAL_CLUSTER || Enable::FEMC_CLUSTER || Enable::EEMC_CLUSTER) - eval->set_do_CLUSTERS(true); + if (Enable::FHCAL_CLUSTER || Enable::FEMC_CLUSTER || Enable::EEMC_CLUSTER) eval->set_do_CLUSTERS(true); + if (Enable::DRCALO_CLUSTER) eval->set_do_DRCALO(true); + if (Enable::LFHCAL_CLUSTER) eval->set_do_LFHCAL(true); + if (Enable::BECAL) eval->set_do_BECAL(true); eval->set_do_MCPARTICLES(true); eval->set_do_HEPMC(true); diff --git a/common/G4_FST_EIC.C b/common/G4_FST_EIC.C index cb891dd0..77bc4284 100644 --- a/common/G4_FST_EIC.C +++ b/common/G4_FST_EIC.C @@ -59,18 +59,17 @@ void FSTSetup(PHG4Reco *g4Reco, const double min_eta = 1.245) const double um = 1e-3 * mm; //Design from Xuan Li @LANL - make_LANL_FST_station("FST_0", g4Reco, 35, 4, 22, 35 * um, 20e-4); //cm - make_LANL_FST_station("FST_1", g4Reco, 57.5, 4.5, 42, 35 * um, 20e-4); - make_LANL_FST_station("FST_2", g4Reco, 80, 6, 43.5, 35 * um, 20e-4); - make_LANL_FST_station("FST_3", g4Reco, 115, 9.3, 46.8, 85 * um, 36.4e-4); - make_LANL_FST_station("FST_4", g4Reco, 125, 9.6, 47.1, 85 * um, 36.4e-4); - - //mirror for e-going FST - make_LANL_FST_station("EFST_0", g4Reco, -35, 4, 22, 35 * um, 20e-4); //cm - make_LANL_FST_station("EFST_1", g4Reco, -57.5, 4.5, 42, 35 * um, 20e-4); - make_LANL_FST_station("EFST_2", g4Reco, -80, 6, 43.5, 35 * um, 20e-4); - make_LANL_FST_station("EFST_3", g4Reco, -115, 9.3, 46.8, 85 * um, 36.4e-4); - make_LANL_FST_station("EFST_4", g4Reco, -125, 9.6, 47.1, 85 * um, 36.4e-4); + make_LANL_FST_station("FST_0", g4Reco, 35, 4, 22, 35 * um, 10e-4); //cm + make_LANL_FST_station("FST_1", g4Reco, 57.5, 4.5, 42, 35 * um, 10e-4); + make_LANL_FST_station("FST_2", g4Reco, 80, 6, 43.5, 35 * um, 10e-4); + make_LANL_FST_station("FST_3", g4Reco, 115, 9.3, 46.8, 35 * um, 10e-4); + make_LANL_FST_station("FST_4", g4Reco, 125, 9.6, 47.1, 35 * um, 10e-4); + + //mirror for e-going FST, after update with new e-spectrometer + make_LANL_FST_station("EFST_0", g4Reco, -35, 4, 22, 35 * um, 10e-4); //cm + make_LANL_FST_station("EFST_1", g4Reco, -57.5, 4.5, 42, 35 * um, 10e-4); + make_LANL_FST_station("EFST_2", g4Reco, -80, 6, 43.5, 35 * um, 10e-4); + make_LANL_FST_station("EFST_3", g4Reco, -107.1, 9.3, 46.8, 35 * um, 10e-4); if (G4FST::SETTING::SUPPORTCYL) { @@ -140,6 +139,20 @@ int make_LANL_FST_station(string name, PHG4Reco *g4Reco, 50e-4 / sqrt(12.), // const float lonres, *ignored in plane detector* 1, // const float eff, 0); // const float noise + TRACKING::FastKalmanFilterInnerTrack->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, + PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + pitch / sqrt(12.), // const float radres, + pitch / sqrt(12.), // const float phires, + 50e-4 / sqrt(12.), // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + TRACKING::FastKalmanFilterSiliconTrack->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, + PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + pitch / sqrt(12.), // const float radres, + pitch / sqrt(12.), // const float phires, + 50e-4 / sqrt(12.), // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise } return 0; } diff --git a/common/G4_GEM_EIC.C b/common/G4_GEM_EIC.C index 676236f8..dbc0303a 100644 --- a/common/G4_GEM_EIC.C +++ b/common/G4_GEM_EIC.C @@ -48,14 +48,14 @@ void BGEM_Init() void EGEMSetup(PHG4Reco *g4Reco) { - make_GEM_station("EGEM_0", g4Reco, -160.0, -1.65, -3.5); -// make_GEM_station("EGEM_1", g4Reco, -190.0, -1.85, -3.6); // replaced by LGAD-TTL + make_GEM_station("EGEM_0", g4Reco, -120.0, -1.68, -3.5); + // make_GEM_station("EGEM_1", g4Reco, -190.0, -1.85, -3.6); // replaced by LGAD-TTL } void FGEMSetup(PHG4Reco *g4Reco, const int N_Sector = 8) { make_GEM_station("FGEM_0", g4Reco, 160.0, 1.65, 3.5, N_Sector); -// make_GEM_station("FGEM_1", g4Reco, 285.0, 1.2, 3.5, N_Sector);// replaced by LGAD-TTL + // make_GEM_station("FGEM_1", g4Reco, 285.0, 1.2, 3.5, N_Sector);// replaced by LGAD-TTL } //! Add drift layers to mini TPC @@ -169,11 +169,20 @@ int make_GEM_station(string name, PHG4Reco *g4Reco, double zpos, double etamin, TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, 1. / sqrt(12.), // const float radres, - 70e-4, // const float phires, + 50e-4, // const float phires, 100e-4, // const float lonres, 1, // const float eff, 0); // const float noise + if (TRACKING::FastKalmanFilterInnerTrack) + TRACKING::FastKalmanFilterInnerTrack->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, + PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + 1. / sqrt(12.), // const float radres, + 50e-4, // const float phires, + 100e-4, // const float lonres, + 1, // const float eff, + 0); // const float noise + return 0; } #endif diff --git a/common/G4_Input.C b/common/G4_Input.C index 436f3cf9..71dd8559 100644 --- a/common/G4_Input.C +++ b/common/G4_Input.C @@ -121,51 +121,57 @@ namespace Input exit(1); } + HepMCGen->PHHepMCGenHelper_Verbosity(10); + //25mrad x-ing as in EIC CDR const double EIC_hadron_crossing_angle = 25e-3; + // beta* for 275*x18 collisions + // Table 4 of + // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf + const double beta_star_p_h = 80; + const double beta_star_p_v = 7.1; + const double beta_star_e_h = 59; + const double beta_star_e_v = 5.7; + // Table 1-2 of + // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf + const double beta_crab_p = 1300e2; + const double beta_crab_e = 150e2; HepMCGen->set_beam_direction_theta_phi( EIC_hadron_crossing_angle, // beamA_theta - M_PI, // beamA_phi + M_PI, // beamA_phi M_PI, // beamB_theta 0 // beamB_phi ); + // Table 4 of + // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf HepMCGen->set_beam_angular_divergence_hv( - 119e-6, 119e-6, // proton beam divergence horizontal & vertical, as in EIC CDR Table 1.1 - 211e-6, 152e-6 // electron beam divergence horizontal & vertical, as in EIC CDR Table 1.1 + 150e-6, 150e-6, // proton beam divergence horizontal & vertical + 202e-6, 187e-6 // electron beam divergence horizontal & vertical ); + // vertex shape from beam_bunch_sim + HepMCGen->use_beam_bunch_sim(true); + // angular kick within a bunch as result of crab cavity - // using an naive assumption of transfer matrix from the cavity to IP, - // which is NOT yet validated with accelerator optics simulations! - const double z_hadron_cavity = 52e2; // CDR Fig 3.3 - const double z_e_cavity = 38e2; // CDR Fig 3.2 + // Eq 5 of + // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf HepMCGen->set_beam_angular_z_coefficient_hv( - -EIC_hadron_crossing_angle / 2. / z_hadron_cavity, 0, - -EIC_hadron_crossing_angle / 2. / z_e_cavity, 0); + -EIC_hadron_crossing_angle / 2. / sqrt(beta_star_p_h * beta_crab_p), 0, + -EIC_hadron_crossing_angle / 2. / sqrt(beta_star_e_h * beta_crab_e), 0); - // calculate beam sigma width at IP as in EIC CDR table 1.1 - const double sigma_p_h = sqrt(80 * 11.3e-7); - const double sigma_p_v = sqrt(7.2 * 1.0e-7); + // Table 4 of + // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf + const double sigma_p_h = sqrt(beta_star_p_h * 18e-7); + const double sigma_p_v = sqrt(beta_star_p_v * 1.6e-7); const double sigma_p_l = 6; - const double sigma_e_h = sqrt(45 * 20.0e-7); - const double sigma_e_v = sqrt(5.6 * 1.3e-7); - const double sigma_e_l = 2; - - // combine two beam gives the collision sigma in z - const double collision_sigma_z = sqrt(sigma_p_l * sigma_p_l + sigma_e_l * sigma_e_l) / 2; - const double collision_sigma_t = collision_sigma_z / 29.9792; // speed of light in cm/ns + const double sigma_e_h = sqrt(beta_star_e_h * 24e-7); + const double sigma_e_v = sqrt(beta_star_e_v * 2.0e-7); + const double sigma_e_l = 0.9; - HepMCGen->set_vertex_distribution_width( - sigma_p_h * sigma_e_h / sqrt(sigma_p_h * sigma_p_h + sigma_e_h * sigma_e_h), //x - sigma_p_v * sigma_e_v / sqrt(sigma_p_v * sigma_p_v + sigma_e_v * sigma_e_v), //y - collision_sigma_z, //z - collision_sigma_t); //t - HepMCGen->set_vertex_distribution_function( - PHHepMCGenHelper::Gaus, //x - PHHepMCGenHelper::Gaus, //y - PHHepMCGenHelper::Gaus, //z - PHHepMCGenHelper::Gaus); //t + HepMCGen->set_beam_bunch_width( + std::vector{sigma_p_h, sigma_p_v, sigma_p_l}, + std::vector{sigma_e_h, sigma_e_v, sigma_e_l}); } //! apply EIC beam parameter to any HepMC generator following EIC CDR, diff --git a/common/G4_LFHCAL.C b/common/G4_LFHCAL.C new file mode 100644 index 00000000..b6af8bb3 --- /dev/null +++ b/common/G4_LFHCAL.C @@ -0,0 +1,211 @@ +#ifndef MACRO_G4LFHCAL_C +#define MACRO_G4LFHCAL_C + +#include + +#include + +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +R__LOAD_LIBRARY(libcalo_reco.so) +R__LOAD_LIBRARY(libg4calo.so) +R__LOAD_LIBRARY(libg4detectors.so) +R__LOAD_LIBRARY(libg4eval.so) + +namespace Enable +{ + bool LFHCAL = false; + bool LFHCAL_ABSORBER = false; + bool LFHCAL_CELL = false; + bool LFHCAL_TOWER = false; + bool LFHCAL_CLUSTER = false; + bool LFHCAL_EVAL = false; + bool LFHCAL_OVERLAPCHECK = false; + int LFHCAL_VERBOSITY = 0; +} // namespace Enable + +namespace G4LFHCAL +{ + // from LFHcal/mapping/towerMap_LFHCAL_v005.txt + double Gz0 = 400.; + double Gdz = 100.; + double outer_radius = 262.; + enum enu_FHcal_clusterizer + { + kFHcalGraphClusterizer, + kFHcalTemplateClusterizer + }; + //template clusterizer, as developed by Sasha Bazilevsky + enu_FHcal_clusterizer FHcal_clusterizer = kFHcalTemplateClusterizer; + // graph clusterizer + //enu_FHcal_clusterizer FHcal_clusterizer = kFHcalGraphClusterizer; + namespace SETTING + { + bool FullEtaAcc = false; + bool HC2x = false; + bool asymmetric = false; + bool wDR = false; + bool FwdSquare = false; + bool longer = false; + } // namespace SETTING +} // namespace G4LFHCAL + +TString GetMappingFile(){ + TString mappingFileName = getenv("CALIBRATIONROOT"); + if (G4LFHCAL::SETTING::HC2x ) + { + if (G4LFHCAL::SETTING::longer) + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_2x-long.txt"; + else + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_2x.txt"; + } + // HCal Fe-Scint surrounding dual readout calorimeter R>50cm + else if (G4LFHCAL::SETTING::wDR) + { + if (G4LFHCAL::SETTING::longer) + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_wDR-long.txt"; + else + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_wDR.txt"; + } + // HCal Fe-Scint surrounding dual readout calorimeter R>50cm + else if (G4LFHCAL::SETTING::FwdSquare) + { + if (G4LFHCAL::SETTING::longer) + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_FwdSquare-long.txt"; + else + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_FwdSquare.txt"; + } + // full HCal Fe-Scint with asymmetric centering around beampipe + else if (G4LFHCAL::SETTING::asymmetric) + { + if (G4LFHCAL::SETTING::longer) + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_asymmetric-long.txt"; + else + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_asymmetric.txt"; + + } + // PSD like HCal Fe-Scint with enlarged beam pipe opening for Mar 2020 beam pipe + else + { + if (G4LFHCAL::SETTING::longer) + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_default-long.txt"; + else + mappingFileName += "/LFHcal/mapping/towerMap_LFHCAL_default.txt"; + } + + return mappingFileName; + +} + + +void LFHCALInit() +{ + if (G4LFHCAL::SETTING::longer){ + G4LFHCAL::Gz0 = 420; + G4LFHCAL::Gdz = 140; + } + // simple way to check if only 1 of the settings is true + if ((G4LFHCAL::SETTING::FullEtaAcc ? 1 : 0) + (G4LFHCAL::SETTING::HC2x ? 1 : 0) > 1) + { + cout << "use only G4LFHCAL::SETTING::FullEtaAcc=true or G4LFHCAL::SETTING::HC2x=true or G4LFHCAL::SETTING::HC4x=true" << endl; + gSystem->Exit(1); + } + + BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, G4LFHCAL::outer_radius); + BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, G4LFHCAL::Gz0 + G4LFHCAL::Gdz / 2.); +} + +void LFHCALSetup(PHG4Reco *g4Reco) +{ + const bool AbsorberActive = Enable::ABSORBER || Enable::LFHCAL_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::LFHCAL_OVERLAPCHECK; + Fun4AllServer *se = Fun4AllServer::instance(); + + /** Use dedicated LFHCAL module */ + PHG4LFHcalSubsystem *fhcal = new PHG4LFHcalSubsystem("LFHCAL"); + + TString mapping_fhcal = GetMappingFile(); + cout << "LFHCAL: "<< mapping_fhcal.Data() << endl; + ostringstream mapping_fhcal_s; + mapping_fhcal_s << mapping_fhcal.Data(); + + fhcal->SetTowerMappingFile(mapping_fhcal_s.str()); + fhcal->OverlapCheck(OverlapCheck); + fhcal->SetActive(); + fhcal->SetDetailed(true); + fhcal->SuperDetector("LFHCAL"); + if (AbsorberActive) fhcal->SetAbsorberActive(); + + g4Reco->registerSubsystem(fhcal); +} + +void LFHCAL_Cells(int verbosity = 0) +{ + return; +} + +void LFHCAL_Towers() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::LFHCAL_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + // Switch to desired calo setup; + // PSD like HCal Fe-Scint with doubled granularity + TString mapping_fhcal = GetMappingFile(); + ostringstream mapping_fhcal_s; + mapping_fhcal_s << mapping_fhcal.Data(); + + RawTowerBuilderByHitIndexLHCal *tower_LFHCAL = new RawTowerBuilderByHitIndexLHCal("TowerBuilder_LFHCAL"); + tower_LFHCAL->Detector("LFHCAL"); + tower_LFHCAL->set_sim_tower_node_prefix("SIM"); + tower_LFHCAL->GeometryTableFile(mapping_fhcal_s.str()); + + se->registerSubsystem(tower_LFHCAL); + + cout << "def: using default for LFHCAL towers" << endl; + RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("LFHCALRawTowerDigitizer"); + TowerDigitizer->Detector("LFHCAL"); + TowerDigitizer->Verbosity(verbosity); + TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kNo_digitization); + se->registerSubsystem(TowerDigitizer); + + + RawTowerCalibration *TowerCalibration = new RawTowerCalibration("LFHCALRawTowerCalibration"); + TowerCalibration->Detector("LFHCAL"); + TowerCalibration->Verbosity(verbosity); + TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + TowerCalibration->set_calib_const_GeV_ADC(1. / (0.03898*0.93)); // temporary factor 0.93 to fix calibration for new tower design + TowerCalibration->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration); +} + +void LFHCAL_Clusters() +{ + + return; +} + +void LFHCAL_Eval(const std::string &outputfile) +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::LFHCAL_VERBOSITY); + Fun4AllServer *se = Fun4AllServer::instance(); + + CaloEvaluator *eval = new CaloEvaluator("LFHCALEVALUATOR", "LFHCAL", outputfile.c_str()); + eval->Verbosity(verbosity); + se->registerSubsystem(eval); + + return; +} +#endif diff --git a/common/G4_Magnet.C b/common/G4_Magnet.C index 970cdc30..1cf7bcf9 100644 --- a/common/G4_Magnet.C +++ b/common/G4_Magnet.C @@ -33,6 +33,7 @@ void MagnetFieldInit() if (G4MAGNET::magfield.empty()) { G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root"); +cout << "Using magnetic field map: " << G4MAGNET::magfield << endl; } } diff --git a/common/G4_Pipe_EIC.C b/common/G4_Pipe_EIC.C index 456f1c11..a6a6599d 100644 --- a/common/G4_Pipe_EIC.C +++ b/common/G4_Pipe_EIC.C @@ -28,6 +28,7 @@ namespace G4PIPE // directly implimenting the central Be section in G4 cylinder for max speed simulation in the detector region. // The jointer lip structure of the pipe R = 3.2cm x L=5mm is ignored here double be_pipe_radius = 3.1000; + double Au_coating_thickness = 2e-4; // 2um Au coating double be_pipe_thickness = 3.1762 - be_pipe_radius; // 760 um for sPHENIX double be_pipe_length_plus = 66.8; // +z beam pipe extend. double be_pipe_length_neg = -79.8; // -z beam pipe extend. @@ -78,9 +79,21 @@ double Pipe(PHG4Reco* g4Reco, cyl->set_double_param("length", be_pipe_length); cyl->set_double_param("place_z", be_pipe_center); cyl->set_string_param("material", "G4_Galactic"); - cyl->set_double_param("thickness", G4PIPE::be_pipe_radius); + cyl->set_double_param("thickness", G4PIPE::be_pipe_radius- G4PIPE::Au_coating_thickness); cyl->SuperDetector("PIPE"); - cyl->OverlapCheck(OverlapCheck); + cyl->OverlapCheck(true); + g4Reco->registerSubsystem(cyl); + + // 2um Au coating for X-ray absorption + cyl = new PHG4CylinderSubsystem("Au_PIPE", 2); + cyl->set_double_param("radius", G4PIPE::be_pipe_radius - G4PIPE::Au_coating_thickness); + cyl->set_int_param("lengthviarapidity", 0); + cyl->set_double_param("length", be_pipe_length); + cyl->set_double_param("place_z", be_pipe_center); + cyl->set_string_param("material", "G4_Au"); + cyl->set_double_param("thickness", G4PIPE::Au_coating_thickness); + cyl->SuperDetector("PIPE"); + cyl->OverlapCheck(true); if (AbsorberActive) cyl->SetActive(); g4Reco->registerSubsystem(cyl); @@ -92,7 +105,7 @@ double Pipe(PHG4Reco* g4Reco, cyl->set_string_param("material", "G4_Be"); cyl->set_double_param("thickness", G4PIPE::be_pipe_thickness); cyl->SuperDetector("PIPE"); - cyl->OverlapCheck(OverlapCheck); + cyl->OverlapCheck(true); if (AbsorberActive) cyl->SetActive(); g4Reco->registerSubsystem(cyl); diff --git a/common/G4_TTL_EIC.C b/common/G4_TTL_EIC.C index 2bcd8be7..206dca8d 100644 --- a/common/G4_TTL_EIC.C +++ b/common/G4_TTL_EIC.C @@ -3,9 +3,8 @@ #include "GlobalVariables.C" +#include #include -#include -#include #include @@ -14,9 +13,13 @@ R__LOAD_LIBRARY(libg4detectors.so) int make_forward_station(string name, PHG4Reco *g4Reco, double zpos, double Rmin, - double Rmax, double tSilicon); -int make_barrel_layer(string name, PHG4Reco *g4Reco, - double radius, double z_start, double z_end, double tSilicon); + double Rmax,double tSilicon, double xoffset=0); +int make_forward_station_basic(string name, PHG4Reco *g4Reco, double zpos, double Rmin, + double Rmax,double tSilicon); +int make_barrel_layer_basic(string name, PHG4Reco *g4Reco, + double radius, double halflength, double tSilicon, double zOffset); +int make_barrel_layer(string name, PHG4Reco *g4Reco, + double radius, double halflength, double tSilicon, double zOffset); //-----------------------------------------------------------------------------------// namespace Enable @@ -24,11 +27,26 @@ namespace Enable bool FTTL = false; bool ETTL = false; bool CTTL = false; - bool TTL_OVERLAPCHECK = false; -} // namespace Enable +} -namespace TTL +namespace G4TTL { + int layer[3] = { 2, 1, 2}; + double positionToVtx[3][3] = { {-169., -172., -309.5}, {80., 114.7, 0. }, { 287., 289., 340.} }; + double minExtension[3][3] = { {8, 8, 15.3}, {218, 180, 0 }, {11.62, 11.7, 13.8 } }; + double maxExtension[3][3] = { {61., 61. , 200}, {-40, 0, 0 }, {170., 170., 250 } }; + namespace SETTING + { + bool optionCEMC = true; + bool optionEEMCH = true; + bool optionBasicGeo = false; + int optionDR = 0; + int optionGeo = 1; + int optionGran = 1; + } // namespace SETTING + + + // 2, LGAD based ToF by Wei Li: // Present the recent simulation studies and detector layout for the LGAD based ToF, // which can be used as an outer tracker. The current detector can reach around 30 micron spatial @@ -37,62 +55,195 @@ namespace TTL // reconstruction within high pseudorapidity regions. const double PositionResolution(30e-4); -} // namespace TTL + +} // namespace G4TTL + //-----------------------------------------------------------------------------------// void TTL_Init() { - BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 200.); + BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 250.); BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 350.); + + if (!G4TTL::SETTING::optionEEMCH){ + G4TTL::positionToVtx[0][0] = -155.5; + G4TTL::positionToVtx[0][1] = -158.5; + } + + if (G4TTL::SETTING::optionCEMC){ + G4TTL::maxExtension[0][0] = 80.; + G4TTL::maxExtension[0][1] = 80.; + G4TTL::positionToVtx[1][0] = 92.; + if(!G4TTL::SETTING::optionBasicGeo) G4TTL::positionToVtx[1][0] = 89.; + } + + if(G4TTL::SETTING::optionBasicGeo){ + G4TTL::minExtension[1][0] = 180; + G4TTL::maxExtension[1][0] = -25; + } + + if (G4TTL::SETTING::optionGeo == 1){ + cout << "TTL setup infront of ECals with 2 layers fwd/bwd & 1 layer barrel" << endl; + } if (G4TTL::SETTING::optionGeo == 2){ + cout << "TTL setup infront of ECals with 2 layers fwd/bwd & 1 layer barrel, lower barrel layer" << endl; + G4TTL::positionToVtx[1][0] = 50.; + if(!G4TTL::SETTING::optionBasicGeo) G4TTL::positionToVtx[1][0] = 65.; + G4TTL::minExtension[1][0] = 100.; + G4TTL::maxExtension[1][0] = 0.; + } if (G4TTL::SETTING::optionGeo == 3){ + cout << "TTL setup infront of ECals with 1 layers fwd/bwd & 1 layer barrel, lower barrel layer" << endl; + G4TTL::positionToVtx[1][0] = 50.; + if(!G4TTL::SETTING::optionBasicGeo) G4TTL::positionToVtx[1][0] = 65.; + G4TTL::minExtension[1][0] = 100.; + G4TTL::maxExtension[1][0] = 0.; + G4TTL::layer[0] = 1; + G4TTL::layer[2] = 1; + G4TTL::positionToVtx[0][0] = G4TTL::positionToVtx[0][1]; + G4TTL::positionToVtx[2][0] = G4TTL::positionToVtx[2][1]; + G4TTL::minExtension[0][0] = G4TTL::minExtension[0][1]; + G4TTL::minExtension[2][0] = G4TTL::minExtension[2][1]; + G4TTL::maxExtension[0][0] = G4TTL::maxExtension[0][1]; + G4TTL::maxExtension[2][0] = G4TTL::maxExtension[2][1]; + } if (G4TTL::SETTING::optionGeo == 4){ + cout << "TTL setup infront of ECals with 2 layers fwd/bwd & 1 layer barrel, 1 layer before HCals everywhere" << endl; + G4TTL::layer[0] = 3; + G4TTL::layer[1] = 2; + if(!G4TTL::SETTING::optionBasicGeo) G4TTL::layer[1] = 1; + G4TTL::layer[2] = 3; + } + + if (G4TTL::SETTING::optionDR == 2 && G4TTL::SETTING::optionGeo == 4 ){ + cout << "TTL setup infront of ECals with 2 layers fwd/bwd & 1 layer barrel, 1 layer before HCals everywhere choosen!" << endl; + cout << "conflicting in forward region with DR calo, reducing by 1 layer!" << endl; + G4TTL::layer[2] = 2; + } else if (G4TTL::SETTING::optionDR == 1 && G4TTL::SETTING::optionGeo == 4 ){ + cout << "TTL setup infront of ECals with 2 layers fwd/bwd & 1 layer barrel, 1 layer before HCals everywhere choosen!" << endl; + cout << "conflicting in forward region with DR calo, reducing adding larger cutout!" << endl; + G4TTL::minExtension[2][2] = 2.5; + } + } //-----------------------------------------------------------------------------------// -void FTTLSetup(PHG4Reco *g4Reco) +void FTTLSetup(PHG4Reco *g4Reco, TString fttloption = "") { const double cm = PHG4Sector::Sector_Geometry::Unit_cm(); const double mm = .1 * cm; const double um = 1e-3 * mm; - make_forward_station("FTTL_0", g4Reco, 287, 3.5, 1.3, 85 * um); - make_forward_station("FTTL_1", g4Reco, 289, 3.5, 1.3, 85 * um); + for (Int_t i = 0; i < G4TTL::layer[2]; i++){ + cout << G4TTL::positionToVtx[2][i] << "\t" << G4TTL::minExtension[2][i] << "\t" << G4TTL::maxExtension[2][i] << endl; + if(!G4TTL::SETTING::optionBasicGeo){ + make_forward_station(Form("FTTL_%d", i), g4Reco, G4TTL::positionToVtx[2][i], G4TTL::minExtension[2][i], G4TTL::maxExtension[2][i], 85*um, -6.0); + } else { + make_forward_station_basic(Form("FTTL_%d", i), g4Reco, G4TTL::positionToVtx[2][i], G4TTL::minExtension[2][i], G4TTL::maxExtension[2][i], 85*um); + } + } } + //-----------------------------------------------------------------------------------// -void ETTLSetup(PHG4Reco *g4Reco) +void ETTLSetup(PHG4Reco *g4Reco, TString ettloption = "") { const double cm = PHG4Sector::Sector_Geometry::Unit_cm(); const double mm = .1 * cm; const double um = 1e-3 * mm; - - make_forward_station("ETTL_0", g4Reco, -190, -1.85, -3.5, 85 * um); // define wit eta - make_forward_station("ETTL_1", g4Reco, -192, -1.85, -3.5, 85 * um); // define wit eta + for (Int_t i = 0; i < G4TTL::layer[0]; i++){ + cout << G4TTL::positionToVtx[0][i] << "\t" << G4TTL::minExtension[0][i] << "\t" << G4TTL::maxExtension[0][i] << endl; + if(!G4TTL::SETTING::optionBasicGeo){ + make_forward_station(Form("ETTL_%d", i), g4Reco, G4TTL::positionToVtx[0][i], G4TTL::minExtension[0][i], G4TTL::maxExtension[0][i], 85*um, 0); + } else { + make_forward_station_basic(Form("ETTL_%d", i), g4Reco, G4TTL::positionToVtx[0][i], G4TTL::minExtension[0][i], G4TTL::maxExtension[0][i], 85*um); + } + } } //-----------------------------------------------------------------------------------// -void CTTLSetup(PHG4Reco *g4Reco) +void CTTLSetup(PHG4Reco *g4Reco, TString cttloption = "") { const double cm = PHG4Sector::Sector_Geometry::Unit_cm(); const double mm = .1 * cm; const double um = 1e-3 * mm; - make_barrel_layer("CTTL_0", g4Reco, 80, -250, 180, 85 * um); + for (Int_t i = 0; i < G4TTL::layer[1]; i++){ + cout << "Radius: " << G4TTL::positionToVtx[1][i] << "\tLength: " << G4TTL::minExtension[1][i] << "\tz-Offset: " << G4TTL::maxExtension[1][i] << endl; + if(G4TTL::SETTING::optionBasicGeo){ + make_barrel_layer_basic(Form("CTTL_%d",i), g4Reco, G4TTL::positionToVtx[1][i], G4TTL::minExtension[1][i], 85*um, G4TTL::maxExtension[1][i]); + } else { + make_barrel_layer(Form("CTTL_%d",i), g4Reco, G4TTL::positionToVtx[1][i], G4TTL::minExtension[1][i], 85*um, G4TTL::maxExtension[1][i]); + } + } } + //-----------------------------------------------------------------------------------// int make_forward_station(string name, PHG4Reco *g4Reco, - double zpos, double etamin, double etamax, - double tSilicon) //silicon thickness + double zpos, double rMin, double rMax, + double tSilicon, //silicon thickness + double xoffset = 0 ) { - bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TTL_OVERLAPCHECK; - + cout << "r min: " << rMin << "\t r max: " << rMax << "\t z: " << zpos << endl; + // always facing the interaction point double polar_angle = 0; - if (zpos < 0) - { + if (zpos < 0){ zpos = -zpos; polar_angle = M_PI; } - if (etamax < etamin) + PHG4TTLSubsystem *ttl; + ttl = new PHG4TTLSubsystem(name); + ttl->SetDetailed(false); + ttl->SuperDetector(name); + ttl->set_double_param("polar_angle", polar_angle); // + ttl->set_int_param("isForward", 1); // + ttl->set_double_param("place_z", zpos * cm); // + ttl->set_double_param("rMin", rMin * cm); // + ttl->set_double_param("rMax", rMax * cm); // + ttl->set_double_param("offset_x", xoffset * cm); // + ttl->set_double_param("tSilicon", tSilicon); // +// ttl->OverlapCheck(true); + ttl->OverlapCheck(false); + + g4Reco->registerSubsystem(ttl); + + + if (TRACKING::FastKalmanFilter) { + TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, + PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4TTL::PositionResolution, // const float radres, + G4TTL::PositionResolution, // const float phires, + tSilicon / sqrt(12.), // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + TRACKING::FastKalmanFilter->add_zplane_state(name, zpos); + TRACKING::ProjectionNames.insert(name); + } + + return 0; +} + + +//-----------------------------------------------------------------------------------// +int make_forward_station_basic(string name, PHG4Reco *g4Reco, + double zpos, double rmin, double rmax, + double tSilicon) //silicon thickness +{ + // cout + // << "make_GEM_station - GEM construction with PHG4SectorSubsystem - make_GEM_station_EdgeReadout of " + // << name << endl; + + // always facing the interaction point + double etamin = -TMath::Log(TMath::Tan(rmin/TMath::Abs(zpos)/2)); + double etamax = -TMath::Log(TMath::Tan(rmax/TMath::Abs(zpos)/2)); + + double polar_angle = 0; + if (zpos < 0){ + zpos = -zpos; + polar_angle = M_PI; + etamin = -etamin; + etamax = -etamax; + } + if (etamax < etamin){ double t = etamax; etamax = etamin; etamin = t; @@ -111,8 +262,8 @@ int make_forward_station(string name, PHG4Reco *g4Reco, ttl->get_geometry().set_min_polar_edge(PHG4Sector::Sector_Geometry::ConeEdge()); ttl->get_geometry().set_N_Sector(1); ttl->get_geometry().set_material("G4_AIR"); - ttl->OverlapCheck(OverlapCheck); - + ttl->OverlapCheck(true); + const double cm = PHG4Sector::Sector_Geometry::Unit_cm(); const double mm = .1 * cm; const double um = 1e-3 * mm; @@ -128,81 +279,120 @@ int make_forward_station(string name, PHG4Reco *g4Reco, g4Reco->registerSubsystem(ttl); + if (TRACKING::FastKalmanFilter) { TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, PHG4TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, - TTL::PositionResolution, // const float radres, - TTL::PositionResolution, // const float phires, + G4TTL::PositionResolution, // const float radres, + G4TTL::PositionResolution, // const float phires, tSilicon / sqrt(12.), // const float lonres, *ignored in plane detector* 1, // const float eff, 0); // const float noise TRACKING::FastKalmanFilter->add_zplane_state(name, zpos); TRACKING::ProjectionNames.insert(name); } - return 0; } //-----------------------------------------------------------------------------------// -int make_barrel_layer(string name, PHG4Reco *g4Reco, - double radius, double z_start, double z_end, double tSilicon) +int make_barrel_layer(string name, PHG4Reco *g4Reco, + double radius, double halflength, double tSilicon, double zOffset ) { - bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TTL_OVERLAPCHECK; + // cout << "r min: " << rMin << "\t r max: " << rMax << "\t z: " << zpos << endl; + if(G4TTL::SETTING::optionCEMC){ + cout << "The improved barrel TTL layer is placed at a fixed radius of rMin = 80 * cm and is meant to only be used with the BECAL!" << endl; + } + PHG4TTLSubsystem *ttl; + ttl = new PHG4TTLSubsystem(name); + ttl->SetDetailed(false); + ttl->SuperDetector(name); + ttl->set_int_param("isForward", 0); // + ttl->set_double_param("place_z", zOffset * cm); // + ttl->set_double_param("rMin", radius * cm); // + ttl->set_double_param("length", 2.0 * halflength * cm); + ttl->set_double_param("tSilicon", tSilicon); // + ttl->OverlapCheck(false); + + g4Reco->registerSubsystem(ttl); + + + if (TRACKING::FastKalmanFilter) + { + TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, + PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype, + tSilicon / sqrt(12.), // const float radres, + G4TTL::PositionResolution, // const float phires, + G4TTL::PositionResolution, // const float lonres, + 1, // const float eff, + 0); // const float noise + TRACKING::FastKalmanFilter->add_cylinder_state(name, radius); + + TRACKING::ProjectionNames.insert(name); + } + + return 0; +} + + + + +//-----------------------------------------------------------------------------------// +int make_barrel_layer_basic(string name, PHG4Reco *g4Reco, + double radius, double halflength, double tSilicon, double zOffset){ + //--------------------------------- //build barrel layer //--------------------------------- const int nSubLayer = 7; - const double cm = PHG4Sector::Sector_Geometry::Unit_cm(); - const double mm = .1 * cm; - const double um = 1e-3 * mm; - - const double halflength = 0.5 * (z_end - z_start); - const double z_center = 0.5 * (z_end + z_start); string layerName[nSubLayer] = {"SiliconSensor", "Metalconnection", "HDI", "Cooling", "Support1", "Support_Gap", "Support2"}; string material[nSubLayer] = {"G4_Si", "G4_Al", "G4_KAPTON", "G4_WATER", "G4_GRAPHITE", "G4_AIR", "G4_GRAPHITE"}; - double thickness[nSubLayer] = {tSilicon, 15 * um, 20 * um, 100 * um, + double thickness[nSubLayer] = {tSilicon , 15 * um, 20 * um, 100 * um, 50 * um, 1, 50 * um}; double max_bh_radius = 0.; - PHG4CylinderSubsystem *cyl; - // cout << "started to create cylinder layer: " << name << endl; - + PHG4CylinderSubsystem* cyl; +// cout << "started to create cylinder layer: " << name << endl; + double currRadius = radius; - // cout << currRadius << endl; - for (int l = 0; l < nSubLayer; l++) - { - // cout << name <<"_"<< layerName[l] << endl; - cyl = new PHG4CylinderSubsystem(name + "_" + layerName[l], l); +// cout << currRadius << endl; + for (int l = 0; l < nSubLayer; l++) { +// cout << name <<"_"<< layerName[l] << endl; + cyl = new PHG4CylinderSubsystem(name + "_" + layerName[l],l); cyl->SuperDetector(name); cyl->set_double_param("radius", currRadius); cyl->set_double_param("length", 2.0 * halflength); - cyl->set_double_param("place_z", z_center); cyl->set_string_param("material", material[l]); cyl->set_double_param("thickness", thickness[l]); + cyl->set_double_param("place_x", 0.); + cyl->set_double_param("place_y", 0.); + cyl->set_double_param("place_z", zOffset); + if (l == 0) cyl->SetActive(); //only the Silicon Sensor is active - cyl->OverlapCheck(OverlapCheck); + cyl->OverlapCheck(true); g4Reco->registerSubsystem(cyl); - currRadius = currRadius + thickness[l]; - // cout << currRadius << endl; + currRadius = currRadius+thickness[l]; +// cout << currRadius << endl; } + if (TRACKING::FastKalmanFilter) { TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + name, // const std::string& phg4hitsNames, PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype, tSilicon / sqrt(12.), // const float radres, - TTL::PositionResolution, // const float phires, - TTL::PositionResolution, // const float lonres, + G4TTL::PositionResolution, // const float phires, + G4TTL::PositionResolution, // const float lonres, 1, // const float eff, 0); // const float noise TRACKING::FastKalmanFilter->add_cylinder_state(name, radius); TRACKING::ProjectionNames.insert(name); } + return 0; } diff --git a/common/G4_TrackingSupport.C b/common/G4_TrackingSupport.C new file mode 100644 index 00000000..7013c7ff --- /dev/null +++ b/common/G4_TrackingSupport.C @@ -0,0 +1,253 @@ +#ifndef MACRO_G4TrackingService_C +#define MACRO_G4TrackingService_C + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +//ECCE Tracking Services +//Should be 5 barrels and 4 cones + +using namespace std; + +class ServiceProperties +{ + public: + ServiceProperties(); + + explicit ServiceProperties(const string &name, + const double &rad_len_copper, + const double &rad_len_aluminum, + const double &rad_len_water, + const double &rad_len_plastic, + const double &rad_len_carbon, + const double &rad_len_iron, + const double &z_south, + const double &z_north, + const double &r_south, + const double &r_north); + + virtual ~ServiceProperties(){}; + + const string get_name(); + const double get_rad_len_copper(); + const double get_rad_len_aluminum(); + const double get_rad_len_water(); + const double get_rad_len_plastic(); + const double get_rad_len_carbon(); + const double get_rad_len_iron(); + const double get_z_south(); + const double get_z_north(); + const double get_r_south(); + const double get_r_north(); + + private: + const string m_name = "service"; + const double m_rad_len_copper = 0.0; + const double m_rad_len_aluminum = 0.0; + const double m_rad_len_water = 0.0; + const double m_rad_len_plastic = 0.0; + const double m_rad_len_carbon = 0.0; + const double m_rad_len_iron = 0.0; + const double m_z_south = 0.0; + const double m_z_north = 0.0; + const double m_r_south = 0.0; + const double m_r_north = 0.0; +}; + +ServiceProperties::ServiceProperties(const string &name, + const double &rad_len_copper, + const double &rad_len_aluminum, + const double &rad_len_water, + const double &rad_len_plastic, + const double &rad_len_carbon, + const double &rad_len_iron, + const double &z_south, + const double &z_north, + const double &r_south, + const double &r_north) + : m_name(name) + , m_rad_len_copper(rad_len_copper) + , m_rad_len_aluminum(rad_len_aluminum) + , m_rad_len_water(rad_len_water) + , m_rad_len_plastic(rad_len_plastic) + , m_rad_len_carbon(rad_len_carbon) + , m_rad_len_iron(rad_len_iron) + , m_z_south(z_south) + , m_z_north(z_north) + , m_r_south(r_south) + , m_r_north(r_north) +{} + +const string ServiceProperties::get_name() { return m_name; } +const double ServiceProperties::get_rad_len_copper() { return m_rad_len_copper; } +const double ServiceProperties::get_rad_len_aluminum() { return m_rad_len_aluminum; } +const double ServiceProperties::get_rad_len_water() { return m_rad_len_water; } +const double ServiceProperties::get_rad_len_plastic() { return m_rad_len_plastic; } +const double ServiceProperties::get_rad_len_carbon() { return m_rad_len_carbon; } +const double ServiceProperties::get_rad_len_iron() { return m_rad_len_iron; } +const double ServiceProperties::get_z_south() { return m_z_south; } +const double ServiceProperties::get_z_north() { return m_z_north; } +const double ServiceProperties::get_r_south() { return m_r_south; } +const double ServiceProperties::get_r_north() { return m_r_north; } + +namespace Enable +{ + bool TrackingService = false; + bool TrackingService_ABSORBER = false; + bool TrackingService_OVERLAPCHECK = false; + int TrackingService_VERBOSITY = 0; + +} // namespace Enable + +namespace G4TrackingService +{ //List materials and radiation length in cm + const int nMaterials = 6; + pair materials[nMaterials] = { make_pair("G4_Cu", 1.436) + , make_pair("G4_Al", 8.897) + , make_pair("G4_WATER", 36.08) + , make_pair("G4_POLYETHYLENE", 50.31) + , make_pair("PEEK", 30.00) + , make_pair("G4_Fe", 1.757) }; + + double GlobalOffset = 0.0; + double ShellThickness = 0.3; //Thickness in cm + int subsysID = 0; +} // namespace G4TrackingService + +vector get_thickness(ServiceProperties *object) +{ + vector thickness = {(object->get_rad_len_copper()/100)*G4TrackingService::materials[0].second + ,(object->get_rad_len_aluminum()/100)*G4TrackingService::materials[1].second + ,(object->get_rad_len_water()/100)*G4TrackingService::materials[2].second + ,(object->get_rad_len_plastic()/100)*G4TrackingService::materials[3].second + ,(object->get_rad_len_carbon()/100)*G4TrackingService::materials[4].second + ,(object->get_rad_len_iron()/100)*G4TrackingService::materials[5].second}; + return thickness; +} + +void TrackingServiceInit() +{ +} + +double TrackingServiceCone(ServiceProperties *object, PHG4Reco* g4Reco, double radius) +{ + bool AbsorberActive = Enable::ABSORBER || Enable::TrackingService_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TrackingService_OVERLAPCHECK; + int verbosity = max(Enable::VERBOSITY, Enable::TrackingService_VERBOSITY); + + PHG4ConeSubsystem* cone; + + double innerRadiusSouth = object->get_r_south(); + double innerRadiusNorth = object->get_r_north(); + double length = abs(object->get_z_north() - object->get_z_south()); + vector thickness = get_thickness(object); + + for (int i = 0; i < G4TrackingService::nMaterials; ++i) + { + if (thickness[i] == 0) continue; + cone = new PHG4ConeSubsystem(object->get_name(), G4TrackingService::subsysID); + cone->Verbosity(verbosity); + cone->SetR1(innerRadiusSouth, innerRadiusSouth + thickness[i]); + cone->SetR2(innerRadiusNorth, innerRadiusNorth + thickness[i]); + cone->SetPlaceZ(object->get_z_south() + length/2 + G4TrackingService::GlobalOffset); + cone->SetZlength(length/2); + cone->SetMaterial(G4TrackingService::materials[i].first); + cone->SuperDetector("TrackingService"); + if (AbsorberActive) cone->SetActive(); + cone->OverlapCheck(OverlapCheck); + g4Reco->registerSubsystem(cone); + ++G4TrackingService::subsysID; + innerRadiusSouth += thickness[i]; + innerRadiusNorth += thickness[i]; + } + radius = max(innerRadiusSouth, innerRadiusNorth); + + return radius; +} + +double TrackingServiceCylinder(ServiceProperties *object, PHG4Reco* g4Reco, double radius) +{ + bool AbsorberActive = Enable::ABSORBER || Enable::TrackingService_ABSORBER; + bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TrackingService_OVERLAPCHECK; + int verbosity = max(Enable::VERBOSITY, Enable::TrackingService_VERBOSITY); + + PHG4CylinderSubsystem* cyl; + + double innerRadius = object->get_r_south(); + double length = abs(object->get_z_north() - object->get_z_south()); + vector thickness = get_thickness(object); + + for (int i = 0; i < G4TrackingService::nMaterials; ++i) + { + if (thickness[i] == 0) continue; + cyl = new PHG4CylinderSubsystem(object->get_name(), G4TrackingService::subsysID); + cyl->Verbosity(verbosity); + cyl->set_double_param("place_z", object->get_z_south() + length/2 + G4TrackingService::GlobalOffset); + cyl->set_double_param("radius", innerRadius); + cyl->set_double_param("length", length); + cyl->set_string_param("material", G4TrackingService::materials[i].first); + cyl->set_double_param("thickness", thickness[i]); + cyl->SuperDetector("TrackingService"); + if (AbsorberActive) cyl->SetActive(); + cyl->OverlapCheck(OverlapCheck); + g4Reco->registerSubsystem(cyl); + ++G4TrackingService::subsysID; + innerRadius += thickness[i]; + } + radius = innerRadius; + + return radius; +} + +double TrackingService(PHG4Reco* g4Reco, double radius) +{ + vector cylinders, cones; + + double shellX0 = 100*G4TrackingService::ShellThickness/G4TrackingService::materials[4].second; + + cylinders.push_back(new ServiceProperties("ETrackingCylinderService_1", 9, 0, 0.42, 0.32, shellX0, 0, -400, -310, 270, 0)); + cones.push_back(new ServiceProperties("ETrackingConeService_1", 9, 0, 0.42, 0.32, shellX0, 0, -310, -300, 270, 68)); + cylinders.push_back(new ServiceProperties("ETrackingCylinderService_2", 17, 0, 0.56, 0.64, shellX0, 0, -300, -200, 68, 0)); + cylinders.push_back(new ServiceProperties("ETrackingCylinderService_3", 15, 0, 0.56, 0.56, shellX0, 0, -200, -152, 68, 0)); + cones.push_back(new ServiceProperties("ETrackingConeService_2", 13, 0, 0.56, 0.48, shellX0, 0, -152, -132.1, 68, 63.2)); + cones.push_back(new ServiceProperties("ETrackingConeService_3", 13, 0, 0.56, 0.48, shellX0, 0, -132.1, -119.9, 63.2, 50.2)); + cones.push_back(new ServiceProperties("ETrackingConeService_4", 11, 0, 0.56, 0.40, shellX0, 0, -119.9, -107, 50.2, 47.1)); + cones.push_back(new ServiceProperties("ETrackingConeService_5", 9, 0, 0.56, 0.32, shellX0, 0, -107, -80, 47.1, 44)); + cones.push_back(new ServiceProperties("ETrackingConeService_6", 7.5, 0, 0.42, 0.24, shellX0, 0, -80, -56.9, 44, 42.7)); + cones.push_back(new ServiceProperties("ETrackingConeService_7", 6, 0, 0.14, 0.159, shellX0, 0, -56.9, -34, 42.7, 24)); + cones.push_back(new ServiceProperties("ETrackingConeService_8", 3, 0, 0, 0.079, shellX0, 0, -34, -18, 24, 7.5)); + + cylinders.push_back(new ServiceProperties("BTrackingCylinderService_1", 0, 0, 0, 0, 0.1, 0, -56.8, 56.8, 48, 0)); + cylinders.push_back(new ServiceProperties("BTrackingCylinderService_2", 0, 0, 0, 0, 0.1, 0, -33.6, 33.6, 24, 0)); + cylinders.push_back(new ServiceProperties("BTrackingCylinderService_3", 0, 0, 0, 0, shellX0, 0, -17.8, 17.8, 7.5, 0)); + + cones.push_back(new ServiceProperties("HTrackingConeService_1", 3, 0, 0, 0.079, shellX0, 0, 18, 34, 7.5, 24)); + cones.push_back(new ServiceProperties("HTrackingConeService_2", 6, 0, 0.14, 0.159, shellX0, 0, 34, 56.9, 24, 42.7)); + cones.push_back(new ServiceProperties("HTrackingConeService_3", 7.5, 0, 0.28, 0.24, shellX0, 0, 56.9, 79.9, 42.7, 44.1)); + cones.push_back(new ServiceProperties("HTrackingConeService_4", 9, 0, 0.42, 0.32, shellX0, 0, 79.9, 114.9, 44.1, 47.2)); + cones.push_back(new ServiceProperties("HTrackingConeService_5", 11, 0, 0.56, 0.40, shellX0, 0, 114.9, 124.9, 47.2, 47.5)); + cones.push_back(new ServiceProperties("HTrackingConeService_6", 13, 0, 0.70, 0.48, shellX0, 0, 124.9, 157.9, 47.5, 69.2)); + cylinders.push_back(new ServiceProperties("HTrackingCylinderService_1", 15, 0, 0.84, 0.56, shellX0, 0, 157.9, 173, 69.2, 0)); + cones.push_back(new ServiceProperties("HTrackingConeService_7", 13, 0, 0.70, 0.48, shellX0, 0, 173, 180, 69.2, 85)); + cones.push_back(new ServiceProperties("HTrackingConeService_8", 13, 0, 0.70, 0.48, shellX0, 0, 180, 195, 85, 100)); + + cylinders.push_back(new ServiceProperties("EEMCalSupport", 0, 0, 0, 0, 0, 171, -200, -197, 62, 0)); + + for (ServiceProperties *cylinder : cylinders) radius += TrackingServiceCylinder(cylinder, g4Reco, radius); + for (ServiceProperties *cone : cones) radius += TrackingServiceCone(cone, g4Reco, radius); + + return radius; +} + +#endif diff --git a/common/G4_Tracking_EIC.C b/common/G4_Tracking_EIC.C index bc05feb2..ffb2b3a5 100644 --- a/common/G4_Tracking_EIC.C +++ b/common/G4_Tracking_EIC.C @@ -27,16 +27,47 @@ namespace G4TRACKING { bool DISPLACED_VERTEX = false; bool PROJECTION_EEMC = false; + bool PROJECTION_EHCAL = false; bool PROJECTION_CEMC = false; + bool PROJECTION_BECAL = false; + bool PROJECTION_HCALOUT = false; bool PROJECTION_FEMC = false; bool PROJECTION_FHCAL = false; + bool PROJECTION_LFHCAL = false; } // namespace G4TRACKING //-----------------------------------------------------------------------------// void TrackingInit() { TRACKING::FastKalmanFilter = new PHG4TrackFastSim("PHG4TrackFastSim"); + TRACKING::FastKalmanFilterSiliconTrack = new PHG4TrackFastSim("FastKalmanFilterSiliconTrack"); + TRACKING::FastKalmanFilterInnerTrack = new PHG4TrackFastSim("FastKalmanFilterInnerTrack"); } + +void InitFastKalmanFilter(PHG4TrackFastSim *kalman_filter) +{ + // kalman_filter->Smearing(false); + if (G4TRACKING::DISPLACED_VERTEX) + { + // do not use truth vertex in the track fitting, + // which would lead to worse momentum resolution for prompt tracks + // but this allows displaced track analysis including DCA and vertex finding + kalman_filter->set_use_vertex_in_fitting(false); + kalman_filter->set_vertex_xy_resolution(0); // do not smear the vertex used in the built-in DCA calculation + kalman_filter->set_vertex_z_resolution(0); // do not smear the vertex used in the built-in DCA calculation + kalman_filter->enable_vertexing(true); // enable vertex finding and fitting + } + else + { + // constraint to a primary vertex and use it as part of the fitting level arm + kalman_filter->set_use_vertex_in_fitting(true); + kalman_filter->set_vertex_xy_resolution(50e-4); + kalman_filter->set_vertex_z_resolution(50e-4); + } + + kalman_filter->set_sub_top_node_name("TRACKS"); +} + //-----------------------------------------------------------------------------// void Tracking_Reco() { @@ -53,27 +84,8 @@ void Tracking_Reco() exit(1); } + InitFastKalmanFilter(TRACKING::FastKalmanFilter); TRACKING::FastKalmanFilter->Verbosity(verbosity); - // TRACKING::FastKalmanFilter->Smearing(false); - if (G4TRACKING::DISPLACED_VERTEX) - { - // do not use truth vertex in the track fitting, - // which would lead to worse momentum resolution for prompt tracks - // but this allows displaced track analysis including DCA and vertex finding - TRACKING::FastKalmanFilter->set_use_vertex_in_fitting(false); - TRACKING::FastKalmanFilter->set_vertex_xy_resolution(0); // do not smear the vertex used in the built-in DCA calculation - TRACKING::FastKalmanFilter->set_vertex_z_resolution(0); // do not smear the vertex used in the built-in DCA calculation - TRACKING::FastKalmanFilter->enable_vertexing(true); // enable vertex finding and fitting - } - else - { - // constraint to a primary vertex and use it as part of the fitting level arm - TRACKING::FastKalmanFilter->set_use_vertex_in_fitting(true); - TRACKING::FastKalmanFilter->set_vertex_xy_resolution(50e-4); - TRACKING::FastKalmanFilter->set_vertex_z_resolution(50e-4); - } - - TRACKING::FastKalmanFilter->set_sub_top_node_name("TRACKS"); TRACKING::FastKalmanFilter->set_trackmap_out_name(TRACKING::TrackNodeName); //------------------------- @@ -95,6 +107,14 @@ void Tracking_Reco() TRACKING::ProjectionNames.insert("FHCAL"); } //------------------------- + // LFHCAL + //------------------------- + if (Enable::LFHCAL && G4TRACKING::PROJECTION_LFHCAL) + { + TRACKING::FastKalmanFilter->add_state_name("LFHCAL"); + TRACKING::ProjectionNames.insert("LFHCAL"); + } + //------------------------- // CEMC //------------------------- if (Enable::CEMC && G4TRACKING::PROJECTION_CEMC) @@ -103,6 +123,22 @@ void Tracking_Reco() TRACKING::ProjectionNames.insert("CEMC"); } //------------------------- + // BECAL + //------------------------- + if (Enable::BECAL && G4TRACKING::PROJECTION_BECAL) + { + TRACKING::FastKalmanFilter->add_state_name("BECAL"); + TRACKING::ProjectionNames.insert("BECAL"); + } + //------------------------- + // HCALOUT + //------------------------- + if (Enable::HCALOUT && G4TRACKING::PROJECTION_HCALOUT) + { + TRACKING::FastKalmanFilter->add_state_name("HCALOUT"); + TRACKING::ProjectionNames.insert("HCALOUT"); + } + //------------------------- // EEMC //------------------------- if (Enable::EEMC && G4TRACKING::PROJECTION_EEMC) @@ -110,8 +146,39 @@ void Tracking_Reco() TRACKING::FastKalmanFilter->add_state_name("EEMC"); TRACKING::ProjectionNames.insert("EEMC"); } + //------------------------- + // EHCAL + //------------------------- + if (Enable::EHCAL && G4TRACKING::PROJECTION_EHCAL) + { + TRACKING::FastKalmanFilter->add_state_name("EHCAL"); + TRACKING::ProjectionNames.insert("EHCAL"); + } se->registerSubsystem(TRACKING::FastKalmanFilter); + + // next, tracks with partial usage of the tracker stack + if (TRACKING::FastKalmanFilterInnerTrack == nullptr) + { + cout << __PRETTY_FUNCTION__ << " : missing the expected initialization for TRACKING::FastKalmanFilterInnerTrack." << endl; + exit(1); + } + InitFastKalmanFilter(TRACKING::FastKalmanFilterInnerTrack); + TRACKING::FastKalmanFilterInnerTrack->Verbosity(verbosity); + TRACKING::FastKalmanFilterInnerTrack->set_trackmap_out_name("InnerTrackMap"); + TRACKING::FastKalmanFilterInnerTrack->enable_vertexing(false); + se->registerSubsystem(TRACKING::FastKalmanFilterInnerTrack); + + if (TRACKING::FastKalmanFilterSiliconTrack == nullptr) + { + cout << __PRETTY_FUNCTION__ << " : missing the expected initialization for TRACKING::FastKalmanFilterSiliconTrack." << endl; + exit(1); + } + InitFastKalmanFilter(TRACKING::FastKalmanFilterSiliconTrack); + TRACKING::FastKalmanFilterSiliconTrack->Verbosity(verbosity); + TRACKING::FastKalmanFilterSiliconTrack->set_trackmap_out_name("SiliconTrackMap"); + TRACKING::FastKalmanFilterSiliconTrack->enable_vertexing(false); + se->registerSubsystem(TRACKING::FastKalmanFilterSiliconTrack); return; } @@ -154,5 +221,18 @@ void Tracking_Eval(const std::string &outputfile) cout << "};" << endl; // override the TRACKING::ProjectionNames in eval macros se->registerSubsystem(fast_sim_eval); + + // now partial track fits + fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval_InnerTrackMap"); + fast_sim_eval->set_trackmapname("InnerTrackMap"); + fast_sim_eval->set_filename(outputfile + ".InnerTrackMap.root"); + fast_sim_eval->Verbosity(verbosity); + se->registerSubsystem(fast_sim_eval); + + fast_sim_eval = new PHG4TrackFastSimEval("FastTrackingEval_SiliconTrackMap"); + fast_sim_eval->set_trackmapname("SiliconTrackMap"); + fast_sim_eval->set_filename(outputfile + ".SiliconTrackMap.root"); + fast_sim_eval->Verbosity(verbosity); + se->registerSubsystem(fast_sim_eval); } #endif diff --git a/common/G4_dRICH.C b/common/G4_dRICH.C index d5158d0c..fe6119ca 100644 --- a/common/G4_dRICH.C +++ b/common/G4_dRICH.C @@ -10,17 +10,19 @@ #include -#include +#include #include #include R__LOAD_LIBRARY(libg4detectors.so) +R__LOAD_LIBRARY(libEICG4dRICH.so) namespace Enable { bool RICH = false; bool RICH_OVERLAPCHECK = false; + int RICH_VERBOSITY = 0; } // namespace Enable void RICHInit() @@ -37,53 +39,19 @@ void RICHInit() void RICHSetup(PHG4Reco* g4Reco) { bool OverlapCheck = Enable::OVERLAPCHECK || Enable::RICH_OVERLAPCHECK; + int verbosity = std::max(Enable::VERBOSITY, Enable::RICH_VERBOSITY); - double z = 185; - double dz = 0; + double z = 185; //Start of dRICH + double dz = 100; //Length of dRICH - dz = 20; - PHG4ConeSubsystem* coneAeroGel = new PHG4ConeSubsystem("dRICHconeAeroGel", 0); - coneAeroGel->SetR1(11, 110); - coneAeroGel->SetR2(11, 120); - coneAeroGel->SetZlength(dz * 0.5); - coneAeroGel->SetPlaceZ(z + dz * 0.5); - coneAeroGel->set_string_param("material", "ePHENIX_AeroGel"); - coneAeroGel->set_color(0.5, 0.5, 0.8); - coneAeroGel->SetActive(); - coneAeroGel->SuperDetector("RICH"); - coneAeroGel->OverlapCheck(OverlapCheck); - g4Reco->registerSubsystem(coneAeroGel); - - z += dz; - dz = 5; - PHG4ConeSubsystem* coneGas = new PHG4ConeSubsystem("dRICHconeGas", 1); - coneGas->SetR1(11, 120); - coneGas->SetR2(11, 170); // <- reduce from 210cm to avoid overlap with the HCal - coneGas->SetZlength(dz * 0.5); - coneGas->SetPlaceZ(z + dz * 0.5); - coneGas->set_string_param("material", "C4F10"); - coneGas->set_color(0.5, 0.5, 0.5); - coneGas->SetActive(); - coneGas->SuperDetector("RICH"); - coneGas->OverlapCheck(OverlapCheck); - g4Reco->registerSubsystem(coneGas); - - z += dz; - dz = 75; - coneGas = new PHG4ConeSubsystem("dRICHconeGas", 2); - coneGas->SetR1(11, 170); - coneGas->SetR2(15, 170); - coneGas->SetZlength(dz * 0.5); - coneGas->SetPlaceZ(z + dz * 0.5); - coneGas->set_string_param("material", "C4F10"); - coneGas->set_color(0.5, 0.5, 0.5); - coneGas->SetActive(); - coneGas->SuperDetector("RICH"); - coneGas->OverlapCheck(OverlapCheck); - g4Reco->registerSubsystem(coneGas); - - z += dz; + EICG4dRICHSubsystem *drichSubsys = new EICG4dRICHSubsystem("dRICh"); + drichSubsys->SetGeometryFile(string(getenv("CALIBRATIONROOT")) + "/dRICH/mapping/drich-g4model_v2.txt"); + drichSubsys->set_double_param("place_z", z + dz*0.5);// relative position to mother vol. + drichSubsys->OverlapCheck(OverlapCheck); + drichSubsys->Verbosity(verbosity); + drichSubsys->SetActive(); + g4Reco->registerSubsystem(drichSubsys); if (TRACKING::FastKalmanFilter) { diff --git a/common/G4_hFarFwdBeamLine_EIC.C b/common/G4_hFarFwdBeamLine_EIC.C index cb171bf4..53023726 100644 --- a/common/G4_hFarFwdBeamLine_EIC.C +++ b/common/G4_hFarFwdBeamLine_EIC.C @@ -12,10 +12,14 @@ #include #include +#include + #include #include +#include + R__LOAD_LIBRARY(libg4detectors.so) float PosFlip(float pos); @@ -42,6 +46,9 @@ namespace Enable bool HFARFWD_VIRTUAL_DETECTORS_IP8 = false; float HFARFWD_ION_ENERGY = 0; + + bool FFR_EVAL = false; + } // namespace Enable namespace hFarFwdBeamLine @@ -81,7 +88,7 @@ void hFarFwdBeamLineInit() { hFarFwdBeamLine::enclosure_z_max = 4500.; BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, hFarFwdBeamLine::starting_z); - hFarFwdBeamLine::enclosure_r_max = 200.; + hFarFwdBeamLine::enclosure_r_max = 200.; } hFarFwdBeamLine::enclosure_center = 0.5 * (hFarFwdBeamLine::starting_z + hFarFwdBeamLine::enclosure_z_max); @@ -103,13 +110,17 @@ void hFarFwdDefineMagnets(PHG4Reco *g4Reco) hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_string_param("material", "G4_Galactic"); hFarFwdBeamLine::hFarFwdBeamLineEnclosure->set_color(.5, .5, .5, 0.2); hFarFwdBeamLine::hFarFwdBeamLineEnclosure->OverlapCheck(overlapCheck); - if (verbosity) - hFarFwdBeamLine::hFarFwdBeamLineEnclosure->Verbosity(verbosity); + if (verbosity) hFarFwdBeamLine::hFarFwdBeamLineEnclosure->Verbosity(verbosity); g4Reco->registerSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); + if (verbosity > 0) + { + std::cout << "hFarFwdBeamLine::hFarFwdBeamLineEnclosure CanBeMotherSubsystem = " << hFarFwdBeamLine::hFarFwdBeamLineEnclosure->CanBeMotherSubsystem() << std::endl; + } + string magFile; if (Enable::HFARFWD_MAGNETS_IP6) - magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farFwdBeamLineMagnets.dat"; + magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farFwdBeamLineMagnets_v2.0.dat"; else if (Enable::HFARFWD_MAGNETS_IP8) magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip8_35mrad_h_farFwdBeamLineMagnets.dat"; else @@ -153,7 +164,7 @@ void hFarFwdDefineMagnets(PHG4Reco *g4Reco) double fieldgradient; if (!(iss >> magname >> x >> y >> z >> inner_radius_zin >> inner_radius_zout >> outer_magnet_diameter >> length >> angle >> dipole_field_x >> fieldgradient)) { - cout << "coud not decode " << line << endl; + cout << "could not decode " << line << endl; gSystem->Exit(1); } else @@ -247,8 +258,7 @@ void hFarFwdDefineMagnets(PHG4Reco *g4Reco) } bl->OverlapCheck(overlapCheck); bl->SuperDetector("BEAMLINEMAGNET"); - if (verbosity) - bl->Verbosity(verbosity); + if (verbosity) bl->Verbosity(verbosity); g4Reco->registerSubsystem(bl); // rag the B0 magnet @@ -281,10 +291,10 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) auto *detZDCsurrogate = new PHG4BlockSubsystem("zdcTruth"); const double detZDCsurrogate_size_z = 0.1; detZDCsurrogate->SuperDetector("ZDCsurrogate"); - detZDCsurrogate->set_double_param("place_x", PosFlip(96.24)); + detZDCsurrogate->set_double_param("place_x", PosFlip(-96.24)); detZDCsurrogate->set_double_param("place_y", 0); - detZDCsurrogate->set_double_param("place_z", 3750 - hFarFwdBeamLine::enclosure_center); - detZDCsurrogate->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg())); + detZDCsurrogate->set_double_param("place_z", 3700 - hFarFwdBeamLine::enclosure_center); + detZDCsurrogate->set_double_param("rot_y", AngleFlip(0.025 * TMath::RadToDeg())); detZDCsurrogate->set_double_param("size_x", 60); detZDCsurrogate->set_double_param("size_y", 60); detZDCsurrogate->set_double_param("size_z", detZDCsurrogate_size_z); @@ -293,8 +303,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detZDCsurrogate->set_color(1, 0, 0, 0.5); detZDCsurrogate->OverlapCheck(overlapCheck); if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole(); - if (verbosity) - detZDCsurrogate->Verbosity(verbosity); + if (verbosity) detZDCsurrogate->Verbosity(verbosity); detZDCsurrogate->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); g4Reco->registerSubsystem(detZDCsurrogate); @@ -302,9 +311,9 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) { EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC"); detZDC->SetActive(); - detZDC->set_double_param("place_z", 3750. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center); - detZDC->set_double_param("place_x", PosFlip(96.24)); - detZDC->set_double_param("rot_y", AngleFlip(-0.025)); + detZDC->set_double_param("place_z", 3700. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center); + detZDC->set_double_param("place_x", PosFlip(-96.24)); + detZDC->set_double_param("rot_y", AngleFlip(0.025)); detZDC->OverlapCheck(overlapCheck); detZDC->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); g4Reco->registerSubsystem(detZDC); @@ -312,7 +321,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) const int offMomDetNr = 2; const double om_zCent[offMomDetNr] = {3450, 3650}; - const double om_xCent[offMomDetNr] = {162, 171}; + const double om_xCent[offMomDetNr] = {-162, -171}; for (int i = 0; i < offMomDetNr; i++) { auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i), i); @@ -320,14 +329,13 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detOM->set_double_param("place_x", PosFlip(om_xCent[i])); detOM->set_double_param("place_y", 0); detOM->set_double_param("place_z", om_zCent[i] - hFarFwdBeamLine::enclosure_center); - detOM->set_double_param("rot_y", AngleFlip(-0.045 * TMath::RadToDeg())); + detOM->set_double_param("rot_y", AngleFlip(0.045 * TMath::RadToDeg())); detOM->set_double_param("size_x", 50); detOM->set_double_param("size_y", 35); detOM->set_double_param("size_z", 0.03); detOM->set_string_param("material", "G4_Si"); detOM->SetActive(); - if (verbosity) - detOM->Verbosity(verbosity); + if (verbosity) detOM->Verbosity(verbosity); detOM->OverlapCheck(overlapCheck); detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); g4Reco->registerSubsystem(detOM); @@ -335,7 +343,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) const int rpDetNr = 2; const double rp_zCent[rpDetNr] = {2600, 2800}; - const double rp_xCent[rpDetNr] = {84.49, 93.59}; + const double rp_xCent[rpDetNr] = {-83.22, -92.20}; for (int i = 0; i < rpDetNr; i++) { ////********************* @@ -358,12 +366,12 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) //// Disk design //// 50 cm in x - auto *detRP = new PHG4CylinderSubsystem(Form("rpTruth_%d", i), i); + auto *detRP = new PHG4CylinderSubsystem(Form("rpTruth_%d", i),i); detRP->SuperDetector("rpTruth"); detRP->set_double_param("place_x", PosFlip(rp_xCent[i])); detRP->set_double_param("place_y", 0); detRP->set_double_param("place_z", rp_zCent[i] - hFarFwdBeamLine::enclosure_center); - detRP->set_double_param("rot_y", AngleFlip(-0.025 * TMath::RadToDeg())); + detRP->set_double_param("rot_y", AngleFlip(0.047 * TMath::RadToDeg())); detRP->set_double_param("radius", 0); detRP->set_double_param("thickness", 25); // This is intentionally made large 25cm radius detRP->set_double_param("length", 0.03); @@ -372,28 +380,30 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); detRP->SetActive(); - if (verbosity) - detRP->Verbosity(verbosity); + if (verbosity) detRP->Verbosity(verbosity); g4Reco->registerSubsystem(detRP); } const int b0DetNr = 4; const double b0Mag_zCent = 590; const double b0Mag_zLen = 120; + + for (int i = 0; i < b0DetNr; i++) { - auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), i); + auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i),i); detB0->SuperDetector("b0Truth"); detB0->set_double_param("radius", 0); detB0->set_double_param("thickness", 20); detB0->set_double_param("length", 0.1); detB0->set_string_param("material", "G4_Si"); - detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); // relative to B0 magnet + detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); detB0->SetActive(true); - if (verbosity) - detB0->Verbosity(verbosity); + if (verbosity) detB0->Verbosity(verbosity); detB0->OverlapCheck(overlapCheck); + detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); + g4Reco->registerSubsystem(detB0); } } @@ -415,31 +425,31 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY); -// const int offMomDetNr = 3; -// const double om_zCent[offMomDetNr] = {4250, 4400, 4550}; -// const double om_xCent[offMomDetNr] = {100, 100, 100}; - + //const int offMomDetNr = 3; + //const double om_xCent[offMomDetNr] = {100, 100, 100}; + //const double om_zCent[offMomDetNr] = {4250, 4400, 4550}; + const int offMomDetNr = 2; const double om_xCent[offMomDetNr] = {46, 49}; const double om_zCent[offMomDetNr] = {3250, 3450}; for (int i = 0; i < offMomDetNr; i++) { - auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i)); + auto *detOM = new PHG4BlockSubsystem(Form("offMomTruth_%d", i),i); + detOM->SuperDetector(Form("SDoffMomTruth_%d",i)); detOM->set_double_param("place_x", PosFlip(om_xCent[i])); detOM->set_double_param("place_y", 0); - detOM->set_double_param("place_z", om_zCent[i]); + detOM->set_double_param("place_z", om_zCent[i] - hFarFwdBeamLine::enclosure_center); detOM->set_double_param("rot_y", AngleFlip(-0.045 * TMath::RadToDeg())); detOM->set_double_param("size_x", 40); // Original design specification detOM->set_double_param("size_y", 35); // Original design specification -// detOM->set_double_param("size_x", 100); -// detOM->set_double_param("size_y", 100); detOM->set_double_param("size_z", 0.03); detOM->set_string_param("material", "G4_Si"); + detOM->OverlapCheck(overlapCheck); + detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); detOM->SetActive(); detOM->set_color(0, 0, 1, 0.5); - if (verbosity) - detOM->Verbosity(verbosity); + if (verbosity) detOM->Verbosity(verbosity); g4Reco->registerSubsystem(detOM); } @@ -448,7 +458,7 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) detZDCsurrogate->SuperDetector("ZDCsurrogate"); detZDCsurrogate->set_double_param("place_x", PosFlip(120)); detZDCsurrogate->set_double_param("place_y", 0); - detZDCsurrogate->set_double_param("place_z", 3350); + detZDCsurrogate->set_double_param("place_z", 3350 - hFarFwdBeamLine::enclosure_center); detZDCsurrogate->set_double_param("rot_y", AngleFlip(-0.035 * TMath::RadToDeg())); detZDCsurrogate->set_double_param("size_x", 60); detZDCsurrogate->set_double_param("size_y", 60); @@ -458,45 +468,50 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) detZDCsurrogate->OverlapCheck(overlapCheck); detZDCsurrogate->set_color(1, 0, 0, 0.5); if (!Enable::ZDC_DISABLE_BLACKHOLE) detZDCsurrogate->BlackHole(); - if (verbosity) - detZDCsurrogate->Verbosity(verbosity); + if (verbosity) detZDCsurrogate->Verbosity(verbosity); g4Reco->registerSubsystem(detZDCsurrogate); - EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC"); - detZDC->SetActive(); - detZDC->set_double_param("place_z", 3350. + detZDCsurrogate_size_z); - detZDC->set_double_param("place_x", PosFlip(120.0)); - detZDC->set_double_param("rot_y", AngleFlip(-0.035)); - detZDC->OverlapCheck(overlapCheck); - g4Reco->registerSubsystem(detZDC); + + + if (Enable::ZDC_DISABLE_BLACKHOLE) + { + + EICG4ZDCSubsystem *detZDC = new EICG4ZDCSubsystem("EICG4ZDC"); + detZDC->SetActive(); + detZDC->set_double_param("place_z", 3350. + detZDCsurrogate_size_z - hFarFwdBeamLine::enclosure_center); + detZDC->set_double_param("place_x", PosFlip(120.0)); + detZDC->set_double_param("rot_y", AngleFlip(-0.035)); + detZDC->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); + detZDC->OverlapCheck(overlapCheck); + g4Reco->registerSubsystem(detZDC); + + } //------------------ // Roman pot set #1 const int rpDetNr = 2; -// const double rp_zCent[rpDetNr] = {2200, 2500, 2800, 3100}; -// const double rp_xCent[rpDetNr] = {75, 75, 75, 75}; + //const double rp_xCent[rpDetNr] = {75, 75, 75, 75}; + //const double rp_zCent[rpDetNr] = {2200, 2500, 2800, 3100}; const double rp_xCent[rpDetNr] = {75.6, 78.15}; const double rp_zCent[rpDetNr] = {2600, 2800}; for (int i = 0; i < rpDetNr; i++) { - auto *detRP = new PHG4BlockSubsystem(Form("rpTruth_%d", i)); - // detRP->SuperDetector("RomanPots"); - detRP->SuperDetector(Form("RomanPots_%d", i)); + auto *detRP = new PHG4BlockSubsystem(Form("RomanPots_%d", i),i); + detRP->SuperDetector(Form("SDRomanPots_%d", i)); detRP->set_double_param("place_x", PosFlip(rp_xCent[i])); detRP->set_double_param("place_y", 0); - detRP->set_double_param("place_z", rp_zCent[i]); + detRP->set_double_param("place_z", rp_zCent[i] - hFarFwdBeamLine::enclosure_center); detRP->set_double_param("rot_y", AngleFlip(-0.035 * TMath::RadToDeg())); detRP->set_double_param("size_x", 25); // Original design specification detRP->set_double_param("size_y", 20); // Original design specification -// detRP->set_double_param("size_x", 100); -// detRP->set_double_param("size_y", 100); detRP->set_double_param("size_z", 0.03); detRP->set_string_param("material", "G4_Si"); + detRP->OverlapCheck(overlapCheck); + detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); detRP->SetActive(); - if (verbosity) - detRP->Verbosity(verbosity); + if (verbosity) detRP->Verbosity(verbosity); g4Reco->registerSubsystem(detRP); } @@ -509,42 +524,45 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) for (int i = 0; i < rp2ndDetNr; i++) { - auto *detRP_2nd = new PHG4BlockSubsystem(Form("rp2ndTruth_%d", i)); - // detRP_2nd->SuperDetector("RomanPots"); - detRP_2nd->SuperDetector(Form("RomanPots_2nd_%d", i)); + auto *detRP_2nd = new PHG4BlockSubsystem(Form("RomanPots_2nd_%d", i),i); + detRP_2nd->SuperDetector(Form("SDRomanPots_2nd_%d", i)); detRP_2nd->set_double_param("place_x", PosFlip(rp_2nd_xCent[i])); detRP_2nd->set_double_param("place_y", 0); - detRP_2nd->set_double_param("place_z", rp_2nd_zCent[i]); + detRP_2nd->set_double_param("place_z", rp_2nd_zCent[i] - hFarFwdBeamLine::enclosure_center); detRP_2nd->set_double_param("rot_y", AngleFlip(-0.029 * TMath::RadToDeg())); -// detRP_2nd->set_double_param("size_x", 10); // Original design specification -// detRP_2nd->set_double_param("size_y", 5); // Original design specification detRP_2nd->set_double_param("size_x", 25); detRP_2nd->set_double_param("size_y", 20); detRP_2nd->set_double_param("size_z", 0.03); detRP_2nd->set_string_param("material", "G4_Si"); + detRP_2nd->OverlapCheck(overlapCheck); + detRP_2nd->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); detRP_2nd->SetActive(); - if (verbosity) - detRP_2nd->Verbosity(verbosity); + if (verbosity) detRP_2nd->Verbosity(verbosity); g4Reco->registerSubsystem(detRP_2nd); } + if (verbosity > 0) + { + std::cout << "B0Magnet can be mother = " << hFarFwdBeamLine::B0Magnet->CanBeMotherSubsystem() << std::endl; + } + const int b0DetNr = 4; const double b0Mag_zCent = 610; const double b0Mag_zLen = 120; for (int i = 0; i < b0DetNr; i++) { - auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i), 0); - //detB0->SuperDetector("B0detectors"); + auto *detB0 = new PHG4CylinderSubsystem(Form("b0Truth_%d", i),i); + detB0->SuperDetector("b0truth"); detB0->set_double_param("radius", 0); detB0->set_double_param("thickness", 20); detB0->set_double_param("length", 0.1); detB0->set_string_param("material", "G4_Si"); - detB0->set_double_param("place_x", PosFlip(21.2)); detB0->set_double_param("place_y", 0); - detB0->set_double_param("place_z", (b0Mag_zCent - b0Mag_zLen / 2) + b0Mag_zLen / (b0DetNr - 1) * i); + detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); + detB0->OverlapCheck(overlapCheck); + detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); detB0->SetActive(true); - if (verbosity) - detB0->Verbosity(verbosity); + if (verbosity) detB0->Verbosity(verbosity); g4Reco->registerSubsystem(detB0); } } @@ -600,7 +618,6 @@ void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco) pipe->set_double_param("place_y", qyC[i]); pipe->set_double_param("place_z", qzC[i]); pipe->SetActive(false); - // pipe->SetActive(true); g4Reco->registerSubsystem(pipe); } @@ -628,7 +645,6 @@ void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco) pipeZDC->set_double_param("place_y", 0); pipeZDC->set_double_param("place_z", 2041.59); pipeZDC->SetActive(false); - // pipeZDC->SetActive(true); g4Reco->registerSubsystem(pipeZDC); //Roman Pot pipe @@ -660,7 +676,7 @@ void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco) float PosFlip(float pos) { if(Enable::HFARFWD_MAGNETS_IP6) { - return -pos; + return pos; } else { return pos; } @@ -668,7 +684,7 @@ float PosFlip(float pos) { float AngleFlip(float angle){ if(Enable::HFARFWD_MAGNETS_IP6) { - return -angle; + return angle; } else { return angle; } @@ -676,13 +692,38 @@ float AngleFlip(float angle){ float MagFieldFlip(float Bfield){ if(Enable::HFARFWD_MAGNETS_IP6) { - return -Bfield; + return Bfield; } else { return Bfield; } } +//------------------------------------------ + +void FFR_Eval(const std::string &outputfile) +{ + + string ip_str; + + + if(Enable::IP6) { + ip_str = "IP6"; + } else { + ip_str = "IP8"; + } + + int verbosity = std::max(Enable::VERBOSITY, Enable::EEMC_VERBOSITY); + + Fun4AllServer *se = Fun4AllServer::instance(); + + FarForwardEvaluator *eval = new FarForwardEvaluator("FARFORWARDEVALUATOR", "FFR", outputfile.c_str(), ip_str); + + eval->Verbosity(verbosity); + se->registerSubsystem(eval); + + return; +} diff --git a/common/G4_mRICH.C b/common/G4_mRICH.C index 56c4a6e9..2928ed76 100644 --- a/common/G4_mRICH.C +++ b/common/G4_mRICH.C @@ -43,7 +43,7 @@ void mRICHSetup(PHG4Reco *g4Reco, const int detectorSetup = 1, //1: full setup; mRICH->set_int_param("detectorSetup", detectorSetup); mRICH->set_int_param("subsystemSetup", mRICHsystemSetup); mRICH->UseCalibFiles(PHG4DetectorSubsystem::xml); - mRICH->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/mRICH/Geometry/")); + mRICH->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/mRICH/GeometryV2/")); mRICH->OverlapCheck(OverlapCheck); // mRICH->Verbosity(5); diff --git a/common/G4_mRwell_EIC.C b/common/G4_mRwell_EIC.C index 154471c9..c0502e3a 100644 --- a/common/G4_mRwell_EIC.C +++ b/common/G4_mRwell_EIC.C @@ -29,10 +29,10 @@ namespace RWELL // const double nom_driftgap[RWELL::n_layer] = {0.4, 0.4}; // const double nom_length[RWELL::n_layer] = {300.0, 300.0}; - const int n_layer = 1; //tracker layers - const double nom_radius[RWELL::n_layer] = {69 - 1.6 - 2.6}; - const double nom_driftgap[RWELL::n_layer] = {0.4}; - const double nom_length[RWELL::n_layer] = {300.0}; + const int n_layer = 3; //tracker layers + const double nom_radius[RWELL::n_layer] = {44.2, 47.4, 67.4}; + const double nom_driftgap[RWELL::n_layer] = {0.4, 0.4, 0.4}; + const double nom_length[RWELL::n_layer] = {140, 150, 280.0}; } //namespace RWELL void RWellInit(int verbosity = 0) @@ -379,6 +379,23 @@ double RWellSetup(PHG4Reco* g4Reco, 55e-4, // const float lonres, 1, // const float eff, 0); // const float noise + if (TRACKING::FastKalmanFilterInnerTrack) + TRACKING::FastKalmanFilterInnerTrack->add_phg4hits(string("G4HIT_") + string(Form("RWELL_%d", ilyr)), // const std::string& phg4hitsNames, + PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype, + 1. / sqrt(12.), // const float radres, + 55e-4, // const float phires, + 55e-4, // const float lonres, + 1, // const float eff, + 0); // const float noise + // only for layers that is close to the silicon tracker system, use in FastKalmanFilterSiliconTrack + if (TRACKING::FastKalmanFilterSiliconTrack and RWELL::nom_radius[ilyr] < 50) + TRACKING::FastKalmanFilterSiliconTrack->add_phg4hits(string("G4HIT_") + string(Form("RWELL_%d", ilyr)), // const std::string& phg4hitsNames, + PHG4TrackFastSim::Cylinder, // const DETECTOR_TYPE phg4dettype, + 1. / sqrt(12.), // const float radres, + 55e-4, // const float phires, + 55e-4, // const float lonres, + 1, // const float eff, + 0); // const float noise } return radius; //cm } diff --git a/common/GlobalVariables.C b/common/GlobalVariables.C index 66f068ba..064d68cc 100644 --- a/common/GlobalVariables.C +++ b/common/GlobalVariables.C @@ -64,6 +64,10 @@ namespace TRACKING PHG4TrackFastSim * FastKalmanFilter(nullptr); + PHG4TrackFastSim * FastKalmanFilterSiliconTrack(nullptr); + + PHG4TrackFastSim * FastKalmanFilterInnerTrack(nullptr); + std::set ProjectionNames; } diff --git a/detectors/EICDetector/Fun4All_G4_EICDetector.C b/detectors/EICDetector/Fun4All_G4_EICDetector.C index e36b564c..f7b81eb9 100644 --- a/detectors/EICDetector/Fun4All_G4_EICDetector.C +++ b/detectors/EICDetector/Fun4All_G4_EICDetector.C @@ -265,7 +265,7 @@ int Fun4All_G4_EICDetector( // whether to simulate the Be section of the beam pipe Enable::PIPE = true; // If need to disable EIC beam pipe extension beyond the Be-section: - // G4PIPE::use_forward_pipes = false; + G4PIPE::use_forward_pipes = false; //EIC hadron far forward magnets and detectors. IP6 and IP8 are incompatible (pick either or); Enable::HFARFWD_MAGNETS = true; Enable::HFARFWD_VIRTUAL_DETECTORS = true; @@ -279,10 +279,12 @@ int Fun4All_G4_EICDetector( // Enable::BGEM = true; // not yet defined in this model Enable::RWELL = true; // barrel tracker + Enable::TrackingService = true; + // Enable::TrackingService_VERBOSITY = INT_MAX - 10; Enable::BARREL = true; // fst Enable::FST = true; - // G4FST::SETTING::SUPPORTCYL = false; // if want to disable support + G4FST::SETTING::SUPPORTCYL = false; // if want to disable support // TOFs Enable::FTTL = true; @@ -290,22 +292,25 @@ int Fun4All_G4_EICDetector( Enable::CTTL = true; Enable::TRACKING = true; - Enable::TRACKING_EVAL = Enable::TRACKING && false; - G4TRACKING::DISPLACED_VERTEX = false; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes + Enable::TRACKING_EVAL = Enable::TRACKING && true; + G4TRACKING::DISPLACED_VERTEX = true; // this option exclude vertex in the track fitting and use RAVE to reconstruct primary and 2ndary vertexes // projections to calorimeters G4TRACKING::PROJECTION_EEMC = true; + G4TRACKING::PROJECTION_BECAL = true; + G4TRACKING::PROJECTION_EHCAL = true; G4TRACKING::PROJECTION_CEMC = true; + G4TRACKING::PROJECTION_HCALOUT = true; G4TRACKING::PROJECTION_FEMC = true; G4TRACKING::PROJECTION_FHCAL = true; + G4TRACKING::PROJECTION_LFHCAL = true; - Enable::CEMC = true; - // Enable::CEMC_ABSORBER = true; - Enable::CEMC_CELL = Enable::CEMC && true; - Enable::CEMC_TOWER = Enable::CEMC_CELL && true; - Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true; - Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && false; + Enable::BECAL = true; + Enable::BECAL_CELL = Enable::BECAL && true; + Enable::BECAL_TOWER = Enable::BECAL_CELL && true; + Enable::BECAL_CLUSTER = Enable::BECAL_TOWER && true; + Enable::BECAL_EVAL = Enable::BECAL_CLUSTER && false; - Enable::HCALIN = true; + Enable::HCALIN = false; // Enable::HCALIN_ABSORBER = true; Enable::HCALIN_CELL = Enable::HCALIN && true; Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true; @@ -322,7 +327,7 @@ int Fun4All_G4_EICDetector( Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && false; // EICDetector geometry - barrel - Enable::DIRC = true; + Enable::DIRC = false; // EICDetector geometry - 'hadron' direction Enable::RICH = true; @@ -336,17 +341,32 @@ int Fun4All_G4_EICDetector( Enable::FEMC_CLUSTER = Enable::FEMC_TOWER && true; Enable::FEMC_EVAL = Enable::FEMC_CLUSTER && false; - Enable::FHCAL = true; - // Enable::FHCAL_ABSORBER = true; - Enable::FHCAL_TOWER = Enable::FHCAL && true; - Enable::FHCAL_CLUSTER = Enable::FHCAL_TOWER && true; - Enable::FHCAL_EVAL = Enable::FHCAL_CLUSTER && false; + //Enable::DRCALO = false; + Enable::DRCALO_CELL = Enable::DRCALO && true; + Enable::DRCALO_TOWER = Enable::DRCALO_CELL && true; + Enable::DRCALO_CLUSTER = Enable::DRCALO_TOWER && true; + Enable::DRCALO_EVAL = Enable::DRCALO_CLUSTER && false; + G4TTL::SETTING::optionDR = 1; + + Enable::LFHCAL = false; + G4LFHCAL::SETTING::longer = true; + G4LFHCAL::SETTING::asymmetric = true; + Enable::LFHCAL_ABSORBER = false; + Enable::LFHCAL_CELL = Enable::LFHCAL && true; + Enable::LFHCAL_TOWER = Enable::LFHCAL_CELL && true; + Enable::LFHCAL_CLUSTER = Enable::LFHCAL_TOWER && true; + Enable::LFHCAL_EVAL = Enable::LFHCAL_CLUSTER && false; // EICDetector geometry - 'electron' direction - Enable::EEMC = true; - Enable::EEMC_TOWER = Enable::EEMC && true; - Enable::EEMC_CLUSTER = Enable::EEMC_TOWER && true; - Enable::EEMC_EVAL = Enable::EEMC_CLUSTER && false; + Enable::EEMCH = true; + G4EEMCH::SETTING::USECEMCGeo = false; + G4EEMCH::SETTING::USEHYBRID = false; + Enable::EEMCH_TOWER = Enable::EEMCH && true; + Enable::EEMCH_CLUSTER = Enable::EEMCH_TOWER && true; + Enable::EEMCH_EVAL = Enable::EEMCH_CLUSTER && true; + G4TTL::SETTING::optionEEMCH = Enable::EEMCH && true; + G4TTL::SETTING::optionCEMC = false; + G4TTL::SETTING::optionGeo = 1; Enable::EHCAL = true; Enable::EHCAL_CELL = Enable::EHCAL && true; @@ -354,6 +374,8 @@ int Fun4All_G4_EICDetector( Enable::EHCAL_CLUSTER = Enable::EHCAL_TOWER && true; Enable::EHCAL_EVAL = Enable::EHCAL_CLUSTER && false; + Enable::FFR_EVAL = Enable::HFARFWD_MAGNETS && Enable::HFARFWD_VIRTUAL_DETECTORS && true; + Enable::PLUGDOOR = false; // Other options @@ -368,6 +390,10 @@ int Fun4All_G4_EICDetector( Enable::BLACKHOLE = true; //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits //BlackHoleGeometry::visible = true; + + // ZDC + // Enable::ZDC = true; + // Enable::ZDC_DISABLE_BLACKHOLE = true; // Enabling the event evaluator? Enable::EVENT_EVAL = false; @@ -447,12 +473,21 @@ int Fun4All_G4_EICDetector( if (Enable::FHCAL_TOWER) FHCAL_Towers(); if (Enable::FHCAL_CLUSTER) FHCAL_Clusters(); + if (Enable::DRCALO_TOWER) DRCALO_Towers(); + if (Enable::DRCALO_CLUSTER) DRCALO_Clusters(); + + if (Enable::LFHCAL_TOWER) LFHCAL_Towers(); + if (Enable::LFHCAL_CLUSTER) LFHCAL_Clusters(); + if (Enable::EEMC_TOWER) EEMC_Towers(); if (Enable::EEMC_CLUSTER) EEMC_Clusters(); if (Enable::EHCAL_TOWER) EHCAL_Towers(); if (Enable::EHCAL_CLUSTER) EHCAL_Clusters(); + if (Enable::BECAL_TOWER) BECAL_Towers(); + if (Enable::BECAL_CLUSTER) BECAL_Clusters(); + if (Enable::DSTOUT_COMPRESS) ShowerCompress(); //-------------- @@ -510,10 +545,13 @@ int Fun4All_G4_EICDetector( if (Enable::EEMC_EVAL) EEMC_Eval(outputroot + "_g4eemc_eval.root"); + if (Enable::FFR_EVAL) FFR_Eval(outputroot + "_g4ffr_eval.root"); + if (Enable::FWDJETS_EVAL) Jet_FwdEval(); if (Enable::USER) UserAnalysisInit(); + //-------------- // Set up Input Managers //-------------- diff --git a/detectors/EICDetector/G4Setup_EICDetector.C b/detectors/EICDetector/G4Setup_EICDetector.C index 4dcc11b4..bb896f24 100644 --- a/detectors/EICDetector/G4Setup_EICDetector.C +++ b/detectors/EICDetector/G4Setup_EICDetector.C @@ -5,33 +5,38 @@ #include +//#include +#include +#include #include #include +#include +#include #include +#include #include #include #include +#include #include #include +#include +#include #include +#include +#include #include +#include +#include #include #include #include -#include - -#include -#include -#include -#include -#include #include +#include +#include #include #include -#include -#include -#include #include @@ -65,32 +70,53 @@ void G4Init() gSystem->Exit(1); } + if(Enable::EEMC and Enable::EEMCH) + { + cout << "Can not enable EEMC and EEMCH at the same time!" << endl; + gSystem->Exit(1); + } + // load detector/material macros and execute Init() function if (Enable::PIPE) PipeInit(); - if (Enable::HFARFWD_MAGNETS) hFarFwdBeamLineInit(); - if (Enable::HFARFWD_MAGNETS) hFarBwdBeamLineInit(); if (Enable::PLUGDOOR) PlugDoorInit(); - if (Enable::TRACKING) TrackingInit(); - if (Enable::EGEM) EGEM_Init(); - if (Enable::FGEM) FGEM_Init(); - if (Enable::RWELL) RWellInit(); - if (Enable::FST) FST_Init(); - if (Enable::FTTL || Enable::ETTL || Enable::CTTL) TTL_Init(); - if (Enable::BARREL) BarrelInit(); + //Farforward/backward + if (Enable::HFARFWD_MAGNETS) hFarBwdBeamLineInit(); //Shouldnt this be far backward enables + if (Enable::HFARFWD_MAGNETS) hFarFwdBeamLineInit(); + + //Barrel + if (Enable::TrackingService) TrackingServiceInit(); + if (Enable::BARREL) BarrelInit(); + if (Enable::RWELL) RWellInit(); if (Enable::CEMC) CEmcInit(72); // make it 2*2*2*3*3 so we can try other combinations + if (Enable::BECAL) BECALInit(); if (Enable::HCALIN) HCalInnerInit(1); if (Enable::MAGNET) MagnetInit(); MagnetFieldInit(); // We want the field - even if the magnet volume is disabled if (Enable::HCALOUT) HCalOuterInit(); + if (Enable::DIRC) DIRCInit(); + + //Forward + if (Enable::FGEM) FGEM_Init(); if (Enable::FEMC) FEMCInit(); + if (Enable::DRCALO) DRCALOInit(); if (Enable::FHCAL) FHCALInit(); + if (Enable::LFHCAL) LFHCALInit(); + if (Enable::RICH) RICHInit(); + + //Backward + if (Enable::EGEM) EGEM_Init(); if (Enable::EEMC) EEMCInit(); + if (Enable::EEMCH) EEMCHInit(); if (Enable::EHCAL) EHCALInit(); - if (Enable::DIRC) DIRCInit(); - if (Enable::RICH) RICHInit(); if (Enable::mRICH) mRICHInit(); + + //Combined + if (Enable::FST) FST_Init(); + if (Enable::FTTL || Enable::ETTL || Enable::CTTL) TTL_Init(); + + if (Enable::USER) UserInit(); if (Enable::BLACKHOLE) BlackHoleInit(); } @@ -143,6 +169,7 @@ int G4Setup() double radius = 0.; if (Enable::PIPE) radius = Pipe(g4Reco, radius); + // Far Forward Region if (Enable::HFARFWD_MAGNETS_IP6 || Enable::HFARFWD_MAGNETS_IP8) hFarFwdDefineMagnets(g4Reco); if (Enable::HFARFWD_VIRTUAL_DETECTORS_IP6) hFarFwdDefineDetectorsIP6(g4Reco); @@ -153,28 +180,35 @@ int G4Setup() if (Enable::HFARBWD_VIRTUAL_DETECTORS_IP6) hFarBwdDefineDetectorsIP6(g4Reco); if (Enable::HFARBWD_VIRTUAL_DETECTORS_IP8) hFarBwdDefineDetectorsIP8(g4Reco); - if (Enable::EGEM) EGEMSetup(g4Reco); - if (Enable::FGEM) FGEMSetup(g4Reco); + //Barrel + if (Enable::TrackingService) TrackingService(g4Reco, radius); + if (Enable::RWELL) RWellSetup(g4Reco); if (Enable::FST) FSTSetup(g4Reco); - if (Enable::FTTL) FTTLSetup(g4Reco); - if (Enable::ETTL) ETTLSetup(g4Reco); if (Enable::CTTL) CTTLSetup(g4Reco); if (Enable::BARREL) Barrel(g4Reco); if (Enable::CEMC) radius = CEmc(g4Reco, radius); + if (Enable::BECAL) BECALSetup(g4Reco); if (Enable::HCALIN) radius = HCalInner(g4Reco, radius, 4); if (Enable::MAGNET) radius = Magnet(g4Reco, radius); if (Enable::HCALOUT) radius = HCalOuter(g4Reco, radius, 4); + if (Enable::DIRC) DIRCSetup(g4Reco); + + //Forward + if (Enable::FGEM) FGEMSetup(g4Reco); + if (Enable::FTTL) FTTLSetup(g4Reco); if (Enable::FEMC) FEMCSetup(g4Reco); + if (Enable::DRCALO) DRCALOSetup(g4Reco); if (Enable::FHCAL) FHCALSetup(g4Reco); + if (Enable::LFHCAL) LFHCALSetup(g4Reco); + if (Enable::RICH) RICHSetup(g4Reco); + + //Backward + if (Enable::ETTL) ETTLSetup(g4Reco); + if (Enable::EGEM) EGEMSetup(g4Reco); if (Enable::EEMC) EEMCSetup(g4Reco); + if (Enable::EEMCH) EEMCHSetup(g4Reco); if (Enable::EHCAL) EHCALSetup(g4Reco); - - //---------------------------------------- - // PID - - if (Enable::DIRC) DIRCSetup(g4Reco); - if (Enable::RICH) RICHSetup(g4Reco); if (Enable::mRICH) mRICHSetup(g4Reco); //---------------------------------------- @@ -211,53 +245,83 @@ void ShowerCompress() // compress->AddHitContainer("G4HIT_RomanPots"); // compress->AddHitContainer("G4HIT_B0detector"); compress->AddHitContainer("G4HIT_FIELDCAGE"); + compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS"); compress->AddHitContainer("G4HIT_CEMC"); compress->AddHitContainer("G4HIT_ABSORBER_CEMC"); compress->AddHitContainer("G4HIT_CEMC_SPT"); - compress->AddHitContainer("G4HIT_ABSORBER_HCALIN"); - compress->AddHitContainer("G4HIT_HCALIN"); - compress->AddHitContainer("G4HIT_HCALIN_SPT"); - compress->AddHitContainer("G4HIT_MAGNET"); - compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT"); - compress->AddHitContainer("G4HIT_HCALOUT"); - compress->AddHitContainer("G4HIT_BH_1"); - compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS"); - compress->AddHitContainer("G4HIT_BH_FORWARD_NEG"); compress->AddCellContainer("G4CELL_CEMC"); - compress->AddCellContainer("G4CELL_HCALIN"); - compress->AddCellContainer("G4CELL_HCALOUT"); compress->AddTowerContainer("TOWER_SIM_CEMC"); compress->AddTowerContainer("TOWER_RAW_CEMC"); compress->AddTowerContainer("TOWER_CALIB_CEMC"); + + compress->AddHitContainer("G4HIT_BECAL"); + compress->AddHitContainer("G4HIT_ABSORBER_BECAL"); + compress->AddCellContainer("G4CELL_BECAL"); + compress->AddTowerContainer("TOWER_SIM_BECAL"); + compress->AddTowerContainer("TOWER_RAW_BECAL"); + compress->AddTowerContainer("TOWER_CALIB_BECAL"); + + compress->AddHitContainer("G4HIT_HCALIN"); + compress->AddHitContainer("G4HIT_ABSORBER_HCALIN"); + compress->AddHitContainer("G4HIT_HCALIN_SPT"); + compress->AddCellContainer("G4CELL_HCALIN"); compress->AddTowerContainer("TOWER_SIM_HCALIN"); compress->AddTowerContainer("TOWER_RAW_HCALIN"); compress->AddTowerContainer("TOWER_CALIB_HCALIN"); + + compress->AddHitContainer("G4HIT_MAGNET"); + + compress->AddHitContainer("G4HIT_HCALOUT"); + compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT"); + compress->AddCellContainer("G4CELL_HCALOUT"); compress->AddTowerContainer("TOWER_SIM_HCALOUT"); compress->AddTowerContainer("TOWER_RAW_HCALOUT"); compress->AddTowerContainer("TOWER_CALIB_HCALOUT"); + compress->AddHitContainer("G4HIT_BH_1"); + compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS"); + compress->AddHitContainer("G4HIT_BH_FORWARD_NEG"); + compress->AddHitContainer("G4HIT_FEMC"); compress->AddHitContainer("G4HIT_ABSORBER_FEMC"); - compress->AddHitContainer("G4HIT_FHCAL"); - compress->AddHitContainer("G4HIT_ABSORBER_FHCAL"); - compress->AddHitContainer("G4HIT_EHCAL"); - compress->AddHitContainer("G4HIT_ABSORBER_EHCAL"); compress->AddCellContainer("G4CELL_FEMC"); - compress->AddCellContainer("G4CELL_FHCAL"); compress->AddTowerContainer("TOWER_SIM_FEMC"); compress->AddTowerContainer("TOWER_RAW_FEMC"); compress->AddTowerContainer("TOWER_CALIB_FEMC"); + + compress->AddHitContainer("G4HIT_DRCALO"); + compress->AddHitContainer("G4HIT_ABSORBER_DRCALO"); + compress->AddCellContainer("G4CELL_DRCALO"); + compress->AddTowerContainer("TOWER_SIM_DRCALO"); + compress->AddTowerContainer("TOWER_RAW_DRCALO"); + compress->AddTowerContainer("TOWER_CALIB_DRCALO"); + + + compress->AddHitContainer("G4HIT_FHCAL"); + compress->AddHitContainer("G4HIT_ABSORBER_FHCAL"); + compress->AddCellContainer("G4CELL_FHCAL"); compress->AddTowerContainer("TOWER_SIM_FHCAL"); compress->AddTowerContainer("TOWER_RAW_FHCAL"); compress->AddTowerContainer("TOWER_CALIB_FHCAL"); + compress->AddHitContainer("G4HIT_LFHCAL"); + compress->AddHitContainer("G4HIT_ABSORBER_LFHCAL"); + compress->AddCellContainer("G4CELL_LFHCAL"); + compress->AddTowerContainer("TOWER_SIM_LFHCAL"); + compress->AddTowerContainer("TOWER_RAW_LFHCAL"); + compress->AddTowerContainer("TOWER_CALIB_LFHCAL"); + compress->AddHitContainer("G4HIT_EEMC"); + compress->AddHitContainer("G4HIT_EEMC_glass"); compress->AddHitContainer("G4HIT_ABSORBER_EEMC"); compress->AddCellContainer("G4CELL_EEMC"); compress->AddTowerContainer("TOWER_SIM_EEMC"); compress->AddTowerContainer("TOWER_RAW_EEMC"); compress->AddTowerContainer("TOWER_CALIB_EEMC"); + + compress->AddHitContainer("G4HIT_EHCAL"); + compress->AddHitContainer("G4HIT_ABSORBER_EHCAL"); compress->AddTowerContainer("TOWER_SIM_EHCAL"); compress->AddTowerContainer("TOWER_RAW_EHCAL"); compress->AddTowerContainer("TOWER_CALIB_EHCAL"); @@ -285,27 +349,36 @@ void DstCompress(Fun4AllDstOutputManager *out) out->StripNode("G4HIT_CEMC"); out->StripNode("G4HIT_ABSORBER_CEMC"); out->StripNode("G4HIT_CEMC_SPT"); + out->StripNode("G4CELL_CEMC"); + out->StripNode("G4HIT_BECAL"); + out->StripNode("G4CELL_BECAL"); + out->StripNode("G4HIT_ABSORBER_BECAL"); out->StripNode("G4HIT_ABSORBER_HCALIN"); out->StripNode("G4HIT_HCALIN"); out->StripNode("G4HIT_HCALIN_SPT"); + out->StripNode("G4CELL_HCALIN"); out->StripNode("G4HIT_MAGNET"); - out->StripNode("G4HIT_ABSORBER_HCALOUT"); out->StripNode("G4HIT_HCALOUT"); + out->StripNode("G4HIT_ABSORBER_HCALOUT"); + out->StripNode("G4CELL_HCALOUT"); out->StripNode("G4HIT_BH_1"); out->StripNode("G4HIT_BH_FORWARD_PLUS"); out->StripNode("G4HIT_BH_FORWARD_NEG"); - out->StripNode("G4CELL_CEMC"); - out->StripNode("G4CELL_HCALIN"); - out->StripNode("G4CELL_HCALOUT"); out->StripNode("G4HIT_FEMC"); out->StripNode("G4HIT_ABSORBER_FEMC"); out->StripNode("G4HIT_FHCAL"); out->StripNode("G4HIT_ABSORBER_FHCAL"); out->StripNode("G4CELL_FEMC"); + out->StripNode("G4HIT_DRCALO"); + out->StripNode("G4HIT_ABSORBER_DRCALO"); + out->StripNode("G4CELL_DRCALO"); out->StripNode("G4CELL_FHCAL"); - + out->StripNode("G4HIT_LFHCAL"); + out->StripNode("G4HIT_ABSORBER_LFHCAL"); + out->StripNode("G4CELL_LFHCAL"); out->StripNode("G4HIT_EEMC"); + out->StripNode("G4HIT_EEMC_glass"); out->StripNode("G4HIT_ABSORBER_EEMC"); out->StripNode("G4CELL_EEMC"); out->StripNode("G4HIT_EHCAL");