Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Data Structures | File List | Namespace Members | Data Fields | Related Pages

ResponseCurve.cpp

00001 #include "ResponseCurve.h"
00002 #include "ChipConfiguration.h"
00003 
00004 #include "TF1.h"
00005 #include <iostream>
00006 #include <sstream>
00007 using namespace std;
00008 
00009 namespace SctData {
00010     ResponseCurve& ResponseCurve::operator=(const ResponseCurve& r) throw(){
00011     shared_ptr<TF1> f=r.getFunction();
00012     getFunction()->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
00013     return *this;
00014     }
00015     
00016     bool ResponseCurve::operator==(const ResponseCurve& r) const throw(){
00017     shared_ptr<TF1> f1=getFunction();
00018         shared_ptr<TF1> f2=r.getFunction(); 
00019     if (r.getCurveName() != this->getCurveName()) return false;
00020     for (int i=0; i<3; ++i){
00021         if (f1->GetParameter(i)!=f2->GetParameter(i)) return false;
00022     }
00023     return true;
00024     }
00025 
00026     ResponseCurve::ResponseCurve(auto_ptr<TF1> f) throw(LogicError) : 
00027     ptr_function(f), ptr_inverse() {
00028     if (ptr_function.get()==0) throw InvariantViolatedError("ResponseCurve constructor couldn't make TF1", __FILE__, __LINE__);
00029     }
00030 
00031     ResponseCurveMap& ResponseCurveMap::getMap() throw() {
00032     static ResponseCurveMap* bertha = new ResponseCurveMap();
00033     return *bertha;
00034     }
00035     
00036     bool ResponseCurveMap::addToMap(const string& curveName, const int index, ResponseCurve& mode) throw() {
00037     //check that this className isn't already in the map
00038     if (nameMap.find(curveName) != nameMap.end()) return false;
00039     //if the index is already taken, just return
00040     if (indexMap.find(index) != indexMap.end()) return false;
00041     
00042     nameMap[curveName] = index;
00043     indexMap[index] = &mode;
00044     return true;
00045     }
00046     
00047     int ResponseCurveMap::getIndex(const string& curvename) const throw(LogicError){
00048     if (nameMap.find(curvename)==nameMap.end()) {
00049         ostringstream os; os <<"ResponseCurveMap::get couldn't find curve `"<< curvename<<"'";
00050         cerr << os.str() << endl;
00051         cerr << "There are " << nameMap.size() << " known response curve names: ";
00052         for (map<string, int>::const_iterator it = nameMap.begin();
00053          it != nameMap.end() ; 
00054          ++it)
00055         {
00056         cerr << (*it).first << endl;
00057         }
00058         throw InvalidArgumentError(os.str(), __FILE__, __LINE__);
00059     }
00060     return (*nameMap.find(curvename)).second;
00061     }
00062     
00063     auto_ptr<ResponseCurve> ResponseCurveMap::get(const string& curvename) const throw(LogicError) {
00064     return get(getIndex(curvename));
00065     }
00066     
00067     auto_ptr<ResponseCurve> ResponseCurveMap::get(const int index) const throw(LogicError) {
00068     if (indexMap.find(index)==indexMap.end()) {
00069         ostringstream os; os <<"ResponseCurveMap::get couldn't find index `"<< index<<"'";
00070         cerr << os.str() << endl;
00071         cerr << "There are " << indexMap.size() << " known response curve indices: ";
00072         for (map<int, ResponseCurve*>::const_iterator it =indexMap.begin();
00073          it != indexMap.end() ; 
00074          ++it)
00075         {
00076         cerr << (*it).first << endl;
00077         }
00078         throw InvalidArgumentError(os.str(), __FILE__, __LINE__);
00079     }
00080     return (*indexMap.find(index)).second->create();
00081     }
00082 
00083     auto_ptr<ResponseCurve> ResponseCurveMap::get(const ChipConfiguration& chip) const throw(LogicError) {
00084     auto_ptr<ResponseCurve> rc = ResponseCurveMap::getMap().get(chip.getRcFunctionIndex()) ;
00085     for (unsigned ipar=0; ipar<3; ++ipar){
00086         rc->getFunction()->SetParameter(ipar,chip.getRcParam(ipar));
00087     }
00088     return rc;
00089     }
00090 
00091 //Linear
00092     LinearResponseCurve::LinearResponseCurve() throw(LogicError) : 
00093     ResponseCurve(auto_ptr<TF1>(new TF1("LinearResponseCurve", linFn, 0, 1, 3))) {
00094     }
00095 
00096     
00097     shared_ptr<TF1> LinearResponseCurve::getInverseFunction() const throw(LogicError) {
00098     if (ptr_inverse.get()==0) {
00099         shared_ptr<TF1> inv(new TF1("InverseLinearResponseCurve",invLinFn, 0, 1, 3));
00100         ptr_inverse=inv;
00101         if (ptr_inverse.get()==0) throw InvariantViolatedError("LinearResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
00102     }
00103     shared_ptr<TF1> f=getFunction();
00104     ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
00105     return ptr_inverse;
00106     }
00107     
00108     bool LinearResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("LinearResponseCurve",4, *new LinearResponseCurve());
00109     
00110     auto_ptr<ResponseCurve> LinearResponseCurve::create() const throw() {
00111     return auto_ptr<ResponseCurve> ( new LinearResponseCurve() );
00112     }
00113     
00114     double LinearResponseCurve::getGain(const double charge) const throw() {
00115     return getFunction()->GetParameter(1);
00116     }
00117 
00118     double LinearResponseCurve::linFn(double *x, double *par) throw() {
00119     return par[0] + par[1]*x[0];
00120     }
00121     double LinearResponseCurve::invLinFn(double *x, double *par) throw() {
00122     return (x[0]-par[0])/par[1];
00123     }
00124 
00125 //Exponential
00126     ExponentialResponseCurve::ExponentialResponseCurve() throw(LogicError): 
00127     ResponseCurve(auto_ptr<TF1>( new TF1("ExponentialResponseCurve", expFn, 0, 1, 3)))
00128     {
00129     getFunction()->SetParameters(1800,8,-800);
00130     }
00131     
00132 
00133     shared_ptr<TF1> ExponentialResponseCurve::getInverseFunction() const throw(LogicError){
00134     if (ptr_inverse.get()==0) {
00135         shared_ptr<TF1> inv( new TF1("InverseExponentialResponseCurve", invExpFn, 0, 1, 3));
00136         ptr_inverse=inv;
00137         if (ptr_inverse.get()==0) throw InvariantViolatedError("ExponentialResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
00138     }
00139     shared_ptr<TF1> f=getFunction();
00140     ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
00141     return ptr_inverse;
00142     }
00143 
00144     bool ExponentialResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("ExponentialResponseCurve", 3, *new ExponentialResponseCurve());
00145     
00146     auto_ptr<ResponseCurve> ExponentialResponseCurve::create() const throw() {
00147     return auto_ptr<ResponseCurve> (new ExponentialResponseCurve() );
00148     }
00149     
00150     double ExponentialResponseCurve::getGain(const double charge) const throw() {
00151     if (getFunction()->GetParameter(1)==0)
00152         return 0.;
00153     double temp = exp(-charge/getFunction()->GetParameter(1));
00154     return getFunction()->GetParameter(0) * pow(1+temp,-2) * temp / getFunction()->GetParameter(1);
00155     }
00156     
00157     double ExponentialResponseCurve::expFn(double *x, double *par) throw() {
00158     return par[2] + par[0]/(1. + exp(-x[0]/par[1]));
00159     }
00160     double ExponentialResponseCurve::invExpFn(double *x, double *par) throw() {
00161     return -par[1] * log(par[0]/(x[0]-par[2])-1);
00162     }
00163     
00164 //Grillo
00165     GrilloResponseCurve::GrilloResponseCurve() throw (LogicError) :
00166     ResponseCurve(auto_ptr<TF1>(new TF1("GrilloResponseCurve", grilloFn, 0, 1, 3))){
00167     getFunction()->SetParameters(70, 40, 0.067);
00168     }
00169     
00170     shared_ptr<TF1> GrilloResponseCurve::getInverseFunction() const throw(LogicError){
00171     if (ptr_inverse.get()==0) {
00172         shared_ptr<TF1> inv(new TF1("InverseGrilloResponseCurve",invGrilloFn, 0, 1, 3));
00173         ptr_inverse=inv;    
00174         if (ptr_inverse.get()==0) throw InvariantViolatedError("GrilloResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
00175     }
00176     shared_ptr<TF1> f=getFunction();
00177     ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
00178     return ptr_inverse;
00179     }
00180 
00181     bool GrilloResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("GrilloResponseCurve", 2, *new GrilloResponseCurve());
00182 
00183     auto_ptr<ResponseCurve> GrilloResponseCurve::create() const throw() {
00184     return auto_ptr<ResponseCurve> (new GrilloResponseCurve()) ;
00185     }
00186 
00187     double GrilloResponseCurve::getGain(const double charge) const throw() {
00188     return getFunction()->GetParameter(1) * pow(1 + pow(getFunction()->GetParameter(2) * charge, 2), -3./2.);
00189     }
00190     double GrilloResponseCurve::grilloFn(double *x, double *par) throw() {
00191     double numerator=x[0] * par[1]*par[2];
00192     if (numerator==0) return 0.;
00193     double sq=x[0]*par[1]; sq=sq*sq;
00194     return par[0] + numerator/sqrt(sq+par[2]*par[2]);
00195     }
00196     double GrilloResponseCurve::invGrilloFn(double *x, double *par) throw() {
00197     double a1=x[0]-par[0]; a1=a1*a1;
00198     double a2=par[2]*par[2];
00199     return 1./par[1] * sqrt( (a1*a2) / (a1+a2) );
00200     }
00201 
00202 //Quadratic
00203     QuadraticResponseCurve::QuadraticResponseCurve() throw(LogicError) : 
00204     ResponseCurve( auto_ptr<TF1>(new TF1("QuadraticResponseCurve", quadFn, 0, 1, 3))) {
00205     getFunction()->SetParameters(70,50,0);
00206     }
00207     
00208     shared_ptr<TF1> QuadraticResponseCurve::getInverseFunction() const throw(LogicError){
00209     if (ptr_inverse.get()==0) {
00210         shared_ptr<TF1> inv( new TF1("InverseQuadraticResponseCurve",invQuadFn, 0, 1, 3));
00211         ptr_inverse=inv;
00212         if (ptr_inverse.get()==0) throw InvariantViolatedError("QuadraticResponseCurve::getInverseFunction() failed", __FILE__, __LINE__);
00213     }
00214     shared_ptr<TF1> f=getFunction();
00215     ptr_inverse->SetParameters(f->GetParameter(0), f->GetParameter(1), f->GetParameter(2));
00216     return ptr_inverse;
00217     }
00218 
00219     bool QuadraticResponseCurve::inMap = ResponseCurveMap::getMap().addToMap("QuadraticResponseCurve", 1, *new QuadraticResponseCurve());
00220     
00221     auto_ptr<ResponseCurve> QuadraticResponseCurve::create() const throw() {
00222     return auto_ptr<ResponseCurve> (new QuadraticResponseCurve()) ;
00223     }
00224 
00225     double QuadraticResponseCurve::getGain(const double charge) const throw() {
00226     return getFunction()->GetParameter(1) + 2.0*charge*getFunction()->GetParameter(2);
00227     }
00228 
00229     double QuadraticResponseCurve::quadFn(double *x, double *par) throw() {
00230     return par[0] + x[0] * (par[1] + par[2]*x[0]);
00231     }
00232     double QuadraticResponseCurve::invQuadFn(double *x, double *par) throw() {
00233     return ( -par[1] - sqrt( par[1]*par[1] - 4 * par[2] * (par[0]-x[0]))) / (2 * par[2]);
00234     }
00235 }

Generated on Thu Jul 15 09:50:50 2004 for SCT DAQ/DCS Software - C++ by doxygen 1.3.5