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

#ifndef REVE_MCHelixLine_H
#define REVE_MCHelixLine_H

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

namespace Reve {

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


  struct MCStruct
  {
    TrackRnrStyle*      fRnrMod;  
    std::vector<MCVertex>* track_points;
    MCVertex               v;
    Float_t                fVelocity; // size of particle velocity
  
    MCStruct(TrackRnrStyle* rs, MCVertex* v0 , Float_t vel, std::vector<MCVertex>* tpv) 
    { 
      fRnrMod = rs;
      v = *v0;
      fVelocity = vel;
      track_points = tpv;
      track_points->push_back(v);
    }
  };

  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;
    };

    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::Pi()/180;
      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 - v.z)/(fLam*TMath::Abs(fR*fPhiStep));
      } else {
	nz = (-fRnrMod->fMaxZ - v.z)/(fLam*TMath::Abs(fR*fPhiStep));
      }
      //  printf("steps in helix line %d nz   %f vz %f\n", NMax, nz, v.z);
      if (nz < NMax) NMax = Int_t(nz);
    
      // check steps if circles intersect
      if(TMath::Sqrt(v.x*v.x+v.y*v.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*/) 
    {
      v.t += fTimeStep;
      v.x += (px*sin - py*(1 - cos))/fA + x_off;
      v.y += (py*sin + px*(1 - cos))/fA + y_off;
      v.z += fLam*TMath::Abs(fR*fPhiStep);
      track_points->push_back(v);
      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 fnsteps = (ez - v.z)/zs;
      Int_t   nsteps  = Int_t((ez - v.z)/zs);
      Float_t sinf = TMath::Sin(fnsteps*fPhiStep);
      Float_t cosf = TMath::Cos(fnsteps*fPhiStep);

      { 
	track_points->push_back(v); 
	if(nsteps > 0){
	  Float_t xf  = v.x + (px*sinf - py*(1 - cosf))/fA;  
	  Float_t yf =  v.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  = v.x + (px*sin - py*(1 - cos))/fA + x_off;
	    yforw =  v.y + (py*sin + px*(1 - cos))/fA + y_off;
	    zforw =  v.z + fLam*TMath::Abs(fR*fPhiStep);
	    if ((xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ) ) {
	      return false;
	    }
	    Step(px,py,pz); 
	  }
	 
	}
	// set time to the end point
	v.t += TMath::Sqrt((v.x-ex)*(v.x-ex)+(v.y-ey)*(v.y-ey) +(v.z-ez)*(v.z-ez))/fVelocity;
	v.x = ex; v.y = ey; v.z = ez;      
	track_points->push_back(v);
      }
      
      { // 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)/TMath::Sqrt(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");
	track_points->push_back(v);
	Float_t xforw,yforw,zforw;
	while(fN < NMax){
	  xforw = v.x + (px*sin - py*(1 - cos))/fA + x_off;
	  yforw = v.y + (py*sin + px*(1 - cos))/fA + y_off;
	  zforw =  v.z + fLam*TMath::Abs(fR*fPhiStep);
	
	  if ((crosR && (xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR)) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ)) {
	    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)
    {
      track_points->push_back(v);
      v.t += TMath::Sqrt((v.x-x1)*(v.x-x1)+(v.y-y1)*(v.y-y1)+(v.z-z1)*(v.z-z1))/fVelocity;
      v.x=x1; v.y=y1; v.z=z1;
      track_points->push_back(v);
    }
  

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