Distribution.cpp 188 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 57 58
#define COMPLAINONFAILURE(rc)     if(rc != H5_SUCCESS)                  \
        ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);

gsell's avatar
gsell committed
59 60 61 62
//
// Class Distribution
// ------------------------------------------------------------------------

63 64 65
namespace AttributesT
{
    enum AttributesT {
66 67 68 69 70 71 72 73 74 75 76 77 78 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
        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
169 170 171
    };
}

172 173 174
namespace LegacyAttributesT
{
    enum LegacyAttributesT {
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
        // DESCRIPTION OF THE DISTRIBUTION:
        DEBIN,
        SBIN,
        TEMISSION,
        SIGLASER,
        AG,
        SIGMAPT,
        TRANSVCUTOFF,
        CUTOFF,
        Z,
        T,
        PT,
        ALPHAX,
        ALPHAY,
        BETAX,
        BETAY,
        DX,
        DDX,
        DY,
        DDY,
        SIZE
    };
197 198
}

gsell's avatar
gsell committed
199
Distribution::Distribution():
200 201 202
    Definition(AttributesT::SIZE + LegacyAttributesT::SIZE, "DISTRIBUTION",
               "The DISTRIBUTION statement defines data for the 6D particle distribution."),
    distrTypeT_m(DistrTypeT::NODIST),
203
    numberOfDistributions_m(1),
204 205 206 207 208 209 210 211 212 213 214 215 216
    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),
217 218 219 220 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
                                                                                    pTotThermal_m(0.0),
                                                                                    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),
                                                                                    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)
253
{
254 255 256 257
    SetAttributes();

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

    try {
260
        OpalData::getInstance()->define(defaultDistribution);
gsell's avatar
gsell committed
261
    } catch(...) {
262 263 264 265
        delete defaultDistribution;
    }

    SetFieldEmissionParameters();
gsell's avatar
gsell committed
266 267 268 269 270 271 272
}
/**
 *
 *
 * @param name
 * @param parent
 */
273
Distribution::Distribution(const std::string &name, Distribution *parent):
gsell's avatar
gsell committed
274
    Definition(name, parent),
275 276
    distT_m(parent->distT_m),
    distrTypeT_m(DistrTypeT::NODIST),
277
    numberOfDistributions_m(parent->numberOfDistributions_m),
278 279 280 281 282 283
    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
284
    tEmission_m(parent->tEmission_m),
285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311
    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),
    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),
312
    avrgpz_m(parent->avrgpz_m),
313 314 315 316 317 318 319 320 321 322 323 324 325
    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),
    distCorr_m(parent->distCorr_m),
    laserProfileFileName_m(parent->laserProfileFileName_m),
    laserImageName_m(parent->laserImageName_m),
    laserIntensityCut_m(parent->laserIntensityCut_m),
    laserProfile_m(NULL),
gsell's avatar
gsell committed
326 327 328 329 330 331 332 333 334 335 336 337
    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),
338 339 340 341 342 343 344
    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),
345 346
    cutoff_m(parent->cutoff_m),
    I_m(parent->I_m),
adelmann's avatar
adelmann committed
347 348 349
    E_m(parent->E_m),
    bega_m(parent->bega_m),
    M_m(parent->M_m)
350
{
gsell's avatar
gsell committed
351 352 353 354 355 356 357
}

Distribution::~Distribution() {

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

358 359 360
    if (energyBins_m != NULL) {
        delete energyBins_m;
        energyBins_m = NULL;
gsell's avatar
gsell committed
361 362
    }

363 364 365 366 367 368 369 370 371 372 373 374 375 376
    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
377 378

}
379
/*
380
  void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
381
  for(int i=0; i<M.size1(); ++i) {
382 383 384 385 386
  for(int j=0; j<M.size2(); ++j) {
  *gmsg  << M(i,j) << " ";
  }
  *gmsg << endl;
  }
387 388
  }
*/
gsell's avatar
gsell committed
389 390

/**
391 392
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
393 394 395
 * @param
 * @param
 * @param
396
 */
397
void Distribution::WriteToFile() {
398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
    /*
      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();
      }
      }
    */
414 415
}

416 417 418
/// Distribution can only be replaced by another distribution.
bool Distribution::canReplaceBy(Object *object) {
    return dynamic_cast<Distribution *>(object) != 0;
419 420
}

421 422 423
Distribution *Distribution::clone(const std::string &name) {
    return new Distribution(name, this);
}
424

425 426
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
427

428 429
void Distribution::update() {
}
adelmann's avatar
adelmann committed
430

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

433
    SetFieldEmissionParameters();
adelmann's avatar
adelmann committed
434

435
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
436

437
    case DistrTypeT::MATCHEDGAUSS:
438
        CreateMatchedGaussDistribution(numberOfParticles, massIneV);
439
        break;
440 441
    case DistrTypeT::FROMFILE:
        CreateDistributionFromFile(numberOfParticles, massIneV);
gsell's avatar
gsell committed
442
        break;
443 444
    case DistrTypeT::GAUSS:
        CreateDistributionGauss(numberOfParticles, massIneV);
gsell's avatar
gsell committed
445
        break;
446 447
    case DistrTypeT::BINOMIAL:
        CreateDistributionBinomial(numberOfParticles, massIneV);
gsell's avatar
gsell committed
448
        break;
449 450
    case DistrTypeT::FLATTOP:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
451
        break;
452 453
    case DistrTypeT::GUNGAUSSFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
454
        break;
455 456
    case DistrTypeT::ASTRAFLATTOPTH:
        CreateDistributionFlattop(numberOfParticles, massIneV);
gsell's avatar
gsell committed
457
        break;
458 459 460 461
    default:
        INFOMSG("Distribution unknown." << endl;);
        break;
    }
gsell's avatar
gsell committed
462

463
    // AAA Scale and shift coordinates according to distribution input.
464 465 466
    ScaleDistCoordinates();
    ShiftDistCoordinates(massIneV);
}
gsell's avatar
gsell committed
467

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

470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498
    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
499

500 501 502 503 504 505 506
                        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
507
                        } else {
508 509
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
510
                        }
511 512 513 514 515 516 517 518 519
                        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
520 521
                    }
                }
522 523 524
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
525
            }
526 527 528
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
529
        }
530
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
531

532
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
533

534 535 536 537 538
        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
539

540 541 542 543
        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
544

545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
                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
562 563

                    }
564

gsell's avatar
gsell committed
565
                }
566 567 568
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
569

570
            }
571

572 573 574 575 576 577
        }
        bg.clearCooridinateArray();
        beam->boundp();//fixme if bg.getN()==0?
    }
    *gmsg << *beam << endl;
}
578

579 580 581
void Distribution::DoRestartOpalT(PartBunch &beam, size_t Np, int restartStep) {
    h5_file_t *H5file;
    std::string fn;
582

583
    IpplTimings::startTimer(beam.distrReload_m);
584

585 586
    //        beam.setTEmission(Attributes::getReal(itsAttr[AttributesT::SIZE + LegacyAttributesT::TEMISSION]));
    fn = OpalData::getInstance()->getInputBasename() + std::string(".h5");
gsell's avatar
gsell committed
587 588 589 590

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

594
    if(H5file == (void*)H5_ERR) {
595 596
	throw OpalException("Distribution::DoRestartOpalT",
                            "could not open file '" + fn + "'");
gsell's avatar
gsell committed
597 598 599
    }

    if(restartStep == -1) {
600
        restartStep = H5GetNumSteps(H5file) - 1 ;
gsell's avatar
gsell committed
601 602 603
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
604 605
            throw OpalException("Distribution::DoRestartOpalT",
                                "can't append to the file '" + fn + "'");
gsell's avatar
gsell committed
606 607 608
        }
    }

609 610
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));

gsell's avatar
gsell committed
611 612
    int N = (int)H5PartGetNumParticles(H5file);

613 614 615 616 617 618 619 620 621
    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
622

623
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
gsell's avatar
gsell committed
624
    N = (int)H5PartGetNumParticles(H5file);
625
    assert(N >= 0);
gsell's avatar
gsell committed
626 627

    double actualT;
628
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
gsell's avatar
gsell committed
629
    beam.setT(actualT);
630

gsell's avatar
gsell committed
631
    double dPhiGlobal;
632
    COMPLAINONFAILURE(H5ReadFileAttribFloat64(H5file, "dPhiGlobal", &dPhiGlobal));
gsell's avatar
gsell committed
633 634
    OpalData::getInstance()->setGlobalPhaseShift(dPhiGlobal);

635
    h5_int64_t ltstep;
636
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
637 638 639
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
640
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep", &gtstep));
641
    beam.setGlobalTrackStep((long long)gtstep);
642

643 644 645 646 647 648 649 650
    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;

651
    std::unique_ptr<char[]> varray(new char[(N)*sizeof(double)]);
652
    h5_float64_t *farray = reinterpret_cast<h5_float64_t *>(varray.get());
653
    h5_int64_t *larray = reinterpret_cast<h5_int64_t *>(varray.get());
gsell's avatar
gsell committed
654

655 656
    beam.create(N);

657
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
gsell's avatar
gsell committed
658
    for(unsigned long int n = 0; n < (unsigned int) N; ++n) {
659 660
        beam.R[n](0) = farray[n];
        beam.Bin[n] = 0; // not initialized
gsell's avatar
gsell committed
661
    }
662
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
gsell's avatar
gsell committed
663
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
664
        beam.R[n](1) = farray[n];
gsell's avatar
gsell committed
665

666
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
gsell's avatar
gsell committed
667
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
668
        beam.R[n](2) = farray[n];
gsell's avatar
gsell committed
669

670
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
gsell's avatar
gsell committed
671
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
672
        beam.P[n](0) = farray[n];
gsell's avatar
gsell committed
673

674
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "py", farray));
gsell's avatar
gsell committed
675
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
676
        beam.P[n](1) = farray[n];
gsell's avatar
gsell committed
677

678
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "pz", farray));
gsell's avatar
gsell committed
679
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
680
        beam.P[n](2) = farray[n];
gsell's avatar
gsell committed
681

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

686
    COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "lastsection", larray));
gsell's avatar
gsell committed
687 688 689 690
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
        beam.LastSection[n] = (short) larray[n];

    Ippl::Comm->barrier();
691
    COMPLAINONFAILURE(H5CloseFile(H5file));
692 693 694

    beam.boundp();

gsell's avatar
gsell committed
695 696
    IpplTimings::stopTimer(beam.distrReload_m);

697 698 699 700
    *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
701 702
}

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

705 706 707 708 709 710 711 712 713 714 715 716 717
    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

718
    if(H5file == (void*)H5_ERR) {
719 720
        throw OpalException("Distribution::DoRestartOpalCycl",
                            "could not open file '" + fn + "'");
721 722 723 724 725 726 727
    }

    if(restartStep == -1) {
        restartStep = H5GetNumSteps(H5file) - 1;
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
728 729
            throw OpalException("Distribution::DoRestartOpalCycl",
                                "can't append to the file '" + fn + "'");
730 731 732
        }
    }

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

736
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
737 738 739 740 741 742 743 744 745 746 747
    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
748
        endi = starti + numberOfParticlesPerNode - 1;
kraus's avatar
kraus committed
749

750
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
751 752 753 754
    const int localN = (int)H5PartGetNumParticles(H5file);
    assert(localN >= 0);

    h5_int64_t ltstep;
755
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
756 757 758
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
759
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep", &gtstep));
760 761 762
    beam.setGlobalTrackStep((long long)gtstep);

    char opalFlavour[128];
763
    COMPLAINONFAILURE(H5ReadStepAttribString(H5file, "OPAL_flavour", opalFlavour));
764

765
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPR",&referencePr_m));
766

767
    referencePt_m = 0.0;
768
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPT",&referencePt_m));
769

770
    referencePz_m = 0.0;
771
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPZ",&referencePz_m));
772

773
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFR",&referenceR_m));
774

775
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFTHETA",&referenceTheta_m));
776

777
    referenceZ_m = 0.0;
778
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFZ",&referenceZ_m));
779

780
    double meanE;
781
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ENERGY", &meanE));
782

adelmann's avatar
adelmann committed
783
    double pathLength;;
784
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "SPOS", &pathLength));
adelmann's avatar
adelmann committed
785 786 787

    beam.setLPath(pathLength);

788 789 790
    *gmsg << "* Restart Energy = " << meanE << " (MeV), Path lenght = " << pathLength << " (m)" <<  endl;

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

794
    bega_m = be * ga;
795

796
    *gmsg << "* Gamma = " << ga << ", Beta = " << be << endl;
797

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

800
    double *farray = reinterpret_cast<double *>(varray.get());
801

802 803 804 805 806
    h5_int64_t *larray = reinterpret_cast<h5_int64_t *>(varray.get());

    beam.create(localN);

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

808 809 810 811
        *gmsg << "Restart from hdf5 file generated by OPAL-t" << endl;

        // force the initial time to zero
        beam.setT(0.0);
812 813
        beam.setLocalTrackStep((long long) 0 );
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
814 815 816
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](0) = -farray[n];

817
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
818 819 820
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](2) = farray[n];

821
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
822 823 824
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](1) = farray[n];

825
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
826 827 828
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](0) = -farray[n];

829
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "py", farray));
830 831 832
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](2) = farray[n];

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

837
        COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "id", larray));
838 839 840 841
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.ID[n] = larray[n];

    } else {
842

843 844
        *gmsg << "Restart from hdf5 file generated by OPAL-cycl" << endl;

845
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "AZIMUTH", &phi_m));
846

847
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ELEVATION", &psi_m));
848 849 850

        h5_int64_t localDump = 0;
        previousH5Local_m = false;
851

852 853 854 855 856 857 858 859 860 861 862
        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 {
863 864

            if (localDump == 1) previousH5Local_m = true;
865 866
        }

867
        double actualT;
868
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
869 870 871
        beam.setT(actualT);

        h5_int64_t SteptoLastInj;
872
        COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "SteptoLastInj", &SteptoLastInj));
873 874 875 876
        beam.setSteptoLastInj((int)SteptoLastInj);
        *gmsg << "Tracking Step since last bunch injection is " << SteptoLastInj << endl;

        h5_int64_t numBunch;
877
        COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "NumBunch", &numBunch));
878 879 880
        beam.setNumBunch((int)numBunch);
        *gmsg << numBunch << " Bunches(bins) exist in this file" << endl;

881
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
882 883 884
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](0) = farray[n];

885
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
886 887 888
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](1) = farray[n];

889
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));
890 891 892
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.R[n](2) = farray[n];

893
        COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "px", farray));
894 895 896
        for(unsigned long int n = 0; n < (unsigned int) localN; ++n)
            beam.P[n](0) = farray[n];

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

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

905
        COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "id", larray));
906 907 908 909 910 911 912 913 914 915 916
        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();
917
    COMPLAINONFAILURE(H5CloseFile(H5file));
918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966
    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

967
    if(H5file == (void*)H5_ERR) {
968
        throw OpalException("Distribution::DoRestartOpalE",
969
                            "could not open file '" + fn + "'");
970 971 972 973 974 975 976
    }

    if(restartStep == -1) {
        restartStep = H5GetNumSteps(H5file) - 1;
        OpalData::getInstance()->setRestartStep(restartStep);
    } else {
        if(restartStep != H5GetNumSteps(H5file) - 1 && !OpalData::getInstance()->hasRestartFile()) {
977 978
            throw OpalException("Distribution::DoRestartOpalE",
                                "can't append to the file '" + fn + "'");
979 980 981
        }
    }

982
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
983 984 985 986 987 988 989 990 991 992
    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();

993
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
994 995 996 997
    N = (int)H5PartGetNumParticles(H5file);
    assert(N >= 0 && (unsigned int) N != beam.numMySlices());

    double actualT;
998
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
999 1000 1001

    beam.setT(actualT);
    double dPhiGlobal;
1002
    COMPLAINONFAILURE(H5ReadFileAttribFloat64(H5file, "dPhiGlobal", &dPhiGlobal));
1003 1004 1005
    OpalData::getInstance()->setGlobalPhaseShift(dPhiGlobal);

    h5_int64_t ltstep;
1006
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
1007 1008 1009
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
1010
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep", &gtstep));
1011 1012 1013 1014 1015 1016
    beam.setGlobalTrackStep((long long)gtstep);

    std::unique_ptr<char[]> varray(new char[(N)*sizeof(double)]);
    double *farray = reinterpret_cast<double *>(varray.get());
    h5_int64_t *larray = reinterpret_cast<h5_int64_t *>(varray.get());

1017
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "x", farray));
1018 1019 1020
    for(unsigned long int n = 0; n < (unsigned int) N; ++n) {
        beam.setX(n, farray[n]);
    }
1021
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "y", farray));
1022 1023 1024
    for(unsigned long int n = 0; n < (unsigned int) N; ++n)
        beam.setY(n, farray[n]);

1025
    COMPLAINONFAILURE(H5PartReadDataFloat64(H5file, "z", farray));