Distribution.h 16.1 KB
Newer Older
gsell's avatar
gsell committed
1 2 3
#ifndef OPAL_Distribution_HH
#define OPAL_Distribution_HH

4
// Distribution class
gsell's avatar
gsell committed
5
//
6 7 8
// Copyright (c) 2008-2020
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
gsell's avatar
gsell committed
9
//
10 11
// OPAL is licensed under GNU GPL version 3.

gsell's avatar
gsell committed
12 13
#include <fstream>
#include <string>
14
#include <vector>
gsell's avatar
gsell committed
15

16 17
#include "AbstractObjects/Definition.h"
#include "Algorithms/PartData.h"
18
#include "Algorithms/Vektor.h"
19
#include "AppTypes/SymTenzor.h"
20
#include "Attributes/Attributes.h"
gsell's avatar
gsell committed
21

22 23
#include "Distribution/SigmaGenerator.h"

24
#include <gsl/gsl_histogram.h>
25 26
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
27

28 29 30
#ifdef WITH_UNIT_TESTS
#include <gtest/gtest_prod.h>
#endif
31

32
class Beam;
33
class Beamline;
frey_m's avatar
frey_m committed
34 35 36 37

template <class T, unsigned Dim>
class PartBunchBase;

38 39
class PartBins;
class LaserProfile;
40
class H5PartWrapper;
gsell's avatar
gsell committed
41

42 43 44
namespace DistrTypeT
{
    enum DistrTypeT {NODIST,
45 46 47 48
                     FROMFILE,
                     GAUSS,
                     BINOMIAL,
                     FLATTOP,
49
                     MULTIGAUSS,
50 51 52
                     GUNGAUSSFLATTOPTH,
                     ASTRAFLATTOPTH,
                     MATCHEDGAUSS
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
                    };
}

namespace EmissionModelT
{
    enum EmissionModelT {NONE,
                         ASTRA,
                         NONEQUIL
                        };
}

namespace InputMomentumUnitsT
{
    enum InputMomentumUnitsT {NONE,
                              EV
                              };
}

71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
namespace Attrib
{
    namespace Distribution
    {
        enum AttributesT {
            TYPE,
            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,
104
            OFFSETP,
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
            SIGMAX,
            SIGMAY,
            SIGMAR,
            SIGMAZ,
            SIGMAT,
            TPULSEFWHM,
            TRISE,
            TFALL,
            SIGMAPX,
            SIGMAPY,
            SIGMAPZ,
            MX,
            MY,
            MZ,
            MT,
            CUTOFFX,
            CUTOFFY,
            CUTOFFR,
            CUTOFFLONG,
            CUTOFFPX,
            CUTOFFPY,
            CUTOFFPZ,
            FTOSCAMPLITUDE,
            FTOSCPERIODS,
129 130
            SEPPEAKS,
            NPEAKS,
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
            R,                          // the correlation matrix (a la transport)
            CORRX,
            CORRY,
            CORRZ,
            CORRT,
            R51,
            R52,
            R61,
            R62,
            LASERPROFFN,
            IMAGENAME,
            INTENSITYCUT,
            FLIPX,
            FLIPY,
            ROTATE90,
            ROTATE180,
            ROTATE270,
            EX,                         // below is for the matched distribution
            EY,
            ET,
151
            SECTOR,                     // Matched-Gauss: single sector or full machine
152
            NSECTORS,                   // Matched-Gauss: number of sectors to average field
153
            NSTEPS,                     // Matched-Gauss: number of steps for closed orbit finder
154 155
            RGUESS,                     // Matched-Gauss: guess for closed orbit finder
            DENERGY,                    // Matched-Gauss: energy step for closed orbit finder
156 157 158 159 160
            LINE,
            RESIDUUM,
            MAXSTEPSCO,
            MAXSTEPSSI,
            ORDERMAPS,
161
            //            E2,
162 163 164 165 166 167 168 169 170 171 172 173 174
            ID1,                       // special particle that the user can set
            ID2,                       // special particle that the user can set
            SCALABLE,
            SIZE
        };
    }

    namespace Legacy
    {
        namespace Distribution
        {
            enum LegacyAttributesT {
                // DESCRIPTION OF THE DISTRIBUTION:
175
                DISTRIBUTION = Attrib::Distribution::SIZE,
176
                // DEBIN,
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
                SBIN,
                SIGMAPT,
                CUTOFF,
                T,
                PT,
                // ALPHAX,
                // ALPHAY,
                // BETAX,
                // BETAY,
                // DX,
                // DDX,
                // DY,
                // DDY,
                SIZE
            };
        }
    }
}

196 197 198 199 200
/*
 * Class Distribution
 *
 * Defines the initial beam that is injected or emitted into the simulation.
 */
gsell's avatar
gsell committed
201 202 203 204 205 206 207 208 209

class Distribution: public Definition {

public:

    Distribution();
    virtual ~Distribution();

    virtual bool canReplaceBy(Object *object);
210
    virtual Distribution *clone(const std::string &name);
gsell's avatar
gsell committed
211 212
    virtual void execute();
    virtual void update();
213
    size_t getNumOfLocalParticlesToCreate(size_t n);
frey_m's avatar
frey_m committed
214
    void createOpalCycl(PartBunchBase<double, 3> *beam,
215
                        size_t numberOfParticles,
kraus's avatar
kraus committed
216
                        double current, const Beamline &bl);
frey_m's avatar
frey_m committed
217
    void createOpalT(PartBunchBase<double, 3> *beam,
218
                     std::vector<Distribution *> addedDistributions,
219 220
                     size_t &numberOfParticles);
    void createOpalT(PartBunchBase<double, 3> *beam, size_t &numberOfParticles);
frey_m's avatar
frey_m committed
221 222
    void doRestartOpalT(PartBunchBase<double, 3> *p, size_t Np, int restartStep, H5PartWrapper *h5wrapper);
    void doRestartOpalCycl(PartBunchBase<double, 3> *p, size_t Np, int restartStep,
223
                        const int specifiedNumBunch, H5PartWrapper *h5wrapper);
frey_m's avatar
frey_m committed
224
    size_t emitParticles(PartBunchBase<double, 3> *beam, double eZ);
225
    double getPercentageEmitted() const;
226
    static Distribution *find(const std::string &name);
gsell's avatar
gsell committed
227

228 229 230 231 232 233 234 235 236
    bool getIfDistEmitting();
    int getLastEmittedEnergyBin();
    double getMaxTOrZ();
    double getMinTOrZ();
    size_t getNumberOfEmissionSteps();
    int getNumberOfEnergyBins();
    double getEmissionDeltaT();
    double getEnergyBinDeltaT();
    double getWeight();
gsell's avatar
gsell committed
237

238
    double getTEmission();
239

kraus's avatar
kraus committed
240
    Vector_t get_pmean() const;
241 242 243

    std::string getTypeofDistribution();

gsell's avatar
gsell committed
244 245
    Inform &printInfo(Inform &os) const;

246
    bool Rebin();
247 248
    void setDistToEmitted(bool emitted);
    void setDistType();
249 250
    void setSigmaR_m();
    void setSigmaP_m(double massIneV);
251 252
    void shiftBeam(double &maxTOrZ, double &minTOrZ);
    double getEmissionTimeShift() const;
gsell's avatar
gsell committed
253

254
    void setNumberOfDistributions(unsigned int n) { numberOfDistributions_m = n; }
255 256

    DistrTypeT::DistrTypeT getType() const;
257
private:
258
#ifdef WITH_UNIT_TESTS
259 260
    FRIEND_TEST(GaussTest, FullSigmaTest1);
    FRIEND_TEST(GaussTest, FullSigmaTest2);
261 262
    FRIEND_TEST(BinomialTest, FullSigmaTest1);
    FRIEND_TEST(BinomialTest, FullSigmaTest2);
263
#endif
264

265
    Distribution(const std::string &name, Distribution *parent);
266

gsell's avatar
gsell committed
267
    // Not implemented.
268 269
    Distribution(const Distribution &) = delete;
    void operator=(const Distribution &) = delete;
270

271 272 273 274 275 276 277 278 279 280 281 282 283
    std::vector<double>& getXDist();
    std::vector<double>& getBGxDist();
    std::vector<double>& getYDist();
    std::vector<double>& getBGyDist();
    std::vector<double>& getTOrZDist();
    std::vector<double>& getBGzDist();
    void eraseXDist();
    void eraseBGxDist();
    void eraseYDist();
    void eraseBGyDist();
    void eraseTOrZDist();
    void eraseBGzDist();

284
    void addDistributions();
285 286
    void applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
    void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
287
    void applyEmissModelNone(double &pz);
288
    void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
289 290 291 292 293 294 295
    void create(size_t &numberOfParticles, double massIneV);
    void calcPartPerDist(size_t numberOfParticles);
    void checkEmissionParameters();
    void checkIfEmitted();
    void checkParticleNumber(size_t &numberOfParticles);
    void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits);
    double converteVToBetaGamma(double valueIneV, double massIneV);
kraus's avatar
kraus committed
296
    size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324

    class BinomialBehaviorSplitter {
    public:
        virtual ~BinomialBehaviorSplitter()
        { }

        virtual double get(double rand) = 0;
    };

    class MDependentBehavior: public BinomialBehaviorSplitter {
    public:
        MDependentBehavior(const MDependentBehavior &rhs):
            ami_m(rhs.ami_m)
        {}

        MDependentBehavior(double a)
        { ami_m = 1.0 / a; }

        virtual double get(double rand);
    private:
        double ami_m;
    };

    class GaussianLikeBehavior: public BinomialBehaviorSplitter {
    public:
        virtual double get(double rand);
    };

325 326
    void createDistributionBinomial(size_t numberOfParticles, double massIneV);
    void createDistributionFlattop(size_t numberOfParticles, double massIneV);
327
    void createDistributionMultiGauss(size_t numberOfParticles, double massIneV);
328 329 330
    void createDistributionFromFile(size_t numberOfParticles, double massIneV);
    void createDistributionGauss(size_t numberOfParticles, double massIneV);
    void createMatchedGaussDistribution(size_t numberOfParticles, double massIneV);
331
    void sampleUniformDisk(gsl_qrng* quasiRandGen2D, double& x1, double& x2);
332 333 334 335 336 337 338 339 340
    void fillEBinHistogram();
    void fillParticleBins();
    size_t findEBin(double tOrZ);
    void generateAstraFlattopT(size_t numberOfParticles);
    void generateBinomial(size_t numberOfParticles);
    void generateFlattopLaserProfile(size_t numberOfParticles);
    void generateFlattopT(size_t numberOfParticles);
    void generateFlattopZ(size_t numberOfParticles);
    void generateGaussZ(size_t numberOfParticles);
341 342
    void generateMatchedGauss(const SigmaGenerator::matrix_t&,
                              size_t numberOfParticles, double massIneV);
343 344
    void generateLongFlattopT(size_t numberOfParticles);
    void generateTransverseGauss(size_t numberOfParticles);
frey_m's avatar
frey_m committed
345 346
    void initializeBeam(PartBunchBase<double, 3> *beam);
    void injectBeam(PartBunchBase<double, 3> *beam);
347 348 349
    void printDist(Inform &os, size_t numberOfParticles) const;
    void printDistBinomial(Inform &os) const;
    void printDistFlattop(Inform &os) const;
350
    void printDistMultiGauss(Inform &os) const;
351 352 353 354 355 356 357 358
    void printDistFromFile(Inform &os) const;
    void printDistGauss(Inform &os) const;
    void printDistMatchedGauss(Inform &os) const;
    void printEmissionModel(Inform &os) const;
    void printEmissionModelAstra(Inform &os) const;
    void printEmissionModelNone(Inform &os) const;
    void printEmissionModelNonEquil(Inform &os) const;
    void printEnergyBins(Inform &os) const;
359
    void adjustPhaseSpace(double massIneV);
360 361
    void reflectDistribution(size_t &numberOfParticles);
    void scaleDistCoordinates();
362 363
    /// Select and allocate gsl random number generator
    gsl_qrng* selectRandomGenerator(std::string, unsigned int dimension);
364 365 366
    void setAttributes();
    void setDistParametersBinomial(double massIneV);
    void setDistParametersFlattop(double massIneV);
367
    void setDistParametersMultiGauss(double massIneV);
368 369
    void setDistParametersGauss(double massIneV);
    void setEmissionTime(double &maxT, double &minT);
frey_m's avatar
frey_m committed
370 371 372
    void setupEmissionModel(PartBunchBase<double, 3> *beam);
    void setupEmissionModelAstra(PartBunchBase<double, 3> *beam);
    void setupEmissionModelNone(PartBunchBase<double, 3> *beam);
373 374
    void setupEmissionModelNonEquil();
    void setupEnergyBins(double maxTOrZ, double minTOrZ);
frey_m's avatar
frey_m committed
375
    void setupParticleBins(double massIneV, PartBunchBase<double, 3> *beam);
376 377 378 379
    void shiftDistCoordinates(double massIneV);
    void writeOutFileHeader();
    void writeOutFileEmission();
    void writeOutFileInjection();
380 381 382

    std::string distT_m;                 /// Distribution type. Declared as string
    DistrTypeT::DistrTypeT distrTypeT_m; /// and list type for switch statements.
gsell's avatar
gsell committed
383

384 385
    unsigned int numberOfDistributions_m;

386 387
    bool emitting_m;                     /// Distribution is an emitted, and is currently
                                         /// emitting, rather than an injected, beam.
gsell's avatar
gsell committed
388

389 390
    PartData particleRefData_m;          /// Reference data for particle type (charge,
                                         /// mass etc.)
gsell's avatar
gsell committed
391

392 393 394
    /// Vector of distributions to be added to this distribution.
    std::vector<Distribution *> addedDistributions_m;
    std::vector<size_t> particlesPerDist_m;
gsell's avatar
gsell committed
395

396 397
    /// Emission Model.
    EmissionModelT::EmissionModelT emissionModel_m;
gsell's avatar
gsell committed
398

399
    /// Emission parameters.
400
    double tEmission_m;
401 402 403 404
    static const double percentTEmission_m; /// Increase tEmission_m by twice this percentage
                                            /// to ensure that no particles fall on the leading edge of
                                            /// the first emission time step or the trailing edge of the last
                                            /// emission time step.
405 406 407 408 409 410 411 412 413 414 415 416
    double tBin_m;
    double currentEmissionTime_m;
    int currentEnergyBin_m;
    int currentSampleBin_m;
    int numberOfEnergyBins_m;       /// Number of energy bins the distribution
                                    /// is broken into. Used for an emitted beam.
    int numberOfSampleBins_m;       /// Number of samples to use per energy bin
                                    /// when emitting beam.
    PartBins *energyBins_m;         /// Distribution energy bins.
    gsl_histogram *energyBinHist_m; /// GSL histogram used to define energy bin
                                    /// structure.

417
    gsl_rng *randGen_m;             /// Random number generator
418 419 420

    // ASTRA and NONE photo emission model.
    double pTotThermal_m;           /// Total thermal momentum.
kraus's avatar
kraus committed
421
    Vector_t pmean_m;
422 423 424 425 426 427 428 429

    // NONEQUIL photo emission model.
    double cathodeWorkFunc_m;       /// Cathode material work function (eV).
    double laserEnergy_m;           /// Laser photon energy (eV).
    double cathodeFermiEnergy_m;    /// Cathode material Fermi energy (eV).
    double cathodeTemp_m;           /// Cathode temperature (K).
    double emitEnergyUpperLimit_m;  /// Upper limit on emission energy distribution (eV).

kraus's avatar
kraus committed
430
    std::vector<std::vector<double> > additionalRNs_m;
431

432 433 434
    size_t totalNumberParticles_m;
    size_t totalNumberEmittedParticles_m;

435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
    // Beam coordinate containers.
    std::vector<double> xDist_m;
    std::vector<double> pxDist_m;
    std::vector<double> yDist_m;
    std::vector<double> pyDist_m;
    std::vector<double> tOrZDist_m;
    std::vector<double> pzDist_m;

    // Initial coordinates for file write.
    std::vector<double> xWrite_m;
    std::vector<double> pxWrite_m;
    std::vector<double> yWrite_m;
    std::vector<double> pyWrite_m;
    std::vector<double> tOrZWrite_m;
    std::vector<double> pzWrite_m;
    std::vector<size_t> binWrite_m;

452
    // for compatibility reasons
453
    double avrgpz_m;
454

455 456 457 458 459 460 461 462 463 464
    //Distribution parameters.
    InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits_m;
    double sigmaTRise_m;
    double sigmaTFall_m;
    double tPulseLengthFWHM_m;
    Vector_t sigmaR_m;
    Vector_t sigmaP_m;
    Vector_t cutoffR_m;
    Vector_t cutoffP_m;
    Vector_t mBinomial_m;
465
    SymTenzor<double, 6> correlationMatrix_m;
466

467 468 469
    // Parameters multiGauss distribution.
    double sepPeaks_m;
    unsigned nPeaks_m;
470

471 472 473 474 475 476
    // Laser profile.
    std::string laserProfileFileName_m;
    std::string laserImageName_m;
    double laserIntensityCut_m;
    LaserProfile *laserProfile_m;

477 478 479
    // AAA This is for the matched distribution
    double I_m;
    double E_m;
480

gsell's avatar
gsell committed
481 482 483 484 485 486 487 488 489 490 491 492
    /// time binned distribution with thermal energy
    double tRise_m;
    double tFall_m;
    double sigmaRise_m;
    double sigmaFall_m;
    double cutoff_m;
};

inline Inform &operator<<(Inform &os, const Distribution &d) {
    return d.printInfo(os);
}

kraus's avatar
kraus committed
493 494 495 496 497
inline
Vector_t Distribution::get_pmean() const {
    return pmean_m;
}

498 499 500 501 502
inline
DistrTypeT::DistrTypeT Distribution::getType() const {
    return distrTypeT_m;
}

503 504
inline
double Distribution::getPercentageEmitted() const {
505 506 507
    return (double)totalNumberEmittedParticles_m / (double)totalNumberParticles_m;
}

508 509 510 511 512
inline
std::string Distribution::getTypeofDistribution() {
    return (std::string) Attributes::getString(itsAttr[Attrib::Distribution::TYPE]);
}

kraus's avatar
kraus committed
513
#endif // OPAL_Distribution_HH