SdMonteCarloPropagator.cc

Go to the documentation of this file.

#include <cmath>
#include <iostream>
#include <iomanip>

#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/Cholesky>

#include <time.h>

#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <fwk/VModule.h>
#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>
#include <utl/Is.h>

// #include <fwk/RandomEngineRegistry.h>
#include <fwk/LocalCoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
// #include <CLHEP/Random/RandGauss.h>
// #include <CLHEP/Random/Randomize.h>

#include <boost/random.hpp>
#include <boost/random/normal_distribution.hpp>

#include "SdMonteCarloPropagator.h"

using namespace sevt;
using namespace evt;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace std;
using namespace Eigen;

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "


namespace SdMonteCarloPropagatorKVI {

  // typedef CLHEP::RandGauss RandGauss;


  enum EPrimaryParameters {
    eBeta = 0,
    eGamma = 1,
    eXCore = 2,
    eYCore = 3,
    eCurvature = 4,
    eShowerSize = 5,
    eUv = 6,
    eVv = 7,
    eCoreTime = 8, // Last!
  };


  const int kNumParameters = eCoreTime + 1;


  typedef Matrix<double, kNumParameters, kNumParameters> MatrixmNxNd;
  typedef Matrix<double, kNumParameters, 1> VectorNd;
  typedef Transpositions<kNumParameters, kNumParameters, int> TranspositionmNxNd;


  VModule::ResultFlag
  SdMonteCarloPropagator::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdMonteCarloPropagator");

    topB.GetChild("infoLevel").GetData(fInfoLevel);
    topB.GetChild("randomSeed").GetData(fRandomSeed);
    // fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

    return eSuccess;
  }


  const int trillion = 1000000000; // Workaround for nanoseconds

  evt::Event* originalEvent = nullptr;


  VModule::ResultFlag
  SdMonteCarloPropagator::Run(evt::Event& event)
  {
    INFO("SdMonteCarloPropagator::Run()");

    // if the old event was stored then replace the event with the old event
    // this should happen on every second call of the SdMonteCarloPropagator
    if (originalEvent) {
      OUT(1) << "Restoring the original event." << endl;
      event = *originalEvent;
      delete originalEvent;
      originalEvent = nullptr;
      return eSuccess;
    }

    // Change the event ids in order to keep track of the original and the varied parameters
    int sEventId = event.GetSEvent().GetHeader().GetId();
    auto eventId = event.GetHeader().GetId();
    bool finishEarly = false;
    if (Is("_original").NotIn(eventId) && Is("_varied").NotIn(eventId)) {
      eventId += string("_original");
      finishEarly = true;
    } else {
      const auto originalStart = eventId.find("_original");
      unsigned int newNum = 0;
      if (originalStart == string::npos) {
        newNum = atoi(eventId.substr(eventId.length()-4).c_str()) + 1;
        sEventId += 1;
      } else {
        eventId = eventId.substr(0, originalStart) + "_varied0000";
        sEventId *= 100;
      }
      eventId = eventId.substr(0, eventId.length()-4);
      ostringstream numOut;
      numOut << eventId << setfill('0') << setw(4) << newNum;
      eventId = numOut.str();
    }
    event.GetHeader().SetId(eventId);
    event.GetSEvent().GetHeader().SetId(sEventId);
    ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();

    // Saving the old event
    originalEvent = new evt::Event(event);

    if (finishEarly) {
      // Do nothing but re-seed and return in order to keep the original parameters as well
      if (fRandomSeed) {
        OUT(2) << "Seeding the random generator with " << sEventId*fRandomSeed << endl;
        fRandgen.seed(sEventId*fRandomSeed);
      }
      return eSuccess;
    }

    // Not sure about how to use the random engine.
    // Can't get the wrapped registry stuff to work if I want to
    // have the option to seed it.
    // CLHEP::HepRandomEngine &randomEngine = fRandomEngine->GetEngine();
    // if (fRandomSeed) {
    //   randomEngine.setSeed(sEventId*fRandomSeed, 0);
    // }
    // Instead I will use the Mersenne twister from boost (for now)

    MatrixmNxNd covarianceMatrix; // Will contain the variances and covariances

    // Setting all elements to zero
    for (int i = 0; i < kNumParameters; ++i) {
      for (int j = 0; j < kNumParameters; ++j) {
        covarianceMatrix(i, j) = 0;
      }
    }

    VectorNd origParameters; // Will store the originally reconstructed primary shower parameters
    VectorNd variedParameters; // Will store the varied primary shower parameters
    VectorNd randomShiftInParameters; // Will contain random gaussian values

    const CoordinateSystemPtr gBaryCS = LocalCoordinateSystem::Create(recShower.GetBarycenter());

    if (!recShower.HasParameter(eShowerAxisX) and !recShower.HasParameter(eShowerAxisY)) {
      WARNING("There is no shower axis available! Sorry, I can't proceed with this event...");
      return eBreakLoop;
    }

    // Store original primary parameters into a vector
    origParameters(eBeta) = recShower.GetBeta();
    origParameters(eGamma) = recShower.GetGamma();
    const auto& core = recShower.GetCorePosition();
    origParameters(eXCore) = core.GetX(gBaryCS);
    origParameters(eYCore) = core.GetY(gBaryCS);
    origParameters(eCurvature) = recShower.GetCurvature();
    origParameters(eShowerSize) = recShower.GetShowerSize();
    origParameters(eUv) = recShower.GetParameter(eShowerAxisX);
    origParameters(eVv) = recShower.GetParameter(eShowerAxisY);
    // It is wrong that the time is stored in a non-floating point data-type. We'll try to work around this.
    origParameters(eCoreTime) = recShower.GetCoreTime().GetGPSNanoSecond();
    const unsigned coreTimeSeconds = recShower.GetCoreTime().GetGPSSecond();

    // Storing variances into the covariance matrix
    covarianceMatrix(eBeta, eBeta) = Sqr(recShower.GetBetaError());
    covarianceMatrix(eGamma, eGamma) = Sqr(recShower.GetGammaError());
    const auto& coreError = recShower.GetCoreError();
    covarianceMatrix(eXCore, eXCore) = Sqr(coreError.GetX(gBaryCS));
    covarianceMatrix(eYCore, eYCore) = Sqr(coreError.GetY(gBaryCS));
    covarianceMatrix(eCurvature, eCurvature) = Sqr(recShower.GetCurvatureError());
    covarianceMatrix(eShowerSize, eShowerSize) = Sqr(recShower.GetShowerSizeError());
    // The new ParameterStorage class does not use a coordinate system.
    covarianceMatrix(eUv, eUv) = recShower.GetParameterCovariance(eShowerAxisX, eShowerAxisX);
    covarianceMatrix(eVv, eVv) = recShower.GetParameterCovariance(eShowerAxisY, eShowerAxisY);
    // Workaround for core time, converting from int to float
    covarianceMatrix(eCoreTime, eCoreTime) = Sqr(recShower.GetCoreTimeError().GetInterval()/ns);

    // Storing covariances into the covariance matrix is done below
    // N.B. Not all covariances are implemented. E.g. the covariances between the core time and the core
    // position are not available). This is wrong. Future versions of offline should contain all relevant
    // covariances and this module should implement these as well.

    // The core error is stored in a datatype that is coordinate system free
    const double xCoreErr = recShower.GetCoreError().GetX(gBaryCS);
    const double yCoreErr = recShower.GetCoreError().GetY(gBaryCS);
    // The correlation of the error however is stored with an implicit coordinate system.
    // This is not in accordance with the 'philosophy' of OffLine
    // Carefull! Implicit assumption that the coordinate system is gBaryCS!!!!!
    const double xyCoreCov = recShower.GetCorrelationXY() * xCoreErr * yCoreErr;
    // The same comments hold here for the covariances between U and V, these are stored with the implicit gBaryCS in mind!
    // Carefull! Implicit assumption that the coordinate system is gBaryCS!!!!!
    const double uvCov = recShower.GetParameterCovariance(eShowerAxisX, eShowerAxisY);
    covarianceMatrix(eUv, eVv) = covarianceMatrix(eVv, eUv) = uvCov;
    covarianceMatrix(eXCore, eYCore) = covarianceMatrix(eYCore, eXCore) = xyCoreCov;

    OUT(2) << "Testing errorpropagator " << sEventId << endl;
    OUT(2) << "Original parameters:" << endl;
    OUT(2) << origParameters << endl;
    OUT(2) << "Covariance matrix:" << endl;
    OUT(2) << covarianceMatrix << endl;

    Eigen::LDLT<MatrixmNxNd> ldlt(covarianceMatrix);
    // Testing wether the covariance matrix is positive semidefinite
    // The isPositive() function is a bit misleading, isNonNegative() would have been better (or simply isPositiveSemiDefinite())
    if (!ldlt.isPositive()) {
      ERROR("Covariance matrix not positive (semi)definite! "
            "You are the weakest link, good bye! "
            "(by the way this should definitely not be happening. "
            "Beter make sure that the SD reconstruction is ok...)");
      return eBreakLoop;
    }

    // The variables below will have to satisfy
    // covarianceMatrix == P^T L D L^T P (see eigen3 documentation on LDLT decomposition)
    MatrixmNxNd cML(ldlt.matrixL());                // lower triangular matrix L
    TranspositionmNxNd cMP(ldlt.transpositionsP()); // permutation matrix P
    VectorNd cMD(ldlt.vectorD());                   // diagonal matrix D (represented as a vector)

    // Calculating the conventional square roots of the diagonals and putting them into a real matrix
    MatrixmNxNd cMDSqrt;
    cMDSqrt.setZero();
    for (int i = 0; i < kNumParameters; ++i)
      cMDSqrt(i, i) = sqrt(cMD(i));

    // Calculating the 'square root' of the covariance matrix
    MatrixmNxNd covarianceMatrixSqrt = cMP.transpose() * cML * cMDSqrt;

    // Testing wether covarianceMatrix == covarianceMatrixSqrt * covarianceMatrixSqrt.transpose()
    OUT(2) << "Testing wether the 'square root' of the covariance matrix is correctly computed... " << endl;
    if (!(covarianceMatrixSqrt * covarianceMatrixSqrt.transpose()).isApprox(covarianceMatrix, 1e-12) ) {
      ERROR("The test for the sanity of the covariance matrix failed! This should not be happening.");
      return eBreakLoop;
    }

    VectorNd randomGaussianValues; // Will contain random gaussian values with mu=0 and sigma=1
    for (int i = 0; i < kNumParameters; ++i)
      randomGaussianValues(i) = fRandomGenerator();

    randomShiftInParameters = covarianceMatrixSqrt * randomGaussianValues;
    variedParameters = origParameters + randomShiftInParameters;

    // Put the varied data back into the shower parameters...
    recShower.SetBeta(variedParameters(eBeta), recShower.GetBetaError());    // beta  Is stored into varied data structure
    recShower.SetGamma(variedParameters(eGamma), recShower.GetGammaError()); // gamma Is stored into varied data structure
    const double variedXCore = variedParameters(eXCore);                     // put XCore in separate variable
    const double variedYCore = variedParameters(eYCore);                     // put YCore in separate variable
    const double variedZcore = core.GetZ(gBaryCS);                           // store unaltered ZCore in separate variable
    recShower.SetCorePosition(utl::Point(variedXCore, variedYCore, variedZcore, gBaryCS));  // core position is stored into varied data structure
    recShower.SetCurvature(variedParameters(eCurvature), recShower.GetCurvatureError());    // curvature is stored into varied data structure
    recShower.SetShowerSize(variedParameters(eShowerSize), recShower.GetShowerSizeError()); // ShowerSize is stored into varied data structure
    recShower.SetParameter(eShowerAxisX, variedParameters(eShowerAxisX));                   // U is stored into varied data structure
    recShower.SetParameter(eShowerAxisY, variedParameters(eShowerAxisY));                   // V is stored into varied data structure
    // Originally the next lines were a workaround for faulty implemented nanoseconds in the core time
    // Now that the nanoseconds have become floats we have to see wether this teeds to be removed or changed.
    unsigned int variedSeconds = coreTimeSeconds;
    unsigned int variedNanoSeconds = trillion + variedParameters(eCoreTime);
    variedSeconds += variedNanoSeconds/trillion - 1;
    variedNanoSeconds %= trillion;
    const TimeStamp variedCoreTime(variedSeconds, variedNanoSeconds);
    recShower.SetCoreTime(variedCoreTime, recShower.GetCoreTimeError());                 // CoreTime is stored into varied data structure

    // Secondary shower parameters
    // Re-calculate the varied axis based on U and V and the gBaryCS
    const double axisU = variedParameters(eUv);
    const double axisV = variedParameters(eVv);
    const double axisW = sqrt(1 - Sqr(axisU) - Sqr(axisV));
    recShower.SetAxis(utl::Vector(axisU, axisV, axisW, gBaryCS));

    // Re-calculate the energy based on the energy calculation of the LDFFinderKG::EnergyConversion::GetEnergy()
    // The beautiful modularity of offline does not allow access to this function itself. Here we are forced to
    // copy the relevant features of this function into this line. If the error calculation in this GetEnergy() function
    // changes then this line below may give inconsistent results.
    const double energyError = recShower.GetEnergyError();
    recShower.SetEnergy(recShower.GetEnergy() * variedParameters(eShowerSize) / origParameters(eShowerSize), energyError);

    OUT(3) << "Varied parameters:\n" << variedParameters << endl;
    OUT(1) << "Finished varying parametrs for EventID " << eventId << ", SEventID " << sEventId << endl;

    return eSuccess;
  }

}