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 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
    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),
    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)
267
{
268 269 270 271
    SetAttributes();

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

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

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

Distribution::~Distribution() {

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

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

377 378 379 380 381 382 383 384 385 386 387 388 389 390
    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
391 392

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

/**
405 406
 * At the moment only write the header into the file dist.dat
 * PartBunch will then append (very uggly)
407 408 409
 * @param
 * @param
 * @param
410
 */
411
void Distribution::WriteToFile() {
412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427
    /*
      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();
      }
      }
    */
428 429
}

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

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

439 440
void Distribution::execute() {
}
adelmann's avatar
adelmann committed
441

442 443
void Distribution::update() {
}
adelmann's avatar
adelmann committed
444

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

447
    SetFieldEmissionParameters();
adelmann's avatar
adelmann committed
448

449
    switch (distrTypeT_m) {
adelmann's avatar
adelmann committed
450

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

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

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

484 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
    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
513

514 515 516 517 518 519 520
                        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
521
                        } else {
522 523
                            beam->R[lowMark + count] = bg.getCooridinate(count);
                            beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
gsell's avatar
gsell committed
524
                        }
525 526 527 528 529 530 531 532 533
                        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
534 535
                    }
                }
536 537 538
                pc++;
                if(pc == Ippl::getNodes())
                    pc = 0;
gsell's avatar
gsell committed
539
            }
540 541 542
            bg.clearCooridinateArray();
            bg.clearMomentaArray();
            beam->boundp();
gsell's avatar
gsell committed
543
        }
544
        *gmsg << *beam << endl;
gsell's avatar
gsell committed
545

546
    } else {// Normal procedure to create primary particles
gsell's avatar
gsell committed
547

548 549 550 551 552
        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
553

554 555 556 557
        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
558

559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575
                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
576 577

                    }
578

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

584
            }
585

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

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

597
    IpplTimings::startTimer(beam.distrReload_m);
598

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

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

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

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

623 624
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));

gsell's avatar
gsell committed
625 626
    int N = (int)H5PartGetNumParticles(H5file);

627 628 629 630 631 632 633 634 635
    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
636

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

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

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

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

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

657 658 659 660 661 662 663 664
    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;

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

669 670
    beam.create(N);

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

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

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

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

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

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

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

    Ippl::Comm->barrier();
705
    COMPLAINONFAILURE(H5CloseFile(H5file));
706 707 708

    beam.boundp();

gsell's avatar
gsell committed
709 710
    IpplTimings::stopTimer(beam.distrReload_m);

711 712 713 714
    *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
715 716
}

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

719 720 721 722 723 724 725 726 727 728 729 730 731
    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

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

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

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

750
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
751 752 753 754 755 756 757 758 759 760 761
    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
762
        endi = starti + numberOfParticlesPerNode - 1;
kraus's avatar
kraus committed
763

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

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

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

    char opalFlavour[128];
777
    COMPLAINONFAILURE(H5ReadStepAttribString(H5file, "OPAL_flavour", opalFlavour));
778

779
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPR",&referencePr_m));
780

781
    referencePt_m = 0.0;
782
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPT",&referencePt_m));
783

784
    referencePz_m = 0.0;
785
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFPZ",&referencePz_m));
786

787
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFR",&referenceR_m));
788

789
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFTHETA",&referenceTheta_m));
790

791
    referenceZ_m = 0.0;
792
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "REFZ",&referenceZ_m));
793

794
    double meanE;
795
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ENERGY", &meanE));
796

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

    beam.setLPath(pathLength);

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

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

808
    bega_m = be * ga;
809

810
    *gmsg << "* Gamma = " << ga << ", Beta = " << be << endl;
811

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

814
    double *farray = reinterpret_cast<double *>(varray.get());
815

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

    beam.create(localN);

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

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

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

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

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

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

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

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

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

    } else {
856

857 858
        *gmsg << "Restart from hdf5 file generated by OPAL-cycl" << endl;

859
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "AZIMUTH", &phi_m));
860

861
        COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "ELEVATION", &psi_m));
862 863 864

        h5_int64_t localDump = 0;
        previousH5Local_m = false;
865

866 867 868 869 870 871 872 873 874 875 876
        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 {
877 878

            if (localDump == 1) previousH5Local_m = true;
879 880
        }

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

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

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

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

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

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

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

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

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

919
        COMPLAINONFAILURE(H5PartReadDataInt64(H5file, "id", larray));
920 921 922 923 924 925 926 927 928 929 930
        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();
931
    COMPLAINONFAILURE(H5CloseFile(H5file));
932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980
    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

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

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

996
    COMPLAINONFAILURE(H5SetStep(H5file, restartStep));
997 998 999 1000 1001 1002 1003 1004 1005 1006
    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();

1007
    COMPLAINONFAILURE(H5PartSetView(H5file, starti, endi));
1008 1009 1010 1011
    N = (int)H5PartGetNumParticles(H5file);
    assert(N >= 0 && (unsigned int) N != beam.numMySlices());

    double actualT;
1012
    COMPLAINONFAILURE(H5ReadStepAttribFloat64(H5file, "TIME", &actualT));
1013 1014 1015

    beam.setT(actualT);
    double dPhiGlobal;
1016
    COMPLAINONFAILURE(H5ReadFileAttribFloat64(H5file, "dPhiGlobal", &dPhiGlobal));
1017 1018 1019
    OpalData::getInstance()->setGlobalPhaseShift(dPhiGlobal);

    h5_int64_t ltstep;
1020
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "LocalTrackStep", &ltstep));
1021 1022 1023
    beam.setLocalTrackStep((long long)ltstep);

    h5_int64_t gtstep;
1024
    COMPLAINONFAILURE(H5ReadStepAttribInt64(H5file, "GlobalTrackStep",