RdLDFMultiFitter.cc

Go to the documentation of this file.

#include "RdLDFMultiFitter.h"


/* ################################################################################################

      Removed support for this module from framework and adst with revision 34470 (01.11.21)

#################################################################################################*/

#include <cmath>
#include <sstream>

#include <fwk/CentralConfig.h>

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/BasicVector.h>
#include <utl/GeometryUtilities.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>

#include <utl/AxialVector.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerRRecData.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <det/VManager.h>

#include "TGraphErrors.h"
#include "TF1.h"
#include "TVector3.h"
#include "TMath.h"
#include "TColor.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"

#include <cstddef>
#include <string>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <stdlib.h>
#include <sys/stat.h>

using namespace revt;
using namespace fwk;
using namespace utl;
using namespace det;
using std::string;
using std::vector;

namespace RdLDFMultiFitter {

  RdLDFMultiFitter::RdLDFMultiFitter() :
      fInfoLevel(eNone),
      fParA200(0),
      fParA201(0),
      fParA210(0),
      fParA211(0),
      fEnergyEstimatorDist(0),
      fEnergyEstimatorPower(0),
      fEnergyEstimatorCoeff(0)
  {
  }

  RdLDFMultiFitter::~RdLDFMultiFitter()
  {
  }

  VModule::ResultFlag RdLDFMultiFitter::Init()
  {
    INFO("RdLDFMultiFitter::Init()");

    utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdLDFMultiFitter");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("parA200").GetData(fParA200);
    topBranch.GetChild("parA201").GetData(fParA201);
    topBranch.GetChild("parA210").GetData(fParA210);
    topBranch.GetChild("parA211").GetData(fParA211);
    topBranch.GetChild("energyEstimatorDist").GetData(fEnergyEstimatorDist);
    topBranch.GetChild("energyEstimatorPower").GetData(fEnergyEstimatorPower);
    topBranch.GetChild("energyEstimatorCoeff").GetData(fEnergyEstimatorCoeff);

    Branch models = topBranch.GetChild("Models");

    for (Branch curmodel = models.GetFirstChild(); curmodel; curmodel = curmodel.GetNextSibling()) {
       if (curmodel.GetName() == "Model") {
         Model modcontainer;
         modcontainer.name = VManager::FindComponent<string>("name", curmodel.GetAttributes());
         for(Branch curp = curmodel.GetFirstChild(); curp; curp = curp.GetNextSibling()) {
           if(curp.GetName() == "formula") {
             curp.GetData(modcontainer.formula);
           }
           if(curp.GetName() == "color") {
             curp.GetData(modcontainer.color);
           }
           if(curp.GetName() == "parameter") {
             parameter p;
             p.id = VManager::FindComponent<int>("id", curp.GetAttributes());
             p.name = VManager::FindComponent<string>("name", curp.GetAttributes());
             p.minval = VManager::FindComponent<double>("minval", curp.GetAttributes());
             p.maxval = VManager::FindComponent<double>("maxval", curp.GetAttributes());
             p.zenithCorrection = VManager::FindComponent<double>("zenithCorrection", curp.GetAttributes());
             curp.GetData(p.startvalues);
             modcontainer.parameters.push_back(p);
           }
         }
         fModels.push_back(modcontainer);
       }
    }


    return eSuccess;
  }

  VModule::ResultFlag RdLDFMultiFitter::Run(evt::Event& event)
  {

    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eSuccess;
    }

    REvent& rEvent = event.GetREvent();
    evt::ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    // Warning / TODO: need functionality for non-radio cores
    if(!showerrrec.HasCorePosition() || !showerrrec.HasAxis()) {
      WARNING("No core or axis. LDF Fit can not be performed.");
      return eSuccess;
    }

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    const utl::Vector showerAxis = showerrrec.GetAxis();
    const utl::Point coordinateOrigin = showerrrec.GetCoordinateOrigin();
    const utl::Point core = showerrrec.GetCorePosition();

    int Nant = showerrrec.GetParameter(revt::eNumberOfStationsWithPulseFound);
    int curAnt = 0;
    double distance[Nant], fieldStr[Nant], distanceEr[Nant], fieldStrEr[Nant];
    double distancetoshoweraxis = 0.0;

    // Warning / TODO: #warning Iterator has to be changed to stations with pulse
    for(REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin(); sIt != rEvent.CandidateStationsEnd(); ++sIt) {

      Station& station = *sIt;

      const StationRecData& sRec = station.GetRecData();
      const utl::Point sPos = rDetector.GetStation(station).GetPosition();

      if(!sRec.GetPulseFound())
        continue;

      distancetoshoweraxis = utl::RadioGeometryUtilities::GetDistanceToAxis(showerAxis, core, sPos);
      distance[curAnt] = distancetoshoweraxis;
      distanceEr[curAnt] = 0;
      // Warning / TODO: #warning need to define error of distance, currently put to zero (case of Tunka-Rex)
      fieldStr[curAnt] = sRec.GetParameter(revt::eSignal);
      fieldStrEr[curAnt] = sRec.GetParameterError(revt::eSignal);

      ++curAnt;
    }

    TGraphErrors *LateralDistribution = new TGraphErrors (Nant, distance,fieldStr,distanceEr,fieldStrEr);

    TF1 *energy_estimator;
    energy_estimator = new TF1("exp_formula", "[0]*exp(-[1]*(x-[2]))", 0, 5000);
    // Warning / TODO: #warning domain is fixed to (0, 5000 m)
    energy_estimator->SetParLimits(0,0.,1.);
    energy_estimator->SetParameter(0,fieldStr[0]);
    energy_estimator->FixParameter(2,fEnergyEstimatorDist);
    LateralDistribution->Fit(energy_estimator, "RS");

    double estimated_energy = fEnergyEstimatorCoeff * pow(energy_estimator->GetParameter(0), fEnergyEstimatorPower);

    vector<TF1> ldfResult;
    vector<TFitResult> ldfFitResult;

    for (vector<Model>::iterator m_it = fModels.begin() ; m_it != fModels.end(); ++m_it) {

      TF1 *fitfuncExp;

      fitfuncExp = new TF1(m_it->name.c_str(),m_it->formula.c_str(), 0, 5000);
      // Warning / TODO: #warning: domain is fixed to (0, 5000 m)
      fitfuncExp->SetLineColor(TColor::GetColor(m_it->color.c_str()));

      for(vector<parameter>::iterator p_it = m_it->parameters.begin(); p_it != m_it->parameters.end(); ++p_it) {

        fitfuncExp->SetParName(p_it->id,p_it->name.c_str());
        fitfuncExp->SetParLimits(p_it->id,p_it->minval,p_it->maxval);

        if(p_it->zenithCorrection) {
          fitfuncExp->FixParameter(p_it->id,ZenithCorrection(showerrrec.GetZenith(), estimated_energy));
      }

    }

      TFitResultPtr fitresult = LateralDistribution->Fit(fitfuncExp, "RS");
      ldfResult.push_back(*fitfuncExp);
      ldfFitResult.push_back(*fitresult);

      fitfuncExp->Delete();
    }

    showerrrec.SetRdLDF(ldfResult);
    showerrrec.SetRdLDFFitResult(ldfFitResult);

    return eSuccess;
  }

  double RdLDFMultiFitter::ZenithCorrection(double zenith, double energy) {
    return p0(energy) + p1(energy)*cos(zenith);
  }

  double RdLDFMultiFitter::p0(double energy) {
    return fParA200 + fParA201*energy;
  }

  double RdLDFMultiFitter::p1(double energy) {
    return fParA210 + fParA211*energy;
  }

  VModule::ResultFlag RdLDFMultiFitter::Finish()
  {
    INFO("RdLDFMultiFitter::Finish()");

    return eSuccess;
  }

}