Distribution.cpp 182 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
        SIZE
gsell's avatar
gsell committed
186 187 188
    };
}

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

gsell's avatar
gsell committed
211
Distribution::Distribution():
212 213
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
214
    distrTypeT_m(DistrTypeT::NODIST),
215
    numberOfDistributions_m(1),
216 217 218 219 220 221
    emitting_m(false),
    scan_m(false),
    emissionModel_m(EmissionModelT::NONE),
    tEmission_m(0.0),
    tBin_m(0.0),
    currentEmissionTime_m(0.0),
kraus's avatar
kraus committed
222
    currentEnergyBin_m(1),
223
    currentSampleBin_m(0),
224 225 226 227
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
228
    randGen_m(NULL),
229
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
230
    pmean_m(0.0),
231 232 233 234 235
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
236 237
    totalNumberParticles_m(0),
    totalNumberEmittedParticles_m(0),
238
    avrgpz_m(0.0),
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
    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),
265 266 267 268
    M_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    bega_m(0.0)
269
{
270
    setAttributes();
271 272 273

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

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

281 282 283 284 285
    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
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 376
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
    gsl_rng_set(randGen_m, Options::seed);
gsell's avatar
gsell committed
377 378 379 380 381 382 383
}

Distribution::~Distribution() {

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

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

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

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

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

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

416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
/**
 * Calculate the local number of particles evenly and adjust node 0
 * such that n is matched exactly. 
 * @param n total number of particles 
 * @return n / #cores 
 * @param
 */

size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {
    
    size_t locNumber = n / Ippl::getNodes();

    // make sure the total number is exact                                                                               
    size_t reminder  = n % Ippl::getNodes();

    if (Ippl::myNode() == 0)
        locNumber += reminder;   

    return locNumber;
}



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

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

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

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

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

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

482
    setFieldEmissionParameters();
adelmann's avatar
adelmann committed
483

484 485
    size_t locNumber = getNumOfLocalParticlesToCreate(numberOfParticles);

486
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
487

488
    case DistrTypeT::MATCHEDGAUSS:
489
        createMatchedGaussDistribution(locNumber, massIneV);
490
        break;
491
    case DistrTypeT::FROMFILE:
492
        createDistributionFromFile(locNumber, massIneV);
gsell's avatar
gsell committed
493
        break;
494
    case DistrTypeT::GAUSS:
495
        createDistributionGauss(locNumber, massIneV);
gsell's avatar
gsell committed
496
        break;
497
    case DistrTypeT::BINOMIAL:
498
        createDistributionBinomial(locNumber, massIneV);
gsell's avatar
gsell committed
499
        break;
500
    case DistrTypeT::FLATTOP:
501
        createDistributionFlattop(locNumber, massIneV);
gsell's avatar
gsell committed
502
        break;
503
    case DistrTypeT::GUNGAUSSFLATTOPTH:
504
        createDistributionFlattop(locNumber, massIneV);
gsell's avatar
gsell committed
505
        break;
506
    case DistrTypeT::ASTRAFLATTOPTH:
507
        createDistributionFlattop(locNumber, massIneV);
gsell's avatar
gsell committed
508
        break;
509 510 511 512
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
513 514

    if (emitting_m) {
515

kraus's avatar
kraus committed
516
        unsigned int numAdditionalRNsPerParticle = 0;
517
        if (emissionModel_m == EmissionModelT::ASTRA ||
518 519 520
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {

kraus's avatar
kraus committed
521
            numAdditionalRNsPerParticle = 2;
522 523
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
            numAdditionalRNsPerParticle = 20;
kraus's avatar
kraus committed
524 525 526 527 528
        }

        if (Options::cZero) {
            numAdditionalRNsPerParticle *= 2;
        }
529 530 531 532 533 534
        
        for (size_t partIndex = 0; partIndex < locNumber; ++ partIndex) {
            std::vector<double> rns;
            for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
                double x = gsl_rng_uniform(randGen_m);
                rns.push_back(x);
kraus's avatar
kraus committed
535
            }
536
            additionalRNs_m.push_back(rns);
kraus's avatar
kraus committed
537 538 539
        }
    }

540 541
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
542
        
543
    Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
544
}
gsell's avatar
gsell committed
545

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

548 549 550
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
551 552
        double vw = this->getVw();
        double vt = this->getvVThermal();
553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576
        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
577

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

609
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
610

611 612 613 614 615
        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
616

617 618 619 620
        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
621

622 623 624 625 626 627 628 629 630
                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
631
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
632 633 634 635 636 637
                        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
638 639

                    }
640

gsell's avatar
gsell committed
641
                }
642 643 644
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
645

646
            }
647

648 649 650 651 652 653
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
654

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

657
    IpplTimings::startTimer(beam.distrReload_m);
658

659 660
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
661

662 663 664 665
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
666

667 668
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
669

670
    beam.create(numParticles);
671

672 673
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
674 675 676

    beam.boundp();

677 678 679 680
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
681 682
    IpplTimings::stopTimer(beam.distrReload_m);

683 684 685 686
    *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
687 688
}

689
void Distribution::doRestartOpalCycl(PartBunch &beam,
690 691 692 693
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
694

695
    // h5_int64_t rc;
696
    IpplTimings::startTimer(beam.distrReload_m);
697
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
698

699 700
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
701

702 703 704 705
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
706

707 708
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
709

710
    beam.create(numParticles);
711

712 713
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
714

715
    beam.Q = beam.getChargePerParticle();
716

717
    beam.boundp();
718

719
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
720

721 722 723 724 725 726
    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);
727

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

731
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
732

733 734
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
735 736
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
737 738
            energyBins_m = new PartBinsCyc(1000, beam.getNumBunch());
            beam.setPBins(energyBins_m);
739 740
        }

741 742
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
743 744 745 746 747 748 749 750 751 752 753 754 755 756

        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);
757
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
758 759 760 761 762 763 764

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

765
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
766 767 768
    IpplTimings::stopTimer(beam.distrReload_m);
}

769
void Distribution::doRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
770
                                  H5PartWrapper *dataSource) {
771
    IpplTimings::startTimer(beam.distrReload_m);
772 773
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
774 775 776 777 778 779

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

780 781
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
782 783 784 785 786 787 788 789 790 791 792 793 794 795 796

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

797
double Distribution::getTEmission() {
798 799 800 801
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

802
    setDistType();
803 804

    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[AttributesT::TPULSEFWHM]);
805
    cutoff_m = Attributes::getReal(itsAttr[LegacyAttributesT::CUTOFF]);
806 807 808 809 810 811 812
    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) {
813 814 815 816 817 818 819 820 821 822 823 824 825 826
    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;
827
        }
828 829 830 831 832 833 834 835 836
        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;
837 838 839 840
    }
    return tEmission_m;
}

841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859
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]);}

860
std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::TYPE]);}
861 862 863 864 865 866 867 868 869 870 871

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;
872 873 874

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

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

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

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

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

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

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

    return os;
}

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

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

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

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

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

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

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

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

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

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

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
988
        applyEmissModelNone(pz);