PartBunch.h 29.5 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 22
#ifndef OPAL_PartBunch_HH
#define OPAL_PartBunch_HH

// ------------------------------------------------------------------------
// $RCSfile: PartBunch.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class PartBunch
//
// ------------------------------------------------------------------------
// Class category: Algorithms
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:33 $
// $Author: Andreas Adelmann  and Co. $
//
// ------------------------------------------------------------------------

23
#include "Ippl.h"
kraus's avatar
kraus committed
24
#include "Algorithms/PBunchDefs.h"
gsell's avatar
gsell committed
25 26 27
#include "Algorithms/Particle.h"
#include "FixedAlgebra/FMatrix.h"
#include "FixedAlgebra/FVector.h"
28
#include "Algorithms/PartBins.h"
gsell's avatar
gsell committed
29 30
#include "Algorithms/PartBinsCyc.h"
#include "Algorithms/PartData.h"
31 32
#include "Utilities/SwitcherError.h"
#include "Physics/Physics.h"
gsell's avatar
gsell committed
33 34 35 36 37 38 39

#include <iosfwd>
#include <vector>


class Distribution;
class LossDataSink;
40
class FieldSolver;
gsell's avatar
gsell committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
class ListElem;

template <class T, int, int> class FMatrix;
template <class T, int> class FVector;

// Class PartBunch.
// ------------------------------------------------------------------------
/// Particle Bunch.
//  A representation of a particle bunch as a vector of particles.

// class PartBunch: public std::vector<Particle>, public ParticleBase< ParticleSpatialLayout<double, 3> > {
class PartBunch: public ParticleBase< ParticleSpatialLayout<double, 3> > {

public:
    /// Default constructor.
    //  Construct empty bunch.
    PartBunch(const PartData *ref);

    /// Conversion.
    PartBunch(const std::vector<Particle> &, const PartData *ref);

    PartBunch(const PartBunch &);
    ~PartBunch();

65 66 67 68 69 70 71 72
    bool GetIfBeamEmitting();
    int GetLastEmittedEnergyBin();
    size_t GetNumberOfEmissionSteps();
    int GetNumberOfEnergyBins();
    void Rebin();
    void SetEnergyBins(int numberOfEnergyBins);
    bool WeHaveEnergyBins();

73 74 75
    // helpers to store and restore a PartBunch
    void stash();
    void pop();
76
    Vector_t getStashIniP() const;
77 78

    enum UnitState_t { units = 0, unitless = 1 };
79
    UnitState_t getUnitState() const;
80 81

    //FIXME: unify methods, use convention that all particles have own dt
82
    void switchToUnitlessPositions(bool use_dt_per_particle = false);
83 84

    //FIXME: unify methods, use convention that all particles have own dt
85
    void switchOffUnitlessPositions(bool use_dt_per_particle = false);
86

gsell's avatar
gsell committed
87 88 89 90 91 92 93 94 95 96
    void makHistograms();

    /** \brief After each Schottky scan we delete
    all the particles.

    */
    void cleanUpParticles();

    void resetIfScan();

97
    double getRho(NDIndex<3> e);
gsell's avatar
gsell committed
98

99
    double getRho(int x, int y, int z);
gsell's avatar
gsell committed
100

101
    bool itIsMyTurn(int *n);
gsell's avatar
gsell committed
102 103 104

    bool hasZeroNLP();

105
    void do_binaryRepart();
gsell's avatar
gsell committed
106

kraus's avatar
kraus committed
107

108

gsell's avatar
gsell committed
109
    /// per default the MT value of the field solver is used
110
    void set_nBinsLineDensity(int n);
gsell's avatar
gsell committed
111 112 113 114 115

    void calcLineDensity();
    void fillArray(double *lineDensity, const std::list<ListElem> &l);
    void getLineDensity(std::vector<double> &lineDensity);

116 117 118 119
    void setDistribution(Distribution *d,
                         std::vector<Distribution *> addedDistributions,
                         size_t &np,
                         bool scan);
gsell's avatar
gsell committed
120 121


122 123
    void setGridIsFixed();
    bool isGridFixed();
gsell's avatar
gsell committed
124 125 126 127 128 129 130

    /*

    Energy bins related functions

    */

131 132 133
    void setTEmission(double t);
    double getTEmission();
    bool doEmission();
gsell's avatar
gsell committed
134

135
    bool weHaveBins() const;
gsell's avatar
gsell committed
136

137
    double getRebinEnergy();
gsell's avatar
gsell committed
138

139
    void weHaveNOBins();
gsell's avatar
gsell committed
140 141 142 143 144 145 146 147

    void setPBins(PartBins *pbin);
    void setPBins(PartBinsCyc *pbin);

    /** \brief Emit particles in the given bin
        i.e. copy the particles from the bin structure into the
        particle container
    */
148
    size_t EmitParticles(double eZ);
gsell's avatar
gsell committed
149 150 151 152 153

    double calcTimeDelay(const double &jifactor);
    void moveBunchToCathode(double &t);
    void printBinHist();

154
    void rebin();
gsell's avatar
gsell committed
155

156
    int getNumBins();
gsell's avatar
gsell committed
157

158
    int getLastemittedBin();
gsell's avatar
gsell committed
159

160
    void updatePartInBin(size_t countLost[]);
gsell's avatar
gsell committed
161 162 163 164 165 166 167 168 169 170 171 172 173

    /** \brief we need this because only node 0 is emitting */
    void updateBinStructure();

    /** \brief report if any particle has been lost */
    void reportParticleLoss();

    /** \brief Compute the gammas of all bins */
    void calcGammas();

    void calcGammas_cycl();

    /** \brief Get gamma of one bin */
174
    double getBinGamma(int bin);
gsell's avatar
gsell committed
175 176

    /** \brief Set the charge of one bin to the value of q and all other to zero */
177
    void setBinCharge(int bin, double q);
gsell's avatar
gsell committed
178 179

    /** \brief Set the charge of all other the ones in bin to zero */
180
    void setBinCharge(int bin);
gsell's avatar
gsell committed
181 182 183 184 185

    /** \brief gets back the maximum dE of all the bins */
    double getMaxdEBins();

    /** \brief calculates back the max/min of the efield on the grid */
186
    std::pair<Vector_t, Vector_t> getEExtrema();
gsell's avatar
gsell committed
187 188 189 190 191 192 193

    /*

    Mesh and Field Layout related functions

    */

194
    const Mesh_t &getMesh() const;
gsell's avatar
gsell committed
195

196
    Mesh_t &getMesh();
gsell's avatar
gsell committed
197

198
    FieldLayout_t &getFieldLayout();
gsell's avatar
gsell committed
199

200
    void setBCAllOpen();
gsell's avatar
gsell committed
201

202 203
    void setBCForDCBeam();

gsell's avatar
gsell committed
204 205 206 207
    void boundp();

    /** delete particles which are too far away from the center of beam*/
    void boundp_destroy();
208

gsell's avatar
gsell committed
209 210 211 212 213 214 215 216 217 218
    /** This is only temporary in order to get the collimator and pepperpot workinh */
    size_t boundp_destroyT();

    /* Definiere virtuelle Funktionen, um die Koordinaten auszulesen
     *     */

    virtual double getPx(int i);
    virtual double getPy(int i);
    virtual double getPz(int i);

219 220
    virtual double getPx0(int i);
    virtual double getPy0(int i);
gsell's avatar
gsell committed
221 222 223 224 225 226 227 228

    virtual double getX(int i);
    virtual double getY(int i);
    virtual double getZ(int i);

    virtual double getX0(int i);
    virtual double getY0(int i);

229
    virtual void setZ(int i, double zcoo);
gsell's avatar
gsell committed
230

231
    virtual void setKR(Vector_t value, int i);
gsell's avatar
gsell committed
232

233
    virtual void setKT(Vector_t value, int i);
gsell's avatar
gsell committed
234

235
    virtual void BetOut(FILE *dat, FILE *sli);
gsell's avatar
gsell committed
236

237
    virtual void plotR();
gsell's avatar
gsell committed
238

239
    void get_bounds(Vector_t &rmin, Vector_t &rmax);
gsell's avatar
gsell committed
240 241 242 243 244 245

    /*
      Compatibility function push_back

    */

246
    void push_back(Particle p);
gsell's avatar
gsell committed
247

248
    void set_part(FVector<double, 6> z, int ii);
gsell's avatar
gsell committed
249

250
    void set_part(Particle p, int ii);
gsell's avatar
gsell committed
251

252
    Particle get_part(int ii);
gsell's avatar
gsell committed
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269

    /// Return bunch distribution.
    //  Return the bunch centroid in [b]centroid[/b],
    //  and the second moments in [b]moment[/b].
    void beamEllipsoid(FVector<double, 6>   &centroid,
                       FMatrix<double, 6, 6> &moment);

    /// Return maximum amplitudes.
    //  The matrix [b]D[/b] is used to normalise the first two modes.
    //  The maximum normalised amplitudes for these modes are stored
    //  in [b]axmax[/b] and [b]aymax[/b].
    void maximumAmplitudes(const FMatrix<double, 6, 6> &D,
                           double &axmax, double &aymax);

    Inform &print(Inform &os);


270 271
    void   setdT(double dt);
    double getdT() const;
gsell's avatar
gsell committed
272

273 274
    void   setT(double t);
    double getT() const;
gsell's avatar
gsell committed
275 276 277
    void   computeSelfFields();

    /** /brief used for self fields with binned distribution */
278
    void   computeSelfFields(int b);
gsell's avatar
gsell committed
279

280 281
    //void computeSelfFields_cycl(double gamma);
    //void computeSelfFields_cycl(int b);
gsell's avatar
gsell committed
282

kraus's avatar
kraus committed
283
    // Replaced computeSelfFields_cycl() with versions that have meanR and the quaternion of the
284
    // rotation of the particle bunch in order to take into account the rotation
285
    // when finding the boundary conditions for the fieldsolver. -DW
kraus's avatar
kraus committed
286
    void computeSelfFields_cycl(double gamma, Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
287 288
				              Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));

kraus's avatar
kraus committed
289
    void computeSelfFields_cycl(int b,        Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
290
                                              Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));
gsell's avatar
gsell committed
291

292 293
    void resetInterpolationCache(bool clearCache = false);

gsell's avatar
gsell committed
294 295 296 297 298 299 300 301 302 303 304 305
    void calcWeightedAverages(Vector_t &CentroidPosition, Vector_t &CentroidMomentum) const;

    /** EulerAngle[0] = rotation about the y-axis, EulerAngle[1] = rotation about x-axis
     *  EulerAngle[2] = rotation about the y'-axis */

    void rotateAbout(const Vector_t &Center, const Vector_t &EulerAngles);

    void moveBy(const Vector_t &Center);

    void ResetLocalCoordinateSystem(const int &i, const Vector_t &Orientation, const double &origin);

    // OLD stuff should go away
306 307
    bool isZPeriodic() const;
    double getGaBeLa() const;
gsell's avatar
gsell committed
308 309 310 311 312 313 314

    /**
     * get the spos of the primary beam
     *
     * @param none
     *
     */
315
    double get_sPos();
gsell's avatar
gsell committed
316

317 318 319 320 321 322
    /// Get average z position from local "lab frame" coordinates, X.
    double getZPos();

    /// Get bounds of local "lab frame" coordinates, X.
    void getXBounds(Vector_t &xMin, Vector_t &xMax);

323 324
    double   get_phase() const;
    double   get_gamma() const;
gsell's avatar
gsell committed
325

326 327 328 329 330 331 332 333
    double get_meanEnergy() const;
    //  double* get_energy() {return energy_m; }
    Vector_t get_origin() const;
    Vector_t get_maxExtend() const;
    Vector_t get_centroid() const;
    Vector_t get_rrms() const;
    Vector_t get_rmean() const;
    Vector_t get_prms() const;
adelmann's avatar
adelmann committed
334
    Vector_t get_rprrms() const;
335 336 337 338
    Vector_t get_pmean() const;
    Vector_t get_emit() const;
    Vector_t get_norm_emit() const;
    Vector_t get_hr() const;
gsell's avatar
gsell committed
339

340 341
    double get_Dx() const;
    double get_Dy() const;
gsell's avatar
gsell committed
342

343 344
    double get_DDx() const;
    double get_DDy() const;
gsell's avatar
gsell committed
345

346 347
    void set_meshEnlargement(double dh);
    double get_meshEnlargement() const;
gsell's avatar
gsell committed
348 349

    void gatherLoadBalanceStatistics();
350
    size_t getLoadBalance(int p) const;
gsell's avatar
gsell committed
351

352
    void get_PBounds(Vector_t &min, Vector_t &max) const;
gsell's avatar
gsell committed
353 354

    void calcBeamParameters();
355
    void calcBeamParametersLight();   // used in autophase and avoides communication
gsell's avatar
gsell committed
356 357 358
    void calcBeamParametersInitial(); // Calculate initial beam parameters before emission.
    void calcBeamParameters_cycl();

359 360
    double getCouplingConstant() const;
    void setCouplingConstant(double c);
gsell's avatar
gsell committed
361 362

    // set the charge per simulation particle
363
    void setCharge(double q);
gsell's avatar
gsell committed
364
    // set the charge per simulation particle when total particle number equals 0
365
    void setChargeZeroPart(double q);
gsell's avatar
gsell committed
366 367

    // set the mass per simulation particle
368
    void setMass(double mass);
gsell's avatar
gsell committed
369 370 371 372 373 374 375 376 377 378 379 380


    /// \brief Need Ek for the Schottky effect calculation (eV)
    double getEkin() const;

    /// Need the work function for the Schottky effect calculation (eV)
    double getWorkFunctionRf() const;

    /// Need the laser energy for the Schottky effect calculation (eV)
    double getLaserEnergy() const;

    /// get the total charge per simulation particle
381
    double getCharge() const;
gsell's avatar
gsell committed
382 383

    /// get the macro particle charge
384
    double getChargePerParticle() const;
gsell's avatar
gsell committed
385 386 387 388 389

    void setSolver(FieldSolver *fs);

    bool hasFieldSolver();

390 391
    void setLPath(double s);
    double getLPath() const;
gsell's avatar
gsell committed
392

393 394
    void setStepsPerTurn(int n);
    int getStepsPerTurn() const;
gsell's avatar
gsell committed
395

396 397 398
    /// step in multiple TRACK commands
    inline void setGlobalTrackStep(long long n) {globalTrackStep_m = n;}
    inline long long getGlobalTrackStep() const {return globalTrackStep_m;}
gsell's avatar
gsell committed
399

400 401
    /// step in a TRACK command
    inline void setLocalTrackStep(long long n) {localTrackStep_m = n;}
402
    inline void incTrackSteps() {globalTrackStep_m++; localTrackStep_m++;}
403 404 405 406
    inline long long getLocalTrackStep() const {return localTrackStep_m;}

    inline void setNumBunch(int n);
    inline int getNumBunch() const;
gsell's avatar
gsell committed
407

408 409
    inline Vector_t getGlobalMeanR() {return globalMeanR_m;}
    inline Vektor<double, 4> getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
410

411 412
    void setSteptoLastInj(int n);
    int getSteptoLastInj();
gsell's avatar
gsell committed
413 414 415 416 417 418 419 420 421 422 423 424

    /// calculate average angle of longitudinal direction of bins
    double calcMeanPhi();

    size_t getNumPartInBin(int BinID) const;

    /// reset Bin[] for each particle
    bool resetPartBinID();
    /// reset Bin[] for each particle according to the method given in paper PAST-AB(064402) by  G. Fubiani et al.
    bool resetPartBinID2(const double eta);


425 426 427 428
    double getQ() const;
    double getM() const;
    double getP() const;
    double getE() const;
gsell's avatar
gsell committed
429

430 431
    void resetQ(double q);
    void resetM(double m);
gsell's avatar
gsell committed
432

433 434 435 436 437 438 439
    double getdE();
    double getBeta() const;
    double getGamma() const;
    virtual double getGamma(int i);
    virtual double getBeta(int i);
    virtual void actT();
    const PartData *getReference() const;
gsell's avatar
gsell committed
440 441

    double getTBin();
442
    double GetEmissionDeltaT();
443
    bool isDcBeam();
444 445

    void iterateEmittedBin(int binNumber);
gsell's avatar
gsell committed
446

447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
    // Particle container attributes
    ParticleAttrib< Vector_t > X;      // local 'lab frame' coordinates;
    ParticleAttrib< Vector_t > P;      // particle momentum //  ParticleSpatialLayout<double, 3>::ParticlePos_t P;
    ParticleAttrib< double >   Q;      // charge per simulation particle, unit: C.
    ParticleAttrib< double >   M;      // mass per simulation particle, for multi-species particle tracking, unit:GeV/c^2.
    ParticleAttrib< Vector_t > Ef;     // e field vector
    ParticleAttrib< Vector_t > Eftmp;  // e field vector for gun simulations

    ParticleAttrib< Vector_t > Bf;   // b field vector
    ParticleAttrib< int >      Bin;   // holds the bin in which the particle is in, if zero particle is marked for deletion
    ParticleAttrib< double >   dt;   // holds the dt timestep for particle
    ParticleAttrib< long >     LastSection; // last em-field section


    ParticleAttrib< short >    PType; // we can distinguish dark current particles from primary particle
    ParticleAttrib< int >      TriID; // holds the ID of triangle that the particle hit. Only for BoundaryGeometry case.

    Vector_t RefPart_R;
    Vector_t RefPart_P;

    /// scalar potential
    Field_t rho_m;

    /// scalar fields for projecttion i.e. line densities
    Field_t tmpFieldZ_m;

    /// vector field on the grid
    VField_t  eg_m;

    /// avoid calls to Ippl::myNode()
    int myNode_m;

    /// avoid calls to Ippl::getNodes()
    int nodes_m;

    /// if the grid does not have to adapt
    bool fixed_grid;

    // The structure for particle binning
    PartBins *pbin_m;

    std::unique_ptr<LossDataSink> lossDs_m;

    // save particles in case of one core
    std::unique_ptr<Inform> pmsg_m;
    std::unique_ptr<std::ofstream> f_stream;

    /// timer for IC, can not be in Distribution.h
    IpplTimings::TimerRef distrReload_m;
    IpplTimings::TimerRef distrCreate_m;

498 499 500
    // For AMTS integrator in OPAL-T
    double dtScInit_m, deltaTau_m;

adelmann's avatar
adelmann committed
501 502 503 504
    /// for the Courant Shnider parameters
    Vector_t csBeta_m;
    Vector_t csAlpha_m;

505 506 507 508 509 510 511 512 513
protected:
    /// timer for selfField calculation
    IpplTimings::TimerRef selfFieldTimer_m;
    IpplTimings::TimerRef compPotenTimer_m;
    IpplTimings::TimerRef boundpTimer_m;
    IpplTimings::TimerRef statParamTimer_m;

    IpplTimings::TimerRef histoTimer_m;

gsell's avatar
gsell committed
514 515 516 517 518 519 520 521 522
private:

    const PartData *reference;
    void calcMoments();    // Calculates bunch moments using only emitted particles.
    void calcMomentsInitial(); // Calcualtes bunch moments by summing over bins (not accurate when any particles have been emitted).

    double calculateAngle(double x, double y);
    double calculateAngle2(double x, double y);

523
 public:
524
    void calcEMean(); // update eKin_m;
525 526
    void correctEnergy(double avrgp);
 private:
527

gsell's avatar
gsell committed
528 529 530 531
    /*
      Member variables starts here
    */

532 533 534 535
    // unit state of PartBunch
    UnitState_t unit_state_;
    UnitState_t stateOfLastBoundP_;

gsell's avatar
gsell committed
536
    /// hold the line-density
537
    std::unique_ptr<double[]> lineDensity_m;
gsell's avatar
gsell committed
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
    /// how many bins the line-density has
    unsigned int nBinsLineDensity_m;

    /// holds the centroid of the beam
    double centroid_m[2 * Dim];

    /// resize mesh to geometry specified
    void resizeMesh();

    /// 6x6 matrix of the moments of the beam
    FMatrix<double, 2 * Dim, 2 * Dim> moments_m;

    /// holds the timestep in seconds
    double dt_m;
    /// holds the actual time of the integration
    double t_m;
    /// mean energy of the bunch (MeV)
    double eKin_m;
    /// energy of the bunch
    double *energy_m;
    /// energy spread of the beam in keV
    double dE_m;

561 562
    /// Initialize the translation vector and rotation quaternion
    /// here. Cyclotron tracker will reset these values each timestep
kraus's avatar
kraus committed
563
    /// TTracker can just use 0 translation and 0 rotation (quat[1 0 0 0]).
564 565
    Vector_t globalMeanR_m = Vector_t(0.0, 0.0, 0.0);
    Quaternion_t globalToLocalQuaternion_m = Quaternion_t(1.0, 0.0, 0.0, 0.0);
566

gsell's avatar
gsell committed
567 568 569 570 571 572 573 574 575 576 577 578 579
    /// maximal extend of particles
    Vector_t rmax_m;
    /// minimal extend of particles
    Vector_t rmin_m;

    /// rms beam size (m)
    Vector_t rrms_m;
    /// rms momenta
    Vector_t prms_m;
    /// mean position (m)
    Vector_t rmean_m;
    /// mean momenta
    Vector_t pmean_m;
580

gsell's avatar
gsell committed
581 582
    /// rms emittance (not normalized)
    Vector_t eps_m;
583

gsell's avatar
gsell committed
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
    /// rms normalized emittance
    Vector_t eps_norm_m;
    /// rms correlation
    Vector_t rprms_m;

    /// dispersion x & y
    double Dx_m;
    double Dy_m;

    /// derivative of the dispersion
    double DDx_m;
    double DDy_m;

    /// meshspacing of cartesian mesh
    Vector_t hr_m;
    /// meshsize of cartesian mesh
    Vektor<int, 3> nr_m;

    /// for defining the boundary conditions
    BConds<double, 3, Mesh_t, Center_t> bc_m;
    BConds<Vector_t, 3, Mesh_t, Center_t> vbc_m;

    /// stores the used field solver
    FieldSolver *fs_m;

    double couplingConstant_m;

    double qi_m;

613 614
    bool interpolationCacheSet_m;

615
    ParticleAttrib<CacheDataCIC<double, 3U> > interpolationCache_m;
616

gsell's avatar
gsell committed
617 618 619 620 621 622 623 624
    /// counter to store the distributin dump
    int distDump_m;

    /*
      For saving particles on a collimator or
      dead particles (absobed)
    */

625 626 627 628 629 630 631 632 633 634 635
    // variables for stashing a bunch
    unsigned int stash_Nloc_m;
    Vector_t stash_iniR_m;
    Vector_t stash_iniP_m;
    PID_t stash_id_m;
    Ppos_t stash_r_m, stash_p_m, stash_x_m;
    ParticleAttrib<double> stash_q_m, stash_dt_m;
    ParticleAttrib<int> stash_bin_m;
    ParticleAttrib<long> stash_ls_m;
    ParticleAttrib<short> stash_ptype_m;
    bool bunchStashed_m;
gsell's avatar
gsell committed
636

637
    PartBunch &operator=(const PartBunch &) = delete;
638

gsell's avatar
gsell committed
639 640 641 642 643 644 645 646 647 648
    ///
    int fieldDBGStep_m;

    /// Mesh enlargement
    double dh_m; /// in % how much the mesh is enlarged

    /// if larger than 0, emitt particles for tEmission_m [s]
    double tEmission_m;

    /// holds the gamma of the bin
649
    std::unique_ptr<double[]> bingamma_m;
gsell's avatar
gsell committed
650 651 652 653

    //FIXME: this should go into the Bin class!
    // holds number of emitted particles of the bin
    // jjyang: opal-cycl use *nBin_m of pbin_m
654
    std::unique_ptr<size_t[]> binemitted_m;
gsell's avatar
gsell committed
655 656 657 658 659 660 661

    /// path length from the start point
    double lPath_m;

    /// steps per turn for OPAL-cycl
    int stepsPerTurn_m;

662 663 664 665 666
    /// step in a TRACK command
    long long localTrackStep_m;

    /// if multiple TRACK commands
    long long globalTrackStep_m;
gsell's avatar
gsell committed
667 668 669 670 671 672 673 674 675 676 677 678 679

    /// current bunch number
    int numBunch_m;

    /// this parameter records the current steps since last bunch injection
    /// it helps to inject new bunches correctly in the restart run of OPAL-cycl
    /// it is stored during phase space dump.
    int SteptoLastInj_m;

    /*
      Data structure for particle load balance information
    */

680 681
    std::unique_ptr<double[]> partPerNode_m;
    std::unique_ptr<double[]> globalPartPerNode_m;
gsell's avatar
gsell committed
682 683 684

    Distribution *dist_m;

685 686 687 688
 // flag to tell if we are a DC-beam
    bool dcBeam_m;


gsell's avatar
gsell committed
689 690
};

kraus's avatar
kraus committed
691
inline
692 693
bool PartBunch::isDcBeam() { return dcBeam_m;}

694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872
inline
Vector_t PartBunch::getStashIniP() const {
    return stash_iniP_m;
}

inline
PartBunch::UnitState_t PartBunch::getUnitState() const {
    return unit_state_;
}

inline
void PartBunch::switchToUnitlessPositions(bool use_dt_per_particle) {

    if(unit_state_ == unitless)
        throw SwitcherError("PartBunch::switchToUnitlessPositions",
                            "Cannot make a unitless PartBunch unitless");

    bool hasToReset = false;
    if(!R.isDirty()) hasToReset = true;

    for(size_t i = 0; i < getLocalNum(); i++) {
        double dt = getdT();
        if(use_dt_per_particle)
            dt = this->dt[i];

        R[i] /= Vector_t(Physics::c * dt);
        X[i] /= Vector_t(Physics::c * dt);
    }

    unit_state_ = unitless;

    if(hasToReset) R.resetDirtyFlag();
}

//FIXME: unify methods, use convention that all particles have own dt
inline
void PartBunch::switchOffUnitlessPositions(bool use_dt_per_particle) {

    if(unit_state_ == units)
        throw SwitcherError("PartBunch::switchOffUnitlessPositions",
                            "Cannot apply units twice to PartBunch");

    bool hasToReset = false;
    if(!R.isDirty()) hasToReset = true;

    for(size_t i = 0; i < getLocalNum(); i++) {
        double dt = getdT();
        if(use_dt_per_particle)
            dt = this->dt[i];

        R[i] *= Vector_t(Physics::c * dt);
        X[i] *= Vector_t(Physics::c * dt);
    }

    unit_state_ = units;

    if(hasToReset) R.resetDirtyFlag();
}


inline
double PartBunch::getRho(NDIndex<3> e) {
    return rho_m.localElement(e);
}

inline
double PartBunch::getRho(int x, int y, int z) {
    return rho_m[x][y][z].get();
}

inline
bool PartBunch::itIsMyTurn(int *n) {
    bool res = (*n == myNode_m);
    n++;
    if(*n == nodes_m) *n = 0;
    return res;
}

inline
void PartBunch::do_binaryRepart() {
    get_bounds(rmin_m, rmax_m);
    BinaryRepartition(*this);
    update();
    get_bounds(rmin_m, rmax_m);
    boundp();
}

inline
void PartBunch::set_nBinsLineDensity(int n) {
    nBinsLineDensity_m = n;
}

inline
void PartBunch::setGridIsFixed() {
    fixed_grid = true;
}

inline
bool PartBunch::isGridFixed() {
    return fixed_grid;
}

inline
void   PartBunch::setTEmission(double t) {
    tEmission_m = t;
}

inline
double PartBunch::getTEmission() {
    return tEmission_m;
}

inline
double PartBunch::getRebinEnergy() {
    return pbin_m->getRebinEnergy();
}

inline
void PartBunch::weHaveNOBins() {
    if(pbin_m != NULL)
        delete pbin_m;
    pbin_m = NULL;
}

inline
void PartBunch::rebin() {
    this->Bin = 0;
    pbin_m->resetBins();
    // delete pbin_m; we did not allocate it!
    pbin_m = NULL;
}

inline
int PartBunch::getNumBins() {
    if(pbin_m != NULL)
        return pbin_m->getNBins();
    else
        return 0;
}

inline
int PartBunch::getLastemittedBin() {
    if(pbin_m != NULL)
        return pbin_m->getLastemittedBin();
    else
        return 0;
}

inline
void PartBunch::updatePartInBin(size_t countLost[]) {
    if(pbin_m != NULL)
        pbin_m->updatePartInBin(countLost);
}

inline
double PartBunch::getBinGamma(int bin) {
    return bingamma_m[bin];
}

inline
void PartBunch::setBinCharge(int bin, double q) {
    this->Q = where(eq(this->Bin, bin), q, 0.0);
}

inline
void PartBunch::setBinCharge(int bin) {
    this->Q = where(eq(this->Bin, bin), this->Q, 0.0);
}

inline
std::pair<Vector_t, Vector_t> PartBunch::getEExtrema() {
    const Vector_t maxE = max(eg_m);
    //      const double maxL = max(dot(eg_m,eg_m));
    const Vector_t minE = min(eg_m);
    // INFOMSG("MaxE= " << maxE << " MinE= " << minE << endl);
    return std::pair<Vector_t, Vector_t>(maxE, minE);
}

inline
873
const Mesh_t &PartBunch::getMesh() const {
874 875 876 877
    return getLayout().getLayout().getMesh();
}

inline
878
Mesh_t &PartBunch::getMesh() {
879 880 881 882
    return getLayout().getLayout().getMesh();
}

inline
883
FieldLayout_t &PartBunch::getFieldLayout() {
884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987
    return dynamic_cast<FieldLayout_t &>(getLayout().getLayout().getFieldLayout());
}

inline
double PartBunch::getPx0(int i) {
    return 0;
}

inline
double PartBunch::getPy0(int i) {
    return 0;
}

inline
void PartBunch::setZ(int i, double zcoo) {};

inline
void PartBunch::setKR(Vector_t value, int i) {};

inline
void PartBunch::setKT(Vector_t value, int i) {};

inline
void PartBunch::BetOut(FILE *dat, FILE *sli) {};

inline
void PartBunch::plotR() {};

inline
void PartBunch::get_bounds(Vector_t &rmin, Vector_t &rmax) {
    bounds(this->R, rmin, rmax);
}

inline
void PartBunch::push_back(Particle p) {
    Inform msg("PartBunch ");

    create(1);
    size_t i = getTotalNum();

    R[i](0) = p[0];
    R[i](1) = p[1];
    R[i](2) = p[2];

    P[i](0) = p[3];
    P[i](1) = p[4];
    P[i](2) = p[5];

    update();
    msg << "Created one particle i= " << i << endl;
}

inline
void PartBunch::set_part(FVector<double, 6> z, int ii) {
    R[ii](0) = z[0];
    P[ii](0) = z[1];
    R[ii](1) = z[2];
    P[ii](1) = z[3];
    R[ii](2) = z[4];
    P[ii](2) = z[5];
}

inline
void PartBunch::set_part(Particle p, int ii) {
    R[ii](0) = p[0];
    P[ii](0) = p[1];
    R[ii](1) = p[2];
    P[ii](1) = p[3];
    R[ii](2) = p[4];
    P[ii](2) = p[5];
}

inline
Particle PartBunch::get_part(int ii) {
    Particle part;
    part[0] = R[ii](0);
    part[1] = P[ii](0);
    part[2] = R[ii](1);
    part[3] = P[ii](1);
    part[4] = R[ii](2);
    part[5] = P[ii](2);
    return part;
}

inline
void   PartBunch::setdT(double dt) {
    dt_m = dt;
}

inline
double PartBunch::getdT() const {
    return dt_m;
}

inline
void   PartBunch::setT(double t) {
    t_m = t;
}

inline
double PartBunch::getT() const {
    return t_m;
}

988 989 990
inline
void PartBunch::resetInterpolationCache(bool clearCache) {
    interpolationCacheSet_m = false;
991
    if(clearCache) {
992 993 994 995
        interpolationCache_m.destroy(interpolationCache_m.size(), 0, true);
    }
}

996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045
inline
bool PartBunch::isZPeriodic() const { // used in Distribution
    return false;
}

inline
double PartBunch::getGaBeLa() const { // used in Distribution
    return 1.0;
}

inline
double PartBunch::get_sPos() {
    if(sum(PType != 0)) {
        const size_t n = getLocalNum();
        size_t numPrimBeamParts = 0;
        double z = 0.0;
        if(n != 0) {
            for(size_t i = 0; i < n; i++) {
                if(PType[i] == 0) {
                    z += R[i](2);
                    numPrimBeamParts++;
                }
            }
            if(numPrimBeamParts != 0)
                z = z / numPrimBeamParts;
        }
        reduce(z, z, OpAddAssign());

        z = z / Ippl::getNodes();

        return z;
    } else {
        const size_t n = getTotalNum();
        if(n > 0)
            return sum(R(2)) / getTotalNum();
        else
            return 0.0;
    }
}

inline
double   PartBunch::get_phase() const {
    return 1.0;
}

inline
double   PartBunch::get_gamma() const {
    return 1.0;
}

1046

1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076
inline
double PartBunch::get_meanEnergy() const {
    return eKin_m;
}

// inline
// double* PartBunch::get_energy() {
//     return energy_m;
//}

inline
Vector_t PartBunch::get_origin() const {
    return rmin_m;
}

inline
Vector_t PartBunch::get_maxExtend() const {
    return rmax_m;
}

inline
Vector_t PartBunch::get_centroid() const {
    return rmean_m;
}

inline
Vector_t PartBunch::get_rrms() const {
    return rrms_m;
}

adelmann's avatar
adelmann committed
1077 1078 1079 1080 1081
inline
Vector_t PartBunch::get_rprrms() const {
    return rprms_m;
}

1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295
inline
Vector_t PartBunch::get_rmean() const {
    return rmean_m;
}

inline
Vector_t PartBunch::get_prms() const {
    return prms_m;
}

inline
Vector_t PartBunch::get_pmean() const {
    return pmean_m;
}

inline
Vector_t PartBunch::get_emit() const {
    return eps_m;
}

inline
Vector_t PartBunch::get_norm_emit() const {
    return eps_norm_m;
}

inline
Vector_t PartBunch::get_hr() const {
    return hr_m;
}

inline
double PartBunch::get_Dx() const {
    return Dx_m;
}

inline
double PartBunch::get_Dy() const {
    return Dy_m;
}

inline
double PartBunch::get_DDx() const {
    return DDx_m;
}

inline
double PartBunch::get_DDy() const {
    return DDy_m;
}

inline
void PartBunch::set_meshEnlargement(double dh) {
    dh_m = dh;
}

inline
double PartBunch::get_meshEnlargement() const {
    return dh_m;
}

inline
size_t PartBunch::getLoadBalance(int p) const {
    return globalPartPerNode_m[p];
}

inline
void PartBunch::get_PBounds(Vector_t &min, Vector_t &max) const {
    bounds(this->P, min, max);
}

inline
double PartBunch::getCouplingConstant() const {
    return couplingConstant_m;
}

inline
void PartBunch::setCouplingConstant(double c) {
    couplingConstant_m = c;
}

inline
void PartBunch::setCharge(double q) {
    qi_m = q;
    if(getTotalNum() != 0)
        Q = qi_m;
    else
        WARNMSG("Could not set total charge in PartBunch::setCharge based on getTotalNum" << endl);
}

inline
void PartBunch::setChargeZeroPart(double q) {
    qi_m = q;
}

inline
void PartBunch::setMass(double mass) {
    M = mass;
}

inline
double PartBunch::getCharge() const {
    return sum(Q);
}

inline
double PartBunch::getChargePerParticle() const {
    return qi_m;
}

inline
void PartBunch::setLPath(double s) {
    lPath_m = s;
}

inline
double PartBunch::getLPath() const {
    return lPath_m;
}

inline
void PartBunch::setStepsPerTurn(int n) {
    stepsPerTurn_m = n;
}

inline
int PartBunch::getStepsPerTurn() const {
    return stepsPerTurn_m;
}

inline
void PartBunch::setNumBunch(int n) {
    numBunch_m = n;
}

inline
int PartBunch::getNumBunch() const {
    return numBunch_m;
}

inline
void PartBunch::setSteptoLastInj(int n) {
    SteptoLastInj_m = n;
}

inline
int PartBunch::getSteptoLastInj() {
    return SteptoLastInj_m;
}

inline
double PartBunch::getQ() const {
    return reference->getQ();
}

inline
double PartBunch::getM() const {
    return reference->getM();
}

inline
double PartBunch::getP() const {
    return reference->getP();
}

inline
double PartBunch::getE() const {
    return reference->getE();
}

inline
void PartBunch::resetQ(double q)  {
    const_cast<PartData *>(reference)->setQ(q);
}

inline
void PartBunch::resetM(double m)  {
    const_cast<PartData *>(reference)->setM(m);
}

inline
double PartBunch::getdE() {
    return dE_m;
}

inline
double PartBunch::getBeta() const {
    return reference->getBeta();
}

inline
double PartBunch::getGamma() const {
    return reference->getGamma();
}

inline
double PartBunch::getGamma(int i) {
    return 0;
}

inline
double PartBunch::getBeta(int i) {
    return 0;
}

inline
void PartBunch::actT() {};

inline
const PartData *PartBunch::getReference() const {
    return reference;
}

inline
Inform &operator<<(Inform &os, PartBunch &p) {
gsell's avatar
gsell committed
1296 1297 1298 1299 1300
    return p.print(os);
}


#endif // OPAL_PartBunch_HH