Changes in / [56af73f:90132e0] in git

Files:

• DelphesEnv.sh (modified) (1 diff)
• Makefile (modified) (1 diff)
• README (modified) (1 diff)
• README.md (modified) (2 diffs)
• classes/DelphesHepMC3Reader.cc (modified) (14 diffs)
• classes/DelphesHepMC3Reader.h (modified) (4 diffs)
• modules/Delphes.cc (modified) (2 diffs)
• modules/ParticlePropagator.cc (modified) (5 diffs)

DelphesEnv.sh

@@ -1 @@
-export PYTHONPATH=`pwd`/python:$PYTHONPATH
-export LD_LIBRARY_PATH=`pwd`:$LD_LIBRARY_PATH
+DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )"
 
+export DELPHES_HOME="$DIR"
+export PYTHONPATH="$DIR/python:${PYTHONPATH}"
+export LD_LIBRARY_PATH="$DIR:${LD_LIBRARY_PATH}"
+export LIBRARY_PATH="$DIR:${LIBRARY_PATH}"

Makefile

@@ -232 +232 @@
         classes/DelphesClasses.h \
         classes/DelphesFactory.h \
-        classes/DelphesHepMC3Reader.h \
         modules/Delphes.h \
         external/ExRootAnalysis/ExRootProgressBar.h \

README

@@ -16 +16 @@
 Finally, we can run Delphes:
 
-   ./DelphesHepMC
+   ./DelphesHepMC3
 
 Command line parameters:
 
-   ./DelphesHepMC config_file output_file [input_file(s)]
+   ./DelphesHepMC3 config_file output_file [input_file(s)]
      config_file - configuration file in Tcl format
      output_file - output file in ROOT format,

README.md

@@ -30 +30 @@
 
 ```
-   ./DelphesHepMC
+   ./DelphesHepMC3
 ```
 
@@ -36 +36 @@
 
 ```
-   ./DelphesHepMC config_file output_file [input_file(s)]
+   ./DelphesHepMC3 config_file output_file [input_file(s)]
      config_file - configuration file in Tcl format
      output_file - output file in ROOT format,

classes/DelphesHepMC3Reader.cc

@@ -56 +56 @@
 DelphesHepMC3Reader::DelphesHepMC3Reader() :
   fInputFile(0), fBuffer(0), fPDG(0),
-  fVertexCounter(-1), fParticleCounter(-1)
+  fVertexCounter(-2), fParticleCounter(-1)
 {
   fBuffer = new char[kBufferSize];
@@ -81 +81 @@
 void DelphesHepMC3Reader::Clear()
 {
-  fWeight.clear();
+  fWeights.clear();
   fMomentumCoefficient = 1.0;
   fPositionCoefficient = 1.0;
-  fVertexCounter = -1;
+  fVertexCounter = -2;
   fParticleCounter = -1;
-  fVertexMap.clear();
+  fVertices.clear();
+  fParticles.clear();
+  fInVertexMap.clear();
+  fOutVertexMap.clear();
+  fMotherMap.clear();
   fDaughterMap.clear();
 }
@@ -94 +98 @@
 bool DelphesHepMC3Reader::EventReady()
 {
-  return (fParticleCounter == 0);
+  return (fVertexCounter == -1) && (fParticleCounter == 0);
 }
 
@@ -119 +123 @@
     Clear();
 
-    fX = 0.0;
-    fY = 0.0;
-    fZ = 0.0;
-    fT = 0.0;
-
     rc = bufferStream.ReadInt(fEventNumber)
       && bufferStream.ReadInt(fVertexCounter)
@@ -168 +167 @@
     while(bufferStream.ReadDbl(weight))
     {
-      fWeight.push_back(weight);
+      fWeights.push_back(weight);
     }
   }
@@ -257 +256 @@
   else if(key == 'V')
   {
+    fParticles.clear();
+
     fX = 0.0;
     fY = 0.0;
@@ -262 +263 @@
     fT = 0.0;
 
-    fM1 = 0;
-    fM2 = 0;
-
     rc = bufferStream.ReadInt(fVertexCode)
       && bufferStream.ReadInt(fVertexStatus);
@@ -286 +284 @@
     while(bufferStream.ReadInt(code))
     {
-      if(code < fM1 || fM1 == 0) fM1 = code;
-      if(code > fM2) fM2 = code;
-      fVertexMap[code] = fVertexCode;
-    }
-
-    if(fM1 == fM2) fM2 = 0;
+      fParticles.push_back(code);
+      bufferStream.FindChr(',');
+    }
 
     if(bufferStream.FindChr('@'))
@@ -307 +302 @@
       }
     }
+
+    AnalyzeVertex(factory, fVertexCode);
   }
   else if(key == 'P' && fParticleCounter > 0)
@@ -329 +326 @@
     }
 
-    itDaughterMap = fDaughterMap.find(fOutVertexCode);
-    if(itDaughterMap == fDaughterMap.end())
-    {
-      fDaughterMap[fOutVertexCode] = make_pair(fParticleCode, fParticleCode);
-    }
-    else
-    {
-      itDaughterMap->second.second = fParticleCode;
-    }
-
-    AnalyzeParticle(factory, allParticleOutputArray,
-      stableParticleOutputArray, partonOutputArray);
+    AnalyzeParticle(factory);
   }
 
   if(EventReady())
   {
-    FinalizeParticles(allParticleOutputArray);
+    FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
   }
 
@@ -363 +349 @@
   element->ProcessID = fProcessID;
   element->MPI = fMPI;
-  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
+  element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
   element->CrossSection = fCrossSection;
   element->CrossSectionError = fCrossSectionError;
@@ -387 +373 @@
 {
   Weight *element;
-  vector<double>::const_iterator itWeight;
-
-  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
+  vector<double>::const_iterator itWeights;
+
+  for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
   {
     element = static_cast<Weight *>(branch->NewEntry());
 
-    element->Weight = *itWeight;
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory,
-  TObjArray *allParticleOutputArray,
+    element->Weight = *itWeights;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate)
+{
+  int index;
+  TLorentzVector *position;
+  TObjArray *array;
+  vector<int>::iterator itParticle;
+  map<int, int>::iterator itVertexMap;
+
+  itVertexMap = fOutVertexMap.find(code);
+  if(itVertexMap == fOutVertexMap.end())
+  {
+    --fVertexCounter;
+
+    index = fVertices.size();
+    fOutVertexMap[code] = index;
+    if(candidate && code > 0) fInVertexMap[code] = index;
+
+    position = factory->New<TLorentzVector>();
+    array = factory->NewArray();
+    position->SetXYZT(0.0, 0.0, 0.0, 0.0);
+    fVertices.push_back(make_pair(position, array));
+  }
+  else
+  {
+    index = itVertexMap->second;
+    position = fVertices[index].first;
+    array = fVertices[index].second;
+  }
+
+  if(candidate)
+  {
+    array->Add(candidate);
+  }
+  else
+  {
+    position->SetXYZT(fX, fY, fZ, fT);
+    for(itParticle = fParticles.begin(); itParticle != fParticles.end(); ++itParticle)
+    {
+      fInVertexMap[*itParticle] = index;
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory)
+{
+  Candidate *candidate;
+
+  candidate = factory->NewCandidate();
+
+  candidate->PID = fPID;
+
+  candidate->Status = fParticleStatus;
+
+  candidate->Mass = fMass;
+
+  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
+
+  candidate->D1 = fParticleCode;
+
+  AnalyzeVertex(factory, fOutVertexCode, candidate);
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray,
   TObjArray *stableParticleOutputArray,
   TObjArray *partonOutputArray)
 {
+  TLorentzVector *position;
+  TObjArray *array;
   Candidate *candidate;
   TParticlePDG *pdgParticle;
   int pdgCode;
-
-  candidate = factory->NewCandidate();
-
-  candidate->PID = fPID;
-  pdgCode = TMath::Abs(candidate->PID);
-
-  candidate->Status = fParticleStatus;
-
-  pdgParticle = fPDG->GetParticle(fPID);
-  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
-  candidate->Mass = fMass;
-
-  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
-  if(fMomentumCoefficient != 1.0)
-  {
-    candidate->Momentum *= fMomentumCoefficient;
-  }
-
-  candidate->Position.SetXYZT(fX, fY, fZ, fT);
-  if(fPositionCoefficient != 1.0)
-  {
-    candidate->Position *= fPositionCoefficient;
-  }
-
-  candidate->D1 = -1;
-  candidate->D2 = -1;
-
-  if(fOutVertexCode < 0)
-  {
-    candidate->M1 = fM1 - 1;
-    candidate->M2 = fM2 - 1;
-  }
-  else
-  {
-    candidate->M1 = fOutVertexCode - 1;
-    candidate->M2 = -1;
-  }
-
-  allParticleOutputArray->Add(candidate);
-
-  if(!pdgParticle) return;
-
-  if(fParticleStatus == 1)
-  {
-    stableParticleOutputArray->Add(candidate);
-  }
-  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
-  {
-    partonOutputArray->Add(candidate);
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray)
-{
-  Candidate *candidate;
   map<int, int >::iterator itVertexMap;
+  map<int, pair<int, int> >::iterator itMotherMap;
   map<int, pair<int, int> >::iterator itDaughterMap;
-  int i, index;
+  int i, j, code, counter;
+
+  counter = 0;
+  for(i = 0; i < fVertices.size(); ++i)
+  {
+    position = fVertices[i].first;
+    array = fVertices[i].second;
+
+    for(j = 0; j < array->GetEntriesFast(); ++j)
+    {
+      candidate = static_cast<Candidate *>(array->At(j));
+
+      candidate->Position = *position;
+      if(fPositionCoefficient != 1.0)
+      {
+        candidate->Position *= fPositionCoefficient;
+      }
+
+      if(fMomentumCoefficient != 1.0)
+      {
+        candidate->Momentum *= fMomentumCoefficient;
+      }
+
+      candidate->M1 = i;
+
+      itDaughterMap = fDaughterMap.find(i);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        fDaughterMap[i] = make_pair(counter, counter);
+      }
+      else
+      {
+        itDaughterMap->second.second = counter;
+      }
+
+      code = candidate->D1;
+
+      itVertexMap = fInVertexMap.find(code);
+      if(itVertexMap == fInVertexMap.end())
+      {
+        candidate->D1 = -1;
+      }
+      else
+      {
+        code = itVertexMap->second;
+
+        candidate->D1 = code;
+
+        itMotherMap = fMotherMap.find(code);
+        if(itMotherMap == fMotherMap.end())
+        {
+          fMotherMap[code] = make_pair(counter, -1);
+        }
+        else
+        {
+          itMotherMap->second.second = counter;
+        }
+      }
+
+      allParticleOutputArray->Add(candidate);
+
+      ++counter;
+
+      pdgParticle = fPDG->GetParticle(candidate->PID);
+
+      candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
+
+      if(!pdgParticle) continue;
+
+      pdgCode = TMath::Abs(candidate->PID);
+
+      if(candidate->Status == 1)
+      {
+        stableParticleOutputArray->Add(candidate);
+      }
+      else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+      {
+        partonOutputArray->Add(candidate);
+      }
+    }
+  }
 
   for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
@@ -472 +549 @@
     candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
 
-    index = i + 1;
-
-    itVertexMap = fVertexMap.find(index);
-    if(itVertexMap != fVertexMap.end())
-    {
-      index = itVertexMap->second;
-    }
-
-    itDaughterMap = fDaughterMap.find(index);
-    if(itDaughterMap == fDaughterMap.end())
+    itMotherMap = fMotherMap.find(candidate->M1);
+    if(itMotherMap == fMotherMap.end())
+    {
+      candidate->M1 = -1;
+      candidate->M2 = -1;
+    }
+    else
+    {
+      candidate->M1 = itMotherMap->second.first;
+      candidate->M2 = itMotherMap->second.second;
+    }
+
+    if(candidate->D1 < 0)
     {
       candidate->D1 = -1;
@@ -488 +568 @@
     else
     {
-      candidate->D1 = itDaughterMap->second.first - 1;
-      candidate->D2 = itDaughterMap->second.second - 1;
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
+      itDaughterMap = fDaughterMap.find(candidate->D1);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        candidate->D1 = -1;
+        candidate->D2 = -1;
+      }
+      else
+      {
+        candidate->D1 = itDaughterMap->second.first;
+        candidate->D2 = itDaughterMap->second.second;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------

classes/DelphesHepMC3Reader.h

@@ -36 +36 @@
 class TStopwatch;
 class TDatabasePDG;
+class TLorentzVector;
 class ExRootTreeBranch;
 class DelphesFactory;
+class Candidate;
 
 class DelphesHepMC3Reader
@@ -61 +63 @@
 
 private:
-  void AnalyzeParticle(DelphesFactory *factory,
-    TObjArray *allParticleOutputArray,
+  void AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate = 0);
+
+  void AnalyzeParticle(DelphesFactory *factory);
+
+  void FinalizeParticles(TObjArray *allParticleOutputArray,
     TObjArray *stableParticleOutputArray,
     TObjArray *partonOutputArray);
-
-  void FinalizeParticles(TObjArray *allParticleOutputArray);
 
   FILE *fInputFile;
@@ -79 +82 @@
   double fMomentumCoefficient, fPositionCoefficient;
 
-  std::vector<double> fWeight;
+  std::vector<double> fWeights;
 
   double fCrossSection, fCrossSectionError;
@@ -92 +95 @@
   double fPx, fPy, fPz, fE, fMass;
 
-  int fM1, fM2;
+  std::vector<std::pair<TLorentzVector *, TObjArray *> > fVertices;
+  std::vector<int> fParticles;
 
-  std::map<int, int> fVertexMap;
+  std::map<int, int> fInVertexMap;
+  std::map<int, int> fOutVertexMap;
+
+  std::map<int, std::pair<int, int> > fMotherMap;
   std::map<int, std::pair<int, int> > fDaughterMap;
 };

modules/Delphes.cc

@@ -61 +61 @@
   fFactory(0)
 {
-  TFolder *folder = new TFolder(name, "");
+  TFolder *folder;
+
   fFactory = new DelphesFactory("ObjectFactory");
+
+  folder = new TFolder(name, "");
 
   SetName(name);
@@ -80 +83 @@
   if(folder)
   {
+    gROOT->GetListOfBrowsables()->Remove(folder);
     folder->Clear();
     delete folder;

modules/ParticlePropagator.cc

@@ -126 +126 @@
   Double_t px, py, pz, pt, pt2, e, q;
   Double_t x, y, z, t, r;
-  Double_t x_c, y_c, r_c, phi_c, phi_0;
-  Double_t x_t, y_t, z_t, r_t;
-  Double_t t_z, t_r;
-  Double_t discr;
+  Double_t x_c, y_c, r_c, phi_0;
+  Double_t x_t, y_t, z_t, r_t, phi_t;
+  Double_t t_r, t_z;
+  Double_t tmp;
   Double_t gammam, omega;
   Double_t xd, yd, zd;
   Double_t l, d0, dz, ctgTheta, alpha;
   Double_t bsx, bsy, bsz;
-  Double_t rxp, rdp, t_R;
-  Double_t td, pio, phid, sign_pz, vz;
+  Double_t td, pio, phid, vz;
 
   const Double_t c_light = 2.99792458E8;
 
   if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
+  {
     beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  }
   else
   {
@@ -161 +162 @@
     particlePosition = particle->Position;
     particleMomentum = particle->Momentum;
-
-    // Constants
 
     x = particlePosition.X() * 1.0E-3;
@@ -208 +207 @@
     else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
     {
-
-      rxp = x*py - y*px;
-      rdp = x*px + y*py;
-
-      discr = fRadius*fRadius*pt*pt - rxp*rxp;
-
-      t_R = e * (sqrt(discr) - rdp) / (c_light * pt * pt);
-      t_z = e * (TMath::Sign(fHalfLengthMax, pz) - z) / ( c_light * pz);
-
-      t = TMath::Min(t_R, t_z);
-
-      x_t = x + px*t*c_light/e;
-      y_t = y + py*t*c_light/e;
-      z_t = z + pz*t*c_light/e;
-      r_t = TMath::Hypot(x_t, y_t);
-
-      l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
+      // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
+      tmp = px * y - py * x;
+      t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
+
+      t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
+
+      t = TMath::Min(t_r, t_z);
+
+      x_t = x + px * t;
+      y_t = y + py * t;
+      z_t = z + pz * t;
+
+      l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
 
       mother = candidate;
-      candidate = static_cast<Candidate*>(candidate->Clone());
+      candidate = static_cast<Candidate *>(candidate->Clone());
 
       candidate->InitialPosition = particlePosition;
-      candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
-      candidate->L = l*1.0E3;
+      candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
+      candidate->L = l * 1.0E3;
 
       candidate->Momentum = particleMomentum;
@@ -260 +255 @@
     {
 
-      // 1.  initial transverse momentum p_{T0}: Part->pt
-      //     initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
-      //     relativistic gamma: gamma = E/mc^2; gammam = gamma * m
-      //     gyration frequency omega = q/(gamma m) fBz
-      //     helix radius r = p_{T0} / (omega gamma m)
-
-      gammam = e*1.0E9 / (c_light*c_light);      // gammam in [eV/c^2]
-      omega = q * fBz / (gammam);                // omega is here in [89875518/s]
-      r = pt / (q * fBz) * 1.0E9/c_light;        // in [m]
+      // 1. initial transverse momentum p_{T0}: Part->pt
+      //    initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
+      //    relativistic gamma: gamma = E / mc^2; gammam = gamma * m
+      //    gyration frequency omega = q * Bz / (gammam)
+      //    helix radius r = p_{T0} / (omega * gammam)
+
+      gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
+      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]
 
       // 2. helix axis coordinates
-      x_c = x + r*TMath::Sin(phi_0);
-      y_c = y - r*TMath::Cos(phi_0);
+      x_c = x + r * TMath::Sin(phi_0);
+      y_c = y - r * TMath::Cos(phi_0);
       r_c = TMath::Hypot(x_c, y_c);
-      phi_c = TMath::ATan(y_c/x_c);
-      if(x_c < 0.0) phi_c -= TMath::Sign(1., phi_c)*TMath::Pi();
-
-      //Find the time of closest approach
-      td = (phi_0 - TMath::ATan(-x_c/y_c))/omega;
-
-      //Remove all the modulo pi that might have come from the atan
-      pio = fabs(TMath::Pi()/omega);
-      while(fabs(td) > 0.5*pio)
-      {
-        td -= TMath::Sign(1., td)*pio;
-      }
-
-      //Compute the coordinate of closed approach to z axis
-      //if wants wtr beamline need to be changedto re-center with a traslation of the z axis
-      phid = phi_0 - omega*td;
-      xd = x_c - r*TMath::Sin(phid);
-      yd = y_c + r*TMath::Cos(phid);
-      zd = z + c_light*(pz/e)*td;
-
-      //Compute momentum at closest approach (perigee??)
-      px = pt*TMath::Cos(phid);
-      py = pt*TMath::Sin(phid);
+
+      // time of closest approach
+      td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
+
+      // remove all the modulo pi that might have come from the atan
+      pio = TMath::Abs(TMath::Pi() / omega);
+      while(TMath::Abs(td) > 0.5 * pio)
+      {
+        td -= TMath::Sign(1.0, td) * pio;
+      }
+
+      vz = pz * c_light / e;
+
+      // calculate coordinates of closest approach to z axis
+      phid = phi_0 - omega * td;
+      xd = x_c - r * TMath::Sin(phid);
+      yd = y_c + r * TMath::Cos(phid);
+      zd = z + vz * td;
+
+      // momentum at closest approach
+      px = pt * TMath::Cos(phid);
+      py = pt * TMath::Sin(phid);
 
       particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
@@ -305 +299 @@
       d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
       dz = zd - bsz;
-      ctgTheta  = 1.0 / TMath::Tan (particleMomentum.Theta());
+      ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
 
       // 3. time evaluation t = TMath::Min(t_r, t_z)
       //    t_r : time to exit from the sides
       //    t_z : time to exit from the front or the back
-      t = 0;
-      t_z = 0;
-      sign_pz = (pz > 0.0) ? 1 : -1;
-      if(pz == 0.0) t_z = 1.0E99;
-      else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
-
-      if(r_c + TMath::Abs(r)  < fRadius)   // helix does not cross the cylinder sides
-      {
+      t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
+
+      if(r_c + TMath::Abs(r) < fRadius)
+      {
+        // helix does not cross the cylinder sides
         t = t_z;
       }
       else
       {
-        alpha = -(fRadius*fRadius - r*r - r_c*r_c)/(2*fabs(r)*r_c);
-        alpha = fabs(TMath::ACos(alpha));
-        t_r = td + alpha/fabs(omega);
+        alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
+        t_r = td + TMath::Abs(alpha / omega);
 
         t = TMath::Min(t_r, t_z);
       }
 
-      x_t = x_c - r*TMath::Sin(phi_0 - omega*t);
-      y_t = y_c + r*TMath::Cos(phi_0 - omega*t);
-      z_t = z + c_light*t*pz/e;
-      r_t =  TMath::Hypot(x_t, y_t);
-
-      // compute path length for an helix
-      vz = pz*1.0E9 / c_light / gammam;
-      //lenght of the path from production to tracker
-      l = t * TMath::Sqrt(vz*vz + r*r*omega*omega);
+      // 4. position in terms of x(t), y(t), z(t)
+      phi_t = phi_0 - omega * t;
+      x_t = x_c - r * TMath::Sin(phi_t);
+      y_t = y_c + r * TMath::Cos(phi_t);
+      z_t = z + vz * t;
+      r_t = TMath::Hypot(x_t, y_t);
+
+      // lenght of the path from production to tracker
+      l = t * TMath::Hypot(vz, r * omega);
 
       if(r_t > 0.0)