Commit 90600494 authored by Daniel Winklehner's avatar Daniel Winklehner
Browse files

Added the option to define a reference Z and PZ to OPAL cyclotron -DW

parent 1cc69499
......@@ -55,6 +55,8 @@ Cyclotron::Cyclotron(const Cyclotron &right):
rinit_m(right.rinit_m),
prinit_m(right.prinit_m),
phiinit_m(right.phiinit_m),
zinit_m(right.zinit_m),
pzinit_m(right.pzinit_m),
type_m(right.type_m),
harm_m(right.harm_m),
bscale_m(right.bscale_m),
......@@ -92,7 +94,6 @@ double Cyclotron::getRinit() const {
return rinit_m;
}
void Cyclotron::setPRinit(double prinit) {
prinit_m = prinit;
}
......@@ -109,6 +110,22 @@ double Cyclotron::getPHIinit() const {
return phiinit_m;
}
void Cyclotron::setZinit(double zinit){
zinit_m = zinit;
}
double Cyclotron::getZinit() const {
return zinit_m;
}
void Cyclotron::setPZinit(double pzinit){
pzinit_m = pzinit;
}
double Cyclotron::getPZinit() const {
return pzinit_m;
}
void Cyclotron::setFieldMapFN(string f) {
fmapfn_m = f;
}
......@@ -352,7 +369,7 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &centroid, const double
} else {
/*
With t his we have B-field AND this is far more
With this we have B-field AND this is far more
intuitive for me ....
*/
r1t1 = idx(ir, it);
......@@ -748,8 +765,8 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) {
BP.Bfact = scaleFactor;
if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
ERRORMSG("Error in Cyclotron::getFieldFromFile()!" << endl);
ERRORMSG(" Cannot open file, please check if it really exists." << endl);
ERRORMSG("* Error in Cyclotron::getFieldFromFile()!" << endl);
ERRORMSG("* Cannot open file, please check if it really exists." << endl);
exit(1);
}
......@@ -813,8 +830,8 @@ void Cyclotron::getFieldFromFile(const double &scaleFactor) {
Bfield.dbtt.resize(Bfield.ntot);
Bfield.dbttt.resize(Bfield.ntot);
*gmsg << "* read -in loop one block per radius" << endl;
*gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;
*gmsg << "* Read-in loop one block per radius" << endl;
*gmsg << "* Rescaling of the fields with factor: " << BP.Bfact << endl;
for(int i = 0; i < Bfield.nrad; i++) {
if(i > 0) {
......@@ -885,12 +902,12 @@ void Cyclotron::initialise(PartBunch *bunch, const int &fieldflag, const double
getFieldFromFile_CYCIAE(scaleFactor);
} else if(fieldflag == 4) {
*gmsg << "* Read AVFEQ data (Riken) use bfield scale factor bs= " << getBScale() << endl;
*gmsg << "* Read AVFEQ data (Riken) use bfield scale factor bs = " << getBScale() << endl;
myBFieldType_m = AVFEQBF;
getFieldFromFile_AVFEQ(scaleFactor);
} else if(fieldflag == 5) {
*gmsg << "Read FFAG data MSU/FNAL " << getBScale() << endl;
*gmsg << "* Read FFAG data MSU/FNAL " << getBScale() << endl;
myBFieldType_m = FFAGBF;
getFieldFromFile_FFAG(scaleFactor);
......@@ -900,7 +917,7 @@ void Cyclotron::initialise(PartBunch *bunch, const int &fieldflag, const double
getFieldFromFile_BandRF(scaleFactor);
} else
ERRORMSG("The field reading function of this TYPE of CYCLOTRON has not implemented yet.!" << endl);
ERRORMSG("* The field reading function of this TYPE of CYCLOTRON has not implemented yet!" << endl);
// calculate the radii of initial grid.
initR(BP.rmin, BP.delr, Bfield.nrad);
......@@ -1001,7 +1018,7 @@ void Cyclotron::getFieldFromFile_FFAG(const double &scaleFactor) {
Bfield.dbtt.resize(Bfield.ntot);
Bfield.dbttt.resize(Bfield.ntot);
*gmsg << "* rescaling of the fields with factor: " << BP.Bfact << endl;
*gmsg << "* Rescaling of the fields with factor: " << BP.Bfact << endl;
int count = 0;
if((Ippl::getNodes()) == 1 && Options::info) {
......@@ -1043,8 +1060,8 @@ void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) {
BP.Bfact = scaleFactor / 1000.;
if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
ERRORMSG("Error in Cyclotron::getFieldFromFile_AVFEQ()!" << endl);
ERRORMSG(" Cannot open file, please check if it really exists." << endl);
ERRORMSG("* Error in Cyclotron::getFieldFromFile_AVFEQ()!" << endl);
ERRORMSG("* Cannot open file, please check if it really exists." << endl);
exit(1);
}
......@@ -1125,7 +1142,7 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) {
if((f = fopen(fmapfn_m.c_str(), "r")) == NULL) {
ERRORMSG("* Error in Cyclotron::getFieldFromFile_Carbon()!" << endl);
ERRORMSG(" Cannot open file, please check if it really exists." << endl);
ERRORMSG("* Cannot open file, please check if it really exists." << endl);
exit(1);
}
......@@ -1136,12 +1153,12 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) {
*gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
*gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;
*gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.dtet));
//if the value is nagtive, the actual value is its reciprocal.
if(BP.dtet < 0.0) BP.dtet = 1.0 / (-BP.dtet);
*gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg.]" << endl;
*gmsg << "* Stepsize in azimuth direction: " << BP.dtet << " [deg]" << endl;
CHECK_CYC_FSCANF_EOF(fscanf(f, "%d", &Bfield.ntet));
*gmsg << "* Index in azimuthal direction: " << Bfield.ntet << endl;
......
......@@ -155,6 +155,12 @@ public:
void setPHIinit(double phiinit);
virtual double getPHIinit() const;
void setZinit(double zinit);
virtual double getZinit() const;
void setPZinit(double zinit);
virtual double getPZinit() const;
void setBScale(double bs);
virtual double getBScale() const;
......@@ -217,6 +223,8 @@ private:
double rinit_m;
double prinit_m;
double phiinit_m;
double zinit_m;
double pzinit_m;
std::string type_m; /* what type of field we use */
double harm_m;
......
......@@ -325,8 +325,10 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
if(!OpalData::getInstance()->inRestartRun()) {
// get values from cyclotron command
referenceR = elptr->getRinit();
referencePr = elptr->getPRinit();
referenceTheta = elptr->getPHIinit();
referenceZ = elptr->getZinit();
referencePr = elptr->getPRinit();
referencePz = elptr->getPZinit();
//msg << "PRINIT= " << pri << " [CU]" << endl;
if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
......@@ -344,8 +346,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
// TEMP for testing if we can inject a bunch vertically for inflector -DW
// This will probably not work for long bunches, also I don't know what happens for a restart
Vector_t pmean = itsBunch->get_pmean();
referencePz = pmean[2];
//Vector_t pmean = itsBunch->get_pmean();
//referencePz = pmean[2];
//referencePr = sqrt(pmean[0] * pmean[0] + pmean[1] * pmean[1]);
float insqrt = referencePtot * referencePtot - referencePr * referencePr - referencePz * referencePz;
if(insqrt < 0) {
......@@ -368,23 +370,33 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
if(referencePtot < 0.0) referencePt *= -1.0;
sinRefTheta_m = sin(referenceTheta / 180.0 * pi);
cosRefTheta_m = cos(referenceTheta / 180.0 * pi);
cosRefTheta_m = cos(referenceTheta / 180.0 * pi);
*gmsg << "" << endl;
*gmsg << "*** Bunch global starting position: ***" << endl;
*gmsg << "* RINIT = " << referenceR << " [mm]" << endl;
*gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
*gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
*gmsg << "" << endl;
*gmsg << "*** Bunch global starting momenta: ***" << endl;
*gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
*gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
*gmsg << "* Total reference momentum = " << referencePtot * 1000.0 << " [MCU]" << endl;
*gmsg << "* Total reference momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference azimuthal momentum = " << referencePt * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference radial momentum = " << referencePr * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference axial momentum = " << referencePz * 1000.0 << " [MCU]" << endl;
*gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
double sym = elptr->getSymmetry();
*gmsg << "* " << sym << "-fold field symmerty " << endl;
......@@ -2412,6 +2424,7 @@ void ParallelCyclotronTracker::Tracker_RK4() {
double E = itsBunch->get_meanEnergy();
itsBunch->R /= Vector_t(1000.0); // mm --> m
itsDataSink->writeStatData(*itsBunch, FDext_m , 0.0, 0.0, 0.0, E);
lastDumpedStep_m = itsDataSink->writePhaseSpace_cycl(*itsBunch, FDext_m, E, referencePr, referenceR, referenceTheta);
itsBunch->R *= Vector_t(1000.0); // m --> mm
*gmsg << "* Phase space dump " << lastDumpedStep_m << " (global frame) at integration step "
......@@ -3497,7 +3510,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
// Initialize global R
itsBunch->R *= Vector_t(1000.0); // m --> mm
Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m, referenceR * sinRefTheta_m, 0.0); // [referenceR] == mm
Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m, referenceR * sinRefTheta_m, referenceZ); // [referenceR] == mm
localToGlobal(itsBunch->R, initialReferenceTheta, initMeanR);
......@@ -3505,6 +3518,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
for(size_t i = 0; i < initialLocalNum_m; ++i) {
itsBunch->P[i](0) += referencePr;
itsBunch->P[i](1) += referencePt;
itsBunch->P[i](2) += referencePz;
}
// Initialize the bin number of the first bunch to 0
......@@ -3517,17 +3531,18 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
double E = itsBunch->get_meanEnergy();
itsBunch->R *= Vector_t(0.001); // mm --> m
itsBunch->calcBeamParameters_cycl();
lastDumpedStep_m = itsDataSink->writePhaseSpace_cycl(*itsBunch, FDext_m, E, referencePr, referenceR, referenceTheta);
itsBunch->R *= Vector_t(1000.0); // m --> mm
*gmsg << "* Phase space dump " << lastDumpedStep_m << " (global frame) at integration step 0 T= 0 [ns]" << endl;
*gmsg << "* Phase space dump " << lastDumpedStep_m << " (global frame) at integration step 0 T = 0 [ns]" << endl;
// Initial dump (if requested in local frame)
} else {
// Initial dump (if requested in local frame)
itsBunch->R *= Vector_t(0.001); // mm --> m
itsBunch->calcBeamParameters_cycl();
double E = itsBunch->get_meanEnergy();
lastDumpedStep_m = itsDataSink->writePhaseSpace_cycl(*itsBunch, FDext_m, E, referencePr, referenceR, referenceTheta);
itsDataSink->writeStatData(*itsBunch, FDext_m, 0, 0, 0, E);
*gmsg << "* Phase space dump " << lastDumpedStep_m << " (local frame) at integration step 0 T= 0 [ns]" << endl;
*gmsg << "* Phase space dump " << lastDumpedStep_m << " (local frame) at integration step 0 T = 0 [ns]" << endl;
itsBunch->R *= Vector_t(1000.0); // m --> mm
}
......@@ -3559,13 +3574,14 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
// Initialize global R
itsBunch->R *= Vector_t(1000.0); // m --> mm
Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m, referenceR * sinRefTheta_m, 0.0); // [referenceR] == mm
Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m, referenceR * sinRefTheta_m, referenceZ); // [referenceR] == mm
localToGlobal(itsBunch->R, initialReferenceTheta, initMeanR);
// Initialize global P
for(size_t i = 0; i < initialLocalNum_m; ++i) {
itsBunch->P[i](0) += referencePr;
itsBunch->P[i](1) += referencePt;
itsBunch->P[i](2) += referencePz;
}
localToGlobal(itsBunch->P, initialReferenceTheta);
......
......@@ -32,27 +32,29 @@ public:
/// The attributes of class OpalCyclotron.
enum {
TYPE,
GEOMETRY, // geometry of boundary
CYHARMON, // The harmonic number of the cyclotron
SYMMETRY, // The symetry of the field
RINIT, // The initial radius [m]
PRINIT, // The initial radial momenta [pr/p0] []
PHIINIT, // The initial phase [deg]
RFFREQ, // First hamonic of the RF system
FMAPFN, // The filename of the fieldmap
RFMAPFN, // The filename of the fieldmap
BSCALE, // A scalar to scale the B-field
ESCALE, // A scalar to scale the RF field
TCR1, //trim coil r1 (mm)
TCR2, //trim coil r2 (mm)
MBTC, //max bfield of trim coil (kG)
SLPTC, //slope of the rising edge
RFPHI, // the initial phase of RF field
GEOMETRY, // geometry of boundary
CYHARMON, // The harmonic number of the cyclotron
SYMMETRY, // The symetry of the field
RINIT, // The initial radius [m]
PRINIT, // The initial radial momentum [pr/p0] []
PHIINIT, // The initial phase [deg]
ZINIT, // The initial z coordinate [m]
PZINIT, // The initial vertical momentum [pz/p0] []
RFFREQ, // First hamonic of the RF system
FMAPFN, // The filename of the fieldmap
RFMAPFN, // The filename of the fieldmap
BSCALE, // A scalar to scale the B-field
ESCALE, // A scalar to scale the RF field
TCR1, //trim coil r1 (mm)
TCR2, //trim coil r2 (mm)
MBTC, //max bfield of trim coil (kG)
SLPTC, //slope of the rising edge
RFPHI, // the initial phase of RF field
SUPERPOSE, // whether the electric field map are superposed or not
MINZ, // minimal vertical extend of the machine
MAXZ, // maximal vertical extend of the machine
MINR, // minimal radial extend of the machine
MAXR, // maximal radial extend of the machine
MINZ, // minimal vertical extend of the machine
MAXZ, // maximal vertical extend of the machine
MINR, // minimal radial extend of the machine
MAXR, // maximal radial extend of the machine
SIZE
};
......
......@@ -927,7 +927,8 @@ void DataSink::writePhaseSpace(PartBunch &beam, Vector_t FDext[], double sposHea
int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double meanEnergy, double refPr, double refR, double refTheta) {
if (!doHDF5_m) return -1;
//if (beam.getLocalNum() == 0) return -1; //TEMP for testing -DW
h5_int64_t rc;
IpplTimings::startTimer(H5PartTimer_m);
......@@ -1123,7 +1124,6 @@ int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double mea
/*
Attributes originally not in OPAL-cycl
*/
double sposHead = 0.0;
......@@ -1141,10 +1141,8 @@ int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double mea
/*
Done attributes originally not in OPAL-cycl
*/
setOPALcycl();
for(size_t i = 0; i < nLoc; i++)
......@@ -1208,7 +1206,7 @@ int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double mea
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
// Need this because rc = H5root does not yet work with
// a vaiable number of data
// avaiable number of data
if(Options::ebDump) {
for(size_t i = 0; i < nLoc; i++)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment