Distribution.cpp 191 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12
// ------------------------------------------------------------------------
// $RCSfile: Distribution.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.3.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Distribution
//   The class for the OPAL Distribution command.
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
13 14

#include "Distribution/Distribution.h"
15
#include "Distribution/SigmaGenerator.h"
16
#include "AbsBeamline/SpecificElementVisitor.h"
17

18 19 20 21 22 23
#include <cmath>
#include <cfloat>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
24
#include <numeric>
25

26 27
#include "AbstractObjects/Expressions.h"
#include "Attributes/Attributes.h"
28
#include "Utilities/Options.h"
29 30 31 32
#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBunch.h"
#include "Algorithms/PartBins.h"
#include "Algorithms/bet/EnvelopeBunch.h"
33
#include "Structure/Beam.h"
34
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
35 36 37
#include "Algorithms/PartBinsCyc.h"
#include "BasicActions/Option.h"
#include "Distribution/LaserProfile.h"
38
#include "Elements/OpalBeamline.h"
39
#include "AbstractObjects/BeamSequence.h"
40 41
#include "Structure/H5PartWrapper.h"
#include "Structure/H5PartWrapperForPC.h"
42
#include "Utilities/Util.h"
43 44 45 46

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

50
#include <sys/time.h>
51

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

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

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

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

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

Distribution::Distribution():
73
    Definition( Attrib::Legacy::Distribution::SIZE, "DISTRIBUTION",
74
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
75
    distrTypeT_m(DistrTypeT::NODIST),
76
    numberOfDistributions_m(1),
77 78 79 80 81 82
    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
83
    currentEnergyBin_m(1),
84
    currentSampleBin_m(0),
85 86 87 88
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
89
    randGen_m(NULL),
90
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
91
    pmean_m(0.0),
92 93 94 95 96
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
97 98
    totalNumberParticles_m(0),
    totalNumberEmittedParticles_m(0),
99
    avrgpz_m(0.0),
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
    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),
126 127 128 129
    M_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    bega_m(0.0)
130
{
131
    setAttributes();
132 133 134

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

    try {
137
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
138
    } catch(...) {
139 140 141
        delete defaultDistribution;
    }

142
    // setFieldEmissionParameters();
143 144 145

    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
gsell's avatar
gsell committed
146 147 148 149 150 151 152
}
/**
 *
 *
 * @param name
 * @param parent
 */
153
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
154
    Definition(name, parent),
155 156
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
157
    numberOfDistributions_m(parent->numberOfDistributions_m),
158 159 160 161 162 163
    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
164
    tEmission_m(parent->tEmission_m),
165 166 167 168 169 170 171 172
    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),
173
    randGen_m(NULL),
174
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
175
    pmean_m(parent->pmean_m),
176 177 178 179 180
    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),
181 182
    totalNumberParticles_m(parent->totalNumberParticles_m),
    totalNumberEmittedParticles_m(parent->totalNumberEmittedParticles_m),
183 184 185 186 187 188 189 190 191 192 193 194
    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),
195
    avrgpz_m(parent->avrgpz_m),
196 197 198 199 200 201 202 203
    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),
204
    correlationMatrix_m(parent->correlationMatrix_m),
205 206 207 208
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
209 210 211 212 213 214 215 216 217 218 219 220
    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),
221 222 223
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
224 225 226
    I_m(parent->I_m),
    E_m(parent->E_m),
    M_m(parent->M_m),
227 228 229 230
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
231
    cutoff_m(parent->cutoff_m),
232
    bega_m(parent->bega_m)
233
{
234 235
    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
gsell's avatar
gsell committed
236 237 238 239 240 241 242
}

Distribution::~Distribution() {

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

243 244 245
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
246 247
    }

248 249 250 251 252
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

253 254 255
    if (randGen_m != NULL) {
        gsl_rng_free(randGen_m);
        randGen_m = NULL;
256 257 258 259 260 261
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
262 263

}
264
/*
265
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
266
  for(int i=0; i<M.size1(); ++i) {
267 268 269 270 271
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
272 273
  }
*/
gsell's avatar
gsell committed
274

275 276
/**
 * Calculate the local number of particles evenly and adjust node 0
277 278 279
 * such that n is matched exactly.
 * @param n total number of particles
 * @return n / #cores
280 281 282 283 284
 * @param
 */

size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {

285
    size_t locNumber = n / Ippl::getNodes();
286

287 288 289 290 291
    // make sure the total number is exact
    size_t remainder  = n % Ippl::getNodes();
    size_t myNode = Ippl::myNode();
    if (myNode < remainder)
        ++ locNumber;
292 293 294 295 296 297

    return locNumber;
}



gsell's avatar
gsell committed
298
/**
299 300
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
301 302 303
 * @param
 * @param
 * @param
304
 */
305
void Distribution::writeToFile() {
306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
    /*
      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();
      }
      }
    */
322 323
}

324 325 326
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
327 328
}

329 330 331
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
332

333 334
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
335

336 337
void Distribution::update() {
}
adelmann's avatar
adelmann committed
338

kraus's avatar
kraus committed
339 340 341 342 343 344 345 346 347 348 349 350
void Distribution::create(size_t &numberOfParticles, double massIneV) {

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

352
    size_t mySeed = Options::seed;
adelmann's avatar
adelmann committed
353

354
    if (Options::seed == -1) {
Andreas Adelmann's avatar
Andreas Adelmann committed
355 356
        struct timeval tv;
        gettimeofday(&tv,0);
357
        mySeed = tv.tv_sec + tv.tv_usec;
Andreas Adelmann's avatar
Andreas Adelmann committed
358
    }
359

360
    if (Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE])) {
kraus's avatar
kraus committed
361
        numberOfLocalParticles = getNumOfLocalParticlesToCreate(numberOfLocalParticles);
362 363 364 365 366 367
        *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << " + core_id\n"
              << "* is scalable with number of particles and cores." << endl;
        mySeed += Ippl::myNode();
    } else {
        *gmsg << level2 << "* Generation of distribution with seed = " << mySeed << "\n"
              << "* isn't scalable with number of particles and cores." << endl;
Andreas Adelmann's avatar
Andreas Adelmann committed
368 369
    }

370 371 372
    gsl_rng_set(randGen_m, mySeed);

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

374
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
375

376
    case DistrTypeT::MATCHEDGAUSS:
kraus's avatar
kraus committed
377
        createMatchedGaussDistribution(numberOfLocalParticles, massIneV);
378
        break;
379
    case DistrTypeT::FROMFILE:
kraus's avatar
kraus committed
380
        createDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
381
        break;
382
    case DistrTypeT::GAUSS:
kraus's avatar
kraus committed
383
        createDistributionGauss(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
384
        break;
385
    case DistrTypeT::BINOMIAL:
kraus's avatar
kraus committed
386
        createDistributionBinomial(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
387
        break;
388
    case DistrTypeT::FLATTOP:
kraus's avatar
kraus committed
389
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
390
        break;
391
    case DistrTypeT::GUNGAUSSFLATTOPTH:
kraus's avatar
kraus committed
392
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
393
        break;
394
    case DistrTypeT::ASTRAFLATTOPTH:
kraus's avatar
kraus committed
395
        createDistributionFlattop(numberOfLocalParticles, massIneV);
gsell's avatar
gsell committed
396
        break;
397 398 399 400
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
kraus's avatar
kraus committed
401 402

    if (emitting_m) {
403

kraus's avatar
kraus committed
404
        unsigned int numAdditionalRNsPerParticle = 0;
405
        if (emissionModel_m == EmissionModelT::ASTRA ||
406 407
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
408

kraus's avatar
kraus committed
409
            numAdditionalRNsPerParticle = 2;
410
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
411 412 413 414 415
            if (Options::cZero) {
                numAdditionalRNsPerParticle = 40;
            } else {
                numAdditionalRNsPerParticle = 20;
            }
kraus's avatar
kraus committed
416 417
        }

418 419 420
        int saveProcessor = -1;
        const int myNode = Ippl::myNode();
        const int numNodes = Ippl::getNodes();
421
        const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
422

kraus's avatar
kraus committed
423
        for (size_t partIndex = 0; partIndex < numberOfLocalParticles; ++ partIndex) {
424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440

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

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

445 446
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
447

kraus's avatar
kraus committed
448 449 450
    // Now reflect particles if Options::cZero is true
    reflectDistribution(numberOfLocalParticles);

451
    if (Options::seed != -1)
Andreas Adelmann's avatar
Andreas Adelmann committed
452
        Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
kraus's avatar
kraus committed
453

454 455 456 457 458
    if (particlesPerDist_m.size() == 0) {
        particlesPerDist_m.push_back(tOrZDist_m.size());
    } else {
        particlesPerDist_m[0] = tOrZDist_m.size();
    }
459
}
gsell's avatar
gsell committed
460

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

463 464 465
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
466 467
        double vw = this->getVw();
        double vt = this->getvVThermal();
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
        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
492

493 494 495 496 497 498 499
                        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
500
                        } else {
501 502
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
503
                        }
504
                        beam->Bin[lowMark + count] = 0;
kraus's avatar
kraus committed
505
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
506 507 508 509 510 511
                        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
512 513
                    }
                }
514 515 516
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
517
            }
518 519 520
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
521
        }
522
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
523

524
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
525

526 527 528 529 530
        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
531

532 533 534 535
        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
536

537 538 539 540 541 542 543 544 545
                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
546
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
547 548 549 550 551 552
                        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
553 554

                    }
555

gsell's avatar
gsell committed
556
                }
557 558 559
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
560

561
            }
562

563 564 565 566 567 568
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
569

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

572
    IpplTimings::startTimer(beam.distrReload_m);
573

574 575
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
gsell's avatar
gsell committed
576

577 578 579 580
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
581

582 583
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
584

585
    beam.create(numParticles);
586

587 588
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
589 590 591

    beam.boundp();

592 593 594 595
    double actualT = beam.getT();
    long long ltstep = beam.getLocalTrackStep();
    long long gtstep = beam.getGlobalTrackStep();

gsell's avatar
gsell committed
596 597
    IpplTimings::stopTimer(beam.distrReload_m);

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

604
void Distribution::doRestartOpalCycl(PartBunch &beam,
605 606 607 608
                                     size_t Np,
                                     int restartStep,
                                     const int specifiedNumBunch,
                                     H5PartWrapper *dataSource) {
adelmann's avatar
adelmann committed
609

610
    // h5_int64_t rc;
611
    IpplTimings::startTimer(beam.distrReload_m);
612
    INFOMSG("---------------- Start reading hdf5 file----------------" << endl);
613

614 615
    long numParticles = dataSource->getNumParticles();
    size_t numParticlesPerNode = numParticles / Ippl::getNodes();
616

617 618 619 620
    size_t firstParticle = numParticlesPerNode * Ippl::myNode();
    size_t lastParticle = firstParticle + numParticlesPerNode - 1;
    if (Ippl::myNode() == Ippl::getNodes() - 1)
        lastParticle = numParticles - 1;
621

622 623
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
624

625
    beam.create(numParticles);
626

627 628
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
629

630
    beam.Q = beam.getChargePerParticle();
631

632
    beam.boundp();
633

634
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
635

636 637 638 639 640 641
    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);
642

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

646
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
647

648 649
    if(dataSource->predecessorIsSameFlavour()) {
        INFOMSG("Restart from hdf5 file generated by OPAL-cycl" << endl);
650 651
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
652 653
            energyBins_m = new PartBinsCyc(1000, beam.getNumBunch());
            beam.setPBins(energyBins_m);
654 655
        }

656 657
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
658 659 660 661 662 663 664 665 666 667 668 669 670 671

        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);
672
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
673 674 675 676 677 678 679

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

680
    INFOMSG("---------------Finished reading hdf5 file---------------" << endl);
681 682 683
    IpplTimings::stopTimer(beam.distrReload_m);
}

684
void Distribution::doRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep,
685
                                  H5PartWrapper *dataSource) {
686
    IpplTimings::startTimer(beam.distrReload_m);
687 688
    int N = dataSource->getNumParticles();
    *gmsg << "total number of slices = " << N << endl;
689 690 691 692 693 694

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

695 696
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
697 698 699 700 701 702 703 704 705 706 707 708 709 710 711

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

712
double Distribution::getTEmission() {
713 714 715 716
    if(tEmission_m > 0.0) {
        return tEmission_m;
    }

717
    setDistType();
718

719 720 721 722
    tPulseLengthFWHM_m = Attributes::getReal(itsAttr[Attrib::Distribution::TPULSEFWHM]);
    cutoff_m = Attributes::getReal(itsAttr[Attrib::Legacy::Distribution::CUTOFF]);
    tRise_m = Attributes::getReal(itsAttr[Attrib::Distribution::TRISE]);
    tFall_m = Attributes::getReal(itsAttr[Attrib::Distribution::TFALL]);
723 724 725 726 727
    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) {
728 729 730 731 732 733 734 735 736 737 738 739 740 741
    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;
742
        }
743 744 745 746 747 748 749 750 751
        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;
752 753 754 755
    }
    return tEmission_m;
}

756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786
double Distribution::getEkin() const {return Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]);}
double Distribution::getLaserEnergy() const {return Attributes::getReal(itsAttr[Attrib::Distribution::ELASER]);}
double Distribution::getWorkFunctionRf() const {return Attributes::getReal(itsAttr[Attrib::Distribution::W]);}

size_t Distribution::getNumberOfDarkCurrentParticles() { return (size_t) Attributes::getReal(itsAttr[Attrib::Distribution::NPDARKCUR]);}
double Distribution::getDarkCurrentParticlesInwardMargin() { return Attributes::getReal(itsAttr[Attrib::Distribution::INWARDMARGIN]);}
double Distribution::getEInitThreshold() { return Attributes::getReal(itsAttr[Attrib::Distribution::EINITHR]);}
double Distribution::getWorkFunction() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNPHIW]); }
double Distribution::getFieldEnhancement() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNBETA]); }
size_t Distribution::getMaxFNemissionPartPerTri() { return (size_t) Attributes::getReal(itsAttr[Attrib::Distribution::FNMAXEMI]);}
double Distribution::getFieldFNThreshold() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNFIELDTHR]);}
double Distribution::getFNParameterA() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNA]);}
double Distribution::getFNParameterB() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNB]);}
double Distribution::getFNParameterY() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNY]);}
double Distribution::getFNParameterVYZero() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNVYZERO]);}
double Distribution::getFNParameterVYSecond() { return Attributes::getReal(itsAttr[Attrib::Distribution::FNVYSECOND]);}
int    Distribution::getSecondaryEmissionFlag() { return Attributes::getReal(itsAttr[Attrib::Distribution::SECONDARYFLAG]);}
bool   Distribution::getEmissionMode() { return Attributes::getBool(itsAttr[Attrib::Distribution::NEMISSIONMODE]);}

std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[Attrib::Distribution::TYPE]);}

double Distribution::getvSeyZero() {return Attributes::getReal(itsAttr[Attrib::Distribution::VSEYZERO]);}// return sey_0 in Vaughan's model
double Distribution::getvEZero() {return Attributes::getReal(itsAttr[Attrib::Distribution::VEZERO]);}// return the energy related to sey_0 in Vaughan's model
double Distribution::getvSeyMax() {return Attributes::getReal(itsAttr[Attrib::Distribution::VSEYMAX]);}// return sey max in Vaughan's model
double Distribution::getvEmax() {return Attributes::getReal(itsAttr[Attrib::Distribution::VEMAX]);}// return Emax in Vaughan's model
double Distribution::getvKenergy() {return Attributes::getReal(itsAttr[Attrib::Distribution::VKENERGY]);}// return fitting parameter denotes the roughness of surface for impact energy in Vaughan's model
double Distribution::getvKtheta() {return Attributes::getReal(itsAttr[Attrib::Distribution::VKTHETA]);}// return fitting parameter denotes the roughness of surface for impact angle in Vaughan's model
double Distribution::getvVThermal() {return Attributes::getReal(itsAttr[Attrib::Distribution::VVTHERMAL]);}// thermal velocity of Maxwellian distribution of secondaries in Vaughan's model
double Distribution::getVw() {return Attributes::getReal(itsAttr[Attrib::Distribution::VW]);}// velocity scalar for parallel plate benchmark;

int Distribution::getSurfMaterial() {return (int)Attributes::getReal(itsAttr[Attrib::Distribution::SURFMATERIAL]);}// Surface material number for Furman-Pivi's Model;
787 788 789

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

790 791
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
792
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
793
        os << "* In restart. Distribution read in from .h5 file." << endl;
794 795
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
796
            os << "* Main Distribution" << endl
797
               << "-----------------" << endl;
798 799

        if (particlesPerDist_m.empty())
800
            printDist(os, 0);
801
        else
802
            printDist(os, particlesPerDist_m.at(0));
803 804 805

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
806
            os << "* " << endl;
adelmann's avatar
adelmann committed
807 808
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
809
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
810 811 812
            distCount++;
        }

813
        os << "* " << endl;
814
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
815
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
816

817
            //            if (numberOfEnergyBins_m > 1)
818
            //    printEnergyBins(os);
819 820 821
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
822 823
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
824 825
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
826
            os << "* " << endl;
827
            printEmissionModel(os);
828
            os << "* " << endl;
829
        } else
adelmann's avatar
adelmann committed
830
            os << "* Distribution is injected." << endl;
831
    }
832 833
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
834 835 836 837

    return os;
}

838
const PartData &Distribution::getReference() const {
839 840 841 842 843
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

844
void Distribution::addDistributions() {
845 846 847
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
848

kraus's avatar
kraus committed
849
    size_t idx = 1;
850 851
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
kraus's avatar
kraus committed
852 853 854
         addedDistIt != addedDistributions_m.end(); ++ addedDistIt, ++ idx) {

        particlesPerDist_m[idx] = (*addedDistIt)->tOrZDist_m.size();
855 856

        std::vector<double>::iterator distIt;
857 858
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
859 860 861
             distIt++) {
            xDist_m.push_back(*distIt);
        }
862
        (*addedDistIt)->eraseXDist();
863

864 865
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
866 867 868
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
869
        (*addedDistIt)->eraseBGxDist();
870

871 872
        for (distIt = (*addedDistIt)->getYDist().begin();
             distIt != (*addedDistIt)->getYDist().end();
873 874 875
             distIt++) {
            yDist_m.push_back(*distIt);
        }
876
        (*addedDistIt)->eraseYDist();
877

878 879
        for (distIt = (*addedDistIt)->getBGyDist().begin();
             distIt != (*addedDistIt)->getBGyDist().end();
880 881 882
             distIt++) {
            pyDist_m.push_back(*distIt);
        }
883
        (*addedDistIt)->eraseBGyDist();
884

885 886
        for (distIt = (*addedDistIt)->getTOrZDist().begin();
             distIt != (*addedDistIt)->getTOrZDist().end();
887 888 889
             distIt++) {
            tOrZDist_m.push_back(*distIt);
        }
890
        (*addedDistIt)->eraseTOrZDist();
891

892 893
        for (distIt = (*addedDistIt)->getBGzDist().begin();
             distIt != (*addedDistIt)->getBGzDist().end();
894 895 896
             distIt++) {
            pzDist_m.push_back(*distIt);
        }
897
        (*addedDistIt)->eraseBGzDist();
898 899 900
    }
}

901
void Distribution::applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs) {
902 903 904 905

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
906
        applyEmissModelNone(pz);
907 908
        break;
    case EmissionModelT::ASTRA:
909
        applyEmissModelAstra(px, py, pz, additionalRNs);
910 911
        break;
    case EmissionModelT::NONEQUIL:
912
        applyEmissModelNonEquil(lowEnergyLimit, px, py, pz, additionalRNs);
913 914 915 916 917 918
        break;
    default:
        break;
    }
}