|
|
|
#include <TFParams.h>
Public Member Functions | |
| double | computePulseWidth (int, double, double) |
| void | copie_colonne_mat (matrice, matrice, int) |
| void | diff_mat (matrice, matrice, matrice) |
| double | f3deg (int, double parom[dimout], double mask[dimmat], double adcpj[dimmat], double errpj[dimmat][dimmat]) |
| double | fitpj (double **, double *, double **, double noise_val, int debug) |
| double | inv3x3 (double a[3][3], double b[3][3]) |
| double | inverpj (int, double g[dimmat][dimmat], double ginv[dimmat][dimmat]) |
| void | inverse_mat (matrice, matrice) |
| double | lastShape (Double_t *, Double_t *) |
| double | lastShape2 (Double_t *, Double_t *) |
| double | mixShape (Double_t *, Double_t *) |
| double | parab (double *, Int_t, Int_t, double *) |
| Double_t | polfit (Int_t ns, Int_t imax, Double_t par3d[dimout], Double_t errpj[dimmat][dimmat], double *) |
| void | print_mat (matrice) |
| void | print_mat_nk (matrice, int) |
| void | produit_mat (matrice, matrice, matrice) |
| void | produit_mat_int (matrice, matrice, matrice) |
| double | pulseShapepj (Double_t *, Double_t *) |
| double | pulseShapepj2 (Double_t *, Double_t *) |
| void | set_const (int, int, int, double, double, int) |
| void | somme_mat_int (matrice, matrice) |
| void | somme_mat_int_scale (matrice, matrice, double) |
| TFParams (int size=SDIM2, int size_sh=PLSHDIM) | |
| void | transpose_mat (matrice, matrice) |
| void | zero_mat (matrice) |
| void | zero_mat_nk (matrice, int) |
| ~TFParams () | |
Public Attributes | |
| char | name_mat [10] |
Private Attributes | |
| double | a1ini |
| double | a2ini |
| double | a3ini |
| double | adclu [26] |
| int | METHODE |
| int | nevtmax |
| int | ns |
| int | nsmax |
| int | nsmin |
| double | step_shape |
| double | weight_matrix [10][10] |
Definition at line 47 of file TFParams.h.
| TFParams::TFParams | ( | int | size = SDIM2, |
| int | size_sh = PLSHDIM |
||
| ) |
| TFParams::~TFParams | ( | ) | [inline] |
Definition at line 66 of file TFParams.h.
{};
| double TFParams::computePulseWidth | ( | int | methode, |
| double | alpha_here, | ||
| double | beta_here | ||
| ) |
Definition at line 1252 of file TFParams.cc.
References dt, testEve_cfg::level, evf::evtn::offset(), and tablePrinter::width.
{
// level of amplitude where we calculate the width ( level = 0.5 if at 50 % )
// (level = 0.3 if at 30 % )
double level = 0.30 ;
// fixed parameters
double amplitude = 1.00 ;
double offset = 7.00;
double amp_max = amplitude;
// steps in time
double t_min = offset-4.50;
double t_max = offset+12.50;
int t_step_max = 3000;
double delta_t = (double)((t_max-t_min)/t_step_max);
// Loop over time ( Loop 2 --> get width )
int t_amp_half_flag = 0;
double t_amp_half_min = 999.;
double t_amp_half_max = -999.;
for (int t_step=0; t_step<t_step_max; t_step++){
double t_val = t_min + (double)t_step*delta_t;
double albet = alpha_here*beta_here ;
double dt = t_val-offset ;
double amp =0;
if( methode == 2 ) { // electronic function
if( (t_val-offset) > -albet) {
amp = amplitude*TMath::Power( ( 1 + ( dt / (alpha_here*beta_here) ) ) , alpha_here ) * TMath::Exp(-1.0*(dt/beta_here));
} else {
amp = 1. ;
}
}
if( amp > (amp_max*level) && t_amp_half_flag == 0) {
t_amp_half_flag = 1;
t_amp_half_min = t_val;
}
if( amp < (amp_max*level) && t_amp_half_flag == 1) {
t_amp_half_flag = 2;
t_amp_half_max = t_val;
}
}
// Compute Width
double width = (t_amp_half_max - t_amp_half_min);
return width;
}
Definition at line 655 of file TFParams.cc.
References matrice::coeff, i, j, gen::k, matrice::nb_colonnes, and matrice::nb_lignes.
{
int i,j ;
int k ;
/* resultat de la copie de A dans un vecteur colonne M */
k = 0 ;
for(i=0 ; i< A.nb_lignes; i++) {
for(j=0 ; j < A.nb_colonnes ; j++) {
M.coeff[nk][k] = A.coeff[i][j] ;
printf(" copie nk %d i %d j %d k %d A %e M %e \n ",nk,i,j,k,A.coeff[i][j],
M.coeff[nk][k]);
k++ ;
}
}
return ;
}
Definition at line 636 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int i,j ;
//resultat de la difference A-B = M
if(A.nb_lignes != B.nb_lignes) {
printf( " Erreur : difference de matrices de tailles incompatibles \n ");
M.coeff = NULL ;
return ;
}
M.nb_lignes = A.nb_lignes ;
M.nb_colonnes = A.nb_colonnes ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j < M.nb_colonnes ; j++) {
M.coeff[i][j] = A.coeff[i][j] - B.coeff[i][j] ;
}
}
return ;
}
| double TFParams::f3deg | ( | int | nmxu, |
| double | parom[dimout], | ||
| double | mask[dimmat], | ||
| double | adcpj[dimmat], | ||
| double | errpj[dimmat][dimmat] | ||
| ) |
Definition at line 893 of file TFParams.cc.
References delta, dimmat, f, h, i, gen::k, prof2calltree::l, asciidump::s, mathSSE::sqrt(), matplotRender::t, and tmp.
{
/* */
/* fit 3rd degree polynomial */
/* nmxu = nb of samples in sample data array adcpj[]
parom values of parameters
errpj inverse of the error matrix
fplo3dg uses only the diagonal terms of errpj[][]
*/
int i , k , l , iworst ;
double h , t2 , tm , delta , tmp ;
double xki2 , dif , difmx , deglib ;
double t[dimmat] , f[dimmat][4] ;
double cov[dimmat][dimmat] , bv[4] , invcov[dimmat][dimmat] , s , deter ;
deglib=(double)nmxu - 4. ;
for ( i=0 ; i<nmxu ; i++ ) {
t[i]=i ;
f[i][0]=1. ;
f[i][1]=t[i] ;
f[i][2]=t[i]*t[i] ;
f[i][3]=f[i][2]*t[i] ;
}
/* computation of covariance matrix */
for ( k=0 ; k<4 ; k++ ) {
for ( l=0 ; l<4 ; l++ ) {
s=0. ;
for (i=0 ; i<nmxu ; i++ ) {
s=s+f[i][k]*f[i][l]*errpj[i][i]*mask[i] ;
}
cov[k][l]=s ;
}
s=0. ;
for (i=0 ; i<nmxu ; i++ ) {
s=s+f[i][k]*adcpj[i]*errpj[i][i]*mask[i] ;
}
bv[k]=s ;
}
/* parameters */
deter = inverpj ( 4 , cov , invcov );
for ( k=0 ; k<4 ; k++ ) {
s=0. ;
for ( l=0 ; l<4 ; l++ ) {
s=s+bv[l]*invcov[l][k] ;
}
parom[k]=s ;
}
if( parom[3] == 0. ){
parom[4] = -1000.;
parom[5] = -1000.;
parom[6] = -1000.;
return 1000000.;
}
/* worst hit and ki2 */
xki2=0. ;
difmx=0. ;
for (i=0 ; i<nmxu ; i++ ){
t2=t[i]*t[i] ;
h= parom[0]+parom[1]*t[i]+parom[2]*t2+parom[3]*t2*t[i] ;
dif=(adcpj[i]-h)*mask[i] ;
if(dif > difmx) {
iworst=i ;
difmx=dif ;
}
}
if(deglib > 0.5) xki2=xki2/deglib ;
/* amplitude and maximum position */
delta=parom[2]*parom[2]-3.*parom[3]*parom[1] ;
if(delta > 0.){
delta=sqrt(delta) ;
tm=-(delta+parom[2])/(3.*parom[3]) ;
tmp=(delta-parom[2])/(3.*parom[3]) ;
}
else{
parom[4] = -1000.;
parom[5] = -1000.;
parom[6] = -1000.;
return xki2 ;
}
parom[4]= tm ;
parom[5]= parom[0]+parom[1]*tm+parom[2]*tm*tm+parom[3]*tm*tm*tm ;
parom[6]= tmp ;
// printf("par --------> %f %f %f %f \n",parom[3],parom[2],parom[1],parom[0]);
return xki2 ;
}
| double TFParams::fitpj | ( | double ** | adcval, |
| double * | parout, | ||
| double ** | db_i, | ||
| double | noise_val, | ||
| int | debug | ||
| ) |
Definition at line 39 of file TFParams.cc.
References alpha, beta, funct::C, matrice::coeff, gather_cfg::cout, cree_mat(), cree_mat_prod(), funct::D, delta, dimmat, dt, funct::exp(), i, j, gen::k, funct::log(), nsamp, ntrack, Gflash::par, funct::pow(), matplotRender::t, X, and Gflash::Z.
Referenced by HcalSiPMShape::computeShape().
{
#define dimn 10
#define dimin 10
#define plshdim 300
#define nsamp 10
#define ntrack 500
//#define debug debug1
// ******************************************************************
// * Definitions of variables used in the routine
// ******************************************************************
double a1,a2,a3,a1i,a2i,a3i,b1,b2;
int iter,nevt;
double errpj[dimmat][dimmat] ;
double bi[ntrack][2],dbi[ntrack][2];
double zi[ntrack][2] ;
double par3degre[3] ;
int ioktk[ntrack],iokchi2[ntrack],nk,nborn_min=0,nborn_max=0;
double cti[ntrack][6],dm1i[ntrack][4];
double aiter[10][3];
double par[4],tsig[1];
double amp,delta[nsamp],delta2,fun;
double num_fit_min[ntrack],num_fit_max[ntrack] ;
int i,j,k,imax[ntrack];
double noise_initialvalue,one_over_noisesq ;
double ampmax[ntrack],dt,t;
double chi2, chi2s, da1[nsamp], da2[nsamp], db1[nsamp], db2[nsamp] ;
double chi2tot;
double fact1,fact2,expo;
double albet,dtsbeta,variab,alpha,beta,puiss ;
double unsurs1 ,unsurs2 ;
double fit3 ;
int numb_a,numb_b,numb_ab,numb_b2,numb_x,ndimi,ndimj ;
fun=0; chi2s=0; chi2tot=0;
matrice DA,DAT,BK,DB,DBT,C,CT,D,DM1,CDM1,CDM1CT,Z,CDM1Z,YK,Y,B,X,XINV,RES2 ;
matrice A_CROISS,ZMCT ;
double *amplu ;
amplu = new double[nsamp] ;
parout[0] = 0. ;
parout[1] = 0. ;
parout[2] = 0. ;
//
// Define Error Matrix
//
noise_initialvalue = noise_val ;
one_over_noisesq=1./(noise_initialvalue * noise_initialvalue) ;
for ( i=0 ; i<dimmat ; i++ ) {
for ( j=0 ; j<dimmat ; j++ ) {
errpj[i][j] = 0. ;
}
errpj[i][i]=one_over_noisesq ;
}
//
// Initialisation of fit parameters
//
a1 = a1ini ;
a2 = a2ini ;
a3 = a3ini ;
if( METHODE==2) {
a2 = a3ini ; // for lastshape BETA is the third parameter ( ... ! )
}
a1i = a1 ;
a2i = a2 ;
a3i = a3 ;
if (debug==1){
printf(" ------> __> valeurs de a1 %f a2 %f a3 %f\n",a1,a2,a3) ;
}
for (i=0 ; i<ntrack ; i++) {
iokchi2[i]=1 ;
for (j=0 ; j<2 ; j++ ) {
bi[i][j] = (double)0. ;
dbi[i][j] = (double)0. ;
zi[i][j]=(double)0. ;
cti[i][j]=(double)0. ;
dm1i[i][j]=(double)0. ;
}
}
ndimi = 2 ;
ndimj = 2 ;
numb_a = 2 ;
//
// Matrices initialisation
//
numb_x = 1 ;
numb_b = 2 ;
numb_ab = numb_a*numb_b ;
numb_b2 = numb_b*numb_b ;
DA = cree_mat(numb_a,numb_x) ;
DAT = cree_mat(numb_x,numb_a) ;
BK = cree_mat_prod(DA,DAT) ;
DB = cree_mat(numb_b,numb_x) ;
DBT = cree_mat(numb_x,numb_b) ;
C = cree_mat(numb_a,numb_b) ;
CT = cree_mat(numb_b,numb_a) ;
D = cree_mat_prod(DB,DBT) ;
DM1 = cree_mat_prod(DB,DBT) ;
CDM1 = cree_mat_prod(C,DM1) ;
CDM1CT = cree_mat_prod(CDM1,CT) ;
Z = cree_mat(numb_b,numb_x) ;
CDM1Z =cree_mat_prod(CDM1,Z) ;
YK =cree_mat(numb_a,numb_x) ;
Y =cree_mat(numb_a,numb_x) ;
B = cree_mat_prod(DA,DAT) ;
X = cree_mat_prod(DA,DAT) ;
XINV = cree_mat_prod(DA,DAT) ;
RES2=cree_mat(numb_a,numb_x) ;
A_CROISS = cree_mat(numb_a,numb_x) ;
ZMCT = cree_mat(numb_b,numb_x) ;
// First Loop on iterations //
for (iter=0 ; iter < 6 ; iter++) {
chi2tot=0;
//
// Set zeros for general matrices
//
if (debug==1){
printf(" Debut de l'iteration numero %d \n",iter) ;
}
zero_mat(CDM1Z) ;
zero_mat(Y) ;
zero_mat(CDM1CT) ;
zero_mat(B) ;
zero_mat(X) ;
zero_mat(CDM1) ;
nk = -1 ;
aiter[iter][0] = a1 ;
aiter[iter][1] = a2 ;
if (debug==1){
printf( " resultats injectes a iterations %d \n",iter) ;
printf( " parametre a1 = %f \n",a1) ;
printf( " parametre a2 = %f \n",a2) ;
printf( " chi2 du fit chi2s = %f \n",chi2s) ;
printf(" value de nevtmax _______________ %d \n",nevtmax) ;
}
// Loop on events //
for (nevt=0 ; nevt < nevtmax ; nevt++) {
// B1 = BI[nk,1] est la normalisation du signal
// B2 = BI[nk,2] ewst le dephasage par rapport a une
// fonction centree en zero
// Nous choisissons ici de demarrer avec les resultats
// de l'ajustement parabolique mais il faudra bien
// entendu verifier que cela ne biaise pas le resultat !
// mise a zero des matrices utilisees dans la boucle
zero_mat(Z) ;
zero_mat(YK) ;
zero_mat(BK) ;
zero_mat(C) ;
zero_mat(D) ;
nk=nevt ;
ampmax[nk] = 0. ;
imax[nk] = 0 ;
for ( k = 0 ; k < 10 ; k++) {
amplu[k]=adcval[nevt][k] ;
if (amplu[k] > ampmax[nk] ) {
ampmax[nk] = amplu[k] ;
imax[nk] = k ;
}
}
if( iter == 0 ) {
// start with degree 3 polynomial ....
//fit3 =polfit(ns ,imax[nk] ,par3degre ,errpj ,amplu ) ;
// std::cout << "Poly Fit Param :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl;
// start with parabol
//fit3 = parab(amplu,4,12,par3degre) ;
fit3 = parab(amplu,2,9,par3degre) ;
//std::cout << "Parab Fit Param :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl;
// start with basic initial values
//par3degre[0]= ampmax+ampmax/10. ;
//par3degre[1]= (double)imax[nk]+0.1 ;
//bi[nk][0] = ampmax[nk] ;
//bi[nk][1] = (double)imax[nk] ;
num_fit_min[nevt] = (double)imax[nk] - (double)nsmin ;
num_fit_max[nevt] = (double)imax[nk] + (double)nsmax ;
bi[nk][0] = par3degre[0] ;
bi[nk][1] = par3degre[1] ;
if (debug==1){
printf("---------> depart ampmax[%d]=%f maximum %f tim %f \n",
nk,ampmax[nk],bi[nk][0],bi[nk][1]);
}
} else {
// in other iterations :
// increment bi[][] parameters with bdi[][]
// calculated in previous
// iteration
bi[nk][0] += dbi[nk][0] ;
bi[nk][1] += dbi[nk][1] ;
if (debug==1){
printf("iter %d valeur de max %f et norma %f poly 3 \n",iter,bi[nk][1],bi[nk][0]) ;
}
}
b1 = bi[nk][0] ;
b2 = bi[nk][1] ;
// Loop on samples //
chi2 = 0. ;
ioktk[nk] = 0 ;
ns =nborn_max-nborn_min+1 ;
for (k=0 ; k < 10 ; k++){
//
// calculation of fonction used to fit
//
dt =(double)k - b2 ;
t = (double)k ;
amp = amplu[k] ;
if (debug==1){
printf(" CHECK sample %f ampli %f \n",t,amp) ;
}
//unsurs1 = 1./sig_err ;
//unsurs2 = 1./(sig_err*sig_err) ;
//unsurs1 = 0.1 ;
//unsurs2 = 0.01 ;
unsurs1=1./noise_val;
unsurs2=(1./noise_val)*(1./noise_val);
//
// Pulse shape function used: pulseShapepj
//
nborn_min = (int)num_fit_min[nevt] ;
nborn_max = (int)num_fit_max[nevt] ;
if(k < nborn_min || k > nborn_max ) continue ;
tsig[0] =(double)k ;
if(METHODE==2) {
par[0]= b1 ;
par[1] = b2 ;
par[2] = a1 ;
par[3] = a2 ;
fun = pulseShapepj( tsig , par) ;
}
if (debug==1){
printf(" valeur ampli %f et function %f min %d max %d \n",amp,fun,nsmin,nsmax) ;
printf("min %f max %f \n",num_fit_min[nevt],num_fit_max[nevt]) ;
}
// we need to determine a1,a2 which are global parameters
// and b1, b2 which are parameters for each individual signal:
// b1, b2 = amplitude and time event by event
// a1, a2 = alpha and beta global
// we first begin to calculate the derivatives used in the following calculation
if(METHODE==2){
alpha = a1 ;
beta = a2 ;
albet=alpha*beta;
if(dt > -albet) {
variab = (double)1. + dt/albet ;
dtsbeta = dt/beta ;
expo = exp(-dtsbeta) ;
puiss = pow(variab,alpha) ;
fact1 = puiss*expo ;
db1[k] = unsurs1*fun/b1 ;
fact2 = fun ;
db2[k] = unsurs1*fact2*dtsbeta/(albet*variab) ;
da1[k] = unsurs1*fact2*(log(variab)-dtsbeta/(alpha*variab)) ;
da2[k] = unsurs1*fact2*dtsbeta*dtsbeta/(albet*variab) ;
}
}
delta[k] = (amp - fun)*unsurs1 ;
if (debug==1){
printf(" ------->iter %d valeur de k %d amp %f fun %f delta %f \n",
iter,k,amp,fun,delta[k]) ;
printf(" -----> valeur de k %d delta %f da1 %f da2 %f \n",
k,delta[k],da1[k],da2[k]) ;
}
chi2 = chi2 + delta[k]*delta[k] ;
if (debug==1){
printf(" CHECK chi2 %f deltachi2 %f sample %d iter %d \n",chi2,delta[k]*delta[k],k, iter) ;
}
}
// End Loop on samples //
double wk1wk2 ;
// Start Loop on samples //
for(int k1=nborn_min ; k1<nborn_max+1 ; k1++) {
wk1wk2 = 1. ;
int k2 = k1 ;
DA.coeff[0][0] = da1[k1]*wk1wk2 ;
DA.coeff[1][0] = da2[k1]*wk1wk2 ;
DAT.coeff[0][0]= da1[k2] ;
DAT.coeff[0][1]= da2[k2] ;
DB.coeff[0][0] = db1[k1]*wk1wk2 ;
DB.coeff[1][0] = db2[k1]*wk1wk2 ;
DBT.coeff[0][0]= db1[k2] ;
DBT.coeff[0][1]= db2[k2] ;
// Compute derivative matrix : matrix b[2][2]
produit_mat_int(DA,DAT,BK) ;
// Compute matrix c[2][2]
produit_mat_int(DA,DBT,C) ;
// Compute matrix d[2][2]
produit_mat_int(DB,DBT,D) ;
// Compute matrix y[3] and z[2] depending of delta (amp-fun)
delta2 = delta[k2] ;
somme_mat_int_scale(DA,YK,delta2) ;
somme_mat_int_scale(DB,Z,delta2) ;
ioktk[nk]++ ;
}
// End Loop on samples //
// Remove events with a bad shape
if(ioktk[nk] < 4 ) {
printf(" event rejected because npamp_used = %d \n",ioktk[nk]);
continue ;
}
chi2s = chi2/(2. + (double)ns + 2.) ;
chi2tot+=chi2s;
if (debug==1){
if (nevt==198 || nevt==199){
std::cout << "adc123 pour l'evt " << nevt <<" = "<<adcval[nevt][nborn_min]<<" = "<<adcval[nevt][imax[nevt]]<<" = "<<adcval[nevt][nborn_max]<<std::endl;
std::cout << "chi2s pour l'evt " << nevt <<" = "<< chi2s<<" "<< chi2<<" "<< ns<<" "<<iter<<std::endl;
std::cout << "chi2tot " << nevt <<" = "<< chi2tot<<" "<<iter<<std::endl;
}
}
// Transpose matrix C ---> CT
transpose_mat(C,CT) ;
// Calculate DM1 (inverse of D matrix 2x2)
inverse_mat(D,DM1) ;
// Set matrix product c*d in memory in order to compute variations
// of parameters B at the end of the iteration loop
// the variations of parameters b are dependant of the variations of
// parameters da[1],da[2]
cti[nk][0] = CT.coeff[0][0] ;
cti[nk][1] = CT.coeff[0][1] ;
cti[nk][2] = CT.coeff[1][0] ;
cti[nk][3] = CT.coeff[1][1] ;
dm1i[nk][0] = DM1.coeff[0][0] ;
dm1i[nk][1] = DM1.coeff[0][1] ;
dm1i[nk][2] = DM1.coeff[1][0] ;
dm1i[nk][3] = DM1.coeff[1][1] ;
zi[nk][0] = Z.coeff[0][0] ;
zi[nk][1] = Z.coeff[1][0] ;
// Sum the matrix b and y after every event
for( k=0 ; k< numb_a ; k++) {
Y.coeff[k][0] += YK.coeff[k][0] ;
}
somme_mat_int(BK,B) ;
// Calculate c(d-1)
produit_mat(C,DM1,CDM1) ;
// Compute c(d-1)ct
produit_mat_int(CDM1,CT,CDM1CT);
// Compute c(d-1)z
produit_mat_int(CDM1,Z,CDM1Z) ;
}
// End Loop on events //
// Compute b-cdm1ct
diff_mat(B,CDM1CT,X) ;
inverse_mat(X,XINV) ;
diff_mat(Y,CDM1Z,RES2) ;
// Calculation is now easy for da[0] da[1]
produit_mat(XINV,RES2,A_CROISS) ;
// A la fin, on peut iterer en mesurant l'accroissement a apporter
// des parametres globaux par la formule db[i] = dm1(z-ct*da[i])
for( k=0 ; k< nk+1 ; k++) {
if(METHODE ==2 ) {
ZMCT.coeff[0][0] = zi[k][0] - (cti[k][0]*A_CROISS.coeff[0][0]+
cti[k][1]*A_CROISS.coeff[1][0]) ;
ZMCT.coeff[1][0] = zi[k][1] - (cti[k][2]*A_CROISS.coeff[0][0]+
cti[k][3]*A_CROISS.coeff[1][0]) ;
}
dbi[k][0] = dm1i[k][0]*ZMCT.coeff[0][0] + dm1i[k][1]*ZMCT.coeff[1][0] ;
dbi[k][1] = dm1i[k][2]*ZMCT.coeff[0][0] + dm1i[k][3]*ZMCT.coeff[1][0] ;
if (debug==1){
if( k < 100 ){
printf(" variations de b1= %f et b2= %f \n",dbi[k][0],dbi[k][1]) ;
}
}
db_i[k][0] = bi[k][0]+ dbi[k][0] ;
db_i[k][1] = bi[k][1]+ dbi[k][1] ;
}
// dbi[0] et dbi[1] mesurent les variations a apporter aux
// parametres du signal
a1 += A_CROISS.coeff[0][0] ;
a2 += A_CROISS.coeff[1][0] ;
if (debug==1){
printf(" CHECK croiss coef0: %f croiss coef1: %f iter %d \n",fabs(A_CROISS.coeff[0][0]),fabs(A_CROISS.coeff[1][0]), iter);
}
if(fabs(A_CROISS.coeff[0][0]) < 0.001 && fabs(A_CROISS.coeff[1][0]) < 0.001)
break;
}
// End Loop on iterations //
parout[0] = a1 ;
parout[1] = a2 ;
parout[2] = a3 ;
if (debug==1){
printf( " resultats trouves au bout de %d iterations \n",iter) ;
printf( " parametre a1 = %f \n",a1) ;
printf( " parametre a2 = %f \n",a2) ;
}
if (debug==1){
std::cout << " Final chi2 / NDOF : "<< chi2tot/nevtmax << std::endl;
std::cout << " Final (alpha,beta) : ("<< a1<<","<<a2<<")"<< std::endl;
}
return chi2tot/nevtmax ;
}
| double TFParams::inv3x3 | ( | double | a[3][3], |
| double | b[3][3] | ||
| ) |
Definition at line 1043 of file TFParams.cc.
{
/* a[ligne][colonne] b[ligne][colonne] */
int i , j ;
double deter=0. ;
b[0][0]=a[1][1]*a[2][2]-a[2][1]*a[1][2] ;
b[1][1]=a[0][0]*a[2][2]-a[2][0]*a[0][2] ;
b[2][2]=a[0][0]*a[1][1]-a[0][1]*a[1][0] ;
printf("a[x][x] %e %e %e %e %e %e %e \n",a[0][0],a[1][1],a[0][1],a[1][0],
a[0][0]*a[1][1],a[0][1]*a[1][0],b[2][2]);
b[0][1]=a[2][1]*a[0][2]-a[0][1]*a[2][2] ;
b[0][2]=a[0][1]*a[1][2]-a[1][1]*a[0][2] ;
b[1][0]=a[1][2]*a[2][0]-a[1][0]*a[2][2] ;
b[1][2]=a[1][0]*a[0][2]-a[0][0]*a[1][2] ;
b[2][0]=a[1][0]*a[2][1]-a[1][1]*a[2][0] ;
b[2][1]=a[0][1]*a[2][0]-a[0][0]*a[2][1] ;
deter=a[0][0]*b[0][0] + a[1][0]*b[0][1] + a[2][0]*b[0][2] ;
printf(" deter = %e \n",deter) ;
for ( i=0 ; i<3 ; i++ ) {
for ( j=0 ; j<3 ; j++ ) {
printf(" avant division a[3][3] %d %d %e \n",i,j,a[i][j]) ;
printf(" avant division b[3][3] %d %d %e %e \n",i,j,b[i][j],deter) ;
b[i][j] = b[i][j]/deter ;
printf(" valeur de b[3][3] apres division %d %d %e %e \n",i,j,b[i][j],deter) ;
}
}
return deter ;
}
| double TFParams::inverpj | ( | int | n, |
| double | g[dimmat][dimmat], | ||
| double | ginv[dimmat][dimmat] | ||
| ) |
Definition at line 981 of file TFParams.cc.
References dimmat, g, i, j, gen::k, n, csvReporter::r, asciidump::s, and mathSSE::sqrt().
{
/* */
/* inversion d une matrice symetrique definie positive de taille n */
/* J.P. Pansart Novembre 99 */
/* */
int i , j , k , jj ;
double r , s ;
double deter=0 ;
double al[dimmat][dimmat] , be[dimmat][dimmat] ;
/* initialisation */
for( i=0 ; i<n ; i++ ) {
for ( j=0 ; j<n ; j++ ) {
al[i][j] = 0. ;
be[i][j] = 0. ;
}
}
/* decomposition en vecteurs sur une base orthonormee */
al[0][0] = sqrt( g[0][0] ) ;
for ( i=1 ; i<n ; i++ ) {
al[i][0] = g[0][i] / al[0][0] ;
for ( j=1 ; j<=i ; j++ ) {
s=0. ;
for ( k=0 ; k<=j-1 ; k++ ) {
s=s+ al[i][k] * al[j][k] ;
}
r= g[i][j] - s ;
if( j < i ) al[i][j] = r/al[j][j] ;
if( j == i ) al[i][j] = sqrt ( r) ;
}
}
/* inversion de la matrice al */
be[0][0] = 1./al[0][0] ;
for ( i=1 ; i<n ; i++ ) {
be[i][i] = 1. / al[i][i] ;
for ( j=0 ; j<i ; j++ ) {
jj=i-j-1 ;
s=0. ;
for ( k=jj+1 ; k<=i ; k++ ) {
s=s+ be[i][k] * al[k][jj] ;
}
be[i][jj]=-s/al[jj][jj] ;
}
}
/* calcul de la matrice ginv */
for ( i=0 ; i<n ; i++ ) {
for ( j=0 ; j<n ; j++ ) {
s=0. ;
for ( k=0 ; k<n ; k++ ) {
s=s+ be[k][i] * be[k][j] ;
}
ginv[i][j]=s ;
// if (debug==1){
//printf("valeur de la matrice %d %d %f \n",i,j,ginv[i][j]) ;
//}
}
}
return deter ;
}
Definition at line 803 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, and matrice::nb_lignes.
{
/* A[ligne][colonne] B[ligne][colonne] */
int i , j ;
double deter=0. ;
/* M est la matrice inverse de A */
if(A.nb_lignes != A.nb_colonnes) {
printf( " attention matrice non inversible !!!! %d lignes %d colonnes \n",
A.nb_lignes,A.nb_colonnes) ;
return ;
}
zero_mat(M) ;
if(A.nb_lignes == 2) {
deter = A.coeff[0][0]*A.coeff[1][1] - A.coeff[0][1]*A.coeff[1][0] ;
M.coeff[0][0] = A.coeff[1][1]/deter ;
M.coeff[0][1] = -A.coeff[0][1]/deter ;
M.coeff[1][0] = -A.coeff[1][0]/deter ;
M.coeff[1][1] = A.coeff[0][0]/deter ;
}
else if(A.nb_lignes == 3) {
M.coeff[0][0]=A.coeff[1][1]*A.coeff[2][2]-A.coeff[2][1]*A.coeff[1][2] ;
M.coeff[1][1]=A.coeff[0][0]*A.coeff[2][2]-A.coeff[2][0]*A.coeff[0][2] ;
M.coeff[2][2]=A.coeff[0][0]*A.coeff[1][1]-A.coeff[0][1]*A.coeff[1][0] ;
M.coeff[0][1]=A.coeff[2][1]*A.coeff[0][2]-A.coeff[0][1]*A.coeff[2][2] ;
M.coeff[0][2]=A.coeff[0][1]*A.coeff[1][2]-A.coeff[1][1]*A.coeff[0][2] ;
M.coeff[1][0]=A.coeff[1][2]*A.coeff[2][0]-A.coeff[1][0]*A.coeff[2][2] ;
M.coeff[1][2]=A.coeff[1][0]*A.coeff[0][2]-A.coeff[0][0]*A.coeff[1][2] ;
M.coeff[2][0]=A.coeff[1][0]*A.coeff[2][1]-A.coeff[1][1]*A.coeff[2][0] ;
M.coeff[2][1]=A.coeff[0][1]*A.coeff[2][0]-A.coeff[0][0]*A.coeff[2][1] ;
deter=A.coeff[0][0]*M.coeff[0][0]+A.coeff[1][0]*M.coeff[0][1]
+A.coeff[2][0]*M.coeff[0][2] ;
for ( i=0 ; i<3 ; i++ ) {
for ( j=0 ; j<3 ; j++ ) M.coeff[i][j] = M.coeff[i][j]/deter ;
}
}
else {
printf(" Attention , on ne peut inverser la MATRICE %d \n",A.nb_lignes) ;
return ;
}
return ;
}
| double TFParams::lastShape | ( | Double_t * | x, |
| Double_t * | par | ||
| ) |
Definition at line 1109 of file TFParams.cc.
References alpha, beta, dt, and funct::exp().
| double TFParams::lastShape2 | ( | Double_t * | x, |
| Double_t * | par | ||
| ) |
Definition at line 1126 of file TFParams.cc.
References alpha, beta, dt, and funct::exp().
| double TFParams::mixShape | ( | Double_t * | x, |
| Double_t * | par | ||
| ) |
Definition at line 1219 of file TFParams.cc.
References alpha, beta, dt, funct::exp(), funct::pow(), and matplotRender::t.
{
Double_t fitval0,fitval;
Double_t alpha,beta,fact,puiss;
Double_t dt,alpha2dt,exponent ;
Double_t b1,b2,alpha2,t ;
b1 = par[0] ;
b2 = par[1] ;
alpha = par[2] ;
alpha2 = par[3] ;
beta = par[4] ;
//
t = x[0] ;
dt= x[0]-b2 ;
//
if(t>0.) {
fact = t/b2 ;
puiss = pow(fact,alpha) ;
fitval0 = puiss*exp(-alpha*dt/b2) ;
}
else
{
fitval0=1. ;
}
dt = x[0] - b2 ;
alpha2dt = dt*alpha2 ;
exponent = -(alpha2dt+(exp(-alpha2dt)-1.))/beta ;
fitval = b1*fitval0*exp(exponent) ;
return fitval;
}
| double TFParams::parab | ( | double * | , |
| Int_t | , | ||
| Int_t | , | ||
| double * | |||
| ) |
Definition at line 1180 of file TFParams.cc.
{
/* Now we calculate the parabolic adjustement in order to get */
/* maximum and time max */
double denom,dt,amp1,amp2,amp3 ;
double ampmax = 0. ;
int imax = 0 ;
int k ;
/*
*/
for ( k = nmin ; k < nmax ; k++) {
if (ampl[k] > ampmax ) {
ampmax = ampl[k] ;
imax = k ;
}
}
amp1 = ampl[imax-1] ;
amp2 = ampl[imax] ;
amp3 = ampl[imax+1] ;
denom=2.*amp2-amp1-amp3 ;
/* */
if(denom>0.){
dt =0.5*(amp3-amp1)/denom ;
}
else {
//printf("denom =%f\n",denom) ;
dt=0.5 ;
}
/* */
/* ampmax correspond au maximum d'amplitude parabolique et dt */
/* decalage en temps par rapport au sample maximum soit k + dt */
parout[0] =amp2+(amp3-amp1)*dt*0.25 ;
parout[1] = (double)imax + dt ;
parout[2] = (double)imax ;
return denom ;
}
| Double_t TFParams::polfit | ( | Int_t | ns, |
| Int_t | imax, | ||
| Double_t | par3d[dimout], | ||
| Double_t | errpj[dimmat][dimmat], | ||
| double * | adcpj | ||
| ) |
Definition at line 847 of file TFParams.cc.
References dimmat, dimout, h, and i.
{
double val , val2 , val3 , adfmx[dimmat] , parfp3[dimout] ;
double ius[dimmat], maskp3[dimmat] ;
double deglib,fit3,tm,h,xki2 ;
int i ,nus ,ilow,isup ;
val=adcpj[imax] ;
val2=val/2. ;
val3=val/3. ;
ilow=0 ;
isup=ns ;
deglib=-4. ;
for (i=0 ; i<ns ; i++ ){
deglib=deglib+1. ;
ius[i] = 1. ;
if((adcpj[i] < val3) && (i<imax) ){
ilow=i ;
}
if(adcpj[i] > val2 ){
isup=i ;
}
}
ilow=ilow+1 ;
if(ilow == imax )ilow=ilow-1 ;
if(isup-ilow < 3) isup=ilow+3 ;
nus=0 ;
for(i=ilow ; i<=isup ; i++){
adfmx[nus]=adcpj[i] ;
maskp3[nus] =0. ;
if(ius[i] == 1) {
maskp3[nus]=1. ;
nus=nus+1 ;
}
}
if(nus < 4) return 10000. ;
xki2 = f3deg ( nus , parfp3 , maskp3 , adfmx , errpj ) ;
tm= parfp3[4] ;
h=parfp3[5] ;
tm=tm+(double)ilow ;
par3d[0] = h ;
par3d[1] = tm ;
fit3 = xki2 ;
return fit3 ;
}
| void TFParams::print_mat | ( | matrice | M | ) |
Definition at line 765 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int i,j ;
if( M.coeff == NULL)
printf(" erreur : affichage d'une matrice vide \n") ;
printf(" m_nli %d M_ncol %d \n",M.nb_lignes,M.nb_colonnes) ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++)
printf(" MATRICE i= %d j= %d ---> %e \n",i,j,M.coeff[i][j]) ;
}
//printf(" apres passage d'impression \n") ;
return ;
}
| void TFParams::print_mat_nk | ( | matrice | M, |
| int | nk | ||
| ) |
Definition at line 792 of file TFParams.cc.
References matrice::coeff, j, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int j ;
if( M.coeff == NULL)
printf(" erreur : affichage d'une matrice vide \n") ;
printf(" nk = %d m_nli %d M_ncol %d \n",nk,M.nb_lignes,M.nb_colonnes) ;
for(j=0 ; j< M.nb_colonnes ; j++)
printf(" MATRICE nk= %d j= %d ---> %e \n",nk,j,M.coeff[nk][j]) ;
printf(" apres passage d'impression \n") ;
return ;
}
Definition at line 595 of file TFParams.cc.
References matrice::coeff, i, j, gen::k, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int i,j,k ;
// resultat du produit A*B = M
if(A.nb_colonnes != B.nb_lignes) {
printf( " Erreur : produit de matrices de tailles incompatibles \n ");
M.coeff = NULL ;
return ;
}
M.nb_lignes = A.nb_lignes ;
M.nb_colonnes = B.nb_colonnes ;
zero_mat(M) ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++) {
for(k=0 ; k< A.nb_colonnes; k++){
M.coeff[i][j] += A.coeff[i][k]*B.coeff[k][j] ;
}
}
}
return ;
}
Definition at line 617 of file TFParams.cc.
References matrice::coeff, i, j, gen::k, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int i,j,k ;
if(A.nb_colonnes != B.nb_lignes) {
printf( " Erreur : produit de matrices de tailles incompatibles \n ");
M.coeff = NULL ;
return ;
}
M.nb_lignes = A.nb_lignes ;
M.nb_colonnes = B.nb_colonnes ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++) {
for(k=0 ; k< A.nb_colonnes; k++){
M.coeff[i][j] += A.coeff[i][k]*B.coeff[k][j] ;
}
}
}
return ;
}
| double TFParams::pulseShapepj | ( | Double_t * | x, |
| Double_t * | par | ||
| ) |
Definition at line 1072 of file TFParams.cc.
References alpha, beta, dt, funct::exp(), h, and funct::pow().
{
Double_t fitfun;
Double_t ped, h, tm, alpha, beta;
Double_t dt, dtsbeta, albet, variab, puiss ;
Double_t b1,b2,a1,a2 ;
b1 = par[0] ;
b2 = par[1] ;
a1 = par[2] ;
a2 = par[3] ;
ped = 0. ;
h = b1 ;
tm = b2 ;
alpha = a1 ;
beta = a2 ;
dt= x[0] - tm ;
//printf(" par %f %f %f %f dt = %f albet = %f",b1,b2,a1,a2,dt,albet) ;
albet = alpha*beta ;
if( albet <= 0 )return( (Double_t)0. );
if(dt > -albet) {
dtsbeta=dt/beta ;
variab=1.+dt/albet ;
puiss=pow(variab,alpha);
fitfun=h*puiss*exp(-dtsbeta) + ped;
//printf(" dt = %f h = %f puiss = %f exp(-dtsbeta) %f \n",dt,h,puiss,
// exp(-dtsbeta)) ;
}
else {
fitfun = ped;
}
return fitfun;
}
| Double_t TFParams::pulseShapepj2 | ( | Double_t * | x, |
| Double_t * | par | ||
| ) |
Definition at line 1145 of file TFParams.cc.
References alpha, beta, dt, funct::exp(), h, and funct::pow().
{
Double_t fitfun;
Double_t ped, h, tm, alpha, beta;
Double_t dt, dtsbeta, albet, variab, puiss;
Double_t b1,b2,a1,a2 ;
b1 = par[0] ;
b2 = par[1] ;
a1 = par[2] ;
a2 = par[3] ;
ped = 0. ;
h = b1 ;
tm = b2 ;
alpha = a1 ;
beta = a2 ;
dt= x[0] ;
albet = alpha*beta ;
if( albet <= 0 )return( (Double_t)0. );
if(dt > -albet) {
dtsbeta=dt/beta ;
variab=1.+dt/albet ;
puiss=pow(variab,alpha);
fitfun=h*puiss*exp(-dtsbeta) + ped;
}
else {
fitfun = ped;
}
/* printf( "fitfun %f %f %f %f, %f %f %f\n", ped, h, tm, alpha, beta, *x, fitfun ); */
return fitfun;
}
| void TFParams::set_const | ( | int | n_samples, |
| int | sample_min, | ||
| int | sample_max, | ||
| double | alpha, | ||
| double | beta, | ||
| int | nevtmaximum | ||
| ) |
Definition at line 581 of file TFParams.cc.
References alpha, beta, and SDIM2.
Referenced by HcalSiPMShape::computeShape().
{
/*------------------------------------------------------------------------*/
ns = n_samples;
nsmin = sample_min ;
nsmax = sample_max ;
nevtmax = nevtmaximum ;
a1ini = alpha ;
a2ini = 0.0 ;
a3ini = beta ;
step_shape = .04;
METHODE = 2;
if(ns > SDIM2) printf("warning: NbOfsamples exceed maximum\n");
}
Definition at line 672 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, matrice::nb_lignes, and NULL.
{
int i,j;
/* resultat de la somme integree M += A */
if(A.nb_lignes != M.nb_lignes) {
printf( " Erreur : somme de matrices de tailles incompatibles \n ");
M.coeff = NULL ;
return ;
}
M.nb_lignes = A.nb_lignes ;
M.nb_colonnes = A.nb_colonnes ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++)
M.coeff[i][j] += A.coeff[i][j] ;
}
return ;
}
Definition at line 689 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, and matrice::nb_lignes.
{
int i,j ;
M.nb_lignes = A.nb_lignes ;
M.nb_colonnes = A.nb_colonnes ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[i][j] += A.coeff[i][j]*delta ;
}
return ;
}
Definition at line 699 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, and matrice::nb_lignes.
| void TFParams::zero_mat | ( | matrice | M | ) |
Definition at line 778 of file TFParams.cc.
References matrice::coeff, i, j, matrice::nb_colonnes, and matrice::nb_lignes.
{
int i,j ;
for(i=0 ; i< M.nb_lignes; i++) {
for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[i][j]=0. ;
}
return ;
}
| void TFParams::zero_mat_nk | ( | matrice | M, |
| int | nk | ||
| ) |
Definition at line 786 of file TFParams.cc.
References matrice::coeff, j, and matrice::nb_colonnes.
{
int j ;
for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[nk][j]=0. ;
return ;
}
double TFParams::a1ini [private] |
Definition at line 55 of file TFParams.h.
double TFParams::a2ini [private] |
Definition at line 56 of file TFParams.h.
double TFParams::a3ini [private] |
Definition at line 57 of file TFParams.h.
double TFParams::adclu[26] [private] |
Definition at line 59 of file TFParams.h.
int TFParams::METHODE [private] |
Definition at line 61 of file TFParams.h.
| char TFParams::name_mat[10] |
Definition at line 79 of file TFParams.h.
int TFParams::nevtmax [private] |
Definition at line 54 of file TFParams.h.
int TFParams::ns [private] |
Definition at line 51 of file TFParams.h.
int TFParams::nsmax [private] |
Definition at line 53 of file TFParams.h.
int TFParams::nsmin [private] |
Definition at line 52 of file TFParams.h.
double TFParams::step_shape [private] |
Definition at line 58 of file TFParams.h.
double TFParams::weight_matrix[10][10] [private] |
Definition at line 60 of file TFParams.h.
1.7.3