From 527500f3f87298d14a346595d60b856f396c6029 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sun, 12 Jun 2022 09:06:39 -0700
Subject: [PATCH] Handle tau polarization correctly

---
 .../Pythia8Decayer/CMakeLists.txt             |  2 +-
 .../Pythia8Decayer/src/Pythia8Decayer.cxx     | 23 ++++++++++++++-----
 2 files changed, 18 insertions(+), 7 deletions(-)

diff --git a/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
index 901e0e3e..57dcf64b 100644
--- a/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
+++ b/Simulation/G4Extensions/Pythia8Decayer/CMakeLists.txt
@@ -18,7 +18,7 @@ atlas_subdir( Pythia8Decayer )
                        src/*.cxx
                        NO_PUBLIC_HEADERS
                        PRIVATE_INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
-                       PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib Pythia8_iLib )
+                       PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} FaserMCTruth GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib Pythia8_iLib )
 
     atlas_add_component( Pythia8Decayer
                          src/components/*.cxx
diff --git a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx b/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
index 7047c554..a6f9ea09 100644
--- a/Simulation/G4Extensions/Pythia8Decayer/src/Pythia8Decayer.cxx
+++ b/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();
    
-- 
GitLab