Distribution.h 17.4 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
    size_t getNumOfLocalParticlesToCreate(size_t n);

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

126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
    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
147

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

kraus's avatar
kraus committed
152
    Vector_t get_pmean() const;
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 181 182
    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
183 184 185

    Inform &printInfo(Inform &os) const;

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

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

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

205
    void setNumberOfDistributions(unsigned int n) { numberOfDistributions_m = n; }
206 207

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

214
    Distribution(const std::string &name, Distribution *parent);
215

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

220 221 222

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

223
    void addDistributions();
224 225
    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);
226
    void applyEmissModelNone(double &pz);
227
    void applyEmissModelNonEquil(double eZ, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
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 286 287
    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();
288 289 290 291 292

    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
293

294 295
    unsigned int numberOfDistributions_m;

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

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

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

306 307 308
    /// 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
309

310 311
    /// Emission Model.
    EmissionModelT::EmissionModelT emissionModel_m;
gsell's avatar
gsell committed
312

313 314 315 316 317 318 319 320 321 322 323 324 325 326
    /// 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.

327
    gsl_rng *randGen_m;             /// Random number generator
328 329 330

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

    // 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
340
    std::vector<std::vector<double> > additionalRNs_m;
341

342 343 344
    size_t totalNumberParticles_m;
    size_t totalNumberEmittedParticles_m;

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

362
    // for compatibility reasons
363
    double avrgpz_m;
364 365 366



367 368 369 370 371 372 373 374 375 376
    //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;
377
    SymTenzor<double, 6> correlationMatrix_m;
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 418 419

    // 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.


420 421 422 423 424 425 426
    // 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
427 428
    double bg_m;                      /// beta gamma
    double M_m;                       /// mass in terms of proton mass
429 430
    std::string bfieldfn_m;           /// only temporarly

431 432


433 434


435 436
    // Some legacy members that need to be cleaned up.

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

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

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

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

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

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

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

477
#endif // OPAL_Distribution_HH