OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
PartBins.cpp
Go to the documentation of this file.
1//
2// Class PartBins
3// Defines a structure to hold energy bins and their associated data
4//
5// Copyright (c) 2007-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
6// All rights reserved
7//
8// This file is part of OPAL.
9//
10// OPAL is free software: you can redistribute it and/or modify
11// it under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// You should have received a copy of the GNU General Public License
16// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
17//
18
19#include "Algorithms/PartBins.h"
21#include "Physics/Physics.h"
22#include "Physics/Units.h"
23
24#include "Utility/Inform.h"
25
26#include <cfloat>
27#include <limits>
28#include <vector>
29#include <cmath>
30
31PartBins::PartBins(int bins, int sbins) :
32 gamma_m(1.0),
33 bins_m(bins),
34 sBins_m(sbins),
35 xmin_m(0.0),
36 xmax_m(0.0),
37 hBin_m(0.0),
38 nemittedBins_m(0) {
39
40 // number of particles in the bins on the local node
41 nBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]);
42 xbinmin_m = std::unique_ptr<double[]>(new double[bins_m]);
43 xbinmax_m = std::unique_ptr<double[]>(new double[bins_m]);
44
45 // flag whether the bin contain particles or not
46 binsEmitted_m = std::unique_ptr<bool[]>(new bool[bins_m]);
47
48 nDelBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]);
49
50 for(int i = 0; i < bins_m; i++) {
51 nDelBin_m[i] = nBin_m[i] = 0;
53 xbinmax_m[i] = -xbinmin_m[i];
54 binsEmitted_m[i] = false;
55 }
56}
57
58
60 size_t s = 0;
61 size_t sd = 0;
62 size_t gs = 0;
63
64 for(int i = 0; i < getLastemittedBin(); i++) {
65 s += nBin_m[i];
66 sd += nDelBin_m[i];
67 }
68 gs = s - sd;
69 reduce(gs, gs, OpAddAssign());
70 return gs;
71}
72
74 size_t s = 0;
75 s = nBin_m[b];
76 reduce(s, s, OpAddAssign());
77 return s;
78}
79
80
81void PartBins::updatePartInBin_cyc(size_t countLost[]) {
82
83 for(int ii = 0; ii < nemittedBins_m; ii++) {
84 if(countLost[ii] > 0)
85 nBin_m[ii] -= countLost[ii];
86 }
87}
88
89void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
90 reduce(maxbinIndex, maxbinIndex, OpMaxAssign());
91 nemittedBins_m = maxbinIndex + 1;
92
93 for(int ii = 0; ii < nemittedBins_m; ii++) {
94 nBin_m[ii] = newPartNum[ii]; // only count particles on the local node
95 setBinEmitted(ii); // set true for this bin
96 }
97}
98
99
101 tmppart_m.clear();
102 isEmitted_m.clear();
103}
104
105
106bool PartBins::getPart(size_t n, int bin, std::vector<double> &p) {
107
108 if(tmppart_m[n][6] == bin) {
109 p = tmppart_m[n];
110 return true;
111 } else
112 return false;
113}
114
122 double sshift = std::sqrt(1. - (1. / (gamma_m * gamma_m))) * Physics::c * 0.1 * Units::ps2s;
123 std::sort(tmppart_m.begin(), tmppart_m.end(), DescendingLocationSort(2));
124 xmax_m = tmppart_m[0][2];
125 xmin_m = tmppart_m.back()[2];
126
127 for(unsigned int n = 0; n < tmppart_m.size(); n++)
128 tmppart_m[n][2] -= xmax_m + sshift; /* push particles back */
129
130 xmin_m -= xmax_m + 0.0001 * (xmax_m - xmin_m) + sshift; /* lower the limits */
131 xmax_m = -sshift;
132
135
137 calcHBins();
138 for(int n = 0; n < bins_m; n++)
139 if(nBin_m[n] == 0) setBinEmitted(n);
140}
141
142
144
145 for(unsigned int n = 0; n < tmppart_m.size(); n++)
146 tmppart_m[n][6] = getBin(tmppart_m[n][2]);
147 calcExtrema();
148}
149
151 for(unsigned int n = 0; n < tmppart_m.size(); n++) {
152 if(xbinmin_m[(int)tmppart_m[n][6]] >= tmppart_m[n][2])
153 xbinmin_m[(int)tmppart_m[n][6]] = tmppart_m[n][2];
154
155 if(xbinmax_m[(int)tmppart_m[n][6]] <= tmppart_m[n][2])
156 xbinmax_m[(int)tmppart_m[n][6]] = tmppart_m[n][2];
157 }
158}
159
161
162 os << "-----------------------------------------" << endl;
163 os << " CREATE BINNED GAUSS DISTRIBUTION DONE " << endl;
164
165 os << "Bins= " << bins_m << " hBin= " << hBin_m << " Particle vector length " << tmppart_m.size() << endl;
166
167 //for(int i = 0; i < gsl_histogram_bins(h_m); i++)
168 //os << "Bin # " << i << " val " << gsl_histogram_get(h_m, i) << endl;
169 for(int i = 0; i < bins_m; i++) {
170 size_t msum = 0;
171 for(int j = 0; j < sBins_m; j++)
172 msum += gsl_histogram_get(h_m.get(), i * sBins_m + j);
173 os << "Bin # " << i << " val " << msum << endl;
174 }
175
176 if(getLastemittedBin() >= 0)
177 os << "Last emitted bin is " << getLastemittedBin() << endl;
178 else
179 os << "No bin is emitted !" << endl;
180 return os;
181}
182
183int PartBins::getBin(double x) {
188 int b = (int) std::floor(std::abs(xmax_m - x) / hBin_m);
189 nBin_m[b]++;
190 return b;
191}
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
Definition: GlobalComm.hpp:55
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:733
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
constexpr double ps2s
Definition: Units.h:53
virtual ~PartBins()
Definition: PartBins.cpp:100
void calcHBins()
Definition: PartBins.cpp:143
std::vector< std::vector< double > > tmppart_m
Definition: PartBins.h:121
int sBins_m
Definition: PartBins.h:107
std::unique_ptr< double[]> xbinmax_m
Definition: PartBins.h:115
std::unique_ptr< double[]> xbinmin_m
Definition: PartBins.h:114
double xmin_m
Definition: PartBins.h:110
double gamma_m
Definition: PartBins.h:99
PartBins(int bins, int sbins)
Definition: PartBins.cpp:31
std::unique_ptr< size_t[]> nDelBin_m
Definition: PartBins.h:163
size_t getTotalNumPerBin(int b)
How many particles are in the bin b.
Definition: PartBins.cpp:73
void resetPartInBin_cyc(size_t newPartNum[], int binID)
Definition: PartBins.cpp:89
std::vector< bool > isEmitted_m
Definition: PartBins.h:122
void setBinEmitted(int bin)
Definition: PartBins.h:86
void calcExtrema()
Definition: PartBins.cpp:150
int getBin(double x)
Definition: PartBins.cpp:183
void updatePartInBin_cyc(size_t countLost[])
Definition: PartBins.cpp:81
void sortArray()
Definition: PartBins.cpp:116
Inform & print(Inform &os)
Definition: PartBins.cpp:160
std::unique_ptr< bool[]> binsEmitted_m
Definition: PartBins.h:124
double hBin_m
Definition: PartBins.h:118
bool getPart(size_t n, int bin, std::vector< double > &p)
Definition: PartBins.cpp:106
int nemittedBins_m
Definition: PartBins.h:157
size_t getTotalNum()
How many particles are in all the bins.
Definition: PartBins.cpp:59
double xmax_m
Definition: PartBins.h:111
int getLastemittedBin()
Definition: PartBins.h:136
int bins_m
Definition: PartBins.h:106
std::unique_ptr< size_t[]> nBin_m
Definition: PartBins.h:160
std::unique_ptr< gsl_histogram > h_m
Definition: PartBins.h:165
Definition: Inform.h:42