Commit 68b56040 authored by kraus's avatar kraus
Browse files

Squashed commit of the following:

commit 50e2a63e
Author: Christof Kraus <christof.j.kraus@gmail.com>
Date:   Sun Feb 28 08:42:19 2021 +0100

    add total charge and mass to the statistics

commit 8792c971
Author: Christof Kraus <christof.j.kraus@gmail.com>
Date:   Sat Feb 27 22:06:57 2021 +0100

    adding more information to loss data sink

commit 93c63b24
Author: Christof Kraus <christof.j.kraus@gmail.com>
Date:   Sun Oct 18 10:14:42 2020 +0200

    initial version for moments computer

commit 8631137e
Author: Christof Kraus <christof.j.kraus@gmail.com>
Date:   Sat Oct 10 23:07:39 2020 +0200

    add initial version
parent c972256f
......@@ -203,9 +203,11 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i],
dtime, intecoords, triId);
if(res >= 0) {
lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i],
itsBunch_m->ID[i], itsBunch_m->getT()*1e9,
turnnumber_m, itsBunch_m->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(itsBunch_m->ID[i],
itsBunch_m->R[i], itsBunch_m->P[i],
itsBunch_m->getT()*1e9,
itsBunch_m->Q[i], itsBunch_m->M[i]),
std::make_pair(turnnumber_m, itsBunch_m->bunchNum[i]));
itsBunch_m->Bin[i] = -1;
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
......@@ -478,7 +480,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << std::boolalpha << "* Superpose electric field maps -> " << superpose << endl;
}
// Read in cyclotron field maps
// Read in cyclotron field maps
cycl_m->initialise(itsBunch_m, cycl_m->getBScale());
double BcParameter[8] = {};
......@@ -1049,7 +1051,7 @@ void ParallelCyclotronTracker::visitVacuum(const Vacuum &vac) {
*gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;
elptr->initialise(itsBunch_m, elptr->getPScale());
BcParameter[1] = temperature;
BcParameter[2] = stop;
......@@ -3524,4 +3526,4 @@ void ParallelCyclotronTracker::initPathLength() {
// we need to reset the path length of each bunch
itsDataSink->setMultiBunchInitialPathLengh(mbHandler_m.get());
}
}
}
\ No newline at end of file
......@@ -35,7 +35,7 @@
#include "Beamlines/Beamline.h"
#include "Beamlines/FlaggedBeamline.h"
#include "Solvers/CSRWakeFunction.hh"
#include "Solvers/CSRWakeFunction.h"
#include "AbstractObjects/OpalData.h"
......@@ -47,7 +47,7 @@
#include "ValueDefinitions/RealVariable.h"
#include "Utilities/Timer.h"
#include "Utilities/OpalException.h"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/ParticleMatterInteractionHandler.h"
#include "Structure/BoundaryGeometry.h"
#include "AbsBeamline/Monitor.h"
......@@ -473,15 +473,14 @@ void ParallelTTracker::emitParticles(long long step) {
void ParallelTTracker::computeSpaceChargeFields(unsigned long long step) {
if (numParticlesInSimulation_m <= minBinEmitted_m) return;
if (!itsBunch_m->hasFieldSolver()) return;
if (numParticlesInSimulation_m <= minBinEmitted_m || !itsBunch_m->hasFieldSolver()) {
return;
}
itsBunch_m->calcBeamParameters();
Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t(0, 0, 1));
CoordinateSystemTrafo beamToReferenceCSTrafo(Vector_t(0, 0, pathLength_m), alignment.conjugate());
CoordinateSystemTrafo referenceToBeamCSTrafo = beamToReferenceCSTrafo.inverted();
const unsigned int localNum1 = itsBunch_m->getLocalNum();
for (unsigned int i = 0; i < localNum1; ++ i) {
itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
......@@ -1408,4 +1407,4 @@ void ParallelTTracker::evenlyDistributeParticles() {
if (requests.size() > 0) {
MPI_Waitall(requests.size(), &(requests[0]), MPI_STATUSES_IGNORE);
}
}
}
\ No newline at end of file
......@@ -60,7 +60,7 @@
#include "Beamlines/Beamline.h"
#include "Elements/OpalBeamline.h"
#include "Solvers/WakeFunction.hh"
#include "Solvers/WakeFunction.h"
#include <list>
#include <vector>
......
......@@ -20,7 +20,7 @@
#include "AbsBeamline/BeamlineVisitor.h"
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/ParticleMatterInteractionHandler.h"
#include "Structure/LossDataSink.h"
#include <cmath>
......@@ -75,7 +75,7 @@ bool CCollimator::doPreCheck(PartBunchBase<double, 3> *bunch) {
// rectangle collimators in cyclotron cylindrical coordinates
// when there is no particlematterinteraction, the particle hitting collimator is deleted directly
bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double /*t*/, const double /*tstep*/) {
bool flagNeedUpdate = false;
size_t tempnum = bunch->getLocalNum();
......@@ -88,7 +88,10 @@ bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumbe
/// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
if ((pflag != 0) && (bunch->Bin[i] != -1)) {
if (!parmatint_m)
lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
bunch->R[i], bunch->P[i],
bunch->getT(), bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
bunch->Bin[i] = -1;
flagNeedUpdate = true;
}
......
......@@ -402,8 +402,9 @@ bool Cyclotron::apply(const size_t& id, const double& t, Vector_t& E, Vector_t&
}
if (flagNeedUpdate) {
lossDs_m->addParticle(RefPartBunch_m->R[id], RefPartBunch_m->P[id],
id, t, 0, RefPartBunch_m->bunchNum[id]);
lossDs_m->addParticle(OpalParticle(id, RefPartBunch_m->R[id], RefPartBunch_m->P[id],
t, RefPartBunch_m->Q[id], RefPartBunch_m->M[id]),
std::make_pair(0, RefPartBunch_m->bunchNum[id]));
RefPartBunch_m->Bin[id] = -1;
}
......@@ -432,23 +433,23 @@ bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/,
// Necessary for gap phase output -DW
if (0 <= tet && tet <= 45) waiting_for_gap = 1;
// dB_{z}/dr, dB_{z}/dtheta, B_{z}
double brint = 0.0, btint = 0.0, bzint = 0.0;
if ( this->interpolate(rad, tet_rad, brint, btint, bzint) ) {
/* Br */
double br = - brint * R[2];
/* Btheta */
double bt = - btint / rad * R[2];
/* Bz */
double bz = - bzint;
this->applyTrimCoil(rad, R[2], tet_rad, br, bz);
/* Br Btheta -> Bx By */
B[0] = br * std::cos(tet_rad) - bt * std::sin(tet_rad);
B[1] = br * std::sin(tet_rad) + bt * std::cos(tet_rad);
......@@ -689,11 +690,11 @@ bool Cyclotron::interpolate(const double& rad,
const double wr1 = xir - (double)ir;
// wr2 : the relative distance to the outer path radius
const double wr2 = 1.0 - wr1;
// the corresponding angle on the field map
// Note: this does not work if the start point of field map does not equal zero.
double tet_map = std::fmod(tet_rad * Physics::rad2deg, 360.0 / symmetry_m);
double xit = tet_map / BP.dtet;
int it = (int) xit;
......@@ -1537,4 +1538,4 @@ void Cyclotron::getFieldFromFile_Synchrocyclotron(const double& scaleFactor) {
void Cyclotron::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const
{ }
#undef CHECK_CYC_FSCANF_EOF
#undef CHECK_CYC_FSCANF_EOF
\ No newline at end of file
......@@ -23,7 +23,7 @@
#include "AbsBeamline/BeamlineVisitor.h"
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/ParticleMatterInteractionHandler.h"
#include "Physics/Physics.h"
#include <memory>
#include <string>
......
......@@ -65,8 +65,8 @@
#include <string>
#include "Structure/BoundaryGeometry.h" // OPAL file
#include "Solvers/WakeFunction.hh"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/WakeFunction.h"
#include "Solvers/ParticleMatterInteractionHandler.h"
ElementBase::ElementBase():
......
......@@ -26,7 +26,7 @@
#include "Fields/Fieldmap.h"
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/ParticleMatterInteractionHandler.h"
#include "Utilities/Util.h"
#include <memory>
......@@ -113,9 +113,10 @@ bool FlexibleCollimator::apply(const size_t &i, const double &t, Vector_t &/*E*/
if (pdead) {
if (lossDs_m) {
double frac = -R(2) / P(2) * recpgamma;
lossDs_m->addParticle(R, P,
RefPartBunch_m->ID[i],
t + frac * dt, 0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R, P,
t + frac * dt,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i]));
}
++ losses_m;
}
......
......@@ -57,7 +57,7 @@ Monitor::Monitor(const std::string &name):
Component(name),
filename_m(""),
plane_m(OFF),
type_m(SPATIAL),
type_m(CollectionType::SPATIAL),
numPassages_m(0)
{}
......@@ -75,16 +75,17 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t
const Vector_t &P = RefPartBunch_m->P[i];
const double &dt = RefPartBunch_m->dt[i];
const Vector_t singleStep = Physics::c * dt * Util::getBeta(P);
if (online_m && type_m == SPATIAL) {
if (online_m && type_m == CollectionType::SPATIAL) {
if (dt * R(2) < 0.0 &&
dt * (R(2) + singleStep(2)) > 0.0) {
double frac = R(2) / singleStep(2);
lossDs_m->addParticle(R + frac * singleStep,
P,
RefPartBunch_m->ID[i],
t + frac * dt,
0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R + frac * singleStep,
P,
t + frac * dt,
RefPartBunch_m->Q[i],
RefPartBunch_m->M[i]));
}
}
......@@ -113,17 +114,18 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
RefPartBunch_m->get_sPos() + ds,
RefPartBunch_m->getGlobalTrackStep());
if (type_m == TEMPORAL) {
if (type_m == CollectionType::TEMPORAL) {
const unsigned int localNum = RefPartBunch_m->getLocalNum();
for (unsigned int i = 0; i < localNum; ++ i) {
Vector_t shift = ((frac - 0.5) * cdt * Util::getBeta(RefPartBunch_m->P[i])
- singleStep);
lossDs_m->addParticle(RefPartBunch_m->R[i] + shift,
RefPartBunch_m->P[i],
RefPartBunch_m->ID[i],
time,
0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
RefPartBunch_m->R[i] + shift,
RefPartBunch_m->P[i],
time,
RefPartBunch_m->Q[i],
RefPartBunch_m->M[i]));
}
OpalData::OPENMODE openMode;
if (numPassages_m > 0) {
......@@ -170,7 +172,7 @@ void Monitor::initialise(PartBunchBase<double, 3> *bunch, double &startField, do
}
}
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m, !Options::asciidump, getType()));
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m, !Options::asciidump, type_m));
}
void Monitor::finalise() {
......@@ -187,7 +189,7 @@ void Monitor::goOffline() {
statFileEntries_sm.insert(std::make_pair(stat.spos_m, stat));
}
if (type_m != TEMPORAL) {
if (type_m != CollectionType::TEMPORAL) {
lossDs_m->save(numPassages_m);
}
}
......
......@@ -54,11 +54,6 @@ public:
XY
};
enum Type {
TEMPORAL,
SPATIAL
};
/// Constructor with given name.
explicit Monitor(const std::string &name);
......@@ -102,7 +97,7 @@ public:
void setOutputFN(std::string fn);
void setType(Type type);
void setCollectionType(CollectionType type);
static void writeStatistics();
......@@ -115,7 +110,7 @@ private:
void operator=(const Monitor &);
std::string filename_m; /**< The name of the outputfile*/
Plane plane_m;
Type type_m;
CollectionType type_m;
unsigned int numPassages_m;
std::unique_ptr<LossDataSink> lossDs_m;
......@@ -125,7 +120,7 @@ private:
};
inline
void Monitor::setType(Monitor::Type type) {
void Monitor::setCollectionType(CollectionType type) {
type_m = type;
}
......
......@@ -108,8 +108,9 @@ bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const
// peak finder uses millimetre not metre
peakfinder_m->addParticle(probepoint * 1e3);
lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i], probepoint, bunch->P[i], t+dt,
bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
}
peakfinder_m->evaluate(turnnumber);
......
......@@ -106,7 +106,10 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
gmsgALL << getName() << ": particle " << id
<< " at " << refPartBunch_m->R[id]
<< " m out of the field map boundary" << endl;
lossDS_m->addParticle(refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id], id);
lossDS_m->addParticle(OpalParticle(id,
refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id],
t,
refPartBunch_m->Q[id], refPartBunch_m->M[id]));
lossDS_m->save();
refPartBunch_m->Bin[id] = -1;
......
......@@ -78,7 +78,7 @@ bool Septum::doPreCheck(PartBunchBase<double, 3> *bunch) {
return false;
}
bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double /*t*/, const double /*tstep*/) {
bool flag = false;
const double slope = (yend_m - ystart_m) / (xend_m - xstart_m);
......@@ -100,7 +100,11 @@ bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/,
R(1) > ystart_m &&
R(1) < yend_m) {
lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
R, bunch->P[i],
bunch->getT(),
bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
bunch->Bin[i] = -1;
flag = true;
}
......@@ -110,4 +114,4 @@ bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/,
ElementBase::ElementType Septum::getType() const {
return SEPTUM;
}
\ No newline at end of file
}
......@@ -54,8 +54,10 @@ bool Source::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t &
if (online_m && R(2) <= 0.0 && P(2) < 0.0) {
double frac = -R(2) / (P(2) * recpgamma);
lossDs_m->addParticle(R + frac * recpgamma * P,
P, RefPartBunch_m->ID[i], t + frac * dt, 0);
lossDs_m->addParticle(OpalParticle(RefPartBunch_m->ID[i],
R + frac * recpgamma * P, P,
t + frac * dt,
RefPartBunch_m->Q[i], RefPartBunch_m->M[i]));
return true;
}
......@@ -70,8 +72,7 @@ void Source::initialise(PartBunchBase<double, 3> *bunch, double &startField, dou
std::string filename = getName();
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename,
!Options::asciidump,
getType()));
!Options::asciidump));
}
......
......@@ -131,8 +131,11 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
strippoint(0) = (B_m * B_m * bunch->R[i](0) - A_m * B_m* bunch->R[i](1) - A_m * C_m) / (R_m * R_m);
strippoint(1) = (A_m * A_m * bunch->R[i](1) - A_m * B_m* bunch->R[i](0) - B_m * C_m) / (R_m * R_m);
strippoint(2) = bunch->R[i](2);
lossDs_m->addParticle(strippoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
lossDs_m->addParticle(OpalParticle(bunch->ID[i],
strippoint, bunch->P[i],
t+dt,
bunch->Q[i], bunch->M[i]),
std::make_pair(turnnumber, bunch->bunchNum[i]));
flagNeedUpdate = true;
if (stop_m) {
......
......@@ -27,8 +27,8 @@
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
#include "Physics/Physics.h"
#include "Solvers/BeamStrippingPhysics.hh"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Solvers/BeamStrippingPhysics.h"
#include "Solvers/ParticleMatterInteractionHandler.h"
#include "Utilities/LogicalError.h"
#include "Utilities/GeneralClassicException.h"
......
......@@ -3,6 +3,7 @@ set (_SRCS
AbstractTracker.cpp
CoordinateSystemTrafo.cpp
DefaultVisitor.cpp
DistributionMoments.cpp
Flagger.cpp
PartBunch.cpp
PartBins.cpp
......@@ -30,6 +31,7 @@ set (HDRS
AbstractTracker.h
CoordinateSystemTrafo.h
DefaultVisitor.h
DistributionMoments.h
Flagger.h
ListElem.h
PartBinsCyc.h
......
//
// Class DistributionMoments
// Computes the statistics of particle distributions.
//
// Copyright (c) 2021, Christof Metzger-Kraus
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "DistributionMoments.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "Message/GlobalComm.h"
#include "Utility/Inform.h"
#include "OpalParticle.h"
#include "PartBunchBase.h"
extern Inform* gmsg;
DistributionMoments::DistributionMoments()
{
reset();
}
void DistributionMoments::compute(PartBunchBase<double, 3> const& bunch)
{
computeStatistics(bunch.begin(), bunch.end());
}
void DistributionMoments::compute(const std::vector<OpalParticle>::const_iterator &first,
const std::vector<OpalParticle>::const_iterator &last)
{
computeStatistics(first, last);
}
/* 2 * Dim centroids + Dim * ( 2 * Dim + 1 ) 2nd moments + 2 * Dim (3rd and 4th order moments)
* --> 1st order moments: 0, ..., 2 * Dim - 1
* --> 2nd order moments: 2 * Dim, ..., Dim * ( 2 * Dim + 1 )
* --> 3rd order moments: Dim * ( 2 * Dim + 1 ) + 1, ..., Dim * ( 2 * Dim + 1 ) + Dim
* (only, <x^3>, <y^3> and <z^3>)
* --> 4th order moments: Dim * ( 2 * Dim + 1 ) + Dim + 1, ..., Dim * ( 2 * Dim + 1 ) + 2 * Dim
*
* For a 6x6 matrix we have each 2nd order moment (except diagonal
* entries) twice. We only store the upper half of the matrix.
*/
template<class InputIt>
void DistributionMoments::computeStatistics(const InputIt &first, const InputIt &last)
{
reset();
unsigned int localNum = last - first;
std::vector<double> localStatistics(41);
for (InputIt it = first; it != last; ++ it) {
OpalParticle const& particle = *it;
//FIXME After issue 287 is resolved this shouldn't be necessary anymore
if (!Options::amr
&& OpalData::getInstance()->isInOPALCyclMode()
&& particle.getId() == 0) {
-- localNum;
continue;
}
unsigned int l = 6;
for (unsigned int i = 0; i < 6; ++ i) {
localStatistics[i] += particle[i];
for (unsigned int j = 0; j <= i; ++ j, ++ l) {
localStatistics[l] += particle[i] * particle[j];
}
}
localStatistics[l++] += particle.getTime();
localStatistics[l++] += std::pow(particle.getTime(), 2);
for (unsigned int i = 0; i < 3; ++ i, l += 2) {
double r2 = std::pow(particle[i], 2);
localStatistics[l] += r2 * particle[i];
localStatistics[l + 1] += r2 * r2;
}
double gamma = Util::getGamma(particle.getP());
double eKin = (gamma - 1.0) * particle.getMass();
localStatistics[l++] += eKin;
localStatistics[l++] += std::pow(eKin, 2);
localStatistics[l++] += gamma;
localStatistics[l++] += particle.getCharge();
localStatistics[l++] += particle.getMass();
}
localStatistics.back() = localNum;
allreduce(localStatistics.data(), localStatistics.size(), std::plus<double>());
totalNumParticles_m = localStatistics.back();
fillMembers(localStatistics);
}
void DistributionMoments::computeMeanKineticEnergy(PartBunchBase<double, 3> const& bunch)
{
double data[] = {0.0, 0.0};
for (OpalParticle const& particle: bunch) {
data[0] += Util::getKineticEnergy(particle.getP(), particle.getMass());
}
data[1] = bunch.getLocalNum();
allreduce(data, 2, std::plus<double>());
meanKineticEnergy_m = data[0] / data[1];
}
void DistributionMoments::fillMembers(std::vector<double> const& localMoments) {
Vector_t squaredEps, fac, squaredSumR, squaredSumP, sumRP;
double perParticle = 1.0 / totalNumParticles_m;
unsigned int l = 0;
for (; l < 6; ++ l) {
centroid_m[l] = localMoments[l];
}
for (unsigned int i = 0; i < 6; ++ i) {
for (unsigned int j = 0; j <= i; ++ j, ++ l) {
moments_m(i, j) = localMoments[l];
moments_m(j, i) = moments_m(i, j);
}
}
meanTime_m = localMoments[l++] * perParticle;
stdTime_m = std::sqrt(localMoments[l++] * perParticle - std::pow(meanTime_m, 2));