Cleanup and bug fixing after changing cyclotron tracker internal units to m -DW

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -222,9 +222,8 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
 
     Vector_t intecoords = 0.0;
 
-    // This has to match the dT in the rk4 pusher! -DW
-    //double dtime = 0.5 * itsBunch->getdT();  // Old
-    double dtime = itsBunch->getdT() * getHarmonicNumber();  // New
+    // This has to match the dT in the rk4 pusher
+    double dtime = itsBunch->getdT() * getHarmonicNumber();
 
     int triId = 0;
     for(size_t i = 0; i < itsBunch->getLocalNum(); i++) {
@@ -523,9 +522,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
     } else //(type == "RING")
         fieldflag = 1;
 
-    // read field map on the  middle plane of cyclotron.
-    // currently scalefactor is set to 1.0
-    // TEMP changed 1.0 to getBScale() to test if we can scale the midplane field -DW
+    // Read in cyclotron field maps (midplane + 3D fields if desired).
     elptr->initialise(itsBunch, fieldflag, elptr->getBScale());
 
     double BcParameter[8];
@@ -3376,6 +3373,9 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
 }
 
 void ParallelCyclotronTracker::applyPluginElements(const double dt) {
+    // Plugin Elements are all defined in mm, change beam to mm before applying
+
+    itsBunch->R *= Vector_t(1000.0);
 
     for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); sindex++) {
         if(((*sindex)->first) == ElementBase::SEPTUM)    {
@@ -3401,6 +3401,9 @@ void ParallelCyclotronTracker::applyPluginElements(const double dt) {
             collim->checkCollimator(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt);
         }
     }
+
+    itsBunch->R *= Vector_t(0.001);
+
 }
 
 bool ParallelCyclotronTracker::deleteParticle(){

src/Classic/AbsBeamline/Cyclotron.cpp

@@ -1141,6 +1141,7 @@ void Cyclotron::getFieldFromFile_FFAG(const double &scaleFactor) {
     for(int i = 0; i < num_of_header_lines; ++i)
         file_to_read.ignore(max_num_of_char_in_a_line, '\n');
 
+/*
     while(!file_to_read.eof()) {
         double r, th, x, y, bz;
         file_to_read >> r >> th >> x >> y >> bz;
@@ -1152,6 +1153,20 @@ void Cyclotron::getFieldFromFile_FFAG(const double &scaleFactor) {
             bzv.push_back(bz);
         }
     }
+*/
+
+    // TEMP for OPAL 2.0 changing this to m -DW
+    while(!file_to_read.eof()) {
+        double r, th, x, y, bz;
+        file_to_read >> r >> th >> x >> y >> bz;
+        if((int)th != 360) {
+            rv.push_back(r);
+            thv.push_back(th);
+            xv.push_back(x);
+            yv.push_back(y);
+            bzv.push_back(bz);
+        }
+    }
 
     double maxtheta = 360.0;
     BP.dtet = thv[1] - thv[0];
@@ -1167,9 +1182,9 @@ void Cyclotron::getFieldFromFile_FFAG(const double &scaleFactor) {
     Bfield.ntet = (int)((maxtheta - thv[0]) / BP.dtet);
     Bfield.nrad  = (int)(rmax - BP.rmin) / BP.delr + 1;
     Bfield.ntetS  = Bfield.ntet + 1;
-    *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
-    *gmsg << "* Maximal radius of measured field map: " << rmax << " [mm]" << endl;
-    *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
+    *gmsg << "* Minimal radius of measured field map: " << 1000.0 * BP.rmin << " [mm]" << endl;
+    *gmsg << "* Maximal radius of measured field map: " << 1000.0 * rmax << " [mm]" << endl;
+    *gmsg << "* Stepsize in radial direction: " << 1000.0 * BP.delr << " [mm]" << endl;
     *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;
     *gmsg << "* Maximal angle of measured field map: " << maxtheta << " [deg.]" << endl;
 
@@ -1235,13 +1250,16 @@ void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) {
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
     *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
+    BP.rmin *= 0.001;  // mm --> m
 
     double rmax;
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &rmax));
     *gmsg << "* Maximal radius of measured field map: " << rmax << " [mm]" << endl;
+    rmax *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
     *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
+    BP.delr *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
     *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;
@@ -1261,7 +1279,6 @@ void Cyclotron::getFieldFromFile_AVFEQ(const double &scaleFactor) {
 
     Bfield.nrad = (int)(rmax - BP.rmin) / BP.delr;
 
-
     int ntotidx = idx(Bfield.nrad, Bfield.ntetS) + 1;
 
     Bfield.ntot = Bfield.ntetS * Bfield.nrad;
@@ -1315,11 +1332,13 @@ void Cyclotron::getFieldFromFile_Carbon(const double &scaleFactor) {
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
     *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
+    BP.rmin *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
     //if the value is negative, the actual value is its reciprocal.
     if(BP.delr < 0.0) BP.delr = 1.0 / (-BP.delr);
     *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
+    BP.delr *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
     *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg]" << endl;
@@ -1413,9 +1432,11 @@ void Cyclotron::getFieldFromFile_CYCIAE(const double &scaleFactor) {
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.rmin));
     *gmsg << "* Minimal radius of measured field map: " << BP.rmin << " [mm]" << endl;
+    BP.rmin *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.delr));
     *gmsg << "* Stepsize in radial direction: " << BP.delr << " [mm]" << endl;
+    BP.delr *= 0.001;  // mm --> m
 
     CHECK_CYC_FSCANF_EOF(fscanf(f, "%lf", &BP.tetmin));
     *gmsg << "* Minimal angle of measured field map: " << BP.tetmin << " [deg.]" << endl;

src/Classic/AbsBeamline/Ring.cpp

@@ -101,8 +101,8 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
         Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
         gmsgALL << getName() << ": particle " << id
                 << " at " << refPartBunch_m->R[id]
-                << " out of the field map boundary" << endl;
-        lossDS_m->addParticle(refPartBunch_m->R[id], refPartBunch_m->P[id], id);
+                << " m out of the field map boundary" << endl;
+        lossDS_m->addParticle(refPartBunch_m->R[id] * Vector_t(1000.0), refPartBunch_m->P[id], id);
         lossDS_m->save();
 
         refPartBunch_m->Bin[id] = -1;
@@ -116,13 +116,14 @@ bool Ring::apply(const Vector_t &R, const Vector_t &P,
     B = Vector_t(0.0, 0.0, 0.0);
     E = Vector_t(0.0, 0.0, 0.0);
 
-    std::vector<RingSection*> sections = getSectionsAt(R);
+    std::vector<RingSection*> sections = getSectionsAt(R); // I think this doesn't actually use R -DW
     bool outOfBounds = true;
     // assume field maps don't set B, E to 0...
     for (size_t i = 0; i < sections.size(); ++i) {
         Vector_t B_temp(0.0, 0.0, 0.0);
         Vector_t E_temp(0.0, 0.0, 0.0);
-        outOfBounds &= sections[i]->getFieldValue(R, refPartBunch_m->get_centroid(), t, E_temp, B_temp);
+        // Super-TEMP! cyclotron tracker now uses m internally, have to change to mm here to match old field limits -DW
+        outOfBounds &= sections[i]->getFieldValue(R * Vector_t(1000.0), refPartBunch_m->get_centroid() * Vector_t(1000.0), t, E_temp, B_temp);
         B += (scale_m * B_temp);
         E += (scale_m * E_temp);
     }