ResponseCurve.cpp

#include "ResponseCurve.h"
#include "ChipConfiguration.h"

#include "TF1.h"
#include <iostream>
#include <sstream>
using namespace std;

namespace SctData {
    ResponseCurve& ResponseCurve::operator=(const ResponseCurve& r) throw(){
    shared_ptr<TF1> f=r.getFunction();
    getFunction()->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
    return *this;
    }
    
    bool ResponseCurve::operator==(const ResponseCurve& r) const throw(){
    shared_ptr<TF1> f1=getFunction();
        shared_ptr<TF1> f2=r.getFunction(); 
    if (r.getCurveName() != this->getCurveName()) return false;
    for (int i=0; i<3; ++i){
        if (f1->GetParameter(i)!=f2->GetParameter(i)) return false;
    }
    return true;
    }
 
   ResponseCurve::~ResponseCurve() throw() {
        //std::cout << "In ~ResponseCurve() " << __FILE__ << " " << __LINE__ << std::endl;
    this->ptr_inverse = shared_ptr<TF1>();
    //std::cout << "I have just reset the ptr_inverse.   Now I am going to throw!" << std::endl;
    this->ptr_function = shared_ptr<TF1>();
    //std::cout << "I have just reset the ptr_function" << std::endl;
   }

    ResponseCurve::ResponseCurve(auto_ptr<TF1> f) throw(LogicError) : 
    ptr_function(f), ptr_inverse() {
    if (ptr_function.get()==0) throw InvariantViolatedError("ResponseCurve constructor couldn't make TF1", __FILE__, __LINE__);
    }

    ResponseCurveMap& ResponseCurveMap::getMap() throw() {
    static ResponseCurveMap* bertha = new ResponseCurveMap();
    return *bertha;
    }
    
    bool ResponseCurveMap::addToMap(const string& curveName, const int index, ResponseCurve& mode) throw() {
    //check that this className isn't already in the map
    if (nameMap.find(curveName) != nameMap.end()) return false;
    //if the index is already taken, just return
    if (indexMap.find(index) != indexMap.end()) return false;
    
    nameMap[curveName] = index;
    indexMap[index] = &mode;
    return true;
    }
    
    int ResponseCurveMap::getIndex(const string& curvename) const throw(LogicError){
    if (nameMap.find(curvename)==nameMap.end()) {
        ostringstream os; os <<"ResponseCurveMap::get couldn't find curve `"<< curvename<<"'";
        cerr << os.str() << endl;
        cerr << "There are " << nameMap.size() << " known response curve names: ";
        for (map<string, int>::const_iterator it = nameMap.begin();
         it != nameMap.end() ; 
         ++it)
        {
        cerr << (*it).first << endl;
        }
        throw InvalidArgumentError(os.str(), __FILE__, __LINE__);
    }
    return (*nameMap.find(curvename)).second;
    }
    
    auto_ptr<ResponseCurve> ResponseCurveMap::get(const string& curvename) const throw(LogicError) {
    return get(getIndex(curvename));
    }
    
    auto_ptr<ResponseCurve> ResponseCurveMap::get(const int index) const throw(LogicError) {
    if (indexMap.find(index)==indexMap.end()) {
        ostringstream os; os <<"ResponseCurveMap::get couldn't find index `"<< index<<"'";
        cerr << os.str() << endl;
        cerr << "There are " << indexMap.size() << " known response curve indices: ";
        for (map<int, ResponseCurve*>::const_iterator it =indexMap.begin();
         it != indexMap.end() ; 
         ++it)
        {
        cerr << (*it).first << endl;
        }
        throw InvalidArgumentError(os.str(), __FILE__, __LINE__);
    }
    return (*indexMap.find(index)).second->create();
    }

    auto_ptr<ResponseCurve> ResponseCurveMap::get(const ChipConfiguration& chip) const throw(LogicError) {
    auto_ptr<ResponseCurve> rc = ResponseCurveMap::getMap().get(chip.getRcFunctionIndex()) ;
    for (unsigned ipar=0; ipar<3; ++ipar){
        rc->getFunction()->SetParameter(ipar,chip.getRcParam(ipar));
    }
    return rc;
    }

  std::string notTheSameName(const char * part) {
    std::ostringstream s;
    s << part;
    static unsigned int ReneBrunMustDie = 0;
    s << (ReneBrunMustDie++);
    return s.str();
  };

//Linear
    LinearResponseCurve::LinearResponseCurve() throw(LogicError) : 
    ResponseCurve(auto_ptr<TF1>(new TF1(notTheSameName("LinearResponseCurve").c_str(), linFn, 0, 1, 3))) {
    m_deb = rand();
    //std::cout << "In LinearResponseCurve() for m_deb="<<m_deb<<std::endl;
    }

    LinearResponseCurve::~LinearResponseCurve() throw() {
    //std::cout << "In ~LinearResponseCurve() for m_deb="<<m_deb<<std::endl;
    };
    
    shared_ptr<TF1> LinearResponseCurve::getInverseFunction() const throw(LogicError) {
    if (ptr_inverse.get()==0) {
        shared_ptr<TF1> inv(new TF1("InverseLinearResponseCurve",invLinFn, 0, 1, 3));
        ptr_inverse=inv;
        if (ptr_inverse.get()==0) throw InvariantViolatedError("LinearResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
    }
    shared_ptr<TF1> f=getFunction();
    ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
    return ptr_inverse;
    }
    
    bool LinearResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("LinearResponseCurve",4, *new LinearResponseCurve());
    
    auto_ptr<ResponseCurve> LinearResponseCurve::create() const throw() {
    return auto_ptr<ResponseCurve> ( new LinearResponseCurve() );
    }
    
    double LinearResponseCurve::getGain(const double charge) const throw() {
    return getFunction()->GetParameter(1);
    }

    double LinearResponseCurve::linFn(double *x, double *par) throw() {
    return par[0] + par[1]*x[0];
    }
    double LinearResponseCurve::invLinFn(double *x, double *par) throw() {
    return (x[0]-par[0])/par[1];
    }

//Exponential
    ExponentialResponseCurve::ExponentialResponseCurve() throw(LogicError): 
    ResponseCurve(auto_ptr<TF1>( new TF1(notTheSameName("ExponentialResponseCurve").c_str(), expFn, 0, 1, 3)))
    {
    getFunction()->SetParameters(1800,8,-800);
    }
    

    shared_ptr<TF1> ExponentialResponseCurve::getInverseFunction() const throw(LogicError){
    if (ptr_inverse.get()==0) {
        shared_ptr<TF1> inv( new TF1("InverseExponentialResponseCurve", invExpFn, 0, 1, 3));
        ptr_inverse=inv;
        if (ptr_inverse.get()==0) throw InvariantViolatedError("ExponentialResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
    }
    shared_ptr<TF1> f=getFunction();
    ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
    return ptr_inverse;
    }

    bool ExponentialResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("ExponentialResponseCurve", 3, *new ExponentialResponseCurve());
    
    auto_ptr<ResponseCurve> ExponentialResponseCurve::create() const throw() {
    return auto_ptr<ResponseCurve> (new ExponentialResponseCurve() );
    }
    
    double ExponentialResponseCurve::getGain(const double charge) const throw() {
    if (getFunction()->GetParameter(1)==0)
        return 0.;
    double temp = exp(-charge/getFunction()->GetParameter(1));
    return getFunction()->GetParameter(0) * pow(1+temp,-2) * temp / getFunction()->GetParameter(1);
    }
    
    double ExponentialResponseCurve::expFn(double *x, double *par) throw() {
    return par[2] + par[0]/(1. + exp(-x[0]/par[1]));
    }
    double ExponentialResponseCurve::invExpFn(double *x, double *par) throw() {
    return -par[1] * log(par[0]/(x[0]-par[2])-1);
    }
    
//Grillo
    GrilloResponseCurve::GrilloResponseCurve() throw (LogicError) :
    ResponseCurve(auto_ptr<TF1>(new TF1(notTheSameName("GrilloResponseCurve").c_str(), grilloFn, 0, 1, 3))){
    getFunction()->SetParameters(70, 40, 0.067);
    }
    
    shared_ptr<TF1> GrilloResponseCurve::getInverseFunction() const throw(LogicError){
    if (ptr_inverse.get()==0) {
        shared_ptr<TF1> inv(new TF1("InverseGrilloResponseCurve",invGrilloFn, 0, 1, 3));
        ptr_inverse=inv;    
        if (ptr_inverse.get()==0) throw InvariantViolatedError("GrilloResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
    }
    shared_ptr<TF1> f=getFunction();
    ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
    return ptr_inverse;
    }

    bool GrilloResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("GrilloResponseCurve", 2, *new GrilloResponseCurve());

    auto_ptr<ResponseCurve> GrilloResponseCurve::create() const throw() {
    return auto_ptr<ResponseCurve> (new GrilloResponseCurve()) ;
    }

    double GrilloResponseCurve::getGain(const double charge) const throw() {
    return getFunction()->GetParameter(1) * pow(1 + pow(getFunction()->GetParameter(2) * charge, 2), -3./2.);
    }
    double GrilloResponseCurve::grilloFn(double *x, double *par) throw() {
    double numerator=x[0] * par[1]*par[2];
    if (numerator==0) return 0.;
    double sq=x[0]*par[1]; sq=sq*sq;
    return par[0] + numerator/sqrt(sq+par[2]*par[2]);
    }
    double GrilloResponseCurve::invGrilloFn(double *x, double *par) throw() {
    double a1=x[0]-par[0]; a1=a1*a1;
    double a2=par[2]*par[2];
    return 1./par[1] * sqrt( (a1*a2) / (a1+a2) );
    }

//Quadratic
    QuadraticResponseCurve::QuadraticResponseCurve() throw(LogicError) : 
    ResponseCurve( auto_ptr<TF1>(new TF1(notTheSameName("QuadraticResponseCurve").c_str(), quadFn, 0, 1, 3))) {
    getFunction()->SetParameters(70,50,0);
    }
    
    shared_ptr<TF1> QuadraticResponseCurve::getInverseFunction() const throw(LogicError){
    if (ptr_inverse.get()==0) {
        shared_ptr<TF1> inv( new TF1("InverseQuadraticResponseCurve",invQuadFn, 0, 1, 3));
        ptr_inverse=inv;
        if (ptr_inverse.get()==0) throw InvariantViolatedError("QuadraticResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
    }
    shared_ptr<TF1> f=getFunction();
    ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
    return ptr_inverse;
    }

    bool QuadraticResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("QuadraticResponseCurve", 1, *new QuadraticResponseCurve());
    
    auto_ptr<ResponseCurve> QuadraticResponseCurve::create() const throw() {
    return auto_ptr<ResponseCurve> (new QuadraticResponseCurve()) ;
    }

    double QuadraticResponseCurve::getGain(const double charge) const throw() {
    return getFunction()->GetParameter(1) + 2.0*charge*getFunction()->GetParameter(2);
    }

    double QuadraticResponseCurve::quadFn(double *x, double *par) throw() {
    return par[0] + x[0] * (par[1] + par[2]*x[0]);
    }
    double QuadraticResponseCurve::invQuadFn(double *x, double *par) throw() {
    return ( -par[1] - sqrt( par[1]*par[1] - 4 * par[2] * (par[0]-x[0]))) / (2 * par[2]);
    }
}

---