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

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

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

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

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

50

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

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

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

58 59 60 61 62 63 64 65
SymTenzor<double, 6> getUnit6x6() {
    SymTenzor<double, 6> unit6x6;
    for (unsigned int i = 0; i < 6u; ++ i) {
        unit6x6(i,i) = 1.0;
    }
    return unit6x6;
}

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

70 71 72
namespace AttributesT
{
    enum AttributesT {
73
        TYPE,
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
        DISTRIBUTION,
        FNAME,
        WRITETOFILE,
        WEIGHT,
        INPUTMOUNITS,
        EMITTED,
        EMISSIONSTEPS,
        EMISSIONMODEL,
        EKIN,
        ELASER,
        W,
        FE,
        CATHTEMP,
        NBIN,
        XMULT,
        YMULT,
        ZMULT,
        TMULT,
        PXMULT,
        PYMULT,
        PZMULT,
        OFFSETX,
        OFFSETY,
        OFFSETZ,
        OFFSETT,
        OFFSETPX,
        OFFSETPY,
        OFFSETPZ,
        SIGMAX,
        SIGMAY,
        SIGMAR,
        SIGMAZ,
        SIGMAT,
        TPULSEFWHM,
        TRISE,
        TFALL,
        SIGMAPX,
        SIGMAPY,
        SIGMAPZ,
        MX,
        MY,
        MZ,
        MT,
        CUTOFFX,
        CUTOFFY,
        CUTOFFR,
        CUTOFFLONG,
        CUTOFFPX,
        CUTOFFPY,
        CUTOFFPZ,
        FTOSCAMPLITUDE,
        FTOSCPERIODS,
        R,                          // the correlation matrix (a la transport)
        CORRX,
        CORRY,
        CORRZ,
        CORRT,
        R51,
        R52,
        R61,
        R62,
        LASERPROFFN,
        IMAGENAME,
        INTENSITYCUT,
138 139 140 141 142
        FLIPX,
        FLIPY,
        ROTATE90,
        ROTATE180,
        ROTATE270,
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
        NPDARKCUR,
        INWARDMARGIN,
        EINITHR,
        FNA,
        FNB,
        FNY,
        FNVYZERO,
        FNVYSECOND,
        FNPHIW,
        FNBETA,
        FNFIELDTHR,
        FNMAXEMI,
        SECONDARYFLAG,
        NEMISSIONMODE,
        VSEYZERO,                   // sey_0 in Vaughn's model.
        VEZERO,                     // Energy related to sey_0 in Vaughan's model.
        VSEYMAX,                    // sey max in Vaughan's model.
        VEMAX,                      // Emax in Vaughan's model.
        VKENERGY,                   // Fitting parameter denotes the roughness of
        // surface for impact energy in Vaughn's model.
        VKTHETA,                    // Fitting parameter denotes the roughness of
        // surface for impact angle in Vaughn's model.
        VVTHERMAL,                  // Thermal velocity of Maxwellian distribution
        // of secondaries in Vaughan's model.
        VW,
        SURFMATERIAL,               // Add material type, currently 0 for copper
        // and 1 for stainless steel.
        EX,                         // below is for the matched distribution
        EY,
        ET,
        MAGSYM,                     // number of sector magnets
        LINE,
        FMAPFN,
frey_m's avatar
frey_m committed
176
        FMTYPE,                     // field map type used in matched gauss distribution
177 178 179 180 181
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
182
        RGUESS,
183 184
        ID1,                       // special particle that the user can set
        ID2,                       // special particle that the user can set
185
        SIZE
gsell's avatar
gsell committed
186 187 188
    };
}

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

gsell's avatar
gsell committed
216
Distribution::Distribution():
217 218
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
219
    distrTypeT_m(DistrTypeT::NODIST),
220
    numberOfDistributions_m(1),
221 222 223 224 225 226
    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
227
    currentEnergyBin_m(1),
228
    currentSampleBin_m(0),
229 230 231 232
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
233
    randGen_m(NULL),
234
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
235
    pmean_m(0.0),
236 237 238 239 240
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
241 242
    totalNumberParticles_m(0),
    totalNumberEmittedParticles_m(0),
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
    inputMoUnits_m(InputMomentumUnitsT::NONE),
    sigmaTRise_m(0.0),
    sigmaTFall_m(0.0),
    tPulseLengthFWHM_m(0.0),
    correlationMatrix_m(getUnit6x6()),
    laserProfileFileName_m(""),
    laserImageName_m(""),
    laserIntensityCut_m(0.0),
    laserProfile_m(NULL),
    darkCurrentParts_m(0),
    darkInwardMargin_m(0.0),
    eInitThreshold_m(0.0),
    workFunction_m(0.0),
    fieldEnhancement_m(0.0),
    fieldThrFN_m(0.0),
    maxFN_m(0),
    paraFNA_m(0.0),
    paraFNB_m(0.0),
    paraFNY_m(0.0),
    paraFNVYSe_m(0.0),
    paraFNVYZe_m(0.0),
    secondaryFlag_m(0),
    ppVw_m(0.0),
    vVThermal_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    avrgpz_m(0.0),
    I_m(0.0),
    E_m(0.0),
    bega_m(0.0),
    M_m(0.0)
274
{
275
    setAttributes();
276 277 278

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

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

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

Distribution::~Distribution() {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if (emitting_m) {
495

kraus's avatar
kraus committed
496
        unsigned int numAdditionalRNsPerParticle = 0;
497
        if (emissionModel_m == EmissionModelT::ASTRA ||
498 499 500
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {

kraus's avatar
kraus committed
501
            numAdditionalRNsPerParticle = 2;
502 503
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
            numAdditionalRNsPerParticle = 20;
kraus's avatar
kraus committed
504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521
        }

        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) {
522
                    double x = gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
523 524 525 526 527
                    rns.push_back(x);
                }
                additionalRNs_m.push_back(rns);
            } else {
                for (unsigned int i = 0; i < numAdditionalRNsPerParticle; ++ i) {
528
                    gsl_rng_uniform(randGen_m);
kraus's avatar
kraus committed
529 530 531 532 533
                }
            }
        }
    }

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

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

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

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

602
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
603

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

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

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

                    }
633

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

639
            }
640

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

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

650
    IpplTimings::startTimer(beam.distrReload_m);
651

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

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

660 661
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
662

663
    beam.create(numParticles);
664

665 666
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
667 668 669

    beam.boundp();

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

gsell's avatar
gsell committed
674 675
    IpplTimings::stopTimer(beam.distrReload_m);

676 677 678 679
    *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
680 681
}

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

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

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

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

700 701
    numParticles = lastParticle - firstParticle + 1;
    PAssert(numParticles >= 0);
702

703
    beam.create(numParticles);
704

705 706
    dataSource->readHeader();
    dataSource->readStep(beam, firstParticle, lastParticle);
707

708
    beam.Q = beam.getChargePerParticle();
709

710
    beam.boundp();
711

712
    double meanE = static_cast<H5PartWrapperForPC*>(dataSource)->getMeanKineticEnergy();
713

714 715 716 717 718 719
    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);
720

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

724
    INFOMSG("* Gamma = " << gamma << ", Beta = " << beta << endl);
725

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

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

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

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

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

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

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

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

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

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

795
    setDistType();
796 797

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

834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852
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]);}

853
std::string Distribution::getTypeofDistribution() { return (std::string) Attributes::getString(itsAttr[AttributesT::TYPE]);}
854 855 856 857 858 859 860 861 862 863 864

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;
865 866 867

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

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

        if (particlesPerDist_m.empty())
878
            printDist(os, 0);
879
        else
880
            printDist(os, particlesPerDist_m.at(0));
881 882 883

        size_t distCount = 1;
        for (unsigned distIndex = 0; distIndex < addedDistributions_m.size(); distIndex++) {
884
            os << "* " << endl;
adelmann's avatar
adelmann committed
885 886
            os << "* Added Distribution #" << distCount << endl;
            os << "* ----------------------" << endl;
887
            addedDistributions_m.at(distIndex)->printDist(os, particlesPerDist_m.at(distCount));
888 889 890
            distCount++;
        }

891
        os << "* " << endl;
892
        if (numberOfEnergyBins_m > 0) {
adelmann's avatar
adelmann committed
893
            os << "* Number of energy bins    = " << numberOfEnergyBins_m << endl;
894

895
            //            if (numberOfEnergyBins_m > 1)
896
            //    printEnergyBins(os);
897 898 899
        }

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

    return os;
}

916
const PartData &Distribution::getReference() const {
917 918 919 920 921
    // Cast away const, to allow logically constant Distribution to update.
    const_cast<Distribution *>(this)->update();
    return particleRefData_m;
}

922
void Distribution::addDistributions() {
923 924 925
    /*
     * Move particle coordinates from added distributions to main distribution.
     */
adelmann's avatar
doc  
adelmann committed
926

927 928 929 930 931
    std::vector<Distribution *>::iterator addedDistIt;
    for (addedDistIt = addedDistributions_m.begin();
         addedDistIt != addedDistributions_m.end(); addedDistIt++) {

        std::vector<double>::iterator distIt;
932 933
        for (distIt = (*addedDistIt)->getXDist().begin();
             distIt != (*addedDistIt)->getXDist().end();
934 935 936
             distIt++) {
            xDist_m.push_back(*distIt);
        }
937
        (*addedDistIt)->eraseXDist();
938

939 940
        for (distIt = (*addedDistIt)->getBGxDist().begin();
             distIt != (*addedDistIt)->getBGxDist().end();
941 942 943
             distIt++) {
            pxDist_m.push_back(*distIt);
        }
944
        (*addedDistIt)->eraseBGxDist();
945

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

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

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

967 968
        for (distIt = (*addedDistIt)->getBGzDist().begin();
             distIt != (*addedDistIt)->getBGzDist().end();
969 970 971
             distIt++) {
            pzDist_m.push_back(*distIt);
        }
972
        (*addedDistIt)->eraseBGzDist();
973 974 975
    }
}

976
void Distribution::applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs) {
977 978 979 980

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
981
        applyEmissModelNone(pz);
982 983
        break;
    case EmissionModelT::ASTRA:
984
        applyEmissModelAstra(px, py, pz, additionalRNs);
985 986
        break;
    case EmissionModelT::NONEQUIL:
987
        applyEmissModelNonEquil(lowEnergyLimit, px, py, pz, additionalRNs);
988 989 990 991 992 993
        break;
    default:
        break;
    }
}

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