Distribution.cpp 186 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12
// ------------------------------------------------------------------------
// $RCSfile: Distribution.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.3.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Distribution
//   The class for the OPAL Distribution command.
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
13 14

#include "Distribution/Distribution.h"
15
#include "Distribution/SigmaGenerator.h"
16
#include "AbsBeamline/SpecificElementVisitor.h"
17

18 19 20 21 22 23
#include <cmath>
#include <cfloat>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
24
#include <numeric>
25

26 27
#include "AbstractObjects/Expressions.h"
#include "Attributes/Attributes.h"
kraus's avatar
kraus committed
28
#include "Utilities/OpalOptions.h"
29 30 31 32 33 34
#include "Utilities/Options.h"
#include "halton1d_sequence.hh"
#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBunch.h"
#include "Algorithms/PartBins.h"
#include "Algorithms/bet/EnvelopeBunch.h"
35
#include "Structure/Beam.h"
36
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
37 38 39
#include "Algorithms/PartBinsCyc.h"
#include "BasicActions/Option.h"
#include "Distribution/LaserProfile.h"
40
#include "Elements/OpalBeamline.h"
41
#include "AbstractObjects/BeamSequence.h"
42 43 44 45

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_erf.h>
46 47
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
gsell's avatar
gsell committed
48

49 50
#include "MagneticField.h"

gsell's avatar
gsell committed
51 52
extern Inform *gmsg;

kraus's avatar
kraus committed
53
#define DISTDBG1
adelmann's avatar
Cleanup  
adelmann committed
54
#define noDISTDBG2
55

56
#define COMPLAINONFAILURE(rc)     complainOnFailure(rc, __FILE__, __LINE__);
57

58 59 60 61 62
inline
void complainOnFailure(int rc, const std::string &file, int line) {
    if(rc != H5_SUCCESS)
        ERRORMSG("H5 rc= " << rc << " in " << file << " @ line " << line << endl);
}
63

64 65 66 67 68 69 70 71
SymTenzor<double, 6> getUnit6x6() {
    SymTenzor<double, 6> unit6x6;
    for (unsigned int i = 0; i < 6u; ++ i) {
        unit6x6(i,i) = 1.0;
    }
    return unit6x6;
}

gsell's avatar
gsell committed
72 73 74 75
//
// Class Distribution
// ------------------------------------------------------------------------

76 77 78
namespace AttributesT
{
    enum AttributesT {
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
        DISTRIBUTION,
        FNAME,
        WRITETOFILE,
        WEIGHT,
        INPUTMOUNITS,
        EMITTED,
        EMISSIONSTEPS,
        EMISSIONMODEL,
        EKIN,
        ELASER,
        W,
        FE,
        CATHTEMP,
        NBIN,
        XMULT,
        YMULT,
        ZMULT,
        TMULT,
        PXMULT,
        PYMULT,
        PZMULT,
        OFFSETX,
        OFFSETY,
        OFFSETZ,
        OFFSETT,
        OFFSETPX,
        OFFSETPY,
        OFFSETPZ,
        SIGMAX,
        SIGMAY,
        SIGMAR,
        SIGMAZ,
        SIGMAT,
        TPULSEFWHM,
        TRISE,
        TFALL,
        SIGMAPX,
        SIGMAPY,
        SIGMAPZ,
        MX,
        MY,
        MZ,
        MT,
        CUTOFFX,
        CUTOFFY,
        CUTOFFR,
        CUTOFFLONG,
        CUTOFFPX,
        CUTOFFPY,
        CUTOFFPZ,
        FTOSCAMPLITUDE,
        FTOSCPERIODS,
        R,                          // the correlation matrix (a la transport)
        CORRX,
        CORRY,
        CORRZ,
        CORRT,
        R51,
        R52,
        R61,
        R62,
        LASERPROFFN,
        IMAGENAME,
        INTENSITYCUT,
        NPDARKCUR,
        INWARDMARGIN,
        EINITHR,
        FNA,
        FNB,
        FNY,
        FNVYZERO,
        FNVYSECOND,
        FNPHIW,
        FNBETA,
        FNFIELDTHR,
        FNMAXEMI,
        SECONDARYFLAG,
        NEMISSIONMODE,
        VSEYZERO,                   // sey_0 in Vaughn's model.
        VEZERO,                     // Energy related to sey_0 in Vaughan's model.
        VSEYMAX,                    // sey max in Vaughan's model.
        VEMAX,                      // Emax in Vaughan's model.
        VKENERGY,                   // Fitting parameter denotes the roughness of
        // surface for impact energy in Vaughn's model.
        VKTHETA,                    // Fitting parameter denotes the roughness of
        // surface for impact angle in Vaughn's model.
        VVTHERMAL,                  // Thermal velocity of Maxwellian distribution
        // of secondaries in Vaughan's model.
        VW,
        SURFMATERIAL,               // Add material type, currently 0 for copper
        // and 1 for stainless steel.
        EX,                         // below is for the matched distribution
        EY,
        ET,
        MAGSYM,                     // number of sector magnets
        LINE,
        FMAPFN,
        RESIDUUM,
        MAXSTEPSCO,
        MAXSTEPSSI,
        ORDERMAPS,
        E2,
        SIZE
gsell's avatar
gsell committed
182 183 184
    };
}

185 186 187
namespace LegacyAttributesT
{
    enum LegacyAttributesT {
188
        // DESCRIPTION OF THE DISTRIBUTION:
189
        DEBIN = AttributesT::SIZE,
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
        SBIN,
        TEMISSION,
        SIGLASER,
        AG,
        SIGMAPT,
        TRANSVCUTOFF,
        CUTOFF,
        Z,
        T,
        PT,
        ALPHAX,
        ALPHAY,
        BETAX,
        BETAY,
        DX,
        DDX,
        DY,
        DDY,
        SIZE
    };
210 211
}

gsell's avatar
gsell committed
212
Distribution::Distribution():
213 214
    Definition( LegacyAttributesT::SIZE, "DISTRIBUTION",
                "The DISTRIBUTION statement defines data for the 6D particle distribution."),
215
    distrTypeT_m(DistrTypeT::NODIST),
216
    numberOfDistributions_m(1),
217 218 219 220 221 222 223 224 225 226 227 228 229
    emitting_m(false),
    scan_m(false),
    emissionModel_m(EmissionModelT::NONE),
    tEmission_m(0.0),
    tBin_m(0.0),
    currentEmissionTime_m(0.0),
    currentEnergyBin_m(0.0),
    currentSampleBin_m(0.0),
    numberOfEnergyBins_m(0),
    numberOfSampleBins_m(0),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
    randGenEmit_m(NULL),
230
    pTotThermal_m(0.0),
kraus's avatar
kraus committed
231
    pmean_m(0.0),
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
    cathodeWorkFunc_m(0.0),
    laserEnergy_m(0.0),
    cathodeFermiEnergy_m(0.0),
    cathodeTemp_m(0.0),
    emitEnergyUpperLimit_m(0.0),
    inputMoUnits_m(InputMomentumUnitsT::NONE),
    sigmaTRise_m(0.0),
    sigmaTFall_m(0.0),
    tPulseLengthFWHM_m(0.0),
    correlationMatrix_m(getUnit6x6()),
    laserProfileFileName_m(""),
    laserImageName_m(""),
    laserIntensityCut_m(0.0),
    laserProfile_m(NULL),
    darkCurrentParts_m(0),
    darkInwardMargin_m(0.0),
    eInitThreshold_m(0.0),
    workFunction_m(0.0),
    fieldEnhancement_m(0.0),
    fieldThrFN_m(0.0),
    maxFN_m(0),
    paraFNA_m(0.0),
    paraFNB_m(0.0),
    paraFNY_m(0.0),
    paraFNVYSe_m(0.0),
    paraFNVYZe_m(0.0),
    secondaryFlag_m(0),
    ppVw_m(0.0),
    vVThermal_m(0.0),
    referencePz_m(0.0),
    referenceZ_m(0.0),
    avrgpz_m(0.0),
    I_m(0.0),
    E_m(0.0),
    bega_m(0.0),
    M_m(0.0)
268
{
269 270 271 272
    SetAttributes();

    Distribution *defaultDistribution = clone("UNNAMED_Distribution");
    defaultDistribution->builtin = true;
gsell's avatar
gsell committed
273 274

    try {
275
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
276
    } catch(...) {
277 278 279 280
        delete defaultDistribution;
    }

    SetFieldEmissionParameters();
gsell's avatar
gsell committed
281 282 283 284 285 286 287
}
/**
 *
 *
 * @param name
 * @param parent
 */
288
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
289
    Definition(name, parent),
290 291
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
292
    numberOfDistributions_m(parent->numberOfDistributions_m),
293 294 295 296 297 298
    emitting_m(parent->emitting_m),
    scan_m(parent->scan_m),
    particleRefData_m(parent->particleRefData_m),
    addedDistributions_m(parent->addedDistributions_m),
    particlesPerDist_m(parent->particlesPerDist_m),
    emissionModel_m(parent->emissionModel_m),
gsell's avatar
gsell committed
299
    tEmission_m(parent->tEmission_m),
300 301 302 303 304 305 306 307 308 309
    tBin_m(parent->tBin_m),
    currentEmissionTime_m(parent->currentEmissionTime_m),
    currentEnergyBin_m(parent->currentEmissionTime_m),
    currentSampleBin_m(parent->currentSampleBin_m),
    numberOfEnergyBins_m(parent->numberOfEnergyBins_m),
    numberOfSampleBins_m(parent->numberOfSampleBins_m),
    energyBins_m(NULL),
    energyBinHist_m(NULL),
    randGenEmit_m(parent->randGenEmit_m),
    pTotThermal_m(parent->pTotThermal_m),
kraus's avatar
kraus committed
310
    pmean_m(parent->pmean_m),
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
    cathodeWorkFunc_m(parent->cathodeWorkFunc_m),
    laserEnergy_m(parent->laserEnergy_m),
    cathodeFermiEnergy_m(parent->cathodeFermiEnergy_m),
    cathodeTemp_m(parent->cathodeTemp_m),
    emitEnergyUpperLimit_m(parent->emitEnergyUpperLimit_m),
    xDist_m(parent->xDist_m),
    pxDist_m(parent->pxDist_m),
    yDist_m(parent->yDist_m),
    pyDist_m(parent->pyDist_m),
    tOrZDist_m(parent->tOrZDist_m),
    pzDist_m(parent->pzDist_m),
    xWrite_m(parent->xWrite_m),
    pxWrite_m(parent->pxWrite_m),
    yWrite_m(parent->yWrite_m),
    pyWrite_m(parent->pyWrite_m),
    tOrZWrite_m(parent->tOrZWrite_m),
    pzWrite_m(parent->pzWrite_m),
328
    avrgpz_m(parent->avrgpz_m),
329 330 331 332 333 334 335 336
    inputMoUnits_m(parent->inputMoUnits_m),
    sigmaTRise_m(parent->sigmaTRise_m),
    sigmaTFall_m(parent->sigmaTFall_m),
    tPulseLengthFWHM_m(parent->tPulseLengthFWHM_m),
    sigmaR_m(parent->sigmaR_m),
    sigmaP_m(parent->sigmaP_m),
    cutoffR_m(parent->cutoffR_m),
    cutoffP_m(parent->cutoffP_m),
337
    correlationMatrix_m(parent->correlationMatrix_m),
338 339 340 341
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
342 343 344 345 346 347 348 349 350 351 352 353
    darkCurrentParts_m(parent->darkCurrentParts_m),
    darkInwardMargin_m(parent->darkInwardMargin_m),
    eInitThreshold_m(parent->eInitThreshold_m),
    workFunction_m(parent->workFunction_m),
    fieldEnhancement_m(parent->fieldEnhancement_m),
    fieldThrFN_m(parent->fieldThrFN_m),
    maxFN_m(parent->maxFN_m),
    paraFNA_m(parent-> paraFNA_m),
    paraFNB_m(parent-> paraFNB_m),
    paraFNY_m(parent-> paraFNY_m),
    paraFNVYSe_m(parent-> paraFNVYSe_m),
    paraFNVYZe_m(parent-> paraFNVYZe_m),
354 355 356 357 358 359 360
    secondaryFlag_m(parent->secondaryFlag_m),
    ppVw_m(parent->ppVw_m),
    vVThermal_m(parent->vVThermal_m),
    tRise_m(parent->tRise_m),
    tFall_m(parent->tFall_m),
    sigmaRise_m(parent->sigmaRise_m),
    sigmaFall_m(parent->sigmaFall_m),
361 362
    cutoff_m(parent->cutoff_m),
    I_m(parent->I_m),
adelmann's avatar
adelmann committed
363 364 365
    E_m(parent->E_m),
    bega_m(parent->bega_m),
    M_m(parent->M_m)
366
{
gsell's avatar
gsell committed
367 368 369 370 371 372 373
}

Distribution::~Distribution() {

    if((Ippl::getNodes() == 1) && (os_m.is_open()))
        os_m.close();

374 375 376
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
377 378
    }

379 380 381 382 383 384 385 386 387 388 389 390 391 392
    if (energyBinHist_m != NULL) {
        gsl_histogram_free(energyBinHist_m);
        energyBinHist_m = NULL;
    }

    if (randGenEmit_m != NULL) {
        delete randGenEmit_m;
        randGenEmit_m = NULL;
    }

    if(laserProfile_m) {
        delete laserProfile_m;
        laserProfile_m = NULL;
    }
gsell's avatar
gsell committed
393 394

}
395
/*
396
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
397
  for(int i=0; i<M.size1(); ++i) {
398 399 400 401 402
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
403 404
  }
*/
gsell's avatar
gsell committed
405 406

/**
407 408
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
409 410 411
 * @param
 * @param
 * @param
412
 */
413
void Distribution::WriteToFile() {
414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429
    /*
      if(Ippl::getNodes() == 1) {
      if(os_m.is_open()) {
      ;
      } else {
      *gmsg << " Write distribution to file data/dist.dat" << endl;
      std::string file("data/dist.dat");
      os_m.open(file.c_str());
      if(os_m.bad()) {
      *gmsg << "Unable to open output file " <<  file << endl;
      }
      os_m << "# x y ti px py pz "  << std::endl;
      os_m.close();
      }
      }
    */
430 431
}

432 433 434
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
435 436
}

437 438 439
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
440

441 442
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
443

444 445
void Distribution::update() {
}
adelmann's avatar
adelmann committed
446

447
void Distribution::Create(size_t &numberOfParticles, double massIneV) {
adelmann's avatar
adelmann committed
448

449
    SetFieldEmissionParameters();
adelmann's avatar
adelmann committed
450

451
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
452

453
    case DistrTypeT::MATCHEDGAUSS:
454
        CreateMatchedGaussDistribution(numberOfParticles, massIneV);
455
        break;
456 457
    case DistrTypeT::FROMFILE:
        CreateDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
458
        break;
459 460
    case DistrTypeT::GAUSS:
        CreateDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
461
        break;
462 463
    case DistrTypeT::BINOMIAL:
        CreateDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
464
        break;
465 466
    case DistrTypeT::FLATTOP:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
467
        break;
468 469
    case DistrTypeT::GUNGAUSSFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
470
        break;
471 472
    case DistrTypeT::ASTRAFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
473
        break;
474 475 476 477
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
gsell's avatar
gsell committed
478

479
    // AAA Scale and shift coordinates according to distribution input.
480 481 482
    ScaleDistCoordinates();
    ShiftDistCoordinates(massIneV);
}
gsell's avatar
gsell committed
483

484
void  Distribution::CreatePriPart(PartBunch *beam, BoundaryGeometry &bg) {
gsell's avatar
gsell committed
485

486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
    if(Options::ppdebug) {  // This is Parallel Plate Benchmark.
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
        double vw = this->GetVw();
        double vt = this->GetvVThermal();
        double f_max = vw / vt * exp(-0.5);
        double test_a = vt / vw;
        double test_asq = test_a * test_a;
        size_t count = 0;
        size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
        size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
        if(Ippl::myNode() == 0)
            N_mean += N_extra;
        if(bg.getN() != 0) {
            for(size_t i = 0; i < bg.getN(); i++) {
                if(pc == Ippl::myNode()) {
                    if(count < N_mean) {
                        /*==============Parallel Plate Benchmark=====================================*/
                        double test_s = 1;
                        double f_x = 0;
                        double test_x = 0;
                        while(test_s > f_x) {
                            test_s = IpplRandom();
                            test_s *= f_max;
                            test_x = IpplRandom();
                            test_x *= 10 * test_a; //range for normalized emission speed(0,10*test_a);
                            f_x = test_x / test_asq * exp(-test_x * test_x / 2 / test_asq);
                        }
                        double v_emi = test_x * vw;
gsell's avatar
gsell committed
515

516 517 518 519 520 521 522
                        double betaemit = v_emi / Physics::c;
                        double betagamma = betaemit / sqrt(1 - betaemit * betaemit);
                        /*============================================================================ */
                        beam->create(1);
                        if(pc != 0) {
                            beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(Ippl::myNode() * N_mean + count);
gsell's avatar
gsell committed
523
                        } else {
524 525
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
526
                        }
527 528 529 530 531 532 533 534 535
                        beam->Bin[lowMark + count] = 0;
                        beam->PType[lowMark + count] = 0; // create primary particle bunch;
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        beam->LastSection[lowMark + count] = 0;
                        beam->Ef[lowMark + count] = Vector_t(0.0);
                        beam->Bf[lowMark + count] = Vector_t(0.0);
                        beam->dt[lowMark + count] = beam->getdT();
                        count ++;
gsell's avatar
gsell committed
536 537
                    }
                }
538 539 540
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
541
            }
542 543 544
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
545
        }
546
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
547

548
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
549

550 551 552 553 554
        int pc = 0;
        size_t lowMark = beam->getLocalNum();
        size_t count = 0;
        size_t N_mean = static_cast<size_t>(floor(bg.getN() / Ippl::getNodes()));
        size_t N_extra = static_cast<size_t>(bg.getN() - N_mean * Ippl::getNodes());
gsell's avatar
gsell committed
555

556 557 558 559
        if(Ippl::myNode() == 0)
            N_mean += N_extra;
        if(bg.getN() != 0) {
            for(size_t i = 0; i < bg.getN(); i++) {
gsell's avatar
gsell committed
560

561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577
                if(pc == Ippl::myNode()) {
                    if(count < N_mean) {
                        beam->create(1);
                        if(pc != 0)
                            beam->R[lowMark + count] = bg.getCooridinate(Ippl::myNode() * N_mean + count + N_extra); // node 0 will emit the particle with coordinate ID from 0 to N_mean+N_extra, so other nodes should shift to node_number*N_mean+N_extra
                        else
                            beam->R[lowMark + count] = bg.getCooridinate(count); // for node0 the particle number N_mean =  N_mean + N_extra
                        beam->P[lowMark + count] = Vector_t(0.0);
                        beam->Bin[lowMark + count] = 0;
                        beam->PType[lowMark + count] = 0; // create primary particle bunch.
                        beam->TriID[lowMark + count] = 0;
                        beam->Q[lowMark + count] = beam->getChargePerParticle();
                        beam->LastSection[lowMark + count] = 0;
                        beam->Ef[lowMark + count] = Vector_t(0.0);
                        beam->Bf[lowMark + count] = Vector_t(0.0);
                        beam->dt[lowMark + count] = beam->getdT();
                        count++;
gsell's avatar
gsell committed
578 579

                    }
580

gsell's avatar
gsell committed
581
                }
582 583 584
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
585

586
            }
587

588 589 590 591 592 593
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
594

595 596 597
void Distribution::DoRestartOpalT(PartBunch &beam, size_t Np, int restartStep) {
    h5_file_t *H5file;
    std::string fn;
598

599
    IpplTimings::startTimer(beam.distrReload_m);
600

601
    //        beam.setTEmission(Attributes::getReal(itsAttr[ LegacyAttributesT::TEMISSION]));
602
    fn = OpalData::getInstance()->getInputBasename() + std::string(".h5");
gsell's avatar
gsell committed
603 604 605 606

#ifdef PARALLEL_IO
    H5file = H5OpenFile(fn.c_str(), H5_O_RDONLY, Ippl::getComm());
#else
607
    H5file = H5OpenFile(fn.c_str(), H5_O_RDONLY, 0);
gsell's avatar
gsell committed
608 609
#endif

610
    if(H5file == (void*)H5_ERR) {
611 612
	throw OpalException("Distribution::DoRestartOpalT",
                            "could not open file '" + fn + "'");
gsell's avatar
gsell committed
613 614 615
    }

    if(restartStep == -1) {
616
        restartStep = H5GetNumSteps(H5file) - 1 ;
gsell's avatar
gsell committed
617 618 619
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
620 621
            throw OpalException("Distribution::DoRestartOpalT",
                                "can't append to the file '" + fn + "'");
gsell's avatar
gsell committed
622 623 624
        }
    }

625 626
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));

gsell's avatar
gsell committed
627 628
    int N = (int)H5PartGetNumParticles(H5file);

629 630 631 632 633 634 635 636 637
    int numberOfParticlesPerNode = (int) floor((double) N / Ippl::getNodes());
    long long starti = Ippl::myNode() * numberOfParticlesPerNode;
    long long endi = starti + numberOfParticlesPerNode - 1;

    // In case we miss some particles we add them at the end on the last core
    if(Ippl::myNode() == (Ippl::getNodes() - 1)) {
        if(Ippl::getNodes()*numberOfParticlesPerNode < N)
            endi += (N - (Ippl::getNodes() * numberOfParticlesPerNode));
    }
gsell's avatar
gsell committed
638

639
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
gsell's avatar
gsell committed
640
    N = (int)H5PartGetNumParticles(H5file);
641
    assert(N >= 0);
gsell's avatar
gsell committed
642 643

    double actualT;
644
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
gsell's avatar
gsell committed
645
    beam.setT(actualT);
646

gsell's avatar
gsell committed
647
    double dPhiGlobal;
648
    COMPLAINONFAILURE(H5ReadFileAttribFloat64(H5file, "dPhiGlobal", &dPhiGlobal));
gsell's avatar
gsell committed
649 650
    OpalData::getInstance()->setGlobalPhaseShift(dPhiGlobal);

651
    h5_int64_t ltstep;
652
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
653 654 655
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
656
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep", &gtstep));
657
    beam.setGlobalTrackStep((long long)gtstep);
658

659 660 661 662 663 664 665 666
    Vektor<double, 3> RefPartR;
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "RefPartR", (h5_float64_t *)&RefPartR));
    beam.RefPart_R = RefPartR;

    Vektor<double, 3> RefPartP;
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "RefPartP", (h5_float64_t *)&RefPartP));
    beam.RefPart_P = RefPartP;

667
    std::unique_ptr<char[]> varray(new char[(N)*sizeof(double)]);
668
    h5_float64_t *farray = reinterpret_cast<h5_float64_t *>(varray.get());
669
    h5_int64_t *larray = reinterpret_cast<h5_int64_t *>(varray.get());
gsell's avatar
gsell committed
670

671 672
    beam.create(N);

673
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
gsell's avatar
gsell committed
674
    for(unsigned long int n = 0; n < (unsigned int) N; ++n) {
675 676
        beam.R[n](0) = farray[n];
        beam.Bin[n] = 0; // not initialized
gsell's avatar
gsell committed
677
    }
678
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
gsell's avatar
gsell committed
679
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
680
        beam.R[n](1) = farray[n];
gsell's avatar
gsell committed
681

682
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
gsell's avatar
gsell committed
683
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
684
        beam.R[n](2) = farray[n];
gsell's avatar
gsell committed
685

686
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
gsell's avatar
gsell committed
687
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
688
        beam.P[n](0) = farray[n];
gsell's avatar
gsell committed
689

690
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "py", farray));
gsell's avatar
gsell committed
691
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
692
        beam.P[n](1) = farray[n];
gsell's avatar
gsell committed
693

694
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "pz", farray));
gsell's avatar
gsell committed
695
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
696
        beam.P[n](2) = farray[n];
gsell's avatar
gsell committed
697

698
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "q", farray));
gsell's avatar
gsell committed
699
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
700
        beam.Q[n] = farray[n];
gsell's avatar
gsell committed
701

702
    COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "lastsection", larray));
gsell's avatar
gsell committed
703 704 705 706
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
        beam.LastSection[n] = (short) larray[n];

    Ippl::Comm->barrier();
707
    COMPLAINONFAILURE(H5CloseFile(H5file));
708 709 710

    beam.boundp();

gsell's avatar
gsell committed
711 712
    IpplTimings::stopTimer(beam.distrReload_m);

713 714 715 716
    *gmsg << "Total number of particles in the h5 file = " << N << " NPerBunch= " << beam.getTotalNum()
          << " Global step " << gtstep << " Local step " << ltstep << endl
          << " restart step= " << restartStep << " time of restart = " << actualT
          << " phishift= " << OpalData::getInstance()->getGlobalPhaseShift() << endl;
gsell's avatar
gsell committed
717 718
}

719
void Distribution::DoRestartOpalCycl(PartBunch &beam, size_t Np, int restartStep, const int specifiedNumBunch) {
adelmann's avatar
adelmann committed
720

721 722 723 724 725 726 727 728 729 730 731 732 733
    h5_int64_t rc;
    IpplTimings::startTimer(beam.distrReload_m);
    *gmsg << "---------------- Start reading hdf5 file----------------" << endl;
    h5_file_t *H5file;

    std::string fn = OpalData::getInstance()->getInputBasename() + std::string(".h5");

#ifdef PARALLEL_IO
    H5file = H5OpenFile(fn.c_str(), H5_O_RDONLY, Ippl::getComm());
#else
    H5file = H5OpenFile(fn.c_str(), H5_O_RDONLY, 0);
#endif

734
    if(H5file == (void*)H5_ERR) {
735 736
        throw OpalException("Distribution::DoRestartOpalCycl",
                            "could not open file '" + fn + "'");
737 738 739 740 741 742 743
    }

    if(restartStep == -1) {
        restartStep = H5GetNumSteps(H5file) - 1;
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
744 745
            throw OpalException("Distribution::DoRestartOpalCycl",
                                "can't append to the file '" + fn + "'");
746 747 748
        }
    }

kraus's avatar
kraus committed
749
    *gmsg << "Restart from hdf5 format file " << fn
750
          << ", read phase space data of DumpStep " << (int)restartStep << endl;
751

752
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
753 754 755 756 757 758 759 760 761 762 763
    const int globalN = (int)H5PartGetNumParticles(H5file);

    *gmsg << "total number of particles = " << globalN << endl;

    int numberOfParticlesPerNode = (int) floor((double) globalN / Ippl::getNodes());
    long long starti = Ippl::myNode() * numberOfParticlesPerNode;
    long long endi = 0;

    if(Ippl::myNode() == Ippl::getNodes() - 1)
        endi = -1;
    else
764
        endi = starti + numberOfParticlesPerNode - 1;
kraus's avatar
kraus committed
765

766
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
767 768 769 770
    const int localN = (int)H5PartGetNumParticles(H5file);
    assert(localN >= 0);

    h5_int64_t ltstep;
771
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
772 773 774
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
775
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep", &gtstep));
776 777 778
    beam.setGlobalTrackStep((long long)gtstep);

    char opalFlavour[128];
779
    COMPLAINONFAILURE(H5ReadStepAttribString(H5file, "OPAL_flavour", opalFlavour));
780

781
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPR",&referencePr_m));
782

783
    referencePt_m = 0.0;
784
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPT",&referencePt_m));
785

786
    referencePz_m = 0.0;
787
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPZ",&referencePz_m));
788

789
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFR",&referenceR_m));
790

791
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFTHETA",&referenceTheta_m));
792

793
    referenceZ_m = 0.0;
794
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFZ",&referenceZ_m));
795

796
    double meanE;
797
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ENERGY", &meanE));
798

adelmann's avatar
adelmann committed
799
    double pathLength;;
800
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "SPOS", &pathLength));
adelmann's avatar
adelmann committed
801 802 803

    beam.setLPath(pathLength);

804 805 806
    *gmsg << "* Restart Energy = " << meanE << " (MeV), Path lenght = " << pathLength << " (m)" <<  endl;

    // beam.getM() in eV, meanE in MeV had to change 1e3 to 1e6 -DW
807 808
    double ga = 1 + meanE / beam.getM() * 1.0e6;
    double be = sqrt(1.0 - (1.0 / (ga * ga)));
809

810
    bega_m = be * ga;
811

812
    *gmsg << "* Gamma = " << ga << ", Beta = " << be << endl;
813

814
    std::unique_ptr<char[]> varray(new char[(localN)*sizeof(double)]);
815

816
    double *farray = reinterpret_cast<double *>(varray.get());
817

818 819 820 821 822
    h5_int64_t *larray = reinterpret_cast<h5_int64_t *>(varray.get());

    beam.create(localN);

    if(strcmp(opalFlavour, "opal-t") == 0) {
823

824 825 826 827
        *gmsg << "Restart from hdf5 file generated by OPAL-t" << endl;

        // force the initial time to zero
        beam.setT(0.0);
828 829
        beam.setLocalTrackStep((long long) 0 );
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
830 831 832
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](0) = -farray[n];

833
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
834 835 836
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](2) = farray[n];

837
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
838 839 840
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](1) = farray[n];

841
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
842 843 844
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](0) = -farray[n];

845
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "py", farray));
846 847 848
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](2) = farray[n];

849
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "pz", farray));
850 851 852
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](1) = farray[n];

853
        COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "id", larray));
854 855 856 857
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.ID[n] = larray[n];

    } else {
858

859 860
        *gmsg << "Restart from hdf5 file generated by OPAL-cycl" << endl;

861
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "AZIMUTH", &phi_m));
862

863
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ELEVATION", &psi_m));
864 865 866

        h5_int64_t localDump = 0;
        previousH5Local_m = false;
867

868 869 870 871 872 873 874 875 876 877 878
        rc = H5ReadStepAttribInt64(H5file, "LOCAL", &localDump);
        if(rc != H5_SUCCESS) {

            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);

            throw OpalException("Error during restart for Cyclotron Tracker:",
                                "You are trying to restart from a legacy file that doesn't contain\
  information on local/global frame. We are working on legacy support, but for now you have to use\
  OPAL 1.3.0!");

        } else {
879 880

            if (localDump == 1) previousH5Local_m = true;
881 882
        }

883
        double actualT;
884
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
885 886 887
        beam.setT(actualT);

        h5_int64_t SteptoLastInj;
888
        COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "SteptoLastInj", &SteptoLastInj));
889 890 891 892
        beam.setSteptoLastInj((int)SteptoLastInj);
        *gmsg << "Tracking Step since last bunch injection is " << SteptoLastInj << endl;

        h5_int64_t numBunch;
893
        COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "NumBunch", &numBunch));
894 895 896
        beam.setNumBunch((int)numBunch);
        *gmsg << numBunch << " Bunches(bins) exist in this file" << endl;

897
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
898 899 900
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](0) = farray[n];

901
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
902 903 904
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](1) = farray[n];

905
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
906 907 908
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](2) = farray[n];

909
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
910 911 912
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](0) = farray[n];

913
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "py", farray));
914 915 916
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](1) = farray[n];

917
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "pz", farray));
918 919 920
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](2) = farray[n];

921
        COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "id", larray));
922 923 924 925 926 927 928 929 930 931 932
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.ID[n] = larray[n];

        // only for multi-bunch mode
        if(specifiedNumBunch > 1) {
            // the allowed maximal bin number is set to 1000
            beam.setPBins(new PartBinsCyc(1000, numBunch));
        }
    }

    Ippl::Comm->barrier();
933
    COMPLAINONFAILURE(H5CloseFile(H5file));
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
    beam.boundp();
    beam.Q = beam.getChargePerParticle();

    if(strcmp(opalFlavour, "opal-t") == 0) {
        Vector_t meanR(0.0, 0.0, 0.0);
        Vector_t meanP(0.0, 0.0, 0.0);
        unsigned long int newLocalN = beam.getLocalNum();
        for(unsigned int i = 0; i < newLocalN; ++i) {
            for(int d = 0; d < 3; ++d) {
                meanR(d) += beam.R[i](d);
                meanP(d) += beam.P[i](d);
            }
        }
        reduce(meanR, meanR, OpAddAssign());
        meanR /= Vector_t(globalN);
        reduce(meanP, meanP, OpAddAssign());
        meanP /= Vector_t(globalN);
        *gmsg << "Rmean = " << meanR << "[m], Pmean=" << meanP << endl;

        for(unsigned int i = 0; i < newLocalN; ++i) {
            beam.R[i] -= meanR;
            beam.P[i] -= meanP;
        }
    }

    *gmsg << "----------------Finish reading hdf5 file----------------" << endl;
    IpplTimings::stopTimer(beam.distrReload_m);
}

void Distribution::DoRestartOpalE(EnvelopeBunch &beam, size_t Np, int restartStep) {
    h5_file_t *H5file;
    std::string fn;

    IpplTimings::startTimer(beam.distrReload_m);

    if(OpalData::getInstance()->hasRestartFile()) {
        fn = OpalData::getInstance()->getRestartFileName();
        *gmsg << "Restart from a specified file:" << fn << endl;

    } else {
        fn = OpalData::getInstance()->getInputBasename() + std::string(".h5");
    }

#ifdef PARALLEL_IO
    H5file = H5OpenFile(fn.c_str(), H5_O_RDONLY, Ippl::getComm());
#else
    H5file = H5PartOpenFile(fn.c_str(), H5_O_RDONLY, 0);
#endif

983
    if(H5file == (void*)H5_ERR) {
984
        throw OpalException("Distribution::DoRestartOpalE",
985
                            "could not open file '" + fn + "'");
986 987 988 989 990 991 992
    }

    if(restartStep == -1) {
        restartStep = H5GetNumSteps(H5file) - 1;
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
993 994
            throw OpalException("Distribution::DoRestartOpalE",
                                "can't append to the file '" + fn + "'");
995 996 997
        }
    }

998
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
999 1000 1001 1002 1003 1004 1005 1006 1007 1008
    int N = (int)H5PartGetNumParticles(H5file);

    h5_int64_t totalSteps = H5GetNumSteps(H5file);
    *gmsg << "total number of slices = " << N << " total steps " << totalSteps << endl;

    beam.distributeSlices(N);
    beam.createBunch();
    long long starti = beam.mySliceStartOffset();
    long long endi = beam.mySliceEndOffset();

1009
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
1010
    N = (int)H5PartGetNumParticles(H5file);
kraus's avatar
kraus committed
1011
    assert(N >= 0 && (unsigned int) N == beam.numMySlices());
1012 1013

    double actualT;
1014
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
1015 1016 1017

    beam.setT(actualT);
    double dPhiGlobal;
1018
    COMPLAINONFAILURE(H5ReadFileAttribFloat64(H5file, "dPhiGlobal", &dPhiGlobal));