Changeset 341014c in git for modules/ParticlePropagator.cc
- Timestamp:
- Feb 12, 2019, 9:29:17 PM (6 years ago)
- Branches:
- ImprovedOutputFile, Timing, llp, master
- Children:
- 6455202
- Parents:
- 45e58be
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/ParticlePropagator.cc
r45e58be r341014c 17 17 */ 18 18 19 20 19 /** \class ParticlePropagator 21 20 * … … 35 34 #include "classes/DelphesFormula.h" 36 35 36 #include "ExRootAnalysis/ExRootClassifier.h" 37 #include "ExRootAnalysis/ExRootFilter.h" 37 38 #include "ExRootAnalysis/ExRootResult.h" 38 #include "ExRootAnalysis/ExRootFilter.h" 39 #include "ExRootAnalysis/ExRootClassifier.h" 40 39 40 #include "TDatabasePDG.h" 41 #include "TFormula.h" 42 #include "TLorentzVector.h" 41 43 #include "TMath.h" 44 #include "TObjArray.h" 45 #include "TRandom3.h" 42 46 #include "TString.h" 43 #include "TFormula.h"44 #include "TRandom3.h"45 #include "TObjArray.h"46 #include "TDatabasePDG.h"47 #include "TLorentzVector.h"48 47 49 48 #include <algorithm> 50 #include <stdexcept>51 49 #include <iostream> 52 50 #include <sstream> 51 #include <stdexcept> 53 52 54 53 using namespace std; … … 67 66 } 68 67 69 70 68 //------------------------------------------------------------------------------ 71 69 … … 73 71 { 74 72 fRadius = GetDouble("Radius", 1.0); 75 fRadius2 = fRadius *fRadius;73 fRadius2 = fRadius * fRadius; 76 74 fHalfLength = GetDouble("HalfLength", 3.0); 77 75 fBz = GetDouble("Bz", 0.0); … … 140 138 const Double_t c_light = 2.99792458E8; 141 139 142 if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)140 if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0) 143 141 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0); 144 142 else 145 143 { 146 Candidate &beamSpotCandidate = *((Candidate *) 144 Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0)); 147 145 beamSpotPosition = beamSpotCandidate.Position; 148 146 } 149 147 150 148 fItInputArray->Reset(); 151 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))149 while((candidate = static_cast<Candidate *>(fItInputArray->Next()))) 152 150 { 153 151 if(candidate->GetCandidates()->GetEntriesFast() == 0) … … 157 155 else 158 156 { 159 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));157 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); 160 158 } 161 159 162 160 particlePosition = particle->Position; 163 161 particleMomentum = particle->Momentum; 164 x = particlePosition.X() *1.0E-3;165 y = particlePosition.Y() *1.0E-3;166 z = particlePosition.Z() *1.0E-3;167 168 bsx = beamSpotPosition.X() *1.0E-3;169 bsy = beamSpotPosition.Y() *1.0E-3;170 bsz = beamSpotPosition.Z() *1.0E-3;162 x = particlePosition.X() * 1.0E-3; 163 y = particlePosition.Y() * 1.0E-3; 164 z = particlePosition.Z() * 1.0E-3; 165 166 bsx = beamSpotPosition.X() * 1.0E-3; 167 bsy = beamSpotPosition.Y() * 1.0E-3; 168 bsz = beamSpotPosition.Z() * 1.0E-3; 171 169 172 170 q = particle->Charge; … … 193 191 { 194 192 mother = candidate; 195 candidate = static_cast<Candidate *>(candidate->Clone());193 candidate = static_cast<Candidate *>(candidate->Clone()); 196 194 197 195 candidate->InitialPosition = particlePosition; … … 207 205 { 208 206 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0 209 tmp = px *y - py*x;210 discr2 = pt2 *fRadius2 - tmp*tmp;207 tmp = px * y - py * x; 208 discr2 = pt2 * fRadius2 - tmp * tmp; 211 209 212 210 if(discr2 < 0.0) … … 216 214 } 217 215 218 tmp = px *x + py*y;216 tmp = px * x + py * y; 219 217 discr = TMath::Sqrt(discr2); 220 t1 = (-tmp + discr) /pt2;221 t2 = (-tmp - discr) /pt2;218 t1 = (-tmp + discr) / pt2; 219 t2 = (-tmp - discr) / pt2; 222 220 t = (t1 < 0.0) ? t2 : t1; 223 221 224 z_t = z + pz *t;222 z_t = z + pz * t; 225 223 if(TMath::Abs(z_t) > fHalfLength) 226 224 { … … 230 228 } 231 229 232 x_t = x + px *t;233 y_t = y + py *t;234 z_t = z + pz *t;235 236 l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));230 x_t = x + px * t; 231 y_t = y + py * t; 232 z_t = z + pz * t; 233 234 l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z)); 237 235 238 236 mother = candidate; 239 candidate = static_cast<Candidate *>(candidate->Clone());237 candidate = static_cast<Candidate *>(candidate->Clone()); 240 238 241 239 candidate->InitialPosition = particlePosition; 242 candidate->Position.SetXYZT(x_t *1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*e*1.0E3);243 candidate->L = l *1.0E3;240 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3); 241 candidate->L = l * 1.0E3; 244 242 245 243 candidate->Momentum = particleMomentum; … … 251 249 switch(TMath::Abs(candidate->PID)) 252 250 { 253 254 255 256 257 258 259 260 251 case 11: 252 fElectronOutputArray->Add(candidate); 253 break; 254 case 13: 255 fMuonOutputArray->Add(candidate); 256 break; 257 default: 258 fChargedHadronOutputArray->Add(candidate); 261 259 } 262 260 } 263 261 else 264 262 { 265 263 fNeutralOutputArray->Add(candidate); 266 264 } 267 265 } … … 275 273 // helix radius r = p_{T0} / (omega gamma m) 276 274 277 gammam = e *1.0E9 / (c_light*c_light);// gammam in [eV/c^2]278 omega = q * fBz / (gammam); 279 r = pt / (q * fBz) * 1.0E9 /c_light;// in [m]275 gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2] 276 omega = q * fBz / (gammam); // omega is here in [89875518/s] 277 r = pt / (q * fBz) * 1.0E9 / c_light; // in [m] 280 278 281 279 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi] 282 280 283 281 // 2. helix axis coordinates 284 x_c = x + r *TMath::Sin(phi_0);285 y_c = y - r *TMath::Cos(phi_0);282 x_c = x + r * TMath::Sin(phi_0); 283 y_c = y - r * TMath::Cos(phi_0); 286 284 r_c = TMath::Hypot(x_c, y_c); 287 285 phi_c = TMath::ATan2(y_c, x_c); … … 290 288 291 289 rcu = TMath::Abs(r); 292 rc2 = r_c *r_c;290 rc2 = r_c * r_c; 293 291 294 292 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd 295 xd = x_c *x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;293 xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c; 296 294 xd = (rc2 > 0.0) ? xd / rc2 : -999; 297 yd = y_c *(-rcu*r_c + rc2);295 yd = y_c * (-rcu * r_c + rc2); 298 296 yd = (rc2 > 0.0) ? yd / rc2 : -999; 299 zd = z + (TMath::Sqrt(xd *xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;297 zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt; 300 298 301 299 // use perigee momentum rather than original particle … … 311 309 // calculate additional track parameters (correct for beamspot position) 312 310 313 d0 = ((x - bsx) * py - (y - bsy) * px) / pt; 314 dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt); 315 p = particleMomentum.P(); 316 ctgTheta = 1.0 / TMath::Tan (particleMomentum.Theta()); 317 311 d0 = ((x - bsx) * py - (y - bsy) * px) / pt; 312 dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt); 313 p = particleMomentum.P(); 314 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta()); 318 315 319 316 // 3. time evaluation t = TMath::Min(t_r, t_z) … … 322 319 t_r = 0.0; // in [ns] 323 320 int sign_pz = (pz > 0.0) ? 1 : -1; 324 if(pz == 0.0) t_z = 1.0E99; 325 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz); 326 327 if(r_c + TMath::Abs(r) < fRadius) 321 if(pz == 0.0) 322 t_z = 1.0E99; 323 else 324 t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz); 325 326 if(r_c + TMath::Abs(r) < fRadius) 328 327 { 329 328 // helix does not cross the cylinder sides … … 332 331 else 333 332 { 334 asinrho = TMath::ASin((fRadius *fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c));333 asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c)); 335 334 delta = phi_0 - phi; 336 if(delta < -TMath::Pi()) delta += 2*TMath::Pi();337 if(delta > TMath::Pi()) delta -= 2 *TMath::Pi();335 if(delta < -TMath::Pi()) delta += 2 * TMath::Pi(); 336 if(delta > TMath::Pi()) delta -= 2 * TMath::Pi(); 338 337 t1 = (delta + asinrho) / omega; 339 338 t2 = (delta + TMath::Pi() - asinrho) / omega; … … 359 358 x_t = x_c + r * TMath::Sin(omega * t - phi_0); 360 359 y_t = y_c + r * TMath::Cos(omega * t - phi_0); 361 z_t = z + pz *1.0E9 / c_light / gammam * t;360 z_t = z + pz * 1.0E9 / c_light / gammam * t; 362 361 r_t = TMath::Hypot(x_t, y_t); 363 362 364 365 363 // compute path length for an helix 366 364 367 alpha = pz *1.0E9 / c_light / gammam;368 l = t * TMath::Sqrt(alpha *alpha + r*r*omega*omega);365 alpha = pz * 1.0E9 / c_light / gammam; 366 l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega); 369 367 370 368 if(r_t > 0.0) … … 374 372 if(particle == candidate) 375 373 { 376 particle->D0 = d0 *1.0E3;377 particle->DZ = dz *1.0E3;374 particle->D0 = d0 * 1.0E3; 375 particle->DZ = dz * 1.0E3; 378 376 particle->P = p; 379 377 particle->PT = pt; … … 383 381 384 382 mother = candidate; 385 candidate = static_cast<Candidate *>(candidate->Clone());383 candidate = static_cast<Candidate *>(candidate->Clone()); 386 384 387 385 candidate->InitialPosition = particlePosition; 388 candidate->Position.SetXYZT(x_t *1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);386 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * c_light * 1.0E3); 389 387 390 388 candidate->Momentum = particleMomentum; 391 389 392 candidate->L = l *1.0E3;393 394 candidate->Xd = xd *1.0E3;395 candidate->Yd = yd *1.0E3;396 candidate->Zd = zd *1.0E3;390 candidate->L = l * 1.0E3; 391 392 candidate->Xd = xd * 1.0E3; 393 candidate->Yd = yd * 1.0E3; 394 candidate->Zd = zd * 1.0E3; 397 395 398 396 candidate->AddCandidate(mother); … … 401 399 switch(TMath::Abs(candidate->PID)) 402 400 { 403 404 405 406 407 408 409 410 401 case 11: 402 fElectronOutputArray->Add(candidate); 403 break; 404 case 13: 405 fMuonOutputArray->Add(candidate); 406 break; 407 default: 408 fChargedHadronOutputArray->Add(candidate); 411 409 } 412 410 } … … 416 414 417 415 //------------------------------------------------------------------------------ 418
Note:
See TracChangeset
for help on using the changeset viewer.