FlashHypothesis.cxx

Go to the documentation of this file.

#include "FlashHypothesis.h"

#include <limits>
#include <iostream>

void opdet::FlashHypothesis::SetHypothesisVectorAndErrorVector( std::vector<float> v , std::vector<float> err)
{
  if(err.size()!=0 && err.size()!=v.size())
    throw std::runtime_error("ERROR in FlashHypothesisVectorSetter: Vector sizes not equal");

  _NPEs_Vector = v;
  if(err.size()==0){
    _NPEs_Vector = v;
    _NPEs_ErrorVector.resize(v.size());
    for(size_t i=0; i<_NPEs_Vector.size(); i++)
      _NPEs_ErrorVector[i] = std::sqrt(_NPEs_Vector[i]);
  }
  else
    _NPEs_ErrorVector = err;

}

void opdet::FlashHypothesis::SetHypothesisAndError( size_t i_opdet, float pe , float err)
{
  SetHypothesis(i_opdet,pe);
  if(err>0) SetHypothesisError(i_opdet,err);
  else SetHypothesisError(i_opdet,std::sqrt(pe));
}

void opdet::FlashHypothesis::Normalize(float const& totalPE_target){

  if( GetTotalPEs() < std::numeric_limits<float>::epsilon() ) return;

  const float PE_ratio = totalPE_target/GetTotalPEs();
  for(size_t i_opdet=0; i_opdet<_NPEs_Vector.size(); i_opdet++){
    _NPEs_Vector[i_opdet] *= PE_ratio;
    _NPEs_ErrorVector[i_opdet] = std::sqrt(_NPEs_Vector[i_opdet]);
  }

}

void opdet::FlashHypothesis::Print()
{
  std::cout << "TotalPEs: " << GetTotalPEs() << " +/- " << GetTotalPEsError() << std::endl;
  std::cout << "Vector size: " << GetVectorSize() << std::endl;
  for(size_t i=0; i<GetVectorSize(); i++)
    std::cout << "\t" << i << ": " << GetHypothesis(i) << " +/- " << GetHypothesisError(i) << std::endl;
}

void opdet::FlashHypothesisCollection::SetTotalHypAndPromptFraction(const FlashHypothesis& total, float frac)
{
  CheckFrac(frac);
  _total_hyp = total;
  _prompt_frac = frac;
  _prompt_hyp = total; _prompt_hyp.Normalize( frac*total.GetTotalPEs() );
  _late_hyp = total; _late_hyp.Normalize( (1.-frac)*total.GetTotalPEs() );
}

void opdet::FlashHypothesisCollection::SetPromptHypAndPromptFraction(const FlashHypothesis& prompt, float frac)
{
  CheckFrac(frac);
  _prompt_hyp = prompt;
  _prompt_frac = frac;
  _late_hyp = prompt; _late_hyp.Normalize( (1/frac - 1.)*prompt.GetTotalPEs() );
  _total_hyp = _prompt_hyp + _late_hyp;
}

void opdet::FlashHypothesisCollection::Normalize(float totalPE_target){
  _prompt_hyp.Normalize(totalPE_target*_prompt_frac);
  _late_hyp.Normalize(totalPE_target*(1.-_prompt_frac));
  UpdateTotalHyp();
}

void opdet::FlashHypothesisCollection::CheckFrac(float f)
{
  if ( std::abs(f-0.0) < std::numeric_limits<float>::epsilon() ||
       std::abs(f-1.0) < std::numeric_limits<float>::epsilon() ||
       (f>0.0 && f<1.0) )
    return;

  throw std::runtime_error("ERROR in FlashHypothesisCollection: Input fraction is not in valid range.");
}

void opdet::FlashHypothesisCollection::UpdateTotalHyp()
{
  _total_hyp = _prompt_hyp + _late_hyp;
  const float total_pe = _total_hyp.GetTotalPEs();
  if(total_pe > std::numeric_limits<float>::epsilon())
    _prompt_frac = _prompt_hyp.GetTotalPEs() / total_pe;
  else
    _prompt_frac = 1.;
}

void opdet::FlashHypothesisCollection::Print()
{
  std::cout << "PromptHyp:" << std::endl;
  _prompt_hyp.Print();
  std::cout << "LateHyp:" << std::endl;
  _late_hyp.Print();
  std::cout << "TotalHyp:" << std::endl;
  _total_hyp.Print();
  std::cout << "PromptFraction: " << _prompt_frac << std::endl;

}