Changeset ae93700 in git
- Timestamp:
- May 27, 2021, 8:36:27 AM (3 years ago)
- Branches:
- master
- Children:
- 716e186, dfdad8b
- Parents:
- c6f9311
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/ParticlePropagator.cc
rc6f9311 rae93700 126 126 Double_t px, py, pz, pt, pt2, e, q; 127 127 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; 132 132 Double_t gammam, omega; 133 133 Double_t xd, yd, zd; 134 134 Double_t l, d0, dz, ctgTheta, alpha; 135 135 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; 138 137 139 138 const Double_t c_light = 2.99792458E8; 140 139 141 140 if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0) 141 { 142 142 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0); 143 } 143 144 else 144 145 { … … 161 162 particlePosition = particle->Position; 162 163 particleMomentum = particle->Momentum; 163 164 // Constants165 164 166 165 x = particlePosition.X() * 1.0E-3; … … 208 207 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9) 209 208 { 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)); 227 222 228 223 mother = candidate; 229 candidate = static_cast<Candidate *>(candidate->Clone());224 candidate = static_cast<Candidate *>(candidate->Clone()); 230 225 231 226 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; 234 229 235 230 candidate->Momentum = particleMomentum; … … 260 255 { 261 256 262 // 1. 263 // initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)264 // relativistic gamma: gamma = E/mc^2; gammam = gamma * m265 // gyration frequency omega = q/(gamma m) fBz266 // helix radius r = p_{T0} / (omega gammam)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] 271 266 272 267 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi] 273 268 274 269 // 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); 277 272 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); 301 295 302 296 particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E()); … … 305 299 d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt; 306 300 dz = zd - bsz; 307 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());301 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta()); 308 302 309 303 // 3. time evaluation t = TMath::Min(t_r, t_z) 310 304 // t_r : time to exit from the sides 311 305 // 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 320 311 t = t_z; 321 312 } 322 313 else 323 314 { 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); 327 317 328 318 t = TMath::Min(t_r, t_z); 329 319 } 330 320 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 helix337 vz = pz*1.0E9 / c_light / gammam; 338 // lenght of the path from production to tracker339 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); 340 330 341 331 if(r_t > 0.0)
Note:
See TracChangeset
for help on using the changeset viewer.