TPNFit.cc

Go to the documentation of this file.

/* 
 *  \class TPNFit
 *
 *  \author: Patrice Verrecchia - CEA/Saclay
 */
 
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
 
#include <iostream>
#include "TVectorD.h"
 
//ClassImp(TPNFit)
 
// Constructor...
TPNFit::TPNFit() {
  fNsamples = 0;
  fNum_samp_bef_max = 0;
  fNum_samp_after_max = 0;
}
 
// Destructor
TPNFit::~TPNFit() {}
 
void TPNFit::init(int nsamples, int firstsample, int lastsample) {
  fNsamples = nsamples;
  fNum_samp_bef_max = firstsample;
  fNum_samp_after_max = lastsample;
  //printf("nsamples=%d firstsample=%d lastsample=%d\n",nsamples,firstsample,lastsample);
 
  if (fNsamples > NMAXSAMP2)
    printf("number of PN samples exceed maximum\n");
 
  for (int k = 0; k < NMAXSAMP2; k++)
    t[k] = (double)k;
 
  return;
}
 
double TPNFit::doFit(int maxsample, double *adc) {
  double chi2 = 0.;
  ampl = 0.;
  timeatmax = 0.;
 
  //printf("maxsample=%d\n",maxsample);
  firstsample = maxsample - fNum_samp_bef_max;
  lastsample = maxsample + fNum_samp_after_max;
 
  if (firstsample <= 0)
    return 101;
  if (lastsample >= fNsamples)
    lastsample = fNsamples - 1;
  if (lastsample - firstsample < 1)
    return 102;
  int nval = lastsample - firstsample + 1;
  //printf("firstsample=%d lastsample=%d nval=%d\n",
  //                        firstsample,lastsample,nval);
  int testneg = 0;
  for (int kn = firstsample; kn <= lastsample; kn++) {
    //printf("adc[%d]=%f\n",kn,adc[kn]);
    if (adc[kn] < 0.)
      testneg = 1;
  }
  if (testneg == 1)
    return 103;
 
  for (int i = firstsample; i <= lastsample; i++) {
    val[i - firstsample] = adc[i];
    fv1[i - firstsample] = 1.;
    fv2[i - firstsample] = (double)(i);
    fv3[i - firstsample] = ((double)(i)) * ((double)(i));
  }
 
  TVectorD y(nval, val);
  //y.Print();
  TVectorD f1(nval, fv1);
  //f1.Print();
  TVectorD f2(nval, fv2);
  //f2.Print();
  TVectorD f3(nval, fv3);
  //f3.Print();
 
  double bj[3];
  bj[0] = f1 * y;
  bj[1] = f2 * y;
  bj[2] = f3 * y;
  TVectorD b(3, bj);
  //b.Print();
 
  double aij[9];
  aij[0] = f1 * f1;
  aij[1] = f1 * f2;
  aij[2] = f1 * f3;
  aij[3] = f2 * f1;
  aij[4] = f2 * f2;
  aij[5] = f2 * f3;
  aij[6] = f3 * f1;
  aij[7] = f3 * f2;
  aij[8] = f3 * f3;
  TMatrixD a(3, 3, aij);
  //a.Print();
 
  double det1;
  a.InvertFast(&det1);
  //a.Print();
 
  TVectorD c = a * b;
  //c.Print();
 
  double *par = c.GetMatrixArray();
  //printf("par0=%f par1=%f par2=%f\n",par[0],par[1],par[2]);
  if (par[2] != 0.) {
    timeatmax = -par[1] / (2. * par[2]);
    ampl = par[0] - par[2] * (timeatmax * timeatmax);
  }
  //printf("amp=%f time=%f\n",amp_max,timeatmax);
 
  if (ampl <= 0.) {
    ampl = adc[maxsample];
    return 1.;
  }
 
  if ((int)timeatmax > lastsample)
    return 103;
  if ((int)timeatmax < firstsample)
    return 103;
 
  return chi2;
}