EVE/Reve/MCHelixLine.hi

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