Fork me on GitHub

Changeset d41ba4a in git


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/ParticlePropagator.cc

    r8555f6d rd41ba4a  
    33 *
    44 *  Propagates charged and neutral particles
    5  *  from a given vertex to a cylinder defined by its radius, 
     5 *  from a given vertex to a cylinder defined by its radius,
    66 *  its half-length, centered at (0,0,0) and with its axis
    77 *  oriented along the z-axis.
     
    3333#include "TLorentzVector.h"
    3434
    35 #include <algorithm> 
     35#include <algorithm>
    3636#include <stdexcept>
    3737#include <iostream>
     
    6262  fBz = GetDouble("Bz", 0.0);
    6363  if(fRadius < 1.0E-2)
    64   { 
     64  {
    6565    cout << "ERROR: magnetic field radius is too low\n";
    6666    return;
     
    106106  Double_t tmp, discr, discr2;
    107107  Double_t delta, gammam, omega, asinrho;
    108  
     108
    109109  const Double_t c_light = 2.99792458E8;
    110    
     110
    111111  fItInputArray->Reset();
    112112  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     
    142142      tmp = px*y - py*x;
    143143      discr2 = pt2*fRadius2 - tmp*tmp;
    144      
     144
    145145      if(discr2 < 0)
    146146      {
     
    153153      t1 = (-tmp + discr)/pt2;
    154154      t2 = (-tmp - discr)/pt2;
    155       t = (t1 < 0) ? t2 : t1; 
     155      t = (t1 < 0) ? t2 : t1;
    156156
    157157      z_t = z + pz*t;
     
    160160        t3 = (+fHalfLength - z) / pz;
    161161        t4 = (-fHalfLength - z) / pz;
    162         t = (t3 < 0) ? t4 : t3; 
     162        t = (t3 < 0) ? t4 : t3;
    163163      }
    164164
     
    170170      candidate = static_cast<Candidate*>(candidate->Clone());
    171171
    172       candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t);
     172      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
    173173
    174174      candidate->Momentum = candidateMomentum;
    175175      candidate->AddCandidate(mother);
    176      
     176
    177177      fOutputArray->Add(candidate);
    178       if(TMath::Abs(q) > 1.0E-9) 
     178      if(TMath::Abs(q) > 1.0E-9)
    179179      {
    180180        switch(TMath::Abs(candidate->PID))
     
    195195
    196196      // 1.  initial transverse momentum p_{T0} : Part->pt
    197       //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) 
     197      //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
    198198      //     relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
    199199      //     giration frequency \omega = q/(gamma m) fBz
     
    201201
    202202      gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c²]
    203       omega = q * fBz / (gammam);                // omega is here in [ 89875518 / s] 
    204       r = pt / (q * fBz) * 1.0E9/c_light;        // in [m] 
     203      omega = q * fBz / (gammam);                // omega is here in [ 89875518 / s]
     204      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
    205205
    206206      phi_0 = TMath::ATan2(py, px); // [rad] in [-pi; pi]
     
    250250        t_rb = TMath::Min(t4, TMath::Min(t5, t6));
    251251        t_r = TMath::Min(t_ra, t_rb);
    252         t = TMath::Min(t_r, t_z); 
     252        t = TMath::Min(t_r, t_z);
    253253      }
    254254
     
    264264        candidate = static_cast<Candidate*>(candidate->Clone());
    265265
    266         candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t);
     266        candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
    267267
    268268        candidate->Momentum = candidateMomentum;
Note: See TracChangeset for help on using the changeset viewer.