#include "CbmTrdModuleRecT.h"
#include "CbmTrdDigi.h"
#include "CbmTrdCluster.h"
#include "CbmTrdHit.h"
#include "CbmTrdParModDigi.h"
#include "CbmTrdFASP.h"

#include <FairLogger.h>

#include <TGeoPhysicalNode.h>

#include <TClonesArray.h>
#include <TGraphErrors.h>
#include <TF1.h>

#include <iostream>
#include <vector>

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

Double_t CbmTrdModuleRecT::fgDT[] = {4.181e-6, 1586, 24};
TGraphErrors* CbmTrdModuleRecT::fgEdep = NULL;
TGraphErrors* CbmTrdModuleRecT::fgT = NULL;
TF1* CbmTrdModuleRecT::fgPRF = NULL;
//_______________________________________________________________________________
CbmTrdModuleRecT::CbmTrdModuleRecT()
  : CbmTrdModuleRec()
  ,fConfigMap(0)
  ,fT0(0)
  ,fBuffer()
  ,vs()
  ,vse()
  ,vt()
  ,vx()
  ,vxe()
{
    
}

//_______________________________________________________________________________
CbmTrdModuleRecT::CbmTrdModuleRecT(Int_t mod, Int_t ly, Int_t rot)
  : CbmTrdModuleRec(mod, ly, rot)
  ,fConfigMap(0)
  ,fT0(0)
  ,fBuffer()
  ,vs()
  ,vse()
  ,vt()
  ,vx()
  ,vxe()
{
  //printf("AddModuleT @ %d\n", mod); Config(1,0);
  SetNameTitle(Form("TrdRecT%d", mod), "Reconstructor for triangular read-out.");
}

//_______________________________________________________________________________
CbmTrdModuleRecT::~CbmTrdModuleRecT()
{
}

//_______________________________________________________________________________
Bool_t CbmTrdModuleRecT::AddDigi(const CbmTrdDigi *d, Int_t id)
{
/* Add digi to cluster fragments. At first clusters are ordered on pad rows and time. 
 * No channel ordering is assumed. The time condition for a digi to enter a cluster 
 * chunk is to have abs(dt)<5 wrt cluster t0 
 */

  if(CWRITE()) cout<<"add @"<<id<<" "<<d->ToString();

  Int_t ch  = d->GetAddressChannel(), col,
        row = GetPadRowCol(ch, col),      
        dtime;
  Double_t t, r = d->GetCharge(t, dtime);
  Int_t tm  = d->GetTimeDAQ()-fT0, terminator(0);
  if(dtime<0) tm+=dtime; // correct for the difference between tilt and rect
  if(r<1) terminator=1;
  else if(t<1) terminator=-1;
  
  if(CWRITE()) printf("row[%2d] col[%2d] tm[%2d] terminator[%d]\n", row, col, tm, terminator);
  CbmTrdCluster *cl(NULL);
  
  // get the link to saved clusters
  std::map<Int_t,std::list<CbmTrdCluster*>>::iterator it=fBuffer.find(row);
  // check for saved 
  if(it!=fBuffer.end()){
    Bool_t kINSERT(kFALSE);
    std::list<CbmTrdCluster*>::iterator itc=fBuffer[row].begin();
    for(; itc!=fBuffer[row].end(); itc++){
      if(CWRITE()) cout<<(*itc)->ToString();
      
      UShort_t tc = (*itc)->GetStartTime();
      Int_t dt = tc - tm;

      if(dt<-5) continue;
      else if(dt<5) {
        if(CWRITE()) printf("match time dt=%d\n", dt);
        if((*itc)->IsChannelInRange(ch)==0){ 
          if(!(*itc)->AddDigi(id, ch, terminator, dt)) break;
          kINSERT=kTRUE;
          if(CWRITE()) cout<<"          => Cl (Ad) "<<(*itc)->ToString();        
          break;
        } 
      } else {
        if(CWRITE()) printf("break for time dt=%d\n", dt);
        break;
      }
    }
    
    if(!kINSERT){
      if(itc != fBuffer[row].end() && itc != fBuffer[row].begin()){
        itc--;
        fBuffer[row].insert(itc, cl = new CbmTrdCluster(fModAddress, id, ch, row, tm));
        if(CWRITE()) cout<<"          => Cl (I) "<<cl->ToString();        
      } else {
        fBuffer[row].push_back(cl = new CbmTrdCluster(fModAddress, id, ch, row, tm));
        if(CWRITE()) cout<<"          => Cl (Pb) "<<cl->ToString();                
      }
      cl->SetTrianglePads();
      if(terminator<0) cl->SetProfileStart();
      else if(terminator>0) cl->SetProfileStop();
    }
  } else {
    fBuffer[row].push_back(cl = new CbmTrdCluster(fModAddress, id, ch, row, tm));
    cl->SetTrianglePads();
    if(terminator<0) cl->SetProfileStart();
    else if(terminator>0) cl->SetProfileStop();
    if(CWRITE()) cout<<"          => Cl (Nw) "<<cl->ToString();                
  }
  
  return kTRUE;
}

//_______________________________________________________________________________
Int_t CbmTrdModuleRecT::GetOverThreshold() const
{
  Int_t nch(0);
  for(std::map<Int_t, std::list<CbmTrdCluster*>>::const_iterator ir=fBuffer.cbegin(); ir!=fBuffer.cend(); ir++){
    for(std::list<CbmTrdCluster*>::const_iterator itc=(*ir).second.cbegin(); itc!=(*ir).second.cend(); itc++) nch+=(*itc)->GetNCols();
  }
  return nch;
}

//_______________________________________________________________________________
Int_t CbmTrdModuleRecT::FindClusters()
{
  CbmTrdCluster *cl(NULL);
  
  // get the link to saved clusters
  Int_t ncl(0);
  std::list<CbmTrdCluster*>::iterator itc0, itc1;
  for(std::map<Int_t, std::list<CbmTrdCluster*>>::iterator ir=fBuffer.begin(); ir!=fBuffer.end(); ir++){
    itc0=(*ir).second.begin();
    while((*ir).second.size()>1 && itc0!=(*ir).second.end()){ // try merge clusters
      itc1=itc0; itc1++;

      Bool_t kMERGE=kFALSE;
      while(itc1!=(*ir).second.end()){
        if(CWRITE()) cout<<"  base cl[0] : "<<(*itc0)->ToString()
                         <<"     + cl[1] : "<<(*itc1)->ToString();
        if((*itc0)->Merge((*itc1))){
          if(CWRITE()) cout<<"  SUM        : "<<(*itc0)->ToString();
          kMERGE = kTRUE;
          delete (*itc1); 
          itc1=(*ir).second.erase(itc1);
          if(itc1==(*ir).second.end()) break;
        } else itc1++;
      }
      if(!kMERGE) itc0++;
    }

    for(itc0=(*ir).second.begin(); itc0!=(*ir).second.end(); itc0++){ 
      cl = (*itc0);
      cl = new((*fClusters)[ncl++]) CbmTrdCluster(*cl);
      cl->SetTrianglePads();
      delete (*itc0); 
    }
  }
  fBuffer.clear();
  
  //printf("fClusters[%p] nCl[%d]\n", (void*)fClusters, fClusters->GetEntriesFast());
  //LOG(info) << GetName() << "::FindClusters : " << ncl;
  return ncl;
}

//_______________________________________________________________________________
Bool_t CbmTrdModuleRecT::MakeHits()
{
  return kTRUE;  
}

//_______________________________________________________________________________
Bool_t CbmTrdModuleRecT::Finalize()
{
/*  Steering routine for classifying hits and apply further analysis
 * -> hit deconvolution (see Deconvolute)
 * -> hit merging for row-cross (see RowCross)
 */

//   Int_t nhits=fHits->GetEntriesFast();
//   //if(CWRITE()) 
//     LOG(info) << "\n"<<GetName()<<"::Finalize("<<nhits<<")";
//   CbmTrdHit *hit(NULL), *hitp(NULL);
//   for(Int_t ih(0); ih<nhits; ih++){
//     hit = (CbmTrdHit*)((*fHits)[ih]); 
//     for(Int_t jh(ih+1); jh<nhits; jh++){
//       hitp = (CbmTrdHit*)((*fHits)[jh]); 
//       //if(CWRITE()) 
//       cout<<ih<<" "<<hit->ToString();
//       cout<<"-> "<<jh<<"  "<<hitp->ToString();
//     }
//   }  
  return kTRUE;
}

#include <TCanvas.h>
#include <TH1.h>
#define NANODE 9
//_______________________________________________________________________________
CbmTrdHit* CbmTrdModuleRecT::MakeHit(Int_t ic, const CbmTrdCluster *cl, std::vector<const CbmTrdDigi*> *digis)
{
  if(!fgEdep){ // first use initialization of PRF helppers
    LOG(info) << GetName()<<"::MakeHit: Init static helpers. ";
    fgEdep = new TGraphErrors;
    fgEdep->SetLineColor(kBlue);fgEdep->SetLineWidth(2);
    fgT = new TGraphErrors;
    fgT->SetLineColor(kBlue);fgT->SetLineWidth(2);
    fgPRF=new TF1("prf", "[0]*exp(-0.5*((x-[1])/[2])**2)", -10, 10);
    fgPRF->SetLineColor(kRed);
    fgPRF->SetParNames("E", "x", "prf");
  }
  TH1 *hf(NULL); TCanvas *cvs(NULL);
  if(CDRAW()){
    cvs = new TCanvas("c", "TRD Anode Hypothesis",10,600,1000,500); cvs->Divide(2,1,1.e-5,1.e-5);
    TVirtualPad *vp = cvs->cd(1);
    vp->SetLeftMargin(0.06904908);
    vp->SetRightMargin(0.009969325);
    vp->SetTopMargin(0.01402806);
    vp = cvs->cd(2);
    vp->SetLeftMargin(0.06904908);
    vp->SetRightMargin(0.009969325);
    vp->SetTopMargin(0.01402806);
    hf = new TH1I("hf", ";x [pw];s [ADC chs]", 100,-3, 3);
    hf->GetYaxis()->SetRangeUser(-50, 4500);
    hf->Draw("p");
  }

  if(CWRITE()) cout<<cl->ToString();
  
  // check cluster integrity and do digi merging if needed
  vector<Bool_t>      vmask(digis->size(), 0);  // masks in case of merged Digis
  vector<CbmTrdDigi*> vdgM; 
  if(cl->GetNCols() != digis->size() && !MergeDigis(digis, &vdgM, &vmask)){
    cout<<cl->ToString();  
    for(vector<const CbmTrdDigi*>::iterator i=digis->begin(); i!= digis->end(); i++) cout<<(*i)->ToString(); 
    return NULL;
  }
  
  // Read in all digis information()
  ULong64_t t0=fT0+cl->GetStartTime();        // absolute hit time (prompt signal)
  Int_t n0(0), ovf(0), cM;
  if(!(n0 = LoadDigis(digis, &vdgM, &vmask, t0, cM))){
    cout<<cl->ToString();  
    for(vector<const CbmTrdDigi*>::iterator i=digis->begin(); i!= digis->end(); i++) cout<<(*i)->ToString(); 
    return NULL;
  }
  if(n0<0){ovf=1; n0*=-1;}
  
  // analyse digis topology; no of signal types, maximum, etc
  Bool_t tM(kTRUE);           // maximum type tilt=1; rect=0
  Int_t nL(0);                // signal index for the max signal
  Int_t col, row= GetPadRowCol(cl->GetStartCh(), col);
  Double_t  max(0.),          // maximum signal
            LS(0.),           // left side sum of signals
            S(0.);            // sum of signals
  Int_t nr(0), nt(0);
  for(Int_t is(1); is<=n0; is++){ 
    if(vxe[is]>0) {       // select tilted coupling
      nt++; S+=vs[is]; 
      if(vs[is]>max){
        max=vs[is]; tM=kTRUE; nL=is; LS+=vs[is];
      }
    } else          {       // select rectangular coupling
      nr++; S+=vs[is];
      if(vs[is]>max){
        max=vs[is]; tM=kFALSE; nL=is; LS+=vs[is]; 
      }
    } 
  }
  col+=cM; S-=LS; LS-=max;
  // evaluate asymmetry 
  Int_t  lr(0),           // max signal left-right asymmetry wrt central pad
         tr(0),           // left-right innequality for symmetric clusters
         nR(n0+1-nL);     // no of signals to the right of maximum 
  if(nL<nR) lr=1;
  else if(nL>nR) lr=-1;
  if(!lr&&(n0%2)) tr=(LS<S?-1:1);
  // compute x and y offset from center pad
  Double_t dx(0.), dy(0.),
           edx(0.21650635), edy(0.77942286); // fixed error parametrization 
  switch(n0){
    case 1: 
      if(nt)  {dx = -0.5; dy=0;}    // T
      else    {dx =  0.5; dy=0;}    // R
      break;
    case 2:
      if(cl->HasOpenStart() && 
         cl->HasOpenStop()) { dx = cM?-1.:0.; dy = -0.5; }  // RT
      else                  { dx = 0.;        dy =  0.5; }  // TR
      break;
    case 3: 
      if(tM&&lr)        {dx = cM?-1.:0.; dy = GetYoffset(n0); }    // TRT asymm
      else if(!tM&&!lr) {dx = 0.;        dy = GetYoffset(n0);}     // TRT symm
      else if(tM&&!lr)  {dx = GetXoffset(n0); dy = 0.;}            // RTR symm
      else if(!tM&&lr)  {dx = GetXoffset(n0); dy = cM?0.5:-0.5;}   // RTR asymm
      break;
    default: 
      dx = GetXoffset(n0); dy = GetYoffset(n0);
      break;
  }
  if(dx<-0.5 && cM>0) { // shift graph representation to fit dx[pw] in [-0.5, 0.5]
    dx += 1.; col-=1; 
    for(UInt_t idx(0); idx<vx.size(); idx++) vx[idx]+=1;
  }
  if(dx>0.5) { // dirty fix for compound clusters TODO
    Int_t ishift = Int_t(dx-0.5)+1;
    dx -= ishift; col+=ishift; 
    for(UInt_t idx(0); idx<vx.size(); idx++) vx[idx]-=ishift;
  }
  dy = dx - dy; // only on natural scalling !
  // go to cm scale
  dx*=fDigiPar->GetPadSizeX(0);
  dy*=fDigiPar->GetPadSizeY(0);

  // apply systematic correction for x (MC derived)
  Int_t typ=0;  // [0] center hit type
                // [1] side hit type
  if((n0==5&&((tM&&!lr)||(!tM&&lr))) || // TRTRT symm/asymm 
      n0==4 ||
     (n0==3&&((tM&&!lr)||(!tM&&lr!=0)))) typ=1; // RTR symm/asymm
  Double_t xcorr(0.);
  Int_t nbins((NBINSCORRX-1)>>1),
        ii = nbins+TMath::Nint(dx/fgCorrXdx),
        i0, i1;
  if(ii<0||ii>=NBINSCORRX) LOG(warn)<<GetName() << "::MakeHit : Idx "<<ii<<" outside range for displacement "<<dx<<".";
  else{
    if(dx<fgCorrXdx*ii) {i0=TMath::Max(0,ii-1); i1=ii;}
    else                {i0=ii; i1=TMath::Min(NBINSCORRX-1,ii+1);}
    Double_t DDx = (fgCorrXval[typ][i1] - fgCorrXval[typ][i0]),
             a   = DDx/fgCorrXdx,  
             b   = fgCorrXval[typ][i0] - DDx*(i0-nbins);
    xcorr = 0.1*(b+a*dx);
  }
  dx+=xcorr; dy+=xcorr; 
  if(dx>0.5*fDigiPar->GetPadSizeX(0)) dx = 0.5*fDigiPar->GetPadSizeX(0); 
  else if(dx<-0.5*fDigiPar->GetPadSizeX(0)) dx = -0.5*fDigiPar->GetPadSizeX(0); 
  
  if(dy>0.5 * fDigiPar->GetPadSizeY(0)) { // 
    //printf("BEFORE dy[%+6.4f] dx[%+6.4f] {n[%d] max[%c] lr[%+d]}\n", dy, dx, n0, tM?'T':'R', lr);
    dy-=fDigiPar->GetPadSizeY(0);
  }
  if(dy<-0.5 * fDigiPar->GetPadSizeY(0)) { // 
    //printf("(BEFORE) dy[%+6.4f] dx[%+6.4f] {n[%d] max[%c] lr[%+d]}\n", dy, dx, n0, tM?'T':'R', lr);
    dy+=fDigiPar->GetPadSizeY(0);
  }
  
  // process y offset
  // apply systematic correction for y (MC derived)
  Float_t fdy(1.), yoff(0.);
  switch(n0){
    case 3:
      fdy = fgCorrYval[0][0]; yoff = fgCorrYval[0][1]; 
      if(tM&&!lr) dy-=tr*0.5*fDigiPar->GetPadSizeY(0);
      else if(lr)  dy-= 0.5*fDigiPar->GetPadSizeY(0);;
      break;
    case 4:
      fdy = fgCorrYval[1][0]; yoff = fgCorrYval[1][1]; 
      if((!tM&&lr==1)||(tM&&lr==-1)) yoff*=-1;
      break;
    case 5:
    case 7:
    case 9:
    case 11:
      fdy = fgCorrYval[2][0]; yoff = fgCorrYval[2][1]; 
      break;
    case 6:
    case 8:
    case 10:
      fdy = fgCorrYval[3][0]; yoff = fgCorrYval[3][1]; 
      if((!tM&&lr==1)||(tM&&lr==-1)) yoff*=-1;
      break;
  }
  dy*=fdy; dy+=yoff; 
  if(dy>0.5 * fDigiPar->GetPadSizeY(0)) dy = 0.5 * fDigiPar->GetPadSizeY(0); 
  else if(dy<-0.5 * fDigiPar->GetPadSizeY(0)) dy = -0.5 * fDigiPar->GetPadSizeY(0); 
  
  // get anode wire offset
  Int_t ia(0); Float_t ya; //  anode position in local pad coordinates
  for(; ia<=NANODE; ia++) {
    ya=-1.35+ia*0.3;
    if(dy > ya + 0.15) continue;
    break;
  }
  
  // Error parametrization X : parabolic model on cl size
  Double_t parX[] = {0.713724, -0.318667, 0.0366036};
  Int_t nex = TMath::Min(n0, 7);
  edx = parX[0] + parX[1] * nex + parX[2] * nex * nex;
  Double_t parY[] = {0.0886413, 0., 0.0435141}; 
  edy = parY[0] + parY[2] * dy * dy;
  
  if(CWRITE()) {
    printf("row[%2d] col[%2d] sz[%d%c] M[%d%c] dx[mm]=%6.3f dy[mm]=%6.3f t0[clk]=%llu OVF[%c]\n", row, col, n0, (lr?(lr<0?'L':'R'):'C'), cM, (tM?'T':'R'), 10*dx, 10*dy, t0, (ovf?'y':'n'));
    for(UInt_t idx(0); idx<vs.size(); idx++){    
      printf("%2d%cdt[%2d] s[ADC]=%6.1f+-%6.1f x[pw]=%5.2f+-%5.2f\n", idx, (UInt_t(nL)==idx?'*':' '), vt[idx], vs[idx], vse[idx], vx[idx], vxe[idx]);
    }
  }
  
  // compute energy
  for(UInt_t idx(0); idx<vs.size(); idx++){
    if(vxe[idx]>0){
      fgEdep->SetPoint(idx, vx[idx]+dy/fDigiPar->GetPadSizeY(0), vs[idx]);
      fgEdep->SetPointError(idx, vxe[idx], vse[idx]);
    } else {
      fgEdep->SetPoint(idx, vx[idx], vs[idx]);
      fgEdep->SetPointError(idx, vxe[idx], vse[idx]);
    }
  }
  Double_t xc=vx[n0+2];
  for(Int_t ip(vs.size()); ip<fgEdep->GetN(); ip++){
    //fgEdep->RemovePoint(ip);
    xc+=0.5; fgEdep->SetPoint(ip, xc, 0); fgEdep->SetPointError(ip, 0., 300);
  }
  if(CWRITE()) fgEdep->Print();
  
  Double_t e(0.), chi(100),
           xlo(*vx.begin()), xhi(*vx.rbegin());
  fgPRF->SetParameter(0, max);
  fgPRF->FixParameter(1, dx/fDigiPar->GetPadSizeX(0));
  fgPRF->SetParameter(2, 0.65); fgPRF->SetParLimits(2, 0.45, 10.5);
  fgEdep->Fit(fgPRF, "QBN", "goff", xlo-0.5, xhi+0.5);
  if(!fgPRF->GetNDF()) return NULL;
  chi = fgPRF->GetChisquare()/fgPRF->GetNDF();
  e   = fgPRF->Integral(xlo-0.5, xhi+0.5);

  // apply MC correction
  Float_t gain0 = 210.21387; //(XeCO2 @ 1900V) 
//   Float_t grel[3] = {1., 0.98547803, 0.93164071},
//           goff[3][3] = {
//             {0.05714, -0.09, -0.09},
//             {0., -0.14742, -0.14742},
//             {0., -0.29, -0.393}
//           };
//   Int_t ian=0;
//   if(TMath::Abs(dy)<=0.3) ian=0;
//   else if(TMath::Abs(dy)<=0.6) ian=1;
//   else if(TMath::Abs(dy)<=0.9) ian=2;
//   Int_t isize=0;
//   if(n0<=3) isize=0;
//   else if(n0<=4) isize=1;
//   else isize=2;
  Float_t gain = gain0;//*grel[ian];
  e/=gain;              // apply effective gain
  //e+=goff[ian][isize];  // apply missing energy offset
  
  TVector3 local_pad, local_pad_err;
  Int_t srow, sector= fDigiPar->GetSectorRow(row, srow);
  fDigiPar->GetPadPosition(sector, col, srow, local_pad, local_pad_err);
  //printf("r[%2d] c[%2d] max[%d] lr[%d] n0[%d] cM[%d] nM[%d] dx[%7.4f] dy[%7.4f] loc[%6.2f %6.2f %6.2f] err[%6.2f %6.2f %6.2f] e[%f] chi[%f]\n", row, col, mtyp, lr, n0, cM, nM, dx, dy, local_pad[0], local_pad[1], local_pad[2], local_pad_err[0], local_pad_err[1], local_pad_err[2], e, chi);
  Double_t 
    local[3] = {
      local_pad[0]+dx,
      local_pad[1]+dy,
      local_pad[2]}, 
    global[3], 
    globalErr[3] = {edx, edy, 0.};
  LocalToMaster(local, global);
  
  // process time profile
  for(Int_t idx(1); idx<=n0; idx++){
    Double_t dtFEE = fgDT[0]*(vs[idx]-fgDT[1])*(vs[idx]-fgDT[1])/CbmTrdDigi::Clk(CbmTrdDigi::kFASP);
    if(vxe[idx]>0) vx[idx] += dy/fDigiPar->GetPadSizeY(0);
    fgT->SetPoint(idx-1, vx[idx], vt[idx] - dtFEE);
  }
  xc=vx[n0+2];
  for(Int_t ip(n0); ip<fgT->GetN(); ip++){ fgT->SetPoint(ip, xc, 0); xc+=0.5;}
  
  Double_t time(-21.), 
           edt(26.33),  // should be parametrized as function of da TODO
           tdrift(100); // should depend on Ua
  if(n0>1 && (fgT->Fit("pol1", "QC", "goff") == 0)) { 
    TF1 *f = fgT->GetFunction("pol1");
    time = f->GetParameter(0)-fgDT[2];
    if(TMath::IsNaN(time)) time = -21;
    //dtime += TMath::Abs(f->GetParameter(1)*(vx[n0+1] - vx[1]));
  }
  
  if(CDRAW()){
    Double_t rangex(vx[0]-0.25), rangeX(vx[n0+2]+0.25);
    cvs->cd(1);
    hf->Draw("p");     hf->GetXaxis()->SetRangeUser(rangex, rangeX);
    hf->SetTitle(Form("%d[%d] row[%d] col[%2d] an[%+d] m[%+4.2f] s[%4.2f] E[%7.2f] chi2[%7.2f]", ic, Int_t(vs.size()), row, col, ia, fgPRF->GetParameter(1), fgPRF->GetParameter(2), e, chi));
    hf->GetXaxis()->SetRangeUser(rangex, rangeX);
    hf->GetYaxis()->SetRangeUser(-100., 4500);
    fgEdep->Draw("pl"); 
    fgPRF->Draw("same");
    
    cvs->cd(2);
    hf = (TH1*)hf->DrawClone("p");     
    hf->GetXaxis()->SetRangeUser(rangex, rangeX);
    hf->GetYaxis()->SetRangeUser(-10, 20);
    //hf->SetTitle(Form("%d row[%d] col[%2d] an[%+d] m[%+4.2f] s[%4.2f] E[%7.2f] chi2[%7.2f]", ic, row, col, ia, fgPRF->GetParameter(1), fgPRF->GetParameter(2), fgPRF->Integral(xlo, xhi), fgPRF->GetChisquare()/fgPRF->GetNDF()));
    //       hf->GetXaxis()->SetRangeUser(xlo-0.25, xhi+0.25);
    //       //hf->GetYaxis()->SetRangeUser(0., 4500);
    fgT->Draw("pl");
    cvs->Modified(); cvs->Update();
    cvs->SaveAs(Form("cl_%02d_A%d.gif", ic, ia));
  }
  
  Int_t nofHits = fHits->GetEntriesFast();
  CbmTrdHit *hit = new ((*fHits)[nofHits]) 
          CbmTrdHit(fModAddress, 
                    global, globalErr, 
                    0., // sxy chi, 
                    ic, 
                    e, // energy
                    CbmTrdDigi::Clk(CbmTrdDigi::kFASP)*(t0+time)- tdrift+30.29, 
                    edt);
  hit->SetClassType();
  hit->SetMaxType(tM);
  if(ovf) hit->SetOverFlow();
  if(CWRITE()) cout<<hit->ToString();
  return hit;
}

//_______________________________________________________________________________
Double_t CbmTrdModuleRecT::GetXoffset(Int_t n0) const
{
  Double_t dx(0.), R(0.);
  for(Int_t ir(1); ir<=n0; ir++){ 
    if(vxe[ir]>0) continue; // select rectangular coupling
    R +=vs[ir];
    dx+=vs[ir]*vx[ir];    
  }    
  if(TMath::Abs(R)>0) return dx/R;
  LOG(warn) << GetName() << "::GetDx : Unexpected null sum.";
  return 0. ;
}

//_______________________________________________________________________________
Double_t CbmTrdModuleRecT::GetYoffset(Int_t n0) const
{
  Double_t dy(0.), T(0.);
  for(Int_t it(1); it<=n0; it++) {
    if(vxe[it]>0){ // select tilted coupling
      T +=vs[it];
      dy+=vs[it]*vx[it];
    }
  }
  if(TMath::Abs(T)>0) return dy/T;
  LOG(warn) << GetName() << "::GetDy : Unexpected null sum.";
  return 0.;  
}

//_______________________________________________________________________________
Int_t CbmTrdModuleRecT::LoadDigis(vector<const CbmTrdDigi*> *digis, vector<CbmTrdDigi*> *vdgM, vector<Bool_t> *vmask, ULong64_t &t0, Int_t &cM)
{
/* Load digis information in working vectors. 
 * The digis as provided by the digitizer are replaced by the merged one 
 * according to the vmask map. Digis are represented in the normal coordinate system of
 * (pad width [pw], DAQ time [clk], signal [ADC chs]) with rectangular signals at integer 
 * positions.  
 */
  vs.clear();vse.clear();
  vx.clear();vxe.clear();
  vt.clear();
  
  cM=0;                   // relative position of maximum signal
  Int_t     n0(0),        // number of measured signals
            ovf(1),       // over-flow flag for at least one of the digis
            dt;
  Char_t    ddt,          // signal time offset wrt prompt  
            dt0(0);       // cluster time offset wrt arbitrary t0
  Double_t  r, re(100.),  // rect signal
            t, te(100.),  // tilt signal
            err,          // final error parametrization for signal
            xc(-0.5),     // running signal-pad position
            max(0.);      // max signal
  Int_t j(0), col0(-1), col1(0);
  const CbmTrdDigi *dg(NULL); 
  vector<CbmTrdDigi*>::iterator idgM=vdgM->begin();
  for(vector<const CbmTrdDigi*>::iterator i=digis->begin(); i!= digis->end(); i++, j++){   
    dg = (*i);
    if((*vmask)[j]) {dg = (*idgM); idgM++;}
    if(CWRITE()) cout<<dg->ToString();
    r = dg->GetCharge(t, dt);
    if(t0==0) t0 = dg->GetTimeDAQ(); // set arbitrary t0 to avoid double digis loop
    if(col0<0) GetPadRowCol(dg->GetAddressChannel(), col0); //  initialilze
    ddt=dg->GetTimeDAQ()-t0; if(ddt < dt0) dt0 = ddt;
    
    // check column wise organization
    GetPadRowCol(dg->GetAddressChannel(), col1);
    if(col0 + j != col1) {
      LOG(error) << GetName() << "::LoadDigis : digis in cluster not in increasing order !";
      return 0;
    }
    
    // process tilt signal
    if(t>0) {
      err = te;
      if(t>4094){ovf=-1; err=150.;}
      t-=CbmTrdFASP::GetBaselineCorr();
      n0++;
      if(t>max){ max=t; cM=j;}
    } else err=300.;
    vt.push_back(ddt);
    vs.push_back(t); 
    vse.push_back(err);
    vx.push_back(xc); 
    vxe.push_back(0.035); 
    xc+=0.5;
    
    // process rect signal
    ddt+=dt; if(ddt < dt0) dt0 = ddt;
    if(r>0){ 
      err = re;
      if(r>4094){ovf=-1; err=150.;}
      r-=CbmTrdFASP::GetBaselineCorr();
      n0++; 
      if(r>max){ max=r; cM=j;}
    } else err=300.;
    vt.push_back(ddt);
    vs.push_back(r);
    vse.push_back(err);
    vx.push_back(xc); 
    vxe.push_back(0.); 
    xc+=0.5;
  } 
  
  // remove merged digis if they were created
  if(idgM != vdgM->end()) LOG(warn) << GetName() << "::LoadDigis : not all merged digis have been consumed !";
  for(idgM=vdgM->begin(); idgM!=vdgM->end(); idgM++) delete (*idgM);
  
  // add front and back anchor points if needed
  if(TMath::Abs(vs[0])>1.e-3) {
    xc = vx[0]; ddt=vt[0]; 
    vs.insert(vs.begin(), 0);
    vse.insert(vse.begin(), 300);
    vt.insert(vt.begin(), ddt);
    vx.insert(vx.begin(), xc-0.5);
    vxe.insert(vxe.begin(), 0);
  }
  Int_t n(vs.size()-1);
  if(TMath::Abs(vs[n])>1.e-3) {
    xc = vx[n]+0.5; ddt=vt[n];
    vs.push_back(0);
    vse.push_back(300);
    vt.push_back(ddt);
    vx.push_back(xc);
    vxe.push_back(0.035);
  }
  // recenter time and space profile
  t0+=dt0;
  for(UInt_t idx(0); idx<vx.size(); idx++){ vt[idx]-=dt0; vx[idx]-=cM;}
  return ovf*n0;
}

//_______________________________________________________________________________
Bool_t CbmTrdModuleRecT::MergeDigis(vector<const CbmTrdDigi*> *digis, vector<CbmTrdDigi*> *vdgM, vector<Bool_t> *vmask)
{
/* Merge digis in the cluster if their topology within it allows it although cuts in the 
 * digi merger procedure (CbmTrdFASP::WriteDigi()) were not fulfilled. 
 * Normally this are boundary signals with large time delays wrt neighbors
 */
  
  CbmTrdDigi *dgM(NULL);
  if(digis->size()<2){            // sanity check ... just in case
    LOG(warn) << GetName() << "::MergeDigis : Bad digi config for cluster :";
    return kFALSE;
  }
  
  Double_t r,t;
  Int_t colR, colT,
        dt,
        contor(0);
  Bool_t kFOUND(0);
  for(vector<const CbmTrdDigi*>::iterator idig = digis->begin(), jdig=idig+1; jdig!=digis->end(); idig++, jdig++, contor++){
    const CbmTrdDigi *dgT((*idig)); // tilt signal digi
    const CbmTrdDigi *dgR((*jdig)); // rect signal digi 
    GetPadRowCol(dgR->GetAddressChannel(), colR);
    GetPadRowCol(dgT->GetAddressChannel(), colT);
    
    if(colR!=colT) continue;
    
    dgM = new CbmTrdDigi(*dgT);
    r = dgR->GetCharge(t,dt); dgT->GetCharge(t,dt);
    dt = dgR->GetTimeDAQ() - dgT->GetTimeDAQ();
    dgM->SetCharge(t,r,dt);        
    Int_t rtrg(dgR->GetTriggerType()&2),
          ttrg(dgT->GetTriggerType()&1);      
    dgM->SetTriggerType(rtrg|ttrg); //merge the triggers
    if(CWRITE()){
      cout<<"MERGE"<<endl;
      cout<<dgT->ToString();
      cout<<dgR->ToString();
      cout<<"TO"<<endl;
      cout<<dgM->ToString();
      cout<<"..."<<endl;
    }
    kFOUND=1;
    
    vdgM->push_back(dgM);
    (*vmask)[contor]=1; 
    jdig = digis->erase(jdig);
    if(jdig==digis->end()) break;
  }
  
  if(!kFOUND){
    LOG(warn) << GetName() << "::MergeDigis : Digi to pads matching failed for cluster :";
    return kFALSE;        
  }
  return kTRUE;
}

Float_t CbmTrdModuleRecT::fgCorrXdx=0.005;
Float_t CbmTrdModuleRecT::fgCorrXval[2][NBINSCORRX] = {{
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000, -0.144, -0.091, -0.134, -0.185, -0.120, -0.115, -0.125, -0.125, -0.124, 
 -0.124, -0.122, -0.120, -0.119, -0.116, -0.114, -0.113, -0.111, -0.109, -0.107, 
 -0.105, -0.102, -0.100, -0.098, -0.097, -0.093, -0.091, -0.089, -0.087, -0.084, 
 -0.082, -0.079, -0.077, -0.074, -0.072, -0.068, -0.065, -0.062, -0.059, -0.056, 
 -0.053, -0.049, -0.046, -0.043, -0.039, -0.036, -0.032, -0.029, -0.025, -0.022, 
 -0.018, -0.015, -0.011, -0.007, -0.003,  0.000,  0.003,  0.007,  0.011,  0.014, 
  0.018,  0.022,  0.025,  0.029,  0.032,  0.036,  0.039,  0.043,  0.046,  0.049, 
  0.053,  0.056,  0.059,  0.062,  0.065,  0.068,  0.071,  0.074,  0.077,  0.080, 
  0.082,  0.084,  0.087,  0.090,  0.091,  0.094,  0.096,  0.098,  0.100,  0.103, 
  0.104,  0.107,  0.108,  0.110,  0.113,  0.115,  0.116,  0.120,  0.121,  0.121, 
  0.123,  0.125,  0.124,  0.127,  0.140,  0.119,  0.114,  0.028,  0.049,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000 },
{ 0.003,  0.013,  0.026,  0.039,  0.052,  0.065,  0.078,  0.091,  0.104,  0.118, 
  0.132,  0.145,  0.160,  0.174,  0.189,  0.203,  0.219,  0.234,  0.250,  0.267, 
  0.283,  0.301,  0.319,  0.338,  0.357,  0.375,  0.398,  0.419,  0.440,  0.464, 
  0.491,  0.514,  0.541,  0.569,  0.598,  0.623,  0.660,  0.696,  0.728,  0.763, 
  0.804,  0.847,  0.888,  0.930,  0.988,  1.015,  1.076,  1.128,  1.167,  1.228, 
  1.297,  1.374,  1.443,  1.511,  1.564,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000,  0.000, 
  0.000, -1.992, -1.884, -1.765, -1.703, -1.609, -1.572, -1.493, -1.426, -1.356, 
 -1.309, -1.243, -1.202, -1.109, -1.069, -1.026, -0.970, -0.933, -0.881, -0.844, 
 -0.803, -0.767, -0.721, -0.691, -0.659, -0.629, -0.596, -0.569, -0.541, -0.514, 
 -0.490, -0.466, -0.441, -0.419, -0.397, -0.377, -0.357, -0.337, -0.319, -0.301, 
 -0.283, -0.267, -0.250, -0.234, -0.218, -0.203, -0.189, -0.174, -0.160, -0.145, 
 -0.131, -0.119, -0.104, -0.091, -0.078, -0.064, -0.052, -0.039, -0.026, -0.013, 
 -0.002}
};
Float_t CbmTrdModuleRecT::fgCorrYval[NBINSCORRY][2] = {
  {2.421729, 0.},
  {1.537359, 0.483472},
  {1.1752,   0.},
  {1.154183,-0.090229}
};

ClassImp(CbmTrdModuleRecT)