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
00038 if (nameMap.find(curveName) != nameMap.end()) return false;
00039
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
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
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
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
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 }