ThickTracker.cpp 17.8 KB
Newer Older
gsell's avatar
gsell committed
1
//
2 3
// Source file of the ThickTracker class
//   Tracks using thick-lens algorithm.
gsell's avatar
gsell committed
4
//
5 6
// Copyright (c) 2018, Philippe Ganz, ETH Zürich
// All rights reserved
gsell's avatar
gsell committed
7
//
8
// OPAL is licensed under GNU GPL version 3.
9

snuverink_j's avatar
snuverink_j committed
10
#include <cfloat>
11 12
#include <fstream>
#include <typeinfo>
13

snuverink_j's avatar
snuverink_j committed
14
#include "Algorithms/ThickTracker.h"
kraus's avatar
kraus committed
15

gsell's avatar
gsell committed
16
#include "Beamlines/Beamline.h"
17
#include "Beamlines/FlaggedBeamline.h"
gsell's avatar
gsell committed
18

snuverink_j's avatar
snuverink_j committed
19
#include "Classic/Algorithms/PartData.h"
20

Philippe Ganz's avatar
Philippe Ganz committed
21 22
#include "Utilities/Options.h"
#include "Utilities/Util.h"
23
#include "Utilities/Timer.h"
gsell's avatar
gsell committed
24 25

#include "Physics/Physics.h"
26

27
//
gsell's avatar
gsell committed
28 29
// Class ThickTracker
// ------------------------------------------------------------------------
30
//
gsell's avatar
gsell committed
31 32
ThickTracker::ThickTracker(const Beamline &beamline,
                           const PartData &reference,
snuverink_j's avatar
snuverink_j committed
33 34 35 36 37 38
                           bool revBeam, bool revTrack)
    : Tracker(beamline, reference, revBeam, revTrack)
    , hamiltonian_m(1)
    , RefPartR_m(0.0)
    , RefPartP_m(0.0)
    , itsDataSink_m(nullptr)
39
    , itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection())
snuverink_j's avatar
snuverink_j committed
40 41 42 43 44 45 46
    , zstart_m(0.0)
    , zstop_m(0.0)
    , threshold_m(1.0e-6)
    , truncOrder_m(1) // linear
    , mapCreation_m(   IpplTimings::getTimer("mapCreation"))
    , mapCombination_m(IpplTimings::getTimer("mapCombination"))
    , mapTracking_m(   IpplTimings::getTimer("mapTracking"))
47
{
snuverink_j's avatar
snuverink_j committed
48
    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
49
                                   beamline.getInitialDirection());
snuverink_j's avatar
snuverink_j committed
50
    referenceToLabCSTrafo_m = labToRef.inverted();
51
}
gsell's avatar
gsell committed
52 53 54


ThickTracker::ThickTracker(const Beamline &beamline,
snuverink_j's avatar
snuverink_j committed
55
                           PartBunchBase<double, 3> *bunch,
56
                           Beam &/*beam*/,
snuverink_j's avatar
snuverink_j committed
57 58 59
                           DataSink &ds,
                           const PartData &reference,
                           bool revBeam, bool revTrack,
60
                           const std::vector<unsigned long long> &/*maxSteps*/,
snuverink_j's avatar
snuverink_j committed
61 62
                           double zstart,
                           const std::vector<double> &zstop,
63
                           const std::vector<double> &/*dt*/,
snuverink_j's avatar
snuverink_j committed
64 65 66 67 68 69
                           const int& truncOrder)
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , hamiltonian_m(truncOrder)
    , RefPartR_m(0.0)
    , RefPartP_m(0.0)
    , itsDataSink_m(&ds)
70
    , itsOpalBeamline_m(beamline.getOrigin3D(), beamline.getInitialDirection())
snuverink_j's avatar
snuverink_j committed
71 72 73 74 75 76 77
    , zstart_m(zstart)
    , zstop_m(zstop[0])
    , threshold_m(1.0e-6)
    , truncOrder_m(truncOrder)
    , mapCreation_m(   IpplTimings::getTimer("mapCreation"))
    , mapCombination_m(IpplTimings::getTimer("mapCombination"))
    , mapTracking_m(   IpplTimings::getTimer("mapTracking"))
78
{
snuverink_j's avatar
snuverink_j committed
79 80 81 82 83 84
    if ( zstop.size() > 1 )
        throw OpalException("ThickTracker::ThickTracker()",
                            "Multiple tracks not yet supported.");


    CoordinateSystemTrafo labToRef(beamline.getOrigin3D(),
85
                                   beamline.getInitialDirection());
snuverink_j's avatar
snuverink_j committed
86
    referenceToLabCSTrafo_m = labToRef.inverted();
87
}
gsell's avatar
gsell committed
88 89 90 91 92 93


ThickTracker::~ThickTracker()
{}


94 95

void ThickTracker::visitBeamline(const Beamline &bl) {
96

97 98
    const FlaggedBeamline* fbl = static_cast<const FlaggedBeamline*>(&bl);
    if (fbl->getRelativeFlag()) {
99
        *gmsg << " do stuff" << endl;
100
        OpalBeamline stash(fbl->getOrigin3D(), fbl->getInitialDirection());
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
        stash.swap(itsOpalBeamline_m);
        fbl->iterate(*this, false);
        itsOpalBeamline_m.prepareSections();
        itsOpalBeamline_m.compute3DLattice();
        stash.merge(itsOpalBeamline_m);
        stash.swap(itsOpalBeamline_m);
    } else {
        fbl->iterate(*this, false);
    }
}


void ThickTracker::prepareSections() {
    itsBeamline_m.accept(*this);
    itsOpalBeamline_m.prepareSections();
}




snuverink_j's avatar
snuverink_j committed
121 122 123 124
/*
//TODO complete and test fringe fields
void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
        double& beta0, double& gamma0, double& P0, double& q, std::array<double,2>& paramFringe, std::string e){
125

126

snuverink_j's avatar
snuverink_j committed
127 128 129 130
    series_t H;
    map_t tempMap, fieldMap;
    lstruct_t::iterator mBLit;
    double lenFringe = paramFringe[0] + paramFringe[1];
131

snuverink_j's avatar
snuverink_j committed
132 133
    //double intHeight = 0.1; //:TODO: change or argument for that value
    double stepSize = 0.01; //:TODO: change or argument for that value maybe minStepIfSC?
134 135


snuverink_j's avatar
snuverink_j committed
136
    //entrFringe = std::fabs(intHeight * std::tan(pSBend->getEntranceAngle()));
137

snuverink_j's avatar
snuverink_j committed
138 139 140
    //get proper slices for the fringe field
    double nSlices = std::ceil(2 * lenFringe / stepSize);
    stepSize = 2 * lenFringe / nSlices;
141

snuverink_j's avatar
snuverink_j committed
142 143 144
    //Create structMapTracking and push in BeamLine
    double K0= pSBend ->getB()*(Physics::c/itsReference.getP());
    K0= std::round(K0*1e6)/1e6 *q*(Physics::c/P0);
145

snuverink_j's avatar
snuverink_j committed
146 147
    double h = 1./ pSBend ->getBendRadius();               //inverse bending radius [1/m]
    K0<0 ? h = -1*h : h = 1*h;
148

snuverink_j's avatar
snuverink_j committed
149 150 151
    series_t az = - K0 / (2 * lenFringe) * (hamiltonian_m.y*hamiltonian_m.y - hamiltonian_m.x*hamiltonian_m.x) * std::tan(pSBend->getEntranceAngle());
    if (e == "out") {
            az = -az;
152
        }
snuverink_j's avatar
snuverink_j committed
153 154 155 156 157 158 159
    std::cout << az.getMaxOrder() << std::endl;
    for ( int i = 0 ; i < nSlices; i++ ){
        double l = stepSize * i;
        series_t ax = 0.5 * K0 / (2 * lenFringe) * (l*l - hamiltonian_m.y*hamiltonian_m.y);
        H = hamiltonian_m.bendFringe(beta0, gamma0, h, K0, ax, az);
        tempMap = ExpMap(-H * stepSize ,truncOrder_m);
        fieldMap = (tempMap * fieldMap).truncate(truncOrder_m);
160 161
    }

snuverink_j's avatar
snuverink_j committed
162 163 164 165 166 167
    structMapTracking fringeField;
    fringeField.elementName=pSBend->getName() + "accumulated Fringe Field Map";
    if (e == "in"){
        fringeField.elementPos=(pSBend->getElementPosition() - paramFringe[0]);
    } else if (e == "out") {
        fringeField.elementPos=(pSBend->getElementPosition() + pSBend->getArcLength() - paramFringe[1]);
168
    }
snuverink_j's avatar
snuverink_j committed
169 170 171 172
    fringeField.stepSize= lenFringe;
    fringeField.nSlices=1;
    fringeField.elementMap=fieldMap;
    mBL.push_back(fringeField);
173
}
174 175 176 177 178
*/


/**
 * @brief Algorithm for Thick Map-Tracking
snuverink_j's avatar
snuverink_j committed
179
 *
180 181 182
 */
void ThickTracker::execute() {

snuverink_j's avatar
snuverink_j committed
183 184 185
    /*
     * global settings
     */
frey_m's avatar
frey_m committed
186
    Inform msg("ThickTracker", *gmsg);
187 188 189 190 191

    OpalData::getInstance()->setInPrepState(true);

    OpalData::getInstance()->setGlobalPhaseShift(0.0);

snuverink_j's avatar
snuverink_j committed
192 193 194
    if ( OpalData::getInstance()->hasPriorTrack() ||
         OpalData::getInstance()->inRestartRun() )
    {
195
        OpalData::getInstance()->setOpenMode(OpalData::OPENMODE::APPEND);
196 197 198 199
    }

    prepareSections();

snuverink_j's avatar
snuverink_j committed
200
    msg << "Truncation order: " << this->truncOrder_m << endl;
201

Philippe Ganz's avatar
Philippe Ganz committed
202

snuverink_j's avatar
snuverink_j committed
203
    this->checkElementOrder_m();
204

snuverink_j's avatar
snuverink_j committed
205
    this->fillGaps_m();
Philippe Ganz's avatar
Philippe Ganz committed
206

snuverink_j's avatar
snuverink_j committed
207 208 209
    //FIXME in local system only
    RefPartR_m = Vector_t(0.0);
    RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
Philippe Ganz's avatar
Philippe Ganz committed
210

snuverink_j's avatar
snuverink_j committed
211
    itsBunch_m->set_sPos(zstart_m);
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
212

snuverink_j's avatar
snuverink_j committed
213
//     IpplTimings::startTimer(mapCreation_m);
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
214

snuverink_j's avatar
snuverink_j committed
215 216 217 218
    // TODO need a way to implement Initial dispersion Values
    //dispInitialVal[0][0]=  0.5069938765;
    //dispInitialVal[0][1]= -0.1681363086;
    OpalData::getInstance()->setInPrepState(false);
Philippe Ganz's avatar
Philippe Ganz committed
219

snuverink_j's avatar
snuverink_j committed
220
    track_m();
221

snuverink_j's avatar
snuverink_j committed
222 223 224
    // fMatrix_t sigMatrix = itsBunch_m->getSigmaMatrix();
    // linSigAnalyze(sigMatrix);
}
Philippe Ganz's avatar
Philippe Ganz committed
225

226

snuverink_j's avatar
snuverink_j committed
227
void ThickTracker::checkElementOrder_m() {
228

snuverink_j's avatar
snuverink_j committed
229 230 231
    // check order of beam line
    FieldList elements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
    beamline_t::const_iterator el = elements_m.cbegin();
232

snuverink_j's avatar
snuverink_j committed
233 234 235 236 237
    double currentEnd = zstart_m;
    for (FieldList::iterator it = elements.begin(); it != elements.end(); ++it) {
        double pos = it->getElement()->getElementPosition();
        // we have to take this length due to dipole --> arclength
        if ( currentEnd - pos > threshold_m ) {
238

snuverink_j's avatar
snuverink_j committed
239 240 241 242 243 244
            throw OpalException("ThickTracker::checkOrder_m()",
                                std::string("Elements are not in ascending order or overlap.") +
                                " Element Name: " + it->getElement()->getName() +
                                " ELEMEDGE position: " + std::to_string(pos) +
                                " Beamline end " + std::to_string(currentEnd) +
                                " element length " + std::to_string(std::get<2>(*el))) ;
frey_m's avatar
frey_m committed
245
        }
snuverink_j's avatar
snuverink_j committed
246 247 248 249 250
        currentEnd = pos + std::get<2>(*el);
        ++el;
        while (std::get<2>(*el) == 0) ++el; // skip zero-length elements (aka fringe fields)
    }
}
251

Philippe Ganz's avatar
Backup  
Philippe Ganz committed
252

snuverink_j's avatar
snuverink_j committed
253 254 255 256 257
///Fills undefined beam path with a Drift Space
/** \f[H_{Drift}= \frac{\delta}{\beta_0} -
* \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } \f]
*/
void ThickTracker::fillGaps_m() {
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
258

snuverink_j's avatar
snuverink_j committed
259
    beamline_t tmp;
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
260

snuverink_j's avatar
snuverink_j committed
261 262
    FieldList elements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
    beamline_t::const_iterator el = elements_m.cbegin();
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
263

snuverink_j's avatar
snuverink_j committed
264 265
    double currentEnd = zstart_m;
    for (FieldList::iterator it = elements.begin(); it != elements.end(); ++it) {
266

snuverink_j's avatar
snuverink_j committed
267
        tmp.push_back( *el );
268

snuverink_j's avatar
snuverink_j committed
269
        double pos = it->getElement()->getElementPosition();
Philippe Ganz's avatar
Philippe Ganz committed
270

snuverink_j's avatar
snuverink_j committed
271 272 273 274 275
        double length = std::abs(pos - currentEnd);
        if ( length  > threshold_m ) {
            //FIXME if gamma changes this is an error!!!
            double gamma = itsReference.getGamma();
            tmp.push_back(std::make_tuple(hamiltonian_m.drift(gamma), 1, length));
frey_m's avatar
frey_m committed
276
        }
277

snuverink_j's avatar
snuverink_j committed
278 279
        currentEnd = pos + std::get<2>(*el);
        ++el;
280
    }
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
281

snuverink_j's avatar
snuverink_j committed
282 283 284 285 286 287
    double length = std::abs(zstop_m - currentEnd);
    if ( length > threshold_m ) {
        //FIXME if gamma changes this is an error!!!
        double gamma = itsReference.getGamma();
        tmp.push_back(std::make_tuple(hamiltonian_m.drift(gamma), 1, length));
    }
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
288

snuverink_j's avatar
snuverink_j committed
289 290
    elements_m.swap(tmp);
}
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
291 292


snuverink_j's avatar
snuverink_j committed
293 294 295
void ThickTracker::track_m()
{
    IpplTimings::startTimer(mapTracking_m);
296

Philippe Ganz's avatar
Backup  
Philippe Ganz committed
297
    fMatrix_t tFMatrix;
snuverink_j's avatar
snuverink_j committed
298 299
    for (int i=0; i<6; i++){
        tFMatrix[i][i]=1.;
frey_m's avatar
frey_m committed
300
    }
301

snuverink_j's avatar
snuverink_j committed
302 303 304 305
    FMatrix<double, 1, 4> dispInitialVal;
    for (int i=2; i<4; i++){
        dispInitialVal[0][i]=0;
    }
306

snuverink_j's avatar
snuverink_j committed
307
    map_t transferMap;
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
308

snuverink_j's avatar
snuverink_j committed
309 310
    fMatrix_t refSigma;
    refSigma = (itsBunch_m->getSigmaMatrix());
311

snuverink_j's avatar
snuverink_j committed
312
    double spos = zstart_m;
313 314


snuverink_j's avatar
snuverink_j committed
315
    this->advanceDispersion_m(tFMatrix, dispInitialVal, spos);
316

317

snuverink_j's avatar
snuverink_j committed
318
    std::size_t step = 0;
319

snuverink_j's avatar
snuverink_j committed
320
    this->dump_m();
321

snuverink_j's avatar
snuverink_j committed
322 323 324
    //(1) Loop Beamline
    for(beamline_t::const_iterator it = elements_m.cbegin(); it != elements_m.end(); ++it) {
        //(2) Loop Slices
325

snuverink_j's avatar
snuverink_j committed
326 327 328
        const series_t& H          = std::get<0>(*it);
        const std::size_t& nSlices = std::get<1>(*it);
        const double& length       = std::get<2>(*it);
329

snuverink_j's avatar
snuverink_j committed
330
        const double ds = length / double(nSlices);
331

snuverink_j's avatar
snuverink_j committed
332
        map_t map = ExpMap(- H * ds, truncOrder_m);
333

snuverink_j's avatar
snuverink_j committed
334 335
        // global truncation order is truncOrder_m + 1
        map = map.truncate(truncOrder_m);
336

snuverink_j's avatar
snuverink_j committed
337
        for (std::size_t slice = 0; slice < nSlices; ++slice) {
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
338

snuverink_j's avatar
snuverink_j committed
339
            this->advanceParticles_m(map);
340

snuverink_j's avatar
snuverink_j committed
341
            this->concatenateMaps_m(map, transferMap);
342

snuverink_j's avatar
snuverink_j committed
343 344
            tFMatrix= transferMap.linearTerms();
            this->advanceDispersion_m(tFMatrix, dispInitialVal, spos);
345

snuverink_j's avatar
snuverink_j committed
346 347
            spos += ds;
            ++step;
348

snuverink_j's avatar
snuverink_j committed
349
            this->update_m(spos, step);
350

snuverink_j's avatar
snuverink_j committed
351
            this->dump_m();
352

snuverink_j's avatar
snuverink_j committed
353
            refSigma = itsBunch_m->getSigmaMatrix();
354

snuverink_j's avatar
snuverink_j committed
355 356 357
            //printPhaseShift(refSigma ,mapBeamLineit->elementMap.linearTerms(), N);
        }
    }
358 359


snuverink_j's avatar
snuverink_j committed
360
    this->write_m(transferMap);
361

snuverink_j's avatar
snuverink_j committed
362
    IpplTimings::stopTimer(mapTracking_m);
frey_m's avatar
frey_m committed
363
}
gsell's avatar
gsell committed
364 365


snuverink_j's avatar
snuverink_j committed
366 367 368 369 370 371 372 373
ThickTracker::particle_t
ThickTracker::particleToVector_m(const Vector_t& R,
                                 const Vector_t& P) const
{
    particle_t particle;
    for (int d = 0; d < 3; ++d) {
        particle[2 * d] = R(d);
        particle[2 *d + 1] = P(d);
374
    }
snuverink_j's avatar
snuverink_j committed
375
    return particle;
Philippe Ganz's avatar
Philippe Ganz committed
376
}
377

378

snuverink_j's avatar
snuverink_j committed
379 380 381 382 383 384 385 386
void ThickTracker::vectorToParticle_m(const particle_t& particle,
                                      Vector_t& R,
                                      Vector_t& P) const
{
    for (int d = 0; d < 3; ++d) {
        R(d) = particle[2 * d];
        P(d) = particle[2 *d + 1];
    }
387 388 389
}


snuverink_j's avatar
snuverink_j committed
390
void ThickTracker::write_m(const map_t& map) {
391

snuverink_j's avatar
snuverink_j committed
392
    if ( Ippl::myNode() == 0 ) {
393

snuverink_j's avatar
snuverink_j committed
394
        static bool first = true;
Philippe Ganz's avatar
Backup  
Philippe Ganz committed
395

snuverink_j's avatar
snuverink_j committed
396
        std::string fn = OpalData::getInstance()->getInputBasename() + ".map";
397

snuverink_j's avatar
snuverink_j committed
398
        std::ofstream out;
399

snuverink_j's avatar
snuverink_j committed
400 401 402 403 404 405
        if ( first ) {
            first = false;
            out.open(fn, std::ios::out);
        } else {
            out.open(fn, std::ios::app);
        }
406

snuverink_j's avatar
snuverink_j committed
407 408 409
        out << std::setprecision(16)
            << map
            << std::endl;
410

snuverink_j's avatar
snuverink_j committed
411 412 413
        out.close();
    }
}
414 415


snuverink_j's avatar
snuverink_j committed
416 417
void ThickTracker::advanceParticles_m(const map_t& map) {
    for (std::size_t ip = 0; ip < itsBunch_m->getLocalNum(); ++ip) {
418

snuverink_j's avatar
snuverink_j committed
419 420
        particle_t particle = this->particleToVector_m(itsBunch_m->R[ip],
                                                       itsBunch_m->P[ip]);
421

snuverink_j's avatar
snuverink_j committed
422
        this->updateParticle_m(particle, map);
423

snuverink_j's avatar
snuverink_j committed
424 425 426
        this->vectorToParticle_m(particle,
                                 itsBunch_m->R[ip],
                                 itsBunch_m->P[ip]);
427 428
    }

snuverink_j's avatar
snuverink_j committed
429 430 431
    // update reference particle
    particle_t particle = this->particleToVector_m(RefPartR_m,
                                                   RefPartP_m);
432

snuverink_j's avatar
snuverink_j committed
433
    this->updateParticle_m(particle, map);
434

snuverink_j's avatar
snuverink_j committed
435 436 437 438
    this->vectorToParticle_m(particle,
                             RefPartR_m,
                             RefPartP_m);
}
439 440


snuverink_j's avatar
snuverink_j committed
441 442 443 444
void ThickTracker::updateParticle_m(particle_t& particle,
                                    const map_t& map)
{
    double betagamma = itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma();
445

snuverink_j's avatar
snuverink_j committed
446 447 448 449
    //Units
    double pParticle= std::sqrt( particle[1] * particle[1] +
                                 particle[3] * particle[3] +
                                 particle[5] * particle[5]); // [beta gamma]
450

snuverink_j's avatar
snuverink_j committed
451 452
    particle[1] /= betagamma;
    particle[3] /= betagamma;
453

snuverink_j's avatar
snuverink_j committed
454 455
    particle[5] = std::sqrt( 1.0 + pParticle * pParticle ) / betagamma
                    - 1.0 / itsBunch_m->getInitialBeta();
456

snuverink_j's avatar
snuverink_j committed
457 458
    //Apply element map
    particle = map * particle;
459

snuverink_j's avatar
snuverink_j committed
460 461 462
    //Units back
    particle[1] *= betagamma;
    particle[3] *= betagamma;
463

snuverink_j's avatar
snuverink_j committed
464 465
    double tempGamma = (particle[5] + 1.0 / itsBunch_m->getInitialBeta())
            * betagamma ;
466

snuverink_j's avatar
snuverink_j committed
467
    pParticle =  std::sqrt( tempGamma * tempGamma -1.0);
468

snuverink_j's avatar
snuverink_j committed
469 470 471 472
    particle[5] = std::sqrt(pParticle * pParticle -
                            particle[1] * particle[1] -
                            particle[3] * particle[3] );
}
473 474


snuverink_j's avatar
snuverink_j committed
475 476 477 478 479
void ThickTracker::concatenateMaps_m(const map_t& x, map_t& y) {
    IpplTimings::startTimer(mapCombination_m);
    y = x * y;
    y = y.truncate(truncOrder_m);
    IpplTimings::stopTimer(mapCombination_m);
480 481 482
}


snuverink_j's avatar
snuverink_j committed
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
///Writes Dispersion in X and Y plane
/** used formula:
 * \f[\begin{pmatrix}\eta_{x} \\ \eta_{p_x}\end{pmatrix}_{s_1} =
 *   \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \end{pmatrix}
 *  \cdot
 *  \begin{pmatrix} \eta_{x} \\ \eta_{p_x} \end{pmatrix}_{s_0}
 *  +
 *  \begin{pmatrix} R_{16} \\ R_{26} \end{pmatrix}\f]
 */
void ThickTracker::advanceDispersion_m(fMatrix_t tempMatrix,
                                       FMatrix<double, 1, 4> initialVal,
                                       double pos)
{
    if ( Ippl::myNode() == 0 ) {
        static bool first = true;
498

snuverink_j's avatar
snuverink_j committed
499 500
        std::ofstream out;
        std::string fn = OpalData::getInstance()->getInputBasename() + ".dispersion";
501

snuverink_j's avatar
snuverink_j committed
502 503 504 505 506 507
        if ( first ) {
            first = false;
            out.open(fn, std::ios::out);
        } else {
            out.open(fn, std::ios::app);
        }
508

snuverink_j's avatar
snuverink_j committed
509
        out << std::setprecision(16);
510 511


snuverink_j's avatar
snuverink_j committed
512 513 514
        FMatrix<double, 2, 2> subx, suby;
        FMatrix<double, 2, 1> subdx, subdy;
        FMatrix<double, 2, 1> dxi, dyi, dx, dy;
515

snuverink_j's avatar
snuverink_j committed
516 517 518
        for (int i=0; i<2; i++){
            dxi[i][0]= initialVal[0][i]; //
            dyi[i][0]= initialVal[0][i+2];
519

snuverink_j's avatar
snuverink_j committed
520 521
            subdx[i][0]=tempMatrix[i][5]*itsBunch_m->getInitialBeta();
            subdy[i][0]=tempMatrix[i+2][5]*itsBunch_m->getInitialBeta();
522

snuverink_j's avatar
snuverink_j committed
523
            for (int j=0; j<2; j++){
524

snuverink_j's avatar
snuverink_j committed
525 526
                subx[i][j]=tempMatrix[i][j];
                suby[i][j]=tempMatrix[i+2][j+2];
527

snuverink_j's avatar
snuverink_j committed
528 529
            }
        }
530

snuverink_j's avatar
snuverink_j committed
531 532
        dx= subx*dxi + subdx;
        dy= suby*dyi + subdy;
533 534


Philippe Ganz's avatar
Philippe Ganz committed
535

snuverink_j's avatar
snuverink_j committed
536 537 538
        out <<pos<< "\t" << dx[0][0] << "\t" << dx[0][1] << "\t" << dy[0][0] << "\t" << dy[0][1] << std::endl;
    }
}
Philippe Ganz's avatar
Philippe Ganz committed
539 540


snuverink_j's avatar
snuverink_j committed
541
void ThickTracker::dump_m() {
Philippe Ganz's avatar
Philippe Ganz committed
542

snuverink_j's avatar
snuverink_j committed
543
    if ( itsBunch_m->getTotalNum() == 0 )
Philippe Ganz's avatar
Philippe Ganz committed
544
        return;
545

snuverink_j's avatar
snuverink_j committed
546
    Inform msg("ThickTracker", *gmsg);
Philippe Ganz's avatar
Philippe Ganz committed
547

snuverink_j's avatar
snuverink_j committed
548
    msg << *itsBunch_m << endl;
Philippe Ganz's avatar
Philippe Ganz committed
549

snuverink_j's avatar
snuverink_j committed
550
    const std::size_t step = itsBunch_m->getGlobalTrackStep();
Philippe Ganz's avatar
Philippe Ganz committed
551

snuverink_j's avatar
snuverink_j committed
552 553
    bool psDump   = ((step + 1) % Options::psDumpFreq == 0);
    bool statDump = ((step + 1) % Options::statDumpFreq == 0);
Philippe Ganz's avatar
Philippe Ganz committed
554 555 556 557

    // Sample fields at (xmin, ymin, zmin), (xmax, ymax, zmax) and the centroid location. We
    // are sampling the electric and magnetic fields at the back, front and
    // center of the beam.
snuverink_j's avatar
snuverink_j committed
558 559

    Vector_t FDext[2];  // FDext = {BHead, EHead, BRef, ERef, BTail, ETail}.
Philippe Ganz's avatar
Philippe Ganz committed
560 561

    if (psDump || statDump) {
snuverink_j's avatar
snuverink_j committed
562 563 564 565 566 567

        Vector_t externalE, externalB;

        Vector_t rmin, rmax;
        itsBunch_m->get_bounds(rmin, rmax);

Philippe Ganz's avatar
Philippe Ganz committed
568 569 570 571 572 573 574 575 576 577 578 579
        externalB = Vector_t(0.0);
        externalE = Vector_t(0.0);
        itsOpalBeamline_m.getFieldAt(referenceToLabCSTrafo_m.transformTo(RefPartR_m),
                                     referenceToLabCSTrafo_m.rotateTo(RefPartP_m),
                                     itsBunch_m->getT() - 0.5 * itsBunch_m->getdT(),
                                     externalE,
                                     externalB);
        FDext[0] = referenceToLabCSTrafo_m.rotateFrom(externalB);
        FDext[1] = referenceToLabCSTrafo_m.rotateFrom(externalE * 1e-6);
    }


snuverink_j's avatar
snuverink_j committed
580
    if ( psDump ) {
frey_m's avatar
frey_m committed
581
        itsDataSink_m->dumpH5(itsBunch_m, FDext);
snuverink_j's avatar
snuverink_j committed
582
    }
583

snuverink_j's avatar
snuverink_j committed
584 585
    if ( statDump ) {
        std::vector<std::pair<std::string, unsigned int> > collimatorLosses;
frey_m's avatar
frey_m committed
586
        itsDataSink_m->dumpSDDS(itsBunch_m, FDext, collimatorLosses);
Philippe Ganz's avatar
Philippe Ganz committed
587
    }
snuverink_j's avatar
snuverink_j committed
588
}
Philippe Ganz's avatar
Philippe Ganz committed
589 590


snuverink_j's avatar
snuverink_j committed
591 592 593 594 595 596 597 598 599 600
void ThickTracker::update_m(const double& spos,
                            const std::size_t& step)
{
    // update dt and t
    double ds = spos - itsBunch_m->get_sPos();
    double gamma = itsBunch_m->get_gamma();
    double beta  = std::sqrt(1.0 - 1.0 / (gamma * gamma) );
    double dt = ds / Physics::c / beta;
    itsBunch_m->setdT(dt);
    itsBunch_m->setT(itsBunch_m->getT() + dt);
Philippe Ganz's avatar
Philippe Ganz committed
601

602 603
    const unsigned int localNum = itsBunch_m->getLocalNum();
    for (unsigned int i = 0; i < localNum; ++i) {
snuverink_j's avatar
snuverink_j committed
604
        itsBunch_m->dt[i] = dt;
605
    }
snuverink_j's avatar
snuverink_j committed
606 607 608 609 610 611


    itsBunch_m->set_sPos(spos);
    itsBunch_m->setGlobalTrackStep(step);
    itsBunch_m->calcBeamParameters();
    itsBunch_m->calcEMean();
612
}