Distribution.cpp 191 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"
29 30 31 32
#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBunch.h"
#include "Algorithms/PartBins.h"
#include "Algorithms/bet/EnvelopeBunch.h"
33
#include "Structure/Beam.h"
34
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
35 36 37
#include "Algorithms/PartBinsCyc.h"
#include "BasicActions/Option.h"
#include "Distribution/LaserProfile.h"
38
#include "Elements/OpalBeamline.h"
39
#include "AbstractObjects/BeamSequence.h"
40 41
#include "Structure/H5PartWrapper.h"
#include "Structure/H5PartWrapperForPC.h"
42
#include "Utilities/Util.h"
43 44 45 46

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

50
#include <sys/time.h>
51

kraus's avatar
kraus committed
52 53
#include <boost/regex.hpp>

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

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

snuverink_j's avatar
snuverink_j committed
58 59 60 61 62 63 64
namespace {
    SymTenzor<double, 6> getUnit6x6() {
        SymTenzor<double, 6> unit6x6;
        for (unsigned int i = 0; i < 6u; ++ i) {
            unit6x6(i,i) = 1.0;
        }
        return unit6x6;
65 66 67
    }
}

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

72 73 74
namespace AttributesT
{
    enum AttributesT {
75
        TYPE,
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 138 139
        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,
140 141 142 143 144
        FLIPX,
        FLIPY,
        ROTATE90,
        ROTATE180,
        ROTATE270,
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
        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
178
        FMTYPE,                     // field map type used in matched gauss distribution
179 180 181 182 183
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
184
        RGUESS,
185 186
        ID1,                       // special particle that the user can set
        ID2,                       // special particle that the user can set
187
        SCALABLE,
188
        SIZE
gsell's avatar
gsell committed
189 190 191
    };
}

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

gsell's avatar
gsell committed
214
Distribution::Distribution():
215 216
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
217
    distrTypeT_m(DistrTypeT::NODIST),
218
    numberOfDistributions_m(1),
219 220 221 222 223 224
    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
225
    currentEnergyBin_m(1),
226
    currentSampleBin_m(0),
227 228 229 230
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
231
    randGen_m(NULL),
232
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
233
    pmean_m(0.0),
234 235 236 237 238
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
239 240
    totalNumberParticles_m(0),
    totalNumberEmittedParticles_m(0),
241
    avrgpz_m(0.0),
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
    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),
    I_m(0.0),
    E_m(0.0),
268 269 270 271
    M_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    bega_m(0.0)
272
{
273
    setAttributes();
274 275 276

    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
        delete defaultDistribution;
    }

284
    setFieldEmissionParameters();
285 286 287

    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
gsell's avatar
gsell committed
288 289 290 291 292 293 294
}
/**
 *
 *
 * @param name
 * @param parent
 */
295
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
296
    Definition(name, parent),
297 298
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
299
    numberOfDistributions_m(parent->numberOfDistributions_m),
300 301 302 303 304 305
    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
306
    tEmission_m(parent->tEmission_m),
307 308 309 310 311 312 313 314
    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),
315
    randGen_m(NULL),
316
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
317
    pmean_m(parent->pmean_m),
318 319 320 321 322
    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),
323 324
    totalNumberParticles_m(parent->totalNumberParticles_m),
    totalNumberEmittedParticles_m(parent->totalNumberEmittedParticles_m),
325 326 327 328 329 330 331 332 333 334 335 336
    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),
337
    avrgpz_m(parent->avrgpz_m),
338 339 340 341 342 343 344 345
    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),
346
    correlationMatrix_m(parent->correlationMatrix_m),
347 348 349 350
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
351 352 353 354 355 356 357 358 359 360 361 362
    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),
363 364 365
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
366 367 368
    I_m(parent->I_m),
    E_m(parent->E_m),
    M_m(parent->M_m),
369 370 371 372
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
373
    cutoff_m(parent->cutoff_m),
374
    bega_m(parent->bega_m)
375
{
376 377
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
gsell's avatar
gsell committed
378 379 380 381 382 383 384
}

Distribution::~Distribution() {

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

385 386 387
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
388 389
    }

390 391 392 393 394
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

395 396 397
    if (randGen_m != NULL) {
        gsl_rng_free(randGen_m);
        randGen_m = NULL;
398 399 400 401 402 403
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
404 405

}
406
/*
407
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
408
  for(int i=0; i<M.size1(); ++i) {
409 410 411 412 413
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
414 415
  }
*/
gsell's avatar
gsell committed
416

417 418
/**
 * Calculate the local number of particles evenly and adjust node 0
419 420 421
 * such that n is matched exactly.
 * @param n total number of particles
 * @return n / #cores
422 423 424 425 426
 * @param
 */

size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {

427
    size_t locNumber = n / Ippl::getNodes();
428

429 430 431 432 433
    // make sure the total number is exact
    size_t remainder  = n % Ippl::getNodes();
    size_t myNode = Ippl::myNode();
    if (myNode < remainder)
        ++ locNumber;
434 435 436 437 438 439

    return locNumber;
}



gsell's avatar
gsell committed
440
/**
441 442
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
443 444 445
 * @param
 * @param
 * @param
446
 */
447
void Distribution::writeToFile() {
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463
    /*
      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();
      }
      }
    */
464 465
}

466 467 468
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
469 470
}

471 472 473
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
474

475 476
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
477

478 479
void Distribution::update() {
}
adelmann's avatar
adelmann committed
480

kraus's avatar
kraus committed
481 482 483 484 485 486 487 488 489 490 491 492
void Distribution::create(size_t &numberOfParticles, double massIneV) {

    /*
     * If Options::cZero is true than we reflect generated distribution
     * to ensure that the transverse averages are 0.0.
     *
     * For now we just cut the number of generated particles in half.
     */
    size_t numberOfLocalParticles = numberOfParticles;
    if (Options::cZero && distrTypeT_m != DistrTypeT::FROMFILE) {
        numberOfLocalParticles = (numberOfParticles + 1) / 2;
    }
adelmann's avatar
adelmann committed
493

494
    size_t mySeed = Options::seed;
adelmann's avatar
adelmann committed
495

496
    if (Options::seed == -1) {
Andreas Adelmann's avatar
Andreas Adelmann committed
497 498
        struct timeval tv;
        gettimeofday(&tv,0);
499
        mySeed = tv.tv_sec + tv.tv_usec;
Andreas Adelmann's avatar
Andreas Adelmann committed
500
    }
501 502

    if (Attributes::getBool(itsAttr[AttributesT::SCALABLE])) {
kraus's avatar
kraus committed
503
        numberOfLocalParticles = getNumOfLocalParticlesToCreate(numberOfLocalParticles);
504 505 506 507 508 509
        *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << " + core_id\n"
              << "* is scalable with number of particles and cores." << endl;
        mySeed += Ippl::myNode();
    } else {
        *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << "\n"
              << "* isn't scalable with number of particles and cores." << endl;
Andreas Adelmann's avatar
Andreas Adelmann committed
510 511
    }

512 513 514
    gsl_rng_set(randGen_m, mySeed);

    setFieldEmissionParameters();
Andreas Adelmann's avatar
Andreas Adelmann committed
515

516
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
517

518
    case DistrTypeT::MATCHEDGAUSS:
kraus's avatar
kraus committed
519
        createMatchedGaussDistribution(numberOfLocalParticles, massIneV);
520
        break;
521
    case DistrTypeT::FROMFILE:
kraus's avatar
kraus committed
522
        createDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
523
        break;
524
    case DistrTypeT::GAUSS:
kraus's avatar
kraus committed
525
        createDistributionGauss(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
526
        break;
527
    case DistrTypeT::BINOMIAL:
kraus's avatar
kraus committed
528
        createDistributionBinomial(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
529
        break;
530
    case DistrTypeT::FLATTOP:
kraus's avatar
kraus committed
531
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
532
        break;
533
    case DistrTypeT::GUNGAUSSFLATTOPTH:
kraus's avatar
kraus committed
534
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
535
        break;
536
    case DistrTypeT::ASTRAFLATTOPTH:
kraus's avatar
kraus committed
537
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
538
        break;
539 540 541 542
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
543 544

    if (emitting_m) {
545

kraus's avatar
kraus committed
546
        unsigned int numAdditionalRNsPerParticle = 0;
547
        if (emissionModel_m == EmissionModelT::ASTRA ||
548 549
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
550

kraus's avatar
kraus committed
551
            numAdditionalRNsPerParticle = 2;
552
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
553 554 555 556 557
            if (Options::cZero) {
                numAdditionalRNsPerParticle = 40;
            } else {
                numAdditionalRNsPerParticle = 20;
            }
kraus's avatar
kraus committed
558 559
        }

560 561 562 563 564
        int saveProcessor = -1;
        const int myNode = Ippl::myNode();
        const int numNodes = Ippl::getNodes();
        const bool scalable = Attributes::getBool(itsAttr[AttributesT::SCALABLE]);

kraus's avatar
kraus committed
565
        for (size_t partIndex = 0; partIndex < numberOfLocalParticles; ++ partIndex) {
566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582

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

            if (scalable || myNode == saveProcessor) {
                std::vector<double> rns;
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
                    double x = gsl_rng_uniform(randGen_m);
                    rns.push_back(x);
                }
                additionalRNs_m.push_back(rns);
            } else {
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
                    gsl_rng_uniform(randGen_m);
                }
kraus's avatar
kraus committed
583 584 585 586
            }
        }
    }

587 588
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
589

kraus's avatar
kraus committed
590 591 592
    // Now reflect particles if Options::cZero is true
    reflectDistribution(numberOfLocalParticles);

593
    if (Options::seed != -1)
Andreas Adelmann's avatar
Andreas Adelmann committed
594
        Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
kraus's avatar
kraus committed
595 596

    particlesPerDist_m[0] = tOrZDist_m.size();
597
}
gsell's avatar
gsell committed
598

599
void  Distribution::createPriPart(PartBunch *beam, BoundaryGeometry &bg) {
gsell's avatar
gsell committed
600

601 602 603
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
604 605
        double vw = this->getVw();
        double vt = this->getvVThermal();
606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
        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
630

631 632 633 634 635 636 637
                        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
638
                        } else {
639 640
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
641
                        }
642
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
643
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
644 645 646 647 648 649
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        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
650 651
                    }
                }
652 653 654
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
655
            }
656 657 658
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
659
        }
660
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
661

662
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
663

664 665 666 667 668
        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
669

670 671 672 673
        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
674

675 676 677 678 679 680 681 682 683
                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
684
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
685 686 687 688 689 690
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        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
691 692

                    }
693

gsell's avatar
gsell committed
694
                }
695 696 697
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
698

699
            }
700

701 702 703 704 705 706
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
707

708
void Distribution::doRestartOpalT(PartBunch &beam, size_t Np, int restartStep, H5PartWrapper *dataSource) {
709

710
    IpplTimings::startTimer(beam.distrReload_m);
711

712 713
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
714

715 716 717 718
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
719

720 721
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
722

723
    beam.create(numParticles);
724

725 726
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
727 728 729

    beam.boundp();

730 731 732 733
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
734 735
    IpplTimings::stopTimer(beam.distrReload_m);

736 737 738 739
    *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
740 741
}

742
void Distribution::doRestartOpalCycl(PartBunch &beam,
743 744 745 746
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
747

748
    // h5_int64_t rc;
749
    IpplTimings::startTimer(beam.distrReload_m);
750
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
751

752 753
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
754

755 756 757 758
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
759

760 761
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
762

763
    beam.create(numParticles);
764

765 766
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
767

768
    beam.Q = beam.getChargePerParticle();
769

770
    beam.boundp();
771

772
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
773

774 775 776 777 778 779
    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);
780

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

784
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
785

786 787
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
788 789
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
790 791
            energyBins_m = new PartBinsCyc(1000, beam.getNumBunch());
            beam.setPBins(energyBins_m);
792 793
        }

794 795
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
796 797 798 799 800 801 802 803 804 805 806 807 808 809

        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);
810
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
811 812 813 814 815 816 817

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

818
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
819 820 821
    IpplTimings::stopTimer(beam.distrReload_m);
}

822
void Distribution::doRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
823
                                  H5PartWrapper *dataSource) {
824
    IpplTimings::startTimer(beam.distrReload_m);
825 826
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
827 828 829 830 831 832

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

833 834
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
835 836 837 838 839 840 841 842 843 844 845 846 847 848 849

    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;
}

850
double Distribution::getTEmission() {
851 852 853 854
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

855
    setDistType();
856 857

    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[AttributesT::TPULSEFWHM]);
858
    cutoff_m = Attributes::getReal(itsAttr[LegacyAttributesT::CUTOFF]);
859 860 861 862 863 864 865
    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) {
866 867 868 869 870 871 872 873 874 875 876 877 878 879
    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;
880
        }
881 882 883 884 885 886 887 888 889
        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;
890 891 892 893
    }
    return tEmission_m;
}

894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912
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]);}

913
std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::TYPE]);}
914 915 916 917 918 919 920 921 922 923 924

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;
925 926 927

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

928 929
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
930
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
931
        os << "* In restart. Distribution read in from .h5 file." << endl;
932 933
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
934
            os << "* Main Distribution" << endl
935
               << "-----------------" << endl;
936 937

        if (particlesPerDist_m.empty())
938
            printDist(os, 0);
939
        else
940
            printDist(os, particlesPerDist_m.at(0));
941 942 943

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
944
            os << "* " << endl;
adelmann's avatar
adelmann committed
945 946
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
947
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
948 949 950
            distCount++;
        }

951
        os << "* " << endl;
952
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
953
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
954

955
            //            if (numberOfEnergyBins_m > 1)
956
            //    printEnergyBins(os);
957 958 959
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
960 961
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
962 963
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
964
            os << "* " << endl;
965
            printEmissionModel(os);
966
            os << "* " << endl;
967
        } else
adelmann's avatar
adelmann committed
968
            os << "* Distribution is injected." << endl;
969
    }
970 971
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
972 973 974 975

    return os;
}

976
const PartData &Distribution::getReference() const {
977 978 979 980 981
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

982
void Distribution::addDistributions() {
983 984 985
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
986

kraus's avatar
kraus committed
987
    size_t idx = 1;
988 989
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
kraus's avatar
kraus committed
990 991 992
         addedDistIt != addedDistributions_m.end(); ++ addedDistIt, ++ idx) {

        particlesPerDist_m[idx] = (*addedDistIt)->tOrZDist_m.size();
993 994

        std::vector<double>::iterator distIt;
995 996
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();