FlashHypothesis.cxx

Go to the documentation of this file.

#include "FlashHypothesis.h"
#include <limits>

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;
  
}