Distribution.h 17.3 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
#ifndef OPAL_Distribution_HH
#define OPAL_Distribution_HH

// ------------------------------------------------------------------------
// $RCSfile: Distribution.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Distribution
//
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:33:44 $
// $Author: Andreas Adelmann $
//
// ------------------------------------------------------------------------
#include <iosfwd>
#include <fstream>
22
#include <forward_list>
gsell's avatar
gsell committed
23 24 25
#include <string>
#include <map>

26 27
#include "AbstractObjects/Definition.h"
#include "Algorithms/PartData.h"
28

29
#include "Algorithms/Vektor.h"
30
#include "Beamlines/Beamline.h"
gsell's avatar
gsell committed
31

32 33 34 35 36 37 38 39
#include "Ippl.h"

#include "H5hut.h"

#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_qrng.h>

40 41 42
#ifdef WITH_UNIT_TESTS
#include <gtest/gtest_prod.h>
#endif
43

44
class Beam;
45 46 47 48 49
class PartBunch;
class PartBins;
class EnvelopeBunch;
class BoundaryGeometry;
class LaserProfile;
50
class H5PartWrapper;
gsell's avatar
gsell committed
51

52 53 54 55 56 57 58 59 60 61
namespace DistrTypeT
{
    enum DistrTypeT {NODIST,
                    FROMFILE,
                    GAUSS,
                    BINOMIAL,
                    FLATTOP,
                    SURFACEEMISSION,
                    SURFACERANDCREATE,
                    GUNGAUSSFLATTOPTH,
62 63
	            ASTRAFLATTOPTH,
		    MATCHEDGAUSS
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
                    };
}

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

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

/*
 * Class Distribution
 *
 * Defines the initial beam that is injected or emitted into the simulation.
 */
gsell's avatar
gsell committed
87 88 89 90 91 92 93 94 95

class Distribution: public Definition {

public:

    Distribution();
    virtual ~Distribution();

    virtual bool canReplaceBy(Object *object);
96
    virtual Distribution *clone(const std::string &name);
gsell's avatar
gsell committed
97 98 99
    virtual void execute();
    virtual void update();

100 101
    void createBoundaryGeometry(PartBunch &p, BoundaryGeometry &bg);
    void createOpalCycl(PartBunch &beam,
102
                        size_t numberOfParticles,
103
			double current, const Beamline &bl,
104
                        bool scan);
105
    void createOpalE(Beam *beam,
106 107 108 109
                     std::vector<Distribution *> addedDistributions,
                     EnvelopeBunch *envelopeBunch,
                     double distCenter,
                     double Bz0);
110
    void createOpalT(PartBunch &beam,
111 112 113
                     std::vector<Distribution *> addedDistributions,
                     size_t &numberOfParticles,
                     bool scan);
114 115 116 117
    void createOpalT(PartBunch &beam, size_t &numberOfParticles, bool scan);
    void createPriPart(PartBunch *beam, BoundaryGeometry &bg);
    void doRestartOpalT(PartBunch &p, size_t Np, int restartStep, H5PartWrapper *h5wrapper);
    void doRestartOpalCycl(PartBunch &p, size_t Np, int restartStep,
118
                        const int specifiedNumBunch, H5PartWrapper *h5wrapper);
119 120 121
    void doRestartOpalE(EnvelopeBunch &p, size_t Np, int restartStep, H5PartWrapper *h5wrapper);
    size_t emitParticles(PartBunch &beam, double eZ);
    double getPercentageEmitted() const;
122
    static Distribution *find(const std::string &name);
gsell's avatar
gsell committed
123

124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
    void eraseXDist();
    void eraseBGxDist();
    void eraseYDist();
    void eraseBGyDist();
    void eraseTOrZDist();
    void eraseBGzDist();
    bool getIfDistEmitting();
    int getLastEmittedEnergyBin();
    double getMaxTOrZ();
    double getMinTOrZ();
    size_t getNumberOfEmissionSteps();
    int getNumberOfEnergyBins();
    double getEmissionDeltaT();
    double getEnergyBinDeltaT();
    double getWeight();
    std::vector<double>& getXDist();
    std::vector<double>& getBGxDist();
    std::vector<double>& getYDist();
    std::vector<double>& getBGyDist();
    std::vector<double>& getTOrZDist();
    std::vector<double>& getBGzDist();
gsell's avatar
gsell committed
145

146
    /// Return the embedded CLASSIC PartData.
147 148
    const PartData &getReference() const;
    double getTEmission();
149

kraus's avatar
kraus committed
150
    Vector_t get_pmean() const;
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 176 177 178 179 180
    double getEkin() const;
    double getLaserEnergy() const;
    double getWorkFunctionRf() const;

    size_t getNumberOfDarkCurrentParticles();
    double getDarkCurrentParticlesInwardMargin();
    double getEInitThreshold();
    double getWorkFunction();
    double getFieldEnhancement();
    size_t getMaxFNemissionPartPerTri();
    double getFieldFNThreshold();
    double getFNParameterA();
    double getFNParameterB();
    double getFNParameterY();
    double getFNParameterVYZero();
    double getFNParameterVYSecond();
    int    getSecondaryEmissionFlag();
    bool   getEmissionMode() ;

    std::string getTypeofDistribution();

    double getvSeyZero();//return sey_0 in Vaughan's model
    double getvEZero();//return the energy related to sey_0 in Vaughan's model
    double getvSeyMax();//return sey max in Vaughan's model
    double getvEmax();//return Emax in Vaughan's model
    double getvKenergy();//return fitting parameter denotes the roughness of surface for impact energy in Vaughan's model
    double getvKtheta();//return fitting parameter denotes the roughness of surface for impact angle in Vaughan's model
    double getvVThermal();//return thermal velocity of Maxwellian distribution of secondaries in Vaughan's model
    double getVw();//return velocity scalar for parallel plate benchmark;
    int getSurfMaterial();//material type for Furman-Pivi's model 0 for copper, 1 for stainless steel
gsell's avatar
gsell committed
181 182 183

    Inform &printInfo(Inform &os) const;

184
    bool Rebin();
185 186 187 188
    void setDistToEmitted(bool emitted);
    void setDistType();
    void shiftBeam(double &maxTOrZ, double &minTOrZ);
    double getEmissionTimeShift() const;
gsell's avatar
gsell committed
189

190
    // in case if OPAL-cycl in restart mode
191
    double GetBeGa() {return bega_m;}
192
    double GetPr() {return referencePr_m;}
193 194
    double GetPt() {return referencePt_m;}
    double GetPz() {return referencePz_m;}
195 196
    double GetR() {return referenceR_m;}
    double GetTheta() {return referenceTheta_m;}
197
    double GetZ() {return referenceZ_m;}
gsell's avatar
gsell committed
198

199 200 201 202
    double GetPhi() {return phi_m;}
    double GetPsi() {return psi_m;}
    bool GetPreviousH5Local() {return previousH5Local_m;}

203
    void setNumberOfDistributions(unsigned int n) { numberOfDistributions_m = n; }
204 205

    DistrTypeT::DistrTypeT getType() const;
206
private:
207
#ifdef WITH_UNIT_TESTS
208 209
    FRIEND_TEST(GaussTest, FullSigmaTest1);
    FRIEND_TEST(GaussTest, FullSigmaTest2);
210
#endif
211

212
    Distribution(const std::string &name, Distribution *parent);
213

gsell's avatar
gsell committed
214 215 216 217
    // Not implemented.
    Distribution(const Distribution &);
    void operator=(const Distribution &);

218 219 220

    //    void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);

221
    void addDistributions();
222 223
    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);
224
    void applyEmissModelNone(double &pz);
225
    void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
226 227 228 229 230 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 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
    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 convertBetaGammaToeV(double valueInbega, double mass);
    double converteVToBetaGamma(double valueIneV, double massIneV);
    double convertMeVPerCToBetaGamma(double valueInMeVPerC, double massIneV);
    void createDistributionBinomial(size_t numberOfParticles, double massIneV);
    void createDistributionFlattop(size_t numberOfParticles, double massIneV);
    void createDistributionFromFile(size_t numberOfParticles, double massIneV);
    void createDistributionGauss(size_t numberOfParticles, double massIneV);
    void createMatchedGaussDistribution(size_t numberOfParticles, double massIneV);
    void destroyBeam(PartBunch &beam);
    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);
    void generateLongFlattopT(size_t numberOfParticles);
    void generateTransverseGauss(size_t numberOfParticles);
    void initializeBeam(PartBunch &beam);
    void injectBeam(PartBunch &beam);
    void printDist(Inform &os, size_t numberOfParticles) const;
    void printDistBinomial(Inform &os) const;
    void printDistFlattop(Inform &os) const;
    void printDistFromFile(Inform &os) const;
    void printDistGauss(Inform &os) const;
    void printDistMatchedGauss(Inform &os) const;
    void printDistSurfEmission(Inform &os) const;
    void printDistSurfAndCreate(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;
    void reflectDistribution(size_t &numberOfParticles);
    void scaleDistCoordinates();
    void setAttributes();
    void setDistParametersBinomial(double massIneV);
    void setDistParametersFlattop(double massIneV);
    void setDistParametersGauss(double massIneV);
    void setEmissionTime(double &maxT, double &minT);
    void setFieldEmissionParameters();
    void setupEmissionModel(PartBunch &beam);
    void setupEmissionModelAstra(PartBunch &beam);
    void setupEmissionModelNone(PartBunch &beam);
    void setupEmissionModelNonEquil();
    void setupEnergyBins(double maxTOrZ, double minTOrZ);
    void setupParticleBins(double massIneV, PartBunch &beam);
    void shiftDistCoordinates(double massIneV);
    void writeOutFileHeader();
    void writeOutFileEmission();
    void writeOutFileInjection();
    void writeToFile();
286 287 288 289 290

    std::string distT_m;                 /// Distribution type. Declared as string
    DistrTypeT::DistrTypeT distrTypeT_m; /// and list type for switch statements.
    std::ofstream os_m;                  /// Output file to write distribution.
    void writeToFileCycl(PartBunch &beam, size_t Np);
gsell's avatar
gsell committed
291

292 293
    unsigned int numberOfDistributions_m;

294 295
    bool emitting_m;                     /// Distribution is an emitted, and is currently
                                         /// emitting, rather than an injected, beam.
gsell's avatar
gsell committed
296

297 298 299
    bool scan_m;                         /// If we are doing a scan, we need to
                                         /// destroy existing particles before
                                         /// each run.
gsell's avatar
gsell committed
300

301 302
    PartData particleRefData_m;          /// Reference data for particle type (charge,
                                         /// mass etc.)
gsell's avatar
gsell committed
303

304 305 306
    /// 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
307

308 309
    /// Emission Model.
    EmissionModelT::EmissionModelT emissionModel_m;
gsell's avatar
gsell committed
310

311 312 313 314 315 316 317 318 319 320 321 322 323 324
    /// Emission parameters.
    double tEmission_m;
    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.

325
    gsl_rng *randGen_m;             /// Random number generator
326 327 328

    // ASTRA and NONE photo emission model.
    double pTotThermal_m;           /// Total thermal momentum.
kraus's avatar
kraus committed
329
    Vector_t pmean_m;
330 331 332 333 334 335 336 337

    // 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
338
    std::vector<std::vector<double> > additionalRNs_m;
339

340 341 342
    size_t totalNumberParticles_m;
    size_t totalNumberEmittedParticles_m;

343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
    // 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;

360
    // for compatibility reasons
361
    double avrgpz_m;
362 363 364



365 366 367 368 369 370 371 372 373 374
    //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;
375
    SymTenzor<double, 6> correlationMatrix_m;
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417

    // Laser profile.
    std::string laserProfileFileName_m;
    std::string laserImageName_m;
    double laserIntensityCut_m;
    LaserProfile *laserProfile_m;

    /*
     * Dark current calculation parameters.
     */
    size_t darkCurrentParts_m;      /// Number of dark current particles.
    double darkInwardMargin_m;      /// Dark current particle initialization position.
                                    /// Inward along the triangle normal, positive.
                                    /// Inside the geometry.
    double eInitThreshold_m;        /// Field threshold (MV/m) beyond which particles
                                    /// will be initialized.
    double workFunction_m;          /// Work function of surface material (eV).
    double fieldEnhancement_m;      /// Field enhancement factor beta for Fowler-
                                    /// Nordheim emission.
    double fieldThrFN_m;            /// Field threshold for Fowler-Nordheim
                                    /// emission (MV/m).
    size_t maxFN_m;                 /// Max. number of electrons emitted from a
                                    /// single triangle for Fowler-Nordheim emission.
    double paraFNA_m;               /// Empirical constant A in Fowler-Nordheim
                                    /// emission model.
    double paraFNB_m;               /// Empirical constant B in Fowler-Nordheim
                                    /// emission model.
    double paraFNY_m;               /// Constant for image charge effect parameter y(E)
                                    /// in Fowler-Nordheim emission model.
    double paraFNVYSe_m;            /// Second order constant for v(y) function in
                                    /// Fowler-Nordheim emission model.
    double paraFNVYZe_m;            /// Zero order constant for v(y) function in
                                    /// Fowler-Nordheim emission model.
    int    secondaryFlag_m;         /// Select the secondary model type:
                                    ///     0           ==> no secondary emission.
                                    ///     1           ==> Furman-Pivi
                                    ///     2 or larger ==> Vaughan's model
    double ppVw_m;                  /// Velocity scalar for parallel plate benchmark.
    double vVThermal_m;             /// Thermal velocity of Maxwellian distribution
                                    /// of secondaries in Vaughan's model.


418 419 420 421 422 423 424
    // AAA This is for the matched distribution
    double ex_m;
    double ey_m;
    double et_m;

    double I_m;
    double E_m;
adelmann's avatar
adelmann committed
425 426
    double bg_m;                      /// beta gamma
    double M_m;                       /// mass in terms of proton mass
427 428 429 430 431 432
    std::string bfieldfn_m;           /// only temporarly





433 434
    // Some legacy members that need to be cleaned up.

gsell's avatar
gsell committed
435 436 437 438 439 440 441
    /// time binned distribution with thermal energy
    double tRise_m;
    double tFall_m;
    double sigmaRise_m;
    double sigmaFall_m;
    double cutoff_m;

442
    // Cyclotron stuff
443
    double referencePr_m;
444
    double referencePt_m;
adelmann's avatar
adelmann committed
445
    double referencePz_m;
446 447
    double referenceR_m;
    double referenceTheta_m;
adelmann's avatar
adelmann committed
448
    double referenceZ_m;
449
    double bega_m;
450

451 452 453 454
    // Cyclotron for restart in local mode
    double phi_m;
    double psi_m;
    bool previousH5Local_m;
gsell's avatar
gsell committed
455 456 457 458 459 460
};

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

kraus's avatar
kraus committed
461 462 463 464 465
inline
Vector_t Distribution::get_pmean() const {
    return pmean_m;
}

466 467 468 469 470 471 472 473 474
inline
DistrTypeT::DistrTypeT Distribution::getType() const {
    return distrTypeT_m;
}

inline double Distribution::getPercentageEmitted() const {
    return (double)totalNumberEmittedParticles_m / (double)totalNumberParticles_m;
}

475
#endif // OPAL_Distribution_HH