BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
SamplerAnalysis.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2022.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "SamplerAnalysis.hh"
20#include "rebdsim.hh"
21
22#include <cmath>
23#include <vector>
24
25double SamplerAnalysis::particleMass = 0;
26
27SamplerAnalysis::SamplerAnalysis():
28 s(nullptr),
29 S(0),
30 debug(false)
31{
32 CommonCtor();
33}
34
35#ifndef __ROOTDOUBLE__
36SamplerAnalysis::SamplerAnalysis(BDSOutputROOTEventSampler<float> *samplerIn,
37 bool debugIn):
38#else
40 bool debugIn):
41#endif
42 s(samplerIn),
43 S(0),
44 debug(debugIn)
45{
46 CommonCtor();
47}
48
50{
51 int id = s->s->partID[0];
52 switch (std::abs(id))
53 {
54 case 11:
55 case -11:
56 {
57 std::cout << "Primary particle: e-" << std::endl;
58 particleMass = 0.000510999; break;}
59 case 2212:
60 {
61 std::cout << "Primary particle: proton" << std::endl;
62 particleMass = 0.938272; break;}
63 default:
64 {
65 std::cout << "Primary particle: unknown -> m=0" << std::endl;
66 particleMass = 0; break;}
67 }
68}
69
70void SamplerAnalysis::UpdateMass(const std::string& particleName)
71{
72 if ((particleName == "e-") || (particleName == "e+"))
73 {
74 std::cout << "Primary particle: e-" << std::endl;
75 particleMass = 0.000510999;;
76 }
77 else if (particleName == "proton")
78 {
79 std::cout << "Primary particle: proton" << std::endl;
80 particleMass = 0.938272;
81 }
82 else
83 {
84 std::cout << "Primary particle: unknown -> m=0" << std::endl;
85 particleMass = 0;
86 }
87}
88
90{
91 npart = 0;
92
93 // initialise a vector to store the coordinates of every event in a sampler
94 coordinates.resize(6, 0);
95
96 // initialise a vector to store the first values in a sampler for assumed mean subtraction
97 offsets.resize(6, 0);
98
99 optical.resize(3); // resize to 3 entries initialised to 0
100 varOptical.resize(3);
101 for(int i=0;i<3;++i)
102 {
103 optical[i].resize(25, 0); // 12 for central values and 12 for errors and 1 for xy correlation
104 varOptical[i].resize(12, 0);
105 }
106
107 covMats.resize(3);
108 for(int i=0;i<3;++i)
109 {
110 covMats[i].resize(3);
111 for(int j=0;j<3;++j)
112 {covMats[i][j].resize(3, 0);}
113 }
114
115 derivMats.resize(3);
116 for(int i=0;i<3;++i)
117 {
118 derivMats[i].resize(12);
119 for(int j=0;j<12;++j)
120 {derivMats[i][j].resize(3, 0);}
121 }
122
123 powSums.resize(6);
124 cenMoms.resize(6);
125
126 powSumsFirst.resize(6); // first sampler
127 cenMomsFirst.resize(6);
128
129 // (x,xp,y,yp,E,t) (x,xp,y,yp,E,t) v1pow, v2pow
130 for(int i=0;i<6;++i)
131 {
132 powSums[i].resize(6);
133 cenMoms[i].resize(6);
134
135 powSumsFirst[i].resize(6);
136 cenMomsFirst[i].resize(6);
137
138 for(int j=0;j<6;++j)
139 {
140 powSums[i][j].resize(5);
141 cenMoms[i][j].resize(5);
142
143 powSumsFirst[i][j].resize(5);
144 cenMomsFirst[i][j].resize(5);
145
146 for(int k=0;k<=4;++k)
147 {
148 powSums[i][j][k].resize(5, 0);
149 cenMoms[i][j][k].resize(5, 0);
150
151 powSumsFirst[i][j][k].resize(5, 0);
152 cenMomsFirst[i][j][k].resize(5, 0);
153 }
154 }
155 }
156}
157
158SamplerAnalysis::~SamplerAnalysis()
159{;}
160
162{
163 npart = 0;
164}
165
166void SamplerAnalysis::Process(bool firstTime)
167{
168 if(debug)
169 {std::cout << __METHOD_NAME__ << "\"" << s->samplerName << "\" with " << s->n << " entries" << std::endl;}
170
171 double m2 = std::pow(particleMass,2);
172
173 // loop over all entries
174 for(int i=0;i<s->n;++i)
175 {
176 if (i == 0)
177 {S = s->S;} // update sampler on first round - do here so not to load default data
178 if (s->parentID[i] != 0)
179 {continue;} // select only primary particles
180 if (s->turnNumber[i] > 1)
181 {continue;} // only use first turn particles
182 if (s->zp[i] <= 0)
183 {continue;} // only forward going particles - sampler can intercept backwards particles
184
185 coordinates[0] = s->x[i];
186 coordinates[1] = s->xp[i];
187 coordinates[2] = s->y[i];
188 coordinates[3] = s->yp[i];
189 coordinates[4] = std::sqrt(std::pow(s->energy[i],2) - m2); // p = sqrt(E^2 - M^2)
190 coordinates[5] = s->T[i];
191
192 if (firstTime)
193 {offsets = coordinates;}
194
195 // power sums
196 for(int a=0;a<6;++a)
197 {
198 for(int b=0;b<6;++b)
199 {
200 for (int j = 0; j <= 4; ++j)
201 {
202 for (int k = 0; k <= 4; ++k)
203 {
204 powSums[a][b][j][k] += std::pow(coordinates[a]-offsets[a],j)*std::pow(coordinates[b]-offsets[b],k);
205 }
206 }
207 }
208 }
209 npart++;
210 }
211}
212
213std::vector<double> SamplerAnalysis::Terminate(std::vector<double> emittance,
214 bool useEmittanceFromFirstSampler)
215{
216 if(debug)
217 {std::cout << " " << __METHOD_NAME__ << s->modelID << " " << npart << std::flush;}
218
219 // determine whether the input emittance is non-zero
220 bool nonZeroEmittanceIn = !std::all_of(emittance.begin(), emittance.end(), [](double l) { return l==0; });
221
222 // central moments
223 for(int a=0;a<6;++a)
224 {
225 for(int b=0;b<6;++b)
226 {
227 for (int j = 0; j <= 4; ++j)
228 {
229 for (int k = 0; k <= 4; ++k)
230 {
231 cenMoms[a][b][j][k] = powSumToCentralMoment(powSums, npart, a, b, j, k);
232 }
233 }
234 }
235 }
236
237 if (debug)
238 {printBeamCorrelationMatrix(cenMoms);}
239
240 // optical function calculation
241 for(int i=0;i<3;++i)
242 {
243 int j = 0;
244 if(i == 1) j = 2;
245 if(i == 2) j = 4;
246
247 // note: optical functions vector not populated in sequential order in order
248 // to apply dispersion correction to lattice funcs.
249
250 optical[i][6] = cenMoms[j][j+1][1][0]; // mean spatial (transv.)/ mean P (longit.)
251 optical[i][7] = cenMoms[j][j+1][0][1]; // mean divergence (transv.)/ mean t (longit.)
252 optical[i][8] = sqrt(cenMoms[j][j+1][2][0]); // sigma spatial (transv.)/ sigma P (longit.)
253 optical[i][9] = sqrt(cenMoms[j][j+1][0][2]); // sigma divergence (transv.)/ sigma t (longit.)
254
255 // transverse optical function calculation skipped for longitudinal plane,
256 // only mean and sigma of longit. coordinates recorded
257 if (i==2)
258 {continue;}
259
260 optical[i][4] = (cenMoms[4][4][1][0]*cenMoms[j][4][1][1])/cenMoms[4][4][2][0]; // eta
261 optical[i][5] = (cenMoms[4][4][1][0]*cenMoms[j+1][4][1][1])/cenMoms[4][4][2][0]; // eta prime
262
263 // check for nans (expected if dP=0) and set disp=0 if found
264 if (!std::isfinite(optical[i][4]))
265 {optical[i][4] = 0;}
266 if (!std::isfinite(optical[i][5]))
267 {optical[i][5] = 0;}
268
269 // temp variables to store dispersion corrected moments, aka moments without dispersion contribution
270 double corrCentMom_2_0 = 0.0, corrCentMom_1_1 = 0.0;
271 double corrCentMom_0_2 = 0.0;
272
273 //Correction example: var[x]_corrected = var[x] - eta_x^2*sigma[dP/P]^2
274 corrCentMom_2_0 =
275 cenMoms[j][j+1][2][0] - std::pow(optical[i][4],2)*cenMoms[4][4][2][0]/std::pow(cenMoms[4][4][1][0],2);
276
277 corrCentMom_0_2 =
278 cenMoms[j][j+1][0][2] - std::pow(optical[i][5],2)*cenMoms[4][4][2][0]/std::pow(cenMoms[4][4][1][0],2);
279
280 corrCentMom_1_1 =
281 cenMoms[j][j+1][1][1] - (optical[i][4]*optical[i][5])*cenMoms[4][4][2][0]/std::pow(cenMoms[4][4][1][0],2);
282
283 if(useEmittanceFromFirstSampler && nonZeroEmittanceIn)
284 {optical[i][0] = emittance[i];}
285 else
286 {optical[i][0] = sqrt(corrCentMom_2_0 * corrCentMom_0_2 - std::pow(corrCentMom_1_1, 2));} //emittance
287
288 optical[i][1] = -corrCentMom_1_1/optical[i][0]; // alpha
289 optical[i][2] = corrCentMom_2_0/optical[i][0]; // beta
290 optical[i][3] = (1+std::pow(optical[i][1],2))/optical[i][2]; // gamma
291
292 optical[i][10] = S;
293 optical[i][11] = npart;
294 }
295
296 // statistical error calculation
297
298 // covariance matrix of central moments for optical functions.
299 for(int i=0;i<3;++i)
300 {
301 for(int j=0;j<3;++j)
302 {
303 for(int k=0;k<3;++k)
304 {
305 covMats[i][j][k]=centMomToCovariance(cenMoms, npart, i, j, k);
306 }
307 }
308 }
309
310 // derivative matrix of parameters for optical functions. Each entry is
311 // a first order derivative w.r.t central moments.
312 for(int i=0;i<3;++i)
313 {
314 for(int j=0;j<12;++j) //loop over optical functions.
315 {
316 for(int k=0;k<3;++k) //loop over derivative indices
317 {
318 derivMats[i][j][k]=centMomToDerivative(cenMoms, i, j, k);
319 }
320 }
321 }
322
323 // compute variances of optical functions
324 for(int i=0;i<3;++i)
325 {
326 for(int j=0;j<12;++j)
327 {
328 for(int k=0;k<3;++k)
329 {
330 for(int l=0;l<3;++l)
331 {
332 varOptical[i][j] += derivMats[i][j][k]*derivMats[i][j][l]*covMats[i][k][l];
333 }
334 }
335 }
336 }
337
338 // compute the sigmas of optical functions
339 for(int i=0;i<3;++i)
340 {
341 for(int j=0;j<12;++j)
342 {
343 if(j==6 || j==7)
344 {optical[i][j+12]=optical[i][j+2]/sqrt(npart);} //errors of the means
345 else if(j == 10 || j == 11)
346 {optical[i][j+12]= 0.0;} //no errors on S and Npart
347 else if (j==0)
348 {
349 if(useEmittanceFromFirstSampler && nonZeroEmittanceIn)
350 {
351 optical[i][j+12]=emittance[j+2];
352 }
353 else
354 {
355 optical[i][j+12]=sqrt(varOptical[i][j]);
356 }
357 }
358 else
359 {optical[i][j+12]=sqrt(varOptical[i][j]);}
360 }
361 }
362
363 // convert meanP back to meanE
364 double m2 = std::pow(particleMass,2);
365 double meanP = optical[2][6];
366 optical[2][6] = std::sqrt(std::pow(meanP,2) + m2);
367
368 //convert sigmaP back to SigmaE
369 double gamma = optical[2][6] / particleMass;
370 double beta = std::sqrt(1.0 - 1.0/std::pow(gamma,2));
371 double sigmaP = optical[2][8];
372 optical[2][8] = (sigmaP/meanP)*std::pow(beta,2);
373
374 // write out the correlation x-y coefficient to the output as a matrix of horizontal-vertical coupling
375 // writen only to the x vector, but 0 is added to y and z vectors to keep all vector sizes the same
376 optical[0][24]=cenMoms[0][2][1][1]/std::sqrt(cenMoms[0][2][2][0]*cenMoms[0][2][0][2]);
377
378 // emitt_x, emitt_y, err_emitt_x, err_emitt_y
379 std::vector<double> emittanceOut = {optical[0][0],
380 optical[1][0],
381 optical[0][12],
382 optical[1][12]};
383
384 return emittanceOut;
385}
386
387double SamplerAnalysis::powSumToCentralMoment(fourDArray& powSumsIn,
388 long long int npartIn,
389 int a,
390 int b,
391 int m,
392 int n)
393{
394 double moment = 0.0;
395
396 // Explicitly multiply out the number of particles to achieve the appropriate power.
397 // Typecast to a double is done to avoid integer overflow for large number of particles.
398 double npartPow1 = (double)npartIn;
399
400 // Done because of potential problems with the pow() funtion in C++ where some
401 // arguments cause NaNs
402 double npartPow2 = npartPow1 * npartPow1;
403 double npartPow3 = npartPow2 * npartPow1;
404 double npartPow4 = npartPow3 * npartPow1;
405
406 if((m == 1 && n == 0) || (m == 0 && n == 1))
407 {
408 double s_1_0 = powSumsIn[a][b][m][n];
409 int k = m > n ? a : b;
410
411 moment = s_1_0/npartIn+offsets[k];
412 }
413
414 else if((n == 2 && m == 0) || (n == 0 && m == 2))
415 {
416 double s_1_0 = 0.0, s_2_0 = 0.0;
417 if(m == 2)
418 {
419 s_1_0 = powSumsIn[a][b][m-1][n];
420 s_2_0 = powSumsIn[a][b][m][n];
421 }
422 else if(n == 2)
423 {
424 s_1_0 = powSumsIn[a][b][m][n-1];
425 s_2_0 = powSumsIn[a][b][m][n];
426 }
427
428 moment = (npartPow1*s_2_0 - std::pow(std::abs(s_1_0),2))/(npartPow1*(npartPow1-1));
429 }
430
431 else if(n == 1 && m == 1)
432 {
433 double s_1_0 = 0.0, s_0_1 = 0.0, s_1_1 = 0.0;
434
435 s_1_0 = powSumsIn[a][b][m][n-1];
436 s_0_1 = powSumsIn[a][b][m-1][n];
437 s_1_1 = powSumsIn[a][b][m][n];
438
439 moment = (npartPow1*s_1_1 - s_0_1*s_1_0)/(npartPow1*(npartPow1-1));
440 }
441
442 else if((n == 4 && m == 0) || (n == 0 && m == 4))
443 {
444 double s_1_0 = 0.0, s_2_0 = 0.0, s_3_0 = 0.0, s_4_0 = 0.0;
445 if(m == 4)
446 {
447 s_1_0 = powSumsIn[a][b][m-3][n];
448 s_2_0 = powSumsIn[a][b][m-2][n];
449 s_3_0 = powSumsIn[a][b][m-1][n];
450 s_4_0 = powSumsIn[a][b][m][n];
451 }
452 else if( n == 4)
453 {
454 s_1_0 = powSumsIn[a][b][m][n-3];
455 s_2_0 = powSumsIn[a][b][m][n-2];
456 s_3_0 = powSumsIn[a][b][m][n-1];
457 s_4_0 = powSumsIn[a][b][m][n];
458 }
459
460 moment = - (3*std::pow(s_1_0,4))/npartPow4 + (6*std::pow(s_1_0,2)*s_2_0)/npartPow3
461 - (4*s_1_0*s_3_0)/npartPow2 + s_4_0/npartPow1;
462 }
463
464 else if((m == 3 && n == 1) || (m == 1 && n ==3))
465 {
466 double s_1_0 = 0.0, s_0_1 = 0.0, s_1_1 = 0.0, s_2_0 = 0.0, s_2_1 = 0.0, s_3_0 = 0.0, s_3_1 = 0.0;
467
468 if(m == 3)
469 {
470 s_1_0 = powSumsIn[a][b][m-2][n-1];
471 s_0_1 = powSumsIn[a][b][m-3][n];
472 s_1_1 = powSumsIn[a][b][m-2][n];
473 s_2_0 = powSumsIn[a][b][m-1][n-1];
474 s_2_1 = powSumsIn[a][b][m-1][n];
475 s_3_0 = powSumsIn[a][b][m][n-1];
476 s_3_1 = powSumsIn[a][b][m][n];
477 }
478 else if(n == 3)
479 {
480 s_1_0 = powSumsIn[a][b][m-1][n-2];
481 s_0_1 = powSumsIn[a][b][m][n-3];
482 s_1_1 = powSumsIn[a][b][m][n-2];
483 s_2_0 = powSumsIn[a][b][m-1][n-1];
484 s_2_1 = powSumsIn[a][b][m][n-1];
485 s_3_0 = powSumsIn[a][b][m-1][n];
486 s_3_1 = powSumsIn[a][b][m][n];
487 }
488
489 moment = - (3*s_0_1*std::pow(s_1_0,3))/npartPow4 + (3*s_1_0*s_1_0*s_1_1)/npartPow3
490 + (3*s_0_1*s_1_0*s_2_0)/npartPow3 - (3*s_1_0*s_2_1)/npartPow2
491 - (s_0_1*s_3_0)/npartPow2 + s_3_1/npartPow1;
492 }
493
494 else if(m == 2 && n == 2)
495 {
496 double s_1_0 = 0.0, s_0_1 = 0.0, s_1_1 = 0.0, s_2_0 = 0.0, s_0_2 = 0.0, s_1_2 = 0.0, s_2_1 = 0.0, s_2_2 = 0.0;
497
498 s_1_0 = powSumsIn[a][b][m-1][n-2];
499 s_0_1 = powSumsIn[a][b][m-2][n-1];
500 s_1_1 = powSumsIn[a][b][m-1][n-1];
501 s_2_0 = powSumsIn[a][b][m][n-2];
502 s_0_2 = powSumsIn[a][b][m-2][n];
503 s_1_2 = powSumsIn[a][b][m-1][n];
504 s_2_1 = powSumsIn[a][b][m][n-1];
505 s_2_2 = powSumsIn[a][b][m][n];
506
507 moment = - (3*std::pow(s_0_1,2)*std::pow(s_1_0,2))/npartPow4 + (s_0_2*std::pow(s_1_0,2))/npartPow3
508 + (4*s_0_1*s_1_0*s_1_1)/npartPow3 - (2*s_1_0*s_1_2)/npartPow2
509 + (std::pow(s_0_1,2)*s_2_0)/npartPow3 - (2*s_0_1*s_2_1)/npartPow2 + s_2_2/npartPow1;
510 }
511
512 return moment;
513}
514
515double SamplerAnalysis::centMomToCovariance(fourDArray& centMoms,
516 long long int npartIn,
517 int k,
518 int i,
519 int j)
520{
521 double cov = 0.0;
522
523 int a = 0;
524 if(k == 1) {a=2;}
525 if(k == 2) {a=4;}
526
527 if((i == 0 && j == 0) || (i == 1 && j == 1))
528 {
529 double m_4_0 = 0.0, m_2_0 = 0.0;
530
531 if(i == 0)
532 {
533 m_4_0 = centMoms[a][a+1][4][0];
534 m_2_0 = centMoms[a][a+1][2][0];
535 }
536
537 else if(i == 1)
538 {
539 m_4_0 = centMoms[a][a+1][0][4];
540 m_2_0 = centMoms[a][a+1][0][2];
541 }
542
543 cov = -((npartIn-3)*std::pow(m_2_0,2))/(npartIn*(npartIn-1)) + m_4_0/npartIn;
544 }
545
546 else if(i == 2 && j == 2)
547 {
548 double m_1_1 = 0.0, m_2_0 = 0.0, m_0_2 = 0.0, m_2_2 = 0.0;
549
550 m_1_1 = centMoms[a][a+1][1][1];
551 m_2_0 = centMoms[a][a+1][2][0];
552 m_0_2 = centMoms[a][a+1][0][2];
553 m_2_2 = centMoms[a][a+1][2][2];
554
555 cov = -((npartIn-2)*std::pow(m_1_1,2))/(npartIn*(npartIn-1)) + (m_0_2*m_2_0)/(npartIn*(npartIn-1)) + m_2_2/npartIn;
556 }
557
558 else if((i == 0 && j == 2) || (i == 1 && j == 2) || (i == 2 && j == 0) || (i == 2 && j == 1))
559 {
560 double m_1_1 = 0.0, m_2_0 = 0.0, m_3_1 = 0.0;
561
562 if((i == 0 && j == 2) || (i == 2 && j == 0 ))
563 {
564 m_1_1 = centMoms[a][a+1][1][1];
565 m_2_0 = centMoms[a][a+1][2][0];
566 m_3_1 = centMoms[a][a+1][3][1];
567 }
568 else if((i == 1 && j == 2) || (i == 2 && j == 1))
569 {
570 m_1_1 = centMoms[a][a+1][1][1];
571 m_2_0 = centMoms[a][a+1][0][2];
572 m_3_1 = centMoms[a][a+1][1][3];
573 }
574
575 cov = -((npartIn-3)*m_1_1*m_2_0)/(npartIn*(npartIn-1)) + m_3_1/npartIn;
576 }
577 else if((i == 0 && j == 1) || (i == 1 && j == 0) )
578 {
579 double m_1_1 = 0.0, m_2_0 = 0.0, m_0_2 = 0.0, m_2_2 = 0.0;
580
581 m_1_1 = centMoms[a][a+1][1][1];
582 m_2_0 = centMoms[a][a+1][0][2];
583 m_0_2 = centMoms[a][a+1][2][0];
584 m_2_2 = centMoms[a][a+1][2][2];
585
586 cov = 2*std::pow(m_1_1,2)/(npartIn*(npartIn-1)) - m_2_0*m_0_2/npartIn + m_2_2/npartIn;
587 }
588
589 return cov;
590}
591
592double SamplerAnalysis::centMomToDerivative(fourDArray& centMoms,
593 int k,
594 int t,
595 int i)
596
597{
598 double deriv = 0.0;
599
600 int a = 0;
601 if(k == 1) {a=2;}
602 if(k == 2) {a=4;}
603
604 switch(t)
605 {
606 case 0:
607 //emittance
608 if(i == 0 && k < 2) // k<2 check selects transverse planes, longitudinal parameters are not calculated
609 {
610 deriv = centMoms[a][a+1][0][2]/(2*sqrt(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2)));
611 }
612 else if(i == 1 && k < 2)
613 {
614 deriv = centMoms[a][a+1][2][0]/(2*sqrt(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2)));
615 }
616 else if(i == 2 && k < 2)
617 {
618 deriv = -centMoms[a][a+1][1][1]/(sqrt(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2)));
619 }
620 else {deriv=0;}
621
622 return deriv;
623 break;
624
625 case 1:
626 //alpha
627 if(i == 0 && k < 2)
628 {
629 deriv = centMoms[a][a+1][0][2]*centMoms[a][a+1][1][1]/(2*std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.));
630 }
631 else if(i == 1 && k < 2)
632 {
633 deriv = centMoms[a][a+1][2][0]*centMoms[a][a+1][1][1]/(2*std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.));
634 }
635 else if(i == 2 && k < 2)
636 {
637 deriv = - centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]/std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.);
638 }
639 else {deriv=0;}
640
641 return deriv;
642 break;
643
644 case 2:
645 //beta
646 if(i == 0 && k < 2)
647 {
648 deriv = (centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-2*std::pow(centMoms[a][a+1][1][1],2))/(2*std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.));
649 }
650 else if(i == 1 && k < 2)
651 {
652 deriv = - std::pow(centMoms[a][a+1][2][0],2)/(2*std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.));
653 }
654 else if(i == 2 && k < 2)
655 {
656 deriv = centMoms[a][a+1][2][0]*centMoms[a][a+1][1][1]/std::pow(centMoms[a][a+1][2][0]*centMoms[a][a+1][0][2]-std::pow(centMoms[a][a+1][1][1],2),3./2.);
657 }
658 else {deriv=0;}
659
660 return deriv;
661 break;
662
663
664 case 3:
665 //gamma
666 return 0; // TODO
667 break;
668
669 case 4:
670 //eta
671 if(i == 0 && k < 2)
672 {
673 deriv = 0;
674 }
675 else if(i == 1 && k < 2)
676 {
677 deriv = -centMoms[a][4][1][1]/std::pow(centMoms[4][4][2][0],2);
678 }
679 else if(i == 2 && k < 2)
680 {
681 deriv = 1/centMoms[4][4][2][0];
682 }
683 else {deriv=0;}
684
685 return deriv;
686 break;
687
688 case 5:
689 //eta prime
690 if(i == 0 && k < 2)
691 {
692 deriv = 0;
693 }
694 else if(i == 1 && k < 2)
695 {
696 deriv = -centMoms[a+1][4][1][1]/std::pow(centMoms[4][4][2][0],2);
697 }
698 else if(i == 2 && k < 2)
699 {
700 deriv = 1/centMoms[4][4][2][0];
701 }
702 else {deriv=0;}
703
704 return deriv;
705 break;
706
707 case 6: // mean derivatives are all 0 as they don't depend on second order moments
708 //mean spatial
709 return 0;
710 break;
711
712 case 7:
713 //mean divergence
714 return 0;
715 break;
716
717 case 8:
718 //sigma spatial
719 if(i == 0)
720 {
721 deriv = 1/(2*sqrt(centMoms[a][a+1][2][0]));
722 }
723
724 else {deriv=0;}
725
726 return deriv;
727 break;
728
729 case 9:
730 //sigma divergence
731 if(i == 1)
732 {
733 deriv = 1/(2*sqrt(centMoms[a][a+1][0][2]));
734 }
735
736 else {deriv=0;}
737
738 return deriv;
739 break;
740
741 default:
742 return 0;
743 break;
744 }
745}
746
747
748void SamplerAnalysis::printBeamCorrelationMatrix(fourDArray& centMoms)
749{
750 std::cout<<"\nCorrelation matrix for sampler: "<<s->samplerName<<std::endl;
751 double corr = 0.0;
752 for (int i=0; i<6; i++)
753 {
754 for(int j=0; j<6; j++)
755 {
756 corr = centMoms[i][j][1][1]/std::sqrt(centMoms[i][j][2][0]*centMoms[i][j][0][2]);
757 std::printf("%- *.6e ", 12, corr);
758 }
759 std::printf("\n");
760 }
761
762}
Analysis routines for an individual sampler.
twoDArray optical
emt, alf, bta, gma, eta, etapr, mean, sigma
twoDArray varOptical
variances of optical functions
double centMomToDerivative(fourDArray &centMoms, int k, int t, int i)
void Initialise()
Initialise variables.
static void UpdateMass(SamplerAnalysis *s)
Set primary particle mass for optical functions from sampler data.
void CommonCtor()
Initialisation of arrays for optical function calculations.
void Process(bool firstTime=false)
Loop over all entries in the sampler and accumulate power sums over variuos moments.
double centMomToCovariance(fourDArray &centMoms, long long int npartIn, int k, int i, int j)
double powSumToCentralMoment(fourDArray &powSum, long long int npartIn, int i, int j, int m, int n)
std::vector< double > Terminate(std::vector< double > emittance, bool useEmittanceFromFirstSampler=true)
Calculate optical functions based on combinations of moments already accumulated.