diff --git a/common/G4_B0Tracking_EIC.C b/common/G4_B0Tracking_EIC.C new file mode 100644 index 00000000..06e46ba6 --- /dev/null +++ b/common/G4_B0Tracking_EIC.C @@ -0,0 +1,155 @@ +#ifndef MACRO_G4B0TRACKINGEIC_C +#define MACRO_G4B0TRACKINGEIC_C + +#include + +#include + +#include +#include + +#include + +#include + +R__LOAD_LIBRARY(libtrack_reco.so) +R__LOAD_LIBRARY(libEICG4B0.so)// + +namespace Enable +{ + bool B0TRACKING = false; + bool B0TRACKING_EVAL = false; + int B0TRACKING_VERBOSITY = 0; +} // namespace Enable + +namespace G4B0TRACKING +{ + bool DISPLACED_VERTEX = false; + bool PROJECTION_B0 = false; + const double PositionResolution(30e-4); +} // namespace G4TRACKING + +//-----------------------------------------------------------------------------// +void B0TrackingInit() +{ + B0TRACKING::FastKalmanFilter = new B0TrackFastSim("B0TrackFastSim"); + B0TRACKING::FastKalmanFilterB0Track = new B0TrackFastSim("FastKalmanFilterB0Track"); //b0Truth +} + +void InitFastKalmanFilter(B0TrackFastSim *kalman_filter) +{ + // kalman_filter->Smearing(false); + if (G4B0TRACKING::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(G4B0TRACKING::PositionResolution); + kalman_filter->set_vertex_z_resolution(G4B0TRACKING::PositionResolution); + } + + kalman_filter->set_sub_top_node_name("TRACKS"); +} + +//-----------------------------------------------------------------------------// +void B0Tracking_Reco() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::B0TRACKING_VERBOSITY); + //--------------- + // Fun4All server + //--------------- + + Fun4AllServer *se = Fun4AllServer::instance(); + + if (B0TRACKING::FastKalmanFilter == nullptr) + { + cout << __PRETTY_FUNCTION__ << " : missing the expected initialization for B0TRACKING::FastKalmanFilter." << endl; + exit(1); + } + + InitFastKalmanFilter(B0TRACKING::FastKalmanFilter); + B0TRACKING::FastKalmanFilter->Verbosity(verbosity); + B0TRACKING::FastKalmanFilter->set_trackmap_out_name(B0TRACKING::TrackNodeName); + + //------------------------- + // B0 + //------------------------- + if (Enable::B0TRACKING && G4B0TRACKING::PROJECTION_B0) + { +// TRACKING::FastKalmanFilter->add_state_name("b0Truth"); +// TRACKING::ProjectionNames.insert("b0Truth"); + } + + se->registerSubsystem(B0TRACKING::FastKalmanFilter); + + // next, tracks with partial usage of the tracker stack + if (B0TRACKING::FastKalmanFilterB0Track == nullptr) + { + cout << __PRETTY_FUNCTION__ << " : missing the expected initialization for TRACKING::FastKalmanFilterB0Track." << endl; + exit(1); + } + InitFastKalmanFilter(B0TRACKING::FastKalmanFilterB0Track); + B0TRACKING::FastKalmanFilterB0Track->Verbosity(verbosity); + B0TRACKING::FastKalmanFilterB0Track->set_trackmap_out_name("B0TrackMap"); + B0TRACKING::FastKalmanFilterB0Track->enable_vertexing(false); + se->registerSubsystem(B0TRACKING::FastKalmanFilterB0Track); + return; +} + +//-----------------------------------------------------------------------------// + +void B0Tracking_Eval(const std::string &outputfile) +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::B0TRACKING_VERBOSITY); + //--------------- + // Fun4All server + //--------------- + + Fun4AllServer *se = Fun4AllServer::instance(); + + //---------------- + // Fast Tracking evaluation + //---------------- + + B0TrackFastSimEval *b0_fast_sim_eval = new B0TrackFastSimEval("B0FastTrackingEval"); + b0_fast_sim_eval->set_trackmapname(B0TRACKING::TrackNodeName); + b0_fast_sim_eval->set_filename(outputfile); + b0_fast_sim_eval->Verbosity(verbosity); + //------------------------- + // FEMC + //------------------------- + + cout << "Tracking_Eval(): configuration of track projections in B0TrackFastSimEval" << endl; + cout << "std::set B0TRACKING::ProjectionNames = {"; + bool first = true; + for (const string &proj : B0TRACKING::B0ProjectionNames) + { + if (first) + first = false; + else + cout << ", "; + cout << "\"" << proj << "\""; + + b0_fast_sim_eval->AddProjection(proj); + } + cout << "};" << endl; // override the TRACKING::ProjectionNames in eval macros + + se->registerSubsystem(b0_fast_sim_eval); + // now partial track fits + //B0TrackFastSimEval *b0_fast_sim_eval = new B0TrackFastSimEval("FastTrackingEval_B0TrackMap"); + b0_fast_sim_eval = new B0TrackFastSimEval("FastTrackingEval_B0TrackMap"); + b0_fast_sim_eval->set_trackmapname("B0TrackMap"); + b0_fast_sim_eval->set_filename(outputfile + ".B0TrackMap_debug.root"); + b0_fast_sim_eval->Verbosity(verbosity); + se->registerSubsystem(b0_fast_sim_eval); +} +#endif diff --git a/common/G4_BWD.C b/common/G4_BWD.C new file mode 100644 index 00000000..9b3fa0ff --- /dev/null +++ b/common/G4_BWD.C @@ -0,0 +1,150 @@ +#ifndef MACRO_G4BWD_C +#define MACRO_G4BWD_C + +#include + +//include our own Bwd Raw Tower Builder +#include + +#include + +#include +// Use Forward Cal Cell Reco . + +#include +// Include our Subsystem + +//Standard RawTowerDefs.h is modified + +#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 BWD = false; + bool BWDN[5]={true,false,false,false,false}; + bool BWD_ABSORBER = false; + bool BWD_CELL = false; + bool BWD_TOWER = false; + bool BWD_CLUSTER = false; + bool BWD_EVAL = false; + bool BWD_OVERLAPCHECK = false; + int BWD_VERBOSITY = 0; +} // namespace Enable + + +namespace G4BWD +{ + + string mapname[5]={"BWD_mapping_v1.txt","BWD_mapping_v2.txt","BWD_mapping_v3.txt","BWD_mapping_v4.txt","BWD_mapping_v5.txt"}; + double minz = -5000; + double maxz = -300; + double radius = 100; + + // Default set to B0 Ecal position at IP6 + + // Digitization (default photon digi): + RawTowerDigitizer::enu_digi_algorithm TowerDigi = RawTowerDigitizer::kNo_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 + +} // namespace G4B0ECAL + +void BWDInit() +{ +} + +void BWDSetup(PHG4Reco *g4Reco) +{ +//Done in G4_hFarBwdBeamLine.C +} + +void BWD_Cells(int verbosity = 0) +{ + return; +} + +void BWD_Towers() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::BWD_VERBOSITY); + + +for (int i = 0; i < 5; i++){ + if (!Enable::BWDN[i])continue; + Fun4AllServer *se = Fun4AllServer::instance(); + ostringstream mapping_bwd; + mapping_bwd << getenv("CALIBRATIONROOT") << "/BWD/mapping/"<Detector(Form("BWD_%d", i)); + tower_BWD->set_sim_tower_node_prefix("SIM"); + tower_BWD->GeometryTableFile(mapping_bwd.str()); + + se->registerSubsystem(tower_BWD); + + + RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer(Form("BWDRawTowerDigitizer_%d",i)); + TowerDigitizer->Detector(Form("BWD_%d", i)); + TowerDigitizer->Verbosity(verbosity); + TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kNo_digitization); + se->registerSubsystem(TowerDigitizer); + + RawTowerCalibration *TowerCalibration = new RawTowerCalibration(Form("BWDRawTowerCalibration_%d",i)); + TowerCalibration->Detector(Form("BWD_%d", i)); + TowerCalibration->Verbosity(verbosity); + TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration); + TowerCalibration->set_calib_const_GeV_ADC(1. ); + TowerCalibration->set_pedstal_ADC(0); + se->registerSubsystem(TowerCalibration); +} +} + +void BWD_Clusters() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::BWD_VERBOSITY); + +for (int i = 0; i < 5; i++){ + if (!Enable::BWDN[i])continue; + Fun4AllServer *se = Fun4AllServer::instance(); + RawClusterBuilderFwd *ClusterBuilder = new RawClusterBuilderFwd(Form("BWDRawClusterBuilderFwd_%d",i)); + ClusterBuilder->Detector(Form("BWD_%d", i)); + ClusterBuilder->Verbosity(verbosity); + ClusterBuilder->set_threshold_energy(0.100); + se->registerSubsystem(ClusterBuilder); +} + return; +} + +void BWD_Eval(const std::string &outputfile) +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::BWD_VERBOSITY); + string filename=outputfile.c_str(); +for (int i = 0; i < 5; i++){ + if (!Enable::BWDN[i])continue; + Fun4AllServer *se = Fun4AllServer::instance(); + CaloEvaluator *eval = new CaloEvaluator(Form("BWDEVALUATOR_%d",i), Form("BWD_%d", i), (filename+Form("_%d.root",i)) ); + eval->Verbosity(verbosity); + se->registerSubsystem(eval); +} + return; +} +#endif diff --git a/common/G4_Input.C b/common/G4_Input.C index 66f733da..13a8eb53 100644 --- a/common/G4_Input.C +++ b/common/G4_Input.C @@ -121,6 +121,152 @@ namespace Input exit(1); } + // --------------------------------------- + +// cout << Enable::HFARFWD_ION_ENERGY << " " << Enable::HFARBWD_E_ENERGY << endl; + + float ION_Energy = Enable::HFARFWD_ION_ENERGY; + float ELECTRON_Energy = Enable::HFARBWD_E_ENERGY; + + TString beam_setting_str; + beam_setting_str.Form("%.0fx%.0f", ION_Energy, ELECTRON_Energy); + + cout << "Beam scattering setting: " << beam_setting_str << endl; + + TString beam_opt; +// beam_opt = "ep-high-acceptance"; + beam_opt = Enable::BEAM_COLLISION_SETTING; +// beam_opt = "ep-high-divergence"; +// beam_opt = "eA"; +// cout << Enable::HFARFWD_ION_ENERGY << endl;; +// cout << Enable::IP6 << endl;; + + string beamFile; + + if (beam_opt == "ep-high-acceptance") { + beamFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_ep_high_acceptance_parameter.dat"; + } else if (beam_opt == "ep-high-divergence") { + beamFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_ep_high_divergence_parameter.dat"; + } else if (beam_opt == "eA") { + beamFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_eAu_parameter.dat"; + } else { + cout << "No beam scattering configuration file was identified." << endl; + gSystem->Exit(1); + } + +// cout << << endl; +// beamFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_ep_high_acceptance_parameter.dat"; +// beamFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_eAu_parameter.dat"; + + string settingname; + + double beta_star_p_h, beta_star_p_v, beta_star_e_h, beta_star_e_v; + double emit_p_h, emit_p_v, emit_e_h, emit_e_v; + + double beam_angular_divergence_p_h; + double beam_angular_divergence_p_v; + double beam_angular_divergence_e_h; + double beam_angular_divergence_e_v; + + double sigma_e_l, sigma_p_l; + + bool setting_found = false; + + float ION_Energy_Setting = 275; + float ION_Energy_Setting_diff = 275; + + std::ifstream infile(beamFile); + + if (infile.is_open()) + { + double biggest_z = 0.; + int imagnet = 0; + std::string line; + while (std::getline(infile, line)) + { + + if (line.find("#")!=std::string::npos) { + continue; + } + + std::istringstream iss(line); + + if (!(iss >> settingname >> beta_star_p_h >> beta_star_p_v >> beta_star_e_h >> beta_star_e_v >> emit_p_h >> emit_p_v >> emit_e_h >> emit_e_v >> beam_angular_divergence_p_h >> beam_angular_divergence_p_v >> beam_angular_divergence_e_h >> beam_angular_divergence_e_v >> sigma_p_l >> sigma_e_l)) + { + cout << "could not decode " << line << endl; + gSystem->Exit(1); + + } else { + + cout << line << endl; + + if (settingname==beam_setting_str) { + setting_found = true; + } + +// cout << "aaaaaa "<< settingname.find("x") << " " << settingname.substr(0, settingname.find("x")) << endl; + + float hadron_setting = stof(settingname.substr(0, settingname.find("x"))); + +// cout << hadron_setting - ION_Energy << endl; + + if (fabs(hadron_setting - ION_Energy) < ION_Energy_Setting_diff ) { + ION_Energy_Setting = hadron_setting; + ION_Energy_Setting_diff = fabs(hadron_setting - ION_Energy); + } + } + } + + // Reseting the pointer to the infile + infile.clear(); + infile.seekg(0,std::ios::beg); + // cout << infile.getline() << endl;; + + if(!setting_found) { + beam_setting_str.Form("%.0fx%.0f", ION_Energy_Setting, ELECTRON_Energy); + } + + while (std::getline(infile, line)) + { + + if (line.find("#")!=std::string::npos) { + continue; + } + + std::istringstream iss(line); + + if (!(iss >> settingname >> beta_star_p_h >> beta_star_p_v >> beta_star_e_h >> beta_star_e_v >> emit_p_h >> emit_p_v >> emit_e_h >> emit_e_v >> beam_angular_divergence_p_h >> beam_angular_divergence_p_v >> beam_angular_divergence_e_h >> beam_angular_divergence_e_v >> sigma_p_l >> sigma_e_l)) + { + cout << "could not decode " << line << endl; + gSystem->Exit(1); + + } else { + + cout << line << endl; + + if (settingname == beam_setting_str) { + +// cout << beta_star_p_h << " "<< beta_star_p_v << endl; +// cout << "BBBbBBBB" << endl; + + setting_found = true; + + break; + } + + } + } + + infile.close(); + } + + if (!setting_found) { + cout << "Could not find the specifed beam collision energy setting!" << endl; + gSystem->Exit(1); + } + + // --------------------------------------- + HepMCGen->PHHepMCGenHelper_Verbosity(VERBOSITY); //25mrad x-ing as in EIC CDR @@ -128,10 +274,16 @@ namespace Input // 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; +// 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; + +// beta_star_p_h = 80; +// beta_star_p_v = 7.1; +// beta_star_e_h = 59; +// 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; @@ -145,11 +297,17 @@ namespace Input ); // Table 4 of // https://github.com/eic/documents/blob/master/reports/general/Note-Simulations-BeamEffects.pdf +// HepMCGen->set_beam_angular_divergence_hv( +// 150e-6, 150e-6, // proton beam divergence horizontal & vertical +// 202e-6, 187e-6 // electron beam divergence horizontal & vertical +// ); + HepMCGen->set_beam_angular_divergence_hv( - 150e-6, 150e-6, // proton beam divergence horizontal & vertical - 202e-6, 187e-6 // electron beam divergence horizontal & vertical + beam_angular_divergence_p_h*1e-6, beam_angular_divergence_p_v*1e-6, // proton beam divergence horizontal & vertical + beam_angular_divergence_e_h*1e-6, beam_angular_divergence_e_v*1e-6 // electron beam divergence horizontal & vertical ); + // vertex shape from beam_bunch_sim HepMCGen->use_beam_bunch_sim(true); @@ -162,12 +320,27 @@ namespace Input // 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(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; + +// 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; +//// sigma_p_l = 6; +// 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; +//// sigma_e_l = 0.9; + + const double sigma_p_h = sqrt(beta_star_p_h * emit_p_h * 1e-7); + const double sigma_p_v = sqrt(beta_star_p_v * emit_p_v * 1e-7); + +// const double sigma_p_l = 6; +// sigma_p_l = 6; + + const double sigma_e_h = sqrt(beta_star_e_h * emit_e_h * 1e-7); + const double sigma_e_v = sqrt(beta_star_e_v * emit_e_v * 1e-7); + +// const double sigma_e_l = 0.9; +// sigma_e_l = 0.9; HepMCGen->set_beam_bunch_width( std::vector{sigma_p_h, sigma_p_v, sigma_p_l}, @@ -570,7 +743,7 @@ void InputRegister() } // here are the various utility modules which read particles and // put them onto the G4 particle stack - if (Input::HEPMC or Input::PYTHIA8 or Input::PYTHIA6 or Input::READEIC) + if (Input::HEPMC or Input::PYTHIA8 or Input::PYTHIA6 or Input::SARTRE or Input::READEIC) { if (Input::HEPMC) { diff --git a/common/G4_hFarBwdBeamLine_EIC.C b/common/G4_hFarBwdBeamLine_EIC.C index 6d9ab585..ec107314 100644 --- a/common/G4_hFarBwdBeamLine_EIC.C +++ b/common/G4_hFarBwdBeamLine_EIC.C @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -115,9 +116,9 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) string magFile; if (Enable::HFARFWD_MAGNETS_IP6) - magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farFwdBeamLineMagnets.dat"; + magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip6_h_farBwdBeamLineMagnets.dat"; else if (Enable::HFARFWD_MAGNETS_IP8) - magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip8_35mrad_h_farFwdBeamLineMagnets.dat"; + magFile = string(getenv("CALIBRATIONROOT")) + "/Beam/ip8_35mrad_h_farBwdBeamLineMagnets.dat"; else { cout << " You have to enable either the IP6 or IP8 Magnet configuration to define magnets! " << endl; @@ -125,7 +126,9 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) } // make magnet active volume if you want to study the hits - bool magnet_active = false; +// bool magnet_active = false; + bool magnet_active = true; + int absorberactive = 0; // if you insert numbers it only displays those magnets, do not comment out the set declaration @@ -141,9 +144,9 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) std::string line; while (std::getline(infile, line)) { - if (!line.compare(0, 1, "B") || - !line.compare(0, 1, "Q") || - !line.compare(0, 1, "S")) + if (!line.compare(0, 3, "eDB") || + !line.compare(0, 2, "eQ") || + !line.compare(0, 2, "eS")) { std::istringstream iss(line); string magname; @@ -194,15 +197,15 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) cout << "dipole_field_x: " << dipole_field_x << endl; cout << "fieldgradient: " << fieldgradient << endl; } - if (!magname.compare(0, 1, "B")) + if (!magname.compare(0, 3, "eDB")) { magtype = "DIPOLE"; } - else if (!magname.compare(0, 1, "Q")) + else if (!magname.compare(0, 2, "eQ")) { magtype = "QUADRUPOLE"; } - else if (!magname.compare(0, 1, "S")) + else if (!magname.compare(0, 2, "eS")) { magtype = "SEXTUPOLE"; } @@ -222,9 +225,10 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) //------------------------ // Linearly scaling down the magnetic field for lower energy proton - if( Enable::HFARFWD_ION_ENERGY != 275 ) { - float scaleFactor = Enable::HFARFWD_ION_ENERGY / 275. ; + if( Enable::HFARBWD_E_ENERGY != 18 ) { + float scaleFactor = Enable::HFARBWD_E_ENERGY / 18. ; dipole_field_x = dipole_field_x*scaleFactor; + fieldgradient = fieldgradient * scaleFactor; } if (magnetlist.empty() || magnetlist.find(imagnet) != magnetlist.end()) @@ -277,6 +281,44 @@ void hFarBwdDefineMagnets(PHG4Reco *g4Reco) void hFarBwdDefineDetectorsIP6(PHG4Reco *g4Reco) { + +// bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARBWD_OVERLAPCHECK; +// if (Enable::HFARBWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARBWD_VIRTUAL_DETECTORS_IP8) +// { +// cout << "You cannot have detectors enabled for both IP6 and IP8 ON at the same time" << endl; +// gSystem->Exit(1); +// } +// +// int verbosity = std::max(Enable::VERBOSITY, Enable::HFARBWD_VERBOSITY); +// +// auto *detBackward = new PHG4CylinderSubsystem("detBackward"); +// detBackward->SuperDetector("backTruth"); +// detBackward->set_double_param("place_x", 0); +// detBackward->set_double_param("place_y", 0); +//// detBackward->set_double_param("place_z", -500); +// detBackward->set_double_param("place_z", -500 - hFarBwdBeamLine::enclosure_center); +// detBackward->set_double_param("rot_y", 0); +// detBackward->set_double_param("radius", 0); +// detBackward->set_double_param("thickness", 30); // This is intentionally made large 25cm radius +// detBackward->set_double_param("length", 0.03); +// detBackward->set_string_param("material", "G4_Si"); +// +// detBackward->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); +// +// +// detBackward->SetActive(); +// detBackward->OverlapCheck(overlapCheck); +// detBackward->set_color(1, 0, 0, 0.5); +// +// detBackward->BlackHole(); +// if (verbosity) detBackward->Verbosity(verbosity); +// g4Reco->registerSubsystem(detBackward); +// +//// int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY); + + // ********************************************** + // Luminosity monitor + bool overlapCheck = Enable::OVERLAPCHECK || Enable::HFARBWD_OVERLAPCHECK; if (Enable::HFARBWD_VIRTUAL_DETECTORS_IP6 && Enable::HFARBWD_VIRTUAL_DETECTORS_IP8) { @@ -286,29 +328,198 @@ void hFarBwdDefineDetectorsIP6(PHG4Reco *g4Reco) int verbosity = std::max(Enable::VERBOSITY, Enable::HFARBWD_VERBOSITY); - auto *detBackward = new PHG4CylinderSubsystem("detBackward"); - detBackward->SuperDetector("backTruth"); - detBackward->set_double_param("place_x", 0); - detBackward->set_double_param("place_y", 0); -// detBackward->set_double_param("place_z", -500); - detBackward->set_double_param("place_z", -500 - hFarBwdBeamLine::enclosure_center); - detBackward->set_double_param("rot_y", 0); - detBackward->set_double_param("radius", 0); - detBackward->set_double_param("thickness", 30); // This is intentionally made large 25cm radius - detBackward->set_double_param("length", 0.03); - detBackward->set_string_param("material", "G4_Si"); - - detBackward->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); - - - detBackward->SetActive(); - detBackward->OverlapCheck(overlapCheck); - detBackward->set_color(1, 0, 0, 0.5); - - detBackward->BlackHole(); - if (verbosity) detBackward->Verbosity(verbosity); - g4Reco->registerSubsystem(detBackward); - +// auto *detLumi = new PHG4CylinderSubsystem("detLumi"); +// detLumi->SuperDetector("backLumi"); +// detLumi->set_double_param("place_x", 0); +// detLumi->set_double_param("place_y", 0); +// detLumi->set_double_param("place_z", -3400 - hFarBwdBeamLine::enclosure_center); +// detLumi->set_double_param("rot_y", 0); +// detLumi->set_double_param("radius", 0); +// detLumi->set_double_param("thickness", 20); // This is intentionally made large 25cm radius +// detLumi->set_double_param("length", 0.03); +// detLumi->set_string_param("material", "G4_Si"); +// detLumi->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + +// detLumi->SetActive(); +// detLumi->OverlapCheck(overlapCheck); +// detLumi->set_color(0, 1, 0, 0.5); + +// detLumi->BlackHole(); +// if (verbosity) detLumi->Verbosity(verbosity); +// g4Reco->registerSubsystem(detLumi); + + // ********************************************** + // Low Q2 Tagger + // There are two set of Q2 taggers: 1st at 24m; 2nd at 37m. + // + // 1st Set of Low Q2 tagger, z location 24m + +// auto *detLowQ2Tag_1 = new PHG4CylinderSubsystem("detLowQ2Tag_1"); + +// auto *detLowQ2Tag_1= new PHG4BlockSubsystem("detLowQ2Tag_1"); + +// detLowQ2Tag_1->SuperDetector("backLowQ2Tag_1"); +// detLowQ2Tag_1->set_double_param("place_x", -50); +// detLowQ2Tag_1->set_double_param("place_y", 0); +// detLowQ2Tag_1->set_double_param("place_z", -2400 - hFarBwdBeamLine::enclosure_center); +// detLowQ2Tag_1->set_double_param("rot_y", 0); +// detLowQ2Tag_1->set_double_param("radius", 0); +// detLowQ2Tag_1->set_double_param("thickness", 30); // This is intentionally made large 25cm radius +// detLowQ2Tag_1->set_double_param("size_x", 50); // This is intentionally made large 25cm radius +// detLowQ2Tag_1->set_double_param("size_y", 35); // This is intentionally made large 25cm radius +// detLowQ2Tag_1->set_double_param("size_z", 0.03); // This is intentionally made large 25cm radius +// detLowQ2Tag_1->set_double_param("length", 0.03); +// detLowQ2Tag_1->set_string_param("material", "G4_Si"); +// detLowQ2Tag_1->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + +// detLowQ2Tag_1->SetActive(); +// detLowQ2Tag_1->OverlapCheck(overlapCheck); +// detLowQ2Tag_1->set_color(1, 0, 0, 0.5); + +// detLowQ2Tag_1->BlackHole(); +// if (verbosity) detLowQ2Tag_1->Verbosity(verbosity); +// g4Reco->registerSubsystem(detLowQ2Tag_1); + + + // 2nd Set of Low Q2 tagger, z location 37m + +// auto *detLowQ2Tag_2 = new PHG4CylinderSubsystem("detLowQ2Tag_2"); +// auto *detLowQ2Tag_2 = new PHG4BlockSubsystem("detLowQ2Tag_2"); +// detLowQ2Tag_2->SuperDetector("backLowQ2Tag_2"); +// detLowQ2Tag_2->set_double_param("place_x", -80); +// detLowQ2Tag_2->set_double_param("place_y", 0); +// detLowQ2Tag_2->set_double_param("place_z", -3700 - hFarBwdBeamLine::enclosure_center); +// detLowQ2Tag_2->set_double_param("rot_y", 0); +// detLowQ2Tag_2->set_double_param("radius", 0); +// detLowQ2Tag_2->set_double_param("thickness", 30); // This is intentionally made large 25cm radius +// detLowQ2Tag_2->set_double_param("size_x", 50); // This is intentionally made large 25cm radius +// detLowQ2Tag_2->set_double_param("size_y", 35); // This is intentionally made large 25cm radius +// detLowQ2Tag_2->set_double_param("size_z", 0.03); // This is intentionally made large 25cm radius +// detLowQ2Tag_2->set_double_param("length", 0.03); +// detLowQ2Tag_2->set_string_param("material", "G4_Si"); +// detLowQ2Tag_2->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + +// detLowQ2Tag_2->SetActive(); +// detLowQ2Tag_2->OverlapCheck(overlapCheck); +// detLowQ2Tag_2->set_color(1, 0, 0, 0.5); + +// detLowQ2Tag_2->BlackHole(); +// if (verbosity) detLowQ2Tag_2->Verbosity(verbosity); +// g4Reco->registerSubsystem(detLowQ2Tag_2); + + if (Enable::BWD) { + int DetNr = 5; //number of Backward Detectors + string bwddetname[5]={" is the first Q2 tagger", " is the second Q2 tagger", " is Lumi 0"," is Lumi +"," is Lumi -"}; +// string mapname[5]={"BWD_mapping_v1.txt","BWD_mapping_v2.txt","BWD_mapping_v3.txt","BWD_mapping_v3.txt","BWD_mapping_v3.txt"}; + float placex[5] = {-50,-80,0,0,0}; + //float placex[5] = {0,0,0,0,0}; + float placey[5] = {0,0,0,30,-30}; + float placez[5] = {-118.5, - 818.5, -818.5,-718.5,-718.5}; + float length = 20; + float Si_length = .1; + float Cu_length = .2; + float C_length = 40; + float width[5] = {40.5, 30, 16,16,16}; + float height[5] = {40.5, 21, 16,8,8}; + int deti=0; + for (int i = 0; i < DetNr; i ++){ + if (!Enable::BWDN[i])continue; + cout <<"Detector "<SetTowerMappingFile(mapping_bwd.str()); + Bwd->SuperDetector(Form("BWD_%d", i)); + Bwd->set_double_param("place_x", placex[i]); + Bwd->set_double_param("place_y", placey[i]); + Bwd->set_double_param("place_z", placez[i]-length/2); + Bwd->set_double_param("length", length); + Bwd->set_double_param("width", width[i]); + Bwd->set_double_param("height", height[i]); + Bwd->set_string_param("material", "G4_PbWO4"); + Bwd->set_double_param("detid",deti); + Bwd->set_double_param("global_x", placex[i]); + Bwd->set_double_param("global_y", placey[i]); + Bwd->set_double_param("global_z", placez[i]-length/2 + hFarBwdBeamLine::enclosure_center); + Bwd->set_int_param("lightyield",1); //Note additional parameter for storing Light Yield in B0 Ecal + Bwd->SetActive(true); + if (verbosity) + Bwd->Verbosity(verbosity); + Bwd->OverlapCheck(overlapCheck); + Bwd->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + g4Reco->registerSubsystem(Bwd); + deti++; + if (i != 2){ + auto *Bwdt = new EICG4BwdSubsystem("BWD"); + Bwdt->SuperDetector(Form("BWDt_%d", i)); + Bwdt->set_double_param("place_x", placex[i]); + Bwdt->set_double_param("place_y", placey[i]); + Bwdt->set_double_param("place_z", placez[i]+Si_length/2+Cu_length/2+1); + Bwdt->set_double_param("length", Si_length); + Bwdt->set_double_param("width", width[i]); + Bwdt->set_double_param("height", height[i]); + Bwdt->set_string_param("material", "G4_Si"); + Bwdt->set_double_param("detid",deti); + Bwdt->set_double_param("global_x", placex[i]); + Bwdt->set_double_param("global_y", placey[i]); + Bwdt->set_double_param("global_z", placez[i]+Si_length/2 +Cu_length/2+1+ hFarBwdBeamLine::enclosure_center); + Bwdt->set_int_param("lightyield",0); //Note additional parameter for storing Light Yield in B0 Ecal + Bwdt->SetActive(true); + if (verbosity) + Bwdt->Verbosity(verbosity); + Bwdt->OverlapCheck(overlapCheck); + Bwdt->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + g4Reco->registerSubsystem(Bwdt); + deti++; + auto *Bwdd = new EICG4BwdSubsystem("BWD"); + Bwdd->SuperDetector(Form("BWDd_%d", i)); + Bwdd->set_double_param("place_x", placex[i]); + Bwdd->set_double_param("place_y", placey[i]); + Bwdd->set_double_param("place_z", placez[i]+Cu_length/2+1); + Bwdd->set_double_param("length", Cu_length); + Bwdd->set_double_param("width", width[i]); + Bwdd->set_double_param("height", height[i]); + Bwdd->set_string_param("material", "G4_Cu"); + Bwdd->set_double_param("detid",deti); + Bwdd->set_double_param("global_x", placex[i]); + Bwdd->set_double_param("global_y", placey[i]); + Bwdd->set_double_param("global_z", placez[i]+Cu_length/2 +1+ hFarBwdBeamLine::enclosure_center); + Bwdd->set_int_param("lightyield",0); //Note additional parameter for storing Light Yield in B0 Ecal + Bwdd->SetActive(false); + if (verbosity) + Bwdd->Verbosity(verbosity); + Bwdd->OverlapCheck(overlapCheck); + Bwdd->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + g4Reco->registerSubsystem(Bwdd); + deti++; + } + else{ + auto *Bwdd = new EICG4BwdSubsystem("BWD"); + Bwdd->SuperDetector(Form("BWDd_%d", i)); + Bwdd->set_double_param("place_x", placex[i]); + Bwdd->set_double_param("place_y", placey[i]); + Bwdd->set_double_param("place_z", placez[i]+C_length/2+1); + Bwdd->set_double_param("length", C_length); + Bwdd->set_double_param("width", width[i]); + Bwdd->set_double_param("height", height[i]); + Bwdd->set_string_param("material", "G4_C"); + Bwdd->set_double_param("detid",deti); + Bwdd->set_double_param("global_x", placex[i]); + Bwdd->set_double_param("global_y", placey[i]); + Bwdd->set_double_param("global_z", placez[i]+C_length/2 +1+ hFarBwdBeamLine::enclosure_center); + Bwdd->set_int_param("lightyield",0); //Note additional parameter for storing Light Yield in B0 Ecal + Bwdd->SetActive(false); + if (verbosity) + Bwdd->Verbosity(verbosity); + Bwdd->OverlapCheck(overlapCheck); + Bwdd->SetMotherSubsystem(hFarBwdBeamLine::hFarBwdBeamLineEnclosure); + g4Reco->registerSubsystem(Bwdd); + deti++; + + } + } + } // int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY); } diff --git a/common/G4_hFarFwdBeamLine_EIC.C b/common/G4_hFarFwdBeamLine_EIC.C index 46e3e384..fc168e71 100644 --- a/common/G4_hFarFwdBeamLine_EIC.C +++ b/common/G4_hFarFwdBeamLine_EIC.C @@ -30,6 +30,7 @@ float PosFlip(float pos); float AngleFlip(float angle); float MagFieldFlip(float Bfield); + // This creates the Enable Flag to be used in the main steering macro namespace Enable { @@ -44,10 +45,9 @@ namespace Enable bool ZDC_DISABLE_BLACKHOLE = false; bool B0_DISABLE_HITPLANE = false; bool B0_FULLHITPLANE = false; + bool B0_VAR_PIPE_HOLE = false; + bool B0_CIRCLE_PIPE_HOLE = false; bool RP_DISABLE_HITPLANE = false; - bool RP_FULLHITPLANE = false; - bool RP2nd_DISABLE_HITPLANE = false; - bool RP2nd_FULLHITPLANE = false; bool B0ECALTOWERS = true; //Set to 'false' for nice PackMan views. Set 'true' for physics studies. //enabled automatically in hFarFwdBeamLineInit(), unless overridden by user @@ -58,8 +58,6 @@ namespace Enable bool HFARFWD_VIRTUAL_DETECTORS_IP6 = false; bool HFARFWD_VIRTUAL_DETECTORS_IP8 = false; - float HFARFWD_ION_ENERGY = 0; - bool FFR_EVAL = false; } // namespace Enable @@ -364,187 +362,92 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detOM->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); g4Reco->registerSubsystem(detOM); } - ////********************* - // RP - // Three choices: 1. Realistic detector; 2. Circulat or square plane; 3. hit plane with realistic detector goemetry - - - const int rpDetNr = 2; - const double rp_zCent[rpDetNr] = {2600, 2800}; - const double rp_xCent[rpDetNr] = {-83.22, -92.20}; - - if (Enable::RP_DISABLE_HITPLANE) + + + + //---------------------- + // Roman Pots + //---------------------- + + if( ! Enable::RP_DISABLE_HITPLANE ) { - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 30.0; //detector width - const double rp_y = 10.0; //detector height - const double rot_y = 0.047; //rotation angle - for (int i = 0; i < rpDetNr; i++) - { - auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i), 2 * 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(rot_y * TMath::RadToDeg())); - detRP->set_double_param("rp_x", rp_x); - detRP->set_double_param("rp_y", rp_y); - detRP->set_double_param("hole_x", hole_x); - detRP->set_double_param("hole_y", hole_y); - detRP->set_double_param("length", rpSi_zLen); - detRP->set_string_param("material", "G4_Si"); - detRP->set_double_param("detid", 2 * i); - detRP->set_double_param("pipe_x", rppipe_x); - detRP->set_double_param("pipe_y", 0); - detRP->set_double_param("pipe_z", 0); - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP->SetActive(true); - if (verbosity) detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - - auto *detRPe = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i + 1), 2 * i + 1); - detRPe->SuperDetector("rpTruth"); - detRPe->set_double_param("place_x", PosFlip(rp_xCent[i])); - detRPe->set_double_param("place_y", 0); - detRPe->set_double_param("place_z", rp_zCent[i] - hFarFwdBeamLine::enclosure_center + (rpSi_zLen + rpCu_zLen) / 2); - detRPe->set_double_param("rot_y", AngleFlip(rot_y * TMath::RadToDeg())); - detRPe->set_double_param("rp_x", rp_x); - detRPe->set_double_param("rp_y", rp_y); - detRPe->set_double_param("hole_x", hole_x); - detRPe->set_double_param("hole_y", hole_y); - detRPe->set_double_param("length", rpCu_zLen); - detRPe->set_string_param("material", "G4_Cu"); - detRPe->set_double_param("detid", 2 * i + 1); - detRPe->set_double_param("pipe_x", rppipe_x); - detRPe->set_double_param("pipe_y", 0); - detRPe->set_double_param("pipe_z", 0); - detRPe->OverlapCheck(overlapCheck); - detRPe->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRPe->SetActive(true); - if (verbosity) detRPe->Verbosity(verbosity); - g4Reco->registerSubsystem(detRPe); - } - } - else - { - if (Enable::RP_FULLHITPLANE) - { - for (int i = 0; i < rpDetNr; i++) - { - ////********************* - //// Square design - //// 25 cm in x - // - // auto *detRP = new PHG4BlockSubsystem(Form("rpTruth_%d",i)); - //// detRP->SuperDetector("RomanPots"); - // detRP->SuperDetector(Form("RomanPots_%d",i)); - // detRP->set_double_param("place_x",rp_xCent[i]); - // detRP->set_double_param("place_y",0); - // detRP->set_double_param("place_z",rp_zCent[i]); - // detRP->set_double_param("rot_y",-0.025*TMath::RadToDeg()); - // detRP->set_double_param("size_x",25); - // detRP->set_double_param("size_y",10); - // detRP->set_double_param("size_z",0.03); - // detRP->set_string_param("material","G4_Si"); - - ////********************* - //// Disk design - //// 50 cm in x - - 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.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); - detRP->set_string_param("material", "G4_Si"); - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - - detRP->SetActive(); - if (verbosity) detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - } - } - else - { - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 30.0; //detector width - const double rp_y = 10.0; //detector height - const double rot_y = 0.047; //rotation angle - for (int i = 0; i < rpDetNr; i++) - { - auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i), 2 * 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(rot_y * TMath::RadToDeg())); - detRP->set_double_param("rp_x", rp_x); - detRP->set_double_param("rp_y", rp_y); - detRP->set_double_param("hole_x", hole_x); - detRP->set_double_param("hole_y", hole_y); - detRP->set_double_param("length", rpSi_zLen); - detRP->set_string_param("material", "G4_Si"); - detRP->set_double_param("detid", 2 * i); - detRP->set_double_param("pipe_x", rppipe_x); - detRP->set_double_param("pipe_y", 0); - detRP->set_double_param("pipe_z", 0); - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP->SetActive(true); - if (verbosity) detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - } - } + string paramFile = string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP6.dat"; + int Nlayers = GetParameterFromFile (paramFile, "Number_layers"); + + for( int layer = 0; layer < Nlayers; layer++ ) { + auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", layer), layer); + detRP->SuperDetector("rpTruth"); + detRP->SetParameterFile( paramFile ); + detRP->set_double_param("FFenclosure_center", hFarFwdBeamLine::enclosure_center ); + detRP->set_int_param("layerNumber", layer + 1); + + detRP->OverlapCheck(overlapCheck); + detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); + detRP->SetActive(true); + if (verbosity) detRP->Verbosity(verbosity); + g4Reco->registerSubsystem(detRP); + } } + //--------------------------------- // B0 implementation // Three choices: 1. Realistic detector; 2. Circulat plane; 3. hit plane with realistic detector goemetry double b0tr_z = 0; //Subsystem position relative to B0 magnet (for iterator) - - if (Enable::B0_DISABLE_HITPLANE) { - - // Choice 1 realistic detector - cout << "Realistic B0"<SuperDetector(Form("b0Truth_%d", i)); detB0->set_double_param("place_x", 0); detB0->set_double_param("place_y", 0); // detB0->set_int_param("ispipe", 0); //for future pipe implementation detB0->set_double_param("pipe_hole", pipe_hole); + detB0->set_double_param("cable_hole", cable_hole); detB0->set_double_param("outer_radius", b0_radius); detB0->set_double_param("d_radius", d_radius); detB0->set_double_param("length", b0Si_zLen); @@ -555,6 +458,10 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detB0->set_double_param("pipe_x", pipe_x); detB0->set_double_param("pipe_y", 0); detB0->set_double_param("pipe_z", 0); + detB0->set_double_param("pipe_hole_r", pipe_hole_r); + detB0->set_double_param("cable_x", cable_x); + detB0->set_double_param("cable_y", 0); + detB0->set_double_param("cable_z", 0); detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet detB0->SetActive(true); if (verbosity) @@ -563,7 +470,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); // For B0 Tracking Implementation -/* if (Enable::B0TRACKING){ + if (Enable::B0TRACKING){ if (B0TRACKING::FastKalmanFilter) { B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, @@ -585,7 +492,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); } } -*/ + auto *detB0e = new EICG4B0Subsystem(Form("b0Dead_%d", i), i); detB0e->SuperDetector("b0Dead"); // detB0e->set_int_param("ispipe", 0); //for future pipe implementation @@ -596,13 +503,17 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detB0e->set_double_param("pipe_x", pipe_x); detB0e->set_double_param("pipe_y", 0); detB0e->set_double_param("pipe_z", 0); + detB0e->set_double_param("pipe_hole_r", pipe_hole_r); + detB0e->set_double_param("cable_x", cable_x); + detB0e->set_double_param("cable_y", 0); + detB0e->set_double_param("cable_z", 0); detB0e->set_double_param("outer_radius", b0_radius); detB0e->set_double_param("length", b0Cu_zLen); detB0e->set_string_param("material", "G4_Cu"); detB0e->set_double_param("detid",i); detB0e->set_double_param("startAngle",start_angle); detB0e->set_double_param("spanningAngle",spanning_angle); - detB0e->set_double_param("place_z", (b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2) +(b0Cu_zLen+b0Si_zLen)/2) ); // relative to B0 magnet + detB0e->set_double_param("place_z", b0tr_z +(b0Cu_zLen+b0Si_zLen)/2) ; // relative to B0 magnet detB0e->SetActive(false); if (verbosity) detB0e->Verbosity(verbosity); @@ -612,11 +523,21 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) } if (Enable::B0ECAL) { + pipe_hole = b0Mag_zLen*cross_angle; + pipe_x = - cross_angle*b0Mag_zCent - hFarFwdBeamLine::B0Magnet_x; + if (Enable::B0_CIRCLE_PIPE_HOLE){ + pipe_hole = 0.1; + pipe_hole_r = pipe_hole_r + b0Mag_zLen*cross_angle/2; + } + cout <<"Starting B0 ECAL "<SetTowerMappingFile(mapping_b0ecal.str()); @@ -628,6 +549,10 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0Ecal->set_double_param("pipe_x", pipe_x); B0Ecal->set_double_param("pipe_y", 0); B0Ecal->set_double_param("pipe_z", 0); + B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecal->set_double_param("cable_x", cable_x); + B0Ecal->set_double_param("cable_y", 0); + B0Ecal->set_double_param("cable_z", 0); B0Ecal->set_double_param("length", b0Ecal_zLen); B0Ecal->set_double_param("outer_radius", b0_radius); B0Ecal->set_double_param("d_radius", d_radius); @@ -656,6 +581,10 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0Ecal->set_double_param("pipe_x", pipe_x); B0Ecal->set_double_param("pipe_y", 0); B0Ecal->set_double_param("pipe_z", 0); + B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecal->set_double_param("cable_x", cable_x); + B0Ecal->set_double_param("cable_y", 0); + B0Ecal->set_double_param("cable_z", 0); B0Ecal->set_double_param("length", b0Ecal_zLen); B0Ecal->set_double_param("outer_radius", b0_radius); B0Ecal->set_double_param("d_radius", d_radius); @@ -681,6 +610,10 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0Ecale->set_double_param("pipe_x", pipe_x); B0Ecale->set_double_param("pipe_y", 0); B0Ecale->set_double_param("pipe_z", 0); + B0Ecale->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecale->set_double_param("cable_x", cable_x); + B0Ecale->set_double_param("cable_y", 0); + B0Ecale->set_double_param("cable_z", 0); B0Ecale->set_double_param("length", b0Cu_zLen); B0Ecale->set_double_param("d_radius", d_radius); B0Ecale->set_double_param("outer_radius", b0_radius); @@ -703,21 +636,12 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) // Choice 2 circular hit planes cout << "Circular hit planes"<SuperDetector(Form("b0Truth_%d", i)); + detB0->SuperDetector("b0Truth"); + //detB0->SuperDetector(Form("b0Truth_%d", i)); detB0->set_double_param("radius", 0); detB0->set_double_param("thickness", 20); detB0->set_double_param("length", 0.1); @@ -730,7 +654,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); -/* if (Enable::B0TRACKING){ + if (Enable::B0TRACKING){ if (B0TRACKING::FastKalmanFilter) { B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, @@ -752,7 +676,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); } } -*/ + } } else { @@ -760,47 +684,53 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) /// Fun4All default B0 planes /// Choice 3 Hit planes with real detector geometry cout << "Realistic hit planes"<SuperDetector(Form("b0Truth_%d", i)); - detB0->set_double_param("place_x", 0); - detB0->set_double_param("place_y", 0); - // detB0->set_int_param("ispipe", 0); //for future pipe implementation - detB0->set_double_param("pipe_hole", pipe_hole); - detB0->set_double_param("outer_radius", b0_radius); - detB0->set_double_param("d_radius", d_radius); - detB0->set_double_param("length", b0Si_zLen); - detB0->set_string_param("material", "G4_Si"); - detB0->set_double_param("detid",i); - detB0->set_double_param("startAngle",start_angle); - detB0->set_double_param("spanningAngle",spanning_angle); - detB0->set_double_param("pipe_x", pipe_x); - detB0->set_double_param("pipe_y", 0); - detB0->set_double_param("pipe_z", 0); - detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet + if (Enable::B0_VAR_PIPE_HOLE){ + pipe_hole = b0tr[i]*cross_angle; + pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else if (Enable::B0_CIRCLE_PIPE_HOLE){ + pipe_hole = 0.1; + pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2; + pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else { + pipe_hole = b0tr[b0DetNr-1]*cross_angle; + pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + cout <<"Starting B0 Tracker layer "<SuperDetector(Form("b0Truth_%d", i)); + detB0->set_double_param("place_x", 0); + detB0->set_double_param("place_y", 0); + // detB0->set_int_param("ispipe", 0); //for future pipe implementation + detB0->set_double_param("pipe_hole", pipe_hole); + detB0->set_double_param("cable_hole", cable_hole); + detB0->set_double_param("outer_radius", b0_radius); + detB0->set_double_param("d_radius", d_radius); + detB0->set_double_param("length", b0Si_zLen); + detB0->set_string_param("material", "G4_Si"); + detB0->set_double_param("startAngle",start_angle); + detB0->set_double_param("spanningAngle",spanning_angle); + detB0->set_double_param("detid",i); + detB0->set_double_param("pipe_x", pipe_x); + detB0->set_double_param("pipe_y", 0); + detB0->set_double_param("pipe_z", 0); + detB0->set_double_param("pipe_hole_r", pipe_hole_r); + detB0->set_double_param("cable_x", cable_x); + detB0->set_double_param("cable_y", 0); + detB0->set_double_param("cable_z", 0); + detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet detB0->SetActive(true); if (verbosity) detB0->Verbosity(verbosity); detB0->OverlapCheck(overlapCheck); detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); -/* if (Enable::B0TRACKING){ + if (Enable::B0TRACKING){ if (B0TRACKING::FastKalmanFilter) { B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, @@ -822,7 +752,7 @@ void hFarFwdDefineDetectorsIP6(PHG4Reco *g4Reco) B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); } } -*/ } + } } } @@ -902,285 +832,29 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) } - //------------------ - // Roman pot set #1 - ////********************* - // RP - // Three choices: 1. Realistic detector; 2. Circulat or square plane; 3. hit plane with realistic detector goemetry + //---------------------- + // Roman Pots: Both sets before and near the secondary focus + //---------------------- - const int rpDetNr = 2; - const double rp_zCent[rpDetNr] = {2600, 2800}; - const double rp_xCent[rpDetNr] = {75.6, 78.15}; - - if (Enable::RP_DISABLE_HITPLANE) - { - // Circular disk design (16cm in) - auto *detRP = new PHG4CylinderSubsystem(Form("rpTruth_%d", 0), 0); - detRP->SuperDetector("rpTruth"); - detRP->set_double_param("place_x", PosFlip(rp_xCent[0])); - detRP->set_double_param("place_y", 0); - detRP->set_double_param("place_z", rp_zCent[0] - hFarFwdBeamLine::enclosure_center); - detRP->set_double_param("rot_y", AngleFlip(-0.035 * TMath::RadToDeg())); - detRP->set_double_param("radius", 5); - detRP->set_double_param("thickness", 10); // 16 cm circulr to cover 25cm x20cm square (IR design) - detRP->set_double_param("length", 0.03); - detRP->set_string_param("material", "G4_Si"); - - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 25.0; //detector width - const double rp_y = 20.0; //detector height - const double rot_y = -0.035; //rotation angle - for (int i = 0; i < rpDetNr; i++) - { - auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i + 1), 2 * i + 1); - 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(rot_y * TMath::RadToDeg())); - detRP->set_double_param("rp_x", rp_x); - detRP->set_double_param("rp_y", rp_y); - detRP->set_double_param("hole_x", hole_x); - detRP->set_double_param("hole_y", hole_y); - detRP->set_double_param("length", rpSi_zLen); - detRP->set_string_param("material", "G4_Si"); - detRP->set_double_param("detid", 2 * i); - detRP->set_double_param("pipe_x", rppipe_x); - detRP->set_double_param("pipe_y", 0); - detRP->set_double_param("pipe_z", 0); - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP->SetActive(true); - if (verbosity) detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - - auto *detRPe = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i + 2), 2 * i + 2); - detRPe->SuperDetector("rpTruth"); - detRPe->set_double_param("place_x", PosFlip(rp_xCent[i])); - detRPe->set_double_param("place_y", 0); - detRPe->set_double_param("place_z", rp_zCent[i] - hFarFwdBeamLine::enclosure_center + (rpSi_zLen + rpCu_zLen) / 2); - detRPe->set_double_param("rot_y", AngleFlip(rot_y * TMath::RadToDeg())); - detRPe->set_double_param("rp_x", rp_x); - detRPe->set_double_param("rp_y", rp_y); - detRPe->set_double_param("hole_x", hole_x); - detRPe->set_double_param("hole_y", hole_y); - detRPe->set_double_param("length", rpCu_zLen); - detRPe->set_string_param("material", "G4_Cu"); - detRPe->set_double_param("detid", 2 * i + 1); - detRPe->set_double_param("pipe_x", rppipe_x); - detRPe->set_double_param("pipe_y", 0); - detRPe->set_double_param("pipe_z", 0); - detRPe->OverlapCheck(overlapCheck); - detRPe->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRPe->SetActive(true); - if (verbosity) detRPe->Verbosity(verbosity); - g4Reco->registerSubsystem(detRPe); - } - } - else + if( ! Enable::RP_DISABLE_HITPLANE ) { - if (Enable::RP_FULLHITPLANE) - { - for (int i = 0; i < rpDetNr; i++) - { - // Circular disk design (16cm in) - 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.035 * TMath::RadToDeg())); - detRP->set_double_param("radius", 5); - detRP->set_double_param("thickness", 10); // 16 cm circulr to cover 25cm x20cm square (IR design) - detRP->set_double_param("length", 0.03); - detRP->set_string_param("material", "G4_Si"); - - // //------------------------------------ - // /// Square Design - // auto *detRP = new PHG4BlockSubsystem(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.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_z", 0.03); - // detRP->set_string_param("material", "G4_Si"); - - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP->SetActive(); - if (verbosity) - detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - } - } - else - { - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 25.0; //detector width - const double rp_y = 20.0; //detector height - const double rot_y = -0.029; //rotation angle - for (int i = 0; i < rpDetNr; i++) - { - auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", 2 * i), 2 * 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(rot_y * TMath::RadToDeg())); - detRP->set_double_param("rp_x", rp_x); - detRP->set_double_param("rp_y", rp_y); - detRP->set_double_param("hole_x", hole_x); - detRP->set_double_param("hole_y", hole_y); - detRP->set_double_param("length", rpSi_zLen); - detRP->set_string_param("material", "G4_Si"); - detRP->set_double_param("detid", 2 * i); - detRP->set_double_param("pipe_x", rppipe_x); - detRP->set_double_param("pipe_y", 0); - detRP->set_double_param("pipe_z", 0); - detRP->OverlapCheck(overlapCheck); - detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP->SetActive(true); - if (verbosity) detRP->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP); - } - } - } - - //------------------ - // Roman pot set #2 before and after the secondary focus - // Three choices: 1. Realistic detector; 2. Circulat or square plane; 3. hit plane with realistic detector goemetry - const int rp2ndDetNr = 2; - const double rp_2nd_xCent[rp2ndDetNr] = {101.94, 106.94}; - const double rp_2nd_zCent[rp2ndDetNr] = {4300, 4450}; - - if (Enable::RP_DISABLE_HITPLANE) - { - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 25.0; //detector width - const double rp_y = 20.0; //detector height - const double rot_y = -0.029; //rotation angle - for (int i = 0; i < rp2ndDetNr; i++) - { - auto *detRP_2nd = new EICG4RPSubsystem(Form("rpTruth2_%d", 2 * i), 2 * i); - detRP_2nd->SuperDetector("rpTruth2"); - 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] - hFarFwdBeamLine::enclosure_center); - detRP_2nd->set_double_param("rot_y", AngleFlip(rot_y * TMath::RadToDeg())); - detRP_2nd->set_double_param("rp_x", rp_x); - detRP_2nd->set_double_param("rp_y", rp_y); - detRP_2nd->set_double_param("hole_x", hole_x); - detRP_2nd->set_double_param("hole_y", hole_y); - detRP_2nd->set_double_param("length", rpSi_zLen); - detRP_2nd->set_string_param("material", "G4_Si"); - detRP_2nd->set_double_param("detid", 2 * i); - detRP_2nd->set_double_param("pipe_x", rppipe_x); - detRP_2nd->set_double_param("pipe_y", 0); - detRP_2nd->set_double_param("pipe_z", 0); - detRP_2nd->OverlapCheck(overlapCheck); - detRP_2nd->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP_2nd->SetActive(true); - if (verbosity) detRP_2nd->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP_2nd); - - auto *detRPe_2nd = new EICG4RPSubsystem(Form("rpTruth2_%d", 2 * i + 1), 2 * i + 1); - detRPe_2nd->SuperDetector("rpTruth2"); - detRPe_2nd->set_double_param("place_x", PosFlip(rp_2nd_xCent[i])); - detRPe_2nd->set_double_param("place_y", 0); - detRPe_2nd->set_double_param("place_z", rp_2nd_zCent[i] - hFarFwdBeamLine::enclosure_center + (rpSi_zLen + rpCu_zLen) / 2); - detRPe_2nd->set_double_param("rot_y", AngleFlip(rot_y * TMath::RadToDeg())); - detRPe_2nd->set_double_param("rp_x", rp_x); - detRPe_2nd->set_double_param("rp_y", rp_y); - detRPe_2nd->set_double_param("hole_x", hole_x); - detRPe_2nd->set_double_param("hole_y", hole_y); - detRPe_2nd->set_double_param("length", rpCu_zLen); - detRPe_2nd->set_string_param("material", "G4_Cu"); - detRPe_2nd->set_double_param("detid", 2 * i + 1); - detRPe_2nd->set_double_param("pipe_x", rppipe_x); - detRPe_2nd->set_double_param("pipe_y", 0); - detRPe_2nd->set_double_param("pipe_z", 0); - detRPe_2nd->OverlapCheck(overlapCheck); - detRPe_2nd->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRPe_2nd->SetActive(true); - if (verbosity) detRPe_2nd->Verbosity(verbosity); - g4Reco->registerSubsystem(detRPe_2nd); - } - } - else - { - if (Enable::RP_FULLHITPLANE) - { - for (int i = 0; i < rp2ndDetNr; i++) - { - auto *detRP_2nd = new PHG4BlockSubsystem(Form("rpTruth2_%d", i), i); - detRP_2nd->SuperDetector("rpTruth2"); - 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] - hFarFwdBeamLine::enclosure_center); - detRP_2nd->set_double_param("rot_y", AngleFlip(-0.029 * TMath::RadToDeg())); - 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); - g4Reco->registerSubsystem(detRP_2nd); - } - } - else - { - const double rpCu_zLen = .2; //B0 dead material length - const double rpSi_zLen = .03; //B0 Si length - const double hole_x = 10.0; //detector cut off for beam pipe - const double rppipe_x = 0.0; //detector cut off for beam pipe position - const double hole_y = 3.0; //detector cut off for beam pipe - const double rp_x = 25.0; //detector width - const double rp_y = 20.0; //detector height - const double rot_y = -0.029; //rotation angle - for (int i = 0; i < rp2ndDetNr; i++) - { - auto *detRP_2nd = new EICG4RPSubsystem(Form("rpTruth2_%d", 2 * i), 2 * i); - detRP_2nd->SuperDetector("rpTruth2"); - 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] - hFarFwdBeamLine::enclosure_center); - detRP_2nd->set_double_param("rot_y", AngleFlip(rot_y * TMath::RadToDeg())); - detRP_2nd->set_double_param("rp_x", rp_x); - detRP_2nd->set_double_param("rp_y", rp_y); - detRP_2nd->set_double_param("hole_x", hole_x); - detRP_2nd->set_double_param("hole_y", hole_y); - detRP_2nd->set_double_param("length", rpSi_zLen); - detRP_2nd->set_string_param("material", "G4_Si"); - detRP_2nd->set_double_param("detid", 2 * i); - detRP_2nd->set_double_param("pipe_x", rppipe_x); - detRP_2nd->set_double_param("pipe_y", 0); - detRP_2nd->set_double_param("pipe_z", 0); - detRP_2nd->OverlapCheck(overlapCheck); - detRP_2nd->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); - detRP_2nd->SetActive(true); - if (verbosity) detRP_2nd->Verbosity(verbosity); - g4Reco->registerSubsystem(detRP_2nd); - } - } + string paramFile = string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP8.dat"; + int Nlayers = GetParameterFromFile (paramFile, "Number_layers"); + + for( int layer = 0; layer < Nlayers; layer++ ) { + auto *detRP = new EICG4RPSubsystem(Form("rpTruth_%d", layer), layer); + detRP->SuperDetector("rpTruth"); + detRP->SetParameterFile( paramFile ); + detRP->set_double_param("FFenclosure_center", hFarFwdBeamLine::enclosure_center ); + detRP->set_int_param("layerNumber", layer + 1); + + detRP->OverlapCheck(overlapCheck); + detRP->SetMotherSubsystem(hFarFwdBeamLine::hFarFwdBeamLineEnclosure); + detRP->SetActive(true); + if (verbosity) detRP->Verbosity(verbosity); + g4Reco->registerSubsystem(detRP); + } } if (verbosity > 0) @@ -1211,55 +885,106 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) //--------------------------------- // B0 implementation // Three choices: 1. Realistic detector; 2. Circulat plane; 3. hit plane with realistic detector goemetry - - - if (Enable::B0_DISABLE_HITPLANE) { - - // Choice 1 realistic detector - + double b0tr_z = 0; //Subsystem position relative to B0 magnet (for iterator) const int b0DetNr = 4; const double b0Mag_zCent = 610; const double b0Mag_zLen = 120; + const double b0tr[4]={10,40,70,100}; const double b0Cu_zLen = .2; //B0 dead material length const double b0Si_zLen = .1; //B0 Si length - const double b0Ecal_zLen = 20.; //B0 Ecal length - const double pipe_hole = 6.0; //detector cut off for beam pipe //- different from IP6 to account for the 35mrad angle - const double pipe_x = -1; //pipe hole position //- different from IP6 to account for the 35mrad angle - const double d_radius = 9.0; //detector cut off Packman //- different from IP6 to account for the 35mrad angle - const double b0_radius = 24.5; //outer radius of B0-detector //- different from IP6 + const double b0Ecal_zLen = 10; //B0 Ecal length + double pipe_hole_r = 3.5; //detector cut off for beam pipe + double pipe_hole = 2.5; + const double cable_hole = 2.0; + const double cable_x = 21.5; + double pipe_x = -1.; //pipe hole position + const double d_radius = 7.0; //detector cut off Packman + const double b0_radius = 23.5; //outer radius of B0-detector + const double b0_magradius = 24.5; //inner radius of B0-magnet const double spanning_angle = 240; //spanning angle Packman - const double b0Ecal_z = 48; - const double tower_size = 2; //Tower size for B0 ECal in cm - double start_angle = -120; //start angle Packman //- mirrored wrt to IP6 - + const double b0Ecal_z = 48;//B0 ECal position (relative to the B0-magnet) + double start_angle = -120; //start angle Packman + const double cross_angle = 0.035; + + if (Enable::B0_DISABLE_HITPLANE) { + + // Choice 1 realistic detector +// const double b0tr[4]={10,45,80,115}; + //const double b0tr[4]={0,30,60,90}; + //const double b0tr[5]={0,25,50,75,100}; + cout << "Realistic B0"<SuperDetector("b0Truth"); + if (Enable::B0_VAR_PIPE_HOLE){ + pipe_hole = b0tr[i]*cross_angle; + pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else if (Enable::B0_CIRCLE_PIPE_HOLE){ + pipe_hole = 0.1; + pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2; + pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else { + pipe_hole = b0tr[b0DetNr-1]*cross_angle; + pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + cout <<"Starting B0 Tracker layer "<SuperDetector(Form("b0Truth_%d", i)); detB0->set_double_param("place_x", 0); detB0->set_double_param("place_y", 0); // detB0->set_int_param("ispipe", 0); //for future pipe implementation detB0->set_double_param("pipe_hole", pipe_hole); + detB0->set_double_param("cable_hole", cable_hole); detB0->set_double_param("outer_radius", b0_radius); detB0->set_double_param("d_radius", d_radius); detB0->set_double_param("length", b0Si_zLen); detB0->set_string_param("material", "G4_Si"); - detB0->set_double_param("detid",2*i); detB0->set_double_param("startAngle",start_angle); detB0->set_double_param("spanningAngle",spanning_angle); + detB0->set_double_param("detid",i); detB0->set_double_param("pipe_x", pipe_x); detB0->set_double_param("pipe_y", 0); detB0->set_double_param("pipe_z", 0); - detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); // relative to B0 magnet + detB0->set_double_param("pipe_hole_r", pipe_hole_r); + detB0->set_double_param("cable_x", cable_x); + detB0->set_double_param("cable_y", 0); + detB0->set_double_param("cable_z", 0); + detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet detB0->SetActive(true); if (verbosity) - detB0->Verbosity(verbosity); + detB0->Verbosity(verbosity); detB0->OverlapCheck(overlapCheck); detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); - - auto *detB0e = new EICG4B0Subsystem(Form("b0Truth_%d", 2*i+1), 2*i+1); - detB0e->SuperDetector("b0Truth"); +// For B0 Tracking Implementation + if (Enable::B0TRACKING){ + if (B0TRACKING::FastKalmanFilter) + { + B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); + } + } + + auto *detB0e = new EICG4B0Subsystem(Form("b0Dead_%d", i), i); + detB0e->SuperDetector("b0Dead"); // detB0e->set_int_param("ispipe", 0); //for future pipe implementation detB0e->set_double_param("pipe_hole", pipe_hole); detB0e->set_double_param("place_x", 0); @@ -1268,26 +993,41 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) detB0e->set_double_param("pipe_x", pipe_x); detB0e->set_double_param("pipe_y", 0); detB0e->set_double_param("pipe_z", 0); + detB0e->set_double_param("pipe_hole_r", pipe_hole_r); + detB0e->set_double_param("cable_x", cable_x); + detB0e->set_double_param("cable_y", 0); + detB0e->set_double_param("cable_z", 0); detB0e->set_double_param("outer_radius", b0_radius); detB0e->set_double_param("length", b0Cu_zLen); detB0e->set_string_param("material", "G4_Cu"); - detB0e->set_double_param("detid",2*i+1); + detB0e->set_double_param("detid",i); detB0e->set_double_param("startAngle",start_angle); detB0e->set_double_param("spanningAngle",spanning_angle); - detB0e->set_double_param("place_z", (b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2) +(b0Cu_zLen+b0Si_zLen)/2) ); // relative to B0 magnet - detB0e->SetActive(true); + detB0e->set_double_param("place_z", b0tr_z +(b0Cu_zLen+b0Si_zLen)/2) ; // relative to B0 magnet + detB0e->SetActive(false); if (verbosity) detB0e->Verbosity(verbosity); detB0e->OverlapCheck(overlapCheck); detB0e->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0e); - - + } + if (Enable::B0ECAL) { + pipe_hole = b0Mag_zLen*cross_angle; + pipe_x = cross_angle*b0Mag_zCent - hFarFwdBeamLine::B0Magnet_x; + if (Enable::B0_CIRCLE_PIPE_HOLE){ + pipe_hole = 0.1; + pipe_hole_r = pipe_hole_r + b0Mag_zLen*cross_angle/2; + } + cout <<"Starting B0 ECAL "<SetTowerMappingFile(mapping_b0ecal.str()); @@ -1299,14 +1039,20 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) B0Ecal->set_double_param("pipe_x", pipe_x); B0Ecal->set_double_param("pipe_y", 0); B0Ecal->set_double_param("pipe_z", 0); + B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecal->set_double_param("cable_x", cable_x); + B0Ecal->set_double_param("cable_y", 0); + B0Ecal->set_double_param("cable_z", 0); B0Ecal->set_double_param("length", b0Ecal_zLen); B0Ecal->set_double_param("outer_radius", b0_radius); B0Ecal->set_double_param("d_radius", d_radius); B0Ecal->set_string_param("material", "G4_PbWO4"); B0Ecal->set_double_param("startAngle",start_angle); B0Ecal->set_double_param("spanningAngle",spanning_angle); - B0Ecal->set_double_param("tower_size",tower_size); B0Ecal->set_double_param("detid",0); + B0Ecal->set_double_param("global_x",hFarFwdBeamLine::B0Magnet_x); + B0Ecal->set_double_param("global_y",hFarFwdBeamLine::B0Magnet_y); + B0Ecal->set_double_param("global_z",hFarFwdBeamLine::B0Magnet_z); B0Ecal->set_int_param("lightyield",1); //Note additional parameter for storing Light Yield in B0 Ecal B0Ecal->SetActive(true); if (verbosity) @@ -1325,6 +1071,10 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) B0Ecal->set_double_param("pipe_x", pipe_x); B0Ecal->set_double_param("pipe_y", 0); B0Ecal->set_double_param("pipe_z", 0); + B0Ecal->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecal->set_double_param("cable_x", cable_x); + B0Ecal->set_double_param("cable_y", 0); + B0Ecal->set_double_param("cable_z", 0); B0Ecal->set_double_param("length", b0Ecal_zLen); B0Ecal->set_double_param("outer_radius", b0_radius); B0Ecal->set_double_param("d_radius", d_radius); @@ -1340,8 +1090,8 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) g4Reco->registerSubsystem(B0Ecal); } - auto *B0Ecale = new EICG4B0Subsystem(Form("b0Truth_%d", 2*b0DetNr+1), 2*b0DetNr+1); //B0 ECal dead layer is the same subsystem as other four dead layers - B0Ecale->SuperDetector("b0Truth"); + auto *B0Ecale = new EICG4B0Subsystem(Form("b0Dead_%d", b0DetNr), b0DetNr); //B0 ECal dead layer is the same subsystem as other four dead layers + B0Ecale->SuperDetector("b0Dead"); // B0Ecale->set_int_param("ispipe", 0); //for future pipe implementation B0Ecale->set_double_param("pipe_hole", pipe_hole); B0Ecale->set_double_param("place_x", 0); @@ -1350,47 +1100,43 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) B0Ecale->set_double_param("pipe_x", pipe_x); B0Ecale->set_double_param("pipe_y", 0); B0Ecale->set_double_param("pipe_z", 0); + B0Ecale->set_double_param("pipe_hole_r", pipe_hole_r); + B0Ecale->set_double_param("cable_x", cable_x); + B0Ecale->set_double_param("cable_y", 0); + B0Ecale->set_double_param("cable_z", 0); B0Ecale->set_double_param("length", b0Cu_zLen); B0Ecale->set_double_param("d_radius", d_radius); B0Ecale->set_double_param("outer_radius", b0_radius); B0Ecale->set_string_param("material", "G4_Cu"); B0Ecale->set_double_param("startAngle",start_angle); B0Ecale->set_double_param("spanningAngle",spanning_angle); - B0Ecale->set_double_param("detid",2*b0DetNr+1); - B0Ecale->SetActive(true); + B0Ecale->set_double_param("detid",b0DetNr+1); + //B0Ecale->SetActive(true); + B0Ecale->SetActive(false); if (verbosity) B0Ecale->Verbosity(verbosity); B0Ecale->OverlapCheck(overlapCheck); B0Ecale->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(B0Ecale); - } - } - + } } else { if (Enable::B0_FULLHITPLANE) { // Choice 2 circular hit planes + cout << "Circular hit planes"<SuperDetector("b0Truth"); + detB0->SuperDetector("b0Truth"); + //detB0->SuperDetector(Form("b0Truth_%d", i)); 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", b0tr_z); // relative to B0 magnet detB0->SetActive(true); if (verbosity) detB0->Verbosity(verbosity); detB0->OverlapCheck(overlapCheck); @@ -1398,6 +1144,28 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); + if (Enable::B0TRACKING){ + if (B0TRACKING::FastKalmanFilter) + { + B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); + } + } } @@ -1405,52 +1173,82 @@ void hFarFwdDefineDetectorsIP8(PHG4Reco *g4Reco) /// Fun4All default B0 planes /// Choice 3 Hit planes with real detector geometry - - const int b0DetNr = 4; - const double b0Mag_zCent = 610; - const double b0Mag_zLen = 120; - const double b0Cu_zLen = .2; //B0 dead material length - const double b0Si_zLen = .1; //B0 Si length - const double b0Ecal_zLen = 20.; //B0 Ecal length - const double pipe_hole = 6.0; //detector cut off for beam pipe - const double pipe_x = -5.1; //pipe hole position - const double d_radius = 9.0; //detector cut off Packman - const double b0_radius = 24.5; //outer radius of B0-detector - const double spanning_angle = 240; //spanning angle Packman - const double b0Ecal_z = 48; - double start_angle = -120; //start angle Packman + cout << "Realistic hit planes"<SuperDetector("b0Truth"); - detB0->set_double_param("place_x", 0); - detB0->set_double_param("place_y", 0); - // detB0->set_int_param("ispipe", 0); //for future pipe implementation - detB0->set_double_param("pipe_hole", pipe_hole); - detB0->set_double_param("outer_radius", b0_radius); - detB0->set_double_param("d_radius", d_radius); - detB0->set_double_param("length", b0Si_zLen); - detB0->set_string_param("material", "G4_Si"); - detB0->set_double_param("detid",2*i); - detB0->set_double_param("startAngle",start_angle); - detB0->set_double_param("spanningAngle",spanning_angle); - detB0->set_double_param("pipe_x", pipe_x); - detB0->set_double_param("pipe_y", 0); - detB0->set_double_param("pipe_z", 0); - detB0->set_double_param("place_z", b0Mag_zLen / (b0DetNr + 1) * (i - b0DetNr / 2)); // relative to B0 magnet + if (Enable::B0_VAR_PIPE_HOLE){ + pipe_hole = b0tr[i]*cross_angle; + pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[i]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else if (Enable::B0_CIRCLE_PIPE_HOLE){ + pipe_hole = 0.1; + pipe_hole_r = pipe_hole_r + b0tr[b0DetNr-1]*cross_angle/2; + pipe_x = - cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + else { + pipe_hole = b0tr[b0DetNr-1]*cross_angle; + pipe_x = cross_angle*(b0Mag_zCent - b0Mag_zLen/2 + b0tr[b0DetNr-1]/2) - hFarFwdBeamLine::B0Magnet_x; + } + cout <<"Starting B0 Tracker layer "<SuperDetector(Form("b0Truth_%d", i)); + detB0->set_double_param("place_x", 0); + detB0->set_double_param("place_y", 0); + // detB0->set_int_param("ispipe", 0); //for future pipe implementation + detB0->set_double_param("pipe_hole", pipe_hole); + detB0->set_double_param("cable_hole", cable_hole); + detB0->set_double_param("outer_radius", b0_radius); + detB0->set_double_param("d_radius", d_radius); + detB0->set_double_param("length", b0Si_zLen); + detB0->set_string_param("material", "G4_Si"); + detB0->set_double_param("startAngle",start_angle); + detB0->set_double_param("spanningAngle",spanning_angle); + detB0->set_double_param("detid",i); + detB0->set_double_param("pipe_x", pipe_x); + detB0->set_double_param("pipe_y", 0); + detB0->set_double_param("pipe_z", 0); + detB0->set_double_param("pipe_hole_r", pipe_hole_r); + detB0->set_double_param("cable_x", cable_x); + detB0->set_double_param("cable_y", 0); + detB0->set_double_param("cable_z", 0); + detB0->set_double_param("place_z", b0tr_z); // relative to B0 magnet detB0->SetActive(true); if (verbosity) detB0->Verbosity(verbosity); detB0->OverlapCheck(overlapCheck); detB0->SetMotherSubsystem(hFarFwdBeamLine::B0Magnet); g4Reco->registerSubsystem(detB0); + if (Enable::B0TRACKING){ + if (B0TRACKING::FastKalmanFilter) + { + B0TRACKING::FastKalmanFilter->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilter->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::FastKalmanFilterB0Track->add_phg4hits(string("G4HIT_") + Form("b0Truth_%d", i) , // const std::string& phg4hitsNames, + B0TrackFastSim::Vertical_Plane, // const DETECTOR_TYPE phg4dettype, + G4B0TRACKING::PositionResolution, // const float radres, + G4B0TRACKING::PositionResolution, // const float phires, + 0, // const float lonres, *ignored in plane detector* + 1, // const float eff, + 0); // const float noise + B0TRACKING::FastKalmanFilterB0Track->add_zplane_state(Form("b0Truth_%d", i), b0Mag_zCent+b0tr_z); + B0TRACKING::B0ProjectionNames.insert(Form("b0Truth_%d", i)); + } + } } } } } + void hFarFwdDefineBeamPipe(PHG4Reco *g4Reco) { int verbosity = std::max(Enable::VERBOSITY, Enable::HFARFWD_VERBOSITY); @@ -1610,5 +1408,4 @@ void FFR_Eval(const std::string &outputfile) } - #endif diff --git a/common/GlobalVariables.C b/common/GlobalVariables.C index 064d68cc..ad012a4b 100644 --- a/common/GlobalVariables.C +++ b/common/GlobalVariables.C @@ -37,6 +37,11 @@ namespace Enable // IP selection require explicit choice in the main macros bool IP6 = false; bool IP8 = false; + + float HFARFWD_ION_ENERGY = 0; + float HFARBWD_E_ENERGY = 0; + TString BEAM_COLLISION_SETTING; + } // namespace Enable // every G4 subsystem needs to implement this @@ -71,6 +76,19 @@ namespace TRACKING std::set ProjectionNames; } +//For B0 Tracking +class B0TrackFastSim; +namespace B0TRACKING +{ + string TrackNodeName = "TrackMap"; + + B0TrackFastSim * FastKalmanFilter(nullptr); + + B0TrackFastSim * FastKalmanFilterB0Track(nullptr); + + std::set B0ProjectionNames; +} + namespace G4MAGNET { // initialize to garbage values - the override is done in the respective diff --git a/detectors/EICDetector/Fun4All_G4_EICDetector.C b/detectors/EICDetector/Fun4All_G4_EICDetector.C index 76bd272e..6a4b1f11 100644 --- a/detectors/EICDetector/Fun4All_G4_EICDetector.C +++ b/detectors/EICDetector/Fun4All_G4_EICDetector.C @@ -67,11 +67,38 @@ int Fun4All_G4_EICDetector( // switching IPs by comment/uncommenting the following lines // used for both beamline setting and for the event generator crossing boost Enable::IP6 = true; - // Enable::IP8 = true; + //Enable::IP8 = true; + + + //=============== + // The following Ion energy and electron energy setting needs to be speficied + // The setting options for e-p high divergence setting (p energy x e energy): + // Option: 275x18, 275x10, 100x10, 100x5, 41x5 + // + // The setting options for e-p high divergence setting (p energy x e energy): + // Option: 275x18, 275x10, 100x10, 100x5, 41x5 + // + // The setting options for e-p high divergence setting (A energy x e energy): + // Option: 110x18, 110x10, 110x5, 41x5 // Setting proton beam pipe energy. If you don't know what to set here, leave it at 275 Enable::HFARFWD_ION_ENERGY = 275; + // Setting electron beam pipe energy. If you don't know what to set here, leave it at 18 + Enable::HFARBWD_E_ENERGY = 18; + + // Beam Scattering configuration setting specified by CDR + // + // Option 1: ep-high-acceptance + // Option 2: ep-high-divergence + // Option 3: eA + // + // Enable::BEAM_COLLISION_SETTING = "ep-high-divergence"; + // If you don't know what to put here, set it to ep-high-divergence + // + // Enable::BEAM_COLLISION_SETTING = "eA"; + Enable::BEAM_COLLISION_SETTING = "ep-high-divergence"; + // Either: // read previously generated g4-hits files, in this case it opens a DST and skips // the simulations step completely. The G4Setup macro is only loaded to get information @@ -111,6 +138,9 @@ int Fun4All_G4_EICDetector( // Input::GUN_NUMBER = 3; // if you need 3 of them // Input::GUN_VERBOSITY = 0; + // Particle ion gun + // Input::IONGUN = true; + // Upsilon generator // Input::UPSILON = true; // Input::UPSILON_NUMBER = 3; // if you need 3 of them @@ -183,6 +213,20 @@ int Fun4All_G4_EICDetector( INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0); INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0); } + + if (Input::IONGUN) + { + float theta = -25e-3; + + INPUTGENERATOR::IonGun[0]->SetA(197); + INPUTGENERATOR::IonGun[0]->SetZ(79); + INPUTGENERATOR::IonGun[0]->SetCharge(79); + INPUTGENERATOR::IonGun[0]->SetMom(sin(theta)*110*197, 0,cos(theta)*110*197); // -25 mrad + + INPUTGENERATOR::IonGun[0]->Print(); + + } + // pythia6 if (Input::PYTHIA6) { @@ -251,7 +295,7 @@ int Fun4All_G4_EICDetector( // Enable::DSTREADER = true; // turn the display on (default off) - // Enable::DISPLAY = true; + // Enable::DISPLAY = true; //====================== // What to run @@ -374,9 +418,17 @@ int Fun4All_G4_EICDetector( // Enable::ZDC_DISABLE_BLACKHOLE = true; // B0 + // B0 shape // Enable::B0_DISABLE_HITPLANE = true; // Enable::B0_FULLHITPLANE = true; + // Enable::B0_VAR_PIPE_HOLE = true; + // Enable::B0_CIRCLE_PIPE_HOLE = false; + + // B0 Tracking + // Enable::B0TRACKING = false; // For B0 Tracking + // Enable::B0TRACKING_EVAL = Enable::B0TRACKING && true; //For B0 Tracking + // B0 ECAL // Enable::B0ECALTOWERS = true; //To Construct Towers of B0ECal instead of one single volume // Enable::B0ECAL = Enable::B0_DISABLE_HITPLANE && true; // Enable::B0ECAL_CELL = Enable::B0ECAL && true; @@ -386,12 +438,19 @@ int Fun4All_G4_EICDetector( // RP // Enable::RP_DISABLE_HITPLANE = true; - // Enable::RP_FULLHITPLANE = true; + + //Far Backward detectors +// Enable::BWD = true; +// Enable::BWDN[0]=true; // 1st Q2 tagger +// Enable::BWDN[1]=true; // 2nd Q2 tagger +// Enable::BWDN[2]=true; // 1st Lumi +// Enable::BWDN[3]=true; // 2nd Lumi (+) +// Enable::BWDN[4]=true; // 3rd Lumi (-) +// Enable::BWD_CELL = Enable::BWD && true; +// Enable::BWD_TOWER = Enable::BWD_CELL && true; +// Enable::BWD_CLUSTER = Enable::BWD_TOWER && true; +// Enable::BWD_EVAL = Enable::BWD_CLUSTER && true; - // RP after 2nd focus for IP8 only - // Enable::RP2nd_DISABLE_HITPLANE = true; - // Enable::RP2nd_FULLHITPLANE = true; - //************************************************************ // details for calos: cells, towers, clusters //************************************************************ @@ -450,6 +509,7 @@ int Fun4All_G4_EICDetector( //--------------- // G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo // G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans + // G4WORLD::WorldMaterial = "G4_Galactic"; // set to G4_GALACTIC for material scans //--------------- // Magnet Settings @@ -536,6 +596,9 @@ int Fun4All_G4_EICDetector( if (Enable::B0ECAL_TOWER) B0ECAL_Towers(); // For B0Ecal if (Enable::B0ECAL_CLUSTER) B0ECAL_Clusters(); //For B0Ecal + + if (Enable::BWD_TOWER) BWD_Towers(); // For Bwd + if (Enable::BWD_CLUSTER) BWD_Clusters(); //For Bwd if (Enable::DSTOUT_COMPRESS) ShowerCompress(); @@ -604,6 +667,8 @@ int Fun4All_G4_EICDetector( if (Enable::B0ECAL_EVAL) B0ECAL_Eval(outputroot + "_g4b0ecal_eval_test.root"); // For B0Ecal + if (Enable::BWD_EVAL) BWD_Eval(outputroot + "_g4bwd_eval_e0100_debug"); // For Bwd + //if (Enable::FWDJETS_EVAL) Jet_FwdEval(); if (Enable::USER) UserAnalysisInit(); diff --git a/detectors/EICDetector/G4Setup_EICDetector.C b/detectors/EICDetector/G4Setup_EICDetector.C index a39e12bb..f2f00634 100644 --- a/detectors/EICDetector/G4Setup_EICDetector.C +++ b/detectors/EICDetector/G4Setup_EICDetector.C @@ -21,6 +21,7 @@ #include #include #include //for B0 ECAL +#include //for Far Backward Detectors #include #include #include @@ -35,7 +36,7 @@ #include #include #include -//#include for B0 Tracking +#include //for B0 Tracking #include #include #include @@ -53,6 +54,7 @@ #include #include #include +#include #include @@ -128,7 +130,7 @@ void G4Init() if (Enable::PIPE) PipeInit(); if (Enable::PLUGDOOR) PlugDoorInit(); if (Enable::TRACKING) TrackingInit(); -// if (Enable::B0TRACKING) B0TrackingInit(); + if (Enable::B0TRACKING) B0TrackingInit(); //Farforward/backward if (Enable::HFARFWD_MAGNETS) hFarBwdBeamLineInit(); //Shouldnt this be far backward enables @@ -166,6 +168,7 @@ void G4Init() if (Enable::EHCAL) EHCALInit(); if (Enable::mRICH) mRICHInit(); if (Enable::ETOF) ETOFInit(); + if (Enable::BWD) BWDInit(); //Combined if (Enable::FST) FST_Init(); @@ -394,6 +397,13 @@ void ShowerCompress() compress->AddTowerContainer("TOWER_RAW_B0ECAL"); compress->AddTowerContainer("TOWER_CALIB_B0ECAL"); + compress->AddHitContainer("G4HIT_BWD"); + compress->AddHitContainer("G4HIT_ABSORBER_BWD"); + compress->AddCellContainer("G4CELL_BWD"); + compress->AddTowerContainer("TOWER_SIM_BWD"); + compress->AddTowerContainer("TOWER_RAW_BWD"); + compress->AddTowerContainer("TOWER_CALIB_BWD"); + se->registerSubsystem(compress); return; @@ -412,6 +422,7 @@ void DstCompress(Fun4AllDstOutputManager *out) // out->StripNode("G4HIT_ZDC"); // out->StripNode("G4HIT_RomanPots"); out->StripNode("G4HIT_b0Truth"); + out->StripNode("G4HIT_BWD"); out->StripNode("G4HIT_SVTXSUPPORT"); out->StripNode("G4HIT_CEMC_ELECTRONICS"); out->StripNode("G4HIT_CEMC"); @@ -453,6 +464,8 @@ void DstCompress(Fun4AllDstOutputManager *out) out->StripNode("G4HIT_ABSORBER_EHCAL"); out->StripNode("G4CELL_EHCAL"); out->StripNode("G4CELL_B0ECAL"); + out->StripNode("G4CELL_BWD"); + } } #endif