diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp index fcfd363f3e0b74fac0c6d12ee1ae85455ed29a16..983ed6cbcd6041b837f7ef88f467fcf9528d6362 100644 --- a/src/Algorithms/ParallelCyclotronTracker.cpp +++ b/src/Algorithms/ParallelCyclotronTracker.cpp @@ -396,8 +396,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { // useful information spiral_flag = cycl_m->getSpiralFlag(); - if(spiral_flag) { - + if (spiral_flag) { *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl; *gmsg << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl; *gmsg << "* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl; @@ -416,7 +415,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { } // Fresh run (no restart): - if(!OpalData::getInstance()->inRestartRun()) { + if (!OpalData::getInstance()->inRestartRun()) { // Get reference values from cyclotron element // For now, these are still stored in mm. should be the only ones. -DW @@ -426,7 +425,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { referencePr = cycl_m->getPRinit(); referencePz = cycl_m->getPZinit(); - if(referenceTheta <= -180.0 || referenceTheta > 180.0) { + if (referenceTheta <= -180.0 || referenceTheta > 180.0) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "PHIINIT is out of [-180, 180)!"); } @@ -437,25 +436,22 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { float insqrt = referencePtot * referencePtot - \ referencePr * referencePr - referencePz * referencePz; - if(insqrt < 0) { - - if(insqrt > -1.0e-10) { + if (insqrt < 0) { + if (insqrt > -1.0e-10) { referencePt = 0.0; - } else { - throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "Pt imaginary!"); } } else { - referencePt = std::sqrt(insqrt); } - if(referencePtot < 0.0) + if (referencePtot < 0.0) { referencePt *= -1.0; + } // End calculate referencePt // Restart a run: @@ -463,15 +459,15 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { // If the user wants to save the restarted run in local frame, // make sure the previous h5 file was local too - if (Options::psDumpFrame != Options::GLOBAL) { - if (!previousH5Local) { + if (Options::psDumpFrame != Options::GLOBAL) { + if (!previousH5Local) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "You are trying a local restart from a global h5 file!"); - } + } // Else, if the user wants to save the restarted run in global frame, // make sure the previous h5 file was global too - } else { - if (previousH5Local) { + } else { + if (previousH5Local) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "You are trying a global restart from a local h5 file!"); } @@ -481,7 +477,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { referencePhi *= Physics::deg2rad; referencePsi *= Physics::deg2rad; referencePtot = bega; - if(referenceTheta <= -180.0 || referenceTheta > 180.0) { + if (referenceTheta <= -180.0 || referenceTheta > 180.0) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "PHIINIT is out of [-180, 180)!"); } @@ -539,21 +535,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { *gmsg << std::boolalpha << "* Superpose electric field maps -> " << superpose << endl; } - /** - * To ease the initialise() function, set a integral parameter fieldflag internally. - * Its value is by the option "TYPE" of the element "CYCLOTRON" - * fieldflag = 1, readin PSI format measured field file (default) - * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL - * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE - * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons - * fieldflag = 5, readin FFA format file for MSU/FNAL FFA - * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron - * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations - */ - int fieldflag = cycl_m->getFieldFlag(type); - - // Read in cyclotron field maps (midplane + 3D fields if desired). - cycl_m->initialise(itsBunch_m, fieldflag, cycl_m->getBScale()); + // Read in cyclotron field maps + cycl_m->initialise(itsBunch_m, cycl_m->getBScale()); double BcParameter[8] = {}; @@ -3520,4 +3503,4 @@ void ParallelCyclotronTracker::initPathLength() { // we need to reset the path length of each bunch itsDataSink->setMultiBunchInitialPathLengh(mbHandler_m.get()); } -} \ No newline at end of file +} diff --git a/src/Classic/AbsBeamline/Cyclotron.cpp b/src/Classic/AbsBeamline/Cyclotron.cpp index 157bf14b024bc313366e54702b72c2b38a998a11..6a42c8102c9c4111ceac5dae47724636fcdf422a 100644 --- a/src/Classic/AbsBeamline/Cyclotron.cpp +++ b/src/Classic/AbsBeamline/Cyclotron.cpp @@ -41,22 +41,20 @@ #include <cstdio> #include <cmath> +#include <boost/filesystem.hpp> + #define CHECK_CYC_FSCANF_EOF(arg) if (arg == EOF)\ throw GeneralClassicException("Cyclotron::getFieldFromFile",\ "fscanf returned EOF at " #arg); -extern Inform *gmsg; - +extern Inform* gmsg; -// Class Cyclotron -// ------------------------------------------------------------------------ Cyclotron::Cyclotron(): Component() { } - -Cyclotron::Cyclotron(const Cyclotron &right): +Cyclotron::Cyclotron(const Cyclotron& right): Component(right), fmapfn_m(right.fmapfn_m), rffrequ_m(right.rffrequ_m), @@ -71,7 +69,7 @@ Cyclotron::Cyclotron(const Cyclotron &right): pzinit_m(right.pzinit_m), spiral_flag_m(right.spiral_flag_m), trimCoilThreshold_m(right.trimCoilThreshold_m), - type_m(right.type_m), + typeName_m(right.typeName_m), harm_m(right.harm_m), bscale_m(right.bscale_m), trimcoils_m(right.trimcoils_m), @@ -86,19 +84,17 @@ Cyclotron::Cyclotron(const Cyclotron &right): RFVCoeff_fn_m(right.RFVCoeff_fn_m) { } - -Cyclotron::Cyclotron(const std::string &name): +Cyclotron::Cyclotron(const std::string& name): Component(name) { } - Cyclotron::~Cyclotron() { } void Cyclotron::applyTrimCoil_m(const double r, const double z, const double tet_rad, - double *br, double *bz) { + double* br, double* bz) { for (auto trimcoil : trimcoils_m) { trimcoil->applyField(r,tet_rad,z,br,bz); } @@ -108,10 +104,9 @@ void Cyclotron::applyTrimCoil(const double r, const double z, const double tet_rad, double& br, double& bz) { //Changed from > to >= to include case where bz == 0 and trimCoilThreshold_m == 0 -DW - if (std::abs(bz) >= trimCoilThreshold_m) { + if (std::abs(bz) >= trimCoilThreshold_m) { applyTrimCoil_m(r, z, tet_rad, &br, &bz); - } - else { + } else { // make sure to have a smooth transition double tmp_bz = 0.0; applyTrimCoil_m(r, z, tet_rad, &br, &tmp_bz); @@ -184,7 +179,12 @@ void Cyclotron::setFieldMapFN(std::string f) { } std::string Cyclotron::getFieldMapFN() const { - return fmapfn_m; + if (boost::filesystem::exists(fmapfn_m)) { + return fmapfn_m; + } else { + throw GeneralClassicException("Cyclotron::getFieldMapFN", + "Failed to open file '" + fmapfn_m + "', please check if it exists"); + } } void Cyclotron::setRfFieldMapFN(std::vector<std::string> f) { @@ -199,45 +199,14 @@ void Cyclotron::setRFVCoeffFN(std::vector<std::string> f) { RFVCoeff_fn_m = f; } -int Cyclotron::getFieldFlag(const std::string& type) const { - /** - * To ease the initialise() function, set a integral parameter fieldflag internally. - * Its value is by the option "TYPE" of the element "CYCLOTRON" - * fieldflag = 1, read in PSI format measured field file (default) - * fieldflag = 2, read in carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL - * fieldflag = 3, read in ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE - * fieldflag = 4, read in AVFEQ format file for Riken cyclotrons - * fieldflag = 5, read in FFA format file for MSU/FNAL FFA - * fieldflag = 6, read in both median plane B field map and 3D E field map - * of RF cavity for compact cyclotron - * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations - */ - int fieldflag; - if (type == std::string("CARBONCYCL")) { - fieldflag = 2; - } else if (type == std::string("CYCIAE")) { - fieldflag = 3; - } else if (type == std::string("AVFEQ")) { - fieldflag = 4; - } else if (type == std::string("FFA")) { - fieldflag = 5; - } else if (type == std::string("BANDRF")) { - fieldflag = 6; - } else if (type == std::string("SYNCHROCYCLOTRON")) { - fieldflag = 7; - } else //(type == "RING") - fieldflag = 1; - return fieldflag; -} - void Cyclotron::setRfPhi(std::vector<double> f) { rfphi_m = f; } double Cyclotron::getRfPhi(unsigned int i) const { - if (i < rfphi_m.size()) + if (i < rfphi_m.size()) { return rfphi_m[i]; - else { + } else { throw GeneralClassicException("Cyclotron::getRfPhi", "RFPHI not defined for CYCLOTRON!"); } @@ -248,9 +217,9 @@ void Cyclotron::setRfFrequ(std::vector<double> f) { } double Cyclotron::getRfFrequ(unsigned int i) const { - if (i < rffrequ_m.size()) + if (i < rffrequ_m.size()) { return rffrequ_m[i]; - else { + } else { throw GeneralClassicException("Cyclotron::getRfFrequ", "RFFREQ not defined for CYCLOTRON!"); } @@ -261,9 +230,9 @@ void Cyclotron::setSuperpose(std::vector<bool> flag) { } bool Cyclotron::getSuperpose(unsigned int i) const { - if (i < superpose_m.size()) + if (i < superpose_m.size()) { return superpose_m[i]; - else { + } else { throw GeneralClassicException("Cyclotron::getSuperpose", "SUPERPOSE not defined for CYCLOTRON!"); } @@ -278,12 +247,12 @@ double Cyclotron::getSymmetry() const { } -void Cyclotron::setType(std::string t) { - type_m = t; +void Cyclotron::setCyclotronType(std::string t) { + typeName_m = t; } const std::string &Cyclotron::getCyclotronType() const { - return type_m; + return typeName_m; } ElementBase::ElementType Cyclotron::getType() const { @@ -307,10 +276,10 @@ void Cyclotron::setEScale(std::vector<double> s) { } double Cyclotron::getEScale(unsigned int i) const { - if (i < escale_m.size()) + if (i < escale_m.size()) { return escale_m[i]; - else { - throw GeneralClassicException("Cyclotron::EScale", + } else { + throw GeneralClassicException("Cyclotron::getEScale", "EScale not defined for CYCLOTRON!"); } } @@ -327,12 +296,10 @@ double Cyclotron::getRmin() const { return BP.rmin; } - double Cyclotron::getRmax() const { return BP.rmin + (Bfield.nrad - 1) * BP.delr; } - void Cyclotron::setMinR(double r) { // DW: This is to let the user keep using mm in the input file for now // while switching internally to m @@ -344,6 +311,7 @@ void Cyclotron::setMaxR(double r) { // while switching internally to m maxr_m = 0.001 * r; } + double Cyclotron::getMinR() const { return minr_m; } @@ -357,19 +325,22 @@ void Cyclotron::setMinZ(double z) { // while switching internally to m minz_m = 0.001 * z; } + double Cyclotron::getMinZ() const { return minz_m; } + void Cyclotron::setMaxZ(double z) { // DW: This is to let the user keep using mm in the input file for now // while switching internally to m maxz_m = 0.001 * z; } + double Cyclotron::getMaxZ() const { return maxz_m; } -void Cyclotron::setTrimCoils(const std::vector<TrimCoil*> &trimcoils) { +void Cyclotron::setTrimCoils(const std::vector<TrimCoil*>& trimcoils) { trimcoils_m = trimcoils; } @@ -379,8 +350,31 @@ double Cyclotron::getFMLowE() const { return fmLowE_m;} void Cyclotron::setFMHighE(double e) { fmHighE_m = e;} double Cyclotron::getFMHighE() const { return fmHighE_m;} +void Cyclotron::setBFieldType() { + if (typeName_m.empty()) { + throw GeneralClassicException("Cyclotron::setBFieldType", + "TYPE is not defined for CYCLOTRON!"); + } else if (typeName_m == std::string("RING")) { + fieldType_m = BFieldType::PSIBF; + } else if (typeName_m == std::string("CARBONCYCL")) { + fieldType_m = BFieldType::CARBONBF; + } else if (typeName_m == std::string("CYCIAE")) { + fieldType_m = BFieldType::ANSYSBF; + } else if (typeName_m == std::string("AVFEQ")) { + fieldType_m = BFieldType::AVFEQBF; + } else if (typeName_m == std::string("FFA")) { + fieldType_m = BFieldType::FFABF; + } else if (typeName_m == std::string("BANDRF")) { + fieldType_m = BFieldType::BANDRF; + } else if (typeName_m == std::string("SYNCHROCYCLOTRON")) { + fieldType_m = BFieldType::SYNCHRO; + } else { + throw GeneralClassicException("Cyclotron::setBFieldType", + "TYPE " + typeName_m + " field reading function of CYCLOTRON is not defined!"); + } +} -bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B) { +bool Cyclotron::apply(const size_t& id, const double& t, Vector_t& E, Vector_t& B) { bool flagNeedUpdate = false; @@ -388,7 +382,7 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t & const double zpos = RefPartBunch_m->R[id](2); Inform gmsgALL("OPAL", INFORM_ALL_NODES); - if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m){ + if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m) { flagNeedUpdate = true; gmsgALL << level4 << getName() << ": Particle " << id @@ -396,9 +390,9 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t & gmsgALL << level4 << getName() << ": Coords: "<< RefPartBunch_m->R[id] << " m" << endl; - } else{ + } else { flagNeedUpdate = apply(RefPartBunch_m->R[id], RefPartBunch_m->P[id], t, E, B); - if (flagNeedUpdate){ + if (flagNeedUpdate) { gmsgALL << level4 << getName() << ": Particle "<< id << " out of the field map boundary!" << endl; @@ -416,8 +410,8 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t & return flagNeedUpdate; } -bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, - const double &t, Vector_t &E, Vector_t &B) { +bool Cyclotron::apply(const Vector_t& R, const Vector_t& /*P*/, + const double& t, Vector_t& E, Vector_t& B) { const double rad = std::hypot(R[0],R[1]); const double tempv = std::atan(R[1] / R[0]); @@ -463,7 +457,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, return true; } - if (myBFieldType_m != SYNCHRO && myBFieldType_m != BANDRF) { + if (fieldType_m != BFieldType::SYNCHRO && fieldType_m != BFieldType::BANDRF) { return false; } @@ -500,7 +494,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, double frequency = (*rffi); // frequency in MHz double ebscale = (*escali); // E and B field scaling - if (myBFieldType_m == SYNCHRO) { + if (fieldType_m == BFieldType::SYNCHRO) { double powert = 1; for (const double fcoef : (*rffci)) { powert *= (t * 1e-9); @@ -518,7 +512,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, E += ebscale * std::cos(phase) * tmpE; B -= ebscale * std::sin(phase) * tmpB; - if (myBFieldType_m != SYNCHRO) + if (fieldType_m != BFieldType::SYNCHRO) continue; // Some phase output -DW @@ -532,8 +526,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, *gmsg << "RF Frequency = " << frequency << " MHz" << endl; waiting_for_gap = 2; - } - else if (tet >= 270.0 && waiting_for_gap == 2) { + } else if (tet >= 270.0 && waiting_for_gap == 2) { phase_print = std::fmod(phase_print, 360) - 360.0; *gmsg << endl << "Gap 2 phase = " << phase_print << " Deg" << endl; @@ -542,7 +535,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &/*P*/, *gmsg << "RF Frequency = " << frequency << " MHz" << endl; waiting_for_gap = 0; } - if (myBFieldType_m == SYNCHRO) { + if (fieldType_m == BFieldType::SYNCHRO) { ++rffci, ++rfvci; } } @@ -569,9 +562,9 @@ bool Cyclotron::bends() const { } // calculate derivatives with 5-point lagrange's formula. -double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr) +double Cyclotron::gutdf5d(double* f, double dx, const int kor, + const int krl, const int lpr) { -{ double C[5][5][3], FAC[3]; double result; int j; @@ -640,7 +633,6 @@ double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, co C[3][4][1] = -104; C[4][4][1] = 35.0; - /* COEFFICIENTS FOR THE 3RD DERIVATIVE: */ C[0][0][2] = -10.0; C[1][0][2] = 36.0; @@ -686,8 +678,8 @@ bool Cyclotron::interpolate(const double& rad, const double& tet_rad, double& brint, double& btint, - double& bzint) -{ + double& bzint) { + const double xir = (rad - BP.rmin) / BP.delr; // ir : the number of path whose radius is less than the 4 points of cell which surround the particle. @@ -720,7 +712,7 @@ bool Cyclotron::interpolate(const double& rad, // r1t1 : the index of the "min angle, min radius" point in the 2D field array. // considering the array start with index of zero, minus 1. - if (myBFieldType_m != FFABF) { + if (fieldType_m != BFieldType::FFABF) { /* For FFA this does not work */ @@ -741,7 +733,6 @@ bool Cyclotron::interpolate(const double& rad, } if ((it >= 0) && (ir >= 0) && (it < Bfield.ntetS) && (ir < Bfield.nrad)) { - // B_{z} double bzf = Bfield.bfld[r1t1] * wr2 * wt2 + Bfield.bfld[r2t1] * wr1 * wt2 + @@ -769,48 +760,48 @@ bool Cyclotron::interpolate(const double& rad, } -void Cyclotron::read(const int &fieldflag, const double &scaleFactor) { - // PSIBF, AVFEQBF, ANSYSBF, FFABF - // for your own format field, you should add your own getFieldFromFile() function by yourself. - - if (fieldflag == 1) { - //*gmsg<<"Read field data from PSI format field map file."<<endl; - myBFieldType_m = PSIBF; - getFieldFromFile(scaleFactor); - - } else if (fieldflag == 2) { - // *gmsg<<"Read data from 450MeV Carbon cyclotron field file"<<endl - myBFieldType_m = CARBONBF; - getFieldFromFile_Carbon(scaleFactor); - - } else if (fieldflag == 3) { - // *gmsg<<"Read data from 100MeV H- cyclotron CYCIAE-100 field file"<<endl; - myBFieldType_m = ANSYSBF; - getFieldFromFile_CYCIAE(scaleFactor); - - } else if (fieldflag == 4) { - *gmsg << "* Read AVFEQ data (Riken) use bfield scale factor bs = " << getBScale() << endl; - myBFieldType_m = AVFEQBF; - getFieldFromFile_AVFEQ(scaleFactor); - - } else if (fieldflag == 5) { - *gmsg << "* Read FFA data MSU/FNAL " << getBScale() << endl; - myBFieldType_m = FFABF; - getFieldFromFile_FFA(scaleFactor); - - } else if (fieldflag == 6) { - *gmsg << "* Read both median plane B field map and 3D E field map of RF cavity for compact cyclotron" << endl; - myBFieldType_m = BANDRF; - getFieldFromFile_BandRF(scaleFactor); - - } else if (fieldflag == 7) { - *gmsg << "* Read midplane B-field, 3D RF fieldmaps, and text files with RF frequency/Voltage coefficients for Synchrocyclotron. (Midplane scaling = " << getBScale() << ")" << endl; - myBFieldType_m = SYNCHRO; - getFieldFromFile_Synchrocyclotron(scaleFactor); +void Cyclotron::read(const double& scaleFactor) { - } else - ERRORMSG("* The field reading function of this TYPE of CYCLOTRON has not implemented yet!" << endl); + setBFieldType(); + switch (fieldType_m) { + case BFieldType::PSIBF: { + *gmsg << "* Read field data from PSI format field map file" << endl; + getFieldFromFile_Ring(scaleFactor); + break; + } + case BFieldType::CARBONBF: { + *gmsg << "* Read data from 450MeV Carbon cyclotron field file" << endl; + getFieldFromFile_Carbon(scaleFactor); + break; + } + case BFieldType::ANSYSBF: { + *gmsg << "* Read data from 100MeV H- cyclotron CYCIAE-100 field file" << endl; + getFieldFromFile_CYCIAE(scaleFactor); + break; + } + case BFieldType::AVFEQBF: { + *gmsg << "* Read AVFEQ data (Riken)" << endl; + getFieldFromFile_AVFEQ(scaleFactor); + break; + } + case BFieldType::FFABF: { + *gmsg << "* Read FFA data MSU/FNAL" << endl; + getFieldFromFile_FFA(scaleFactor); + break; + } + case BFieldType::BANDRF: { + *gmsg << "* Read both median plane B field map and 3D E field map of RF cavity for compact cyclotron" << endl; + getFieldFromFile_BandRF(scaleFactor); + break; + } + case BFieldType::SYNCHRO: { + *gmsg << "* Read midplane B-field, 3D RF fieldmaps, and text files with RF frequency/Voltage coefficients for Synchrocyclotron" << endl; + getFieldFromFile_Synchrocyclotron(scaleFactor); + break; + } + } + // calculate the radii of initial grid. initR(BP.rmin, BP.delr, Bfield.nrad); @@ -834,7 +825,6 @@ void Cyclotron::getdiffs() { Bfield.g3.resize(Bfield.ntot); for (int i = 0; i < Bfield.nrad; i++) { - for (int k = 0; k < Bfield.ntet; k++) { double dtheta = Physics::deg2rad * BP.dtet; @@ -848,15 +838,12 @@ void Cyclotron::getdiffs() { int index = idx(i, k); int indexkEdge = idx(i, kEdge); - Bfield.dbt[index] = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 0, dkFromEdge, 1); Bfield.dbtt[index] = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 1, dkFromEdge, 1); Bfield.dbttt[index] = gutdf5d(&Bfield.bfld[indexkEdge], dtheta, 2, dkFromEdge, 1); } } - - for (int k = 0; k < Bfield.ntet; k++) { // inner loop varies R for (int i = 0; i < Bfield.nrad; i++) { @@ -869,7 +856,6 @@ void Cyclotron::getdiffs() { int index = idx(i, k); int indexredg = idx(iredg, k); - Bfield.dbr[index] = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 0, irtak, Bfield.ntetS); Bfield.dbrr[index] = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 1, irtak, Bfield.ntetS); Bfield.dbrrr[index] = gutdf5d(&Bfield.bfld[indexredg], BP.delr, 2, irtak, Bfield.ntetS); @@ -915,7 +901,6 @@ void Cyclotron::getdiffs() { Bfield.f2[iend] = Bfield.f2[istart]; Bfield.f3[iend] = Bfield.f3[istart]; Bfield.g3[iend] = Bfield.g3[istart]; - } /* debug @@ -931,8 +916,32 @@ void Cyclotron::getdiffs() { */ } -// read field map from external file. -void Cyclotron::getFieldFromFile(const double &scaleFactor) { + +// Calculates Radii of initial grid. +// dimensions in [m]! +void Cyclotron::initR(double rmin, double dr, int nrad) { + BP.rarr.resize(nrad); + for (int i = 0; i < nrad; i++) { + BP.rarr[i] = rmin + i * dr; + } + BP.delr = dr; +} + +void Cyclotron::initialise(PartBunchBase<double, 3>* bunch, double& /*startField*/, double& /*endField*/) { + RefPartBunch_m = bunch; + online_m = true; +} + +void Cyclotron::initialise(PartBunchBase<double, 3>* bunch, const double& scaleFactor) { + RefPartBunch_m = bunch; + lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump)); + + this->read(scaleFactor); +} + + +// Read field map from external file. +void Cyclotron::getFieldFromFile_Ring(const double& scaleFactor) { FILE *f = NULL; int lpar; @@ -946,11 +955,8 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) { BP.Bfact = scaleFactor; - if ((f = std::fopen(fmapfn_m.c_str(), "r")) == NULL) { - throw GeneralClassicException("Cyclotron::getFieldFromField", - "failed to open file '" + fmapfn_m + "', please check if it exists"); - } - + f = std::fopen(fmapfn_m.c_str(), "r"); + CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%lf", &BP.rmin)); *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl; BP.rmin *= 0.001; // mm --> m @@ -969,8 +975,9 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) { if (BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet); *gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl; - for (int i = 0; i < 13; i++) + for (int i = 0; i < 13; i++) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%s", fout)); + } CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%d", &Bfield.nrad)); *gmsg << "* Index in radial direction: " << Bfield.nrad << endl; @@ -1007,7 +1014,6 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%s", fout)); } Bfield.ntot = idx(Bfield.nrad - 1, Bfield.ntet) + 1; - //jjyang *gmsg << "* Total stored grid point number ( ntetS * nrad ) : " << Bfield.ntot << endl; Bfield.bfld.resize(Bfield.ntot); @@ -1018,7 +1024,6 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) { *gmsg << "* Read-in loop one block per radius" << endl; *gmsg << "* Rescaling of the magnetic fields with factor: " << BP.Bfact << endl; for (int i = 0; i < Bfield.nrad; i++) { - if (i > 0) { for (int dummy = 0; dummy < 6; dummy++) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%s", fout)); // INFO-LINE @@ -1043,36 +1048,11 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) { } std::fclose(f); - *gmsg << "* Field Map read successfully!" << endl << endl; } - -// Calculates Radii of initial grid. -// dimensions in [m]! -void Cyclotron::initR(double rmin, double dr, int nrad) { - BP.rarr.resize(nrad); - for (int i = 0; i < nrad; i++) { - BP.rarr[i] = rmin + i * dr; - } - BP.delr = dr; -} - -void Cyclotron::initialise(PartBunchBase<double, 3> *bunch, double &/*startField*/, double &/*endField*/) { - RefPartBunch_m = bunch; - online_m = true; -} - -void Cyclotron::initialise(PartBunchBase<double, 3> *bunch, const int &fieldflag, const double &scaleFactor) { - RefPartBunch_m = bunch; - lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump)); - - this->read(fieldflag, scaleFactor); -} - - -void Cyclotron::getFieldFromFile_FFA(const double &/*scaleFactor*/) { +void Cyclotron::getFieldFromFile_FFA(const double& /*scaleFactor*/) { /* Field is read in from ascii file (COSY output) in the order: @@ -1101,7 +1081,7 @@ void Cyclotron::getFieldFromFile_FFA(const double &/*scaleFactor*/) { std::vector<double>::iterator vit; *gmsg << "* ----------------------------------------------" << endl; - *gmsg << "* READ IN FFA FIELD MAP " << endl; + *gmsg << "* READ IN FFA FIELD MAP " << endl; *gmsg << "* ----------------------------------------------" << endl; BP.Bfact = -10.0; // T->kG and H- for the current FNAL FFA @@ -1187,11 +1167,12 @@ void Cyclotron::getFieldFromFile_FFA(const double &/*scaleFactor*/) { *gmsg << "* Field Map read successfully nelem= " << count << endl << endl; } -void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) { + +void Cyclotron::getFieldFromFile_AVFEQ(const double& scaleFactor) { FILE *f = NULL; *gmsg << "* ----------------------------------------------" << endl; - *gmsg << "* READ IN AVFEQ CYCLOTRON FIELD MAP " << endl; + *gmsg << "* READ IN AVFEQ CYCLOTRON FIELD MAP " << endl; *gmsg << "* ----------------------------------------------" << endl; /* From Hiroki-san @@ -1210,11 +1191,7 @@ void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) { BP.Bfact = scaleFactor / 1000.; - if ((f = std::fopen(fmapfn_m.c_str(), "r")) == NULL) { - throw GeneralClassicException( - "Cyclotron::getFieldFromFile_AVFEQ", - "failed to open file '" + fmapfn_m + "', please check if it exists"); - } + f = std::fopen(fmapfn_m.c_str(), "r"); CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%lf", &BP.rmin)); *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl; @@ -1291,8 +1268,7 @@ void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) { } -// read field map from external file. -void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) { +void Cyclotron::getFieldFromFile_Carbon(const double& scaleFactor) { FILE *f = NULL; *gmsg << "* ----------------------------------------------" << endl; @@ -1301,10 +1277,7 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) { BP.Bfact = scaleFactor; - if ((f = std::fopen(fmapfn_m.c_str(), "r")) == NULL) { - throw GeneralClassicException("Cyclotron::getFieldFromFile_Carbon", - "failed to open file '" + fmapfn_m + "', please check if it exists"); - } + f = std::fopen(fmapfn_m.c_str(), "r"); CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%lf", &BP.rmin)); *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl; @@ -1330,7 +1303,7 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%d", &Bfield.nrad)); *gmsg << "* Grid points along radius (nrad): " << Bfield.nrad << endl; -// Bfield.ntetS = Bfield.ntet; + //Bfield.ntetS = Bfield.ntet; Bfield.ntetS = Bfield.ntet + 1; //*gmsg << "* Accordingly, total grid point along azimuth: " << Bfield.ntetS << endl; @@ -1394,8 +1367,7 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) { } -// read field map from external file. -void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) { +void Cyclotron::getFieldFromFile_CYCIAE(const double& scaleFactor) { FILE *f = NULL; char fout[100]; @@ -1407,10 +1379,7 @@ void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) { BP.Bfact = scaleFactor; - if ((f = std::fopen(fmapfn_m.c_str(), "r")) == NULL) { - throw GeneralClassicException("Cyclotron::getFieldFromFile_CYCIAE", - "failed to open file '" + fmapfn_m + "', please check if it exists"); - } + f = std::fopen(fmapfn_m.c_str(), "r"); CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%lf", &BP.rmin)); *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl; @@ -1450,10 +1419,9 @@ void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) { int nHalfPoints = Bfield.ntet / 2.0 + 1; for (int i = 0; i < Bfield.nrad; i++) { - - for (int ii = 0; ii < 13; ii++) + for (int ii = 0; ii < 13; ii++) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%s", fout)); - + } for (int k = 0; k < nHalfPoints; k++) { CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%d", &dtmp)); CHECK_CYC_FSCANF_EOF(std::fscanf(f, "%d", &dtmp)); @@ -1462,7 +1430,6 @@ void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) { // T --> kGs, minus for minus hydrogen Bfield.bfld[idx(i, k)] = Bfield.bfld[idx(i, k)] * (-10.0); } - for (int k = nHalfPoints; k < Bfield.ntet; k++) { Bfield.bfld[idx(i, k)] = Bfield.bfld[idx(i, Bfield.ntet-k)]; } @@ -1473,17 +1440,13 @@ void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) { *gmsg << "* Field Map read successfully!" << endl << endl; } -void Cyclotron::getFieldFromFile_BandRF(const double &scaleFactor) { + +void Cyclotron::getFieldFromFile_BandRF(const double& scaleFactor) { // read 3D E&B field data file // loop over all field maps and superpose fields for (auto& fm: RFfilename_m) { Fieldmap *f = Fieldmap::getFieldmap(fm, false); - if (f == NULL) { - throw GeneralClassicException( - "Cyclotron::getFieldFromFile_BandRF", - "failed to open file '" + fm + "', please check if it exists"); - } *gmsg << "* Reading " << fm << endl; f->readMap(); RFfields_m.push_back(f); @@ -1492,7 +1455,8 @@ void Cyclotron::getFieldFromFile_BandRF(const double &scaleFactor) { getFieldFromFile_Carbon(scaleFactor); } -void Cyclotron::getFieldFromFile_Synchrocyclotron(const double &scaleFactor) { + +void Cyclotron::getFieldFromFile_Synchrocyclotron(const double& scaleFactor) { // read 3D E&B field data file std::vector<std::string>::const_iterator fm = RFfilename_m.begin(); @@ -1510,11 +1474,6 @@ void Cyclotron::getFieldFromFile_Synchrocyclotron(const double &scaleFactor) { for (; fm != RFfilename_m.end(); ++fm, ++rffcfni, ++rfvcfni, ++fcount) { Fieldmap *f = Fieldmap::getFieldmap(*fm, false); - if (f == NULL) { - throw GeneralClassicException( - "Cyclotron::getFieldFromFile_Synchrocyclotron", - "failed to open file '" + *fm + "', please check if it exists"); - } f->readMap(); // if (IPPL::Comm->getOutputLevel() != 0) f->getInfo(gmsg); RFfields_m.push_back(f); @@ -1575,7 +1534,7 @@ void Cyclotron::getFieldFromFile_Synchrocyclotron(const double &scaleFactor) { getFieldFromFile_Carbon(scaleFactor); } -void Cyclotron::getDimensions(double &/*zBegin*/, double &/*zEnd*/) const +void Cyclotron::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const { } #undef CHECK_CYC_FSCANF_EOF diff --git a/src/Classic/AbsBeamline/Cyclotron.h b/src/Classic/AbsBeamline/Cyclotron.h index 3561bfb7afc942d24881cf269c516434ae0700cf..1d574efc062febc4fa25e8528bd984cff9fe982a 100644 --- a/src/Classic/AbsBeamline/Cyclotron.h +++ b/src/Classic/AbsBeamline/Cyclotron.h @@ -35,7 +35,6 @@ class Fieldmap; class LossDataSink; class TrimCoil; -enum BFieldType {PSIBF,CARBONBF,ANSYSBF,AVFEQBF,FFABF,BANDRF,SYNCHRO}; struct BfieldData { std::string filename; @@ -60,19 +59,12 @@ struct BfieldData { std::vector<double> g3; // for Btheta // Grid-Size - //need to be read from inputfile. - int nrad, ntet; - - // one more grid line is stored in azimuthal direction: - int ntetS; - - // total grid points number. - int ntot; + int nrad, ntet; //need to be read from inputfile. + int ntetS; // one more grid line is stored in azimuthal direction + int ntot; // total grid points number. // Mean and Maximas double bacc, dbtmx, dbttmx, dbtttmx; - - }; struct BPositions { @@ -88,16 +80,22 @@ struct BPositions { }; -// Class Cyclotron -// ------------------------------------------------------------------------ -/// Interface for a Cyclotron. -// This class defines the abstract interface for a Cyclotron. - class Cyclotron: public Component { + public: + enum class BFieldType { + PSIBF, + CARBONBF, + ANSYSBF, + AVFEQBF, + FFABF, + BANDRF, + SYNCHRO + }; + /// Constructor with given name. - explicit Cyclotron(const std::string &name); + explicit Cyclotron(const std::string& name); Cyclotron(); Cyclotron(const Cyclotron &); @@ -122,13 +120,11 @@ public: void setRFFCoeffFN(std::vector<std::string> rff_coeff_fn); void setRFVCoeffFN(std::vector<std::string> rfv_coeff_fn); - int getFieldFlag(const std::string& type) const; - - void setType(std::string t); + void setCyclotronType(std::string t); const std::string &getCyclotronType() const; virtual ElementBase::ElementType getType() const; - virtual void getDimensions(double &zBegin, double &zEnd) const; + virtual void getDimensions(double& zBegin, double& zEnd) const; unsigned int getNumberOfTrimcoils() const; @@ -165,7 +161,7 @@ public: void setEScale(std::vector<double> bs); virtual double getEScale(unsigned int i) const; - void setTrimCoils(const std::vector<TrimCoil*> &trimcoils); + void setTrimCoils(const std::vector<TrimCoil*>& trimcoils); void setSuperpose(std::vector<bool> flag); virtual bool getSuperpose(unsigned int i) const; @@ -191,17 +187,17 @@ public: void setSpiralFlag(bool spiral_flag); virtual bool getSpiralFlag() const; - virtual bool apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B); + virtual bool apply(const size_t& id, const double& t, Vector_t& E, Vector_t& B); - virtual bool apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B); + virtual bool apply(const Vector_t& R, const Vector_t& P, const double& t, Vector_t& E, Vector_t& B); virtual void apply(const double& rad, const double& z, const double& tet_rad, double& br, double& bt, double& bz); - virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField); + virtual void initialise(PartBunchBase<double, 3>* bunch, double& startField, double& endField); - virtual void initialise(PartBunchBase<double, 3> *bunch, const int &fieldflag, const double &scaleFactor); + virtual void initialise(PartBunchBase<double, 3>* bunch, const double& scaleFactor); virtual void finalise(); @@ -217,34 +213,38 @@ public: double& bt, double& bz); - void read(const int &fieldflag, const double &scaleFactor); - + void read(const double& scaleFactor); + + void setBFieldType(); + + BFieldType fieldType_m; + private: /// Apply trim coils (calculate field contributions) with smooth field transition void applyTrimCoil (const double r, const double z, const double tet_rad, double& br, double& bz); /// Apply trim coils (calculate field contributions) - void applyTrimCoil_m(const double r, const double z, const double tet_rad, double *br, double *bz); + void applyTrimCoil_m(const double r, const double z, const double tet_rad, double* br, double* bz); + protected: - - + void getdiffs(); - double gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr); + double gutdf5d(double* f, double dx, const int kor, const int krl, const int lpr); void initR(double rmin, double dr, int nrad); - void getFieldFromFile(const double &scaleFactor); - void getFieldFromFile_Carbon(const double &scaleFactor); - void getFieldFromFile_CYCIAE(const double &scaleFactor); - void getFieldFromFile_AVFEQ(const double &scaleFactor); - void getFieldFromFile_FFA(const double &scaleFactor); - void getFieldFromFile_BandRF(const double &scaleFactor); - void getFieldFromFile_Synchrocyclotron(const double &scaleFactor); + void getFieldFromFile_Ring(const double& scaleFactor); + void getFieldFromFile_Carbon(const double& scaleFactor); + void getFieldFromFile_CYCIAE(const double& scaleFactor); + void getFieldFromFile_AVFEQ(const double& scaleFactor); + void getFieldFromFile_FFA(const double& scaleFactor); + void getFieldFromFile_BandRF(const double& scaleFactor); + void getFieldFromFile_Synchrocyclotron(const double& scaleFactor); inline int idx(int irad, int ktet) {return (ktet + Bfield.ntetS * irad);} - + private: std::string fmapfn_m; /* stores the filename of the fieldmap */ @@ -267,7 +267,7 @@ private: bool spiral_flag_m; double trimCoilThreshold_m; ///< B-field threshold for applying trim coil - std::string type_m; /* what type of field we use */ + std::string typeName_m; // name of the TYPE parameter in cyclotron double harm_m; double bscale_m; // a scale factor for the B-field @@ -287,9 +287,6 @@ private: // Not implemented. void operator=(const Cyclotron &) = delete; - - BFieldType myBFieldType_m; - // RF field map handler // Fieldmap *RFfield; std::vector<Fieldmap *> RFfields_m; diff --git a/src/Distribution/ClosedOrbitFinder.h b/src/Distribution/ClosedOrbitFinder.h index 27771a60b7ceaa8a02918a163571a410cd0a8f9b..f4b27192add67d98e6837bac269980aeab3307f9 100644 --- a/src/Distribution/ClosedOrbitFinder.h +++ b/src/Distribution/ClosedOrbitFinder.h @@ -290,8 +290,7 @@ ClosedOrbitFinder<Value_type, N_m /= cycl_m->getSymmetry(); } - cycl_m->read(cycl_m->getFieldFlag(cycl_m->getCyclotronType()), - cycl_m->getBScale()); + cycl_m->read(cycl_m->getBScale()); // reserve storage for the orbit and momentum (--> size = 0, capacity = N_m+1) /* @@ -925,4 +924,4 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, cons } -#endif \ No newline at end of file +#endif diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp index b1fc947c6c98eb6a820559bb818c89c5e95cb486..d09917bd5f0d9f9959a49ab25becad7d6d3dc2b4 100644 --- a/src/Distribution/Distribution.cpp +++ b/src/Distribution/Distribution.cpp @@ -1225,9 +1225,10 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, else *gmsg << "* SECTOR: " << "match using single sector" << endl; - *gmsg << "* NSTEPS = " << Nint << endl - << "* HN= " << CyclotronElement->getCyclHarm() - << " PHIINIT= " << CyclotronElement->getPHIinit() << endl + *gmsg << "* NSTEPS = " << Nint << endl + << "* HN = " << CyclotronElement->getCyclHarm() + << " PHIINIT = " << CyclotronElement->getPHIinit() << endl + << "* FIELD MAP = " << CyclotronElement->getFieldMapFN() << endl << "* ----------------------------------------------------" << endl; if ( CyclotronElement->getFMLowE() < 0 || diff --git a/src/Elements/OpalCyclotron.cpp b/src/Elements/OpalCyclotron.cpp index 7f7f09458f528bb976a73381049326cfce2c8135..3fe0d1f473fcea10c3d39406a5ad8e799a76db72 100644 --- a/src/Elements/OpalCyclotron.cpp +++ b/src/Elements/OpalCyclotron.cpp @@ -170,7 +170,7 @@ void OpalCyclotron::update() { cycl->setBScale(bscale); - cycl->setType(type); + cycl->setCyclotronType(type); cycl->setCyclHarm(harmnum); cycl->setMinR(minr); @@ -218,13 +218,13 @@ void OpalCyclotron::update() { cycl->setRfFrequ(rff_str); cycl->setSuperpose(superpose_str); - if(itsAttr[GEOMETRY] && obgeo_m == nullptr) { + if (itsAttr[GEOMETRY] && obgeo_m == nullptr) { obgeo_m = BoundaryGeometry::find(Attributes::getString(itsAttr[GEOMETRY])); - if(obgeo_m) { + if (obgeo_m) { cycl->setBoundaryGeometry(obgeo_m); } } // Transmit "unknown" attributes. OpalElement::updateUnknown(cycl); -} \ No newline at end of file +}