ROOT logo
//*-- Author  : A.Ivashkin <ivashkin@inr.ru>
//*-- Modified: A.Sadovsky <sadovski@fz-rossendorf.de> (04.11.2004 Hydra implementation)
//*-- Last modified : 14/08/2005 by Ilse Koenig

using namespace std;
#include <iostream>
#include <fstream>
#include "hrungekutta.h"
#include "hmdcsizescells.h"
#include <math.h>
#include "TMatrixD.h"
#include "TDecompLU.h"
ClassImp(HRungeKutta)
		
HRungeKutta::HRungeKutta() {

  // constructor
  clear();

  // Eventually these parameters should be stored in a parameter container!
  initialStepSize=50.;
  stepSizeDec=0.75;
  stepSizeInc=1.25;
  maxStepSize=200.;
  minStepSize=10.;
  minPrecision=2.e-04;
  maxPrecision=2.e-05;
  maxNumSteps=1000;
  maxDist=2.5;
  dpar[0]=0.01; dpar[1]=0.001; dpar[2]=0.001; dpar[3]=0.01; dpar[4]=0.01;
  resolution[0]=0.140*2.00;  resolution[1]=0.140*1.00; // MDC1 intrinsic resolution
  resolution[2]=0.140*2.00;  resolution[3]=0.140*1.00; // MDC2
  resolution[4]=0.140*2.00;  resolution[5]=0.140*1.00; // MDC3
  resolution[6]=0.140*2.00;  resolution[7]=0.140*1.00; // MDC4
  for(Int_t i=0;i<8;i++) {sig[i]=resolution[i];}
  mTol=DBL_EPSILON;
}


void HRungeKutta::clear() {
  // clears magnetic field definition and vertex position and direction
  fieldMap=0;
  fieldScalFact=0.;
  for(Int_t m=0; m<4; m++) {
    for(Int_t s=0; s<6; s++) {
      mdcInstalled[m][s]=kFALSE;
      mdcPos[m][s][0]=0.;
      mdcPos[m][s][1]=0.;
      mdcPos[m][s][2]=0.;
      normVecMdc[m][s][0]=0.;
      normVecMdc[m][s][1]=0.;
      normVecMdc[m][s][2]=0.;
    }       
  }       
}


void HRungeKutta::setField(HMdcTrackGField* pField, Float_t fpol) {
  // sets the field map and the scaling factor
  fieldMap=pField;
  fieldScalFact=fpol;
}


void HRungeKutta::setMdcPosition(Int_t s, Int_t m, HGeomTransform& gtrMod2Sec) {
  // sets the MDC module position in the sector coordinate system
  // calculates the normal vector on the planes
  mdcInstalled[m][s]=kTRUE;
  HGeomVector aTransl=gtrMod2Sec.getTransVector();
  mdcPos[m][s][0] = aTransl.getX();
  mdcPos[m][s][1] = aTransl.getY();
  mdcPos[m][s][2] = aTransl.getZ();
  HGeomVector r0_mod(0.0, 0.0, 0.0);
  HGeomVector rz_mod(0.0, 0.0, 1.0);
  HGeomVector norm= gtrMod2Sec.transFrom(rz_mod) - gtrMod2Sec.transFrom(r0_mod);
  norm/=norm.length();  // to be shure it is normalized after transformation
  normVecMdc[m][s][0]=norm.getX();
  normVecMdc[m][s][1]=norm.getY();
  normVecMdc[m][s][2]=norm.getZ();
  if (m==0 || m==1) {
    dzfacMdc[m][s]=tan(acos(norm.getZ()));
    HGeomVector v(0.0, 0.0, -500.0);
    if (m==1 && mdcInstalled[0][s]==kFALSE) v.setZ(-650.);
    HGeomVector vs=gtrMod2Sec.transFrom(v);
    pointForVertexRec[s][0]=vs.getX();   
    pointForVertexRec[s][1]=vs.getY();   
    pointForVertexRec[s][2]=vs.getZ();
  }
}


Bool_t HRungeKutta::fit4Hits(Double_t* u1pos, Double_t* u1dir,
                             Double_t* u2pos, Double_t* u2dir,
                             Float_t* sigMult, Int_t sec, Double_t pGuess) {
  // Input:  position and direction of the two segments u1 and u2,
  //         array with multiplication factors for errors
  //         sector number, initial guess for momentum
  // From the segments the hit positions are calculated.
  // Calls fitMdc for Runge-Kutta track propagation and fitting
  nMaxMod=4;
  ndf = 3; //{x,y}*4mdc=8 - {p,x,y,dx,dy}=5 == 3
  sector=sec;
  p0Guess = pGuess;
  findMdcIntersectionPoint(u1pos,u1dir,0,&hit[0][0]);  // hit = hit position in modules
  findMdcIntersectionPoint(u1pos,u1dir,1,&hit[1][0]);
  findMdcIntersectionPoint(u2pos,u2dir,2,&hit[2][0]);
  findMdcIntersectionPoint(u2pos,u2dir,3,&hit[3][0]);
  for(Int_t i=0;i<4;i++) {
    mdcModules[i]=i;
    mdcLookup[i]=i;
    sig[2*i]=sigMult[2*i]*resolution[2*i];
    sig[2*i+1]=sigMult[2*i+1]*resolution[2*i+1];
  }
  return fitMdc();
}


Bool_t HRungeKutta::fit3Hits(Double_t* u1pos, Double_t* u1dir,
                            Double_t* u2pos, Double_t* u2dir,
                            Float_t* sigMult, Int_t sec, Double_t pGuess) {
  // Input:  position and direction of the two segments u1 and u2
  //         array with multiplication factors for errors
  //         sector number, initial guess for momentum
  // From the segments the hit positions are calculated.
  // If the multiplication factors for errors of a hit -1, the hit is discarded.
  // Calls fitMdc for Runge-Kutta track propagation and fitting
  nMaxMod=3;
  ndf = 1; //{x,y}*3mdc=6 - {p,x,y,dx,dy}=5 == 1
  sector=sec;
  p0Guess = pGuess;
  Int_t m=0;
  // hit = hit position in modules
  for (Int_t i=0;i<2;i++) {
    if (sigMult[2*i]>0) {
      findMdcIntersectionPoint(u1pos,u1dir,i,&hit[m][0]);
      mdcModules[m]=i;
      mdcLookup[i]=m;
      sig[2*m]=sigMult[2*i]*resolution[2*i];
      sig[2*m+1]=sigMult[2*i+1]*resolution[2*i+1];
      m++;
    } else mdcLookup[i]=-1;      
  }
  for (Int_t i=2;i<4;i++) {
    if (sigMult[2*i]>0) {
      findMdcIntersectionPoint(u2pos,u2dir,i,&hit[m][0]);
      mdcModules[m]=i;
      mdcLookup[i]=m;
      sig[2*m]=sigMult[2*i]*resolution[2*i];
      sig[2*m+1]=sigMult[2*i+1]*resolution[2*i+1];
      m++;
    } else mdcLookup[i]=-1;      
  }
  mdcModules[3]=-1;
  return fitMdc();
}


Bool_t HRungeKutta::fitMdc() {
  initMom(); // sets pfit
  if (fabs(pfit)<10.) return kFALSE;	

  Bool_t rc=kTRUE;
  Float_t chi2old=0.F;
  for(Int_t i=0;i<3;i++) fit[0][i]=hit[0][i];  // fit = fitted position in first module
  Float_t dist=distance(&hit[0][0],&hit[1][0]);
  if (dist==0.) { return kFALSE; }
  Float_t uh=(hit[1][0]-hit[0][0])/dist;       // cos(phi)
  Float_t uv=(hit[1][1]-hit[0][1])/dist;       // cos(theta)
  for(Int_t iter=0;iter<11;iter++) {
    polbending=fieldScalFact*pfit;
    chi2=0.;
    parSetNew0(uh,uv,iter);
    if (jstep>=maxNumSteps) { return kFALSE; }
    for(Int_t m=0;m<nMaxMod;m++) {
      for(Int_t i=0;i<2;i++) {
        Float_t chi=(hit[m][i]-fit[m][i])/sig[(m*2)+i];
        chi2+=chi*chi;
      }
    }
    if (finite(pfit)==0||finite(chi2)==0) {
      Warning("HRungeKuttta::fitMdc()","NaN or Inf value received for pFit or chi2...skipped!");
      return kFALSE;
    }
    if (chi2==0.) chi2=0.001;
    if ((TMath::Abs(chi2-chi2old)/chi2)<0.005||iter==10||(iter>5&&chi2old>1.e+06&&chi2>chi2old)) break;

    chi2old=chi2;
    parSetNewVar(uh,uv);
    if (jstep>=maxNumSteps) { return kFALSE; }

    derive(&fitdp[0][0],0);
    derive(&fitdh[0][0],1);	
    derive(&fitdv[0][0],2);
    derive(&fit1h[0][0],3);
    derive(&fit1v[0][0],4);

    for(Int_t i=0;i<nMaxMod;i++) {
      for(Int_t k=0;k<2;k++) {
        ydata[i*2+k]=fit[i][k]-hit[i][k];
      }
    }

    if (!(rc=linSys())) break;
    Double_t dpp=fitpar[0]/pfit;
    if ( dpp>=0) {
      pfit+=fitpar[0]*(1+dpp);
    } else {
      pfit+=fitpar[0]/(1-dpp);
    }
    uh+=fitpar[1];
    uv+=fitpar[2];
    fit[0][0]+=fitpar[3];
    fit[0][1]+=fitpar[4];
    fit[0][2]-=fitpar[4]*dzfacMdc[(mdcModules[0])][sector];
    if (pfit==0.) pfit=10.;
  }
  return rc;
}


void HRungeKutta::findMdcIntersectionPoint(Double_t* upos,Double_t* udir,Int_t mdcN,Double_t* xyz) {
  // calculates the hit points in the MDC module (index mdcN 0..3)
  // input:  posline == pos from track piece 
  // input:  dirline == pos from segment 
  //         mdcN    == module number 0..3
  // output: xyz     == hit position in module
  Double_t* normplane=&normVecMdc[mdcN][sector][0];
  Double_t* pplane=&mdcPos[mdcN][sector][0];
  findIntersectionPoint(upos,udir,pplane,normplane,xyz);
}


void HRungeKutta::initMom() {
  // Here we set initial approximation momentum value
  if (p0Guess!=999999999.) {
    // User supplied momentum guess, so we use it as 0-approximation
    pfit=p0Guess;
  } else {
    // User did not supply a momentum guess. We start from simple-minded kickplane,
    // which needs the tracklet directions on input.
    Double_t dir1[3],dir2[3];
    for(Int_t i=0;i<3;i++) {
      dir1[i]=hit[1][i]-hit[0][i];
      dir2[i]=hit[3][i]-hit[2][i];
    }
    Float_t temp=sqrt((dir1[1]*dir1[1]+dir1[2]*dir1[2])*(dir2[1]*dir2[1]+dir2[2]*dir2[2]));
    if (temp==0.) temp=1.;	
    Float_t cosxi=(dir1[1]*dir2[1]+dir1[2]*dir2[2])/temp;
    if (cosxi> 1.0) cosxi=1;
    if (cosxi<-1.0) cosxi=-1;
    Float_t xi=acos(cosxi);
    Float_t vprod=(dir1[1]*dir2[2]-dir1[2]*dir2[1]);
    temp=sin(xi/2.);
    if (temp==0.) temp=0.001;
    Float_t temp1=TMath::Abs(vprod);
    if (temp1==0.) temp1=0.001;
    pfit=75./(2*temp)*vprod/temp1;
  }
}
		

Float_t HRungeKutta::distance(Double_t* p1,Double_t* p2) {
  // calculates the distance between two points
  Float_t dist=sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+
                    (p1[1]-p2[1])*(p1[1]-p2[1])+
                    (p1[2]-p2[2])*(p1[2]-p2[2]));
  return dist;
}		


void HRungeKutta::parSetNew0(Float_t uh,Float_t uv,Int_t iter) {
  // RK propagation with initial or last fit parameters
  // uv - vertical direction in MDC1
  // uh - horizontal direction
  Double_t upos[3], udir[3];
  for(Int_t i=0;i<3;i++) {
    upos[i]=fit[0][i];     // position of track on MDC1
  }
  cosines(uh,uv,udir);
  if (iter>0) {
    // start with point on plane 
    findMdcIntersectionPoint(upos,udir,mdcModules[0],&fit[0][0]);
    for(Int_t i=0;i<3;i++) { upos[i]=fit[0][i]; }
  }
  for(Int_t i=0;i<3;i++) {
    dirAtFirstMDC[i]=udir[i];  // direction of track on first MDC
  }
  trackLength=0.;
  gentrkNew(pfit,&upos[0],&udir[0],&fit[0][0],1);
}


void HRungeKutta::parSetNewVar(Float_t uh,Float_t uv) {
  // calculates the derivatives (RK propagation with parameter variations)
  // uv - vertical direction
  // uh - horizontal direction in MDC1
  Double_t upos[3], udir[3];
  upos[0]=fit[0][0];
  upos[1]=fit[0][1];
  upos[2]=fit[0][2];
  cosines(uh,uv,udir);
  gentrkNew(pfit+dpar[0],&upos[0],&udir[0],&fitdp[0][0],2); // momentum additive

  upos[0]=fit[0][0];
  upos[1]=fit[0][1];
  upos[2]=fit[0][2];
  cosines(uh+dpar[1],uv,udir);         // x direction additive
  gentrkNew(pfit,&upos[0],&udir[0],&fitdh[0][0],3);

  upos[0]=fit[0][0];
  upos[1]=fit[0][1];
  upos[2]=fit[0][2];
  cosines(uh,uv+dpar[2],udir);         // y direction additive
  gentrkNew(pfit,&upos[0],&udir[0],&fitdv[0][0],4);

  upos[0]=fit[0][0]+dpar[3];           // x position additive
  upos[1]=fit[0][1];
  upos[2]=fit[0][2];
  cosines(uh,uv,udir);
  gentrkNew(pfit,&upos[0],&udir[0],&fit1h[0][0],5);

  upos[0]=fit[0][0];
  upos[1]=fit[0][1]+dpar[4];                   // y position additive
  upos[2]=fit[0][2]-dpar[4]*dzfacMdc[(mdcModules[0])][sector]; // change in the first MDC plane
  cosines(uh,uv,udir);
  gentrkNew(pfit,&upos[0],&udir[0],&fit1v[0][0],6);
}


void HRungeKutta::cosines(Float_t uh, Float_t uv, Double_t* dir) {
  // calculates the direction
  Float_t prod=sqrt(uh*uh+uv*uv);
  if (prod==0.) prod=0.0001;
  if (prod<1.) { 
    dir[0]=uh;
    dir[1]=uv;
    dir[2]=sqrt(1.-prod*prod);
  } else {
    dir[0]=uh/prod;
    dir[1]=uv/prod;
    dir[2]=0.;
  }	
}


void HRungeKutta::derive(Double_t* hitpar, Int_t kind) {
  // calculates the derivative (hit-fit)/variation
  Double_t dp=dpar[kind];
  for(Int_t i=0;i<nMaxMod;i++) {
    Double_t* h=&hitpar[i*3];
    for(Int_t k=0;k<2;k++) {
      dydpar[kind][i*2+k]=(h[k]-fit[i][k])/dp;
    }
  }
}		


Bool_t HRungeKutta::linSys() {
  // solves the system of linear equations and sets the new fitting parameters
  Double_t data[5][5], elements[5];
  Int_t kMax=nMaxMod*2;
  for(Int_t i=0;i<5;i++) {	  
    for(Int_t j=0;j<5;j++) {
      data[j][i]=0.;
      for(Int_t k=0;k<kMax;k++) {
        data[j][i]+=dydpar[j][k]*dydpar[i][k]/(sig[k]*sig[k]);
      } 	 
    }	  	
  }
  for(Int_t i=0;i<5;i++) {
    elements[i]=0.;
    for(Int_t k=0;k<kMax;k++) {
      elements[i]-= ydata[k]*dydpar[i][k]/(sig[k]*sig[k]);
    }
  }
		

  //-------------------------------------------
  // test if matrix is singular
  TMatrixD test(5,5,(Double_t*)data,NULL);
  Double_t sign;
  Int_t nrZeros;
  if(!decompose(test,sign,mTol,nrZeros)||nrZeros>0){return kFALSE;}
  //-------------------------------------------

  TMatrixD m(5,5,(Double_t*)data,NULL);

  TMatrixD b(5,1,(Double_t*)elements,NULL);
  TMatrixD m1(5,1);
  m.Invert(0);
  m1.Mult(m,b);
  
  for(Int_t i=0;i<5;i++) {
    fitpar[i]=m1(i,0);
  }

  return kTRUE;
}


void HRungeKutta::gentrkNew(Float_t p, Double_t* upos, Double_t* udir, Double_t* hitxyz, Int_t kind ) {
  // propages the track through the MDC system
  // p      == momentum
  // upos   == initial position
  // udir   == initial direction
  // hitxyz == hits in chambers
  // kind   == type of variation
  Double_t* pHit=&hitxyz[0];
  for(Int_t i=0;i<3;i++) pHit[i]=upos[i];
  jstep=0;
  Double_t d=0.;	
  Float_t step=initialStepSize;
  Float_t nextStep=step;
  Float_t stepFac=1.;
  for(Int_t nMod=1;nMod<nMaxMod;nMod++) {
    pHit=&hitxyz[nMod*3];
    do {	
      stepFac=rkgtrk(step,p,upos,udir,kind);
      findMdcIntersectionPoint(upos,udir,mdcModules[nMod],pHit);
      d=distance(upos,pHit);
      if (step<nextStep) {  // last step was rest step to chamber
        if (stepFac<1.) nextStep=step*stepFac;
      } else nextStep*=stepFac;
      if (d>nextStep||d<maxDist) step=nextStep;
      else step=d;
    } while (d>=maxDist&&((polbending>=0&&upos[2]<pHit[2])||(polbending<0&&upos[1]<pHit[1]))&&jstep<maxNumSteps);
  }
  if (kind==1) {
    // we store the track position after the last RK step for subsequent propagation
    // to Meta (this step ends near the chamber, while the position on the chamber is
    // stored in fit[3])
    for(Int_t i=0;i<3;i++) {
      posNearLastMDC[i] = upos[i];
      dirAtLastMDC[i] = udir[i]; // track direction
    }
    stepSizeAtLastMDC=nextStep;
  }
}


Float_t HRungeKutta::rkgtrk(Float_t totStep, Float_t p, Double_t* upos, Double_t* udir, Int_t kind) {
  // one step of track tracing from the position upos in the direction of udir
  Float_t  k1[3],k2[3],k3[3],k4[3],hstep,qstep;
  Double_t u1pos[3], u1dir[3], u2pos[3], u2dir[3], ucm[3], b[3];
  Double_t est=0.;
  Float_t step=totStep;
  Float_t stepFac=1.;
  for(Int_t i=0;i<3;i++) ucm[i]=upos[i]/10.;		
  fieldMap->calcField(ucm,b,fieldScalFact);
  rkeqfw(b,p,udir,k1);

  do {
    jstep++;
    hstep=step*0.5;
    qstep=step*step*0.125;

    for(Int_t i=0;i<3;i++) {
      u1pos[i]=upos[i]+hstep*udir[i]+qstep*k1[i];
      u1dir[i]=udir[i]+hstep*k1[i];
    }
    for(Int_t i=0;i<3;i++) ucm[i]=u1pos[i]/10.;		
    fieldMap->calcField(ucm,b,fieldScalFact);
    rkeqfw(b,p,u1dir,k2);

    for(Int_t i=0;i<3;i++) {
      u1dir[i] = udir[i] + hstep*k2[i];
    }
    rkeqfw(b,p,u1dir,k3);	  

    for(Int_t i=0;i<3;i++) {
      u2pos[i]=upos[i]+step*udir[i]+2.*hstep*hstep*k3[i];
      u2dir[i]=udir[i]+step*k3[i];
    }
    for(Int_t i=0;i<3;i++) ucm[i]=u2pos[i]/10.;		
    fieldMap->calcField(ucm,b,fieldScalFact);
    rkeqfw(b,p,u2dir,k4);

    est=0.;
    for(Int_t i=0;i<3;i++) {
      est+=fabs(k1[i]+k4[i]-k2[i]-k3[i])*hstep;
    }
    if (est<minPrecision||step<=minStepSize) {
      break; 
    } else {
      stepFac*=stepSizeDec;
      step*=stepSizeDec;
    } 
  } while (jstep<maxNumSteps);
  for(Int_t i=0;i<3;i++) {
    upos[i]=upos[i]+step*udir[i]+step*step*(k1[i]+k2[i]+k3[i])/6.;
    udir[i]=udir[i]+step*(k1[i]+2.*k2[i]+2.*k3[i]+k4[i])/6.;
  }
  if (kind==1) {
    trackLength+=step;  // calculate track length
  }
  if (est<maxPrecision&&step<maxStepSize) {
    stepFac*=stepSizeInc;
  }
  return stepFac;
}


void HRungeKutta::rkeqfw(Double_t b[3], Float_t p, Double_t* udir, Float_t* duds) {
  // calculates 2nd derivative
  Double_t  con=0.2997925;
  duds[0]=con*(udir[1]*b[2]-udir[2]*b[1])/p;
  duds[1]=con*(udir[2]*b[0]-udir[0]*b[2])/p;
  duds[2]=con*(udir[0]*b[1]-udir[1]*b[0])/p;	
}	
				

void HRungeKutta::traceToVertex(HMdcSizesCells* pMSizesCells) {
  // propogates the track backwards from first MDC till Z=zVertex plane
  // should be called only after RK-minimisation is done
  Float_t step=initialStepSize;
  Float_t stepFac=1.;
  jstep=0;
  Double_t dist=0.;
  Double_t* normplane=&normVecMdc[(mdcModules[0])][sector][0];
  Double_t* pplane=&pointForVertexRec[sector][0];
  Double_t upos[3], udir[3], xyz[3];
  for(Int_t i=0;i<3;i++) {
    upos[i]=fit[0][i];     // track position at first MDC
    udir[i]=-dirAtFirstMDC[i]; // direction (reversed to go back from first MDC to target)
  }
  do {
    findIntersectionPoint(upos,udir,pplane,normplane,&xyz[0]);
    dist=distance(upos,xyz);
    if (dist<maxDist||upos[2]<=xyz[2]) {
      break;
    }
    if (dist<step) { step=dist; }
    stepFac=rkgtrk(step,-pfit,upos,udir,1); // with reversed momentum for correct track curvature
    step*=stepFac;
  } while (jstep<maxNumSteps);
  HMdcSizesCells::calcMdcSeg(upos[0],upos[1],upos[2],
        upos[0]-udir[0],upos[1]-udir[1],upos[2]-udir[2],
        zSeg1,rSeg1,thetaSeg1,phiSeg1);
  HMdcSizesCellsSec& fSCSec=(*pMSizesCells)[sector];
  const HGeomVector& target=fSCSec.getTargetMiddlePoint(); 
  xyz[0]=target.getX();
  xyz[1]=target.getY();
  xyz[2]=target.getZ();
  dist=distance(xyz,upos);
  trackLength+=dist;
  trackLengthAtLastMDC=trackLength;
}


void HRungeKutta::traceToMETA(HGeomVector& metaHit,HGeomVector& metaNormVec,
                              HGeomVector* pathCorrPos,HGeomVector* pathCorrNorm) {
  // propagates the track from the last fitted MDC position till track intersection
  // with META plane (is defined by META-hit and normal vector to rod's plane
  // should be called only after RK-minimisation is done
  trackLength=trackLengthAtLastMDC;
  Double_t  pointOnMETA[3], normVectorOfMETA[3];
  pointOnMETA[0]=metaHit.getX();
  pointOnMETA[1]=metaHit.getY();
  pointOnMETA[2]=metaHit.getZ();
  normVectorOfMETA[0] = metaNormVec.getX();
  normVectorOfMETA[1] = metaNormVec.getY();
  normVectorOfMETA[2] = metaNormVec.getZ();
   
  Float_t step = stepSizeAtLastMDC;
  Float_t stepFac=1.;
  jstep = 0;
  Double_t upos[3], udir[3];
  for(Int_t i=0;i<3;i++) {
    upos[i]=posNearLastMDC[i];   // tracking position near last MDC
    udir[i]=dirAtLastMDC[i];     // and direction
  }
  Double_t dist=0.;
  Bool_t rc=kTRUE;
  do {
    rc=findIntersectionPoint(upos,udir,pointOnMETA,normVectorOfMETA,&trackPosOnMETA[0]);
    if (!rc) break;
    dist=distance(upos,trackPosOnMETA);
    if ((polbending>=0&&upos[2]>trackPosOnMETA[2])||(polbending<0&&upos[1]>trackPosOnMETA[1])) {
      trackLength-=dist;
      break;
    } else if (dist<maxDist) {
      trackLength+=dist;
      break;
    } else if (upos[1]>=2475.||upos[1]<=0.) {  
      trackLength=0.;
      break;
    }
    if (dist<step) { step=dist; }
    stepFac=rkgtrk(step,pfit,upos,udir,1); // we use pfit as value for momentum to propagate
    step*=stepFac;
  } while (jstep<maxNumSteps);
  if (!rc&&jstep>=maxNumSteps) {
    trackLength=0.;
  }
  HMdcSizesCells::calcMdcSeg(upos[0],upos[1],upos[2],
        upos[0]+udir[0],upos[1]+udir[1],upos[2]+udir[2],
        zSeg2,rSeg2,thetaSeg2,phiSeg2);
  if (pathCorrPos!=0&&pathCorrNorm!=0&&trackLength>0.) {
    Double_t p1[3], pn[3], p2[3];
    p1[0]=pathCorrPos->getX();
    p1[1]=pathCorrPos->getY();
    p1[2]=pathCorrPos->getZ();
    pn[0]=pathCorrNorm->getX();
    pn[1]=pathCorrNorm->getY();
    pn[2]=pathCorrNorm->getZ();
    findIntersectionPoint(trackPosOnMETA,udir,p1,pn,&p2[0]);
    trackLength-=distance(trackPosOnMETA,p2);
  }
}


Bool_t HRungeKutta::findIntersectionPoint(Double_t* upos, Double_t* udir,
                                          Double_t* planeR0, Double_t* planeNorm,
                                          Double_t* pointIntersect) {
  // upos, udir - "moving point" current position and movement direction vector
  // planeR0[3] some point on the plane to which intersection is supposed to happen
  // planeNorm[3] perpendicular vector to the plane
  // pointIntersect[3] - output - the intersection point of "moving point" with plane
  Double_t denom = (planeNorm[0]*udir[0] + planeNorm[1]*udir[1] + planeNorm[2]*udir[2]);
  if (denom!=0.0) {
    Double_t t = ((planeR0[0]-upos[0])*planeNorm[0] + 
                  (planeR0[1]-upos[1])*planeNorm[1] + 
                  (planeR0[2]-upos[2])*planeNorm[2]) / denom;
    for(Int_t i=0;i<3;i++) pointIntersect[i] = upos[i] + udir[i]*t;  
    return kTRUE;
  } else {
    Warning("findIntersectionPoint","No intersection point found : (plane || track)");
    return kFALSE;
  }
}


Bool_t HRungeKutta::decompose(TMatrixD &lu,Double_t &sign,Double_t tol,Int_t &nrZeros)
{
    // THis part is copid from ROOT TDecompLU::Decompose() to get
    // rid of error message!
    // Crout/Doolittle algorithm of LU decomposing a square matrix, with implicit partial
    // pivoting.  The decomposition is stored in fLU: U is explicit in the upper triag
    // and L is in multiplier form in the subdiagionals .
    // Row permutations are mapped out in fIndex. fSign, used for calculating the
    // determinant, is +/- 1 for even/odd row permutations. .

  // copied part from  TDecompLU constructor-------

  Int_t  fNIndex = lu.GetNcols();
  Int_t* index = new Int_t[fNIndex];
  memset(index,0,fNIndex*sizeof(Int_t));

  const Int_t kWorkMax = 100; // size of work array's in several routines
  //-----------------------------------------------

  const Int_t     n     = lu.GetNcols();
        Double_t *pLU   = lu.GetMatrixArray();

  Double_t work[kWorkMax];
  Bool_t isAllocated = kFALSE;
  Double_t *scale = work;
  if (n > kWorkMax) {
    isAllocated = kTRUE;
    scale = new Double_t[n];
  }

  sign    = 1.0;
  nrZeros = 0;
  // Find implicit scaling factors for each row
  for (Int_t i = 0; i < n ; i++) {
    const Int_t off_i = i*n;
    Double_t max = 0.0;
    for (Int_t j = 0; j < n; j++) {
      const Double_t tmp = TMath::Abs(pLU[off_i+j]);
      if (tmp > max)
        max = tmp;
    }
    scale[i] = (max == 0.0 ? 0.0 : 1.0/max);
  }

  for (Int_t j = 0; j < n; j++) {
    const Int_t off_j = j*n;
    Int_t i;
    // Run down jth column from top to diag, to form the elements of U.
    for (i = 0; i < j; i++) {
      const Int_t off_i = i*n;
      Double_t r = pLU[off_i+j];
      for (Int_t k = 0; k < i; k++) {
        const Int_t off_k = k*n;
        r -= pLU[off_i+k]*pLU[off_k+j];
      }
      pLU[off_i+j] = r;
    }

    // Run down jth subdiag to form the residuals after the elimination of
    // the first j-1 subdiags.  These residuals divided by the appropriate
    // diagonal term will become the multipliers in the elimination of the jth.
    // subdiag. Find fIndex of largest scaled term in imax.

    Double_t max = 0.0;
    Int_t imax = 0;
    for (i = j; i < n; i++) {
      const Int_t off_i = i*n;
      Double_t r = pLU[off_i+j];
      for (Int_t k = 0; k < j; k++) {
        const Int_t off_k = k*n;
        r -= pLU[off_i+k]*pLU[off_k+j];
      }
      pLU[off_i+j] = r;
      const Double_t tmp = scale[i]*TMath::Abs(r);
      if (tmp >= max) {
        max = tmp;
        imax = i;
      }
    }

    // Permute current row with imax
    if (j != imax) {
      const Int_t off_imax = imax*n;
      for (Int_t k = 0; k < n; k++ ) {
        const Double_t tmp = pLU[off_imax+k];
        pLU[off_imax+k] = pLU[off_j+k];
        pLU[off_j+k]    = tmp;
      }
      sign = -sign;
      scale[imax] = scale[j];
    }
    index[j] = imax;

    // If diag term is not zero divide subdiag to form multipliers.
    if (pLU[off_j+j] != 0.0) {
      if (TMath::Abs(pLU[off_j+j]) < tol)
        nrZeros++;
      if (j != n-1) {
        const Double_t tmp = 1.0/pLU[off_j+j];
        for (Int_t i = j+1; i < n; i++) {
          const Int_t off_i = i*n;
          pLU[off_i+j] *= tmp;
        }
      }
    } else {
      if (isAllocated) delete [] scale;
      delete [] index;
      return kFALSE;
    }
  }

  if (isAllocated)
    delete [] scale;

  delete [] index;
  return kTRUE;
}

 hrungekutta.cc:1
 hrungekutta.cc:2
 hrungekutta.cc:3
 hrungekutta.cc:4
 hrungekutta.cc:5
 hrungekutta.cc:6
 hrungekutta.cc:7
 hrungekutta.cc:8
 hrungekutta.cc:9
 hrungekutta.cc:10
 hrungekutta.cc:11
 hrungekutta.cc:12
 hrungekutta.cc:13
 hrungekutta.cc:14
 hrungekutta.cc:15
 hrungekutta.cc:16
 hrungekutta.cc:17
 hrungekutta.cc:18
 hrungekutta.cc:19
 hrungekutta.cc:20
 hrungekutta.cc:21
 hrungekutta.cc:22
 hrungekutta.cc:23
 hrungekutta.cc:24
 hrungekutta.cc:25
 hrungekutta.cc:26
 hrungekutta.cc:27
 hrungekutta.cc:28
 hrungekutta.cc:29
 hrungekutta.cc:30
 hrungekutta.cc:31
 hrungekutta.cc:32
 hrungekutta.cc:33
 hrungekutta.cc:34
 hrungekutta.cc:35
 hrungekutta.cc:36
 hrungekutta.cc:37
 hrungekutta.cc:38
 hrungekutta.cc:39
 hrungekutta.cc:40
 hrungekutta.cc:41
 hrungekutta.cc:42
 hrungekutta.cc:43
 hrungekutta.cc:44
 hrungekutta.cc:45
 hrungekutta.cc:46
 hrungekutta.cc:47
 hrungekutta.cc:48
 hrungekutta.cc:49
 hrungekutta.cc:50
 hrungekutta.cc:51
 hrungekutta.cc:52
 hrungekutta.cc:53
 hrungekutta.cc:54
 hrungekutta.cc:55
 hrungekutta.cc:56
 hrungekutta.cc:57
 hrungekutta.cc:58
 hrungekutta.cc:59
 hrungekutta.cc:60
 hrungekutta.cc:61
 hrungekutta.cc:62
 hrungekutta.cc:63
 hrungekutta.cc:64
 hrungekutta.cc:65
 hrungekutta.cc:66
 hrungekutta.cc:67
 hrungekutta.cc:68
 hrungekutta.cc:69
 hrungekutta.cc:70
 hrungekutta.cc:71
 hrungekutta.cc:72
 hrungekutta.cc:73
 hrungekutta.cc:74
 hrungekutta.cc:75
 hrungekutta.cc:76
 hrungekutta.cc:77
 hrungekutta.cc:78
 hrungekutta.cc:79
 hrungekutta.cc:80
 hrungekutta.cc:81
 hrungekutta.cc:82
 hrungekutta.cc:83
 hrungekutta.cc:84
 hrungekutta.cc:85
 hrungekutta.cc:86
 hrungekutta.cc:87
 hrungekutta.cc:88
 hrungekutta.cc:89
 hrungekutta.cc:90
 hrungekutta.cc:91
 hrungekutta.cc:92
 hrungekutta.cc:93
 hrungekutta.cc:94
 hrungekutta.cc:95
 hrungekutta.cc:96
 hrungekutta.cc:97
 hrungekutta.cc:98
 hrungekutta.cc:99
 hrungekutta.cc:100
 hrungekutta.cc:101
 hrungekutta.cc:102
 hrungekutta.cc:103
 hrungekutta.cc:104
 hrungekutta.cc:105
 hrungekutta.cc:106
 hrungekutta.cc:107
 hrungekutta.cc:108
 hrungekutta.cc:109
 hrungekutta.cc:110
 hrungekutta.cc:111
 hrungekutta.cc:112
 hrungekutta.cc:113
 hrungekutta.cc:114
 hrungekutta.cc:115
 hrungekutta.cc:116
 hrungekutta.cc:117
 hrungekutta.cc:118
 hrungekutta.cc:119
 hrungekutta.cc:120
 hrungekutta.cc:121
 hrungekutta.cc:122
 hrungekutta.cc:123
 hrungekutta.cc:124
 hrungekutta.cc:125
 hrungekutta.cc:126
 hrungekutta.cc:127
 hrungekutta.cc:128
 hrungekutta.cc:129
 hrungekutta.cc:130
 hrungekutta.cc:131
 hrungekutta.cc:132
 hrungekutta.cc:133
 hrungekutta.cc:134
 hrungekutta.cc:135
 hrungekutta.cc:136
 hrungekutta.cc:137
 hrungekutta.cc:138
 hrungekutta.cc:139
 hrungekutta.cc:140
 hrungekutta.cc:141
 hrungekutta.cc:142
 hrungekutta.cc:143
 hrungekutta.cc:144
 hrungekutta.cc:145
 hrungekutta.cc:146
 hrungekutta.cc:147
 hrungekutta.cc:148
 hrungekutta.cc:149
 hrungekutta.cc:150
 hrungekutta.cc:151
 hrungekutta.cc:152
 hrungekutta.cc:153
 hrungekutta.cc:154
 hrungekutta.cc:155
 hrungekutta.cc:156
 hrungekutta.cc:157
 hrungekutta.cc:158
 hrungekutta.cc:159
 hrungekutta.cc:160
 hrungekutta.cc:161
 hrungekutta.cc:162
 hrungekutta.cc:163
 hrungekutta.cc:164
 hrungekutta.cc:165
 hrungekutta.cc:166
 hrungekutta.cc:167
 hrungekutta.cc:168
 hrungekutta.cc:169
 hrungekutta.cc:170
 hrungekutta.cc:171
 hrungekutta.cc:172
 hrungekutta.cc:173
 hrungekutta.cc:174
 hrungekutta.cc:175
 hrungekutta.cc:176
 hrungekutta.cc:177
 hrungekutta.cc:178
 hrungekutta.cc:179
 hrungekutta.cc:180
 hrungekutta.cc:181
 hrungekutta.cc:182
 hrungekutta.cc:183
 hrungekutta.cc:184
 hrungekutta.cc:185
 hrungekutta.cc:186
 hrungekutta.cc:187
 hrungekutta.cc:188
 hrungekutta.cc:189
 hrungekutta.cc:190
 hrungekutta.cc:191
 hrungekutta.cc:192
 hrungekutta.cc:193
 hrungekutta.cc:194
 hrungekutta.cc:195
 hrungekutta.cc:196
 hrungekutta.cc:197
 hrungekutta.cc:198
 hrungekutta.cc:199
 hrungekutta.cc:200
 hrungekutta.cc:201
 hrungekutta.cc:202
 hrungekutta.cc:203
 hrungekutta.cc:204
 hrungekutta.cc:205
 hrungekutta.cc:206
 hrungekutta.cc:207
 hrungekutta.cc:208
 hrungekutta.cc:209
 hrungekutta.cc:210
 hrungekutta.cc:211
 hrungekutta.cc:212
 hrungekutta.cc:213
 hrungekutta.cc:214
 hrungekutta.cc:215
 hrungekutta.cc:216
 hrungekutta.cc:217
 hrungekutta.cc:218
 hrungekutta.cc:219
 hrungekutta.cc:220
 hrungekutta.cc:221
 hrungekutta.cc:222
 hrungekutta.cc:223
 hrungekutta.cc:224
 hrungekutta.cc:225
 hrungekutta.cc:226
 hrungekutta.cc:227
 hrungekutta.cc:228
 hrungekutta.cc:229
 hrungekutta.cc:230
 hrungekutta.cc:231
 hrungekutta.cc:232
 hrungekutta.cc:233
 hrungekutta.cc:234
 hrungekutta.cc:235
 hrungekutta.cc:236
 hrungekutta.cc:237
 hrungekutta.cc:238
 hrungekutta.cc:239
 hrungekutta.cc:240
 hrungekutta.cc:241
 hrungekutta.cc:242
 hrungekutta.cc:243
 hrungekutta.cc:244
 hrungekutta.cc:245
 hrungekutta.cc:246
 hrungekutta.cc:247
 hrungekutta.cc:248
 hrungekutta.cc:249
 hrungekutta.cc:250
 hrungekutta.cc:251
 hrungekutta.cc:252
 hrungekutta.cc:253
 hrungekutta.cc:254
 hrungekutta.cc:255
 hrungekutta.cc:256
 hrungekutta.cc:257
 hrungekutta.cc:258
 hrungekutta.cc:259
 hrungekutta.cc:260
 hrungekutta.cc:261
 hrungekutta.cc:262
 hrungekutta.cc:263
 hrungekutta.cc:264
 hrungekutta.cc:265
 hrungekutta.cc:266
 hrungekutta.cc:267
 hrungekutta.cc:268
 hrungekutta.cc:269
 hrungekutta.cc:270
 hrungekutta.cc:271
 hrungekutta.cc:272
 hrungekutta.cc:273
 hrungekutta.cc:274
 hrungekutta.cc:275
 hrungekutta.cc:276
 hrungekutta.cc:277
 hrungekutta.cc:278
 hrungekutta.cc:279
 hrungekutta.cc:280
 hrungekutta.cc:281
 hrungekutta.cc:282
 hrungekutta.cc:283
 hrungekutta.cc:284
 hrungekutta.cc:285
 hrungekutta.cc:286
 hrungekutta.cc:287
 hrungekutta.cc:288
 hrungekutta.cc:289
 hrungekutta.cc:290
 hrungekutta.cc:291
 hrungekutta.cc:292
 hrungekutta.cc:293
 hrungekutta.cc:294
 hrungekutta.cc:295
 hrungekutta.cc:296
 hrungekutta.cc:297
 hrungekutta.cc:298
 hrungekutta.cc:299
 hrungekutta.cc:300
 hrungekutta.cc:301
 hrungekutta.cc:302
 hrungekutta.cc:303
 hrungekutta.cc:304
 hrungekutta.cc:305
 hrungekutta.cc:306
 hrungekutta.cc:307
 hrungekutta.cc:308
 hrungekutta.cc:309
 hrungekutta.cc:310
 hrungekutta.cc:311
 hrungekutta.cc:312
 hrungekutta.cc:313
 hrungekutta.cc:314
 hrungekutta.cc:315
 hrungekutta.cc:316
 hrungekutta.cc:317
 hrungekutta.cc:318
 hrungekutta.cc:319
 hrungekutta.cc:320
 hrungekutta.cc:321
 hrungekutta.cc:322
 hrungekutta.cc:323
 hrungekutta.cc:324
 hrungekutta.cc:325
 hrungekutta.cc:326
 hrungekutta.cc:327
 hrungekutta.cc:328
 hrungekutta.cc:329
 hrungekutta.cc:330
 hrungekutta.cc:331
 hrungekutta.cc:332
 hrungekutta.cc:333
 hrungekutta.cc:334
 hrungekutta.cc:335
 hrungekutta.cc:336
 hrungekutta.cc:337
 hrungekutta.cc:338
 hrungekutta.cc:339
 hrungekutta.cc:340
 hrungekutta.cc:341
 hrungekutta.cc:342
 hrungekutta.cc:343
 hrungekutta.cc:344
 hrungekutta.cc:345
 hrungekutta.cc:346
 hrungekutta.cc:347
 hrungekutta.cc:348
 hrungekutta.cc:349
 hrungekutta.cc:350
 hrungekutta.cc:351
 hrungekutta.cc:352
 hrungekutta.cc:353
 hrungekutta.cc:354
 hrungekutta.cc:355
 hrungekutta.cc:356
 hrungekutta.cc:357
 hrungekutta.cc:358
 hrungekutta.cc:359
 hrungekutta.cc:360
 hrungekutta.cc:361
 hrungekutta.cc:362
 hrungekutta.cc:363
 hrungekutta.cc:364
 hrungekutta.cc:365
 hrungekutta.cc:366
 hrungekutta.cc:367
 hrungekutta.cc:368
 hrungekutta.cc:369
 hrungekutta.cc:370
 hrungekutta.cc:371
 hrungekutta.cc:372
 hrungekutta.cc:373
 hrungekutta.cc:374
 hrungekutta.cc:375
 hrungekutta.cc:376
 hrungekutta.cc:377
 hrungekutta.cc:378
 hrungekutta.cc:379
 hrungekutta.cc:380
 hrungekutta.cc:381
 hrungekutta.cc:382
 hrungekutta.cc:383
 hrungekutta.cc:384
 hrungekutta.cc:385
 hrungekutta.cc:386
 hrungekutta.cc:387
 hrungekutta.cc:388
 hrungekutta.cc:389
 hrungekutta.cc:390
 hrungekutta.cc:391
 hrungekutta.cc:392
 hrungekutta.cc:393
 hrungekutta.cc:394
 hrungekutta.cc:395
 hrungekutta.cc:396
 hrungekutta.cc:397
 hrungekutta.cc:398
 hrungekutta.cc:399
 hrungekutta.cc:400
 hrungekutta.cc:401
 hrungekutta.cc:402
 hrungekutta.cc:403
 hrungekutta.cc:404
 hrungekutta.cc:405
 hrungekutta.cc:406
 hrungekutta.cc:407
 hrungekutta.cc:408
 hrungekutta.cc:409
 hrungekutta.cc:410
 hrungekutta.cc:411
 hrungekutta.cc:412
 hrungekutta.cc:413
 hrungekutta.cc:414
 hrungekutta.cc:415
 hrungekutta.cc:416
 hrungekutta.cc:417
 hrungekutta.cc:418
 hrungekutta.cc:419
 hrungekutta.cc:420
 hrungekutta.cc:421
 hrungekutta.cc:422
 hrungekutta.cc:423
 hrungekutta.cc:424
 hrungekutta.cc:425
 hrungekutta.cc:426
 hrungekutta.cc:427
 hrungekutta.cc:428
 hrungekutta.cc:429
 hrungekutta.cc:430
 hrungekutta.cc:431
 hrungekutta.cc:432
 hrungekutta.cc:433
 hrungekutta.cc:434
 hrungekutta.cc:435
 hrungekutta.cc:436
 hrungekutta.cc:437
 hrungekutta.cc:438
 hrungekutta.cc:439
 hrungekutta.cc:440
 hrungekutta.cc:441
 hrungekutta.cc:442
 hrungekutta.cc:443
 hrungekutta.cc:444
 hrungekutta.cc:445
 hrungekutta.cc:446
 hrungekutta.cc:447
 hrungekutta.cc:448
 hrungekutta.cc:449
 hrungekutta.cc:450
 hrungekutta.cc:451
 hrungekutta.cc:452
 hrungekutta.cc:453
 hrungekutta.cc:454
 hrungekutta.cc:455
 hrungekutta.cc:456
 hrungekutta.cc:457
 hrungekutta.cc:458
 hrungekutta.cc:459
 hrungekutta.cc:460
 hrungekutta.cc:461
 hrungekutta.cc:462
 hrungekutta.cc:463
 hrungekutta.cc:464
 hrungekutta.cc:465
 hrungekutta.cc:466
 hrungekutta.cc:467
 hrungekutta.cc:468
 hrungekutta.cc:469
 hrungekutta.cc:470
 hrungekutta.cc:471
 hrungekutta.cc:472
 hrungekutta.cc:473
 hrungekutta.cc:474
 hrungekutta.cc:475
 hrungekutta.cc:476
 hrungekutta.cc:477
 hrungekutta.cc:478
 hrungekutta.cc:479
 hrungekutta.cc:480
 hrungekutta.cc:481
 hrungekutta.cc:482
 hrungekutta.cc:483
 hrungekutta.cc:484
 hrungekutta.cc:485
 hrungekutta.cc:486
 hrungekutta.cc:487
 hrungekutta.cc:488
 hrungekutta.cc:489
 hrungekutta.cc:490
 hrungekutta.cc:491
 hrungekutta.cc:492
 hrungekutta.cc:493
 hrungekutta.cc:494
 hrungekutta.cc:495
 hrungekutta.cc:496
 hrungekutta.cc:497
 hrungekutta.cc:498
 hrungekutta.cc:499
 hrungekutta.cc:500
 hrungekutta.cc:501
 hrungekutta.cc:502
 hrungekutta.cc:503
 hrungekutta.cc:504
 hrungekutta.cc:505
 hrungekutta.cc:506
 hrungekutta.cc:507
 hrungekutta.cc:508
 hrungekutta.cc:509
 hrungekutta.cc:510
 hrungekutta.cc:511
 hrungekutta.cc:512
 hrungekutta.cc:513
 hrungekutta.cc:514
 hrungekutta.cc:515
 hrungekutta.cc:516
 hrungekutta.cc:517
 hrungekutta.cc:518
 hrungekutta.cc:519
 hrungekutta.cc:520
 hrungekutta.cc:521
 hrungekutta.cc:522
 hrungekutta.cc:523
 hrungekutta.cc:524
 hrungekutta.cc:525
 hrungekutta.cc:526
 hrungekutta.cc:527
 hrungekutta.cc:528
 hrungekutta.cc:529
 hrungekutta.cc:530
 hrungekutta.cc:531
 hrungekutta.cc:532
 hrungekutta.cc:533
 hrungekutta.cc:534
 hrungekutta.cc:535
 hrungekutta.cc:536
 hrungekutta.cc:537
 hrungekutta.cc:538
 hrungekutta.cc:539
 hrungekutta.cc:540
 hrungekutta.cc:541
 hrungekutta.cc:542
 hrungekutta.cc:543
 hrungekutta.cc:544
 hrungekutta.cc:545
 hrungekutta.cc:546
 hrungekutta.cc:547
 hrungekutta.cc:548
 hrungekutta.cc:549
 hrungekutta.cc:550
 hrungekutta.cc:551
 hrungekutta.cc:552
 hrungekutta.cc:553
 hrungekutta.cc:554
 hrungekutta.cc:555
 hrungekutta.cc:556
 hrungekutta.cc:557
 hrungekutta.cc:558
 hrungekutta.cc:559
 hrungekutta.cc:560
 hrungekutta.cc:561
 hrungekutta.cc:562
 hrungekutta.cc:563
 hrungekutta.cc:564
 hrungekutta.cc:565
 hrungekutta.cc:566
 hrungekutta.cc:567
 hrungekutta.cc:568
 hrungekutta.cc:569
 hrungekutta.cc:570
 hrungekutta.cc:571
 hrungekutta.cc:572
 hrungekutta.cc:573
 hrungekutta.cc:574
 hrungekutta.cc:575
 hrungekutta.cc:576
 hrungekutta.cc:577
 hrungekutta.cc:578
 hrungekutta.cc:579
 hrungekutta.cc:580
 hrungekutta.cc:581
 hrungekutta.cc:582
 hrungekutta.cc:583
 hrungekutta.cc:584
 hrungekutta.cc:585
 hrungekutta.cc:586
 hrungekutta.cc:587
 hrungekutta.cc:588
 hrungekutta.cc:589
 hrungekutta.cc:590
 hrungekutta.cc:591
 hrungekutta.cc:592
 hrungekutta.cc:593
 hrungekutta.cc:594
 hrungekutta.cc:595
 hrungekutta.cc:596
 hrungekutta.cc:597
 hrungekutta.cc:598
 hrungekutta.cc:599
 hrungekutta.cc:600
 hrungekutta.cc:601
 hrungekutta.cc:602
 hrungekutta.cc:603
 hrungekutta.cc:604
 hrungekutta.cc:605
 hrungekutta.cc:606
 hrungekutta.cc:607
 hrungekutta.cc:608
 hrungekutta.cc:609
 hrungekutta.cc:610
 hrungekutta.cc:611
 hrungekutta.cc:612
 hrungekutta.cc:613
 hrungekutta.cc:614
 hrungekutta.cc:615
 hrungekutta.cc:616
 hrungekutta.cc:617
 hrungekutta.cc:618
 hrungekutta.cc:619
 hrungekutta.cc:620
 hrungekutta.cc:621
 hrungekutta.cc:622
 hrungekutta.cc:623
 hrungekutta.cc:624
 hrungekutta.cc:625
 hrungekutta.cc:626
 hrungekutta.cc:627
 hrungekutta.cc:628
 hrungekutta.cc:629
 hrungekutta.cc:630
 hrungekutta.cc:631
 hrungekutta.cc:632
 hrungekutta.cc:633
 hrungekutta.cc:634
 hrungekutta.cc:635
 hrungekutta.cc:636
 hrungekutta.cc:637
 hrungekutta.cc:638
 hrungekutta.cc:639
 hrungekutta.cc:640
 hrungekutta.cc:641
 hrungekutta.cc:642
 hrungekutta.cc:643
 hrungekutta.cc:644
 hrungekutta.cc:645
 hrungekutta.cc:646
 hrungekutta.cc:647
 hrungekutta.cc:648
 hrungekutta.cc:649
 hrungekutta.cc:650
 hrungekutta.cc:651
 hrungekutta.cc:652
 hrungekutta.cc:653
 hrungekutta.cc:654
 hrungekutta.cc:655
 hrungekutta.cc:656
 hrungekutta.cc:657
 hrungekutta.cc:658
 hrungekutta.cc:659
 hrungekutta.cc:660
 hrungekutta.cc:661
 hrungekutta.cc:662
 hrungekutta.cc:663
 hrungekutta.cc:664
 hrungekutta.cc:665
 hrungekutta.cc:666
 hrungekutta.cc:667
 hrungekutta.cc:668
 hrungekutta.cc:669
 hrungekutta.cc:670
 hrungekutta.cc:671
 hrungekutta.cc:672
 hrungekutta.cc:673
 hrungekutta.cc:674
 hrungekutta.cc:675
 hrungekutta.cc:676
 hrungekutta.cc:677
 hrungekutta.cc:678
 hrungekutta.cc:679
 hrungekutta.cc:680
 hrungekutta.cc:681
 hrungekutta.cc:682
 hrungekutta.cc:683
 hrungekutta.cc:684
 hrungekutta.cc:685
 hrungekutta.cc:686
 hrungekutta.cc:687
 hrungekutta.cc:688
 hrungekutta.cc:689
 hrungekutta.cc:690
 hrungekutta.cc:691
 hrungekutta.cc:692
 hrungekutta.cc:693
 hrungekutta.cc:694
 hrungekutta.cc:695
 hrungekutta.cc:696
 hrungekutta.cc:697
 hrungekutta.cc:698
 hrungekutta.cc:699
 hrungekutta.cc:700
 hrungekutta.cc:701
 hrungekutta.cc:702
 hrungekutta.cc:703
 hrungekutta.cc:704
 hrungekutta.cc:705
 hrungekutta.cc:706
 hrungekutta.cc:707
 hrungekutta.cc:708
 hrungekutta.cc:709
 hrungekutta.cc:710
 hrungekutta.cc:711
 hrungekutta.cc:712
 hrungekutta.cc:713
 hrungekutta.cc:714
 hrungekutta.cc:715
 hrungekutta.cc:716
 hrungekutta.cc:717
 hrungekutta.cc:718
 hrungekutta.cc:719
 hrungekutta.cc:720
 hrungekutta.cc:721
 hrungekutta.cc:722
 hrungekutta.cc:723
 hrungekutta.cc:724
 hrungekutta.cc:725
 hrungekutta.cc:726
 hrungekutta.cc:727
 hrungekutta.cc:728
 hrungekutta.cc:729
 hrungekutta.cc:730
 hrungekutta.cc:731
 hrungekutta.cc:732
 hrungekutta.cc:733
 hrungekutta.cc:734
 hrungekutta.cc:735
 hrungekutta.cc:736
 hrungekutta.cc:737
 hrungekutta.cc:738
 hrungekutta.cc:739
 hrungekutta.cc:740
 hrungekutta.cc:741
 hrungekutta.cc:742
 hrungekutta.cc:743
 hrungekutta.cc:744
 hrungekutta.cc:745
 hrungekutta.cc:746
 hrungekutta.cc:747
 hrungekutta.cc:748
 hrungekutta.cc:749
 hrungekutta.cc:750