Distribution.cpp 188 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

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

kraus's avatar
kraus committed
54 55
#define SMALLESTCUTOFF 1e-12

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

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

70 71 72
namespace AttributesT
{
    enum AttributesT {
73
        TYPE,
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,
182
        RGUESS,
183 184
        ID1,                       // special particle that the user can set
        ID2,                       // special particle that the user can set
185
        SCALABLE,
186
        SIZE
gsell's avatar
gsell committed
187 188 189
    };
}

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

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

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

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

282
    setFieldEmissionParameters();
283 284 285

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

Distribution::~Distribution() {

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

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

388 389 390 391 392
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

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

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
402 403

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

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

size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {

425
    size_t locNumber = n / Ippl::getNodes();
426

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

    return locNumber;
}



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

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

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

473 474
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
475

476 477
void Distribution::update() {
}
adelmann's avatar
adelmann committed
478

479
void Distribution::create(size_t &totalNumberOfParticles, double massIneV) {
adelmann's avatar
adelmann committed
480

481
    size_t mySeed = Options::seed;
adelmann's avatar
adelmann committed
482

483
    if (Options::seed == -1) {
Andreas Adelmann's avatar
Andreas Adelmann committed
484 485
        struct timeval tv;
        gettimeofday(&tv,0);
486
        mySeed = tv.tv_sec + tv.tv_usec;
Andreas Adelmann's avatar
Andreas Adelmann committed
487
    }
488 489 490 491 492 493 494 495 496 497

    size_t numberOfParticles = totalNumberOfParticles;
    if (Attributes::getBool(itsAttr[AttributesT::SCALABLE])) {
        numberOfParticles = getNumOfLocalParticlesToCreate(totalNumberOfParticles);
        *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
498 499
    }

500 501 502
    gsl_rng_set(randGen_m, mySeed);

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

504
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
505

506
    case DistrTypeT::MATCHEDGAUSS:
507
        createMatchedGaussDistribution(numberOfParticles, massIneV);
508
        break;
509
    case DistrTypeT::FROMFILE:
510
        createDistributionFromFile(totalNumberOfParticles, massIneV);
gsell's avatar
gsell committed
511
        break;
512
    case DistrTypeT::GAUSS:
513
        createDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
514
        break;
515
    case DistrTypeT::BINOMIAL:
516
        createDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
517
        break;
518
    case DistrTypeT::FLATTOP:
519
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
520
        break;
521
    case DistrTypeT::GUNGAUSSFLATTOPTH:
522
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
523
        break;
524
    case DistrTypeT::ASTRAFLATTOPTH:
525
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
526
        break;
527 528 529 530
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
531 532

    if (emitting_m) {
533

kraus's avatar
kraus committed
534
        unsigned int numAdditionalRNsPerParticle = 0;
535
        if (emissionModel_m == EmissionModelT::ASTRA ||
536 537
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
538

kraus's avatar
kraus committed
539
            numAdditionalRNsPerParticle = 2;
540
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
541 542 543 544 545
            if (Options::cZero) {
                numAdditionalRNsPerParticle = 40;
            } else {
                numAdditionalRNsPerParticle = 20;
            }
kraus's avatar
kraus committed
546 547
        }

548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
        int saveProcessor = -1;
        const int myNode = Ippl::myNode();
        const int numNodes = Ippl::getNodes();
        const bool scalable = Attributes::getBool(itsAttr[AttributesT::SCALABLE]);

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

            // 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
571 572 573 574
            }
        }
    }

575 576
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
577 578

    if (Options::seed != -1)
Andreas Adelmann's avatar
Andreas Adelmann committed
579
        Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
580
}
gsell's avatar
gsell committed
581

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

584 585 586
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
587 588
        double vw = this->getVw();
        double vt = this->getvVThermal();
589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
        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
613

614 615 616 617 618 619 620
                        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
621
                        } else {
622 623
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
624
                        }
625
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
626
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
627 628 629 630 631 632
                        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
633 634
                    }
                }
635 636 637
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
638
            }
639 640 641
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
642
        }
643
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
644

645
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
646

647 648 649 650 651
        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
652

653 654 655 656
        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
657

658 659 660 661 662 663 664 665 666
                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
667
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
668 669 670 671 672 673
                        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
674 675

                    }
676

gsell's avatar
gsell committed
677
                }
678 679 680
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
681

682
            }
683

684 685 686 687 688 689
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
690

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

693
    IpplTimings::startTimer(beam.distrReload_m);
694

695 696
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
697

698 699 700 701
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
702

703 704
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
705

706
    beam.create(numParticles);
707

708 709
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
710 711 712

    beam.boundp();

713 714 715 716
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
717 718
    IpplTimings::stopTimer(beam.distrReload_m);

719 720 721 722
    *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
723 724
}

725
void Distribution::doRestartOpalCycl(PartBunch &beam,
726 727 728 729
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
730

731
    // h5_int64_t rc;
732
    IpplTimings::startTimer(beam.distrReload_m);
733
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
734

735 736
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
737

738 739 740 741
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
742

743 744
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
745

746
    beam.create(numParticles);
747

748 749
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
750

751
    beam.Q = beam.getChargePerParticle();
752

753
    beam.boundp();
754

755
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
756

757 758 759 760 761 762
    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);
763

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

767
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
768

769 770
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
771 772
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
773 774
            energyBins_m = new PartBinsCyc(1000, beam.getNumBunch());
            beam.setPBins(energyBins_m);
775 776
        }

777 778
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
779 780 781 782 783 784 785 786 787 788 789 790 791 792

        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);
793
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
794 795 796 797 798 799 800

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

801
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
802 803 804
    IpplTimings::stopTimer(beam.distrReload_m);
}

805
void Distribution::doRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
806
                                  H5PartWrapper *dataSource) {
807
    IpplTimings::startTimer(beam.distrReload_m);
808 809
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
810 811 812 813 814 815

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

816 817
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
818 819 820 821 822 823 824 825 826 827 828 829 830 831 832

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

833
double Distribution::getTEmission() {
834 835 836 837
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

838
    setDistType();
839 840

    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[AttributesT::TPULSEFWHM]);
841
    cutoff_m = Attributes::getReal(itsAttr[LegacyAttributesT::CUTOFF]);
842 843 844 845 846 847 848
    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) {
849 850 851 852 853 854 855 856 857 858 859 860 861 862
    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;
863
        }
864 865 866 867 868 869 870 871 872
        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;
873 874 875 876
    }
    return tEmission_m;
}

877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895
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]);}

896
std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::TYPE]);}
897 898 899 900 901 902 903 904 905 906 907

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;
908 909 910

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

911 912
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
913
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
914
        os << "* In restart. Distribution read in from .h5 file." << endl;
915 916
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
917
            os << "* Main Distribution" << endl
918
               << "-----------------" << endl;
919 920

        if (particlesPerDist_m.empty())
921
            printDist(os, 0);
922
        else
923
            printDist(os, particlesPerDist_m.at(0));
924 925 926

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
927
            os << "* " << endl;
adelmann's avatar
adelmann committed
928 929
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
930
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
931 932 933
            distCount++;
        }

934
        os << "* " << endl;
935
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
936
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
937

938
            //            if (numberOfEnergyBins_m > 1)
939
            //    printEnergyBins(os);
940 941 942
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
943 944
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
945 946
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
947
            os << "* " << endl;
948
            printEmissionModel(os);
949
            os << "* " << endl;
950
        } else
adelmann's avatar
adelmann committed
951
            os << "* Distribution is injected." << endl;
952
    }
953 954
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
955 956 957 958

    return os;
}

959
const PartData &Distribution::getReference() const {
960 961 962 963 964
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

965
void Distribution::addDistributions() {
966 967 968
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
969

970 971 972 973 974
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
         addedDistIt != addedDistributions_m.end(); addedDistIt++) {

        std::vector<double>::iterator distIt;
975 976
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
977 978 979
             distIt++) {
            xDist_m.push_back(*distIt);
        }
980
        (*addedDistIt)->eraseXDist();
981

982 983
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
984 985 986
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
987
        (*addedDistIt)->eraseBGxDist();
988

989 990
        for (distIt = (*addedDistIt)->getYDist().begin();
             distIt != (*addedDistIt)->getYDist().end();
991 992 993
             distIt++) {
            yDist_m.push_back(*distIt);
        }
994
        (*addedDistIt)->eraseYDist();