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

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

49
#include <sys/time.h>
50

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

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

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

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

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

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

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

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

141
    // setFieldEmissionParameters();
142 143 144

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

Distribution::~Distribution() {

239
    if ((Ippl::getNodes() == 1) && (os_m.is_open()))
gsell's avatar
gsell committed
240 241
        os_m.close();

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

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

252 253 254
    if (randGen_m != NULL) {
        gsl_rng_free(randGen_m);
        randGen_m = NULL;
255 256
    }

257
    if (laserProfile_m) {
258 259 260
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
261 262

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

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

size_t Distribution::getNumOfLocalParticlesToCreate(size_t n) {

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

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

    return locNumber;
}



gsell's avatar
gsell committed
297
/**
298 299
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
300 301 302
 * @param
 * @param
 * @param
303
 */
304
void Distribution::writeToFile() {
305
    /*
306 307
      if (Ippl::getNodes() == 1) {
      if (os_m.is_open()) {
308 309 310 311 312
      ;
      } else {
      *gmsg << " Write distribution to file data/dist.dat" << endl;
      std::string file("data/dist.dat");
      os_m.open(file.c_str());
313
      if (os_m.bad()) {
314 315 316 317 318 319 320
      *gmsg << "Unable to open output file " <<  file << endl;
      }
      os_m << "# x y ti px py pz "  << std::endl;
      os_m.close();
      }
      }
    */
321 322
}

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

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

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

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

kraus's avatar
kraus committed
338 339 340 341 342 343 344 345 346 347 348 349
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
350

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

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

359
    if (Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE])) {
kraus's avatar
kraus committed
360
        numberOfLocalParticles = getNumOfLocalParticlesToCreate(numberOfLocalParticles);
361 362 363 364 365 366
        *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
367 368
    }

369 370 371
    gsl_rng_set(randGen_m, mySeed);

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

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

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

    if (emitting_m) {
402

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

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

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

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

            // 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
440 441 442 443
            }
        }
    }

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

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

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

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

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

462
    if (Options::ppdebug) {  // This is Parallel Plate Benchmark.
463 464
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
465 466
        double vw = this->getVw();
        double vt = this->getvVThermal();
467 468 469 470 471 472
        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());
473
        if (Ippl::myNode() == 0)
474
            N_mean += N_extra;
475 476 477 478
        if (bg.getN() != 0) {
            for (size_t i = 0; i < bg.getN(); i++) {
                if (pc == Ippl::myNode()) {
                    if (count < N_mean) {
479 480 481 482 483
                        /*==============Parallel Plate Benchmark=====================================*/
                        double test_s = 1;
                        double f_x = 0;
                        double test_x = 0;
                        while(test_s > f_x) {
484
                            test_s = gsl_rng_uniform(randGen_m);
485
                            test_s *= f_max;
486
                            test_x = gsl_rng_uniform(randGen_m);
487 488 489 490
                            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
491

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

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

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

531
        if (Ippl::myNode() == 0)
532
            N_mean += N_extra;
533 534
        if (bg.getN() != 0) {
            for (size_t i = 0; i < bg.getN(); i++) {
gsell's avatar
gsell committed
535

536 537
                if (pc == Ippl::myNode()) {
                    if (count < N_mean) {
538
                        beam->create(1);
539
                        if (pc != 0)
540 541 542 543 544
                            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
545
                        beam->PType[lowMark + count] = ParticleType::REGULAR;
546 547 548 549 550 551
                        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
552 553

                    }
554

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

560
            }
561

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

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

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

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

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

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

584
    beam.create(numParticles);
585

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

    beam.boundp();

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

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

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

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

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

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

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

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

624
    beam.create(numParticles);
625

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

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

631
    beam.boundp();
632

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

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

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

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

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

655 656
    } else {
        INFOMSG("Restart from hdf5 file generated by OPAL-t" << endl);
657 658 659 660

        Vector_t meanR(0.0, 0.0, 0.0);
        Vector_t meanP(0.0, 0.0, 0.0);
        unsigned long int newLocalN = beam.getLocalNum();
661 662
        for (unsigned int i = 0; i < newLocalN; ++i) {
            for (int d = 0; d < 3; ++d) {
663 664 665 666 667 668 669 670
                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);
671
        INFOMSG("Rmean = " << meanR << "[m], Pmean=" << meanP << endl);
672

673
        for (unsigned int i = 0; i < newLocalN; ++i) {
674 675 676 677 678
            beam.R[i] -= meanR;
            beam.P[i] -= meanP;
        }
    }

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

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

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

694 695
    dataSource->readHeader();
    dataSource->readStep(beam, starti, endi);
696 697 698 699 700 701 702 703

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

704
    if (dist == 0) {
705 706 707 708 709 710
        throw OpalException("Distribution::find()", "Distribution \"" + name + "\" not found.");
    }

    return dist;
}

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

716
    setDistType();
717

718 719 720 721
    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]);
722 723 724 725 726
    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) {
727 728 729 730 731 732 733 734
    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;
735
        for (int i = 0; i < 10; ++ i) {
736 737 738 739 740
            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;
741
        }
742 743 744 745 746 747 748 749 750
        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;
751 752 753 754 755 756
    }
    return tEmission_m;
}

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

757 758
    os << "* ************* D I S T R I B U T I O N ********************************************" << endl;
    os << "* " << endl;
759
    if (OpalData::getInstance()->inRestartRun()) {
adelmann's avatar
adelmann committed
760
        os << "* In restart. Distribution read in from .h5 file." << endl;
761 762
    } else {
        if (addedDistributions_m.size() > 0)
adelmann's avatar
adelmann committed
763
            os << "* Main Distribution" << endl
764
               << "-----------------" << endl;
765 766

        if (particlesPerDist_m.empty())
767
            printDist(os, 0);
768
        else
769
            printDist(os, particlesPerDist_m.at(0));
770 771 772

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
773
            os << "* " << endl;
adelmann's avatar
adelmann committed
774 775
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
776
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
777 778 779
            distCount++;
        }

780
        os << "* " << endl;
781
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
782
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
783

784
            //            if (numberOfEnergyBins_m > 1)
785
            //    printEnergyBins(os);
786 787 788
        }

        if (emitting_m) {
adelmann's avatar
adelmann committed
789 790
            os << "* Distribution is emitted. " << endl;
            os << "* Emission time            = " << tEmission_m << " [sec]" << endl;
791 792
            os << "* Time per bin             = " << tEmission_m / numberOfEnergyBins_m << " [sec]" << endl;
            os << "* Delta t during emission  = " << tBin_m / numberOfSampleBins_m << " [sec]" << endl;
793
            os << "* " << endl;
794
            printEmissionModel(os);
795
            os << "* " << endl;
796
        } else
adelmann's avatar
adelmann committed
797
            os << "* Distribution is injected." << endl;
798
    }
799 800
    os << "* " << endl;
    os << "* **********************************************************************************" << endl;
801 802 803 804

    return os;
}

805
const PartData &Distribution::getReference() const {
806 807 808 809 810
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

811
void Distribution::addDistributions() {
812 813 814
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
815

kraus's avatar
kraus committed
816
    size_t idx = 1;
817 818
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
kraus's avatar
kraus committed
819 820 821
         addedDistIt != addedDistributions_m.end(); ++ addedDistIt, ++ idx) {

        particlesPerDist_m[idx] = (*addedDistIt)->tOrZDist_m.size();
822 823

        std::vector<double>::iterator distIt;
824 825
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
826 827 828
             distIt++) {
            xDist_m.push_back(*distIt);
        }
829
        (*addedDistIt)->eraseXDist();
830

831 832
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
833 834 835
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
836
        (*addedDistIt)->eraseBGxDist();
837

838 839
        for (distIt = (*addedDistIt)->getYDist().begin();
             distIt != (*addedDistIt)->getYDist().end();
840 841 842
             distIt++) {
            yDist_m.push_back(*distIt);
        }
843
        (*addedDistIt)->eraseYDist();
844

845 846
        for (distIt = (*addedDistIt)->getBGyDist().begin();
             distIt != (*addedDistIt)->getBGyDist().end();
847 848 849
             distIt++) {
            pyDist_m.push_back(*distIt);
        }
850
        (*addedDistIt)->eraseBGyDist();
851

852 853
        for (distIt = (*addedDistIt)->getTOrZDist().begin();
             distIt != (*addedDistIt)->getTOrZDist().end();
854 855 856
             distIt++) {
            tOrZDist_m.push_back(*distIt);
        }
857
        (*addedDistIt)->eraseTOrZDist();
858

859 860
        for (distIt = (*addedDistIt)->getBGzDist().begin();
             distIt != (*addedDistIt)->getBGzDist().end();