Fork me on GitHub

Changeset ae93700 in git


Ignore:
Timestamp:
May 27, 2021, 8:36:27 AM (3 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
master
Children:
716e186, dfdad8b
Parents:
c6f9311
Message:

improve ParticlePropagator code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/ParticlePropagator.cc

    rc6f9311 rae93700  
    126126  Double_t px, py, pz, pt, pt2, e, q;
    127127  Double_t x, y, z, t, r;
    128   Double_t x_c, y_c, r_c, phi_c, phi_0;
    129   Double_t x_t, y_t, z_t, r_t;
    130   Double_t t_z, t_r;
    131   Double_t discr;
     128  Double_t x_c, y_c, r_c, phi_0;
     129  Double_t x_t, y_t, z_t, r_t, phi_t;
     130  Double_t t_r, t_z;
     131  Double_t tmp;
    132132  Double_t gammam, omega;
    133133  Double_t xd, yd, zd;
    134134  Double_t l, d0, dz, ctgTheta, alpha;
    135135  Double_t bsx, bsy, bsz;
    136   Double_t rxp, rdp, t_R;
    137   Double_t td, pio, phid, sign_pz, vz;
     136  Double_t td, pio, phid, vz;
    138137
    139138  const Double_t c_light = 2.99792458E8;
    140139
    141140  if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
     141  {
    142142    beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
     143  }
    143144  else
    144145  {
     
    161162    particlePosition = particle->Position;
    162163    particleMomentum = particle->Momentum;
    163 
    164     // Constants
    165164
    166165    x = particlePosition.X() * 1.0E-3;
     
    208207    else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
    209208    {
    210 
    211       rxp = x*py - y*px;
    212       rdp = x*px + y*py;
    213 
    214       discr = fRadius*fRadius*pt*pt - rxp*rxp;
    215 
    216       t_R = e * (sqrt(discr) - rdp) / (c_light * pt * pt);
    217       t_z = e * (TMath::Sign(fHalfLengthMax, pz) - z) / ( c_light * pz);
    218 
    219       t = TMath::Min(t_R, t_z);
    220 
    221       x_t = x + px*t*c_light/e;
    222       y_t = y + py*t*c_light/e;
    223       z_t = z + pz*t*c_light/e;
    224       r_t = TMath::Hypot(x_t, y_t);
    225 
    226       l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
     209      // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
     210      tmp = px * y - py * x;
     211      t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
     212
     213      t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
     214
     215      t = TMath::Min(t_r, t_z);
     216
     217      x_t = x + px * t;
     218      y_t = y + py * t;
     219      z_t = z + pz * t;
     220
     221      l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
    227222
    228223      mother = candidate;
    229       candidate = static_cast<Candidate*>(candidate->Clone());
     224      candidate = static_cast<Candidate *>(candidate->Clone());
    230225
    231226      candidate->InitialPosition = particlePosition;
    232       candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
    233       candidate->L = l*1.0E3;
     227      candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
     228      candidate->L = l * 1.0E3;
    234229
    235230      candidate->Momentum = particleMomentum;
     
    260255    {
    261256
    262       // 1.  initial transverse momentum p_{T0}: Part->pt
    263       //     initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
    264       //     relativistic gamma: gamma = E/mc^2; gammam = gamma * m
    265       //     gyration frequency omega = q/(gamma m) fBz
    266       //     helix radius r = p_{T0} / (omega gamma m)
    267 
    268       gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c^2]
    269       omega = q * fBz / (gammam);                // omega is here in [89875518/s]
    270       r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
     257      // 1. initial transverse momentum p_{T0}: Part->pt
     258      //    initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
     259      //    relativistic gamma: gamma = E / mc^2; gammam = gamma * m
     260      //    gyration frequency omega = q * Bz / (gammam)
     261      //    helix radius r = p_{T0} / (omega * gammam)
     262
     263      gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
     264      omega = q * fBz / gammam; // omega is here in [89875518/s]
     265      r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
    271266
    272267      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    273268
    274269      // 2. helix axis coordinates
    275       x_c = x + r*TMath::Sin(phi_0);
    276       y_c = y - r*TMath::Cos(phi_0);
     270      x_c = x + r * TMath::Sin(phi_0);
     271      y_c = y - r * TMath::Cos(phi_0);
    277272      r_c = TMath::Hypot(x_c, y_c);
    278       phi_c = TMath::ATan(y_c/x_c);
    279       if(x_c < 0.0) phi_c -= TMath::Sign(1., phi_c)*TMath::Pi();
    280 
    281       //Find the time of closest approach
    282       td = (phi_0 - TMath::ATan(-x_c/y_c))/omega;
    283 
    284       //Remove all the modulo pi that might have come from the atan
    285       pio = fabs(TMath::Pi()/omega);
    286       while(fabs(td) > 0.5*pio)
    287       {
    288         td -= TMath::Sign(1., td)*pio;
    289       }
    290 
    291       //Compute the coordinate of closed approach to z axis
    292       //if wants wtr beamline need to be changedto re-center with a traslation of the z axis
    293       phid = phi_0 - omega*td;
    294       xd = x_c - r*TMath::Sin(phid);
    295       yd = y_c + r*TMath::Cos(phid);
    296       zd = z + c_light*(pz/e)*td;
    297 
    298       //Compute momentum at closest approach (perigee??)
    299       px = pt*TMath::Cos(phid);
    300       py = pt*TMath::Sin(phid);
     273
     274      // time of closest approach
     275      td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
     276
     277      // remove all the modulo pi that might have come from the atan
     278      pio = TMath::Abs(TMath::Pi() / omega);
     279      while(TMath::Abs(td) > 0.5 * pio)
     280      {
     281        td -= TMath::Sign(1.0, td) * pio;
     282      }
     283
     284      vz = pz * c_light / e;
     285
     286      // calculate coordinates of closest approach to z axis
     287      phid = phi_0 - omega * td;
     288      xd = x_c - r * TMath::Sin(phid);
     289      yd = y_c + r * TMath::Cos(phid);
     290      zd = z + vz * td;
     291
     292      // momentum at closest approach
     293      px = pt * TMath::Cos(phid);
     294      py = pt * TMath::Sin(phid);
    301295
    302296      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
     
    305299      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
    306300      dz = zd - bsz;
    307       ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
     301      ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
    308302
    309303      // 3. time evaluation t = TMath::Min(t_r, t_z)
    310304      //    t_r : time to exit from the sides
    311305      //    t_z : time to exit from the front or the back
    312       t = 0;
    313       t_z = 0;
    314       sign_pz = (pz > 0.0) ? 1 : -1;
    315       if(pz == 0.0) t_z = 1.0E99;
    316       else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
    317 
    318       if(r_c + TMath::Abs(r)  < fRadius)   // helix does not cross the cylinder sides
    319       {
     306      t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
     307
     308      if(r_c + TMath::Abs(r) < fRadius)
     309      {
     310        // helix does not cross the cylinder sides
    320311        t = t_z;
    321312      }
    322313      else
    323314      {
    324         alpha = -(fRadius*fRadius - r*r - r_c*r_c)/(2*fabs(r)*r_c);
    325         alpha = fabs(TMath::ACos(alpha));
    326         t_r = td + alpha/fabs(omega);
     315        alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
     316        t_r = td + TMath::Abs(alpha / omega);
    327317
    328318        t = TMath::Min(t_r, t_z);
    329319      }
    330320
    331       x_t = x_c - r*TMath::Sin(phi_0 - omega*t);
    332       y_t = y_c + r*TMath::Cos(phi_0 - omega*t);
    333       z_t = z + c_light*t*pz/e;
    334       r_t =  TMath::Hypot(x_t, y_t);
    335 
    336       // compute path length for an helix
    337       vz = pz*1.0E9 / c_light / gammam;
    338       //lenght of the path from production to tracker
    339       l = t * TMath::Sqrt(vz*vz + r*r*omega*omega);
     321      // 4. position in terms of x(t), y(t), z(t)
     322      phi_t = phi_0 - omega * t;
     323      x_t = x_c - r * TMath::Sin(phi_t);
     324      y_t = y_c + r * TMath::Cos(phi_t);
     325      z_t = z + vz * t;
     326      r_t = TMath::Hypot(x_t, y_t);
     327
     328      // lenght of the path from production to tracker
     329      l = t * TMath::Hypot(vz, r * omega);
    340330
    341331      if(r_t > 0.0)
Note: See TracChangeset for help on using the changeset viewer.