diff --git a/CMakeLists.txt b/CMakeLists.txt index 28cd2fb3f5dd47b5f28ef6a34dd24486b6c76b8f..5b87a14e8fe1f9958a451015a671a3d016f3880a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,10 +44,10 @@ find_package(fmt REQUIRED) find_package(spdlog REQUIRED) find_package(range-v3 REQUIRED) -#add_subdirectory(src/runplan) add_subdirectory(src/monitor) add_subdirectory(src/calibration) add_subdirectory(src/epics) +#add_subdirectory(src/runplan) #add_subdirectory(replay) # ------------------------- diff --git a/src/runplan/CMakeLists.txt b/src/runplan/CMakeLists.txt index 4b74142a7ebd0ff4ee00e29aa841379b0859a307..b9e7a4447c08918671c40612c47667167471b83e 100644 --- a/src/runplan/CMakeLists.txt +++ b/src/runplan/CMakeLists.txt @@ -1,73 +1,55 @@ cmake_minimum_required (VERSION 3.1) -find_package(InSANE QUIET) +add_library(RunPlan SHARED + src/arrange_tables.cxx + src/RunStatus.cxx + src/RunTable.cxx + src/hallc_settings.cxx + ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict.cxx + ) -#add_library(RunPlan SHARED -# src/arrange_tables.cxx -# #src/arrange_tables2.cxx -# src/RunStatus.cxx -# src/RunTable.cxx -# src/hallc_settings.cxx -# #$<$<BOOL:${InSANE_FOUND}>: src/sidis_settings_rate.cxx > -# $<$<BOOL:${InSANE_FOUND}>: src/sidis_rates.cxx > -# ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict.cxx -# ) -# -#message( "InSANE_FOUND : ${InSANE_FOUND}") -#set(lib_HEADERS -# include/runplan/arrange_tables.h -# #include/runplan/sidis_settings_rate.h -# include/runplan/sidis_rates.h -# include/runplan/RunStatus.h -# include/runplan/RunTable.h -# ) -# -#message("INSANE_FOUND : ${InSANE_INCLUDE_DIRS}") -# -#if(InSANE_FOUND) -#get_target_property(insane_include_dir InSANE::InSANEnew_xsec INTERFACE_INCLUDE_DIRECTORIES) -#endif(InSANE_FOUND) -# -## build root dictionaries -#root_generate_dictionary(RunPlanDict -# #-I${insane_include_dir} -# -I${CMAKE_CURRENT_SOURCE_DIR}/include -# -I${CMAKE_CURRENT_SOURCE_DIR} ${lib_HEADERS} -# LINKDEF include/LinkDef.h -# OPTIONS -p) -#add_custom_target(RunPlan_ROOTDICTS DEPENDS ${lib_HEADERS} ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict.cxx) -#add_dependencies(RunPlan RunPlan_ROOTDICTS) -# -#target_include_directories(RunPlan -# PUBLIC -# $<INSTALL_INTERFACE:include/runplan> -# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include/runplan> -# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include> -# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> -# PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src -# PUBLIC ${EIGEN3_INCLUDE_DIRS} -# PUBLIC $<${InSANE_FOUND}:${insane_include_dir}> -#) -#target_compile_features(RunPlan PUBLIC cxx_std_17) -# -#target_link_libraries(RunPlan -# PUBLIC fmt::fmt -# ${CXX_FILESYSTEM_LIB} -# #$<$<BOOL:${InSANE_FOUND}>:PUBLIC InSANE::InSANEbase InSANE::InSANEnew_xsec> -#) -# -# -#install(DIRECTORY include/runplan -# DESTINATION include -# ) -#install(FILES ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict_rdict.pcm -# DESTINATION ${CMAKE_INSTALL_LIBDIR} -# ) -# -#include(GNUInstallDirs) -#install(TARGETS RunPlan -# EXPORT RunPlannerTargets -# LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} -# ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} -#) +set(lib_HEADERS + include/runplan/arrange_tables.h + include/runplan/RunStatus.h + include/runplan/RunTable.h + ) + +root_generate_dictionary(RunPlanDict + -I${CMAKE_CURRENT_SOURCE_DIR}/include + -I${CMAKE_CURRENT_SOURCE_DIR} ${lib_HEADERS} + LINKDEF include/LinkDef.h + OPTIONS -p) +add_custom_target(RunPlan_ROOTDICTS DEPENDS ${lib_HEADERS} ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict.cxx) +add_dependencies(RunPlan RunPlan_ROOTDICTS) + +target_include_directories(RunPlan + PUBLIC + $<INSTALL_INTERFACE:include/runplan> + $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include/runplan> + $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include> + $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src + PUBLIC ${EIGEN3_INCLUDE_DIRS} +) +target_compile_features(RunPlan PUBLIC cxx_std_17) + +target_link_libraries(RunPlan + PUBLIC fmt::fmt + ${CXX_FILESYSTEM_LIB} +) + + +install(DIRECTORY include/runplan + DESTINATION include + ) +install(FILES ${CMAKE_CURRENT_BINARY_DIR}/RunPlanDict_rdict.pcm + DESTINATION ${CMAKE_INSTALL_LIBDIR} + ) + +include(GNUInstallDirs) +install(TARGETS RunPlan + EXPORT HallcToolsTargets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} +) diff --git a/src/runplan/include/runplan/sidis_rates.h b/src/runplan/include/runplan/sidis_rates.h deleted file mode 100644 index 615e2c3ec581ed303eb3e12fb275fd261a736f4f..0000000000000000000000000000000000000000 --- a/src/runplan/include/runplan/sidis_rates.h +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef hallc_tools_runplan_settings_rates_HH -#define hallc_tools_runplan_settings_rates_HH - -#include <chrono> -#include <cmath> -#include <fstream> -#include <future> -#include <iostream> -#include <string> -#include <thread> -#include <typeinfo> - - -#include "hallc_settings.h" -#include "CSV_settings.h" - -#include "RunTable.h" - - -namespace hallc { - - /** Builds a run-plan table etnry for a single kinematic. - * - */ - RunPlanTableEntry build_table_entry(const Kinematic& kine); - - /** Builds a table of - * - */ - RunTable build_sidis_table(const csv::CSV_Settings& settings); - - /** - * - */ - void save_tables(const RunTable& table_rows, std::string table_name, - const std::ios_base::openmode& mode = std::ios_base::app); - - void save_json_tables(const RunTable& table_rows, std::string table_name, - const std::ios_base::openmode& mode = std::ios_base::app); - - /** - * - */ - std::vector<double> sidis_config_xs_and_rates(RunPlanTableEntry& run_setting); - -} - -#endif diff --git a/src/runplan/src/LH2_settings_rate.cxx b/src/runplan/src/LH2_settings_rate.cxx deleted file mode 100644 index 3bdde90dbc2bc30f26c8a0cca77e0f0b0935c844..0000000000000000000000000000000000000000 --- a/src/runplan/src/LH2_settings_rate.cxx +++ /dev/null @@ -1,698 +0,0 @@ -#include <iostream> -#include <fstream> -#include <string> -#include <cmath> -#include <thread> -#include <future> -#include <iostream> -#include <string> -#include <chrono> - -#include <fmt/core.h> -#include <fmt/ostream.h> -R__LOAD_LIBRARY(libfmt.so) - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" - -#include "ROOT/TFuture.hxx" - -#include "InSANE/MaterialProperties.h" -#include "InSANE/Luminosity.h" -#include "PhaseSpaceVariables.h" -#include "PSSampler.h" -#include "NFoldDifferential.h" -#include "Jacobians.h" -#include "DiffCrossSection.h" -#include "FinalState.h" -#include <typeinfo> - - -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" - -#include "InSANE/Stat2015_UPDFs.h" -#include "DSSFragmentationFunctions.h" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" -#include "TBufferJSON.h" - -#include "ROOT/RDataFrame.hxx" - -#include "hallc_settings.h" - -using namespace ROOT; -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; -using namespace ROOT::Experimental; - - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - -std::vector<double> sidis_config_sigma_rate( RunPlanTableEntry& run_setting); - -/** Calculates a table of rates. - * - */ -void LH2_settings_rate() { - - // The fortran code used for the fragmentation functions - // doesn't allow multi threaded to run without modifying - // the way it uses the common blocks - ROOT::EnableImplicitMT(1); - - std::vector<std::pair<double,std::vector<RunPlanTableEntry>>> all_tables; - std::vector<std::pair<double,std::vector<RunPlanTableEntry>>> all_kaon_tables; - - for(auto& setting_group : csv::LH2_settings) { - - // We only are interested in the LH2 target rates for the two lowest Q2 settings - bool use_LH2 = true; - //if( std::get<0>(setting_group) > 5.5 ) { use_LH2 = false; } - bool only_LH2 = true; - - std::vector<RunPlanTableEntry> run_plan_kaon_table; - std::vector<RunPlanTableEntry> run_plan_table; - for(auto& a_setting: setting_group.second) { - - a_setting.Q2 = std::get<0>(setting_group); - - auto plus_entry = ROOT::Experimental::Async([&a_setting,&use_LH2,&only_LH2](){ - std::vector<RunPlanTableEntry> resulting_entries; - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = 1; - entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - entry0.rates.total_rate = entry0.rates.window_rate + entry0.rates.LH2_rate; - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - //if( ! only_LH2 ) { - resulting_entries.push_back( entry0 ); - //} - //if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto minus_entry = ROOT::Experimental::Async([&a_setting,&use_LH2,&only_LH2](){ - std::vector<RunPlanTableEntry> resulting_entries; - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = -1; - entry0.Ibeam = 50.0; - entry0.counts = 30000; - if( entry0.kinematic.z > 0.55){ - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - entry0.rates.total_rate = entry0.rates.window_rate + entry0.rates.LH2_rate; - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - //if( ! only_LH2 ) { - resulting_entries.push_back( entry0 ); - //} - //if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto table_entry0 = plus_entry.get() ; - auto table_entry1 = minus_entry.get(); - - auto Kplus_entry = ROOT::Experimental::Async([&a_setting,&table_entry0,&use_LH2,&only_LH2](){ - std::vector<RunPlanTableEntry> resulting_entries; - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = 2; - entry0.Ibeam = 25.0; - double scale = 50.0 / entry0.Ibeam; - sidis_config_sigma_rate(entry0); - entry0.rates.total_rate = entry0.rates.window_rate + entry0.rates.LH2_rate; - entry0.time = table_entry0[0].time; //.counts/entry0.rate/(60.0*60.0); - entry0.counts = (60.0 * 60.0) * entry0.time * entry0.rates.total_rate; - entry0.LD2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LD2_rate; - entry0.LH2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LH2_rate; - entry0.window_counts = (60.0 * 60.0) * entry0.time * entry0.rates.window_rate; - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - //if( ! only_LH2 ) { - resulting_entries.push_back( entry0 ); - //} - //if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = table_entry0[0].time;//entry1.counts/entry1.rate/(60.0*60.0); - // entry1.counts = (60.0 * 60.0) * entry1.time * entry1.rate; - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto Kminus_entry = ROOT::Experimental::Async([&a_setting,&table_entry1,&use_LH2,&only_LH2](){ - std::vector<RunPlanTableEntry> resulting_entries; - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = -2; - entry0.Ibeam = 50.0; - double scale = 50.0 / entry0.Ibeam; - sidis_config_sigma_rate(entry0); - entry0.rates.total_rate = entry0.rates.window_rate + entry0.rates.LH2_rate; - entry0.time = table_entry1[0].time; //.counts/entry0.rate/(60.0*60.0); - entry0.counts = (60.0 * 60.0) * entry0.time * entry0.rates.total_rate; - entry0.LD2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LD2_rate; - entry0.LH2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LH2_rate; - entry0.window_counts = (60.0 * 60.0) * entry0.time * entry0.rates.window_rate; - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - //if( ! only_LH2 ) { - resulting_entries.push_back( entry0 ); - //} - //if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = table_entry1[0].time;//entry1.counts/entry1.rate/(60.0*60.0); - // entry1.counts = (60.0 * 60.0) * entry1.time * entry1.rate; - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - //run_plan_table.push_back( entry_plus ); - //run_plan_table.push_back( entry_minus ); - auto table_entry2 = Kplus_entry.get() ; - auto table_entry3 = Kminus_entry.get(); - run_plan_table.insert( run_plan_table.end(), table_entry0.begin(), table_entry0.end() ); - run_plan_table.insert( run_plan_table.end(), table_entry1.begin(), table_entry1.end() ); - run_plan_kaon_table.insert( run_plan_kaon_table.end(), table_entry2.begin(), table_entry2.end() ); - run_plan_kaon_table.insert( run_plan_kaon_table.end(), table_entry3.begin(), table_entry3.end() ); - - std::cout << "done\n"; - //run_plan_table.push_back( table_entry1[0] ); - //run_plan_kaon_table.push_back( Kplus_entry.get()[0] ); - //run_plan_kaon_table.push_back( Kminus_entry.get()[0] ); - - //run_plan_table.push_back( Kplus_entry.get() ); - //run_plan_table.push_back( Kminus_entry.get() ); - } - all_tables.push_back(std::make_pair(setting_group.first, run_plan_table)); - all_kaon_tables.push_back(std::make_pair(setting_group.first, run_plan_kaon_table)); - } - - { - double total_time = 0.0; - std::ofstream ofs("tables/LH2_run_plan_table.txt",std::ios_base::trunc); - std::ofstream wiki_file("tables/LH2_run_plan_table.wiki",std::ios_base::trunc); - for(const auto& atable : all_tables){ - //ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for(const auto& entry : atable.second){ - if( std::abs(entry.polarity) >= 2 ) { - continue; - } - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - - fmt::print(ofs, "{:-<79s}\n", "-"); - fmt::print(wiki_file, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) << "\n"; - wiki_file << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) << "\n"; - } - fmt::print(ofs, "{:=<79s}\n", "="); - ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - fmt::print(wiki_file, "{:=<79s}\n", "="); - wiki_file << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - std::ofstream json_ofile("tables/LH2_run_plan_table.json",std::ios_base::trunc); - json_ofile << TBufferJSON::ToJSON(&all_tables); - json_ofile.close(); - - } - - { - double total_time = 0.0; - std::ofstream ofs("tables/LH2_run_plan_kaon_table.txt",std::ios_base::trunc); - std::ofstream wiki_file("tables/LH2_run_plan_kaon_table.wiki",std::ios_base::trunc); - for(const auto& atable : all_kaon_tables){ - //ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for(const auto& entry : atable.second){ - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - fmt::print(ofs, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0); - ofs << "\n"; - } - fmt::print(ofs, "{:=<79s}\n", "="); - ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - std::ofstream json_ofile("tables/LH2_run_plan_kaon_table.json",std::ios_base::trunc); - json_ofile << TBufferJSON::ToJSON(&all_tables); - json_ofile.close(); - } -} - -/** Compute the cross section and rate for a given setting - * - */ -std::vector<double> sidis_config_sigma_rate( RunPlanTableEntry& run_setting){ - - hallc::HCKinematic& kine_set = run_setting.kinematic; - hallc::HallCSetting& hcSet = run_setting.hcSet; - double x_set = kine_set.x; - double z_set = kine_set.z; - double nu_set = kine_set.nu; - double W_set = kine_set.W; - double Wp_set = kine_set.Wp; - int polarity = run_setting.polarity; - hcSet.HMS_theta = kine_set.th_e * degree; - hcSet.SHMS_theta = kine_set.th_q * degree; - hcSet.HMS_p0 = kine_set.Ee; - hcSet.SHMS_p0 = kine_set.Ppi; - - //SHMS 4 msr dp -10 +20 - //HMS 6 msr dp +-8 - - // ------------------------------------------------------- - // Phase space variables - IPSV phi_SHMS_psv({hcSet.SHMS_phi_min(), hcSet.SHMS_phi_max()}, "phi_shms"); - IPSV phi_HMS_psv({hcSet.HMS_phi_min(), hcSet.HMS_phi_max()}, "phi_hms"); - IPSV theta_HMS_psv({hcSet.HMS_theta_min(), hcSet.HMS_theta_max()}, "theta_hms"); - IPSV theta_SHMS_psv({hcSet.SHMS_theta_min(), hcSet.SHMS_theta_max()}, "theta_shms"); - IPSV P_HMS_psv({hcSet.HMS_P_min(), hcSet.HMS_P_max()}, "P_hms"); - IPSV P_SHMS_psv({hcSet.SHMS_P_min(), hcSet.SHMS_P_max()}, "P_shms"); - - auto diff_spectrometers = make_diff(P_HMS_psv, theta_HMS_psv, phi_HMS_psv, - P_SHMS_psv, theta_SHMS_psv, phi_SHMS_psv); - auto ps_spectrometers = make_phase_space(diff_spectrometers); - //diff_spectrometers.Print(); - - // ------------------------------ - // Prepare the initial state - double P1 = 11.0; - double m1 = 0.000511; - double E1 = std::sqrt(P1*P1+m1*m1); - - double P2 = 0.0; - double m2 = M_p/GeV; - double E2 = std::sqrt(P2*P2+m2*m2); - - InitialState init_state(P1, P2, m2); - - int xs_select = 1; - PDFs pdf; - FragFs ffs; - - /** SIDIS diferential cross section. - * - * d(sigma) - * -------------------------- - * d(x)d(y)d(z)d(phi)d(phi_e) - * - * Note: "phi" is the angle between the leptonic and hadronic planes. - * See http://inspirehep.net/record/677636 for details. - * "phi_e" is the scattered electron's azimuthal angle. - * - */ - auto SIDIS_xs_func = [&](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - double Mh = M_pion/GeV; - - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - if(!(kine.isValid)) return 0.0; - - //return 1.0; - - // Gaussian k_perp and p_perp distributions for PT dependence - double kperp2_mean = 0.28; // GeV^2 - double pperp2_mean = 0.25; // GeV^2 - double PT_mean = kine.z * kine.z * kperp2_mean + pperp2_mean; - double PT2 = kine.P_hPerp * kine.P_hPerp; - auto f_PT = std::exp(-PT2 / PT_mean) / (PT_mean * pi); - auto t1 = kine.SIDIS_xs_term() * f_PT; - double FF = 0.0; - auto pdf_vals = pdf.Calculate(kine.x, kine.Q2); - - if (polarity == 1) { - // pi+ - if( xs_select == 1 || xs_select ==2){ - //proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if( xs_select == 0 || xs_select ==2){ - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -1) { - // pi- - if( xs_select == 1 || xs_select ==2){ - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if( xs_select == 0 || xs_select ==2){ - //neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == 2) { - // pi+ - if( xs_select == 1 || xs_select ==2){ - //proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if( xs_select == 0 || xs_select ==2){ - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -2) { - // pi- - if( xs_select == 1 || xs_select ==2){ - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if( xs_select == 0 || xs_select ==2){ - //neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } - FF *= kine.x; - - /// \f$ F_{UU,T} = x \sum e_a^2 f_1^a(x) D_1^a(z) \f$ - //double res = Elastic_XS({GE2, GM2}, kine); - //res = res * hbarc2_gev_nb * kine.sin_th; - return(t1*FF); - }; - - // ------------------------------------------------------- - // - auto Eprime_and_Ph_func = [](const InitialState& is, const std::array<double, 6>& x) - { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - //SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - - XYZVector l2(Polar3D<double>(kine.P_l2, kine.theta_l2, kine.phi_l2)); - XYZVector p_h(Polar3D<double>(kine.P_p2, kine.theta_p2, kine.phi_p2)); - //XYZVector p_h_lab(r1*r2*p_h); - - auto E_prime = ROOT::Math::XYZTVector(l2.x(), l2.y(), l2.z(), kine.E_l2); - auto P_hadron = ROOT::Math::XYZTVector(p_h.x(), p_h.y(), p_h.z(), kine.E_p2); - return std::array<ROOT::Math::XYZTVector,2>({E_prime,P_hadron}); - }; - - // returns SIDIS_vars using the SPEC_vars (measured by spectrometers) - auto to_sidis_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - //kine.Print(); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto to_spectrometer_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto diff_xs_sidis = make_diff_cross_section(ps_spectrometers, SIDIS_xs_func); - - auto fs_parts = make_final_state_particles({11, 211}, {0.000511, M_pion / GeV}, - Eprime_and_Ph_func); - - auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - - auto p_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - auto n_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - //auto sampler = make_ps_sampler(p_integrated_XS); - //auto sampler2 = make_ps_sampler(n_integrated_XS); - //sampler.SetFoamCells(500); - //sampler.SetFoamSample(100); - //sampler.SetFoamChat(0); - //sampler2.SetFoamCells(500); - //sampler2.SetFoamSample(100); - //sampler2.SetFoamChat(0); - - //auto evgen = make_event_generator(sampler, fs_SIDIS); - xs_select = 1; - auto total_XS_p = p_integrated_XS.CalculateTotalXS(); - //auto total_XS_p = evgen.Init(); - - //auto evgen2 = make_event_generator(sampler2, fs_SIDIS); - xs_select = 0; - auto total_XS_n = n_integrated_XS.CalculateTotalXS(); - //auto total_XS_n = evgen2.Init(); - - double total_XS = total_XS_n + total_XS_p; - - //// -------------------------- - // - //std::cout << "Total cross section = " << total_XS << " nb\n"; - int N_sim = 100000; - - using namespace insane; - double Loop1_LD2_Lumi = materials::ComputeLuminosityPerN({1, 2}, - materials::LD2_density /*g/cm3*/, - 10.0 /*cm*/, - 50.0e-6); - double Loop1_LH2_Lumi = materials::ComputeLuminosityPerN({1, 1}, - materials::LH2_density /*g/cm3*/, - 10.0 /*cm*/, - 50.0e-6); - double Loop1_entrance_window_thickness = 0.0104; - double Loop1_exit_window_thickness = 0.0133; - double Loop3_entrance_window_thickness = 0.0130; - double Loop3_exit_window_thickness = 0.0188; - - double Loop1_windows_Lumi = materials::ComputeLuminosityPerN({13, 27}, - materials::Al_density , - Loop1_entrance_window_thickness + - Loop1_exit_window_thickness , - 50.0e-6); - - TargetRates targ_rates; - - targ_rates.LD2_Lumi = Loop1_LD2_Lumi; - targ_rates.LH2_Lumi = Loop1_LH2_Lumi; - targ_rates.window_Lumi = Loop1_windows_Lumi; - - // double total_cross_section = 2*total_XS_p * nanobarn; - //double luminosity = Lumi / cm2; - targ_rates.p_XS = total_XS_p; - targ_rates.n_XS = total_XS_n; - targ_rates.LD2_XS = total_XS_n + total_XS_p; - targ_rates.LH2_XS = total_XS_p; - targ_rates.window_XS = (total_XS_n * 14 + total_XS_p * 13); - - targ_rates.LD2_rate = targ_rates.LD2_XS * targ_rates.LD2_Lumi * (nanobarn / cm2); - targ_rates.LH2_rate = targ_rates.LH2_XS * targ_rates.LH2_Lumi * (nanobarn / cm2); - targ_rates.window_rate = targ_rates.window_XS * targ_rates.window_Lumi * (nanobarn / cm2); - targ_rates.total_LD2_rate = targ_rates.window_rate + targ_rates.LD2_rate; - targ_rates.total_LH2_rate = targ_rates.window_rate + targ_rates.LH2_rate; - targ_rates.total_rate = targ_rates.window_rate + targ_rates.LD2_rate; - - std::cout << "Cross sections : \n"; - std::cout << " [p] " << targ_rates.p_XS << " nb \n"; - std::cout << " [n] " << targ_rates.n_XS << " nb \n"; - std::cout << " [LD2] " << targ_rates.LD2_XS << " nb \n"; - std::cout << " [LH2] " << targ_rates.LH2_XS << " nb \n"; - std::cout << " [windows] " << targ_rates.window_XS << " nb \n"; - //std::cout << " : " << total_cross_section/microbarn << " ub\n"; - //std::cout << " : " << total_cross_section/barn << " b\n"; - std::cout << "Luminosities : \n"; - std::cout << " [LH2] " << targ_rates.LH2_Lumi << " 1/cm2 s\n"; - std::cout << " [LD2] " << targ_rates.LD2_Lumi << " 1/cm2 s\n"; - std::cout << " [windows] " << targ_rates.window_Lumi << " 1/cm2 s\n"; - std::cout << "Rates : \n"; - std::cout << " [LD2] " << targ_rates.LD2_rate << " 1/s\n"; - std::cout << " [LH2] " << targ_rates.LH2_rate << " 1/s\n"; - std::cout << " [windows] " << targ_rates.window_rate << " 1/s\n"; - std::cout << "Total Rates : \n"; - std::cout << " [LD2] " << targ_rates.total_LD2_rate << " 1/s\n"; - std::cout << " [LH2] " << targ_rates.total_LH2_rate << " 1/s\n"; - std::cout << " " << targ_rates.total_rate << " 1/s\n"; - - run_setting.rates = targ_rates; - //run_setting.LD2_Luminosity = Loop1_LD2_Lumi; - //run_setting.window_Luminosity = Loop1_windows_Lumi; - //run_setting.LD2_sigma = total_xs_LD2; - //run_setting.window_sigma = total_xs_windows; - std::vector<double> res = { - targ_rates.LD2_XS, - targ_rates.LD2_rate, - targ_rates.total_LD2_rate, - targ_rates.total_LH2_rate, - targ_rates.LH2_XS, - targ_rates.LH2_rate - }; - return res; -} - - diff --git a/src/runplan/src/fmt_test.cxx b/src/runplan/src/fmt_test.cxx deleted file mode 100644 index a469cc13eca9abec2b23c673d76a6facf31bdb88..0000000000000000000000000000000000000000 --- a/src/runplan/src/fmt_test.cxx +++ /dev/null @@ -1,18 +0,0 @@ -#include <iostream> - -#include <fmt/core.h> - -#include <fmt/ostream.h> - -R__LOAD_LIBRARY(libfmt.so) - -void fmt_test() { - -fmt::print( - std::cout, "{:>5.3f}, {:>5.3f}\n", 1.0,2.0); - - std::cout << " soihasdfoisdoafij aosoaijweoifjaoefwoijaseofi j" << "oiajsdfoiajsidfoijsadflkj" << - "\n"; - -} - diff --git a/src/runplan/src/json_test.cxx b/src/runplan/src/json_test.cxx deleted file mode 100644 index 917482f78ce68f4530e43ee32eeca95c5105adb0..0000000000000000000000000000000000000000 --- a/src/runplan/src/json_test.cxx +++ /dev/null @@ -1,74 +0,0 @@ -#include <iostream> -#include <fstream> -#include <string> -#include <cmath> -#include <thread> -#include <future> -#include <iostream> -#include <string> -#include <chrono> - -#include <fmt/core.h> -#include <fmt/ostream.h> - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" - -#include "InSANE/MaterialProperties.h" -#include "InSANE/Luminosity.h" -#include "PhaseSpaceVariables.h" -#include "PSSampler.h" -#include "NFoldDifferential.h" -#include "Jacobians.h" -#include "DiffCrossSection.h" -#include "FinalState.h" -#include <typeinfo> - - -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" - -#include "InSANE/Stat2015_UPDFs.h" -#include "DSSFragmentationFunctions.h" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" -#include "TBufferJSON.h" - -#include "ROOT/RDataFrame.hxx" - -#include "csv_settings.h" - -using namespace ROOT; -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - -using FullTable = std::vector<std::pair<double,std::vector<RunPlanTableEntry>>>; - -void json_test() { - std::ifstream t("run_plan_table.json"); - std::stringstream buffer; - buffer << t.rdbuf(); - - FullTable* all_tables = nullptr; - TBufferJSON::FromJSON(all_tables, buffer.str().c_str()); - std::cout << all_tables->size() << "\n"; - //std::cout << json << std::endl; -} - diff --git a/src/runplan/src/sidis_pion_eg.cxx b/src/runplan/src/sidis_pion_eg.cxx deleted file mode 100644 index e29e253f266e2cca0be19b1ad3e7cc45bfc8e547..0000000000000000000000000000000000000000 --- a/src/runplan/src/sidis_pion_eg.cxx +++ /dev/null @@ -1,495 +0,0 @@ -#include <iostream> -#include <fstream> -#include <string> -#include <cmath> - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" - -#include "PhaseSpaceVariables.h" -#include "NFoldDifferential.h" -#include "Jacobians.h" -#include "DiffCrossSection.h" -#include "FinalState.h" -#include <typeinfo> - -#include "PSSampler.h" - -//#include "InSANE/DoubleSpinAsymmetries.h" -//#include "InSANE/AMTFormFactors.h" -//#include "InSANE/KellyFormFactors.h" -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" - -#include "InSANE/Stat2015_UPDFs.h" -#include "DSSFragmentationFunctions.h" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" - -#include "ROOT/RDataFrame.hxx" - -using namespace ROOT; // TDataFrame's namespace - -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; - -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - - -/** Hall C Spectrometer Settings. - * - * SHMS: - * Central Momentum: Momentum Acceptance: Momentum Resolution: - * 2 to 11 GeV/c -10% < δ < +22% 0.03%-0.08% - * Scattering Angle: Horizontal acceptance: Vertical acceptance: - * 5.5 to 40° +- 24mrad +- 40 mrad - * Solid angle acceptance: - * 4.0 msr - * Vertex Length (90deg): Vertex length res (Ytar): - * 30cm 0.1-0.3 cm - * Horizontal resolution (YPtar): Vertical resolution (XPtar): - * 0.5-1.2 mrad 0.3-1.1 mrad - * - * HMS: - * https://www.jlab.org/Hall-C/equipment/HMS/performance.html - * - * - * Coordinate system used: - * z = beam direction - * x = SHMS (beam left) - * y = up (towards hall ceiling) - * - */ -struct HallCSetting { - const double SHMS_dtheta = 0.024; // 24 mrad - const double SHMS_dphi = 0.040; // 40 mrad - const double SHMS_dP_low = 0.1; // -10% - const double SHMS_dP_high = 0.22; // 22% - const double SHMS_solid_angle = 4.0*SHMS_dtheta*SHMS_dphi; // 4 msr - const double SHMS_phi = 0.0; // SHMS sits on the +x side - const double HMS_dtheta = 0.04; // - const double HMS_dphi = 0.04; // - const double HMS_dP_low = 0.09; // -9% - const double HMS_dP_high = 0.09; // 9% - const double HMS_solid_angle = 4.0*HMS_dtheta*HMS_dphi; // 4 msr - const double HMS_phi = insane::units::pi-0.04; - - // HMS detects electron - double HMS_theta = 14.5*degree; - double HMS_p0 = 5.0; // GeV/c - // SHMS detects hadron - double SHMS_theta = 13.5*degree; - double SHMS_p0 = 3.0; // GeV/c - - double HMS_P_min() const { return HMS_p0 * (1.0 - HMS_dP_low); } - double HMS_P_max() const { return HMS_p0 * (1.0 + HMS_dP_high); } - double SHMS_P_min() const { return SHMS_p0 * (1.0 - SHMS_dP_low); } - double SHMS_P_max() const { return SHMS_p0 * (1.0 + SHMS_dP_high); } - - double HMS_phi_min() const { return HMS_phi - HMS_dphi; } - double HMS_phi_max() const { return HMS_phi + HMS_dphi; } - double SHMS_phi_min() const { return SHMS_phi - SHMS_dphi; } - double SHMS_phi_max() const { return SHMS_phi + SHMS_dphi; } - - double HMS_theta_min() const { return HMS_theta - HMS_dtheta; } - double HMS_theta_max() const { return HMS_theta + HMS_dtheta; } - double SHMS_theta_min() const { return SHMS_theta - SHMS_dtheta; } - double SHMS_theta_max() const { return SHMS_theta + SHMS_dtheta; } -}; - -void sidis_pion_eg(){ - - HallCSetting hc_setting; - - //IPSV phi_HMS_psv( {-pi, pi }, "phi_hms"); - //IPSV phi_SHMS_psv({-pi, pi }, "phi_shms"); - ////IPSV phi_SHMS_psv({hc_setting.SHMS_phi_min(), hc_setting.SHMS_phi_max()}, "phi_shms"); - //IPSV theta_HMS_psv( {hc_setting.HMS_theta_min(), hc_setting.HMS_theta_max()}, "theta_hms"); - //IPSV theta_SHMS_psv({hc_setting.SHMS_theta_min(), hc_setting.SHMS_theta_max()}, "theta_shms"); - ////IPSV P_HMS_psv( {hc_setting.HMS_P_min(), hc_setting.HMS_P_max()}, "P_hms"); - ////IPSV P_SHMS_psv({hc_setting.SHMS_P_min(), hc_setting.SHMS_P_max()}, "P_shms"); - //IPSV P_HMS_psv( {0.5,11.0}, "P_hms"); - //IPSV P_SHMS_psv({0.5,11.0}, "P_shms"); - - //auto diff_spectrometers = make_diff(P_HMS_psv, theta_HMS_psv, phi_HMS_psv, - // P_SHMS_psv, theta_SHMS_psv, phi_SHMS_psv); - //auto ps_spectrometers = make_phase_space(diff_spectrometers); - //diff_spectrometers.Print(); - - // ------------------------------------------------------- - // Phase space variables - - PSV<Invariant> x_psv ({0.01,0.99}, Invariant::x, "x"); - PSV<Invariant> y_psv ({0.01,0.99}, Invariant::y, "y"); - PSV<Invariant> z_psv ({0.01,0.99}, Invariant::z, "z"); - PSV<Invariant> Q2_psv({0.5, 10.0}, Invariant::Q2, "Q2"); - IPSV phi_psv( {-1.0*pi, pi}, "phi"); - IPSV phi_LAB_psv( {-1.0*pi, pi}, "phi_e"); - IPSV PT_LAB_psv( {0.0, 1.0}, "PT"); - - - //SHMS 4 msr dp -10 +20 - //HMS 6 msr dp +-8 - - // ------------------------------ - // Prepare the initial state - double P1 = 11.0; - double m1 = 0.000511; - double E1 = std::sqrt(P1*P1+m1*m1); - - double P2 = 0.0; - double m2 = M_p/GeV; - double E2 = std::sqrt(P2*P2+m2*m2); - - InitialState init_state(P1, P2, m2); - - auto SIDIS_xs_func = [&](const InitialState& is, const std::array<double,6>& x){ - // SIDIS diferential cross section - // - // d(sigma) - // -------------------------- - // d(x)d(y)d(z)d(phi)d(phi_e) - // - // Note: "phi" is the angle between the leptonic and hadronic planes. - // See http://inspirehep.net/record/677636 for details. - // "phi_e" is the scattered electron's azimuthal angle. - // - using namespace insane::units; - static PDFs pdf; - static FragFs ffs; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - - if(!(kine.isValid)) return 0.0; - return 1.0; - - // Gaussian k_perp and p_perp distributions for PT dependence - double kperp2_mean = 0.28; //GeV^2 - double pperp2_mean = 0.25; //GeV^2 - double PT_mean = kine.z*kine.z *kperp2_mean + pperp2_mean; - double PT2 = kine.P_hPerp*kine.P_hPerp; - auto f_PT = std::exp(-PT2/PT_mean)/(PT_mean*pi); - auto t1 = kine.SIDIS_xs_term()*f_PT; - - double FF = 0.0; - - auto pdf_vals = pdf.Calculate(kine.x, kine.Q2); - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF *= kine.x; - - /// \f$ F_{UU,T} = x \sum e_a^2 f_1^a(x) D_1^a(z) \f$ - //double res = Elastic_XS({GE2, GM2}, kine); - //res = res * hbarc2_gev_nb * kine.sin_th; - return(t1*FF*hbarc2_gev_nb); - }; - - // ------------------------------------------------------- - // - auto Eprime_and_Ph_func = [](const InitialState& is, const std::array<double, 6>& x) - { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - //SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - - XYZVector l2(Polar3D<double>(kine.P_l2, kine.theta_l2, kine.phi_l2)); - XYZVector p_h(Polar3D<double>(kine.P_p2, kine.theta_p2, kine.phi_p2)); - //XYZVector p_h_lab(r1*r2*p_h); - - auto E_prime = ROOT::Math::XYZTVector(l2.x(), l2.y(), l2.z(), kine.E_l2); - auto P_hadron = ROOT::Math::XYZTVector(p_h.x(), p_h.y(), p_h.z(), kine.E_p2); - return std::array<ROOT::Math::XYZTVector,2>({E_prime,P_hadron}); - }; - - // returns SIDIS_vars using the SPEC_vars (measured by spectrometers) - auto to_sidis_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - //kine.Print(); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto to_spectrometer_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - //auto to_sidis_tup = std::make_tuple( - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<0>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<1>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<2>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<3>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<4>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<5>(to_sidis_variables(is,x));} - // ); - - - //IPSV phi_SHMS_psv({hc_setting.SHMS_phi_min(), hc_setting.SHMS_phi_max()}, "phi_shms"); - //IPSV theta_HMS_psv( {hc_setting.HMS_theta_min(), hc_setting.HMS_theta_max()}, "theta_hms"); - //IPSV theta_SHMS_psv({hc_setting.SHMS_theta_min(), hc_setting.SHMS_theta_max()}, "theta_shms"); - //IPSV P_HMS_psv( {hc_setting.HMS_P_min(), hc_setting.HMS_P_max()}, "P_hms"); - //IPSV P_SHMS_psv({hc_setting.SHMS_P_min(), hc_setting.SHMS_P_max()}, "P_shms"); - // --------------------------------------------------- - // cross section differential in x y z PT phi_h phi_q - auto diff_dxyzPTphiphie = make_diff(x_psv, y_psv, z_psv, PT_LAB_psv, phi_psv, phi_LAB_psv); - auto ps_dxyzPTphiphie = make_phase_space( - diff_dxyzPTphiphie, - make_PSV("P_hms", {-pi, pi},//hc_setting.HMS_P_min(), hc_setting.HMS_P_max()}, - [&](const InitialState& is, const std::array<double, 6>& x) { - return std::get<0>(to_spectrometer_variables(is, x)); - }), - make_PSV("theta_hms", {hc_setting.HMS_theta_min(), hc_setting.HMS_theta_max()}, - [&](const InitialState& is, const std::array<double, 6>& x) { - return std::get<1>(to_spectrometer_variables(is, x)); - })); - //make_PSV("phi_hms", {hc_setting.HMS_phi_min(), hc_setting.HMS_phi_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<2>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("P_shms", {hc_setting.SHMS_P_min(), hc_setting.SHMS_P_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<3>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("theta_shms", {hc_setting.SHMS_theta_min(), hc_setting.SHMS_theta_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<4>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("phi_shms", {hc_setting.SHMS_phi_min(), hc_setting.SHMS_phi_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<5>(to_spectrometer_variables(is, x)); - // })); - ps_dxyzPTphiphie.Print(); - auto diff_xs_sidis = make_diff_cross_section(ps_dxyzPTphiphie, SIDIS_xs_func); - - // Transform to a function of the 2 measured momenta - //auto jm_spec_dxyzPTphiphie = make_jacobian(diff_dxyzPTphiphie, - // diff_spectrometers, - // to_sidis_tup); - - //auto xs_spectrometers = make_jacobian_transformed_xs(diff_xs_sidis, - // init_state, - // jm_spec_dxyzPTphiphie, - // ps_spectrometers); - - //// ---------------------------- - // - auto sampler = make_ps_sampler(make_integrated_cross_section(init_state, diff_xs_sidis)); - //auto sampler = make_ps_sampler(make_integrated_cross_section(init_state, xs_spectrometers)); - - auto fs_parts = make_final_state_particles({11, 211}, {0.000511, M_pion / GeV}, - Eprime_and_Ph_func); - - auto fs_SIDIS = make_final_state(ps_dxyzPTphiphie, fs_parts); - //auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - - auto evgen = make_event_generator(sampler, fs_SIDIS); - auto total_XS = evgen.Init(); - - //// -------------------------- - // - std::cout << "Total cross section = " << total_XS/1000.0 << " ub\n"; - int N_sim = 100000; - - // -------------------------- - // - TFile* f = new TFile("eg_output_sidis_pion.root","RECREATE"); - RDataFrame d1(N_sim); - - auto gen_func = [&](){ - return evgen(); - }; - using gen_func_t = insane::helpers::result_t<decltype(gen_func)>; - - auto gen_func2 = [&](gen_func_t data){ - auto [vars, parts, pdgs] = data; - auto res = to_vector(parts); - //std::transform(res.begin(), res.end(), res.begin(),[&](auto v){return boost_to_lab(v);}); - return res; - }; - using gen_func2_t = typename insane::helpers::result_t<decltype(gen_func2)>; - - auto gen_func_nucleon_rest = [&](gen_func_t data){ - auto [vars, parts, pdgs] = data; - auto res = to_vector(parts); - return res; - }; - using gen_func_nucleon_rest_t = typename insane::helpers::result_t<decltype(gen_func_nucleon_rest)>; - - auto gen_func3 = [&](gen_func_t data){ - auto [vars, parts, pdgs] = data; - return to_vector(vars); - }; - using gen_func3_t = typename insane::helpers::result_t<decltype(gen_func3)>; - - // Generate the event returning the tuple containing - // [thrown variables, final_state, pdg codes] - auto d2 = d1.Define("EvGen", gen_func , {}) - .Define("Pgen", gen_func2 , {"EvGen"}) - .Define("Prest_frame", gen_func_nucleon_rest , {"EvGen"}); - - //// Create leaf for thrown variables - //auto d4 = d3.Define("vars", gen_func3 , {"EvGen"}) - //.Define("Q2", [](const gen_func3_t& vars){ return vars[0];}, {"vars"}) - //.Define("phi", [](const gen_func3_t& vars){ return vars[1];}, {"vars"}) - //.Define("Q20", - // [&](const gen_func2_t& parts ){ - // auto q = init_state_EIC.p1()-parts.at(0); - // double Q2_ev = -1.0*(q.Dot(q)); - // return Q2_ev ; - // }, - // {"Pgen"}) - //.Define("Delta", - // [&](const gen_func2_t& parts ){ - // auto Delta = init_state_EIC.p2()-parts.at(1); - // return Delta ; - // }, {"Pgen"}) - //.Define("u", - // [&](const gen_func2_t& parts ){ - // auto delta = init_state_EIC.p1()-parts.at(1); - // return delta.Dot(delta) ; - // }, {"Pgen"}) - //.Define("Delta_Long", - // [&](const gen_func2_t& parts ){ - // auto Delta = init_state_EIC.p2()-parts.at(1); - // return Delta.pz()/(init_state_EIC.p2().pz()); - // }, {"Pgen"}) - //.Define("Delta_p_over_p", - // [&](const gen_func2_t& parts ){ - // auto Delta = init_state_EIC.p2()-parts.at(1); - // return Delta.P()/(init_state_EIC.p2().P()); - // }, {"Pgen"}) - //.Define("t", - // [&](const gen_func2_t& parts ){ - // auto Delta = init_state_EIC.p2()-parts.at(1); - // double t_ev = (Delta.Dot(Delta)); - // return t_ev ; - // }, {"Pgen"}) - //.Define("theta_e", - // [&](const gen_func2_t& parts ){return parts.at(0).Theta()/degree; }, - // {"Pgen"}) - //.Define("P_e", - // [&](const gen_func2_t& parts ){return parts.at(0).P(); }, - // {"Pgen"}) - //.Define("theta_p", - // [&](const gen_func2_t& parts ){return parts.at(1).Theta()/degree; }, - // {"Pgen"}) - //.Define("P_p", - // [&](const gen_func2_t& parts ){return parts.at(1).P(); }, - // {"Pgen"}) - //.Define("tau", [&](const double& Qsq ){return Qsq/(4.0*0.938*0.938); }, {"Q2"}) - //.Define("theta_e_Rest", - // [&](const gen_func_nucleon_rest_t& P_rest_frame ){ - // double theta_e = pi - P_rest_frame.at(0).Theta(); - // return theta_e/degree; }, - // {"Prest_frame"}) - //.Define("epsilon", - // [&](const double& tau, const gen_func_nucleon_rest_t& P_rest_frame ){ - // double theta_e = pi - P_rest_frame.at(0).Theta(); - // double den = 1.0+2.0*(1.0+tau)*std::pow(std::tan(theta_e/2.0),2) ; - // return 1.0/den; - // }, {"tau","Prest_frame"}) - //.Define("MinusEpsilon", - // [&](const double& tau, const gen_func_nucleon_rest_t& P_rest_frame ){ - // double theta_e = pi - P_rest_frame.at(0).Theta(); - // double den = 1.0+2.0*(1.0+tau)*std::pow(std::tan(theta_e/2.0),2) ; - // return 1.0-1.0/den; - // }, {"tau","Prest_frame"}) - //.Define("Aperp", - // [&](const double& tau, const double& theta_e, const double& eps){ - // double GEoverGM = (1.0/2.79); - // double num = -2.0*std::sqrt(tau*(1.0+tau))*std::tan(theta_e*degree/2.0)*GEoverGM; - // double den = GEoverGM*GEoverGM + tau/eps; - // return num/den; - // }, - // {"tau","theta_e_Rest","epsilon"}) - //.Define("Asym", - // [&](const gen_func2_t& parts, const double& th_e, const double& Q2){ - // static insane::physics::AMTFormFactors ffs; - // XYZVector Ptarg{1,0,0}; - // auto n_0 = init_state_REST.p1().Vect().Cross(parts.at(0).Vect()).Unit(); - // auto q = init_state_REST.p1()-parts.at(0); - // auto q3 = q.Vect(); - // auto z_q = q.Vect().Unit(); - // auto u_x = n_0.Cross(z_q).Unit(); - // auto Pxy = Ptarg - (Ptarg.Dot(z_q))*z_q; - // auto th_star = VectorUtil::Angle(Ptarg, z_q); - // auto phi_star = VectorUtil::Angle(Pxy, u_x); - // return insane::physics::A_elastic(ffs, Q2, th_e*degree, th_star, phi_star); - // }, {"Prest_frame","theta_e_Rest","Q2"}) - //.Define("r_p", - // [&](const double& th_p ){return z_detector*std::sin(th_p*degree);}, - // {"theta_p"}); - - //auto d5 = d4.Filter( - // [](double A){ - // if ( std::abs(A) >0.1 ){ - // return true; - // } - // return false; - // },{"Asym"}); - - //auto d6 = d4.Filter( - // [](double pr){ - // if ( std::abs(pr - 0.01) < 0.005 ){ - // return true; - // } - // return false; - // }, {"r_p"}); - - //// -------------------------- - //d4.Snapshot("elastic_ep", "eg_output_elastic_ep.root"); - d2.Snapshot("pion_sidis", "eg_output_sidis_pion.root"); - - - ////auto tot_count_wCut = d6.Count(); - //auto tot_count_wCut = d5.Count(); - //double total_detected = *tot_count_wCut; - //double day_integrated_lum = 100.0/116.0; // 1/fb per day - //double count_per_day = total_XS*day_integrated_lum; - //double count_per_second = count_per_day/(24.0*60.0*60.0); - //double simulated_time = double(N_sim)/count_per_second; - //double meas_rate_in_bin = total_detected/simulated_time; - - //std::cout << " count_per_second = " << count_per_second << " Hz\n"; - //std::cout << " total_detected = " << total_detected << "\n"; - //std::cout << " stat_unc for 1s = " << std::sqrt(count_per_second)/count_per_second << " \%\n"; - //std::cout << " simulated_time = " << simulated_time << " s\n"; - //std::cout << " meas_rate_in_bin = " << meas_rate_in_bin << " Hz\n"; - //std::cout << " meas_unc for 1s = " << std::sqrt(meas_rate_in_bin)/meas_rate_in_bin << " \%\n"; - //std::cout << " meas_unc for 1m = " << std::sqrt(60.0*meas_rate_in_bin)/(60.0*meas_rate_in_bin) << " \%\n"; - //std::cout << " z_detector = " << z_detector << " m\n"; - - f->Close(); -} - diff --git a/src/runplan/src/sidis_pion_eg2.cxx b/src/runplan/src/sidis_pion_eg2.cxx deleted file mode 100644 index 6ab3425fb0d7f7942a02d96959dc4b634c0d4ea4..0000000000000000000000000000000000000000 --- a/src/runplan/src/sidis_pion_eg2.cxx +++ /dev/null @@ -1,526 +0,0 @@ -#include <iostream> -#include <fstream> -#include <string> -#include <cmath> - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" - -#include "InSANE/MaterialProperties.h" -#include "InSANE/Luminosity.h" -#include "PhaseSpaceVariables.h" -#include "NFoldDifferential.h" -#include "Jacobians.h" -#include "DiffCrossSection.h" -#include "FinalState.h" -#include <typeinfo> - -#include "PSSampler.h" - -//#include "InSANE/DoubleSpinAsymmetries.h" -//#include "InSANE/AMTFormFactors.h" -//#include "InSANE/KellyFormFactors.h" -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" - -#include "InSANE/Stat2015_UPDFs.h" -#include "DSSFragmentationFunctions.h" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" - -#include "ROOT/RDataFrame.hxx" - -#include "csv_settings.h" - -using namespace ROOT; // TDataFrame's namespace - -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; - -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - - -/** Hall C Spectrometer Settings. - * - * SHMS: - * Central Momentum: Momentum Acceptance: Momentum Resolution: - * 2 to 11 GeV/c -10% < δ < +22% 0.03%-0.08% - * Scattering Angle: Horizontal acceptance: Vertical acceptance: - * 5.5 to 40° +- 24mrad +- 40 mrad - * Solid angle acceptance: - * 4.0 msr - * Vertex Length (90deg): Vertex length res (Ytar): - * 30cm 0.1-0.3 cm - * Horizontal resolution (YPtar): Vertical resolution (XPtar): - * 0.5-1.2 mrad 0.3-1.1 mrad - * - * HMS: - * https://www.jlab.org/Hall-C/equipment/HMS/performance.html - * - * - * Coordinate system used: - * z = beam direction - * x = SHMS (beam left) - * y = up (towards hall ceiling) - * - */ -struct HallCSetting { - const double SHMS_dtheta = 0.024; // 24 mrad - const double SHMS_dphi = 0.040; // 40 mrad - const double SHMS_dP_low = 0.1; // -10% - const double SHMS_dP_high = 0.22; // 22% - const double SHMS_solid_angle = 4.0*SHMS_dtheta*SHMS_dphi; // 4 msr - const double SHMS_phi = 0.0; // SHMS sits on the +x side - const double HMS_dtheta = 0.04; // - const double HMS_dphi = 0.04; // - const double HMS_dP_low = 0.09; // -9% - const double HMS_dP_high = 0.09; // 9% - const double HMS_solid_angle = 4.0*HMS_dtheta*HMS_dphi; // 4 msr - const double HMS_phi = insane::units::pi-0.04; - - // HMS detects electron - double HMS_theta = 14.5*degree; - double HMS_p0 = 5.0; // GeV/c - // SHMS detects hadron - double SHMS_theta = 13.5*degree; - double SHMS_p0 = 3.0; // GeV/c - - double HMS_P_min() const { return HMS_p0 * (1.0 - HMS_dP_low); } - double HMS_P_max() const { return HMS_p0 * (1.0 + HMS_dP_high); } - double SHMS_P_min() const { return SHMS_p0 * (1.0 - SHMS_dP_low); } - double SHMS_P_max() const { return SHMS_p0 * (1.0 + SHMS_dP_high); } - - double HMS_phi_min() const { return HMS_phi - HMS_dphi; } - double HMS_phi_max() const { return HMS_phi + HMS_dphi; } - double SHMS_phi_min() const { return SHMS_phi - SHMS_dphi; } - double SHMS_phi_max() const { return SHMS_phi + SHMS_dphi; } - - double HMS_theta_min() const { return HMS_theta - HMS_dtheta; } - double HMS_theta_max() const { return HMS_theta + HMS_dtheta; } - double SHMS_theta_min() const { return SHMS_theta - SHMS_dtheta; } - double SHMS_theta_max() const { return SHMS_theta + SHMS_dtheta; } -}; - -std::vector<double> sidis_pion_eg2(csv::CSV_kinematic kine_set){ - - HallCSetting hcSet; - double x_set = kine_set[0]; - double z_set = kine_set[1]; - double nu_set = kine_set[2]; - double W_set = kine_set[3]; - double Wp_set = kine_set[4]; - hcSet.HMS_theta = kine_set[5]; - hcSet.SHMS_theta = kine_set[6]; - hcSet.HMS_p0 = kine_set[7]; - hcSet.SHMS_p0 = kine_set[8]; - - //SHMS 4 msr dp -10 +20 - //HMS 6 msr dp +-8 - - // ------------------------------------------------------- - // Phase space variables - IPSV phi_SHMS_psv({hcSet.SHMS_phi_min(), hcSet.SHMS_phi_max()}, "phi_shms"); - IPSV phi_HMS_psv({hcSet.HMS_phi_min(), hcSet.HMS_phi_max()}, "phi_hms"); - IPSV theta_HMS_psv({hcSet.HMS_theta_min(), hcSet.HMS_theta_max()}, "theta_hms"); - IPSV theta_SHMS_psv({hcSet.SHMS_theta_min(), hcSet.SHMS_theta_max()}, "theta_shms"); - IPSV P_HMS_psv({hcSet.HMS_P_min(), hcSet.HMS_P_max()}, "P_hms"); - IPSV P_SHMS_psv({hcSet.SHMS_P_min(), hcSet.SHMS_P_max()}, "P_shms"); - - auto diff_spectrometers = make_diff(P_HMS_psv, theta_HMS_psv, phi_HMS_psv, - P_SHMS_psv, theta_SHMS_psv, phi_SHMS_psv); - auto ps_spectrometers = make_phase_space(diff_spectrometers); - diff_spectrometers.Print(); - - // ------------------------------ - // Prepare the initial state - double P1 = 11.0; - double m1 = 0.000511; - double E1 = std::sqrt(P1*P1+m1*m1); - - double P2 = 0.0; - double m2 = M_p/GeV; - double E2 = std::sqrt(P2*P2+m2*m2); - - InitialState init_state(P1, P2, m2); - - /** SIDIS diferential cross section. - * - * d(sigma) - * -------------------------- - * d(x)d(y)d(z)d(phi)d(phi_e) - * - * Note: "phi" is the angle between the leptonic and hadronic planes. - * See http://inspirehep.net/record/677636 for details. - * "phi_e" is the scattered electron's azimuthal angle. - * - */ - auto SIDIS_xs_func = [&](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - static PDFs pdf; - static FragFs ffs; - static double Mh = M_pion/GeV; - - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - if(!(kine.isValid)) return 0.0; - //return 1.0; - - // Gaussian k_perp and p_perp distributions for PT dependence - double kperp2_mean = 0.28; // GeV^2 - double pperp2_mean = 0.25; // GeV^2 - double PT_mean = kine.z * kine.z * kperp2_mean + pperp2_mean; - double PT2 = kine.P_hPerp * kine.P_hPerp; - auto f_PT = std::exp(-PT2 / PT_mean) / (PT_mean * pi); - auto t1 = kine.SIDIS_xs_term() * f_PT; - double FF = 0.0; - auto pdf_vals = pdf.Calculate(kine.x, kine.Q2); - - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF *= kine.x; - - /// \f$ F_{UU,T} = x \sum e_a^2 f_1^a(x) D_1^a(z) \f$ - //double res = Elastic_XS({GE2, GM2}, kine); - //res = res * hbarc2_gev_nb * kine.sin_th; - return(t1*FF); - }; - - // ------------------------------------------------------- - // - auto Eprime_and_Ph_func = [](const InitialState& is, const std::array<double, 6>& x) - { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - //SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - - XYZVector l2(Polar3D<double>(kine.P_l2, kine.theta_l2, kine.phi_l2)); - XYZVector p_h(Polar3D<double>(kine.P_p2, kine.theta_p2, kine.phi_p2)); - //XYZVector p_h_lab(r1*r2*p_h); - - auto E_prime = ROOT::Math::XYZTVector(l2.x(), l2.y(), l2.z(), kine.E_l2); - auto P_hadron = ROOT::Math::XYZTVector(p_h.x(), p_h.y(), p_h.z(), kine.E_p2); - return std::array<ROOT::Math::XYZTVector,2>({E_prime,P_hadron}); - }; - - // returns SIDIS_vars using the SPEC_vars (measured by spectrometers) - auto to_sidis_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - //kine.Print(); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto to_spectrometer_variables = [=](const InitialState& is, const std::array<double,6>& x){ - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion/GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - std::array<double,6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - //auto to_sidis_tup = std::make_tuple( - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<0>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<1>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<2>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<3>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<4>(to_sidis_variables(is,x));}, - // [&](const InitialState& is, const std::array<double,6>& x){ return std::get<5>(to_sidis_variables(is,x));} - // ); - - - //IPSV phi_SHMS_psv({hc_setting.SHMS_phi_min(), hc_setting.SHMS_phi_max()}, "phi_shms"); - //IPSV theta_HMS_psv( {hc_setting.HMS_theta_min(), hc_setting.HMS_theta_max()}, "theta_hms"); - //IPSV theta_SHMS_psv({hc_setting.SHMS_theta_min(), hc_setting.SHMS_theta_max()}, "theta_shms"); - //IPSV P_HMS_psv( {hc_setting.HMS_P_min(), hc_setting.HMS_P_max()}, "P_hms"); - //IPSV P_SHMS_psv({hc_setting.SHMS_P_min(), hc_setting.SHMS_P_max()}, "P_shms"); - // --------------------------------------------------- - // cross section differential in x y z PT phi_h phi_q - //auto ps_spectrometer = make_phase_space( diff_spectrometer );//, - //make_PSV("P_hms", {-pi, pi},//hc_setting.HMS_P_min(), hc_setting.HMS_P_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<0>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("theta_hms", {hc_setting.HMS_theta_min(), hc_setting.HMS_theta_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<1>(to_spectrometer_variables(is, x)); - // })); - //make_PSV("phi_hms", {hc_setting.HMS_phi_min(), hc_setting.HMS_phi_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<2>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("P_shms", {hc_setting.SHMS_P_min(), hc_setting.SHMS_P_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<3>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("theta_shms", {hc_setting.SHMS_theta_min(), hc_setting.SHMS_theta_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<4>(to_spectrometer_variables(is, x)); - // }), - //make_PSV("phi_shms", {hc_setting.SHMS_phi_min(), hc_setting.SHMS_phi_max()}, - // [&](const InitialState& is, const std::array<double, 6>& x) { - // return std::get<5>(to_spectrometer_variables(is, x)); - // })); - //ps_dxyzPTphiphie.Print(); - // - auto diff_xs_sidis = make_diff_cross_section(ps_spectrometers, SIDIS_xs_func); - - // Transform to a function of the 2 measured momenta - //auto jm_spec_dxyzPTphiphie = make_jacobian(diff_dxyzPTphiphie, - // diff_spectrometers, - // to_sidis_tup); - - //auto xs_spectrometers = make_jacobian_transformed_xs(diff_xs_sidis, - // init_state, - // jm_spec_dxyzPTphiphie, - // ps_spectrometers); - - //// ---------------------------- - // - auto sampler = make_ps_sampler(make_integrated_cross_section(init_state, diff_xs_sidis)); - - auto fs_parts = make_final_state_particles({11, 211}, {0.000511, M_pion / GeV}, - Eprime_and_Ph_func); - - auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - //auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - - auto evgen = make_event_generator(sampler, fs_SIDIS); - auto total_XS = evgen.Init(); - - //// -------------------------- - // - std::cout << "Total cross section = " << total_XS << " nb\n"; - int N_sim = 100000; - - using namespace insane; - double Loop1_LD2_Lumi = materials::ComputeLuminosityPerN({1, 2}, - materials::LD2_density /*g/cm3*/, - 10.0 /*cm*/, - 50.0e-6); - double Loop1_entrance_window_thickness = 0.104; - double Loop1_exit_window_thickness = 0.133; - double Loop3_entrance_window_thickness = 0.130; - double Loop3_exit_window_thickness = 0.188; - - double Loop1_windows_Lumi = materials::ComputeLuminosityPerN({13, 27}, - materials::Al_density , - Loop1_entrance_window_thickness + - Loop1_exit_window_thickness , - 50.0e-6); - - double Lumi = Loop1_LD2_Lumi + Loop1_windows_Lumi; - - double total_cross_section = total_XS * nanobarn; - double luminosity = Lumi / cm2; - double rate = total_cross_section * luminosity; - std::cout << "Total XS : " << total_XS << " nb\n"; - std::cout << " : " << total_cross_section/microbarn << " ub\n"; - std::cout << " : " << total_cross_section/barn << " b\n"; - std::cout << "Luminosity : " << Lumi << " 1/cm2 s\n"; - std::cout << " [LD2] " << Loop1_LD2_Lumi << " 1/cm2 s\n"; - std::cout << " [windows] " << Loop1_windows_Lumi << " 1/cm2 s\n"; - std::cout << "Rate : " << rate << " 1/s\n"; - - std::vector<double> res = {total_XS, rate}; - - return res; - - // -------------------------- - // - //TFile* f = new TFile("eg_output_sidis_pion.root","RECREATE"); - //RDataFrame d1(N_sim); - - //auto gen_func = [&](){ - // return evgen(); - //}; - //using gen_func_t = insane::helpers::result_t<decltype(gen_func)>; - - //auto gen_func2 = [&](gen_func_t data){ - // auto [vars, parts, pdgs] = data; - // auto res = to_vector(parts); - // //std::transform(res.begin(), res.end(), res.begin(),[&](auto v){return boost_to_lab(v);}); - // return res; - //}; - //using gen_func2_t = typename insane::helpers::result_t<decltype(gen_func2)>; - - //auto gen_func_nucleon_rest = [&](gen_func_t data){ - // auto [vars, parts, pdgs] = data; - // auto res = to_vector(parts); - // return res; - //}; - //using gen_func_nucleon_rest_t = typename insane::helpers::result_t<decltype(gen_func_nucleon_rest)>; - - //auto gen_func3 = [&](gen_func_t data){ - // auto [vars, parts, pdgs] = data; - // return to_vector(vars); - //}; - //using gen_func3_t = typename insane::helpers::result_t<decltype(gen_func3)>; - - //// Generate the event returning the tuple containing - //// [thrown variables, final_state, pdg codes] - //auto d2 = d1.Define("EvGen", gen_func , {}) - // .Define("Pgen", gen_func2 , {"EvGen"}); - // //.Define("Prest_frame", gen_func_nucleon_rest , {"EvGen"}); - - ////// Create leaf for thrown variables - ////auto d4 = d3.Define("vars", gen_func3 , {"EvGen"}) - ////.Define("Q2", [](const gen_func3_t& vars){ return vars[0];}, {"vars"}) - ////.Define("phi", [](const gen_func3_t& vars){ return vars[1];}, {"vars"}) - ////.Define("Q20", - //// [&](const gen_func2_t& parts ){ - //// auto q = init_state_EIC.p1()-parts.at(0); - //// double Q2_ev = -1.0*(q.Dot(q)); - //// return Q2_ev ; - //// }, - //// {"Pgen"}) - ////.Define("Delta", - //// [&](const gen_func2_t& parts ){ - //// auto Delta = init_state_EIC.p2()-parts.at(1); - //// return Delta ; - //// }, {"Pgen"}) - ////.Define("u", - //// [&](const gen_func2_t& parts ){ - //// auto delta = init_state_EIC.p1()-parts.at(1); - //// return delta.Dot(delta) ; - //// }, {"Pgen"}) - ////.Define("Delta_Long", - //// [&](const gen_func2_t& parts ){ - //// auto Delta = init_state_EIC.p2()-parts.at(1); - //// return Delta.pz()/(init_state_EIC.p2().pz()); - //// }, {"Pgen"}) - ////.Define("Delta_p_over_p", - //// [&](const gen_func2_t& parts ){ - //// auto Delta = init_state_EIC.p2()-parts.at(1); - //// return Delta.P()/(init_state_EIC.p2().P()); - //// }, {"Pgen"}) - ////.Define("t", - //// [&](const gen_func2_t& parts ){ - //// auto Delta = init_state_EIC.p2()-parts.at(1); - //// double t_ev = (Delta.Dot(Delta)); - //// return t_ev ; - //// }, {"Pgen"}) - ////.Define("theta_e", - //// [&](const gen_func2_t& parts ){return parts.at(0).Theta()/degree; }, - //// {"Pgen"}) - ////.Define("P_e", - //// [&](const gen_func2_t& parts ){return parts.at(0).P(); }, - //// {"Pgen"}) - ////.Define("theta_p", - //// [&](const gen_func2_t& parts ){return parts.at(1).Theta()/degree; }, - //// {"Pgen"}) - ////.Define("P_p", - //// [&](const gen_func2_t& parts ){return parts.at(1).P(); }, - //// {"Pgen"}) - ////.Define("tau", [&](const double& Qsq ){return Qsq/(4.0*0.938*0.938); }, {"Q2"}) - ////.Define("theta_e_Rest", - //// [&](const gen_func_nucleon_rest_t& P_rest_frame ){ - //// double theta_e = pi - P_rest_frame.at(0).Theta(); - //// return theta_e/degree; }, - //// {"Prest_frame"}) - ////.Define("epsilon", - //// [&](const double& tau, const gen_func_nucleon_rest_t& P_rest_frame ){ - //// double theta_e = pi - P_rest_frame.at(0).Theta(); - //// double den = 1.0+2.0*(1.0+tau)*std::pow(std::tan(theta_e/2.0),2) ; - //// return 1.0/den; - //// }, {"tau","Prest_frame"}) - ////.Define("MinusEpsilon", - //// [&](const double& tau, const gen_func_nucleon_rest_t& P_rest_frame ){ - //// double theta_e = pi - P_rest_frame.at(0).Theta(); - //// double den = 1.0+2.0*(1.0+tau)*std::pow(std::tan(theta_e/2.0),2) ; - //// return 1.0-1.0/den; - //// }, {"tau","Prest_frame"}) - ////.Define("Aperp", - //// [&](const double& tau, const double& theta_e, const double& eps){ - //// double GEoverGM = (1.0/2.79); - //// double num = -2.0*std::sqrt(tau*(1.0+tau))*std::tan(theta_e*degree/2.0)*GEoverGM; - //// double den = GEoverGM*GEoverGM + tau/eps; - //// return num/den; - //// }, - //// {"tau","theta_e_Rest","epsilon"}) - ////.Define("Asym", - //// [&](const gen_func2_t& parts, const double& th_e, const double& Q2){ - //// static insane::physics::AMTFormFactors ffs; - //// XYZVector Ptarg{1,0,0}; - //// auto n_0 = init_state_REST.p1().Vect().Cross(parts.at(0).Vect()).Unit(); - //// auto q = init_state_REST.p1()-parts.at(0); - //// auto q3 = q.Vect(); - //// auto z_q = q.Vect().Unit(); - //// auto u_x = n_0.Cross(z_q).Unit(); - //// auto Pxy = Ptarg - (Ptarg.Dot(z_q))*z_q; - //// auto th_star = VectorUtil::Angle(Ptarg, z_q); - //// auto phi_star = VectorUtil::Angle(Pxy, u_x); - //// return insane::physics::A_elastic(ffs, Q2, th_e*degree, th_star, phi_star); - //// }, {"Prest_frame","theta_e_Rest","Q2"}) - ////.Define("r_p", - //// [&](const double& th_p ){return z_detector*std::sin(th_p*degree);}, - //// {"theta_p"}); - - ////auto d5 = d4.Filter( - //// [](double A){ - //// if ( std::abs(A) >0.1 ){ - //// return true; - //// } - //// return false; - //// },{"Asym"}); - - ////auto d6 = d4.Filter( - //// [](double pr){ - //// if ( std::abs(pr - 0.01) < 0.005 ){ - //// return true; - //// } - //// return false; - //// }, {"r_p"}); - - ////// -------------------------- - ////d4.Snapshot("elastic_ep", "eg_output_elastic_ep.root"); - //d2.Snapshot("pion_sidis", "eg_output_sidis_pion.root"); - - - //////auto tot_count_wCut = d6.Count(); - ////auto tot_count_wCut = d5.Count(); - ////double total_detected = *tot_count_wCut; - ////double day_integrated_lum = 100.0/116.0; // 1/fb per day - ////double count_per_day = total_XS*day_integrated_lum; - ////double count_per_second = count_per_day/(24.0*60.0*60.0); - ////double simulated_time = double(N_sim)/count_per_second; - ////double meas_rate_in_bin = total_detected/simulated_time; - - ////std::cout << " count_per_second = " << count_per_second << " Hz\n"; - ////std::cout << " total_detected = " << total_detected << "\n"; - ////std::cout << " stat_unc for 1s = " << std::sqrt(count_per_second)/count_per_second << " \%\n"; - ////std::cout << " simulated_time = " << simulated_time << " s\n"; - ////std::cout << " meas_rate_in_bin = " << meas_rate_in_bin << " Hz\n"; - ////std::cout << " meas_unc for 1s = " << std::sqrt(meas_rate_in_bin)/meas_rate_in_bin << " \%\n"; - ////std::cout << " meas_unc for 1m = " << std::sqrt(60.0*meas_rate_in_bin)/(60.0*meas_rate_in_bin) << " \%\n"; - ////std::cout << " z_detector = " << z_detector << " m\n"; - - //f->Close(); -} - diff --git a/src/runplan/src/sidis_rates.cxx b/src/runplan/src/sidis_rates.cxx deleted file mode 100644 index d51928ecb4f31c374ad42428d44558ac87cc8a0e..0000000000000000000000000000000000000000 --- a/src/runplan/src/sidis_rates.cxx +++ /dev/null @@ -1,644 +0,0 @@ -#include "sidis_rates.h" -#include "runplan/CSV_settings.h" -#include <iostream> - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" -#include "ROOT/TFuture.hxx" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" -#include "TBufferJSON.h" -#include "ROOT/RDataFrame.hxx" - -// insane libraries -#include "InSANE/MaterialProperties.h" -#include "InSANE/Luminosity.h" -#include "InSANE/PhaseSpaceVariables.h" -#include "InSANE/PSSampler.h" -#include "InSANE/NFoldDifferential.h" -#include "InSANE/Jacobians.h" -#include "InSANE/DiffCrossSection.h" -#include "InSANE/FinalState.h" -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" -#include "InSANE/Stat2015_UPDFs.h" -#include "InSANE/DSSFragmentationFunctions.h" - -#include <fmt/core.h> -#include <fmt/ostream.h> -R__LOAD_LIBRARY(libfmt.so) - -using namespace ROOT; -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; -using namespace ROOT::Experimental; - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - -namespace hallc { - - /** Calculates the table row entry. - * - * Computes cross section, rate, and time based on the desired statistics. - */ - RunPlanTableEntry build_table_entry(const Kinematic& kine) { - RunPlanTableEntry entry0; - entry0.kinematic = kine; - entry0.polarity = 1; - if(kine.pid_had < 0 ) { - entry0.polarity = -1; - } - //entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_xs_and_rates(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - // if( use_LH2 ) { - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - //} - return entry0; - } - - RunTable build_sidis_table(const csv::CSV_Settings& settings) { - // The fortran code used for the fragmentation functions - // doesn't allow multi threaded to run without modifying - // the way it uses the common blocks - ROOT::EnableImplicitMT(1); - //{ std::ofstream json_ofile("tables/pion_rates.json", std::ios_base::trunc); } - - // Copy the settings - auto csv_settings = settings; - RunTable run_plan_table; - - for (auto& a_kinematic : csv_settings) { - - a_kinematic.Compute(); - - auto plus_entry = ROOT::Experimental::Async([&a_kinematic]() { - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_kinematic; - entry0.polarity = 1; - entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_xs_and_rates(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - return entry0; - }); - - auto minus_entry = ROOT::Experimental::Async([&a_kinematic]() { - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_kinematic; - entry0.polarity = -1; - entry0.Ibeam = 50.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_xs_and_rates(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - // if( ! only_LH2 ) { - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return entry0; - }); - - run_plan_table.add(plus_entry.get()); - run_plan_table.add(minus_entry.get()); - } - // all_tables.push_back(std::make_pair(setting_group.first, run_plan_table)); - // all_kaon_tables.push_back(std::make_pair(setting_group.first, run_plan_kaon_table)); - //} - return run_plan_table; - } - - void save_json_tables(const RunTable& table_rows, std::string table_name, - const std::ios_base::openmode& mode) { - double total_time = 0.0; - std::ofstream ofs("tables/run_plan_table.txt", mode); - // for(const auto& atable : all_tables){ - // ofs << "Q2 = " << atable.first << "\n"; - double Q2 = table_rows.GetVector().at(0).kinematic.Q2; - double Q2group_total_time = 0.0; - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", Q2)); - RunPlanTableEntry::PrintHeader(ofs); - for (const auto& entry : table_rows.GetVector()) { - entry.Print(ofs); - total_time += entry.time; - Q2group_total_time += entry.time; - } - fmt::print(ofs, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", Q2) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0); - ofs << "\n"; - //} - fmt::print(ofs, "{:=<79s}\n", "="); - ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - - std::ofstream json_ofile("tables/run_plan_table.json", std::ios_base::trunc); - //json_ofile << TBufferJSON::ToJSON(&all_tables); - json_ofile.close(); - } - - void save_tables(const RunTable& table_rows, std::string table_name, - const std::ios_base::openmode& mode) { - { - double total_time = 0.0; - std::ofstream ofs(table_name + ".txt", std::ios_base::app); - std::ofstream wiki_file(table_name + ".wiki", std::ios_base::app); - // for (const auto& atable : all_tables) { - // ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - - double Q2 = table_rows.GetVector().at(0).kinematic.Q2; - - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", Q2)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", Q2)); - - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for (const auto& entry : table_rows.GetVector()) { - if (std::abs(entry.polarity) >= 2) { - continue; - } - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - - fmt::print(ofs, "{:-<79s}\n", "-"); - fmt::print(wiki_file, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", Q2) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - wiki_file << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", Q2) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - //} - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " - // days\n"; fmt::print(wiki_file, "{:=<79s}\n", "="); wiki_file << " total time: " << - // total_time << " hours or " << total_time / 24.0 - // << " days\n"; - - // std::ofstream json_ofile("tables/LD2_run_plan_table.json", std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - } - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/LD2_run_plan_kaon_table.txt", std::ios_base::trunc); - // std::ofstream wiki_file("tables/LD2_run_plan_kaon_table.wiki", std::ios_base::trunc); - // for (const auto& atable : all_kaon_tables) { - // // ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - // fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // RunPlanTableEntry::PrintWikiHeader(wiki_file); - // for (const auto& entry : atable.second) { - // entry.Print(ofs); - // entry.PrintWiki(wiki_file); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // RunPlanTableEntry::PrintWikiFooter(wiki_file); - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - - // std::ofstream json_ofile("tables/LD2_run_plan_kaon_table.json", std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/run_plan_kaon_table.txt",std::ios_base::trunc); - // for(const auto& atable : all_kaon_tables){ - // //ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // for(const auto& entry : atable.second){ - // entry.Print(ofs); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - // std::ofstream json_ofile("tables/run_plan_kaon_table.json",std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - } - - /** Compute the cross section and rate for a given setting - * - */ - std::vector<double> sidis_config_xs_and_rates(RunPlanTableEntry& run_setting) { - - hallc::Kinematic& kine_set = run_setting.kinematic; - hallc::HallCSetting& hcSet = run_setting.hcSet; - double x_set = kine_set.x; - double z_set = kine_set.z; - double nu_set = kine_set.nu; - double W_set = kine_set.W; - double Wp_set = kine_set.Wp; - int polarity = run_setting.polarity; - hcSet.HMS_theta = kine_set.th_e * degree; - hcSet.SHMS_theta = kine_set.th_q * degree; - hcSet.HMS_p0 = kine_set.Ee; - hcSet.SHMS_p0 = kine_set.Ppi; - - // SHMS 4 msr dp -10 +20 - // HMS 6 msr dp +-8 - // ------------------------------------------------------- - // Phase space variables - //hcSet.SHMS_dP_high = -0.05; - IPSV lowdp_phi_SHMS_psv({hcSet.SHMS_phi_min(), hcSet.SHMS_phi_max()}, "phi_shms"); - IPSV lowdp_phi_HMS_psv({hcSet.HMS_phi_min(), hcSet.HMS_phi_max()}, "phi_hms"); - IPSV lowdp_theta_HMS_psv({hcSet.HMS_theta_min(), hcSet.HMS_theta_max()}, "theta_hms"); - IPSV lowdp_theta_SHMS_psv({hcSet.SHMS_theta_min(), hcSet.SHMS_theta_max()}, "theta_shms"); - IPSV lowdp_P_HMS_psv({hcSet.HMS_P_min(), hcSet.HMS_P_max()}, "P_hms"); - IPSV lowdp_P_SHMS_psv({hcSet.SHMS_P_min(), hcSet.SHMS_P_max()}, "P_shms"); - - - //hcSet.SHMS_dP_high = hallc::shms::SHMS_dP_high; - //hcSet.SHMS_dP_low = -0.05; - IPSV phi_SHMS_psv({hcSet.SHMS_phi_min(), hcSet.SHMS_phi_max()}, "phi_shms"); - IPSV phi_HMS_psv({hcSet.HMS_phi_min(), hcSet.HMS_phi_max()}, "phi_hms"); - IPSV theta_HMS_psv({hcSet.HMS_theta_min(), hcSet.HMS_theta_max()}, "theta_hms"); - IPSV theta_SHMS_psv({hcSet.SHMS_theta_min(), hcSet.SHMS_theta_max()}, "theta_shms"); - IPSV P_HMS_psv({hcSet.HMS_P_min(), hcSet.HMS_P_max()}, "P_hms"); - IPSV P_SHMS_psv({hcSet.SHMS_P_min(), hcSet.SHMS_P_max()}, "P_shms"); - - hcSet.SHMS_dP_low = hallc::shms::SHMS_dP_low; - - auto diff_spectrometers = - make_diff(P_HMS_psv, theta_HMS_psv, phi_HMS_psv, P_SHMS_psv, theta_SHMS_psv, phi_SHMS_psv); - auto ps_spectrometers = make_phase_space(diff_spectrometers); - - auto diff_spectrometers_lowdp = - make_diff(lowdp_P_HMS_psv, lowdp_theta_HMS_psv, lowdp_phi_HMS_psv, lowdp_P_SHMS_psv, lowdp_theta_SHMS_psv, lowdp_phi_SHMS_psv); - auto ps_spectrometers_lowdp = make_phase_space(diff_spectrometers_lowdp); - //diff_spectrometers.Print(); - //diff_spectrometers_lowdp.Print(); - - // ------------------------------ - // Prepare the initial state - double P1 = kine_set.E0; - double m1 = 0.000511; - double E1 = std::sqrt(P1 * P1 + m1 * m1); - - double P2 = 0.0; - double m2 = M_p / GeV; - double E2 = std::sqrt(P2 * P2 + m2 * m2); - - InitialState init_state(P1, P2, m2); - - int xs_select = 1; - PDFs pdf; - FragFs ffs; - - /** SIDIS diferential cross section. - * - * d(sigma) - * -------------------------- - * d(x)d(y)d(z)d(phi)d(phi_e) - * - * Note: "phi" is the angle between the leptonic and hadronic planes. - * See http://inspirehep.net/record/677636 for details. - * "phi_e" is the scattered electron's azimuthal angle. - * - */ - auto SIDIS_xs_func = [&](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - double Mh = M_pion / GeV; - - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - if (!(kine.isValid)) - return 0.0; - - // Gaussian k_perp and p_perp distributions for PT dependence - double kperp2_mean = 0.28; // GeV^2 - double pperp2_mean = 0.25; // GeV^2 - double PT_mean = kine.z * kine.z * kperp2_mean + pperp2_mean; - double PT2 = kine.P_hPerp * kine.P_hPerp; - auto f_PT = std::exp(-PT2 / PT_mean) / (PT_mean * pi); - auto t1 = kine.SIDIS_xs_term() * f_PT; - double FF = 0.0; - auto pdf_vals = pdf.Calculate(kine.x, kine.Q2); - - if (polarity == 1) { - // pi+ - if (xs_select == 1 || xs_select == 2) { - // proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -1) { - // pi- - if (xs_select == 1 || xs_select == 2) { - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == 2) { - // pi+ - if (xs_select == 1 || xs_select == 2) { - // proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -2) { - // pi- - if (xs_select == 1 || xs_select == 2) { - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } - FF *= kine.x; - - /// \f$ F_{UU,T} = x \sum e_a^2 f_1^a(x) D_1^a(z) \f$ - // double res = Elastic_XS({GE2, GM2}, kine); - // res = res * hbarc2_gev_nb * kine.sin_th; - return (t1 * FF); - }; - - // ------------------------------------------------------- - // - auto Eprime_and_Ph_func = [](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - XYZVector l2(Polar3D<double>(kine.P_l2, kine.theta_l2, kine.phi_l2)); - XYZVector p_h(Polar3D<double>(kine.P_p2, kine.theta_p2, kine.phi_p2)); - auto E_prime = ROOT::Math::XYZTVector(l2.x(), l2.y(), l2.z(), kine.E_l2); - auto P_hadron = ROOT::Math::XYZTVector(p_h.x(), p_h.y(), p_h.z(), kine.E_p2); - return std::array<ROOT::Math::XYZTVector, 2>({E_prime, P_hadron}); - }; - - // returns SIDIS_vars using the SPEC_vars (measured by spectrometers) - auto to_sidis_variables = [=](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - // kine.Print(); - std::array<double, 6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto to_spectrometer_variables = [=](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - std::array<double, 6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto diff_xs_sidis = make_diff_cross_section(ps_spectrometers, SIDIS_xs_func); - auto diff_xs_sidis_lowdp = make_diff_cross_section(ps_spectrometers_lowdp, SIDIS_xs_func); - - auto fs_parts = - make_final_state_particles({11, 211}, {0.000511, M_pion / GeV}, Eprime_and_Ph_func); - - auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - - auto p_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - auto n_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - - //auto p_integrated_XS_lowdp = make_integrated_cross_section(init_state, diff_xs_sidis_lowdp); - //auto n_integrated_XS_lowdp = make_integrated_cross_section(init_state, diff_xs_sidis_lowdp); - // auto sampler = make_ps_sampler(p_integrated_XS); - // auto sampler2 = make_ps_sampler(n_integrated_XS); - // sampler.SetFoamCells(500); - // sampler.SetFoamSample(100); - // sampler.SetFoamChat(0); - // sampler2.SetFoamCells(500); - // sampler2.SetFoamSample(100); - // sampler2.SetFoamChat(0); - - // auto evgen = make_event_generator(sampler, fs_SIDIS); - xs_select = 1; - auto total_XS_p = p_integrated_XS.CalculateTotalXS();// + p_integrated_XS_lowdp.CalculateTotalXS(); - // auto total_XS_p = evgen.Init(); - - // auto evgen2 = make_event_generator(sampler2, fs_SIDIS); - xs_select = 0; - auto total_XS_n = n_integrated_XS.CalculateTotalXS();// + n_integrated_XS_lowdp.CalculateTotalXS(); - // auto total_XS_n = evgen2.Init(); - - double total_XS = total_XS_n + total_XS_p; - - //// -------------------------- - // - // std::cout << "Total cross section = " << total_XS << " nb\n"; - int N_sim = 100000; - - using namespace insane; - double Loop1_LD2_Lumi = materials::ComputeLuminosityPerN( - {1, 2}, materials::LD2_density /*g/cm3*/, 10.0 /*cm*/, 50.0e-6); - double Loop1_LH2_Lumi = materials::ComputeLuminosityPerN( - {1, 1}, materials::LH2_density /*g/cm3*/, 10.0 /*cm*/, 50.0e-6); - double Loop1_entrance_window_thickness = 0.0104; - double Loop1_exit_window_thickness = 0.0133; - double Loop3_entrance_window_thickness = 0.0130; - double Loop3_exit_window_thickness = 0.0188; - - double Loop1_windows_Lumi = materials::ComputeLuminosityPerN( - {13, 27}, materials::Al_density, - Loop1_entrance_window_thickness + Loop1_exit_window_thickness, 50.0e-6); - - // double Lumi = Loop1_LD2_Lumi + Loop1_windows_Lumi; - TargetRates targ_rates; - - targ_rates.LD2_Lumi = Loop1_LD2_Lumi; - targ_rates.LH2_Lumi = Loop1_LH2_Lumi; - targ_rates.window_Lumi = Loop1_windows_Lumi; - - // double total_cross_section = 2*total_XS_p * nanobarn; - // double luminosity = Lumi / cm2; - targ_rates.p_XS = total_XS_p; - targ_rates.n_XS = total_XS_n; - targ_rates.LD2_XS = total_XS_n + total_XS_p; - targ_rates.LH2_XS = total_XS_p; - targ_rates.window_XS = (total_XS_n * 14 + total_XS_p * 13); - - targ_rates.LD2_rate = targ_rates.LD2_XS * targ_rates.LD2_Lumi * (nanobarn / cm2); - targ_rates.LH2_rate = targ_rates.LH2_XS * targ_rates.LH2_Lumi * (nanobarn / cm2); - targ_rates.window_rate = targ_rates.window_XS * targ_rates.window_Lumi * (nanobarn / cm2); - targ_rates.total_LD2_rate = targ_rates.window_rate + targ_rates.LD2_rate; - targ_rates.total_LH2_rate = targ_rates.window_rate + targ_rates.LH2_rate; - targ_rates.total_rate = targ_rates.window_rate + targ_rates.LD2_rate; - - run_setting.rates = targ_rates; - // run_setting.LD2_Luminosity = Loop1_LD2_Lumi; - // run_setting.window_Luminosity = Loop1_windows_Lumi; - // run_setting.LD2_sigma = total_xs_LD2; - // run_setting.window_sigma = total_xs_windows; - std::vector<double> res = {targ_rates.LD2_XS, targ_rates.LD2_rate, - targ_rates.total_LD2_rate, targ_rates.total_LH2_rate, - targ_rates.LH2_XS, targ_rates.LH2_rate}; - return res; - } - -} // namespace hallc diff --git a/src/runplan/src/sidis_settings_rate.cxx b/src/runplan/src/sidis_settings_rate.cxx deleted file mode 100644 index a1b551a675c851534fdeed9c75e0a1573c258bd4..0000000000000000000000000000000000000000 --- a/src/runplan/src/sidis_settings_rate.cxx +++ /dev/null @@ -1,1012 +0,0 @@ -#include "sidis_settings_rate.h" - -#include "Math/Vector3Dfwd.h" -#include "Math/Vector4Dfwd.h" -#include "Math/Transform3D.h" -#include "Math/Rotation3D.h" -#include "Math/RotationY.h" -#include "Math/RotationZ.h" -#include "ROOT/TFuture.hxx" - -#include "TH1F.h" -#include "TFile.h" -#include "TTree.h" -#include "TBufferJSON.h" -#include "ROOT/RDataFrame.hxx" - -// insane libraries -#include "InSANE/MaterialProperties.h" -#include "InSANE/Luminosity.h" -#include "InSANE/PhaseSpaceVariables.h" -#include "InSANE/PSSampler.h" -#include "InSANE/NFoldDifferential.h" -#include "InSANE/Jacobians.h" -#include "InSANE/DiffCrossSection.h" -#include "InSANE/FinalState.h" -#include "InSANE/Helpers.h" -#include "InSANE/Physics.h" -#include "InSANE/Kinematics.h" -#include "InSANE/Stat2015_UPDFs.h" -#include "InSANE/DSSFragmentationFunctions.h" - -#include <fmt/core.h> -#include <fmt/ostream.h> -R__LOAD_LIBRARY(libfmt.so) - -using namespace ROOT; -using namespace insane::physics; -using namespace insane::kinematics; -using namespace insane::units; -using namespace insane::helpers; -using namespace ROOT::Math; -using namespace ROOT::Experimental; - -using PDFs = insane::physics::Stat2015_UPDFs; -using FragFs = insane::physics::DSSFragmentationFunctions; -using SIDIS_vars = insane::kinematics::variables::SIDIS_x_y_z_phih_phie; -using SPEC_vars = insane::kinematics::variables::SIDIS_eprime_ph; -using Q2_variable = insane::kinematics::variables::MomentumTransfer; -using Theta_variable = insane::kinematics::variables::Theta; - -namespace hallc { - - /** Calculates the table row entry. - * Computes cross section, rate, and time based on the desired statistics. - */ - RunPlanTableEntry build_table_entry(const Kinematic& kine) { - RunPlanTableEntry entry0; - entry0.kinematic = kine; - entry0.polarity = 1; - if(kine.pid_had < 0 ) { - entry0.polarity = -1; - } - //entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - //if( ! only_LH2 ) { - // resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return entry0; - } - - std::vector<RunPlanTableEntry> build_sidis_table(const csv::CSV_Settings& settings, - double Q2_new) { - - // The fortran code used for the fragmentation functions - // doesn't allow multi threaded to run without modifying - // the way it uses the common blocks - ROOT::EnableImplicitMT(1); - - //{ std::ofstream json_ofile("tables/pion_rates.json", std::ios_base::trunc); } - - // Copy the settings - auto csv_settings = settings; - //for (auto& setting_group : csv_all_settings) { - - // We only are interested in the LH2 target rates for the two lowest Q2 settings - bool use_LH2 = false; - // if( std::get<0>(setting_group) > 5.5 ) { use_LH2 = false; } - bool only_LH2 = false; - // if( !use_LH2) { - // continue; - //} - - std::vector<RunPlanTableEntry> run_plan_table; - - for (auto& a_setting : csv_settings) { - - a_setting = RecomputeKinematic_SIDIS(a_setting, Q2_new); - //std::cout << a_setting.Q2 << " \n"; - // a_setting.Q2 = std::get<0>(setting_group); - - auto plus_entry = ROOT::Experimental::Async([&a_setting, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = 1; - entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto minus_entry = ROOT::Experimental::Async([&a_setting, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = -1; - entry0.Ibeam = 50.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto table_entry0 = plus_entry.get(); - auto table_entry1 = minus_entry.get(); - - run_plan_table.insert(run_plan_table.end(), table_entry0.begin(), table_entry0.end()); - run_plan_table.insert(run_plan_table.end(), table_entry1.begin(), table_entry1.end()); - - } - //all_tables.push_back(std::make_pair(setting_group.first, run_plan_table)); - //all_kaon_tables.push_back(std::make_pair(setting_group.first, run_plan_kaon_table)); - //} - return run_plan_table; - } - - /** Calculates a table of rates. - * - */ - void sidis_settings_rate() { - - // The fortran code used for the fragmentation functions - // doesn't allow multi threaded to run without modifying - // the way it uses the common blocks - ROOT::EnableImplicitMT(1); - - std::vector<std::pair<double, std::vector<RunPlanTableEntry>>> all_tables; - std::vector<std::pair<double, std::vector<RunPlanTableEntry>>> all_kaon_tables; - - //{ std::ofstream json_ofile("tables/pion_rates.json", std::ios_base::trunc); } - - // Copy the settings - auto csv_all_settings = csv::all_settings; - - for (auto& setting_group : csv_all_settings) { - - // We only are interested in the LH2 target rates for the two lowest Q2 settings - bool use_LH2 = false; - // if( std::get<0>(setting_group) > 5.5 ) { use_LH2 = false; } - bool only_LH2 = false; - // if( !use_LH2) { - // continue; - //} - - std::vector<RunPlanTableEntry> run_plan_kaon_table; - std::vector<RunPlanTableEntry> run_plan_table; - for (auto& a_setting : setting_group.second) { - - a_setting = RecomputeKinematic_SIDIS(a_setting, std::get<0>(setting_group)); - std::cout << a_setting.Q2 << " \n"; - // a_setting.Q2 = std::get<0>(setting_group); - - auto plus_entry = ROOT::Experimental::Async([&a_setting, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = 1; - entry0.Ibeam = 25.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto minus_entry = ROOT::Experimental::Async([&a_setting, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = -1; - entry0.Ibeam = 50.0; - entry0.counts = 30000; - if (entry0.kinematic.z > 0.55) { - entry0.counts = 20000; - } - sidis_config_sigma_rate(entry0); - double scale = 50.0 / entry0.Ibeam; - entry0.time = scale * (entry0.counts / entry0.rates.total_rate) / (60.0 * 60.0); - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = entry1.counts/entry1.rate/(60.0*60.0); - // entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto table_entry0 = plus_entry.get(); - auto table_entry1 = minus_entry.get(); - - { - std::ofstream json_ofile("tables/pion_rates.json", std::ios_base::app); - json_ofile << TBufferJSON::ToJSON(&table_entry0); - json_ofile << TBufferJSON::ToJSON(&table_entry1); - } - - auto Kplus_entry = - ROOT::Experimental::Async([&a_setting, &table_entry0, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi plus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = 2; - entry0.Ibeam = 25.0; - double scale = 50.0 / entry0.Ibeam; - sidis_config_sigma_rate(entry0); - entry0.time = table_entry0[0].time; //.counts/entry0.rate/(60.0*60.0); - entry0.counts = (60.0 * 60.0) * entry0.time * entry0.rates.total_rate; - entry0.LD2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LD2_rate; - entry0.LH2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LH2_rate; - entry0.window_counts = (60.0 * 60.0) * entry0.time * entry0.rates.window_rate; - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = - // table_entry0[0].time;//entry1.counts/entry1.rate/(60.0*60.0); entry1.counts = - // (60.0 * 60.0) * entry1.time * entry1.rate; entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - auto Kminus_entry = - ROOT::Experimental::Async([&a_setting, &table_entry1, &use_LH2, &only_LH2]() { - std::vector<RunPlanTableEntry> resulting_entries; - // pi minus - RunPlanTableEntry entry0; - entry0.kinematic = a_setting; - entry0.polarity = -2; - entry0.Ibeam = 50.0; - double scale = 50.0 / entry0.Ibeam; - sidis_config_sigma_rate(entry0); - entry0.time = table_entry1[0].time; //.counts/entry0.rate/(60.0*60.0); - entry0.counts = (60.0 * 60.0) * entry0.time * entry0.rates.total_rate; - entry0.LD2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LD2_rate; - entry0.LH2_counts = (60.0 * 60.0) * entry0.time * entry0.rates.LH2_rate; - entry0.window_counts = (60.0 * 60.0) * entry0.time * entry0.rates.window_rate; - RunPlanTableEntry::PrintHeader(); - entry0.Print(); - - // if( ! only_LH2 ) { - resulting_entries.push_back(entry0); - //} - // if( use_LH2 ) { - // RunPlanTableEntry entry1 = entry0; - // entry1.use_LH2_target = true; - // entry1.rate = res_rate[3]*(entry1.Ibeam/50.0); - // entry1.sigma = res_rate[4]; - // entry1.time = - // table_entry1[0].time;//entry1.counts/entry1.rate/(60.0*60.0); entry1.counts = - // (60.0 * 60.0) * entry1.time * entry1.rate; entry1.Print(); - // resulting_entries.push_back( entry1 ); - //} - return resulting_entries; - }); - - // run_plan_table.push_back( entry_plus ); - // run_plan_table.push_back( entry_minus ); - auto table_entry2 = Kplus_entry.get(); - auto table_entry3 = Kminus_entry.get(); - - run_plan_table.insert(run_plan_table.end(), table_entry0.begin(), table_entry0.end()); - run_plan_table.insert(run_plan_table.end(), table_entry1.begin(), table_entry1.end()); - - run_plan_kaon_table.insert(run_plan_kaon_table.end(), table_entry2.begin(), - table_entry2.end()); - run_plan_kaon_table.insert(run_plan_kaon_table.end(), table_entry3.begin(), - table_entry3.end()); - - std::cout << "done\n"; - // run_plan_table.push_back( table_entry1[0] ); - // run_plan_kaon_table.push_back( Kplus_entry.get()[0] ); - // run_plan_kaon_table.push_back( Kminus_entry.get()[0] ); - - // run_plan_table.push_back( Kplus_entry.get() ); - // run_plan_table.push_back( Kminus_entry.get() ); - } - all_tables.push_back(std::make_pair(setting_group.first, run_plan_table)); - all_kaon_tables.push_back(std::make_pair(setting_group.first, run_plan_kaon_table)); - } - - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/run_plan_table.txt",std::ios_base::trunc); - // for(const auto& atable : all_tables){ - // //ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // for(const auto& entry : atable.second){ - // entry.Print(ofs); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - // std::ofstream json_ofile("tables/run_plan_table.json",std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - { - double total_time = 0.0; - std::ofstream ofs("tables/LD2_run_plan_table.txt", std::ios_base::trunc); - std::ofstream wiki_file("tables/LD2_run_plan_table.wiki", std::ios_base::trunc); - for (const auto& atable : all_tables) { - // ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for (const auto& entry : atable.second) { - if (std::abs(entry.polarity) >= 2) { - continue; - } - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - - fmt::print(ofs, "{:-<79s}\n", "-"); - fmt::print(wiki_file, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - wiki_file << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - } - fmt::print(ofs, "{:=<79s}\n", "="); - ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - fmt::print(wiki_file, "{:=<79s}\n", "="); - wiki_file << " total time: " << total_time << " hours or " << total_time / 24.0 - << " days\n"; - - std::ofstream json_ofile("tables/LD2_run_plan_table.json", std::ios_base::trunc); - json_ofile << TBufferJSON::ToJSON(&all_tables); - json_ofile.close(); - } - { - double total_time = 0.0; - std::ofstream ofs("tables/LD2_run_plan_kaon_table.txt", std::ios_base::trunc); - std::ofstream wiki_file("tables/LD2_run_plan_kaon_table.wiki", std::ios_base::trunc); - for (const auto& atable : all_kaon_tables) { - // ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for (const auto& entry : atable.second) { - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - fmt::print(ofs, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0); - ofs << "\n"; - } - fmt::print(ofs, "{:=<79s}\n", "="); - ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - - std::ofstream json_ofile("tables/LD2_run_plan_kaon_table.json", std::ios_base::trunc); - json_ofile << TBufferJSON::ToJSON(&all_tables); - json_ofile.close(); - } - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/run_plan_kaon_table.txt",std::ios_base::trunc); - // for(const auto& atable : all_kaon_tables){ - // //ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // for(const auto& entry : atable.second){ - // entry.Print(ofs); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - // std::ofstream json_ofile("tables/run_plan_kaon_table.json",std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - } - - void save_tables(const std::vector<RunPlanTableEntry>& table_rows, std::string table_name){ - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/run_plan_table.txt",std::ios_base::trunc); - // for(const auto& atable : all_tables){ - // //ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // for(const auto& entry : atable.second){ - // entry.Print(ofs); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - // std::ofstream json_ofile("tables/run_plan_table.json",std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - { - double total_time = 0.0; - std::ofstream ofs(table_name+".txt", std::ios_base::app); - std::ofstream wiki_file(table_name+".wiki", std::ios_base::app); - //for (const auto& atable : all_tables) { - // ofs << "Q2 = " << atable.first << "\n"; - double Q2group_total_time = 0.0; - - double Q2 = table_rows.at(0).kinematic.Q2; - - fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", Q2)); - fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", Q2)); - - RunPlanTableEntry::PrintHeader(ofs); - RunPlanTableEntry::PrintWikiHeader(wiki_file); - for (const auto& entry : table_rows) { - if (std::abs(entry.polarity) >= 2) { - continue; - } - entry.Print(ofs); - entry.PrintWiki(wiki_file); - total_time += entry.time; - Q2group_total_time += entry.time; - } - RunPlanTableEntry::PrintWikiFooter(wiki_file); - - fmt::print(ofs, "{:-<79s}\n", "-"); - fmt::print(wiki_file, "{:-<79s}\n", "-"); - ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", Q2) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - wiki_file << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", Q2) - << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - Q2group_total_time / 24.0) - << "\n"; - //} - //fmt::print(ofs, "{:=<79s}\n", "="); - //ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - //fmt::print(wiki_file, "{:=<79s}\n", "="); - //wiki_file << " total time: " << total_time << " hours or " << total_time / 24.0 - // << " days\n"; - - //std::ofstream json_ofile("tables/LD2_run_plan_table.json", std::ios_base::trunc); - //json_ofile << TBufferJSON::ToJSON(&all_tables); - //json_ofile.close(); - } - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/LD2_run_plan_kaon_table.txt", std::ios_base::trunc); - // std::ofstream wiki_file("tables/LD2_run_plan_kaon_table.wiki", std::ios_base::trunc); - // for (const auto& atable : all_kaon_tables) { - // // ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - // fmt::print(wiki_file, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2", atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // RunPlanTableEntry::PrintWikiHeader(wiki_file); - // for (const auto& entry : atable.second) { - // entry.Print(ofs); - // entry.PrintWiki(wiki_file); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // RunPlanTableEntry::PrintWikiFooter(wiki_file); - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time / 24.0 << " days\n"; - - // std::ofstream json_ofile("tables/LD2_run_plan_kaon_table.json", std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - //{ - // double total_time = 0.0; - // std::ofstream ofs("tables/run_plan_kaon_table.txt",std::ios_base::trunc); - // for(const auto& atable : all_kaon_tables){ - // //ofs << "Q2 = " << atable.first << "\n"; - // double Q2group_total_time = 0.0; - // fmt::print(ofs, "{:^80}\n", fmt::format("Q2 = {:3.2f} GeV2",atable.first)); - // RunPlanTableEntry::PrintHeader(ofs); - // for(const auto& entry : atable.second){ - // entry.Print(ofs); - // total_time += entry.time; - // Q2group_total_time += entry.time; - // } - // fmt::print(ofs, "{:-<79s}\n", "-"); - // ofs << fmt::format(" Time for Q2 = {:3.2f} GeV2 setting : ", atable.first) - // << fmt::format("{:6.3f} hrs or {:6.3f} days\n", Q2group_total_time, - // Q2group_total_time / 24.0); - // ofs << "\n"; - // } - // fmt::print(ofs, "{:=<79s}\n", "="); - // ofs << " total time: " << total_time << " hours or " << total_time/24.0 << " days\n"; - - // std::ofstream json_ofile("tables/run_plan_kaon_table.json",std::ios_base::trunc); - // json_ofile << TBufferJSON::ToJSON(&all_tables); - // json_ofile.close(); - //} - } - - /** Compute the cross section and rate for a given setting - * - */ - std::vector<double> sidis_config_sigma_rate(RunPlanTableEntry& run_setting) { - - hallc::Kinematic& kine_set = run_setting.kinematic; - hallc::HallCSetting& hcSet = run_setting.hcSet; - double x_set = kine_set.x; - double z_set = kine_set.z; - double nu_set = kine_set.nu; - double W_set = kine_set.W; - double Wp_set = kine_set.Wp; - int polarity = run_setting.polarity; - hcSet.HMS_theta = kine_set.th_e * degree; - hcSet.SHMS_theta = kine_set.th_q * degree; - hcSet.HMS_p0 = kine_set.Ee; - hcSet.SHMS_p0 = kine_set.Ppi; - - // SHMS 4 msr dp -10 +20 - // HMS 6 msr dp +-8 - // ------------------------------------------------------- - // Phase space variables - IPSV phi_SHMS_psv({hcSet.SHMS_phi_min(), hcSet.SHMS_phi_max()}, "phi_shms"); - IPSV phi_HMS_psv({hcSet.HMS_phi_min(), hcSet.HMS_phi_max()}, "phi_hms"); - IPSV theta_HMS_psv({hcSet.HMS_theta_min(), hcSet.HMS_theta_max()}, "theta_hms"); - IPSV theta_SHMS_psv({hcSet.SHMS_theta_min(), hcSet.SHMS_theta_max()}, "theta_shms"); - IPSV P_HMS_psv({hcSet.HMS_P_min(), hcSet.HMS_P_max()}, "P_hms"); - IPSV P_SHMS_psv({hcSet.SHMS_P_min(), hcSet.SHMS_P_max()}, "P_shms"); - - auto diff_spectrometers = - make_diff(P_HMS_psv, theta_HMS_psv, phi_HMS_psv, P_SHMS_psv, theta_SHMS_psv, phi_SHMS_psv); - auto ps_spectrometers = make_phase_space(diff_spectrometers); - // diff_spectrometers.Print(); - - // ------------------------------ - // Prepare the initial state - double P1 = kine_set.E0; - double m1 = 0.000511; - double E1 = std::sqrt(P1 * P1 + m1 * m1); - - double P2 = 0.0; - double m2 = M_p / GeV; - double E2 = std::sqrt(P2 * P2 + m2 * m2); - - InitialState init_state(P1, P2, m2); - - int xs_select = 1; - PDFs pdf; - FragFs ffs; - - /** SIDIS diferential cross section. - * - * d(sigma) - * -------------------------- - * d(x)d(y)d(z)d(phi)d(phi_e) - * - * Note: "phi" is the angle between the leptonic and hadronic planes. - * See http://inspirehep.net/record/677636 for details. - * "phi_e" is the scattered electron's azimuthal angle. - * - */ - auto SIDIS_xs_func = [&](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - double Mh = M_pion / GeV; - - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - if (!(kine.isValid)) - return 0.0; - - // Gaussian k_perp and p_perp distributions for PT dependence - double kperp2_mean = 0.28; // GeV^2 - double pperp2_mean = 0.25; // GeV^2 - double PT_mean = kine.z * kine.z * kperp2_mean + pperp2_mean; - double PT2 = kine.P_hPerp * kine.P_hPerp; - auto f_PT = std::exp(-PT2 / PT_mean) / (PT_mean * pi); - auto t1 = kine.SIDIS_xs_term() * f_PT; - double FF = 0.0; - auto pdf_vals = pdf.Calculate(kine.x, kine.Q2); - - if (polarity == 1) { - // pi+ - if (xs_select == 1 || xs_select == 2) { - // proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -1) { - // pi- - if (xs_select == 1 || xs_select == 2) { - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_piminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == 2) { - // pi+ - if (xs_select == 1 || xs_select == 2) { - // proton - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kplus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } else if (polarity == -2) { - // pi- - if (xs_select == 1 || xs_select == 2) { - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - if (xs_select == 0 || xs_select == 2) { - // neutron - FF += PartonCharge2[q_id(Parton::u)] * ffs.D_u_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::d)]; - FF += PartonCharge2[q_id(Parton::ubar)] * ffs.D_ubar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::dbar)]; - FF += PartonCharge2[q_id(Parton::d)] * ffs.D_d_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::u)]; - FF += PartonCharge2[q_id(Parton::dbar)] * ffs.D_dbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::ubar)]; - FF += PartonCharge2[q_id(Parton::s)] * ffs.D_s_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::s)]; - FF += PartonCharge2[q_id(Parton::sbar)] * ffs.D_sbar_Kminus(kine.z, kine.Q2) * - pdf_vals[q_id(Parton::sbar)]; - } - } - FF *= kine.x; - - /// \f$ F_{UU,T} = x \sum e_a^2 f_1^a(x) D_1^a(z) \f$ - // double res = Elastic_XS({GE2, GM2}, kine); - // res = res * hbarc2_gev_nb * kine.sin_th; - return (t1 * FF); - }; - - // ------------------------------------------------------- - // - auto Eprime_and_Ph_func = [](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - XYZVector l2(Polar3D<double>(kine.P_l2, kine.theta_l2, kine.phi_l2)); - XYZVector p_h(Polar3D<double>(kine.P_p2, kine.theta_p2, kine.phi_p2)); - auto E_prime = ROOT::Math::XYZTVector(l2.x(), l2.y(), l2.z(), kine.E_l2); - auto P_hadron = ROOT::Math::XYZTVector(p_h.x(), p_h.y(), p_h.z(), kine.E_p2); - return std::array<ROOT::Math::XYZTVector, 2>({E_prime, P_hadron}); - }; - - // returns SIDIS_vars using the SPEC_vars (measured by spectrometers) - auto to_sidis_variables = [=](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SIDIS_vars()); - // kine.Print(); - std::array<double, 6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto to_spectrometer_variables = [=](const InitialState& is, const std::array<double, 6>& x) { - using namespace insane::units; - using namespace ROOT::Math; - double Mh = M_pion / GeV; - SIDISKinematics kine(is.p1().E(), Mh, x, SPEC_vars()); - std::array<double, 6> res = {kine.x, kine.y, kine.z, kine.P_hPerp, kine.phi_h, kine.phi_q}; - return res; - }; - - auto diff_xs_sidis = make_diff_cross_section(ps_spectrometers, SIDIS_xs_func); - - auto fs_parts = - make_final_state_particles({11, 211}, {0.000511, M_pion / GeV}, Eprime_and_Ph_func); - - auto fs_SIDIS = make_final_state(ps_spectrometers, fs_parts); - - auto p_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - auto n_integrated_XS = make_integrated_cross_section(init_state, diff_xs_sidis); - // auto sampler = make_ps_sampler(p_integrated_XS); - // auto sampler2 = make_ps_sampler(n_integrated_XS); - // sampler.SetFoamCells(500); - // sampler.SetFoamSample(100); - // sampler.SetFoamChat(0); - // sampler2.SetFoamCells(500); - // sampler2.SetFoamSample(100); - // sampler2.SetFoamChat(0); - - // auto evgen = make_event_generator(sampler, fs_SIDIS); - xs_select = 1; - auto total_XS_p = p_integrated_XS.CalculateTotalXS(); - // auto total_XS_p = evgen.Init(); - - // auto evgen2 = make_event_generator(sampler2, fs_SIDIS); - xs_select = 0; - auto total_XS_n = n_integrated_XS.CalculateTotalXS(); - // auto total_XS_n = evgen2.Init(); - - double total_XS = total_XS_n + total_XS_p; - - //// -------------------------- - // - // std::cout << "Total cross section = " << total_XS << " nb\n"; - int N_sim = 100000; - - using namespace insane; - double Loop1_LD2_Lumi = materials::ComputeLuminosityPerN( - {1, 2}, materials::LD2_density /*g/cm3*/, 10.0 /*cm*/, 50.0e-6); - double Loop1_LH2_Lumi = materials::ComputeLuminosityPerN( - {1, 1}, materials::LH2_density /*g/cm3*/, 10.0 /*cm*/, 50.0e-6); - double Loop1_entrance_window_thickness = 0.0104; - double Loop1_exit_window_thickness = 0.0133; - double Loop3_entrance_window_thickness = 0.0130; - double Loop3_exit_window_thickness = 0.0188; - - double Loop1_windows_Lumi = materials::ComputeLuminosityPerN( - {13, 27}, materials::Al_density, - Loop1_entrance_window_thickness + Loop1_exit_window_thickness, 50.0e-6); - - // double Lumi = Loop1_LD2_Lumi + Loop1_windows_Lumi; - - TargetRates targ_rates; - - targ_rates.LD2_Lumi = Loop1_LD2_Lumi; - targ_rates.LH2_Lumi = Loop1_LH2_Lumi; - targ_rates.window_Lumi = Loop1_windows_Lumi; - - // double total_cross_section = 2*total_XS_p * nanobarn; - // double luminosity = Lumi / cm2; - targ_rates.p_XS = total_XS_p; - targ_rates.n_XS = total_XS_n; - targ_rates.LD2_XS = total_XS_n + total_XS_p; - targ_rates.LH2_XS = total_XS_p; - targ_rates.window_XS = (total_XS_n * 14 + total_XS_p * 13); - - targ_rates.LD2_rate = targ_rates.LD2_XS * targ_rates.LD2_Lumi * (nanobarn / cm2); - targ_rates.LH2_rate = targ_rates.LH2_XS * targ_rates.LH2_Lumi * (nanobarn / cm2); - targ_rates.window_rate = targ_rates.window_XS * targ_rates.window_Lumi * (nanobarn / cm2); - targ_rates.total_LD2_rate = targ_rates.window_rate + targ_rates.LD2_rate; - targ_rates.total_LH2_rate = targ_rates.window_rate + targ_rates.LH2_rate; - targ_rates.total_rate = targ_rates.window_rate + targ_rates.LD2_rate; - - //std::cout << "Cross sections : \n"; - //std::cout << " [p] " << targ_rates.p_XS << " nb \n"; - //std::cout << " [n] " << targ_rates.n_XS << " nb \n"; - //std::cout << " [LD2] " << targ_rates.LD2_XS << " nb \n"; - //std::cout << " [LH2] " << targ_rates.LH2_XS << " nb \n"; - //std::cout << " [windows] " << targ_rates.window_XS << " nb \n"; - //// std::cout << " : " << total_cross_section/microbarn << " ub\n"; - //// std::cout << " : " << total_cross_section/barn << " b\n"; - //std::cout << "Luminosities : \n"; - //std::cout << " [LH2] " << targ_rates.LH2_Lumi << " 1/cm2 s\n"; - //std::cout << " [LD2] " << targ_rates.LD2_Lumi << " 1/cm2 s\n"; - //std::cout << " [windows] " << targ_rates.window_Lumi << " 1/cm2 s\n"; - //std::cout << "Rates : \n"; - //std::cout << " [LD2] " << targ_rates.LD2_rate << " 1/s\n"; - //std::cout << " [LH2] " << targ_rates.LH2_rate << " 1/s\n"; - //std::cout << " [windows] " << targ_rates.window_rate << " 1/s\n"; - //std::cout << "Total Rates : \n"; - //std::cout << " [LD2] " << targ_rates.total_LD2_rate << " 1/s\n"; - //std::cout << " [LH2] " << targ_rates.total_LH2_rate << " 1/s\n"; - //std::cout << " " << targ_rates.total_rate << " 1/s\n"; - - run_setting.rates = targ_rates; - // run_setting.LD2_Luminosity = Loop1_LD2_Lumi; - // run_setting.window_Luminosity = Loop1_windows_Lumi; - // run_setting.LD2_sigma = total_xs_LD2; - // run_setting.window_sigma = total_xs_windows; - std::vector<double> res = {targ_rates.LD2_XS, targ_rates.LD2_rate, - targ_rates.total_LD2_rate, targ_rates.total_LH2_rate, - targ_rates.LH2_XS, targ_rates.LH2_rate}; - return res; - } - -} // namespace hallc diff --git a/src/runplan/include/runplan/CSV_settings.h b/src/runplan/tests/CSV_settings.h similarity index 100% rename from src/runplan/include/runplan/CSV_settings.h rename to src/runplan/tests/CSV_settings.h