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 54
#define SMALLESTCUTOFF 1e-12

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

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

69 70 71
namespace AttributesT
{
    enum AttributesT {
72
        TYPE,
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
        SBIN,
        SIGMAPT,
        CUTOFF,
        T,
        PT,
198 199 200 201 202 203 204 205
        // ALPHAX,
        // ALPHAY,
        // BETAX,
        // BETAY,
        // DX,
        // DDX,
        // DY,
        // DDY,
206 207
        SIZE
    };
208 209
}

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

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

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

280 281 282 283 284
    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
285 286 287 288 289 290 291
}
/**
 *
 *
 * @param name
 * @param parent
 */
292
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
293
    Definition(name, parent),
294 295
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
296
    numberOfDistributions_m(parent->numberOfDistributions_m),
297 298 299 300 301 302
    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
303
    tEmission_m(parent->tEmission_m),
304 305 306 307 308 309 310 311
    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),
312
    randGen_m(NULL),
313
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
314
    pmean_m(parent->pmean_m),
315 316 317 318 319
    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),
320 321
    totalNumberParticles_m(parent->totalNumberParticles_m),
    totalNumberEmittedParticles_m(parent->totalNumberEmittedParticles_m),
322 323 324 325 326 327 328 329 330 331 332 333
    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),
334
    avrgpz_m(parent->avrgpz_m),
335 336 337 338 339 340 341 342
    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),
343
    correlationMatrix_m(parent->correlationMatrix_m),
344 345 346 347
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
348 349 350 351 352 353 354 355 356 357 358 359
    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),
360 361 362
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
363 364 365
    I_m(parent->I_m),
    E_m(parent->E_m),
    M_m(parent->M_m),
366 367 368 369
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
370
    cutoff_m(parent->cutoff_m),
371
    bega_m(parent->bega_m)
372
{
373 374 375
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
    gsl_rng_set(randGen_m, Options::seed);
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 417
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
418 419 420
 * @param
 * @param
 * @param
421
 */
422
void Distribution::writeToFile() {
423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
    /*
      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();
      }
      }
    */
439 440
}

441 442 443
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
444 445
}

446 447 448
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
449

450 451
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
452

453 454
void Distribution::update() {
}
adelmann's avatar
adelmann committed
455

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

458
    setFieldEmissionParameters();
adelmann's avatar
adelmann committed
459

460
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
461

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

    if (emitting_m) {
489

kraus's avatar
kraus committed
490
        unsigned int numAdditionalRNsPerParticle = 0;
491
        if (emissionModel_m == EmissionModelT::ASTRA ||
492 493 494
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {

kraus's avatar
kraus committed
495
            numAdditionalRNsPerParticle = 2;
496 497
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
            numAdditionalRNsPerParticle = 20;
kraus's avatar
kraus committed
498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515
        }

        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) {
516
                    double x = gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
517 518 519 520 521
                    rns.push_back(x);
                }
                additionalRNs_m.push_back(rns);
            } else {
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
522
                    gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
523 524 525 526 527
                }
            }
        }
    }

528 529 530
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
    Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
531
}
gsell's avatar
gsell committed
532

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

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

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

596
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
597

598 599 600 601 602
        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
603

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

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

                    }
627

gsell's avatar
gsell committed
628
                }
629 630 631
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
632

633
            }
634

635 636 637 638 639 640
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
641

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

644
    IpplTimings::startTimer(beam.distrReload_m);
645

646 647
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
648

649 650 651 652
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
653

654 655
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
656

657
    beam.create(numParticles);
658

659 660
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
661 662 663

    beam.boundp();

664 665 666 667
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
668 669
    IpplTimings::stopTimer(beam.distrReload_m);

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

676
void Distribution::doRestartOpalCycl(PartBunch &beam,
677 678 679 680
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
681

682
    // h5_int64_t rc;
683
    IpplTimings::startTimer(beam.distrReload_m);
684
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
685

686 687
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
688

689 690 691 692
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
693

694 695
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
696

697
    beam.create(numParticles);
698

699 700
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
701

702
    beam.Q = beam.getChargePerParticle();
703

704
    beam.boundp();
705

706
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
707

708 709 710 711 712 713
    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);
714

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

718
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
719

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

728 729
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
730 731 732 733 734 735 736 737 738 739 740 741 742 743

        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);
744
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
745 746 747 748 749 750 751

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

752
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
753 754 755
    IpplTimings::stopTimer(beam.distrReload_m);
}

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

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

767 768
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
769 770 771 772 773 774 775 776 777 778 779 780 781 782 783

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

784
double Distribution::getTEmission() {
785 786 787 788
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

789
    setDistType();
790 791

    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[AttributesT::TPULSEFWHM]);
792
    cutoff_m = Attributes::getReal(itsAttr[LegacyAttributesT::CUTOFF]);
793 794 795 796 797 798 799
    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) {
800 801 802 803 804 805 806 807 808 809 810 811 812 813
    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;
814
        }
815 816 817 818 819 820 821 822 823
        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;
824 825 826 827
    }
    return tEmission_m;
}

828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846
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]);}

847
std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::TYPE]);}
848 849 850 851 852 853 854 855 856 857 858

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;
859 860 861

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

862 863
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
864
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
865
        os << "* In restart. Distribution read in from .h5 file." << endl;
866 867
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
868
            os << "* Main Distribution" << endl
869
               << "-----------------" << endl;
870 871

        if (particlesPerDist_m.empty())
872
            printDist(os, 0);
873
        else
874
            printDist(os, particlesPerDist_m.at(0));
875 876 877

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
878
            os << "* " << endl;
adelmann's avatar
adelmann committed
879 880
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
881
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
882 883 884
            distCount++;
        }

885
        os << "* " << endl;
886
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
887
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
888

889
            //            if (numberOfEnergyBins_m > 1)
890
            //    printEnergyBins(os);
891 892 893
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
894 895
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
896 897
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
898
            os << "* " << endl;
899
            printEmissionModel(os);
900
            os << "* " << endl;
901
        } else
adelmann's avatar
adelmann committed
902
            os << "* Distribution is injected." << endl;
903
    }
904 905
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
906 907 908 909

    return os;
}

910
const PartData &Distribution::getReference() const {
911 912 913 914 915
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

916
void Distribution::addDistributions() {
917 918 919
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
920

921 922 923 924 925
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
         addedDistIt != addedDistributions_m.end(); addedDistIt++) {

        std::vector<double>::iterator distIt;
926 927
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
928 929 930
             distIt++) {
            xDist_m.push_back(*distIt);
        }
931
        (*addedDistIt)->eraseXDist();
932

933 934
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
935 936 937
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
938
        (*addedDistIt)->eraseBGxDist();
939

940 941
        for (distIt = (*addedDistIt)->getYDist().begin();
             distIt != (*addedDistIt)->getYDist().end();
942 943 944
             distIt++) {
            yDist_m.push_back(*distIt);
        }
945
        (*addedDistIt)->eraseYDist();
946

947 948
        for (distIt = (*addedDistIt)->getBGyDist().begin();
             distIt != (*addedDistIt)->getBGyDist().end();
949 950 951
             distIt++) {
            pyDist_m.push_back(*distIt);
        }
952
        (*addedDistIt)->eraseBGyDist();
953

954 955
        for (distIt = (*addedDistIt)->getTOrZDist().begin();
             distIt != (*addedDistIt)->getTOrZDist().end();
956 957 958
             distIt++) {
            tOrZDist_m.push_back(*distIt);
        }
959
        (*addedDistIt)->eraseTOrZDist();
960

961 962
        for (distIt = (*addedDistIt)->getBGzDist().begin();
             distIt != (*addedDistIt)->getBGzDist().end();
963 964 965
             distIt++) {
            pzDist_m.push_back(*distIt);
        }
966
        (*addedDistIt)->eraseBGzDist();
967 968 969
    }
}

970
void Distribution::applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs) {
971 972 973 974

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
975
        applyEmissModelNone(pz);
976 977
        break;
    case EmissionModelT::ASTRA:
978
        applyEmissModelAstra(px, py, pz, additionalRNs);
979 980
        break;
    case EmissionModelT::NONEQUIL:
981
        applyEmissModelNonEquil(lowEnergyLimit, px, py, pz, additionalRNs);
982 983 984 985 986 987
        break;
    default:
        break;
    }
}

kraus's avatar