//
// ============================================================================
//
// Straw tube response simulation
// Corrections to StrawCharge and Eject for delta electrons (dec 2007)
//
//  version with one threshold only Julich signal 
//
//  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
// 
// ------------------------------------------------------------------------------
//
// Full simulation: calls at each track for MC applications
//
//  define particle and its in-out hits
//
//  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 TimnsToDiscm(Double_t time);  from time (ns)
//                                         to the reconstructed distance (cm)
//
//  Double_t PartToADC();  ADC charge
//                        signal corresponding to the energy loss of StrawCharge
//                        for energy loss simulation        
//
//
//-----------------------------------------------------------------------------
//
// Fast simulation: calls at each track for MC applications
//
//  define particle and its in-out hits
//
//  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 FastRec(Double_t TrueDcm, Int_t flag);  fast approximation: 
//                                                   from the true distance   (cm)
//                                                   to the reconstructed one (cm)
//                                                   flag=0 simulation
//                                                   flag=1 Julich exp curve COSY-TOF (old)
//                                                   flag=2 Julich exp curve COSY-TOF (Feb 2011)
//                                                   flag=3 Julich exp curve from prototype
//                                                   (active with press = 2 in TConst)
//
//  void TInit(Double_t xPMass, Double_t xPMom, Double_t InOut[]); constant initialization
//
//  Double_t FastPartToADC();  ADC charge
//                        signal corresponding to the energy loss of StrawCharge
//                        with the Urban model for energy loss simulation        
//
// ----------------------------------------------------------------------------
// where:
//
// time: measured drift time in ns
//
// 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
//
//
//
// For other applications see the routines listed in PndSttSingleStraw.h
// 
// Author: A. Rotondi (november 2005, revised and extended in  August 2006)
//
// ================================================================================
//

#include <iostream>
#include <string>
#include "TRandom.h"
#include "TMath.h"
#include "PndSttSingleStraw.h"
#include "TImage.h"
#include "TVector3.h"
#include "TVectorD.h"

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


using namespace std;

 ClassImp(PndSttSingleStraw);


// ============================================================
PndSttSingleStraw::PndSttSingleStraw() {

  CNumb=0;
  memset(CumClus,0,sizeof(CumClus));
  memset(CH4Clus,0,sizeof(CH4Clus));
  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 PndSttSingleStraw::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 = CO2 percentage   ---> C02Perc
  //

  Radius = R1;
  pSTP   = P1;
  Field  = 1;   // magnetic field as a default (to change use PutField)
  // 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=100000.;
  // Ar/CO2 90/10 2 bar (NTP) 1800 V experimental Juelich value
  GasGain = 50000.;

  // argon ----------------------------------------------------
  AAr = 39.948;     // Argon (39.948)
  ZAr = 18.0;       // Argon (18)
  RhoAr = pSTP*1.662*1.e-03;  // g/cm3  (1.662 mg/cm3)
  IAr  =  188*1.e-09;         // ionization potential (GeV) (188 eV)
  WiAr =  27.0;     // energy to create an ion pair (standard 26.7 eV)
  NclAr = 25.;      // cluster/cm in Argon
  //  CO2 -----------------------------------------------------
  ACO2 = 44;        // CO2 
  ZCO2 = 22.;       // CO2 
  RhoCO2 = pSTP*1.842*1.e-03;  // g/cm3  CO2 (1.842 mg/cm3)
  ICO2   = 95.8*1.e-09;        // ionization potential (GeV) (96 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) (45 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 densities ----------------------------------------------------
  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.4;
  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++) {
    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 PndSttSingleStraw::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 PndSttSingleStraw::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  = 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  = 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);
  // number of clusters
  Double_t nAr  = ArWPerc  * RhoMixCO2/AAr;
  Double_t nCO2 = CO2WPerc * RhoMixCO2/ACO2;
  //Ncl=pSTP*((ArPerc*NclAr)+(CO2Perc*NclCO2))*Emed/Emin; // <--- Cluster/cm
  Ncl=pSTP*((nAr*NclAr+nCO2*NclCO2)/(nAr+nCO2))*Emed/Emin; // <--- Cluster/cm
 
  Lcl=1./Ncl;           // mean free path between clusters (cm)
  Ecl = 2.8;            // mean number of electrons per clusters (def 2.8)
  Cutoff = Ncl*Ecl*40;  // limit to the number of primary electrons (*20)
  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 = 0.5*Ncl*Ecl* 0.05;  // default 0.5*Ncl*Ecl* 0.05;
  Thresh2 = 0.5*Ncl*Ecl* 0.10;  // default 0.5*Ncl*Ecl* 0.10;
  OffT = 10.; // 10.;   // time offset in ns for trheshold calibration;

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

  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 PndSttSingleStraw::STEloss() {
  return  gRandom->Landau(Emp*Dx,Csi*Dx);  // in GeV
}

// =============================================================
Double_t PndSttSingleStraw::STUrban() {
 
//
// energy loss in GeV from the Urbam distribution NIM A 362(1995)416
// Author: A. Rotondi December 2007
//

  Double_t Sig1Ar, Sig2Ar, Sig3Ar; 
  Double_t Sig1Mix, Sig2Mix, Sig3Mix;
  Double_t Sig1CO2, Sig2CO2, Sig3CO2, IMix, nAr, nCO2;
  Double_t f1, f2, e1, e2, e2f, ru, Cu, Cuc;
  Double_t N1Mix, N2Mix, N3Mix, E1Mix, E2Mix, E3Mix, ZMix;
  Double_t Sig21, Sig22, Sig23, Medd;
  
  Double_t eVGeV = 1.e-09;

  // ru= 0.6; 
  ru= 0.55; 

  Cuc = 2.*PMass*beta*beta*gamma*gamma;

  // calculation for Argon
  Cu = EmedAr;
  f2 = 2./ZAr;
  f1 = 1. -f2;
  e2 = 10.*ZAr*ZAr*eVGeV;
  e2f = TMath::Power(e2,f2);
  e1 = TMath::Power(IAr/e2f,1./f1);
  Sig1Ar = (1.-ru)*Cu*(f1/e1)*(TMath::Log(Cuc/e1)  - beta*beta)/
                              (TMath::Log(Cuc/IAr) - beta*beta);
  Sig2Ar = (1.-ru)*Cu*(f2/e2)*(TMath::Log(Cuc/e2)  - beta*beta)/
                              (TMath::Log(Cuc/IAr) - beta*beta);
  Sig3Ar = Cu*Emax*ru / (IAr*(Emax+IAr)*TMath::Log((Emax+IAr)/IAr));

  // Calculation for CO2
  Cu = EmedCO2;
  f2 = 2./ZCO2;
  f1 = 1. -f2;
  e2 = 10.*ZCO2*ZCO2*eVGeV;
  e2f = TMath::Power(e2,f2);
  e1 = TMath::Power(ICO2/e2f,1./f1);
 
  Sig1CO2 = (1.-ru)*Cu*(f1/e1)*(TMath::Log(Cuc/e1)   - beta*beta)/
                               (TMath::Log(Cuc/ICO2) - beta*beta);
  Sig2CO2 = (1.-ru)*Cu*(f2/e2)*(TMath::Log(Cuc/e2)   - beta*beta)/
                               (TMath::Log(Cuc/ICO2) - beta*beta);
  Sig3CO2 = Cu*Emax*ru / (ICO2*(Emax+ICO2)*TMath::Log((Emax+ICO2)/ICO2));

  // calculation for the mixture
  nAr  = ArWPerc  * RhoMixCO2/AAr;
  nCO2 = CO2WPerc * RhoMixCO2/ACO2;
  IMix = TMath::Exp((nAr*ZAr*TMath::Log(IAr) + nCO2*ZCO2*TMath::Log(ICO2))/
		    (nAr*ZAr + nCO2*ZCO2));

  ZMix = ArWPerc*ZAr + CO2WPerc*ZCO2;
  f2 = 2./ZMix;
  f1 = 1. -f2;
  E2Mix = 10.*ZMix*ZMix*eVGeV;
  e2f = TMath::Power(E2Mix,f2);
  E1Mix = TMath::Power(IMix/e2f,1./f1);

  Sig1Mix = Dx*RhoMixCO2*(ArWPerc*Sig1Ar/RhoAr + CO2WPerc*Sig1CO2/RhoCO2);
  Sig2Mix = Dx*RhoMixCO2*(ArWPerc*Sig2Ar/RhoAr + CO2WPerc*Sig2CO2/RhoCO2);
  Sig3Mix = Dx*RhoMixCO2*(ArWPerc*Sig3Ar/RhoAr + CO2WPerc*Sig3CO2/RhoCO2); 

  NUrban = Sig1Mix+Sig2Mix+Sig3Mix;  // total number of collisions (mean value)

  // Urban distribution calculation
  N1Mix = gRandom->Poisson(Sig1Mix);
  N2Mix = gRandom->Poisson(Sig2Mix);
  N3Mix = gRandom->Poisson(Sig3Mix);

  Int_t N3 = TMath::Nint(N3Mix+gRandom->Uniform());

  E3Mix = 0.;
  for(Int_t j=1;j<=N3;j++){
    E3Mix += ((Emax+IMix)*IMix)/(IMix + Emax*(1.-gRandom->Uniform()));
  }

  // variance calculation
  
  Eup = IMix/(1.-0.98*Emax/(Emax+IMix));  // 98% of the delta electron area (GeV)
  Sig21 = IMix*(Emax+IMix)/Emax;
  Medd  = Sig21*TMath::Log(Eup/IMix);
  Sig22 = Sig21*Eup - Medd*Medd;
  Sig23 =  Sig1Mix*E1Mix*E1Mix + Sig2Mix*E2Mix*E2Mix 
         + Sig3Mix*Medd*Medd  + Sig3Mix*Sig22;

  AvUrb = Sig1Mix*E1Mix + Sig2Mix*E2Mix + Sig3Mix*Medd;
  
  SigUrb = TMath::Sqrt(Sig23);

  // return the energy lost in keV
  return  N1Mix*E1Mix + N2Mix*E2Mix + E3Mix;  // in GeV
}


// =============================================================
Double_t PndSttSingleStraw::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)
  // 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;
    }

    // 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 PndSttSingleStraw::Cluster() {

  // calculate the number of clusters CNumb
  // their distance CDist and 
  // 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 PndSttSingleStraw::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 {
    // long delta electron tail
    return Inelect = (Int_t)(20./(1.03-gRandom->Uniform()));
  }
}


// =============================================================
Double_t  PndSttSingleStraw::RRise(Double_t gamma2) {

  // 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(gamma2);

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

  return Rise;
}


// ==========================================================================
void PndSttSingleStraw::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 PndSttSingleStraw::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  PndSttSingleStraw::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 PndSttSingleStraw::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 PndSttSingleStraw::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);
	etime*=1000.;   // from micros to ns
      }
      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);
	etime*=1000.;   // from micros to ns
      }
    }
    else{
      // 2.2 bar 2150 V, 5 mm, 90/10 and/or 80/20 wire 20 micron 
      // (may 2009 data)
      if(Field==1 &&  CO2Perc <0.15){
	etime= -1.569123e-03 + 2.228707e-01*dst -2.100147e-01*pow(dst,2)  
	  +1.240429*pow(dst,3)  -5.929662e-01*pow(dst,4);
	etime=etime*1000.;  // from micros to ns
      }
      else if(Field==1 &&  CO2Perc >0.15){
	etime = -3.022419e-04 + 1.490671e-01*dst -1.947376e-01*pow(dst,2)
	  +2.764081e+00*pow(dst,3) -2.129739e+00*pow(dst,4);
	etime=etime*1000.;  // from micros to ns
      }
      //no magnetic field 90/10
      else{
	etime= -1.472305e-03 + 2.139223e-01*dst -5.835399e-02*pow(dst,2) 
	  -6.727544e-02*pow(dst,3) +6.334636e-01*pow(dst,4);
	etime=etime*1000.;  // from micros to ns
      }
    }

    TeleTime.push_back(etime);    
  }

  return TeleTime.size();
}


// =====================================================================
Double_t PndSttSingleStraw::Signal(Double_t t, Double_t t0) {

  // electric signal at time  t of a cluster arriving at t0
  Double_t elesig;
  //Double_t A = 1.03e-03; //[R.K. 01/2017] unused variable?
  //Double_t B = 3.95; //[R.K. 01/2017] unused variable?
  //Double_t C = 0.228; //[R.K. 01/2017] unused variable?
  //Double_t D = -3.839e-02; //[R.K. 01/2017] unused variable?
  //Double_t E = 1.148e-03; //[R.K. 01/2017] unused variable?

  Double_t x= t - t0; 
  //if(x>0) elesig = A*exp(B*log(x)-C*x)*(1+D*x+E*x*x);  // Sokolov Signal
  if(x>0) elesig =pow(2.*x/10.,2)*exp(-2.*x/10.);  // Wirtz signal
  else elesig = 0.;

  return elesig;
}


// =====================================================================
Int_t  PndSttSingleStraw::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 (constant fraction)
  // Thresh1=0.05*Pmax;
  // Thresh2=0.15*Pmax;

  return neltot;
}


// =====================================================================
Int_t  PndSttSingleStraw::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=1;
  
  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)-OffT;  // time in ns;
 
  return ind;
} 


// =====================================================================
Int_t  PndSttSingleStraw::StrawTot() {

  // simulate the discrimination of the straw signal 
  // and give the time window over threshold of the signal
  // return the integer difference between the two memory locations
  // and set the two times PulseTime and PulseTime1


  PulseTime=0.;
  PulseTime1=0.;

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

    //  find the second threshold
    if(Pulse[k]>Thresh1)ind1=k;
  }
  if(flag1==1 && flag2==1) PulseTime=PulseT.at(ind);
  PulseTime1 = PulseT.at(ind1);

  return ind1-ind;
} 


// =====================================================================
Double_t  PndSttSingleStraw::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  PndSttSingleStraw::TimnsToDiscm(Double_t time) {

  // distance in cm from time in ns for pSTP =1 and pSTP=2 atm 
  // where the time is given by the threshold on the signal
  // utility routine for the track reconstruction
  // last update: A. Rotondi 24-5-2009

  Double_t drift;


  // 1  absolute atm (NTP) 20-7-2006 GARFIELd Ar/CO2 90/10 MAGY
  //  ns -->  mm
 
  if(pSTP < 1.9) {
    if(Radius < 0.5){
 
      drift =  -0.0140 -1.37281e-01  
               +5.13978e-02*time
               +7.65443e-05*pow(time,2)   
               -9.53479e-06*pow(time,3)   
               +1.19432e-07*pow(time,4) 
	       -6.19861e-10*pow(time,5)  
               +1.35458e-12*pow(time,6)
               -1.10933e-15*pow(time,7);   
    } 


    else{
    // 1600 V 5 mm 
      if(time < 120.){

      drift =    0.0300 -1.07377e-01  
                 +3.65134e-02*time  
                 +1.20909e-03*pow(time,2) 
                 -4.56678e-05*pow(time,3) 
                 +6.70207e-07*pow(time,4)  
                 -4.99204e-09*pow(time,5) 
                 +2.19079e-11*pow(time,6)  
                 -8.01791e-14*pow(time,7)  
                 +2.16778e-16*pow(time,8);
      }
      else{    

      drift =    0.0300 -8.91701e-01   
                 +4.68487e-02*time  
                 +1.00902e-03*pow(time,2) 
                 -4.00359e-05*pow(time,3) 
                 +6.23768e-07*pow(time,4)  
                 -5.20556e-09*pow(time,5) 
                 +2.41502e-11*pow(time,6)  
                 -5.85450e-14*pow(time,7)  
                 +5.77250e-17*pow(time,8);  
      }
    }
  }

  else {
    // 2 absolute atm  5 mm  2150 V   90/10 (1 bar overpressure may 2009)
    if(CO2Perc < 0.15){
      //  constant fraction threshold 5%  5 ns offset
      //        drift = -1.343144e-01  
      // 	 +5.101958e-02 *time    
      // 	 +9.116481e-05 *pow(time,2)
      // 	 -5.787444e-06 *pow(time,3)
      // 	 +5.015327e-08 *pow(time,4)
      // 	 -1.624539e-10 *pow(time,5)
      // 	 +1.414059e-13 *pow(time,6);
      
      //  fixed threshold 5% with 10 ns offset 90/10 and 80/20 (may 2009)
      drift = 2.498729e-02    
	-2.788627e-03 *time  
	+3.997619e-03 *pow(time,2)  
	-1.289555e-04 *pow(time,3)
	+2.101670e-06 *pow(time,4)
	-1.959622e-08 *pow(time,5)
	+1.053198e-10 *pow(time,6)
	-3.030411e-13 *pow(time,7)
	+3.609173e-16 *pow(time,8);
    }
    else{
      // fixed threshold 80/20
      drift = -4.843174e-02  
	+4.859828e-03 *time    
	+4.593485e-03 *pow(time,2)     
	-1.677031e-04 *pow(time,3)  
	+2.955270e-06 *pow(time,4)  
	-2.996812e-08 *pow(time,5)  
	+1.832841e-10 *pow(time,6)  
	-6.678677e-13 *pow(time,7)  
	+1.335729e-15 *pow(time,8)  
	-1.128932e-18 *pow(time,9);  
    }
  }

  drift =  0.1*drift;  
  if(drift < 0.) drift=0.;
  return drift;  

}


// --------------------------------------------------------------------------------
Double_t PndSttSingleStraw::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
  Out4 = StrawTot();              // find the time over threshold

  return PulseTime;
}


// --------------------------------------------------------------------------------   
Double_t  PndSttSingleStraw::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  PndSttSingleStraw::FastPartToADC() {

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

  // number of elecrons. Wi is the mean energy lost per free
  // electrn in eV
  Int_t NtotEle = (Int_t) (1.e+09 * STUrban()/Wi);

  for(Int_t j=1; j<= NtotEle; j++){
   
      ADCsignal += bPolya * GasGain * PolyaSamp();

  }

  return ADCsignal;
}


// ----------------------------------------------------------------------------------   
Double_t  PndSttSingleStraw::FastRec(Double_t TrueDcm, Int_t Flag) {

  // having the true distance TrueDcm (cm) as input,
  // return the reconstructed distance (cm), in a fast
  // and approximated way
  // by sampling on the simulated reconstruction curve
  // When Press =2 and Flag=1 one uses the julich experimental data
  // A. Rotondi March 2007
  // When Press = 2 and Flag=2 one uses the juelich experimental data from COSY-TOF (Feb 2011)
  // When Press = 2 and Flag=3 one uses the juelich experimental data from the small setup
  // (Apr 2010)
  

  Double_t  resmic;
  
  // 1 atm pressure
  if(pSTP < 1.9) {
    if(Radius < 0.45){
      if(TrueDcm < 0.38){
	resmic =   1.24506e+02  -1.80117e+02*TrueDcm
	  +3.76905e+03*pow(TrueDcm,2) -4.63251e+04*pow(TrueDcm,3)  
	  +1.80068e+05*pow(TrueDcm,4)   -2.21094e+05*pow(TrueDcm,5);
      }
      else resmic = 57.;
    }
    // radius > 0.4 cm
    else{   
      //      if(TrueDcm < 0.48){
      resmic =     1.53656e+02  
	-5.07749e+03*TrueDcm 
	+1.73707e+05*pow(TrueDcm,2)  
	-2.72285e+06*pow(TrueDcm,3)   
	+2.28719e+07*pow(TrueDcm,4)   
	-1.12921e+08*pow(TrueDcm,5) 
	+3.39427e+08*pow(TrueDcm,6) 
	-6.12741e+08 *pow(TrueDcm,7) 
	+6.12041e+08 *pow(TrueDcm,8)  
	-2.60444e+08*pow(TrueDcm,9) ;  
      //       }
      //       else resmic = 72.;
    }
  }
  
  // 2 atm pressure radius 5 cm
  else {
    if(Flag==0) {
      // data from simulation
      //       if(TrueDcm < 0.48){
      resmic =      +1.06966e+02  
	-4.03073e+03 *TrueDcm 
	+1.60851e+05 *pow(TrueDcm,2)  
	-2.87722e+06 *pow(TrueDcm,3)   
	+2.67581e+07 *pow(TrueDcm,4)   
	-1.43397e+08 *pow(TrueDcm,5) 
	+4.61046e+08 *pow(TrueDcm,6) 
	-8.79170e+08 *pow(TrueDcm,7) 
	+9.17095e+08 *pow(TrueDcm,8)  
	-4.03253e+08 *pow(TrueDcm,9) ;  
      //       }
      //       else resmic=30.;
    }
    else if (Flag==1) {
      // data from COSY-TOF (old)
      //       if(TrueDcm < 0.48){
      resmic =    +1.48048e+02
	-3.35951e+02*TrueDcm 
	-1.87575e+03*pow(TrueDcm,2)  
	+1.92910e+04*pow(TrueDcm,3)   
	-6.90036e+04*pow(TrueDcm,4)   
	+1.07960e+05*pow(TrueDcm,5) 
	-5.90064e+04*pow(TrueDcm,6) ;  
      //     }
      //       else resmic=65.;
    }
    else if (Flag==2) {
      // data from COSY-TOF (Feb 2011)
      // the parametrization comes from mm vs mm:
      // => TrueDcm must be in mm
      TrueDcm *= 10.; // cm -> mm
      // and resmic will be given in mm ...

      // old parametriz
      //     resmic = +0.02152
      // 	+0.6764*TrueDcm
      // 	-1.008*pow(TrueDcm,2)
      // 	+0.7421*pow(TrueDcm,3)
      // 	-0.3036*pow(TrueDcm,4)
      // 	+0.06955*pow(TrueDcm,5)
      // 	-0.008327*pow(TrueDcm,6)
      // 	+0.0004049*pow(TrueDcm,7);

      // pol5 parametriz
      resmic = 0.188119 
	+ 0.00211993 * TrueDcm
	+ 0.00336004 * pow(TrueDcm, 2)
	- 0.0103979  * pow(TrueDcm, 3)
	+ 0.0033461  * pow(TrueDcm, 4)
	-0.000315764 * pow(TrueDcm, 5);
      
   // convert resmic to micron and TrueDcm to cm
      resmic *= 1000.;
      TrueDcm *= 0.1;
    }
    else {
      // data from juelich prototype
      resmic = 4.521331e-01
	-2.087216e-01 *10.*TrueDcm
	+4.911102e-02 *pow(10.*TrueDcm,2)
	-3.934728e-03 *pow(10.*TrueDcm,3);
      resmic = resmic*1000.;
    }
  }
  
  //real distance in cm
  Double_t rsim = gRandom->Gaus(TrueDcm, resmic*0.0001);
  if (rsim<0.) rsim = TrueDcm - TrueDcm*gRandom->Uniform(0.,1.);
  else if (rsim>0.5) rsim = TrueDcm + (0.5-TrueDcm)*gRandom->Uniform(0.,1.);
  
  return rsim;
}


// --------------------------------------------------------------------------------   
Double_t  PndSttSingleStraw::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  PndSttSingleStraw::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  PndSttSingleStraw::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

}