RootCFMatrixOutput.cc

Go to the documentation of this file.

#include "RootCFMatrixOutput.h"
#include "CherenkovFluorescenceMatrix.h"
#include "lowerTriangularMatrix.h"
#include "matrixInversionExceptions.h"
#include "diagonalMatrix.h"
#include "lowerTriangularMatrix.h"
#include "upperTriangularMatrix.h"

#include <TFile.h>
#include <TGraph.h>
#include <TH2D.h>

#include <boost/format.hpp>
#include <evt/Event.h>
#include <evt/ShowerFRecData.h>
#include <evt/GaisserHillas4Parameter.h>
#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/EyeHeader.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/AugerUnits.h>
#include <utl/UTMPoint.h>
#include <utl/Reader.h>
#include <det/Detector.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <fwk/CentralConfig.h>

#include <sstream>
#include <iostream>
#include <string>
#include <cmath>

using namespace FdProfileReconstructorKG;
using namespace oBLAS;
using namespace fevt;
using namespace sevt;
using namespace evt;
using namespace utl;
using namespace atm;
using namespace det;
using namespace fwk;
using namespace boost;
using namespace std;


RootCFMatrixOutput::RootCFMatrixOutput()
{
  Init();
}


void
RootCFMatrixOutput::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  const Branch topB = cc->GetTopBranch("FdProfileReconstructor");
  const Branch rootB = topB.GetChild("rootOutput");
  rootB.GetChild("verbosity").GetData(fVerbosity);
  string outputFileName;
  rootB.GetChild("outputFileName").GetData(outputFileName);
  fOutputFile = new TFile(outputFileName.c_str(),"RECREATE");
  if (fVerbosity >- 1) {
    std::ostringstream info;
    info << "\n\t saving matrix visualization to "
         << outputFileName << endl;
    INFO(info);
  }
}


void
RootCFMatrixOutput::Write(const evt::Event& event,
                          const fevt::Eye& eye,
                          CherenkovFluorescenceMatrix chfl)
{
  if (fVerbosity > -1)
    cout << "\n writing ROOT output" << endl;

  if (!event.HasSEvent())
    return;
  const SEvent& sEvent = event.GetSEvent();

  const int eventId = sEvent.GetHeader().GetId();
  const int eyeId = eye.GetId();

  const string eventDirName = (format("Event_%04d_eye%1d") % eventId % eyeId ).str();

  fOutputFile->mkdir(eventDirName.c_str(), " ");
  fOutputFile->cd(eventDirName.c_str());

  WriteLateral(eye, chfl);
}


void
RootCFMatrixOutput::WriteLateral(const fevt::Eye& eye,
                                 CherenkovFluorescenceMatrix& ChFl)
{
  namespace ublas = boost::numeric::ublas;

  const evt::ShowerFRecData& showerFRecData = eye.GetRecData().GetFRecShower();
  const evt::GaisserHillas4Parameter* const gh4 =
    dynamic_cast<const evt::GaisserHillas4Parameter*>(&showerFRecData.GetGHParameters());

  const double Xmax = gh4->GetXMax();
  const double dEdXmax = gh4->GetNMax();
  const double lambda = gh4->GetShapeParameter(gh::eLambda);
  const double X0 = gh4->GetShapeParameter(gh::eX0);

  const oBLAS::lowerTriangularMatrix& chFlMatrix = *ChFl.GetCherenkovFluorescenceMatrix();
  unsigned int n = chFlMatrix.size1();

  const std::deque<double>& depth = *ChFl.GetDepth();

  ublas::vector<double> dEdX(n);
  for (unsigned int i = 0; i < n; ++i)
    dEdX(i) = GaisserHillas(depth[i], X0, Xmax, dEdXmax, lambda);

  // all light
  double zeta = 90.*degree;
  ChFl.SetZeta(zeta);
  ublas::vector<double> allFluo = ublas::prod(*ChFl.GetFluorescenceMatrix(), dEdX);
  ublas::vector<double> allChMScat =ublas::prod(*ChFl.GetMieScatteredCherenkovMatrix(), dEdX);
  ublas::vector<double> allChRScat =ublas::prod(*ChFl.GetRayleighScatteredCherenkovMatrix(), dEdX);
  ublas::vector<double> allScat = allChMScat + allChRScat;

  vector<double> groundLevelDistance(1, 0.);
  vector<double> groundLevelScattered(1, 0.);
  vector<double> groundLevelFluo(1, 0.);
  vector<double> xmaxLevelDistance(1, 0.);
  vector<double> xmaxLevelScattered(1, 0.);
  vector<double> xmaxLevelFluo(1, 0.);

  int iGround = n - 1;
  int iXmax = -1;
  double minDist = 0;
  for (unsigned int i = 0; i < depth.size(); ++i) {
    if (i == 0 || fabs(Xmax - depth[i]) < minDist) {
      iXmax = i;
      minDist = fabs(Xmax - depth[i]);
    }
  }

  zeta = 0.05*degree;
  while (zeta < 10.*degree) {
    ChFl.SetZeta(zeta);
    ublas::vector<double> vFluo = ublas::prod(*ChFl.GetFluorescenceMatrix(), dEdX);
    ublas::vector<double> vChMScat = ublas::prod(*ChFl.GetMieScatteredCherenkovMatrix(), dEdX);
    ublas::vector<double> vChRScat = ublas::prod(*ChFl.GetRayleighScatteredCherenkovMatrix(), dEdX);

    groundLevelDistance.push_back(ChFl.GetZetaDistance(iGround));
    groundLevelFluo.push_back(vFluo[iGround] / allFluo[iGround]);
    groundLevelScattered.push_back((vChMScat[iGround] + vChRScat[iGround]) / allScat[iGround]);

    if (iXmax >= 0) {
      xmaxLevelDistance.push_back(ChFl.GetZetaDistance(iXmax));
      xmaxLevelFluo.push_back(vFluo[iXmax] / allFluo[iXmax]);
      xmaxLevelScattered.push_back((vChMScat[iXmax] + vChRScat[iXmax]) / allScat[iXmax]);
    }

    if (fVerbosity > -1)
      cout << "calc. zeta " << zeta/degree << endl;

    zeta += 0.1*degree;
  }

  double ageAtGround = floor(ShowerAge(depth[iGround], Xmax)*10. + 0.5) / 10.;
  const string groundAge = (format("s=%3.1f") % ageAtGround).str();

  TGraph groundFluoGraph(groundLevelDistance.size(),
                         &groundLevelDistance.front(),
                         &groundLevelFluo.front());
  groundFluoGraph.SetName("groundFluo");
  groundFluoGraph.SetTitle(groundAge.c_str());
  groundFluoGraph.Write();

  TGraph groundScatteredGraph(groundLevelDistance.size(),
                              &groundLevelDistance.front(),
                              &groundLevelScattered.front());
  groundScatteredGraph.SetName("groundScattered");
  groundScatteredGraph.SetTitle(groundAge.c_str());
  groundScatteredGraph.Write();

  TGraph xmaxFluoGraph(xmaxLevelDistance.size(),
                       &xmaxLevelDistance.front(),
                       &xmaxLevelFluo.front());
  xmaxFluoGraph.SetName("xmaxFluo");
  xmaxFluoGraph.SetTitle("s=1");
  xmaxFluoGraph.Write();

  TGraph xmaxScatteredGraph(xmaxLevelDistance.size(),
                            &xmaxLevelDistance.front(),
                            &xmaxLevelScattered.front());
  xmaxScatteredGraph.SetName("xmaxScattered");
  xmaxScatteredGraph.SetTitle("s=1");
  xmaxScatteredGraph.Write();
}


RootCFMatrixOutput::~RootCFMatrixOutput()
{
  fOutputFile->Write();
  fOutputFile->Close();

  delete fOutputFile;
}