Changeset d41ba4a in git for modules/ParticlePropagator.cc

Timestamp:

Dec 20, 2013, 11:30:49 PM (11 years ago)

Author:

pavel <pavel@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

41f30b4

Parents:

8555f6d

Message:

fix time calculation in ParticlePropagator

File:

• modules/ParticlePropagator.cc (modified) (12 diffs)

modules/ParticlePropagator.cc

@@ -3 +3 @@
  *
  *  Propagates charged and neutral particles
- *  from a given vertex to a cylinder defined by its radius, 
+ *  from a given vertex to a cylinder defined by its radius,
  *  its half-length, centered at (0,0,0) and with its axis
  *  oriented along the z-axis.
@@ -33 +33 @@
 #include "TLorentzVector.h"
 
-#include <algorithm> 
+#include <algorithm>
 #include <stdexcept>
 #include <iostream>
@@ -62 +62 @@
   fBz = GetDouble("Bz", 0.0);
   if(fRadius < 1.0E-2)
-  { 
+  {
     cout << "ERROR: magnetic field radius is too low\n";
     return;
@@ -106 +106 @@
   Double_t tmp, discr, discr2;
   Double_t delta, gammam, omega, asinrho;
-  
+
   const Double_t c_light = 2.99792458E8;
-    
+
   fItInputArray->Reset();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
@@ -142 +142 @@
       tmp = px*y - py*x;
       discr2 = pt2*fRadius2 - tmp*tmp;
-      
+
       if(discr2 < 0)
       {
@@ -153 +153 @@
       t1 = (-tmp + discr)/pt2;
       t2 = (-tmp - discr)/pt2;
-      t = (t1 < 0) ? t2 : t1; 
+      t = (t1 < 0) ? t2 : t1;
 
       z_t = z + pz*t;
@@ -160 +160 @@
         t3 = (+fHalfLength - z) / pz;
         t4 = (-fHalfLength - z) / pz;
-        t = (t3 < 0) ? t4 : t3; 
+        t = (t3 < 0) ? t4 : t3;
       }
 
@@ -170 +170 @@
       candidate = static_cast<Candidate*>(candidate->Clone());
 
-      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t);
+      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
 
       candidate->Momentum = candidateMomentum;
       candidate->AddCandidate(mother);
-      
+
       fOutputArray->Add(candidate);
-      if(TMath::Abs(q) > 1.0E-9) 
+      if(TMath::Abs(q) > 1.0E-9)
       {
         switch(TMath::Abs(candidate->PID))
@@ -195 +195 @@
 
       // 1.  initial transverse momentum p_{T0} : Part->pt
-      //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) 
+      //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
       //     relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
       //     giration frequency \omega = q/(gamma m) fBz
@@ -201 +201 @@
 
       gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c²]
-      omega = q * fBz / (gammam);                // omega is here in [ 89875518 / s]  
-      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]  
+      omega = q * fBz / (gammam);                // omega is here in [ 89875518 / s]
+      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
 
       phi_0 = TMath::ATan2(py, px); // [rad] in [-pi; pi]
@@ -250 +250 @@
         t_rb = TMath::Min(t4, TMath::Min(t5, t6));
         t_r = TMath::Min(t_ra, t_rb);
-        t = TMath::Min(t_r, t_z); 
+        t = TMath::Min(t_r, t_z);
       }
 
@@ -264 +264 @@
         candidate = static_cast<Candidate*>(candidate->Clone());
 
-        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t);
+        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
 
         candidate->Momentum = candidateMomentum;