PartBunchBase.h 20.5 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
//
// Class PartBunchBase
//   Base class for representing particle bunches.
//
// Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
18 19 20
#ifndef PART_BUNCH_BASE_H
#define PART_BUNCH_BASE_H

21 22
#include "Particle/AbstractParticle.h"
#include "Particle/ParticleAttrib.h"
kraus's avatar
kraus committed
23 24
#include "Utility/IpplTimings.h"
#include "Utilities/GeneralClassicException.h"
25

kraus's avatar
kraus committed
26
#include "Algorithms/DistributionMoments.h"
27
#include "Algorithms/CoordinateSystemTrafo.h"
28 29 30 31 32
#include "Algorithms/OpalParticle.h"
#include "Algorithms/PBunchDefs.h"
#include "Algorithms/Quaternion.h"
#include "Algorithms/Vektor.h"

33 34 35
#include "FixedAlgebra/FMatrix.h"
#include "FixedAlgebra/FVector.h"

36 37
#include <memory>
#include <utility>
38 39
#include <vector>

frey_m's avatar
frey_m committed
40
class Distribution;
41 42 43 44 45
class FieldSolver;
class PartBins;
class PartBinsCyc;
class PartData;

46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
enum class ParticleOrigin:short {REGULAR,
                               SECONDARY,
                               STRIPPED};

enum class ParticleType:short {ELECTRON,
                               PROTON,
                               POSITRON,
                               ANTIPROTON,
                               CARBON,
                               HMINUS,
                               URANIUM,
                               MUON,
                               DEUTERON,
                               XENON,
                               H2P,
                               ALPHA,
                               HYDROGEN,
                               H3P,
                               UNNAMED};
65 66

template <class T, unsigned Dim>
kraus's avatar
kraus committed
67
class PartBunchBase : std::enable_shared_from_this<PartBunchBase<T, Dim>>
68
{
69
public:
70
    typedef typename AbstractParticle<T, Dim>::ParticlePos_t   ParticlePos_t;
71
    typedef typename AbstractParticle<T, Dim>::ParticleIndex_t ParticleIndex_t;
72 73
    typedef typename AbstractParticle<T, Dim>::UpdateFlags     UpdateFlags_t;
    typedef typename AbstractParticle<T, Dim>::Position_t      Position_t;
kraus's avatar
kraus committed
74

75
    typedef std::pair<Vector_t, Vector_t> VectorPair_t;
kraus's avatar
kraus committed
76

77
    static const unsigned Dimension = Dim;
kraus's avatar
kraus committed
78

79
    enum UnitState_t { units = 0, unitless = 1 };
kraus's avatar
kraus committed
80

81
public:
frey_m's avatar
frey_m committed
82
    virtual ~PartBunchBase() { }
kraus's avatar
kraus committed
83

84
    PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData* ref);
kraus's avatar
kraus committed
85

86
    PartBunchBase(const PartBunchBase& rhs) = delete; // implement if needed
kraus's avatar
kraus committed
87

88 89 90
    /*
     * Bunch common member functions
     */
kraus's avatar
kraus committed
91

92
    // This is required since we initialize the Layout and the RegionLayout with default constructor
93
    virtual void initialize(FieldLayout_t* fLayout) = 0;
kraus's avatar
kraus committed
94

95
    bool getIfBeamEmitting();
kraus's avatar
kraus committed
96

97
    int getLastEmittedEnergyBin();
kraus's avatar
kraus committed
98

99
    size_t getNumberOfEmissionSteps();
kraus's avatar
kraus committed
100

101
    int getNumberOfEnergyBins();
kraus's avatar
kraus committed
102

103
    void Rebin();
kraus's avatar
kraus committed
104

105
    void setEnergyBins(int numberOfEnergyBins);
kraus's avatar
kraus committed
106

107
    bool weHaveEnergyBins();
kraus's avatar
kraus committed
108

109
    //FIXME: unify methods, use convention that all particles have own dt
110
    void switchToUnitlessPositions(bool use_dt_per_particle = false);
kraus's avatar
kraus committed
111

112
    //FIXME: unify methods, use convention that all particles have own dt
113
    void switchOffUnitlessPositions(bool use_dt_per_particle = false);
kraus's avatar
kraus committed
114

115 116 117
    void setDistribution(Distribution* d,
                         std::vector<Distribution*> addedDistributions,
                         size_t& np);
kraus's avatar
kraus committed
118

119
    bool isGridFixed() const;
120

121
    bool hasBinning() const;
kraus's avatar
kraus committed
122 123


124 125 126
    /*
       Energy bins related functions
     */
kraus's avatar
kraus committed
127

128
    void setTEmission(double t);
kraus's avatar
kraus committed
129

130
    double getTEmission();
kraus's avatar
kraus committed
131

132
    bool doEmission();
kraus's avatar
kraus committed
133

134
    bool weHaveBins() const;
kraus's avatar
kraus committed
135

136
    void setPBins(PartBins* pbin);
kraus's avatar
kraus committed
137

138
    void setPBins(PartBinsCyc* pbin);
kraus's avatar
kraus committed
139

140 141 142 143
    /** \brief Emit particles in the given bin
        i.e. copy the particles from the bin structure into the
        particle container
    */
144
    size_t emitParticles(double eZ);
kraus's avatar
kraus committed
145

146
    void updateNumTotal();
kraus's avatar
kraus committed
147

148
    void rebin();
kraus's avatar
kraus committed
149

150
    int getLastemittedBin();
kraus's avatar
kraus committed
151

152 153
    void setLocalBinCount(size_t num, int bin);

154
    /** \brief Compute the gammas of all bins */
155
    void calcGammas();
kraus's avatar
kraus committed
156

157
    void calcGammas_cycl();
kraus's avatar
kraus committed
158

159 160
    /** \brief Get gamma of one bin */
    double getBinGamma(int bin);
kraus's avatar
kraus committed
161

162
    /** \brief Set the charge of one bin to the value of q and all other to zero */
163
    virtual void setBinCharge(int bin, double q);
kraus's avatar
kraus committed
164

165
    /** \brief Set the charge of all other the ones in bin to zero */
166
    virtual void setBinCharge(int bin);
kraus's avatar
kraus committed
167

168
    /** \brief returns the number of particles outside of a box defined by x */
169
    size_t calcNumPartsOutside(Vector_t x);
kraus's avatar
kraus committed
170

171 172
    void calcLineDensity(unsigned int nBins, std::vector<double>& lineDensity,
                         std::pair<double, double>& meshInfo);
kraus's avatar
kraus committed
173

174 175
    void setBeamFrequency(double v);

176 177 178
    /*
       Mesh and Field Layout related functions
     */
kraus's avatar
kraus committed
179

180
    virtual void boundp();
kraus's avatar
kraus committed
181

frey_m's avatar
frey_m committed
182
    /** delete particles which are too far away from the center of beam*/
183
    void boundp_destroyCycl();
kraus's avatar
kraus committed
184

185
    /** This is only temporary in order to get the collimator and pepperpot working */
frey_m's avatar
frey_m committed
186
    size_t boundp_destroyT();
kraus's avatar
kraus committed
187

frey_m's avatar
frey_m committed
188
    size_t destroyT();
kraus's avatar
kraus committed
189

190 191 192 193 194
    /*
       Read out coordinates
     */
    virtual double getPx(int i);
    virtual double getPy(int i);
195 196 197
    virtual double getPz(int i);

    virtual double getPx0(int i);
198
    virtual double getPy0(int i);
199 200 201 202 203 204 205 206 207

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

    virtual void setZ(int i, double zcoo);
kraus's avatar
kraus committed
208

209
    void get_bounds(Vector_t& rmin, Vector_t& rmax);
kraus's avatar
kraus committed
210

211
    void getLocalBounds(Vector_t& rmin, Vector_t& rmax);
kraus's avatar
kraus committed
212

213
    std::pair<Vector_t, double> getBoundingSphere();
kraus's avatar
kraus committed
214

215
    std::pair<Vector_t, double> getLocalBoundingSphere();
kraus's avatar
kraus committed
216 217


218 219 220
    /*
       Compatibility function push_back
     */
kraus's avatar
kraus committed
221 222 223 224 225 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 286 287 288 289 290 291 292 293
    void push_back(OpalParticle const& p);

    void setParticle(FVector<double, 6> z, int ii);

    void setParticle(OpalParticle const& p, int ii);

    OpalParticle getParticle(int ii);

    class ConstIterator {
        friend class PartBunchBase<T, Dim>;

    public:
        ConstIterator():
            bunch_m(nullptr),
            index_m(0)
        {}
        ConstIterator(PartBunchBase const* bunch, unsigned int i):
            bunch_m(bunch),
            index_m(i)
        {}

        ~ConstIterator()
        {}

        bool operator == (ConstIterator const& rhs) const
        {
            return bunch_m == rhs.bunch_m && index_m == rhs.index_m;
        }

        bool operator != (ConstIterator const& rhs) const
        {
            return bunch_m != rhs.bunch_m || index_m != rhs.index_m;
        }

        OpalParticle operator*() const
        {
            if (index_m >= bunch_m->getLocalNum()) {
                throw GeneralClassicException("PartBunchBase::ConstIterator::operator*", "out of bounds");
            }
            return OpalParticle(bunch_m->ID[index_m],
                                bunch_m->R[index_m],
                                bunch_m->P[index_m],
                                bunch_m->getT(),
                                bunch_m->Q[index_m],
                                bunch_m->getM() * 1e-6);
        }

        ConstIterator operator++()
        {
            ++index_m;
            return *this;
        }

        ConstIterator operator++(int)
        {
            ConstIterator it = *this;
            ++index_m;

            return it;
        }

        int operator-(const ConstIterator& other) const
        {
            return index_m - other.index_m;
        }
    private:
        PartBunchBase const* bunch_m;
        unsigned int index_m;
    };

    ConstIterator begin() const {
        return ConstIterator(this, 0);
    }
kraus's avatar
kraus committed
294

kraus's avatar
kraus committed
295 296 297
    ConstIterator end() const {
        return ConstIterator(this, getLocalNum());
    }
kraus's avatar
kraus committed
298

299 300 301 302
    /// 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].
303 304
    void maximumAmplitudes(const FMatrix<double, 6, 6>& D,
                           double& axmax, double& aymax);
kraus's avatar
kraus committed
305

306
    void   setdT(double dt);
307 308 309
    double getdT() const;

    void   setT(double t);
310
    void   incrementT();
311
    double getT() const;
kraus's avatar
kraus committed
312

313 314 315 316 317 318
    /**
     * get the spos of the primary beam
     *
     * @param none
     *
     */
319
    double get_sPos() const;
kraus's avatar
kraus committed
320

321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
    void set_sPos(double s);

    double get_gamma() const;
    double get_meanKineticEnergy() const;
    Vector_t get_origin() const;
    Vector_t get_maxExtent() const;
    Vector_t get_centroid() const;
    Vector_t get_rrms() const;
    Vector_t get_rprms() const;
    Vector_t get_rmean() const;
    Vector_t get_prms() const;
    Vector_t get_pmean() const;
    Vector_t get_pmean_Distribution() const;
    Vector_t get_emit() const;
    Vector_t get_norm_emit() const;
336
    Vector_t get_halo() const;
frey_m's avatar
frey_m committed
337
    virtual Vector_t get_hr() const;
338 339 340 341 342

    double get_Dx() const;
    double get_Dy() const;
    double get_DDx() const;
    double get_DDy() const;
kraus's avatar
kraus committed
343

frey_m's avatar
frey_m committed
344
    virtual void set_meshEnlargement(double dh);
345 346 347

    void gatherLoadBalanceStatistics();
    size_t getLoadBalance(int p) const;
kraus's avatar
kraus committed
348

349
    void get_PBounds(Vector_t &min, Vector_t &max) const;
kraus's avatar
kraus committed
350

351 352
    void calcBeamParameters();
    void calcBeamParametersInitial(); // Calculate initial beam parameters before emission.
kraus's avatar
kraus committed
353

354 355
    double getCouplingConstant() const;
    void setCouplingConstant(double c);
kraus's avatar
kraus committed
356

357 358 359 360 361 362 363
    // set the charge per simulation particle
    void setCharge(double q);
    // set the charge per simulation particle when total particle number equals 0
    void setChargeZeroPart(double q);

    // set the mass per simulation particle
    void setMass(double mass);
kraus's avatar
kraus committed
364
    void setMassZeroPart(double mass);
kraus's avatar
kraus committed
365

366 367 368 369 370
    /// get the total charge per simulation particle
    double getCharge() const;

    /// get the macro particle charge
    double getChargePerParticle() const;
kraus's avatar
kraus committed
371

kraus's avatar
kraus committed
372 373
    double getMassPerParticle() const;

374
    virtual void setSolver(FieldSolver *fs);
kraus's avatar
kraus committed
375

376
    bool hasFieldSolver();
kraus's avatar
kraus committed
377

378
    std::string getFieldSolverType() const;
kraus's avatar
kraus committed
379

380 381 382 383 384 385 386 387 388 389 390 391
    void setStepsPerTurn(int n);
    int getStepsPerTurn() const;

    /// step in multiple TRACK commands
    void setGlobalTrackStep(long long n);
    long long getGlobalTrackStep() const;

    /// step in a TRACK command
    void setLocalTrackStep(long long n);
    void incTrackSteps();
    long long getLocalTrackStep() const;

frey_m's avatar
frey_m committed
392 393 394 395 396 397 398 399 400 401 402 403 404 405
    void setNumBunch(short n);
    short getNumBunch() const;

    // used in ParallelCyclotronTracker for multi-bunch mode
    void setTotalNumPerBunch(size_t numpart, short n);
    size_t getTotalNumPerBunch(short n) const;

    void setLocalNumPerBunch(size_t numpart, short n);
    size_t getLocalNumPerBunch(short n) const;

    /* used in initializeTracking_m of ParallelCyclotronTracker
     * for multi-bunch mode
     */
    void countTotalNumPerBunch();
406 407 408 409 410

    void setGlobalMeanR(Vector_t globalMeanR);
    Vector_t getGlobalMeanR();
    void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion);
    Quaternion_t getGlobalToLocalQuaternion();
kraus's avatar
kraus committed
411

412
    void setSteptoLastInj(int n);
413
    int getSteptoLastInj() const;
kraus's avatar
kraus committed
414

415 416
    /// calculate average angle of longitudinal direction of bins
    double calcMeanPhi();
kraus's avatar
kraus committed
417

418 419
    /// 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);
kraus's avatar
kraus committed
420

frey_m's avatar
frey_m committed
421 422
    bool resetPartBinBunch();

423
    ///@{ Access to reference data
424 425 426 427
    double getQ() const;
    double getM() const;
    double getP() const;
    double getE() const;
428 429 430
    ParticleOrigin getPOrigin() const;
    ParticleType getPType() const;
    std::string getPTypeString();
431 432 433 434
    double getInitialBeta() const;
    double getInitialGamma() const;
    ///@}
    ///@{ Set reference data
435 436
    void resetQ(double q);
    void resetM(double m);
437 438
    void setPOrigin(ParticleOrigin);
    void setPType(std::string type);
439
    ///@}
440
    double getdE() const;
441 442 443
    virtual double getGamma(int i);
    virtual double getBeta(int i);
    virtual void actT();
kraus's avatar
kraus committed
444

445
    const PartData* getReference() const;
kraus's avatar
kraus committed
446

447
    double getEmissionDeltaT();
kraus's avatar
kraus committed
448

449 450 451 452 453 454
    Quaternion_t getQKs3D();
    void         setQKs3D(Quaternion_t q);
    Vector_t     getKs3DRefr();
    void         setKs3DRefr(Vector_t r);
    Vector_t     getKs3DRefp();
    void         setKs3DRefp(Vector_t p);
kraus's avatar
kraus committed
455

456
    void iterateEmittedBin(int binNumber);
kraus's avatar
kraus committed
457

458
    void calcEMean();
kraus's avatar
kraus committed
459

460
    Inform& print(Inform& os);
kraus's avatar
kraus committed
461

462
    /*
kraus's avatar
kraus committed
463
     * (Pure) virtual member functions
464
     */
kraus's avatar
kraus committed
465

466
    virtual void runTests();
kraus's avatar
kraus committed
467

frey_m's avatar
frey_m committed
468
    virtual void do_binaryRepart();
kraus's avatar
kraus committed
469

470
    virtual void resetInterpolationCache(bool clearCache = false);
kraus's avatar
kraus committed
471

472
    //brief calculates back the max/min of the efield on the grid
473
    virtual VectorPair_t getEExtrema() = 0;
kraus's avatar
kraus committed
474

475
    virtual double getRho(int x, int y, int z) = 0;
kraus's avatar
kraus committed
476

477
    virtual void computeSelfFields() = 0;
478

479
    //brief used for self fields with binned distribution
480
    virtual void computeSelfFields(int bin) = 0;
481

482
    virtual void computeSelfFields_cycl(double gamma) = 0;
483
    virtual void computeSelfFields_cycl(int bin) = 0;
kraus's avatar
kraus committed
484

frey_m's avatar
frey_m committed
485
    virtual void swap(unsigned int i, unsigned int j);
kraus's avatar
kraus committed
486

frey_m's avatar
frey_m committed
487 488 489
    /*
       Mesh and Field Layout related functions
     */
kraus's avatar
kraus committed
490

frey_m's avatar
frey_m committed
491 492 493 494
    virtual void setBCAllPeriodic();
    virtual void setBCAllOpen();

    virtual void setBCForDCBeam();
kraus's avatar
kraus committed
495 496


497
//     virtual void setMesh(Mesh_t* mesh) = 0;
498
//     virtual Mesh_t &getMesh() = 0;
kraus's avatar
kraus committed
499

500
//     virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
501
    virtual FieldLayout_t& getFieldLayout() = 0;
kraus's avatar
kraus committed
502

503 504
    virtual void resizeMesh() { };

505 506 507
    /*
     * Wrapped member functions of IpplParticleBase
     */
kraus's avatar
kraus committed
508

509 510 511 512
    size_t getTotalNum() const;
    void setTotalNum(size_t n);
    void setLocalNum(size_t n);
    size_t getLocalNum() const;
kraus's avatar
kraus committed
513

514 515
    size_t getDestroyNum() const;
    size_t getGhostNum() const;
kraus's avatar
kraus committed
516

517
    ParticleLayout<T, Dim>& getLayout();
518
    const ParticleLayout<T, Dim>& getLayout() const;
kraus's avatar
kraus committed
519

520 521
    bool getUpdateFlag(UpdateFlags_t f) const;
    void setUpdateFlag(UpdateFlags_t f, bool val);
kraus's avatar
kraus committed
522 523


524
    ParticleBConds<Position_t, Dimension>& getBConds() {
kraus's avatar
kraus committed
525
        return pbase_m->getBConds();
526
    }
kraus's avatar
kraus committed
527

528
    void setBConds(const ParticleBConds<Position_t, Dimension>& bc) {
kraus's avatar
kraus committed
529
        pbase_m->setBConds(bc);
530
    }
kraus's avatar
kraus committed
531

532
    bool singleInitNode() const;
kraus's avatar
kraus committed
533

534
    void resetID();
kraus's avatar
kraus committed
535

536 537
    void update();
    void update(const ParticleAttrib<char>& canSwap);
kraus's avatar
kraus committed
538

539 540 541
    void createWithID(unsigned id);
    void create(size_t M);
    void globalCreate(size_t np);
kraus's avatar
kraus committed
542

543 544 545
    void destroy(size_t M, size_t I, bool doNow = false);
    void performDestroy(bool updateLocalNum = false);
    void ghostDestroy(size_t M, size_t I);
kraus's avatar
kraus committed
546

547
protected:
548
    size_t calcMoments();    // Calculates bunch moments using only emitted particles.
kraus's avatar
kraus committed
549

snuverink_j's avatar
snuverink_j committed
550
    /* Calculates bunch moments by summing over bins
551 552 553
     * (not accurate when any particles have been emitted).
     */
    void calcMomentsInitial();
snuverink_j's avatar
snuverink_j committed
554
    /// angle range [0~2PI) degree
555
    double calculateAngle(double x, double y);
kraus's avatar
kraus committed
556 557


frey_m's avatar
frey_m committed
558
private:
frey_m's avatar
frey_m committed
559
    virtual void updateDomainLength(Vektor<int, 3>& grid) = 0;
kraus's avatar
kraus committed
560

frey_m's avatar
frey_m committed
561
    virtual void updateFields(const Vector_t& hr, const Vector_t& origin);
kraus's avatar
kraus committed
562

frey_m's avatar
frey_m committed
563
    void setup(AbstractParticle<T, Dim>* pb);
kraus's avatar
kraus committed
564

565 566 567 568 569 570
public:
    /*
     * Bunch attributes
     */
    ParticlePos_t& R;
    ParticleIndex_t& ID;
kraus's avatar
kraus committed
571

572
    // Particle container attributes
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
    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< double >       Phi;    // the electric potential
    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< ParticleType > PType;  // particle names
    ParticleAttrib< ParticleOrigin > POrigin;  // 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.
    ParticleAttrib< short >        cavityGapCrossed; // particle just crossed cavity gap (for ParallelCyclotronTracker)
    ParticleAttrib< short >        bunchNum; // bunch number to which particle belongs (multi-bunch mode)
kraus's avatar
kraus committed
588

589 590
    Vector_t RefPartR_m;
    Vector_t RefPartP_m;
591

592
    CoordinateSystemTrafo toLabTrafo_m;
kraus's avatar
kraus committed
593

594 595 596
    ParticleOrigin refPOrigin_m;
    ParticleType refPType_m;

597
    // The structure for particle binning
598
    PartBins* pbin_m;
599 600 601 602

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

604 605 606
    // For AMTS integrator in OPAL-T
    double dtScInit_m, deltaTau_m;

607 608 609
    // get 2nd order momentum matrix
    FMatrix<double, 2 * Dim, 2 * Dim> getSigmaMatrix();

610 611 612 613 614 615 616
private:
    // save particles in case of one core
    std::unique_ptr<Inform> pmsg_m;
    std::unique_ptr<std::ofstream> f_stream;
    /// if the grid does not have to adapt
    bool fixed_grid;

frey_m's avatar
frey_m committed
617
protected:
618
    IpplTimings::TimerRef boundpTimer_m;
619 620
    IpplTimings::TimerRef boundpBoundsTimer_m;
    IpplTimings::TimerRef boundpUpdateTimer_m;
621 622 623
    IpplTimings::TimerRef statParamTimer_m;

    IpplTimings::TimerRef histoTimer_m;
624 625
    /// timer for selfField calculation
    IpplTimings::TimerRef selfFieldTimer_m;
kraus's avatar
kraus committed
626

627
    const PartData* reference;
kraus's avatar
kraus committed
628

629 630 631
    /*
       Member variables starts here
     */
632 633 634 635 636 637 638

    // unit state of PartBunch
    UnitState_t unit_state_;
    UnitState_t stateOfLastBoundP_;

    /// holds the centroid of the beam
    double centroid_m[2 * Dim];
kraus's avatar
kraus committed
639

640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668
    /// 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;
    /// the position along design trajectory
    double spos_m;

    /// Initialize the translation vector and rotation quaternion
    /// here. Cyclotron tracker will reset these values each timestep
    /// TTracker can just use 0 translation and 0 rotation (quat[1 0 0 0]).
    //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);
    Vector_t globalMeanR_m;
    Quaternion_t globalToLocalQuaternion_m;

    /// maximal extend of particles
    Vector_t rmax_m;
    /// minimal extend of particles
    Vector_t rmin_m;

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

    /// stores the used field solver
669
    FieldSolver* fs_m;
670 671 672 673

    double couplingConstant_m;

    double qi_m;
kraus's avatar
kraus committed
674
    double massPerParticle_m;
kraus's avatar
kraus committed
675

676
    /// counter to store the distribution dump
677
    int distDump_m;
kraus's avatar
kraus committed
678

679
    /// Mesh enlargement
680
    double dh_m; /// relative enlargement of the mesh
681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696

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

    /// holds the gamma of the bin
    std::unique_ptr<double[]> bingamma_m;

    //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
    std::unique_ptr<size_t[]> binemitted_m;

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

    /// step in a TRACK command
697
    long long localTrackStep_m;
698 699 700 701 702

    /// if multiple TRACK commands
    long long globalTrackStep_m;

    /// current bunch number
frey_m's avatar
frey_m committed
703 704 705 706 707
    short numBunch_m;

    /// number of particles per bunch
    std::vector<size_t> bunchTotalNum_m;
    std::vector<size_t> bunchLocalNum_m;
708 709 710 711 712 713 714 715 716

    /// 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
    */
kraus's avatar
kraus committed
717

718 719
    std::unique_ptr<size_t[]> globalPartPerNode_m;

kraus's avatar
kraus committed
720 721
    Distribution *dist_m;
    DistributionMoments momentsComputer_m;
722 723 724

    // flag to tell if we are a DC-beam
    bool dcBeam_m;
725
    double periodLength_m;
kraus's avatar
kraus committed
726
    std::shared_ptr<AbstractParticle<T, Dim> > pbase_m;
727 728
};

kraus's avatar
kraus committed
729 730 731 732 733 734 735 736 737 738
template<class T, unsigned Dim>
typename PartBunchBase<T, Dim>::ConstIterator begin(PartBunchBase<T, Dim> const& bunch) {
    return bunch.begin();
}

template<class T, unsigned Dim>
typename PartBunchBase<T, Dim>::ConstIterator end(PartBunchBase<T, Dim> const& bunch) {
    return bunch.end();
}

739 740
#include "PartBunchBase.hpp"

741
#endif