#include "AddOns/Analysis/Observables/Normalized_Observable.H" #include "ATOOLS/Org/Exception.H" using namespace ANALYSIS; using namespace ATOOLS; Normalized_Observable::Normalized_Observable(): Primitive_Observable_Base() { } Normalized_Observable:: Normalized_Observable(int type,double xmin,double xmax,int nbins, const std::string &name,const int mode): Primitive_Observable_Base(type,xmin,xmax,nbins), m_mode(mode) { if (!(type&4)) THROW(fatal_error,"Must be initialized with type X+4"); p_obs = new Histogram(m_type,m_xmin,m_xmax,m_nbins,name); p_norm = new Histogram(m_type,m_xmin,m_xmax,m_nbins,name); } Normalized_Observable:: Normalized_Observable(const Normalized_Observable & old): Primitive_Observable_Base(old) { if (old.p_histo) { p_obs = new Histogram(m_type,m_xmin,m_xmax,m_nbins,old.Name()); p_norm = new Histogram(m_type,m_xmin,m_xmax,m_nbins,old.Name()); } } Normalized_Observable::~Normalized_Observable() { if (p_histo) { delete p_obs; delete p_norm; } } void Normalized_Observable::Reset() { p_obs->Reset(); p_norm->Reset(); } Primitive_Observable_Base & Normalized_Observable::operator+=(const Primitive_Observable_Base &obs) { const Normalized_Observable *nobs((const Normalized_Observable*)&obs); if (p_histo) { (*p_obs)+=(*nobs->p_obs); (*p_norm)+=(*nobs->p_norm); if (!m_copied) { double n=ATOOLS::Max(1.0,double(p_obs->Fills())); for (int i=0;iBin(i)==0.0?0.0: p_obs->Bin(i)/p_norm->Bin(i); p_histo->SetBin(i,val); double w2mo(n*p_obs->Bin2(i)), wmo(p_obs->Bin(i)); double sigo(sqrt((w2mo-wmo*wmo)/Max(1.0,n-1.0))); double w2mi(n*p_norm->Bin2(i)), wmi(p_norm->Bin(i)); double sigi(sqrt((w2mi-wmi*wmi)/Max(1.0,n-1.0))); double relerr=p_obs->Bin(i)==0.0?0.0: (sigo/dabs(p_obs->Bin(i))+sigi/dabs(p_norm->Bin(i))); p_histo->SetBin2(i,sqr(val*relerr)); p_histo->SetBinPS(i,p_norm->Bin(i)); p_histo->SetBinPS2(i,sqr(sigi)); } p_histo->SetFills((long int)p_obs->Fills()); } } else { THROW(critical_error,Name()+" has not overloaded the operator"); } return *this; } void Normalized_Observable::EndEvaluation(double scale) { p_obs->MPISync(); p_norm->MPISync(); double n=ATOOLS::Max(1.0,double(p_obs->Fills())); if (m_mode==0) p_obs->Scale(scale*m_nbins/(m_xmax-m_xmin)/n); else p_obs->Scale(scale/n); p_norm->Scale(scale/n); if (!m_copied) { for (int i=0;iBin(i)==0.0?0.0: p_obs->Bin(i)/p_norm->Bin(i); p_histo->SetBin(i,val); double w2mo(n*p_obs->Bin2(i)), wmo(p_obs->Bin(i)); double sigo(sqrt((w2mo-wmo*wmo)/Max(1.0,n-1.0))); double w2mi(n*p_norm->Bin2(i)), wmi(p_norm->Bin(i)); double sigi(sqrt((w2mi-wmi*wmi)/Max(1.0,n-1.0))); double relerr=p_obs->Bin(i)==0.0?0.0: (sigo/dabs(p_obs->Bin(i))+sigi/dabs(p_norm->Bin(i))); p_histo->SetBin2(i,sqr(val*relerr)); p_histo->SetBinPS(i,p_norm->Bin(i)); p_histo->SetBinPS2(i,sqr(sigi)); } p_histo->SetFills((long int)p_obs->Fills()); } } void Normalized_Observable::Restore(double scale) { double n=ATOOLS::Max(1.0,double(p_obs->Fills())); if (m_mode==0) p_obs->Scale(scale*n*(m_xmax-m_xmin)/m_nbins); else p_obs->Scale(scale*n); p_norm->Scale(scale*n); } void Normalized_Observable::Fill (const double &x,const double &y, const double &weight,const double &ntrial) { p_obs->Insert(x,y*weight,ntrial); p_norm->Insert(x,weight,ntrial); } void Normalized_Observable::FillMCB (const double &x,const double &y, const double &weight,const double &ntrial) { p_obs->InsertMCB(x,y*weight,ntrial); p_norm->InsertMCB(x,weight,ntrial); } void Normalized_Observable::FinishMCB() { p_obs->FinishMCB(); p_norm->FinishMCB(); } Primitive_Observable_Base *Normalized_Observable::Copy() const { return new Normalized_Observable(m_type,m_xmin,m_xmax,m_nbins,Name()); }