Distribution.cpp 179 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12
// ------------------------------------------------------------------------
// $RCSfile: Distribution.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.3.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Distribution
//   The class for the OPAL Distribution command.
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
13 14

#include "Distribution/Distribution.h"
15
#include "Distribution/SigmaGenerator.h"
16
#include "AbsBeamline/SpecificElementVisitor.h"
17

18 19 20 21 22 23
#include <cmath>
#include <cfloat>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
24
#include <numeric>
25

26 27
#include "AbstractObjects/Expressions.h"
#include "Attributes/Attributes.h"
28
#include "Utilities/Options.h"
adelmann's avatar
adelmann committed
29
#include "Utilities/OpalException.h"
30 31 32 33 34
#include "halton1d_sequence.hh"
#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBunch.h"
#include "Algorithms/PartBins.h"
#include "Algorithms/bet/EnvelopeBunch.h"
35
#include "Structure/Beam.h"
36
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
37 38 39
#include "Algorithms/PartBinsCyc.h"
#include "BasicActions/Option.h"
#include "Distribution/LaserProfile.h"
40
#include "Elements/OpalBeamline.h"
41
#include "AbstractObjects/BeamSequence.h"
42 43
#include "Structure/H5PartWrapper.h"
#include "Structure/H5PartWrapperForPC.h"
44 45 46 47

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_erf.h>
48 49
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
gsell's avatar
gsell committed
50

51

gsell's avatar
gsell committed
52 53
extern Inform *gmsg;

kraus's avatar
kraus committed
54
#define DISTDBG1
adelmann's avatar
Cleanup  
adelmann committed
55
#define noDISTDBG2
56

kraus's avatar
kraus committed
57 58
#define SMALLESTCUTOFF 1e-12

59 60 61 62 63 64 65 66
SymTenzor<double, 6> getUnit6x6() {
    SymTenzor<double, 6> unit6x6;
    for (unsigned int i = 0; i < 6u; ++ i) {
        unit6x6(i,i) = 1.0;
    }
    return unit6x6;
}

gsell's avatar
gsell committed
67 68 69 70
//
// Class Distribution
// ------------------------------------------------------------------------

71 72 73
namespace AttributesT
{
    enum AttributesT {
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
        DISTRIBUTION,
        FNAME,
        WRITETOFILE,
        WEIGHT,
        INPUTMOUNITS,
        EMITTED,
        EMISSIONSTEPS,
        EMISSIONMODEL,
        EKIN,
        ELASER,
        W,
        FE,
        CATHTEMP,
        NBIN,
        XMULT,
        YMULT,
        ZMULT,
        TMULT,
        PXMULT,
        PYMULT,
        PZMULT,
        OFFSETX,
        OFFSETY,
        OFFSETZ,
        OFFSETT,
        OFFSETPX,
        OFFSETPY,
        OFFSETPZ,
        SIGMAX,
        SIGMAY,
        SIGMAR,
        SIGMAZ,
        SIGMAT,
        TPULSEFWHM,
        TRISE,
        TFALL,
        SIGMAPX,
        SIGMAPY,
        SIGMAPZ,
        MX,
        MY,
        MZ,
        MT,
        CUTOFFX,
        CUTOFFY,
        CUTOFFR,
        CUTOFFLONG,
        CUTOFFPX,
        CUTOFFPY,
        CUTOFFPZ,
        FTOSCAMPLITUDE,
        FTOSCPERIODS,
        R,                          // the correlation matrix (a la transport)
        CORRX,
        CORRY,
        CORRZ,
        CORRT,
        R51,
        R52,
        R61,
        R62,
        LASERPROFFN,
        IMAGENAME,
        INTENSITYCUT,
138 139 140 141 142
        FLIPX,
        FLIPY,
        ROTATE90,
        ROTATE180,
        ROTATE270,
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
        NPDARKCUR,
        INWARDMARGIN,
        EINITHR,
        FNA,
        FNB,
        FNY,
        FNVYZERO,
        FNVYSECOND,
        FNPHIW,
        FNBETA,
        FNFIELDTHR,
        FNMAXEMI,
        SECONDARYFLAG,
        NEMISSIONMODE,
        VSEYZERO,                   // sey_0 in Vaughn's model.
        VEZERO,                     // Energy related to sey_0 in Vaughan's model.
        VSEYMAX,                    // sey max in Vaughan's model.
        VEMAX,                      // Emax in Vaughan's model.
        VKENERGY,                   // Fitting parameter denotes the roughness of
        // surface for impact energy in Vaughn's model.
        VKTHETA,                    // Fitting parameter denotes the roughness of
        // surface for impact angle in Vaughn's model.
        VVTHERMAL,                  // Thermal velocity of Maxwellian distribution
        // of secondaries in Vaughan's model.
        VW,
        SURFMATERIAL,               // Add material type, currently 0 for copper
        // and 1 for stainless steel.
        EX,                         // below is for the matched distribution
        EY,
        ET,
        MAGSYM,                     // number of sector magnets
        LINE,
        FMAPFN,
frey_m's avatar
frey_m committed
176
        FMTYPE,                     // field map type used in matched gauss distribution
177 178 179 180 181
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
Andreas Adelmann's avatar
Andreas Adelmann committed
182
	RGUESS,
183 184
        ID1,                       // special particle that the user can set
        ID2,                       // special particle that the user can set
185
        SIZE
gsell's avatar
gsell committed
186 187 188
    };
}

189 190 191
namespace LegacyAttributesT
{
    enum LegacyAttributesT {
192
        // DESCRIPTION OF THE DISTRIBUTION:
193
        DEBIN = AttributesT::SIZE,
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
        SBIN,
        TEMISSION,
        SIGLASER,
        AG,
        SIGMAPT,
        TRANSVCUTOFF,
        CUTOFF,
        Z,
        T,
        PT,
        ALPHAX,
        ALPHAY,
        BETAX,
        BETAY,
        DX,
        DDX,
        DY,
        DDY,
        SIZE
    };
214 215
}

gsell's avatar
gsell committed
216
Distribution::Distribution():
217 218
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
219
    distrTypeT_m(DistrTypeT::NODIST),
220
    numberOfDistributions_m(1),
221 222 223 224 225 226
    emitting_m(false),
    scan_m(false),
    emissionModel_m(EmissionModelT::NONE),
    tEmission_m(0.0),
    tBin_m(0.0),
    currentEmissionTime_m(0.0),
kraus's avatar
kraus committed
227
    currentEnergyBin_m(1),
228
    currentSampleBin_m(0),
229 230 231 232 233
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
    randGenEmit_m(NULL),
234
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
235
    pmean_m(0.0),
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
    inputMoUnits_m(InputMomentumUnitsT::NONE),
    sigmaTRise_m(0.0),
    sigmaTFall_m(0.0),
    tPulseLengthFWHM_m(0.0),
    correlationMatrix_m(getUnit6x6()),
    laserProfileFileName_m(""),
    laserImageName_m(""),
    laserIntensityCut_m(0.0),
    laserProfile_m(NULL),
    darkCurrentParts_m(0),
    darkInwardMargin_m(0.0),
    eInitThreshold_m(0.0),
    workFunction_m(0.0),
    fieldEnhancement_m(0.0),
    fieldThrFN_m(0.0),
    maxFN_m(0),
    paraFNA_m(0.0),
    paraFNB_m(0.0),
    paraFNY_m(0.0),
    paraFNVYSe_m(0.0),
    paraFNVYZe_m(0.0),
    secondaryFlag_m(0),
    ppVw_m(0.0),
    vVThermal_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    avrgpz_m(0.0),
    I_m(0.0),
    E_m(0.0),
    bega_m(0.0),
    M_m(0.0)
272
{
273 274 275 276
    SetAttributes();

    Distribution *defaultDistribution = clone("UNNAMED_Distribution");
    defaultDistribution->builtin = true;
gsell's avatar
gsell committed
277 278

    try {
279
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
280
    } catch(...) {
281 282 283 284
        delete defaultDistribution;
    }

    SetFieldEmissionParameters();
gsell's avatar
gsell committed
285 286 287 288 289 290 291
}
/**
 *
 *
 * @param name
 * @param parent
 */
292
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
293
    Definition(name, parent),
294 295
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
296
    numberOfDistributions_m(parent->numberOfDistributions_m),
297 298 299 300 301 302
    emitting_m(parent->emitting_m),
    scan_m(parent->scan_m),
    particleRefData_m(parent->particleRefData_m),
    addedDistributions_m(parent->addedDistributions_m),
    particlesPerDist_m(parent->particlesPerDist_m),
    emissionModel_m(parent->emissionModel_m),
gsell's avatar
gsell committed
303
    tEmission_m(parent->tEmission_m),
304 305 306 307 308 309 310 311 312 313
    tBin_m(parent->tBin_m),
    currentEmissionTime_m(parent->currentEmissionTime_m),
    currentEnergyBin_m(parent->currentEmissionTime_m),
    currentSampleBin_m(parent->currentSampleBin_m),
    numberOfEnergyBins_m(parent->numberOfEnergyBins_m),
    numberOfSampleBins_m(parent->numberOfSampleBins_m),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
    randGenEmit_m(parent->randGenEmit_m),
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
314
    pmean_m(parent->pmean_m),
315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
    cathodeWorkFunc_m(parent->cathodeWorkFunc_m),
    laserEnergy_m(parent->laserEnergy_m),
    cathodeFermiEnergy_m(parent->cathodeFermiEnergy_m),
    cathodeTemp_m(parent->cathodeTemp_m),
    emitEnergyUpperLimit_m(parent->emitEnergyUpperLimit_m),
    xDist_m(parent->xDist_m),
    pxDist_m(parent->pxDist_m),
    yDist_m(parent->yDist_m),
    pyDist_m(parent->pyDist_m),
    tOrZDist_m(parent->tOrZDist_m),
    pzDist_m(parent->pzDist_m),
    xWrite_m(parent->xWrite_m),
    pxWrite_m(parent->pxWrite_m),
    yWrite_m(parent->yWrite_m),
    pyWrite_m(parent->pyWrite_m),
    tOrZWrite_m(parent->tOrZWrite_m),
    pzWrite_m(parent->pzWrite_m),
332
    avrgpz_m(parent->avrgpz_m),
333 334 335 336 337 338 339 340
    inputMoUnits_m(parent->inputMoUnits_m),
    sigmaTRise_m(parent->sigmaTRise_m),
    sigmaTFall_m(parent->sigmaTFall_m),
    tPulseLengthFWHM_m(parent->tPulseLengthFWHM_m),
    sigmaR_m(parent->sigmaR_m),
    sigmaP_m(parent->sigmaP_m),
    cutoffR_m(parent->cutoffR_m),
    cutoffP_m(parent->cutoffP_m),
341
    correlationMatrix_m(parent->correlationMatrix_m),
342 343 344 345
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
346 347 348 349 350 351 352 353 354 355 356 357
    darkCurrentParts_m(parent->darkCurrentParts_m),
    darkInwardMargin_m(parent->darkInwardMargin_m),
    eInitThreshold_m(parent->eInitThreshold_m),
    workFunction_m(parent->workFunction_m),
    fieldEnhancement_m(parent->fieldEnhancement_m),
    fieldThrFN_m(parent->fieldThrFN_m),
    maxFN_m(parent->maxFN_m),
    paraFNA_m(parent-> paraFNA_m),
    paraFNB_m(parent-> paraFNB_m),
    paraFNY_m(parent-> paraFNY_m),
    paraFNVYSe_m(parent-> paraFNVYSe_m),
    paraFNVYZe_m(parent-> paraFNVYZe_m),
358 359 360 361 362 363 364
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
365 366
    cutoff_m(parent->cutoff_m),
    I_m(parent->I_m),
adelmann's avatar
adelmann committed
367 368 369
    E_m(parent->E_m),
    bega_m(parent->bega_m),
    M_m(parent->M_m)
370
{
gsell's avatar
gsell committed
371 372 373 374 375 376 377
}

Distribution::~Distribution() {

    if((Ippl::getNodes() == 1) && (os_m.is_open()))
        os_m.close();

378 379 380
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
381 382
    }

383 384 385 386 387 388
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

    if (randGenEmit_m != NULL) {
389
        gsl_rng_free(randGenEmit_m);
390 391 392 393 394 395 396
        randGenEmit_m = NULL;
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
397 398

}
399
/*
400
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
401
  for(int i=0; i<M.size1(); ++i) {
402 403 404 405 406
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
407 408
  }
*/
gsell's avatar
gsell committed
409 410

/**
411 412
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
413 414 415
 * @param
 * @param
 * @param
416
 */
417
void Distribution::WriteToFile() {
418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
    /*
      if(Ippl::getNodes() == 1) {
      if(os_m.is_open()) {
      ;
      } else {
      *gmsg << " Write distribution to file data/dist.dat" << endl;
      std::string file("data/dist.dat");
      os_m.open(file.c_str());
      if(os_m.bad()) {
      *gmsg << "Unable to open output file " <<  file << endl;
      }
      os_m << "# x y ti px py pz "  << std::endl;
      os_m.close();
      }
      }
    */
434 435
}

436 437 438
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
439 440
}

441 442 443
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
444

445 446
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
447

448 449
void Distribution::update() {
}
adelmann's avatar
adelmann committed
450

451
void Distribution::Create(size_t &numberOfParticles, double massIneV) {
adelmann's avatar
adelmann committed
452

453
    SetFieldEmissionParameters();
adelmann's avatar
adelmann committed
454

455
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
456

457
    case DistrTypeT::MATCHEDGAUSS:
458
        CreateMatchedGaussDistribution(numberOfParticles, massIneV);
459
        break;
460 461
    case DistrTypeT::FROMFILE:
        CreateDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
462
        break;
463 464
    case DistrTypeT::GAUSS:
        CreateDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
465
        break;
466 467
    case DistrTypeT::BINOMIAL:
        CreateDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
468
        break;
469 470
    case DistrTypeT::FLATTOP:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
471
        break;
472 473
    case DistrTypeT::GUNGAUSSFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
474
        break;
475 476
    case DistrTypeT::ASTRAFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
477
        break;
478 479 480 481
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516

    if (emitting_m) {
        unsigned int numAdditionalRNsPerParticle = 0;
        if (emissionModel_m == EmissionModelT::ASTRA) {
            numAdditionalRNsPerParticle = 2;
        }

        if (Options::cZero) {
            numAdditionalRNsPerParticle *= 2;
        }

        int saveProcessor = -1;

        for (size_t partIndex = 0; partIndex < numberOfParticles; ++ partIndex) {

            // Save to each processor in turn.
            ++ saveProcessor;
            if (saveProcessor >= Ippl::getNodes())
                saveProcessor = 0;

            if (Ippl::myNode() == saveProcessor) {
                std::vector<double> rns;
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
                    double x = gsl_rng_uniform(randGenEmit_m);
                    rns.push_back(x);
                }
                additionalRNs_m.push_back(rns);
            } else {
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
                    gsl_rng_uniform(randGenEmit_m);
                }
            }
        }
    }

517
    // AAA Scale and shift coordinates according to distribution input.
Andreas Adelmann's avatar
Andreas Adelmann committed
518
    ScaleDistCoordinates();
519
}
gsell's avatar
gsell committed
520

521
void  Distribution::CreatePriPart(PartBunch *beam, BoundaryGeometry &bg) {
gsell's avatar
gsell committed
522

523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
        double vw = this->GetVw();
        double vt = this->GetvVThermal();
        double f_max = vw / vt * exp(-0.5);
        double test_a = vt / vw;
        double test_asq = test_a * test_a;
        size_t count = 0;
        size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
        size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
        if(Ippl::myNode() == 0)
            N_mean += N_extra;
        if(bg.getN() != 0) {
            for(size_t i = 0; i < bg.getN(); i++) {
                if(pc == Ippl::myNode()) {
                    if(count < N_mean) {
                        /*==============Parallel Plate Benchmark=====================================*/
                        double test_s = 1;
                        double f_x = 0;
                        double test_x = 0;
                        while(test_s > f_x) {
                            test_s = IpplRandom();
                            test_s *= f_max;
                            test_x = IpplRandom();
                            test_x *= 10 * test_a; //range for normalized emission speed(0,10*test_a);
                            f_x = test_x / test_asq * exp(-test_x * test_x / 2 / test_asq);
                        }
                        double v_emi = test_x * vw;
gsell's avatar
gsell committed
552

553 554 555 556 557 558 559
                        double betaemit = v_emi / Physics::c;
                        double betagamma = betaemit / sqrt(1 - betaemit * betaemit);
                        /*============================================================================ */
                        beam->create(1);
                        if(pc != 0) {
                            beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(Ippl::myNode() * N_mean + count);
gsell's avatar
gsell committed
560
                        } else {
561 562
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
563
                        }
564
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
565
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
566 567 568 569 570 571 572
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        beam->LastSection[lowMark + count] = 0;
                        beam->Ef[lowMark + count] = Vector_t(0.0);
                        beam->Bf[lowMark + count] = Vector_t(0.0);
                        beam->dt[lowMark + count] = beam->getdT();
                        count ++;
gsell's avatar
gsell committed
573 574
                    }
                }
575 576 577
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
578
            }
579 580 581
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
582
        }
583
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
584

585
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
586

587 588 589 590 591
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
        size_t count = 0;
        size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
        size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
gsell's avatar
gsell committed
592

593 594 595 596
        if(Ippl::myNode() == 0)
            N_mean += N_extra;
        if(bg.getN() != 0) {
            for(size_t i = 0; i < bg.getN(); i++) {
gsell's avatar
gsell committed
597

598 599 600 601 602 603 604 605 606
                if(pc == Ippl::myNode()) {
                    if(count < N_mean) {
                        beam->create(1);
                        if(pc != 0)
                            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
                        else
                            beam->R[lowMark + count] = bg.getCooridinate(count); // for node0 the particle number N_mean =  N_mean + N_extra
                        beam->P[lowMark + count] = Vector_t(0.0);
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
607
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
608 609 610 611 612 613 614
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        beam->LastSection[lowMark + count] = 0;
                        beam->Ef[lowMark + count] = Vector_t(0.0);
                        beam->Bf[lowMark + count] = Vector_t(0.0);
                        beam->dt[lowMark + count] = beam->getdT();
                        count++;
gsell's avatar
gsell committed
615 616

                    }
617

gsell's avatar
gsell committed
618
                }
619 620 621
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
622

623
            }
624

625 626 627 628 629 630
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
631

632
void Distribution::DoRestartOpalT(PartBunch &beam, size_t Np, int restartStep, H5PartWrapper *dataSource) {
633

634
    IpplTimings::startTimer(beam.distrReload_m);
635

636 637
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
638

639 640 641 642
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
643

644 645
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
646

647
    beam.create(numParticles);
648

649 650
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
651 652 653

    beam.boundp();

654 655 656 657
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
658 659
    IpplTimings::stopTimer(beam.distrReload_m);

660 661 662 663
    *gmsg << "Total number of particles in the h5 file= " << beam.getTotalNum() << "\n"
          << "Global step= " << gtstep << "; Local step= " << ltstep << "; "
          << "restart step= " << restartStep << "\n"
          << "time of restart= " << actualT << "; phishift= " << OpalData::getInstance()->getGlobalPhaseShift() << endl;
gsell's avatar
gsell committed
664 665
}

666 667 668 669 670
void Distribution::DoRestartOpalCycl(PartBunch &beam,
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
671

672
    // h5_int64_t rc;
673
    IpplTimings::startTimer(beam.distrReload_m);
674
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
675

676 677
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
678

679 680 681 682
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
683

684 685
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
686

687
    beam.create(numParticles);
688

689 690
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
691

692
    beam.Q = beam.getChargePerParticle();
693

694
    beam.boundp();
695

696
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
697

698 699 700 701 702 703
    const int globalN = beam.getTotalNum();
    INFOMSG("Restart from hdf5 format file " << OpalData::getInstance()->getInputBasename() << ".h5" << endl);
    INFOMSG("total number of particles = " << globalN << endl);
    INFOMSG("* Restart Energy = " << meanE << " (MeV), Path lenght = " << beam.getLPath() << " (m)" <<  endl);
    INFOMSG("Tracking Step since last bunch injection is " << beam.getSteptoLastInj() << endl);
    INFOMSG(beam.getNumBunch() << " Bunches(bins) exist in this file" << endl);
704

705 706
    double gamma = 1 + meanE / beam.getM() * 1.0e6;
    double beta = sqrt(1.0 - (1.0 / std::pow(gamma, 2.0)));
707

708
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
709

710 711
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
712 713
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
714
            beam.setPBins(new PartBinsCyc(1000, beam.getNumBunch()));
715 716
        }

717 718
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
719 720 721 722 723 724 725 726 727 728 729 730 731 732

        Vector_t meanR(0.0, 0.0, 0.0);
        Vector_t meanP(0.0, 0.0, 0.0);
        unsigned long int newLocalN = beam.getLocalNum();
        for(unsigned int i = 0; i < newLocalN; ++i) {
            for(int d = 0; d < 3; ++d) {
                meanR(d) += beam.R[i](d);
                meanP(d) += beam.P[i](d);
            }
        }
        reduce(meanR, meanR, OpAddAssign());
        meanR /= Vector_t(globalN);
        reduce(meanP, meanP, OpAddAssign());
        meanP /= Vector_t(globalN);
733
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
734 735 736 737 738 739 740

        for(unsigned int i = 0; i < newLocalN; ++i) {
            beam.R[i] -= meanR;
            beam.P[i] -= meanP;
        }
    }

741
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
742 743 744
    IpplTimings::stopTimer(beam.distrReload_m);
}

745 746
void Distribution::DoRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
                                  H5PartWrapper *dataSource) {
747
    IpplTimings::startTimer(beam.distrReload_m);
748 749
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
750 751 752 753 754 755

    beam.distributeSlices(N);
    beam.createBunch();
    long long starti = beam.mySliceStartOffset();
    long long endi = beam.mySliceEndOffset();

756 757
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788

    beam.setCharge(beam.getChargePerParticle());
    IpplTimings::stopTimer(beam.distrReload_m);
}

Distribution *Distribution::find(const std::string &name) {
    Distribution *dist = dynamic_cast<Distribution *>(OpalData::getInstance()->find(name));

    if(dist == 0) {
        throw OpalException("Distribution::find()", "Distribution \"" + name + "\" not found.");
    }

    return dist;
}

double Distribution::GetTEmission() {
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

    distT_m = Attributes::getString(itsAttr[AttributesT::DISTRIBUTION]);
    if(distT_m == "GAUSS")
        distrTypeT_m = DistrTypeT::GAUSS;
    else if(distT_m == "GUNGAUSSFLATTOPTH")
        distrTypeT_m = DistrTypeT::GUNGAUSSFLATTOPTH;
    else if(distT_m == "FROMFILE")
        distrTypeT_m = DistrTypeT::FROMFILE;
    else if(distT_m == "BINOMIAL")
        distrTypeT_m = DistrTypeT::BINOMIAL;

    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[AttributesT::TPULSEFWHM]);
789
    cutoff_m = Attributes::getReal(itsAttr[ LegacyAttributesT::CUTOFF]);
790 791 792 793 794 795 796
    tRise_m = Attributes::getReal(itsAttr[AttributesT::TRISE]);
    tFall_m = Attributes::getReal(itsAttr[AttributesT::TFALL]);
    double tratio = sqrt(2.0 * log(10.0)) - sqrt(2.0 * log(10.0 / 9.0));
    sigmaRise_m = tRise_m / tratio;
    sigmaFall_m = tFall_m / tratio;

    switch(distrTypeT_m) {
797 798 799 800 801 802 803 804 805 806 807 808 809 810
    case DistrTypeT::ASTRAFLATTOPTH: {
        double a = tPulseLengthFWHM_m / 2;
        double sig = tRise_m / 2;
        double inv_erf08 = 0.906193802436823; // erfinv(0.8)
        double sqr2 = sqrt(2.);
        double t = a - sqr2 * sig * inv_erf08;
        double tmps = sig;
        double tmpt = t;
        for(int i = 0; i < 10; ++ i) {
            sig = (t + tRise_m - a) / (sqr2 * inv_erf08);
            t = a - 0.5 * sqr2 * (sig + tmps) * inv_erf08;
            sig = (0.5 * (t + tmpt) + tRise_m - a) / (sqr2 * inv_erf08);
            tmps = sig;
            tmpt = t;
811
        }
812 813 814 815 816 817 818 819 820
        tEmission_m = tPulseLengthFWHM_m + 10 * sig;
        break;
    }
    case DistrTypeT::GUNGAUSSFLATTOPTH: {
        tEmission_m = tPulseLengthFWHM_m + (cutoff_m - sqrt(2.0 * log(2.0))) * (sigmaRise_m + sigmaFall_m);
        break;
    }
    default:
        tEmission_m = 0.0;
821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858
    }
    return tEmission_m;
}

double Distribution::GetEkin() const {return Attributes::getReal(itsAttr[AttributesT::EKIN]);}
double Distribution::GetLaserEnergy() const {return Attributes::getReal(itsAttr[AttributesT::ELASER]);}
double Distribution::GetWorkFunctionRf() const {return Attributes::getReal(itsAttr[AttributesT::W]);}

size_t Distribution::GetNumberOfDarkCurrentParticles() { return (size_t) Attributes::getReal(itsAttr[AttributesT::NPDARKCUR]);}
double Distribution::GetDarkCurrentParticlesInwardMargin() { return Attributes::getReal(itsAttr[AttributesT::INWARDMARGIN]);}
double Distribution::GetEInitThreshold() { return Attributes::getReal(itsAttr[AttributesT::EINITHR]);}
double Distribution::GetWorkFunction() { return Attributes::getReal(itsAttr[AttributesT::FNPHIW]); }
double Distribution::GetFieldEnhancement() { return Attributes::getReal(itsAttr[AttributesT::FNBETA]); }
size_t Distribution::GetMaxFNemissionPartPerTri() { return (size_t) Attributes::getReal(itsAttr[AttributesT::FNMAXEMI]);}
double Distribution::GetFieldFNThreshold() { return Attributes::getReal(itsAttr[AttributesT::FNFIELDTHR]);}
double Distribution::GetFNParameterA() { return Attributes::getReal(itsAttr[AttributesT::FNA]);}
double Distribution::GetFNParameterB() { return Attributes::getReal(itsAttr[AttributesT::FNB]);}
double Distribution::GetFNParameterY() { return Attributes::getReal(itsAttr[AttributesT::FNY]);}
double Distribution::GetFNParameterVYZero() { return Attributes::getReal(itsAttr[AttributesT::FNVYZERO]);}
double Distribution::GetFNParameterVYSecond() { return Attributes::getReal(itsAttr[AttributesT::FNVYSECOND]);}
int    Distribution::GetSecondaryEmissionFlag() { return Attributes::getReal(itsAttr[AttributesT::SECONDARYFLAG]);}
bool   Distribution::GetEmissionMode() { return Attributes::getBool(itsAttr[AttributesT::NEMISSIONMODE]);}

std::string Distribution::GetTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::DISTRIBUTION]);}

double Distribution::GetvSeyZero() {return Attributes::getReal(itsAttr[AttributesT::VSEYZERO]);}// return sey_0 in Vaughan's model
double Distribution::GetvEZero() {return Attributes::getReal(itsAttr[AttributesT::VEZERO]);}// return the energy related to sey_0 in Vaughan's model
double Distribution::GetvSeyMax() {return Attributes::getReal(itsAttr[AttributesT::VSEYMAX]);}// return sey max in Vaughan's model
double Distribution::GetvEmax() {return Attributes::getReal(itsAttr[AttributesT::VEMAX]);}// return Emax in Vaughan's model
double Distribution::GetvKenergy() {return Attributes::getReal(itsAttr[AttributesT::VKENERGY]);}// return fitting parameter denotes the roughness of surface for impact energy in Vaughan's model
double Distribution::GetvKtheta() {return Attributes::getReal(itsAttr[AttributesT::VKTHETA]);}// return fitting parameter denotes the roughness of surface for impact angle in Vaughan's model
double Distribution::GetvVThermal() {return Attributes::getReal(itsAttr[AttributesT::VVTHERMAL]);}// thermal velocity of Maxwellian distribution of secondaries in Vaughan's model
double Distribution::GetVw() {return Attributes::getReal(itsAttr[AttributesT::VW]);}// velocity scalar for parallel plate benchmark;

int Distribution::GetSurfMaterial() {return (int)Attributes::getReal(itsAttr[AttributesT::SURFMATERIAL]);}// Surface material number for Furman-Pivi's Model;

Inform &Distribution::printInfo(Inform &os) const {

859 860
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
861
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
862
        os << "* In restart. Distribution read in from .h5 file." << endl;
863 864
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
865
            os << "* Main Distribution" << endl
866
               << "-----------------" << endl;
867 868 869 870 871 872 873 874

        if (particlesPerDist_m.empty())
            PrintDist(os, 0);
        else
            PrintDist(os, particlesPerDist_m.at(0));

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
875
            os << "* " << endl;
adelmann's avatar
adelmann committed
876 877
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
878 879 880 881
            addedDistributions_m.at(distIndex)->PrintDist(os, particlesPerDist_m.at(distCount));
            distCount++;
        }

882
        os << "* " << endl;
883
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
884
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
885

886
            //            if (numberOfEnergyBins_m > 1)
adelmann's avatar
adelmann committed
887
            //    PrintEnergyBins(os);
888 889 890
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
891 892
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
893 894
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
895
            os << "* " << endl;
896
            PrintEmissionModel(os);
897
            os << "* " << endl;
898
        } else
adelmann's avatar
adelmann committed
899
            os << "* Distribution is injected." << endl;
900
    }
901 902
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
903 904 905 906 907 908 909 910 911 912 913 914 915 916

    return os;
}

const PartData &Distribution::GetReference() const {
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

void Distribution::AddDistributions() {
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
917

918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
         addedDistIt != addedDistributions_m.end(); addedDistIt++) {

        std::vector<double>::iterator distIt;
        for (distIt = (*addedDistIt)->GetXDist().begin();
             distIt != (*addedDistIt)->GetXDist().end();
             distIt++) {
            xDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseXDist();

        for (distIt = (*addedDistIt)->GetBGxDist().begin();
             distIt != (*addedDistIt)->GetBGxDist().end();
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseBGxDist();

        for (distIt = (*addedDistIt)->GetYDist().begin();
             distIt != (*addedDistIt)->GetYDist().end();
             distIt++) {
            yDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseYDist();

        for (distIt = (*addedDistIt)->GetBGyDist().begin();
             distIt != (*addedDistIt)->GetBGyDist().end();
             distIt++) {
            pyDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseBGyDist();

        for (distIt = (*addedDistIt)->GetTOrZDist().begin();
             distIt != (*addedDistIt)->GetTOrZDist().end();
             distIt++) {
            tOrZDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseTOrZDist();

        for (distIt = (*addedDistIt)->GetBGzDist().begin();
             distIt != (*addedDistIt)->GetBGzDist().end();
             distIt++) {
            pzDist_m.push_back(*distIt);
        }
        (*addedDistIt)->EraseBGzDist();
    }
}

kraus's avatar
kraus committed
967
void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
968 969 970 971 972 973 974

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
        ApplyEmissModelNone(pz);
        break;
    case EmissionModelT::ASTRA:
kraus's avatar
kraus committed
975
        ApplyEmissModelAstra(px, py, pz, additionalRNs);
976 977 978 979 980 981 982 983 984
        break;
    case EmissionModelT::NONEQUIL:
        ApplyEmissModelNonEquil(eZ, px, py, pz);
        break;
    default:
        break;
    }
}

kraus's avatar
kraus committed
985
void Distribution::ApplyEmissModelAstra(double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
986

kraus's avatar
kraus committed
987 988
    double phi = 2.0 * acos(sqrt(additionalRNs[0]));
    double theta = 2.0 * Physics::pi * additionalRNs[1];
989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005

    px = pTotThermal_m * sin(phi) * cos(theta);
    py = pTotThermal_m * sin(phi) * sin(theta);
    pz = pTotThermal_m * std::abs(cos(phi));

}

void Distribution::ApplyEmissModelNone(double &pz) {
    pz += pTotThermal_m;
}

void Distribution::ApplyEmissModelNonEquil(double eZ,
                                           double &bgx,
                                           double &bgy,
                                           double &bgz) {

    double phiEffective = cathodeWorkFunc_m
1006 1007
        - sqrt(Physics::q_e * eZ /
               (4.0 * Physics::pi * Physics::epsilon_0));
1008 1009 1010 1011 1012 1013
    double lowEnergyLimit = cathodeFermiEnergy_m + phiEffective - laserEnergy_m;

    // Generate emission energy.
    double energy = 0.0;
    bool allow = false;
    while (!allow) {
1014 1015
        energy = lowEnergyLimit + (gsl_rng_uniform(randGenEmit_m)*emitEnergyUpperLimit_m);
        double randFuncValue = gsl_rng_uniform(randGenEmit_m);
1016 1017 1018 1019 1020
        double funcValue = (1.0
                            - 1.0
                            / (1.0
                               + exp((energy + laserEnergy_m - cathodeFermiEnergy_m)
                                     / (Physics::kB * cathodeTemp_m))))
1021 1022 1023
            / (1.0
               + exp((energy - cathodeFermiEnergy_m)
                     / (Physics::kB * cathodeTemp_m)));
1024 1025 1026 1027 1028 1029 1030
        if (randFuncValue <= funcValue)
            allow = true;
    }

    // Compute emission angles.
    double energyInternal = energy + laserEnergy_m;
    double energyExternal = energy + laserEnergy_m
1031
        - cathodeFermiEnergy_m - phiEffective;
1032 1033 1034

    double thetaInMax = acos(sqrt((cathodeFermiEnergy_m + phiEffective)
                                  / (energy + laserEnergy_m)));
1035
    double thetaIn = gsl_rng_uniform(randGenEmit_m)*thetaInMax;
1036
    double sinThetaOut = sin(thetaIn) * sqrt(energyInternal / energyExternal);
1037
    double phi = Physics::two_pi * gsl_rng_uniform(randGenEmit_m);
1038 1039 1040

    // Compute emission momenta.
    double betaGammaExternal
1041
        = sqrt(pow(energyExternal / (Physics::m_e * 1.0e9) + 1.0, 2.0) - 1.0);
1042 1043 1044 1045 1046 1047 1048 1049 1050

    bgx = betaGammaExternal * sinThetaOut * cos(phi);
    bgy = betaGammaExternal * sinThetaOut * sin(phi);
    bgz = betaGammaExternal * sqrt(1.0 - pow(sinThetaOut, 2.0));

}

void Distribution::CalcPartPerDist(size_t numberOfParticles) {

1051
    typedef std::vector<Distribution *>::iterator iterator;
1052

1053
    if (numberOfDistributions_m == 1)
1054 1055
        particlesPerDist_m.push_back(numberOfParticles);
    else {
1056 1057 1058 1059
        double totalWeight = GetWeight();
        for (iterator it = addedDistributions_m.begin(); it != addedDistributions_m.end(); it++) {
            totalWeight += (*it)->GetWeight();
        }
1060 1061 1062

        particlesPerDist_m.push_back(0);
        size_t numberOfCommittedPart = 0;
1063 1064 1065 1066
        for (iterator it = addedDistributions_m.begin(); it != addedDistributions_m.end(); it++) {
            size_t particlesCurrentDist = numberOfParticles * (*it)->GetWeight() / totalWeight;
            particlesPerDist_m.push_back(particlesCurrentDist);
            numberOfCommittedPart += particlesCurrentDist;
1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095
        }

        // Remaining particles go into main distribution.
        particlesPerDist_m.at(0) = numberOfParticles - numberOfCommittedPart;

    }

}

void Distribution::CheckEmissionParameters() {

    // There must be at least on energy bin for an emitted beam.
    numberOfEnergyBins_m
        = std::abs(static_cast<int> (Attributes::getReal(itsAttr[AttributesT::NBIN])));
    if (numberOfEnergyBins_m == 0)
        numberOfEnergyBins_m = 1;

    int emissionSteps = std::abs(static_cast<int> (Attributes::getReal(itsAttr[AttributesT::EMISSIONSTEPS])));

    // There must be at least one emission step.
    if (emissionSteps == 0)
        emissionSteps = 1;

    // Set number of sample bins per energy bin from the number of emission steps.
    numberOfSampleBins_m = static_cast<int> (std::ceil(emissionSteps / numberOfEnergyBins_m));
    if (numberOfSampleBins_m <= 0)
        numberOfSampleBins_m = 1;

    // Initialize emission counters.
kraus's avatar
kraus committed
1096
    currentEnergyBin_m = 1;
1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119
    currentSampleBin_m = 0;

}

void Distribution::CheckIfEmitted() {

    emitting_m = Attributes::getBool(itsAttr[AttributesT::EMITTED]);

    switch (distrTypeT_m) {

    case DistrTypeT::ASTRAFLATTOPTH:
        emitting_m = true;
        break;
    case DistrTypeT::GUNGAUSSFLATTOPTH:
        emitting_m = true;
        break;
    default:
        break;
    }
}

void Distribution::CheckParticleNumber(size_t &numberOfParticles) {

1120 1121
    size_t numberOfDistParticles = tOrZDist_m.size();
    reduce(numberOfDistParticles, numberOfDistParticles, OpAddAssign());
adelmann's avatar
adelmann committed
1122
    
1123
    if (numberOfDistParticles != numberOfParticles) {
1124
        *gmsg << "\n--------------------------------------------------" << endl
Andreas Adelmann's avatar
Andreas Adelmann committed
1125
              << "Error!! The number of particles in the initial" << endl
1126
              << "distribution is " << numberOfDistParticles << "." << endl << endl
1127
              << "This is different from the number of particles" << endl
1128
              << "defined by the BEAM command: " << numberOfParticles << endl << endl
1129 1130 1131 1132 1133
              << "This is often happens when using a FROMFILE type"