UnivTimeKGTester.cc

Go to the documentation of this file.

#include "UnivTimeKGTester.h"

#include <utl/ErrorLogger.h>

#include <tls/Atmosphere.h>
#include <tls/UnivParam.h>
#include <tls/UnivTimeKGLogNormal.h>
#include <tls/UnivTimeKGGamma.h>

#include <iostream>

using namespace fwk;
using namespace utl;
using namespace std;

namespace UnivTimeKGTester {


  VModule::ResultFlag
  UnivTimeKGTester::Init()
  {
    return eSuccess;
  }


  VModule::ResultFlag
  UnivTimeKGTester::Run(evt::Event& /*event*/)
  {
    // initialize time models and fill parameters from xml files
    vector<UnivTimeKG::TimeModel*> timeModels;

    for (int icomp = 0; icomp < 4; ++icomp) {
      timeModels.push_back(new UnivTimeKG::GammaTimeModel(icomp));
      //timeModels.push_back(new UnivTimeKG::LogNormalTimeModel(icomp));
    }

    UnivParamNS::UnivParam fpar(0); // 0: Model for Water Cherenkov Detector
    AtmosphereNS::Atmosphere fatm("monthly");
    fatm.SetCurrentAtmosphere(1);

    const double gpercm2 = utl::gram / utl::cm2;

    double density = 1.04e-3 * gram / cm3;
    double xmax_eg = 800 * gpercm2;
    //double xmax_mu = 500 * gpercm2;
    double Nmu = 1.0;
    double logE = 19.0;
    double theta = 36 * degree;

    double height = 1400 * meter;
    int month = 4;
    int sigComp = 1;  // 1: pure electromagnetic signal

    double tracebinning = 25 * nanosecond;

    double rho = 1000 * meter;
    double psi = 90 * degree;

    double DX = fatm.Get_DX_DL_i(rho / cm, psi, xmax_eg / gpercm2, theta, height / cm, false, true) * gpercm2;

    //double integratedSignal = 30.;
    double integratedSignal = fpar.GetSignal(rho / cm, xmax_eg / gpercm2, logE,
                              theta, psi, density / (gram / cm3), height / cm,
                              Nmu, sigComp, month);

    // time referred to plane front arrival
    double planeFrontTime = 400 * nanosecond;
    // arrival time of first particle (difference between curved and plane front, depending on X0)
    double fpTime = 50 * nanosecond;

    timeModels[sigComp]->setShapeParameters(DX / gpercm2, rho, psi, theta, logE);
    // total signal of this component in the station
    double S = timeModels[sigComp]->pdf(planeFrontTime - fpTime) * integratedSignal * tracebinning;

    cout << "S " << S << endl;

    double Scdf = timeModels[sigComp]->cdf(planeFrontTime - fpTime);

    cout << "cdf S (normed)" << Scdf << endl;

    double tinv = fpTime + timeModels[sigComp]->invcdf(Scdf);

    cout << "t " << tinv << endl;

    double pftimeit = fpTime;
    double cdf = 0;

    while (cdf < 0.99) {
      cdf = timeModels[sigComp]->cdf(pftimeit - fpTime);
      cout << "rel. time [ns]: " << pftimeit << ", CDF (normed): " << cdf << endl;
      pftimeit += 25;
    }

    double startTimeProb = timeModels[0]->firstParticlePdfSmeared(planeFrontTime - fpTime, 20);

    cout << "start time probability: " << startTimeProb << endl;

    return eSuccess;
  }


  VModule::ResultFlag
  UnivTimeKGTester::Finish()
  {
    return eSuccess;
  }

}