Handle tau polarization correctly for Genie events

Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx

@@ -5,8 +5,7 @@
 // Header for my class
 #include "Pythia8Decayer.h"
 
-// The actual decayers.  Singleton classes, not toolhandles
-// #include "Pythia8Instance.h"
+#include "FaserMCTruth/FaserTrackInformation.h"
 
 // For passing things around
 #include "CLHEP/Vector/LorentzVector.h"
@@ -14,6 +13,7 @@
 #include "G4DynamicParticle.hh"
 #include "G4ParticleTable.hh"
 #include "G4DecayProducts.hh"
+#include "G4VProcess.hh"
 
 Pythia8Decayer::Pythia8Decayer( const std::string s )
  : G4VExtDecayer(s)
@@ -84,15 +84,26 @@ G4DecayProducts* Pythia8Decayer::ImportDecayProducts(const G4Track& aTrack){
 
    // specify polarization, if any
    
+   // special logic for primary taus(anti-taus), assumed to have polarization -1(+1), respectively
+   // verified from the polarization info in Genie output
+   double spinup;
+   FaserTrackInformation* info = dynamic_cast<FaserTrackInformation*>(aTrack.GetUserInformation());
+   if (info != nullptr && abs(pdgid) == 15 && (info->GetClassification() == TrackClassification::Primary || info->GetClassification() == TrackClassification::RegeneratedPrimary))
+   {
+       G4cout << "Primary tau decay identified." << G4endl;
+       spinup = (pdgid > 0 ? -1 : +1);
+   }
+   else
    // NOTE: while in Py8 polarization is a double variable , 
    //       in reality it's expected to be -1, 0., or 1 in case of "external" tau's, 
    //       similar to LHA SPINUP; see Particle Decays, Hadron and Tau Decays in docs at
    //       https://pythia.org/manuals/pythia8305/Welcome.html
    //       so it's not able to handle anything like 0.99, thus we're rounding off    
-   m_decayer->event.back().pol( 
-      round( std::cos( aTrack.GetPolarization().angle( aTrack.GetMomentumDirection() ) )
-      ) 
-   );
+   {
+       spinup = round( std::cos( aTrack.GetPolarization().angle( aTrack.GetMomentumDirection() ) ) );
+   }
+   G4cout << "Using " << aTrack.GetParticleDefinition()->GetParticleName() << " helicity " << spinup << " for Pythia8 decay." << G4endl;
+   m_decayer->event.back().pol( spinup );
 
    int npart_before_decay = m_decayer->event.size();