Add new component MDIReader

CMakeLists.txt

@@ -38,7 +38,8 @@ endif()
 
 set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
 
-#---------------------------------------------------------------
+set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE CACHE BOOL "RPATH USE LINK PATH")
+
 add_subdirectory(k4Gen)
 
 #---------------------------------------------------------------

k4Gen/CMakeLists.txt

@@ -10,25 +10,15 @@ find_package(EvtGen)
 file(GLOB k4gen_plugin_sources src/components/*.cpp)
 gaudi_add_module(k4Gen
                  SOURCES ${k4gen_plugin_sources}
-                 LINK Gaudi::GaudiKernel
-                      ${HEPMC3_LIBRARIES} 
-                      Gaudi::GaudiAlgLib
-                      k4FWCore::k4FWCore
-                      ${EVTGEN_LIBRARIES}
-                      EDM4HEP::edm4hep
-                      EDM4HEP::edm4hepDict
-                      HepPDT::heppdt
-                  )
-
-target_include_directories(k4Gen PUBLIC 
-                           ${PYTHIA8_INCLUDE_DIRS} 
-                           ${HEPMC3_INCLUDE_DIR}
-                           ${EVTGEN_INCLUDE_DIR}
-                           $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/include>
-                           $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
-                          )
-
-target_link_options(k4Gen PRIVATE "-Wl,--enable-new-dtags")
+                 LINK Gaudi::GaudiKernel ${HEPMC3_LIBRARIES} Gaudi::GaudiAlgLib k4FWCore::k4FWCore ${HEPPDT_LIBRARIES} ${EVTGEN_LIBRARIES} EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${HEPMC_LIBRARIES})
+
+target_include_directories(k4Gen PUBLIC ${PYTHIA8_INCLUDE_DIRS} 
+                                        ${HEPMC3_INCLUDE_DIR}
+                                        ${HEPPDT_INCLUDE_DIR}
+                                        ${EVTGEN_INCLUDE_DIR}
+					${HEPMC_INCLUDE_DIR}
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/include>
+  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
 
 
 install(TARGETS k4Gen
@@ -79,6 +69,12 @@ set_tests_properties(Pythia8ExtraSettings PROPERTIES
   )
 set_test_env(Pythia8ExtraSettings)
 
+add_test(NAME MDIreader
+	      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
+	      COMMAND  k4run ${CMAKE_CURRENT_LIST_DIR}/options/mdireader_test.py
+	      )
+set_test_env(MDIreader)
+
 #--- Install the example options to the directory where the spack installation
 #--- points the $K4GEN environment variable
 install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/options

k4Gen/options/mdireader_test.py

@@ -0,0 +1,68 @@
+from Gaudi.Configuration import *
+import os
+
+# Data service
+from Configurables import FCCDataSvc
+podioevent = FCCDataSvc("EventDataSvc")
+
+from Configurables import MDIReader
+#mdi_converter = MDIReader("Reader",MDIFilename="k4Gen/options/mdireader_testparticles.dat")
+mdi_converter = MDIReader("Reader",MDIFilename=(os.path.join(os.environ["K4GEN"],'../options/mdireader_testparticles.dat')))
+mdi_converter.GenParticles.Path = "allGenParticles"
+mdi_converter.CrossingAngle = 0.015
+mdi_converter.LongitudinalCut = 0
+mdi_converter.InputType = "xtrack"
+mdi_converter.BeamEnergy = 45.6
+
+## DD4hep geometry service
+## Parses the given xml file
+#from Configurables import GeoSvc
+#geoservice = GeoSvc("GeoSvc", detectors=[
+ #                                        #os.path.join(os.environ["FCCDETECTORS"], 'Detector/DetFCCeeCLD/compact/FCCee_DectMaster.xml'),
+  #                                       #os.path.join(os.environ["FCCDETECTORS"], 'Detector/DetFCCeeCLD/compact/FCCee_o2_v02/FCCee_o2_v02.xml'),
+   #                                      os.path.join(os.environ["FCCDETECTORS"], 'Detector/DetCommon/compact/Box.xml'),
+    #                                     ])
+
+## Geant4 service
+## Configures the Geant simulation: geometry, physics list and user actions
+#from Configurables import SimG4Svc, SimG4FullSimActions, SimG4UserLimitPhysicsList, SimG4UserLimitRegion
+#from GaudiKernel.SystemOfUnits import mm
+
+
+## giving the names of tools will initialize the tools of that type
+#geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist="SimG4FtfpBert"
+ #                       )
+
+
+## Geant4 algorithm
+## Translates EDM to G4Event, passes the event to G4, writes out outputs via tools
+#from Configurables import SimG4Alg, SimG4PrimariesFromEdmTool
+
+#particle_converter = SimG4PrimariesFromEdmTool("EdmConverter")
+#particle_converter.GenParticles.Path = "allGenParticles"
+
+
+#geantsim = SimG4Alg("SimG4Alg",
+ #                   outputs= [
+  #      ],
+   #                 eventProvider=particle_converter)
+
+
+# PODIO algorithm
+from Configurables import PodioOutput
+out = PodioOutput("out",
+                  filename="mdireader_test_out.root",
+                  OutputLevel=INFO)
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr( #TopAlg = [mdi_converter, geantsim, out],
+                TopAlg = [mdi_converter, out], 
+                EvtSel = 'NONE',
+                EvtMax   = 1,
+                # order is important, as GeoSvc is needed by SimG4Svc
+                #ExtSvc = [podioevent, geoservice, geantservice],
+                ExtSvc = [podioevent],
+                OutputLevel=INFO
+               )

k4Gen/options/mdireader_testparticles.dat

@@ -0,0 +1,3 @@
+0.0		     0.0		 0.0		    0.0			 0.0		    0.0		       0.0
+0.5		     -20.0e-3		 20.0e-3	    0.5			 -0.5		    0.0		       -0.001
+-0.5		     20.0e-3		 -20.0e-3	    -0.5		 0.5		    0.0		       0.001

k4Gen/src/components/MDIReader.cpp

@@ -0,0 +1,216 @@
+
+#include "MDIReader.h"
+
+#include "GaudiKernel/IIncidentSvc.h"
+#include "GaudiKernel/Incident.h"
+#include "GaudiKernel/IEventProcessor.h"
+
+#include "edm4hep/MCParticleCollection.h"
+
+#include <math.h>
+#include <algorithm>
+#include <iostream>
+
+DECLARE_COMPONENT(MDIReader)
+
+MDIReader::MDIReader(const std::string& name, ISvcLocator* svcLoc):
+GaudiAlgorithm(name, svcLoc),
+  m_filename() {
+
+  declareProperty("MDIFilename", m_filename="", "Name of the MDI file to read");
+  declareProperty("GenParticles", m_genphandle, "Generated particles collection (output)");
+  declareProperty("CrossingAngle",xing,"Half the crossing angle beam in [rad]");
+  declareProperty("LongitudinalCut",cut_z,"the value for cut_z used in GP++ in [um]");
+  declareProperty("InputType",input_type,"string: guineapig, xtrack");
+  declareProperty("BeamEnergy",beam_energy,"beam energy [GeV], necessary for xtrack type");
+}
+
+StatusCode MDIReader::initialize()
+{
+  StatusCode sc = GaudiAlgorithm::initialize();
+
+  m_input.open(m_filename.c_str(), std::ifstream::in);
+
+  if ( !m_input.good() )   {
+    error() << "Failed to open input stream:"+m_filename << endmsg;
+    return StatusCode::FAILURE;
+  }
+  return StatusCode::SUCCESS;
+}
+
+StatusCode MDIReader::execute()
+{
+  // First check the input file status
+  if(m_input.eof())
+    {
+      error() << "End of file reached" << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+  // Check the input type flag
+  if(input_type!="guineapig" && input_type!="xtrack")
+    {
+      error() << "Input type flag - wrong definition: "<< input_type << endmsg;
+      return StatusCode::FAILURE;
+    }
+  else{
+    debug() <<"Selected input type : "<<input_type<< endmsg;
+  }
+
+  //  Loop over particles
+  int ISTHEP = 1;   // status code
+  int IDHEP;    // PDG code
+  int CHARGE;   // charge
+  //int JMOHEP1;  // first mother
+  //int JMOHEP2;  // last mother
+  //int JDAHEP1;  // first daughter
+  //int JDAHEP2;  // last daughter
+  double PHEP1; // px in GeV/c
+  double PHEP2; // py in GeV/c
+  double PHEP3; // pz in GeV/c
+  double PHEP4; // energy in GeV
+  double PHEP5; // mass in GeV/c**2
+  double VHEP1; // x vertex position in mm
+  double VHEP2; // y vertex position in mm
+  double VHEP3; // z vertex position in mm
+  double VHEP4 = 0; // production time in mm/c
+
+  //guineapig IPP variables which do not end up in the EDM format
+  double process; // 1=Breit-Wheeler 2=Bethe-Heitler 3=Landau-Lifshitz
+  double trash;   // trash variable
+  double id_ee;   // same id means they are a pair
+  double temp_x,temp_y,temp_z,temp_px, temp_py,temp_pz,temp_e;
+
+  debug() <<"The crossing angle is "<<xing<<" [rad]"<< endmsg;
+  //std::cout <<"The crossing angle is "<<xing<<" [rad]"<< endmsg;
+  edm4hep::MCParticleCollection* particles = new edm4hep::MCParticleCollection();
+
+
+  PHEP5 = 5.11e-4;
+  while(m_input.good())
+    {
+      if(input_type=="guineapig"){
+	m_input >> PHEP4
+		>> PHEP1 >> PHEP2 >> PHEP3
+		>> VHEP1 >> VHEP2 >> VHEP3
+		>> process >> trash >> id_ee;
+
+
+	//std::cout<<PHEP4<<" "<<sqrt(PHEP1*PHEP1 + PHEP2*PHEP2 + PHEP3*PHEP3)<<" "<<sqrt((PHEP1*PHEP1 + PHEP2*PHEP2 + PHEP3*PHEP3)*PHEP4*PHEP4 + PHEP5*PHEP5)<<std::endl;
+
+	if(m_input.eof())break;
+	else if(!m_input.good())
+	  {
+	    debug() << "End of file reached before reading all the hits" << endmsg;
+	    error() << "End of file reached before reading all the hits" << endmsg;
+	    return StatusCode::FAILURE;
+	  }
+
+	IDHEP = 11;
+	CHARGE = -1;
+	if(PHEP4<0){
+	  IDHEP=-11;
+	  CHARGE=1;
+	}
+
+	VHEP1 *=1e-9; //convert from nm to m
+	VHEP2 *=1e-9; //convert from nm to m
+	VHEP3 *=1e-9; //convert from nm to m
+
+	//---boost section
+	PHEP4 = abs(PHEP4);
+	temp_x  = cut_z*1e-6*tan(xing) + VHEP1*sqrt(1+pow(tan(xing),2));
+	temp_px = PHEP4*tan(xing) + PHEP1*PHEP4*sqrt(1+pow(tan(xing),2));
+	temp_e  = abs( PHEP4*sqrt(1+pow(tan(xing),2)) + PHEP1*PHEP4*tan(xing) );
+
+	VHEP1 = temp_x;
+	PHEP1 = temp_px;
+	PHEP2 = PHEP2*PHEP4;
+	PHEP3 = PHEP3*PHEP4;
+	PHEP4 = temp_e;
+
+	//---end boost section
+
+	VHEP1 *=1e3; //convert from m to mm
+	VHEP2 *=1e3; //convert from m to mm 
+	VHEP3 *=1e3; //convert from m to mm
+
+	//boost to lab frame
+
+	VHEP1=0;
+	VHEP2=0;
+	VHEP3=0;
+      }//end if guineapig
+
+      else if(input_type=="xtrack"){
+	m_input >> VHEP3
+		>> VHEP1 >> VHEP2
+		>> PHEP1 >> PHEP2
+		>> temp_z >> PHEP4 ;
+
+	if(m_input.eof())break;
+	else if(!m_input.good())
+	  {
+	    debug() << "End of file reached before reading all the hits" << endmsg;
+	    error() << "End of file reached before reading all the hits" << endmsg;
+	    return StatusCode::FAILURE;
+	  }
+
+	//VHEP3 = VHEP3 + temp_z; //is ct necessary? prob not
+
+	//---position rotation in CLD frame
+	temp_z =  VHEP3*cos(-xing) + VHEP1*sin(-xing);
+	temp_x = -VHEP3*sin(-xing) + VHEP1*cos(-xing);
+	VHEP1 = temp_x;
+	VHEP3 = temp_z;
+
+	VHEP1 *=1e3; //convert from m to mm
+	VHEP2 *=1e3; //convert from m to mm 
+	VHEP3 *=1e3; //convert from m to mm
+
+	PHEP4 = (1.+PHEP4)*beam_energy;
+	PHEP1 = PHEP1*beam_energy;
+	PHEP2 = PHEP2*beam_energy;
+	PHEP3 = sqrt(pow(PHEP4,2) - pow(PHEP1,2) - pow(PHEP2,2) );
+
+	//---momentum rotation in CLD frame
+	temp_pz =  PHEP3*cos(-xing) + PHEP1*sin(-xing);
+	temp_px = -PHEP3*sin(-xing) + PHEP1*cos(-xing);
+	PHEP3 = temp_pz;
+	PHEP1 = temp_px;
+
+	IDHEP = 11;
+	CHARGE = -1;
+
+      }	//end if xtrack
+
+
+      edm4hep::MutableMCParticle particle = particles->create();
+
+      particle.setPDG(IDHEP);
+      particle.setCharge(CHARGE);
+      particle.setGeneratorStatus(ISTHEP); 
+      particle.setMomentum({
+	  (float) (PHEP1),
+	  (float) (PHEP2),
+	  (float) (PHEP3),
+	});
+      particle.setMass(PHEP5);
+      particle.setVertex({
+                          VHEP1,
+                          VHEP2,
+                          VHEP3,
+                        });
+      particle.setTime(VHEP4);
+    }
+
+    m_genphandle.put(particles);
+  return StatusCode::SUCCESS;
+}
+
+StatusCode MDIReader::finalize() 
+{
+  m_input.close();
+  debug() << "MDIReader finalization" << endmsg;
+  return GaudiAlgorithm::finalize();
+}