RdGeoCeLDFFitter.h

Go to the documentation of this file.

#ifndef _RdGeoCeLDFFitter_RdGeoCeLDFFitter_h_
#define _RdGeoCeLDFFitter_RdGeoCeLDFFitter_h_

#include <fwk/VModule.h>

#include <utl/CoordinateSystemPtr.h>

#include <revt/REvent.h>
#include <evt/ShowerRRecData.h>

#include <rdet/RDetector.h>

#include <string>
#include <vector>
#include <algorithm>
#include <cmath>

#include <TF2.h>
#include <TRandom3.h>


namespace evt {
  class Event;
}

namespace revt {
  class Station;
}


namespace RdGeoCeLDFFitter {

  class CorRand {

    TRandom3* rnd = nullptr;
    Double_t RT[2][2] = { { 0 } };
    Double_t myMu[2] = { 0 };
    TF2* gausFkt = nullptr;

  public:
    CorRand(const Double_t mu[2], const Double_t sig[2], const Double_t rho, unsigned int seed)
    {
      myMu[0] = mu[0];
      myMu[1] = mu[1];
      rnd = new TRandom3(seed);
      RT[0][0] = sig[0];
      RT[1][0] = 0;
      RT[0][1] = rho * sig[1];
      RT[1][1] = sig[1] * sqrt(1 - rho * rho);
      gausFkt = new TF2(
        "2Dim Gauss",
        "exp( -0.5 / (1.-[4]**2) * ( ((x-[0])/[1])**2 + ((y-[2])/[3])**2 - 2.*[4]*((x-[0])/[1])*((y-[2])/[3])) ) / ( 2.*TMath::Pi()*[1]*[3]*sqrt(1.-[4]**2) )",
        mu[0] - 10 * sig[0], mu[0] + 10 * sig[0], mu[1] - 10 * sig[1], mu[1] + 10 * sig[1]
      );
      gausFkt->SetParameters(mu[0], sig[0], mu[1], sig[1], rho);
    }

    void
    GetRnd(double* const a, double* const b)
    {
      const Double_t normRnd[2] = { rnd->Gaus(0, 1), rnd->Gaus(0, 1) };
      const Double_t corRnd[2] = {
        myMu[0] + RT[0][0] * normRnd[0],  // RT[1][0]=0!
        myMu[1] + RT[0][1] * normRnd[0] + RT[1][1] * normRnd[1]
      };
      *a = corRnd[0];
      *b = corRnd[1];
    }

    TF2* GetFkt() { return gausFkt; }

    ~CorRand()
    {
      delete rnd;
      delete gausFkt;
    }

  };

  class RdGeoCeLDFFitter : public fwk::VModule {

  public:
    fwk::VModule::ResultFlag Init() override;
    fwk::VModule::ResultFlag Run(evt::Event& event) override;
    fwk::VModule::ResultFlag Finish() override;

  private:
    bool fFixCore = false;
    bool fUserSubThresholdStations = true;
    int fMinNuberOfStationsForCoreFit = 5;
    unsigned int fRandomSeed = 0;

    REGISTER_MODULE("RdGeoCeLDFFitter", RdGeoCeLDFFitter);
  };

}

#endif