19#include "SamplerAnalysis.hh"
25double SamplerAnalysis::particleMass = 0;
27SamplerAnalysis::SamplerAnalysis():
51 int id = s->s->partID[0];
57 std::cout <<
"Primary particle: e-" << std::endl;
58 particleMass = 0.000510999;
break;}
61 std::cout <<
"Primary particle: proton" << std::endl;
62 particleMass = 0.938272;
break;}
65 std::cout <<
"Primary particle: unknown -> m=0" << std::endl;
66 particleMass = 0;
break;}
72 if ((particleName ==
"e-") || (particleName ==
"e+"))
74 std::cout <<
"Primary particle: e-" << std::endl;
75 particleMass = 0.000510999;;
77 else if (particleName ==
"proton")
79 std::cout <<
"Primary particle: proton" << std::endl;
80 particleMass = 0.938272;
84 std::cout <<
"Primary particle: unknown -> m=0" << std::endl;
94 coordinates.resize(6, 0);
110 covMats[i].resize(3);
112 {covMats[i][j].resize(3, 0);}
118 derivMats[i].resize(12);
119 for(
int j=0;j<12;++j)
120 {derivMats[i][j].resize(3, 0);}
126 powSumsFirst.resize(6);
127 cenMomsFirst.resize(6);
132 powSums[i].resize(6);
133 cenMoms[i].resize(6);
135 powSumsFirst[i].resize(6);
136 cenMomsFirst[i].resize(6);
140 powSums[i][j].resize(5);
141 cenMoms[i][j].resize(5);
143 powSumsFirst[i][j].resize(5);
144 cenMomsFirst[i][j].resize(5);
146 for(
int k=0;k<=4;++k)
148 powSums[i][j][k].resize(5, 0);
149 cenMoms[i][j][k].resize(5, 0);
151 powSumsFirst[i][j][k].resize(5, 0);
152 cenMomsFirst[i][j][k].resize(5, 0);
158SamplerAnalysis::~SamplerAnalysis()
169 {std::cout << __METHOD_NAME__ <<
"\"" << s->samplerName <<
"\" with " << s->n <<
" entries" << std::endl;}
171 double m2 = std::pow(particleMass,2);
174 for(
int i=0;i<s->n;++i)
178 if (s->parentID[i] != 0)
180 if (s->turnNumber[i] > 1)
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);
190 coordinates[5] = s->T[i];
193 {offsets = coordinates;}
200 for (
int j = 0; j <= 4; ++j)
202 for (
int k = 0; k <= 4; ++k)
204 powSums[a][b][j][k] += std::pow(coordinates[a]-offsets[a],j)*std::pow(coordinates[b]-offsets[b],k);
214 bool useEmittanceFromFirstSampler)
217 {std::cout <<
" " << __METHOD_NAME__ << s->modelID <<
" " << npart << std::flush;}
220 bool nonZeroEmittanceIn = !std::all_of(emittance.begin(), emittance.end(), [](
double l) { return l==0; });
227 for (
int j = 0; j <= 4; ++j)
229 for (
int k = 0; k <= 4; ++k)
238 {printBeamCorrelationMatrix(cenMoms);}
250 optical[i][6] = cenMoms[j][j+1][1][0];
251 optical[i][7] = cenMoms[j][j+1][0][1];
252 optical[i][8] = sqrt(cenMoms[j][j+1][2][0]);
253 optical[i][9] = sqrt(cenMoms[j][j+1][0][2]);
260 optical[i][4] = (cenMoms[4][4][1][0]*cenMoms[j][4][1][1])/cenMoms[4][4][2][0];
261 optical[i][5] = (cenMoms[4][4][1][0]*cenMoms[j+1][4][1][1])/cenMoms[4][4][2][0];
264 if (!std::isfinite(
optical[i][4]))
266 if (!std::isfinite(
optical[i][5]))
270 double corrCentMom_2_0 = 0.0, corrCentMom_1_1 = 0.0;
271 double corrCentMom_0_2 = 0.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);
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);
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);
283 if(useEmittanceFromFirstSampler && nonZeroEmittanceIn)
284 {
optical[i][0] = emittance[i];}
286 {
optical[i][0] = sqrt(corrCentMom_2_0 * corrCentMom_0_2 - std::pow(corrCentMom_1_1, 2));}
314 for(
int j=0;j<12;++j)
326 for(
int j=0;j<12;++j)
332 varOptical[i][j] += derivMats[i][j][k]*derivMats[i][j][l]*covMats[i][k][l];
341 for(
int j=0;j<12;++j)
345 else if(j == 10 || j == 11)
349 if(useEmittanceFromFirstSampler && nonZeroEmittanceIn)
351 optical[i][j+12]=emittance[j+2];
364 double m2 = std::pow(particleMass,2);
366 optical[2][6] = std::sqrt(std::pow(meanP,2) + m2);
369 double gamma =
optical[2][6] / particleMass;
370 double beta = std::sqrt(1.0 - 1.0/std::pow(gamma,2));
372 optical[2][8] = (sigmaP/meanP)*std::pow(beta,2);
376 optical[0][24]=cenMoms[0][2][1][1]/std::sqrt(cenMoms[0][2][2][0]*cenMoms[0][2][0][2]);
379 std::vector<double> emittanceOut = {
optical[0][0],
388 long long int npartIn,
398 double npartPow1 = (double)npartIn;
402 double npartPow2 = npartPow1 * npartPow1;
403 double npartPow3 = npartPow2 * npartPow1;
404 double npartPow4 = npartPow3 * npartPow1;
406 if((m == 1 && n == 0) || (m == 0 && n == 1))
408 double s_1_0 = powSumsIn[a][b][m][n];
409 int k = m > n ? a : b;
411 moment = s_1_0/npartIn+offsets[k];
414 else if((n == 2 && m == 0) || (n == 0 && m == 2))
416 double s_1_0 = 0.0, s_2_0 = 0.0;
419 s_1_0 = powSumsIn[a][b][m-1][n];
420 s_2_0 = powSumsIn[a][b][m][n];
424 s_1_0 = powSumsIn[a][b][m][n-1];
425 s_2_0 = powSumsIn[a][b][m][n];
428 moment = (npartPow1*s_2_0 - std::pow(std::abs(s_1_0),2))/(npartPow1*(npartPow1-1));
431 else if(n == 1 && m == 1)
433 double s_1_0 = 0.0, s_0_1 = 0.0, s_1_1 = 0.0;
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];
439 moment = (npartPow1*s_1_1 - s_0_1*s_1_0)/(npartPow1*(npartPow1-1));
442 else if((n == 4 && m == 0) || (n == 0 && m == 4))
444 double s_1_0 = 0.0, s_2_0 = 0.0, s_3_0 = 0.0, s_4_0 = 0.0;
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];
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];
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;
464 else if((m == 3 && n == 1) || (m == 1 && n ==3))
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;
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];
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];
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;
494 else if(m == 2 && n == 2)
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;
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];
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;
516 long long int npartIn,
527 if((i == 0 && j == 0) || (i == 1 && j == 1))
529 double m_4_0 = 0.0, m_2_0 = 0.0;
533 m_4_0 = centMoms[a][a+1][4][0];
534 m_2_0 = centMoms[a][a+1][2][0];
539 m_4_0 = centMoms[a][a+1][0][4];
540 m_2_0 = centMoms[a][a+1][0][2];
543 cov = -((npartIn-3)*std::pow(m_2_0,2))/(npartIn*(npartIn-1)) + m_4_0/npartIn;
546 else if(i == 2 && j == 2)
548 double m_1_1 = 0.0, m_2_0 = 0.0, m_0_2 = 0.0, m_2_2 = 0.0;
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];
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;
558 else if((i == 0 && j == 2) || (i == 1 && j == 2) || (i == 2 && j == 0) || (i == 2 && j == 1))
560 double m_1_1 = 0.0, m_2_0 = 0.0, m_3_1 = 0.0;
562 if((i == 0 && j == 2) || (i == 2 && j == 0 ))
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];
568 else if((i == 1 && j == 2) || (i == 2 && j == 1))
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];
575 cov = -((npartIn-3)*m_1_1*m_2_0)/(npartIn*(npartIn-1)) + m_3_1/npartIn;
577 else if((i == 0 && j == 1) || (i == 1 && j == 0) )
579 double m_1_1 = 0.0, m_2_0 = 0.0, m_0_2 = 0.0, m_2_2 = 0.0;
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];
586 cov = 2*std::pow(m_1_1,2)/(npartIn*(npartIn-1)) - m_2_0*m_0_2/npartIn + m_2_2/npartIn;
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)));
612 else if(i == 1 && k < 2)
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)));
616 else if(i == 2 && k < 2)
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)));
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.));
631 else if(i == 1 && k < 2)
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.));
635 else if(i == 2 && k < 2)
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.);
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.));
650 else if(i == 1 && k < 2)
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.));
654 else if(i == 2 && k < 2)
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.);
675 else if(i == 1 && k < 2)
677 deriv = -centMoms[a][4][1][1]/std::pow(centMoms[4][4][2][0],2);
679 else if(i == 2 && k < 2)
681 deriv = 1/centMoms[4][4][2][0];
694 else if(i == 1 && k < 2)
696 deriv = -centMoms[a+1][4][1][1]/std::pow(centMoms[4][4][2][0],2);
698 else if(i == 2 && k < 2)
700 deriv = 1/centMoms[4][4][2][0];
721 deriv = 1/(2*sqrt(centMoms[a][a+1][2][0]));
733 deriv = 1/(2*sqrt(centMoms[a][a+1][0][2]));
748void SamplerAnalysis::printBeamCorrelationMatrix(fourDArray& centMoms)
750 std::cout<<
"\nCorrelation matrix for sampler: "<<s->samplerName<<std::endl;
752 for (
int i=0; i<6; i++)
754 for(
int j=0; j<6; j++)
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);
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 ¢Moms, 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 ¢Moms, 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.