////////////////////////////////////
// KRATTA Event analyser
// for the Asy-Eos experiment
// TKratAnaStep1 DECLARATION
// Nov 2012
// revison 11/2012
// E.d.F ver 1.0 
// sebastian.kupny@uj.edu.pl
// jerzy.lukasik@ifj.edu.pl
// Changes:
////////////////////////////////////



#include <string>
#include <iostream>
    


#include "TClonesArray.h"
#include "TMath.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TGeoManager.h"

#include "TKratAnaStep1.h"


using std::cout;
using std::endl;


//______________________________________________________________________
void FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{

  Double_t chi2=0.;

  Double_t cu,eu,fu,fsum;
  Double_t *zik=0;
  Double_t *pl0=0;
  Int_t i, bin,binx,biny,binz;
  Axis_t binlow, binup, binsize;

  TFumili *fitter = (TFumili*)TVirtualFitter::GetFitter();

  TH1D *hist = (TH1D*)fitter->GetObjectFit();
  TF1 *f1 = (TF1*)fitter->GetUserFunc();
  int nbx = hist->GetNbinsX();
  //printf("(1) %f %f %f %f\n",par[0],par[1],par[2]);
  Double_t x[3]={0.,0.,0.};
  npar = f1->GetNpar();
  fitter->SetParNumber(npar);
  if(iflag == 9) return;
  zik = fitter->GetZ();
  pl0 = fitter->GetPL0();

//  Double_t *df=new Double_t[npar];
Double_t df[12]={0.};
  f1->InitArgs(x,par);
  f = 0;
  //printf("(2) %f %f %f %f\n",par[0],par[1],par[2]);

  // printf("%s\n",hist->GetName());
  // printf("npar %d %d\n",npar,iflag);
  //   printf("%7.4f %8.3f %11.3f %11.3f %e\n",par[0],par[1],par[2],par[3],par[4]);

  int np=0;
  for(int ix=1;ix<=nbx;ix++){
    x[0]=(Double_t)hist->GetXaxis()->GetBinCenter(ix);
    cu=(Double_t)hist->GetBinContent(ix);
    eu=(Double_t)hist->GetBinError(ix);
    //  eu=1;
    //  eu=sqrt(fabs(cu));
    //   if(cu>0 && eu>0){
    if(eu>0){
    fu=f1->EvalPar(x,par);
    //     if(fu==-1.){
    //       f=1.e9;
    //       return;
    //     }
    np++;
    fitter->Derivatives(df,x);
    Int_t N = 0;
    fsum = (fu-cu)/eu;
    if(iflag!=1){
      for(i=0;i<npar;i++) 
      if(pl0[i]>0){
      df[N] = df[i]/eu; 
      // left only non-fixed param derivatives / by Sigma
      gin[i] += df[N]*fsum;
      N++;
      }
      Int_t L = 0;
      for (i=0;i<N;i++)
      for (Int_t j=0;j<=i;j++) 
      zik[L++] += df[i]*df[j];
    }
    //     f += .5*fsum*fsum;
    f += fsum*fsum;
    }
  }

  f1->SetNumberFitPoints(np);
//  delete[] df;

  chi2=f;

  // take care of the fixed parameters...
  Int_t  ndof=np-1-(npar-1);
  if(iflag==3){
  //printf("%c\n",7);
  //printf("%c\n",7);

  //   double pr_chi2 = 1.-TMath::Gamma((Double_t)(0.5*ndof),(Double_t)(0.5*chi2));  
  //   double ch2 = (ndof)? chi2/ndof : chi2;
  //   pr_chi2 = 1.-TMath::Gamma((Double_t)(0.5*ndof),(Double_t)(0.5*chi2));  
  //   ch2 = (ndof>0)? chi2/ndof : chi2;
  //printf("chi2=%f chi2/ndf=%f pr_chi2=%f\n",chi2,ch2,pr_chi2);

  //   printf("%7.1f %5.3f %5.3f %7.1f %e\n",par[0],par[1],par[2],par[3],par[4]);
  //   printf("%7.1f %5.3f %5.3f %7.1f %e\n",
  //   fitter->fParamError[0],
  //   fitter->fParamError[1],
  //   fitter->fParamError[2],
  //   fitter->fParamError[3],
  //   fitter->fParamError[4]
  //   );

  // printf("%8.3f %8.3f  %7.3f %7.3f  %7.3f %7.3f  %8.3f %8.3f  %7.3f\n",
  // par[0],fitter->fParamError[0],
  // par[1],fitter->fParamError[1],
  // par[2],fitter->fParamError[2],
  // par[3],fitter->fParamError[3],ch2);

  //  fitter->PrintResults(4,f);

  //   for(int i=0;i<5;i++){
  //     fitter->GetParameter(i, &name[0], value, verr, vlow, vhigh);
  //     fitter->GetErrors(i, eplus, eminus, eparab, globcc);
  //     ppar[i]=value;
  //     epar[i]=eparab;
  //     printf("%d %s %f %f %f %f %f\n",i,name,value,eplus, eminus, eparab, globcc);
  //   } 
  }
  //printf("chi2/ndf=%f\n",chi2/ndof);

}


//______________________________________________________________________
double RC4DEXP2TR(Double_t *xx, Double_t *par){

  //rozne czasy narastania dla si i csi

  double x=xx[0];// 

  double Q0     = par[0]; // 
  double Q1     = par[1]; // 
  double Q2     = par[2]; // 
  double Q3     = par[3]; // 
  double x0     = par[4]; // 
  double t0     = par[5]; // decay si0
  double t1     = par[6]; // decay si1
  double t2     = par[7]; // decay csi0
  double t3     = par[8]; // decay csi1
  double tr1    = par[9]; // common si rise time
  double tr2    = par[10]; // common csi rise time
  double rc     = par[11]; // 

  // !!!!  Q0 is in fact Q0/C

  double fun = 0.; // 
  double exp_rc = exp(-(x-x0)/rc); // 
  if(x<=x0) fun=0;
  else {
  if(t0!=tr1){
    fun+=Q0*rc*(exp_rc*rc/((rc-t0)*(rc-tr1))+exp(-(x-x0)/t0)*t0/((t0-rc)*(t0-tr1))+exp(-(x-x0)/tr1)*tr1/((tr1-rc)*(tr1-t0)));
  }
  else{
    fun+=Q0*rc*(exp_rc*rc*t0-exp(-(x-x0)/t0)*(rc*(x-x0+t0)-(x-x0)*t0))/(t0*pow(rc-t0,2));
  }
  if(t1!=tr1){
    fun+=Q1*rc*(exp_rc*rc/((rc-t1)*(rc-tr1))+exp(-(x-x0)/t1)*t1/((t1-rc)*(t1-tr1))+exp(-(x-x0)/tr1)*tr1/((tr1-rc)*(tr1-t1)));
  }
  else{
    fun+=Q1*rc*(exp_rc*rc*t1-exp(-(x-x0)/t1)*(rc*(x-x0+t1)-(x-x0)*t1))/(t1*pow(rc-t1,2));
  }
  if(t2!=tr2){
    fun+=Q2*rc*(exp_rc*rc/((rc-t2)*(rc-tr2))+exp(-(x-x0)/t2)*t2/((t2-rc)*(t2-tr2))+exp(-(x-x0)/tr2)*tr2/((tr2-rc)*(tr2-t2)));
  }
  else{
    fun+=Q2*rc*(exp_rc*rc*t2-exp(-(x-x0)/t2)*(rc*(x-x0+t2)-(x-x0)*t2))/(t2*pow(rc-t2,2));
  }
  if(t3!=tr2){
    fun+=Q3*rc*(exp_rc*rc/((rc-t3)*(rc-tr2))+exp(-(x-x0)/t3)*t3/((t3-rc)*(t3-tr2))+exp(-(x-x0)/tr2)*tr2/((tr2-rc)*(tr2-t3)));
  }
  else{
    fun+=Q3*rc*(exp_rc*rc*t3-exp(-(x-x0)/t3)*(rc*(x-x0+t3)-(x-x0)*t3))/(t3*pow(rc-t3,2));
  }
  }

  return fun;
}

//______________________________________________________________________

double RC4DEXP2TR1(Double_t *xx, Double_t *par){

  double x=xx[0];// 

  double Q0     = par[0]; // 
  double x0     = par[1]; // 
  double t0     = par[2]; // 
  double tr     = par[3]; // 
  double rc     = par[4]; // 

  // !!!!  Q0 is in fact Q0/C

  double fun = 0.; // 
  double exp_rc = exp(-(x-x0)/rc); // 
  if(x<=x0) fun=0;
  else {
    if(t0!=tr){
      fun+=Q0*rc*(exp_rc*rc/((rc-t0)*(rc-tr))+exp(-(x-x0)/t0)*t0/((t0-rc)*(t0-tr))+exp(-(x-x0)/tr)*tr/((tr-rc)*(tr-t0)));
    }
    else{
      fun+=Q0*rc*(exp_rc*rc*t0-exp(-(x-x0)/t0)*(rc*(x-x0+t0)-(x-x0)*t0))/(t0*pow(rc-t0,2));
    }
  }
  return fun;
}

//______________________________________________________________________
TKratAnaStep1::TKratAnaStep1() :
  FairTask("[TKratAnaStep1:] created ")
  {
  fVerbose = 1;
  
     //ievt_tot = 0;
     //ievt_tot_good = 0;

     //mod = -1;
     //isig = -1;

     tls=2000;
     tlo=2000;
     tm=0;
     am=0;
     at=0;

     p0=0;
     p1=0;
     p2=0;
     p3=0;
     p4=0;
     p5=0;
     p6=0;
     p7=0;
     p8=0;
     p9=0;
     p10=0;
     p11=0;
     chi2=0;
    fEventCounter = 0;
    fCurrentRun = 0;
    fSaveOutputToTree   = kFALSE;
    fKratContOutName    = "KRATTA_EVENT_PARAM_CLONE"; 
    fKratContInName     = "KRATTARAWEVENT";  
    fMbstsContInName    = "MBSTSCLONE"; 
    fMbsinfoContInName  = "MBSINFO"; 
    
    
}


TKratAnaStep1::TKratAnaStep1(Int_t verbose) :
  FairTask("[TKratAnaStep1:] created ", verbose)
{
  fVerbose = verbose;
}

//______________________________________________________________________
TKratAnaStep1::~TKratAnaStep1() {
}

//______________________________________________________________________
void TKratAnaStep1::SetParContainers() {

  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  if ( ! rtdb ) Fatal("SetParContainers", "No runtime database");

}



//______________________________________________________________________
InitStatus TKratAnaStep1::Init() {

  /// Get input data containers:
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No FairRootManager");
  
  fKrattaRawEventCopy = (TClonesArray*)  ioman->GetObject( fKratContInName );
  fKrattaRawEvent     = (TKRATRawEvent*) fKrattaRawEventCopy;
  
  fRootMBSTSCopy =  (TClonesArray*) ioman->GetObject( fMbstsContInName );
  fMBSTS = (TRootTS*)fRootMBSTSCopy;
  
  fMbsInfoClone =  (TClonesArray*) ioman->GetObject( fMbsinfoContInName );
  fMbsInfo  = (TMbsInfo*)fMbsInfoClone;
  
  /// Check input data containers:
  if (fKrattaRawEventCopy==NULL){
    cout << "[TKratAnaStep1::Init:] Error: Couldn't find necessary input ";
    cout << "data container with name: " << fKratContInName << " (of class TKRATRawEvent).";
    cout << " It will breaks the analysis." << endl;
    getchar();
  }
  if (fRootMBSTSCopy==NULL){
    cout << "[TKratAnaStep1::Init:] Error: Couldn't find necessary input ";
    cout << "data container with name: " << fMbstsContInName << " (of class TRootTS).";
    cout << " It will breaks the analysis." << endl;
    getchar();
  }
  if (fMbsInfoClone==NULL){
    cout << "[TKratAnaStep1::Init:] Error: Couldn't find necessary input ";
    cout << "data container with name: " << fMbsinfoContInName << " (of class TMbsInfo).";
    cout << " It will breaks the analysis." << endl;
    getchar();
  }
  
  /// Alocating memory for used objects:
  fASYEvent = new ASYEvent();

  ///for (int index_modules = 0;index_modules < 35; index_modules++){
  ///  fASYEvent-> AddPeak();
  ///}
  fPeakData = new PeakData();
  
  fASYEventCopy = new  TClonesArray("ASYEvent");
  fASYEventCopy = (TClonesArray*) fASYEvent;
  
  ///ioman->Register("KRATTA_EVENT_PARAM", "Kratta ASYEvent Tree", fASYEvent, kFALSE); /// Commented 2014-03-20 by SK TODO: check if  necessary
  ioman->Register(fKratContOutName , "Kratta ASYEvent Tree CLONE",fASYEventCopy, fSaveOutputToTree);
  
  if ( ! gGeoManager )cout << "Init " << " No GeoManager"<<endl;
  
  
  /// Load parameters from files: GEO, MAPPING, OTHER..
  string dir= gSystem->Getenv("VMCWORKDIR");
  
  string fGeoDir ( dir + "/KRATTA/krattaTasks/TKratAnaStep1/mktreedst/geo/"); 
  string fGeoMap ( fGeoDir + "geo_map.geo"); 
  string fGeoPar ( fGeoDir + "parRC3GAM2DEXP13.dat");
  
  cout << "[TKratAnaStep1:] fGeoMap = " << fGeoMap << endl;
  cout << "[TKratAnaStep1:] fGeoPar = " << fGeoPar << endl;
  
  gg = new ASYgeo();
  gg -> ReadGeo( fGeoMap.c_str() );
  gg -> LoadParRC3GAM2DEXP13(fCurrentRun, fGeoPar.c_str() );
  
  for(int i=0;i<35;i++){
    printf("%2d %4.2f %4.2f %4.2f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f %5.1f\n",
    i,
    gg->P4[i][0],gg->P4[i][1],gg->P4[i][2],
    gg->P5[i][0],gg->P5[i][1],gg->P5[i][2],
    gg->P6[i][0],gg->P6[i][1],gg->P6[i][2],
    gg->P7[i][0],gg->P7[i][1],gg->P7[i][2],
    gg->P8[i][0],gg->P8[i][1],gg->P8[i][2],
    gg->P9[i][0],gg->P9[i][1],gg->P9[i][2],
    gg->P10[i][0]/gg->P4[i][0],gg->P10[i][1]/gg->P4[i][1],gg->P10[i][2]/gg->P4[i][2]);
  }
  
  
  /// Define fitter
  TVirtualFitter::SetDefaultFitter("Fumili");
  fitter=TVirtualFitter::Fitter(0,12);


  /// Allocate histograms
  hfit_512 = new TH1D("hfit_512","hfit_512",512,0.,512.);
  funn_512  = new TF1("funn_512",RC4DEXP2TR,0.,512.,12);
  funn0_512 = new TF1("funn0_512",RC4DEXP2TR1,0,512.,5);
  funn1_512 = new TF1("funn1_512",RC4DEXP2TR1,0,512.,5);
  funn2_512 = new TF1("funn2_512",RC4DEXP2TR1,0,512.,5);
  funn3_512 = new TF1("funn3_512",RC4DEXP2TR1,0,512.,5);

  hfit_1024 = new TH1D("hfit_1024","hfit_1024",1024,0.,1024.);
  funn_1024  = new TF1("funn_1024",RC4DEXP2TR,0.,1024.,12);
  funn0_1024 = new TF1("funn0_1024",RC4DEXP2TR1,0,1024.,5);
  funn1_1024 = new TF1("funn1_1024",RC4DEXP2TR1,0,1024.,5);
  funn2_1024 = new TF1("funn2_1024",RC4DEXP2TR1,0,1024.,5);
  funn3_1024 = new TF1("funn3_1024",RC4DEXP2TR1,0,1024.,5);

  fp2    = new TF1("fp2","pol3",0.,1024.);
  mrk    = new TMarker(0,1,29);
  mrk0   = new TMarker(0,1,29);
  marker = new TMarker(0,1,29);
  li1e   = new TLine(0,0,0,16000);
  fexp   = new TF1("fexp","exp(-x/[0]+[1])",0,1024);
  gr     = new TGraph();
  mrkp   = new TMarker(0,1,29);
  fp21   = new TF1("fp21","pol2",0,1);
  fp23   = new TF1("fp23","[0]*x+[1]*x*x+[2]*x*x*x",0,1);
  
  /// Done
  cout << "[TKratAnaStep1:] Init ............................ done" << endl;     
  cout << "[TKratAnaStep1:] Press Enter to continue..." << endl;     
  getchar();
  
  return kSUCCESS;

}


//______________________________________________________________________
void TKratAnaStep1::Exec(Option_t* opt)
{
  
   ///cout<<"[TKratAnaStep1:] New event..............................."  << (fKrattaRawEvent->ToString()).Data() <<endl;
   
   ///Int_t krattaMult = fKrattaRawEvent->GetKRATTAMultiplicity();
   ///cout<<"[TKratAnaStep1:] Multiplicity = " <<  krattaMult << endl;
   ///Int_t MBS_TimeStamp=    fMBSTS->TSM*65535+fMBSTS->TSS;   
   ///cout<<"[TKratAnaStep1:] MBS_TimeStamp = " <<  MBS_TimeStamp << endl;
   
   fEventCounter++;
   Int_t FiredMoudlesCounter = 0;
   ///TODO: ZEROWAC fEventCounter na poczatku kazdego run'u.
   
   Bool_t CurrentModuleFired  = kFALSE;
   
   
   fASYEvent->Clear(); ///TODO Has to be cleared 
   
   for (int imodule = 0; imodule < 35; imodule++ ){ /// Loop over modules
     
     CurrentModuleFired = kFALSE;
     fPeakData->Clear();
     for (int ph_signal = 0; ph_signal < 3; ph_signal++ ){ /// Loop over photodiodes
     
       if ( fKrattaRawEvent-> ContainSignalFromPh ( imodule, ph_signal)  )
       {
         CurrentModuleFired = kTRUE;
         
         ifadc    = (Int_t *) fKrattaRawEvent->GetSignalArray ( imodule, ph_signal  );
         n_sample = fKrattaRawEvent->GetSignalArraySize ( imodule, ph_signal );
         GetPulseParameters( imodule, ph_signal);
         
         ///cout << "[TKratAnaStep1:] FOUND SIGNAL FROM module=" << imodule << ", signal="<<  ph << endl;
         ///fPhotodiodePtr = fKrattaRawEvent->GetPhotodiode ( imodule, ph_signal );
         ///cout <<  fPhotodiodePtr->ToString() << endl;
         
         ///Int_t photodiodeID = fPhotodiodePtr->GetId();
         ///Int_t module = 7 * fPhotodiodePtr->GetRow_Int() + fPhotodiodePtr->GetColumn_Int();
         ///Int_t signal = fPhotodiodePtr->GetPhotodiode_Int();
         ///cout << "[TKratAnaStep1:] FOUND SIGNAL FROM module=" <<module << ", signal="<<  ph << endl;
         
       }
     }
     
     if ( CurrentModuleFired )
     {
       
       
       fPeakData->mod = imodule;
       
       fASYEvent-> AddPeak();
       FillASYEvent( FiredMoudlesCounter );
       
       FiredMoudlesCounter++;
     }
     
   } ///End loop over modules
   
   fASYEvent->npeak = FiredMoudlesCounter; ///Or size of the ClonesArray
   fASYEvent->evt = fMbsInfo->GetCamacEventNumber();
   fASYEvent->run = fMbsInfo->GetRunNumber();
   fASYEvent->tms = fMBSTS->TSM*65535+fMBSTS->TSS; 
   /**
   for (int ph = 0; ph < 105; ph++ ){
     
     if (fKrattaRawEvent-> ContainSignalFromPh(ph)  )
     {
       fPhotodiodePtr = fKrattaRawEvent->GetPhotodiode (ph );
       
       ifadc    = (Int_t *) fKrattaRawEvent->GetSignalArray ( ph );
       n_sample = fKrattaRawEvent->GetSignalArraySize ( ph );
       
       Int_t photodiodeID = fPhotodiodePtr->GetId();
       Int_t module = 7 * fPhotodiodePtr->GetRow_Int() + fPhotodiodePtr->GetColumn_Int();
       
       Int_t signal = fPhotodiodePtr->GetPhotodiode_Int();
       
       cout << "PHotodiode [" << module << "." << signal <<  "]";
       cout << " with id= " << photodiodeID << " stored pulse len = "<< n_sample << endl;
       cout <<  "  [";
       //for (int ii = 0; ii < n_sample; ii++ )
       //{
       //  cout << ", "  << ifadc[ ii ];
       //}
       //cout << "]" << endl;
       
       
       mod  = module;
       isig = signal;
       if ( module == 7 ) {
         GetPulseParameters();
       }
     }
     else{
      // cout << " NO photodiode " << ph << endl;
     }
    // getchar();
    
     /// MOD SYGN SIZE ARRAY
   }
   */
}
//______________________________________________________________________
void TKratAnaStep1::FillASYEvent( Int_t A_FiredModuleCounter )
{
  
  ASYFadcPeak * ptrASYFadcPeak;
  ptrASYFadcPeak = fASYEvent->Peak( A_FiredModuleCounter );

  ///cout << "[TKratAnaStep1:] IM INSIDE FILLING FIELD" << endl;

  ptrASYFadcPeak->mod = fPeakData->mod;
  
  ///ALERT A_Module == ptrASYFadcPeak->mod 
 
  ptrASYFadcPeak->pds0 =fPeakData->pds[0];
  ptrASYFadcPeak->pds1 =fPeakData->pds[1];
  ptrASYFadcPeak->pds2 =fPeakData->pds[2];

  //cout << " fPeakData->pds[0] = " << fPeakData->pds[0] << endl;
  //cout << " fPeakData->pds[1] = " << fPeakData->pds[1] << endl;
  //cout << " fPeakData->pds[2] = " <<  fPeakData->pds[2] << endl;

  //cout << " ptrASYFadcPeak->pds0 = " << ptrASYFadcPeak->pds0 << endl;
  //cout << " ptrASYFadcPeak->pds1 = " << ptrASYFadcPeak->pds1 << endl;
  //cout << " ptrASYFadcPeak->pds2 = " <<  ptrASYFadcPeak->pds2<< endl;
  
  ptrASYFadcPeak->pdr0 =fPeakData->pdr[0];
  ptrASYFadcPeak->pdr1 =fPeakData->pdr[1];
  ptrASYFadcPeak->pdr2 =fPeakData->pdr[2];

  ptrASYFadcPeak->tls0  =fPeakData->tls[0];
  ptrASYFadcPeak->tls1  =fPeakData->tls[1];
  ptrASYFadcPeak->tls2  =fPeakData->tls[2];

  ptrASYFadcPeak->tlo0  =fPeakData->tlo[0];
  ptrASYFadcPeak->tlo1  =fPeakData->tlo[1];
  ptrASYFadcPeak->tlo2  =fPeakData->tlo[2];

  ptrASYFadcPeak->tm0  =fPeakData->tm[0];
  ptrASYFadcPeak->tm1  =fPeakData->tm[1];
  ptrASYFadcPeak->tm2  =fPeakData->tm[2];

  ptrASYFadcPeak->am0  =fPeakData->am[0];
  ptrASYFadcPeak->am1  =fPeakData->am[1];
  ptrASYFadcPeak->am2  =fPeakData->am[2];

  ptrASYFadcPeak->at0  =fPeakData->at[0];
  ptrASYFadcPeak->at1  =fPeakData->at[1];
  ptrASYFadcPeak->at2  =fPeakData->at[2];

  ptrASYFadcPeak->sam0  =fPeakData->sam[0];
  ptrASYFadcPeak->sam1  =fPeakData->sam[1];
  ptrASYFadcPeak->sam2  =fPeakData->sam[2];

  ptrASYFadcPeak->sat0  =fPeakData->sat[0];
  ptrASYFadcPeak->sat1  =fPeakData->sat[1];
  ptrASYFadcPeak->sat2  =fPeakData->sat[2];

  ptrASYFadcPeak->s01  =fPeakData->s1[0];
  ptrASYFadcPeak->s02  =fPeakData->s2[0];
  ptrASYFadcPeak->s03  =fPeakData->s3[0];
  ptrASYFadcPeak->s04  =fPeakData->s4[0];
  ptrASYFadcPeak->s05  =fPeakData->s5[0];
  ptrASYFadcPeak->s06  =fPeakData->s6[0];
  ptrASYFadcPeak->s07  =fPeakData->s7[0];
  ptrASYFadcPeak->s08  =fPeakData->s8[0];
  ptrASYFadcPeak->s09  =fPeakData->s9[0];
  ptrASYFadcPeak->s09f =fPeakData->s9f[0];

  ptrASYFadcPeak->s11  =fPeakData->s1[1];
  ptrASYFadcPeak->s12  =fPeakData->s2[1];
  ptrASYFadcPeak->s13  =fPeakData->s3[1];
  ptrASYFadcPeak->s14  =fPeakData->s4[1];
  ptrASYFadcPeak->s15  =fPeakData->s5[1];
  ptrASYFadcPeak->s16  =fPeakData->s6[1];
  ptrASYFadcPeak->s17  =fPeakData->s7[1];
  ptrASYFadcPeak->s18  =fPeakData->s8[1];
  ptrASYFadcPeak->s19  =fPeakData->s9[1];
  ptrASYFadcPeak->s19f =fPeakData->s9f[1];

  ptrASYFadcPeak->s21  =fPeakData->s1[2];
  ptrASYFadcPeak->s22  =fPeakData->s2[2];
  ptrASYFadcPeak->s23  =fPeakData->s3[2];
  ptrASYFadcPeak->s24  =fPeakData->s4[2];
  ptrASYFadcPeak->s25  =fPeakData->s5[2];
  ptrASYFadcPeak->s26  =fPeakData->s6[2];
  ptrASYFadcPeak->s27  =fPeakData->s7[2];
  ptrASYFadcPeak->s28  =fPeakData->s8[2];
  ptrASYFadcPeak->s29  =fPeakData->s9[2];
  ptrASYFadcPeak->s29f =fPeakData->s9f[2];
  
  /// falling edge:
  ptrASYFadcPeak->sf00 =fPeakData->sf0[0];
  ptrASYFadcPeak->sf01 =fPeakData->sf1[0];
  ptrASYFadcPeak->sf02 =fPeakData->sf2[0];
  ptrASYFadcPeak->sf03 =fPeakData->sf3[0];
  ptrASYFadcPeak->sf04 =fPeakData->sf4[0];
  ptrASYFadcPeak->sf05 =fPeakData->sf5[0];
  ptrASYFadcPeak->sf06 =fPeakData->sf6[0];
  ptrASYFadcPeak->sf07 =fPeakData->sf7[0];
  ptrASYFadcPeak->sf08 =fPeakData->sf8[0];
  ptrASYFadcPeak->sf09 =fPeakData->sf9[0];

  ptrASYFadcPeak->sf10 =fPeakData->sf0[1];
  ptrASYFadcPeak->sf11 =fPeakData->sf1[1];
  ptrASYFadcPeak->sf12 =fPeakData->sf2[1];
  ptrASYFadcPeak->sf13 =fPeakData->sf3[1];
  ptrASYFadcPeak->sf14 =fPeakData->sf4[1];
  ptrASYFadcPeak->sf15 =fPeakData->sf5[1];
  ptrASYFadcPeak->sf16 =fPeakData->sf6[1];
  ptrASYFadcPeak->sf17 =fPeakData->sf7[1];
  ptrASYFadcPeak->sf18 =fPeakData->sf8[1];
  ptrASYFadcPeak->sf19 =fPeakData->sf9[1];

  ptrASYFadcPeak->sf20 =fPeakData->sf0[2];
  ptrASYFadcPeak->sf21 =fPeakData->sf1[2];
  ptrASYFadcPeak->sf22 =fPeakData->sf2[2];
  ptrASYFadcPeak->sf23 =fPeakData->sf3[2];
  ptrASYFadcPeak->sf24 =fPeakData->sf4[2];
  ptrASYFadcPeak->sf25 =fPeakData->sf5[2];
  ptrASYFadcPeak->sf26 =fPeakData->sf6[2];
  ptrASYFadcPeak->sf27 =fPeakData->sf7[2];
  ptrASYFadcPeak->sf28 =fPeakData->sf8[2];
  ptrASYFadcPeak->sf29 =fPeakData->sf9[2];
  
  ///moments:
  ptrASYFadcPeak->m00  =fPeakData->m0[0];
  ptrASYFadcPeak->m01  =fPeakData->m1[0];
  ptrASYFadcPeak->m02  =fPeakData->m2[0];
  ptrASYFadcPeak->m03  =fPeakData->m3[0];
  ptrASYFadcPeak->m04  =fPeakData->m4[0];
  ptrASYFadcPeak->m05  =fPeakData->m5[0];
  ptrASYFadcPeak->m06  =fPeakData->m6[0];
  ptrASYFadcPeak->m07  =fPeakData->m7[0];
  ptrASYFadcPeak->m08  =fPeakData->m8[0];
  ptrASYFadcPeak->m09  =fPeakData->m9[0];

  ptrASYFadcPeak->m10  =fPeakData->m0[1];
  ptrASYFadcPeak->m11  =fPeakData->m1[1];
  ptrASYFadcPeak->m12  =fPeakData->m2[1];
  ptrASYFadcPeak->m13  =fPeakData->m3[1];
  ptrASYFadcPeak->m14  =fPeakData->m4[1];
  ptrASYFadcPeak->m15  =fPeakData->m5[1];
  ptrASYFadcPeak->m16  =fPeakData->m6[1];
  ptrASYFadcPeak->m17  =fPeakData->m7[1];
  ptrASYFadcPeak->m18  =fPeakData->m8[1];
  ptrASYFadcPeak->m19  =fPeakData->m9[1];

  ptrASYFadcPeak->m20  =fPeakData->m0[2];
  ptrASYFadcPeak->m21  =fPeakData->m1[2];
  ptrASYFadcPeak->m22  =fPeakData->m2[2];
  ptrASYFadcPeak->m23  =fPeakData->m3[2];
  ptrASYFadcPeak->m24  =fPeakData->m4[2];
  ptrASYFadcPeak->m25  =fPeakData->m5[2];
  ptrASYFadcPeak->m26  =fPeakData->m6[2];
  ptrASYFadcPeak->m27  =fPeakData->m7[2];
  ptrASYFadcPeak->m28  =fPeakData->m8[2];
  ptrASYFadcPeak->m29  =fPeakData->m9[2];

  ptrASYFadcPeak->mx0  =fPeakData->mx[0];
  ptrASYFadcPeak->mx1  =fPeakData->mx[1];
  ptrASYFadcPeak->mx2  =fPeakData->mx[2];

  ptrASYFadcPeak->my0  =fPeakData->my[0];
  ptrASYFadcPeak->my1  =fPeakData->my[1];
  ptrASYFadcPeak->my2  =fPeakData->my[2];

  ptrASYFadcPeak->mz0  =fPeakData->mz[0];
  ptrASYFadcPeak->mz1  =fPeakData->mz[1];
  ptrASYFadcPeak->mz2  =fPeakData->mz[2];

  ptrASYFadcPeak->ma0  =fPeakData->ma[0];
  ptrASYFadcPeak->ma1  =fPeakData->ma[1];
  ptrASYFadcPeak->ma2  =fPeakData->ma[2];

  ptrASYFadcPeak->mb0  =fPeakData->mb[0];
  ptrASYFadcPeak->mb1  =fPeakData->mb[1];
  ptrASYFadcPeak->mb2  =fPeakData->mb[2];

  ptrASYFadcPeak->mc0  =fPeakData->mc[0];
  ptrASYFadcPeak->mc1  =fPeakData->mc[1];
  ptrASYFadcPeak->mc2  =fPeakData->mc[2];

  ptrASYFadcPeak->p00  =fPeakData->p0[0];
  ptrASYFadcPeak->p04  =fPeakData->p4[0];
  ptrASYFadcPeak->p05  =fPeakData->p5[0];
  ptrASYFadcPeak->p09  =fPeakData->p9[0];
  ptrASYFadcPeak->rc0  =fPeakData->rc[0];
  ptrASYFadcPeak->chi20  =fPeakData->chi2[0];

  ptrASYFadcPeak->p10  =fPeakData->p0[1];
  ptrASYFadcPeak->p12  =fPeakData->p2[1]; 
  ptrASYFadcPeak->p13  =fPeakData->p3[1];
  ptrASYFadcPeak->p14  =fPeakData->p4[1];
  ptrASYFadcPeak->p15  =fPeakData->p5[1];
  ptrASYFadcPeak->p17  =fPeakData->p7[1];
  ptrASYFadcPeak->p18  =fPeakData->p8[1];
  ptrASYFadcPeak->p19  =fPeakData->p9[1];
  ptrASYFadcPeak->p110 =fPeakData->p10[1];
  ptrASYFadcPeak->rc1  =fPeakData->rc[1];
  ptrASYFadcPeak->chi21  =fPeakData->chi2[1];

  ptrASYFadcPeak->p20  =fPeakData->p0[2];
  ptrASYFadcPeak->p22  =fPeakData->p2[2];
  ptrASYFadcPeak->p23  =fPeakData->p3[2];
  ptrASYFadcPeak->p24  =fPeakData->p4[2];
  ptrASYFadcPeak->p25  =fPeakData->p5[2];
  ptrASYFadcPeak->p27  =fPeakData->p7[2];
  ptrASYFadcPeak->p28  =fPeakData->p8[2];
  ptrASYFadcPeak->p29  =fPeakData->p9[2];
  ptrASYFadcPeak->p210 =fPeakData->p10[2];
  ptrASYFadcPeak->rc2  =fPeakData->rc[2];
  ptrASYFadcPeak->chi22  =fPeakData->chi2[2];
}

//______________________________________________________________________
///Parameters:  A_Module - module number [0..34], A_Signal - signal id from photodiode [0..2]
void TKratAnaStep1::GetPulseParameters( Int_t A_Module, Int_t A_Signal )
{

  Int_t mod  = A_Module;
  Int_t isig = A_Signal;

  //find max/min
  int chmin = 0;
  int chmax = 0;
  int vmin =  1.e9;
  int vmax = -1.e9;
  
  for(int is = 0;is<n_sample;is++){
    if(ifadc[is]<vmin){
      vmin=ifadc[is];
      chmin = is;
    }
    if(ifadc[is]>vmax){
      vmax=ifadc[is];
      chmax = is;
    }
  }
  //find max/min

  //pds 0-approx	    
  double pds0   = 0;
  double pds00   = 0;
  double pds    = 0;
  double pdsrms = 0;
  double pdsrms00 = 0;
  double pds1m  = 0;
  double pds1r  = 0; 
  double rpds1m = 0; 
  double rpds1r = 0;
  double pds1m2 = 0; 

  int np = n_sample;
  double fnp = 0;
  int np0 = 10;
  for(int i=0; i<np0; i++){
    fnp   +=1;
    pds1m += double(ifadc[i]); 
    pds1r += pow(double(ifadc[i]),2); 
  }

  pds0    = pds1m/fnp;
  pdsrms = sqrt(fabs(pds1r/fnp-pow(pds0,2)));

  for(int i=np0; i<n_sample; i++){
    fnp   +=1;
    pds1m += double(ifadc[i]); 
    pds1r += pow(double(ifadc[i]),2); 
  }
  pds00    = pds1m/fnp;
  pdsrms00 = sqrt(fabs(pds1r/fnp-pow(pds00,2)));
  //pds 0-approx	    

  // printf("#     %9d %d   -------------------------------------\n",ievt,ich);

  if(pds00<pds0) pds0=pds00;
  // treat events without flat part at the beginning
  if(ifadc[n_sample-1]<ifadc[0]-10) pds0=double(ifadc[n_sample-1]);//!!!!!?????
  //if(ifadc[n_sample-1]<ifadc[0]-20) pds0=double(ifadc[n_sample-1]);//!!!!!?????

//   printf("# pds0 %9.2f\n", pds0);


  double t00=0;
  double t000=0;
  double a000=0;
  double b000=0;
  double t035=0;
  double t025=0;
  double t005=0;
  double y035=0;
  double y025=0;
  double y005=0;
  double dnpt=0;
  double hmaxi = 0;
  int icount=0;
  int imin = (chmax>=50) ? chmax-50 : 0;
  int imax = (chmax<np-50) ? chmax+50 : np-1;
  
  for(int i=imin; i<=imax; i++){
    hmaxi += double(ifadc[i]);
    icount++;
  }
  hmaxi /= float(icount);
  for(int i=chmax; i>=1; i--){
    //   if((double(ifadc[i])-pds0)<0.35*(hmaxi-pds0) && i<400){
    if((double(ifadc[i])-pds0)<0.35*(hmaxi-pds0)){
      t035 = double(i);
      y035 = double(ifadc[i]);
      break;
    }
  }
  for(int i=int(t035); i>=1; i--){
    if((double(ifadc[i])-pds0)<=0.25*(hmaxi-pds0)){
      t025 = double(i);
      y025 = double(ifadc[i]);
      break;
    }
  }
  for(int i=int(t025); i>=1; i--){
    if((double(ifadc[i])-pds0)<0.05*(hmaxi-pds0)){
      t005 = double(i);
      y005 = double(ifadc[i]);
      break;
    }
  }

  // 2 lines
  Double_t As00=0.,Bs00=0.;
  Double_t As20=0.,Bs20=0.;
  Double_t fT10 = 0;
  int T10 = 0;
  int Tmin =int(t005)-20;
  if(Tmin<1) Tmin=1;
  //Tmin = 1; //!!!!!!!!!!!!1
  int Tmax =int(t035);
  {

    double Dmin=1.e10;

    Double_t As0=0.,Bs0=0.;
    Double_t As2=0.,Bs2=0.;
    Double_t sw1, swy1;
    Double_t sw2, swx2, swy2, swxx2, swxy2;
    Double_t xi,yi;
    int iii1 = 0;
    int iii2 = 0;

    sw1=swy1=0.;
    sw2=swxy2=swx2=swy2=swxx2=0.;
    for(int T1=1;T1<Tmin;T1++){
      xi=double(T1)-0.5;
      yi=double(ifadc[T1-1]);
      sw1	+= 1.;
      swy1  += yi;
      iii1++;
    }
    for(int T1=Tmin;T1<=Tmax;T1++){
      xi=double(T1)-0.5;
      yi=double(ifadc[T1-1]);
      sw2	+= 1.;
      swx2  += xi;
      swy2  += yi;
      swxx2 += xi*xi;
      swxy2 += xi*yi;
      iii2++;
    }
    for(int T1=Tmin;T1<=Tmax-3;T1++){
      xi=double(T1)-0.5;
      yi=double(ifadc[T1-1]);
      sw1	+= 1.;
      swy1  += yi;
      iii1++;

      sw2	-= 1.;
      swx2  -= xi;
      swy2  -= yi;
      swxx2 -= xi*xi;
      swxy2 -= xi*yi;
      iii2--;

      if(sw1==0){printf("!!!!! sw1==0\n");exit(1);}
      As0 = 0;
      Bs0 = swy1/sw1;

      if(sw2==0){printf("!!!!! sw2==0\n");exit(1);}
      double deno = sw2*swxx2-swx2*swx2;
      if(deno==0){printf("!!!!! sw2*swxx2-swx2*swx2==0\n");exit(1);}
      As2 = (sw2*swxy2-swx2*swy2)/(sw2*swxx2-swx2*swx2);
      Bs2 = (swy2-As2*swx2)/sw2;

      double Dmin12=0;
      for(int it=1;it<=T1;it++){
        xi=double(it)-0.5;
        yi=double(ifadc[it-1]);
        Dmin12 += pow(As0*xi+Bs0-yi,2);
      }
      for(int it=T1+1;it<=Tmax;it++){
        xi=double(it)-0.5;
        yi=double(ifadc[it-1]);
        Dmin12 += pow(As2*xi+Bs2-yi,2);
      }


      //   printf("--> %d %f %f | %f %f %f %f %d\n",
      //   T1,Dmin12,Dmin,(Bs0-Bs2)/As2,fabs((Bs0-Bs2)/As2-T1-1),Bs0,As2,iii2);
      if(Dmin12<Dmin && As2>0 && fabs((Bs0-Bs2)/As2-T1-1)<10){// 
        As00=As0;Bs00=Bs0;
        As20=As2;Bs20=Bs2;
        Dmin=Dmin12;
        T10 = T1;
      }
    }
    if(As20!=0){
      fT10 = (Bs00-Bs20)/As20;
    }
    else{
      fT10 = T10;
      Bs00 = pds0;
      //printf("!!!!!!!!!!!!As20==0\n");
    }
    //fT10 -=1;
    //printf("T10 %3d %f |  %f %f %f\n",T10,fT10,Bs00,As20,Bs20);
  }
  // double tt00 = t005-2.*(As20*t005+Bs20-Bs00)/As20; //styczna w t005 i min w tt00
  // printf("tt00 %f\n",tt00);
  // 2 lines

  t00=fT10;
  if(t00<0) t00=0; // bez tego segm.viol.
  pds=Bs00;

//  printf("# pds1 %f %f %f\n", pds,t00,fT10);

  if(n_sample==512){
    hfit = hfit_512;
    funn = funn_512;
    funn0 = funn0_512;
    funn1 = funn1_512;
    funn2 = funn2_512;
    funn3 = funn3_512;
    //fexp = fexp_512;
  }
  else  {
    hfit = hfit_1024;
    funn = funn_1024;
    funn0 = funn0_1024;
    funn1 = funn1_1024;
    funn2 = funn2_1024;
    funn3 = funn3_1024;
    //fexp = fexp_1024;
  }
  
  
  hfit->Reset();    
  am=-1.e10;
  for(int i=0; i<n_sample; i++){
  
    hfit->SetBinContent( i+1,    double( ifadc[i]) - pds );
    hfit->SetBinError  ( i+1, sqrt(fabs( double(ifadc[i])  - pds )) );
    if( (ifadc[i] - pds) >am) am = ifadc[i] - pds;
  }
  
  
  //// fit tail ////


  double sxi=0;
  double syi=0;
  double sxxi=0;
  double sxyi=0;
  double si=0;
  int ibeg=hfit->GetMaximumBin()+150;
  if(n_sample>512) ibeg+=150;
  if(n_sample-ibeg<10) ibeg=n_sample-100;
  //printf("ibeg %d\n",ibeg);
  for(int i=ibeg; i<n_sample; i++){
    double xi = hfit->GetBinCenter(i+1);
    double yi = hfit->GetBinContent(i+1);
    if(yi>0){
      yi=log(yi);
      sxi+=xi;
      syi+=yi;
      sxxi+=xi*xi;
      sxyi+=xi*yi;
      si+=1;
    }
  }
  double den = si*sxyi-sxi*syi;
  double tail = 0.1;
  if(den!=0) tail = -(si*sxxi-sxi*sxi)/den;
  double btail = (syi+sxi/tail)/si;
  //// fit tail ////
  
  
  
  double Q0     = 1;
  double Q1     = 1;
  double Q2     = 1;
  double x0     = 0;
  //double t0     = 6.7;
  double t0     = gg->P4[mod][isig];
  double t1     = 43;
  double t9     = gg->P4[mod][isig]*2.813;//Si risetime (9.2)
  //double t9     = gg->P4[mod][isig]*2.887;//Si risetime (9.44)
  //double t9     = 9.2;
  //double t9     = gg->P4[mod][isig]*3.06;//Si risetime (10)
  //double t9     = gg->P4[mod][isig]*2.45;//Si risetime (8)
  double t2     = 0.8*gg->P8[mod][isig]; //CsI  fast decay
  double t3     = 0.85*gg->P9[mod][isig];//CsI slow decay
  //if(isig==1 && (*iter)->am[2]>0) t3 *= 1.15; //increase for punch through
  // double t10    = 16;
  //double t10    = gg->P4[mod][isig]*4.893;//CsI risetime (16)
  //double t10    = gg->P10[mod][isig]*1.575;//CsI risetime (16)
  double t10    = gg->P10[mod][isig]*1.378;//CsI risetime (14)
  if(mod==30) t10 *= 1.30;
  if(mod== 3) t10 *= 0.95;
  //if(MOD== 9) t10 *= 0.8;
  //if(MOD==3 ||MOD==9 ||MOD==26) t10 *= 1.05;
  //double t10    = t9;
  double rc     = gg->P7[mod][isig];
  //if(isig==1 && (*iter)->am[2]>0) t10 = 10; //increase for punch through


  tls = tail;
  tlo = btail;

  double rat = (tail-rc)/rc;

  // for(int i=0; i<n_sample; i++){
  // //   double yi = exp(-(double(i)+0.5)/tail+btail); // extrapolate using exp fit
  //   double yi = exp(-(double(i)+0.5-t00)/rc); // extrapolate using exp fit
  //    hfit->SetBinContent(i+1,(*iter)->dfadc[isig][i]/yi);
  // 
  //   //hfit->SetBinError(i+1,sqrt(fabs((*iter)->dfadc[isig][i])));
  // }


  double parRCGAM12SI[20];
  double parRCGAM12SI0[20];
  double parRCGAM12SI1[20];
  double parRCGAM12SI2[20];
  double parRCGAM12SI3[20];
  double xx[2];

  double PARPD1[5]={0};


  double tx0=1;
  double tx1=2;
  double ty0=1;
  double ty1=2;
  double Qy0=0;
  double Qy1=0;
  double Qx0=0;
  double Qx1=0;
  double tz0=1.5;
  double tz1=10.4;
  double Qz0=am;
  double Qz1=am;
  double xi0=0;
  double xi00=t00;
  double DIFF0=1.e10;
  int    dzi=0;

  //double rc1 = tail;
  double rc1 = rc;

  double ts0=1.5;
  double ts1=t9;
  double Qs0=am;


  double xpar[5]={0};
  double ypar[5]={0};
  int    ipar=0;
  icount =0;
  //bool Done=false;
  bool Done=true;
  double xdif[5]={0,-0.2,0.2,0,0};
  //do{

  int it01 = 0;
  int it010 = int(t00-0.5);
  //int it01 = 188;
  for(int j=it010-10;j<=it010+5;j++){
  it01 = j;
  //int nbins = hfit->GetMaximumBin()-it01;
  //int nbins = n_sample-it01-1;
  // int nbins = (n_sample-it01)/3;
  // int nbins = (hfit->GetMaximumBin()-it01)/2;
  //int nbins = 15;
  int nbins = 18;
  if(isig==0) nbins = (n_sample-it01)/2;
  xi0 = hfit->GetBinCenter(it01+1);
  // if(ipar==0) xi0 = hfit->GetBinCenter(it01+1)+xdif[ipar];
  // else        xi0 = xpar[0]+xdif[ipar];

  double M0=0;
  double M1=0;
  double M2=0;
  for(int i=it01; i<it01+nbins; i++){
    double xi = double(i-it01);
    double yi = hfit->GetBinContent(i+1);
    M0+=yi;
    M1+=xi*yi;
    M2+=xi*xi*yi;
  }
  // double TMAX=hfit->GetBinCenter(it01+nbins)-xi0;
  double TMAX=hfit->GetBinCenter(it01+nbins+1)-xi0;
  // printf("OO %d %d\n",it01,n_sample);
  // printf("OOO %f %f %f %f\n",xi0,TMAX,hfit->GetBinCenter(it01),am);

  double ERC = (1.-exp(-TMAX/rc1))*rc1*rc1;


  double e1 = ((1.-exp(-(TMAX/ts1)))*pow(ts1,2))/(rc1-ts1);

  ts0=(ERC*(M2-2*M1*(rc1+TMAX)+M0*TMAX*(2*rc1+TMAX))-(rc1-ts1)*
  (e1*(M2-2*M1*(ts1+TMAX)+M0*TMAX*(2*ts1+TMAX))+
  TMAX*(M2+2*M0*(rc1+ts1)*TMAX-M1*(2*(rc1+ts1)+TMAX))))/
  (2*ERC*(M1-M0*(rc1+TMAX))+(rc1-ts1)*(2*e1*(-M1+M0*(ts1+TMAX))+
  TMAX*(-2*M1+M0*(2*(rc1+ts1)+TMAX))));

  Qs0=((rc1-ts1)*(2*pow(M1,2)-2*M0*M1*TMAX+M0*(-M2+M0*pow(TMAX,2))))/
  (rc1*(2*ERC*(M1-M0*(rc1+TMAX))+(rc1-ts1)*(2*e1*(-M1+M0*(ts1+TMAX))+
  TMAX*(-2*M1+M0*(2*(rc1+ts1)+TMAX)))));


  //if(ty0<0 || ty1<0) continue;

  // printf("M0 -> %f,\n",M0);
  // printf("M1 -> %f,\n",M1);
  // printf("M2 -> %f,\n",M2);
  // printf("M3 -> %f,\n",M3);
  // printf("M4 -> %f,\n",M4);
  // printf("rc -> %f,\n",rc1);
  // printf("ERC -> %f,\n",ERC);
  // printf("TMAX -> %f,\n",TMAX);
  // printf("Q0 -> %f\n",Qy0);

  //printf("MM %g %g\n",tx0,tx1);
  //printf("Mx %7.4f %7.4f %7.4f\n",tx0,tx1,Qy0);
  // printf("MM %7.4f %7.4f %7.4f\n",ty0,ty1,Qy0);
  // printf("Ms %7.4f         %7.4f\n",ts0,Qs0);


  Qy0=Qs0;
  ty0=ts0;
  ty1=ts1;

  funn1->SetParameter(0,Qy0);
  funn1->SetParameter(1,xi0);
  funn1->SetParameter(2,ty0);
  funn1->SetParameter(3,ty1);
  funn1->SetParameter(4,rc1);

  //     double DIFF2=0;
  //     double DIFF =0;
  //     for(int i=0; i<n_sample; i++){
  //       double xi = double(i)+0.5;
  //       double yhi = hfit->GetBinContent(i+1);
  //       double yfi = funn1->Eval(xi);
  //       DIFF = yhi-yfi;
  //       DIFF2 += DIFF*DIFF;
  //     }
  // 
  // DIFF2/=(n_sample-5);
  // printf("DIFF2 %f %f\n",xi0,DIFF2);

  double DIFF2=0;
  double DIFF =0;
  for(int i=it01; i<it01+nbins; i++){
    double xi  = hfit->GetBinCenter(i+1);
    double yhi = hfit->GetBinContent(i+1);
    double yfi = funn1->Eval(xi);
    DIFF = yhi-yfi;
    DIFF2 += DIFF*DIFF;
  }

  DIFF2/=(n_sample-5);
  //printf("DIFF2 %f %f\n",xi0,DIFF2);

  //printf("MM %3d %7.4f %7.4f %7.4f %10.4f\n",it01,ty0,ty1,Qy0,DIFF2);



  if(isig>=1){
    if(ty0>0 && Qy0<DIFF0){
      DIFF0=Qy0;
      tz0=ty0;
      tz1=ty1;
      Qz0=Qy0;
      xi00=xi0;
      //printf("MIN\n");
    }
  }
  else{
    if(ty0>0 && DIFF2<DIFF0){
      DIFF0=DIFF2;
      tz0=ty0;
      tz1=ty1;
      Qz0=Qy0;
      xi00=xi0;
      //printf("MIN\n");
    }
  }


  }
  if(icount>10) Done=true;
  //}
  //while (!Done);

  ty0=tz0;
  ty1=tz1;
  Qy0=Qz0;
  xi0=xi00;
  //  printf("MM0 %7.4f %7.4f %7.4f %g %g\n",ty0,ty1,Qy0,xi0,DIFF0);

  hfit->SetTitle(Form("ev %d mod %2d sig %d",fEventCounter,mod,isig));
  //if(0){
  // if(isig>0){

  fitter->SetObjectFit(hfit);
  fitter->SetUserFunc(funn);
  fitter->SetFitMethod("chisquare");
  fitter->SetFCN(FCN);
  for(int ip=0;ip<12;ip++){ fitter->ReleaseParameter(ip);}
  fitter->Clear(); //!!!!! this is it !!!!!

  Double_t arglist[100];
  arglist[0] = 2;
  fitter->ExecuteCommand("SET PRINT", arglist, 1); // no print
  //  arglist[0] = 0;
  //  fitter->ExecuteCommand("SET NOWAR", arglist, 1); // no warnings


  double q0st = 1-2*rat;
  if(q0st<0.05) q0st=0.05;
  if(q0st>1)    q0st=1.;
  fitter->SetParameter(0,"Q0",am*q0st,.1,0.,2.*am);//Q0
  //  fitter->SetParameter(0,"Q0",Qy0,.1,0.,2.*am);//Q0

  //    fitter->SetParameter(1,"Q1",0.04*am*q0st     ,.1,0.,2.*am);//Q1
  fitter->SetParameter(1,"Q1",0.     ,.1,0.,2.*am);//Q1
  fitter->FixParameter(1);//Q1

  double q2st = 1.6*rat;
  if(q2st<0.05) q2st=0.01;
  if(q2st>0.9)  q2st=0.9;
  fitter->SetParameter(2,"Q2",am*q2st,.1,0.,2.*am);//Q2
  double q3st = 1.9*rat;
  if(q3st<0.05) q3st=0.01;
  if(q3st>0.9)  q3st=0.9;
  fitter->SetParameter(3,"Q3",am*q3st,.1,0.,2.*am);//Q3

  //printf(" start Q %f %f %f |  %f %f %f | rat %f tail %f\n",am*q0st,am*q2st,am*q3st,q0st,q2st,q3st,rat,tail);

  //   fitter->SetParameter(2,"Q2",0.5*am,.1,0.,2.*am);//Q2
  //   fitter->SetParameter(3,"Q3",0.25*am,.1,0.,2.*am);//Q3
  //   if(fabs(tail-rc)<30){
  //   printf("...........fabs(tail-rc)<30\n");
  //   fitter->SetParameter(2,"Q2",0.05*am,.1,0.,2.*am);//Q2
  //   fitter->SetParameter(3,"Q3",0.,.1,0.,2.*am);//Q3
  //   fitter->FixParameter(3);//Q3
  //   }


  fitter->SetParameter(4,"x0",t00,.1,-n_sample/10, n_sample);//x0

  //cout << "DEBUG YYYYYYYY ty0 = "<< ty0 << endl;
  fitter->SetParameter(5,"t0",ty0,.001,.01,30);//T0
  //fitter->FixParameter(5);//T0

  fitter->SetParameter(6,"t1",25,.001,90,350);//T1
  fitter->FixParameter(6);//T1

  fitter->SetParameter(7,"t2",t2,.01,40, 100);//T1 -> decay1
  //if(MOD==30&& isig==1)fitter->SetParameter(7,"t2",t2/2.,.01,20, 100);//T1 -> decay1


  fitter->SetParameter(8,"t3",t3,.01,150, 1000);//T3 -> decay2
  fitter->FixParameter(8);//T3 -> decay2

  fitter->SetParameter(9,"t9",t9,.01,.1, 30);//tr -> common si rise time
  //if(MOD==30&& isig==1)fitter->SetParameter(9,"t9",t9/2.,.01,.1, 30);//tr -> common si rise time
  //fitter->FixParameter(9);//tr -> common si rise time

  fitter->SetParameter(10,"t10",t10,0.01,0.1,30);//tr -> common csi rise time
  //if(MOD==30&& isig==1)fitter->SetParameter(10,"t10",t10/2.,0.01,0.1,30);//tr -> common csi rise time
  fitter->FixParameter(10);//tr -> common csi rise time

  //cout << "DEBUG XXXXXXXXXXX rat = " << rat << endl;
  if(isig==0){
    fitter->SetParameter(5,"t0",ty0,.001,.01,30);//T0
    fitter->SetParameter(6,"t1",30,.001,.01,2000);//T1
    fitter->FixParameter(6);//T1

    fitter->SetParameter(2,"Q2",0.,.1,0.,2.*am);//Q2
    fitter->FixParameter(2);//rc
    fitter->SetParameter(3,"Q3",0.,.1,0.,2.*am);//Q3
    fitter->FixParameter(3);//rc
    fitter->SetParameter(7,"t2",10,.01,40, 250);//T1 -> decay1
    fitter->FixParameter(7);//rc
    fitter->SetParameter(8,"t3",11,.01,250, 2000);//T3 -> decay2
    fitter->FixParameter(8);//rc
  }
  else{
    fitter->SetParameter(1,"Q1",0.     ,.1,0.,2.*am);//Q1
    fitter->FixParameter(1);//Q1
    fitter->SetParameter(5,"t0",ty0,.001,.01,30);//T0
    fitter->SetParameter(9,"t9",t9,.001,.1, 50);//tr -> common si rise time
    fitter->FixParameter(9);//Q0


    if(rat>0.05){
      fitter->FixParameter(5);//T0
    }
    //   if(MOD==30 || MOD==3 || MOD==9 || MOD==26) fitter->FixParameter(5);//T0
    if(mod==30) fitter->FixParameter(5);//T0


    //   else{
    //     printf("rat %f\n",rat);
    //     //plot=1;
    //   }

    //   if(ratexp<rat){
    //     printf("rat>ratexp %f %f\n",rat,ratexp);
    //     fitter->FixParameter(5);//T0
    //   }
    //   else{
    //     printf("rat<ratexp %f %f\n",rat,ratexp);
    //     //plot=1;
    //   }
  }

  fitter->SetParameter(11,"rc",rc1,0.1,0.1,3000);//tr -> common csi rise time
  fitter->FixParameter(11);//rc


  arglist[0] = 100;
  arglist[1] = 0.001;

  fitter->ExecuteCommand("MIGRAD", arglist, 2);
  //fitter->Minimize();

  //fitter->PrintResults(4, 0);

  Double_t eplus, eminus, eparab, globcc;
  Double_t value, verr, vlow, vhigh;
  Double_t ppar[12], epar[12];
  char name[80];
  for(int i=0;i<12;i++){
    fitter->GetParameter(i, &name[0], value, verr, vlow, vhigh);
    fitter->GetErrors(i, eplus, eminus, eparab, globcc);
    ppar[i]=value;
    epar[i]=eparab;
    //     printf("%d %s %f %f %f %f %f\n",i,name,value,eplus, eminus, eparab, globcc);
    funn->SetParameter(i,value);
  }

  //invert par[5] par[6]
  if(funn->GetParameter(5)>funn->GetParameter(6)){
    double ddum = funn->GetParameter(6);
    funn->SetParameter(6,funn->GetParameter(5));
    funn->SetParameter(5,ddum);
    ddum = funn->GetParameter(1);
    funn->SetParameter(1,funn->GetParameter(0));
    funn->SetParameter(0,ddum);
    funn->Update();
  }




  funn0->SetParameter(0,funn->GetParameter(0));
  funn0->SetParameter(1,funn->GetParameter(4));
  funn0->SetParameter(2,funn->GetParameter(5));
  funn0->SetParameter(3,funn->GetParameter(9));
  funn0->SetParameter(4,funn->GetParameter(11));
  funn1->SetParameter(0,funn->GetParameter(1));
  funn1->SetParameter(1,funn->GetParameter(4));
  funn1->SetParameter(2,funn->GetParameter(6));
  funn1->SetParameter(3,funn->GetParameter(9));
  funn1->SetParameter(4,funn->GetParameter(11));
  funn2->SetParameter(0,funn->GetParameter(2));
  funn2->SetParameter(1,funn->GetParameter(4));
  funn2->SetParameter(2,funn->GetParameter(7));
  funn2->SetParameter(3,funn->GetParameter(10));
  funn2->SetParameter(4,funn->GetParameter(11));
  funn3->SetParameter(0,funn->GetParameter(3));
  funn3->SetParameter(1,funn->GetParameter(4));
  funn3->SetParameter(2,funn->GetParameter(8));
  funn3->SetParameter(3,funn->GetParameter(10));
  funn3->SetParameter(4,funn->GetParameter(11));

  //printf("XX %7.4f %7.4f %7.4f\n",funn->GetParameter(5),funn->GetParameter(9),funn->GetParameter(0));


  int npt = 0;
  double DIFF2=0;
  double DIFF =0;
  for(int i=0; i<n_sample; i++){
    double xi = double(i)+0.5;
    double yhi = hfit->GetBinContent(i+1);
    double yfi = funn->Eval(xi);
    double ei  = fabs(yhi);
    if(ei){
      DIFF += pow(yhi-yfi,2)/ei;
      npt++;
    }
  }
  DIFF2 = DIFF/(npt-5);

  //DIFF2/=(n_sample-8);
  //printf("DIFF3 %f %f\n",funn->GetParameter(4),DIFF2);

  Int_t plot = 0;
  if(plot){
    printf("\n%d\n",fMbsInfo->GetCamacEventNumber());
    printf("%8d %2d %d  ", fEventCounter ,  mod, isig);
    // printf("%7.1f %7.1f %7.1f %7.1f %7.1f | %7.1f %7.1f | %7.1f %7.1f | %7.1f  %7.1f %7.1f | %7.3f\n", 
    printf("%.3f %.3f %.3f %.3f %.3f | %.3f %.3f | %.3f %.3f | %.3f  %.3f %.3f | %7.3f\n", 
    funn->GetParameter(0),
    funn->GetParameter(1),
    funn->GetParameter(2),
    funn->GetParameter(3),
    funn->GetParameter(4),
    funn->GetParameter(5),
    funn->GetParameter(6),
    funn->GetParameter(7),
    funn->GetParameter(8),
    funn->GetParameter(9),
    funn->GetParameter(10),
    //funn->GetParameter(11),
    tail,
    // funn->GetChisquare()/funn->GetNDF()
    //sqrt(DIFF2)
    DIFF2
    );
    getchar();
  }
  ////// MONDAY 2012-11-26

  //printf("%d %d %d\n",indx[0],indx[1],indx[2]);
  //funn->Update();
  //funn->EvalPar(xx);


  parRCGAM12SI0[0]=-1;
  parRCGAM12SI0[1]=-1;
  parRCGAM12SI0[2]=-1;
  parRCGAM12SI0[3]=-1;
  parRCGAM12SI0[4]=-1;
  parRCGAM12SI1[0]=-1;
  parRCGAM12SI1[1]=-1;
  parRCGAM12SI1[2]=-1;
  parRCGAM12SI1[3]=-1;
  parRCGAM12SI1[4]=-1;
  parRCGAM12SI2[0]=-1;
  parRCGAM12SI2[1]=-1;
  parRCGAM12SI2[2]=-1;
  parRCGAM12SI2[3]=-1;
  parRCGAM12SI2[4]=-1;
  parRCGAM12SI3[0]=-1;
  parRCGAM12SI3[1]=-1;
  parRCGAM12SI3[2]=-1;
  parRCGAM12SI3[3]=-1;
  parRCGAM12SI3[4]=-1;


  int it00=int(t00);


  parRCGAM12SI[0]=funn->GetParameter(0);
  parRCGAM12SI[1]=funn->GetParameter(1);
  parRCGAM12SI[2]=funn->GetParameter(2);
  parRCGAM12SI[3]=funn->GetParameter(3);
  parRCGAM12SI[4]=funn->GetParameter(4);
  parRCGAM12SI[5]=funn->GetParameter(5);
  parRCGAM12SI[6]=funn->GetParameter(6);
  parRCGAM12SI[7]=funn->GetParameter(7);
  parRCGAM12SI[8]=funn->GetParameter(8);
  parRCGAM12SI[9]=funn->GetParameter(9);
  parRCGAM12SI[10]=funn->GetParameter(10);
  parRCGAM12SI[11]=funn->GetParameter(11);
  parRCGAM12SI0[0]=funn->GetParameter(0);
  parRCGAM12SI1[0]=funn->GetParameter(1);
  parRCGAM12SI2[0]=funn->GetParameter(2);
  parRCGAM12SI3[0]=funn->GetParameter(3);
  parRCGAM12SI0[1]=funn->GetParameter(4);
  parRCGAM12SI1[1]=funn->GetParameter(4);
  parRCGAM12SI2[1]=funn->GetParameter(4);
  parRCGAM12SI3[1]=funn->GetParameter(4);

  p0=funn->GetParameter(0);
  p1=funn->GetParameter(1);
  p2=funn->GetParameter(2);
  p3=funn->GetParameter(3);
  p4=funn->GetParameter(4);
  p5=funn->GetParameter(5);
  p6=funn->GetParameter(6);
  p7=funn->GetParameter(7);
  p8=funn->GetParameter(8);
  p9=funn->GetParameter(9);
  p10=funn->GetParameter(10);
  p11=funn->GetParameter(11);
  //chi2=funn->GetChisquare()/funn->GetNDF();
  chi2=DIFF2;



  //interpolation and moments
  // POL3
  int bmax = hfit->GetMaximumBin();

  // int bl  = (bmax-p4)/4.; //bins left
  // int br  = (bmax-p4)/2.; //bins right
  int bl  = (bmax-(int)p4)/3.; //bins left
  int br  = (bmax-(int)p4)/1.; //bins right
  if(bmax-(int)p4<150) {bl*=2;br*=0.75;}
  // int bl  =  50; //bins left
  // int br  =  70; //bins right
  double  SX6,SX5,SX4,SX3,SX2,SX1,SX0,SYX3,SYX2,SYX1,SYX0;
  double  XI,YI,X0,X1,X2,X3,X4,X5,X6;
  double  a,b,c,d,xmax,ymax;
  SX6=SX5=SX4=SX3=SX2=SX1=SX0=SYX3=SYX2=SYX1=SYX0=0;
  int bl0 = bmax-bl;
  if(bl0<1) bl0=1;
  int br0 = bmax+br;
  if(br0>hfit->GetNbinsX()) br0=hfit->GetNbinsX();
  for(int i=bl0;i<br0;i++){
    XI=hfit->GetBinCenter(i);
    YI=hfit->GetBinContent(i);
    X0=1.;
    X1=XI;
    X2=X1*XI;
    X3=X2*XI;
    X4=X3*XI;
    X5=X4*XI;
    X6=X5*XI;

    SX6+=X6;
    SX5+=X5;
    SX4+=X4;
    SX3+=X3;
    SX2+=X2;
    SX1+=X1;
    SX0+=X0;
    SYX3+=YI*X3;
    SYX2+=YI*X2;
    SYX1+=YI*X1;
    SYX0+=YI;
  }
  const int npar = 4;
  double DF[npar+1][npar];
  DF[0][0]=SX6; DF[1][0]=SX5; DF[2][0]=SX4; DF[3][0]=SX3; DF[4][0]=-SYX3;
  DF[0][1]=SX5; DF[1][1]=SX4; DF[2][1]=SX3; DF[3][1]=SX2; DF[4][1]=-SYX2;
  DF[0][2]=SX4; DF[1][2]=SX3; DF[2][2]=SX2; DF[3][2]=SX1; DF[4][2]=-SYX1;
  DF[0][3]=SX3; DF[1][3]=SX2; DF[2][3]=SX1; DF[3][3]=SX0; DF[4][3]=-SYX0;

  //GAUSS ELIMINATION
  double ROW[npar+1];
  for(int iy0=1;iy0<npar;iy0++){
    for(int iy=iy0;iy<npar;iy++){
      for(int ix=iy0-1;ix<npar+1;ix++) ROW[ix] = DF[ix][iy0-1]*DF[iy0-1][iy]/DF[iy0-1][iy0-1];
      for(int ix=iy0-1;ix<npar+1;ix++){
        DF[ix][iy] -= ROW[ix];
      }
    }
  }


  //BACK SUBSTITUTION
  double dx[npar]={0};
  for(int iy=npar-1;iy>=0;iy--){
    dx[iy] = -DF[npar][iy];
    for(int ix=iy+1;ix<npar;ix++){
      dx[iy] -= dx[ix]*DF[ix][iy];
    }
    dx[iy]/=DF[iy][iy];
  }

  a  = dx[0];
  b  = dx[1];
  c  = dx[2];
  d  = dx[3];

  double del = 4.*b*b-12.*a*c;
  xmax = bmax;
  ymax = hfit->GetMaximum();
  if(del>=0 && a!=0){
    double sdel = sqrt(del);
    double x1 = (-2.*b-sdel)/(6.*a);
    double x2 = (-2.*b+sdel)/(6.*a);
    //printf("x1 x2  %f %f \n",x1,x2);
    xmax = (fabs(x1-bmax)<fabs(x2-bmax)) ? x1 : x2;
    ymax = a*xmax*xmax*xmax+b*xmax*xmax+c*xmax+d;
  }
  //printf("max %f %f %f %f\n",a,b,c,d);
  //printf("max %f %f\n",xmax,ymax);

  ///(*iter)->sam[isig]  = ymax; ///Not necessary here
  ///(*iter)->sat[isig]  = xmax; ///Not necessary here


  vector <Double_t > vect = getPercents(hfit,p4,xmax,ymax);

  ///if(vect.size()){
  ///(*iter)->s1[isig]   = vect[ 0];
  ///(*iter)->s2[isig]   = vect[ 2];
  ///(*iter)->s3[isig]   = vect[ 4];
  ///(*iter)->s4[isig]   = vect[ 6];
  ///(*iter)->s5[isig]   = vect[ 8];
  ///(*iter)->s6[isig]   = vect[10];
  ///(*iter)->s7[isig]   = vect[12];
  ///(*iter)->s8[isig]   = vect[14];
  ///(*iter)->s9[isig]   = vect[16];
  ///(*iter)->s9f[isig]  = vect[18];
  ///}
  //  printf("p4,xmax,ymax : %.3f %.3f %.3f\n",p4,xmax,ymax); 

  
  /// FALLING EDGE PROOBING: Wyliczenie probkowania sygnalu na zboczu narastajacym 
  vector <Double_t > vectSamplesFromTheFallingEdge;
  vectSamplesFromTheFallingEdge = GetSamplesFromTheFallingEdge(hfit,p4,xmax,ymax );

  //moments 
  const int NMT = 48;
  double NNN[NMT]={0,1,2,3,4,5,6,7,8,9,10,12,15,19,24,30,37,45,54,64,75,87,100,
  114,129,145,162,180,200,220,240,260,280,300,320,340,370,400,450,500,550,600,
  650,700,750,800,850,900};
  double YYY[NMT];
  double MT[NMT]={0};
  it01 = int(p4+0.5); ///removed here int (doubled declaration)
  // int nbmax = int(rc+0.5);
  int nbmax = 800;
  double TMAX = double(nbmax);
  if(TMAX<=0) TMAX=1.;
  for(int i=it01; i<it01+nbmax+1; i++){
    //double xi = hfit->GetBinCenter(i+1)-hfit->GetBinCenter(it01+1);
    double xi = double(i-it01);
    double yi = 0;
    if(i<n_sample){
      yi = hfit->GetBinContent(i+1); // take measured values
    }
    else{
      yi = exp(-(double(i)+0.5)/tail+btail); // extrapolate using exp fit
    }
    xi /= TMAX;
    for(int in=0;in<NMT;in++){
      MT[in] += yi*pow(xi,NNN[in]);
    }
  }

  // printf("00, %f %6.3f,\n",MT[ 0],1.);
  //printf("  0, %f %7.3f %7.3f,\n",1.,1.,0.);
  if(MT[0]){
    for(int in=1;in<NMT;in++){
      MT[in] /= MT[0];
      //  printf("%3d, %f %7.3f %7.3f,\n",int(NNN[in]),
      //  MT[in],1/MT[in],1/MT[in]-(1+NNN[in]));
      if(MT[in]){
        YYY[in]=1./MT[in]-(1+NNN[in]);
        //  YYY[in]=1/MT[in];
      }
    }
  }
  //MT[0] = 1;

  ///(*iter)->m0[isig]   = MT[0];
  ///(*iter)->m1[isig]   = MT[1];
  ///(*iter)->m2[isig]   = MT[2];
  ///(*iter)->m3[isig]   = MT[3];
  ///(*iter)->m4[isig]   = MT[4];
  ///(*iter)->m5[isig]   = MT[5];
  ///(*iter)->m6[isig]   = MT[6];
  ///(*iter)->m7[isig]   = MT[7];
  ///(*iter)->m8[isig]   = MT[8];
  ///(*iter)->m9[isig]   = MT[9];



  double ymx=1;
  double xmx=0;
  int    imx=0;
  for(int i=0;i<NMT;i++){
    if(YYY[i]>ymx){
      ymx = YYY[i];
      xmx = NNN[i];
      imx = i;
    }
  }
  int bl1  = imx-3; //bins left
  int br1  = imx+3; //bins right
  if(bl1<0)   bl1=0;
  if(br1>NMT) br1=NMT;
  double  A,B,C,D,E,F,G,H;
  double  XXI,YYI,XX0,XX1,XX2,XX3,XX4,XX5,XX6;
  double  W,Wa,Wb,Wc,a1,b1,c1,xmax1,ymax1;
  A=B=C=D=E=F=G=H=0;
  a1=b1=c1=0;
  xmax1=ymax1=0;
  for(int i=bl1;i<br1;i++){
    XXI=NNN[i];
    YYI=YYY[i];
    XX0=1.;
    XX1=XXI;
    XX2=XX1*XXI;
    XX3=XX2*XXI;
    XX4=XX3*XXI;

    A+=XX4;
    B+=XX3;
    C+=XX2;
    D+=XX1;
    E+=YYI*XX2;
    F+=YYI*XX1;
    G+=YYI;
    H+=XX0;
  }
  W  = A*C*H+2*B*D*C-C*C*C-A*D*D-H*B*B;
  Wa = E*C*H+F*D*C+B*D*G-G*C*C-D*D*E-F*B*H;
  Wb = A*F*H+B*G*C+E*D*C-F*C*C-A*G*D-H*B*E;
  Wc = A*C*G+B*D*E+B*F*C-E*C*C-A*D*F-G*B*B;
  //printf("W %f %f %f %f\n",W,Wa,Wb,Wc);

  if(W){ 
  a1  = Wa/W;
  b1  = Wb/W;
  c1  = Wc/W;
  } 
  //double Sa = sqrt((C*H-D*D)/W);
  //double Sb = sqrt((A*H-C*C)/W);
  //double Sc = sqrt((A*C-B*B)/W);

  if(a1){
    xmax1 = -b1/(2.*a1);
    ymax1 = a1*xmax1*xmax1+b1*xmax1+c1;
    //  printf("max1 %f %f %f\n",a1,b1,c1);
    //  printf("max1 %f %f\n",xmax1,ymax1);
  }

  ///(*iter)->mx[isig]   = xmax1;
  ///(*iter)->my[isig]   = ymax1;
  ///(*iter)->mz[isig]   = c1;

  double  a0,b0,c0;
  A=B=C=D=E=F=G=H=0;
  a0=b0=c0=0;
  for(int i=0;i<8;i++){
    XXI=NNN[i];
    YYI=YYY[i];

    XX0=1.;
    XX1=XXI;
    XX2=XX1*XXI;
    XX3=XX2*XXI;
    XX4=XX3*XXI;
    XX5=XX4*XXI;
    XX6=XX5*XXI;

    A+=XX6;
    B+=XX5;
    C+=XX4;
    D+=XX3;
    H+=XX2;
    E+=YYI*XX3;
    F+=YYI*XX2;
    G+=YYI*XX1;
  }
  W  = A*C*H+2*B*D*C-C*C*C-A*D*D-H*B*B;
  Wa = E*C*H+F*D*C+B*D*G-G*C*C-D*D*E-F*B*H;
  Wb = A*F*H+B*G*C+E*D*C-F*C*C-A*G*D-H*B*E;
  Wc = A*C*G+B*D*E+B*F*C-E*C*C-A*D*F-G*B*B;
  //printf("W %f %f %f %f\n",W,Wa,Wb,Wc);

  if(W){ 
    a0  = Wa/W;
    b0  = Wb/W;
    c0  = Wc/W;
  } 

  ///(*iter)->ma[isig]   = a0; ///Here it is not necessary
  ///(*iter)->mb[isig]   = b0; ///Here it is not necessary
  ///(*iter)->mc[isig]   = c0; ///Here it is not necessary

  //printf("a0 b0 c0 %f %f %f\n",a0,b0,c0); 
  //interpolation and moments


  fPeakData->sam[isig]  = ymax;
  fPeakData->sat[isig]  = xmax;

  fPeakData->pds[isig]  = pds;
  fPeakData->pdr[isig]  = pdsrms;
  fPeakData->tm[isig]   = t00;
  fPeakData->am[isig]   = am;
  
  ///TODO remove those commented lines:
  ///ppd->pds[isig]  = pds;
  ///ppd->pdr[isig]  = pdsrms;
  ///ppd->tm[isig]   = t00;
  ///ppd->evt[isig]  = ievt;
    
  ///am=-1.e10;
  ///for(int i=0; i<n_sample; i++){
  /// ppd->dfadc[isig][i] = double(ifadc[i])-pds;
  /// if(ppd->dfadc[isig][i]>am) am=ppd->dfadc[isig][i];
  ///}
  ///ppd->am[isig]  = am;

  fPeakData->s1[isig]   = vect[ 0];
  fPeakData->s2[isig]   = vect[ 2];
  fPeakData->s3[isig]   = vect[ 4];
  fPeakData->s4[isig]   = vect[ 6];
  fPeakData->s5[isig]   = vect[ 8];
  fPeakData->s6[isig]   = vect[10];
  fPeakData->s7[isig]   = vect[12];
  fPeakData->s8[isig]   = vect[14];
  fPeakData->s9[isig]   = vect[16];
  fPeakData->s9f[isig]  = vect[18];

  if( vectSamplesFromTheFallingEdge.size() == 20){ /// because not always is possible proobe the falling slope
    fPeakData->sf0[isig]   = vectSamplesFromTheFallingEdge[ 1]; ///vectSamplesFromTheFallingEdge=[x10,y10,x20,y20,...]
    fPeakData->sf1[isig]   = vectSamplesFromTheFallingEdge[ 3];
    fPeakData->sf2[isig]   = vectSamplesFromTheFallingEdge[ 5];
    fPeakData->sf3[isig]   = vectSamplesFromTheFallingEdge[ 7];
    fPeakData->sf4[isig]   = vectSamplesFromTheFallingEdge[ 9];
    fPeakData->sf5[isig]   = vectSamplesFromTheFallingEdge[11];
    fPeakData->sf6[isig]   = vectSamplesFromTheFallingEdge[13];
    fPeakData->sf7[isig]   = vectSamplesFromTheFallingEdge[15];
    fPeakData->sf8[isig]   = vectSamplesFromTheFallingEdge[17];
    fPeakData->sf9[isig]   = vectSamplesFromTheFallingEdge[19];
  }/// else: members will remain empty
  
  fPeakData->m0[isig]   = MT[0];
  fPeakData->m1[isig]   = MT[1];
  fPeakData->m2[isig]   = MT[2];
  fPeakData->m3[isig]   = MT[3];
  fPeakData->m4[isig]   = MT[4];
  fPeakData->m5[isig]   = MT[5];
  fPeakData->m6[isig]   = MT[6];
  fPeakData->m7[isig]   = MT[7];
  fPeakData->m8[isig]   = MT[8];
  fPeakData->m9[isig]   = MT[9];

  fPeakData->mx[isig]   = xmax1;
  fPeakData->my[isig]   = ymax1;
  fPeakData->mz[isig]   = c1;

  fPeakData->ma[isig]   = a0;
  fPeakData->mb[isig]   = b0;
  fPeakData->mc[isig]   = c0;

  fPeakData->tls[isig] = tls;
  fPeakData->tlo[isig] = tlo;
  fPeakData->at[isig]  = hfit->GetMaximumBin();

  fPeakData->p0[isig]  = p0;
  fPeakData->p1[isig]  = p1;
  fPeakData->p2[isig]  = p2;
  fPeakData->p3[isig]  = p3;
  fPeakData->p4[isig]  = p4;
  fPeakData->p5[isig]  = p5;
  fPeakData->p6[isig]  = p6;
  fPeakData->p7[isig]  = p7;
  fPeakData->p8[isig]  = p8;
  fPeakData->p9[isig]  = p9;
  fPeakData->p10[isig] = p10;
  fPeakData->rc[isig]  = p11;
  fPeakData->chi2[isig]  = chi2;

  return;
} /// End of the function GetPulseParameters()

//______________________________________________________________________
void TKratAnaStep1::Reset(){
 
}   

//______________________________________________________________________
void TKratAnaStep1::Finish()
{
   if ( fVerbose > 2 ) cout << "[TKratAnaStep1::Finish()] fEventCounter=" << fEventCounter << endl;
}

	      

ClassImp(TKratAnaStep1)