Issue #268 (and issue #269):

modified:   Algorithms/ParallelCyclotronTracker.cpp
modified:   Algorithms/ParallelCyclotronTracker.h
modified:   Classic/Algorithms/PartBunchBase.h
modified:   Classic/Algorithms/PartBunchBase.hpp
modified:   Structure/DataSink.cpp
modified:   Structure/DataSink.h
modified:   Track/TrackRun.cpp

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -132,6 +132,7 @@ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
     Tracker(beamline, reference, revBeam, revTrack),
     itsDataSink(nullptr),
     bgf_m(nullptr),
+    numBunch_m(1),
     lastDumpedStep_m(0),
     eta_m(0.01),
     myNode_m(Ippl::myNode()),
@@ -162,10 +163,12 @@ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
                                                    DataSink &ds,
                                                    const PartData &reference,
                                                    bool revBeam, bool revTrack,
-                                                   int maxSTEPS, int timeIntegrator):
+                                                   int maxSTEPS, int timeIntegrator,
+                                                   int numBunch):
     Tracker(beamline, bunch, reference, revBeam, revTrack),
     bgf_m(nullptr),
     maxSteps_m(maxSTEPS),
+    numBunch_m(numBunch),
     lastDumpedStep_m(0),
     eta_m(0.01),
     myNode_m(Ippl::myNode()),
@@ -2618,8 +2621,8 @@ void ParallelCyclotronTracker::singleParticleDump() {
 void ParallelCyclotronTracker::bunchDumpStatData(){
 
     // don't dump stat file in case of multi-bunch mode
-    if ( multiBunchMode_m != MB_MODE::NONE )
-        return;
+//     if ( multiBunchMode_m != MB_MODE::NONE )
+//         return;
 
     IpplTimings::startTimer(DumpTimer_m);
 

src/Algorithms/ParallelCyclotronTracker.h

@@ -87,7 +87,9 @@ public:
     //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
     //  If [b]revTrack[/b] is true, we track against the beam.
     explicit ParallelCyclotronTracker(const Beamline &bl, PartBunchBase<double, 3> *bunch, DataSink &ds,
-                                      const PartData &data, bool revBeam, bool revTrack, int maxSTEPS, int timeIntegrator);
+                                      const PartData &data, bool revBeam,
+                                      bool revTrack, int maxSTEPS,
+                                      int timeIntegrator, int numBunch);
 
     virtual ~ParallelCyclotronTracker();
 

src/Classic/Algorithms/PartBunchBase.h

@@ -232,6 +232,7 @@ public:
     Vector_t get_pmean_Distribution() const;
     Vector_t get_emit() const;
     Vector_t get_norm_emit() const;
+    Vector_t get_halo(std::size_t bin = 0) const;
     virtual Vector_t get_hr() const;
 
     double get_Dx() const;
@@ -590,6 +591,9 @@ protected:
 
     /// rms normalized emittance
     Vector_t eps_norm_m;
+
+    std::vector<Vector_t> halo_m;
+
     /// rms correlation
     Vector_t rprms_m;
 

src/Classic/Algorithms/PartBunchBase.hpp

@@ -43,6 +43,7 @@ PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb)
       pmean_m(0.0),
       eps_m(0.0),
       eps_norm_m(0.0),
+      halo_m(1, Vector_t(0.0, 0.0, 0.0)),
       rprms_m(0.0),
       Dx_m(0.0),
       Dy_m(0.0),
@@ -1183,6 +1184,14 @@ Vector_t PartBunchBase<T, Dim>::get_norm_emit() const {
 }
 
 
+template <class T, unsigned Dim>
+Vector_t PartBunchBase<T, Dim>::get_halo(std::size_t bin) const {
+    if ( bin >= halo_m.size() )
+        throw OpalException("PartBunchBase<T, Dim>::get_halo() ", "Out of range.");
+    return halo_m[bin];
+}
+
+
 template <class T, unsigned Dim>
 Vector_t PartBunchBase<T, Dim>::get_hr() const {
     return hr_m;
@@ -1914,8 +1923,18 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
      */
     std::vector<double> loc_moments(2 * Dim + Dim * ( 2 * Dim + 1 ));
 
+    /* Stored as
+     * bin 0: <x^2>, <x^4>, <y^2>, <y^4>, <z^2>, <z^4>
+     * bin 1: <x^2>, <x^4>, <y^2>, <y^4>, <z^2>, <z^4>
+     * ...
+     */
+    std::vector<double> loc_halo_moments_per_bin;
+
     long int totalNum = this->getTotalNum();
     if (OpalData::getInstance()->isInOPALCyclMode()) {
+
+        loc_halo_moments_per_bin.resize(2 * Dim * numBunch_m, 0.0);
+
         for(unsigned long k = 0; k < localNum; ++ k) {
             if (ID[k] == 0) {
                 part[1] = P[k](0);
@@ -1932,6 +1951,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
                         loc_moments[l++] -= part[i] * part[j];
                     }
                 }
+
                 --totalNum;
                 break;
             }
@@ -1955,6 +1975,16 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
                 loc_moments[l++] += part[i] * part[j];
             }
         }
+
+        if ( OpalData::getInstance()->isInOPALCyclMode() && ID[k] != 0) {
+            for (unsigned int i = 0; i < Dim; ++i) {
+                // <x^2>, <y^2>, <z^2>
+                loc_halo_moments_per_bin[2 * Dim * Bin[k] + 2 * i] += R[k](i) * R[k](i);
+                // <x^4>, <y^4>, <z^4>
+                loc_halo_moments_per_bin[2 * Dim * Bin[k] + 2 * i + 1] += R[k](i) * R[k](i) *
+                                                                          R[k](i) * R[k](i);
+            }
+        }
     }
 
     allreduce(loc_moments.data(), loc_moments.size(), std::plus<double>());
@@ -1971,6 +2001,26 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
         }
     }
 
+    if ( OpalData::getInstance()->isInOPALCyclMode() ) {
+        allreduce(loc_halo_moments_per_bin.data(),
+                  loc_halo_moments_per_bin.size(),
+                  std::plus<double>());
+
+        if ( (int)halo_m.size() < numBunch_m ) {
+            halo_m.resize(numBunch_m);
+        }
+
+        for (unsigned int i = 0; i < halo_m.size(); ++i) {
+            for (unsigned int j = 0; j < Dim; ++j) {
+                double w4 = loc_halo_moments_per_bin[2 * Dim * i + 2 * j + 1];
+                double w2 = loc_halo_moments_per_bin[2 * Dim * i + 2 * j];
+
+                if ( w2 > 0.0 )
+                    halo_m[i](j) =  w4 / (w2 * w2);
+            }
+        }
+    }
+
     return totalNum;
 }
 

src/Structure/DataSink.cpp

@@ -31,6 +31,7 @@
 extern Inform *gmsg;
 
 DataSink::DataSink() :
+    nMaxBunches_m(1),
     H5call_m(0),
     lossWrCounter_m(0),
     doHDF5_m(true),
@@ -39,6 +40,7 @@ DataSink::DataSink() :
 
 DataSink::DataSink(H5PartWrapper *h5wrapper, int restartStep):
     mode_m(std::ios::out),
+    nMaxBunches_m(1),
     H5call_m(0),
     lossWrCounter_m(0),
     h5wrapper_m(h5wrapper)
@@ -121,6 +123,7 @@ DataSink::DataSink(H5PartWrapper *h5wrapper, int restartStep):
 
 DataSink::DataSink(H5PartWrapper *h5wrapper):
     mode_m(std::ios::out),
+    nMaxBunches_m(1),
     H5call_m(0),
     lossWrCounter_m(0),
     h5wrapper_m(h5wrapper)
@@ -159,6 +162,10 @@ void DataSink::storeCavityInformation() {
     h5wrapper_m->storeCavityInformation();
 }
 
+void DataSink::setMaxNumBunches(int nBunches) {
+    nMaxBunches_m = nBunches;
+}
+
 void DataSink::writePhaseSpace(PartBunchBase<double, 3> *beam, Vector_t FDext[]) {
 
     if (!doHDF5_m) return;
@@ -434,6 +441,23 @@ void DataSink::doWriteStatData(PartBunchBase<double, 3> *beam, Vector_t FDext[],
             os_statData << beam->P[0](2) << std::setw(pwi) << "\t";         // 48 P0_z
         }
 
+        if (OpalData::getInstance()->isInOPALCyclMode()) {
+
+            int nBunch = beam->getNumBunch();
+
+            for (int bin = 0; bin < nBunch; ++bin) {
+                Vector_t halo = beam->get_halo(bin);
+                for (int i = 0; i < 3; ++i)
+                    os_statData << halo(i) << std::setw(pwi) << "\t";
+            }
+
+            for ( int bin = nBunch; bin < nMaxBunches_m; ++bin ) {
+                Vector_t halo = Vector_t(0.0, 0.0, 0.0);
+                for (int i = 0; i < 3; ++i)
+                    os_statData << halo(i) << std::setw(pwi) << "\t";
+            }
+        }
+
 
         for(size_t i = 0; i < losses.size(); ++ i) {
             os_statData << losses[i].second << std::setw(pwi) << "\t";
@@ -920,6 +944,29 @@ void DataSink::writeSDDSHeader(std::ofstream &outputFile,
         columnStart = 49;
     }
 
+    if (OpalData::getInstance()->isInOPALCyclMode()) {
+        char dir[] = { 'x', 'y', 'z' };
+
+        for (int bin = 0; bin < nMaxBunches_m; ++bin) {
+            for (int i = 0; i < 3; ++i) {
+                std::stringstream ss;
+                ss << "&column\n" << indent << "name=halo_" << dir[i];
+                
+                if ( nMaxBunches_m > 1 ) {
+                    ss << " (energy bin " << bin << "),\n";
+                } else {
+                    ss << ",\n";
+                }
+                ss << indent << "type=double,\n"
+                   << indent << "units=1,\n"
+                   << indent << "description=\"" << columnStart++ << " Halo in "
+                   << dir[i] << "\"\n"
+                   << "&end\n";
+                outputFile << ss.str();
+            }
+        }
+    }
+
     for (size_t i = 0; i < losses.size(); ++ i) {
         outputFile << "&column\n"
                    << indent << "name=" << losses[i].first << ",\n"

src/Structure/DataSink.h

@@ -71,6 +71,8 @@ public:
      */
     void storeCavityInformation();
 
+    void setMaxNumBunches(int nBunches);
+
     /** \brief Write statistical data.
      *
      * Writes statistical beam data to proper output file. This is information such as RMS beam parameters
@@ -258,6 +260,8 @@ private:
      */
     std::ios_base::openmode mode_m;
 
+    int nMaxBunches_m;
+
     /** \brief First write to the H5 surface loss file.
      *
      * If true, file name will be assigned and file will be prepared to write.

src/Track/TrackRun.cpp

@@ -757,13 +757,13 @@ void TrackRun::setupCyclotronTracker(){
     *gmsg << "* Phase space dump frequency = " << Options::psDumpFreq << endl;
     *gmsg << "* Statistics dump frequency  = " << Options::statDumpFreq << " w.r.t. the time step." << endl;
     *gmsg << "* ********************************************************************************** " << endl;
+    
+    ds->setMaxNumBunches(specifiedNumBunch);
 
     itsTracker = new ParallelCyclotronTracker(*Track::block->use->fetchLine(),
                                               Track::block->bunch, *ds, Track::block->reference,
                                               false, false, Track::block->localTimeSteps.front(),
-                                              Track::block->timeIntegrator);
-
-    itsTracker->setNumBunch(specifiedNumBunch);
+                                              Track::block->timeIntegrator, specifiedNumBunch);
 
     if(opal->inRestartRun()) {