weight_sum.cpp 7.59 KB
Newer Older
1
/**
Tobias WEBER's avatar
Tobias WEBER committed
2
 * Calculates the integrated weights/energies for setup 3
3
4
5
6
7
8
9
 * @author Tobias Weber <tweber@ill.fr>
 * @date jun-20
 * @license GPLv2 (see 'LICENSE' file)
 */

#include "core/heli.h"
#include "core/skx.h"
Tobias WEBER's avatar
Tobias WEBER committed
10
#include "tlibs2/libs/phys.h"
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#include <fstream>

#include <boost/histogram.hpp>
namespace hist = boost::histogram;

using t_real = double;
using t_cplx = std::complex<t_real>;
using t_vec = ublas::vector<t_real>;
const auto g_j = t_cplx(0,1);

#include "core/skx_default_gs.cxx"


Tobias WEBER's avatar
Tobias WEBER committed
25
26
27
const t_real g_T = 28.5;
const t_real g_eps = 1e-5;

28
#define E_BINS 200
Tobias WEBER's avatar
Tobias WEBER committed
29
30
31
#define NUM_ANGLES 512
#define COL_SIZE 15

32
33
34
35
36

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,
Tobias WEBER's avatar
Tobias WEBER committed
37
	t_real q, int iProj=1)
38
39
{
	Skx<t_real, t_cplx, DEF_SKX_ORDER> skx;
Tobias WEBER's avatar
Tobias WEBER committed
40
	Heli<t_real, t_cplx, DEF_HELI_ORDER> heli;
41
42
43
44
45
46
47
48

	std::vector<ublas::vector<t_cplx>> 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<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));

	skx.SetFourier(fourier_skx);
Tobias WEBER's avatar
Tobias WEBER committed
49

50
	skx.SetProjNeutron(iProj!=0);
Tobias WEBER's avatar
Tobias WEBER committed
51
52
	heli.SetProjNeutron(iProj!=0);

53
	skx.SetT(-1000.);
Tobias WEBER's avatar
Tobias WEBER committed
54
55
	heli.SetT(g_T);

Tobias WEBER's avatar
Tobias WEBER committed
56
	skx.SetB(25.);
Tobias WEBER's avatar
Tobias WEBER committed
57
58
	heli.SetB(0.17);

59
	skx.GenFullFourier();
Tobias WEBER's avatar
Tobias WEBER committed
60

61
 	skx.SetFilterZeroWeight(1);
Tobias WEBER's avatar
Tobias WEBER committed
62
63
	heli.SetFilterZeroWeight(1);

64
	skx.SetWeightEps(1e-6);
Tobias WEBER's avatar
Tobias WEBER committed
65
	heli.SetWeightEps(1e-6);
66
67

	skx.SetCoords(Bx,By,Bz, Px,Py,Pz);
Tobias WEBER's avatar
Tobias WEBER committed
68
69
	heli.SetCoords(Bx,By,Bz);

70
	skx.SetG(Gx, Gy, Gz);
Tobias WEBER's avatar
Tobias WEBER committed
71
72
73
	heli.SetG(Gx, Gy, Gz);


74
75
	t_real Erange = 0.1;

Tobias WEBER's avatar
Tobias WEBER committed
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
102
103
104
	t_real angle_offs = 90/180.*M_PI;
	t_real angle_begin = -135/180.*M_PI + angle_offs;
	t_real angle_end = 135/180.*M_PI + angle_offs;
	t_real angle_delta = 2*M_PI/t_real(NUM_ANGLES);

	auto histWeightsNSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
	auto histWeightsSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));

	auto histWeightsHeliNSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
	auto histWeightsHeliSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));


	std::ofstream ofstr_raw("weightsum_skx.dat");
	std::ofstream ofstr_raw_heli("weightsum_heli.dat");

    // write file header
    for(std::ostream* ostr : {&ofstr_raw, &ofstr_raw_heli})
    {
        ostr->precision(8);
        (*ostr)  << std::left << std::setw(COL_SIZE) << "# angle"
            << " " << std::left << std::setw(COL_SIZE) << "qh"
            << " " << std::left << std::setw(COL_SIZE) << "qk"
            << " " << std::left << std::setw(COL_SIZE) << "ql"
            << " " << std::left << std::setw(COL_SIZE) << "E"
            << " " << std::left << std::setw(COL_SIZE) << "wSF1"
            << " " << std::left << std::setw(COL_SIZE) << "wSF2"
            << " " << std::left << std::setw(COL_SIZE) << "wNSF"
            << "\n";
    }
105
106
107
108
109
110
111
112


	for(t_real angle=angle_begin; angle<angle_end; angle+=angle_delta)
	{
		t_real Qx = q * cos(angle) + Gx;
		t_real Qy = q * sin(angle) + Gy;
		t_real Qz = Gz;

Tobias WEBER's avatar
Tobias WEBER committed
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
		std::cout << "# angle: " << angle/M_PI*180. << ", Q = (" << Qx << ", " << Qy << ", " << Qz << ")" << std::endl;

		{
			auto [Es, wsUnpol, wsSF1, wsSF2, wsNSF] = skx.GetDisp(Qx, Qy, Qz, -Erange, Erange);
			for(std::size_t i=0; i<Es.size(); ++i)
			{
				histWeightsNSF(Es[i], hist::weight(wsNSF[i]*0.5));
				histWeightsSF(Es[i], hist::weight(wsSF1[i]));

                ofstr_raw << std::left << std::setw(COL_SIZE) << angle
                    << " " << std::left << std::setw(COL_SIZE) << (Qx-Gx)
                    << " " << std::left << std::setw(COL_SIZE) << (Qy-Gy)
                    << " " << std::left << std::setw(COL_SIZE) << (Qz-Gz)
                    << " " << std::left << std::setw(COL_SIZE) << Es[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsSF1[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsSF2[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsNSF[i]
                    << std::endl;
			}
		}

		{
			auto [EsH, wsUnpolH, wsSF1H, wsSF2H, wsNSFH] = heli.GetDisp(Qx, Qy, Qz, -Erange, Erange);
			for(std::size_t i=0; i<EsH.size(); ++i)
			{
				histWeightsHeliNSF(EsH[i], hist::weight(wsNSFH[i]*0.5));
				histWeightsHeliSF(EsH[i], hist::weight(wsSF1H[i]));

                ofstr_raw_heli << std::left << std::setw(COL_SIZE) << angle
                    << " " << std::left << std::setw(COL_SIZE) << (Qx-Gx)
                    << " " << std::left << std::setw(COL_SIZE) << (Qy-Gy)
                    << " " << std::left << std::setw(COL_SIZE) << (Qz-Gz)
                    << " " << std::left << std::setw(COL_SIZE) << EsH[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsSF1H[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsSF2H[i]
                    << " " << std::left << std::setw(COL_SIZE) << wsNSFH[i]
                    << std::endl;
			}
		}
152
153
154
	}


Tobias WEBER's avatar
Tobias WEBER committed
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
	std::ofstream ofstrBinned("../data/weightbin.dat");
	ofstrBinned.precision(8);

	ofstrBinned << std::left << std::setw(COL_SIZE) << "# E (meV)"
		<< " " << std::left << std::setw(COL_SIZE) << "bose"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_skx_sf"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_skx_nsf"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_skx_sum"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_skx_sum_bose"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_heli_sf"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_heli_nsf"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_heli_sum"
		<< " " << std::left << std::setw(COL_SIZE) << "weight_heli_sum_bose"
		<< "\n";

	auto iterHeliNSF = boost::histogram::indexed(histWeightsHeliNSF).begin();
	auto iterHeliSF = boost::histogram::indexed(histWeightsHeliSF).begin();
	auto iterNSF = boost::histogram::indexed(histWeightsNSF).begin();
	for(const auto& val : boost::histogram::indexed(histWeightsSF))
174
175
	{
		t_real E = val.bin().lower() + 0.5*(val.bin().upper() - val.bin().lower());
Tobias WEBER's avatar
Tobias WEBER committed
176
		t_real w = *val / t_real{E_BINS};
177

Tobias WEBER's avatar
Tobias WEBER committed
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
		t_real E_nsf = iterNSF->bin().lower() + 0.5*(iterNSF->bin().upper() - iterNSF->bin().lower());
		t_real w_nsf = **iterNSF / t_real{E_BINS};

		t_real E_h = iterHeliSF->bin().lower() + 0.5*(iterHeliSF->bin().upper() - iterHeliSF->bin().lower());
		t_real w_h = **iterHeliSF / t_real{E_BINS};

		t_real E_h_nsf = iterHeliNSF->bin().lower() + 0.5*(iterHeliNSF->bin().upper() - iterHeliNSF->bin().lower());
		t_real w_h_nsf = **iterHeliNSF / t_real{E_BINS};

		if(!tl2::float_equal<t_real>(E, E_h, g_eps) || !tl2::float_equal<t_real>(E, E_h_nsf, g_eps) ||
			!tl2::float_equal<t_real>(E, E_nsf, g_eps))

		{
			std::cerr << "Energy binning mismatch: " << E << " != " << E_h << std::endl;
			break;
		}

		t_real bose = tl2::bose(E, g_T);

		ofstrBinned << std::left << std::setw(COL_SIZE) << E
			<< " " << std::left << std::setw(COL_SIZE) << bose
			<< " " << std::left << std::setw(COL_SIZE) << w
			<< " " << std::left << std::setw(COL_SIZE) << w_nsf
			<< " " << std::left << std::setw(COL_SIZE) << (w + w_nsf)
			<< " " << std::left << std::setw(COL_SIZE) << (w + w_nsf)*bose
			<< " " << std::left << std::setw(COL_SIZE) << w_h
			<< " " << std::left << std::setw(COL_SIZE) << w_h_nsf
			<< " " << std::left << std::setw(COL_SIZE) << (w_h + w_h_nsf)
			<< " " << std::left << std::setw(COL_SIZE) << (w_h + w_h_nsf)*bose
			<< std::endl;

		++iterHeliSF;
		++iterHeliNSF;
		++iterNSF;
212
213
214
215
216
217
	}
}


int main()
{
Tobias WEBER's avatar
Tobias WEBER committed
218
    // equivalent to the used setup around (000)
Tobias WEBER's avatar
Tobias WEBER committed
219
	t_real Gx = 1., Gy = 1., Gz = 0.;
220
221
222
223
224
225
226
	t_real Bx = 0., By = 0., Bz = 1.;
	t_real q = 0.0123;
	int proj = 1;

	calc_disp(Gx,Gy,Gz, Bx,By,Bz, 1,0,0, q, proj);
	return 0;
}