Add time for decay vertices

src/erhic/TreeToHepMC.cxx

@@ -213,8 +213,10 @@ Long64_t TreeToHepMC(const std::string& inputFileName,
     mcTree->GetEntry(i);
 
     // Construct new empty HepMC3 event and fill it.
-    // Using GeV and cm (!)
-    GenEvent hepmc3evt( HepMC3::Units::GEV, HepMC3::Units::CM );
+    // Using GeV and mm (!) as that's what afterburner for beam effects expect
+    // It has been confirmed that BeAGLE uses MM
+    // Need to check for other models as well
+    GenEvent hepmc3evt( HepMC3::Units::GEV, HepMC3::Units::MM );
     hepmc3evt.set_event_number(i);
     hepmc3evt.weights().clear();
     hepmc3evt.weights().push_back(1.0);
@@ -781,7 +783,28 @@ Long64_t TreeToHepMC(const std::string& inputFileName,
       // update prod vertex?
       if ( momindex > 1){
         auto vnew = inParticle->GetVertex();
-        momend->set_position( FourVector( vnew.x(), vnew.y(), vnew.z(), 0));
+
+        // recalcuate vertex time using: dt = (SV-PV).Mag()/(pMother/eMother)
+        // where PV is the mother vertex
+        // add mother decay time recursively to account for cascading decays
+
+        int daughter_index = inParticle->GetIndex();
+        int mother_index = momindex;
+        double dt = 0;
+        while (mother_index > 1 && mother_index != beagle_final_index){
+          const Particle* daughter = inEvent->GetTrack(daughter_index-1);
+          const Particle* mother = inEvent->GetTrack(mother_index-1);
+          double decay_length = (daughter->GetVertex()-mother->GetVertex()).Mag();
+	  if(decay_length>1e-3){
+	    double pMother = mother->GetP();
+	    double eMother = mother->GetE();
+	    dt += decay_length * eMother / pMother; // E/p = 1/beta; this order of operations is numerically more stable
+          }
+
+          daughter_index = mother_index;
+          mother_index = mother->GetParentIndex();
+        }
+        momend->set_position( FourVector( vnew.x(), vnew.y(), vnew.z(), dt));
       }
       // file->write_event(hepmc3evt);
     }