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