weight_sum.cpp 2.58 KB
 Tobias WEBER committed Jun 26, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 ``````/** * Calculates the integrated weights/energies for the reseda experiment * @author Tobias Weber * @date jun-20 * @license GPLv2 (see 'LICENSE' file) */ #include "core/heli.h" #include "core/skx.h" #include #include #include #include namespace hist = boost::histogram; using t_real = double; using t_cplx = std::complex; using t_vec = ublas::vector; const auto g_j = t_cplx(0,1); #include "core/skx_default_gs.cxx" #define E_BINS 200 void calc_disp( t_real Gx, t_real Gy, t_real Gz, t_real Bx, t_real By, t_real Bz, t_real Px, t_real Py, t_real Pz, t_real q, int iProj=1) { Skx skx; std::vector> fourier_skx; fourier_skx.reserve(_skxgs_allcomps.size()/3); for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3) fourier_skx.push_back(tl2::make_vec>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]})); skx.SetFourier(fourier_skx); skx.SetProjNeutron(iProj!=0); skx.SetT(-1000.); skx.SetB(25.); // BC2 = 40.3425 skx.GenFullFourier(); skx.SetFilterZeroWeight(1); skx.SetWeightEps(1e-6); skx.SetCoords(Bx,By,Bz, Px,Py,Pz); skx.SetG(Gx, Gy, Gz); t_real Erange = 0.1; t_real angle_begin = -135/180.*tl2::pi; t_real angle_end = 135/180.*tl2::pi; t_real angle_delta = 2*tl2::pi/100.; auto histWeights = hist::make_histogram(hist::axis::regular(E_BINS, -Erange, Erange, "E")); for(t_real angle=angle_begin; angle*180. << ", Q = (" << Qx << ", " << Qy << ", " << Qz << ")\n"; auto [Es, wsUnpol, wsSF1, wsSF2, wsNSF] = skx.GetDisp(Qx, Qy, Qz, -Erange, Erange); for(std::size_t i=0; i