diff --git a/GeneratorInterface/TauolaInterface/interface/TauolappInterface.h b/GeneratorInterface/TauolaInterface/interface/TauolappInterface.h
index 21a439db97c2c..513e45d534ada 100644
--- a/GeneratorInterface/TauolaInterface/interface/TauolappInterface.h
+++ b/GeneratorInterface/TauolaInterface/interface/TauolappInterface.h
@@ -1,70 +1,55 @@
 #ifndef gen_TauolaInterface_TauolappInterface_h
 #define gen_TauolaInterface_TauolappInterface_h
-
 #include "HepPDT/ParticleDataTable.hh"
-
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "FWCore/Framework/interface/ESHandle.h"
 #include "FWCore/Framework/interface/EventSetup.h"
 #include "GeneratorInterface/TauolaInterface/interface/TauolaInterfaceBase.h"
-
-namespace HepMC 
+namespace HepMC
 {
-class GenEvent;
+  class GenEvent;
 }
-
 namespace CLHEP
 {
-class HepRandomEngine;
+  class HepRandomEngine;
 }
-
 namespace gen {
   extern "C" {
     void ranmar_( float *rvec, int *lenv );
     void rmarin_( int*, int*, int* );
   }
-
-   class TauolappInterface : public TauolaInterfaceBase {
-      public:
-      
-      // ctor & dtor
-     TauolappInterface( const edm::ParameterSet& );
-     ~TauolappInterface();
-     
-      void enablePolarization()  { fPolarization = true; return; }
-      void disablePolarization() { fPolarization = false; return; }
-      void init( const edm::EventSetup& );
-      const std::vector<int>& operatesOnParticles() { return fPDGs; }
-      HepMC::GenEvent* decay( HepMC::GenEvent* );
-      void statistics() ;
-
-      void setRandomEngine(CLHEP::HepRandomEngine* v) { fRandomEngine = v; }
-      static double flat();
-      
-      private: 
-      // member function(s)
-      void decodeMDTAU( int );
-      void selectDecayByMDTAU();
-      int selectLeptonic();
-      int selectHadronic();
-      
-      //
-      static CLHEP::HepRandomEngine*           fRandomEngine;            
-      std::vector<int>                         fPDGs;
-      bool                                     fPolarization;      
-      edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable ;
-      edm::ParameterSet*                       fPSet;
-      bool                                     fIsInitialized;
-      
-      int                                      fMDTAU;
-      bool                                     fSelectDecayByEvent;
-      std::vector<int>                         fLeptonModes;
-      std::vector<int>                         fHadronModes;
-      std::vector<double>                      fScaledLeptonBrRatios;
-      std::vector<double>                      fScaledHadronBrRatios;
-      
-   };
-
+  class TauolappInterface : public TauolaInterfaceBase {
+  public:
+    // ctor & dtor
+    TauolappInterface( const edm::ParameterSet& );
+    ~TauolappInterface();
+    void enablePolarization() { fPolarization = true; return; }
+    void disablePolarization() { fPolarization = false; return; }
+    void init( const edm::EventSetup& );
+    const std::vector<int>& operatesOnParticles() { return fPDGs; }
+    HepMC::GenEvent* decay( HepMC::GenEvent* );
+    void statistics() ;
+    void setRandomEngine(CLHEP::HepRandomEngine* v) { fRandomEngine = v; }
+    static double flat();
+  private:
+    // member function(s)
+    void decodeMDTAU( int );
+    void selectDecayByMDTAU();
+    int selectLeptonic();
+    int selectHadronic();
+    //
+    static CLHEP::HepRandomEngine* fRandomEngine;
+    std::vector<int> fPDGs;
+    bool fPolarization;
+    edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable ;
+    edm::ParameterSet* fPSet;
+    bool fIsInitialized;
+    int fMDTAU;
+    bool fSelectDecayByEvent;
+    std::vector<int> fLeptonModes;
+    std::vector<int> fHadronModes;
+    std::vector<double> fScaledLeptonBrRatios;
+    std::vector<double> fScaledHadronBrRatios;
+  };
 }
-
 #endif
diff --git a/GeneratorInterface/TauolaInterface/plugins/Tauolapp/TauolappInterface.cc b/GeneratorInterface/TauolaInterface/plugins/Tauolapp/TauolappInterface.cc
index f11ed57c9692d..cb028e4e5cc37 100644
--- a/GeneratorInterface/TauolaInterface/plugins/Tauolapp/TauolappInterface.cc
+++ b/GeneratorInterface/TauolaInterface/plugins/Tauolapp/TauolappInterface.cc
@@ -1,48 +1,33 @@
 /* this is the code for the new Tauola++ */
-
 #include <iostream>
-
 #include "GeneratorInterface/TauolaInterface/interface/TauolappInterface.h"
-
 #include "Tauola/Tauola.h"
 #include "Tauola/TauolaHepMCEvent.h"
 #include "Tauola/Log.h"
-
 #include "FWCore/ServiceRegistry/interface/Service.h"
 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
 #include "FWCore/Utilities/interface/Exception.h"
-
 #include "CLHEP/Random/RandomEngine.h"
-
 #include "HepMC/GenEvent.h"
 #include "HepMC/IO_HEPEVT.h"
 #include "HepMC/HEPEVT_Wrapper.h"
-
 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
-
-
 using namespace gen;
 using namespace edm;
 using namespace std;
-
-CLHEP::HepRandomEngine*  TauolappInterface::fRandomEngine = nullptr;
-
+CLHEP::HepRandomEngine* TauolappInterface::fRandomEngine = nullptr;
 extern "C" {
-
   void gen::ranmar_( float *rvec, int *lenv )
   {
     for(int i = 0; i < *lenv; i++)
       *rvec++ = TauolappInterface::flat();
     return;
   }
-  
   void gen::rmarin_( int*, int*, int* )
   {
     return;
   }
-  
 }
-
 TauolappInterface::TauolappInterface( const edm::ParameterSet& pset):
   fPolarization(false),
   fPSet(0),
@@ -53,257 +38,188 @@ TauolappInterface::TauolappInterface( const edm::ParameterSet& pset):
   if ( fPSet != 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to override Tauola an existing ParameterSet\n" << std::endl;
   fPSet = new ParameterSet(pset);
 }
-
 TauolappInterface::~TauolappInterface(){if ( fPSet != 0 ) delete fPSet;}
-
 void TauolappInterface::init( const edm::EventSetup& es ){
   if ( fIsInitialized ) return; // do init only once
-  if ( fPSet == 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n" << std::endl;   
-  
+  if ( fPSet == 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n" << std::endl;
   fIsInitialized = true;
-      
-   es.getData( fPDGTable ) ;
-
-   Tauolapp::Tauola::setDecayingParticle(15);
-   // --> ??? Tauola::setRadiation(false);
-
-   // polarization switch 
-   //
-   // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
-   fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
-   
-   // read tau decay mode switches
-   //
-   ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
-   
-   fMDTAU = cards.getParameter< int >( "mdtau" );
-
-   if ( fMDTAU == 0 || fMDTAU == 1 )
-   {
+  es.getData( fPDGTable ) ;
+  Tauolapp::Tauola::setDecayingParticle(15);
+  // --> ??? Tauola::setRadiation(false);
+  // polarization switch
+  //
+  // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
+  fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
+  // read tau decay mode switches
+  //
+  ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
+  fMDTAU = cards.getParameter< int >( "mdtau" );
+  if ( fMDTAU == 0 || fMDTAU == 1 )
+    {
       Tauolapp::Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
-   }
-
-   Tauolapp::Tauola::setTauLifetime(0.0);
-   Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
-
-   // some more options, copied over from an example 
-   // - maybe will use later...
-   //
-   //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off. 
-   //
-
-//
-//   const HepPDT::ParticleData* 
-//         PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauola::getDecayingParticle()) )) ;
-//   double lifetime = PData->lifetime().value();
-//   Tauola::setTauLifetime( lifetime );
-
-   fPDGs.push_back( Tauolapp::Tauola::getDecayingParticle() );
-
-   Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);         
-   Tauolapp::Tauola::initialize();
-
-   Tauolapp::Tauola::spin_correlation.setAll(fPolarization);// Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
-
-   if (fPSet->exists("parameterSets")){
-     std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
-     for (unsigned int ip=0; ip<par.size(); ++ip ){
-       std::string curSet = par[ip];
-       
-       if(curSet=="setNewCurrents") Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
-       if(curSet=="spinCorrelationSetAll") Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
-       if(curSet=="spinCorrelationGAMMA") Tauolapp::Tauola::spin_correlation.GAMMA=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationZ0") Tauolapp::Tauola::spin_correlation.Z0=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationHIGGS") Tauolapp::Tauola::spin_correlation.HIGGS=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationHIGGSH") Tauolapp::Tauola::spin_correlation.HIGGS_H=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationHIGGSA") Tauolapp::Tauola::spin_correlation.HIGGS_A=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationHIGGSPLUS") Tauolapp::Tauola::spin_correlation.HIGGS_PLUS=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationHIGGSMINUS") Tauolapp::Tauola::spin_correlation.HIGGS_MINUS=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationWPLUS") Tauolapp::Tauola::spin_correlation.W_PLUS=fPSet->getParameter<bool>(curSet);
-       if(curSet=="spinCorrelationWMINUS") Tauolapp::Tauola::spin_correlation.W_MINUS=fPSet->getParameter<bool>(curSet);
-       
-       if(curSet=="setHiggsScalarPseudoscalarPDG") Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
-       if(curSet=="setHiggsScalarPseudoscalarMixingAngle") Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
-       
-       if(curSet=="setRadiation") Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
-       if(curSet=="setRadiationCutOff") Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
-       
-       if(curSet=="setEtaK0sPi"){
-	 std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
-	 if(vpar.size()==3) Tauolapp::Tauola::setEtaK0sPi(vpar[0],vpar[1],vpar[2]);
-	 else {std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;}
-       }
-       
-       if(curSet=="setTaukle"){
-	 std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
-	 if(vpar.size()==4) Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
-	 else {std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;}
-       }
-       
-       if(curSet=="setTauBr"){
-	 edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
-	 std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
-	 std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
-	 if(vJAK.size() == vBR.size()){
-	   for(unsigned int i=0;i<vJAK.size();i++) Tauolapp::Tauola::setTauBr(vJAK[i],vBR[i]);
-	 }
-	 else {std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;}
-       }
-     }
-   }
-
-   // override decay modes if needs be
-   //
-   // we have to do it AFTER init because otherwises branching ratios are NOT filled in
-   //
-   if ( fMDTAU != 0 && fMDTAU != 1 )
-   {
+    }
+  Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
+  // some more options, copied over from an example
+  // Default values
+  //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
+  const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauolapp::Tauola::getDecayingParticle()) )) ;
+  double lifetime = PData->lifetime().value();
+  Tauolapp::Tauola::setTauLifetime( lifetime );
+  //std::cout << "Setting lifetime: ctau=" << lifetime << "m or tau=" << lifetime/(299792458.0*1000.0) << "s" << std::endl;
+  fPDGs.push_back( Tauolapp::Tauola::getDecayingParticle() );
+  Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);
+  Tauolapp::Tauola::initialize();
+  Tauolapp::Tauola::spin_correlation.setAll(fPolarization);// Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
+  if (fPSet->exists("parameterSets")){
+    std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
+    for (unsigned int ip=0; ip<par.size(); ++ip ){
+      std::string curSet = par[ip];
+      if(curSet=="setNewCurrents") Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
+      if(curSet=="spinCorrelationSetAll") Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
+      if(curSet=="spinCorrelationGAMMA") Tauolapp::Tauola::spin_correlation.GAMMA=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationZ0") Tauolapp::Tauola::spin_correlation.Z0=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationHIGGS") Tauolapp::Tauola::spin_correlation.HIGGS=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationHIGGSH") Tauolapp::Tauola::spin_correlation.HIGGS_H=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationHIGGSA") Tauolapp::Tauola::spin_correlation.HIGGS_A=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationHIGGSPLUS") Tauolapp::Tauola::spin_correlation.HIGGS_PLUS=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationHIGGSMINUS") Tauolapp::Tauola::spin_correlation.HIGGS_MINUS=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationWPLUS") Tauolapp::Tauola::spin_correlation.W_PLUS=fPSet->getParameter<bool>(curSet);
+      if(curSet=="spinCorrelationWMINUS") Tauolapp::Tauola::spin_correlation.W_MINUS=fPSet->getParameter<bool>(curSet);
+      if(curSet=="setHiggsScalarPseudoscalarPDG") Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
+      if(curSet=="setHiggsScalarPseudoscalarMixingAngle") Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
+      if(curSet=="setRadiation") Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
+      if(curSet=="setRadiationCutOff") Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
+      if(curSet=="setEtaK0sPi"){
+	std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
+	if(vpar.size()==3) Tauolapp::Tauola::setEtaK0sPi(vpar[0],vpar[1],vpar[2]);
+	else {std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;}
+      }
+      if(curSet=="setTaukle"){
+	std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
+	if(vpar.size()==4) Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
+	else {std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;}
+      }
+      if(curSet=="setTauBr"){
+	edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
+	std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
+	std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
+	if(vJAK.size() == vBR.size()){
+	  for(unsigned int i=0;i<vJAK.size();i++) Tauolapp::Tauola::setTauBr(vJAK[i],vBR[i]);
+	}
+	else {std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;}
+      }
+    }
+  }
+  // override decay modes if needs be
+  //
+  // we have to do it AFTER init because otherwises branching ratios are NOT filled in
+  //
+  if ( fMDTAU != 0 && fMDTAU != 1 )
+    {
       decodeMDTAU( fMDTAU );
-   }
-
-   Tauolapp::Log::LogWarning(false);
-   
-   return;
+    }
+  Tauolapp::Log::LogWarning(false);
+  return;
 }
-
 double TauolappInterface::flat()
 {
-   if ( !fRandomEngine ) {
-     throw cms::Exception("LogicError")
-       << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
-       << "This might mean that the code was modified to generate a random number outside the\n"
-       << "event and beginLuminosityBlock methods, which is not allowed.\n";
-   }
-   return fRandomEngine->flat();
+  if ( !fRandomEngine ) {
+    throw cms::Exception("LogicError")
+      << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
+      << "This might mean that the code was modified to generate a random number outside the\n"
+      << "event and beginLuminosityBlock methods, which is not allowed.\n";
+  }
+  return fRandomEngine->flat();
 }
-
 HepMC::GenEvent* TauolappInterface::decay( HepMC::GenEvent* evt ){
-   if ( !fIsInitialized ) return evt;
-   Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);   // rest tauola++ random number incase other modules use tauola++
-   int NPartBefore = evt->particles_size();
-   int NVtxBefore  = evt->vertices_size();
-   
-   // what do we do if Hep::GenEvent size is larger than 10K ???
-   // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
-   // and in case of CMS, it's only up to 4K !!!
-   //
-   // if ( NPartBefore > 10000 ) return evt;
-   //
-   
-   // override decay mode if needs be
-   if ( fSelectDecayByEvent )
-   {
+  if ( !fIsInitialized ) return evt;
+  Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat); // rest tauola++ random number incase other modules use tauola++
+  int NPartBefore = evt->particles_size();
+  int NVtxBefore = evt->vertices_size();
+  // what do we do if Hep::GenEvent size is larger than 10K ???
+  // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
+  // and in case of CMS, it's only up to 4K !!!
+  //
+  // if ( NPartBefore > 10000 ) return evt;
+  //
+  // override decay mode if needs be
+  if ( fSelectDecayByEvent )
+    {
       selectDecayByMDTAU();
-   }
-   
-    //construct tmp TAUOLA event
-    //
-    auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
-   
-    // another option: if one lets Pythia or another master gen to decay taus, 
-    //                 we have to undecay them first
-    // t_event->undecayTaus();
-    
-    // run Tauola on the tmp event - HepMC::GenEvernt will be MODIFIED !!!
-    //
-    t_event->decayTaus();
-    
-    // delet tmp Tauola event
-    //
-    delete t_event; 
-    
-    // do we also need to apply the lifetime and vtx position shift ??? 
-    // (see TauolappInterface, for example)
-    //
-    // NOTE: the procedure ASSYMES that vertex barcoding is COUNTIUOUS/SEQUENTIAL,
-    // and that the abs(barcode) corresponds to vertex "plain indexing"
-    //
-    for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ )
-    {
-       HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
-       HepMC::GenParticle* GenPart = *(GenVtx->particles_in_const_begin());
-       HepMC::GenVertex* ProdVtx = GenPart->production_vertex();
-       HepMC::FourVector PMom = GenPart->momentum();
-       double mass = GenPart->generated_mass();
-       const HepPDT::ParticleData* 
-             PData = fPDGTable->particle(HepPDT::ParticleID(abs(GenPart->pdg_id()))) ;
-       double lifetime = PData->lifetime().value();
-       float prob = 0.;
-       int length=1;
-       ranmar_(&prob,&length);
-       double ct = -lifetime * std::log(prob);
-       double VxDec = GenVtx->position().x();
-       VxDec += ct * (PMom.px()/mass);
-       VxDec += ProdVtx->position().x();
-       double VyDec = GenVtx->position().y();
-       VyDec += ct * (PMom.py()/mass);
-       VyDec += ProdVtx->position().y();
-       double VzDec = GenVtx->position().z();
-       VzDec += ct * (PMom.pz()/mass);
-       VzDec += ProdVtx->position().z();
-       double VtDec = GenVtx->position().t();
-       VtDec += ct * (PMom.e()/mass);
-       VtDec += ProdVtx->position().t();
-       GenVtx->set_position( HepMC::FourVector(VxDec,VyDec,VzDec,VtDec) ); 
-       //
-       // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
-       // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
-       // thus we have to take a 2-step procedure
-       //
-       std::vector<int> BCodes;
-       BCodes.clear();
-       for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
-                                               pitr != GenVtx->particles_end(HepMC::children); ++pitr) 
-       {
+    }
+  //construct tmp TAUOLA event
+  //
+  auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
+  // another option: if one lets Pythia or another master gen to decay taus,
+  // we have to undecay them first
+  // t_event->undecayTaus();
+  // run Tauola on the tmp event - HepMC::GenEvernt will be MODIFIED !!!
+  //
+  t_event->decayTaus();
+  // delet tmp Tauola event
+  //
+  delete t_event;
+  // do we also need to apply the lifetime and vtx position shift ???
+  // (see TauolappInterface, for example)
+  //
+  // NOTE: the procedure ASSYMES that vertex barcoding is COUNTIUOUS/SEQUENTIAL,
+  // and that the abs(barcode) corresponds to vertex "plain indexing"
+  //
+  for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ )
+    {
+      HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
+      //
+      // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
+      // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
+      // thus we have to take a 2-step procedure
+      //
+      std::vector<int> BCodes;
+      BCodes.clear();
+      for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
+	   pitr != GenVtx->particles_end(HepMC::children); ++pitr)
+	{
 	  if ( (*pitr)->barcode() > 10000 )
-	  {
-	     BCodes.push_back( (*pitr)->barcode() );
-	  }
-       }
-       if ( BCodes.size() > 0 )
-       {
-          for ( size_t ibc=0; ibc<BCodes.size(); ibc++ )
-	  {
-	     HepMC::GenParticle* p1 = evt->barcode_to_particle( BCodes[ibc] );
-	     int nbc = p1->barcode() - 10000 + NPartBefore;
-             p1->suggest_barcode( nbc );
-	  }
-       }             
-    }
-        
-    return evt;
-      
+	    {
+	      BCodes.push_back( (*pitr)->barcode() );
+	    }
+	}
+      if ( BCodes.size() > 0 )
+	{
+	  for ( size_t ibc=0; ibc<BCodes.size(); ibc++ )
+	    {
+	      HepMC::GenParticle* p1 = evt->barcode_to_particle( BCodes[ibc] );
+	      int nbc = p1->barcode() - 10000 + NPartBefore;
+	      p1->suggest_barcode( nbc );
+	    }
+	}
+    }
+  return evt;
 }
-
 void TauolappInterface::statistics()
 {
-   return;
+  return;
 }
-
 void TauolappInterface::decodeMDTAU( int mdtau )
 {
-
-   // Note-1:
-   // I have to hack the common block directly because set<...>DecayMode(...)
-   // only changes it in the Tauola++ instance but does NOT passes it over
-   // to the Fortran core - this it does only one, via initialize() stuff...
-   //
-   // So I'll do both ways of settings, just for consistency...
-   // but I probably need to communicate it to the Tauola(++) team...
-   //
-   
-   // Note-2: 
-   // originally, the 1xx settings are meant for tau's from hard event,
-   // and the 2xx settings are for any tau in the event record;
-   //
-   // later one, we'll have to take this into account...
-   // but first I'll have to sort out what happens in the 1xx case
-   // to tau's coming outside of hard event (if any in the record)   
-   //
-   
-   if ( mdtau == 101 || mdtau == 201 )
-   {
+  // Note-1:
+  // I have to hack the common block directly because set<...>DecayMode(...)
+  // only changes it in the Tauola++ instance but does NOT passes it over
+  // to the Fortran core - this it does only one, via initialize() stuff...
+  //
+  // So I'll do both ways of settings, just for consistency...
+  // but I probably need to communicate it to the Tauola(++) team...
+  //
+  // Note-2:
+  // originally, the 1xx settings are meant for tau's from hard event,
+  // and the 2xx settings are for any tau in the event record;
+  //
+  // later one, we'll have to take this into account...
+  // but first I'll have to sort out what happens in the 1xx case
+  // to tau's coming outside of hard event (if any in the record)
+  //
+  if ( mdtau == 101 || mdtau == 201 )
+    {
       // override with electron mode for both tau's
       //
       Tauolapp::jaki_.jak1 = 1;
@@ -311,10 +227,9 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
       return;
-   }
-   
-   if ( mdtau == 102 || mdtau == 202 )
-   {
+    }
+  if ( mdtau == 102 || mdtau == 202 )
+    {
       // override with muon mode for both tau's
       //
       Tauolapp::jaki_.jak1 = 2;
@@ -322,11 +237,10 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
       return;
-   }
-
-   if ( mdtau == 111 || mdtau == 211 )
-   {
-      // override with electron mode for 1st tau 
+    }
+  if ( mdtau == 111 || mdtau == 211 )
+    {
+      // override with electron mode for 1st tau
       // and any mode for 2nd tau
       //
       Tauolapp::jaki_.jak1 = 1;
@@ -334,11 +248,10 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
-   }
-
-   if ( mdtau == 112 || mdtau == 212 )
-   {
-      // override with muon mode for the 1st tau 
+    }
+  if ( mdtau == 112 || mdtau == 212 )
+    {
+      // override with muon mode for the 1st tau
       // and any mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 2;
@@ -346,11 +259,10 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
-   }
-   
-   if ( mdtau == 121 || mdtau == 221 )
-   {
-      // override with any mode for the 1st tau 
+    }
+  if ( mdtau == 121 || mdtau == 221 )
+    {
+      // override with any mode for the 1st tau
       // and electron mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
@@ -358,11 +270,10 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
       return;
-   }
-   
-   if ( mdtau == 122 || mdtau == 222 )
-   {
-      // override with any mode for the 1st tau 
+    }
+  if ( mdtau == 122 || mdtau == 222 )
+    {
+      // override with any mode for the 1st tau
       // and muon mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
@@ -370,22 +281,20 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
       return;
-   }
-
-   if ( mdtau == 140 || mdtau == 240 )
-   {
-      // override with pi+/- nutau mode for both tau's 
+    }
+  if ( mdtau == 140 || mdtau == 240 )
+    {
+      // override with pi+/- nutau mode for both tau's
       //
       Tauolapp::jaki_.jak1 = 3;
       Tauolapp::jaki_.jak2 = 3;
       Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
       return;
-   }
-
-   if ( mdtau == 141 || mdtau == 241 )
-   {
-      // override with pi+/- nutau mode for the 1st tau 
+    }
+  if ( mdtau == 141 || mdtau == 241 )
+    {
+      // override with pi+/- nutau mode for the 1st tau
       // and any mode for the 2nd tau
       //
       Tauolapp::jaki_.jak1 = 3;
@@ -393,11 +302,10 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
       return;
-   }
-
-   if ( mdtau == 142 || mdtau == 242 )
-   {
-      // override with any mode for the 1st tau 
+    }
+  if ( mdtau == 142 || mdtau == 242 )
+    {
+      // override with any mode for the 1st tau
       // and pi+/- nutau mode for 2nd tau
       //
       Tauolapp::jaki_.jak1 = 0;
@@ -405,58 +313,47 @@ void TauolappInterface::decodeMDTAU( int mdtau )
       Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
       Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
       return;
-   }
-   
-   // OK, we come here for semi-inclusive modes
-   //
-   
-   // First of all, leptons and hadron modes sums
-   //
-   // re-scale branching ratios, just in case...
-   //
-   double sumBra = 0;
-   
-   // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
-   // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
-   //
-   
-   for ( int i=0; i<22; i++ )
-   {
+    }
+  // OK, we come here for semi-inclusive modes
+  //
+  // First of all, leptons and hadron modes sums
+  //
+  // re-scale branching ratios, just in case...
+  //
+  double sumBra = 0;
+  // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
+  // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
+  //
+  for ( int i=0; i<22; i++ )
+    {
       sumBra += Tauolapp::taubra_.gamprt[i];
-   }
-   if ( sumBra == 0. ) return ; // perhaps need to throw ?
-   for ( int i=0; i<22; i++ )
-   {
+    }
+  if ( sumBra == 0. ) return ; // perhaps need to throw ?
+  for ( int i=0; i<22; i++ )
+    {
       double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
-      Tauolapp::Tauola::setTauBr( i+1, newBra ); 
-   }
-   sumBra = 1.0;
-   
-   double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
-   double sumHadronBra = sumBra - sumLeptonBra;
-   
-   for ( int i=0; i<2; i++ )
-   {
+      Tauolapp::Tauola::setTauBr( i+1, newBra );
+    }
+  sumBra = 1.0;
+  double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
+  double sumHadronBra = sumBra - sumLeptonBra;
+  for ( int i=0; i<2; i++ )
+    {
       fLeptonModes.push_back( i+1 );
-      fScaledLeptonBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumLeptonBra) );  
-   }
-   for ( int i=2; i<22; i++ )
-   {
+      fScaledLeptonBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumLeptonBra) );
+    }
+  for ( int i=2; i<22; i++ )
+    {
       fHadronModes.push_back( i+1 );
-      fScaledHadronBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumHadronBra) ); 
-   }
-
-   fSelectDecayByEvent = true;
-   return;
-      
+      fScaledHadronBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumHadronBra) );
+    }
+  fSelectDecayByEvent = true;
+  return;
 }
-
 void TauolappInterface::selectDecayByMDTAU()
 {
-
-   
-   if ( fMDTAU == 100 || fMDTAU == 200 )
-   {
+  if ( fMDTAU == 100 || fMDTAU == 200 )
+    {
       int mode = selectLeptonic();
       Tauolapp::jaki_.jak1 = mode;
       Tauolapp::Tauola::setSameParticleDecayMode( mode );
@@ -464,167 +361,140 @@ void TauolappInterface::selectDecayByMDTAU()
       Tauolapp::jaki_.jak2 = mode;
       Tauolapp::Tauola::setOppositeParticleDecayMode( mode );
       return ;
-   }
-   
-   int modeL = selectLeptonic();
-   int modeH = selectHadronic();
-   
-   if ( fMDTAU == 110 || fMDTAU == 210 )
-   {
+    }
+  int modeL = selectLeptonic();
+  int modeH = selectHadronic();
+  if ( fMDTAU == 110 || fMDTAU == 210 )
+    {
       Tauolapp::jaki_.jak1 = modeL;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( modeL );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
       return ;
-   }
-   
-   if ( fMDTAU == 120 || fMDTAU == 22 )
-   {
+    }
+  if ( fMDTAU == 120 || fMDTAU == 22 )
+    {
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = modeL;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
-      return;      
-   }
-   
-   if ( fMDTAU == 114 || fMDTAU == 214 )
-   {
+      return;
+    }
+  if ( fMDTAU == 114 || fMDTAU == 214 )
+    {
       Tauolapp::jaki_.jak1 = modeL;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( modeL );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
-      return;      
-   }
-
-   if ( fMDTAU == 124 || fMDTAU == 224 )
-   {
+      return;
+    }
+  if ( fMDTAU == 124 || fMDTAU == 224 )
+    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = modeL;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
-      return;      
-   }
-
-   if ( fMDTAU == 115 || fMDTAU == 215 )
-   {
+      return;
+    }
+  if ( fMDTAU == 115 || fMDTAU == 215 )
+    {
       Tauolapp::jaki_.jak1 = 1;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 1 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
-      return;      
-   }
-
-   if ( fMDTAU == 125 || fMDTAU == 225 )
-   {
+      return;
+    }
+  if ( fMDTAU == 125 || fMDTAU == 225 )
+    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 1;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 1 );
-      return;      
-   }
-
-   if ( fMDTAU == 116 || fMDTAU == 216 )
-   {
+      return;
+    }
+  if ( fMDTAU == 116 || fMDTAU == 216 )
+    {
       Tauolapp::jaki_.jak1 = 2;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 2 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
-      return;      
-   }
-
-   if ( fMDTAU == 126 || fMDTAU == 226 )
-   {
+      return;
+    }
+  if ( fMDTAU == 126 || fMDTAU == 226 )
+    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 2;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 2 );
-      return;      
-   }
-
-   if ( fMDTAU == 130 || fMDTAU == 230 )
-   {
+      return;
+    }
+  if ( fMDTAU == 130 || fMDTAU == 230 )
+    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = selectHadronic();
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( Tauolapp::jaki_.jak2 );
-      return;      
-   }
-
-   if ( fMDTAU == 131 || fMDTAU == 231 )
-   {
+      return;
+    }
+  if ( fMDTAU == 131 || fMDTAU == 231 )
+    {
       Tauolapp::jaki_.jak1 = modeH;
       Tauolapp::jaki_.jak2 = 0;
       Tauolapp::Tauola::setSameParticleDecayMode( modeH );
       Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
-      return;      
-   }
-
-   if ( fMDTAU == 132 || fMDTAU == 232 )
-   {
+      return;
+    }
+  if ( fMDTAU == 132 || fMDTAU == 232 )
+    {
       Tauolapp::jaki_.jak1 = 0;
       Tauolapp::jaki_.jak2 = modeH;
       Tauolapp::Tauola::setSameParticleDecayMode( 0 );
       Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
-      return;      
-   }
-   
-   // unlikely that we get here on unknown mdtau 
-   // - there's a protection earlier
-   // but if we do, just set defaults
-   // probably need to spit a warning...
-   //
-   Tauolapp::Tauola::setSameParticleDecayMode( 0 );
-   Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
-      
-   return;
-   
-
+      return;
+    }
+  // unlikely that we get here on unknown mdtau
+  // - there's a protection earlier
+  // but if we do, just set defaults
+  // probably need to spit a warning...
+  //
+  Tauolapp::Tauola::setSameParticleDecayMode( 0 );
+  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
+  return;
 }
-
 int TauolappInterface::selectLeptonic()
 {
-   
-   float prob = flat();
-   
-   if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] ) 
-   {
+  float prob = flat();
+  if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] )
+    {
       return 1;
-   }
-   else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
-   {
+    }
+  else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
+    {
       return 2;
-   }
-      
-   return 0;
+    }
+  return 0;
 }
-
 int TauolappInterface::selectHadronic()
 {
-
-   float prob = 0.;
-   int len = 1;
-   ranmar_(&prob,&len);
-   
-   double sumBra = fScaledHadronBrRatios[0];
-   if ( prob > 0. && prob <= sumBra ) 
-   {
+  float prob = 0.;
+  int len = 1;
+  ranmar_(&prob,&len);
+  double sumBra = fScaledHadronBrRatios[0];
+  if ( prob > 0. && prob <= sumBra )
+    {
       return fHadronModes[0];
-   }
-   else
-   {
+    }
+  else
+    {
       int NN = fScaledHadronBrRatios.size();
       for ( int i=1; i<NN; i++ )
-      {
-         if ( prob > sumBra && prob <= (sumBra+fScaledHadronBrRatios[i]) ) 
-	 {
-	    return fHadronModes[i];
-	 }
-	 sumBra += fScaledHadronBrRatios[i];
-      }
-   }
-   
-   return 0;
-
+	{
+	  if ( prob > sumBra && prob <= (sumBra+fScaledHadronBrRatios[i]) )
+	    {
+	      return fHadronModes[i];
+	    }
+	  sumBra += fScaledHadronBrRatios[i];
+	}
+    }
+  return 0;
 }
-
-