Move the WCSimTrackingAction code into HKG4TrackingAction & multiply …

…inherit from Tool & G4UserTrackingAction. Untested (like last couple of commits) since I forgot to add new tools into the toolchain until now, and now I need to toolify more tools to setup G4 user classes in the correct order

HKG4TrackingAction/HKG4TrackingAction.cpp

@@ -1,6 +1,17 @@
 #include "HKG4TrackingAction.h"
 
-#include "WCSimTrackingAction.hh"
+#include "G4ParticleTypes.hh"
+#include "G4TrackingManager.hh"
+#include "G4Track.hh"
+#include "G4ios.hh"
+#include "G4VProcess.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+
+#include "WCSimTrajectory.hh"
+#include "WCSimTrackInformation.hh"
+#include "WCSimTrackingMessenger.hh"
+#include "WCSimPrimaryGeneratorAction.hh"
 
 using namespace HK::Ghost::G4;
 
@@ -19,7 +30,52 @@ bool HKG4TrackingAction::Initialise(std::string configfile, DataModel& data) {
 		m_verbose = 1;
 
 	m_p_UI = G4UImanager::GetUIpointer();
-	m_data->m_p_run_manager->SetUserAction(new WCSimTrackingAction);
+
+	ProcessList.insert("Decay") ;                         // Michel e- from pi+ and mu+
+  ProcessList.insert("conv") ;                         // Products of gamma conversion
+
+  //ProcessList.insert("muMinusCaptureAtRest") ;          // Includes Muon decay from K-shell: for Michel e- from mu0. This dominates/replaces the mu- decay (noticed when switching off this process in PhysicsList)                                                   // TF: IMPORTANT: ONLY USE FROM G4.9.6 onwards because buggy/double counting before.
+  ////////// ToDo: switch ON the above when NuPRISM uses G4 >= 4.9.6
+  ProcessList.insert("nCapture");
+
+//   ProcessList.insert("conv");
+
+  // F. Nova One can check here if the photon comes from WLS
+  ProcessList.insert("OpWLS");
+
+  ParticleList.insert(111); // pi0
+  ParticleList.insert(211); // pion+
+  ParticleList.insert(-211);
+  ParticleList.insert(321);
+  ParticleList.insert(-321); // kaon-
+  ParticleList.insert(311); // kaon0
+  ParticleList.insert(-311); // kaon0 bar
+  //ParticleList.insert(22); // I add photons (B.Q)
+  ParticleList.insert(11); // e-
+  ParticleList.insert(-11); // e+
+  ParticleList.insert(13); // mu-
+  ParticleList.insert(-13); // mu+
+  // don't put gammas there or there'll be too many
+
+  //TF: add protons and neutrons
+  ParticleList.insert(2212);
+  ParticleList.insert(2112);
+
+  percentageOfCherenkovPhotonsToDraw = 0.0;
+#ifndef WCSIM_SAVE_PHOTON_HISTORY
+  SAVE_PHOTON_HISTORY = false;
+#else
+  SAVE_PHOTON_HISTORY = true;
+#endif
+
+  //messenger = new WCSimTrackingMessenger(this);
+	*m_log << ML(0) << "Need to put WCSimTrackingMessenger functionality into the HKG4TrackingAction tool config file" << std::endl;
+
+  // Max time for radioactive decay:
+  fMaxTime    = 1. * CLHEP::second;
+  fTime_birth = 0.;
+
+	m_data->m_p_run_manager->SetUserAction(this);
 
 	return true;
 }
@@ -33,3 +89,163 @@ bool HKG4TrackingAction::Finalise() {
 
 	return true;
 }
+
+
+void HKG4TrackingAction::PreUserTrackingAction(const G4Track* aTrack)
+{
+  //TF: userdefined now
+  //G4double percentageOfCherenkovPhotonsToDraw = 100.0;
+  // if larger than zero, will keep trajectories of many secondaries as well
+  // and store them in output file. Difficult to control them all, so best only
+  // use for visualization, not for storing in ROOT.
+
+  if ( aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()
+       || G4UniformRand() < percentageOfCherenkovPhotonsToDraw/100. )
+    {
+      WCSimTrajectory* thisTrajectory = new WCSimTrajectory(aTrack);
+      thisTrajectory->SetSavePhotonTrack(true); // mark to save photon track in ROOT
+      fpTrackingManager->SetTrajectory(thisTrajectory);
+      fpTrackingManager->SetStoreTrajectory(true);
+    }
+  else if (SAVE_PHOTON_HISTORY)
+    {
+      // Keep the trajectory but not saving in ROOT
+      WCSimTrajectory* thisTrajectory = new WCSimTrajectory(aTrack);
+      thisTrajectory->SetSavePhotonTrack(false);
+      fpTrackingManager->SetTrajectory(thisTrajectory);
+      fpTrackingManager->SetStoreTrajectory(true);
+    }
+  else 
+    fpTrackingManager->SetStoreTrajectory(false);
+
+  WCSimPrimaryGeneratorAction *primaryGenerator = (WCSimPrimaryGeneratorAction *) (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
+  if(!primaryGenerator->IsConversionFound()) {
+    if(aTrack->GetParentID()==0){
+      primaryID = aTrack->GetTrackID();
+    }
+    else if(aTrack->GetParentID() == primaryID) {
+      if (aTrack->GetCreatorProcess()->GetProcessName() == "conv") {
+          primaryGenerator->FoundConversion();
+      }
+      G4EventManager::GetEventManager()->AbortCurrentEvent();
+      G4EventManager::GetEventManager()->GetNonconstCurrentEvent()->SetEventAborted();
+    }
+  }
+
+  // Kill nucleus generated after TrackID 1
+  G4ParticleDefinition* particle = aTrack->GetDefinition();
+  G4String name   = particle->GetParticleName();
+  G4double fCharge = particle->GetPDGCharge();
+
+  G4Track* tr = (G4Track*) aTrack;
+  if ( aTrack->GetTrackID() == 1 ) {
+  	// Re-initialize time
+  	fTime_birth = 0;
+  	// Ask G4 to kill the track when all secondary are done (will exclude other decays)
+  	if ( fCharge > 2. )
+  		tr->SetTrackStatus(fStopButAlive);
+  }
+
+  if ( aTrack->GetTrackID() == 2 ) {
+  	// First track of the decay save time
+  	fTime_birth = aTrack->GetGlobalTime(); 
+  }
+}
+
+void HKG4TrackingAction::PostUserTrackingAction(const G4Track* aTrack)
+{
+  // added by M Fechner
+  const G4VProcess* creatorProcess = aTrack->GetCreatorProcess();
+  //  if ( creatorProcess )
+
+
+  WCSimTrackInformation* anInfo;
+  if (aTrack->GetUserInformation())
+    anInfo = (WCSimTrackInformation*)(aTrack->GetUserInformation());   //eg. propagated to all secondaries blelow.
+  else anInfo = new WCSimTrackInformation();
+
+    // Simplified tracking: track all primaries, particles from the chosen list of particle types or creation processes,
+    // Optionally, all particles producing cherekov hits and their parents, grandparents, etc. will be added later.
+
+    // is it a primary ?
+    // is the process in the set ? 
+    // is the particle in the set ?
+    if( aTrack->GetParentID()==0
+	|| ((creatorProcess!=0) && ProcessList.count(creatorProcess->GetProcessName()))
+	|| (ParticleList.count(aTrack->GetDefinition()->GetPDGEncoding()))
+      )
+    {
+      // if so the track is worth saving
+      anInfo->WillBeSaved(true);
+    }
+  else {
+    anInfo->WillBeSaved(false);
+  }
+
+
+  G4Track* theTrack = (G4Track*)aTrack;
+  theTrack->SetUserInformation(anInfo);
+
+  // pass parent trajectory to children
+  G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
+  WCSimTrajectory *currentTrajectory = (WCSimTrajectory*)fpTrackingManager->GimmeTrajectory();
+  if(currentTrajectory && !anInfo->GetMyTrajectory())
+    anInfo->SetMyTrajectory(currentTrajectory);
+  if(secondaries)
+  {
+    size_t nSeco = secondaries->size();
+    if(nSeco>0)
+    {
+      for(size_t i=0;i<nSeco;i++)
+      {
+        WCSimTrackInformation* infoSec = new WCSimTrackInformation();
+        infoSec->SetParentTrajectory(anInfo->GetMyTrajectory());
+        (*secondaries)[i]->SetUserInformation(infoSec);
+      }
+    } 
+  }
+
+  // If this track produces a hit, traverse back through parent trajectories to flag that the parents produce a hit
+  if (anInfo->GetProducesHit() && saveHitProducingTracks){
+      WCSimTrajectory* parentTrajectory = anInfo->GetParentTrajectory();
+      while(parentTrajectory != 0 && !parentTrajectory->GetProducesHit()){
+          if (parentTrajectory->GetPDGEncoding()==0 && !parentTrajectory->GetSavePhotonTrack()) break; // do not save unwanted photon tracks
+          parentTrajectory->SetProducesHit(true);
+          parentTrajectory = parentTrajectory->GetParentTrajectory();
+      }
+  }
+
+  if (currentTrajectory)
+    //   if (aTrack->GetDefinition()->GetPDGCharge() == 0) 
+  {
+
+    G4ThreeVector currentPosition      = aTrack->GetPosition();
+    G4VPhysicalVolume* currentVolume   = aTrack->GetVolume();
+
+    currentTrajectory->SetStoppingPoint(currentPosition);
+    currentTrajectory->SetStoppingVolume(currentVolume);
+
+    if (aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
+    {
+      currentTrajectory->SetSaveFlag(anInfo->isSaved());// mark it for WCSimEventAction ;
+      // don't save the optical photon tracks themselves simply when they produce a hit, since that info is already saved in WCSimRootCherenkovHitTime
+      // optical photon tracks can still be saved if wanted by explicitly adding appropriate entries to the ParticleList or ProcessList via mac file commands
+      currentTrajectory->SetProducesHit(anInfo->GetProducesHit());
+    }
+    else if (currentTrajectory->GetSavePhotonTrack()) // only save the wanted photon tracks
+      currentTrajectory->SetSaveFlag(anInfo->isSaved()); 
+    else 
+      currentTrajectory->SetSaveFlag(false);
+  }
+
+  WCSimPrimaryGeneratorAction *primaryGenerator = (WCSimPrimaryGeneratorAction *) (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
+  if(!primaryGenerator->IsConversionFound() && 
+     aTrack->GetTrackID() == primaryID &&
+     aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep() &&
+     aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "conv"){
+      for(int i=0; i<2; i++){
+          primaryGenerator->SetConversionProductParticle(i, secondaries->at(i)->GetParticleDefinition());
+          primaryGenerator->SetConversionProductMomentum(i, secondaries->at(i)->GetMomentum());
+      }
+  }
+}

HKG4TrackingAction/HKG4TrackingAction.h

@@ -8,6 +8,7 @@
 #include "Tool.h"
 
 #include "G4UImanager.hh"
+#include "G4UserTrackingAction.hh"
 
 #include "WCSimTrackingAction.hh"
 
@@ -22,21 +23,49 @@
 namespace HK {
 	namespace Ghost {
 		namespace G4 {
-			class HKG4TrackingAction : public Tool {
+			class HKG4TrackingAction : public Tool, public G4UserTrackingAction {
 
 			public:
 
 				HKG4TrackingAction();  ///< Simple constructor
 				bool Initialise(std::string configfile,
-												DataModel& data);  ///< Initialise Function for setting up Tool resources. @param
+												DataModel& data) override;  ///< Initialise Function for setting up Tool resources. @param
 				///< configfile The path and name of the dynamic configuration file
 				///< to read in. @param data A reference to the transient data
 				///< class used to pass information between Tools.
-				bool Execute();                    ///< Execute function used to perform Tool purpose.
-				bool Finalise();                   ///< Finalise funciton used to clean up resources.
+				bool Execute() override;                    ///< Execute function used to perform Tool purpose.
+				bool Finalise() override;                   ///< Finalise funciton used to clean up resources.
 
 			private:
 				G4UImanager* m_p_UI;
+
+			public:
+				void PreUserTrackingAction (const G4Track* aTrack) override;
+				void PostUserTrackingAction(const G4Track*) override;
+
+				void SetFractionChPhotons(G4double fraction){percentageOfCherenkovPhotonsToDraw = fraction;}
+
+				void AddProcess(const G4String &process){ProcessList.insert(process);}
+				void AddParticle(G4int pid){ParticleList.insert(pid);}
+				void SetSaveHitProducingTracks(G4bool save){saveHitProducingTracks = save;}
+				G4bool GetSaveHitProducingTracks() const {return saveHitProducingTracks;}
+
+			private:
+				std::set<G4String> ProcessList;
+				std::set<G4int> ParticleList;
+				G4bool saveHitProducingTracks = true;
+
+				G4double fTime_birth;
+				G4double fMaxTime;
+				// TF: define in macro now
+				G4double percentageOfCherenkovPhotonsToDraw;
+
+				bool SAVE_PHOTON_HISTORY;
+
+				WCSimTrackingMessenger* messenger;
+
+				G4int primaryID;
+
 			};
 		} // namespace G4
 	} // namespace Ghost