Distribution.cpp 175 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"
kraus's avatar
kraus committed
28
#include "Utilities/OpalOptions.h"
29 30 31 32 33 34
#include "Utilities/Options.h"
#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 52
#include "MagneticField.h"

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

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

58 59 60 61 62 63 64 65
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
66 67 68 69
//
// Class Distribution
// ------------------------------------------------------------------------

70 71 72
namespace AttributesT
{
    enum AttributesT {
73 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
        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,
137 138 139 140 141
        FLIPX,
        FLIPY,
        ROTATE90,
        ROTATE180,
        ROTATE270,
142 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 176 177 178 179 180
        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,
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
        SIZE
gsell's avatar
gsell committed
181 182 183
    };
}

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

gsell's avatar
gsell committed
211
Distribution::Distribution():
212 213
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
214
    distrTypeT_m(DistrTypeT::NODIST),
215
    numberOfDistributions_m(1),
216 217 218 219 220 221
    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
222
    currentEnergyBin_m(1),
223
    currentSampleBin_m(0),
224 225 226 227 228
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
    randGenEmit_m(NULL),
229
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
230
    pmean_m(0.0),
231 232 233 234 235 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
    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)
267
{
268 269 270 271
    SetAttributes();

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

    try {
274
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
275
    } catch(...) {
276 277 278 279
        delete defaultDistribution;
    }

    SetFieldEmissionParameters();
gsell's avatar
gsell committed
280 281 282 283 284 285 286
}
/**
 *
 *
 * @param name
 * @param parent
 */
287
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
288
    Definition(name, parent),
289 290
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
291
    numberOfDistributions_m(parent->numberOfDistributions_m),
292 293 294 295 296 297
    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
298
    tEmission_m(parent->tEmission_m),
299 300 301 302 303 304 305 306 307 308
    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
309
    pmean_m(parent->pmean_m),
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
    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),
327
    avrgpz_m(parent->avrgpz_m),
328 329 330 331 332 333 334 335
    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),
336
    correlationMatrix_m(parent->correlationMatrix_m),
337 338 339 340
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
341 342 343 344 345 346 347 348 349 350 351 352
    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),
353 354 355 356 357 358 359
    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),
360 361
    cutoff_m(parent->cutoff_m),
    I_m(parent->I_m),
adelmann's avatar
adelmann committed
362 363 364
    E_m(parent->E_m),
    bega_m(parent->bega_m),
    M_m(parent->M_m)
365
{
gsell's avatar
gsell committed
366 367 368 369 370 371 372
}

Distribution::~Distribution() {

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

373 374 375
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
376 377
    }

378 379 380 381 382 383
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

    if (randGenEmit_m != NULL) {
384
        gsl_rng_free(randGenEmit_m);
385 386 387 388 389 390 391
        randGenEmit_m = NULL;
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
392 393

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

/**
406 407
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
408 409 410
 * @param
 * @param
 * @param
411
 */
412
void Distribution::WriteToFile() {
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
    /*
      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();
      }
      }
    */
429 430
}

431 432 433
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
434 435
}

436 437 438
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
439

440 441
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
442

443 444
void Distribution::update() {
}
adelmann's avatar
adelmann committed
445

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

448
    SetFieldEmissionParameters();
adelmann's avatar
adelmann committed
449

450
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
451

452
    case DistrTypeT::MATCHEDGAUSS:
453
        CreateMatchedGaussDistribution(numberOfParticles, massIneV);
454
        break;
455 456
    case DistrTypeT::FROMFILE:
        CreateDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
457
        break;
458 459
    case DistrTypeT::GAUSS:
        CreateDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
460
        break;
461 462
    case DistrTypeT::BINOMIAL:
        CreateDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
463
        break;
464 465
    case DistrTypeT::FLATTOP:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
466
        break;
467 468
    case DistrTypeT::GUNGAUSSFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
469
        break;
470 471
    case DistrTypeT::ASTRAFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
472
        break;
473 474 475 476
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
gsell's avatar
gsell committed
477

478
    // AAA Scale and shift coordinates according to distribution input.
479 480 481
    ScaleDistCoordinates();
    ShiftDistCoordinates(massIneV);
}
gsell's avatar
gsell committed
482

483
void  Distribution::CreatePriPart(PartBunch *beam, BoundaryGeometry &bg) {
gsell's avatar
gsell committed
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
    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
514

515 516 517 518 519 520 521
                        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
522
                        } else {
523 524
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
525
                        }
526
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
527
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
528 529 530 531 532 533 534
                        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
535 536
                    }
                }
537 538 539
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
540
            }
541 542 543
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
544
        }
545
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
546

547
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
548

549 550 551 552 553
        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
554

555 556 557 558
        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
559

560 561 562 563 564 565 566 567 568
                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
569
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
570 571 572 573 574 575 576
                        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
577 578

                    }
579

gsell's avatar
gsell committed
580
                }
581 582 583
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
584

585
            }
586

587 588 589 590 591 592
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
593

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

596
    IpplTimings::startTimer(beam.distrReload_m);
597

598 599
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
600

601 602 603 604
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
605

606 607
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
608

609
    beam.create(numParticles);
610

611 612
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
613 614 615

    beam.boundp();

616 617 618 619
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
620 621
    IpplTimings::stopTimer(beam.distrReload_m);

622 623 624 625
    *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
626 627
}

628 629 630 631 632
void Distribution::DoRestartOpalCycl(PartBunch &beam,
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
633

634
    // h5_int64_t rc;
635
    IpplTimings::startTimer(beam.distrReload_m);
636
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
637

638 639
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
640

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

646 647
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
648

649
    beam.create(numParticles);
650

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

654
    beam.Q = beam.getChargePerParticle();
655

656
    beam.boundp();
657

658
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
659

660 661 662 663 664 665
    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);
666

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

670
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
671

672 673
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
674 675
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
676
            beam.setPBins(new PartBinsCyc(1000, beam.getNumBunch()));
677 678
        }

679 680
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
681 682 683 684 685 686 687 688 689 690 691 692 693 694

        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);
695
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
696 697 698 699 700 701 702

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

703
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
704 705 706
    IpplTimings::stopTimer(beam.distrReload_m);
}

707 708
void Distribution::DoRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
                                  H5PartWrapper *dataSource) {
709
    IpplTimings::startTimer(beam.distrReload_m);
710 711
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
712 713 714 715 716 717

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

718 719
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750

    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]);
751
    cutoff_m = Attributes::getReal(itsAttr[ LegacyAttributesT::CUTOFF]);
752 753 754 755 756 757 758
    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) {
759 760 761 762 763 764 765 766 767 768 769 770 771 772
    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;
773
        }
774 775 776 777 778 779 780 781 782
        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;
783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820
    }
    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 {

821 822
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
823
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
824
        os << "* In restart. Distribution read in from .h5 file." << endl;
825 826
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
827
            os << "* Main Distribution" << endl
828
               << "-----------------" << endl;
829 830 831 832 833 834 835 836

        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++) {
837
            os << "* " << endl;
adelmann's avatar
adelmann committed
838 839
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
840 841 842 843
            addedDistributions_m.at(distIndex)->PrintDist(os, particlesPerDist_m.at(distCount));
            distCount++;
        }

844
        os << "* " << endl;
845
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
846
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
847

848
            //            if (numberOfEnergyBins_m > 1)
adelmann's avatar
adelmann committed
849
            //    PrintEnergyBins(os);
850 851 852
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
853 854
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
855 856
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
857
            os << "* " << endl;
858
            PrintEmissionModel(os);
859
            os << "* " << endl;
860
        } else
adelmann's avatar
adelmann committed
861
            os << "* Distribution is injected." << endl;
862
    }
863 864
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
865 866 867 868 869 870 871 872 873 874 875 876 877 878

    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
879

880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 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
    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();
    }
}

void Distribution::ApplyEmissionModel(double eZ, double &px, double &py, double &pz) {

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
        ApplyEmissModelNone(pz);
        break;
    case EmissionModelT::ASTRA:
        ApplyEmissModelAstra(px, py, pz);
        break;
    case EmissionModelT::NONEQUIL:
        ApplyEmissModelNonEquil(eZ, px, py, pz);
        break;
    default:
        break;
    }
}

void Distribution::ApplyEmissModelAstra(double &px, double &py, double &pz) {

949 950
    double phi = 2.0 * acos(sqrt(gsl_rng_uniform(randGenEmit_m)));
    double theta = 2.0 * Physics::pi * gsl_rng_uniform(randGenEmit_m);
951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967

    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
968 969
        - sqrt(Physics::q_e * eZ /
               (4.0 * Physics::pi * Physics::epsilon_0));
970 971 972 973 974 975
    double lowEnergyLimit = cathodeFermiEnergy_m + phiEffective - laserEnergy_m;

    // Generate emission energy.
    double energy = 0.0;
    bool allow = false;
    while (!allow) {
976 977
        energy = lowEnergyLimit + (gsl_rng_uniform(randGenEmit_m)*emitEnergyUpperLimit_m);
        double randFuncValue = gsl_rng_uniform(randGenEmit_m);
978 979 980 981 982
        double funcValue = (1.0
                            - 1.0
                            / (1.0
                               + exp((energy + laserEnergy_m - cathodeFermiEnergy_m)
                                     / (Physics::kB * cathodeTemp_m))))
983 984 985
            / (1.0
               + exp((energy - cathodeFermiEnergy_m)
                     / (Physics::kB * cathodeTemp_m)));
986 987 988 989 990 991 992
        if (randFuncValue <= funcValue)
            allow = true;
    }

    // Compute emission angles.
    double energyInternal = energy + laserEnergy_m;
    double energyExternal = energy + laserEnergy_m
993
        - cathodeFermiEnergy_m - phiEffective;
994 995 996

    double thetaInMax = acos(sqrt((cathodeFermiEnergy_m + phiEffective)
                                  / (energy + laserEnergy_m)));
997
    double thetaIn = gsl_rng_uniform(randGenEmit_m)*thetaInMax;
998
    double sinThetaOut = sin(thetaIn) * sqrt(energyInternal / energyExternal);
999
    double phi = Physics::two_pi * gsl_rng_uniform(randGenEmit_m);
1000 1001 1002

    // Compute emission momenta.
    double betaGammaExternal
1003
        = sqrt(pow(energyExternal / (Physics::m_e * 1.0e9) + 1.0, 2.0) - 1.0);
1004 1005 1006 1007 1008 1009 1010 1011 1012

    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) {

1013
    typedef std::vector<Distribution *>::iterator iterator;
1014

1015
    if (numberOfDistributions_m == 1)
1016 1017
        particlesPerDist_m.push_back(numberOfParticles);
    else {
1018 1019 1020 1021
        double totalWeight = GetWeight();
        for (iterator it = addedDistributions_m.begin(); it != addedDistributions_m.end(); it++) {
            totalWeight += (*it)->GetWeight();
        }
1022 1023 1024

        particlesPerDist_m.push_back(0);
        size_t numberOfCommittedPart = 0;
1025 1026 1027 1028
        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;
1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057
        }

        // 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
1058
    currentEnergyBin_m = 1;
1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081
    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) {

1082 1083
    size_t numberOfDistParticles = tOrZDist_m.size();
    reduce(numberOfDistParticles, numberOfDistParticles, OpAddAssign());
1084

1085
    if (numberOfDistParticles != numberOfParticles) {
1086 1087
        *gmsg << "\n--------------------------------------------------" << endl
              << "Warning!! The number of particles in the initial" << endl
1088
              << "distribution is " << numberOfDistParticles << "." << endl << endl
1089
              << "This is different from the number of particles" << endl
1090
              << "defined by the BEAM command: " << numberOfParticles << endl << endl
1091 1092 1093 1094 1095
              << "This is often happens when using a FROMFILE type" << endl
              << "distribution and not matching the number of" << endl
              << "particles in the particles file(s) with the number" << endl
              << "given in the BEAM command." << endl << endl
              << "The number of particles in the initial distribution" << endl
1096
              << "(" << numberOfDistParticles << ") "
1097 1098 1099
              << "will take precedence." << endl
              << "---------------------------------------------------\n" << endl;
    }
1100
    numberOfParticles = numberOfDistParticles;
1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118
}

void Distribution::ChooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits) {

    /*
     * Toggle what units to use for inputing momentum.
     */
    std::string inputUnits = Attributes::getString(itsAttr[AttributesT::INPUTMOUNITS]);
    if (inputUnits == "NONE")
        inputMoUnits_m = InputMomentumUnitsT::NONE;
    else if (inputUnits == "EV")
        inputMoUnits_m = InputMomentumUnitsT::EV;
    else
        inputMoUnits_m = inputMoUnits;

}

double Distribution::ConvertBetaGammaToeV(double valueInBetaGamma, double massIneV) {
1119 1120 1121 1122
    if (valueInBetaGamma < 0)
        return -1.0 * (sqrt(pow(valueInBetaGamma, 2.0) + 1.0) - 1.0) * massIneV;
    else
        return (sqrt(pow(valueInBetaGamma, 2.0) + 1.0) - 1.0) * massIneV;
1123 1124 1125
}

double Distribution::ConverteVToBetaGamma(double valueIneV, double massIneV) {
1126 1127 1128 1129
    if (valueIneV < 0)
        return -1.0 * sqrt( pow( -1.0 * valueIneV / massIneV + 1.0, 2.0) - 1.0);
    else
        return sqrt( pow( valueIneV / massIneV + 1.0, 2.0) - 1.0);
1130 1131 1132 1133 1134 1135 1136 1137
}

double Distribution::ConvertMeVPerCToBetaGamma(double valueInMeVPerC, double massIneV) {
    return sqrt(pow(valueInMeVPerC * 1.0e6 * Physics::c / massIneV + 1.0, 2.0) - 1.0);
}

void Distribution::CreateDistributionBinomial(size_t numberOfParticles, double massIneV) {

kraus's avatar
kraus committed
1138
    SetDistParametersBinomial(massIneV);
1139 1140 1141 1142 1143 1144 1145 1146 1147