//
// ================================================================================
// Straw tube response simulation
//
//
//
//  mandatory call at the beginning
//   once to set constants
//
//  void TConst(Double_t Radius, Double_t pSTP, Double_t ArP, Double_t CO2P);
//
// Pysical (garfield) values:
//   Radius  Press (pSTP)  % Ar (ArP)      % CO2 (CO2P)
//     0.4     1           0.9              0.1
//     0.5     1           0.9              0.1
//     0.5     2           0.8              0.2
// 
//
// initialization at each track for MC applications
//
//  define particle and its in-out hits
//  void TInit(Double_t Mass, Double_t Momentum, Double_t InOut[]) 
//
//  void PutWireXYZ(Double_t w1, Double_t w2, Double_t w3, 
//		  Double_t w4, Double_t w5, Double_t w6);     define  the wire
//
//  Double_t PartToTime(Double_t Mass, Double_t Momentum, 
//                                     Double_t InOut[]   );   straw time
//
//  Double_t PartToADC();  ADC charge
//                        signal corresponding to the energy loss of StrawCharge
//                        for energy loss simulation        
//
// where:
//
// pSTP = pressure in bar
// Radius = tube radius
// ArP and CO2P percentages of Argon and CO2
// w1:w6         = coordinates of the extremes of the straw wire (straw axis)
// Mass, Momentm = mass and momentum of the MC particle
// InOut[0:5]    = straw input and output coordinates traversed from the MC particle
//
// Routine for off line reconstruction
//
//  Double_t TimnsToDiscm(Double_t time);  from time (ns) to 
//                                         radius in cm
// where:
// time: measured drift time in ns
//
//
//
// For other applications see the routines listed in TStraw.h
// 
// Author: A. Rotondi (november 2005, revised and extended in  August 2006)
//
// ================================================================================
//

#include "TRandom.h"
#include "TMath.h"
#include "TStraw.h"
#include "TImage.h"
#include "TVector3.h"
#include "TVectorD.h"

#include "TMatrixD.h"
#include "TMatrixDEigen.h"


 ClassImp(TStraw);

// ============================================================
TStraw::TStraw() {

  CNumb=0;
  memset(CumClus,0,sizeof(CumClus));
  memset(PolyaCum,0,sizeof(PolyaCum));
  memset(Xs,0,sizeof(Xs));

  // class constructor

  // clear
  CNumb=0;
  CDist.clear();  
  CDistC.clear();
  CNele.clear();
  CNeleT.clear();
  TeleTime.clear();
  AmplSig.clear();
  Pulse.clear();
  PulseT.clear();
  WDist.clear();
}

// ----------------------------------------------------------------------
void TStraw::TConst(Double_t R1, Double_t P1, Double_t A1, Double_t C1) {

  // set constants for the simulation
  // Input
  // P1    = tube absolute pressure pSTP
  // R1  = tube radius Radius
  // A1 = argon percentage ArPerc
  // C1 = C=2 percentage C=2Perc
  //

  Radius=R1;
  pSTP=P1;
  // Input for the media (volume percentages)
  ArPerc  = A1;
  CO2Perc = C1;

  // cluster dimensions in Ar and CO2 (experimantal values)

//   Double_t PClus[20] = 
//      {0., .656,   .150,   .064,  .035,  .0225, .0155, .0105, .0081, 
//       .0061,  .0049, .0039, .0030, .0025, .0020, .0016, .0012, 
//       .00095, .00075, .00063}; // Fischle fit

   Double_t PClus[20] = 
    {0., .802,   .0707,   .020,  .013,  .008, .006, .005, .006, 
           .008,  .009, .007, .0050, .0040, .0033, .0029, .0025, 
           .0023, .0022, .002}; // Lapique 1st calculation

// Double_t PClus[20] = 
//     {0., .841,   .0340,   .021,  .013,  .008, .006, .004, .003, 
//            .008,  .013, .008, .0050, .004, .0030, .0028, .0025, 
//            .0023, .0022, .002}; // Lapique 2nd calculation


// Double_t PClus[20] = 
//     {0., .656, .150,   .064,  .035,  .0225, .0155, .0105, .0081, 
//            .0061,  .0049, .0039, .0030, .0025, .0020, .0016, .0012, 
//            .00080, .00059, .00045}; // Fischle empirical 

// PDouble_t Clus[20] = 
//     {0., .656,   .148,   .0649,  .0337,  .0244, .0141, .0078, .0095, 
//         .0063,  .0062, .0042, .0028, .0018, .0023, .0017, .0014, 
// 	   .00060, .00050, .00063};    // Fischle exp

   Double_t CO2Clus[20] = 
    {0., .730, .162,   .038,  .020,  .0110, .0147, .0060, .0084, 
           .0052,  .0020, .0042, .0021, .0025, .0038, .0021, .0009, 
           .00013, .00064, .00048}; // Fischle exp

//    Double_t CH4Clus[20] = 
//     {0., .786, .120, .032,  .013,  .0098, .0055, .0057, .0027, 
//            .0029,  .0020, .0016, .0013, .0010, .0012, .0006, .0005, 
//            .00042, .00037, .00033};  // Fischle exp
  

  CH4Perc = 0.07;

  // -----------------------------------------------------
  // gain of the avalanche 
  // Ar/CO2 90/10 1 bar (NTP) MAGY GARFIELD 20-7-2006 
  GasGain=90000.;

  // argon ----------------------------------------------------
  AAr = 39.948;   // Argon (39.948)
  ZAr= 18.0;       // Argon (18)
  RhoAr = pSTP*1.78*1.e-03;  // g/cm3  (1.78 mg/cm3)
  IAr  = 188*1.e-09;     // ionization potential (GeV) (188 eV)
  WiAr =26.7;   // energy to create an ion pair  (standard 26.7 eV)
  NclAr = 23.;  // cluster/cm in Argon
  //  CO2 -----------------------------------------------------
  ACO2 = 44;    // CO2 
  ZCO2 = 22.;   // CO2 
  RhoCO2 = pSTP*1.98*1.e-03;  // g/cm3  CO2 (1.98 mg/cm3)
  ICO2   = 95.8*1.e-09;     // ionization potential (GeV) (188 eV)
  WiCO2  = 33.0;   // energy to create an ion pair (33 eV)
  NclCO2 = 35.5;   // clusters/cm CO2  35.5
  // Methane CH4 ---------------------------------------------------------
  ACH4   = 16;   // CO2 (39.948)
  ZCH4   = 10.;   // CO2 (18)
  RhoCH4 = pSTP*0.71*1.e-03;  // g/cm3  CO2 (0.71 mg/cm3)
  ICH4   = 40.6*1.e-09;     // ionization potential (GeV) (188 eV)
  WiCH4  = 28.0;   // energy to create an ion pair  
  NclCH4 = 25.0;
  // Input for the media (weight percentages) ----------------------------
  ArWPerc  = ArPerc *AAr /(ArPerc*AAr + CO2Perc*ACO2);
  CO2WPerc = CO2Perc*ACO2/(ArPerc*AAr + CO2Perc*ACO2);
 
  // mixture densiies ----------------------------------------------------
  RhoMixCO2 = 1./((ArWPerc/RhoAr) + (CO2WPerc/RhoCO2)); 
  RhoMixCH4 = 1./((ArWPerc/RhoAr) + (CH4WPerc/RhoCH4)); 

  //----------------------------------------------------------------------
  //  particles  (Gev, energy losses in Kev)

  PZeta = 1;   // projectile charge
  piMass = 0.139;  // particle mass (GeV)
  eMass = 0.511/1000.;  // electron mass (GeV) (0.511 MeV)
  prMass = 0.93827; // proton mass (GeV)

  // ---------------------------------------------------------------------
  // thresholds for the straw tubes  (default values) see TInit for current values
  Thresh1=10;
  Thresh2=30;
  // channels for the signal
  Nchann = 500;  

  // ---------------------------------------------Emin------------------------

  NPolya= 100; // steps for the calculation of the Polya distributions
  Xmax=0.;     // Polya istribution is calculated between o and Xmax (see Polya)
  bPolya = 0.5;
  Polya(bPolya);  // cumulative of the Polya distribution
  // -----------------------------------------------------------------------

  // cumulative for the number of electron per cluster

  Double_t Wnorm = (ArPerc*NclAr + CO2Perc*NclCO2);
  CumClus[0]=(ArPerc*NclAr*PClus[0] + CO2Perc*NclCO2*CO2Clus[0])/Wnorm;

  //for(Int_t i=1; i<=20; i++) --->> out of bounds was changed
  for(Int_t i=1; i<20; i++) {
    CumClus[i]=(ArPerc*NclAr*PClus[i] + CO2Perc*NclCO2*CO2Clus[i])/Wnorm
                                       + CumClus[i-1];
  }
  CumClus[20]=1.;

  Double_t sum=0.;
  for(Int_t i=0; i<=19; i++) {
    sum += PClus[i];
    //cout<<" PClus["<<i<<"] = "<<PClus[i]<<endl;
  }  
  //for(Int_t i=0;i<=20;i++) {cout<<"CumClus["<<i<<"] = "<<CumClus[i]<<endl;}
  //cout<<" Sum of Probabilities = "<<sum<<endl;

}

// ==============================================================
void TStraw::PutWireXYZ(Double_t w1, Double_t w2, Double_t w3, 
			Double_t w4, Double_t w5, Double_t w6){
  // get wire coordinates

  Wx1=w1;
  Wy1=w2;
  Wz1=w3;
 
  Wx2=w4;
  Wy2=w5;
  Wz2=w6;

}
// =============================================================


void TStraw::TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {

 // initialization of the constants for each track and each straw

  // transfer of data 
  // track geometrical quantities
  Xin=InOut[0];   Yin=InOut[1];   Zin=InOut[2];
  Xout=InOut[3];  Yout=InOut[4];  Zout=InOut[5];
  PMass = xPMass;
  PMom = xPMom;
 
  // path into the straw
  Dx = sqrt((Xout-Xin)*(Xout-Xin) +
	    (Yout-Yin)*(Yout-Yin) +
	    (Zout-Zin)*(Zout-Zin));

  // clear
  CNumb=0;
  PulseMax=0;
  PulseTime=0;
  CDist.clear();  
  CDistC.clear();
  CNele.clear();
  CNeleT.clear();
  WDist.clear();
  TeleTime.clear();
  AmplSig.clear();
  Pulse.clear();
  PulseT.clear();

 // ---------------------------------------------------------------------

 // initialization of the constants
 // maximum energy transfer Emax (GeV) and related quantities 
 
  PEn = sqrt(PMom*PMom + PMass*PMass); // particle energy GeV
  beta = PMom/PEn;
  gamma= 1/sqrt(1.-beta*beta);
  Double_t begam = beta*gamma;  // local
  Double_t mratio= eMass/PMass;  // local

  // calculation of the polarization Sternheimer factor  for Argon
  Double_t ms   = 2.80;
  Double_t Xst  = log10(begam*sqrt(pSTP));
  Double_t Cs   =-11.92;
  Double_t X1   = 4;
  Double_t Xo   = 1.96;
  Double_t as   = 0.389;
  Delta=0;

  if(Xo <= Xst && Xst <= X1) Delta =  4.6051702*Xst +Cs +as*pow(X1-Xst,ms);
  else if(X1< Xst)           Delta =  4.6051702*Xst +Cs;
  if(Delta < 0) Delta=0;

  // calculation of other typical quantites

  Emax = (1.022*begam*begam/(1.+2.*gamma*mratio+mratio*mratio))/1000.; // GeV

  CsiAr  = pSTP*0.5*0.3071*PZeta*PZeta*ZAr*RhoAr/(beta*beta*AAr)*1e-03; //GeV
  Double_t fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(IAr*IAr);
  EmedAr = 2.*CsiAr*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
  EminAr  = pSTP*2.70*1.e-06;   // GeV/cm  (2.5 keV)
  // most prob energy (GeV) 
  EmpAr = EmedAr + CsiAr*(0.422784 + beta*beta + log(CsiAr/Emax)); 

  CsiCO2  = pSTP*0.5*0.3071*PZeta*PZeta*ZCO2*RhoCO2/(beta*beta*ACO2)*1e-03; //GeV
  fact = 2.*eMass*beta*beta*gamma*gamma*Emax/(ICO2*ICO2);
  EmedCO2 = 2.*CsiCO2*(0.5*TMath::Log(fact)-beta*beta - 0.5*Delta); // GeV/cm
  EminCO2  = pSTP*3.60*1.e-06;   // GeV/cm  (2.5 keV)
  // most prob energy (GeV) 
  EmpCO2 = EmedCO2 + CsiCO2*(0.422784 + beta*beta + log(CsiCO2/Emax)); 

  Csi  = RhoMixCO2 * ((ArWPerc*CsiAr/RhoAr)  + (CO2WPerc*CsiCO2/RhoCO2));;
  Emed = RhoMixCO2 * ((ArWPerc*EmedAr/RhoAr) + (CO2WPerc*EmedCO2/RhoCO2));
  Emin = RhoMixCO2 * ((ArWPerc*EminAr/RhoAr) + (CO2WPerc*EminCO2/RhoCO2));
  Emp  = RhoMixCO2 * ((ArWPerc*EmpAr/RhoAr)  + (CO2WPerc*EmpCO2/RhoCO2));

  // mean weighted interaction
  Wi = (ArPerc*NclAr*WiAr + CO2Perc*NclCO2*WiCO2)
                              /(ArPerc*NclAr + CO2Perc*NclCO2);
 
  Ncl=pSTP*((ArPerc*NclAr) + (CO2Perc*NclCO2))*Emed/Emin; // <--- Cluster/cm
 
  Lcl=1./Ncl;          // mean free path between clusters (cm)
  Ecl = 2.8;           // mean number of electrons per clusters
  Cutoff = Ncl*Ecl*20;  // limit to the number of primary electrons
  Ntote = Ecl*Ncl;  // total mean number of electrons per cm

  // ---------------------------------------------------------------------
  // thresholds for the straw tubes  (current values)
  // total number of electrons * scaling factor (max of signal x electron=1)

  Thresh1=Ncl*Ecl* 0.05;
  Thresh2=Ncl*Ecl* 0.15;

  // -----------------------------------------------------------------------

  TDirCos();  // director cosine of the track
 
  // control prints

//  cout<<" Dx "<<Dx<<" Csi  "<<Csi*1.e+09<<" Emax MeV "<<
//    Emax*1.e+03<<" PMass  "<<PMass<<" gamma "<<
//    gamma<<endl;
//   cout<<" RRise "<<RRise(gamma)<<endl;
//   cout<<" Thresh1 = "<<Thresh1<<" Thresh2 = "<<Thresh2<<endl;

}


// =============================================================

// energy loss in GeV
Double_t TStraw::STEloss() {


  return  gRandom->Landau(Emp*Dx,Csi*Dx);  // in GeV

}

// =============================================================


Double_t TStraw::StrawCharge() {

  // sampling of the ionization energy loss, number of
  // primary electrons and their wire distances for each track
  // energy loss in GeV
  
  // clear 
  Double_t TotEnEle = 0., ClusEner=0., Ecurr;
  CNeleT.clear();
  WDist.clear();

  // threshold for delta rays (>1.5 KeV from NIM A301(1991)202)
  Double_t thr=Csi*Dx*1.e+09/(Ntote*1500.);

  // total energy released by the electrons
  //Number of cluster in the length
  NNClus = Cluster(); // set vectors CDist CDistC e CNele
  // calculation of total energy released by the electrons
  for(Int_t j=1;j<=NNClus;j++){
    ClusEner = 0.;
    for(Int_t jj=1; jj<=(Int_t)(CNele.at(j-1)); jj++){
      Ecurr =  Wi;
      TotEnEle += Ecurr;
      ClusEner += Ecurr;
      // delta rays simulation over 1.keV. Upper threshol 12 keV
      if(gRandom->Uniform() <= thr) {
	TotEnEle += int(Wi/(1.002- gRandom->Uniform()));
        ClusEner += int(Wi/(1.002- gRandom->Uniform()));
      }
    }

    // effective number of electrons per cluster

    CNeleT.push_back(int(ClusEner/Wi));
  }

  // record the distance from the wire of any electron
  // adding the diffusion effects
  Int_t owfl =0;
  for(Int_t k=1; k<=NNClus; k++) {
    for (Int_t kk=0; kk<(Int_t)CNeleT.at(k-1); kk++) {
      if(owfl > (Int_t) Cutoff) break;
      owfl++;
      WDist.push_back(WDistCalc(CDistC.at(k-1)));
    }
  }
 
  return TotEnEle*1.e-09;    // in Gev
}

// =============================================================


Int_t TStraw::Cluster() {

  // calculate the number of clusters CNumb
  // their distance CDist end 
  // their number of electrons CNele

  CNumb=0;
  CDist.clear();   
  CDistC.clear();
  CNele.clear();

  Double_t DisTot = 0;

  for(Int_t k=1; k<=1000; k++) {
    // distance
   
    Double_t path = -Lcl*log(1-gRandom->Uniform()+0.000001);  
    DisTot+= path;
    if(DisTot>Dx) break;
    CNumb=k;
    CDist.push_back(path);   
    CDistC.push_back(DisTot);   
    CNele.push_back((Double_t)Eject() );
  }
 
  return CNumb;
}


// =============================================================


Int_t TStraw::Eject() {
  // find the number of electrons in a cluster
  Int_t Inelect;
  Double_t nelect=0.;
  
  nelect= TMath::BinarySearch(21,CumClus,gRandom->Uniform());
  if(nelect<19) {
    return Inelect = (Int_t) (nelect) +1;
  }
  else {
    return Inelect = (Int_t)(20./(1.01-gRandom->Uniform()));
  }
}


// =============================================================


Double_t  TStraw::RRise(Double_t gamma_r) {

  // interpolate the relativisic rise of the
  // number of cluster per cm, starting from the one
  // measured at the ionization minimum

  Double_t Rise;
  Double_t lg = log10(gamma_r);

  if(1.<=gamma_r && gamma_r <= 2.2){
    Rise = -2.159*lg +1.7;
  }
  else if(2.2<=gamma_r && gamma_r <= 6.){
    Rise = 1.;
  }
  else if(6.<=gamma_r && gamma_r <= 200.){
    Rise = 0.302*lg + 0.765;
  }
  else if(200.<=gamma_r && gamma_r <= 1000.){
    Rise = 0.1431*lg + 1.131;
  }
  else{
    Rise= 1.54;
  }

  return Rise;
}

// ==========================================================================

void TStraw::Polya(Double_t bpar) {

  // calculate the cumulative of the Polya distribution 
  // for samplig the gain fluctuations

  Double_t eps = 0.0001;
  Double_t Dxx = 0.05, xx=0., x1,x2;
  Double_t PMax;
  
  Double_t k=1./bpar -1.;

  // find Xmax

  PMax = eps * pow(k,k)*exp(-k);
  Double_t value = 1.e+06;
  Xmax =2*k;
  while(value > PMax){
    Xmax +=Dxx;
    value = pow(Xmax,k)*exp(-Xmax);
  }
  Xmax += -0.5*Dxx;


  // calculate the cumulative

  Double_t dx= Xmax/NPolya;
  Xs[0]=0;
  for(Int_t i=1; i<NPolya; i++){
    x1 = xx;
    xx += dx;
    x2= xx;
    Xs[i]=x2;
    if(i>1)
      PolyaCum[i] = PolyaCum[i-1] + 
	0.5*dx*(pow(x1,k)*exp(-x1) + pow(x2,k)*exp(-x2))/TMath::Gamma(k+1);
    else 
      PolyaCum[i] = PolyaCum[i-1] + 
	0.5*dx* pow(x2,k)*exp(-x2)/TMath::Gamma(k+1);
  }

  // adjust the normalization

  for(Int_t ii=0; ii<NPolya; ii++) PolyaCum[ii] /= PolyaCum[NPolya-1]; 
 
}


// =====================================================================

Double_t TStraw::PolyaSamp() {

  // sampling a wire gain fluctuation in the gas

  Double_t xr = gRandom->Uniform();
  Int_t n = TMath::BinarySearch(NPolya,PolyaCum,xr);
  Double_t xsamp = Xmax;
  if(n<NPolya-1){
    xsamp= Xs[n] + 
      ((xr-PolyaCum[n])/(PolyaCum[n+1] - PolyaCum[n])) *(Xs[n+1]-Xs[n]);
  }
  return xsamp;
}

// =====================================================================

TVector3  TStraw::WDistCalc(Double_t DisTot) {

  // calculation of the distance of the cluster from the wire
  // DisTot is the cluster distance from the track entrance
  // diffusion effects are considered

  // wire director cosines

  Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) + 
			   (Wy2-Wy1)*(Wy2-Wy1) + 
			   (Wz2-Wz1)*(Wz2-Wz1) );
  Wp = (Wx2-Wx1)/Wlength;
  Wq = (Wy2-Wy1)/Wlength;
  Wr = (Wz2-Wz1)/Wlength;

  // current track coordinates
  Double_t xcor = Calpha * DisTot + Xin;
  Double_t ycor = Cbeta  * DisTot + Yin;
  Double_t zcor = Cgamma * DisTot + Zin;
 
  // cross product VV1= (wire x track) for the distance
  Double_t XX1 = Wr*(ycor-Wy1) - Wq*(zcor-Wz1);
  Double_t YY1 = Wp*(zcor-Wz1) - Wr*(xcor-Wx1);
  Double_t ZZ1 = Wq*(xcor-Wx1) - Wp*(ycor-Wy1);
 
  // vector for the 3-D distance from the  wire (from wire x VV1)
  Double_t XX = Wq*ZZ1 - Wr*YY1;
  Double_t YY = Wr*XX1 - Wp*ZZ1;
  Double_t ZZ = Wp*YY1 - Wq*XX1;
  //cout<<" XYZ  "<<XX<<" "<<YY<<" "<<ZZ<<endl;
  Double_t DDistcm = sqrt(XX*XX + YY*YY +ZZ*ZZ);

  //director cosines of the distance vector
  Double_t cosx = XX/DDistcm;  
  Double_t cosy = YY/DDistcm;  
  Double_t cosz = ZZ/DDistcm;
 
  // sampling of the diffusion
  //Double_t DDist = DDistcm*10.;  // distance in mm

  // longitudinal coefficient of diffusion (GARFIELD) cm --> micron
  // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
  Double_t SigL = DiffLong(DDistcm);
  // ... in cm
  SigL *= 1.e-04;  // in cm 

  // tranverse coefficient (GARFIELD)  cm--> micron 
  // MAGY Ar/CO2 90/10 20-7-2006 GARFIELD
  Double_t SigT = DiffTran(DDistcm);
  SigT *= 1.e-04;  // in cm 

  // sampling of Longitudinal and Transverse diffusion
  Double_t difL = gRandom->Gaus(0.,SigL);
  Double_t difT = gRandom->Gaus(0.,SigT);
 
  // vector addition to the distance
  // the transverse component has the same dir cos of the wire
  XX += difL*cosx + difT*Wp;
  YY += difL*cosy + difT*Wq;
  ZZ += difL*cosz + difT*Wr;
  //cout<<" XYZ+ dif  "<<XX<<" "<<YY<<" "<<ZZ<<endl;
  //cout<<"  --------------------------------------------------- "<<endl;
  TVector3 TDist(XX,YY,ZZ);

  return TDist;

}


// =====================================================================

void TStraw::TDirCos() {

  // director cosines of the track (called for each track)

  //path into the straw
  Rpath = sqrt((Xout-Xin)*(Xout-Xin) +
	       (Yout-Yin)*(Yout-Yin) +
	       (Zout-Zin)*(Zout-Zin));

  //director cosines
  Calpha = (Xout-Xin)/Rpath;
  Cbeta  = (Yout-Yin)/Rpath;
  Cgamma = (Zout-Zin)/Rpath;
}


// =====================================================================

Int_t TStraw::TimeEle() {

  TeleTime.clear();

  // calculate the arrival times (ns) for each electron in TeleTime
  // from cm to ns  <---------------------------------
  // return the number of electrons
  // Author:A. Rotondi 20-7-2006

  Int_t ido = WDist.size();
  Double_t etime;

  // cutoff the total charge to 20 times the average
  if(ido > (Int_t) Cutoff) ido =  (Int_t) Cutoff;

  for(Int_t k=0; k<ido; k++) {
    Double_t dst = WDist.at(k).Mag();  // distance in cm

    // space to time calculation from GARFIELD

    if(pSTP <1.9){
      if(Radius<0.5){
      // V=1600V diameter 4 mm 90/10
	etime= -1.624e-05 + 0.1258*dst + 0.8079*pow(dst,2) -2.918*pow(dst,3)
	  + 10.33*pow(dst,4) -10.84*pow(dst,5);
      }
      else{
      // V=1600V diameter 5 mm   90/10
	etime= -6.763e-05 + 0.1471*dst +0.3625*pow(dst,2) +0.3876*pow(dst,3)
	  +1.04*pow(dst,4) -1.693*pow(dst,5);
      }
    }
    else{
      // 2 bar 2000V, 5 mm, 80/20
      etime= -0.0001014 + 0.1463*dst -0.1694*pow(dst,2) +2.4248*pow(dst,3)
	-1.793*pow(dst,4);
    }

    etime*=1000.;   // nano seconds 
    TeleTime.push_back(etime);    
  }

  return TeleTime.size();
}

// =====================================================================

Double_t TStraw::Signal(Double_t t, Double_t t0) {

  // electric signal at time  t of a cluster arriving at t0
  Double_t elesig = 0;
  Double_t x= t - t0;

  if (x>0){
      Double_t A = 1.03e-03;
      Double_t B = 3.95;
      Double_t C = 0.228;
      Double_t D = -3.839e-02;
      Double_t E = 1.148e-03;    

      elesig = A*exp(B*log(x)-C*x)*(1+D*x+E*x*x);
  }

  return elesig;
}

// =====================================================================

Int_t  TStraw::StrawSignal(Int_t nsteps) {

  // creation of nstep values of
  // the straw global Pulse (sum on all clusters)
  // return the number of primary electrons

  PulseMax=0;
  Pulse.clear();
  PulseT.clear();
  AmplSig.clear();

  Int_t neltot = TimeEle(); // creation and size of TeleTime (electron times)

  Double_t Tmax = 1.e-25;
  for(Int_t k=0; k< neltot; k++) if(Tmax<TeleTime.at(k)) Tmax=TeleTime.at(k);
  Tmax +=  100.;
  
  Double_t Dt = Tmax/nsteps;  // number of steps of the signal

  //AmplSig is the amplitude of each electron

  for(Int_t j=0; j< neltot; j++){
    AmplSig.push_back(bPolya*PolyaSamp());
  }

  // creation of the signal Pulse(PulseT) time in ns
  for(Int_t j=0; j< nsteps; j++){
    Double_t te = j*Dt;
    Double_t sumele=0.;
    for(Int_t jj=0; jj< neltot; jj++){
      Double_t te0 = TeleTime.at(jj);
      sumele += AmplSig.at(jj)*Signal(te,te0);
    }
    
    Pulse.push_back(sumele);
    PulseT.push_back(te);
  }

  // add a random noise (3% of maximum) plus 
  // a 1% of 200 ns periodic signal

   Double_t Pmax = 1.e-25;
  for(Int_t k=0; k<(Int_t)Pulse.size(); k++) 
                  if(Pmax<Pulse.at(k)) Pmax=Pulse.at(k);

  for(Int_t k=0; k<(Int_t)Pulse.size(); k++){
    Pulse.at(k) += 0.03*Pmax*gRandom->Uniform()
                       + 0.01*Pmax*(sin(6.28*PulseT.at(k)/120.));
      //    if(Pulse[k]<0) Pulse[k] *= -1.;
  }

  PulseMax=Pmax;

  // set variable threshold for the signals
   Thresh1=0.05*Pmax;
   Thresh2=0.15*Pmax;

  return neltot;
}

// =====================================================================

Int_t  TStraw::StrawTime() {

  // simulate the discrimination of the straw signal 
  // and give the time
  // discriminator technique: set 2 threshold, select the first one
  // (the low one) only if the second one is fired (FINUDA system)


  PulseTime=0.;

  Int_t ind=0;
  Int_t flag1=0;
  Int_t flag2=0;

  for(Int_t k=0; k<(Int_t) Pulse.size(); k++){
    if(flag1==0 && Pulse[k]>Thresh1) {
      flag1=1;
      ind=k;
    }

    if(Pulse[k]<Thresh1) {flag1=0;  ind=0;} // reset if signal decreases

    if(flag1==1 && Pulse[k]>Thresh2) {
      flag2=1;
      break;
   }

  }
  if(flag1==1 && flag2==1) PulseTime=PulseT.at(ind);
 
  return ind;
} 


// =====================================================================

Double_t  TStraw::TrueDist(Double_t Point[]) {
  // service routine that finds the distance in cm from the wire
  // by knowing the wire coordinates (class variables)
  // and the input-output points Point[6]


  Double_t truedist = 0;

  // wire director cosines
  Double_t Wlength = sqrt( (Wx2-Wx1)*(Wx2-Wx1) + 
			   (Wy2-Wy1)*(Wy2-Wy1) + 
			   (Wz2-Wz1)*(Wz2-Wz1) );
  Wp = (Wx2-Wx1)/Wlength;
  Wq = (Wy2-Wy1)/Wlength;
  Wr = (Wz2-Wz1)/Wlength;

  // director cosines of the given track
  Double_t Modu = sqrt( (Point[3]*Point[0])*(Point[3]*Point[0]) +
		      (Point[4]*Point[1])*(Point[4]*Point[1]) +
		      (Point[5]*Point[2])*(Point[5]*Point[2]) );
  Double_t dcx = (Point[3]-Point[0])/Modu;
  Double_t dcy = (Point[4]-Point[1])/Modu;
  Double_t dcz = (Point[5]-Point[2])/Modu;
 
  //distance formula
  Double_t p1 = (Point[0]-Wx1)*(dcy*Wr-dcz*Wq);
  Double_t p2 =-(Point[1]-Wy1)*(dcx*Wr-dcz*Wp);  
  Double_t p3 = (Point[2]-Wz1)*(dcx*Wq-dcy*Wp);
  Double_t Det = p1+p2+p3;
  Double_t Disc = sqrt( (dcy*Wr-dcz*Wq)*(dcy*Wr-dcz*Wq) +
			(dcz*Wp-dcx*Wr)*(dcz*Wp-dcx*Wr) +
			(dcx*Wq-dcy*Wp)*(dcx*Wq-dcy*Wp) );
  if(Disc >0) truedist = TMath::Abs(Det/Disc);
 
  return truedist;   // distance in cm
}

// =====================================================================

Double_t  TStraw::TimnsToDiscm(Double_t time) {

  // distance in cm from time in ns for pSTP =1 and pSTP=2 atm 
  // utility routine for the track reconstruction
  //last update: A. Rotondi 3-3-2007

  Double_t drift, truer;


  // 1  absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
  //  ns -->  mm
 
  if(pSTP < 1.9) {
    if(Radius < 0.5){
    // 1600 V 4 mm variable theshold
//     drift =    -5.68738e-02 +7.06505e-03*time 
//              +3.70423e-03*pow(time,2) -1.39795e-04*pow(time,3)
//              +2.59863e-06*pow(time,4) -2.75729e-08*pow(time,5)
//              +1.69315e-10*pow(time,6) -5.59127e-13*pow(time,7)
//              +7.66371e-16*pow(time,8);


      // 1600 V 4 mm fixed threshold

//       drift =      7.50700e-02    -4.58303e-02*time 
//               +6.01476e-03*pow(time,2)  -1.78491e-04*pow(time,3)
//               +2.82377e-06*pow(time,4)  -2.60709e-08*pow(time,5)
//               +1.40082e-10*pow(time,6)  -4.04386e-13*pow(time,7)
//               +4.81504e-16*pow(time,8);

      drift = -0.070 -1.18253e-02  
              -1.65315e-02 *time
              +6.18870e-03 *pow(time,2)   
              -2.50647e-04 *pow(time,3)   
              +5.15779e-06 *pow(time,4) 
	      -6.04826e-08 *pow(time,5)  
              +4.07345e-10 *pow(time,6)
              -1.46437e-12 *pow(time,7)
	      +2.17349e-15 *pow(time,8);   
    } 


    else{
    // 1600 V 5 mm 

      drift =  - 7.43091e-02   
               + 6.34826e-03  *time  
               + 2.70409e-03  *pow(time,2) 
               - 7.96339e-05  *pow(time,3) 
               + 1.12156e-06  *pow(time,4)  
               - 8.76461e-09  *pow(time,5) 
               + 3.87359e-11  *pow(time,6)  
               - 9.05398e-14  *pow(time,7)  
               + 8.69003e-17  *pow(time,8);  
    }
  }

  else {
     // 2 absolute atm  5 mm   80/20 (1 bar overpressure)
    if(time < 100.) {
      drift = + 1.71849e-02   
	- 3.99942e-02 *time 
	+ 8.96815e-03 *pow(time,2) 
	- 3.23438e-04 *pow(time,3) 
	+ 4.88349e-06 *pow(time,4) 
	- 2.27659e-08 *pow(time,5)  
	- 1.00527e-10 *pow(time,6)  
	- 6.18934e-13 *pow(time,7)  
	+ 2.69373e-14 *pow(time,8) 
	- 1.26629e-16 *pow(time,9);
    }  

    else{

      drift =  - 2.65873e-01   
	+ 6.20514e-02 *time 
	+ 3.42564e-05 *pow(time,2)  
	- 1.17871e-05 *pow(time,3) 
	+ 1.59485e-07 *pow(time,4) 
	- 9.61794e-10 *pow(time,5)
	+ 2.79989e-12 *pow(time,6)  
	- 3.20470e-15 *pow(time,7); 
    }


 // correction of the systematic shift for r=0.5 cm  2 atm
//     if( drift<=4.8){      
//       truer=    2.36461e-01  -5.01598e-01*drift + 3.89468*pow(drift,2)
// 	-4.17036*pow(drift,3)  +2.19117*pow(drift,4)
// 	-6.02374e-01*pow(drift,5)  +8.32138e-02*pow(drift,6)   
// 	-4.55919e-03*pow(drift,7);
      
//       drift= TMath::Abs(truer);
//    }
  }
 
  drift =  0.1*drift;  
  if(drift < 0.) drift=0.;
  return drift;  

}

// --------------------------------------------------------------------------------

Double_t TStraw::PartToTime(Double_t xPMass, Double_t xPMom, Double_t InOut[]) {

 // find the time of a particle of mass xPmass, momentum xPMom, with 
 // input-output coordinate InOut[6]
 // Useful for MC pplication after a call to PutWireXYZ

  TInit(xPMass, xPMom, InOut);    // start the event                   
  Out1 = StrawCharge();           // energy loss (GeV) to generate charge
  Out2 = StrawSignal(Nchann);     // generate the straw signal
  Out3 = StrawTime();             // find the straw drift time PulseTime
  
  return PulseTime;
  
}

// --------------------------------------------------------------------------------   
Double_t  TStraw::PartToADC() {

  // return the energy loss in the tube as charge signal
  // taking into account the Polya fluctuations
  Double_t ADCsignal=0.;

  for(Int_t j=1; j<= NNClus; j++){
   
    for(Int_t jc=1; jc<=(Int_t) CNeleT[j-1]; jc++) 
      ADCsignal += bPolya * GasGain * PolyaSamp();
  }

  return ADCsignal;
}
// --------------------------------------------------------------------------------   
Double_t  TStraw::DiffLong(Double_t Distcm) {

  // return the longitudinal diffusion
  // from cm to microns 
  // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP) 20-7-2006
  //
  Double_t DiffMic;

  if(pSTP<1.9){
    if(Radius < 0.5){
      // 4 mm 1600 V   90/10
      DiffMic = 0.896 + 1387.*Distcm - 1.888e+04*pow(Distcm,2)
	+ 1.799e+05*pow(Distcm,3) - 9.848e+05*pow(Distcm,4)
	+ 3.009e+06*pow(Distcm,5) - 4.777e+06*pow(Distcm,6)
	+ 3.074e+06*pow(Distcm,7);
    }
    else{
      // 5 mm 1600 V
      DiffMic = 1.537 + 1246.*Distcm - 1.357e+04*pow(Distcm,2)
	+ 1.049e+05*pow(Distcm,3) - 4.755e+05*pow(Distcm,4)
	+ 1.211e+06*pow(Distcm,5) - 1.6e+06*pow(Distcm,6)
	+ 8.533e+05*pow(Distcm,7);
    }
  }
  
  else{
    // 2000V 5 mm 2 bar 80/20
    DiffMic = 2.135 +818.*Distcm - 1.044e+04*pow(Distcm,2)
      + 8.31e+04*pow(Distcm,3) - 3.492e+05*pow(Distcm,4)
      + 7.959e+05*pow(Distcm,5) - 9.378e+05*pow(Distcm,6)
      + 4.492e+05*pow(Distcm,7);
  }
  return DiffMic;
}
// --------------------------------------------------------------------------------   
Double_t  TStraw::DiffTran(Double_t Distcm) {

  // return the transverse diffusion in microns
  // from cm to microns 
  // MAGY GARFIELD Ar/CO2 90/10 1 bar (NTP)  20-7-2006
  //
  Double_t DiffMic;

  if(pSTP<1.9){
    if(Radius < 0.5){
      // 4 mm 1600 V   90/10
//       DiffMic = + 0.8513 + 1648.*Distcm   - 1.085e+04*pow(Distcm,2)
// 	+ 7.38e+04*pow(Distcm,3) - 3.025e+05*pow(Distcm,4)
// 	+ 6.067e+05*pow(Distcm,5) - 4.643e+04*pow(Distcm,6);
      DiffMic = + 1.482 + 1529.*Distcm   - 6755.*pow(Distcm,2)
	+ 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
	+ 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
    }
    else{
      // 5 mm 1600 V  90/10
      DiffMic = + 1.482 + 1529.*Distcm   - 6755.*pow(Distcm,2)
	+ 2.924e+04*pow(Distcm,3) - 0.9246e+05*pow(Distcm,4)
	+ 1.548e+05*pow(Distcm,5) - 1.002e+05*pow(Distcm,6);
    }
  }
  else{

    // 5 mm   2000 V  2 bar  80/20
    DiffMic = +2.094 + 1138.*Distcm   - 7557.*pow(Distcm,2)
      + 2.968e+04*pow(Distcm,3) - 6.577e+04*pow(Distcm,4)
      + 7.581e+04*pow(Distcm,5) - 3.497e+04*pow(Distcm,6);
  }
  return DiffMic;
}
// --------------------------------------------------------
Double_t  TStraw::DistEle(Double_t time) {

  // dist  in cm  from time in ns for pSTP =1 and pSTP=2 
  // utility routine for the >SINGLE< electron reconstruction
  //last update: A. Rotondi 20-7-2006

  Double_t drift;

  // 1  absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
  //  ns -->  cm

  time *= 0.001; // time in micro s

  if(pSTP < 2.) {
    if(Radius <0.5){
      // 1600 V 4 mm  90/10
      drift =   0.001629 + 6.194*time 
	- 56.55*pow(time,2) + 355.8*pow(time,3)
	- 903.2*pow(time,4);
    }
    else{
      // 1600 V 5 mm  90/10
      drift =   0.003365 + 5.734*time 
	- 41.88*pow(time,2) + 191.2*pow(time,3)
	- 333.4*pow(time,4) ;
    }
  }

  else {
       
    // 2 absolute atm 2000 V 5 mm 80/20 (1 bar overpressure)

    drift =   0.003365 + 7.056*time 
            - 62.28*pow(time,2) + 306.1*pow(time,3)
            - 558.7*pow(time,4);

  }
  return drift;  // in cm

}