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

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

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

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

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

50
#include <sys/time.h>
51

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

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

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

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

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

189 190 191
namespace LegacyAttributesT
{
    enum LegacyAttributesT {
192
        // DESCRIPTION OF THE DISTRIBUTION:
193
        DEBIN = AttributesT::SIZE,
194 195 196 197 198 199 200
        SBIN,
        SIGMAPT,
        CUTOFF,
        T,
        PT,
        SIZE
    };
201 202
}

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

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

    try {
268
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
269
    } catch(...) {
270 271 272
        delete defaultDistribution;
    }

273
    setFieldEmissionParameters();
gsell's avatar
gsell committed
274 275 276 277 278 279 280
}
/**
 *
 *
 * @param name
 * @param parent
 */
281
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
282
    Definition(name, parent),
283 284
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
285
    numberOfDistributions_m(parent->numberOfDistributions_m),
286 287 288 289 290 291
    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
292
    tEmission_m(parent->tEmission_m),
293 294 295 296 297 298 299 300
    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),
301
    randGen_m(NULL),
302
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
303
    pmean_m(parent->pmean_m),
304 305 306 307 308
    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),
309 310
    totalNumberParticles_m(parent->totalNumberParticles_m),
    totalNumberEmittedParticles_m(parent->totalNumberEmittedParticles_m),
311 312 313 314 315 316 317 318 319 320 321 322
    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),
323
    avrgpz_m(parent->avrgpz_m),
324 325 326 327 328 329 330 331
    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),
332
    correlationMatrix_m(parent->correlationMatrix_m),
333 334 335 336
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
337 338 339 340 341 342 343 344 345 346 347 348
    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),
349 350 351
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
352 353 354
    I_m(parent->I_m),
    E_m(parent->E_m),
    M_m(parent->M_m),
355 356 357 358
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
359
    cutoff_m(parent->cutoff_m),
Andreas Adelmann's avatar
Andreas Adelmann committed
360 361
    bega_m(parent->bega_m),
    mySeed_m(parent->mySeed_m)
362
{
Andreas Adelmann's avatar
Andreas Adelmann committed
363

gsell's avatar
gsell committed
364 365 366 367 368 369 370
}

Distribution::~Distribution() {

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

371 372 373
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
374 375
    }

376 377 378 379 380
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

381 382 383
    if (randGen_m != NULL) {
        gsl_rng_free(randGen_m);
        randGen_m = NULL;
384 385 386 387 388 389
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
390 391

}
392
/*
393
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
394
  for(int i=0; i<M.size1(); ++i) {
395 396 397 398 399
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
400 401
  }
*/
gsell's avatar
gsell committed
402

403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425
/**
 * Calculate the local number of particles evenly and adjust node 0
 * such that n is matched exactly. 
 * @param n total number of particles 
 * @return n / #cores 
 * @param
 */

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

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

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

    return locNumber;
}



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

452 453 454
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
455 456
}

457 458 459
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
460

461 462
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
463

464 465
void Distribution::update() {
}
adelmann's avatar
adelmann committed
466

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

Andreas Adelmann's avatar
Andreas Adelmann committed
469
    setFieldEmissionParameters(); /// is also called in the constructor
adelmann's avatar
adelmann committed
470

471 472
    size_t locNumber = getNumOfLocalParticlesToCreate(numberOfParticles);

Andreas Adelmann's avatar
Andreas Adelmann committed
473 474 475 476 477 478 479 480 481 482 483 484 485 486 487
    if (Options::seed ==-1) {
        struct timeval tv;
        gettimeofday(&tv,0);
        mySeed_m = tv.tv_sec + tv.tv_usec + Ippl::myNode();
        *gmsg << "* Distribution generation with non-portable rng" << endl;
    }
    else {
        mySeed_m = Ippl::myNode();
        *gmsg << "* Distribution generation with portable rng" << endl;
    }

    gsl_rng_env_setup();
    randGen_m = gsl_rng_alloc(gsl_rng_default);
    gsl_rng_set(randGen_m, mySeed_m);

488
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
489

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

    if (emitting_m) {
517

kraus's avatar
kraus committed
518
        unsigned int numAdditionalRNsPerParticle = 0;
519
        if (emissionModel_m == EmissionModelT::ASTRA ||
520 521
            distrTypeT_m == DistrTypeT::ASTRAFLATTOPTH ||
            distrTypeT_m == DistrTypeT::GUNGAUSSFLATTOPTH) {
kraus's avatar
kraus committed
522
            numAdditionalRNsPerParticle = 2;
523 524
        } else if (emissionModel_m == EmissionModelT::NONEQUIL) {
            numAdditionalRNsPerParticle = 20;
kraus's avatar
kraus committed
525 526 527 528 529
        }

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

541 542
    // Scale coordinates according to distribution input.
    scaleDistCoordinates();
Andreas Adelmann's avatar
Andreas Adelmann committed
543 544
    if (Options::seed != -1)        
        Options::seed = gsl_rng_uniform_int(randGen_m, gsl_rng_max(randGen_m));
545
}
gsell's avatar
gsell committed
546

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

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

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

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

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

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

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

                    }
641

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

647
            }
648

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

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

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

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

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

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

671
    beam.create(numParticles);
672

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

    beam.boundp();

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

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

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

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

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

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

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

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

711
    beam.create(numParticles);
712

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

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

718
    beam.boundp();
719

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

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

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

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

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

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

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

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

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

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

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

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

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

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

803
    setDistType();
804 805

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

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

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

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

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

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

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

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

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

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

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

    return os;
}

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

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

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

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

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

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

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

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

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

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

    switch (emissionModel_m) {

    case EmissionModelT::NONE:
989
        applyEmissModelNone(pz);
990 991
        break;
    case EmissionModelT::ASTRA: