SiStripDigitizer/test/parametrization.C

0001 #include "TObject.h"
0002 #include "TF1.h"
0003 #include "TMath.h"
0004 #include <iostream>
0005
0006 #include "TkPulseShape.h"
0007
0008 // ROOT script to generate fast and accurate pulse shape functions
0009 // use (in ROOT interpreter) :
0010 //.L parametrization.C+
0011 // parametrizePulse deconv(1)
0012 // parametrizePulse peak(2)
0013 // deconv.generateCode(-30,35,0.1); > deconv.C
0014 // peak.generateCode(-60,400,0.5); > peak.C
0015
0016 class parametrizePulse {
0017  public:
0018   // constructor: takes the mode as input
0019   // mode 1 = deconvolution mode
0020   // mode 2 = peak mode
0021   parametrizePulse(int fitMode = 1);
0022   // destructor
0023   virtual ~parametrizePulse();
0024   // returns a pointer to the theoretical pulse function
0025   TF1* getTheoreticalPulse() const;
0026   // gives the correction factor for a given offset
0027   float getCorrectionFactor(float offset) const;
0028   // returns an array evaluated in a given range with a given step
0029   void getArray(float low, float high, float step);
0030   // generates code to evaluate the function
0031   void generateCode(float low, float high, float step);
0032  private:
0033   // member functions
0034   // members
0035   int fitMode_;
0036 };
0037
0038 parametrizePulse::parametrizePulse(int fitMode)
0039 {
0040   fitMode_ = fitMode; // 1 = deconvolution, 2 = peak
0041 }
0042
0043 parametrizePulse::~parametrizePulse()
0044 {
0045
0046 }
0047
0048 float parametrizePulse::getCorrectionFactor(float offset) const {
0049     TF1* pulse = getTheoreticalPulse();
0050     if(fitMode_ == 1) { return 1./pulse->Eval(offset); }
0051     else { return 1./pulse->Eval(100+offset); }
0052 }
0053
0054 TF1* parametrizePulse::getTheoreticalPulse() const {
0055     TF1* output = NULL;
0056     if(fitMode_ == 1) {
0057       TF1* deconv_fitter = TkPulseShape::GetDeconvFitter();
0058       deconv_fitter->SetParameters(0,-2.816035506,0.066320437,50,20);
0059       output = deconv_fitter;
0060     }
0061     else {
0062       TF1* peak_fitter = TkPulseShape::GetPeakFitter();
0063       peak_fitter->SetParameters(0,-45.90116379,.056621237,50,20);
0064       output = peak_fitter;
0065     }
0066     return output;
0067 }
0068
0069 void parametrizePulse::getArray(float low, float high, float step)
0070 {
0071   TF1* pulse = getTheoreticalPulse();
0072   double end = high+step/2.;
0073   double maximum = pulse->GetMaximum();
0074   for(float val=low;val<end;val+=step) {
0075     std::cout << "(" << val << ", " << pulse->Eval(val)/maximum << ")" << std::endl;
0076   }
0077 }
0078
0079 void parametrizePulse::generateCode(float low, float high, float step)
0080 {
0081   unsigned int size =0;
0082   TF1* pulse = getTheoreticalPulse();
0083   double end = high+step/2.;
0084   double maximum = pulse->GetMaximum();
0085   double maximumX = pulse->GetMaximumX();
0086   for(float val=low;val<end;val+=step) ++size;
0087   std::cout << "float evaluate(float x) {" << std::endl;
0088   std::cout << "  // Automatically generated using parametrizePulse::generateCode(low="
0089             << low << ", high=" << high << ", step="<< step << ")" << std::endl;
0090   std::cout << "  float valuesArray[" << size << "] = { " ;
0091   size=0;
0092   for(float val=low;val<end;val+=step) {
0093     if(size) std::cout << ", " ;
0094     std::cout << pulse->Eval(val+maximumX)/maximum;
0095     if(!((++size)%5)) std::cout << std::endl << "                           ";
0096   }
0097   std::cout << " };" << std::endl;
0098   std::cout << "  if(x<("<<low<<")) return 0;" << std::endl;
0099   std::cout << "  if(x>("<<high<<")) return 0;" << std::endl;
0100   std::cout << "  return valuesArray[unsigned int(((x-("<<low<<"))/("<<step<<"))+0.5)];" << std::endl;
0101   std::cout << "}" << std::endl;
0102 }
0103