Fork me on GitHub

Changeset 38b4e15 in git for modules


Ignore:
Timestamp:
May 17, 2021, 5:05:51 PM (3 years ago)
Author:
michele <michele.selvaggi@…>
Branches:
master
Children:
6216fb2
Parents:
914d1be
Message:

fixed bug in propagation at low pT

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/ParticlePropagator.cc

    r914d1be r38b4e15  
    125125  TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
    126126  Double_t px, py, pz, pt, pt2, e, q;
    127   Double_t x, y, z, t, r, phi;
     127  Double_t x, y, z, t, r;
    128128  Double_t x_c, y_c, r_c, phi_c, phi_0;
    129129  Double_t x_t, y_t, z_t, r_t;
    130   Double_t t1, t2, t3, t4, t5, t6;
    131   Double_t t_z, t_r, t_ra, t_rb;
    132   Double_t tmp, discr, discr2;
    133   Double_t delta, gammam, omega, asinrho;
    134   Double_t rcu, rc2, xd, yd, zd;
    135   Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
     130  Double_t t_z, t_r;
     131  Double_t discr;
     132  Double_t gammam, omega;
     133  Double_t xd, yd, zd;
     134  Double_t l, d0, dz, ctgTheta, alpha;
    136135  Double_t bsx, bsy, bsz;
    137   Double_t s0, s1, sd;
    138  
     136  Double_t rxp, rdp, t_R;
     137  Double_t td, pio, phid, sign_pz, vz;
     138
    139139  const Double_t c_light = 2.99792458E8;
    140140
     
    161161    particlePosition = particle->Position;
    162162    particleMomentum = particle->Momentum;
     163
     164    // Constants
     165
    163166    x = particlePosition.X() * 1.0E-3;
    164167    y = particlePosition.Y() * 1.0E-3;
     
    205208    else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
    206209    {
    207       // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
    208       tmp = px * y - py * x;
    209       discr2 = pt2 * fRadius2 - tmp * tmp;
    210 
    211       if(discr2 < 0.0)
    212       {
    213         // no solutions
    214         continue;
    215       }
    216 
    217       tmp = px * x + py * y;
    218       discr = TMath::Sqrt(discr2);
    219       t1 = (-tmp + discr) / pt2;
    220       t2 = (-tmp - discr) / pt2;
    221       t = (t1 < 0.0) ? t2 : t1;
    222 
    223       z_t = z + pz * t;
    224       if(TMath::Abs(z_t) > fHalfLength)
    225       {
    226         t3 = (+fHalfLength - z) / pz;
    227         t4 = (-fHalfLength - z) / pz;
    228         t = (t3 < 0.0) ? t4 : t3;
    229       }
    230 
    231       x_t = x + px * t;
    232       y_t = y + py * t;
    233       z_t = z + pz * t;
    234 
    235       l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
     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));
    236227
    237228      mother = candidate;
    238       candidate = static_cast<Candidate *>(candidate->Clone());
     229      candidate = static_cast<Candidate*>(candidate->Clone());
    239230
    240231      candidate->InitialPosition = particlePosition;
    241       candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
    242       candidate->L = l * 1.0E3;
     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;
    243234
    244235      candidate->Momentum = particleMomentum;
     
    246237
    247238      fOutputArray->Add(candidate);
     239
    248240      if(TMath::Abs(q) > 1.0E-9)
    249241      {
     
    274266      //     helix radius r = p_{T0} / (omega gamma m)
    275267
    276       gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
    277       omega = q * fBz / (gammam); // omega is here in [89875518/s]
    278       r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
     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]
    279271
    280272      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
    281273
    282274      // 2. helix axis coordinates
    283       x_c = x + r * TMath::Sin(phi_0);
    284       y_c = y - r * TMath::Cos(phi_0);
     275      x_c = x + r*TMath::Sin(phi_0);
     276      y_c = y - r*TMath::Cos(phi_0);
    285277      r_c = TMath::Hypot(x_c, y_c);
    286       phi_c = TMath::ATan2(y_c, x_c);
    287       phi = phi_c;
    288       if(x_c < 0.0) phi += TMath::Pi();
    289 
    290       rcu = TMath::Abs(r);
    291       rc2 = r_c * r_c;
    292 
    293       // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
    294       xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
    295       xd = (rc2 > 0.0) ? xd / rc2 : -999;
    296       yd = y_c * (-rcu * r_c + rc2);
    297       yd = (rc2 > 0.0) ? yd / rc2 : -999;
    298       zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
    299 
    300       // proper calculation of the DCAz coordinate
    301       // s0: track circle parameter at the track origin
    302       // s1: track circle parameter at the closest approach to beam pipe
    303       // sd: s1-s0 signed angular difference
    304       s0 = atan2(y - y_c, x - x_c);
    305       s1 = atan2(yd - y_c, xd - x_c);
    306       sd = atan2(sin(s1 - s0), cos(s1 - s0));
    307       zd = z - r * pz / pt * sd;
    308      
    309       // use perigee momentum rather than original particle
    310       // momentum, since the orignal particle momentum isn't known
    311 
    312       px = TMath::Sign(1.0, r) * pt * (-y_c / r_c);
    313       py = TMath::Sign(1.0, r) * pt * (x_c / r_c);
    314       etap = particleMomentum.Eta();
    315       phip = TMath::ATan2(py, px);
    316 
    317       particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
     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);
     301
     302      particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
    318303
    319304      // calculate additional track parameters (correct for beamspot position)
    320 
    321       d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
    322       dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
    323       p = particleMomentum.P();
    324       ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
     305      d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
     306      dz = zd - bsz;
     307      ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
    325308
    326309      // 3. time evaluation t = TMath::Min(t_r, t_z)
    327310      //    t_r : time to exit from the sides
    328311      //    t_z : time to exit from the front or the back
    329       t_r = 0.0; // in [ns]
    330       int sign_pz = (pz > 0.0) ? 1 : -1;
    331       if(pz == 0.0)
    332         t_z = 1.0E99;
     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      {
     320        t = t_z;
     321      }
    333322      else
    334         t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
    335 
    336       if(r_c + TMath::Abs(r) < fRadius)
    337       {
    338         // helix does not cross the cylinder sides
    339         t = t_z;
    340       }
    341       else
    342       {
    343         asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
    344         delta = phi_0 - phi;
    345         if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
    346         if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
    347         t1 = (delta + asinrho) / omega;
    348         t2 = (delta + TMath::Pi() - asinrho) / omega;
    349         t3 = (delta + TMath::Pi() + asinrho) / omega;
    350         t4 = (delta - asinrho) / omega;
    351         t5 = (delta - TMath::Pi() - asinrho) / omega;
    352         t6 = (delta - TMath::Pi() + asinrho) / omega;
    353 
    354         if(t1 < 0.0) t1 = 1.0E99;
    355         if(t2 < 0.0) t2 = 1.0E99;
    356         if(t3 < 0.0) t3 = 1.0E99;
    357         if(t4 < 0.0) t4 = 1.0E99;
    358         if(t5 < 0.0) t5 = 1.0E99;
    359         if(t6 < 0.0) t6 = 1.0E99;
    360 
    361         t_ra = TMath::Min(t1, TMath::Min(t2, t3));
    362         t_rb = TMath::Min(t4, TMath::Min(t5, t6));
    363         t_r = TMath::Min(t_ra, t_rb);
     323      {
     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);
     327
    364328        t = TMath::Min(t_r, t_z);
    365329      }
    366330
    367       // 4. position in terms of x(t), y(t), z(t)
    368       x_t = x_c + r * TMath::Sin(omega * t - phi_0);
    369       y_t = y_c + r * TMath::Cos(omega * t - phi_0);
    370       z_t = z + pz * 1.0E9 / c_light / gammam * t;
    371       r_t = TMath::Hypot(x_t, y_t);
     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);
    372335
    373336      // compute path length for an helix
    374 
    375       alpha = pz * 1.0E9 / c_light / gammam;
    376       l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
     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);
    377340
    378341      if(r_t > 0.0)
    379342      {
    380 
    381343        // store these variables before cloning
    382344        if(particle == candidate)
     
    384346          particle->D0 = d0 * 1.0E3;
    385347          particle->DZ = dz * 1.0E3;
    386           particle->P = p;
     348          particle->P = particleMomentum.P();
    387349          particle->PT = pt;
    388350          particle->CtgTheta = ctgTheta;
    389           particle->Phi = phip;
     351          particle->Phi = particleMomentum.Phi();
    390352        }
    391353
Note: See TracChangeset for help on using the changeset viewer.