[u/mrichter/AliRoot.git] / EVE / Reve / MCHelixLine.hi

#ifndef REVE_MCHelixLine_H
#define REVE_MCHelixLine_H

#include <Reve/Track.h>
#include <cassert>
#include <vector>

namespace Reve {

struct MCVertex
{
  Float_t x,y,z,t;

  MCVertex() : x(0), y(0), z(0), t(0) {}
  MCVertex(Float_t _x, Float_t _y, Float_t _z, Float_t _t=0) :
    x(_x), y(_y), z(_z), t(_t) {}

  Float_t Mag()  const { return TMath::Sqrt(x*x+y*y+z*z);}
  Float_t Mag2() const { return x*x+y*y+z*z;}

  Float_t Perp()  const { return TMath::Sqrt(x*x+y*y);}
  Float_t Perp2() const { return x*x+y*y;}
  Float_t R()     const { return Perp(); }

  MCVertex operator + (const MCVertex & b)
  { return MCVertex(x + b.x, y + b.y, z + b.z, t + b.t); }

  MCVertex operator - (const MCVertex & b)
  { return MCVertex(x - b.x, y - b.y, z - b.z, t - b.t); }

  MCVertex operator * (Float_t a)
  { return MCVertex(a*x, a*y, a*z, a*t); }

  MCVertex& operator +=(const MCVertex & b)
  { x += b.x; y += b.y; z += b.z; t += b.t; return *this; }
};


struct MCStruct
{
  MCStruct(const MCStruct&);            // Not implemented
  MCStruct& operator=(const MCStruct&); // Not implemented
 
  TrackRnrStyle*         fRnrMod;  
  std::vector<MCVertex>* fPoints;
  MCVertex               fV;
  Float_t                fVelocity; // size of particle velocity
  
  MCStruct(TrackRnrStyle* rs, MCVertex* v0 , Float_t vel, std::vector<MCVertex>* tpv) :
    fRnrMod   (rs),
    fPoints   (tpv),
    fV        (*v0),
    fVelocity (vel)
  {
    fPoints->push_back(fV);
  }
  virtual ~MCStruct() {}
};

/**************************************************************************/
//  HELIX
/**************************************************************************/

struct MCHelix : public MCStruct
{
  // constant
  Float_t fA;       // contains charge and magnetic field data

  //parameters dependend pT and pZ size, set in init function
  Float_t fLam;        // momentum ratio pT/pZ
  Float_t fR;          // a/pT
  Float_t fPhiStep;    // step size in xy projection, dependent of RnrMode and momentum
  Float_t fTimeStep;   
 
  Int_t   fN;           // step number in helix;
  Int_t   NMax;         // max number of points in helix
  Float_t x_off, y_off; // offset for fitting daughters
  Float_t sin, cos;
  Bool_t  crosR;       

  MCHelix(TrackRnrStyle* rs, MCVertex* v0, Float_t vel,
          std::vector<MCVertex>* tpv, Float_t a) :
    MCStruct(rs, v0 , vel, tpv),
    fA   (a),
    fLam (0),   fR (0), fPhiStep (0), fTimeStep (0),
    fN   (0), NMax (0),
    x_off(0), y_off(0),
    sin  (0), cos  (0), crosR (0)
  {}

  void Init(Float_t pT, Float_t pZ)
  {
    fN    = 0;
    crosR = false;
    x_off = 0;
    y_off = 0;
    fLam  = pZ/pT;
    fR    = pT/fA;

    fPhiStep = fRnrMod->fMinAng * TMath::DegToRad();
    if (fRnrMod->fDelta < TMath::Abs(fR))
    {
      Float_t ang  = 2*TMath::ACos(1 - fRnrMod->fDelta/TMath::Abs(fR));
      if (ang < fPhiStep) fPhiStep = ang; 
    }
    if (fA < 0) fPhiStep = -fPhiStep;

    // printf("MCHelix::init (%f/%f) labda %f time step %e phi step %f \n", pT, pZ,fLam,  fTimeStep,fPhiStep);
    fTimeStep = TMath::Abs(fR*fPhiStep)*TMath::Sqrt(1+fLam*fLam)/fVelocity;
    fTimeStep *= 0.01; //cm->m

    sin = TMath::Sin(fPhiStep); 
    cos = TMath::Cos(fPhiStep);
  }

  void SetBounds()
  {
    // check steps for max orbits
    NMax = Int_t(fRnrMod->fMaxOrbs*TMath::TwoPi()/TMath::Abs(fPhiStep));
    // check steps for Z boundaries
    Float_t nz;
    if(fLam > 0) {
      nz = (fRnrMod->fMaxZ - fV.z)/(fLam*TMath::Abs(fR*fPhiStep));
    } else {
      nz = (-fRnrMod->fMaxZ - fV.z)/(fLam*TMath::Abs(fR*fPhiStep));
    }
    //  printf("steps in helix line %d nz   %f vz %f\n", NMax, nz, fV.z);
    if (nz < NMax) NMax = Int_t(nz + 1);
    
    // check steps if circles intersect
    if(TMath::Sqrt(fV.x*fV.x+fV.y*fV.y) < fRnrMod->fMaxR + TMath::Abs(fR))
    {
      crosR = true;
    }
    // printf("end steps in helix line %d \n", NMax);
  }


  void Step(Float_t &px, Float_t &py, Float_t &/*pz*/) 
  {
    fV.t += fTimeStep;
    fV.x += (px*sin - py*(1 - cos))/fA + x_off;
    fV.y += (py*sin + px*(1 - cos))/fA + y_off;
    fV.z += fLam*TMath::Abs(fR*fPhiStep);
    fPoints->push_back(fV);
    Float_t px_t = px*cos - py*sin;
    Float_t py_t = py*cos + px*sin;
    px = px_t;
    py = py_t;
    ++fN;
  }


  Bool_t LoopToVertex(Float_t &px, Float_t &py, Float_t &pz, 
                      Float_t  ex, Float_t  ey, Float_t  ez)
  {
    Float_t p0x = px, p0y = py;
    Float_t zs = fLam*TMath::Abs(fR*fPhiStep);
    Float_t maxrsq  = fRnrMod->fMaxR * fRnrMod->fMaxR;
    Float_t fnsteps = (ez - fV.z)/zs;
    Int_t   nsteps  = Int_t((ez - fV.z)/zs);
    Float_t sinf = TMath::Sin(fnsteps*fPhiStep);
    Float_t cosf = TMath::Cos(fnsteps*fPhiStep);

    { 
      if (nsteps > 0)
      {
        Float_t xf  = fV.x + (px*sinf - py*(1 - cosf))/fA;  
        Float_t yf =  fV.y + (py*sinf + px*(1 - cosf))/fA;
        x_off =  (ex - xf)/fnsteps;
        y_off =  (ey - yf)/fnsteps;
        Float_t xforw, yforw, zforw;
        for (Int_t l=0; l<nsteps; l++)
	{
          xforw  = fV.x + (px*sin - py*(1 - cos))/fA + x_off;
          yforw =  fV.y + (py*sin + px*(1 - cos))/fA + y_off;
          zforw =  fV.z + fLam*TMath::Abs(fR*fPhiStep);
          if (xforw*xforw+yforw*yforw > maxrsq   ||
              TMath::Abs(zforw) > fRnrMod->fMaxZ)
          {
            return false;
          }
          Step(px, py, pz);
	  if (fN > NMax)
	    break;
        }
	 
      }
      // set time to the end point
      fV.t += TMath::Sqrt((fV.x-ex)*(fV.x-ex)+(fV.y-ey)*(fV.y-ey) +(fV.z-ez)*(fV.z-ez))/fVelocity;
      fV.x = ex; fV.y = ey; fV.z = ez;      
      fPoints->push_back(fV);
    }
      
    { // fix momentum in the remaining part
      Float_t cosr = TMath::Cos((fnsteps-nsteps)*fPhiStep); 
      Float_t sinr = TMath::Sin((fnsteps-nsteps)*fPhiStep); 
      Float_t px_t = px*cosr - py*sinr;
      Float_t py_t = py*cosr + px*sinr;
      px = px_t;
      py = py_t;
    }
    { // calculate direction of faked px,py
      Float_t pxf = (p0x*cosf - p0y*sinf)/TMath::Abs(fA) + x_off/fPhiStep;
      Float_t pyf = (p0y*cosf + p0x*sinf)/TMath::Abs(fA) + y_off/fPhiStep;
      Float_t fac = TMath::Sqrt((p0x*p0x + p0y*p0y) / (pxf*pxf + pyf*pyf));
      px = fac*pxf;
      py = fac*pyf;
    }
    return true;
  }
    
  Bool_t LoopToBounds(Float_t &px, Float_t &py, Float_t &pz)
  {
    //    printf("MC helix  loop_to_bounds\n");
    SetBounds();
    if (NMax > 0)
    {
      // printf("NMax MC helix  loop_to_bounds\n");
      MCVertex forw;
      Float_t maxrsq = fRnrMod->fMaxR * fRnrMod->fMaxR;
      while (fN < NMax)
      {
        forw.x = fV.x + (px*sin - py*(1 - cos))/fA + x_off;
        forw.y = fV.y + (py*sin + px*(1 - cos))/fA + y_off;
        forw.z = fV.z + fLam*TMath::Abs(fR*fPhiStep);
        forw.t = fV.t + fTimeStep;
	
        if (crosR && forw.Perp2() > maxrsq)
        {
          Float_t t = (fRnrMod->fMaxR - fV.R()) / (forw.R() - fV.R());
          assert(t >= 0 && t <= 1);
          fPoints->push_back(fV + (forw-fV)*t);
          return false;
        }
        if (TMath::Abs(forw.z) > fRnrMod->fMaxZ)
        {
          Float_t t = (fRnrMod->fMaxZ - TMath::Abs(fV.z)) / TMath::Abs((forw.z - fV.z));
          assert(t >= 0 && t <= 1);
          fPoints->push_back(fV + (forw-fV)*t);
          return false;
        }

        Step(px, py, pz);
      }
      return true;
    }
    return false;
  }
};

/**************************************************************************/
//  LINE
/**************************************************************************/

struct MCLine : public MCStruct
{
  MCLine(TrackRnrStyle* rs, MCVertex* v0 ,Float_t vel, std::vector<MCVertex>* tpv):
    MCStruct(rs, v0 , vel, tpv)
  {}

  Bool_t InBounds(Float_t ex, Float_t ey, Float_t ez)
  {
    if(TMath::Abs(ez) > fRnrMod->fMaxZ ||
       ex*ex + ey*ey  > fRnrMod->fMaxR*fRnrMod->fMaxR)
      return false;
    else
      return true;
  }

  void GotoVertex(Float_t  x1, Float_t  y1, Float_t  z1)
  {
    fV.t += TMath::Sqrt((fV.x-x1)*(fV.x-x1)+(fV.y-y1)*(fV.y-y1)+(fV.z-z1)*(fV.z-z1))/fVelocity;
    fV.x=x1; fV.y=y1; fV.z=z1;
    fPoints->push_back(fV);
  }
  

  void GotoBounds( Float_t px, Float_t py, Float_t pz)
  {
    Float_t tZ = 0,Tb = 0;
    // time where particle intersect +/- fMaxZ
    if (pz > 0) {
      tZ = (fRnrMod->fMaxZ - fV.z)/pz;
    }
    else  if (pz < 0 ) {
      tZ = (-1)*(fRnrMod->fMaxZ + fV.z)/pz;
    }
    // time where particle intersects cylinder
    Float_t tR=0;
    Double_t a = px*px + py*py;
    Double_t b = 2*(fV.x*px + fV.y*py);
    Double_t c = fV.x*fV.x + fV.y*fV.y - fRnrMod->fMaxR*fRnrMod->fMaxR;
    Double_t D = b*b - 4*a*c;
    if(D >= 0) {
      Double_t D_sqrt=TMath::Sqrt(D);
      tR = ( -b - D_sqrt )/(2*a);
      if( tR < 0) {
        tR = ( -b + D_sqrt )/(2*a);
      }

      // compare the two times
      Tb = tR < tZ ? tR : tZ;
    } else {
      Tb = tZ;
    }

    GotoVertex(fV.x+px*Tb, fV.y+py*Tb, fV.z+ pz*Tb);
  }
}; // struct Line


} // namespace Reve

#endif
Commit	Line	Data
5a5a1232	1	#ifndef REVE_MCHelixLine_H
	2	#define REVE_MCHelixLine_H
	3
	4	#include <Reve/Track.h>
32e219c2	5	#include <cassert>
5a5a1232	6	#include <vector>
	7
	8	namespace Reve {
	9
bbfaa1c4	10	struct MCVertex
	11	{
	12	Float_t x,y,z,t;
32e219c2	13
	14	MCVertex() : x(0), y(0), z(0), t(0) {}
	15	MCVertex(Float_t _x, Float_t _y, Float_t _z, Float_t _t=0) :
	16	x(_x), y(_y), z(_z), t(_t) {}
	17
	18	Float_t Mag() const { return TMath::Sqrt(xx+yy+z*z);}
	19	Float_t Mag2() const { return xx+yy+z*z;}
	20
	21	Float_t Perp() const { return TMath::Sqrt(xx+yy);}
	22	Float_t Perp2() const { return xx+yy;}
	23	Float_t R() const { return Perp(); }
	24
	25	MCVertex operator + (const MCVertex & b)
	26	{ return MCVertex(x + b.x, y + b.y, z + b.z, t + b.t); }
	27
	28	MCVertex operator - (const MCVertex & b)
	29	{ return MCVertex(x - b.x, y - b.y, z - b.z, t - b.t); }
	30
	31	MCVertex operator * (Float_t a)
	32	{ return MCVertex(ax, ay, az, at); }
	33
	34	MCVertex& operator +=(const MCVertex & b)
	35	{ x += b.x; y += b.y; z += b.z; t += b.t; return *this; }
bbfaa1c4	36	};
5a5a1232	37
5a5a1232	38
bbfaa1c4	39	struct MCStruct
	40	{
	41	MCStruct(const MCStruct&); // Not implemented
	42	MCStruct& operator=(const MCStruct&); // Not implemented
	43
	44	TrackRnrStyle* fRnrMod;
	45	std::vector<MCVertex>* fPoints;
	46	MCVertex fV;
	47	Float_t fVelocity; // size of particle velocity
5a5a1232	48
bbfaa1c4	49	MCStruct(TrackRnrStyle* rs, MCVertex* v0 , Float_t vel, std::vector<MCVertex>* tpv) :
	50	fRnrMod (rs),
	51	fPoints (tpv),
	52	fV (*v0),
	53	fVelocity (vel)
b97c26da	54	{
bbfaa1c4	55	fPoints->push_back(fV);
	56	}
	57	virtual ~MCStruct() {}
	58	};
5a5a1232	59
bbfaa1c4	60	/**************************************************************************/
	61	// HELIX
	62	/**************************************************************************/
5a5a1232	63
bbfaa1c4	64	struct MCHelix : public MCStruct
	65	{
	66	// constant
	67	Float_t fA; // contains charge and magnetic field data
5a5a1232	68
bbfaa1c4	69	//parameters dependend pT and pZ size, set in init function
	70	Float_t fLam; // momentum ratio pT/pZ
	71	Float_t fR; // a/pT
	72	Float_t fPhiStep; // step size in xy projection, dependent of RnrMode and momentum
	73	Float_t fTimeStep;
	74
	75	Int_t fN; // step number in helix;
	76	Int_t NMax; // max number of points in helix
	77	Float_t x_off, y_off; // offset for fitting daughters
	78	Float_t sin, cos;
	79	Bool_t crosR;
5a5a1232	80
bbfaa1c4	81	MCHelix(TrackRnrStyle* rs, MCVertex* v0, Float_t vel,
	82	std::vector<MCVertex>* tpv, Float_t a) :
	83	MCStruct(rs, v0 , vel, tpv),
	84	fA (a),
	85	fLam (0), fR (0), fPhiStep (0), fTimeStep (0),
	86	fN (0), NMax (0),
	87	x_off(0), y_off(0),
	88	sin (0), cos (0), crosR (0)
	89	{}
5a5a1232	90
bbfaa1c4	91	void Init(Float_t pT, Float_t pZ)
bbfaa1c4	92	{
32e219c2	93	fN = 0;
bbfaa1c4	94	crosR = false;
32e219c2	95	x_off = 0;
	96	y_off = 0;
	97	fLam = pZ/pT;
	98	fR = pT/fA;
5a5a1232	99
32e219c2	100	fPhiStep = fRnrMod->fMinAng * TMath::DegToRad();
	101	if (fRnrMod->fDelta < TMath::Abs(fR))
	102	{
bbfaa1c4	103	Float_t ang = 2*TMath::ACos(1 - fRnrMod->fDelta/TMath::Abs(fR));
bbfaa1c4	104	if (ang < fPhiStep) fPhiStep = ang;
5a5a1232	105	}
32e219c2	106	if (fA < 0) fPhiStep = -fPhiStep;
5a5a1232	107
bbfaa1c4	108	// printf("MCHelix::init (%f/%f) labda %f time step %e phi step %f \n", pT, pZ,fLam, fTimeStep,fPhiStep);
	109	fTimeStep = TMath::Abs(fRfPhiStep)TMath::Sqrt(1+fLam*fLam)/fVelocity;
	110	fTimeStep *= 0.01; //cm->m
	111
	112	sin = TMath::Sin(fPhiStep);
	113	cos = TMath::Cos(fPhiStep);
	114	}
	115
	116	void SetBounds()
	117	{
	118	// check steps for max orbits
	119	NMax = Int_t(fRnrMod->fMaxOrbs*TMath::TwoPi()/TMath::Abs(fPhiStep));
	120	// check steps for Z boundaries
	121	Float_t nz;
	122	if(fLam > 0) {
	123	nz = (fRnrMod->fMaxZ - fV.z)/(fLamTMath::Abs(fRfPhiStep));
	124	} else {
	125	nz = (-fRnrMod->fMaxZ - fV.z)/(fLamTMath::Abs(fRfPhiStep));
	126	}
	127	// printf("steps in helix line %d nz %f vz %f\n", NMax, nz, fV.z);
32e219c2	128	if (nz < NMax) NMax = Int_t(nz + 1);
5a5a1232	129
bbfaa1c4	130	// check steps if circles intersect
32e219c2	131	if(TMath::Sqrt(fV.xfV.x+fV.yfV.y) < fRnrMod->fMaxR + TMath::Abs(fR))
32e219c2	132	{
bbfaa1c4	133	crosR = true;
5a5a1232	134	}
bbfaa1c4	135	// printf("end steps in helix line %d \n", NMax);
bbfaa1c4	136	}
5a5a1232	137
5a5a1232	138
bbfaa1c4	139	void Step(Float_t &px, Float_t &py, Float_t &/pz/)
	140	{
	141	fV.t += fTimeStep;
	142	fV.x += (pxsin - py(1 - cos))/fA + x_off;
	143	fV.y += (pysin + px(1 - cos))/fA + y_off;
	144	fV.z += fLamTMath::Abs(fRfPhiStep);
	145	fPoints->push_back(fV);
32e219c2	146	Float_t px_t = pxcos - pysin;
32e219c2	147	Float_t py_t = pycos + pxsin;
bbfaa1c4	148	px = px_t;
bbfaa1c4	149	py = py_t;
b855ed81	150	++fN;
bbfaa1c4	151	}
5a5a1232	152
5a5a1232	153
bbfaa1c4	154	Bool_t LoopToVertex(Float_t &px, Float_t &py, Float_t &pz,
f3e27855	155	Float_t ex, Float_t ey, Float_t ez)
bbfaa1c4	156	{
	157	Float_t p0x = px, p0y = py;
	158	Float_t zs = fLamTMath::Abs(fRfPhiStep);
32e219c2	159	Float_t maxrsq = fRnrMod->fMaxR * fRnrMod->fMaxR;
bbfaa1c4	160	Float_t fnsteps = (ez - fV.z)/zs;
	161	Int_t nsteps = Int_t((ez - fV.z)/zs);
	162	Float_t sinf = TMath::Sin(fnsteps*fPhiStep);
	163	Float_t cosf = TMath::Cos(fnsteps*fPhiStep);
5a5a1232	164
bbfaa1c4	165	{
b855ed81	166	if (nsteps > 0)
b855ed81	167	{
bbfaa1c4	168	Float_t xf = fV.x + (pxsinf - py(1 - cosf))/fA;
	169	Float_t yf = fV.y + (pysinf + px(1 - cosf))/fA;
	170	x_off = (ex - xf)/fnsteps;
	171	y_off = (ey - yf)/fnsteps;
	172	Float_t xforw, yforw, zforw;
b855ed81	173	for (Int_t l=0; l<nsteps; l++)
b855ed81	174	{
bbfaa1c4	175	xforw = fV.x + (pxsin - py(1 - cos))/fA + x_off;
	176	yforw = fV.y + (pysin + px(1 - cos))/fA + y_off;
	177	zforw = fV.z + fLamTMath::Abs(fRfPhiStep);
32e219c2	178	if (xforwxforw+yforwyforw > maxrsq \|\|
	179	TMath::Abs(zforw) > fRnrMod->fMaxZ)
	180	{
bbfaa1c4	181	return false;
bbfaa1c4	182	}
32e219c2	183	Step(px, py, pz);
b855ed81	184	if (fN > NMax)
b855ed81	185	break;
bbfaa1c4	186	}
5a5a1232	187
5a5a1232	188	}
bbfaa1c4	189	// set time to the end point
	190	fV.t += TMath::Sqrt((fV.x-ex)(fV.x-ex)+(fV.y-ey)(fV.y-ey) +(fV.z-ez)*(fV.z-ez))/fVelocity;
	191	fV.x = ex; fV.y = ey; fV.z = ez;
	192	fPoints->push_back(fV);
	193	}
5a5a1232	194
bbfaa1c4	195	{ // fix momentum in the remaining part
32e219c2	196	Float_t cosr = TMath::Cos((fnsteps-nsteps)*fPhiStep);
	197	Float_t sinr = TMath::Sin((fnsteps-nsteps)*fPhiStep);
	198	Float_t px_t = pxcosr - pysinr;
	199	Float_t py_t = pycosr + pxsinr;
bbfaa1c4	200	px = px_t;
	201	py = py_t;
	202	}
	203	{ // calculate direction of faked px,py
	204	Float_t pxf = (p0xcosf - p0ysinf)/TMath::Abs(fA) + x_off/fPhiStep;
	205	Float_t pyf = (p0ycosf + p0xsinf)/TMath::Abs(fA) + y_off/fPhiStep;
32e219c2	206	Float_t fac = TMath::Sqrt((p0xp0x + p0yp0y) / (pxfpxf + pyfpyf));
bbfaa1c4	207	px = fac*pxf;
bbfaa1c4	208	py = fac*pyf;
5a5a1232	209	}
bbfaa1c4	210	return true;
bbfaa1c4	211	}
5a5a1232	212
bbfaa1c4	213	Bool_t LoopToBounds(Float_t &px, Float_t &py, Float_t &pz)
	214	{
	215	// printf("MC helix loop_to_bounds\n");
	216	SetBounds();
b855ed81	217	if (NMax > 0)
b855ed81	218	{
32e219c2	219	// printf("NMax MC helix loop_to_bounds\n");
	220	MCVertex forw;
	221	Float_t maxrsq = fRnrMod->fMaxR * fRnrMod->fMaxR;
b855ed81	222	while (fN < NMax)
b855ed81	223	{
32e219c2	224	forw.x = fV.x + (pxsin - py(1 - cos))/fA + x_off;
	225	forw.y = fV.y + (pysin + px(1 - cos))/fA + y_off;
	226	forw.z = fV.z + fLamTMath::Abs(fRfPhiStep);
	227	forw.t = fV.t + fTimeStep;
5a5a1232	228
32e219c2	229	if (crosR && forw.Perp2() > maxrsq)
	230	{
	231	Float_t t = (fRnrMod->fMaxR - fV.R()) / (forw.R() - fV.R());
	232	assert(t >= 0 && t <= 1);
	233	fPoints->push_back(fV + (forw-fV)*t);
bbfaa1c4	234	return false;
bbfaa1c4	235	}
32e219c2	236	if (TMath::Abs(forw.z) > fRnrMod->fMaxZ)
	237	{
	238	Float_t t = (fRnrMod->fMaxZ - TMath::Abs(fV.z)) / TMath::Abs((forw.z - fV.z));
	239	assert(t >= 0 && t <= 1);
	240	fPoints->push_back(fV + (forw-fV)*t);
	241	return false;
	242	}
	243
	244	Step(px, py, pz);
5a5a1232	245	}
bbfaa1c4	246	return true;
5a5a1232	247	}
bbfaa1c4	248	return false;
	249	}
	250	};
5a5a1232	251
bbfaa1c4	252	/**************************************************************************/
	253	// LINE
	254	/**************************************************************************/
5a5a1232	255
bbfaa1c4	256	struct MCLine : public MCStruct
	257	{
	258	MCLine(TrackRnrStyle* rs, MCVertex* v0 ,Float_t vel, std::vector<MCVertex>* tpv):
	259	MCStruct(rs, v0 , vel, tpv)
	260	{}
5a5a1232	261
bbfaa1c4	262	Bool_t InBounds(Float_t ex, Float_t ey, Float_t ez)
	263	{
	264	if(TMath::Abs(ez) > fRnrMod->fMaxZ \|\|
	265	exex + eyey > fRnrMod->fMaxR*fRnrMod->fMaxR)
	266	return false;
	267	else
	268	return true;
	269	}
5a5a1232	270
bbfaa1c4	271	void GotoVertex(Float_t x1, Float_t y1, Float_t z1)
bbfaa1c4	272	{
bbfaa1c4	273	fV.t += TMath::Sqrt((fV.x-x1)(fV.x-x1)+(fV.y-y1)(fV.y-y1)+(fV.z-z1)*(fV.z-z1))/fVelocity;
	274	fV.x=x1; fV.y=y1; fV.z=z1;
	275	fPoints->push_back(fV);
	276	}
5a5a1232	277
5a5a1232	278
bbfaa1c4	279	void GotoBounds( Float_t px, Float_t py, Float_t pz)
	280	{
	281	Float_t tZ = 0,Tb = 0;
	282	// time where particle intersect +/- fMaxZ
	283	if (pz > 0) {
	284	tZ = (fRnrMod->fMaxZ - fV.z)/pz;
	285	}
	286	else if (pz < 0 ) {
	287	tZ = (-1)*(fRnrMod->fMaxZ + fV.z)/pz;
	288	}
	289	// time where particle intersects cylinder
	290	Float_t tR=0;
	291	Double_t a = pxpx + pypy;
	292	Double_t b = 2(fV.xpx + fV.y*py);
	293	Double_t c = fV.xfV.x + fV.yfV.y - fRnrMod->fMaxR*fRnrMod->fMaxR;
	294	Double_t D = bb - 4a*c;
	295	if(D >= 0) {
	296	Double_t D_sqrt=TMath::Sqrt(D);
	297	tR = ( -b - D_sqrt )/(2*a);
	298	if( tR < 0) {
	299	tR = ( -b + D_sqrt )/(2*a);
5a5a1232	300	}
5a5a1232	301
bbfaa1c4	302	// compare the two times
	303	Tb = tR < tZ ? tR : tZ;
	304	} else {
	305	Tb = tZ;
5a5a1232	306	}
bbfaa1c4	307
	308	GotoVertex(fV.x+pxTb, fV.y+pyTb, fV.z+ pz*Tb);
	309	}
	310	}; // struct Line
5a5a1232	311
	312
	313	} // namespace Reve
	314
	315	#endif