SamplerAnalysis.hh

/* 
Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway, 
University of London 2001 - 2023.
 
This file is part of BDSIM.
 
BDSIM is free software: you can redistribute it and/or modify 
it under the terms of the GNU General Public License as published 
by the Free Software Foundation version 3 of the License.
 
BDSIM is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with BDSIM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SAMPLERANALYSIS_H
#define SAMPLERANALYSIS_H
 
#include "BDSOutputROOTEventSampler.hh"
 
class SamplerAnalysis
{
public:
  SamplerAnalysis();
 
  // Note compile time float / double templating.
#ifndef __ROOTDOUBLE__
  SamplerAnalysis(BDSOutputROOTEventSampler<float>* samplerIn,
                  bool debugIn = false);
#else 
  SamplerAnalysis(BDSOutputROOTEventSampler<double>* samplerIn,
                  bool debugIn = false);
#endif
  void CommonCtor();
  virtual ~SamplerAnalysis();
 
  void Initialise();
 
  void Process(bool firstTime = false);
 
  std::vector<double>  Terminate(std::vector<double> emittance,
                                 bool useEmittanceFromFirstSampler = true);
 
  std::vector<std::vector<double> > GetOpticalFunctions() {return optical;}
 
  static void UpdateMass(SamplerAnalysis* s);
 
  static void UpdateMass(const std::string& particleName);
 
#ifndef __ROOTDOUBLE__
  BDSOutputROOTEventSampler<float> *s;
#else 
  BDSOutputROOTEventSampler<double> *s;
  BDSOutputROOTEventSampler<double> *p;
#endif
protected:
 
  static double particleMass;
  
  // sums - initialised to zero as that's what they start at
  long long int npart;
  double S;
 
  //6d phase space coordinates for each event
  std::vector<double> coordinates;
 
  // initial values to use for mean subtraction.
  std::vector<double> offsets;
 
  typedef std::vector<std::vector<double>>                           twoDArray;
  typedef std::vector<std::vector<std::vector<double>>>              threeDArray; 
  typedef std::vector<std::vector<std::vector<std::vector<double>>>> fourDArray;
 
  fourDArray    powSums;
  fourDArray    cenMoms;
 
  fourDArray    powSumsFirst;
  fourDArray    cenMomsFirst;
  
  threeDArray   covMats;
  threeDArray   derivMats;
  twoDArray     optical;     
  twoDArray     varOptical;  
 
 
 
  double powSumToCentralMoment(fourDArray& powSum,
                               long long int npartIn,
                               int i,
                               int j,
                               int m,
                               int n);
 
  double centMomToCovariance(fourDArray&   centMoms,
                             long long int npartIn,
                             int k,
                             int i,
                             int j);
 
  double centMomToDerivative(fourDArray &centMoms, int k, int t, int i);
 
 
  //Prints out the beam correlation matrix at that sampler to the console
  //The beam correlation matrix gives good idication of the presence and magnitude of all couplings.
  void printBeamCorrelationMatrix(fourDArray&   centMoms);
 
private:
  bool debug;
};
 
#endif