Muon DIS with Pythia6 run through Geant4

sndFairTasks/CMakeLists.txt

@@ -10,6 +10,9 @@ ${CMAKE_SOURCE_DIR}/sndFairTasks
 ${CMAKE_SOURCE_DIR}/veto
 ${XROOTD_INCLUDE_DIR}
 ${XROOTD_INCLUDE_DIR}/..
+${GEANT4_ROOT}/include/Geant4
+${GEANT4_VMC_ROOT}/include/geant4vmc
+${CMAKE_SOURCE_DIR}/G4Processes
 )
 
 include_directories( ${INCLUDE_DIRECTORIES}  ${VMC_INCLUDE_DIRS}  ${FAIRROOT_INCLUDE_DIR} ${FairLogger_INCDIR} ${FMT_INCLUDE_DIR}
@@ -20,7 +23,7 @@ set(LINK_DIRECTORIES
 ${ROOT_LIBRARY_DIR}
 ${XROOTD_LIBRARY_DIR}
 ${FAIRROOT_LIBRARY_DIR}
-
+${G4Processes_LIBRARY_DIR}
 )
 
 link_directories( ${LINK_DIRECTORIES})
@@ -30,11 +33,12 @@ DigiTaskSND.cxx
 ConvRawData.cxx
 boardMappingParser.cxx
 MCEventBuilder.cxx
+MuonDISProcessInjector.cxx
 )
 
 Set(HEADERS)
 Set(LINKDEF sndFairTasksLinkDef.h)
 Set(LIBRARY_NAME sndFairTasks)
-Set(DEPENDENCIES Base ShipData shipLHC GeoBase ParBase Geom Core)
+Set(DEPENDENCIES Base ShipData shipLHC GeoBase ParBase Geom Core G4Processes)
 
 GENERATE_LIBRARY()

sndFairTasks/MuonDISProcessInjector.cxx

@@ -0,0 +1,86 @@
+#include "MuonDISProcessInjector.h"
+#include "G4MuonMinus.hh"
+#include "G4MuonPlus.hh"
+#include "G4ProcessManager.hh"
+#include "G4MuonDISProcess.hh"
+#include "FairLogger.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include <exception>
+#include <vector>
+
+using std::vector;
+
+// Convert ROOT's TGraph to G4PhysicsFreeVector
+G4PhysicsFreeVector* GraphToFreeVector(TGraph* graph) {
+    auto xsec_data = new G4PhysicsFreeVector();
+    for (int i = 0; i < graph->GetN(); ++i) {
+        double x, y;
+        graph->GetPoint(i, x, y);
+        xsec_data->InsertValues(x, y);
+    }
+    return xsec_data;
+}
+
+MuonDISProcessInjector::MuonDISProcessInjector(char *nucleon, vector<float> x_range, vector<float> y_range,
+                                               vector<float> z_range, char *volume, char *xsec_filename)
+   : FairTask("MuonDISProcessInjector")
+{
+   fNucleon = nucleon;
+   fVolumeName = volume;
+   fXRange = x_range;
+   fYRange = y_range;
+   fZRange = z_range;
+
+   TFile* xsec_file = new TFile(xsec_filename);
+   if(!xsec_file || xsec_file->IsZombie()){
+      LOG(FATAL) << "DIS cross section file not found";
+      exit(0);
+   }
+   TGraph* crsec_muminus = static_cast<TGraph*> (xsec_file->Get(Form("g_13_%s", fNucleon)));
+   TGraph* crsec_muplus = static_cast<TGraph*> (xsec_file->Get(Form("g_-13_%s", fNucleon)));
+   if (!crsec_muminus || !crsec_muplus){
+     LOG(FATAL) << "DIS cross section graphs per " << fNucleon <<" are missing in "<< xsec_filename;
+     exit(0);
+   }
+
+   fXsecTables = std::make_shared<std::map<int, G4PhysicsFreeVector*>>();
+   fXsecTables->insert({13, GraphToFreeVector(crsec_muminus)});
+   fXsecTables->insert({-13, GraphToFreeVector(crsec_muplus)});
+
+   LOG(WARNING) << "MuonDISProcessInjector: Setting xyz ranges[cm] for muon DIS\nx: "
+                << fXRange[0] << ", " << fXRange[1] << "\ny: " << fYRange[0] << ", " << fYRange[1] << "\nz: "
+                << fZRange[0] << ", " << fZRange[1] << "\nand volume name '" << fVolumeName
+                << "' and nucleon type " << fNucleon;
+}
+
+InitStatus MuonDISProcessInjector::Init()
+{
+   muMinus = G4MuonMinus::Definition();
+   muPlus = G4MuonPlus::Definition();
+
+   DISProcess = new G4MuonDISProcess();
+
+   // Set the nucleon type
+   DISProcess->SetNucleonType(fNucleon);
+   // Set DIS xsec lookup tables per set nucleon type
+   DISProcess->SetCrossSectionTables(fXsecTables);
+   // Set the z range of the DIS interaction
+   DISProcess->SetRange(&fXRange, &fYRange, &fZRange);
+   // Set the name of the geo volume for the DIS interaction
+   DISProcess->SetVolume(fVolumeName);
+
+   auto process_man_muminus = muMinus->GetProcessManager();
+   if (process_man_muminus) {
+      G4int id_muminus = process_man_muminus->AddDiscreteProcess(DISProcess);
+   }
+   auto process_man_muplus = muPlus->GetProcessManager();
+   if (process_man_muplus) {
+      G4int id_muplus = process_man_muplus->AddDiscreteProcess(DISProcess);
+   }
+
+   LOG(WARNING) << "Adding the external generation for muon DIS using Pythia6";
+   return kSUCCESS;
+}
+
+ClassImp(MuonDISProcessInjector)

sndFairTasks/MuonDISProcessInjector.h

@@ -0,0 +1,33 @@
+#ifndef MUONDISPROCESSINJECTOR_H
+#define MUONDISPROCESSINJECTOR_H
+
+#include "FairTask.h"
+#include "G4ParticleDefinition.hh"
+#include "G4MuonDISProcess.hh"
+#include <vector>
+#include <map>
+
+class G4PhysicsFreeVector;
+
+class MuonDISProcessInjector : public FairTask {
+public:
+   MuonDISProcessInjector(char *nucleon, std::vector<float> x_range, std::vector<float> y_range,
+                          std::vector<float> z_range, char *volume, char *xsec_filename);
+   virtual ~MuonDISProcessInjector() {}
+
+   virtual InitStatus Init();
+
+   char *fNucleon;             //! nucleon type
+   char *fVolumeName;          //! geometry volume name to place the DIS
+   std::vector<float> fXRange; /// x range to place the DIS
+   std::vector<float> fYRange; /// y range to place the DIS
+   std::vector<float> fZRange; /// z range to place the DIS
+   std::shared_ptr<std::map<int, G4PhysicsFreeVector*>> fXsecTables; /// DIS cross section lookup tables
+   const G4ParticleDefinition *muPlus;
+   const G4ParticleDefinition *muMinus;
+   G4MuonDISProcess *DISProcess;
+
+   ClassDef(MuonDISProcessInjector, 1);
+};
+
+#endif

sndFairTasks/sndFairTasksLinkDef.h

@@ -7,8 +7,5 @@
 #pragma link C++ class DigiTaskSND;
 #pragma link C++ class ConvRawData;
 #pragma link C++ class MCEventBuilder;
+#pragma link C++ class MuonDISProcessInjector;
 #endif
-
-
-
-