diff --git a/Generators/TruthIO/TruthIO/DumpMC.h b/Generators/TruthIO/TruthIO/DumpMC.h
index 8587f83806880008b5f8a93b385f2f6335fea8ee..e9ab377bdba4dffd1ac4f9e15b9c3bd6220a01e7 100644
--- a/Generators/TruthIO/TruthIO/DumpMC.h
+++ b/Generators/TruthIO/TruthIO/DumpMC.h
@@ -20,6 +20,7 @@ public:
bool m_VerboseOutput;
bool m_DeepCopy;
bool m_EtaPhi;
+ bool m_PrintQuasiStableParticles;
};
#endif
diff --git a/Generators/TruthIO/src/DumpMC.cxx b/Generators/TruthIO/src/DumpMC.cxx
index 46fc84b3fd6102dd3dd22f8751bad1505f8974d8..0a5f215389339caed682bc0d50a51fafaf05d489 100644
--- a/Generators/TruthIO/src/DumpMC.cxx
+++ b/Generators/TruthIO/src/DumpMC.cxx
@@ -7,6 +7,7 @@
using namespace Gaudi::Units;
using namespace std;
+#include <CLHEP/Vector/LorentzVector.h>
DumpMC::DumpMC(const string& name, ISvcLocator* pSvcLocator)
: GenBase(name, pSvcLocator)
@@ -15,6 +16,7 @@ DumpMC::DumpMC(const string& name, ISvcLocator* pSvcLocator)
declareProperty("VerboseOutput", m_VerboseOutput=true);
declareProperty("DeepCopy", m_DeepCopy=false);
declareProperty("EtaPhi", m_EtaPhi=false);
+ declareProperty("PrintQuasiStableParticles", m_PrintQuasiStableParticles=false);
}
@@ -139,6 +141,30 @@ StatusCode DumpMC::execute() {
cout << " eta = " << rapid<< " Phi = " << phi << " Et = " <<et/GeV << " Status= " << p_stat << " PDG ID= "<< p_id << endl;
}
}
+ if(m_PrintQuasiStableParticles) {
+ const HepMC::GenEvent* genEvt = (*itr);
+ for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) {
+ int p_stat=(*pitr)->status();
+ if(p_stat==2 && (*pitr)->production_vertex() && (*pitr)->end_vertex()) {
+ const auto& prodVtx = (*pitr)->production_vertex()->position();
+ const auto& endVtx = (*pitr)->end_vertex()->position();
+ const CLHEP::HepLorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
+ const CLHEP::HepLorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
+
+ CLHEP::HepLorentzVector dist4D(lv1);
+ dist4D-=lv0;
+ CLHEP::Hep3Vector dist3D=dist4D.vect();
+ if(dist3D.mag()>1*Gaudi::Units::mm) {
+ const auto& GenMom = (*pitr)->momentum();
+ const CLHEP::HepLorentzVector mom( GenMom.x(), GenMom.y(), GenMom.z(), GenMom.t() );
+ ATH_MSG_INFO("Quasi stable particle "<<**pitr);
+ ATH_MSG_INFO(" Prod V:"<<*((*pitr)->production_vertex ()));
+ ATH_MSG_INFO(" Decay V:"<<*((*pitr)->end_vertex ()));
+ ATH_MSG_INFO(" gamma(Momentum)="<<mom.gamma()<<" gamma(Vertices)="<<dist4D.gamma());
+ }
+ }
+ }
+ }
}
}
return StatusCode::SUCCESS;
diff --git a/Simulation/ISF/ISF_Core/ISF_Services/src/InputConverter.cxx b/Simulation/ISF/ISF_Core/ISF_Services/src/InputConverter.cxx
index 8cf7f05c97eaacd6ed3e0e3d574acd64b22bf46f..c36ca3dc68675ade26aea43fa0dcf13276f70245 100644
--- a/Simulation/ISF/ISF_Core/ISF_Services/src/InputConverter.cxx
+++ b/Simulation/ISF/ISF_Core/ISF_Services/src/InputConverter.cxx
@@ -255,6 +255,27 @@ ISF::InputConverter::convertParticle(HepMC::GenParticle* genPartPtr, int bcid) c
if(std::abs(e-teste)/e>0.01) {
ATH_MSG_WARNING("Difference in energy for: " << genPart<<" Morg="<<pMomentum.m()<<" Mmod="<<pMass<<" Eorg="<<e<<" Emod="<<teste);
}
+ if(genPart.status()==2 && genPart.production_vertex() && genPart.end_vertex()) { //check for possible changes of gamma for quasi stable particles
+ const auto& prodVtx = genPart.production_vertex()->position();
+ const auto& endVtx = genPart.end_vertex()->position();
+ CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
+
+ if(dist3D.mag()>1*Gaudi::Units::mm) {
+ const auto& GenMom = genPart.momentum();
+ CLHEP::HepLorentzVector mom( GenMom.x(), GenMom.y(), GenMom.z(), GenMom.t() );
+ double gamma_org=mom.gamma();
+ mom.setE(teste);
+ double gamma_new=mom.gamma();
+
+ if(std::abs(gamma_new-gamma_org)/(gamma_new+gamma_org)>0.001) {
+ ATH_MSG_WARNING("Difference in boost gamma for Quasi stable particle "<<genPart);
+ ATH_MSG_WARNING(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
+ } else {
+ ATH_MSG_VERBOSE("Quasi stable particle "<<genPart);
+ ATH_MSG_VERBOSE(" gamma(m="<<mom.m()<<")="<<gamma_org<<" gamma(m="<<pMass<<")="<<gamma_new);
+ }
+ }
+ }
}
const int pPdgId = genPart.pdg_id();
@@ -417,15 +438,25 @@ G4PrimaryParticle* ISF::InputConverter::getG4PrimaryParticle(const HepMC::GenPar
// to validate for every generator on earth...
const auto& prodVtx = genpart.production_vertex()->position();
const auto& endVtx = genpart.end_vertex()->position();
- const G4LorentzVector lv0 ( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
- const G4LorentzVector lv1 ( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
- g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
+ //const G4LorentzVector lv0 ( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
+ //const G4LorentzVector lv1 ( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
+ //Old calculation, not taken because vertex information is not sufficiently precise
+ //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
+
+ CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
+ double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
+ pmag2*=pmag2; //magnitude of particle momentum squared
+ double e2=g4particle->GetTotalEnergy(); //energy of particle
+ e2*=e2; //energy of particle squared
+ double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
+ double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
+ g4particle->SetProperTime( std::sqrt(tau2) );
if(m_quasiStableParticlesIncluded) {
ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
ATH_MSG_VERBOSE( "Will add the primary particle set on." );
ATH_MSG_VERBOSE( "Primary Particle: " << genpart );
- ATH_MSG_VERBOSE( "Number of daughters of "<<genpart.barcode()<<": " << genpart.end_vertex()->particles_out_size() );
+ ATH_MSG_VERBOSE( "Number of daughters of "<<genpart.barcode()<<": " << genpart.end_vertex()->particles_out_size()<<" at position "<<*genpart.end_vertex());
}
else {
ATH_MSG_WARNING( "Detected primary particle with end vertex." );
@@ -524,16 +555,39 @@ G4PrimaryParticle* ISF::InputConverter::getG4PrimaryParticle(const ISF::ISFParti
/// that we have to validate for every generator on earth...
const auto& prodVtx = genpart->production_vertex()->position();
const auto& endVtx = genpart->end_vertex()->position();
- const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
- const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
- g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
-
+
+ CLHEP::Hep3Vector dist3D(endVtx.x()-prodVtx.x(), endVtx.y()-prodVtx.y(), endVtx.z()-prodVtx.z());
+
+ double pmag2=g4particle->GetTotalMomentum(); //magnitude of particle momentum
+ pmag2*=pmag2; //magnitude of particle momentum squared
+ double e2=g4particle->GetTotalEnergy(); //energy of particle
+ e2*=e2; //energy of particle squared
+ double beta2=pmag2/e2; //beta^2=v^2/c^2 for particle
+
+ double tau2=dist3D.mag2()*(1/beta2-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
+
+ g4particle->SetProperTime( std::sqrt(tau2) );
+
+ if(msgLvl(MSG::VERBOSE)) {
+ const G4LorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
+ const G4LorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
+ //Old calculation, not taken because vertex information is not sufficiently precise
+ //g4particle->SetProperTime( (lv1-lv0).mag()/Gaudi::Units::c_light );
+ G4LorentzVector dist4D(lv1);
+ dist4D-=lv0;
+
+ G4LorentzVector fourmom(g4particle->GetMomentum(),g4particle->GetTotalEnergy());
+ //double tau2=dist3D.mag2()*(1/(fourmom.beta()*fourmom.beta())-1)/Gaudi::Units::c_light/Gaudi::Units::c_light;
+
+ ATH_MSG_VERBOSE( "gammaVertex="<<dist4D.gamma()<<" gammamom="<<fourmom.gamma()<<" gamma(beta)="<<1/std::sqrt(1-beta2)<<" lifetime tau="<<std::sqrt(tau2));
+ }
+
if(m_quasiStableParticlesIncluded) {
ATH_MSG_VERBOSE( "Detected primary particle with end vertex." );
ATH_MSG_VERBOSE( "Will add the primary particle set on." );
ATH_MSG_VERBOSE( "ISF Particle: " << isp );
ATH_MSG_VERBOSE( "Primary Particle: " << *genpart );
- ATH_MSG_VERBOSE( "Number of daughters of "<<genpart->barcode()<<": " << genpart->end_vertex()->particles_out_size() );
+ ATH_MSG_VERBOSE( "Number of daughters of "<<genpart->barcode()<<": " << genpart->end_vertex()->particles_out_size() << " at position "<< *(genpart->end_vertex()));
}
else {
ATH_MSG_WARNING( "Detected primary particle with end vertex. This should only be the case if" );
diff --git a/Simulation/ISF/ISF_Core/ISF_Services/src/TruthSvc.cxx b/Simulation/ISF/ISF_Core/ISF_Services/src/TruthSvc.cxx
index 23651a5f9246324f9709996cde69bdc667ff1a93..dcdff36d3ca26e1f6bf526526eb654d4d3cfccdf 100644
--- a/Simulation/ISF/ISF_Core/ISF_Services/src/TruthSvc.cxx
+++ b/Simulation/ISF/ISF_Core/ISF_Services/src/TruthSvc.cxx
@@ -452,7 +452,18 @@ HepMC::GenVertex *ISF::TruthSvc::createGenVertexFromTruthIncident( ISF::ITruthIn
}
else {
//oldVertex->suggest_barcode( vtxbcode );
- oldVertex->set_position( ti.position() );
+ const auto& old_pos=oldVertex->position();
+ const auto& new_pos=ti.position();
+ HepMC::ThreeVector diff(new_pos.x()-old_pos.x(),new_pos.y()-old_pos.y(),new_pos.z()-old_pos.z()); //complicated, but HepMC::ThreeVector and FourVector have no + or - operators
+
+ if(diff.r()>1*Gaudi::Units::mm) { //Check for a change of the vertex position by more than 1mm
+ ATH_MSG_WARNING("For particle: " << *parent);
+ ATH_MSG_WARNING(" decay vertex before QS partice sim: " << *oldVertex );
+ oldVertex->set_position( ti.position() );
+ ATH_MSG_WARNING(" decay vertex after QS partice sim: " << *oldVertex );
+ } else {
+ oldVertex->set_position( ti.position() );
+ }
oldVertex->set_id( vtxID );
oldVertex->weights() = weights;
#ifdef DEBUG_TRUTHSVC