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

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

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

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

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
        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
175
        FMTYPE,                     // field map type used in matched gauss distribution
176 177 178 179 180
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
181
        RGUESS,
182 183
        ID1,                       // special particle that the user can set
        ID2,                       // special particle that the user can set
184
        SIZE
gsell's avatar
gsell committed
185 186 187
    };
}

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

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

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

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

285 286 287 288 289
    setFieldEmissionParameters();

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

Distribution::~Distribution() {

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

388 389 390
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
391 392
    }

393 394 395 396 397
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

398 399 400
    if (randGen_m != NULL) {
        gsl_rng_free(randGen_m);
        randGen_m = NULL;
401 402 403 404 405 406
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
407 408

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

/**
421 422
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
423 424 425
 * @param
 * @param
 * @param
426
 */
427
void Distribution::writeToFile() {
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
    /*
      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();
      }
      }
    */
444 445
}

446 447 448
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
449 450
}

451 452 453
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
454

455 456
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
457

458 459
void Distribution::update() {
}
adelmann's avatar
adelmann committed
460

461
void Distribution::create(size_t &numberOfParticles, double massIneV) {
adelmann's avatar
adelmann committed
462

463
    setFieldEmissionParameters();
adelmann's avatar
adelmann committed
464

465
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
466

467
    case DistrTypeT::MATCHEDGAUSS:
468
        createMatchedGaussDistribution(numberOfParticles, massIneV);
469
        break;
470
    case DistrTypeT::FROMFILE:
471
        createDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
472
        break;
473
    case DistrTypeT::GAUSS:
474
        createDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
475
        break;
476
    case DistrTypeT::BINOMIAL:
477
        createDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
478
        break;
479
    case DistrTypeT::FLATTOP:
480
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
481
        break;
482
    case DistrTypeT::GUNGAUSSFLATTOPTH:
483
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
484
        break;
485
    case DistrTypeT::ASTRAFLATTOPTH:
486
        createDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
487
        break;
488 489 490 491
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
492 493

    if (emitting_m) {
494 495

        std::string model = Util::toUpper(Attributes::getString(itsAttr[AttributesT::EMISSIONMODEL]));
kraus's avatar
kraus committed
496
        unsigned int numAdditionalRNsPerParticle = 0;
497 498 499 500
        if (model == "ASTRA" ||
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {

kraus's avatar
kraus committed
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519
            numAdditionalRNsPerParticle = 2;
        }

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

        int saveProcessor = -1;

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

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

            if (Ippl::myNode() == saveProcessor) {
                std::vector<double> rns;
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
520
                    double x = gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
521 522 523 524 525
                    rns.push_back(x);
                }
                additionalRNs_m.push_back(rns);
            } else {
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
526
                    gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
527 528 529 530 531
                }
            }
        }
    }

532 533 534
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
    Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
535
}
gsell's avatar
gsell committed
536

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

539 540 541
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
542 543
        double vw = this->getVw();
        double vt = this->getvVThermal();
544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567
        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
568

569 570 571 572 573 574 575
                        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
576
                        } else {
577 578
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
579
                        }
580
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
581
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
582 583 584 585 586 587
                        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
588 589
                    }
                }
590 591 592
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
593
            }
594 595 596
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
597
        }
598
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
599

600
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
601

602 603 604 605 606
        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
607

608 609 610 611
        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
612

613 614 615 616 617 618 619 620 621
                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
622
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
623 624 625 626 627 628
                        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
629 630

                    }
631

gsell's avatar
gsell committed
632
                }
633 634 635
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
636

637
            }
638

639 640 641 642 643 644
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
645

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

648
    IpplTimings::startTimer(beam.distrReload_m);
649

650 651
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
652

653 654 655 656
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
657

658 659
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
660

661
    beam.create(numParticles);
662

663 664
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
665 666 667

    beam.boundp();

668 669 670 671
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
672 673
    IpplTimings::stopTimer(beam.distrReload_m);

674 675 676 677
    *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
678 679
}

680
void Distribution::doRestartOpalCycl(PartBunch &beam,
681 682 683 684
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
685

686
    // h5_int64_t rc;
687
    IpplTimings::startTimer(beam.distrReload_m);
688
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
689

690 691
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
692

693 694 695 696
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
697

698 699
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
700

701
    beam.create(numParticles);
702

703 704
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
705

706
    beam.Q = beam.getChargePerParticle();
707

708
    beam.boundp();
709

710
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
711

712 713 714 715 716 717
    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);
718

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

722
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
723

724 725
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
726 727
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
728 729
            energyBins_m = new PartBinsCyc(1000, beam.getNumBunch());
            beam.setPBins(energyBins_m);
730 731
        }

732 733
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
734 735 736 737 738 739 740 741 742 743 744 745 746 747

        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);
748
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
749 750 751 752 753 754 755

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

756
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
757 758 759
    IpplTimings::stopTimer(beam.distrReload_m);
}

760
void Distribution::doRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
761
                                  H5PartWrapper *dataSource) {
762
    IpplTimings::startTimer(beam.distrReload_m);
763 764
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
765 766 767 768 769 770

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

771 772
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
773 774 775 776 777 778 779 780 781 782 783 784 785 786 787

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

788
double Distribution::getTEmission() {
789 790 791 792
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

793
    distT_m = Util::toUpper(Attributes::getString(itsAttr[AttributesT::DISTRIBUTION]));
794 795 796 797 798 799 800 801 802 803
    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]);
804
    cutoff_m = Attributes::getReal(itsAttr[ LegacyAttributesT::CUTOFF]);
805 806 807 808 809 810 811
    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) {
812 813 814 815 816 817 818 819 820 821 822 823 824 825
    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;
826
        }
827 828 829 830 831 832 833 834 835
        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;
836 837 838 839
    }
    return tEmission_m;
}

840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
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;
871 872 873

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

874 875
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
876
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
877
        os << "* In restart. Distribution read in from .h5 file." << endl;
878 879
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
880
            os << "* Main Distribution" << endl
881
               << "-----------------" << endl;
882 883

        if (particlesPerDist_m.empty())
884
            printDist(os, 0);
885
        else
886
            printDist(os, particlesPerDist_m.at(0));
887 888 889

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
890
            os << "* " << endl;
adelmann's avatar
adelmann committed
891 892
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
893
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
894 895 896
            distCount++;
        }

897
        os << "* " << endl;
898
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
899
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
900

901
            //            if (numberOfEnergyBins_m > 1)
902
            //    printEnergyBins(os);
903 904 905
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
906 907
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
908 909
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
910
            os << "* " << endl;
911
            printEmissionModel(os);
912
            os << "* " << endl;
913
        } else
adelmann's avatar
adelmann committed
914
            os << "* Distribution is injected." << endl;
915
    }
916 917
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
918 919 920 921

    return os;
}

922
const PartData &Distribution::getReference() const {
923 924 925 926 927
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

928
void Distribution::addDistributions() {
929 930 931
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
932

933 934 935 936 937
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
         addedDistIt != addedDistributions_m.end(); addedDistIt++) {

        std::vector<double>::iterator distIt;
938 939
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
940 941 942
             distIt++) {
            xDist_m.push_back(*distIt);
        }
943
        (*addedDistIt)->eraseXDist();
944

945 946
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
947 948 949
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
950
        (*addedDistIt)->eraseBGxDist();
951

952 953
        for (distIt = (*addedDistIt)->getYDist().begin();
             distIt != (*addedDistIt)->getYDist().end();
954 955 956
             distIt++) {
            yDist_m.push_back(*distIt);
        }
957
        (*addedDistIt)->eraseYDist();
958

959 960
        for (distIt = (*addedDistIt)->getBGyDist().begin();
             distIt != (*addedDistIt)->getBGyDist().end();
961 962 963
             distIt++) {
            pyDist_m.push_back(*distIt);
        }
964
        (*addedDistIt)->eraseBGyDist();
965

966 967
        for (distIt = (*addedDistIt)->getTOrZDist().begin();
             distIt != (*addedDistIt)->getTOrZDist().end();
968 969 970
             distIt++) {
            tOrZDist_m.push_back(*distIt);
        }
971
        (*addedDistIt)->eraseTOrZDist();
972

973 974
        for (distIt = (*addedDistIt)->getBGzDist().begin();
             distIt != (*addedDistIt)->getBGzDist().end();
975 976 977
             distIt++) {
            pzDist_m.push_back(*distIt);
        }
978
        (*addedDistIt)->eraseBGzDist();
979 980 981
    }
}

982
void Distribution::applyEmissionModel(double eZ, double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
983 984 985 986

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
987
        applyEmissModelNone(pz);
988 989
        break;
    case EmissionModelT::ASTRA:
990
        applyEmissModelAstra(px, py, pz, additionalRNs);
991 992
        break;
    case EmissionModelT::NONEQUIL:
993
        applyEmissModelNonEquil(eZ, px, py, pz);
994 995 996 997 998 999
        break;
    default:
        break;
    }
}

1000
void Distribution::applyEmissModelAstra(double &px, double &py, double &pz, const std::vector<double> &additionalRNs) {
1001

kraus's avatar
kraus committed
1002 1003
    double phi = 2.0 * acos(sqrt(additionalRNs[0]));
    double theta = 2.0 * Physics::pi * additionalRNs[1];
1004 1005