Changeset 38b4e15 in git
- Timestamp:
- May 17, 2021, 5:05:51 PM (3 years ago)
- Branches:
- master
- Children:
- 6216fb2
- Parents:
- 914d1be
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/ParticlePropagator.cc
r914d1be r38b4e15 125 125 TLorentzVector particlePosition, particleMomentum, beamSpotPosition; 126 126 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; 128 128 Double_t x_c, y_c, r_c, phi_c, phi_0; 129 129 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; 136 135 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 139 139 const Double_t c_light = 2.99792458E8; 140 140 … … 161 161 particlePosition = particle->Position; 162 162 particleMomentum = particle->Momentum; 163 164 // Constants 165 163 166 x = particlePosition.X() * 1.0E-3; 164 167 y = particlePosition.Y() * 1.0E-3; … … 205 208 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9) 206 209 { 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)); 236 227 237 228 mother = candidate; 238 candidate = static_cast<Candidate 229 candidate = static_cast<Candidate*>(candidate->Clone()); 239 230 240 231 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; 243 234 244 235 candidate->Momentum = particleMomentum; … … 246 237 247 238 fOutputArray->Add(candidate); 239 248 240 if(TMath::Abs(q) > 1.0E-9) 249 241 { … … 274 266 // helix radius r = p_{T0} / (omega gamma m) 275 267 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] 279 271 280 272 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi] 281 273 282 274 // 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); 285 277 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()); 318 303 319 304 // 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()); 325 308 326 309 // 3. time evaluation t = TMath::Min(t_r, t_z) 327 310 // t_r : time to exit from the sides 328 311 // 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 } 333 322 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 364 328 t = TMath::Min(t_r, t_z); 365 329 } 366 330 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); 372 335 373 336 // 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); 377 340 378 341 if(r_t > 0.0) 379 342 { 380 381 343 // store these variables before cloning 382 344 if(particle == candidate) … … 384 346 particle->D0 = d0 * 1.0E3; 385 347 particle->DZ = dz * 1.0E3; 386 particle->P = p ;348 particle->P = particleMomentum.P(); 387 349 particle->PT = pt; 388 350 particle->CtgTheta = ctgTheta; 389 particle->Phi = p hip;351 particle->Phi = particleMomentum.Phi(); 390 352 } 391 353
Note:
See TracChangeset
for help on using the changeset viewer.