EVE/Reve/MCHelixLine.hi

   1 #ifndef REVE_MCHelixLine_H
   2 #define REVE_MCHelixLine_H
   3
   4 #include <Reve/Track.h>
   5 #include <cassert>
   6 #include <vector>
   7
   8 namespace Reve {
   9
  10 struct MCVertex
  11 {
  12   Float_t x,y,z,t;
  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(x*x+y*y+z*z);}
  19   Float_t Mag2() const { return x*x+y*y+z*z;}
  20
  21   Float_t Perp()  const { return TMath::Sqrt(x*x+y*y);}
  22   Float_t Perp2() const { return x*x+y*y;}
  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(a*x, a*y, a*z, a*t); }
  33
  34   MCVertex& operator +=(const MCVertex & b)
  35   { x += b.x; y += b.y; z += b.z; t += b.t; return *this; }
  36 };
  37
  38
  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
  48
  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)
  54   {
  55     fPoints->push_back(fV);
  56   }
  57   virtual ~MCStruct() {}
  58 };
  59
  60 /**************************************************************************/
  61 //  HELIX
  62 /**************************************************************************/
  63
  64 struct MCHelix : public MCStruct
  65 {
  66   // constant
  67   Float_t fA;       // contains charge and magnetic field data
  68
  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;
  80
  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   {}
  90
  91   void Init(Float_t pT, Float_t pZ)
  92   {
  93     fN    = 0;
  94     crosR = false;
  95     x_off = 0;
  96     y_off = 0;
  97     fLam  = pZ/pT;
  98     fR    = pT/fA;
  99
 100     fPhiStep = fRnrMod->fMinAng * TMath::DegToRad();
 101     if (fRnrMod->fDelta < TMath::Abs(fR))
 102     {
 103       Float_t ang  = 2*TMath::ACos(1 - fRnrMod->fDelta/TMath::Abs(fR));
 104       if (ang < fPhiStep) fPhiStep = ang;
 105     }
 106     if (fA < 0) fPhiStep = -fPhiStep;
 107
 108     // printf("MCHelix::init (%f/%f) labda %f time step %e phi step %f \n", pT, pZ,fLam,  fTimeStep,fPhiStep);
 109     fTimeStep = TMath::Abs(fR*fPhiStep)*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)/(fLam*TMath::Abs(fR*fPhiStep));
 124     } else {
 125       nz = (-fRnrMod->fMaxZ - fV.z)/(fLam*TMath::Abs(fR*fPhiStep));
 126     }
 127     //  printf("steps in helix line %d nz   %f vz %f\n", NMax, nz, fV.z);
 128     if (nz < NMax) NMax = Int_t(nz + 1);
 129
 130     // check steps if circles intersect
 131     if(TMath::Sqrt(fV.x*fV.x+fV.y*fV.y) < fRnrMod->fMaxR + TMath::Abs(fR))
 132     {
 133       crosR = true;
 134     }
 135     // printf("end steps in helix line %d \n", NMax);
 136   }
 137
 138
 139   void Step(Float_t &px, Float_t &py, Float_t &/*pz*/)
 140   {
 141     fV.t += fTimeStep;
 142     fV.x += (px*sin - py*(1 - cos))/fA + x_off;
 143     fV.y += (py*sin + px*(1 - cos))/fA + y_off;
 144     fV.z += fLam*TMath::Abs(fR*fPhiStep);
 145     fPoints->push_back(fV);
 146     Float_t px_t = px*cos - py*sin;
 147     Float_t py_t = py*cos + px*sin;
 148     px = px_t;
 149     py = py_t;
 150     ++fN;
 151   }
 152
 153
 154   Bool_t LoopToVertex(Float_t &px, Float_t &py, Float_t &pz,
 155                       Float_t  ex, Float_t  ey, Float_t  ez)
 156   {
 157     Float_t p0x = px, p0y = py;
 158     Float_t zs = fLam*TMath::Abs(fR*fPhiStep);
 159     Float_t maxrsq  = fRnrMod->fMaxR * fRnrMod->fMaxR;
 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);
 164
 165     {
 166       if (nsteps > 0)
 167       {
 168         Float_t xf  = fV.x + (px*sinf - py*(1 - cosf))/fA;
 169         Float_t yf =  fV.y + (py*sinf + px*(1 - cosf))/fA;
 170         x_off =  (ex - xf)/fnsteps;
 171         y_off =  (ey - yf)/fnsteps;
 172         Float_t xforw, yforw, zforw;
 173         for (Int_t l=0; l<nsteps; l++)
 174         {
 175           xforw  = fV.x + (px*sin - py*(1 - cos))/fA + x_off;
 176           yforw =  fV.y + (py*sin + px*(1 - cos))/fA + y_off;
 177           zforw =  fV.z + fLam*TMath::Abs(fR*fPhiStep);
 178           if (xforw*xforw+yforw*yforw > maxrsq   ||
 179               TMath::Abs(zforw) > fRnrMod->fMaxZ)
 180           {
 181             return false;
 182           }
 183           Step(px, py, pz);
 184           if (fN > NMax)
 185             break;
 186         }
 187
 188       }
 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     }
 194
 195     { // fix momentum in the remaining part
 196       Float_t cosr = TMath::Cos((fnsteps-nsteps)*fPhiStep);
 197       Float_t sinr = TMath::Sin((fnsteps-nsteps)*fPhiStep);
 198       Float_t px_t = px*cosr - py*sinr;
 199       Float_t py_t = py*cosr + px*sinr;
 200       px = px_t;
 201       py = py_t;
 202     }
 203     { // calculate direction of faked px,py
 204       Float_t pxf = (p0x*cosf - p0y*sinf)/TMath::Abs(fA) + x_off/fPhiStep;
 205       Float_t pyf = (p0y*cosf + p0x*sinf)/TMath::Abs(fA) + y_off/fPhiStep;
 206       Float_t fac = TMath::Sqrt((p0x*p0x + p0y*p0y) / (pxf*pxf + pyf*pyf));
 207       px = fac*pxf;
 208       py = fac*pyf;
 209     }
 210     return true;
 211   }
 212
 213   Bool_t LoopToBounds(Float_t &px, Float_t &py, Float_t &pz)
 214   {
 215     //    printf("MC helix  loop_to_bounds\n");
 216     SetBounds();
 217     if (NMax > 0)
 218     {
 219       // printf("NMax MC helix  loop_to_bounds\n");
 220       MCVertex forw;
 221       Float_t maxrsq = fRnrMod->fMaxR * fRnrMod->fMaxR;
 222       while (fN < NMax)
 223       {
 224         forw.x = fV.x + (px*sin - py*(1 - cos))/fA + x_off;
 225         forw.y = fV.y + (py*sin + px*(1 - cos))/fA + y_off;
 226         forw.z = fV.z + fLam*TMath::Abs(fR*fPhiStep);
 227         forw.t = fV.t + fTimeStep;
 228
 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);
 234           return false;
 235         }
 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);
 245       }
 246       return true;
 247     }
 248     return false;
 249   }
 250 };
 251
 252 /**************************************************************************/
 253 //  LINE
 254 /**************************************************************************/
 255
 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   {}
 261
 262   Bool_t InBounds(Float_t ex, Float_t ey, Float_t ez)
 263   {
 264     if(TMath::Abs(ez) > fRnrMod->fMaxZ ||
 265        ex*ex + ey*ey  > fRnrMod->fMaxR*fRnrMod->fMaxR)
 266       return false;
 267     else
 268       return true;
 269   }
 270
 271   void GotoVertex(Float_t  x1, Float_t  y1, Float_t  z1)
 272   {
 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   }
 277
 278
 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 = px*px + py*py;
 292     Double_t b = 2*(fV.x*px + fV.y*py);
 293     Double_t c = fV.x*fV.x + fV.y*fV.y - fRnrMod->fMaxR*fRnrMod->fMaxR;
 294     Double_t D = b*b - 4*a*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);
 300       }
 301
 302       // compare the two times
 303       Tb = tR < tZ ? tR : tZ;
 304     } else {
 305       Tb = tZ;
 306     }
 307
 308     GotoVertex(fV.x+px*Tb, fV.y+py*Tb, fV.z+ pz*Tb);
 309   }
 310 }; // struct Line
 311
 312
 313 } // namespace Reve
 314
 315 #endif