OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
Distribution.cpp
Go to the documentation of this file.
1 // ------------------------------------------------------------------------
2 // $RCSfile: Distribution.cpp,v $
3 // ------------------------------------------------------------------------
4 // $Revision: 1.3.4.1 $
5 // ------------------------------------------------------------------------
6 // Copyright: see Copyright.readme
7 // ------------------------------------------------------------------------
8 //
9 // Class: Distribution
10 // The class for the OPAL Distribution command.
11 //
12 // ------------------------------------------------------------------------
13 
14 #include "Distribution/Distribution.h"
15 #include "Distribution/SigmaGenerator.h"
16 #include "AbsBeamline/SpecificElementVisitor.h"
17 
18 #include <cmath>
19 #include <cfloat>
20 #include <iomanip>
21 #include <iostream>
22 #include <string>
23 #include <vector>
24 #include <numeric>
25 #include <limits>
26 
27 #include "AbstractObjects/Expressions.h"
28 #include "Utilities/Options.h"
29 #include "AbstractObjects/OpalData.h"
30 #include "Algorithms/PartBins.h"
31 #include "Algorithms/PartBunchBase.h"
32 #include "Algorithms/bet/EnvelopeBunch.h"
33 #include "Structure/Beam.h"
34 #include "Structure/BoundaryGeometry.h"
35 #include "Algorithms/PartBinsCyc.h"
36 #include "BasicActions/Option.h"
37 #include "Distribution/LaserProfile.h"
38 #include "Elements/OpalBeamline.h"
39 #include "AbstractObjects/BeamSequence.h"
40 #include "Structure/H5PartWrapper.h"
41 #include "Structure/H5PartWrapperForPC.h"
42 #include "Utilities/Util.h"
43 #include "Utilities/EarlyLeaveException.h"
44 
45 #include <gsl/gsl_cdf.h>
46 #include <gsl/gsl_randist.h>
47 #include <gsl/gsl_sf_erf.h>
48 #include <gsl/gsl_linalg.h>
49 #include <gsl/gsl_blas.h>
50 
51 #include <sys/time.h>
52 
53 #include <boost/regex.hpp>
54 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
55 
56 extern Inform *gmsg;
57 
58 constexpr double SMALLESTCUTOFF = 1e-12;
59 
60 namespace {
61  SymTenzor<double, 6> getUnit6x6() {
62  SymTenzor<double, 6> unit6x6;
63  for (unsigned int i = 0; i < 6u; ++ i) {
64  unit6x6(i,i) = 1.0;
65  }
66  return unit6x6;
67  }
68 }
69 
70 //
71 // Class Distribution
72 // ------------------------------------------------------------------------
73 
74 Distribution::Distribution():
75  Definition( Attrib::Legacy::Distribution::SIZE, "DISTRIBUTION",
76  "The DISTRIBUTION statement defines data for the 6D particle distribution."),
77  distrTypeT_m(DistrTypeT::NODIST),
78  numberOfDistributions_m(1),
79  emitting_m(false),
80  emissionModel_m(EmissionModelT::NONE),
81  tEmission_m(0.0),
82  tBin_m(0.0),
83  currentEmissionTime_m(0.0),
84  currentEnergyBin_m(1),
85  currentSampleBin_m(0),
86  numberOfEnergyBins_m(0),
87  numberOfSampleBins_m(0),
88  energyBins_m(NULL),
89  energyBinHist_m(NULL),
90  randGen_m(NULL),
91  pTotThermal_m(0.0),
92  pmean_m(0.0),
93  cathodeWorkFunc_m(0.0),
94  laserEnergy_m(0.0),
95  cathodeFermiEnergy_m(0.0),
96  cathodeTemp_m(0.0),
97  emitEnergyUpperLimit_m(0.0),
98  totalNumberParticles_m(0),
99  totalNumberEmittedParticles_m(0),
100  avrgpz_m(0.0),
101  inputMoUnits_m(InputMomentumUnitsT::NONE),
102  sigmaTRise_m(0.0),
103  sigmaTFall_m(0.0),
104  tPulseLengthFWHM_m(0.0),
105  correlationMatrix_m(getUnit6x6()),
106  laserProfileFileName_m(""),
107  laserImageName_m(""),
108  laserIntensityCut_m(0.0),
109  laserProfile_m(NULL),
110  darkCurrentParts_m(0),
111  darkInwardMargin_m(0.0),
112  eInitThreshold_m(0.0),
113  workFunction_m(0.0),
114  fieldEnhancement_m(0.0),
115  fieldThrFN_m(0.0),
116  maxFN_m(0),
117  paraFNA_m(0.0),
118  paraFNB_m(0.0),
119  paraFNY_m(0.0),
120  paraFNVYSe_m(0.0),
121  paraFNVYZe_m(0.0),
122  secondaryFlag_m(0),
123  ppVw_m(0.0),
124  vVThermal_m(0.0),
125  I_m(0.0),
126  E_m(0.0)
127 {
128  setAttributes();
129 
130  Distribution *defaultDistribution = clone("UNNAMED_Distribution");
131  defaultDistribution->builtin = true;
132 
133  try {
134  OpalData::getInstance()->define(defaultDistribution);
135  } catch(...) {
136  delete defaultDistribution;
137  }
138 
139  // setFieldEmissionParameters();
140 
141  gsl_rng_env_setup();
142  randGen_m = gsl_rng_alloc(gsl_rng_default);
143 }
150 Distribution::Distribution(const std::string &name, Distribution *parent):
151  Definition(name, parent),
152  distT_m(parent->distT_m),
153  distrTypeT_m(DistrTypeT::NODIST),
154  numberOfDistributions_m(parent->numberOfDistributions_m),
155  emitting_m(parent->emitting_m),
156  particleRefData_m(parent->particleRefData_m),
157  addedDistributions_m(parent->addedDistributions_m),
158  particlesPerDist_m(parent->particlesPerDist_m),
159  emissionModel_m(parent->emissionModel_m),
160  tEmission_m(parent->tEmission_m),
161  tBin_m(parent->tBin_m),
162  currentEmissionTime_m(parent->currentEmissionTime_m),
163  currentEnergyBin_m(parent->currentEmissionTime_m),
164  currentSampleBin_m(parent->currentSampleBin_m),
165  numberOfEnergyBins_m(parent->numberOfEnergyBins_m),
166  numberOfSampleBins_m(parent->numberOfSampleBins_m),
167  energyBins_m(NULL),
168  energyBinHist_m(NULL),
169  randGen_m(NULL),
170  pTotThermal_m(parent->pTotThermal_m),
171  pmean_m(parent->pmean_m),
172  cathodeWorkFunc_m(parent->cathodeWorkFunc_m),
173  laserEnergy_m(parent->laserEnergy_m),
174  cathodeFermiEnergy_m(parent->cathodeFermiEnergy_m),
175  cathodeTemp_m(parent->cathodeTemp_m),
176  emitEnergyUpperLimit_m(parent->emitEnergyUpperLimit_m),
177  totalNumberParticles_m(parent->totalNumberParticles_m),
178  totalNumberEmittedParticles_m(parent->totalNumberEmittedParticles_m),
179  xDist_m(parent->xDist_m),
180  pxDist_m(parent->pxDist_m),
181  yDist_m(parent->yDist_m),
182  pyDist_m(parent->pyDist_m),
183  tOrZDist_m(parent->tOrZDist_m),
184  pzDist_m(parent->pzDist_m),
185  xWrite_m(parent->xWrite_m),
186  pxWrite_m(parent->pxWrite_m),
187  yWrite_m(parent->yWrite_m),
188  pyWrite_m(parent->pyWrite_m),
189  tOrZWrite_m(parent->tOrZWrite_m),
190  pzWrite_m(parent->pzWrite_m),
191  avrgpz_m(parent->avrgpz_m),
192  inputMoUnits_m(parent->inputMoUnits_m),
193  sigmaTRise_m(parent->sigmaTRise_m),
194  sigmaTFall_m(parent->sigmaTFall_m),
195  tPulseLengthFWHM_m(parent->tPulseLengthFWHM_m),
196  sigmaR_m(parent->sigmaR_m),
197  sigmaP_m(parent->sigmaP_m),
198  cutoffR_m(parent->cutoffR_m),
199  cutoffP_m(parent->cutoffP_m),
200  correlationMatrix_m(parent->correlationMatrix_m),
201  laserProfileFileName_m(parent->laserProfileFileName_m),
202  laserImageName_m(parent->laserImageName_m),
203  laserIntensityCut_m(parent->laserIntensityCut_m),
204  laserProfile_m(NULL),
205  darkCurrentParts_m(parent->darkCurrentParts_m),
206  darkInwardMargin_m(parent->darkInwardMargin_m),
207  eInitThreshold_m(parent->eInitThreshold_m),
208  workFunction_m(parent->workFunction_m),
209  fieldEnhancement_m(parent->fieldEnhancement_m),
210  fieldThrFN_m(parent->fieldThrFN_m),
211  maxFN_m(parent->maxFN_m),
212  paraFNA_m(parent-> paraFNA_m),
213  paraFNB_m(parent-> paraFNB_m),
214  paraFNY_m(parent-> paraFNY_m),
215  paraFNVYSe_m(parent-> paraFNVYSe_m),
216  paraFNVYZe_m(parent-> paraFNVYZe_m),
217  secondaryFlag_m(parent->secondaryFlag_m),
218  ppVw_m(parent->ppVw_m),
219  vVThermal_m(parent->vVThermal_m),
220  I_m(parent->I_m),
221  E_m(parent->E_m),
222  tRise_m(parent->tRise_m),
223  tFall_m(parent->tFall_m),
224  sigmaRise_m(parent->sigmaRise_m),
225  sigmaFall_m(parent->sigmaFall_m),
226  cutoff_m(parent->cutoff_m)
227 {
228  gsl_rng_env_setup();
229  randGen_m = gsl_rng_alloc(gsl_rng_default);
230 }
231 
232 Distribution::~Distribution() {
233 
234  delete energyBins_m;
235  gsl_histogram_free(energyBinHist_m);
236  gsl_rng_free(randGen_m);
237  delete laserProfile_m;
238 }
239 
240 /*
241  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
242  for (int i=0; i<M.size1(); ++i) {
243  for (int j=0; j<M.size2(); ++j) {
244  *gmsg << M(i,j) << " ";
245  }
246  *gmsg << endl;
247  }
248  }
249 */
250 
259 size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {
260 
261  size_t locNumber = n / Ippl::getNodes();
262 
263  // make sure the total number is exact
264  size_t remainder = n % Ippl::getNodes();
265  size_t myNode = Ippl::myNode();
266  if (myNode < remainder)
267  ++ locNumber;
268 
269  return locNumber;
270 }
271 
273 bool Distribution::canReplaceBy(Object *object) {
274  return dynamic_cast<Distribution *>(object) != 0;
275 }
276 
277 Distribution *Distribution::clone(const std::string &name) {
278  return new Distribution(name, this);
279 }
280 
281 void Distribution::execute() {
282 }
283 
284 void Distribution::update() {
285 }
286 
287 void Distribution::create(size_t &numberOfParticles, double massIneV) {
288 
289  /*
290  * If Options::cZero is true than we reflect generated distribution
291  * to ensure that the transverse averages are 0.0.
292  *
293  * For now we just cut the number of generated particles in half.
294  */
295  size_t numberOfLocalParticles = numberOfParticles;
296  if (Options::cZero && distrTypeT_m != DistrTypeT::FROMFILE) {
297  numberOfLocalParticles = (numberOfParticles + 1) / 2;
298  }
299 
300  size_t mySeed = Options::seed;
301 
302  if (Options::seed == -1) {
303  struct timeval tv;
304  gettimeofday(&tv,0);
305  mySeed = tv.tv_sec + tv.tv_usec;
306  }
307 
308  if (Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE])) {
309  numberOfLocalParticles = getNumOfLocalParticlesToCreate(numberOfLocalParticles);
310  *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << " + core_id\n"
311  << "* is scalable with number of particles and cores." << endl;
312  mySeed += Ippl::myNode();
313  } else {
314  *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << "\n"
315  << "* isn't scalable with number of particles and cores." << endl;
316  }
317 
318  gsl_rng_set(randGen_m, mySeed);
319 
320  setFieldEmissionParameters();
321 
322  switch (distrTypeT_m) {
323 
324  case DistrTypeT::MATCHEDGAUSS:
325  createMatchedGaussDistribution(numberOfLocalParticles, massIneV);
326  break;
327  case DistrTypeT::FROMFILE:
328  createDistributionFromFile(numberOfParticles, massIneV);
329  break;
330  case DistrTypeT::GAUSS:
331  createDistributionGauss(numberOfLocalParticles, massIneV);
332  break;
333  case DistrTypeT::BINOMIAL:
334  createDistributionBinomial(numberOfLocalParticles, massIneV);
335  break;
336  case DistrTypeT::FLATTOP:
337  case DistrTypeT::GUNGAUSSFLATTOPTH:
338  case DistrTypeT::ASTRAFLATTOPTH:
339  createDistributionFlattop(numberOfLocalParticles, massIneV);
340  break;
341  default:
342  INFOMSG("Distribution unknown." << endl;);
343  break;
344  }
345 
346  if (emitting_m) {
347 
348  unsigned int numAdditionalRNsPerParticle = 0;
349  if (emissionModel_m == EmissionModelT::ASTRA ||
350  distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
351  distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
352 
353  numAdditionalRNsPerParticle = 2;
354  } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
355  if (Options::cZero) {
356  numAdditionalRNsPerParticle = 40;
357  } else {
358  numAdditionalRNsPerParticle = 20;
359  }
360  }
361 
362  int saveProcessor = -1;
363  const int myNode = Ippl::myNode();
364  const int numNodes = Ippl::getNodes();
365  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
366 
367  for (size_t partIndex = 0; partIndex < numberOfLocalParticles; ++ partIndex) {
368 
369  // Save to each processor in turn.
370  ++ saveProcessor;
371  if (saveProcessor >= numNodes)
372  saveProcessor = 0;
373 
374  if (scalable || myNode == saveProcessor) {
375  std::vector<double> rns;
376  for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
377  double x = gsl_rng_uniform(randGen_m);
378  rns.push_back(x);
379  }
380  additionalRNs_m.push_back(rns);
381  } else {
382  for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
383  gsl_rng_uniform(randGen_m);
384  }
385  }
386  }
387  }
388 
389  // Scale coordinates according to distribution input.
390  scaleDistCoordinates();
391 
392  // Now reflect particles if Options::cZero is true
393  reflectDistribution(numberOfLocalParticles);
394 
395  adjustPhaseSpace(massIneV);
396 
397  if (Options::seed != -1)
398  Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
399 
400  if (particlesPerDist_m.size() == 0) {
401  particlesPerDist_m.push_back(tOrZDist_m.size());
402  } else {
403  particlesPerDist_m[0] = tOrZDist_m.size();
404  }
405 }
406 
407 void Distribution::createPriPart(PartBunchBase<double, 3> *beam, BoundaryGeometry &bg) {
408 
409  if (Options::ppdebug) { // This is Parallel Plate Benchmark.
410  int pc = 0;
411  size_t lowMark = beam->getLocalNum();
412  double vw = this->getVw();
413  double vt = this->getvVThermal();
414  double f_max = vw / vt * exp(-0.5);
415  double test_a = vt / vw;
416  double test_asq = test_a * test_a;
417  size_t count = 0;
418  size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
419  size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
420  if (Ippl::myNode() == 0)
421  N_mean += N_extra;
422  if (bg.getN() != 0) {
423  for (size_t i = 0; i < bg.getN(); i++) {
424  if (pc == Ippl::myNode()) {
425  if (count < N_mean) {
426  /*==============Parallel Plate Benchmark=====================================*/
427  double test_s = 1;
428  double f_x = 0;
429  double test_x = 0;
430  while(test_s > f_x) {
431  test_s = gsl_rng_uniform(randGen_m);
432  test_s *= f_max;
433  test_x = gsl_rng_uniform(randGen_m);
434  test_x *= 10 * test_a; //range for normalized emission speed(0,10*test_a);
435  f_x = test_x / test_asq * exp(-test_x * test_x / 2 / test_asq);
436  }
437  double v_emi = test_x * vw;
438 
439  double betaemit = v_emi / Physics::c;
440  double betagamma = betaemit / sqrt(1 - betaemit * betaemit);
441  /*============================================================================ */
442  beam->create(1);
443  if (pc != 0) {
444  beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra);
445  beam->P[lowMark + count] = betagamma * bg.getMomenta(Ippl::myNode() * N_mean + count);
446  } else {
447  beam->R[lowMark + count] = bg.getCooridinate(count);
448  beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
449  }
450  beam->Bin[lowMark + count] = 0;
451  beam->PType[lowMark + count] = ParticleType::REGULAR;
452  beam->TriID[lowMark + count] = 0;
453  beam->Q[lowMark + count] = beam->getChargePerParticle();
454  beam->Ef[lowMark + count] = Vector_t(0.0);
455  beam->Bf[lowMark + count] = Vector_t(0.0);
456  beam->dt[lowMark + count] = beam->getdT();
457  count ++;
458  }
459  }
460  pc++;
461  if (pc == Ippl::getNodes())
462  pc = 0;
463  }
464  bg.clearCooridinateArray();
465  bg.clearMomentaArray();
466  beam->boundp();
467  }
468  *gmsg << *beam << endl;
469 
470  } else {// Normal procedure to create primary particles
471 
472  int pc = 0;
473  size_t lowMark = beam->getLocalNum();
474  size_t count = 0;
475  size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
476  size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
477 
478  if (Ippl::myNode() == 0)
479  N_mean += N_extra;
480  if (bg.getN() != 0) {
481  for (size_t i = 0; i < bg.getN(); i++) {
482 
483  if (pc == Ippl::myNode()) {
484  if (count < N_mean) {
485  beam->create(1);
486  if (pc != 0)
487  beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra); // node 0 will emit the particle with coordinate ID from 0 to N_mean+N_extra, so other nodes should shift to node_number*N_mean+N_extra
488  else
489  beam->R[lowMark + count] = bg.getCooridinate(count); // for node0 the particle number N_mean = N_mean + N_extra
490  beam->P[lowMark + count] = Vector_t(0.0);
491  beam->Bin[lowMark + count] = 0;
492  beam->PType[lowMark + count] = ParticleType::REGULAR;
493  beam->TriID[lowMark + count] = 0;
494  beam->Q[lowMark + count] = beam->getChargePerParticle();
495  beam->Ef[lowMark + count] = Vector_t(0.0);
496  beam->Bf[lowMark + count] = Vector_t(0.0);
497  beam->dt[lowMark + count] = beam->getdT();
498  count++;
499 
500  }
501 
502  }
503  pc++;
504  if (pc == Ippl::getNodes())
505  pc = 0;
506 
507  }
508 
509  }
510  bg.clearCooridinateArray();
511  beam->boundp();//fixme if bg.getN()==0?
512  }
513  *gmsg << *beam << endl;
514 }
515 
516 void Distribution::doRestartOpalT(PartBunchBase<double, 3> *beam, size_t Np, int restartStep, H5PartWrapper *dataSource) {
517 
518  IpplTimings::startTimer(beam->distrReload_m);
519 
520  long numParticles = dataSource->getNumParticles();
521  size_t numParticlesPerNode = numParticles / Ippl::getNodes();
522 
523  size_t firstParticle = numParticlesPerNode * Ippl::myNode();
524  size_t lastParticle = firstParticle + numParticlesPerNode - 1;
525  if (Ippl::myNode() == Ippl::getNodes() - 1)
526  lastParticle = numParticles - 1;
527  OpalData::getInstance()->addProblemCharacteristicValue("NP", numParticles);
528 
529  numParticles = lastParticle - firstParticle + 1;
530  PAssert_GE(numParticles, 0);
531 
532  beam->create(numParticles);
533 
534  dataSource->readHeader();
535  dataSource->readStep(beam, firstParticle, lastParticle);
536 
537  beam->boundp();
538 
539  double actualT = beam->getT();
540  long long ltstep = beam->getLocalTrackStep();
541  long long gtstep = beam->getGlobalTrackStep();
542 
543  IpplTimings::stopTimer(beam->distrReload_m);
544 
545  *gmsg << "Total number of particles in the h5 file= " << beam->getTotalNum() << "\n"
546  << "Global step= " << gtstep << "; Local step= " << ltstep << "; "
547  << "restart step= " << restartStep << "\n"
548  << "time of restart= " << actualT << "; phishift= " << OpalData::getInstance()->getGlobalPhaseShift() << endl;
549 }
550 
551 void Distribution::doRestartOpalCycl(PartBunchBase<double, 3> *beam,
552  size_t Np,
553  int restartStep,
554  const int specifiedNumBunch,
555  H5PartWrapper *dataSource) {
556 
557  // h5_int64_t rc;
558  IpplTimings::startTimer(beam->distrReload_m);
559  INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
560 
561  long numParticles = dataSource->getNumParticles();
562  size_t numParticlesPerNode = numParticles / Ippl::getNodes();
563 
564  size_t firstParticle = numParticlesPerNode * Ippl::myNode();
565  size_t lastParticle = firstParticle + numParticlesPerNode - 1;
566  if (Ippl::myNode() == Ippl::getNodes() - 1)
567  lastParticle = numParticles - 1;
568  OpalData::getInstance()->addProblemCharacteristicValue("NP", numParticles);
569 
570  numParticles = lastParticle - firstParticle + 1;
571  PAssert_GE(numParticles, 0);
572 
573  beam->create(numParticles);
574 
575  dataSource->readHeader();
576  dataSource->readStep(beam, firstParticle, lastParticle);
577 
578  beam->Q = beam->getChargePerParticle();
579 
580  beam->boundp();
581 
582  double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
583 
584  const int globalN = beam->getTotalNum();
585  INFOMSG("Restart from hdf5 format file " << OpalData::getInstance()->getRestartFileName() << endl);
586  INFOMSG("total number of particles = " << globalN << endl);
587  INFOMSG("* Restart Energy = " << meanE << " (MeV), Path lenght = " << beam->get_sPos() << " (m)" << endl);
588  INFOMSG("Tracking Step since last bunch injection is " << beam->getSteptoLastInj() << endl);
589  INFOMSG(beam->getNumBunch() << " Bunches(bins) exist in this file" << endl);
590 
591  double gamma = 1 + meanE / beam->getM() * 1.0e6;
592  double beta = sqrt(1.0 - (1.0 / std::pow(gamma, 2.0)));
593 
594  INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
595 
596  if (dataSource->predecessorIsSameFlavour()) {
597  INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
598  if (specifiedNumBunch > 1) {
599  // the allowed maximal bin number is set to 1000
600  energyBins_m = new PartBinsCyc(1000, beam->getNumBunch());
601  beam->setPBins(energyBins_m);
602  }
603 
604  } else {
605  INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
606 
607  Vector_t meanR(0.0, 0.0, 0.0);
608  Vector_t meanP(0.0, 0.0, 0.0);
609  unsigned long int newLocalN = beam->getLocalNum();
610  for (unsigned int i = 0; i < newLocalN; ++i) {
611  for (int d = 0; d < 3; ++d) {
612  meanR(d) += beam->R[i](d);
613  meanP(d) += beam->P[i](d);
614  }
615  }
616  reduce(meanR, meanR, OpAddAssign());
617  meanR /= Vector_t(globalN);
618  reduce(meanP, meanP, OpAddAssign());
619  meanP /= Vector_t(globalN);
620  INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
621 
622  for (unsigned int i = 0; i < newLocalN; ++i) {
623  beam->R[i] -= meanR;
624  beam->P[i] -= meanP;
625  }
626  }
627 
628  INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
629  IpplTimings::stopTimer(beam->distrReload_m);
630 }
631 
632 void Distribution::doRestartOpalE(EnvelopeBunch *beam, size_t Np, int restartStep,
633  H5PartWrapper *dataSource) {
634  IpplTimings::startTimer(beam->distrReload_m);
635  int N = dataSource->getNumParticles();
636  *gmsg << "total number of slices = " << N << endl;
637 
638  beam->distributeSlices(N);
639  beam->createBunch();
640  long long starti = beam->mySliceStartOffset();
641  long long endi = beam->mySliceEndOffset();
642 
643  dataSource->readHeader();
644  dataSource->readStep(beam, starti, endi);
645 
646  beam->setCharge(beam->getChargePerParticle());
647  IpplTimings::stopTimer(beam->distrReload_m);
648 }
649 
650 Distribution *Distribution::find(const std::string &name) {
651  Distribution *dist = dynamic_cast<Distribution *>(OpalData::getInstance()->find(name));
652 
653  if (dist == 0) {
654  throw OpalException("Distribution::find()", "Distribution \"" + name + "\" not found.");
655  }
656 
657  return dist;
658 }
659 
660 double Distribution::getTEmission() {
661  if (tEmission_m > 0.0) {
662  return tEmission_m;
663  }
664 
665  setDistType();
666 
667  tPulseLengthFWHM_m = Attributes::getReal(itsAttr[Attrib::Distribution::TPULSEFWHM]);
668  cutoff_m = Attributes::getReal(itsAttr[Attrib::Legacy::Distribution::CUTOFF]);
669  tRise_m = Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]);
670  tFall_m = Attributes::getReal(itsAttr[Attrib::Distribution::TFALL]);
671  double tratio = sqrt(2.0 * log(10.0)) - sqrt(2.0 * log(10.0 / 9.0));
672  sigmaRise_m = tRise_m / tratio;
673  sigmaFall_m = tFall_m / tratio;
674 
675  switch(distrTypeT_m) {
676  case DistrTypeT::ASTRAFLATTOPTH: {
677  double a = tPulseLengthFWHM_m / 2;
678  double sig = tRise_m / 2;
679  double inv_erf08 = 0.906193802436823; // erfinv(0.8)
680  double sqr2 = sqrt(2.);
681  double t = a - sqr2 * sig * inv_erf08;
682  double tmps = sig;
683  double tmpt = t;
684  for (int i = 0; i < 10; ++ i) {
685  sig = (t + tRise_m - a) / (sqr2 * inv_erf08);
686  t = a - 0.5 * sqr2 * (sig + tmps) * inv_erf08;
687  sig = (0.5 * (t + tmpt) + tRise_m - a) / (sqr2 * inv_erf08);
688  tmps = sig;
689  tmpt = t;
690  }
691  tEmission_m = tPulseLengthFWHM_m + 10 * sig;
692  break;
693  }
694  case DistrTypeT::GUNGAUSSFLATTOPTH: {
695  tEmission_m = tPulseLengthFWHM_m + (cutoff_m - sqrt(2.0 * log(2.0))) * (sigmaRise_m + sigmaFall_m);
696  break;
697  }
698  default:
699  tEmission_m = 0.0;
700  }
701  return tEmission_m;
702 }
703 
704 Inform &Distribution::printInfo(Inform &os) const {
705 
706  os << "\n"
707  << "* ************* D I S T R I B U T I O N ********************************************" << endl;
708  os << "* " << endl;
709  if (OpalData::getInstance()->inRestartRun()) {
710  os << "* In restart. Distribution read in from .h5 file." << endl;
711  } else {
712  if (addedDistributions_m.size() > 0)
713  os << "* Main Distribution" << endl
714  << "-----------------" << endl;
715 
716  if (particlesPerDist_m.empty())
717  printDist(os, 0);
718  else
719  printDist(os, particlesPerDist_m.at(0));
720 
721  size_t distCount = 1;
722  for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
723  os << "* " << endl;
724  os << "* Added Distribution #" << distCount << endl;
725  os << "* ----------------------" << endl;
726  addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
727  distCount++;
728  }
729 
730  os << "* " << endl;
731  if (numberOfEnergyBins_m > 0) {
732  os << "* Number of energy bins = " << numberOfEnergyBins_m << endl;
733 
734  // if (numberOfEnergyBins_m > 1)
735  // printEnergyBins(os);
736  }
737 
738  if (emitting_m) {
739  os << "* Distribution is emitted. " << endl;
740  os << "* Emission time = " << tEmission_m << " [sec]" << endl;
741  os << "* Time per bin = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
742  os << "* Delta t during emission = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
743  os << "* " << endl;
744  printEmissionModel(os);
745  } else {
746  os << "* Distribution is injected." << endl;
747  }
748  }
749  os << "* " << endl;
750  os << "* **********************************************************************************" << endl;
751 
752  return os;
753 }
754 
755 const PartData &Distribution::getReference() const {
756  // Cast away const, to allow logically constant Distribution to update.
757  const_cast<Distribution *>(this)->update();
758  return particleRefData_m;
759 }
760 
761 void Distribution::addDistributions() {
762  /*
763  * Move particle coordinates from added distributions to main distribution.
764  */
765 
766  size_t idx = 1;
767  std::vector<Distribution *>::iterator addedDistIt;
768  for (addedDistIt = addedDistributions_m.begin();
769  addedDistIt != addedDistributions_m.end(); ++ addedDistIt, ++ idx) {
770 
771  particlesPerDist_m[idx] = (*addedDistIt)->tOrZDist_m.size();
772 
773  for (double dist : (*addedDistIt)->getXDist()) {
774  xDist_m.push_back(dist);
775  }
776  (*addedDistIt)->eraseXDist();
777 
778  for (double dist : (*addedDistIt)->getBGxDist()) {
779  pxDist_m.push_back(dist);
780  }
781  (*addedDistIt)->eraseBGxDist();
782 
783  for (double dist : (*addedDistIt)->getYDist()) {
784  yDist_m.push_back(dist);
785  }
786  (*addedDistIt)->eraseYDist();
787 
788  for (double dist : (*addedDistIt)->getBGyDist()) {
789  pyDist_m.push_back(dist);
790  }
791  (*addedDistIt)->eraseBGyDist();
792 
793  for (double dist : (*addedDistIt)->getTOrZDist()) {
794  tOrZDist_m.push_back(dist);
795  }
796  (*addedDistIt)->eraseTOrZDist();
797 
798  for (double dist : (*addedDistIt)->getBGzDist()) {
799  pzDist_m.push_back(dist);
800  }
801  (*addedDistIt)->eraseBGzDist();
802  }
803 }
804 
805 void Distribution::applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs) {
806 
807  switch (emissionModel_m) {
808 
809  case EmissionModelT::NONE:
810  applyEmissModelNone(pz);
811  break;
812  case EmissionModelT::ASTRA:
813  applyEmissModelAstra(px, py, pz, additionalRNs);
814  break;
815  case EmissionModelT::NONEQUIL:
816  applyEmissModelNonEquil(lowEnergyLimit, px, py, pz, additionalRNs);
817  break;
818  default:
819  break;
820  }
821 }
822 
823 void Distribution::applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs) {
824 
825  double phi = 2.0 * acos(sqrt(additionalRNs[0]));
826  double theta = 2.0 * Physics::pi * additionalRNs[1];
827 
828  px = pTotThermal_m * sin(phi) * cos(theta);
829  py = pTotThermal_m * sin(phi) * sin(theta);
830  pz = pTotThermal_m * std::abs(cos(phi));
831 
832 }
833 
834 void Distribution::applyEmissModelNone(double &pz) {
835  pz += pTotThermal_m;
836 }
837 
838 void Distribution::applyEmissModelNonEquil(double lowEnergyLimit,
839  double &bgx,
840  double &bgy,
841  double &bgz,
842  std::vector<double> &additionalRNs) {
843 
844  // Generate emission energy.
845  double energy = 0.0;
846  bool allow = false;
847 
848  const double expRelativeLaserEnergy = exp(laserEnergy_m / cathodeTemp_m);
849  // double energyRange = emitEnergyUpperLimit_m - lowEnergyLimit;
850  unsigned int counter = 0;
851  while (!allow) {
852  energy = lowEnergyLimit + additionalRNs[counter++] * emitEnergyUpperLimit_m;
853  double randFuncValue = additionalRNs[counter++];
854  double expRelativeEnergy = exp((energy - cathodeFermiEnergy_m) / cathodeTemp_m);
855  double funcValue = ((1.0
856  - 1.0 / (1.0 + expRelativeEnergy * expRelativeLaserEnergy)) /
857  (1.0 + expRelativeEnergy));
858 
859  if (randFuncValue <= funcValue)
860  allow = true;
861 
862  if (counter == additionalRNs.size()) {
863  for (unsigned int i = 0; i < counter; ++ i) {
864  additionalRNs[i] = gsl_rng_uniform(randGen_m);
865  }
866 
867  counter = 0;
868  }
869  }
870 
871  while (additionalRNs.size() - counter < 2) {
872  -- counter;
873  additionalRNs[counter] = gsl_rng_uniform(randGen_m);
874  }
875 
876  // Compute emission angles.
877  double energyInternal = energy + laserEnergy_m;
878  double energyExternal = energy - lowEnergyLimit; // uniformly distributed (?) value between 0 and emitEnergyUpperLimit_m
879 
880  double thetaInMax = acos(sqrt((lowEnergyLimit + laserEnergy_m) / (energyInternal)));
881  double thetaIn = additionalRNs[counter++] * thetaInMax;
882  double sinThetaOut = sin(thetaIn) * sqrt(energyInternal / energyExternal);
883  double phi = Physics::two_pi * additionalRNs[counter];
884 
885  // Compute emission momenta.
886  double betaGammaExternal
887  = sqrt(pow(energyExternal / (Physics::m_e * 1.0e9) + 1.0, 2.0) - 1.0);
888 
889  bgx = betaGammaExternal * sinThetaOut * cos(phi);
890  bgy = betaGammaExternal * sinThetaOut * sin(phi);
891  bgz = betaGammaExternal * sqrt(1.0 - pow(sinThetaOut, 2.0));
892 
893 }
894 
895 void Distribution::calcPartPerDist(size_t numberOfParticles) {
896 
897  if (numberOfDistributions_m == 1) {
898  particlesPerDist_m.push_back(numberOfParticles);
899  return;
900  }
901 
902  std::map<unsigned int, size_t> nPartFromFiles;
903  double totalWeight = 0.0;
904  for (unsigned int i = 0; i <= addedDistributions_m.size(); ++ i) {
905  Distribution *currDist = this;
906  if (i > 0)
907  currDist = addedDistributions_m[i - 1];
908 
909  if (currDist->distrTypeT_m == DistrTypeT::FROMFILE) {
910  std::ifstream inputFile;
911  if (Ippl::myNode() == 0) {
912  std::string fileName = Attributes::getString(currDist->itsAttr[Attrib::Distribution::FNAME]);
913  inputFile.open(fileName.c_str());
914  }
915  size_t nPart = getNumberOfParticlesInFile(inputFile);
916  nPartFromFiles.insert(std::make_pair(i, nPart));
917  if (nPart > numberOfParticles) {
918  throw OpalException("Distribution::calcPartPerDist",
919  "Number of particles is too small");
920  }
921 
922  numberOfParticles -= nPart;
923  } else {
924  totalWeight += currDist->getWeight();
925  }
926  }
927 
928  size_t numberOfCommittedPart = 0;
929  for (unsigned int i = 0; i <= addedDistributions_m.size(); ++ i) {
930  Distribution *currDist = this;
931  if (i > 0)
932  currDist = addedDistributions_m[i - 1];
933 
934  if (currDist->distrTypeT_m == DistrTypeT::FROMFILE) {
935  particlesPerDist_m.push_back(nPartFromFiles[i]);
936  } else {
937  size_t particlesCurrentDist = numberOfParticles * currDist->getWeight() / totalWeight;
938  particlesPerDist_m.push_back(particlesCurrentDist);
939  numberOfCommittedPart += particlesCurrentDist;
940  }
941  }
942 
943  // Remaining particles go into first distribution that isn't FROMFILE
944  if (numberOfParticles != numberOfCommittedPart) {
945  size_t diffNumber = numberOfParticles - numberOfCommittedPart;
946  for (unsigned int i = 0; i <= addedDistributions_m.size(); ++ i) {
947  Distribution *currDist = this;
948  if (i > 0)
949  currDist = addedDistributions_m[i - 1];
950 
951  if (currDist->distrTypeT_m != DistrTypeT::FROMFILE) {
952  particlesPerDist_m.at(i) += diffNumber;
953  diffNumber = 0;
954  break;
955  }
956  }
957  if (diffNumber != 0) {
958  throw OpalException("Distribution::calcPartPerDist",
959  "Particles can't be distributed to distributions in array");
960  }
961  }
962 }
963 
964 void Distribution::checkEmissionParameters() {
965 
966  // There must be at least on energy bin for an emitted beam->
967  numberOfEnergyBins_m
968  = std::abs(static_cast<int> (Attributes::getReal(itsAttr[Attrib::Distribution::NBIN])));
969  if (numberOfEnergyBins_m == 0)
970  numberOfEnergyBins_m = 1;
971 
972  int emissionSteps = std::abs(static_cast<int> (Attributes::getReal(itsAttr[Attrib::Distribution::EMISSIONSTEPS])));
973 
974  // There must be at least one emission step.
975  if (emissionSteps == 0)
976  emissionSteps = 1;
977 
978  // Set number of sample bins per energy bin from the number of emission steps.
979  numberOfSampleBins_m = static_cast<int> (std::ceil(emissionSteps / numberOfEnergyBins_m));
980  if (numberOfSampleBins_m <= 0)
981  numberOfSampleBins_m = 1;
982 
983  // Initialize emission counters.
984  currentEnergyBin_m = 1;
985  currentSampleBin_m = 0;
986 
987 }
988 
989 void Distribution::checkIfEmitted() {
990 
991  emitting_m = Attributes::getBool(itsAttr[Attrib::Distribution::EMITTED]);
992 
993  switch (distrTypeT_m) {
994 
995  case DistrTypeT::ASTRAFLATTOPTH:
996  case DistrTypeT::GUNGAUSSFLATTOPTH:
997  emitting_m = true;
998  break;
999  default:
1000  break;
1001  }
1002 }
1003 
1004 void Distribution::checkParticleNumber(size_t &numberOfParticles) {
1005 
1006  size_t numberOfDistParticles = tOrZDist_m.size();
1007  reduce(numberOfDistParticles, numberOfDistParticles, OpAddAssign());
1008 
1009  if (numberOfDistParticles != numberOfParticles) {
1010  *gmsg << "\n--------------------------------------------------" << endl
1011  << "Warning!! The number of particles in the initial" << endl
1012  << "distribution is " << numberOfDistParticles << "." << endl << endl
1013  << "This is different from the number of particles" << endl
1014  << "defined by the BEAM command: " << numberOfParticles << endl << endl
1015  << "This is often happens when using a FROMFILE type" << endl
1016  << "distribution and not matching the number of" << endl
1017  << "particles in the particles file(s) with the number" << endl
1018  << "given in the BEAM command." << endl << endl
1019  << "The number of particles in the initial distribution" << endl
1020  << "(" << numberOfDistParticles << ") "
1021  << "will take precedence." << endl
1022  << "---------------------------------------------------\n" << endl;
1023  }
1024  numberOfParticles = numberOfDistParticles;
1025 }
1026 
1027 void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits) {
1028 
1029  /*
1030  * Toggle what units to use for inputing momentum.
1031  */
1032  std::string inputUnits = Util::toUpper(Attributes::getString(itsAttr[Attrib::Distribution::INPUTMOUNITS]));
1033  if (inputUnits == "NONE")
1034  inputMoUnits_m = InputMomentumUnitsT::NONE;
1035  else if (inputUnits == "EV")
1036  inputMoUnits_m = InputMomentumUnitsT::EV;
1037  else
1038  inputMoUnits_m = inputMoUnits;
1039 
1040 }
1041 
1042 double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
1043  double value = std::copysign(sqrt(pow(std::abs(valueIneV) / massIneV + 1.0, 2.0) - 1.0), valueIneV);
1044  if (std::abs(value) < std::numeric_limits<double>::epsilon())
1045  value = std::copysign(sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
1046 
1047  return value;
1048 }
1049 
1050 void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
1051 
1052  setDistParametersBinomial(massIneV);
1053  generateBinomial(numberOfParticles);
1054 }
1055 
1056 void Distribution::createDistributionFlattop(size_t numberOfParticles, double massIneV) {
1057 
1058  setDistParametersFlattop(massIneV);
1059 
1060  if (emitting_m) {
1061  if (laserProfile_m == NULL)
1062  generateFlattopT(numberOfParticles);
1063  else
1064  generateFlattopLaserProfile(numberOfParticles);
1065  } else
1066  generateFlattopZ(numberOfParticles);
1067 }
1068 
1069 size_t Distribution::getNumberOfParticlesInFile(std::ifstream &inputFile) {
1070 
1071  size_t numberOfParticlesRead = 0;
1072  if (Ippl::myNode() == 0) {
1073  const boost::regex commentExpr("[[:space:]]*#.*");
1074  const std::string repl("");
1075  std::string line;
1076  std::stringstream linestream;
1077  long tempInt = 0;
1078 
1079  do {
1080  std::getline(inputFile, line);
1081  line = boost::regex_replace(line, commentExpr, repl);
1082  } while (line.length() == 0 && !inputFile.fail());
1083 
1084  linestream.str(line);
1085  linestream >> tempInt;
1086  if (tempInt <= 0) {
1087  throw OpalException("Distribution::getNumberOfParticlesInFile",
1088  "The file '" +
1089  Attributes::getString(itsAttr[Attrib::Distribution::FNAME]) +
1090  "' does not seem to be an ASCII file containing a distribution.");
1091  }
1092  numberOfParticlesRead = static_cast<size_t>(tempInt);
1093  }
1094  reduce(numberOfParticlesRead, numberOfParticlesRead, OpAddAssign());
1095 
1096  return numberOfParticlesRead;
1097 }
1098 
1099 void Distribution::createDistributionFromFile(size_t numberOfParticles, double massIneV) {
1100 
1101  *gmsg << level3 << "\n"
1102  << "------------------------------------------------------------------------------------\n";
1103  *gmsg << "READ INITIAL DISTRIBUTION FROM FILE \""
1104  << Attributes::getString(itsAttr[Attrib::Distribution::FNAME])
1105  << "\"\n";
1106  *gmsg << "------------------------------------------------------------------------------------\n" << endl;
1107 
1108  // Data input file is only read by node 0.
1109  std::ifstream inputFile;
1110  std::string fileName = Attributes::getString(itsAttr[Attrib::Distribution::FNAME]);
1111  if (Ippl::myNode() == 0) {
1112  inputFile.open(fileName.c_str());
1113  if (inputFile.fail())
1114  throw OpalException("Distribution::createDistributionFromFile",
1115  "Open file operation failed, please check if \""
1116  + fileName
1117  + "\" really exists.");
1118  }
1119 
1120  size_t numberOfParticlesRead = getNumberOfParticlesInFile(inputFile);
1121  /*
1122  * We read in the particle information using node zero, but save the particle
1123  * data to each processor in turn.
1124  */
1125  unsigned int saveProcessor = 0;
1126 
1127  pmean_m = 0.0;
1128 
1129  size_t numPartsToSend = 0;
1130  unsigned int distributeFrequency = 1000;
1131  size_t singleDataSize = (/*sizeof(int) +*/ 6 * sizeof(double));
1132  unsigned int dataSize = distributeFrequency * singleDataSize;
1133  std::vector<char> data;
1134 
1135  data.reserve(dataSize);
1136 
1137  const char* buffer;
1138  if (Ippl::myNode() == 0) {
1139  char lineBuffer[1024];
1140  unsigned int numParts = 0;
1141  while (!inputFile.eof()) {
1142  inputFile.getline(lineBuffer, 1024);
1143 
1144  Vector_t R(0.0), P(0.0);
1145 
1146  std::istringstream line(lineBuffer);
1147  line >> R(0);
1148  if (line.rdstate()) break;
1149  line >> P(0);
1150  line >> R(1);
1151  line >> P(1);
1152  line >> R(2);
1153  line >> P(2);
1154 
1155  if (saveProcessor >= (unsigned)Ippl::getNodes())
1156  saveProcessor = 0;
1157 
1158  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
1159  P(0) = converteVToBetaGamma(P(0), massIneV);
1160  P(1) = converteVToBetaGamma(P(1), massIneV);
1161  P(2) = converteVToBetaGamma(P(2), massIneV);
1162  }
1163 
1164  pmean_m += P;
1165 
1166  if (saveProcessor > 0u) {
1167  buffer = reinterpret_cast<const char*>(&R[0]);
1168  data.insert(data.end(), buffer, buffer + 3 * sizeof(double));
1169  buffer = reinterpret_cast<const char*>(&P[0]);
1170  data.insert(data.end(), buffer, buffer + 3 * sizeof(double));
1171  ++ numPartsToSend;
1172 
1173  if (numPartsToSend % distributeFrequency == 0) {
1174  MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
1175  MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
1176  numPartsToSend = 0;
1177 
1178  std::vector<char>().swap(data);
1179  data.reserve(dataSize);
1180  }
1181  } else {
1182  xDist_m.push_back(R(0));
1183  yDist_m.push_back(R(1));
1184  tOrZDist_m.push_back(R(2));
1185  pxDist_m.push_back(P(0));
1186  pyDist_m.push_back(P(1));
1187  pzDist_m.push_back(P(2));
1188  }
1189 
1190  ++ numParts;
1191  ++ saveProcessor;
1192  }
1193 
1194  dataSize = (numberOfParticlesRead == numParts? data.size(): std::numeric_limits<unsigned int>::max());
1195  MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
1196  if (numberOfParticlesRead != numParts) {
1197  throw OpalException("Distribution::createDistributionFromFile",
1198  "Found " +
1199  std::to_string(numParts) +
1200  " particles in file '" +
1201  fileName +
1202  "' instead of " +
1203  std::to_string(numberOfParticlesRead));
1204  }
1205  MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
1206 
1207  } else {
1208  do {
1209  MPI_Bcast(&dataSize, 1, MPI_INT, 0, Ippl::getComm());
1210  if (dataSize == std::numeric_limits<unsigned int>::max()) {
1211  throw OpalException("Distribution::createDistributionFromFile",
1212  "Couldn't find " +
1213  std::to_string(numberOfParticlesRead) +
1214  " particles in file '" +
1215  fileName + "'");
1216  }
1217  MPI_Bcast(&data[0], dataSize, MPI_CHAR, 0, Ippl::getComm());
1218 
1219  size_t i = 0;
1220  while (i < dataSize) {
1221 
1222  if (saveProcessor + 1 == (unsigned) Ippl::myNode()) {
1223  const double *tmp = reinterpret_cast<const double*>(&data[i]);
1224  xDist_m.push_back(tmp[0]);
1225  yDist_m.push_back(tmp[1]);
1226  tOrZDist_m.push_back(tmp[2]);
1227  pxDist_m.push_back(tmp[3]);
1228  pyDist_m.push_back(tmp[4]);
1229  pzDist_m.push_back(tmp[5]);
1230  i += 6 * sizeof(double);
1231  } else {
1232  i += singleDataSize;
1233  }
1234 
1235  ++ saveProcessor;
1236  if (saveProcessor + 1 >= (unsigned) Ippl::getNodes()) {
1237  saveProcessor = 0;
1238  }
1239  }
1240  } while (dataSize == distributeFrequency * singleDataSize);
1241  }
1242 
1243  pmean_m /= numberOfParticlesRead;
1244  reduce(pmean_m, pmean_m, OpAddAssign());
1245 
1246  if (Ippl::myNode() == 0)
1247  inputFile.close();
1248 }
1249 
1250 
1251 void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, double massIneV) {
1252 
1253  /*
1254  ToDo:
1255  - eliminate physics and error
1256  */
1257 
1258  std::string LineName = Attributes::getString(itsAttr[Attrib::Distribution::LINE]);
1259  if (LineName == "") return;
1260 
1261  const BeamSequence* LineSequence = BeamSequence::find(LineName);
1262 
1263  if (LineSequence == NULL)
1264  throw OpalException("Distribution::CreateMatchedGaussDistribution",
1265  "didn't find any Cyclotron element in line");
1266 
1267  SpecificElementVisitor<Cyclotron> CyclotronVisitor(*LineSequence->fetchLine());
1268  CyclotronVisitor.execute();
1269  size_t NumberOfCyclotrons = CyclotronVisitor.size();
1270 
1271  if (NumberOfCyclotrons > 1) {
1272  throw OpalException("Distribution::createMatchedGaussDistribution",
1273  "I am confused, found more than one Cyclotron element in line");
1274  }
1275  if (NumberOfCyclotrons == 0) {
1276  throw OpalException("Distribution::createMatchedGaussDistribution",
1277  "didn't find any Cyclotron element in line");
1278  }
1279 
1280  /* FIXME we need to remove const-ness otherwise we can't update the injection radius
1281  * and injection radial momentum. However, there should be a better solution ..
1282  */
1283  Cyclotron* CyclotronElement = const_cast<Cyclotron*>(CyclotronVisitor.front());
1284 
1285  bool full = !Attributes::getBool(itsAttr[Attrib::Distribution::SECTOR]);
1286  int Nint = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSTEPS]);
1287  int Nsectors = (int)Attributes::getReal(itsAttr[Attrib::Distribution::NSECTORS]);
1288 
1289  if ( Nint < 0 )
1290  throw OpalException("Distribution::CreateMatchedGaussDistribution()",
1291  "Negative number of integration steps");
1292 
1293  if ( Nsectors < 0 )
1294  throw OpalException("Distribution::CreateMatchedGaussDistribution()",
1295  "Negative number of sectors");
1296 
1297  if ( Nsectors > 1 && full == false )
1298  throw OpalException("Distribution::CreateMatchedGaussDistribution()",
1299  "Averaging over sectors can only be done with SECTOR=FALSE");
1300 
1301  *gmsg << "* ----------------------------------------------------" << endl;
1302  *gmsg << "* About to find closed orbit and matched distribution " << endl;
1303  *gmsg << "* I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
1304  *gmsg << "* EX= " << Attributes::getReal(itsAttr[Attrib::Distribution::EX])
1305  << " EY= " << Attributes::getReal(itsAttr[Attrib::Distribution::EY])
1306  << " ET= " << Attributes::getReal(itsAttr[Attrib::Distribution::ET]) << endl;
1307  if ( full ) {
1308  *gmsg << "* SECTOR: " << "match using all sectors, with" << endl
1309  << "* NSECTORS = " << Nsectors << " to average the field over" << endl;
1310  }
1311  else
1312  *gmsg << "* SECTOR: " << "match using single sector" << endl;
1313 
1314  *gmsg << "* NSTEPS = " << Nint << endl
1315  << "* HN= " << CyclotronElement->getCyclHarm()
1316  << " PHIINIT= " << CyclotronElement->getPHIinit() << endl
1317  << "* ----------------------------------------------------" << endl;
1318 
1319  if ( CyclotronElement->getFMLowE() < 0 ||
1320  CyclotronElement->getFMHighE() < 0 )
1321  {
1322  throw OpalException("Distribution::CreateMatchedGaussDistribution()",
1323  "Missing attributes 'FMLOWE' and/or 'FMHIGHE' in "
1324  "'CYCLOTRON' definition.");
1325  }
1326 
1327  std::size_t maxitCOF =
1328  Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSCO]);
1329 
1330  double rguess =
1331  Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]);
1332 
1333  double denergy = 1000.0 *
1334  Attributes::getReal(itsAttr[Attrib::Distribution::DENERGY]);
1335 
1336  if ( denergy < 0.0 )
1337  throw OpalException("Distribution:CreateMatchedGaussDistribution()",
1338  "DENERGY < 0");
1339 
1340  double accuracy =
1341  Attributes::getReal(itsAttr[Attrib::Distribution::RESIDUUM]);
1342 
1343  if ( Options::cloTuneOnly ) {
1344  *gmsg << "* Stopping after closed orbit and tune calculation" << endl;
1345  typedef std::vector<double> container_t;
1346  typedef boost::numeric::odeint::runge_kutta4<container_t> rk4_t;
1347  typedef ClosedOrbitFinder<double,unsigned int, rk4_t> cof_t;
1348 
1349  cof_t cof(massIneV*1E-6, Nint, CyclotronElement, false, Nsectors);
1350  cof.findOrbit(accuracy, maxitCOF, E_m*1E-6, denergy, rguess, true);
1351 
1352  throw EarlyLeaveException("Distribution::CreateMatchedGaussDistribution()",
1353  "Do only tune calculation.");
1354  }
1355 
1356  bool writeMap = true;
1357 
1358  typedef SigmaGenerator<double, unsigned int> sGenerator_t;
1359  sGenerator_t *siggen = new sGenerator_t(I_m,
1360  Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
1361  Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
1362  Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
1363  E_m*1E-6,
1364  massIneV*1E-6,
1365  CyclotronElement,
1366  Nint,
1367  Nsectors,
1368  Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
1369  writeMap);
1370 
1371  if (siggen->match(accuracy,
1372  Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
1373  maxitCOF,
1374  CyclotronElement,
1375  denergy,
1376  rguess,
1377  false, full)) {
1378 
1379  std::array<double,3> Emit = siggen->getEmittances();
1380 
1381  if (rguess > 0)
1382  *gmsg << "* RGUESS " << rguess << " (m) " << endl;
1383 
1384  *gmsg << "* Converged (Ex, Ey, Ez) = (" << Emit[0] << ", " << Emit[1] << ", "
1385  << Emit[2] << ") pi mm mrad for E= " << E_m*1E-6 << " (MeV)" << endl;
1386  *gmsg << "* Sigma-Matrix " << endl;
1387 
1388  for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
1389  *gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
1390  for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
1391  if (siggen->getSigma()(i,j) < 10e-12){
1392  *gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
1393  }
1394  else{
1395  *gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
1396  }
1397 
1398  }
1399  *gmsg << " \\\\" << endl;
1400  }
1401 
1402  /*
1403 
1404  Now setup the distribution generator
1405  Units of the Sigma Matrix:
1406 
1407  spatial: mm
1408  moment: rad
1409 
1410  */
1411  double gamma = E_m / massIneV + 1.0;
1412  double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
1413 
1414  auto sigma = siggen->getSigma();
1415  // change units from mm to m
1416  for (unsigned int i = 0; i < 6; ++ i)
1417  for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
1418 
1419  for (unsigned int i = 0; i < 3; ++ i) {
1420  if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
1421  throw OpalException("Distribution::CreateMatchedGaussDistribution()",
1422  "Negative value on the diagonal of the sigma matrix.");
1423  }
1424 
1425  sigmaR_m[0] = std::sqrt(sigma(0, 0));
1426  sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
1427  sigmaR_m[2] = std::sqrt(sigma(2, 2));
1428  sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
1429  sigmaR_m[1] = std::sqrt(sigma(4, 4));
1430 
1431  //p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
1432 
1433  double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma -
1434  sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
1435 
1436  double pl = std::sqrt(pl2);
1437  sigmaP_m[1] = gamma*(pl - beta*gamma);
1438 
1439  correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
1440  correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
1441  correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5)));
1442  correlationMatrix_m(4, 0) = sigma(0, 4) / (sqrt(sigma(0, 0) * sigma(4, 4)));
1443  correlationMatrix_m(4, 1) = sigma(1, 4) / (sqrt(sigma(1, 1) * sigma(4, 4)));
1444  correlationMatrix_m(5, 0) = sigma(0, 5) / (sqrt(sigma(0, 0) * sigma(5, 5)));
1445  correlationMatrix_m(5, 1) = sigma(1, 5) / (sqrt(sigma(1, 1) * sigma(5, 5)));
1446 
1447  createDistributionGauss(numberOfParticles, massIneV);
1448 
1449  // update injection radius and radial momentum
1450  CyclotronElement->setRinit(siggen->getInjectionRadius() * 1.0e3);
1451  CyclotronElement->setPRinit(siggen->getInjectionMomentum());
1452  }
1453  else {
1454  *gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
1455 
1456  delete siggen;
1457 
1458  throw OpalException("Distribution::CreateMatchedGaussDistribution",
1459  "didn't find any matched distribution.");
1460  }
1461 
1462  delete siggen;
1463 }
1464 
1465 void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV) {
1466 
1467  setDistParametersGauss(massIneV);
1468 
1469  if (emitting_m) {
1470  generateTransverseGauss(numberOfParticles);
1471  generateLongFlattopT(numberOfParticles);
1472  } else {
1473  generateGaussZ(numberOfParticles);
1474  }
1475 }
1476 
1477 void Distribution::createBoundaryGeometry(PartBunchBase<double, 3> *beam, BoundaryGeometry &bg) {
1478 
1479  int pc = 0;
1480  size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
1481  size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
1482  if (Ippl::myNode() == 0)
1483  N_mean += N_extra;
1484  size_t count = 0;
1485  size_t lowMark = beam->getLocalNum();
1486  if (bg.getN() != 0) {
1487 
1488  for (size_t i = 0; i < bg.getN(); i++) {
1489  if (pc == Ippl::myNode()) {
1490  if (count < N_mean) {
1491  beam->create(1);
1492  if (pc != 0)
1493  beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra);
1494  else
1495  beam->R[lowMark + count] = bg.getCooridinate(count);
1496  beam->P[lowMark + count] = Vector_t(0.0);
1497  beam->Bin[lowMark + count] = 0;
1498  beam->PType[lowMark + count] = ParticleType::FIELDEMISSION;
1499  beam->TriID[lowMark + count] = 0;
1500  beam->Q[lowMark + count] = beam->getChargePerParticle();
1501  beam->Ef[lowMark + count] = Vector_t(0.0);
1502  beam->Bf[lowMark + count] = Vector_t(0.0);
1503  beam->dt[lowMark + count] = beam->getdT();
1504  count ++;
1505  }
1506  }
1507  pc++;
1508  if (pc == Ippl::getNodes())
1509  pc = 0;
1510  }
1511  }
1512  bg.clearCooridinateArray();
1513  beam->boundp();
1514  *gmsg << &beam << endl;
1515 }
1516 
1517 void Distribution::createOpalCycl(PartBunchBase<double, 3> *beam,
1518  size_t numberOfParticles,
1519  double current, const Beamline &bl) {
1520 
1521  /*
1522  * setup data for matched distribution generation
1523  */
1524 
1525  E_m = (beam->getInitialGamma()-1.0)*beam->getM();
1526  I_m = current;
1527 
1528  /*
1529  Fixme:
1530 
1531  avrgpz_m = beam->getP()/beam->getM();
1532  */
1533  size_t numberOfPartToCreate = numberOfParticles;
1534  totalNumberParticles_m = numberOfParticles;
1535  if (beam->getTotalNum() != 0) {
1536  numberOfPartToCreate = beam->getLocalNum();
1537  }
1538 
1539  // Setup particle bin structure.
1540  setupParticleBins(beam->getM(),beam);
1541 
1542  /*
1543  * Set what units to use for input momentum units. Default is
1544  * eV.
1545  */
1546  chooseInputMomentumUnits(InputMomentumUnitsT::EV);
1547 
1548  /*
1549  * Determine the number of particles for each distribution. For OPAL-cycl
1550  * there are currently no arrays of distributions supported
1551  */
1552  calcPartPerDist(numberOfParticles);
1553 
1554  // Set distribution type.
1555  setDistType();
1556 
1557  // Emitting particles is not supported in OPAL-cycl.
1558  emitting_m = false;
1559 
1560  // Create distribution.
1561  create(numberOfPartToCreate, beam->getM());
1562 
1563  // this variable is needed to be compatible with OPAL-T
1564  particlesPerDist_m.push_back(tOrZDist_m.size());
1565 
1566  shiftDistCoordinates(beam->getM());
1567 
1568  // Setup energy bins.
1569  if (numberOfEnergyBins_m > 0) {
1570  fillParticleBins();
1571  beam->setPBins(energyBins_m);
1572  }
1573 
1574  // Check number of particles in distribution.
1575  checkParticleNumber(numberOfParticles);
1576 
1577  initializeBeam(beam);
1578  writeOutFileHeader();
1579 
1580  injectBeam(beam);
1581 
1582  OpalData::getInstance()->addProblemCharacteristicValue("NP", numberOfParticles);
1583 }
1584 
1585 void Distribution::createOpalE(Beam *beam,
1586  std::vector<Distribution *> addedDistributions,
1587  EnvelopeBunch *envelopeBunch,
1588  double distCenter,
1589  double Bz0) {
1590 
1591  IpplTimings::startTimer(envelopeBunch->distrCreate_m);
1592 
1593  double beamEnergy = beam->getMass() * (beam->getGamma() - 1.0) * 1.0e9;
1594  numberOfEnergyBins_m = static_cast<int>(fabs(Attributes::getReal(itsAttr[Attrib::Distribution::NBIN])));
1595 
1596  /*
1597  * Set what units to use for input momentum units. Default is
1598  * unitless (i.e. BetaXGamma, BetaYGamma, BetaZGamma).
1599  */
1600  chooseInputMomentumUnits(InputMomentumUnitsT::NONE);
1601 
1602  setDistType();
1603 
1604  // Check if this is to be an emitted beam->
1605  checkIfEmitted();
1606 
1607  switch (distrTypeT_m) {
1608 
1609  case DistrTypeT::FLATTOP:
1610  setDistParametersFlattop(beam->getMass());
1611  beamEnergy = Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]);
1612  break;
1613  case DistrTypeT::GAUSS:
1614  setDistParametersGauss(beam->getMass());
1615  break;
1616  case DistrTypeT::GUNGAUSSFLATTOPTH:
1617  setDistParametersFlattop(beam->getMass());
1618  beamEnergy = Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]);
1619  break;
1620  default:
1621  *gmsg << "Only FLATTOP, GAUSS and GUNGAUSSFLATTOPTH distribution types supported " << endl
1622  << "in envelope mode. Assuming FLATTOP." << endl;
1623  distrTypeT_m = DistrTypeT::FLATTOP;
1624  setDistParametersFlattop(beam->getMass());
1625  break;
1626  }
1627 
1628  tEmission_m = tPulseLengthFWHM_m + (cutoffR_m[2] - sqrt(2.0 * log(2.0)))
1629  * (sigmaTRise_m + sigmaTFall_m);
1630  double beamWidth = tEmission_m * Physics::c * sqrt(1.0 - (1.0 / pow(beam->getGamma(), 2)));
1631  double beamCenter = -1.0 * beamWidth / 2.0;
1632 
1633  envelopeBunch->initialize(beam->getNumberOfSlices(),
1634  beam->getCharge(),
1635  beamEnergy,
1636  beamWidth,
1637  tEmission_m,
1638  0.9,
1639  beam->getCurrent(),
1640  beamCenter,
1641  sigmaR_m[0],
1642  sigmaR_m[1],
1643  0.0,
1644  0.0,
1645  Bz0,
1646  numberOfEnergyBins_m);
1647 
1648  IpplTimings::stopTimer(envelopeBunch->distrCreate_m);
1649 }
1650 
1651 void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
1652  std::vector<Distribution *> addedDistributions,
1653  size_t &numberOfParticles) {
1654 
1655  addedDistributions_m = addedDistributions;
1656  createOpalT(beam, numberOfParticles);
1657 }
1658 
1659 void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
1660  size_t &numberOfParticles) {
1661 
1662  IpplTimings::startTimer(beam->distrCreate_m);
1663 
1664  // This is PC from BEAM
1665  double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
1666  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
1667  deltaP = converteVToBetaGamma(deltaP, beam->getM());
1668  }
1669 
1670  avrgpz_m = beam->getP()/beam->getM() + deltaP;
1671 
1672  totalNumberParticles_m = numberOfParticles;
1673 
1674  /*
1675  * Set what units to use for input momentum units. Default is
1676  * unitless (i.e. BetaXGamma, BetaYGamma, BetaZGamma).
1677  */
1678  chooseInputMomentumUnits(InputMomentumUnitsT::NONE);
1679 
1680  // Set distribution type(s).
1681  setDistType();
1682  for (Distribution* addedDist : addedDistributions_m)
1683  addedDist->setDistType();
1684 
1685  /*
1686  * Determine the number of particles for each distribution. Note
1687  * that if a distribution is generated from an input file, then
1688  * the number of particles in that file will override what is
1689  * determined here.
1690  */
1691  calcPartPerDist(numberOfParticles);
1692 
1693  // Check if this is to be an emitted beam->
1694  checkIfEmitted();
1695 
1696  /*
1697  * Force added distributions to the same emission condition as the main
1698  * distribution.
1699  */
1700  for (Distribution* addedDist : addedDistributions_m)
1701  addedDist->setDistToEmitted(emitting_m);
1702 
1703  if (emitting_m)
1704  setupEmissionModel(beam);
1705 
1706  // Create distributions.
1707  create(particlesPerDist_m.at(0), beam->getM());
1708 
1709  size_t distCount = 1;
1710  for (Distribution* addedDist : addedDistributions_m) {
1711  addedDist->create(particlesPerDist_m.at(distCount), beam->getM());
1712  distCount++;
1713  }
1714 
1715  // Move added distribution particles to main distribution.
1716  addDistributions();
1717 
1718  if (emitting_m && emissionModel_m == EmissionModelT::NONE)
1719  setupEmissionModelNone(beam);
1720 
1721  // Check number of particles in distribution.
1722  checkParticleNumber(numberOfParticles);
1723 
1724  if (emitting_m) {
1725  checkEmissionParameters();
1726  } else {
1727  if (distrTypeT_m != DistrTypeT::FROMFILE) {
1728  pmean_m = Vector_t(0, 0, avrgpz_m);
1729  }
1730  }
1731 
1732  /*
1733  * Find max. and min. particle positions in the bunch. For the
1734  * case of an emitted beam these will be in time (seconds). For
1735  * an injected beam in z (meters).
1736  */
1737  double maxTOrZ = getMaxTOrZ();
1738  double minTOrZ = getMinTOrZ();
1739 
1740  /*
1741  * Set emission time and/or shift distributions if necessary.
1742  * For an emitted beam we place all particles at negative time.
1743  * For an injected beam we just ensure that there are no
1744  * particles at z < 0.
1745  */
1746 
1747  if (emitting_m) {
1748  setEmissionTime(maxTOrZ, minTOrZ);
1749  }
1750  shiftBeam(maxTOrZ, minTOrZ);
1751 
1752  shiftDistCoordinates(beam->getM());
1753 
1754  if (numberOfEnergyBins_m > 0) {
1755  setupEnergyBins(maxTOrZ, minTOrZ);
1756  fillEBinHistogram();
1757  }
1758 
1759  initializeBeam(beam);
1760  writeOutFileHeader();
1761 
1762  if (emitting_m && Options::cZero) {
1763  std::vector<std::vector<double> > mirrored;
1764  const auto end = additionalRNs_m.end();
1765 
1766  if (emissionModel_m == EmissionModelT::ASTRA ||
1767  distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
1768  distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
1769 
1770  for (auto it = additionalRNs_m.begin(); it != end; ++ it) {
1771  std::vector<double> tmp;
1772  tmp.push_back((*it).front());
1773  tmp.push_back((*it).back() + 0.5);
1774  mirrored.push_back(tmp);
1775  }
1776  } else {
1777  size_t numAdditionals = additionalRNs_m.front().size() / 2;
1778  for (auto it = additionalRNs_m.begin(); it != end; ++ it) {
1779  std::vector<double> tmp((*it).begin() + numAdditionals, (*it).end());
1780  mirrored.push_back(tmp);
1781  (*it).erase((*it).begin() + numAdditionals, (*it).end());
1782  }
1783  }
1784 
1785  additionalRNs_m.insert(additionalRNs_m.end(), mirrored.begin(), mirrored.end());
1786  }
1787 
1788  /*
1789  * If this is an injected beam, we create particles right away.
1790  * Emitted beams get created during the course of the simulation.
1791  */
1792  if (!emitting_m)
1793  injectBeam(beam);
1794 
1795  OpalData::getInstance()->addProblemCharacteristicValue("NP", numberOfParticles);
1796  IpplTimings::stopTimer(beam->distrCreate_m);
1797 }
1798 
1823 size_t Distribution::emitParticles(PartBunchBase<double, 3> *beam, double eZ) {
1824 
1825  // Number of particles that have already been emitted and are on this processor.
1826  size_t numberOfEmittedParticles = beam->getLocalNum();
1827  size_t oldNumberOfEmittedParticles = numberOfEmittedParticles;
1828 
1829  if (!tOrZDist_m.empty() && emitting_m) {
1830 
1831  /*
1832  * Loop through emission beam coordinate containers and emit particles with
1833  * the appropriate time coordinate. Once emitted, remove particle from the
1834  * "to be emitted" list.
1835  */
1836  std::vector<size_t> particlesToBeErased;
1837  double phiEffective = (cathodeWorkFunc_m
1838  - sqrt(std::max(0.0, (Physics::q_e * beam->getQ() * eZ) /
1839  (4.0 * Physics::pi * Physics::epsilon_0))));
1840  double lowEnergyLimit = cathodeFermiEnergy_m + phiEffective - laserEnergy_m;
1841 
1842  for (size_t particleIndex = 0; particleIndex < tOrZDist_m.size(); particleIndex++) {
1843 
1844  // Advance particle time.
1845  tOrZDist_m.at(particleIndex) += beam->getdT();
1846 
1847  // If particle time is now greater than zero, we emit it.
1848  if (tOrZDist_m.at(particleIndex) >= 0.0) {
1849 
1850  particlesToBeErased.push_back(particleIndex);
1851 
1852  beam->create(1);
1853  double deltaT = tOrZDist_m.at(particleIndex);
1854  beam->dt[numberOfEmittedParticles] = deltaT;
1855 
1856  double oneOverCDt = 1.0 / (Physics::c * deltaT);
1857 
1858  double px = pxDist_m.at(particleIndex);
1859  double py = pyDist_m.at(particleIndex);
1860  double pz = pzDist_m.at(particleIndex);
1861  std::vector<double> additionalRNs;
1862  if (additionalRNs_m.size() > particleIndex) {
1863  additionalRNs = additionalRNs_m[particleIndex];
1864  } else {
1865  throw OpalException("Distribution::emitParticles",
1866  "not enough additional particles");
1867  }
1868  applyEmissionModel(lowEnergyLimit, px, py, pz, additionalRNs);
1869 
1870  double particleGamma
1871  = std::sqrt(1.0
1872  + std::pow(px, 2.0)
1873  + std::pow(py, 2.0)
1874  + std::pow(pz, 2.0));
1875 
1876  beam->R[numberOfEmittedParticles]
1877  = Vector_t(oneOverCDt * (xDist_m.at(particleIndex)
1878  + px * deltaT * Physics::c / (2.0 * particleGamma)),
1879  oneOverCDt * (yDist_m.at(particleIndex)
1880  + py * deltaT * Physics::c / (2.0 * particleGamma)),
1881  pz / (2.0 * particleGamma));
1882  beam->P[numberOfEmittedParticles]
1883  = Vector_t(px, py, pz);
1884  beam->Bin[numberOfEmittedParticles] = currentEnergyBin_m - 1;
1885  beam->Q[numberOfEmittedParticles] = beam->getChargePerParticle();
1886  beam->Ef[numberOfEmittedParticles] = Vector_t(0.0);
1887  beam->Bf[numberOfEmittedParticles] = Vector_t(0.0);
1888  beam->PType[numberOfEmittedParticles] = ParticleType::REGULAR;
1889  beam->TriID[numberOfEmittedParticles] = 0;
1890  numberOfEmittedParticles++;
1891 
1892  beam->iterateEmittedBin(currentEnergyBin_m - 1);
1893 
1894  // Save particles to vectors for writing initial distribution.
1895  xWrite_m.push_back(xDist_m.at(particleIndex));
1896  pxWrite_m.push_back(px);
1897  yWrite_m.push_back(yDist_m.at(particleIndex));
1898  pyWrite_m.push_back(py);
1899  tOrZWrite_m.push_back(-(beam->getdT() - deltaT + currentEmissionTime_m));
1900  pzWrite_m.push_back(pz);
1901  binWrite_m.push_back(currentEnergyBin_m);
1902  }
1903  }
1904 
1905  // Now erase particles that were emitted.
1906  std::vector<size_t>::reverse_iterator ptbErasedIt;
1907  for (ptbErasedIt = particlesToBeErased.rbegin();
1908  ptbErasedIt < particlesToBeErased.rend();
1909  ++ptbErasedIt) {
1910 
1911  /*
1912  * We don't use the vector class function erase because it
1913  * is much slower than doing a swap and popping off the end
1914  * of the vector.
1915  */
1916  std::swap( xDist_m.at(*ptbErasedIt), xDist_m.back());
1917  std::swap(pxDist_m.at(*ptbErasedIt), pxDist_m.back());
1918  std::swap( yDist_m.at(*ptbErasedIt), yDist_m.back());
1919  std::swap(pyDist_m.at(*ptbErasedIt), pyDist_m.back());
1920  std::swap(tOrZDist_m.at(*ptbErasedIt), tOrZDist_m.back());
1921  std::swap(pzDist_m.at(*ptbErasedIt), pzDist_m.back());
1922  if (additionalRNs_m.size() == xDist_m.size()) {
1923  std::swap(additionalRNs_m.at(*ptbErasedIt), additionalRNs_m.back());
1924 
1925  additionalRNs_m.pop_back();
1926  }
1927 
1928  xDist_m.pop_back();
1929  pxDist_m.pop_back();
1930  yDist_m.pop_back();
1931  pyDist_m.pop_back();
1932  tOrZDist_m.pop_back();
1933  pzDist_m.pop_back();
1934 
1935  }
1936 
1937  /*
1938  * Set energy bin to emitted if all particles in the bin have been emitted.
1939  * However, be careful with the last energy bin. We cannot emit it until all
1940  * of the particles have been accounted for. So when on the last bin, keep it
1941  * open for the rest of the beam->
1942  */
1943  currentSampleBin_m++;
1944  if (currentSampleBin_m == numberOfSampleBins_m) {
1945 
1946  INFOMSG(level3 << "* Bin number: "
1947  << currentEnergyBin_m
1948  << " has emitted all particles (new emit)." << endl);
1949  currentSampleBin_m = 0;
1950  currentEnergyBin_m++;
1951 
1952  }
1953 
1954  /*
1955  * Set beam to emitted. Make sure temporary storage is cleared.
1956  */
1957  if (currentEnergyBin_m > numberOfEnergyBins_m || tOrZDist_m.empty()) {
1958  emitting_m = false;
1959 
1960  xDist_m.clear();
1961  pxDist_m.clear();
1962  yDist_m.clear();
1963  pyDist_m.clear();
1964  tOrZDist_m.clear();
1965  pzDist_m.clear();
1966 
1967  currentEnergyBin_m = numberOfEnergyBins_m;
1968  }
1969 
1970  }
1971  currentEmissionTime_m += beam->getdT();
1972 
1973  // Write emitted particles to file.
1974  writeOutFileEmission();
1975 
1976  size_t currentEmittedParticles = numberOfEmittedParticles - oldNumberOfEmittedParticles;
1977  reduce(currentEmittedParticles, currentEmittedParticles, OpAddAssign());
1978  totalNumberEmittedParticles_m += currentEmittedParticles;
1979 
1980  return currentEmittedParticles;
1981 
1982 }
1983 
1984 void Distribution::eraseXDist() {
1985  xDist_m.erase(xDist_m.begin(), xDist_m.end());
1986 }
1987 
1988 void Distribution::eraseBGxDist() {
1989  pxDist_m.erase(pxDist_m.begin(), pxDist_m.end());
1990 }
1991 
1992 void Distribution::eraseYDist() {
1993  yDist_m.erase(yDist_m.begin(), yDist_m.end());
1994 }
1995 
1996 void Distribution::eraseBGyDist() {
1997  pyDist_m.erase(pyDist_m.begin(), pyDist_m.end());
1998 }
1999 
2000 void Distribution::eraseTOrZDist() {
2001  tOrZDist_m.erase(tOrZDist_m.begin(), tOrZDist_m.end());
2002 }
2003 
2004 void Distribution::eraseBGzDist() {
2005  pzDist_m.erase(pzDist_m.begin(), pzDist_m.end());
2006 }
2007 
2008 void Distribution::fillEBinHistogram() {
2009  double upper = 0.0;
2010  double lower = 0.0;
2011  gsl_histogram_get_range(energyBinHist_m,
2012  gsl_histogram_bins(energyBinHist_m) - 1,
2013  &lower, &upper);
2014  const size_t numberOfParticles = tOrZDist_m.size();
2015  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2016 
2017  double &tOrZ = tOrZDist_m.at(partIndex);
2018 
2019  if (gsl_histogram_increment(energyBinHist_m, tOrZ) == GSL_EDOM) {
2020  gsl_histogram_increment(energyBinHist_m, 0.5 * (lower + upper));
2021  }
2022  }
2023 
2024  reduce(energyBinHist_m->bin, energyBinHist_m->bin + energyBinHist_m->n, energyBinHist_m->bin, OpAddAssign());
2025 }
2026 
2027 void Distribution::fillParticleBins() {
2028 
2029  for (size_t particleIndex = 0; particleIndex < tOrZDist_m.size(); particleIndex++) {
2030 
2031  std::vector<double> partCoord;
2032 
2033  partCoord.push_back(xDist_m.at(particleIndex));
2034  partCoord.push_back(yDist_m.at(particleIndex));
2035  partCoord.push_back(tOrZDist_m.at(particleIndex));
2036  partCoord.push_back(pxDist_m.at(particleIndex));
2037  partCoord.push_back(pyDist_m.at(particleIndex));
2038  partCoord.push_back(pzDist_m.at(particleIndex));
2039  partCoord.push_back(0.0);
2040 
2041  energyBins_m->fill(partCoord);
2042 
2043  }
2044 
2045  energyBins_m->sortArray();
2046 }
2047 
2048 size_t Distribution::findEBin(double tOrZ) {
2049 
2050  if (tOrZ <= gsl_histogram_min(energyBinHist_m)) {
2051  return 0;
2052  } else if (tOrZ >= gsl_histogram_max(energyBinHist_m)) {
2053  return numberOfEnergyBins_m - 1;
2054  } else {
2055  size_t binNumber;
2056  gsl_histogram_find(energyBinHist_m, tOrZ, &binNumber);
2057  return binNumber / numberOfSampleBins_m;
2058  }
2059 }
2060 
2061 void Distribution::generateAstraFlattopT(size_t numberOfParticles) {
2062 
2063  /*
2064  * Legacy function to support the ASTRAFLATTOPOTH
2065  * distribution type.
2066  */
2067  checkEmissionParameters();
2068 
2069  gsl_qrng *quasiRandGen = gsl_qrng_alloc(gsl_qrng_halton, 2);
2070 
2071  int numberOfSampleBins
2072  = std::abs(static_cast<int> (Attributes::getReal(itsAttr[Attrib::Legacy::Distribution::SBIN])));
2073  int numberOfEnergyBins
2074  = std::abs(static_cast<int> (Attributes::getReal(itsAttr[Attrib::Distribution::NBIN])));
2075 
2076  int binTotal = numberOfSampleBins * numberOfEnergyBins;
2077 
2078  double *distributionTable = new double[binTotal + 1];
2079 
2080  double a = tPulseLengthFWHM_m / 2.;
2081  double sig = tRise_m / 2.;
2082  double inv_erf08 = 0.906193802436823; // erfinv(0.8)
2083  double sqr2 = sqrt(2.0);
2084  double t = a - sqr2 * sig * inv_erf08;
2085  double tmps = sig;
2086  double tmpt = t;
2087 
2088  for (int i = 0; i < 10; ++ i) {
2089  sig = (t + tRise_m - a) / (sqr2 * inv_erf08);
2090  t = a - 0.5 * sqr2 * (sig + tmps) * inv_erf08;
2091  sig = (0.5 * (t + tmpt) + tRise_m - a) / (sqr2 * inv_erf08);
2092  tmps = sig;
2093  tmpt = t;
2094  }
2095 
2096  /*
2097  * Emission time is set here during distribution particle creation only for
2098  * the ASTRAFLATTOPTH distribution type.
2099  */
2100  tEmission_m = tPulseLengthFWHM_m + 10. * sig;
2101  tBin_m = tEmission_m / numberOfEnergyBins;
2102 
2103  double lo = -tBin_m / 2.0 * numberOfEnergyBins;
2104  double hi = tBin_m / 2.0 * numberOfEnergyBins;
2105  double dx = tBin_m / numberOfSampleBins;
2106  double x = lo;
2107  double tot = 0;
2108  double weight = 2.0;
2109 
2110  // sample the function that describes the histogram of the requested distribution
2111  for (int i = 0; i < binTotal + 1; ++ i, x += dx, weight = 6. - weight) {
2112  distributionTable[i] = gsl_sf_erf((x + a) / (sqrt(2) * sig)) - gsl_sf_erf((x - a) / (sqrt(2) * sig));
2113  tot += distributionTable[i] * weight;
2114  }
2115  tot -= distributionTable[binTotal] * (5. - weight);
2116  tot -= distributionTable[0];
2117 
2118  int saveProcessor = -1;
2119  const int myNode = Ippl::myNode();
2120  const int numNodes = Ippl::getNodes();
2121  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2122  double tCoord = 0.0;
2123 
2124  int effectiveNumParticles = 0;
2125  int largestBin = 0;
2126  std::vector<int> numParticlesInBin(numberOfEnergyBins,0);
2127  for (int k = 0; k < numberOfEnergyBins; ++ k) {
2128  double loc_fraction = -distributionTable[numberOfSampleBins * k] / tot;
2129 
2130  weight = 2.0;
2131  for (int i = numberOfSampleBins * k; i < numberOfSampleBins * (k + 1) + 1;
2132  ++ i, weight = 6. - weight) {
2133  loc_fraction += distributionTable[i] * weight / tot;
2134  }
2135  loc_fraction -= distributionTable[numberOfSampleBins * (k + 1)]
2136  * (5. - weight) / tot;
2137  numParticlesInBin[k] = static_cast<int>(std::floor(loc_fraction * numberOfParticles + 0.5));
2138  effectiveNumParticles += numParticlesInBin[k];
2139  if (numParticlesInBin[k] > numParticlesInBin[largestBin]) largestBin = k;
2140  }
2141 
2142  numParticlesInBin[largestBin] += (numberOfParticles - effectiveNumParticles);
2143 
2144  for (int k = 0; k < numberOfEnergyBins; ++ k) {
2145  gsl_ran_discrete_t *table
2146  = gsl_ran_discrete_preproc(numberOfSampleBins,
2147  &(distributionTable[numberOfSampleBins * k]));
2148 
2149  for (int i = 0; i < numParticlesInBin[k]; i++) {
2150  double xx[2];
2151  gsl_qrng_get(quasiRandGen, xx);
2152  tCoord = hi * (xx[1] + static_cast<int>(gsl_ran_discrete(randGen_m, table))
2153  - binTotal / 2 + k * numberOfSampleBins) / (binTotal / 2);
2154 
2155  saveProcessor++;
2156  if (saveProcessor >= numNodes)
2157  saveProcessor = 0;
2158 
2159  if (scalable || myNode == saveProcessor) {
2160  tOrZDist_m.push_back(tCoord);
2161  pzDist_m.push_back(0.0);
2162  }
2163  }
2164  gsl_ran_discrete_free(table);
2165  }
2166 
2167  gsl_qrng_free(quasiRandGen);
2168  delete [] distributionTable;
2169 
2170 }
2171 
2172 void Distribution::generateBinomial(size_t numberOfParticles) {
2173 
2193  // Calculate emittance and Twiss parameters.
2194  Vector_t emittance;
2195  Vector_t alpha, beta, gamma;
2196  for (unsigned int index = 0; index < 3; index++) {
2197  emittance(index) = sigmaR_m[index] * sigmaP_m[index]
2198  * cos(asin(correlationMatrix_m(2 * index + 1, 2 * index)));
2199 
2200  if (std::abs(emittance(index)) > std::numeric_limits<double>::epsilon()) {
2201  beta(index) = pow(sigmaR_m[index], 2.0) / emittance(index);
2202  gamma(index) = pow(sigmaP_m[index], 2.0) / emittance(index);
2203  } else {
2204  beta(index) = sqrt(std::numeric_limits<double>::max());
2205  gamma(index) = sqrt(std::numeric_limits<double>::max());
2206  }
2207  alpha(index) = -correlationMatrix_m(2 * index + 1, 2 * index)
2208  * sqrt(beta(index) * gamma(index));
2209  }
2210 
2211  Vector_t M, PM, L, PL, X, PX;
2212  Vector_t CHI, COSCHI, SINCHI(0.0);
2213  Vector_t AMI;
2214  Vektor<BinomialBehaviorSplitter*, 3> splitter;
2215  for (unsigned int index = 0; index < 3; index++) {
2216  // gamma(index) *= cutoffR_m[index];
2217  // beta(index) *= cutoffP_m[index];
2218  COSCHI[index] = sqrt(1.0 / (1.0 + pow(alpha(index), 2.0)));
2219  SINCHI[index] = -alpha(index) * COSCHI[index];
2220  CHI[index] = atan2(SINCHI[index], COSCHI[index]);
2221  AMI[index] = 1.0 / mBinomial_m[index];
2222  M[index] = sqrt(emittance(index) * beta(index));
2223  PM[index] = sqrt(emittance(index) * gamma(index));
2224  L[index] = sqrt((mBinomial_m[index] + 1.0) / 2.0) * M[index];
2225  PL[index] = sqrt((mBinomial_m[index] + 1.0) / 2.0) * PM[index];
2226 
2227  if (mBinomial_m[index] < 10000) {
2228  X[index] = L[index];
2229  PX[index] = PL[index];
2230  splitter[index] = new MDependentBehavior(mBinomial_m[index]);
2231  } else {
2232  X[index] = M[index];
2233  PX[index] = PM[index];
2234  splitter[index] = new GaussianLikeBehavior();
2235  }
2236  }
2237 
2238  int saveProcessor = -1;
2239  const int myNode = Ippl::myNode();
2240  const int numNodes = Ippl::getNodes();
2241  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2242 
2243  double temp = 1.0 - std::pow(correlationMatrix_m(1, 0), 2.0);
2244  const double tempa = copysign(sqrt(std::abs(temp)), temp);
2245  const double l32 = (correlationMatrix_m(4, 1) -
2246  correlationMatrix_m(1, 0) * correlationMatrix_m(4, 0)) / tempa;
2247  temp = 1 - std::pow(correlationMatrix_m(4, 0), 2.0) - l32 * l32;
2248  const double l33 = copysign(sqrt(std::abs(temp)), temp);
2249 
2250  const double l42 = (correlationMatrix_m(5, 1) -
2251  correlationMatrix_m(1, 0) * correlationMatrix_m(5, 0)) / tempa;
2252  const double l43 = (correlationMatrix_m(5, 4) -
2253  correlationMatrix_m(4, 0) * correlationMatrix_m(5, 0) - l42 * l32) / l33;
2254  temp = 1 - std::pow(correlationMatrix_m(5, 0), 2.0) - l42 * l42 - l43 * l43;
2255  const double l44 = copysign(sqrt(std::abs(temp)), temp);
2256 
2257  Vector_t x = Vector_t(0.0);
2258  Vector_t p = Vector_t(0.0);
2259  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2260 
2261  double A = 0.0;
2262  double AL = 0.0;
2263  double Ux = 0.0, U = 0.0;
2264  double Vx = 0.0, V = 0.0;
2265 
2266  A = splitter[0]->get(gsl_rng_uniform(randGen_m));
2267  AL = Physics::two_pi * gsl_rng_uniform(randGen_m);
2268  Ux = A * cos(AL);
2269  Vx = A * sin(AL);
2270  x[0] = X[0] * Ux;
2271  p[0] = PX[0] * (Ux * SINCHI[0] + Vx * COSCHI[0]);
2272 
2273  A = splitter[1]->get(gsl_rng_uniform(randGen_m));
2274  AL = Physics::two_pi * gsl_rng_uniform(randGen_m);
2275  U = A * cos(AL);
2276  V = A * sin(AL);
2277  x[1] = X[1] * U;
2278  p[1] = PX[1] * (U * SINCHI[1] + V * COSCHI[1]);
2279 
2280  A = splitter[2]->get(gsl_rng_uniform(randGen_m));
2281  AL = Physics::two_pi * gsl_rng_uniform(randGen_m);
2282  U = A * cos(AL);
2283  V = A * sin(AL);
2284  x[2] = X[2] * (Ux * correlationMatrix_m(4, 0) + Vx * l32 + U * l33);
2285  p[2] = PX[2] * (Ux * correlationMatrix_m(5, 0) + Vx * l42 + U * l43 + V * l44);
2286 
2287  // Save to each processor in turn.
2288  saveProcessor++;
2289  if (saveProcessor >= numNodes)
2290  saveProcessor = 0;
2291 
2292  if (scalable || myNode == saveProcessor) {
2293  xDist_m.push_back(x[0]);
2294  pxDist_m.push_back(p[0]);
2295  yDist_m.push_back(x[1]);
2296  pyDist_m.push_back(p[1]);
2297  tOrZDist_m.push_back(x[2]);
2298  pzDist_m.push_back(avrgpz_m + p[2]);
2299  }
2300  }
2301 
2302  for (unsigned int index = 0; index < 3; index++) {
2303  delete splitter[index];
2304  }
2305 }
2306 
2307 void Distribution::generateFlattopLaserProfile(size_t numberOfParticles) {
2308 
2309  int saveProcessor = -1;
2310  const int myNode = Ippl::myNode();
2311  const int numNodes = Ippl::getNodes();
2312  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2313 
2314  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2315 
2316  double x = 0.0;
2317  double px = 0.0;
2318  double y = 0.0;
2319  double py = 0.0;
2320 
2321  laserProfile_m->getXY(x, y);
2322 
2323  // Save to each processor in turn.
2324  saveProcessor++;
2325  if (saveProcessor >= numNodes)
2326  saveProcessor = 0;
2327 
2328  if (scalable || myNode == saveProcessor) {
2329  xDist_m.push_back(x * sigmaR_m[0]);
2330  pxDist_m.push_back(px);
2331  yDist_m.push_back(y * sigmaR_m[1]);
2332  pyDist_m.push_back(py);
2333  }
2334  }
2335 
2336  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH)
2337  generateAstraFlattopT(numberOfParticles);
2338  else
2339  generateLongFlattopT(numberOfParticles);
2340 
2341 }
2342 
2343 void Distribution::generateFlattopT(size_t numberOfParticles) {
2344 
2345  gsl_qrng *quasiRandGen2D = selectRandomGenerator(Options::rngtype,2);
2346 
2347  int saveProcessor = -1;
2348  const int myNode = Ippl::myNode();
2349  const int numNodes = Ippl::getNodes();
2350  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2351  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2352 
2353  double x = 0.0;
2354  double px = 0.0;
2355  double y = 0.0;
2356  double py = 0.0;
2357 
2358  bool allow = false;
2359  while (!allow) {
2360 
2361  if (quasiRandGen2D != NULL) {
2362  double randNums[2];
2363  gsl_qrng_get(quasiRandGen2D, randNums);
2364  x = -1.0 + 2.0 * randNums[0];
2365  y = -1.0 + 2.0 * randNums[1];
2366  } else {
2367  x = -1.0 + 2.0 * gsl_rng_uniform(randGen_m);
2368  y = -1.0 + 2.0 * gsl_rng_uniform(randGen_m);
2369  }
2370 
2371  allow = (pow(x, 2.0) + pow(y, 2.0) <= 1.0);
2372 
2373  }
2374  x *= sigmaR_m[0];
2375  y *= sigmaR_m[1];
2376 
2377  // Save to each processor in turn.
2378  saveProcessor++;
2379  if (saveProcessor >= numNodes)
2380  saveProcessor = 0;
2381 
2382  if (scalable || myNode == saveProcessor) {
2383  xDist_m.push_back(x);
2384  pxDist_m.push_back(px);
2385  yDist_m.push_back(y);
2386  pyDist_m.push_back(py);
2387  }
2388 
2389  }
2390 
2391  gsl_qrng_free(quasiRandGen2D);
2392 
2393  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH)
2394  generateAstraFlattopT(numberOfParticles);
2395  else
2396  generateLongFlattopT(numberOfParticles);
2397 
2398 }
2399 
2400 void Distribution::generateFlattopZ(size_t numberOfParticles) {
2401 
2402  gsl_qrng *quasiRandGen1D = selectRandomGenerator(Options::rngtype,1);
2403  gsl_qrng *quasiRandGen2D = selectRandomGenerator(Options::rngtype,2);
2404 
2405  int saveProcessor = -1;
2406  const int myNode = Ippl::myNode();
2407  const int numNodes = Ippl::getNodes();
2408  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2409 
2410  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2411 
2412  double x = 0.0;
2413  double px = 0.0;
2414  double y = 0.0;
2415  double py = 0.0;
2416  double z = 0.0;
2417  double pz = 0.0;
2418 
2419  bool allow = false;
2420  while (!allow) {
2421 
2422  if (quasiRandGen2D != NULL) {
2423  double randNums[2];
2424  gsl_qrng_get(quasiRandGen2D, randNums);
2425  x = -1.0 + 2.0 * randNums[0];
2426  y = -1.0 + 2.0 * randNums[1];
2427  } else {
2428  x = -1.0 + 2.0 * gsl_rng_uniform(randGen_m);
2429  y = -1.0 + 2.0 * gsl_rng_uniform(randGen_m);
2430  }
2431 
2432  allow = (pow(x, 2.0) + pow(y, 2.0) <= 1.0);
2433 
2434  }
2435  x *= sigmaR_m[0];
2436  y *= sigmaR_m[1];
2437 
2438  if (quasiRandGen1D != NULL)
2439  gsl_qrng_get(quasiRandGen1D, &z);
2440  else
2441  z = gsl_rng_uniform(randGen_m);
2442 
2443  z = (z - 0.5) * sigmaR_m[2];
2444 
2445  // Save to each processor in turn.
2446  saveProcessor++;
2447  if (saveProcessor >= numNodes)
2448  saveProcessor = 0;
2449 
2450  if (scalable || myNode == saveProcessor) {
2451  xDist_m.push_back(x);
2452  pxDist_m.push_back(px);
2453  yDist_m.push_back(y);
2454  pyDist_m.push_back(py);
2455  tOrZDist_m.push_back(z);
2456  pzDist_m.push_back(avrgpz_m + pz);
2457  }
2458  }
2459 
2460  gsl_qrng_free(quasiRandGen1D);
2461  gsl_qrng_free(quasiRandGen2D);
2462 }
2463 
2464 void Distribution::generateGaussZ(size_t numberOfParticles) {
2465 
2466  gsl_matrix * corMat = gsl_matrix_alloc (6, 6);
2467  gsl_vector * rx = gsl_vector_alloc(6);
2468  gsl_vector * ry = gsl_vector_alloc(6);
2469 
2470  for (unsigned int i = 0; i < 6; ++ i) {
2471  gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
2472  for (unsigned int j = 0; j < i; ++ j) {
2473  gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
2474  gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
2475  }
2476  }
2477 
2478 #define DISTDBG1
2479 #ifdef DISTDBG1
2480  *gmsg << "* m before gsl_linalg_cholesky_decomp" << endl;
2481  for (int i = 0; i < 6; i++) {
2482  for (int j = 0; j < 6; j++) {
2483  if (j==0)
2484  *gmsg << "* r(" << std::setprecision(1) << i << ", : ) = "
2485  << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
2486  else
2487  *gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
2488  }
2489  *gmsg << " \\\\" << endl;
2490  }
2491 #endif
2492 /*
2493  //Sets the GSL error handler off, exception will be handled internally with a renormalization method
2494  gsl_set_error_handler_off();
2495 */
2496  int errcode = gsl_linalg_cholesky_decomp(corMat);
2497 /*
2498  double rn = 1e-12;
2499 
2500  while (errcode == GSL_EDOM) {
2501 
2502  // Resets the correlation matrix
2503  for (unsigned int i = 0; i < 6; ++ i) {
2504  gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
2505  for (unsigned int j = 0; j < i; ++ j) {
2506  gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
2507  gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
2508  }
2509  }
2510  // Applying a renormalization method corMat = corMat + rn*Unitymatrix
2511  // This is the renormalization
2512  for(int i = 0; i < 6; i++){
2513  double corMati = gsl_matrix_get(corMat, i , i);
2514  gsl_matrix_set(corMat, i, i, corMati + rn);
2515  }
2516  //Just to be sure that the renormalization worked
2517  errcode = gsl_linalg_cholesky_decomp(corMat);
2518  if(errcode != GSL_EDOM) *gmsg << "* WARNING: Correlation matrix had to be renormalized"<<endl;
2519  else rn *= 10;
2520  }
2521  //Sets again the standard GSL error handler on
2522  gsl_set_error_handler(NULL);
2523 */
2524  //Just to be sure
2525  if (errcode == GSL_EDOM) {
2526  throw OpalException("Distribution::GenerateGaussZ",
2527  "gsl_linalg_cholesky_decomp failed");
2528  }
2529  // so make the upper part zero.
2530  for (int i = 0; i < 5; ++ i) {
2531  for (int j = i+1; j < 6 ; ++ j) {
2532  gsl_matrix_set (corMat, i, j, 0.0);
2533  }
2534  }
2535 #define DISTDBG2
2536 #ifdef DISTDBG2
2537  *gmsg << "* m after gsl_linalg_cholesky_decomp" << endl;
2538  for (int i = 0; i < 6; i++) {
2539  for (int j = 0; j < 6; j++) {
2540  if (j==0)
2541  *gmsg << "* r(" << std::setprecision(1) << i << ", : ) = "
2542  << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
2543  else
2544  *gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
2545  }
2546  *gmsg << " \\\\" << endl;
2547  }
2548 #endif
2549 
2550  int saveProcessor = -1;
2551  const int myNode = Ippl::myNode();
2552  const int numNodes = Ippl::getNodes();
2553  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2554 
2555  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2556  bool allow = false;
2557 
2558  double x = 0.0;
2559  double px = 0.0;
2560  double y = 0.0;
2561  double py = 0.0;
2562  double z = 0.0;
2563  double pz = 0.0;
2564 
2565  while (!allow) {
2566  gsl_vector_set(rx, 0, gsl_ran_gaussian(randGen_m, 1.0));
2567  gsl_vector_set(rx, 1, gsl_ran_gaussian(randGen_m, 1.0));
2568  gsl_vector_set(rx, 2, gsl_ran_gaussian(randGen_m, 1.0));
2569  gsl_vector_set(rx, 3, gsl_ran_gaussian(randGen_m, 1.0));
2570  gsl_vector_set(rx, 4, gsl_ran_gaussian(randGen_m, 1.0));
2571  gsl_vector_set(rx, 5, gsl_ran_gaussian(randGen_m, 1.0));
2572 
2573  gsl_blas_dgemv(CblasNoTrans, 1.0, corMat, rx, 0.0, ry);
2574 
2575  x = gsl_vector_get(ry, 0);
2576  px = gsl_vector_get(ry, 1);
2577  y = gsl_vector_get(ry, 2);
2578  py = gsl_vector_get(ry, 3);
2579  z = gsl_vector_get(ry, 4);
2580  pz = gsl_vector_get(ry, 5);
2581 
2582  bool xAndYOk = false;
2583  if (cutoffR_m[0] < SMALLESTCUTOFF && cutoffR_m[1] < SMALLESTCUTOFF)
2584  xAndYOk = true;
2585  else if (cutoffR_m[0] < SMALLESTCUTOFF)
2586  xAndYOk = (std::abs(y) <= cutoffR_m[1]);
2587  else if (cutoffR_m[1] < SMALLESTCUTOFF)
2588  xAndYOk = (std::abs(x) <= cutoffR_m[0]);
2589  else
2590  xAndYOk = (pow(x / cutoffR_m[0], 2.0) + pow(y / cutoffR_m[1], 2.0) <= 1.0);
2591 
2592  bool pxAndPyOk = false;
2593  if (cutoffP_m[0] < SMALLESTCUTOFF && cutoffP_m[1] < SMALLESTCUTOFF)
2594  pxAndPyOk = true;
2595  else if (cutoffP_m[0] < SMALLESTCUTOFF)
2596  pxAndPyOk = (std::abs(py) <= cutoffP_m[1]);
2597  else if (cutoffP_m[1] < SMALLESTCUTOFF)
2598  pxAndPyOk = (std::abs(px) <= cutoffP_m[0]);
2599  else
2600  pxAndPyOk = (pow(px / cutoffP_m[0], 2.0) + pow(py / cutoffP_m[1], 2.0) <= 1.0);
2601 
2602  allow = (xAndYOk && pxAndPyOk && std::abs(z) < cutoffR_m[2] && std::abs(pz) < cutoffP_m[2]);
2603  }
2604 
2605  x *= sigmaR_m[0];
2606  y *= sigmaR_m[1];
2607  z *= sigmaR_m[2];
2608  px *= sigmaP_m[0];
2609  py *= sigmaP_m[1];
2610  pz *= sigmaP_m[2];
2611 
2612  // Save to each processor in turn.
2613  saveProcessor++;
2614  if (saveProcessor >= numNodes)
2615  saveProcessor = 0;
2616 
2617  if (scalable || myNode == saveProcessor) {
2618  xDist_m.push_back(x);
2619  pxDist_m.push_back(px);
2620  yDist_m.push_back(y);
2621  pyDist_m.push_back(py);
2622  tOrZDist_m.push_back(z);
2623  pzDist_m.push_back(avrgpz_m + pz);
2624  }
2625  }
2626 
2627  gsl_vector_free(rx);
2628  gsl_vector_free(ry);
2629  gsl_matrix_free(corMat);
2630 }
2631 
2632 void Distribution::generateLongFlattopT(size_t numberOfParticles) {
2633 
2634  double flattopTime = tPulseLengthFWHM_m
2635  - sqrt(2.0 * log(2.0)) * (sigmaTRise_m + sigmaTFall_m);
2636 
2637  if (flattopTime < 0.0)
2638  flattopTime = 0.0;
2639 
2640  double normalizedFlankArea = 0.5 * sqrt(2.0 * Physics::pi) * gsl_sf_erf(cutoffR_m[2] / sqrt(2.0));
2641  double distArea = flattopTime
2642  + (sigmaTRise_m + sigmaTFall_m) * normalizedFlankArea;
2643 
2644  // Find number of particles in rise, fall and flat top.
2645  size_t numRise = numberOfParticles * sigmaTRise_m * normalizedFlankArea / distArea;
2646  size_t numFall = numberOfParticles * sigmaTFall_m * normalizedFlankArea / distArea;
2647  size_t numFlat = numberOfParticles - numRise - numFall;
2648 
2649  // Generate particles in tail.
2650  int saveProcessor = -1;
2651  const int myNode = Ippl::myNode();
2652  const int numNodes = Ippl::getNodes();
2653  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2654 
2655  for (size_t partIndex = 0; partIndex < numFall; partIndex++) {
2656 
2657  double t = 0.0;
2658  double pz = 0.0;
2659 
2660  bool allow = false;
2661  while (!allow) {
2662  t = gsl_ran_gaussian_tail(randGen_m, 0, sigmaTFall_m);
2663  if (t <= sigmaTFall_m * cutoffR_m[2]) {
2664  t = -t + sigmaTFall_m * cutoffR_m[2];
2665  allow = true;
2666  }
2667  }
2668 
2669  // Save to each processor in turn.
2670  saveProcessor++;
2671  if (saveProcessor >= numNodes)
2672  saveProcessor = 0;
2673 
2674  if (scalable || myNode == saveProcessor) {
2675  tOrZDist_m.push_back(t);
2676  pzDist_m.push_back(pz);
2677  }
2678  }
2679 
2680  /*
2681  * Generate particles in flat top. The flat top can also have sinusoidal
2682  * modulations.
2683  */
2684  double modulationAmp = Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCAMPLITUDE]) / 100.0;
2685  if (modulationAmp > 1.0)
2686  modulationAmp = 1.0;
2687  double numModulationPeriods
2688  = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCPERIODS]));
2689  double modulationPeriod = 0.0;
2690  if (numModulationPeriods != 0.0)
2691  modulationPeriod = flattopTime / numModulationPeriods;
2692 
2693  gsl_qrng *quasiRandGen1D = selectRandomGenerator(Options::rngtype,1);
2694  gsl_qrng *quasiRandGen2D = selectRandomGenerator(Options::rngtype,2);
2695 
2696  for (size_t partIndex = 0; partIndex < numFlat; partIndex++) {
2697 
2698  double t = 0.0;
2699  double pz = 0.0;
2700 
2701  if (modulationAmp == 0.0 || numModulationPeriods == 0.0) {
2702 
2703  if (quasiRandGen1D != NULL)
2704  gsl_qrng_get(quasiRandGen1D, &t);
2705  else
2706  t = gsl_rng_uniform(randGen_m);
2707 
2708  t = flattopTime * t;
2709 
2710  } else {
2711 
2712  double funcAmp = 0.0;
2713  bool allow = false;
2714  while (!allow) {
2715  if (quasiRandGen2D != NULL) {
2716  double randNums[2] = {0.0, 0.0};
2717  gsl_qrng_get(quasiRandGen2D, randNums);
2718  t = randNums[0] * flattopTime;
2719  funcAmp = randNums[1];
2720  } else {
2721  t = gsl_rng_uniform(randGen_m)*flattopTime;
2722  funcAmp = gsl_rng_uniform(randGen_m);
2723  }
2724 
2725  double funcValue = (1.0 + modulationAmp
2726  * sin(Physics::two_pi * t / modulationPeriod))
2727  / (1.0 + std::abs(modulationAmp));
2728 
2729  allow = (funcAmp <= funcValue);
2730 
2731  }
2732  }
2733  // Shift the uniform part of distribution behind the fall
2734  t += sigmaTFall_m * cutoffR_m[2];
2735 
2736  // Save to each processor in turn.
2737  saveProcessor++;
2738  if (saveProcessor >= numNodes)
2739  saveProcessor = 0;
2740 
2741  if (scalable || myNode == saveProcessor) {
2742  tOrZDist_m.push_back(t);
2743  pzDist_m.push_back(pz);
2744  }
2745  }
2746 
2747  // Generate particles in rise.
2748  for (size_t partIndex = 0; partIndex < numRise; partIndex++) {
2749 
2750  double t = 0.0;
2751  double pz = 0.0;
2752 
2753  bool allow = false;
2754  while (!allow) {
2755  t = gsl_ran_gaussian_tail(randGen_m, 0, sigmaTRise_m);
2756  if (t <= sigmaTRise_m * cutoffR_m[2]) {
2757  t += sigmaTFall_m * cutoffR_m[2] + flattopTime;
2758  allow = true;
2759  }
2760  }
2761 
2762  // Save to each processor in turn.
2763  saveProcessor++;
2764  if (saveProcessor >= numNodes)
2765  saveProcessor = 0;
2766 
2767  if (scalable || myNode == saveProcessor) {
2768  tOrZDist_m.push_back(t);
2769  pzDist_m.push_back(pz);
2770  }
2771  }
2772 
2773  gsl_qrng_free(quasiRandGen1D);
2774  gsl_qrng_free(quasiRandGen2D);
2775 }
2776 
2777 void Distribution::generateTransverseGauss(size_t numberOfParticles) {
2778 
2779  // Generate coordinates.
2780  gsl_matrix * corMat = gsl_matrix_alloc (4, 4);
2781  gsl_vector * rx = gsl_vector_alloc(4);
2782  gsl_vector * ry = gsl_vector_alloc(4);
2783 
2784  for (unsigned int i = 0; i < 4; ++ i) {
2785  gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
2786  for (unsigned int j = 0; j < i; ++ j) {
2787  gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
2788  gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
2789  }
2790  }
2791 
2792  int errcode = gsl_linalg_cholesky_decomp(corMat);
2793 
2794  if (errcode == GSL_EDOM) {
2795  throw OpalException("Distribution::GenerateTransverseGauss",
2796  "gsl_linalg_cholesky_decomp failed");
2797  }
2798 
2799  for (int i = 0; i < 3; ++ i) {
2800  for (int j = i+1; j < 4 ; ++ j) {
2801  gsl_matrix_set (corMat, i, j, 0.0);
2802  }
2803  }
2804 
2805  int saveProcessor = -1;
2806  const int myNode = Ippl::myNode();
2807  const int numNodes = Ippl::getNodes();
2808  const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
2809 
2810  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2811 
2812  double x = 0.0;
2813  double px = 0.0;
2814  double y = 0.0;
2815  double py = 0.0;
2816 
2817  bool allow = false;
2818  while (!allow) {
2819 
2820  gsl_vector_set(rx, 0, gsl_ran_gaussian (randGen_m,1.0));
2821  gsl_vector_set(rx, 1, gsl_ran_gaussian (randGen_m,1.0));
2822  gsl_vector_set(rx, 2, gsl_ran_gaussian (randGen_m,1.0));
2823  gsl_vector_set(rx, 3, gsl_ran_gaussian (randGen_m,1.0));
2824 
2825  gsl_blas_dgemv(CblasNoTrans, 1.0, corMat, rx, 0.0, ry);
2826  x = gsl_vector_get(ry, 0);
2827  px = gsl_vector_get(ry, 1);
2828  y = gsl_vector_get(ry, 2);
2829  py = gsl_vector_get(ry, 3);
2830 
2831  bool xAndYOk = false;
2832  if (cutoffR_m[0] < SMALLESTCUTOFF && cutoffR_m[1] < SMALLESTCUTOFF)
2833  xAndYOk = true;
2834  else if (cutoffR_m[0] < SMALLESTCUTOFF)
2835  xAndYOk = (std::abs(y) <= cutoffR_m[1]);
2836  else if (cutoffR_m[1] < SMALLESTCUTOFF)
2837  xAndYOk = (std::abs(x) <= cutoffR_m[0]);
2838  else
2839  xAndYOk = (pow(x / cutoffR_m[0], 2.0) + pow(y / cutoffR_m[1], 2.0) <= 1.0);
2840 
2841  bool pxAndPyOk = false;
2842  if (cutoffP_m[0] < SMALLESTCUTOFF && cutoffP_m[1] < SMALLESTCUTOFF)
2843  pxAndPyOk = true;
2844  else if (cutoffP_m[0] < SMALLESTCUTOFF)
2845  pxAndPyOk = (std::abs(py) <= cutoffP_m[1]);
2846  else if (cutoffP_m[1] < SMALLESTCUTOFF)
2847  pxAndPyOk = (std::abs(px) <= cutoffP_m[0]);
2848  else
2849  pxAndPyOk = (pow(px / cutoffP_m[0], 2.0) + pow(py / cutoffP_m[1], 2.0) <= 1.0);
2850 
2851  allow = (xAndYOk && pxAndPyOk);
2852 
2853  }
2854  x *= sigmaR_m[0];
2855  y *= sigmaR_m[1];
2856  px *= sigmaP_m[0];
2857  py *= sigmaP_m[1];
2858 
2859  // Save to each processor in turn.
2860  saveProcessor++;
2861  if (saveProcessor >= numNodes)
2862  saveProcessor = 0;
2863 
2864  if (scalable || myNode == saveProcessor) {
2865  xDist_m.push_back(x);
2866  pxDist_m.push_back(px);
2867  yDist_m.push_back(y);
2868  pyDist_m.push_back(py);
2869  }
2870  }
2871 
2872  gsl_matrix_free(corMat);
2873  gsl_vector_free(rx);
2874  gsl_vector_free(ry);
2875 }
2876 
2877 void Distribution::initializeBeam(PartBunchBase<double, 3> *beam) {
2878 
2879  /*
2880  * Set emission time, the current beam bunch number and
2881  * set the beam energy bin structure, if it has one.
2882  */
2883  beam->setTEmission(tEmission_m);
2884  beam->setNumBunch(1);
2885  if (numberOfEnergyBins_m > 0)
2886  beam->setEnergyBins(numberOfEnergyBins_m);
2887 }
2888 
2889 void Distribution::injectBeam(PartBunchBase<double, 3> *beam) {
2890 
2891  writeOutFileInjection();
2892 
2893  std::vector<double> id1 = Attributes::getRealArray(itsAttr[Attrib::Distribution::ID1]);
2894  std::vector<double> id2 = Attributes::getRealArray(itsAttr[Attrib::Distribution::ID2]);
2895 
2896  bool hasID1 = (id1.size() != 0);
2897  bool hasID2 = (id2.size() != 0);
2898 
2899  if (hasID1 || hasID2)
2900  *gmsg << "* Use special ID1 or ID2 particle in distribution" << endl;
2901 
2902 
2903  size_t numberOfParticles = tOrZDist_m.size();
2904  beam->create(numberOfParticles);
2905  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
2906 
2907  beam->R[partIndex] = Vector_t(xDist_m.at(partIndex),
2908  yDist_m.at(partIndex),
2909  tOrZDist_m.at(partIndex));
2910 
2911  beam->P[partIndex] = Vector_t(pxDist_m.at(partIndex),
2912  pyDist_m.at(partIndex),
2913  pzDist_m.at(partIndex));
2914 
2915  beam->Q[partIndex] = beam->getChargePerParticle();
2916  beam->Ef[partIndex] = Vector_t(0.0);
2917  beam->Bf[partIndex] = Vector_t(0.0);
2918  beam->PType[partIndex] = ParticleType::REGULAR;
2919  beam->TriID[partIndex] = 0;
2920 
2921  if (numberOfEnergyBins_m > 0) {
2922  size_t binNumber = findEBin(tOrZDist_m.at(partIndex));
2923  beam->Bin[partIndex] = binNumber;
2924  beam->iterateEmittedBin(binNumber);
2925  } else
2926  beam->Bin[partIndex] = 0;
2927 
2928  if (hasID1) {
2929  if (beam->ID[partIndex] == 1) {
2930  beam->R[partIndex] = Vector_t(id1[0],id1[1],id1[2]);
2931  beam->P[partIndex] = Vector_t(id1[3],id1[4],id1[5]);
2932  }
2933  }
2934 
2935  if (hasID2) {
2936  if (beam->ID[partIndex] == 2) {
2937  beam->R[partIndex] = Vector_t(id2[0],id2[1],id2[2]);
2938  beam->P[partIndex] = Vector_t(id2[3],id2[4],id2[5]);
2939  }
2940  }
2941  }
2942 
2943  xDist_m.clear();
2944  pxDist_m.clear();
2945  yDist_m.clear();
2946  pyDist_m.clear();
2947  tOrZDist_m.clear();
2948  pzDist_m.clear();
2949 
2950  beam->boundp();
2951  beam->calcEMean();
2952 }
2953 
2954 bool Distribution::getIfDistEmitting() {
2955  return emitting_m;
2956 }
2957 
2958 int Distribution::getLastEmittedEnergyBin() {
2959  return currentEnergyBin_m;
2960 }
2961 
2962 double Distribution::getMaxTOrZ() {
2963 
2964  std::vector<double>::iterator longIt = tOrZDist_m.begin();
2965  double maxTOrZ = *longIt;
2966  for (++longIt; longIt != tOrZDist_m.end(); ++longIt) {
2967  if (maxTOrZ < *longIt)
2968  maxTOrZ = *longIt;
2969  }
2970 
2971  reduce(maxTOrZ, maxTOrZ, OpMaxAssign());
2972 
2973  return maxTOrZ;
2974 }
2975 
2976 double Distribution::getMinTOrZ() {
2977 
2978  std::vector<double>::iterator longIt = tOrZDist_m.begin();
2979  double minTOrZ = *longIt;
2980  for (++longIt; longIt != tOrZDist_m.end(); ++longIt) {
2981  if (minTOrZ > *longIt)
2982  minTOrZ = *longIt;
2983  }
2984 
2985  reduce(minTOrZ, minTOrZ, OpMinAssign());
2986 
2987  return minTOrZ;
2988 }
2989 
2990 size_t Distribution::getNumberOfEmissionSteps() {
2991  return static_cast<size_t> (numberOfEnergyBins_m * numberOfSampleBins_m);
2992 }
2993 
2994 int Distribution::getNumberOfEnergyBins() {
2995  return numberOfEnergyBins_m;
2996 }
2997 
2998 double Distribution::getEmissionDeltaT() {
2999  return tBin_m / numberOfSampleBins_m;
3000 }
3001 
3002 double Distribution::getEnergyBinDeltaT() {
3003  return tBin_m;
3004 }
3005 
3006 double Distribution::getWeight() {
3007  return std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::WEIGHT]));
3008 }
3009 
3010 std::vector<double>& Distribution::getXDist() {
3011  return xDist_m;
3012 }
3013 
3014 std::vector<double>& Distribution::getBGxDist() {
3015  return pxDist_m;
3016 }
3017 
3018 std::vector<double>& Distribution::getYDist() {
3019  return yDist_m;
3020 }
3021 
3022 std::vector<double>& Distribution::getBGyDist() {
3023  return pyDist_m;
3024 }
3025 
3026 std::vector<double>& Distribution::getTOrZDist() {
3027  return tOrZDist_m;
3028 }
3029 
3030 std::vector<double>& Distribution::getBGzDist() {
3031  return pzDist_m;
3032 }
3033 
3034 void Distribution::printDist(Inform &os, size_t numberOfParticles) const {
3035 
3036  if (numberOfParticles > 0) {
3037  size_t np = numberOfParticles * (Options::cZero && !(distrTypeT_m == DistrTypeT::FROMFILE)? 2: 1);
3038  reduce(np, np, OpAddAssign());
3039  os << "* Number of particles: "
3040  << np
3041  << endl
3042  << "* " << endl;
3043  }
3044 
3045  switch (distrTypeT_m) {
3046 
3047  case DistrTypeT::FROMFILE:
3048  printDistFromFile(os);
3049  break;
3050  case DistrTypeT::GAUSS:
3051  printDistGauss(os);
3052  break;
3053  case DistrTypeT::BINOMIAL:
3054  printDistBinomial(os);
3055  break;
3056  case DistrTypeT::SURFACEEMISSION:
3057  printDistSurfEmission(os);
3058  break;
3059  case DistrTypeT::SURFACERANDCREATE:
3060  printDistSurfAndCreate(os);
3061  break;
3062  case DistrTypeT::FLATTOP:
3063  case DistrTypeT::GUNGAUSSFLATTOPTH:
3064  case DistrTypeT::ASTRAFLATTOPTH:
3065  printDistFlattop(os);
3066  break;
3067  case DistrTypeT::MATCHEDGAUSS:
3068  printDistMatchedGauss(os);
3069  break;
3070  default:
3071  INFOMSG("Distribution unknown." << endl;);
3072  break;
3073  }
3074 
3075 }
3076 
3077 void Distribution::printDistBinomial(Inform &os) const {
3078 
3079  os << "* Distribution type: BINOMIAL" << endl;
3080  os << "* " << endl;
3081  os << "* SIGMAX = " << sigmaR_m[0] << " [m]" << endl;
3082  os << "* SIGMAY = " << sigmaR_m[1] << " [m]" << endl;
3083  if (emitting_m)
3084  os << "* SIGMAT = " << sigmaR_m[2] << " [sec]" << endl;
3085  else
3086  os << "* SIGMAZ = " << sigmaR_m[2] << " [m]" << endl;
3087  os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
3088  os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
3089  os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
3090  os << "* MX = " << mBinomial_m[0] << endl;
3091  os << "* MY = " << mBinomial_m[1] << endl;
3092  if (emitting_m)
3093  os << "* MT = " << mBinomial_m[2] << endl;
3094  else
3095  os << "* MZ = " << mBinomial_m[2] << endl;
3096  os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
3097  os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
3098  os << "* CORRZ = " << correlationMatrix_m(5, 4) << endl;
3099  os << "* R61 = " << correlationMatrix_m(5, 0) << endl;
3100  os << "* R62 = " << correlationMatrix_m(5, 1) << endl;
3101  os << "* R51 = " << correlationMatrix_m(4, 0) << endl;
3102  os << "* R52 = " << correlationMatrix_m(4, 1) << endl;
3103 }
3104 
3105 void Distribution::printDistFlattop(Inform &os) const {
3106 
3107  switch (distrTypeT_m) {
3108 
3109  case DistrTypeT::ASTRAFLATTOPTH:
3110  os << "* Distribution type: ASTRAFLATTOPTH" << endl;
3111  break;
3112 
3113  case DistrTypeT::GUNGAUSSFLATTOPTH:
3114  os << "* Distribution type: GUNGAUSSFLATTOPTH" << endl;
3115  break;
3116 
3117  default:
3118  os << "* Distribution type: FLATTOP" << endl;
3119  break;
3120 
3121  }
3122  os << "* " << endl;
3123 
3124  if (laserProfile_m != NULL) {
3125 
3126  os << "* Transverse profile determined by laser image: " << endl
3127  << endl
3128  << "* Laser profile: " << laserProfileFileName_m << endl
3129  << "* Laser image: " << laserImageName_m << endl
3130  << "* Intensity cut: " << laserIntensityCut_m << endl;
3131 
3132  } else {
3133 
3134  os << "* SIGMAX = " << sigmaR_m[0] << " [m]" << endl;
3135  os << "* SIGMAY = " << sigmaR_m[1] << " [m]" << endl;
3136 
3137  }
3138 
3139  if (emitting_m) {
3140 
3141  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH) {
3142 
3143  os << "* Time Rise = " << tRise_m
3144  << " [sec]" << endl;
3145  os << "* TPULSEFWHM = " << tPulseLengthFWHM_m
3146  << " [sec]" << endl;
3147 
3148  } else {
3149  os << "* Sigma Time Rise = " << sigmaTRise_m
3150  << " [sec]" << endl;
3151  os << "* TPULSEFWHM = " << tPulseLengthFWHM_m
3152  << " [sec]" << endl;
3153  os << "* Sigma Time Fall = " << sigmaTFall_m
3154  << " [sec]" << endl;
3155  os << "* Longitudinal cutoff = " << cutoffR_m[2]
3156  << " [units of Sigma Time]" << endl;
3157  os << "* Flat top modulation amplitude = "
3158  << Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCAMPLITUDE])
3159  << " [Percent of distribution amplitude]" << endl;
3160  os << "* Flat top modulation periods = "
3161  << std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCPERIODS]))
3162  << endl;
3163  }
3164 
3165  } else
3166  os << "* SIGMAZ = " << sigmaR_m[2]
3167  << " [m]" << endl;
3168 }
3169 
3170 void Distribution::printDistFromFile(Inform &os) const {
3171  os << "* Distribution type: FROMFILE" << endl;
3172  os << "* " << endl;
3173  os << "* Input file: "
3174  << Attributes::getString(itsAttr[Attrib::Distribution::FNAME]) << endl;
3175 }
3176 
3177 
3178 void Distribution::printDistMatchedGauss(Inform &os) const {
3179  os << "* Distribution type: MATCHEDGAUSS" << endl;
3180  os << "* SIGMAX = " << sigmaR_m[0] << " [m]" << endl;
3181  os << "* SIGMAY = " << sigmaR_m[1] << " [m]" << endl;
3182  os << "* SIGMAZ = " << sigmaR_m[2] << " [m]" << endl;
3183  os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
3184  os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
3185  os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
3186  os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
3187 
3188  os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
3189  os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
3190  os << "* CORRZ = " << correlationMatrix_m(5, 4) << endl;
3191  os << "* R61 = " << correlationMatrix_m(5, 0) << endl;
3192  os << "* R62 = " << correlationMatrix_m(5, 1) << endl;
3193  os << "* R51 = " << correlationMatrix_m(4, 0) << endl;
3194  os << "* R52 = " << correlationMatrix_m(4, 1) << endl;
3195  os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
3196  os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
3197  os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
3198  os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
3199  os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
3200  os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
3201 }
3202 
3203 void Distribution::printDistGauss(Inform &os) const {
3204  os << "* Distribution type: GAUSS" << endl;
3205  os << "* " << endl;
3206  if (emitting_m) {
3207  os << "* Sigma Time Rise = " << sigmaTRise_m
3208  << " [sec]" << endl;
3209  os << "* TPULSEFWHM = " << tPulseLengthFWHM_m
3210  << " [sec]" << endl;
3211  os << "* Sigma Time Fall = " << sigmaTFall_m
3212  << " [sec]" << endl;
3213  os << "* Longitudinal cutoff = " << cutoffR_m[2]
3214  << " [units of Sigma Time]" << endl;
3215  os << "* Flat top modulation amplitude = "
3216  << Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCAMPLITUDE])
3217  << " [Percent of distribution amplitude]" << endl;
3218  os << "* Flat top modulation periods = "
3219  << std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::FTOSCPERIODS]))
3220  << endl;
3221  os << "* SIGMAX = " << sigmaR_m[0] << " [m]" << endl;
3222  os << "* SIGMAY = " << sigmaR_m[1] << " [m]" << endl;
3223  os << "* SIGMAPX = " << sigmaP_m[0]
3224  << " [Beta Gamma]" << endl;
3225  os << "* SIGMAPY = " << sigmaP_m[1]
3226  << " [Beta Gamma]" << endl;
3227  os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
3228  os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
3229  os << "* CUTOFFX = " << cutoffR_m[0]
3230  << " [units of SIGMAX]" << endl;
3231  os << "* CUTOFFY = " << cutoffR_m[1]
3232  << " [units of SIGMAY]" << endl;
3233  os << "* CUTOFFPX = " << cutoffP_m[0]
3234  << " [units of SIGMAPX]" << endl;
3235  os << "* CUTOFFPY = " << cutoffP_m[1]
3236  << " [units of SIGMAPY]" << endl;
3237  } else {
3238  os << "* SIGMAX = " << sigmaR_m[0] << " [m]" << endl;
3239  os << "* SIGMAY = " << sigmaR_m[1] << " [m]" << endl;
3240  os << "* SIGMAZ = " << sigmaR_m[2] << " [m]" << endl;
3241  os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
3242  os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
3243  os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
3244  os << "* AVRGPZ = " << avrgpz_m << " [Beta Gamma]" << endl;
3245 
3246  os << "* CORRX = " << correlationMatrix_m(1, 0) << endl;
3247  os << "* CORRY = " << correlationMatrix_m(3, 2) << endl;
3248  os << "* CORRZ = " << correlationMatrix_m(5, 4) << endl;
3249  os << "* R61 = " << correlationMatrix_m(5, 0) << endl;
3250  os << "* R62 = " << correlationMatrix_m(5, 1) << endl;
3251  os << "* R51 = " << correlationMatrix_m(4, 0) << endl;
3252  os << "* R52 = " << correlationMatrix_m(4, 1) << endl;
3253  os << "* CUTOFFX = " << cutoffR_m[0] << " [units of SIGMAX]" << endl;
3254  os << "* CUTOFFY = " << cutoffR_m[1] << " [units of SIGMAY]" << endl;
3255  os << "* CUTOFFLONG = " << cutoffR_m[2] << " [units of SIGMAZ]" << endl;
3256  os << "* CUTOFFPX = " << cutoffP_m[0] << " [units of SIGMAPX]" << endl;
3257  os << "* CUTOFFPY = " << cutoffP_m[1] << " [units of SIGMAPY]" << endl;
3258  os << "* CUTOFFPZ = " << cutoffP_m[2] << " [units of SIGMAPY]" << endl;
3259  }
3260 }
3261 
3262 void Distribution::printDistSurfEmission(Inform &os) const {
3263 
3264  os << "* Distribution type: SURFACEEMISION" << endl;
3265  os << "* " << endl;
3266  os << "* * Number of electrons for surface emission "
3267  << Attributes::getReal(itsAttr[Attrib::Distribution::NPDARKCUR]) << endl;
3268  os << "* * Initialized electrons inward margin for surface emission "
3269  << Attributes::getReal(itsAttr[Attrib::Distribution::INWARDMARGIN]) << endl;
3270  os << "* * E field threshold (MV), only in position r with E(r)>EINITHR that "
3271  << "particles will be initialized "
3272  << Attributes::getReal(itsAttr[Attrib::Distribution::EINITHR]) << endl;
3273  os << "* * Field enhancement for surface emission "
3274  << Attributes::getReal(itsAttr[Attrib::Distribution::FNBETA]) << endl;
3275  os << "* * Maximum number of electrons emitted from a single triangle for "
3276  << "Fowler-Nordheim emission "
3277  << Attributes::getReal(itsAttr[Attrib::Distribution::FNMAXEMI]) << endl;
3278  os << "* * Field Threshold for Fowler-Nordheim emission (MV/m) "
3279  << Attributes::getReal(itsAttr[Attrib::Distribution::FNFIELDTHR]) << endl;
3280  os << "* * Empirical constant A for Fowler-Nordheim emission model "
3281  << Attributes::getReal(itsAttr[Attrib::Distribution::FNA]) << endl;
3282  os << "* * Empirical constant B for Fowler-Nordheim emission model "
3283  << Attributes::getReal(itsAttr[Attrib::Distribution::FNB]) << endl;
3284  os << "* * Constant for image charge effect parameter y(E) in Fowler-Nordheim "
3285  << "emission model "
3286  << Attributes::getReal(itsAttr[Attrib::Distribution::FNY]) << endl;
3287  os << "* * Zero order constant for image charge effect parameter v(y) in "
3288  << "Fowler-Nordheim emission model "
3289  << Attributes::getReal(itsAttr[Attrib::Distribution::FNVYZERO]) << endl;
3290  os << "* * Second order constant for image charge effect parameter v(y) in "
3291  << "Fowler-Nordheim emission model "
3292  << Attributes::getReal(itsAttr[Attrib::Distribution::FNVYSECOND]) << endl;
3293  os << "* * Select secondary model type(0:no secondary emission; 1:Furman-Pivi; "
3294  << "2 or larger: Vaughan's model "
3295  << Attributes::getReal(itsAttr[Attrib::Distribution::SECONDARYFLAG]) << endl;
3296  os << "* * Secondary emission mode type(true: emit n true secondaries; false: "
3297  << "emit one particle with n times charge "
3298  << Attributes::getBool(itsAttr[Attrib::Distribution::NEMISSIONMODE]) << endl;
3299  os << "* * Sey_0 in Vaughan's model "
3300  << Attributes::getReal(itsAttr[Attrib::Distribution::VSEYZERO]) << endl;
3301  os << "* * Energy related to sey_0 in Vaughan's model in eV "
3302  << Attributes::getReal(itsAttr[Attrib::Distribution::VEZERO]) << endl;
3303  os << "* * Sey max in Vaughan's model "
3304  << Attributes::getReal(itsAttr[Attrib::Distribution::VSEYMAX]) << endl;
3305  os << "* * Emax in Vaughan's model in eV "
3306  << Attributes::getReal(itsAttr[Attrib::Distribution::VEMAX]) << endl;
3307  os << "* * Fitting parameter denotes the roughness of surface for impact "
3308  << "energy in Vaughan's model "
3309  << Attributes::getReal(itsAttr[Attrib::Distribution::VKENERGY]) << endl;
3310  os << "* * Fitting parameter denotes the roughness of surface for impact angle "
3311  << "in Vaughan's model "
3312  << Attributes::getReal(itsAttr[Attrib::Distribution::VKTHETA]) << endl;
3313  os << "* * Thermal velocity of Maxwellian distribution of secondaries in "
3314  << "Vaughan's model "
3315  << Attributes::getReal(itsAttr[Attrib::Distribution::VVTHERMAL]) << endl;
3316  os << "* * VW denote the velocity scalar for Parallel plate benchmark "
3317  << Attributes::getReal(itsAttr[Attrib::Distribution::VW]) << endl;
3318  os << "* * Material type number of the cavity surface for Furman-Pivi's model, "
3319  << "0 for copper, 1 for stainless steel "
3320  << Attributes::getReal(itsAttr[Attrib::Distribution::SURFMATERIAL]) << endl;
3321 }
3322 
3323 void Distribution::printDistSurfAndCreate(Inform &os) const {
3324 
3325  os << "* Distribution type: SURFACERANDCREATE" << endl;
3326  os << "* " << endl;
3327  os << "* * Number of electrons initialized on the surface as primaries "
3328  << Attributes::getReal(itsAttr[Attrib::Distribution::NPDARKCUR]) << endl;
3329  os << "* * Initialized electrons inward margin for surface emission "
3330  << Attributes::getReal(itsAttr[Attrib::Distribution::INWARDMARGIN]) << endl;
3331  os << "* * E field threshold (MV), only in position r with E(r)>EINITHR that "
3332  << "particles will be initialized "
3333  << Attributes::getReal(itsAttr[Attrib::Distribution::EINITHR]) << endl;
3334 }
3335 
3336 void Distribution::printEmissionModel(Inform &os) const {
3337 
3338  os << "* ------------- THERMAL EMITTANCE MODEL --------------------------------------------" << endl;
3339 
3340  switch (emissionModel_m) {
3341 
3342  case EmissionModelT::NONE:
3343  printEmissionModelNone(os);
3344  break;
3345  case EmissionModelT::ASTRA:
3346  printEmissionModelAstra(os);
3347  break;
3348  case EmissionModelT::NONEQUIL:
3349  printEmissionModelNonEquil(os);
3350  break;
3351  default:
3352  break;
3353  }
3354 
3355  os << "* ----------------------------------------------------------------------------------" << endl;
3356 
3357 }
3358 
3359 void Distribution::printEmissionModelAstra(Inform &os) const {
3360  os << "* THERMAL EMITTANCE in ASTRA MODE" << endl;
3361  os << "* Kinetic energy (thermal emittance) = "
3362  << std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]))
3363  << " [eV] " << endl;
3364 }
3365 
3366 void Distribution::printEmissionModelNone(Inform &os) const {
3367  os << "* THERMAL EMITTANCE in NONE MODE" << endl;
3368  os << "* Kinetic energy added to longitudinal momentum = "
3369  << std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]))
3370  << " [eV] " << endl;
3371 }
3372 
3373 void Distribution::printEmissionModelNonEquil(Inform &os) const {
3374  os << "* THERMAL EMITTANCE in NONEQUIL MODE" << endl;
3375  os << "* Cathode work function = " << cathodeWorkFunc_m << " [eV] " << endl;
3376  os << "* Cathode Fermi energy = " << cathodeFermiEnergy_m << " [eV] " << endl;
3377  os << "* Cathode temperature = " << cathodeTemp_m / Physics::kB << " [K] " << endl;
3378  os << "* Photocathode laser energy = " << laserEnergy_m << " [eV] " << endl;
3379 }
3380 
3381 void Distribution::printEnergyBins(Inform &os) const {
3382 
3383  os << "* " << endl;
3384  for (int binIndex = numberOfEnergyBins_m - 1; binIndex >=0; binIndex--) {
3385  size_t sum = 0;
3386  for (int sampleIndex = 0; sampleIndex < numberOfSampleBins_m; sampleIndex++)
3387  sum += gsl_histogram_get(energyBinHist_m,
3388  binIndex * numberOfSampleBins_m + sampleIndex);
3389 
3390  os << "* Energy Bin # " << numberOfEnergyBins_m - binIndex
3391  << " contains " << sum << " particles" << endl;
3392  }
3393  os << "* " << endl;
3394 
3395 }
3396 
3397 bool Distribution::Rebin() {
3398  /*
3399  * Allow a rebin (numberOfEnergyBins_m = 0) if all particles are emitted or
3400  * injected.
3401  */
3402  if (!emitting_m) {
3403  numberOfEnergyBins_m = 0;
3404  return true;
3405  } else {
3406  return false;
3407  }
3408 }
3409 
3410 void Distribution::reflectDistribution(size_t &numberOfParticles) {
3411 
3412  if (!Options::cZero || (distrTypeT_m == DistrTypeT::FROMFILE))
3413  return;
3414 
3415  size_t currentNumPart = tOrZDist_m.size();
3416  for (size_t partIndex = 0; partIndex < currentNumPart; partIndex++) {
3417  xDist_m.push_back(-xDist_m.at(partIndex));
3418  pxDist_m.push_back(-pxDist_m.at(partIndex));
3419  yDist_m.push_back(-yDist_m.at(partIndex));
3420  pyDist_m.push_back(-pyDist_m.at(partIndex));
3421  tOrZDist_m.push_back(tOrZDist_m.at(partIndex));
3422  pzDist_m.push_back(pzDist_m.at(partIndex));
3423  }
3424  numberOfParticles = tOrZDist_m.size();
3425  reduce(numberOfParticles, numberOfParticles, OpAddAssign());
3426 
3427  // if numberOfParticles is odd then delete last particle
3428  if (Ippl::myNode() == 0 &&
3429  (numberOfParticles + 1) / 2 != numberOfParticles / 2) {
3430  xDist_m.pop_back();
3431  pxDist_m.pop_back();
3432  yDist_m.pop_back();
3433  pyDist_m.pop_back();
3434  tOrZDist_m.pop_back();
3435  pzDist_m.pop_back();
3436  }
3437 }
3438 
3439 void Distribution::scaleDistCoordinates() {
3440  // at this point the distributions of an array of distributions are still separated
3441 
3442  const double xmult = Attributes::getReal(itsAttr[Attrib::Distribution::XMULT]);
3443  const double pxmult = Attributes::getReal(itsAttr[Attrib::Distribution::PXMULT]);
3444  const double ymult = Attributes::getReal(itsAttr[Attrib::Distribution::YMULT]);
3445  const double pymult = Attributes::getReal(itsAttr[Attrib::Distribution::PYMULT]);
3446  const double longmult = (emitting_m ?
3447  Attributes::getReal(itsAttr[Attrib::Distribution::TMULT]) :
3448  Attributes::getReal(itsAttr[Attrib::Distribution::ZMULT]));
3449  const double pzmult = Attributes::getReal(itsAttr[Attrib::Distribution::PZMULT]);
3450 
3451  for (size_t particleIndex = 0; particleIndex < tOrZDist_m.size(); ++ particleIndex) {
3452  xDist_m.at(particleIndex) *= xmult;
3453  pxDist_m.at(particleIndex) *= pxmult;
3454  yDist_m.at(particleIndex) *= ymult;
3455  pyDist_m.at(particleIndex) *= pymult;
3456  tOrZDist_m.at(particleIndex) *= longmult;
3457  pzDist_m.at(particleIndex) *= pzmult;
3458  }
3459 }
3460 
3461 gsl_qrng* Distribution::selectRandomGenerator(std::string,unsigned int dimension)
3462 {
3463  gsl_qrng *quasiRandGen = nullptr;
3464  if (Options::rngtype != std::string("RANDOM")) {
3465  INFOMSG("RNGTYPE= " << Options::rngtype << endl);
3466  if (Options::rngtype == std::string("HALTON")) {
3467  quasiRandGen = gsl_qrng_alloc(gsl_qrng_halton, dimension);
3468  } else if (Options::rngtype == std::string("SOBOL")) {
3469  quasiRandGen = gsl_qrng_alloc(gsl_qrng_sobol, dimension);
3470  } else if (Options::rngtype == std::string("NIEDERREITER")) {
3471  quasiRandGen = gsl_qrng_alloc(gsl_qrng_niederreiter_2, dimension);
3472  } else {
3473  INFOMSG("RNGTYPE= " << Options::rngtype << " not known, using HALTON" << endl);
3474  quasiRandGen = gsl_qrng_alloc(gsl_qrng_halton, dimension);
3475  }
3476  }
3477  return quasiRandGen;
3478 }
3479 
3480 void Distribution::setAttributes() {
3481  itsAttr[Attrib::Distribution::TYPE]
3482  = Attributes::makeString("TYPE","Distribution type: "
3483  "FROMFILE, "
3484  "GAUSS, "
3485  "BINOMIAL, "
3486  "FLATTOP, "
3487  "SURFACEEMISSION, "
3488  "SURFACERANDCREATE, "
3489  "GUNGAUSSFLATTOPTH, "
3490  "ASTRAFLATTOPTH, "
3491  "GAUSSMATCHED");
3492  itsAttr[Attrib::Legacy::Distribution::DISTRIBUTION]
3493  = Attributes::makeString("DISTRIBUTION","This attribute isn't supported any more. Use TYPE instead");
3494  itsAttr[Attrib::Distribution::LINE]
3495  = Attributes::makeString("LINE", "Beamline that contains a cyclotron or ring ", "");
3496  itsAttr[Attrib::Distribution::EX]
3497  = Attributes::makeReal("EX", "Projected normalized emittance EX (m-rad), used to create matched distribution ", 1E-6);
3498  itsAttr[Attrib::Distribution::EY]
3499  = Attributes::makeReal("EY", "Projected normalized emittance EY (m-rad) used to create matched distribution ", 1E-6);
3500  itsAttr[Attrib::Distribution::ET]
3501  = Attributes::makeReal("ET", "Projected normalized emittance ET (m-rad) used to create matched distribution ", 1E-6);
3502  // itsAttr[Attrib::Distribution::E2]
3503  // = Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
3504  itsAttr[Attrib::Distribution::RESIDUUM]
3505  = Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generator ", 1e-8);
3506  itsAttr[Attrib::Distribution::MAXSTEPSCO]
3507  = Attributes::makeReal("MAXSTEPSCO", "Maximum steps used to find closed orbit ", 100);
3508  itsAttr[Attrib::Distribution::MAXSTEPSSI]
3509  = Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",500);
3510  itsAttr[Attrib::Distribution::ORDERMAPS]
3511  = Attributes::makeReal("ORDERMAPS", "Order used in the field expansion ", 7);
3512  itsAttr[Attrib::Distribution::SECTOR]
3513  = Attributes::makeBool("SECTOR", "Match using single sector (true)"
3514  " (false: using all sectors) (default: true)",
3515  true);
3516  itsAttr[Attrib::Distribution::NSTEPS]
3517  = Attributes::makeReal("NSTEPS", "Number of integration steps of closed orbit finder (matched gauss)"
3518  " (default: 720)", 720);
3519 
3520  itsAttr[Attrib::Distribution::RGUESS]
3521  = Attributes::makeReal("RGUESS", "Guess value of radius (m) for closed orbit finder ", -1);
3522 
3523  itsAttr[Attrib::Distribution::DENERGY]
3524  = Attributes::makeReal("DENERGY", "Energy step size for closed orbit finder [GeV]", 0.001);
3525 
3526  itsAttr[Attrib::Distribution::NSECTORS]
3527  = Attributes::makeReal("NSECTORS", "Number of sectors to average field, for closed orbit finder ", 1);
3528 
3529  itsAttr[Attrib::Distribution::FNAME]
3530  = Attributes::makeString("FNAME", "File for reading in 6D particle "
3531  "coordinates.", "");
3532 
3533 
3534  itsAttr[Attrib::Distribution::WRITETOFILE]
3535  = Attributes::makeBool("WRITETOFILE", "Write initial distribution to file.",
3536  false);
3537 
3538  itsAttr[Attrib::Distribution::WEIGHT]
3539  = Attributes::makeReal("WEIGHT", "Weight of distribution when used in a "
3540  "distribution list.", 1.0);
3541 
3542  itsAttr[Attrib::Distribution::INPUTMOUNITS]
3543  = Attributes::makeString("INPUTMOUNITS", "Tell OPAL what input units are for momentum."
3544  " Currently \"NONE\" or \"EV\".", "");
3545 
3546  // Attributes for beam emission.
3547  itsAttr[Attrib::Distribution::EMITTED]
3548  = Attributes::makeBool("EMITTED", "Emitted beam, from cathode, as opposed to "
3549  "an injected beam.", false);
3550  itsAttr[Attrib::Distribution::EMISSIONSTEPS]
3551  = Attributes::makeReal("EMISSIONSTEPS", "Number of time steps to use during emission.",
3552  1);
3553  itsAttr[Attrib::Distribution::EMISSIONMODEL]
3554  = Attributes::makeString("EMISSIONMODEL", "Model used to emit electrons from a "
3555  "photocathode.", "None");
3556  itsAttr[Attrib::Distribution::EKIN]
3557  = Attributes::makeReal("EKIN", "Kinetic energy used in ASTRA thermal emittance "
3558  "model (eV). (Thermal energy added in with random "
3559  "number generator.)", 1.0);
3560  itsAttr[Attrib::Distribution::ELASER]
3561  = Attributes::makeReal("ELASER", "Laser energy (eV) for photocathode "
3562  "emission. (Default 255 nm light.)", 4.86);
3563  itsAttr[Attrib::Distribution::W]
3564  = Attributes::makeReal("W", "Workfunction of photocathode material (eV)."
3565  " (Default atomically clean copper.)", 4.31);
3566  itsAttr[Attrib::Distribution::FE]
3567  = Attributes::makeReal("FE", "Fermi energy of photocathode material (eV)."
3568  " (Default atomically clean copper.)", 7.0);
3569  itsAttr[Attrib::Distribution::CATHTEMP]
3570  = Attributes::makeReal("CATHTEMP", "Temperature of photocathode (K)."
3571  " (Default 300 K.)", 300.0);
3572  itsAttr[Attrib::Distribution::NBIN]
3573  = Attributes::makeReal("NBIN", "Number of energy bins to use when doing a space "
3574  "charge solve.", 0.0);
3575 
3576  // Parameters for shifting or scaling initial distribution.
3577  itsAttr[Attrib::Distribution::XMULT] = Attributes::makeReal("XMULT", "Multiplier for x dimension.", 1.0);
3578  itsAttr[Attrib::Distribution::YMULT] = Attributes::makeReal("YMULT", "Multiplier for y dimension.", 1.0);
3579  itsAttr[Attrib::Distribution::ZMULT] = Attributes::makeReal("TMULT", "Multiplier for z dimension.", 1.0);
3580  itsAttr[Attrib::Distribution::TMULT] = Attributes::makeReal("TMULT", "Multiplier for t dimension.", 1.0);
3581 
3582  itsAttr[Attrib::Distribution::PXMULT] = Attributes::makeReal("PXMULT", "Multiplier for px dimension.", 1.0);
3583  itsAttr[Attrib::Distribution::PYMULT] = Attributes::makeReal("PYMULT", "Multiplier for py dimension.", 1.0);
3584  itsAttr[Attrib::Distribution::PZMULT] = Attributes::makeReal("PZMULT", "Multiplier for pz dimension.", 1.0);
3585 
3586  itsAttr[Attrib::Distribution::OFFSETX]
3587  = Attributes::makeReal("OFFSETX", "Average x offset of distribution.", 0.0);
3588  itsAttr[Attrib::Distribution::OFFSETY]
3589  = Attributes::makeReal("OFFSETY", "Average y offset of distribution.", 0.0);
3590  itsAttr[Attrib::Distribution::OFFSETZ]
3591  = Attributes::makeReal("OFFSETZ", "Average z offset of distribution.", 0.0);
3592  itsAttr[Attrib::Distribution::OFFSETT]
3593  = Attributes::makeReal("OFFSETT", "Average t offset of distribution.", 0.0);
3594 
3595  itsAttr[Attrib::Distribution::OFFSETPX]
3596  = Attributes::makeReal("OFFSETPX", "Average px offset of distribution.", 0.0);
3597  itsAttr[Attrib::Distribution::OFFSETPY]
3598  = Attributes::makeReal("OFFSETPY", "Average py offset of distribution.", 0.0);
3599  itsAttr[Attrib::Distribution::OFFSETPZ]
3600  = Attributes::makeReal("OFFSETPZ", "Average pz offset of distribution.", 0.0);
3601  itsAttr[Attrib::Distribution::OFFSETP]
3602  = Attributes::makeReal("OFFSETP", "Momentum shift relative to referenc particle.", 0.0);
3603 
3604  // Parameters for defining an initial distribution.
3605  itsAttr[Attrib::Distribution::SIGMAX] = Attributes::makeReal("SIGMAX", "SIGMAx (m)", 0.0);
3606  itsAttr[Attrib::Distribution::SIGMAY] = Attributes::makeReal("SIGMAY", "SIGMAy (m)", 0.0);
3607  itsAttr[Attrib::Distribution::SIGMAR] = Attributes::makeReal("SIGMAR", "SIGMAr (m)", 0.0);
3608  itsAttr[Attrib::Distribution::SIGMAZ] = Attributes::makeReal("SIGMAZ", "SIGMAz (m)", 0.0);
3609  itsAttr[Attrib::Distribution::SIGMAT] = Attributes::makeReal("SIGMAT", "SIGMAt (m)", 0.0);
3610  itsAttr[Attrib::Distribution::TPULSEFWHM]
3611  = Attributes::makeReal("TPULSEFWHM", "Pulse FWHM for emitted distribution.", 0.0);
3612  itsAttr[Attrib::Distribution::TRISE]
3613  = Attributes::makeReal("TRISE", "Rise time for emitted distribution.", 0.0);
3614  itsAttr[Attrib::Distribution::TFALL]
3615  = Attributes::makeReal("TFALL", "Fall time for emitted distribution.", 0.0);
3616  itsAttr[Attrib::Distribution::SIGMAPX] = Attributes::makeReal("SIGMAPX", "SIGMApx", 0.0);
3617  itsAttr[Attrib::Distribution::SIGMAPY] = Attributes::makeReal("SIGMAPY", "SIGMApy", 0.0);
3618  itsAttr[Attrib::Distribution::SIGMAPZ] = Attributes::makeReal("SIGMAPZ", "SIGMApz", 0.0);
3619 
3620  itsAttr[Attrib::Distribution::MX]
3621  = Attributes::makeReal("MX", "Defines the binomial distribution in x, "
3622  "0.0 ... infinity.", 10001.0);
3623  itsAttr[Attrib::Distribution::MY]
3624  = Attributes::makeReal("MY", "Defines the binomial distribution in y, "
3625  "0.0 ... infinity.", 10001.0);
3626  itsAttr[Attrib::Distribution::MZ]
3627  = Attributes::makeReal("MZ", "Defines the binomial distribution in z, "
3628  "0.0 ... infinity.", 10001.0);
3629  itsAttr[Attrib::Distribution::MT]
3630  = Attributes::makeReal("MT", "Defines the binomial distribution in t, "
3631  "0.0 ... infinity.", 10001.0);
3632 
3633  itsAttr[Attrib::Distribution::CUTOFFX] = Attributes::makeReal("CUTOFFX", "Distribution cutoff x "
3634  "direction in units of sigma.", 3.0);
3635  itsAttr[Attrib::Distribution::CUTOFFY] = Attributes::makeReal("CUTOFFY", "Distribution cutoff r "
3636  "direction in units of sigma.", 3.0);
3637  itsAttr[Attrib::Distribution::CUTOFFR] = Attributes::makeReal("CUTOFFR", "Distribution cutoff radial "
3638  "direction in units of sigma.", 3.0);
3639  itsAttr[Attrib::Distribution::CUTOFFLONG]
3640  = Attributes::makeReal("CUTOFFLONG", "Distribution cutoff z or t direction in "
3641  "units of sigma.", 3.0);
3642  itsAttr[Attrib::Distribution::CUTOFFPX] = Attributes::makeReal("CUTOFFPX", "Distribution cutoff px "
3643  "dimension in units of sigma.", 3.0);
3644  itsAttr[Attrib::Distribution::CUTOFFPY] = Attributes::makeReal("CUTOFFPY", "Distribution cutoff py "
3645  "dimension in units of sigma.", 3.0);
3646  itsAttr[Attrib::Distribution::CUTOFFPZ] = Attributes::makeReal("CUTOFFPZ", "Distribution cutoff pz "
3647  "dimension in units of sigma.", 3.0);
3648 
3649  itsAttr[Attrib::Distribution::FTOSCAMPLITUDE]
3650  = Attributes::makeReal("FTOSCAMPLITUDE", "Amplitude of oscillations superimposed "
3651  "on flat top portion of emitted GAUSS "
3652  "distribtuion (in percent of flat top "
3653  "amplitude)",0.0);
3654  itsAttr[Attrib::Distribution::FTOSCPERIODS]
3655  = Attributes::makeReal("FTOSCPERIODS", "Number of oscillations superimposed on "
3656  "flat top portion of emitted GAUSS "
3657  "distribution", 0.0);
3658 
3659  /*
3660  * TODO: Find out what these correlations really mean and write
3661  * good descriptions for them.
3662  */
3663  itsAttr[Attrib::Distribution::CORRX]
3664  = Attributes::makeReal("CORRX", "x/px correlation, (R12 in transport "
3665  "notation).", 0.0);
3666  itsAttr[Attrib::Distribution::CORRY]
3667  = Attributes::makeReal("CORRY", "y/py correlation, (R34 in transport "
3668  "notation).", 0.0);
3669  itsAttr[Attrib::Distribution::CORRZ]
3670  = Attributes::makeReal("CORRZ", "z/pz correlation, (R56 in transport "
3671  "notation).", 0.0);
3672  itsAttr[Attrib::Distribution::CORRT]
3673  = Attributes::makeReal("CORRT", "t/pt correlation, (R56 in transport "
3674  "notation).", 0.0);
3675 
3676  itsAttr[Attrib::Distribution::R51]
3677  = Attributes::makeReal("R51", "x/z correlation, (R51 in transport "
3678  "notation).", 0.0);
3679  itsAttr[Attrib::Distribution::R52]
3680  = Attributes::makeReal("R52", "xp/z correlation, (R52 in transport "
3681  "notation).", 0.0);
3682 
3683  itsAttr[Attrib::Distribution::R61]
3684  = Attributes::makeReal("R61", "x/pz correlation, (R61 in transport "
3685  "notation).", 0.0);
3686  itsAttr[Attrib::Distribution::R62]
3687  = Attributes::makeReal("R62", "xp/pz correlation, (R62 in transport "
3688  "notation).", 0.0);
3689 
3690  itsAttr[Attrib::Distribution::R]
3691  = Attributes::makeRealArray("R", "r correlation");
3692 
3693  // Parameters for using laser profile to generate a distribution.
3694  itsAttr[Attrib::Distribution::LASERPROFFN]
3695  = Attributes::makeString("LASERPROFFN", "File containing a measured laser image "
3696  "profile (x,y).", "");
3697  itsAttr[Attrib::Distribution::IMAGENAME]
3698  = Attributes::makeString("IMAGENAME", "Name of the laser image.", "");
3699  itsAttr[Attrib::Distribution::INTENSITYCUT]
3700  = Attributes::makeReal("INTENSITYCUT", "For background subtraction of laser "
3701  "image profile, in % of max intensity.",
3702  0.0);
3703  itsAttr[Attrib::Distribution::FLIPX]
3704  = Attributes::makeBool("FLIPX", "Flip laser profile horizontally", false);
3705  itsAttr[Attrib::Distribution::FLIPY]
3706  = Attributes::makeBool("FLIPY", "Flip laser profile vertically", false);
3707  itsAttr[Attrib::Distribution::ROTATE90]
3708  = Attributes::makeBool("ROTATE90", "Rotate laser profile 90 degrees counter clockwise", false);
3709  itsAttr[Attrib::Distribution::ROTATE180]
3710  = Attributes::makeBool("ROTATE180", "Rotate laser profile 180 degrees counter clockwise", false);
3711  itsAttr[Attrib::Distribution::ROTATE270]
3712  = Attributes::makeBool("ROTATE270", "Rotate laser profile 270 degrees counter clockwise", false);
3713 
3714  // Dark current and field emission model parameters.
3715  itsAttr[Attrib::Distribution::NPDARKCUR]
3716  = Attributes::makeReal("NPDARKCUR", "Number of dark current particles.", 1000.0);
3717  itsAttr[Attrib::Distribution::INWARDMARGIN]
3718  = Attributes::makeReal("INWARDMARGIN", "Inward margin of initialized dark "
3719  "current particle positions.", 0.001);
3720  itsAttr[Attrib::Distribution::EINITHR]
3721  = Attributes::makeReal("EINITHR", "E field threshold (MV/m), only in position r "
3722  "with E(r)>EINITHR that particles will be "
3723  "initialized.", 0.0);
3724  itsAttr[Attrib::Distribution::FNA]
3725  = Attributes::makeReal("FNA", "Empirical constant A for Fowler-Nordheim "
3726  "emission model.", 1.54e-6);
3727  itsAttr[Attrib::Distribution::FNB]
3728  = Attributes::makeReal("FNB", "Empirical constant B for Fowler-Nordheim "
3729  "emission model.", 6.83e9);
3730  itsAttr[Attrib::Distribution::FNY]
3731  = Attributes::makeReal("FNY", "Constant for image charge effect parameter y(E) "
3732  "in Fowler-Nordheim emission model.", 3.795e-5);
3733  itsAttr[Attrib::Distribution::FNVYZERO]
3734  = Attributes::makeReal("FNVYZERO", "Zero order constant for v(y) function in "
3735  "Fowler-Nordheim emission model.", 0.9632);
3736  itsAttr[Attrib::Distribution::FNVYSECOND]
3737  = Attributes::makeReal("FNVYSECOND", "Second order constant for v(y) function "
3738  "in Fowler-Nordheim emission model.", 1.065);
3739  itsAttr[Attrib::Distribution::FNPHIW]
3740  = Attributes::makeReal("FNPHIW", "Work function of gun surface material (eV).",
3741  4.65);
3742  itsAttr[Attrib::Distribution::FNBETA]
3743  = Attributes::makeReal("FNBETA", "Field enhancement factor for Fowler-Nordheim "
3744  "emission.", 50.0);
3745  itsAttr[Attrib::Distribution::FNFIELDTHR]
3746  = Attributes::makeReal("FNFIELDTHR", "Field threshold for Fowler-Nordheim "
3747  "emission (MV/m).", 30.0);
3748  itsAttr[Attrib::Distribution::FNMAXEMI]
3749  = Attributes::makeReal("FNMAXEMI", "Maximum number of electrons emitted from a "
3750  "single triangle for Fowler-Nordheim "
3751  "emission.", 20.0);
3752  itsAttr[Attrib::Distribution::SECONDARYFLAG]
3753  = Attributes::makeReal("SECONDARYFLAG", "Select the secondary model type 0:no "
3754  "secondary emission; 1:Furman-Pivi; "
3755  "2 or larger: Vaughan's model.", 0.0);
3756  itsAttr[Attrib::Distribution::NEMISSIONMODE]
3757  = Attributes::makeBool("NEMISSIONMODE", "Secondary emission mode type true: "
3758  "emit n true secondaries; false: emit "
3759  "one particle with n times charge.", true);
3760  itsAttr[Attrib::Distribution::VSEYZERO]
3761  = Attributes::makeReal("VSEYZERO", "Sey_0 in Vaughan's model.", 0.5);
3762  itsAttr[Attrib::Distribution::VEZERO]
3763  = Attributes::makeReal("VEZERO", "Energy related to sey_0 in Vaughan's model "
3764  "in eV.", 12.5);
3765  itsAttr[Attrib::Distribution::VSEYMAX]
3766  = Attributes::makeReal("VSEYMAX", "Sey max in Vaughan's model.", 2.22);
3767  itsAttr[Attrib::Distribution::VEMAX]
3768  = Attributes::makeReal("VEMAX", "Emax in Vaughan's model in eV.", 165.0);
3769  itsAttr[Attrib::Distribution::VKENERGY]
3770  = Attributes::makeReal("VKENERGY", "Fitting parameter denotes the roughness of "
3771  "surface for impact energy in Vaughan's "
3772  "model.", 1.0);
3773  itsAttr[Attrib::Distribution::VKTHETA]
3774  = Attributes::makeReal("VKTHETA", "Fitting parameter denotes the roughness of "
3775  "surface for impact angle in Vaughan's "
3776  "model.", 1.0);
3777  itsAttr[Attrib::Distribution::VVTHERMAL]
3778  = Attributes::makeReal("VVTHERMAL", "Thermal velocity of Maxwellian distribution "
3779  "of secondaries in Vaughan's model.", 7.268929821 * 1e5);
3780  itsAttr[Attrib::Distribution::VW]
3781  = Attributes::makeReal("VW", "VW denote the velocity scalar for parallel plate "
3782  "benchmark.", 1.0);
3783  itsAttr[Attrib::Distribution::SURFMATERIAL]
3784  = Attributes::makeReal("SURFMATERIAL", "Material type number of the cavity "
3785  "surface for Furman-Pivi's model, 0 "
3786  "for copper, 1 for stainless steel.", 0.0);
3787 
3788  /*
3789  * Feature request Issue #14
3790  */
3791 
3792  itsAttr[Attrib::Distribution::ID1]
3793  = Attributes::makeRealArray("ID1", "User defined particle with ID=1");
3794  itsAttr[Attrib::Distribution::ID2]
3795  = Attributes::makeRealArray("ID2", "User defined particle with ID=2");
3796 
3797 
3798  itsAttr[Attrib::Distribution::SCALABLE]
3799  = Attributes::makeBool("SCALABLE", "If true then distribution is scalable with "
3800  "respect of number of particles and number of cores", false);
3801 
3802  /*
3803  * Legacy attributes (or ones that need to be implemented.)
3804  */
3805 
3806  // Parameters for emitting a distribution.
3807  // itsAttr[Attrib::Legacy::Distribution::DEBIN]
3808  // = Attributes::makeReal("DEBIN", "Energy band for a bin in keV that defines "
3809  // "when to combine bins. That is, when the energy "
3810  // "spread of all bins is below this number "
3811  // "combine bins into a single bin.", 1000000.0);
3812  itsAttr[Attrib::Legacy::Distribution::SBIN]
3813  = Attributes::makeReal("SBIN", "Number of sample bins to use per energy bin "
3814  "when emitting a distribution.", 100.0);
3815  /*
3816  * Specific to type GAUSS and GUNGAUSSFLATTOPTH and ASTRAFLATTOPTH. These
3817  * last two distribution will eventually just become a subset of GAUSS.
3818  */
3819  itsAttr[Attrib::Legacy::Distribution::SIGMAPT] = Attributes::makeReal("SIGMAPT", "SIGMApt", 0.0);
3820 
3821  itsAttr[Attrib::Legacy::Distribution::CUTOFF]
3822  = Attributes::makeReal("CUTOFF", "Longitudinal cutoff for Gaussian in units "
3823  "of sigma.", 3.0);
3824 
3825 
3826  // Mixed use attributes (used by more than one distribution type).
3827  itsAttr[Attrib::Legacy::Distribution::T]
3828  = Attributes::makeReal("T", "Not supported anymore");
3829 
3830  itsAttr[Attrib::Legacy::Distribution::PT]
3831  = Attributes::makeReal("PT", "Not supported anymore.");
3832 
3833 
3834  // Attributes that are not yet implemented.
3835  // itsAttr[Attrib::Legacy::Distribution::ALPHAX]
3836  // = Attributes::makeReal("ALPHAX", "Courant Snyder parameter.", 0.0);
3837  // itsAttr[Attrib::Legacy::Distribution::ALPHAY]
3838  // = Attributes::makeReal("ALPHAY", "Courant Snyder parameter.", 0.0);
3839  // itsAttr[Attrib::Legacy::Distribution::BETAX]
3840  // = Attributes::makeReal("BETAX", "Courant Snyder parameter.", 1.0);
3841  // itsAttr[Attrib::Legacy::Distribution::BETAY]
3842  // = Attributes::makeReal("BETAY", "Courant Snyder parameter.", 1.0);
3843 
3844  // itsAttr[Attrib::Legacy::Distribution::DX]
3845  // = Attributes::makeReal("DX", "Dispersion in x (R16 in Transport notation).", 0.0);
3846  // itsAttr[Attrib::Legacy::Distribution::DDX]
3847  // = Attributes::makeReal("DDX", "First derivative of DX.", 0.0);
3848 
3849  // itsAttr[Attrib::Legacy::Distribution::DY]
3850  // = Attributes::makeReal("DY", "Dispersion in y (R36 in Transport notation).", 0.0);
3851  // itsAttr[Attrib::Legacy::Distribution::DDY]
3852  // = Attributes::makeReal("DDY", "First derivative of DY.", 0.0);
3853 
3854  registerOwnership(AttributeHandler::STATEMENT);
3855 }
3856 
3857 void Distribution::setFieldEmissionParameters() {
3858 
3859  darkCurrentParts_m = static_cast<size_t> (Attributes::getReal(itsAttr[Attrib::Distribution::NPDARKCUR]));
3860  darkInwardMargin_m = Attributes::getReal(itsAttr[Attrib::Distribution::INWARDMARGIN]);
3861  eInitThreshold_m = Attributes::getReal(itsAttr[Attrib::Distribution::EINITHR]);
3862  workFunction_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNPHIW]);
3863  fieldEnhancement_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNBETA]);
3864  maxFN_m = static_cast<size_t> (Attributes::getReal(itsAttr[Attrib::Distribution::FNMAXEMI]));
3865  fieldThrFN_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNFIELDTHR]);
3866  paraFNA_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNA]);
3867  paraFNB_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNB]);
3868  paraFNY_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNY]);
3869  paraFNVYZe_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNVYZERO]);
3870  paraFNVYSe_m = Attributes::getReal(itsAttr[Attrib::Distribution::FNVYSECOND]);
3871  secondaryFlag_m = Attributes::getReal(itsAttr[Attrib::Distribution::SECONDARYFLAG]);
3872  ppVw_m = Attributes::getReal(itsAttr[Attrib::Distribution::VW]);
3873  vVThermal_m = Attributes::getReal(itsAttr[Attrib::Distribution::VVTHERMAL]);
3874 
3875 }
3876 
3877 void Distribution::setDistToEmitted(bool emitted) {
3878  emitting_m = emitted;
3879 }
3880 
3881 void Distribution::setDistType() {
3882  if (itsAttr[Attrib::Legacy::Distribution::DISTRIBUTION]) {
3883  throw OpalException("Distribution::setDistType()",
3884  "The attribute DISTRIBUTION isn't supported any more, use TYPE instead");
3885  }
3886 
3887  distT_m = Util::toUpper(Attributes::getString(itsAttr[Attrib::Distribution::TYPE]));
3888  if (distT_m == "FROMFILE")
3889  distrTypeT_m = DistrTypeT::FROMFILE;
3890  else if (distT_m == "GAUSS")
3891  distrTypeT_m = DistrTypeT::GAUSS;
3892  else if (distT_m == "BINOMIAL")
3893  distrTypeT_m = DistrTypeT::BINOMIAL;
3894  else if (distT_m == "FLATTOP")
3895  distrTypeT_m = DistrTypeT::FLATTOP;
3896  else if (distT_m == "GUNGAUSSFLATTOPTH")
3897  distrTypeT_m = DistrTypeT::GUNGAUSSFLATTOPTH;
3898  else if (distT_m == "ASTRAFLATTOPTH")
3899  distrTypeT_m = DistrTypeT::ASTRAFLATTOPTH;
3900  else if (distT_m == "SURFACEEMISSION")
3901  distrTypeT_m = DistrTypeT::SURFACEEMISSION;
3902  else if (distT_m == "SURFACERANDCREATE")
3903  distrTypeT_m = DistrTypeT::SURFACERANDCREATE;
3904  else if (distT_m == "GAUSSMATCHED")
3905  distrTypeT_m = DistrTypeT::MATCHEDGAUSS;
3906  else {
3907  throw OpalException("Distribution::setDistType()",
3908  "The distribution \"" + distT_m + "\" isn't known.\n" +
3909  "Known distributions are:\n"
3910  "FROMFILE\n"
3911  "GAUSS\n"
3912  "BINOMIAL\n"
3913  "FLATTOP\n"
3914  "GUNGAUSSFLATTOPTH\n"
3915  "ASTRAFLATTTOPTH\n"
3916  "SURFACEEMISSION\n"
3917  "SURFACERANDCREATE\n"
3918  "GAUSSMATCHED");
3919  }
3920 }
3921 
3922 void Distribution::setEmissionTime(double &maxT, double &minT) {
3923 
3924  if (addedDistributions_m.size() == 0) {
3925 
3926  switch (distrTypeT_m) {
3927 
3928  case DistrTypeT::FLATTOP:
3929  case DistrTypeT::GAUSS:
3930  case DistrTypeT::GUNGAUSSFLATTOPTH:
3931  tEmission_m = tPulseLengthFWHM_m + (cutoffR_m[2] - sqrt(2.0 * log(2.0)))
3932  * (sigmaTRise_m + sigmaTFall_m);
3933  break;
3934 
3935  case DistrTypeT::ASTRAFLATTOPTH:
3936  /*
3937  * Don't do anything. Emission time is set during the distribution
3938  * creation. Only this distribution type does it this way. This is
3939  * a legacy feature.
3940  */
3941  break;
3942 
3943  default:
3944  /*
3945  * Increase maxT and decrease minT by 0.05% of total time
3946  * to ensure that no particles fall on the leading edge of
3947  * the first emission time step or the trailing edge of the last
3948  * emission time step.
3949  */
3950  double deltaT = maxT - minT;
3951  maxT += deltaT * 0.0005;
3952  minT -= deltaT * 0.0005;
3953  tEmission_m = (maxT - minT);
3954  break;
3955  }
3956 
3957  } else {
3958  /*
3959  * Increase maxT and decrease minT by 0.05% of total time
3960  * to ensure that no particles fall on the leading edge of
3961  * the first emission time step or the trailing edge of the last
3962  * emission time step.
3963  */
3964  double deltaT = maxT - minT;
3965  maxT += deltaT * 0.0005;
3966  minT -= deltaT * 0.0005;
3967  tEmission_m = (maxT - minT);
3968  }
3969  tBin_m = tEmission_m / numberOfEnergyBins_m;
3970 }
3971 
3972 void Distribution::setDistParametersBinomial(double massIneV) {
3973 
3974  /*
3975  * Set Distribution parameters. Do all the necessary checks depending
3976  * on the input attributes.
3977  */
3978  std::vector<double> cr = Attributes::getRealArray(itsAttr[Attrib::Distribution::R]);
3979 
3980  if (cr.size()>0) {
3981  throw OpalException("Distribution::setDistParametersBinomial",
3982  "Attribute R is not supported for binomial distribution\n"
3983  "use CORR[X|Y|Z] and R51, R52, R61, R62 instead");
3984  }
3985 
3986  correlationMatrix_m(1, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRX]);
3987  correlationMatrix_m(3, 2) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRY]);
3988  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRT]);
3989  correlationMatrix_m(4, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::R51]);
3990  correlationMatrix_m(4, 1) = Attributes::getReal(itsAttr[Attrib::Distribution::R52]);
3991  correlationMatrix_m(5, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::R61]);
3992  correlationMatrix_m(5, 1) = Attributes::getReal(itsAttr[Attrib::Distribution::R62]);
3993 
3994  //CORRZ overrides CORRT. We initially use CORRT for legacy compatibility.
3995  if (!itsAttr[Attrib::Distribution::CORRZ].defaultUsed())
3996  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]);
3997 
3998  sigmaR_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAX])),
3999  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAY])),
4000  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT])));
4001 
4002  // SIGMAZ overrides SIGMAT. We initially use SIGMAT for legacy compatibility.
4003  if (Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAZ]) != 0.0)
4004  sigmaR_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAZ]));
4005 
4006  sigmaP_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPX])),
4007  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPY])),
4008  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPZ])));
4009 
4010  // SIGMAPZ overrides SIGMAPT.
4011  if (!itsAttr[Attrib::Legacy::Distribution::SIGMAPT].defaultUsed()) {
4012  WARNMSG("The attribute SIGMAPT may be removed in a future version\n"
4013  << "use SIGMAPZ instead" << endl;)
4014  sigmaP_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Legacy::Distribution::SIGMAPT]));
4015  }
4016  // Check what input units we are using for momentum.
4017  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
4018  sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
4019  sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
4020  sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
4021  }
4022 
4023  mBinomial_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::MX])),
4024  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::MY])),
4025  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::MT])));
4026 
4027  cutoffR_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFX]),
4028  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFY]),
4029  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFLONG]));
4030 
4031  cutoffP_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPX]),
4032  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPY]),
4033  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPZ]));
4034 
4035  if (mBinomial_m[2] == 0.0
4036  || Attributes::getReal(itsAttr[Attrib::Distribution::MZ]) != 0.0)
4037  mBinomial_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::MZ]));
4038 
4039  if (emitting_m) {
4040  sigmaR_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
4041  mBinomial_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::MT]));
4042  correlationMatrix_m(5, 4) = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::CORRT]));
4043  }
4044 }
4045 
4046 void Distribution::setDistParametersFlattop(double massIneV) {
4047 
4048  /*
4049  * Set distribution parameters. Do all the necessary checks depending
4050  * on the input attributes.
4051  */
4052  sigmaP_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPX])),
4053  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPY])),
4054  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPZ])));
4055 
4056  // Check what input units we are using for momentum.
4057  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
4058  sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
4059  sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
4060  sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
4061  }
4062 
4063  cutoffR_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFX]),
4064  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFY]),
4065  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFLONG]));
4066 
4067  correlationMatrix_m(1, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRX]);
4068  correlationMatrix_m(3, 2) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRY]);
4069  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRT]);
4070 
4071  // CORRZ overrides CORRT.
4072  if (Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]) != 0.0)
4073  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]);
4074 
4075 
4076  if (emitting_m) {
4077  sigmaR_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAX])),
4078  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAY])),
4079  0.0);
4080 
4081  sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
4082  sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
4083 
4084  tPulseLengthFWHM_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TPULSEFWHM]));
4085 
4086  /*
4087  * If TRISE and TFALL are defined > 0.0 then these attributes
4088  * override SIGMAT.
4089  */
4090  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0
4091  || std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0) {
4092 
4093  double timeRatio = sqrt(2.0 * log(10.0)) - sqrt(2.0 * log(10.0 / 9.0));
4094 
4095  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0)
4096  sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]))
4097  / timeRatio;
4098 
4099  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0)
4100  sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL]))
4101  / timeRatio;
4102 
4103  }
4104 
4105  // For an emitted beam, the longitudinal cutoff >= 0.
4106  cutoffR_m[2] = std::abs(cutoffR_m[2]);
4107 
4108  } else {
4109 
4110  sigmaR_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAX])),
4111  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAY])),
4112  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAZ])));
4113 
4114  }
4115 
4116  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR])) > 0.0) {
4117  sigmaR_m[0] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR]));
4118  sigmaR_m[1] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR]));
4119  }
4120 
4121  // Set laser profile/
4122  laserProfileFileName_m = Attributes::getString(itsAttr[Attrib::Distribution::LASERPROFFN]);
4123  if (!(laserProfileFileName_m == std::string(""))) {
4124  laserImageName_m = Attributes::getString(itsAttr[Attrib::Distribution::IMAGENAME]);
4125  laserIntensityCut_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::INTENSITYCUT]));
4126  short flags = 0;
4127  if (Attributes::getBool(itsAttr[Attrib::Distribution::FLIPX])) flags |= LaserProfile::FLIPX;
4128  if (Attributes::getBool(itsAttr[Attrib::Distribution::FLIPY])) flags |= LaserProfile::FLIPY;
4129  if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE90])) flags |= LaserProfile::ROTATE90;
4130  if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE180])) flags |= LaserProfile::ROTATE180;
4131  if (Attributes::getBool(itsAttr[Attrib::Distribution::ROTATE270])) flags |= LaserProfile::ROTATE270;
4132 
4133  laserProfile_m = new LaserProfile(laserProfileFileName_m,
4134  laserImageName_m,
4135  laserIntensityCut_m,
4136  flags);
4137  }
4138 
4139  // Legacy for ASTRAFLATTOPTH.
4140  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH)
4141  tRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]));
4142 
4143 }
4144 
4145 void Distribution::setDistParametersGauss(double massIneV) {
4146 
4147  /*
4148  * Set distribution parameters. Do all the necessary checks depending
4149  * on the input attributes.
4150  * In case of DistrTypeT::MATCHEDGAUSS we only need to set the cutoff parameters
4151  */
4152 
4153 
4154  cutoffP_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPX]),
4155  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPY]),
4156  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFPZ]));
4157 
4158 
4159  cutoffR_m = Vector_t(Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFX]),
4160  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFY]),
4161  Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFLONG]));
4162 
4163  if (distrTypeT_m != DistrTypeT::MATCHEDGAUSS) {
4164  sigmaP_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPX])),
4165  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPY])),
4166  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPZ])));
4167 
4168  // SIGMAPZ overrides SIGMAPT. We initially use SIGMAPT for legacy compatibility.
4169  if (Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPZ]) != 0.0)
4170  sigmaP_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAPZ]));
4171 
4172  // Check what input units we are using for momentum.
4173  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
4174  sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
4175  sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
4176  sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
4177  }
4178 
4179  std::vector<double> cr = Attributes::getRealArray(itsAttr[Attrib::Distribution::R]);
4180 
4181  if (cr.size()>0) {
4182  if (cr.size() == 15) {
4183  *gmsg << "* Use r to specify correlations" << endl;
4184  unsigned int k = 0;
4185  for (unsigned int i = 0; i < 5; ++ i) {
4186  for (unsigned int j = i + 1; j < 6; ++ j, ++ k) {
4187  correlationMatrix_m(j, i) = cr.at(k);
4188  }
4189  }
4190 
4191  }
4192  else {
4193  throw OpalException("Distribution::SetDistParametersGauss",
4194  "Inconsistent set of correlations specified, check manual");
4195  }
4196  }
4197  else {
4198  correlationMatrix_m(1, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRX]);
4199  correlationMatrix_m(3, 2) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRY]);
4200  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRT]);
4201  correlationMatrix_m(4, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::R51]);
4202  correlationMatrix_m(4, 1) = Attributes::getReal(itsAttr[Attrib::Distribution::R52]);
4203  correlationMatrix_m(5, 0) = Attributes::getReal(itsAttr[Attrib::Distribution::R61]);
4204  correlationMatrix_m(5, 1) = Attributes::getReal(itsAttr[Attrib::Distribution::R62]);
4205 
4206  // CORRZ overrides CORRT.
4207  if (Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]) != 0.0)
4208  correlationMatrix_m(5, 4) = Attributes::getReal(itsAttr[Attrib::Distribution::CORRZ]);
4209  }
4210  }
4211 
4212  if (emitting_m) {
4213 
4214  sigmaR_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAX])),
4215  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAY])),
4216  0.0);
4217 
4218  sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
4219  sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT]));
4220 
4221  tPulseLengthFWHM_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TPULSEFWHM]));
4222 
4223  /*
4224  * If TRISE and TFALL are defined then these attributes
4225  * override SIGMAT.
4226  */
4227  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0
4228  || std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0) {
4229 
4230  double timeRatio = sqrt(2.0 * log(10.0)) - sqrt(2.0 * log(10.0 / 9.0));
4231 
4232  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE])) > 0.0)
4233  sigmaTRise_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]))
4234  / timeRatio;
4235 
4236  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL])) > 0.0)
4237  sigmaTFall_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::TFALL]))
4238  / timeRatio;
4239 
4240  }
4241 
4242  // For and emitted beam, the longitudinal cutoff >= 0.
4243  cutoffR_m[2] = std::abs(cutoffR_m[2]);
4244 
4245  } else {
4246  if (distrTypeT_m != DistrTypeT::MATCHEDGAUSS) {
4247  sigmaR_m = Vector_t(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAX])),
4248  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAY])),
4249  std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAT])));
4250 
4251  // SIGMAZ overrides SIGMAT. We initially use SIGMAT for legacy compatibility.
4252  if (Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAZ]) != 0.0)
4253  sigmaR_m[2] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAZ]));
4254 
4255  }
4256  }
4257 
4258  if (std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR])) > 0.0) {
4259  sigmaR_m[0] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR]));
4260  sigmaR_m[1] = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::SIGMAR]));
4261  cutoffR_m[0] = Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFR]);
4262  cutoffR_m[1] = Attributes::getReal(itsAttr[Attrib::Distribution::CUTOFFR]);
4263  }
4264 }
4265 
4266 void Distribution::setupEmissionModel(PartBunchBase<double, 3> *beam) {
4267 
4268  std::string model = Util::toUpper(Attributes::getString(itsAttr[Attrib::Distribution::EMISSIONMODEL]));
4269  if (model == "ASTRA")
4270  emissionModel_m = EmissionModelT::ASTRA;
4271  else if (model == "NONEQUIL")
4272  emissionModel_m = EmissionModelT::NONEQUIL;
4273  else
4274  emissionModel_m = EmissionModelT::NONE;
4275 
4276  /*
4277  * The ASTRAFLATTOPTH of GUNGAUSSFLATTOPTH distributions always uses the
4278  * ASTRA emission model.
4279  */
4280  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH
4281  || distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH)
4282  emissionModel_m = EmissionModelT::ASTRA;
4283 
4284  switch (emissionModel_m) {
4285 
4286  case EmissionModelT::ASTRA:
4287  setupEmissionModelAstra(beam);
4288  break;
4289  case EmissionModelT::NONEQUIL:
4290  setupEmissionModelNonEquil();
4291  break;
4292  default:
4293  break;
4294  }
4295 
4296 }
4297 
4298 void Distribution::setupEmissionModelAstra(PartBunchBase<double, 3> *beam) {
4299 
4300  double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
4301  pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
4302  pmean_m = Vector_t(0.0, 0.0, 0.5 * pTotThermal_m);
4303 }
4304 
4305 void Distribution::setupEmissionModelNone(PartBunchBase<double, 3> *beam) {
4306 
4307  double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
4308  pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
4309  double avgPz = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
4310  size_t numParticles = pzDist_m.size();
4311  reduce(avgPz, avgPz, OpAddAssign());
4312  reduce(numParticles, numParticles, OpAddAssign());
4313  avgPz /= numParticles;
4314 
4315  pmean_m = Vector_t(0.0, 0.0, pTotThermal_m + avgPz);
4316 }
4317 
4318 void Distribution::setupEmissionModelNonEquil() {
4319 
4320  cathodeWorkFunc_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::W]));
4321  laserEnergy_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::ELASER]));
4322  cathodeFermiEnergy_m = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::FE]));
4323  cathodeTemp_m = Physics::kB * std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::CATHTEMP]));
4324 
4325  /*
4326  * Upper limit on energy distribution theoretically goes to infinity.
4327  * Practically we limit to a probability of 1 part in a billion.
4328  */
4329  emitEnergyUpperLimit_m = cathodeFermiEnergy_m
4330  + cathodeTemp_m * log(1.0e9 - 1.0);
4331 
4332  // TODO: get better estimate of pmean
4333  pmean_m = Vector_t(0, 0, sqrt(pow(0.5 * emitEnergyUpperLimit_m / (Physics::m_e * 1e9) + 1.0, 2) - 1.0));
4334 }
4335 
4336 void Distribution::setupEnergyBins(double maxTOrZ, double minTOrZ) {
4337 
4338  energyBinHist_m = gsl_histogram_alloc(numberOfSampleBins_m * numberOfEnergyBins_m);
4339 
4340  if (emitting_m)
4341  gsl_histogram_set_ranges_uniform(energyBinHist_m, -tEmission_m, 0.0);
4342  else
4343  gsl_histogram_set_ranges_uniform(energyBinHist_m, minTOrZ, maxTOrZ);
4344 
4345 }
4346 
4347 void Distribution::setupParticleBins(double massIneV, PartBunchBase<double, 3> *beam) {
4348 
4349  numberOfEnergyBins_m
4350  = static_cast<int>(std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::NBIN])));
4351 
4352  if (numberOfEnergyBins_m > 0) {
4353  delete energyBins_m;
4354 
4355  int sampleBins = static_cast<int>(std::abs(Attributes::getReal(itsAttr[Attrib::Legacy::Distribution::SBIN])));
4356  energyBins_m = new PartBins(numberOfEnergyBins_m, sampleBins);
4357 
4358  if (!itsAttr[Attrib::Legacy::Distribution::PT].defaultUsed())
4359  throw OpalException("Distribution::setupParticleBins",
4360  "PT is obsolete. The moments of the beam is defined with OFFSETPZ");
4361 
4362  // we get gamma from PC of the beam
4363  const double pz = beam->getP()/beam->getM();
4364  double gamma = sqrt(pow(pz, 2.0) + 1.0);
4365  energyBins_m->setGamma(gamma);
4366 
4367  } else {
4368  energyBins_m = NULL;
4369  }
4370 }
4371 
4372 void Distribution::shiftBeam(double &maxTOrZ, double &minTOrZ) {
4373 
4374  if (emitting_m) {
4375 
4376  if (addedDistributions_m.size() == 0) {
4377 
4378  if (distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH) {
4379  for (double& tOrZ : tOrZDist_m)
4380  tOrZ -= tEmission_m / 2.0;
4381 
4382  minTOrZ -= tEmission_m / 2.0;
4383  maxTOrZ -= tEmission_m / 2.0;
4384  } else if (distrTypeT_m == DistrTypeT::GAUSS
4385  || distrTypeT_m == DistrTypeT::FLATTOP
4386  || distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
4387  for (double& tOrZ : tOrZDist_m)
4388  tOrZ -= tEmission_m;
4389 
4390  minTOrZ -= tEmission_m;
4391  maxTOrZ -= tEmission_m;
4392  } else {
4393  for (double& tOrZ : tOrZDist_m)
4394  tOrZ -= maxTOrZ;
4395 
4396  minTOrZ -= maxTOrZ;
4397  maxTOrZ -= maxTOrZ;
4398  }
4399 
4400  } else {
4401  for (double& tOrZ : tOrZDist_m)
4402  tOrZ -= maxTOrZ;
4403 
4404  minTOrZ -= maxTOrZ;
4405  maxTOrZ -= maxTOrZ;
4406  }
4407 
4408  } else if (distrTypeT_m != DistrTypeT::FROMFILE) {
4409  double avgZ[] = {0.0, 1.0 * tOrZDist_m.size()};
4410  for (double tOrZ : tOrZDist_m)
4411  avgZ[0] += tOrZ;
4412 
4413  reduce(avgZ, avgZ + 2, avgZ, OpAddAssign());
4414  avgZ[0] /= avgZ[1];
4415 
4416  for (double& tOrZ : tOrZDist_m)
4417  tOrZ -= avgZ[0];
4418  }
4419 }
4420 
4421 double Distribution::getEmissionTimeShift() const {
4422  if (emitting_m)
4423  return Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETT]);
4424 
4425  return 0.0;
4426 }
4427 
4428 void Distribution::shiftDistCoordinates(double massIneV) {
4429 
4430  size_t startIdx = 0;
4431  for (unsigned int i = 0; i <= addedDistributions_m.size(); ++ i) {
4432  Distribution *currDist = this;
4433  if (i > 0)
4434  currDist = addedDistributions_m[i - 1];
4435  double deltaX = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETX]);
4436  double deltaY = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETY]);
4437 
4438  /*
4439  * OFFSETZ overrides T if it is nonzero. We initially use T
4440  * for legacy compatiblity. OFFSETT always overrides T, even
4441  * when zero, for an emitted beam.
4442  */
4443  if (currDist->itsAttr[Attrib::Legacy::Distribution::T]) {
4444  throw OpalException("Distribution::shiftDistCoordinates",
4445  "Attribute T isn't supported anymore; use OFFSETZ instead");
4446  }
4447 
4448  double deltaTOrZ = 0.0;
4449  if (!emitting_m)
4450  if (Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETZ]) != 0.0)
4451  deltaTOrZ = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETZ]);
4452 
4453  double deltaPx = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETPX]);
4454  double deltaPy = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETPY]);
4455  double deltaPz = Attributes::getReal(currDist->itsAttr[Attrib::Distribution::OFFSETPZ]);
4456 
4457  if (Attributes::getReal(currDist->itsAttr[Attrib::Legacy::Distribution::PT])!=0.0)
4458  WARNMSG("PT & PZ are obsolete and will be ignored. The moments of the beam is defined with PC" << endl);
4459 
4460  // Check input momentum units.
4461  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
4462  deltaPx = converteVToBetaGamma(deltaPx, massIneV);
4463  deltaPy = converteVToBetaGamma(deltaPy, massIneV);
4464  deltaPz = converteVToBetaGamma(deltaPz, massIneV);
4465  }
4466 
4467  size_t endIdx = startIdx + particlesPerDist_m[i];
4468  for (size_t particleIndex = startIdx; particleIndex < endIdx; ++ particleIndex) {
4469  xDist_m.at(particleIndex) += deltaX;
4470  pxDist_m.at(particleIndex) += deltaPx;
4471  yDist_m.at(particleIndex) += deltaY;
4472  pyDist_m.at(particleIndex) += deltaPy;
4473  tOrZDist_m.at(particleIndex) += deltaTOrZ;
4474  pzDist_m.at(particleIndex) += deltaPz;
4475  }
4476 
4477  startIdx = endIdx;
4478  }
4479 }
4480 
4481 void Distribution::writeOutFileHeader() {
4482 
4483  if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE]) == false)
4484  return;
4485 
4486  unsigned int totalNum = tOrZDist_m.size();
4487  reduce(totalNum, totalNum, OpAddAssign());
4488  if (Ippl::myNode() != 0)
4489  return;
4490 
4491  std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
4492  *gmsg << "\n"
4493  << std::left << std::setw(84) << std::setfill('*') << "* " << "\n"
4494  << "* Write initial distribution to file \"" << fname << "\"\n"
4495  << std::left << std::setw(84) << std::setfill('*') << "* "
4496  << std::setfill(' ') << endl;
4497 
4498  std::ofstream outputFile(fname);
4499  if (outputFile.bad()) {
4500  *gmsg << "Unable to open output file \"" << fname << "\"" << endl;
4501  } else {
4502  outputFile.setf(std::ios::left);
4503  outputFile << "# ";
4504  outputFile.width(17);
4505  outputFile << "x [m]";
4506  outputFile.width(17);
4507  outputFile << "px [betax gamma]";
4508  outputFile.width(17);
4509  outputFile << "y [m]";
4510  outputFile.width(17);
4511  outputFile << "py [betay gamma]";
4512  if (emitting_m) {
4513  outputFile.width(17);
4514  outputFile << "t [s]";
4515  } else {
4516  outputFile.width(17);
4517  outputFile << "z [m]";
4518  }
4519  outputFile.width(17);
4520  outputFile << "pz [betaz gamma]" ;
4521  if (emitting_m) {
4522  outputFile.width(17);
4523  outputFile << "Bin Number" << std::endl;
4524  } else {
4525  if (numberOfEnergyBins_m > 0) {
4526  outputFile.width(17);
4527  outputFile << "Bin Number";
4528  }
4529  outputFile << std::endl;
4530 
4531  outputFile << "# " << totalNum << std::endl;
4532  }
4533  }
4534  outputFile.close();
4535 }
4536 
4537 void Distribution::writeOutFileEmission() {
4538 
4539  if (!Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE])) {
4540  xWrite_m.clear();
4541  pxWrite_m.clear();
4542  yWrite_m.clear();
4543  pyWrite_m.clear();
4544  tOrZWrite_m.clear();
4545  pzWrite_m.clear();
4546  binWrite_m.clear();
4547 
4548  return;
4549  }
4550 
4551  // Gather particles to be written onto node 0.
4552  std::vector<char> msgbuf;
4553  const size_t bitsPerParticle = (6 * sizeof(double) + sizeof(size_t));
4554  size_t totalSendBits = xWrite_m.size() * bitsPerParticle;
4555 
4556  std::vector<long> numberOfBits(Ippl::getNodes(), 0);
4557  numberOfBits[Ippl::myNode()] = totalSendBits;
4558 
4559  if (Ippl::myNode() == 0) {
4560  MPI_Reduce(MPI_IN_PLACE, &(numberOfBits[0]), Ippl::getNodes(), MPI_UNSIGNED_LONG, MPI_SUM, 0, Ippl::getComm());
4561  } else {
4562  MPI_Reduce(&(numberOfBits[0]), NULL, Ippl::getNodes(), MPI_UNSIGNED_LONG, MPI_SUM, 0, Ippl::getComm());
4563  }
4564 
4565  Ippl::Comm->barrier();
4566  int tag = Ippl::Comm->next_tag(IPPL_APP_TAG2, IPPL_APP_CYCLE);
4567  if (Ippl::myNode() > 0) {
4568  if (totalSendBits > 0) {
4569  msgbuf.reserve(totalSendBits);
4570  const char *buffer;
4571  for (size_t idx = 0; idx < xWrite_m.size(); ++ idx) {
4572  buffer = reinterpret_cast<const char*>(&(xWrite_m[idx]));
4573  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4574  buffer = reinterpret_cast<const char*>(&(pxWrite_m[idx]));
4575  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4576  buffer = reinterpret_cast<const char*>(&(yWrite_m[idx]));
4577  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4578  buffer = reinterpret_cast<const char*>(&(pyWrite_m[idx]));
4579  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4580  buffer = reinterpret_cast<const char*>(&(tOrZWrite_m[idx]));
4581  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4582  buffer = reinterpret_cast<const char*>(&(pzWrite_m[idx]));
4583  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(double));
4584  buffer = reinterpret_cast<const char*>(&(binWrite_m[idx]));
4585  msgbuf.insert(msgbuf.end(), buffer, buffer + sizeof(size_t));
4586  }
4587 
4588  Ippl::Comm->raw_send(&(msgbuf[0]), totalSendBits, 0, tag);
4589  }
4590  } else {
4591  std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
4592  std::ofstream outputFile(fname, std::fstream::app);
4593  if (outputFile.bad()) {
4594  ERRORMSG(level1 << "Unable to write to file \"" << fname << "\"" << endl);
4595  for (int node = 1; node < Ippl::getNodes(); ++ node) {
4596  if (numberOfBits[node] == 0) continue;
4597  char *recvbuf = new char[numberOfBits[node]];
4598  Ippl::Comm->raw_receive(recvbuf, numberOfBits[node], node, tag);
4599  delete[] recvbuf;
4600  }
4601  } else {
4602 
4603  outputFile.precision(9);
4604  outputFile.setf(std::ios::scientific);
4605  outputFile.setf(std::ios::right);
4606 
4607  for (size_t partIndex = 0; partIndex < xWrite_m.size(); partIndex++) {
4608 
4609  outputFile << std::setw(17) << xWrite_m.at(partIndex)
4610  << std::setw(17) << pxWrite_m.at(partIndex)
4611  << std::setw(17) << yWrite_m.at(partIndex)
4612  << std::setw(17) << pyWrite_m.at(partIndex)
4613  << std::setw(17) << tOrZWrite_m.at(partIndex)
4614  << std::setw(17) << pzWrite_m.at(partIndex)
4615  << std::setw(17) << binWrite_m.at(partIndex) << std::endl;
4616  }
4617 
4618  int numSenders = 0;
4619  for (int i = 1; i < Ippl::getNodes(); ++ i) {
4620  if (numberOfBits[i] > 0) numSenders ++;
4621  }
4622 
4623  for (int i = 0; i < numSenders; ++ i) {
4624  int node = Communicate::COMM_ANY_NODE;
4625  char *recvbuf;
4626  const int bufsize = Ippl::Comm->raw_probe_receive(recvbuf, node, tag);
4627 
4628  int j = 0;
4629  while (j < bufsize) {
4630  const double *dbuffer = reinterpret_cast<const double*>(recvbuf + j);
4631  const size_t *sbuffer = reinterpret_cast<const size_t*>(recvbuf + j + 6 * sizeof(double));
4632  outputFile << std::setw(17) << dbuffer[0]
4633  << std::setw(17) << dbuffer[1]
4634  << std::setw(17) << dbuffer[2]
4635  << std::setw(17) << dbuffer[3]
4636  << std::setw(17) << dbuffer[4]
4637  << std::setw(17) << dbuffer[5]
4638  << std::setw(17) << sbuffer[0]
4639  << std::endl;
4640  j += bitsPerParticle;
4641  }
4642 
4643  delete[] recvbuf;
4644 
4645  }
4646  }
4647  outputFile.close();
4648 
4649  }
4650 
4651  // Clear write vectors.
4652  xWrite_m.clear();
4653  pxWrite_m.clear();
4654  yWrite_m.clear();
4655  pyWrite_m.clear();
4656  tOrZWrite_m.clear();
4657  pzWrite_m.clear();
4658  binWrite_m.clear();
4659 
4660 }
4661 
4662 void Distribution::writeOutFileInjection() {
4663 
4664  if (Attributes::getBool(itsAttr[Attrib::Distribution::WRITETOFILE]) == false)
4665  return;
4666 
4667  std::string fname = "data/" + OpalData::getInstance()->getInputBasename() + "_" + getOpalName() + ".dat";
4668  // Nodes take turn writing particles to file.
4669  for (int nodeIndex = 0; nodeIndex < Ippl::getNodes(); nodeIndex++) {
4670 
4671  // Write to file if its our turn.
4672  size_t numberOfParticles = 0;
4673  if (Ippl::myNode() == nodeIndex) {
4674  std::ofstream outputFile(fname, std::fstream::app);
4675  if (outputFile.bad()) {
4676  *gmsg << "Node " << Ippl::myNode() << " unable to write"
4677  << "to file \"" << fname << "\"" << endl;
4678  } else {
4679 
4680  outputFile.precision(9);
4681  outputFile.setf(std::ios::scientific);
4682  outputFile.setf(std::ios::right);
4683 
4684  numberOfParticles = tOrZDist_m.size();
4685  for (size_t partIndex = 0; partIndex < numberOfParticles; partIndex++) {
4686 
4687  outputFile.width(17);
4688  outputFile << xDist_m.at(partIndex);
4689  outputFile.width(17);
4690  outputFile << pxDist_m.at(partIndex);
4691  outputFile.width(17);
4692  outputFile << yDist_m.at(partIndex);
4693  outputFile.width(17);
4694  outputFile << pyDist_m.at(partIndex);
4695  outputFile.width(17);
4696  outputFile << tOrZDist_m.at(partIndex);
4697  outputFile.width(17);
4698  outputFile << pzDist_m.at(partIndex);
4699  if (numberOfEnergyBins_m > 0) {
4700  size_t binNumber = findEBin(tOrZDist_m.at(partIndex));
4701  outputFile.width(17);
4702  outputFile << binNumber;
4703  }
4704  outputFile << std::endl;
4705 
4706  }
4707  }
4708  outputFile.close();
4709  }
4710 
4711  // Wait for writing node before moving on.
4712  reduce(numberOfParticles, numberOfParticles, OpAddAssign());
4713  }
4714 }
4715 
4716 double Distribution::MDependentBehavior::get(double rand) {
4717  return 2.0 * sqrt(1.0 - pow(rand, ami_m));
4718 }
4719 
4720 double Distribution::GaussianLikeBehavior::get(double rand) {
4721  return sqrt(-2.0 * log(rand));
4722 }
4723 
4724 void Distribution::adjustPhaseSpace(double massIneV) {
4725  if (emitting_m || distrTypeT_m == DistrTypeT::FROMFILE || OpalData::getInstance()->isInOPALCyclMode())
4726  return;
4727 
4728  double deltaPx = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPX]);
4729  double deltaPy = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPY]);
4730  // Check input momentum units.
4731  if (inputMoUnits_m == InputMomentumUnitsT::EV) {
4732  deltaPx = converteVToBetaGamma(deltaPx, massIneV);
4733  deltaPy = converteVToBetaGamma(deltaPy, massIneV);
4734  }
4735 
4736  double avrg[6];
4737  avrg[0] = std::accumulate( xDist_m.begin(), xDist_m.end(), 0.0) / totalNumberParticles_m;
4738  avrg[1] = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0) / totalNumberParticles_m;
4739  avrg[2] = std::accumulate( yDist_m.begin(), yDist_m.end(), 0.0) / totalNumberParticles_m;
4740  avrg[3] = std::accumulate(pyDist_m.begin(), pyDist_m.end(), 0.0) / totalNumberParticles_m;
4741  avrg[4] = std::accumulate(tOrZDist_m.begin(), tOrZDist_m.end(), 0.0) / totalNumberParticles_m;
4742  avrg[5] = 0.0;
4743  for (unsigned int i = 0; i < pzDist_m.size(); ++ i) {
4744  // taylor series of sqrt(px*px + py*py + pz*pz) = pz * sqrt(1 + eps*eps) where eps << 1
4745  avrg[5] += (pzDist_m[i] +
4746  (std::pow(pxDist_m[i] + deltaPx, 2) +
4747  std::pow(pyDist_m[i] + deltaPy, 2)) / (2 * pzDist_m[i]));
4748  }
4749  allreduce(&(avrg[0]), 6, std::plus<double>());
4750  avrg[5] /= totalNumberParticles_m;
4751 
4752  // solve
4753  // \sum_{i = 0}^{N-1} \sqrt{(pz_i + \eps)^2 + px_i^2 + py_i^2} = N p
4754  double eps = avrgpz_m - avrg[5];
4755  for (unsigned int i = 0; i < pzDist_m.size(); ++ i) {
4756  xDist_m[i] -= avrg[0];
4757  pxDist_m[i] -= avrg[1];
4758  yDist_m[i] -= avrg[2];
4759  pyDist_m[i] -= avrg[3];
4760  tOrZDist_m[i] -= avrg[4];
4761  pzDist_m[i] += eps;
4762  }
4763 }
4764 
4765 // vi: set et ts=4 sw=4 sts=4:
4766 // Local Variables:
4767 // mode:c++
4768 // c-basic-offset: 4
4769 // indent-tabs-mode:nil
4770 // End:
Options::seed
int seed
The current random seed.
Definition: Options.cpp:41
Distribution::pmean_m
Vector_t pmean_m
Total thermal momentum.
Definition: Distribution.h:490
Distribution::getEmissionDeltaT
double getEmissionDeltaT()
Definition: Distribution.cpp:2998
Distribution::getBGxDist
std::vector< double > & getBGxDist()
Definition: Distribution.cpp:3014
Distribution::doRestartOpalCycl
void doRestartOpalCycl(PartBunchBase< double, 3 > *p, size_t Np, int restartStep, const int specifiedNumBunch, H5PartWrapper *h5wrapper)
Definition: Distribution.cpp:551
Distribution::MDependentBehavior::get
virtual double get(double rand)
Definition: Distribution.cpp:4716
IpplInfo::getNodes
static int getNodes()
Definition: IpplInfo.cpp:773
PartBunchBase::P
ParticleAttrib< Vector_t > P
Definition: PartBunchBase.h:484
Distribution::printDist
void printDist(Inform &os, size_t numberOfParticles) const
Definition: Distribution.cpp:3034
Distribution::emissionModel_m
EmissionModelT::EmissionModelT emissionModel_m
Emission Model.
Definition: Distribution.h:470
Distribution::getWeight
double getWeight()
Definition: Distribution.cpp:3006
Distribution::~Distribution
virtual ~Distribution()
Definition: Distribution.cpp:232
Distribution::applyEmissionModel
void applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector< double > &additionalRNs)
Definition: Distribution.cpp:805
Distribution::chooseInputMomentumUnits
void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits)
Definition: Distribution.cpp:1027
Distribution::getReference
const PartData & getReference() const
Return the embedded CLASSIC PartData.
Definition: Distribution.cpp:755
PartBunchBase::getT
double getT() const
Definition: PartBunchBase.hpp:1094
Distribution::update
virtual void update()
Update this object.
Definition: Distribution.cpp:284
Distribution::yWrite_m
std::vector< double > yWrite_m
Definition: Distribution.h:515
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
LaserProfile::getXY
void getXY(double &x, double &y)
Definition: LaserProfile.cpp:368
Distribution::doRestartOpalE
void doRestartOpalE(EnvelopeBunch *p, size_t Np, int restartStep, H5PartWrapper *h5wrapper)
Definition: Distribution.cpp:632
Distribution::particlesPerDist_m
std::vector< size_t > particlesPerDist_m
Definition: Distribution.h:467
Physics::kB
constexpr double kB
Boltzman&#39;s constant in eV/K.
Definition: Physics.h:67
PartBunchBase::distrCreate_m
IpplTimings::TimerRef distrCreate_m
Definition: PartBunchBase.h:528
acos
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
Definition: PETE.h:808
Physics::e
constexpr double e
The value of .
Definition: Physics.h:40
PartBunchBase::Ef
ParticleAttrib< Vector_t > Ef
Definition: PartBunchBase.h:488
Distribution::darkCurrentParts_m
size_t darkCurrentParts_m
Definition: Distribution.h:547
Cyclotron::setRinit
void setRinit(double rinit)
Definition: Cyclotron.cpp:108
Distribution::tOrZWrite_m
std::vector< double > tOrZWrite_m
Definition: Distribution.h:517
Distribution::adjustPhaseSpace
void adjustPhaseSpace(double massIneV)
Definition: Distribution.cpp:4724
Distribution::initializeBeam
void initializeBeam(PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:2877
OpalData::define
void define(Object *newObject)
Define a new object.
Definition: OpalData.cpp:538
Distribution::cathodeWorkFunc_m
double cathodeWorkFunc_m
Definition: Distribution.h:493
Distribution::xDist_m
std::vector< double > xDist_m
Definition: Distribution.h:505
Definition
The base class for all OPAL definitions.
Definition: Definition.h:30
ClosedOrbitFinder
Finds a closed orbit of a cyclotron for a given energy.
Definition: ClosedOrbitFinder.h:42
Distribution::numberOfSampleBins_m
int numberOfSampleBins_m
Definition: Distribution.h:480
Distribution::pTotThermal_m
double pTotThermal_m
Random number generator.
Definition: Distribution.h:489
PartBunchBase::create
void create(size_t M)
Definition: PartBunchBase.hpp:2319
PartBunchBase::getLocalTrackStep
long long getLocalTrackStep() const
Definition: PartBunchBase.hpp:1540
Distribution::tOrZDist_m
std::vector< double > tOrZDist_m
Definition: Distribution.h:509
EnvelopeBunch
core of the envelope tracker based on Rene Bakkers BET implementation
Definition: EnvelopeBunch.h:60
Distribution::pxWrite_m
std::vector< double > pxWrite_m
Definition: Distribution.h:514
Beam::getMass
double getMass() const
Return Particle&#39;s rest mass in GeV.
Definition: Beam.cpp:197
Distribution::getNumberOfEnergyBins
int getNumberOfEnergyBins()
Definition: Distribution.cpp:2994
Distribution::createOpalCycl
void createOpalCycl(PartBunchBase< double, 3 > *beam, size_t numberOfParticles, double current, const Beamline &bl)
Definition: Distribution.cpp:1517
Cyclotron
Interface for a Cyclotron.
Definition: Cyclotron.h:91
Distribution::canReplaceBy
virtual bool canReplaceBy(Object *object)
Distribution can only be replaced by another distribution.
Definition: Distribution.cpp:273
Distribution::createDistributionFromFile
void createDistributionFromFile(size_t numberOfParticles, double massIneV)
Definition: Distribution.cpp:1099
BoundaryGeometry::getCooridinate
Vector_t getCooridinate(size_t i)
Definition: BoundaryGeometry.h:137
Beam::getNumberOfSlices
size_t getNumberOfSlices()
Return the number of slices.
Definition: Beam.cpp:164
fabs
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
Definition: PETE.h:815
EnvelopeBunch::setCharge
void setCharge(double _Q)
set the charge of the bunch
Definition: EnvelopeBunch.h:203
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:123
Cyclotron::getPHIinit
virtual double getPHIinit() const
Definition: Cyclotron.cpp:128
Distribution::tEmission_m
double tEmission_m
Emission parameters.
Definition: Distribution.h:473
Distribution::darkInwardMargin_m
double darkInwardMargin_m
Number of dark current particles.
Definition: Distribution.h:548
OpalException
The base class for all OPAL exceptions.
Definition: OpalException.h:28
X
Definition: rbendmap.h:8
sin
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
Distribution::GaussianLikeBehavior::get
virtual double get(double rand)
Definition: Distribution.cpp:4720
Distribution::laserProfileFileName_m
std::string laserProfileFileName_m
Definition: Distribution.h:539
Distribution::additionalRNs_m
std::vector< std::vector< double > > additionalRNs_m
Upper limit on emission energy distribution (eV).
Definition: Distribution.h:499
copysign
PETE_TBTree< FnCopysign, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > copysign(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
Definition: IpplExpressions.h:244
Distribution::sigmaFall_m
double sigmaFall_m
Definition: Distribution.h:587
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:399
PartBunchBase::Q
ParticleAttrib< double > Q
Definition: PartBunchBase.h:485
Physics::two_pi
constexpr double two_pi
The value of .
Definition: Physics.h:34
Distribution::shiftBeam
void shiftBeam(double &maxTOrZ, double &minTOrZ)
Definition: Distribution.cpp:4372
Beam::getCurrent
double getCurrent() const
Return the beam current in A.
Definition: Beam.cpp:189
Distribution::calcPartPerDist
void calcPartPerDist(size_t numberOfParticles)
Definition: Distribution.cpp:895
ceil
PETE_TUTree< FnCeil, typename T::PETE_Expr_t > ceil(const PETE_Expr< T > &l)
Definition: PETE.h:811
Distribution::getBGzDist
std::vector< double > & getBGzDist()
Definition: Distribution.cpp:3030
exp
Tps< T > exp(const Tps< T > &x)
Exponential.
Definition: TpsMath.h:165
EnvelopeBunch::initialize
void initialize(int sli, double charge, double energy, double width, double te, double frac, double current, double center, double bX, double bY, double mX, double mY, double Bz, int nbin)
Definition: EnvelopeBunch.cpp:1374
Distribution::generateBinomial
void generateBinomial(size_t numberOfParticles)
Definition: Distribution.cpp:2172
Communicate::raw_probe_receive
virtual int raw_probe_receive(char *&, int &, int &)
Definition: Communicate.h:208
Distribution::inputMoUnits_m
InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits_m
Definition: Distribution.h:527
Distribution::sigmaP_m
Vector_t sigmaP_m
Definition: Distribution.h:532
Distribution::paraFNVYZe_m
double paraFNVYZe_m
Definition: Distribution.h:568
PartData
Particle reference data.
Definition: PartData.h:38
Distribution::currentSampleBin_m
int currentSampleBin_m
Definition: Distribution.h:477
PartBunchBase::getdT
double getdT() const
Definition: PartBunchBase.hpp:1076
gmsg
Inform * gmsg
Definition: Main.cpp:21
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:130
PartBunchBase::getM
double getM() const
Definition: PartBunchBase.hpp:1771
Distribution::emitEnergyUpperLimit_m
double emitEnergyUpperLimit_m
Cathode temperature (K).
Definition: Distribution.h:497
PartBins::sortArray
void sortArray()
Definition: PartBins.cpp:147
Distribution::createDistributionGauss
void createDistributionGauss(size_t numberOfParticles, double massIneV)
Definition: Distribution.cpp:1465
EnvelopeBunch::createBunch
void createBunch()
create and initialize local num slices
Definition: EnvelopeBunch.cpp:491
PartBins::setGamma
void setGamma(double gamma)
Definition: PartBins.h:132
Distribution::getEmissionTimeShift
double getEmissionTimeShift() const
Definition: Distribution.cpp:4421
Distribution::paraFNVYSe_m
double paraFNVYSe_m
Definition: Distribution.h:566
Distribution::find
static Distribution * find(const std::string &name)
Definition: Distribution.cpp:650
Communicate::barrier
void barrier(void)
Definition: Communicate.cpp:440
Distribution::getEnergyBinDeltaT
double getEnergyBinDeltaT()
Definition: Distribution.cpp:3002
Distribution::writeOutFileEmission
void writeOutFileEmission()
Definition: Distribution.cpp:4537
Distribution::setFieldEmissionParameters
void setFieldEmissionParameters()
Definition: Distribution.cpp:3857
IpplInfo::myNode
static int myNode()
Definition: IpplInfo.cpp:794
H5PartWrapper::predecessorIsSameFlavour
virtual bool predecessorIsSameFlavour() const =0
PX
Definition: rbendmap.h:8
Distribution::xWrite_m
std::vector< double > xWrite_m
Definition: Distribution.h:513
IPPL_APP_CYCLE
#define IPPL_APP_CYCLE
Definition: Tags.h:113
Object::itsAttr
std::vector< Attribute > itsAttr
The object attributes (see Attribute.hh).
Definition: Object.h:214
Distribution::createBoundaryGeometry
void createBoundaryGeometry(PartBunchBase< double, 3 > *p, BoundaryGeometry &bg)
Definition: Distribution.cpp:1477
Inform::setf
FmtFlags_t setf(FmtFlags_t setbits, FmtFlags_t field)
Definition: Inform.h:104
Distribution::getYDist
std::vector< double > & getYDist()
Definition: Distribution.cpp:3018
Communicate::raw_send
virtual bool raw_send(void *, int, int, int)
Definition: Communicate.h:192
Distribution::scaleDistCoordinates
void scaleDistCoordinates()
Definition: Distribution.cpp:3439
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1213
Distribution::totalNumberParticles_m
size_t totalNumberParticles_m
Definition: Distribution.h:501
PartBunchBase::getTotalNum
size_t getTotalNum() const
Definition: PartBunchBase.hpp:2225
log
Tps< T > log(const Tps< T > &x)
Natural logarithm.
Definition: TpsMath.h:182
Distribution::addedDistributions_m
std::vector< Distribution * > addedDistributions_m
Vector of distributions to be added to this distribution.
Definition: Distribution.h:466
Distribution::generateFlattopZ
void generateFlattopZ(size_t numberOfParticles)
Definition: Distribution.cpp:2400
Distribution::sigmaTFall_m
double sigmaTFall_m
Definition: Distribution.h:529
TagMaker::next_tag
int next_tag(int t, int s=1000)
Definition: TagMaker.h:43
Distribution::emitParticles
size_t emitParticles(PartBunchBase< double, 3 > *beam, double eZ)
Definition: Distribution.cpp:1823
Options::cZero
bool cZero
if true create symmetric distribution
Definition: Options.cpp:16
Distribution::printInfo
Inform & printInfo(Inform &os) const
Definition: Distribution.cpp:704
Attributes::getBool
bool getBool(const Attribute &attr)
Return logical value.
Definition: Attributes.cpp:66
PartBunchBase::getP
double getP() const
Definition: PartBunchBase.hpp:1777
Physics::m_e
constexpr double m_e
The electron rest mass in GeV.
Definition: Physics.h:85
level2
Inform & level2(Inform &inf)
Definition: Inform.cpp:46
Physics::alpha
constexpr double alpha
The fine structure constant, no dimension.
Definition: Physics.h:79
Distribution::maxFN_m
size_t maxFN_m
Definition: Distribution.h:558
Distribution::fillParticleBins
void fillParticleBins()
Definition: Distribution.cpp:2027
Distribution::setupEmissionModelNonEquil
void setupEmissionModelNonEquil()
Definition: Distribution.cpp:4318
Distribution::setupEnergyBins
void setupEnergyBins(double maxTOrZ, double minTOrZ)
Definition: Distribution.cpp:4336
Cyclotron::getFMHighE
virtual double getFMHighE() const
Definition: Cyclotron.cpp:339
PartBunchBase::distrReload_m
IpplTimings::TimerRef distrReload_m
timer for IC, can not be in Distribution.h
Definition: PartBunchBase.h:527
reduce
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
Definition: GlobalComm.hpp:55
Distribution::setDistParametersBinomial
void setDistParametersBinomial(double massIneV)
Definition: Distribution.cpp:3972
Cyclotron::getCyclHarm
virtual double getCyclHarm() const
Definition: Cyclotron.cpp:281
OpalData::getInstance
static OpalData * getInstance()
Definition: OpalData.cpp:209
Distribution::generateGaussZ
void generateGaussZ(size_t numberOfParticles)
Definition: Distribution.cpp:2464
Distribution::printDistSurfEmission
void printDistSurfEmission(Inform &os) const
Definition: Distribution.cpp:3262
Distribution::numberOfDistributions_m
unsigned int numberOfDistributions_m
and list type for switch statements.
Definition: Distribution.h:457
Distribution::printEmissionModelNonEquil
void printEmissionModelNonEquil(Inform &os) const
Definition: Distribution.cpp:3373
Object::getOpalName
const std::string & getOpalName() const
Return object name.
Definition: Object.cpp:284
Distribution::getNumberOfEmissionSteps
size_t getNumberOfEmissionSteps()
Definition: Distribution.cpp:2990
Beam::getCharge
double getCharge() const
Return the charge number in elementary charge.
Definition: Beam.cpp:193
BeamSequence::find
static BeamSequence * find(const std::string &name)
Find a BeamSequence by name.
Definition: BeamSequence.cpp:36
Distribution::createDistributionFlattop
void createDistributionFlattop(size_t numberOfParticles, double massIneV)
Definition: Distribution.cpp:1056
Beam::getGamma
double getGamma() const
Definition: Beam.cpp:298
PAssert_GE
#define PAssert_GE(a, b)
Definition: PAssert.h:124
Distribution::pzWrite_m
std::vector< double > pzWrite_m
Definition: Distribution.h:518
Distribution::injectBeam
void injectBeam(PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:2889
Distribution::create
void create(size_t &numberOfParticles, double massIneV)
Definition: Distribution.cpp:287
Distribution::createPriPart
void createPriPart(PartBunchBase< double, 3 > *beam, BoundaryGeometry &bg)
Definition: Distribution.cpp:407
PartBunchBase::PType
ParticleAttrib< short > PType
Definition: PartBunchBase.h:495
PartBunchBase::ID
ParticleIndex_t & ID
Definition: PartBunchBase.h:480
Distribution::getVw
double getVw()
Definition: Distribution.h:748
OpalData::getGlobalPhaseShift
double getGlobalPhaseShift()
units: (sec)
Definition: OpalData.cpp:497
Distribution::setupEmissionModelAstra
void setupEmissionModelAstra(PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:4298
Distribution::cutoffR_m
Vector_t cutoffR_m
Definition: Distribution.h:533
Physics::pi
constexpr double pi
The value of .
Definition: Physics.h:31
Distribution::printDistGauss
void printDistGauss(Inform &os) const
Definition: Distribution.cpp:3203
Distribution::numberOfEnergyBins_m
int numberOfEnergyBins_m
Definition: Distribution.h:478
Distribution::printDistFromFile
void printDistFromFile(Inform &os) const
Definition: Distribution.cpp:3170
Distribution::laserProfile_m
LaserProfile * laserProfile_m
Definition: Distribution.h:542
Inform::width
int width() const
Definition: Inform.h:111
Distribution::setAttributes
void setAttributes()
Definition: Distribution.cpp:3480
Distribution::printDistMatchedGauss
void printDistMatchedGauss(Inform &os) const
Definition: Distribution.cpp:3178
Distribution::sigmaTRise_m
double sigmaTRise_m
Definition: Distribution.h:528
PartBins::fill
void fill(std::vector< double > &p)
Add a particle to the temporary container.
Definition: PartBins.h:60
Distribution::execute
virtual void execute()
Execute the command.
Definition: Distribution.cpp:281
Cyclotron::setPRinit
void setPRinit(double prinit)
Definition: Cyclotron.cpp:116
Distribution::printEnergyBins
void printEnergyBins(Inform &os) const
Definition: Distribution.cpp:3381
PartBunchBase::TriID
ParticleAttrib< int > TriID
Definition: PartBunchBase.h:496
PartBins.h
Defines a structure to hold energy bins and their associated data.
Distribution::generateFlattopLaserProfile
void generateFlattopLaserProfile(size_t numberOfParticles)
Definition: Distribution.cpp:2307
allreduce
void allreduce(const T *input, T *output, int count, Op op)
Definition: GlobalComm.hpp:510
Physics::c
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:52
OpalData::addProblemCharacteristicValue
void addProblemCharacteristicValue(const std::string &name, unsigned int value)
Definition: OpalData.cpp:818
Options::rngtype
std::string rngtype
random number generator
Definition: Options.cpp:85
PartBunchBase::setEnergyBins
void setEnergyBins(int numberOfEnergyBins)
Definition: PartBunchBase.hpp:173
Distribution::checkEmissionParameters
void checkEmissionParameters()
Definition: Distribution.cpp:964
Object::registerOwnership
void registerOwnership(const AttributeHandler::OwnerType &itsClass) const
Definition: Object.cpp:194
BeamSequence::fetchLine
virtual Beamline * fetchLine() const =0
Return the embedded CLASSIC beam line.
IpplTimings::startTimer
static void startTimer(TimerRef t)
Definition: IpplTimings.h:187
PartBunchBase::getChargePerParticle
double getChargePerParticle() const
get the macro particle charge
Definition: PartBunchBase.hpp:1462
PartBunchBase::setTEmission
void setTEmission(double t)
Definition: PartBunchBase.hpp:276
BoundaryGeometry::getMomenta
Vector_t getMomenta(size_t i)
Definition: BoundaryGeometry.h:143
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:397
Distribution::findEBin
size_t findEBin(double tOrZ)
Definition: Distribution.cpp:2048
SpecificElementVisitor::execute
virtual void execute()
Execute the algorithm on the attached beam line.
Definition: SpecificElementVisitor.h:272
Distribution::randGen_m
gsl_rng * randGen_m
Definition: Distribution.h:486
Attributes::getRealArray
std::vector< double > getRealArray(const Attribute &attr)
Get array value.
Definition: Attributes.cpp:258
Distribution::generateAstraFlattopT
void generateAstraFlattopT(size_t numberOfParticles)
Definition: Distribution.cpp:2061
SMALLESTCUTOFF
constexpr double SMALLESTCUTOFF
Definition: Distribution.cpp:58
Distribution::particleRefData_m
PartData particleRefData_m
Definition: Distribution.h:462
Distribution::cutoff_m
double cutoff_m
Definition: Distribution.h:588
atan2
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
Definition: IpplExpressions.h:244
EnvelopeBunch::mySliceEndOffset
size_t mySliceEndOffset()
Definition: EnvelopeBunch.h:197
Distribution::setupEmissionModel
void setupEmissionModel(PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:4266
Distribution::sigmaR_m
Vector_t sigmaR_m
Definition: Distribution.h:531
Distribution::getTOrZDist
std::vector< double > & getTOrZDist()
Definition: Distribution.cpp:3026
WARNMSG
#define WARNMSG(msg)
Definition: IpplInfo.h:398
Distribution::cathodeTemp_m
double cathodeTemp_m
Cathode material Fermi energy (eV).
Definition: Distribution.h:496
Distribution::addDistributions
void addDistributions()
Definition: Distribution.cpp:761
Distribution::checkIfEmitted
void checkIfEmitted()
Definition: Distribution.cpp:989
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
PartBunchBase::getInitialGamma
double getInitialGamma() const
Definition: PartBunchBase.hpp:1823
Vector_t
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6
Cyclotron::getFMLowE
virtual double getFMLowE() const
Definition: Cyclotron.cpp:336
Beamline
An abstract sequence of beam line components.
Definition: Beamline.h:37
IpplInfo::getComm
static MPI_Comm getComm()
Definition: IpplInfo.h:178
PartBunchBase::getLocalNum
size_t getLocalNum() const
Definition: PartBunchBase.hpp:2230
Distribution::currentEnergyBin_m
int currentEnergyBin_m
Definition: Distribution.h:476
Inform::precision
int precision() const
Definition: Inform.h:115
Distribution::converteVToBetaGamma
double converteVToBetaGamma(double valueIneV, double massIneV)
Definition: Distribution.cpp:1042
Distribution::laserEnergy_m
double laserEnergy_m
Cathode material work function (eV).
Definition: Distribution.h:494
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
Distribution::getMaxTOrZ
double getMaxTOrZ()
Definition: Distribution.cpp:2962
Distribution::laserImageName_m
std::string laserImageName_m
Definition: Distribution.h:540
Distribution::cutoffP_m
Vector_t cutoffP_m
Definition: Distribution.h:534
Distribution::setupParticleBins
void setupParticleBins(double massIneV, PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:4347
PartBunchBase::dt
ParticleAttrib< double > dt
Definition: PartBunchBase.h:493
Distribution::pyWrite_m
std::vector< double > pyWrite_m
Definition: Distribution.h:516
BeamSequence
The base class for all OPAL beam lines and sequences.
Definition: BeamSequence.h:32
PartBunchBase::getNumBunch
short getNumBunch() const
Definition: PartBunchBase.hpp:1554
Distribution::generateTransverseGauss
void generateTransverseGauss(size_t numberOfParticles)
Definition: Distribution.cpp:2777
Distribution::applyEmissModelNone
void applyEmissModelNone(double &pz)
Definition: Distribution.cpp:834
Distribution::printDistSurfAndCreate
void printDistSurfAndCreate(Inform &os) const
Definition: Distribution.cpp:3323
Distribution::getvVThermal
double getvVThermal()
Definition: Distribution.h:742
Distribution::selectRandomGenerator
gsl_qrng * selectRandomGenerator(std::string, unsigned int dimension)
Select and allocate gsl random number generator.
Definition: Distribution.cpp:3461
Distribution::generateLongFlattopT
void generateLongFlattopT(size_t numberOfParticles)
Definition: Distribution.cpp:2632
Distribution::writeOutFileHeader
void writeOutFileHeader()
Definition: Distribution.cpp:4481
Distribution::printEmissionModelNone
void printEmissionModelNone(Inform &os) const
Definition: Distribution.cpp:3366
Distribution::setupEmissionModelNone
void setupEmissionModelNone(PartBunchBase< double, 3 > *beam)
Definition: Distribution.cpp:4305
Distribution::getNumberOfParticlesInFile
size_t getNumberOfParticlesInFile(std::ifstream &inputFile)
Definition: Distribution.cpp:1069
H5PartWrapper::readStep
virtual void readStep(PartBunchBase< double, 3 > *, h5_ssize_t firstParticle, h5_ssize_t lastParticle)=0
Distribution::eInitThreshold_m
double eInitThreshold_m
Definition: Distribution.h:551
Distribution::currentEmissionTime_m
double currentEmissionTime_m
Definition: Distribution.h:475
OpalData::find
Object * find(const std::string &name)
Find entry.
Definition: OpalData.cpp:618
Object
The base class for all OPAL objects.
Definition: Object.h:48
Distribution::avrgpz_m
double avrgpz_m
Definition: Distribution.h:522
Distribution::paraFNB_m
double paraFNB_m
Definition: Distribution.h:562
Distribution::sigmaRise_m
double sigmaRise_m
Definition: Distribution.h:586
cos
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
Distribution::tRise_m
double tRise_m
time binned distribution with thermal energy
Definition: Distribution.h:584
PartBunchBase::setNumBunch
void setNumBunch(short n)
Definition: PartBunchBase.hpp:1546
Beam
The BEAM definition.
Definition: Beam.h:35
Distribution::createMatchedGaussDistribution
void createMatchedGaussDistribution(size_t numberOfParticles, double massIneV)
Definition: Distribution.cpp:1251
Attributes::makeBool
Attribute makeBool(const std::string &name, const std::string &help)
Make logical attribute.
Definition: Attributes.cpp:56
Distribution::clone
virtual Distribution * clone(const std::string &name)
Return a clone.
Definition: Distribution.cpp:277
Physics::q_e
constexpr double q_e
The elementary charge in As.
Definition: Physics.h:76
PartBunchBase::getQ
double getQ() const
Access to reference data.
Definition: PartBunchBase.hpp:1765
name
const std::string name
Definition: MaxNormRadialPeak.cpp:3
EnvelopeBunch::distributeSlices
void distributeSlices(int nSlice=101)
distributes nSlice amongst processors and initializes slices
Definition: EnvelopeBunch.cpp:464
Distribution::doRestartOpalT
void doRestartOpalT(PartBunchBase< double, 3 > *p, size_t Np, int restartStep, H5PartWrapper *h5wrapper)
Definition: Distribution.cpp:516
Distribution::createOpalE
void createOpalE(Beam *beam, std::vector< Distribution * > addedDistributions, EnvelopeBunch *envelopeBunch, double distCenter, double Bz0)
Definition: Distribution.cpp:1585
Distribution::writeOutFileInjection
void writeOutFileInjection()
Definition: Distribution.cpp:4662
Distribution::paraFNA_m
double paraFNA_m
Definition: Distribution.h:560
EnvelopeBunch::getChargePerParticle
double getChargePerParticle()
returns charge per slice
Definition: EnvelopeBunch.h:190
Distribution::energyBins_m
PartBins * energyBins_m
Definition: Distribution.h:482
Distribution::totalNumberEmittedParticles_m
size_t totalNumberEmittedParticles_m
Definition: Distribution.h:502
Distribution::applyEmissModelAstra
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector< double > &additionalRNs)
Definition: Distribution.cpp:823
PartBunchBase::Bin
ParticleAttrib< int > Bin
Definition: PartBunchBase.h:492
Communicate::raw_receive
virtual int raw_receive(char *, int, int &, int &)
Definition: Communicate.h:200
Distribution::getXDist
std::vector< double > & getXDist()
Definition: Distribution.cpp:3010
PartBunchBase::boundp
virtual void boundp()
Definition: PartBunchBase.hpp:537
Distribution::eraseTOrZDist
void eraseTOrZDist()
Definition: Distribution.cpp:2000
Distribution::shiftDistCoordinates
void shiftDistCoordinates(double massIneV)
Definition: Distribution.cpp:4428
H5PartWrapper::readHeader
virtual void readHeader()=0
Options::cloTuneOnly
bool cloTuneOnly
Do closed orbit and tune calculation only.
Definition: Options.cpp:87
Physics::epsilon_0
constexpr double epsilon_0
The permittivity of vacuum in As/Vm.
Definition: Physics.h:58
Attributes::makeRealArray
Attribute makeRealArray(const std::string &name, const std::string &help)
Create real array attribute.
Definition: Attributes.cpp:253
Object::builtin
bool builtin
Built-in flag.
Definition: Object.h:231
Distribution::workFunction_m
double workFunction_m
Definition: Distribution.h:553
PartBunchBase::Bf
ParticleAttrib< Vector_t > Bf
Definition: PartBunchBase.h:491
Distribution::getMinTOrZ
double getMinTOrZ()
Definition: Distribution.cpp:2976
Distribution::setDistParametersFlattop
void setDistParametersFlattop(double massIneV)
Definition: Distribution.cpp:4046
asin
PETE_TUTree< FnArcSin, typename T::PETE_Expr_t > asin(const PETE_Expr< T > &l)
Definition: PETE.h:809
Distribution::distrTypeT_m
DistrTypeT::DistrTypeT distrTypeT_m
Distribution type. Declared as string.
Definition: Distribution.h:455
PartBunchBase::R
ParticlePos_t & R
Definition: PartBunchBase.h:479
Distribution::cathodeFermiEnergy_m
double cathodeFermiEnergy_m
Laser photon energy (eV).
Definition: Distribution.h:495
Distribution::applyEmissModelNonEquil
void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector< double > &additionalRNs)
Definition: Distribution.cpp:838
Distribution::getLastEmittedEnergyBin
int getLastEmittedEnergyBin()
Definition: Distribution.cpp:2958
PartBunchBase::iterateEmittedBin
void iterateEmittedBin(int binNumber)
Definition: PartBunchBase.hpp:1860
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
mslang::iterator
std::string::iterator iterator
Definition: MSLang.h:16
Distribution::fillEBinHistogram
void fillEBinHistogram()
Definition: Distribution.cpp:2008
Distribution::createDistributionBinomial
void createDistributionBinomial(size_t numberOfParticles, double massIneV)
Definition: Distribution.cpp:1050
Attributes::getReal
double getReal(const Attribute &attr)
Return real value.
Definition: Attributes.cpp:217
Distribution::getBGyDist
std::vector< double > & getBGyDist()
Definition: Distribution.cpp:3022
Distribution::fieldThrFN_m
double fieldThrFN_m
Definition: Distribution.h:556
Distribution::setDistToEmitted
void setDistToEmitted(bool emitted)
Definition: Distribution.cpp:3877
IpplTimings::stopTimer
static void stopTimer(TimerRef t)
Definition: IpplTimings.h:192
Options::ppdebug
bool ppdebug
ppdebug flag.
Definition: Options.cpp:13
OpalData::getInputBasename
std::string getInputBasename()
get input file name without extension
Definition: OpalData.cpp:717
PM
#define PM(a, b, c, d)
Definition: fftpack.cpp:96
Distribution::tPulseLengthFWHM_m
double tPulseLengthFWHM_m
Definition: Distribution.h:530
Distribution::printDistFlattop
void printDistFlattop(Inform &os) const
Definition: Distribution.cpp:3105
EnvelopeBunch::mySliceStartOffset
size_t mySliceStartOffset()
Definition: EnvelopeBunch.h:196
Distribution::pzDist_m
std::vector< double > pzDist_m
Definition: Distribution.h:510
Inform
Definition: Inform.h:41
Attributes::makeString
Attribute makeString(const std::string &name, const std::string &help)
Make string attribute.
Definition: Attributes.cpp:296
Distribution::yDist_m
std::vector< double > yDist_m
Definition: Distribution.h:507
Distribution::mBinomial_m
Vector_t mBinomial_m
Definition: Distribution.h:535
Distribution::setDistParametersGauss
void setDistParametersGauss(double massIneV)
Definition: Distribution.cpp:4145
H5PartWrapper::getNumParticles
size_t getNumParticles() const
Definition: H5PartWrapper.cpp:518
Distribution::getTEmission
double getTEmission()
Definition: Distribution.cpp:660
Distribution::pyDist_m
std::vector< double > pyDist_m
Definition: Distribution.h:508
IpplInfo::Comm
static Communicate * Comm
Definition: IpplInfo.h:93
Distribution::generateFlattopT
void generateFlattopT(size_t numberOfParticles)
Definition: Distribution.cpp:2343
Attributes::makeReal
Attribute makeReal(const std::string &name, const std::string &help)
Make real attribute.
Definition: Attributes.cpp:205
level1
Inform & level1(Inform &inf)
Definition: Inform.cpp:45
Distribution::checkParticleNumber
void checkParticleNumber(size_t &numberOfParticles)
Definition: Distribution.cpp:1004
Distribution::reflectDistribution
void reflectDistribution(size_t &numberOfParticles)
Definition: Distribution.cpp:3410
Distribution::printDistBinomial
void printDistBinomial(Inform &os) const
Definition: Distribution.cpp:3077
Distribution::laserIntensityCut_m
double laserIntensityCut_m
Definition: Distribution.h:541
PartBinsCyc.h
Defines a structure to hold energy bins and their associated data for multi-bunch tracking in cyclotr...
Distribution::createOpalT
void createOpalT(PartBunchBase< double, 3 > *beam, std::vector< Distribution * > addedDistributions, size_t &numberOfParticles)
Definition: Distribution.cpp:1651
Distribution::fieldEnhancement_m
double fieldEnhancement_m
Work function of surface material (eV).
Definition: Distribution.h:554
Distribution::printEmissionModel
void printEmissionModel(Inform &os) const
Definition: Distribution.cpp:3336
Distribution::printEmissionModelAstra
void printEmissionModelAstra(Inform &os) const
Definition: Distribution.cpp:3359
Distribution::setEmissionTime
void setEmissionTime(double &maxT, double &minT)
Definition: Distribution.cpp:3922
PartBunchBase::setPBins
void setPBins(PartBins *pbin)
Definition: PartBunchBase.hpp:307
Distribution::energyBinHist_m
gsl_histogram * energyBinHist_m
Distribution energy bins.
Definition: Distribution.h:483
Distribution::pxDist_m
std::vector< double > pxDist_m
Definition: Distribution.h:506
floor
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:816
Distribution::distT_m
std::string distT_m
Definition: Distribution.h:454
Distribution::vVThermal_m
double vVThermal_m
Velocity scalar for parallel plate benchmark.
Definition: Distribution.h:575
IPPL_APP_TAG2
#define IPPL_APP_TAG2
Definition: Tags.h:105
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
SigmaGenerator
This class computes the matched distribution.
Definition: SigmaGenerator.h:64
Distribution::getIfDistEmitting
bool getIfDistEmitting()
Definition: Distribution.cpp:2954
Distribution::paraFNY_m
double paraFNY_m
Definition: Distribution.h:564
Distribution::correlationMatrix_m
SymTenzor< double, 6 > correlationMatrix_m
Definition: Distribution.h:536
Distribution::tFall_m
double tFall_m
Definition: Distribution.h:585
Attributes::getString
std::string getString(const Attribute &attr)
Get string value.
Definition: Attributes.cpp:307
PartBunchBase::getGlobalTrackStep
long long getGlobalTrackStep() const
Definition: PartBunchBase.hpp:1522
Distribution::getNumOfLocalParticlesToCreate
size_t getNumOfLocalParticlesToCreate(size_t n)
Definition: Distribution.cpp:259
Distribution::binWrite_m
std::vector< size_t > binWrite_m
Definition: Distribution.h:519
Generated on Sat Jan 4 2020 14:48:00 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5