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;
}