diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp index e21c6ec18528bb451183234d243d53f506b9d142..dacc7dc651a9e799c24610c7a67afb439757fca3 100644 --- a/src/Algorithms/ParallelCyclotronTracker.cpp +++ b/src/Algorithms/ParallelCyclotronTracker.cpp @@ -85,16 +85,15 @@ class Beamline; class PartData; -using Physics::c; +//using Physics::c; using Physics::m_p; // GeV using Physics::PMASS; using Physics::PCHARGE; using Physics::pi; using Physics::q_e; -const double c_mmtns = c * 1.0e-6; // m/s --> mm/ns -const double mass_coeff = 1.0e18 * q_e / c / c; // from GeV/c^2 to basic unit: GV*C*s^2/m^2 -const double PIOVER180 = pi / 180.0; +const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns +const double mass_coeff = 1.0e18 * q_e / Physics::c / Physics::c; // from GeV/c^2 to basic unit: GV*C*s^2/m^2 Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0); Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0); @@ -115,8 +114,8 @@ extern Inform *gmsg; * @param revTrack */ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline, - const PartData &reference, - bool revBeam, bool revTrack): + const PartData &reference, + bool revBeam, bool revTrack): Tracker(beamline, reference, revBeam, revTrack), eta_m(0.01), myNode_m(Ippl::myNode()), @@ -190,18 +189,18 @@ ParallelCyclotronTracker::~ParallelCyclotronTracker() { * @param none */ void ParallelCyclotronTracker::initializeBoundaryGeometry() { - for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) { - bgf_m = dynamic_cast<ElementBase *>(*compindex)->getBoundaryGeometry(); - if(!bgf_m) - continue; - else - break; - } - if (bgf_m) { - itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk")); - OpalData::getInstance()->setGlobalGeometry(bgf_m); - *gmsg << "* Boundary geometry initialized " << endl; - } + for(std::list<Component *>::iterator compindex = myElements.begin(); compindex != myElements.end(); compindex++) { + bgf_m = dynamic_cast<ElementBase *>(*compindex)->getBoundaryGeometry(); + if(!bgf_m) + continue; + else + break; + } + if (bgf_m) { + itsDataSink->writeGeomToVtk(*bgf_m, std::string("data/testGeometry-00000.vtk")); + OpalData::getInstance()->setGlobalGeometry(bgf_m); + *gmsg << "* Boundary geometry initialized " << endl; + } } /** @@ -210,29 +209,29 @@ void ParallelCyclotronTracker::initializeBoundaryGeometry() { * @param fn Base file name */ void ParallelCyclotronTracker::bgf_main_collision_test() { - if(!bgf_m) return; + if(!bgf_m) return; - Inform msg("bgf_main_collision_test "); + Inform msg("bgf_main_collision_test "); - /** - *Here we check if a particles is outside the domain, flag it for deletion - */ + /** + *Here we check if a particles is outside the domain, flag it for deletion + */ - Vector_t intecoords = 0.0; + 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! -DW + //double dtime = 0.5 * itsBunch->getdT(); // Old + double dtime = itsBunch->getdT() * getHarmonicNumber(); // New - int triId = 0; - size_t Nimpact = 0; - for(size_t i = 0; i < itsBunch->getLocalNum(); i++) { - int res = bgf_m->PartInside(itsBunch->R[i]*1.0e-3, itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId); - if(res >= 0) { - itsBunch->Bin[i] = -1; - Nimpact++; + int triId = 0; + size_t Nimpact = 0; + for(size_t i = 0; i < itsBunch->getLocalNum(); i++) { + int res = bgf_m->PartInside(itsBunch->R[i]*1.0e-3, itsBunch->P[i], dtime, itsBunch->PType[i], itsBunch->Q[i], intecoords, triId); + if(res >= 0) { + itsBunch->Bin[i] = -1; + Nimpact++; + } } - } } @@ -321,15 +320,15 @@ void ParallelCyclotronTracker::visitRing(const Ring &ring) { if(referencePtot < 0.0) referencePt *= -1.0; - sinRefTheta_m = sin(referenceTheta * PIOVER180); - cosRefTheta_m = cos(referenceTheta * PIOVER180); + sinRefTheta_m = sin(referenceTheta * Physics::deg2rad); + cosRefTheta_m = cos(referenceTheta * Physics::deg2rad); double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - buildupFieldList(BcParameter, "OPALRING", opalRing_m); + buildupFieldList(BcParameter, ElementBase::RING, opalRing_m); // Finally print some diagnostic *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl; @@ -359,21 +358,21 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone()); myElements.push_back(elptr); - // Is this a Spiral Inflector Simulation? If yes, we'll give the user some + // Is this a Spiral Inflector Simulation? If yes, we'll give the user some // useful information spiral_flag = elptr->getSpiralFlag(); 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 << "* 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; *gmsg << "* and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl; - *gmsg << "* 2.) It is strongly recommended to use the SAAMG fieldsolver," << endl; + *gmsg << "* 2.) It is strongly recommended to use the SAAMG fieldsolver," << endl; *gmsg << "* FFT does not give the correct results (boundaty conditions are missing)." << endl; *gmsg << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl; - *gmsg << "* There will be no transformations of the bunch into a local frame und consequently," << endl; - *gmsg << "* the problem will be treated non-relativistically!" << endl; + *gmsg << "* There will be no transformations of the bunch into a local frame und consequently," << endl; + *gmsg << "* the problem will be treated non-relativistically!" << endl; *gmsg << "* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl; *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl; *gmsg << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl; @@ -400,7 +399,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum: float insqrt = referencePtot * referencePtot - \ - referencePr * referencePr - referencePz * referencePz; + referencePr * referencePr - referencePz * referencePz; if(insqrt < 0) { @@ -423,7 +422,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { referencePt *= -1.0; // End calculate referencePt - // Restart a run: + // Restart a run: } else { // If the user wants to save the restarted run in local frame, @@ -435,8 +434,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { 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 the user wants to save the restarted run in global frame, + // make sure the previous h5 file was global too } else { if (previousH5Local) { @@ -447,8 +446,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { } // Adjust some of the reference variables from the h5 file - referencePhi *= PIOVER180; - referencePsi *= PIOVER180; + referencePhi *= Physics::deg2rad; + referencePsi *= Physics::deg2rad; referencePtot = bega; if(referenceTheta <= -180.0 || referenceTheta > 180.0) { @@ -458,8 +457,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { } } - sinRefTheta_m = sin(referenceTheta * PIOVER180); - cosRefTheta_m = cos(referenceTheta * PIOVER180); + sinRefTheta_m = sin(referenceTheta * Physics::deg2rad); + cosRefTheta_m = cos(referenceTheta * Physics::deg2rad); *gmsg << endl; *gmsg << "* Bunch global starting position:" << endl; @@ -498,22 +497,22 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [mm]"<< endl; /** - bool Sflag = elptr->getSuperpose(); - std::string flagsuperposed; - if (Sflag) - flagsuperposed="yes"; - else - flagsuperposed="no"; - *gmsg << "* Electric field maps are superposed? " << flagsuperposed << " " << endl; - */ + bool Sflag = elptr->getSuperpose(); + std::string flagsuperposed; + if (Sflag) + flagsuperposed="yes"; + else + flagsuperposed="no"; + *gmsg << "* Electric field maps are superposed? " << flagsuperposed << " " << endl; + */ double h = elptr->getCyclHarm(); *gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl; *gmsg << "* Harmonic number h = " << h << " " << endl; /** - if (elptr->getSuperpose()) - *gmsg << "* Fields are superposed " << endl; - */ + if (elptr->getSuperpose()) + *gmsg << "* Fields are superposed " << endl; + */ /** * To ease the initialise() function, set a integral parameter fieldflag internally. @@ -549,13 +548,13 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { double BcParameter[8]; for(int i = 0; i < 8; i++) - BcParameter[i] = 0.0; + BcParameter[i] = 0.0; BcParameter[0] = elptr->getRmin(); BcParameter[1] = elptr->getRmax(); // Store inner radius and outer radius of cyclotron field map in the list - buildupFieldList(BcParameter, "CYCLOTRON", elptr); + buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr); } /** @@ -604,13 +603,13 @@ void ParallelCyclotronTracker::visitCollimator(const Collimator &coll) { double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - std::string ElementType = "CCOLLIMATOR"; + BcParameter[0] = xstart ; BcParameter[1] = xend; BcParameter[2] = ystart ; BcParameter[3] = yend; BcParameter[4] = width ; - buildupFieldList(BcParameter, ElementType, elptr); + buildupFieldList(BcParameter, ElementBase::COLLIMATOR, elptr); } /** @@ -668,11 +667,11 @@ void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) { void ParallelCyclotronTracker::visitOffset(const Offset & off) { if (opalRing_m == NULL) throw OpalException( - "ParallelCylcotronTracker::visitOffset", - "Attempt to place an offset when Ring not defined"); + "ParallelCylcotronTracker::visitOffset", + "Attempt to place an offset when Ring not defined"); Offset* offNonConst = const_cast<Offset*>(&off); offNonConst->updateGeometry(opalRing_m->getNextPosition(), - opalRing_m->getNextNormal()); + opalRing_m->getNextNormal()); opalRing_m->appendElement(off); } @@ -741,7 +740,7 @@ void ParallelCyclotronTracker::visitProbe(const Probe &prob) { double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - std::string ElementType = "PROBE"; + BcParameter[0] = xstart ; BcParameter[1] = xend; BcParameter[2] = ystart ; @@ -749,7 +748,7 @@ void ParallelCyclotronTracker::visitProbe(const Probe &prob) { BcParameter[4] = width ; // store probe parameters in the list - buildupFieldList(BcParameter, ElementType, elptr); + buildupFieldList(BcParameter, ElementBase::PROBE, elptr); } /** @@ -768,7 +767,7 @@ void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) { opalRing_m->appendElement(bend); else throw OpalException("ParallelCyclotronTracker::visitSBend3D", - "Need to define a RINGDEFINITION to use SBend3D element"); + "Need to define a RINGDEFINITION to use SBend3D element"); } void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav) { @@ -777,7 +776,7 @@ void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav opalRing_m->appendElement(cav); else throw OpalException("ParallelCyclotronTracker::visitVariableRFCavity", - "Need to define a RINGDEFINITION to use VariableRFCavity element"); + "Need to define a RINGDEFINITION to use VariableRFCavity element"); } /** @@ -825,7 +824,7 @@ void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) { /* Setup time dependence and in case of no timedependence use a polynom with a_0 = 1 and a_k = 0, k = 1,2,3. - */ + */ std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr; std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr; @@ -836,24 +835,24 @@ void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) { unityVec.push_back(0.); unityVec.push_back(0.); unityVec.push_back(0.); - + if (elptr->getFrequencyModelName() != "") { - freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName()); - *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl; + freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName()); + *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl; } else freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec)); if (elptr->getAmplitudeModelName() != "") { - ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName()); - *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl; + ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName()); + *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl; } else ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec)); if (elptr->getPhaseModelName() != "") { - phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName()); - *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl; + phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName()); + *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl; } else phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec)); @@ -861,18 +860,18 @@ void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) { // read cavity voltage profile data from file. elptr->initialise(itsBunch, 1.0, freq_atd, ampl_atd, phase_atd); -// elptr->initialise(itsBunch, 1.0); + // elptr->initialise(itsBunch, 1.0); double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - std::string ElementType = "CAVITY"; + BcParameter[0] = rmin; BcParameter[1] = rmax; BcParameter[2] = pdis; BcParameter[3] = angle; - buildupFieldList(BcParameter, ElementType, elptr); + buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr); } /** @@ -939,7 +938,7 @@ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) { double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - std::string ElementType = "SEPTUM"; + BcParameter[0] = xstart ; BcParameter[1] = xend; BcParameter[2] = ystart ; @@ -947,7 +946,7 @@ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) { BcParameter[4] = width ; // store septum parameters in the list - buildupFieldList(BcParameter, ElementType, elptr); + buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr); } /** @@ -994,7 +993,7 @@ void ParallelCyclotronTracker::visitCyclotronValley(const CyclotronValley &cv) { * @param scale */ void ParallelCyclotronTracker::applyEntranceFringe(double angle, double curve, - const BMultipoleField &field, double scale) { + const BMultipoleField &field, double scale) { } @@ -1037,7 +1036,7 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) { double BcParameter[8]; for(int i = 0; i < 8; i++) BcParameter[i] = 0.0; - std::string ElementType = "STRIPPER"; + BcParameter[0] = xstart ; BcParameter[1] = xend; BcParameter[2] = ystart ; @@ -1046,12 +1045,12 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) { BcParameter[5] = opcharge; BcParameter[6] = opmass; - buildupFieldList(BcParameter, ElementType, elptr); + buildupFieldList(BcParameter, ElementBase::STRIPPER, elptr); } void ParallelCyclotronTracker::applyExitFringe(double angle, double curve, - const BMultipoleField &field, double scale) { + const BMultipoleField &field, double scale) { } @@ -1063,11 +1062,11 @@ void ParallelCyclotronTracker::applyExitFringe(double angle, double curve, * @param ElementType * @param elptr */ -void ParallelCyclotronTracker::buildupFieldList(double BcParameter[], std::string ElementType, Component *elptr) { +void ParallelCyclotronTracker::buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr) { beamline_list::iterator sindex; type_pair *localpair = new type_pair(); - localpair->first = ElementType; + localpair->first = elementType; for(int i = 0; i < 8; i++) *(((localpair->second).first) + i) = *(BcParameter + i); @@ -1075,7 +1074,7 @@ void ParallelCyclotronTracker::buildupFieldList(double BcParameter[], std::strin (localpair->second).second = elptr; // always put cyclotron as the first element in the list. - if(ElementType == "OPALRING" || ElementType == "CYCLOTRON") { + if(elementType == ElementBase::RING || elementType == ElementBase::CYCLOTRON) { sindex = FieldDimensions.begin(); } else { sindex = FieldDimensions.end(); @@ -1131,7 +1130,7 @@ void ParallelCyclotronTracker::execute() { *gmsg << "* The selected Beam line elements are :" << endl; for(beamline_list::iterator sindex = FieldDimensions.begin(); sindex != FieldDimensions.end(); sindex++) - *gmsg << "* -> " << ((*sindex)->first) << endl; + *gmsg << "* -> " << ElementBase::getTypeString((*sindex)->first) << endl; *gmsg << "* -------------------------------------" << endl; @@ -1230,7 +1229,7 @@ void ParallelCyclotronTracker::Tracker_LF() { if(initialTotalNum_m == 1) openFiles(OpalData::getInstance()->getInputBasename()); - double const initialReferenceTheta = referenceTheta * PIOVER180; + double const initialReferenceTheta = referenceTheta * Physics::deg2rad; initDistInGlobalFrame(); @@ -1553,7 +1552,7 @@ void ParallelCyclotronTracker::Tracker_LF() { // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global aperture bool flagNeedUpdate = deleteParticle(); if(itsBunch->weHaveBins() && flagNeedUpdate) - itsBunch->resetPartBinID2(eta_m); + itsBunch->resetPartBinID2(eta_m); // recalculate bingamma and reset the BinID for each particles according to its current gamma if((itsBunch->weHaveBins()) && BunchCount_m > 1 && step_m % resetBinFreq == 0) @@ -1581,7 +1580,7 @@ void ParallelCyclotronTracker::Tracker_LF() { outfThetaEachTurn_m << "#Turn number = " << turnnumber_m << ", Time = " << itsBunch->getT() * 1e9 << " [ns]" << std::endl; outfThetaEachTurn_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -1589,13 +1588,13 @@ void ParallelCyclotronTracker::Tracker_LF() { // FixMe: should be defined elesewhere ! // define 3 special azimuthal angles where dump particle's six parameters at each turn into 3 ASCII files. const double azimuth_angle0 = 0.0; - const double azimuth_angle1 = 22.5 * PIOVER180; - const double azimuth_angle2 = 45.0 * PIOVER180; + const double azimuth_angle1 = 22.5 * Physics::deg2rad; + const double azimuth_angle2 = 45.0 * Physics::deg2rad; if((oldReferenceTheta < azimuth_angle0 - deltaTheta) && (temp_meanTheta >= azimuth_angle0 - deltaTheta)) { outfTheta0_m << "#Turn number = " << turnnumber_m << ", Time = " << itsBunch->getT() * 1e9 << " [ns]" << std::endl; outfTheta0_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -1605,7 +1604,7 @@ void ParallelCyclotronTracker::Tracker_LF() { outfTheta1_m << "#Turn number = " << turnnumber_m << ", Time = " << itsBunch->getT() * 1e9 << " [ns]" << std::endl; outfTheta1_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -1615,7 +1614,7 @@ void ParallelCyclotronTracker::Tracker_LF() { outfTheta2_m << "#Turn number = " << turnnumber_m << ", Time = " << itsBunch->getT() * 1e9 << " [ns]" << std::endl; outfTheta2_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -1662,14 +1661,14 @@ void ParallelCyclotronTracker::Tracker_LF() { } if (!(itsBunch->getTotalNum()>0)) - throw OpalException("ParallelCyclotronTracker::Tracker_LF()", "No particles anymore, give up now!"); + throw OpalException("ParallelCyclotronTracker::Tracker_LF()", "No particles anymore, give up now!"); for(size_t ii = 0; ii < (itsBunch->getLocalNum()); ii++) { if(itsBunch->ID[ii] == 0) { double FinalMomentum2 = pow(itsBunch->P[ii](0), 2.0) + - pow(itsBunch->P[ii](1), 2.0) + - pow(itsBunch->P[ii](2), 2.0); + pow(itsBunch->P[ii](1), 2.0) + + pow(itsBunch->P[ii](2), 2.0); double FinalEnergy = (sqrt(1.0 + FinalMomentum2) - 1.0) * itsBunch->getM() * 1.0e-6; *gmsg << "* Final energy of reference particle = " << FinalEnergy << " MeV" << endl; *gmsg << "* Total phase space dump number including the initial distribution) = " << lastDumpedStep_m + 1 << endl; @@ -1797,24 +1796,24 @@ void ParallelCyclotronTracker::Tracker_RK4() { for(; (step_m < maxSteps_m) && (itsBunch->getTotalNum()>0); step_m++) { //*gmsg << "step_m= " << step_m << endl; /* - Vector_t B(0., 0., 0.), E(0., 0., 0.); - Vector_t Pt = itsBunch->P[0]; - (*FieldDimensions.begin())->second.second->apply(0, itsBunch->getT()*1e9, E, B); - double PTot = sqrt(Pt(0)*Pt(0)+Pt(1)*Pt(1)+Pt(2)*Pt(2)); - std::cerr << step_m << " P: " << PTot << " "; - std::cerr << "R: "; - for (int jj = 0; jj < 3; jj++) - std::cerr << itsBunch->R[0](jj) << " "; - std::cerr << "B: "; - for (int jj = 0; jj < 3; jj++) - std::cerr << B(jj) << " "; - std::cerr << "E: "; - for (int jj = 0; jj < 3; jj++) - std::cerr << E(jj) << " "; - std::cerr << "P: "; - for (int jj = 0; jj < 3; jj++) - std::cerr << itsBunch->P[0](jj) << " "; - std::cerr << std::endl; + Vector_t B(0., 0., 0.), E(0., 0., 0.); + Vector_t Pt = itsBunch->P[0]; + (*FieldDimensions.begin())->second.second->apply(0, itsBunch->getT()*1e9, E, B); + double PTot = sqrt(Pt(0)*Pt(0)+Pt(1)*Pt(1)+Pt(2)*Pt(2)); + std::cerr << step_m << " P: " << PTot << " "; + std::cerr << "R: "; + for (int jj = 0; jj < 3; jj++) + std::cerr << itsBunch->R[0](jj) << " "; + std::cerr << "B: "; + for (int jj = 0; jj < 3; jj++) + std::cerr << B(jj) << " "; + std::cerr << "E: "; + for (int jj = 0; jj < 3; jj++) + std::cerr << E(jj) << " "; + std::cerr << "P: "; + for (int jj = 0; jj < 3; jj++) + std::cerr << itsBunch->P[0](jj) << " "; + std::cerr << std::endl; */ bool dumpEachTurn = false; @@ -1830,86 +1829,86 @@ void ParallelCyclotronTracker::Tracker_RK4() { int tag = Ippl::Comm->next_tag(IPPL_APP_TAG4, IPPL_APP_CYCLE); - // for all nodes, find the location of particle with ID = 0 & 1 in bunch containers - int found[2] = { -1, -1}; - int counter = 0; + // for all nodes, find the location of particle with ID = 0 & 1 in bunch containers + int found[2] = { -1, -1}; + int counter = 0; - for(size_t ii = 0; ii < (itsBunch->getLocalNum()); ii++) { - if(itsBunch->ID[ii] == 0) { - found[counter] = ii; - counter++; - } - if(itsBunch->ID[ii] == 1) { - found[counter] = ii; - counter++; - } + for(size_t ii = 0; ii < (itsBunch->getLocalNum()); ii++) { + if(itsBunch->ID[ii] == 0) { + found[counter] = ii; + counter++; } - // for the regular modes only the space data of particles with ID = 0 and 1 need be transfored - if(myNode_m == 0) { - // for root node - int notReceived = Ippl::getNodes() - 1; - int numberOfPart = 0; - - while(notReceived > 0) { - int node = COMM_ANY_NODE; - Message *rmsg = Ippl::Comm->receive_block(node, tag); - if(rmsg == 0) - ERRORMSG("Could not receive from client nodes in main." << endl); - notReceived--; - rmsg->get(&numberOfPart); - for(int ii = 0; ii < numberOfPart; ii++) { - rmsg->get(&id); - tmpi.push_back(id); - rmsg->get(&x); - tmpr.push_back(x); - rmsg->get(&x); - tmpr.push_back(x); - rmsg->get(&x); - tmpr.push_back(x); - rmsg->get(&x); - tmpr.push_back(x); - rmsg->get(&x); - tmpr.push_back(x); - rmsg->get(&x); - tmpr.push_back(x); - } - delete rmsg; - + if(itsBunch->ID[ii] == 1) { + found[counter] = ii; + counter++; + } + } + // for the regular modes only the space data of particles with ID = 0 and 1 need be transfored + if(myNode_m == 0) { + // for root node + int notReceived = Ippl::getNodes() - 1; + int numberOfPart = 0; + + while(notReceived > 0) { + int node = COMM_ANY_NODE; + Message *rmsg = Ippl::Comm->receive_block(node, tag); + if(rmsg == 0) + ERRORMSG("Could not receive from client nodes in main." << endl); + notReceived--; + rmsg->get(&numberOfPart); + for(int ii = 0; ii < numberOfPart; ii++) { + rmsg->get(&id); + tmpi.push_back(id); + rmsg->get(&x); + tmpr.push_back(x); + rmsg->get(&x); + tmpr.push_back(x); + rmsg->get(&x); + tmpr.push_back(x); + rmsg->get(&x); + tmpr.push_back(x); + rmsg->get(&x); + tmpr.push_back(x); + rmsg->get(&x); + tmpr.push_back(x); } - for(int ii = 0; ii < 1; ii++) { - tmpi.push_back(itsBunch->ID[found[ii]]); - for(int jj = 0; jj < 3; jj++) { - tmpr.push_back(itsBunch->R[found[ii]](jj)); - tmpr.push_back(itsBunch->P[found[ii]](jj)); - } + delete rmsg; + + } + for(int ii = 0; ii < 1; ii++) { + tmpi.push_back(itsBunch->ID[found[ii]]); + for(int jj = 0; jj < 3; jj++) { + tmpr.push_back(itsBunch->R[found[ii]](jj)); + tmpr.push_back(itsBunch->P[found[ii]](jj)); } - std::vector<double>::iterator itParameter = tmpr.begin(); - std::vector<int>::iterator itId = tmpi.begin(); - - for(itId = tmpi.begin(); itId != tmpi.end(); itId++) { - outfTrackOrbit_m << "ID" << *itId; - for(int ii = 0; ii < 12; ii++) { - outfTrackOrbit_m << " " << *itParameter; - itParameter++; - } - outfTrackOrbit_m << std::endl; + } + std::vector<double>::iterator itParameter = tmpr.begin(); + std::vector<int>::iterator itId = tmpi.begin(); + + for(itId = tmpi.begin(); itId != tmpi.end(); itId++) { + outfTrackOrbit_m << "ID" << *itId; + for(int ii = 0; ii < 12; ii++) { + outfTrackOrbit_m << " " << *itParameter; + itParameter++; } - // sample frequency = SinglePartDumpFreq - } else { - // for other nodes - Message *smsg = new Message(); - smsg->put(counter); - for(int ii = 0; ii < counter; ii++) { - smsg->put(itsBunch->ID[found[ii]]); - for(int jj = 0; jj < 3; jj++) { - smsg->put(itsBunch->R[found[ii]](jj)); - smsg->put(itsBunch->P[found[ii]](jj)); - } + outfTrackOrbit_m << std::endl; + } + // sample frequency = SinglePartDumpFreq + } else { + // for other nodes + Message *smsg = new Message(); + smsg->put(counter); + for(int ii = 0; ii < counter; ii++) { + smsg->put(itsBunch->ID[found[ii]]); + for(int jj = 0; jj < 3; jj++) { + smsg->put(itsBunch->R[found[ii]](jj)); + smsg->put(itsBunch->P[found[ii]](jj)); } - bool res = Ippl::Comm->send(smsg, 0, tag); - if(!res) - ERRORMSG("Ippl::Comm->send(smsg, 0, tag) failed " << endl); } + bool res = Ippl::Comm->send(smsg, 0, tag); + if(!res) + ERRORMSG("Ippl::Comm->send(smsg, 0, tag) failed " << endl); + } IpplTimings::stopTimer(DumpTimer_m); } @@ -2293,6 +2292,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { bgf_main_collision_test(); // New // **************************************************************************************** + IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < (itsBunch->getLocalNum()); i++) { flagNoDeletion = true; // change phase space parameters from localframe of bunch (dr,dtheta,dz) to global Cartesian frame (X,Y,Z) @@ -2319,7 +2319,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { bool tag_crossing = false; double DistOld = 0.0; //mm RFCavity * rfcav; - if(((*sindex)->first) == "CAVITY") { + if(((*sindex)->first) == ElementBase::RFCAVITY) { // here check gap cross in the list, if do , set tag_crossing to TRUE for(int j = 0; j < 3; j++) rnew_m[j] = variable_m[j]; @@ -2331,7 +2331,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { double oldBetgam = sqrt(oldMomentum2); double oldGamma = sqrt(1.0 + oldMomentum2); double oldBeta = oldBetgam / oldGamma; - double dt1 = DistOld / (c * oldBeta * 1.0e-6); // ns + double dt1 = DistOld / (Physics::c * oldBeta * 1.0e-6); // ns double dt2 = dt - dt1; // retrack particle from the old postion to cavity gap point @@ -2366,6 +2366,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { } // end if: gap-crossing monentum kicking at certain cavity } //end for: finish checking for all cavities } //end for: finish one step tracking for all particles for initialTotalNum_m != 2 mode + IpplTimings::stopTimer(IntegrationTimer_m); // *** This has to go before the actual timestep -DW ****************************** // apply the plugin elements: probe, collimator, stripper, septum @@ -2379,11 +2380,11 @@ void ParallelCyclotronTracker::Tracker_RK4() { bool flagNeedUpdate = deleteParticle(); if(itsBunch->weHaveBins() && flagNeedUpdate) - itsBunch->resetPartBinID2(eta_m); + itsBunch->resetPartBinID2(eta_m); // recalculate bingamma and reset the BinID for each particles according to its current gamma if((itsBunch->weHaveBins()) && BunchCount_m > 1 && step_m % resetBinFreq == 0) - itsBunch->resetPartBinID2(eta_m); + itsBunch->resetPartBinID2(eta_m); if((step_m > 10) && ((step_m + 1) % stepsPerTurn) == 0) { ++turnnumber_m; @@ -2403,6 +2404,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { // trigger SEO mode (swith off cavity) and calculate betatron osciliation tuning. double r_tuning[2], z_tuning[2] ; + IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < (itsBunch->getLocalNum()); i++) { // change phase space parameters from local frame of bunch (dr,dtheta,dz) to global Cartesian frame (X,Y,Z) @@ -2411,7 +2413,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { if((step_m % SinglePartDumpFreq == 0)) { outfTrackOrbit_m << "ID" << (itsBunch->ID[i]); outfTrackOrbit_m << " " << variable_m[0] << " " << variable_m[3] << " " << variable_m[1] << " " - << variable_m[4] << " " << variable_m[2] << " " << variable_m[5] << std::endl; + << variable_m[4] << " " << variable_m[2] << " " << variable_m[5] << std::endl; } double OldTheta = 0.0; @@ -2429,6 +2431,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { for(int j = 0; j < 3; j++) itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: dimensionless, beta*gama }//end for: finish one step tracking for all particles for initialTotalNum_m != 2 mode + IpplTimings::stopTimer(IntegrationTimer_m); // store dx and dz for future tune calculation if higher precision needed, reduce freqSample. if(step_m % SinglePartDumpFreq == 0) { @@ -2439,167 +2442,169 @@ void ParallelCyclotronTracker::Tracker_RK4() { } } else if(initialTotalNum_m == 1) { - // initialTotalNum_m == 1 trigger single particle mode + // initialTotalNum_m == 1 trigger single particle mode - IpplTimings::startTimer(IntegrationTimer_m); - flagNoDeletion = true; + IpplTimings::startTimer(IntegrationTimer_m); + flagNoDeletion = true; - // *** This was moved here b/c collision should be tested before the ********************** - // *** actual timestep (bgf_main_collision_test() predicts the next step automatically) -DW - // apply the plugin elements: probe, collimator, stripper, septum - applyPluginElements(dt); // New + // *** This was moved here b/c collision should be tested before the ********************** + // *** actual timestep (bgf_main_collision_test() predicts the next step automatically) -DW + // apply the plugin elements: probe, collimator, stripper, septum + applyPluginElements(dt); // New - // check if we loose particles at the boundary - bgf_main_collision_test(); // New - // **************************************************************************************** + // check if we loose particles at the boundary + bgf_main_collision_test(); // New + // **************************************************************************************** - // track for one step - for( unsigned int i = 0; i < itsBunch->getLocalNum(); i++) { + // track for one step + IpplTimings::startTimer(IntegrationTimer_m); + for( unsigned int i = 0; i < itsBunch->getLocalNum(); i++) { - if((step_m % SinglePartDumpFreq == 0)) { - outfTrackOrbit_m << "ID" <<itsBunch->ID[i]; - outfTrackOrbit_m << " " << itsBunch->R[i](0) << " " <<itsBunch->P[i](0) << " " <<itsBunch->R[i](1) - << " " << itsBunch->P[i](1) << " " <<itsBunch->R[i](2) << " " <<itsBunch->P[i](2) << std::endl; - } + if((step_m % SinglePartDumpFreq == 0)) { + outfTrackOrbit_m << "ID" <<itsBunch->ID[i]; + outfTrackOrbit_m << " " << itsBunch->R[i](0) << " " <<itsBunch->P[i](0) << " " <<itsBunch->R[i](1) + << " " << itsBunch->P[i](1) << " " <<itsBunch->R[i](2) << " " <<itsBunch->P[i](2) << std::endl; + } - // change phase space parameters from local frame of bunch (dr,dtheta,dz) to global Cartesian frame (X,Y,Z) - for(int j = 0; j < 3; j++) { - variable_m[j] = itsBunch->R[i](j); //[x,y,z] units: [mm] - variable_m[j+3] = itsBunch->P[i](j); //[px,py,pz] units: dimensionless - rold_m[j] = variable_m[j]; // used for gap cross checking - pold_m[j] = variable_m[j+3]; // used for gap cross - } + // change phase space parameters from local frame of bunch (dr,dtheta,dz) to global Cartesian frame (X,Y,Z) + for(int j = 0; j < 3; j++) { + variable_m[j] = itsBunch->R[i](j); //[x,y,z] units: [mm] + variable_m[j+3] = itsBunch->P[i](j); //[px,py,pz] units: dimensionless + rold_m[j] = variable_m[j]; // used for gap cross checking + pold_m[j] = variable_m[j+3]; // used for gap cross + } - double temp_meanTheta = calculateAngle2(variable_m[0], variable_m[1]);//[ -pi ~ pi ] + double temp_meanTheta = calculateAngle2(variable_m[0], variable_m[1]);//[ -pi ~ pi ] - if((step_m > 10) && ((step_m + 1) % stepsPerTurn) == 0) { - ++turnnumber_m; - dumpEachTurn = true; - *gmsg << "Turn " << turnnumber_m << endl; - - outfThetaEachTurn_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; - outfThetaEachTurn_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) - << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / pi * 180 - << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) - << " " << variable_m[2] - << " " << variable_m[5] << std::endl; - } + if((step_m > 10) && ((step_m + 1) % stepsPerTurn) == 0) { + ++turnnumber_m; + dumpEachTurn = true; + *gmsg << "Turn " << turnnumber_m << endl; - //define 3 special azimuthal angles where dump particle's six parameters at each turn into 3 ASCII files. - const double azimuth_angle0 = 0.0; - const double azimuth_angle1 = 22.5 / 180.0 * pi; - const double azimuth_angle2 = 45.0 / 180.0 * pi; + outfThetaEachTurn_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; + outfThetaEachTurn_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) + << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) + << " " << temp_meanTheta / pi * 180 + << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) + << " " << variable_m[2] + << " " << variable_m[5] << std::endl; + } - if((oldReferenceTheta < azimuth_angle0 - deltaTheta) && (temp_meanTheta >= azimuth_angle0 - deltaTheta)) { - outfTheta0_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; - outfTheta0_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) - << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / pi * 180 - << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) - << " " << variable_m[2] - << " " << variable_m[5] << std::endl; - } + //define 3 special azimuthal angles where dump particle's six parameters at each turn into 3 ASCII files. + const double azimuth_angle0 = 0.0; + const double azimuth_angle1 = 22.5 / 180.0 * pi; + const double azimuth_angle2 = 45.0 / 180.0 * pi; - if((oldReferenceTheta < azimuth_angle1 - deltaTheta) && (temp_meanTheta >= azimuth_angle1 - deltaTheta)) { - outfTheta1_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; - outfTheta1_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) - << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / pi * 180 - << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) - << " " << variable_m[2] - << " " << variable_m[5] << std::endl; - } + if((oldReferenceTheta < azimuth_angle0 - deltaTheta) && (temp_meanTheta >= azimuth_angle0 - deltaTheta)) { + outfTheta0_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; + outfTheta0_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) + << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) + << " " << temp_meanTheta / pi * 180 + << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) + << " " << variable_m[2] + << " " << variable_m[5] << std::endl; + } - if((oldReferenceTheta < azimuth_angle2 - deltaTheta) && (temp_meanTheta >= azimuth_angle2 - deltaTheta)) { - outfTheta2_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; - outfTheta2_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) - << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) - << " " << temp_meanTheta / pi * 180 - << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) - << " " << variable_m[2] - << " " << variable_m[5] << std::endl; - } + if((oldReferenceTheta < azimuth_angle1 - deltaTheta) && (temp_meanTheta >= azimuth_angle1 - deltaTheta)) { + outfTheta1_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; + outfTheta1_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) + << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) + << " " << temp_meanTheta / pi * 180 + << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) + << " " << variable_m[2] + << " " << variable_m[5] << std::endl; + } - oldReferenceTheta = temp_meanTheta; + if((oldReferenceTheta < azimuth_angle2 - deltaTheta) && (temp_meanTheta >= azimuth_angle2 - deltaTheta)) { + outfTheta2_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; + outfTheta2_m << " " << sqrt(variable_m[0]*variable_m[0] + variable_m[1]*variable_m[1]) + << " " << variable_m[3]*cos(temp_meanTheta) + variable_m[4]*sin(temp_meanTheta) + << " " << temp_meanTheta / pi * 180 + << " " << -variable_m[3]*sin(temp_meanTheta) + variable_m[4]*cos(temp_meanTheta) + << " " << variable_m[2] + << " " << variable_m[5] << std::endl; + } - // integrate for one step in the lab Cartesian frame (absolute value). - flagNoDeletion = rk4(variable_m, t, dt, i); + oldReferenceTheta = temp_meanTheta; - if(!flagNoDeletion) { - *gmsg << "particle" << "is lost at " << step_m << "th step!" << endl; - throw OpalException("ParallelCyclotronTracker", "the particle is out of the region of interest."); - } + // integrate for one step in the lab Cartesian frame (absolute value). + flagNoDeletion = rk4(variable_m, t, dt, i); - for(int j = 0; j < 3; j++) itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] - for(int j = 0; j < 3; j++) itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [] beta*gama - - //If gap crossing happens, do momenta kicking - - for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); sindex++) { - bool tag_crossing = false; - double DistOld = 0.0; //mm - RFCavity * rfcav; - if(((*sindex)->first) == "CAVITY") { - // here check gap cross in the list, if do , set tag_crossing to TRUE - for(int j = 0; j < 3; j++) - rnew_m[j] = variable_m[j]; - rfcav = static_cast<RFCavity *>(((*sindex)->second).second); - tag_crossing = checkGapCross(rold_m, rnew_m, rfcav, DistOld); - } - if(tag_crossing) { - double oldMomentum2 = dot(pold_m, pold_m); - double oldBetgam = sqrt(oldMomentum2); - double oldGamma = sqrt(1.0 + oldMomentum2); - double oldBeta = oldBetgam / oldGamma; - double dt1 = DistOld / (c * oldBeta * 1.0e-6); // ns - double dt2 = dt - dt1; - - // retrack particle from the old postion to cavity gap point - // restore the old coordinates and momenta - for(int j = 0; j < 3; j++) { - variable_m[j] = rold_m[j]; - variable_m[j+3] = pold_m[j]; + if(!flagNoDeletion) { + *gmsg << "particle" << "is lost at " << step_m << "th step!" << endl; + throw OpalException("ParallelCyclotronTracker", "the particle is out of the region of interest."); } - if(dt / dt1 < 1.0e9) rk4(variable_m, t, dt1, i); + for(int j = 0; j < 3; j++) itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] + for(int j = 0; j < 3; j++) itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [] beta*gama - for(int j = 0; j < 3; j++) { - itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] - itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [] beta*gama - } + //If gap crossing happens, do momenta kicking - //momentum kick - RFkick(rfcav, t, dt1, i); + for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); sindex++) { + bool tag_crossing = false; + double DistOld = 0.0; //mm + RFCavity * rfcav; + if(((*sindex)->first) == ElementBase::RFCAVITY) { + // here check gap cross in the list, if do , set tag_crossing to TRUE + for(int j = 0; j < 3; j++) + rnew_m[j] = variable_m[j]; + rfcav = static_cast<RFCavity *>(((*sindex)->second).second); + tag_crossing = checkGapCross(rold_m, rnew_m, rfcav, DistOld); + } + if(tag_crossing) { + double oldMomentum2 = dot(pold_m, pold_m); + double oldBetgam = sqrt(oldMomentum2); + double oldGamma = sqrt(1.0 + oldMomentum2); + double oldBeta = oldBetgam / oldGamma; + double dt1 = DistOld / (Physics::c * oldBeta * 1.0e-6); // ns + double dt2 = dt - dt1; - // retrack particle from cavity gap point for the left time to finish the entire timestep - for(int j = 0; j < 3; j++) { - variable_m[j] = itsBunch->R[i](j); //[x,y,z] units: [mm] - variable_m[j+3] = itsBunch->P[i](j); //[px,py,pz] units: [] - } + // retrack particle from the old postion to cavity gap point + // restore the old coordinates and momenta + for(int j = 0; j < 3; j++) { + variable_m[j] = rold_m[j]; + variable_m[j+3] = pold_m[j]; + } - if(dt / dt2 < 1.0e9) rk4(variable_m, t, dt2, i); + if(dt / dt1 < 1.0e9) rk4(variable_m, t, dt1, i); - for(int j = 0; j < 3; j++) { - itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] - itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [], beta*gama - } - }// end if: gap-crossing monentum kicking at certain cavity - }//end for: finish checking for all cavities - } + for(int j = 0; j < 3; j++) { + itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] + itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [] beta*gama + } - // *** This has to go before the actual timestep -DW ****************************** - // apply the plugin elements: probe, collimator, stripper, septum - //applyPluginElements(dt); // Old + //momentum kick + RFkick(rfcav, t, dt1, i); + + // retrack particle from cavity gap point for the left time to finish the entire timestep + for(int j = 0; j < 3; j++) { + variable_m[j] = itsBunch->R[i](j); //[x,y,z] units: [mm] + variable_m[j+3] = itsBunch->P[i](j); //[px,py,pz] units: [] + } - // check if we loose particles at the boundary - //bgf_main_collision_test(); // Old - // ******************************************************************************** + if(dt / dt2 < 1.0e9) rk4(variable_m, t, dt2, i); - // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture - deleteParticle(); + for(int j = 0; j < 3; j++) { + itsBunch->R[i](j) = variable_m[j] ; //[x,y,z] units: [mm] + itsBunch->P[i](j) = variable_m[j+3] ; //[px,py,pz] units: [], beta*gama + } + }// end if: gap-crossing monentum kicking at certain cavity + }//end for: finish checking for all cavities + } + IpplTimings::stopTimer(IntegrationTimer_m); - IpplTimings::stopTimer(IntegrationTimer_m); + // *** This has to go before the actual timestep -DW ****************************** + // apply the plugin elements: probe, collimator, stripper, septum + //applyPluginElements(dt); // Old + + // check if we loose particles at the boundary + //bgf_main_collision_test(); // Old + // ******************************************************************************** + + // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture + deleteParticle(); + + IpplTimings::stopTimer(IntegrationTimer_m); }//end if: finish one step tracking either for initialTotalNum_m==2 || initialTotalNum_m==2 || initialTotalNum_m==1 mode // update bunch and some parameters and output some info. after one time step. @@ -2632,7 +2637,7 @@ void ParallelCyclotronTracker::Tracker_RK4() { } // end for: the integration is DONE after maxSteps_m steps! if (!(itsBunch->getTotalNum()>0)) - throw OpalException("ParallelCyclotronTracker::Tracker_MTS()", "No particles anymore, give up now!"); + throw OpalException("ParallelCyclotronTracker::Tracker_MTS()", "No particles anymore, give up now!"); // some post-integration works *gmsg << endl; @@ -2680,9 +2685,9 @@ void ParallelCyclotronTracker::Tracker_RK4() { #ifdef GENERICTRACKER /** - * - * - */ +* +* +*/ void ParallelCyclotronTracker::Tracker_Generic() { // Generic Tracker that has three modes defined by timeIntegrator_m: // 0 --- RK-4 (default) @@ -2704,7 +2709,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { int scSolveFreq; if (spiral_flag) - scSolveFreq = 1; + scSolveFreq = 1; else scSolveFreq = Options::scSolveFreq; @@ -2720,10 +2725,10 @@ void ParallelCyclotronTracker::Tracker_Generic() { // Define 3 special azimuthal angles where dump particle's six parameters // at each turn into 3 ASCII files. only important in single particle tracking const double azimuth_angle0 = 0.0; - const double azimuth_angle1 = 22.5 * PIOVER180; - const double azimuth_angle2 = 45.0 * PIOVER180; + const double azimuth_angle1 = 22.5 * Physics::deg2rad; + const double azimuth_angle2 = 45.0 * Physics::deg2rad; - const double initialReferenceTheta = referenceTheta * PIOVER180; + const double initialReferenceTheta = referenceTheta * Physics::deg2rad; double oldReferenceTheta = initialReferenceTheta; // init here, reset each step const double deltaTheta = pi / (stepsPerTurn); // half of the average angle per step @@ -2839,13 +2844,13 @@ void ParallelCyclotronTracker::Tracker_Generic() { // If we are using the 2nd order Leap Frog method, push half a step /* TEMP move this right before the second half step for testing -DW - if (timeIntegrator_m == 1) { + if (timeIntegrator_m == 1) { - itsBunch->R *= Vector_t(0.001); // mm --> m - push(0.5 * dt * 1.0e-9); // ns --> s - itsBunch->R *= Vector_t(1000.0); // m --> mm + itsBunch->R *= Vector_t(0.001); // mm --> m + push(0.5 * dt * 1.0e-9); // ns --> s + itsBunch->R *= Vector_t(1000.0); // m --> mm - } + } */ // Find out if we need to do bunch injection @@ -2911,10 +2916,10 @@ void ParallelCyclotronTracker::Tracker_Generic() { // read initial distribution from h5 file if (multiBunchMode_m == 1) { - + if (BunchCount_m == 2) saveOneBunch(); - + readOneBunch(BunchCount_m - 1); if (timeIntegrator_m == 0) itsBunch->resetPartBinID2(eta_m); @@ -3019,60 +3024,60 @@ void ParallelCyclotronTracker::Tracker_Generic() { } else { // --- Single bunch mode --- // - - // If we are doing a spiral inflector simulation and are using the SAAMG solver + + // If we are doing a spiral inflector simulation and are using the SAAMG solver // we don't rotate or translate the bunch and gamma is 1.0 (non-relativistic). if (spiral_flag and itsBunch->getFieldSolverType() == "SAAMG") { - itsBunch->R *= Vector_t(0.001); // mm --> m + itsBunch->R *= Vector_t(0.001); // mm --> m - IpplTimings::stopTimer(TransformTimer_m); + IpplTimings::stopTimer(TransformTimer_m); - itsBunch->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0)); - itsBunch->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0)); + itsBunch->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0)); + itsBunch->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0)); - itsBunch->computeSelfFields_cycl(1.0); + itsBunch->computeSelfFields_cycl(1.0); - IpplTimings::startTimer(TransformTimer_m); + IpplTimings::startTimer(TransformTimer_m); - itsBunch->R *= Vector_t(1000.0); // m --> mm + itsBunch->R *= Vector_t(1000.0); // m --> mm } else { - double temp_meangamma = sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP)); + double temp_meangamma = sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP)); - Quaternion_t quaternionToYAxis; + Quaternion_t quaternionToYAxis; - getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis); + getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis); - globalToLocal(itsBunch->R, quaternionToYAxis, meanR); + globalToLocal(itsBunch->R, quaternionToYAxis, meanR); - itsBunch->R *= Vector_t(0.001); // mm --> m + itsBunch->R *= Vector_t(0.001); // mm --> m - if((step_m + 1) % boundpDestroyFreq == 0) - itsBunch->boundp_destroy(); - else - itsBunch->boundp(); + if((step_m + 1) % boundpDestroyFreq == 0) + itsBunch->boundp_destroy(); + else + itsBunch->boundp(); - IpplTimings::stopTimer(TransformTimer_m); + IpplTimings::stopTimer(TransformTimer_m); - repartition(); + repartition(); - itsBunch->setGlobalMeanR(0.001 * meanR); - itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis); + itsBunch->setGlobalMeanR(0.001 * meanR); + itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis); - itsBunch->computeSelfFields_cycl(temp_meangamma); + itsBunch->computeSelfFields_cycl(temp_meangamma); - IpplTimings::startTimer(TransformTimer_m); + IpplTimings::startTimer(TransformTimer_m); - //scale coordinates back - itsBunch->R *= Vector_t(1000.0); // m --> mm + //scale coordinates back + itsBunch->R *= Vector_t(1000.0); // m --> mm - // Transform coordinates back to global - localToGlobal(itsBunch->R, quaternionToYAxis, meanR); + // Transform coordinates back to global + localToGlobal(itsBunch->R, quaternionToYAxis, meanR); - // Transform self field back to global frame (rotate only) - localToGlobal(itsBunch->Ef, quaternionToYAxis); - localToGlobal(itsBunch->Bf, quaternionToYAxis); + // Transform self field back to global frame (rotate only) + localToGlobal(itsBunch->Ef, quaternionToYAxis); + localToGlobal(itsBunch->Bf, quaternionToYAxis); } } @@ -3126,6 +3131,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { } + IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < (itsBunch->getLocalNum()); i++) { flagNoDeletion = true; @@ -3179,7 +3185,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { double DistOld = 0.0; //mm RFCavity * rfcav; - if(((*sindex)->first) == "CAVITY") { + if(((*sindex)->first) == ElementBase::RFCAVITY) { // here check gap cross in the list, if do , set tag_crossing to TRUE for(int j = 0; j < 3; j++) @@ -3195,7 +3201,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { double oldBetgam = sqrt(oldMomentum2); double oldGamma = sqrt(1.0 + oldMomentum2); double oldBeta = oldBetgam / oldGamma; - double dt1 = DistOld / (c * oldBeta * 1.0e-6); // ns + double dt1 = DistOld / (Physics::c * oldBeta * 1.0e-6); // ns double dt2 = dt - dt1; // retrack particle from the old postion to cavity gap point @@ -3231,6 +3237,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { } // end for: finish checking for all cavities } // end if: which mode are we in LF2 or RK4 } // end for: finish one step tracking for all particles for initialTotalNum_m != 2 mode + IpplTimings::stopTimer(IntegrationTimer_m); if (timeIntegrator_m == 1) { // LF-2 method @@ -3284,6 +3291,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { } + IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < (itsBunch->getLocalNum()); i++) { for(int j = 0; j < 3; j++) variable_m[j] = itsBunch->R[i](j); //[x,y,z] (mm) @@ -3315,7 +3323,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { getFieldsAtPoint(t, i, externalE, externalB); - IpplTimings::startTimer(IntegrationTimer_m); + // IpplTimings::startTimer(IntegrationTimer_m); // Do LF2 momentum kick pusher.kick(itsBunch->R[i], @@ -3326,7 +3334,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { itsBunch->M[i] * 1.0e9, itsBunch->Q[i] / q_e); - IpplTimings::stopTimer(IntegrationTimer_m); + // IpplTimings::stopTimer(IntegrationTimer_m); } @@ -3343,6 +3351,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { if( (i == 0) && (step_m > 10) && ((step_m%stepsPerTurn) == 0)) turnnumber_m++; } // end for: finish one step tracking for all particles for initialTotalNum_m == 2 mode + IpplTimings::stopTimer(IntegrationTimer_m); if (timeIntegrator_m == 1) { // LF-2 method @@ -3392,6 +3401,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { } + IpplTimings::startTimer(IntegrationTimer_m); for( unsigned int i = 0; i < itsBunch->getLocalNum(); i++) { if((step_m % SinglePartDumpFreq == 0)) { @@ -3422,7 +3432,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { outfThetaEachTurn_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; outfThetaEachTurn_m << " " << sqrt(variable_m[0] * variable_m[0] + variable_m[1] * variable_m[1]) << " " << variable_m[3] * cos(temp_meanTheta) + variable_m[4] * sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3] * sin(temp_meanTheta) + variable_m[4] * cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -3432,7 +3442,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { outfTheta0_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; outfTheta0_m << " " << sqrt(variable_m[0] * variable_m[0] + variable_m[1] * variable_m[1]) << " " << variable_m[3] * cos(temp_meanTheta) + variable_m[4] * sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3] * sin(temp_meanTheta) + variable_m[4] * cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -3442,7 +3452,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { outfTheta1_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; outfTheta1_m << " " << sqrt(variable_m[0] * variable_m[0] + variable_m[1] * variable_m[1]) << " " << variable_m[3] * cos(temp_meanTheta) + variable_m[4] * sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3] * sin(temp_meanTheta) + variable_m[4] * cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -3452,7 +3462,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { outfTheta2_m << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl; outfTheta2_m << " " << sqrt(variable_m[0] * variable_m[0] + variable_m[1] * variable_m[1]) << " " << variable_m[3] * cos(temp_meanTheta) + variable_m[4] * sin(temp_meanTheta) - << " " << temp_meanTheta / PIOVER180 + << " " << temp_meanTheta / Physics::deg2rad << " " << -variable_m[3] * sin(temp_meanTheta) + variable_m[4] * cos(temp_meanTheta) << " " << variable_m[2] << " " << variable_m[5] << std::endl; @@ -3469,11 +3479,11 @@ void ParallelCyclotronTracker::Tracker_Generic() { if (outOfBound) flagNoDeletion = false; // Do LF2 momentum kick - IpplTimings::startTimer(IntegrationTimer_m); + // IpplTimings::startTimer(IntegrationTimer_m); pusher.kick(itsBunch->R[i], itsBunch->P[i], externalE, externalB, dt * 1.0e-9, itsBunch->M[i] * 1.0e9, itsBunch->Q[i] / q_e); - IpplTimings::stopTimer(IntegrationTimer_m); + // IpplTimings::stopTimer(IntegrationTimer_m); if(!flagNoDeletion) { @@ -3506,7 +3516,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { double DistOld = 0.0; //mm RFCavity * rfcav; - if (((*sindex)->first) == "CAVITY") { + if (((*sindex)->first) == ElementBase::RFCAVITY) { // here check gap cross in the list, if do , set tag_crossing to TRUE for(int j = 0; j < 3; j++) @@ -3522,7 +3532,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { double oldBetgam = sqrt(oldMomentum2); double oldGamma = sqrt(1.0 + oldMomentum2); double oldBeta = oldBetgam / oldGamma; - double dt1 = DistOld / (c * oldBeta * 1.0e-6); // ns + double dt1 = DistOld / (Physics::c * oldBeta * 1.0e-6); // ns double dt2 = dt - dt1; // retrack particle from the old postion to cavity gap point @@ -3562,6 +3572,7 @@ void ParallelCyclotronTracker::Tracker_Generic() { } // end for: finish checking for all cavities } // end if: LF2 or RK4 method } // end for: track one particle one step + IpplTimings::stopTimer(IntegrationTimer_m); if (timeIntegrator_m == 1) { // LF-2 method @@ -3600,14 +3611,14 @@ void ParallelCyclotronTracker::Tracker_Generic() { if (itsBunch->getTotalNum()>0) { // Only dump last step if we have particles left. // Check separately for phase space (ps) and statistics (stat) data dump frequency if((((step_m + 1) % Options::psDumpFreq == 0) && initialTotalNum_m != 2) - || (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) { + || (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) { // Write phase space data to h5 file bunchDumpPhaseSpaceData(); - } + } if((((step_m + 1) % Options::statDumpFreq == 0) && initialTotalNum_m != 2) - || (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) { + || (doDumpAfterEachTurn && dumpEachTurn && initialTotalNum_m != 2)) { // Write statistics data to stat file bunchDumpStatData(); @@ -3757,18 +3768,18 @@ bool ParallelCyclotronTracker::derivate(double *y, const double &t, double *yp, double tempgamma = sqrt(1 + (y[3] * y[3] + y[4] * y[4] + y[5] * y[5])); /* - d2x/dt2 = q/m * ( E + v x B ) - dx/dt =vx - units: mm, and [] + d2x/dt2 = q/m * ( E + v x B ) + dx/dt =vx + units: mm, and [] */ yp[0] = c_mmtns / tempgamma * y[3]; // [mn/ns] yp[1] = c_mmtns / tempgamma * y[4]; // [mm/ns] yp[2] = c_mmtns / tempgamma * y[5]; // [mm/ns] - yp[3] = (externalE(0) / c + (externalB(2) * y[4] - externalB(1) * y[5]) / tempgamma) * qtom; // [1/ns] - yp[4] = (externalE(1) / c - (externalB(2) * y[3] - externalB(0) * y[5]) / tempgamma) * qtom; // [1/ns]; - yp[5] = (externalE(2) / c + (externalB(1) * y[3] - externalB(0) * y[4]) / tempgamma) * qtom; // [1/ns]; + yp[3] = (externalE(0) / Physics::c + (externalB(2) * y[4] - externalB(1) * y[5]) / tempgamma) * qtom; // [1/ns] + yp[4] = (externalE(1) / Physics::c - (externalB(2) * y[3] - externalB(0) * y[5]) / tempgamma) * qtom; // [1/ns]; + yp[5] = (externalE(2) / Physics::c + (externalB(1) * y[3] - externalB(0) * y[4]) / tempgamma) * qtom; // [1/ns]; return outOfBound; } @@ -3792,8 +3803,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta // tau Step size (usually time step) // Pindex index of particel, not used yet - IpplTimings::startTimer(IntegrationTimer_m); - bool outOfBound = false; double deriv1[PSdim]; double deriv2[PSdim]; @@ -3805,7 +3814,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta outOfBound = derivate(x, t, deriv1 , Pindex); if (outOfBound) { - IpplTimings::stopTimer(IntegrationTimer_m); return false; } @@ -3818,7 +3826,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta outOfBound = derivate(xtemp, t_half, deriv2 , Pindex); if (outOfBound) { - IpplTimings::stopTimer(IntegrationTimer_m); return false; } @@ -3828,7 +3835,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta outOfBound = derivate(xtemp, t_half, deriv3 , Pindex); if (outOfBound) { - IpplTimings::stopTimer(IntegrationTimer_m); return false; } @@ -3839,7 +3845,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta outOfBound = derivate(xtemp, t_full, deriv4 , Pindex); if (outOfBound) { - IpplTimings::stopTimer(IntegrationTimer_m); return false; } @@ -3847,7 +3852,6 @@ bool ParallelCyclotronTracker::rk4(double x[], const double &t, const double &ta for(int i = 0; i < PSdim; i++) x[i] += tau / 6.*(deriv1[i] + deriv4[i] + 2.*(deriv2[i] + deriv3[i])); - IpplTimings::stopTimer(IntegrationTimer_m); return true; } @@ -3876,25 +3880,25 @@ bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew, RFCav } bool ParallelCyclotronTracker::RFkick(RFCavity * rfcavity, const double t, const double dt, const int Pindex){ - double radius = sqrt(pow(itsBunch->R[Pindex](0), 2.0) + pow(itsBunch->R[Pindex](1), 2.0) - - pow(rfcavity->getPerpenDistance() , 2.0)); - double rmin = rfcavity->getRmin(); - double rmax = rfcavity->getRmax(); - double nomalRadius = (radius - rmin) / (rmax - rmin); - double tempP[3]; - if(nomalRadius <= 1.0 && nomalRadius >= 0.0) { - - for(int j = 0; j < 3; j++) - tempP[j] = itsBunch->P[Pindex](j); //[px,py,pz] units: dimensionless - - // here evaluate voltage and conduct momenta kicking; - rfcavity -> getMomentaKick(nomalRadius, tempP, t, dt, itsBunch->ID[Pindex], itsBunch->getM(), itsBunch->getQ()); // t : ns - - for(int j = 0; j < 3; j++) - itsBunch->P[Pindex](j) = tempP[j]; - return true; - } - return false; + double radius = sqrt(pow(itsBunch->R[Pindex](0), 2.0) + pow(itsBunch->R[Pindex](1), 2.0) + - pow(rfcavity->getPerpenDistance() , 2.0)); + double rmin = rfcavity->getRmin(); + double rmax = rfcavity->getRmax(); + double nomalRadius = (radius - rmin) / (rmax - rmin); + double tempP[3]; + if(nomalRadius <= 1.0 && nomalRadius >= 0.0) { + + for(int j = 0; j < 3; j++) + tempP[j] = itsBunch->P[Pindex](j); //[px,py,pz] units: dimensionless + + // here evaluate voltage and conduct momenta kicking; + rfcavity -> getMomentaKick(nomalRadius, tempP, t, dt, itsBunch->ID[Pindex], itsBunch->getM(), itsBunch->getQ()); // t : ns + + for(int j = 0; j < 3; j++) + itsBunch->P[Pindex](j) = tempP[j]; + return true; + } + return false; } @@ -4152,13 +4156,13 @@ double ParallelCyclotronTracker::getHarmonicNumber() const { if (elcycl != NULL) return elcycl->getCyclHarm(); throw OpalException("ParallelCyclotronTracker::Tracker_MTS()", - std::string("The first item in the FieldDimensions list does not ") - +std::string("seem to be an Ring or a Cyclotron element")); + std::string("The first item in the FieldDimensions list does not ") + +std::string("seem to be an Ring or a Cyclotron element")); } void ParallelCyclotronTracker::Tracker_MTS() { - IpplTimings::startTimer(IpplTimings::getTimer("MTS")); - IpplTimings::startTimer(IpplTimings::getTimer("MTS-Various")); + IpplTimings::startTimer(IpplTimings::getTimer("MTS")); + IpplTimings::startTimer(IpplTimings::getTimer("MTS-Various")); const double harm = getHarmonicNumber(); const double dt = itsBunch->getdT() * harm; @@ -4172,7 +4176,7 @@ void ParallelCyclotronTracker::Tracker_MTS() { if(OpalData::getInstance()->inRestartRun()) { restartStep0_m = itsBunch->getLocalTrackStep(); step_m = restartStep0_m; - if (numBunch_m > 1) itsBunch->resetPartBinID2(eta_m); + if (numBunch_m > 1) itsBunch->resetPartBinID2(eta_m); *gmsg << "* Restart at integration step " << restartStep0_m << endl; } if(OpalData::getInstance()->hasBunchAllocated() && Options::scan) { @@ -4267,7 +4271,7 @@ void ParallelCyclotronTracker::Tracker_MTS() { borisExternalFields(dt_inner); } - IpplTimings::startTimer(IpplTimings::getTimer("MTS-Various")); + IpplTimings::startTimer(IpplTimings::getTimer("MTS-Various")); // bunch injection // TODO: Where is correct location for this piece of code? Beginning/end of step? Before field solve? if(numBunch_m > 1) { @@ -4521,8 +4525,8 @@ void ParallelCyclotronTracker::Tracker_MTS() { for(size_t ii = 0; ii < itsBunch->getLocalNum(); ++ii) { if(itsBunch->ID[ii] == 0) { double FinalMomentum2 = pow(itsBunch->P[ii](0), 2.0) + - pow(itsBunch->P[ii](1), 2.0) + - pow(itsBunch->P[ii](2), 2.0); + pow(itsBunch->P[ii](1), 2.0) + + pow(itsBunch->P[ii](2), 2.0); double FinalEnergy = (sqrt(1.0 + FinalMomentum2) - 1.0) * itsBunch->getM() * 1.0e-6; *gmsg << "* Final energy of reference particle = " << FinalEnergy << " [MeV]" << endl; *gmsg << "* Total phase space dump number including the initial distribution) = " << lastDumpedStep_m + 1 << endl; @@ -4579,8 +4583,8 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t> & particle particleVectors -= translationToGlobal; Tenzor<double, 3> const rotation( cos(phi), sin(phi), 0, - -sin(phi), cos(phi), 0, - 0, 0, 1); // clockwise rotation + -sin(phi), cos(phi), 0, + 0, 0, 1); // clockwise rotation for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { @@ -4593,7 +4597,7 @@ void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t> & particle Tenzor<double, 3> const rotation(cos(phi), -sin(phi), 0, sin(phi), cos(phi), 0, - 0, 0, 1); // counter-clockwise rotation + 0, 0, 1); // counter-clockwise rotation for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { @@ -4680,7 +4684,7 @@ inline void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t> & p //if (fabs(psi) > tolerance) rotateAroundX(particleVectors, -psi); - // Rotation to align P_mean with x-axis + // Rotation to align P_mean with x-axis rotateAroundZ(particleVectors, -phi); // Translation from local to global @@ -4699,7 +4703,7 @@ inline void ParallelCyclotronTracker::localToGlobal(Vector_t & myVector, //if (fabs(psi) > tolerance) rotateAroundX(myVector, -psi); - // Rotation to align P_mean with x-axis + // Rotation to align P_mean with x-axis rotateAroundZ(myVector, -phi); // Translation from local to global @@ -4714,9 +4718,9 @@ inline void ParallelCyclotronTracker::rotateWithQuaternion(ParticleAttrib<Vector for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { particleVectors[i] = 2.0f * dot(quaternionVectorComponent, particleVectors[i]) * quaternionVectorComponent + - (quaternionScalarComponent * quaternionScalarComponent - - dot(quaternionVectorComponent, quaternionVectorComponent)) * particleVectors[i] + 2.0f * - quaternionScalarComponent * cross(quaternionVectorComponent, particleVectors[i]); + (quaternionScalarComponent * quaternionScalarComponent - + dot(quaternionVectorComponent, quaternionVectorComponent)) * particleVectors[i] + 2.0f * + quaternionScalarComponent * cross(quaternionVectorComponent, particleVectors[i]); } } @@ -4748,8 +4752,8 @@ inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t> & p // Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis Tenzor<double, 3> const rotation( cos(phi), sin(phi), 0, - -sin(phi), cos(phi), 0, - 0, 0, 1); + -sin(phi), cos(phi), 0, + 0, 0, 1); for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { @@ -4761,8 +4765,8 @@ inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t & myVector, double // Clockwise rotation of single vector 'myVector' by 'phi' around Z axis Tenzor<double, 3> const rotation( cos(phi), sin(phi), 0, - -sin(phi), cos(phi), 0, - 0, 0, 1); + -sin(phi), cos(phi), 0, + 0, 0, 1); myVector = dot(rotation, myVector); } @@ -4841,7 +4845,7 @@ void ParallelCyclotronTracker::push(double h) { std::list<CavityCrossData> cavCrossDatas; for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) { - if(((*sindex)->first) == "CAVITY") { + if(((*sindex)->first) == ElementBase::RFCAVITY) { RFCavity * cav = static_cast<RFCavity *>(((*sindex)->second).second); CavityCrossData ccd = {cav, cav->getSinAzimuth(), cav->getCosAzimuth(), cav->getPerpenDistance() * 0.001}; cavCrossDatas.push_back(ccd); @@ -4850,7 +4854,7 @@ void ParallelCyclotronTracker::push(double h) { for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { Vector_t const oldR = itsBunch->R[i]; double const gamma = sqrt(1.0 + dot(itsBunch->P[i], itsBunch->P[i])); - double const c_gamma = c / gamma; + double const c_gamma = Physics::c / gamma; Vector_t const v = itsBunch->P[i] * c_gamma; itsBunch->R[i] += h * v; for(std::list<CavityCrossData>::const_iterator ccd = cavCrossDatas.begin(); ccd != cavCrossDatas.end(); ++ccd) { @@ -4882,7 +4886,7 @@ void ParallelCyclotronTracker::push(double h) { // Path length update double const dotP = dot(itsBunch->P[0], itsBunch->P[0]); double const gamma = sqrt(1.0 + dotP); - PathLength_m += h * sqrt(dotP) * c / gamma; + PathLength_m += h * sqrt(dotP) * Physics::c / gamma; itsBunch->setLPath(PathLength_m); IpplTimings::stopTimer(IntegrationTimer_m); @@ -4892,8 +4896,8 @@ void ParallelCyclotronTracker::kick(double h) { IpplTimings::startTimer(IntegrationTimer_m); double const q = itsBunch->Q[0] / q_e; // For now all particles have the same charge double const M = itsBunch->M[0] * 1.0e9; // For now all particles have the same rest energy - double const h12Halfqc_M = 0.5 * h * q * c / M; - double const h12Halfqcc_M = h12Halfqc_M * c; + double const h12Halfqc_M = 0.5 * h * q * Physics::c / M; + double const h12Halfqcc_M = h12Halfqc_M * Physics::c; for(unsigned int i = 0; i < itsBunch->getLocalNum(); ++i) { // Half step E @@ -4945,7 +4949,7 @@ void ParallelCyclotronTracker::borisExternalFields(double h) { push(0.5 * h); IpplTimings::stopTimer(IpplTimings::getTimer("MTS-PushAndRFKick")); - IpplTimings::startTimer(IpplTimings::getTimer("MTS-PluginElements")); + IpplTimings::startTimer(IpplTimings::getTimer("MTS-PluginElements")); // apply the plugin elements: probe, collimator, stripper, septum itsBunch->R *= Vector_t(1000.0); // applyPluginElements expects [R] = mm applyPluginElements(h * 1e9); // expects [dt] = ns @@ -4959,80 +4963,80 @@ void ParallelCyclotronTracker::borisExternalFields(double h) { void ParallelCyclotronTracker::applyPluginElements(const double dt) { - for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); sindex++) { - if(((*sindex)->first) == "SEPTUM") { - (static_cast<Septum *>(((*sindex)->second).second))->checkSeptum(*itsBunch); - } + for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); sindex++) { + if(((*sindex)->first) == ElementBase::SEPTUM) { + (static_cast<Septum *>(((*sindex)->second).second))->checkSeptum(*itsBunch); + } - if(((*sindex)->first) == "PROBE") { - (static_cast<Probe *>(((*sindex)->second).second))->checkProbe(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); - } + if(((*sindex)->first) == ElementBase::PROBE) { + (static_cast<Probe *>(((*sindex)->second).second))->checkProbe(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); + } - if(((*sindex)->first) == "STRIPPER") { - bool flag_stripper = (static_cast<Stripper *>(((*sindex)->second).second)) - -> checkStripper(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); - if(flag_stripper) { - itsBunch->boundp(); - *gmsg << "* Total number of particles after stripping = " << itsBunch->getTotalNum() << endl; - } - } + if(((*sindex)->first) == ElementBase::STRIPPER) { + bool flag_stripper = (static_cast<Stripper *>(((*sindex)->second).second)) + -> checkStripper(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); + if(flag_stripper) { + itsBunch->boundp(); + *gmsg << "* Total number of particles after stripping = " << itsBunch->getTotalNum() << endl; + } + } - if(((*sindex)->first) == "CCOLLIMATOR") { - Collimator * collim; - collim = static_cast<Collimator *>(((*sindex)->second).second); - collim->checkCollimator(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); + if(((*sindex)->first) == ElementBase::COLLIMATOR) { + Collimator * collim; + collim = static_cast<Collimator *>(((*sindex)->second).second); + collim->checkCollimator(*itsBunch, turnnumber_m, itsBunch->getT() * 1e9, dt); + } } - } } bool ParallelCyclotronTracker::deleteParticle(){ - // Update immediately if any particles are lost during this step + // Update immediately if any particles are lost during this step - bool flagNeedUpdate = (min(itsBunch->Bin) < 0); - reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign()); + bool flagNeedUpdate = (min(itsBunch->Bin) < 0); + reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign()); - size_t lostParticleNum = 0; + size_t lostParticleNum = 0; - if(flagNeedUpdate) { - Vector_t const meanR = calcMeanR(); - Vector_t const meanP = calcMeanP(); + if(flagNeedUpdate) { + Vector_t const meanR = calcMeanR(); + Vector_t const meanP = calcMeanP(); - // Bunch (local) azimuth at meanR w.r.t. y-axis - double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; + // Bunch (local) azimuth at meanR w.r.t. y-axis + double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; - // Bunch (local) elevation at meanR w.r.t. xy plane - double const psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); + // Bunch (local) elevation at meanR w.r.t. xy plane + double const psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); - // For statistics data, the bunch is transformed into a local coordinate system - // with meanP in direction of y-axis -DW - globalToLocal(itsBunch->R, phi, psi, meanR); - globalToLocal(itsBunch->P, phi, psi, Vector_t(0.0)); // P should be rotated, but not shifted + // For statistics data, the bunch is transformed into a local coordinate system + // with meanP in direction of y-axis -DW + globalToLocal(itsBunch->R, phi, psi, meanR); + globalToLocal(itsBunch->P, phi, psi, Vector_t(0.0)); // P should be rotated, but not shifted - itsBunch->R *= Vector_t(0.001); // mm --> m + itsBunch->R *= Vector_t(0.001); // mm --> m - for(unsigned int i = 0; i < itsBunch->getLocalNum(); i++) { - if(itsBunch->Bin[i] < 0) { - lostParticleNum++; - itsBunch->destroy(1, i); - } - } + for(unsigned int i = 0; i < itsBunch->getLocalNum(); i++) { + if(itsBunch->Bin[i] < 0) { + lostParticleNum++; + itsBunch->destroy(1, i); + } + } - // Now destroy particles and update pertinent parameters in local frame - // Note that update() will be called within boundp() -DW - itsBunch->boundp(); - //itsBunch->update(); + // Now destroy particles and update pertinent parameters in local frame + // Note that update() will be called within boundp() -DW + itsBunch->boundp(); + //itsBunch->update(); - itsBunch->calcBeamParameters(); + itsBunch->calcBeamParameters(); - itsBunch->R *= Vector_t(1000.0); // m --> mm + itsBunch->R *= Vector_t(1000.0); // m --> mm - localToGlobal(itsBunch->R, phi, psi, meanR); - localToGlobal(itsBunch->P, phi, psi, Vector_t(0.0)); + localToGlobal(itsBunch->R, phi, psi, meanR); + localToGlobal(itsBunch->P, phi, psi, Vector_t(0.0)); - reduce(lostParticleNum, lostParticleNum, OpAddAssign()); - INFOMSG("Step " << step_m << ", " << lostParticleNum << " particles lost on stripper, collimator, septum, or out of cyclotron aperture" << endl); - } - return flagNeedUpdate; + reduce(lostParticleNum, lostParticleNum, OpAddAssign()); + INFOMSG("Step " << step_m << ", " << lostParticleNum << " particles lost on stripper, collimator, septum, or out of cyclotron aperture" << endl); + } + return flagNeedUpdate; } void ParallelCyclotronTracker::initTrackOrbitFile() { @@ -5061,7 +5065,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() { if(!OpalData::getInstance()->inRestartRun()) { // Start a new run (no restart) - double const initialReferenceTheta = referenceTheta * PIOVER180; + double const initialReferenceTheta = referenceTheta * Physics::deg2rad; PathLength_m = 0.0; // Force the initial phase space values of the particle with ID = 0 to zero, @@ -5119,7 +5123,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() { } } - // Else: Restart from the distribution in the h5 file + // Else: Restart from the distribution in the h5 file } else { // Do a local frame restart (we have already checked that the old h5 file was saved in local @@ -5146,7 +5150,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() { itsBunch->Bin[i] = 0; } - // Or do a global frame restart (no transformations necessary) + // Or do a global frame restart (no transformations necessary) } else { *gmsg << "* Restart in the global frame" << endl; @@ -5448,7 +5452,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { // ---------------- Re-calculate reference values in format of input values ----------------- // // Position: referenceR = sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1)); - referenceTheta = theta / PIOVER180; + referenceTheta = theta / Physics::deg2rad; referenceZ = meanR(2); // Momentum in Theta-hat, R-hat, Z-hat coordinates at position meanR: @@ -5473,23 +5477,23 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { // -------------- If flag DumpLocalFrame is not set, dump bunch in global frame ------------- // if (Options::psDumpFreq < std::numeric_limits<int>::max() ){ if(!(Options::psDumpLocalFrame) && (Options::psDumpFreq < std::numeric_limits<int>::max())) { - + FDext_m[0] = extB_m * 0.1; // kgauss --> T FDext_m[1] = extE_m; lastDumpedStep_m = itsDataSink->writePhaseSpace_cycl(*itsBunch, // Global and in m - FDext_m, E, - referencePr, - referencePt, - referencePz, - referenceR, - referenceTheta, - referenceZ, - phi / PIOVER180, // P_mean azimuth - // at ref. R/Th/Z - psi / PIOVER180, // P_mean elevation - // at ref. R/Th/Z - false); // Flag localFrame + FDext_m, E, + referencePr, + referencePt, + referencePz, + referenceR, + referenceTheta, + referenceZ, + phi / Physics::deg2rad, // P_mean azimuth + // at ref. R/Th/Z + psi / Physics::deg2rad, // P_mean elevation + // at ref. R/Th/Z + false); // Flag localFrame // Tell user in which mode we are dumping *gmsg << endl << "* Phase space dump " << lastDumpedStep_m @@ -5502,7 +5506,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { // The bunch is transformed into a local coordinate system with meanP in direction of y-axis globalToLocal(itsBunch->R, phi, psi, meanR * Vector_t(0.001)); globalToLocal(itsBunch->P, phi, psi); // P should only be rotated - + globalToLocal(extB_m, phi, psi); globalToLocal(extE_m, phi, psi); @@ -5510,18 +5514,18 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { FDext_m[1] = extE_m; lastDumpedStep_m = itsDataSink->writePhaseSpace_cycl(*itsBunch, // Local and in m - FDext_m, E, - referencePr, - referencePt, - referencePz, - referenceR, - referenceTheta, - referenceZ, - phi / PIOVER180, // P_mean azimuth - // at ref. R/Th/Z - psi / PIOVER180, // P_mean elevation - // at ref. R/Th/Z - true); // Flag localFrame + FDext_m, E, + referencePr, + referencePt, + referencePz, + referenceR, + referenceTheta, + referenceZ, + phi / Physics::deg2rad, // P_mean azimuth + // at ref. R/Th/Z + psi / Physics::deg2rad, // P_mean elevation + // at ref. R/Th/Z + true); // Flag localFrame localToGlobal(itsBunch->R, phi, psi, meanR * Vector_t(0.001)); localToGlobal(itsBunch->P, phi, psi); @@ -5542,8 +5546,8 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { *gmsg << "* Bunch position: R = " << referenceR << " mm" << ", Theta = " << referenceTheta << " Deg" << ", Z = " << referenceZ << " mm" << endl; - *gmsg << "* Local Azimuth = " << phi / PIOVER180 << " Deg" - << ", Local Elevation = " << psi / PIOVER180 << " Deg" << endl; + *gmsg << "* Local Azimuth = " << phi / Physics::deg2rad << " Deg" + << ", Local Elevation = " << psi / Physics::deg2rad << " Deg" << endl; IpplTimings::stopTimer(DumpTimer_m); } @@ -5583,4 +5587,4 @@ void ParallelCyclotronTracker::evaluateSpaceChargeField() { localToGlobal(itsBunch->Bf, phi); localToGlobal(itsBunch->R, phi, meanR); } -} +} \ No newline at end of file diff --git a/src/Algorithms/ParallelCyclotronTracker.h b/src/Algorithms/ParallelCyclotronTracker.h index b6ba03227a9e75b4833bc9677143af2364dceb54..ba25455902cbe7377d7b7321f066ebfcb9994840 100644 --- a/src/Algorithms/ParallelCyclotronTracker.h +++ b/src/Algorithms/ParallelCyclotronTracker.h @@ -21,6 +21,7 @@ #include "Algorithms/Tracker.h" #include "Structure/DataSink.h" +#include "AbsBeamline/ElementBase.h" #include <vector> class BMultipoleField; @@ -47,7 +48,7 @@ class ParallelCyclotronTracker: public Tracker { public: typedef std::pair<double[8], Component *> element_pair; - typedef std::pair<std::string, element_pair> type_pair; + typedef std::pair<ElementBase::ElementType, element_pair> type_pair; typedef std::list<type_pair *> beamline_list; /// Constructor. // The beam line to be tracked is "bl". @@ -268,7 +269,7 @@ private: #endif // GENERICTRACKER /* - Local Variables both used by the integration methods + Local Variables both used by the integration methods */ Vector_t rold_m, pold_m, rnew_m, ptmp_m; @@ -339,7 +340,7 @@ private: void applyExitFringe(double edge, double curve, const BMultipoleField &field, double scale); - void buildupFieldList(double BcParameter[], std::string ElementType, Component *elptr); + void buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr); bool derivate(double *y, const double &t, double *yp, const size_t &Pindex); @@ -503,4 +504,4 @@ inline double ParallelCyclotronTracker::calculateAngle2(double x, double y) { return atan2(y,x); } -#endif // OPAL_ParallelCyclotronTracker_HH +#endif // OPAL_ParallelCyclotronTracker_HH \ No newline at end of file