CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/CalibCalorimetry/EcalLaserAnalyzer/src/TFParams.cc

Go to the documentation of this file.
00001 
00013 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TFParams.h>
00014 #include "TMatrixD.h"
00015 #include "TMath.h"
00016 
00017 #include <iostream>
00018 #include <time.h>
00019 
00020 
00021 //ClassImp(TFParams)
00022 
00023 using namespace std;
00024 
00025 TFParams::TFParams( int size, int size_sh ) {
00026 
00027   //int  sdim = size;
00028   //int plshdim = size_sh;
00029 
00030 for (int i=0 ; i<10 ; i++) {
00031   for (int j=0 ; j<10 ; j++) {
00032       weight_matrix[i][j] = 8.; 
00033     }
00034   } 
00035 
00036 }
00037 
00038 double TFParams::fitpj(double **adcval , double *parout , double **db_i, double noise_val, int debug)
00039 {
00040   
00041 #define dimn 10
00042 #define dimin 10
00043 #define plshdim 300
00044 #define nsamp 10                                            
00045 #define ntrack 500       
00046   //#define debug debug1
00047 
00048   // ******************************************************************
00049   // *  Definitions of variables used in the routine                    
00050   // ******************************************************************
00051   
00052   
00053   double a1,a2,a3,b1,b2;
00054   int iter,nevt;
00055   //double errpj[dimmat][dimmat]  ;
00056   double bi[ntrack][2],dbi[ntrack][2];
00057   double zi[ntrack][2] ;
00058   double par3degre[3] ;
00059   int    ioktk[ntrack],nk,nborn_min=0,nborn_max=0;
00060   double cti[ntrack][6],dm1i[ntrack][4]; 
00061   double par[4],tsig[1];
00062   double amp,delta[nsamp],delta2,fun;
00063   double num_fit_min[ntrack],num_fit_max[ntrack] ;
00064   int i,j,k,imax[ntrack];
00065 
00066   double ampmax[ntrack],dt,t;
00067   double chi2, chi2s, da1[nsamp], da2[nsamp], db1[nsamp], db2[nsamp] ;
00068   double chi2tot;
00069   double fact2;
00070   double albet,dtsbeta,variab,alpha,beta;
00071   double  unsurs1 /*,unsurs2*/ ;
00072 //  double fit3 ;
00073   int numb_a,numb_b,numb_x;
00074 
00075   fun=0; chi2s=0; chi2tot=0;
00076   matrice DA,DAT,BK,DB,DBT,C,CT,D,DM1,CDM1,CDM1CT,Z,CDM1Z,YK,Y,B,X,XINV,RES2 ;
00077   matrice A_CROISS,ZMCT ;
00078 
00079   double *amplu ;
00080   amplu = new double[nsamp] ;
00081   
00082   parout[0] = 0. ;
00083   parout[1] = 0. ;
00084   parout[2] = 0. ;
00085 
00086   //  
00087   //  Initialisation of fit parameters 
00088   //  
00089 
00090   a1 = a1ini ;
00091   a2 = a2ini ;
00092   a3 = a3ini ;
00093   if( METHODE==2) {
00094     a2 = a3ini ;   // for lastshape BETA is the third parameter ( ... ! )
00095   }
00096   if (debug==1){
00097     printf(" ------> __> valeurs de a1 %f a2 %f a3 %f\n",a1,a2,a3) ;
00098   }
00099   for (i=0 ; i<ntrack ; i++) {
00100     for (j=0 ; j<2 ; j++ ) {
00101       bi[i][j] = (double)0. ;
00102       dbi[i][j] = (double)0. ;
00103       zi[i][j]=(double)0. ;
00104       cti[i][j]=(double)0. ;
00105       dm1i[i][j]=(double)0. ;
00106     }
00107   }
00108 
00109   numb_a = 2 ;
00110 
00111 
00112   //
00113   //  Matrices initialisation
00114   //
00115 
00116   numb_x = 1 ;
00117   numb_b = 2 ;
00118   DA = cree_mat(numb_a,numb_x) ;
00119   DAT = cree_mat(numb_x,numb_a) ;
00120   BK = cree_mat_prod(DA,DAT) ;
00121   DB = cree_mat(numb_b,numb_x) ;
00122   DBT = cree_mat(numb_x,numb_b) ;
00123   C = cree_mat(numb_a,numb_b) ;
00124   CT = cree_mat(numb_b,numb_a) ;
00125   D = cree_mat_prod(DB,DBT) ;
00126   DM1 = cree_mat_prod(DB,DBT) ;
00127   CDM1 = cree_mat_prod(C,DM1) ;
00128   CDM1CT = cree_mat_prod(CDM1,CT) ;
00129   Z = cree_mat(numb_b,numb_x) ;
00130   CDM1Z =cree_mat_prod(CDM1,Z) ;
00131   YK =cree_mat(numb_a,numb_x) ;
00132   Y =cree_mat(numb_a,numb_x) ;
00133   B = cree_mat_prod(DA,DAT) ;
00134   X = cree_mat_prod(DA,DAT) ;
00135   XINV = cree_mat_prod(DA,DAT) ;
00136   RES2=cree_mat(numb_a,numb_x) ;
00137   A_CROISS = cree_mat(numb_a,numb_x) ;
00138   ZMCT = cree_mat(numb_b,numb_x) ;
00139 
00140 
00142   // First Loop on iterations //                            
00144  
00145   for (iter=0 ; iter < 6 ; iter++) {
00146 
00147     chi2tot=0;
00148 
00149     //    
00150     //    Set zeros for general matrices
00151     //                                                                      
00152     
00153     if (debug==1){
00154       printf(" Debut de l'iteration numero %d \n",iter) ;
00155     }
00156     zero_mat(CDM1Z) ;
00157     zero_mat(Y) ;
00158     zero_mat(CDM1CT) ;
00159     zero_mat(B) ;
00160     zero_mat(X) ;
00161     zero_mat(CDM1) ;
00162 
00163 
00164     nk = -1 ;   
00165     if (debug==1){
00166       printf( " resultats injectes a iterations %d \n",iter) ;
00167       printf( " parametre a1 = %f \n",a1) ;
00168       printf( " parametre a2 = %f \n",a2) ;
00169       printf( " chi2 du fit chi2s = %f \n",chi2s) ;
00170       
00171       printf(" value de nevtmax _______________ %d \n",nevtmax) ;
00172     }
00173     
00174     
00176     //  Loop on events //
00178     
00179     for (nevt=0 ; nevt < nevtmax ; nevt++) {
00180                  
00181       //       B1 = BI[nk,1] est la normalisation du signal                    
00182       //       B2 = BI[nk,2] ewst le dephasage par rapport a une            
00183       //                 fonction centree en zero                                 
00184       //        Nous choisissons ici de demarrer avec les resultats           
00185       //            de l'ajustement parabolique mais il faudra bien             
00186       //            entendu verifier que cela ne biaise pas le resultat !      
00187       //            mise a zero des matrices utilisees dans la boucle             
00188       
00189       
00190       zero_mat(Z) ;
00191       zero_mat(YK) ;
00192       zero_mat(BK) ;
00193       zero_mat(C) ;
00194       zero_mat(D) ;
00195       
00196       nk=nevt ;  
00197       ampmax[nk] = 0. ;
00198       imax[nk] = 0 ;
00199       for ( k = 0 ; k < 10 ; k++) {     
00200         amplu[k]=adcval[nevt][k] ; 
00201         if (amplu[k] > ampmax[nk] ) {
00202           ampmax[nk] = amplu[k] ;
00203           imax[nk] = k ; 
00204         }
00205       }
00206       
00207       
00208       if( iter == 0 ) {
00209 
00210         // start with degree 3 polynomial .... 
00211         //fit3 =polfit(ns ,imax[nk] ,par3degre ,errpj ,amplu ) ;
00212         //      std::cout << "Poly Fit Param  :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl; 
00213         
00214         // start with parabol
00215         //fit3 = parab(amplu,4,12,par3degre) ;
00216         /*fit3 =*/ parab(amplu,2,9,par3degre) ;
00217         //std::cout << "Parab Fit Param :"<< par3degre[0] <<" "<< par3degre[1]<< std::endl; 
00218 
00219 
00220         // start with basic initial values
00221         //par3degre[0]= ampmax+ampmax/10. ;
00222         //par3degre[1]= (double)imax[nk]+0.1 ;
00223         //bi[nk][0] = ampmax[nk] ;
00224         //bi[nk][1] = (double)imax[nk] ;
00225 
00226         num_fit_min[nevt] = (double)imax[nk] - (double)nsmin ;
00227         num_fit_max[nevt] = (double)imax[nk] + (double)nsmax ;
00228         
00229         
00230         bi[nk][0] = par3degre[0] ;
00231         bi[nk][1] = par3degre[1] ;
00232         
00233         
00234         if (debug==1){
00235           printf("---------> depart ampmax[%d]=%f   maximum %f tim %f \n",
00236                  nk,ampmax[nk],bi[nk][0],bi[nk][1]);
00237         }
00238         
00239       } else {
00240         
00241         
00242         //  in other iterations  :                                              
00243         //   increment bi[][] parameters with bdi[][]                         
00244         //   calculated in previous                                           
00245         //   iteration                                                                      
00246         
00247 
00248         bi[nk][0] +=  dbi[nk][0] ;
00249         bi[nk][1] +=  dbi[nk][1] ;
00250         
00251         if (debug==1){
00252           printf("iter %d valeur de max %f et norma %f poly 3 \n",iter,bi[nk][1],bi[nk][0]) ;
00253         }       
00254       }
00255     
00256       b1 = bi[nk][0] ;
00257       b2 = bi[nk][1] ;
00258 
00259 
00261       // Loop on samples //                         
00263       
00264       chi2 = 0. ;
00265       ioktk[nk] = 0 ;
00266       ns =nborn_max-nborn_min+1 ;
00267        
00268       for (k=0 ; k < 10 ; k++){
00269         
00270         //
00271         //      calculation of fonction used to fit
00272         //       
00273         
00274         dt =(double)k - b2 ;
00275         t = (double)k ;
00276         amp = amplu[k] ;
00277         if (debug==1){
00278           printf(" CHECK sample %f ampli %f \n",t,amp) ;
00279         }
00280         //unsurs1 = 1./sig_err ;
00281         //unsurs2 = 1./(sig_err*sig_err) ;
00282         //unsurs1 = 0.1 ;
00283         //unsurs2 = 0.01 ;
00284         
00285         
00286         unsurs1=1./noise_val;
00287         //unsurs2=(1./noise_val)*(1./noise_val);
00288 
00289         //                                                   
00290         // Pulse shape function used: pulseShapepj
00291         //
00292         
00293         nborn_min = (int)num_fit_min[nevt] ;
00294         nborn_max = (int)num_fit_max[nevt] ;
00295         if(k < nborn_min || k > nborn_max ) continue ;
00296         tsig[0] =(double)k  ;
00297 
00298         
00299         if(METHODE==2) {
00300         par[0]=  b1 ;
00301         par[1] = b2 ;
00302         par[2] = a1 ;
00303         par[3] = a2 ;
00304         fun = pulseShapepj( tsig , par) ;
00305         }
00306         if (debug==1){
00307           printf(" valeur ampli %f et function %f min %d max %d \n",amp,fun,nsmin,nsmax) ;
00308           printf("min %f max %f \n",num_fit_min[nevt],num_fit_max[nevt]) ;
00309         }
00310         
00311         //       we need to determine a1,a2 which are global parameters         
00312         //        and  b1, b2 which are parameters for each individual signal: 
00313         //        b1, b2 = amplitude and time event by event
00314         //        a1, a2 = alpha and beta global      
00315         //        we first begin to calculate the derivatives used in the following calculation                                              
00316         
00317         if(METHODE==2){
00318           alpha = a1 ;
00319           beta  = a2 ;
00320           albet=alpha*beta;
00321           if(dt > -albet)  { 
00322             variab = (double)1. + dt/albet ;
00323             dtsbeta = dt/beta ;
00324             db1[k] = unsurs1*fun/b1 ;
00325             fact2 =  fun ;
00326             db2[k] = unsurs1*fact2*dtsbeta/(albet*variab) ;
00327             da1[k] = unsurs1*fact2*(log(variab)-dtsbeta/(alpha*variab)) ;
00328             da2[k] = unsurs1*fact2*dtsbeta*dtsbeta/(albet*variab) ;
00329           }
00330         }
00331         delta[k] = (amp - fun)*unsurs1 ;
00332         if (debug==1){
00333           printf(" ------->iter %d valeur de k %d amp %f fun %f delta %f \n",
00334                  iter,k,amp,fun,delta[k]) ;
00335           printf(" -----> valeur de k %d delta %f da1 %f da2 %f  \n",
00336                  k,delta[k],da1[k],da2[k]) ;
00337         }
00338 
00339         chi2 = chi2 + delta[k]*delta[k]     ;
00340         
00341         if (debug==1){
00342           printf(" CHECK chi2 %f deltachi2 %f sample %d iter %d \n",chi2,delta[k]*delta[k],k, iter) ;
00343         }
00344 
00345       }
00346     
00347       
00349       // End Loop on samples //                     
00351       
00352       
00353       double wk1wk2 ;
00354 
00356       // Start Loop on samples //                           
00358 
00359       for(int k1=nborn_min ; k1<nborn_max+1 ; k1++) {
00360         wk1wk2 = 1. ;
00361         int k2 = k1 ;
00362         
00363         DA.coeff[0][0] = da1[k1]*wk1wk2 ;
00364         DA.coeff[1][0] = da2[k1]*wk1wk2 ;
00365         DAT.coeff[0][0]= da1[k2] ;
00366         DAT.coeff[0][1]= da2[k2] ;
00367         DB.coeff[0][0] = db1[k1]*wk1wk2 ;
00368         DB.coeff[1][0] = db2[k1]*wk1wk2 ;
00369         DBT.coeff[0][0]= db1[k2] ;
00370         DBT.coeff[0][1]= db2[k2] ;
00371         
00372         //  Compute derivative matrix : matrix b[2][2]  
00373         
00374         produit_mat_int(DA,DAT,BK) ;
00375         
00376         //  Compute matrix c[2][2]                                  
00377         
00378         produit_mat_int(DA,DBT,C) ;
00379         
00380         //  Compute matrix d[2][2]                                   
00381         
00382         produit_mat_int(DB,DBT,D) ;
00383         
00384         //  Compute matrix y[3] and z[2] depending of delta (amp-fun)        
00385         
00386         delta2 = delta[k2] ;
00387         
00388         somme_mat_int_scale(DA,YK,delta2) ;                     
00389         somme_mat_int_scale(DB,Z,delta2) ;
00390         
00391         ioktk[nk]++ ;
00392         
00393       }
00394       
00396       // End Loop on samples //                     
00398       
00399       
00400       //  Remove events with a bad shape 
00401       
00402       if(ioktk[nk] < 4 ) {
00403         printf(" event rejected because npamp_used = %d \n",ioktk[nk]);
00404         continue ;
00405       }
00406       chi2s = chi2/(2. + (double)ns + 2.)  ;      
00407       chi2tot+=chi2s;
00408 
00409       if (debug==1){
00410         if (nevt==198 || nevt==199){
00411           std::cout << "adc123 pour l'evt " << nevt <<" = "<<adcval[nevt][nborn_min]<<" = "<<adcval[nevt][imax[nevt]]<<" = "<<adcval[nevt][nborn_max]<<std::endl;
00412           std::cout << "chi2s  pour l'evt " << nevt <<" = "<< chi2s<<" "<< chi2<<" "<< ns<<"  "<<iter<<std::endl;
00413           std::cout << "chi2tot           " << nevt <<" = "<< chi2tot<<"  "<<iter<<std::endl;
00414         }
00415       }
00416       
00417       //  Transpose matrix C ---> CT                        
00418      
00419       transpose_mat(C,CT) ;
00420       
00421       //  Calculate DM1 (inverse of D matrix 2x2)                 
00422       
00423       inverse_mat(D,DM1) ;
00424 
00425       
00426       //  Set matrix product c*d in memory in order to compute variations    
00427       //   of parameters B at the end of the iteration loop                   
00428       //   the variations of parameters b are dependant of the variations of  
00429       //   parameters da[1],da[2]                                            
00430          
00431       cti[nk][0] = CT.coeff[0][0]  ;
00432       cti[nk][1] = CT.coeff[0][1]  ;
00433       cti[nk][2] = CT.coeff[1][0]  ;
00434       cti[nk][3] = CT.coeff[1][1]  ;
00435       
00436       
00437        dm1i[nk][0] = DM1.coeff[0][0] ;
00438        dm1i[nk][1] = DM1.coeff[0][1] ;
00439        dm1i[nk][2] = DM1.coeff[1][0] ;
00440        dm1i[nk][3] = DM1.coeff[1][1] ; 
00441 
00442        zi[nk][0] = Z.coeff[0][0]  ;
00443        zi[nk][1] = Z.coeff[1][0]  ;
00444 
00445        //   Sum the matrix b and y after every event            
00446       
00447        for( k=0 ; k< numb_a ; k++) {
00448          Y.coeff[k][0] += YK.coeff[k][0] ;
00449        }
00450        somme_mat_int(BK,B) ;
00451        
00452 
00453        //   Calculate c(d-1)                                     
00454 
00455       produit_mat(C,DM1,CDM1) ;
00456 
00457       // Compute c(d-1)ct                                         
00458 
00459       produit_mat_int(CDM1,CT,CDM1CT);
00460                                                                  
00461       // Compute c(d-1)z                                          
00462       
00463       produit_mat_int(CDM1,Z,CDM1Z) ;
00464 
00465 
00466     }
00468       // End Loop on events //                      
00470     
00471     
00472     //  Compute b-cdm1ct
00473        
00474     diff_mat(B,CDM1CT,X) ;
00475     inverse_mat(X,XINV) ;
00476     diff_mat(Y,CDM1Z,RES2) ;
00477 
00478                                                                   
00479     // Calculation is now easy for da[0] da[1]                         
00480                                                                   
00481     produit_mat(XINV,RES2,A_CROISS) ;
00482 
00483     
00484     //  A la fin, on peut iterer en mesurant l'accroissement a apporter
00485     //    des parametres globaux par la formule db[i] = dm1(z-ct*da[i])  
00486          
00487     for( k=0 ; k< nk+1 ; k++) {
00488       
00489       if(METHODE ==2 ) {
00490         ZMCT.coeff[0][0] = zi[k][0] - (cti[k][0]*A_CROISS.coeff[0][0]+
00491                                        cti[k][1]*A_CROISS.coeff[1][0]) ;
00492         ZMCT.coeff[1][0] = zi[k][1] - (cti[k][2]*A_CROISS.coeff[0][0]+
00493                                        cti[k][3]*A_CROISS.coeff[1][0]) ;
00494       }
00495       
00496       dbi[k][0] = dm1i[k][0]*ZMCT.coeff[0][0] + dm1i[k][1]*ZMCT.coeff[1][0] ;
00497       dbi[k][1] = dm1i[k][2]*ZMCT.coeff[0][0] + dm1i[k][3]*ZMCT.coeff[1][0] ;
00498       if (debug==1){
00499         if( k < 100 ){
00500           printf(" variations de b1= %f et b2= %f  \n",dbi[k][0],dbi[k][1]) ;
00501         } 
00502       }
00503       db_i[k][0] = bi[k][0]+ dbi[k][0]   ;
00504       db_i[k][1] = bi[k][1]+ dbi[k][1]   ;
00505     }
00506     
00507     
00508     //   dbi[0] et dbi[1] mesurent les variations a apporter aux       
00509     //   parametres du signal                                          
00510                                                                      
00511     a1 += A_CROISS.coeff[0][0] ;
00512     a2 += A_CROISS.coeff[1][0] ;
00513     
00514 
00515     if (debug==1){
00516       printf(" CHECK croiss coef0: %f  croiss coef1: %f iter %d \n",fabs(A_CROISS.coeff[0][0]),fabs(A_CROISS.coeff[1][0]), iter);
00517     }
00518     if(fabs(A_CROISS.coeff[0][0]) < 0.001 && fabs(A_CROISS.coeff[1][0]) < 0.001)
00519       break;
00520     
00521   }
00522   
00524   //  End Loop on iterations //
00526 
00527   parout[0] = a1 ;
00528   parout[1] = a2 ;
00529   parout[2] = a3 ;
00530   if (debug==1){
00531     printf( " resultats trouves au bout de %d iterations \n",iter) ;
00532     printf( " parametre a1 = %f \n",a1) ;
00533     printf( " parametre a2 = %f \n",a2) ;
00534   }
00535 
00536   if (debug==1){
00537     std::cout << " Final chi2 / NDOF  :  "<< chi2tot/nevtmax << std::endl;
00538     std::cout << " Final (alpha,beta) : ("<< a1<<","<<a2<<")"<< std::endl;
00539   }
00540 
00541   return chi2tot/nevtmax ;
00542 
00543 }
00544 
00546 // End  Fitpj //
00548 
00549 
00550 
00551 /**************************************************************************/
00552 void TFParams::set_const( int n_samples, int sample_min, int sample_max ,
00553          double alpha ,double beta ,int nevtmaximum) {
00554 /*------------------------------------------------------------------------*/
00555   ns      = n_samples;
00556   nsmin   = sample_min ;
00557   nsmax   = sample_max ;
00558   nevtmax = nevtmaximum ;
00559   a1ini = alpha ;
00560   a2ini = 0.0 ;
00561   a3ini = beta ;
00562   step_shape = .04;
00563   METHODE = 2;
00564   if(ns > SDIM2) printf("warning: NbOfsamples exceed maximum\n");
00565 } 
00566 void TFParams::produit_mat(matrice A , matrice B , matrice M)
00567 {
00568   int i,j,k ;
00569 //  resultat du produit A*B = M 
00570   if(A.nb_colonnes != B.nb_lignes) {
00571     printf( " Erreur : produit de matrices de tailles incompatibles \n ");
00572     M.coeff = NULL ;
00573     return ;
00574   }
00575   M.nb_lignes = A.nb_lignes ;
00576   M.nb_colonnes = B.nb_colonnes ;
00577   zero_mat(M) ;
00578   for(i=0 ; i< M.nb_lignes; i++) {
00579     for(j=0 ; j< M.nb_colonnes ; j++) {
00580       for(k=0 ; k< A.nb_colonnes; k++){
00581         M.coeff[i][j] += A.coeff[i][k]*B.coeff[k][j] ;
00582       }
00583     }
00584   }
00585   return  ;
00586 }
00587 
00588 void TFParams::produit_mat_int(matrice A , matrice B, matrice M)
00589 {
00590   int i,j,k ;
00591   if(A.nb_colonnes != B.nb_lignes) {
00592     printf( " Erreur : produit de matrices de tailles incompatibles \n ");
00593     M.coeff = NULL ;
00594     return ;
00595   }
00596   M.nb_lignes = A.nb_lignes ;
00597   M.nb_colonnes = B.nb_colonnes ;
00598   for(i=0 ; i< M.nb_lignes; i++) {
00599     for(j=0 ; j< M.nb_colonnes ; j++) {
00600       for(k=0 ; k< A.nb_colonnes; k++){
00601         M.coeff[i][j] += A.coeff[i][k]*B.coeff[k][j] ;
00602       }
00603     }
00604   }
00605   return  ;
00606 }
00607 void TFParams::diff_mat(matrice A , matrice B , matrice M)
00608 {
00609   int i,j ;
00610 //resultat de la difference A-B = M 
00611   if(A.nb_lignes != B.nb_lignes) {
00612     printf( " Erreur : difference de matrices de tailles incompatibles \n ");
00613     M.coeff = NULL ;
00614     return ;
00615   }
00616   M.nb_lignes = A.nb_lignes ;
00617   M.nb_colonnes = A.nb_colonnes ;
00618   for(i=0 ; i< M.nb_lignes; i++) {
00619     for(j=0 ; j < M.nb_colonnes ; j++) {
00620       M.coeff[i][j] = A.coeff[i][j] - B.coeff[i][j] ;
00621     }
00622   }
00623   return  ;
00624   
00625 }
00626 void TFParams::copie_colonne_mat(matrice A , matrice M , int nk)
00627 {
00628   int i,j ;
00629   int k ;
00630  /* resultat de la copie de A dans un vecteur colonne M */
00631   k = 0 ;
00632   for(i=0 ; i< A.nb_lignes; i++) {
00633     for(j=0 ; j < A.nb_colonnes ; j++) {
00634       M.coeff[nk][k] = A.coeff[i][j] ;
00635    printf(" copie nk %d  i %d j %d k %d A %e M %e \n ",nk,i,j,k,A.coeff[i][j],
00636           M.coeff[nk][k]);      
00637       k++ ;
00638     }
00639   }
00640   return  ;
00641 }
00642 
00643 void TFParams::somme_mat_int(matrice A , matrice M)
00644 {
00645   int i,j;
00646  /* resultat de la somme integree M += A */
00647   if(A.nb_lignes != M.nb_lignes) {
00648     printf( " Erreur : somme de matrices de tailles incompatibles \n ");
00649     M.coeff = NULL ;
00650     return ;
00651   }
00652   M.nb_lignes = A.nb_lignes ;
00653   M.nb_colonnes = A.nb_colonnes ;
00654   for(i=0 ; i< M.nb_lignes; i++) {
00655     for(j=0 ; j< M.nb_colonnes ; j++) 
00656       M.coeff[i][j] += A.coeff[i][j] ;
00657   }
00658   return  ;
00659 }
00660 void TFParams::somme_mat_int_scale(matrice A , matrice M , double delta)
00661 {
00662   int i,j ;
00663   M.nb_lignes = A.nb_lignes ;
00664   M.nb_colonnes = A.nb_colonnes ;
00665   for(i=0 ; i< M.nb_lignes; i++) {
00666     for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[i][j] += A.coeff[i][j]*delta ;
00667     }
00668   return  ;
00669 }
00670 void TFParams::transpose_mat(matrice A , matrice M)
00671 {
00672   int i,j;
00673 // resultat de la transposition = matrice M 
00674   for(i=0 ; i< A.nb_lignes; i++) {
00675     for(j=0 ; j< A.nb_colonnes ; j++) {
00676         M.coeff[j][i] = A.coeff[i][j]  ;
00677     }
00678   }
00679   return  ;
00680 }
00681 matrice cree_mat_prod(matrice A , matrice B)
00682 {
00683   int i,j;
00684   matrice M ; /* resultat de la creation */
00685  
00686   M.nb_lignes = A.nb_lignes ;
00687   M.nb_colonnes = B.nb_colonnes ;
00688   M.coeff = (double**)malloc(M.nb_lignes*sizeof(double*)) ;
00689   for(i=0 ; i< M.nb_lignes; i++) 
00690     M.coeff[i]=(double*)calloc(M.nb_colonnes,sizeof(double));
00691   for(i=0 ; i< M.nb_lignes; i++) {
00692 
00693     for(j=0 ; j< M.nb_colonnes ; j++) {
00694         M.coeff[i][j] = 0.  ;
00695     }
00696   }
00697   //printf(" creation de matrice ---->  nlignes %d ncolonnes %d  \n",
00698 //       M.nb_lignes,M.nb_colonnes) ;
00699   return (M) ;
00700 }
00701 matrice cree_mat(int n_lignes,int n_colonnes)
00702 {
00703   int i,j;
00704   matrice M ; /* resultat de la creation */
00705  
00706   M.nb_lignes = n_lignes ;
00707   M.nb_colonnes = n_colonnes ;
00708   M.coeff = (double**)malloc(M.nb_lignes*sizeof(double*)) ;
00709   for(i=0 ; i< M.nb_lignes; i++) 
00710     M.coeff[i]=(double*)calloc(M.nb_colonnes,sizeof(double));
00711   for(i=0 ; i< M.nb_lignes; i++) {
00712     for(j=0 ; j< M.nb_colonnes ; j++) {
00713         M.coeff[i][j] = 0.  ;
00714     }
00715   }
00716   //printf(" creation de matrice --->  nlignes %d ncolonnes %d  \n",
00717         // M.nb_lignes,M.nb_colonnes) ;
00718   return (M) ;
00719 }
00720 
00721 void fill_mat( matrice A , matrice M)
00722 {
00723   int i,j;
00724   /* on remplit la matrice M avec la matrice A */
00725  
00726   M.nb_lignes = A.nb_lignes ;
00727   M.nb_colonnes = A.nb_colonnes ;
00728   for(i=0 ; i< M.nb_lignes; i++) {
00729     for(j=0 ; j< M.nb_colonnes ; j++) {
00730         M.coeff[i][j] = A.coeff[i][j]  ;
00731         printf("matrice remplie %e \n",M.coeff[i][j]) ;
00732     }
00733   } 
00734   return ;
00735 }
00736 void TFParams::print_mat(matrice M)
00737 {
00738   int i,j ;
00739   if( M.coeff == NULL) 
00740     printf(" erreur : affichage d'une matrice vide \n") ;
00741   printf(" m_nli %d M_ncol %d \n",M.nb_lignes,M.nb_colonnes) ;
00742   for(i=0 ; i< M.nb_lignes; i++) {
00743     for(j=0 ; j< M.nb_colonnes ; j++) 
00744       printf(" MATRICE i= %d j= %d ---> %e \n",i,j,M.coeff[i][j]) ;
00745   }
00746   //printf(" apres passage d'impression \n") ;
00747   return ;
00748 }
00749 void TFParams::zero_mat(matrice M)
00750 {
00751   int i,j ;
00752   for(i=0 ; i< M.nb_lignes; i++) {
00753     for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[i][j]=0. ; 
00754   }
00755   return ;
00756 }
00757 void TFParams::zero_mat_nk(matrice M,int nk)
00758 {
00759   int j ;
00760     for(j=0 ; j< M.nb_colonnes ; j++) M.coeff[nk][j]=0. ;
00761   return ;
00762 }
00763 void TFParams::print_mat_nk(matrice M,int nk)
00764 {
00765   int j ;
00766   if( M.coeff == NULL)
00767     printf(" erreur : affichage d'une matrice vide \n") ;
00768   printf(" nk = %d m_nli %d M_ncol %d \n",nk,M.nb_lignes,M.nb_colonnes) ;
00769     for(j=0 ; j< M.nb_colonnes ; j++) 
00770       printf(" MATRICE nk= %d j= %d  ---> %e \n",nk,j,M.coeff[nk][j]) ;    
00771   printf(" apres passage d'impression \n") ;
00772   return ;
00773 }
00774 void TFParams::inverse_mat( matrice A , matrice M )
00775 {
00776 /*   A[ligne][colonne]   B[ligne][colonne]   */
00777 int i , j   ;
00778 double  deter=0.  ;
00779 /*  M est la matrice inverse de A */
00780  
00781  if(A.nb_lignes != A.nb_colonnes) {
00782    printf( " attention matrice non inversible !!!! %d lignes %d colonnes \n",
00783            A.nb_lignes,A.nb_colonnes) ;
00784    return ;
00785  }
00786   zero_mat(M) ;
00787   if(A.nb_lignes == 2) {
00788     deter = A.coeff[0][0]*A.coeff[1][1] - A.coeff[0][1]*A.coeff[1][0] ;
00789     M.coeff[0][0] = A.coeff[1][1]/deter ;
00790     M.coeff[0][1] = -A.coeff[0][1]/deter ;
00791     M.coeff[1][0] = -A.coeff[1][0]/deter ;
00792     M.coeff[1][1] = A.coeff[0][0]/deter ;
00793   }
00794  else if(A.nb_lignes == 3) {
00795       M.coeff[0][0]=A.coeff[1][1]*A.coeff[2][2]-A.coeff[2][1]*A.coeff[1][2] ;
00796       M.coeff[1][1]=A.coeff[0][0]*A.coeff[2][2]-A.coeff[2][0]*A.coeff[0][2] ;
00797 
00798       M.coeff[2][2]=A.coeff[0][0]*A.coeff[1][1]-A.coeff[0][1]*A.coeff[1][0] ;
00799       M.coeff[0][1]=A.coeff[2][1]*A.coeff[0][2]-A.coeff[0][1]*A.coeff[2][2] ;
00800       M.coeff[0][2]=A.coeff[0][1]*A.coeff[1][2]-A.coeff[1][1]*A.coeff[0][2] ;
00801       M.coeff[1][0]=A.coeff[1][2]*A.coeff[2][0]-A.coeff[1][0]*A.coeff[2][2] ;
00802       M.coeff[1][2]=A.coeff[1][0]*A.coeff[0][2]-A.coeff[0][0]*A.coeff[1][2] ;
00803       M.coeff[2][0]=A.coeff[1][0]*A.coeff[2][1]-A.coeff[1][1]*A.coeff[2][0] ;
00804       M.coeff[2][1]=A.coeff[0][1]*A.coeff[2][0]-A.coeff[0][0]*A.coeff[2][1] ;
00805       deter=A.coeff[0][0]*M.coeff[0][0]+A.coeff[1][0]*M.coeff[0][1]
00806         +A.coeff[2][0]*M.coeff[0][2] ;
00807       for ( i=0 ; i<3 ; i++ ) {
00808         for ( j=0 ; j<3 ; j++ ) M.coeff[i][j] = M.coeff[i][j]/deter  ;
00809       }
00810  }
00811  else {
00812  printf(" Attention , on ne peut inverser la MATRICE %d \n",A.nb_lignes) ;
00813  return ;
00814  }
00815   
00816  return ;
00817 }
00818 Double_t TFParams::polfit(Int_t ns ,Int_t imax , Double_t par3d[dimout] , 
00819                  Double_t errpj[dimmat][dimmat] ,double *adcpj )
00820  {
00821   double val , val2 , val3 , adfmx[dimmat] , parfp3[dimout]  ;
00822   double ius[dimmat], maskp3[dimmat] ;
00823   double deglib,fit3,tm,h,xki2 ;
00824   int i ,nus ,ilow,isup ;
00825   val=adcpj[imax] ;
00826   val2=val/2.  ;
00827   val3=val/3.  ;
00828   ilow=0       ;
00829   isup=ns    ;
00830   deglib=-4.  ;
00831   for (i=0 ; i<ns ; i++ ){
00832     deglib=deglib+1.  ;
00833     ius[i] = 1. ;
00834     if((adcpj[i] < val3) && (i<imax) ){
00835       ilow=i  ;
00836     }
00837     if(adcpj[i] > val2 ){
00838       isup=i  ;
00839     }
00840   }
00841   ilow=ilow+1   ;
00842   if(ilow == imax )ilow=ilow-1 ;
00843   if(isup-ilow < 3) isup=ilow+3 ;
00844   nus=0  ;
00845   for(i=ilow ; i<=isup ; i++){
00846     
00847     adfmx[nus]=adcpj[i]  ;
00848     maskp3[nus] =0. ;
00849     if(ius[i] == 1) {
00850       maskp3[nus]=1. ;
00851       nus=nus+1    ;
00852     }
00853   }
00854   if(nus < 4) return 10000. ;
00855   xki2 =  f3deg (  nus , parfp3 ,  maskp3 , adfmx ,  errpj ) ;
00856   tm= parfp3[4]  ;
00857   h=parfp3[5] ;
00858   tm=tm+(double)ilow  ;
00859   par3d[0] = h ;
00860   par3d[1] = tm ;
00861   fit3 = xki2 ;
00862   return fit3 ;
00863 }
00864 double TFParams::f3deg (  int nmxu ,  double parom[dimout] , double mask[dimmat] , double adcpj[dimmat] , double errpj[dimmat][dimmat] ) {
00865 /*                                                                   */
00866 /*  fit   3rd degree polynomial                                      */
00867 /*  nmxu = nb of samples in sample data array adcpj[]
00868     parom   values of parameters
00869     errpj  inverse of the error matrix
00870     fplo3dg uses only the diagonal terms of errpj[][]
00871 */
00872   int i , k , l  /*,iworst*/ ;
00873   double  h , t2 , tm , delta , tmp ;
00874   double xki2 , dif , difmx , deglib   ;
00875   double t[dimmat] ,  f[dimmat][4]   ;
00876   double cov[dimmat][dimmat] , bv[4] , invcov[dimmat][dimmat] , s /*, deter*/  ;
00877   
00878   deglib=(double)nmxu - 4.  ;
00879   for ( i=0 ; i<nmxu ; i++ ) {
00880     t[i]=i ;
00881     f[i][0]=1. ;
00882     f[i][1]=t[i]  ;
00883     f[i][2]=t[i]*t[i]  ;
00884     f[i][3]=f[i][2]*t[i] ;
00885   }
00886 /*   computation of covariance matrix     */
00887   for ( k=0 ; k<4 ; k++ ) {
00888     for ( l=0 ; l<4 ; l++ ) {
00889       s=0.   ;
00890       for (i=0 ; i<nmxu ; i++ ) {
00891         s=s+f[i][k]*f[i][l]*errpj[i][i]*mask[i]   ;
00892       }
00893       cov[k][l]=s  ;
00894     }
00895     s=0.    ;
00896     for (i=0 ; i<nmxu ; i++ ) {
00897         s=s+f[i][k]*adcpj[i]*errpj[i][i]*mask[i]   ;
00898     }
00899       bv[k]=s  ;
00900   }
00901 /*     parameters                          */
00902   /*deter =*/ inverpj ( 4 , cov , invcov );
00903   for ( k=0 ; k<4 ; k++ ) {
00904     s=0.  ;
00905     for ( l=0 ; l<4 ; l++ ) {
00906       s=s+bv[l]*invcov[l][k]   ;
00907     }
00908     parom[k]=s  ;
00909   }
00910 
00911   if( parom[3] == 0. ){
00912     parom[4] = -1000.;
00913     parom[5] = -1000.;
00914     parom[6] = -1000.;
00915     return 1000000.;
00916   }
00917 /*    worst hit and ki2                    */
00918   xki2=0.    ;
00919   difmx=0.   ;
00920     for (i=0 ; i<nmxu ; i++ ){
00921       t2=t[i]*t[i]  ;
00922       h= parom[0]+parom[1]*t[i]+parom[2]*t2+parom[3]*t2*t[i] ;
00923       dif=(adcpj[i]-h)*mask[i]     ;
00924         if(dif > difmx) {
00925           // iworst=i  ;
00926           difmx=dif ;
00927         }
00928     }
00929     if(deglib > 0.5) xki2=xki2/deglib ;
00930 /*     amplitude and maximum position                    */
00931   delta=parom[2]*parom[2]-3.*parom[3]*parom[1]  ;
00932   if(delta > 0.){
00933     delta=sqrt(delta)  ;
00934     tm=-(delta+parom[2])/(3.*parom[3])  ;
00935     tmp=(delta-parom[2])/(3.*parom[3])  ;
00936   }
00937   else{
00938     parom[4] = -1000.;
00939     parom[5] = -1000.;
00940     parom[6] = -1000.;
00941     return xki2  ;
00942   }
00943   parom[4]= tm  ;
00944   parom[5]= parom[0]+parom[1]*tm+parom[2]*tm*tm+parom[3]*tm*tm*tm ;
00945   parom[6]= tmp ;
00946   // printf("par --------> %f %f %f %f \n",parom[3],parom[2],parom[1],parom[0]);
00947   
00948    return xki2  ;
00949 }
00950 /*------------------------------------------------------------------*/
00951 
00952 double TFParams::inverpj(int n,double g[dimmat][dimmat],double ginv[dimmat][dimmat] )
00953 {
00954 /*                                                                   */
00955 /*  inversion d une matrice symetrique definie positive de taille n  */
00956 /*  J.P. Pansart   Novembre 99                                       */
00957 /*                                                                   */
00958 int i , j , k , jj  ;
00959 double r ,  s  ;
00960 double deter=0  ;
00961 double al[dimmat][dimmat] , be[dimmat][dimmat]  ;
00962 /*   initialisation                                                  */
00963  for( i=0 ; i<n ; i++ ) {
00964    for ( j=0 ; j<n ; j++ ) {
00965     al[i][j] = 0.  ;
00966     be[i][j] = 0.  ;
00967    }
00968  }
00969 /*  decomposition en vecteurs sur une base orthonormee               */
00970  al[0][0] =  sqrt( g[0][0] )  ;
00971  for ( i=1 ; i<n ; i++ ) {
00972  al[i][0] = g[0][i] / al[0][0]  ;
00973    for ( j=1 ; j<=i ; j++ ) {
00974     s=0.   ;
00975     for ( k=0 ; k<=j-1 ; k++ ) {
00976      s=s+ al[i][k] * al[j][k]  ;
00977     }
00978     r= g[i][j] - s   ;
00979    if( j < i ) al[i][j] = r/al[j][j]  ;
00980    if( j == i ) al[i][j] =  sqrt ( r)  ;
00981    }
00982  }
00983 /*  inversion de la matrice al                                       */
00984  be[0][0] = 1./al[0][0]  ;
00985  for ( i=1 ; i<n ; i++ ) {
00986  be[i][i] = 1. / al[i][i]  ;
00987    for ( j=0 ; j<i ; j++ ) {
00988     jj=i-j-1  ;
00989     s=0.   ;
00990     for ( k=jj+1 ; k<=i ; k++ ) {
00991      s=s+ be[i][k] * al[k][jj]  ;
00992     }
00993     be[i][jj]=-s/al[jj][jj]  ;
00994    }
00995  }
00996 /*   calcul de la matrice ginv                                       */
00997  for ( i=0 ; i<n ; i++ ) {
00998    for ( j=0 ; j<n ; j++ ) {
00999     s=0.   ;
01000     for ( k=0 ; k<n ; k++ ) {
01001      s=s+ be[k][i] * be[k][j]  ;
01002     }
01003     ginv[i][j]=s  ;
01004     //    if (debug==1){
01005     //printf("valeur de la matrice %d %d %f \n",i,j,ginv[i][j]) ;
01006     //}
01007    }
01008  }
01009   return deter ;
01010 }
01011 /*                                                                   */
01012 /*  inversion d une matrice 3x3                                      */
01013 /*                                                                   */
01014 double TFParams::inv3x3(double a[3][3] , double b[3][3] )
01015 {
01016 /*   a[ligne][colonne]   b[ligne][colonne]   */
01017 int i , j   ;
01018 double  deter=0.  ;
01019       b[0][0]=a[1][1]*a[2][2]-a[2][1]*a[1][2] ;
01020       b[1][1]=a[0][0]*a[2][2]-a[2][0]*a[0][2] ;
01021       b[2][2]=a[0][0]*a[1][1]-a[0][1]*a[1][0] ;
01022       printf("a[x][x] %e %e %e %e %e %e %e \n",a[0][0],a[1][1],a[0][1],a[1][0],
01023              a[0][0]*a[1][1],a[0][1]*a[1][0],b[2][2]);
01024       b[0][1]=a[2][1]*a[0][2]-a[0][1]*a[2][2] ;
01025       b[0][2]=a[0][1]*a[1][2]-a[1][1]*a[0][2] ;
01026       b[1][0]=a[1][2]*a[2][0]-a[1][0]*a[2][2] ;
01027       b[1][2]=a[1][0]*a[0][2]-a[0][0]*a[1][2] ;
01028       b[2][0]=a[1][0]*a[2][1]-a[1][1]*a[2][0] ;
01029       b[2][1]=a[0][1]*a[2][0]-a[0][0]*a[2][1] ;
01030       deter=a[0][0]*b[0][0] + a[1][0]*b[0][1] + a[2][0]*b[0][2] ;
01031       printf(" deter = %e \n",deter) ;
01032  for ( i=0 ; i<3 ; i++ ) {
01033    for ( j=0 ; j<3 ; j++ ) {
01034      printf(" avant division a[3][3] %d %d  %e \n",i,j,a[i][j]) ;
01035      printf(" avant division b[3][3] %d %d  %e %e \n",i,j,b[i][j],deter) ;
01036      b[i][j] = b[i][j]/deter  ;
01037      printf(" valeur de b[3][3] apres division %d %d  %e %e \n",i,j,b[i][j],deter) ;
01038    }
01039  }
01040   return deter ;
01041 }
01042 
01043 double TFParams::pulseShapepj( Double_t *x, Double_t *par )
01044 {
01045 
01046   Double_t fitfun;
01047   Double_t ped, h, tm, alpha, beta;
01048   Double_t  dt, dtsbeta, albet, variab, puiss;
01049   Double_t b1,b2,a1,a2 ;
01050   b1 = par[0] ;
01051   b2 = par[1] ;
01052   a1 = par[2] ;
01053   a2 = par[3] ;
01054 
01055   ped   =  0. ;
01056   h     =  b1 ;
01057   tm    =  b2 ;
01058   alpha =  a1 ;
01059   beta  =  a2 ;
01060   dt= x[0] - tm  ;
01061   //printf(" par %f %f %f %f dt = %f albet = %f",b1,b2,a1,a2,dt,albet) ;
01062   albet = alpha*beta ;
01063   if( albet <= 0 )return( (Double_t)0. );
01064 
01065   if(dt > -albet)  {
01066     dtsbeta=dt/beta ;
01067     variab=1.+dt/albet ;
01068     puiss=pow(variab,alpha);
01069     fitfun=h*puiss*exp(-dtsbeta) + ped;
01070     //printf(" dt = %f h = %f puiss = %f exp(-dtsbeta) %f \n",dt,h,puiss,
01071     // exp(-dtsbeta)) ;
01072      }
01073   else {
01074       fitfun = ped;
01075      }
01076 
01077   return fitfun;
01078 }
01079 
01080 double TFParams::lastShape( Double_t *x, Double_t *par )
01081 {
01082 
01083   Double_t fitfun;
01084   Double_t alpha, beta;
01085   Double_t dt,alphadt,exponent ;
01086   Double_t b1,b2 ;
01087   b1 = par[0] ;
01088   b2 = par[1] ;
01089   alpha = par[2] ;
01090   beta  = par[3] ;
01091   dt= x[0] - b2  ;
01092   alphadt = alpha*dt ;
01093   exponent = -(alphadt+(exp(-alphadt)-1.))/beta ; 
01094   fitfun = b1*exp(exponent) ; 
01095   return fitfun;
01096 }
01097 double TFParams::lastShape2( Double_t *x, Double_t *par )
01098 {
01099 
01100   Double_t fitfun;
01101   Double_t alpha, beta;
01102   Double_t dt,expo1,dt2,exponent ;
01103   Double_t b1,b2 ;
01104   b1 = par[0] ;
01105   b2 = par[1] ;
01106   alpha = par[2] ;
01107   beta  = par[3] ;
01108   dt= x[0] - b2  ;
01109   expo1 = exp(-beta*dt) ;
01110   dt2 = dt*dt ;
01111   exponent = -(alpha*dt2+(expo1-1.)) ;
01112   fitfun = b1*exp(exponent) ; 
01113   return fitfun;
01114 }
01115 
01116 Double_t TFParams::pulseShapepj2( Double_t *x, Double_t *par )
01117 {
01118 
01119   Double_t fitfun;
01120   Double_t ped, h, /*tm,*/ alpha, beta;
01121   Double_t  dt, dtsbeta, albet, variab, puiss;
01122   Double_t b1,/*b2,*/a1,a2 ;
01123   b1 = par[0] ;
01124   //b2 = par[1] ;
01125   a1 = par[2] ;
01126   a2 = par[3] ;
01127   ped   =  0. ;
01128   h     =  b1 ;
01129   //tm    =  b2 ;
01130   alpha =  a1 ;
01131   beta  =  a2 ;
01132   dt= x[0]  ;
01133   albet = alpha*beta ;
01134   if( albet <= 0 )return( (Double_t)0. );
01135 
01136   if(dt > -albet)  {
01137     dtsbeta=dt/beta ;
01138     variab=1.+dt/albet ;
01139     puiss=pow(variab,alpha);
01140     fitfun=h*puiss*exp(-dtsbeta) + ped;
01141   }
01142   else {
01143     fitfun = ped;
01144   }
01145 
01146   /*  printf( "fitfun %f %f %f %f, %f %f %f\n", ped, h, tm, alpha, beta, *x, fitfun );  */
01147 
01148   return fitfun;
01149 }
01150 
01151 double TFParams::parab(Double_t ampl[nsamp],Int_t nmin,Int_t nmax,Double_t parout[3])
01152 {
01153 /* Now we calculate the parabolic adjustement in order to get        */
01154 /*    maximum and time max                                           */
01155   
01156   double denom,dt,amp1,amp2,amp3 ; 
01157   double ampmax = 0. ;                          
01158   int imax = 0 ;
01159   int k ;
01160 /*
01161                                                                    */     
01162   for ( k = nmin ; k < nmax ; k++) {
01163     if (ampl[k] > ampmax ) {
01164       ampmax = ampl[k] ;
01165       imax = k ;
01166     }
01167   }
01168         amp1 = ampl[imax-1] ;
01169         amp2 = ampl[imax] ;
01170         amp3 = ampl[imax+1] ;
01171         denom=2.*amp2-amp1-amp3  ;
01172 /*                                                                   */       
01173         if(denom>0.){
01174           dt =0.5*(amp3-amp1)/denom  ;
01175         }
01176         else {
01177           //printf("denom =%f\n",denom)  ;
01178           dt=0.5  ;
01179         }
01180 /*                                                                   */        
01181 /* ampmax correspond au maximum d'amplitude parabolique et dt        */
01182 /* decalage en temps par rapport au sample maximum soit k + dt       */
01183                 
01184         parout[0] =amp2+(amp3-amp1)*dt*0.25 ;
01185         parout[1] = (double)imax + dt ;
01186         parout[2] = (double)imax ;
01187 return denom ;
01188 }
01189 
01190 double TFParams::mixShape( Double_t *x, Double_t *par )
01191 {
01192   Double_t fitval0,fitval;
01193   Double_t alpha,beta,fact,puiss;
01194   Double_t dt,alpha2dt,exponent ;
01195   Double_t b1,b2,alpha2,t ;
01196   b1 = par[0] ;
01197   b2 = par[1] ;
01198   alpha  = par[2] ;
01199   alpha2 = par[3] ;
01200   beta   = par[4] ;
01201 //
01202   t = x[0] ;
01203   dt= x[0]-b2  ;
01204 //
01205   if(t>0.) {
01206   fact = t/b2 ;
01207   puiss = pow(fact,alpha) ;
01208   fitval0 = puiss*exp(-alpha*dt/b2) ;
01209   }
01210   else
01211   {
01212   fitval0=1. ;
01213   }
01214   dt = x[0] - b2 ;
01215   alpha2dt = dt*alpha2 ;
01216   exponent = -(alpha2dt+(exp(-alpha2dt)-1.))/beta ;
01217   fitval = b1*fitval0*exp(exponent) ;  
01218   return fitval;
01219 }
01220 //=========================
01221 // Method computePulseWidth
01222 //=========================
01223 double TFParams::computePulseWidth( int methode , double alpha_here , double beta_here){ 
01224 
01225 // level of amplitude where we calculate the width ( level = 0.5 if at 50 % )
01226 //   (level = 0.3 if at 30 % )
01227   double level = 0.30 ;
01228 // fixed parameters
01229   double amplitude   = 1.00 ;
01230   double offset      = 7.00;  
01231   double amp_max     = amplitude;
01232 
01233 // steps in time
01234   double t_min       =  offset-4.50;
01235   double t_max       =  offset+12.50;
01236 
01237   int    t_step_max  = 3000;
01238   double delta_t     =  (double)((t_max-t_min)/t_step_max);
01239       
01240 // Loop over time ( Loop 2  --> get width )
01241   int    t_amp_half_flag =    0;
01242   double t_amp_half_min  =  999.; 
01243   double t_amp_half_max  = -999.; 
01244 
01245   for (int t_step=0; t_step<t_step_max; t_step++){
01246 
01247        double t_val = t_min + (double)t_step*delta_t;
01248        double albet = alpha_here*beta_here ;
01249        double dt = t_val-offset ;
01250        double amp =0;
01251 
01252        if( methode == 2 ) { // electronic function
01253           if( (t_val-offset) > -albet) {
01254 
01255             amp =  amplitude*TMath::Power( ( 1 + ( dt / (alpha_here*beta_here) ) ) , alpha_here ) * TMath::Exp(-1.0*(dt/beta_here));
01256           } else {
01257             
01258             amp = 1. ;
01259           }
01260        } 
01261 
01262        if( amp > (amp_max*level) && t_amp_half_flag == 0) {
01263            t_amp_half_flag = 1;
01264            t_amp_half_min = t_val;
01265        }
01266 
01267        if( amp < (amp_max*level) && t_amp_half_flag == 1) {
01268            t_amp_half_flag = 2;
01269            t_amp_half_max = t_val;
01270        }          
01271 
01272    }
01273     
01274 // Compute Width
01275   double width = (t_amp_half_max - t_amp_half_min);
01276 
01277   return width;
01278 }
01279 
01280