diff --git a/CMakeLists.txt b/CMakeLists.txt index 198f641..716cb5c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,6 +21,7 @@ find_package(ROOT COMPONENTS RIO Tree MathCore) find_package(EDM4HEP) find_package(k4FWCore) find_package(Gaudi) +find_package(GSL) #--------------------------------------------------------------- include(cmake/Key4hepConfig.cmake) diff --git a/DCHdigi/CMakeLists.txt b/DCHdigi/CMakeLists.txt index e378554..d34d7d3 100644 --- a/DCHdigi/CMakeLists.txt +++ b/DCHdigi/CMakeLists.txt @@ -28,7 +28,7 @@ install(FILES install(FILES dataFormatExtension/driftChamberHit.yaml - DESTINATION "${CMAKE_INSTALL_DATADIR}/extension" COMPONENT dev) + DESTINATION "${CMAKE_INSTALL_PREFIX}/extension" COMPONENT dev) install(FILES "${PROJECT_BINARY_DIR}/libextensionDict_rdict.pcm" @@ -42,15 +42,29 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h ) +include(CheckIncludeFileCXX) +set(CMAKE_REQUIRED_LIBRARIES DD4hep::DDRec) +CHECK_INCLUDE_FILE_CXX(DDRec/DCH_info.h DCH_INFO_H_EXIST) +set(CMAKE_REQUIRED_LIBRARIES) +set(FILES_DEPENDINGON_DCH_INFO_H "DCHdigi.cpp" ) +if(NOT DCH_INFO_H_EXIST) + list(FILTER sources EXCLUDE REGEX "${FILES_DEPENDINGON_DCH_INFO_H}") + message(WARNING "Gaudi algorithm defined in ${FILES_DEPENDINGON_DCH_INFO_H} will not be built because header file DDRec/DCH_info.h was not found") +endif() + gaudi_add_module(${PackageName} SOURCES ${sources} LINK + k4FWCore::k4FWCore + k4FWCore::k4Interface Gaudi::GaudiKernel EDM4HEP::edm4hep extensionDict - k4FWCore::k4FWCore - k4FWCore::k4Interface DD4hep::DDRec + DD4hep::DDCore + ROOT::MathCore + ROOT::MathMore + GSL::gsl ) target_include_directories(${PackageName} PUBLIC @@ -71,7 +85,7 @@ install(TARGETS ${PackageName} EXPORT ${CMAKE_PROJECT_NAME}Targets RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib - PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/@{CMAKE_PROJECT_NAME}" COMPONENT dev + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}" COMPONENT dev ) install(FILES ${scripts} DESTINATION test) @@ -83,3 +97,8 @@ set_test_env(${test_name}) SET(test_name "test_DCHsimpleDigitizerExtendedEdm") ADD_TEST(NAME ${test_name} COMMAND k4run test/runDCHsimpleDigitizerExtendedEdm.py) set_test_env(${test_name}) + +SET(test_name "test_runDCHdigiV2") +ADD_TEST(NAME ${test_name} COMMAND sh +x test_DCHdigi.sh ) +set_test_env(${test_name}) +set_tests_properties(${test_name} PROPERTIES WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}/DCHdigi/test/test_DCHdigi") diff --git a/DCHdigi/README.md b/DCHdigi/README.md new file mode 100644 index 0000000..e3faa87 --- /dev/null +++ b/DCHdigi/README.md @@ -0,0 +1,20 @@ +# Drift chamber (DCH) digitizers + +## DCHdigi_v01 + +* Each simulated hit is transformed into a digitized hit. The digitized hit position is the projection of the simulated hit position onto the sense wire (at the center of the cell) +* Smearing of the digitized hit position along the wire and radially is done according to the input parameter values (`zResolution_mm` and `xyResolution_mm`, respectively) +* The digitized hit adds dNdx information if flag `calculate_dndx` is enabled (default not). This information consist on number of clusters and their size, which are derived from precalculated distributions contained in an input file specified by the parameter `fileDataAlg`. The method and distributions corresponds to the option 3 described in F. Cuna et al, arXiv:2105.07064 +* It requires that the cellID contain the layer and number of cell within the layer (nphi). It does not matter if the segmentation comes from geometrical segmentation by using twisted tubes and hyperboloids (and the cellID is created out of volume IDs), or the segmentation is virtual DD4hep segmentation +* New digitized hit class is used as an EDM4hep data extension, to be integrated into EDM4hep +* Debug histograms are created if `create_debug_histograms` option is enabled (output file name can be given) +* Stand alone test run simulation of the drift chamber based on twisted tubes, and then apply the digitizer. Dedicated directory with all the files needed is given in `DCHdigi/test/test_DCHdigi/` +* Random number generator uses the seeds calculated on an event basis by the UID service, from the podio header information (run/event number) +* This digitizer is meant to be used with `DriftChamber_o1_v02` from k4geo and is expected to work for the upcoming `DriftChamber_o1_v03` + +## DCHsimpleDigitizerExtendedEdm + +* Algorithm for creating digitized drift chamber hits (based on edm4hep::TrackerHit3D) from edm4hep::SimTrackerHit. Resolution along z and xy (distance to the wire) has to be specified. The smearing is applied in the wire reference frame, by means of the placement matrix of the wires +* No cluster counting information is added into the digitized output +* It relies on a dedicated data extension, similar to DCHdigi_v01 +* This digitizer is meant to be used with `DriftChamber_o1_v01` from k4geo. Deprecated. diff --git a/DCHdigi/dataFormatExtension/driftChamberHit.yaml b/DCHdigi/dataFormatExtension/driftChamberHit.yaml index 41f2d7e..0edf251 100644 --- a/DCHdigi/dataFormatExtension/driftChamberHit.yaml +++ b/DCHdigi/dataFormatExtension/driftChamberHit.yaml @@ -3,6 +3,7 @@ schema_version: 1 options: # should getters / setters be prefixed with get / set? getSyntax: True + setSyntax: True # should POD members be exposed with getters/setters in classes that have them as members? exposePODMembers: False includeSubfolder: True @@ -41,3 +42,120 @@ datatypes: OneToOneRelations: - extension::DriftChamberDigi digi // reference to the digitized hit - edm4hep::SimTrackerHit sim // reference to the simulated hit + + extension::DriftChamberDigiV2: + Description: "Digitized hit (before tracking) for Drift Chamber v2 (requires data extension)." + Author: "A. Tolosa-Delgado, B. Francois, CERN" + Members: + - uint64_t cellID // ID of the sensor that created this hit + - int32_t type // type of the raw data hit + - int32_t quality // quality bit flag of the hit + - float time // time of the hit [ns] + - float eDep // energy deposited on the hit [GeV] + - float eDepError // error measured on eDep [GeV] + - edm4hep::Vector3d position // point on the sensitive wire (SW) which is closest to the simhit [mm] + - edm4hep::Vector3d directionSW // direction of SW + - float distanceToWire // distance hit-wire [mm] + - uint32_t nCluster // number of clusters associated to this hit + VectorMembers: + - uint16_t nElectrons // number of electrons for each cluster + + extension::MCRecoDriftChamberDigiV2Association: + Description: "Association between a DriftChamberDigi and the corresponding simulated hit" + Author: "B. Francois, CERN" + Members: + - float weight // weight of this association + OneToOneRelations: + - extension::DriftChamberDigiV2 digi // reference to the digitized hit + - edm4hep::SimTrackerHit sim // reference to the simulated hit + + extension::Track: + Description: "Reconstructed track" + Author: "EDM4hep authors" + Members: + - int32_t type // flagword that defines the type of track.Bits 16-31 are used internally + - float chi2 // Chi^2 of the track fit + - int32_t ndf // number of degrees of freedom of the track fit + - float dEdx // dEdx of the track + - float dEdxError // error of dEdx + - float radiusOfInnermostHit // radius of the innermost hit that has been used in the track fit + VectorMembers: + - int32_t subdetectorHitNumbers // number of hits in particular subdetectors + - edm4hep::TrackState trackStates // track states + - edm4hep::Quantity dxQuantities // different measurements of dx quantities + OneToManyRelations: + - extension::TrackerHit trackerHits // hits that have been used to create this track + - extension::Track tracks // tracks (segments) that have been combined to create this track + +# the following is just because interface looks for e.g. extension/TrackerHit3DCollection.h + extension::TrackerHit3D: + Description: "Tracker hit" + Author: "EDM4hep authors" + Members: + - uint64_t cellID // ID of the sensor that created this hit + - int32_t type // type of raw data hit + - int32_t quality // quality bit flag of the hit + - float time [ns] // time of the hit + - float eDep [GeV] // energy deposited on the hit + - float eDepError [GeV] // error measured on EDep + - edm4hep::Vector3d position [mm] // hit position + - edm4hep::CovMatrix3f covMatrix // covariance matrix of the position (x,y,z) + ExtraCode: + includes: "#include " + declaration: " + /// Get the position covariance matrix value for the two passed dimensions\n + float getCovMatrix(edm4hep::Cartesian dimI, edm4hep::Cartesian dimJ) const { return getCovMatrix().getValue(dimI, dimJ); }\n + " + MutableExtraCode: + includes: "#include " + declaration: " + /// Set the position covariance matrix value for the two passed dimensions\n + void setCovMatrix(float value, edm4hep::Cartesian dimI, edm4hep::Cartesian dimJ) { getCovMatrix().setValue(value, dimI, dimJ); }\n + " + extension::TrackerHitPlane: + Description: "Tracker hit plane" + Author: "EDM4hep authors" + Members: + - uint64_t cellID // ID of the sensor that created this hit + - int32_t type // type of raw data hit + - int32_t quality // quality bit flag of the hit + - float time [ns] // time of the hit + - float eDep [GeV] // energy deposited on the hit + - float eDepError [GeV] // error measured on EDep + - edm4hep::Vector2f u // measurement direction vector, u lies in the x-y plane + - edm4hep::Vector2f v // measurement direction vector, v is along z + - float du // measurement error along the direction + - float dv // measurement error along the direction + - edm4hep::Vector3d position [mm] // hit position + - edm4hep::CovMatrix3f covMatrix // covariance of the position (x,y,z) + ExtraCode: + includes: "#include " + declaration: " + /// Get the position covariance matrix value for the two passed dimensions\n + float getCovMatrix(edm4hep::Cartesian dimI, edm4hep::Cartesian dimJ) const { return getCovMatrix().getValue(dimI, dimJ); }\n + " + MutableExtraCode: + includes: "#include " + declaration: " + /// Set the position covariance matrix value for the two passed dimensions\n + void setCovMatrix(float value, edm4hep::Cartesian dimI, edm4hep::Cartesian dimJ) { getCovMatrix().setValue(value, dimI, dimJ); }\n + " + +# end of "the following is just because interface looks for e.g. extension/TrackerHit3DCollection.h" + +interfaces: + extension::TrackerHit: + Description: "Tracker hit interface class" + Author: "Thomas Madlener, DESY" + Members: + - uint64_t cellID // ID of the sensor that created this hit + - int32_t type // type of the raw data hit + - int32_t quality // quality bit flag of the hit + - float time [ns] // time of the hit + - float eDep [GeV] // energy deposited on the hit + - float eDepError [GeV] // error measured on eDep + - edm4hep::Vector3d position [mm] // hit position + Types: + - extension::TrackerHitPlane + - extension::TrackerHit3D + - extension::DriftChamberDigiV2 diff --git a/DCHdigi/include/AlgData.h b/DCHdigi/include/AlgData.h new file mode 100644 index 0000000..d43f2b3 --- /dev/null +++ b/DCHdigi/include/AlgData.h @@ -0,0 +1,537 @@ +#ifndef ALGDATA_H_INCLUDED +#define ALGDATA_H_INCLUDED + +#include + +#include +#include +#include +#include "Math/InterpolationTypes.h" +#include "Math/Interpolator.h" +#include "TCanvas.h" +#include "TChain.h" +#include "TF1.h" +#include "TFile.h" +#include "TGraph.h" +#include "TTree.h" + +class AlgData { +public: + inline AlgData(); + inline ~AlgData(); + +public: + inline void read_file(TString dataAlg); + inline TF1* read_graph(TString dataAlg, TString cvName, TString fitName); + + inline TF1* get_fit(); + inline TFormula* get_formula(); + inline double get_Fitvalue(double betagamma) const; + inline ROOT::Math::Interpolator* Interpvalue(std::vector bg, std::vector Ydata, int Intertype); + inline void Load_file(TString dataAlg); + inline void Load_interp(); + inline double get_MPVExtra(double betagamma) const; + inline double get_SgmExtra(double betagamma) const; + inline double get_MeanExtra1(double betagamma) const; + inline double get_SgmExtra1(double betagamma) const; + inline double get_SlopeExtra1(double betagamma) const; + inline double get_FracExtra1(double betagamma) const; + inline double get_FfracExtra(double betagamma); + inline double get_maxEx0(double betagamma); + inline double get_maxExSlp(); + inline double get_ExSgmlep(); + inline double get_ExSgmhad(); + inline double get_Ffrac(double betagamma); + inline double get_Fmpv1(double betagamma); + inline double get_Fsgm1(double betagamma); + inline double get_Fmpv2(double betagamma); + inline double get_Fsgm2(double betagamma); + inline double get_ClSzCorrInt(double betagamma); + inline double get_ClSzCorrSlp(double betagamma); + inline std::vector get_ClSzCorrpmean(double betagamma); + inline std::vector get_ClSzCorrpsgm(double betagamma); + + inline std::vector get_ClSzCorrdgmean(double betagamma); + inline std::vector get_ClSzCorrdgsgm(double betagamma); + + inline std::vector get_ClSzCorrdglfrac(double betagamma); + inline std::vector get_ClSzCorrdglmpvl(double betagamma); + inline std::vector get_ClSzCorrdglsgml(double betagamma); + inline std::vector get_ClSzCorrdglmeang(double betagamma); + inline std::vector get_ClSzCorrdglsgmg(double betagamma); + +private: + TString DataAlg; + TChain* data1; + TChain* datalep; + TChain* datahad; + TFile* file; + + double bgT, bgTlep, bgThad, tmpmaxEx0Tot, tmpErrmaxEx0Tot, tmpmaxExSlpTot, tmpErrmaxExSlpTot, tmpExSgmTotlep, + tmpExSgmTothad; + double tmptFfracTot, tmptFerrfracTot, tmptFmpv1Tot, tmptFerrmpv1Tot, tmptFsgm1Tot, tmptFerrsgm1Tot, tmptFmpv2Tot, + tmptFerrmpv2Tot, tmptFsgm2Tot, tmptFerrsgm2Tot; + double tmpCorrIntTot, tmpErrCorrIntTot, tmpCorrSlpTot, tmpErrCorrSlpTot; + + std::vector* tmpCorrmeanSliceTot = nullptr; + std::vector* tmpErrCorrmeanSliceTot = nullptr; + std::vector* tmpCorrsgmSliceTot = nullptr; + std::vector* tmpErrCorrsgmSliceTot = nullptr; + + std::vector* tmpCorrdglfracSliceTot = nullptr; + std::vector* tmpErrCorrdglfracSliceTot = nullptr; + std::vector* tmpCorrdglmeangSliceTot = nullptr; + std::vector* tmpErrCorrdglmeangSliceTot = nullptr; + std::vector* tmpCorrdglsgmgSliceTot = nullptr; + std::vector* tmpErrCorrdglsgmgSliceTot = nullptr; + std::vector* tmpCorrdglmpvSliceTot = nullptr; + std::vector* tmpErrCorrdglmpvSliceTot = nullptr; + std::vector* tmpCorrdglsgmlSliceTot = nullptr; + std::vector* tmpErrCorrdglsgmlSliceTot = nullptr; + + std::vector* tmpCorrdgmeanSliceTot = nullptr; + std::vector* tmpErrCorrdgmeanSliceTot = nullptr; + std::vector* tmpCorrdgsgmSliceTot = nullptr; + std::vector* tmpErrCorrdgsgmSliceTot = nullptr; + + //____________________________________________________// + std::vector bgv, bglep, bghad; + std::vector maxEx0Tot, ErrmaxEx0Tot, maxExSlpTot, ErrmaxExSlpTot, ExSgmTotlep, ExSgmTothad; + std::vector tFfracTot, tFerrfracTot, tFmpv1Tot, tFerrmpv1Tot, tFsgm1Tot, tFerrsgm1Tot, tFmpv2Tot, + tFerrmpv2Tot, tFsgm2Tot, tFerrsgm2Tot; + std::vector CorrmeanSliceTot[20], ErrCorrmeanSliceTot[20], CorrsgmSliceTot[20], ErrCorrsgmSliceTot[20]; + std::vector CorrIntTot, ErrCorrIntTot, CorrSlpTot, ErrCorrSlpTot; + std::vector CorrdglfracSliceTot[20], ErrCorrdglfracSliceTot[20], CorrdglmeangSliceTot[20], + ErrCorrdglmeangSliceTot[20], CorrdglsgmgSliceTot[20], ErrCorrdglsgmgSliceTot[20], CorrdglmpvSliceTot[20], + ErrCorrdglmpvSliceTot[20], CorrdglsgmlSliceTot[20], ErrCorrdglsgmlSliceTot[20]; + std::vector CorrdgmeanSliceTot[20], ErrCorrdgmeanSliceTot[20], CorrdgsgmSliceTot[20], + ErrCorrdgsgmSliceTot[20]; + + //__________________________________________// + TF1* Ft; + TF1* FitMPVExtra; + TF1* FitSgmExtra; + TF1* FitMeanExtra1; + TF1* FitSgmExtra1; + TF1* FitSlopeExtra1; + TF1* FitFracExtra1; + + TFormula* FtFormula; + + std::vector Corrmean; + std::vector Corrsgm; + + std::vector Corrdgmean; + std::vector Corrdgsgm; + + std::vector Corrdglfrac; + std::vector Corrdglmeang; + std::vector Corrdglsgmg; + std::vector Corrdglmpv; + std::vector Corrdglsgml; + + std::vector ClSzCorrpmean; + std::vector ClSzCorrpsgm; + std::vector ClSzCorrdgmean; + std::vector ClSzCorrdgsgm; + std::vector ClSzCorrdglfrac; + std::vector ClSzCorrdglmpvl; + std::vector ClSzCorrdglsgml; + std::vector ClSzCorrdglmeang; + std::vector ClSzCorrdglsgmg; + + double maxExSlp; + double ExSgmlep; + double ExSgmhad; + inline void Calc_maxExSlp(); + inline void Calc_ExSgmlep(); + inline void Calc_ExSgmhad(); + + //___________________________________________// + + std::map itpm; + ROOT::Math::Interpolator* itp1; +}; + +//_______________________________________________// + +AlgData::AlgData() { + DataAlg = ""; + data1 = nullptr; + datalep = nullptr; + datahad = nullptr; + file = nullptr; + + tmpCorrmeanSliceTot = nullptr; + tmpErrCorrmeanSliceTot = nullptr; + tmpCorrsgmSliceTot = nullptr; + tmpErrCorrsgmSliceTot = nullptr; + + tmpCorrdglfracSliceTot = nullptr; + tmpErrCorrdglfracSliceTot = nullptr; + tmpCorrdglmeangSliceTot = nullptr; + tmpErrCorrdglmeangSliceTot = nullptr; + tmpCorrdglsgmgSliceTot = nullptr; + tmpErrCorrdglsgmgSliceTot = nullptr; + tmpCorrdglmpvSliceTot = nullptr; + tmpErrCorrdglmpvSliceTot = nullptr; + tmpCorrdglsgmlSliceTot = nullptr; + tmpErrCorrdglsgmlSliceTot = nullptr; + + tmpCorrdgmeanSliceTot = nullptr; + tmpErrCorrdgmeanSliceTot = nullptr; + tmpCorrdgsgmSliceTot = nullptr; + tmpErrCorrdgsgmSliceTot = nullptr; + + Ft = nullptr; + FtFormula = nullptr; + FitMPVExtra = nullptr; + FitSgmExtra = nullptr; + FitMeanExtra1 = nullptr; + FitSgmExtra1 = nullptr; + FitSlopeExtra1 = nullptr; + FitFracExtra1 = nullptr; +} + +AlgData::~AlgData() { + if (file->IsOpen()) + file->Close(); +} + +void AlgData::read_file(TString dataAlg) { + DataAlg = dataAlg; + + std::cout << " ---reading of---" << DataAlg << std::endl; + + data1 = new TChain("DataAlg"); + datalep = new TChain("DataAlglep"); + datahad = new TChain("DataAlghad"); + + if (dataAlg.Contains(".root")) { + data1->Add(dataAlg.Data()); + datalep->Add(dataAlg.Data()); + datahad->Add(dataAlg.Data()); + } + + data1->SetBranchAddress("bgT", &bgT); + data1->SetBranchAddress("tmpmaxEx0Tot", &tmpmaxEx0Tot); + data1->SetBranchAddress("tmpErrmaxEx0Tot", &tmpErrmaxEx0Tot); + data1->SetBranchAddress("tmpmaxExSlpTot", &tmpmaxExSlpTot); + data1->SetBranchAddress("tmpErrmaxExSlpTot", &tmpErrmaxExSlpTot); + data1->SetBranchAddress("tmptFfracTot", &tmptFfracTot); + data1->SetBranchAddress("tFerrfracTot", &tmptFerrfracTot); + data1->SetBranchAddress("tmptFerrfracTot", &tmptFmpv1Tot); + data1->SetBranchAddress("tmptFerrmpv1Tot", &tmptFerrmpv1Tot); + data1->SetBranchAddress("tmptFerrmpv1Tot", &tmptFerrmpv1Tot); + data1->SetBranchAddress("tmptFsgm1Tot", &tmptFsgm1Tot); + data1->SetBranchAddress("tmptFerrsgm1Tot", &tmptFerrsgm1Tot); + data1->SetBranchAddress("tmptFmpv2Tot", &tmptFmpv2Tot); + data1->SetBranchAddress("tmptFerrmpv2Tot", &tmptFerrmpv2Tot); + data1->SetBranchAddress("tmptFsgm2Tot", &tmptFsgm2Tot); + data1->SetBranchAddress("tmptFerrsgm2Tot", &tmptFerrsgm2Tot); + + data1->SetBranchAddress("tmpCorrIntTot", &tmpCorrIntTot); + data1->SetBranchAddress("tmpErrCorrIntTot", &tmpErrCorrIntTot); + data1->SetBranchAddress("tmpCorrSlpTot", &tmpCorrSlpTot); + data1->SetBranchAddress("tmpErrCorrSlpTot", &tmpErrCorrSlpTot); + + data1->SetBranchAddress("tmpCorrmeanSliceTot", &tmpCorrmeanSliceTot); + data1->SetBranchAddress("tmpErrCorrmeanSliceTot", &tmpErrCorrmeanSliceTot); + data1->SetBranchAddress("tmpCorrsgmSliceTot", &tmpCorrsgmSliceTot); + data1->SetBranchAddress("tmpErrCorrsgmSliceTot", &tmpErrCorrsgmSliceTot); + + data1->SetBranchAddress("tmpCorrdgmeanSliceTot", &tmpCorrdgmeanSliceTot); + data1->SetBranchAddress("tmpErrCorrdgmeanSliceTot", &tmpErrCorrdgmeanSliceTot); + data1->SetBranchAddress("tmpCorrdgsgmSliceTot", &tmpCorrdgsgmSliceTot); + data1->SetBranchAddress("tmpErrCorrdgsgmSliceTot", &tmpErrCorrdgsgmSliceTot); + + data1->SetBranchAddress("tmpCorrdglfracSliceTot", &tmpCorrdglfracSliceTot); + data1->SetBranchAddress("tmpErrCorrdglfracSliceTot", &tmpErrCorrdglfracSliceTot); + data1->SetBranchAddress("tmpCorrdglmeangSliceTot", &tmpCorrdglmeangSliceTot); + data1->SetBranchAddress("tmpErrCorrdglmeangSliceTot", &tmpErrCorrdglmeangSliceTot); + data1->SetBranchAddress("tmpCorrdglsgmgSliceTot", &tmpCorrdglsgmgSliceTot); + data1->SetBranchAddress("tmpErrCorrdglsgmgSliceTot", &tmpErrCorrdglsgmgSliceTot); + data1->SetBranchAddress("tmpCorrdglmpvSliceTot", &tmpCorrdglmpvSliceTot); + data1->SetBranchAddress("tmpErrCorrdglmpvSliceTot", &tmpErrCorrdglmpvSliceTot); + data1->SetBranchAddress("tmpCorrdglsgmlSliceTot", &tmpCorrdglsgmlSliceTot); + data1->SetBranchAddress("tmpErrCorrdglsgmlSliceTot", &tmpErrCorrdglsgmlSliceTot); + + datalep->SetBranchAddress("bgTlep", &bgTlep); + datalep->SetBranchAddress("tmpExSgmTotlep", &tmpExSgmTotlep); + + datahad->SetBranchAddress("bgThad", &bgThad); + datahad->SetBranchAddress("tmpExSgmTothad", &tmpExSgmTothad); + + for (int i = 0; i < data1->GetEntries(); i++) { + data1->GetEntry(i); + bgv.push_back(bgT); + // std::cout<size(); ++j) { + CorrmeanSliceTot[j].push_back(tmpCorrmeanSliceTot->at(j)); + ErrCorrmeanSliceTot[j].push_back(tmpErrCorrmeanSliceTot->at(j)); + CorrsgmSliceTot[j].push_back(tmpCorrsgmSliceTot->at(j)); + ErrCorrsgmSliceTot[j].push_back(tmpErrCorrsgmSliceTot->at(j)); + } + + for (unsigned int n = 0; n < tmpCorrdgmeanSliceTot->size(); ++n) { + CorrdgmeanSliceTot[n].push_back(tmpCorrdgmeanSliceTot->at(n)); + ErrCorrdgmeanSliceTot[n].push_back(tmpErrCorrdgmeanSliceTot->at(n)); + CorrdgsgmSliceTot[n].push_back(tmpCorrdgsgmSliceTot->at(n)); + ErrCorrdgsgmSliceTot[n].push_back(tmpErrCorrdgsgmSliceTot->at(n)); + } + + for (unsigned int k = 0; k < tmpCorrdglfracSliceTot->size(); ++k) { + CorrdglfracSliceTot[k].push_back(tmpCorrdglfracSliceTot->at(k)); + ErrCorrdglfracSliceTot[k].push_back(tmpErrCorrdglfracSliceTot->at(k)); + CorrdglmeangSliceTot[k].push_back(tmpCorrdglmeangSliceTot->at(k)); + ErrCorrdglmeangSliceTot[k].push_back(tmpErrCorrdglmeangSliceTot->at(k)); + CorrdglsgmgSliceTot[k].push_back(tmpCorrdglsgmgSliceTot->at(k)); + ErrCorrdglsgmgSliceTot[k].push_back(tmpErrCorrdglsgmgSliceTot->at(k)); + CorrdglmpvSliceTot[k].push_back(tmpCorrdglmpvSliceTot->at(k)); + ErrCorrdglmpvSliceTot[k].push_back(tmpErrCorrdglmpvSliceTot->at(k)); + CorrdglsgmlSliceTot[k].push_back(tmpCorrdglsgmlSliceTot->at(k)); + ErrCorrdglsgmlSliceTot[k].push_back(tmpErrCorrdglsgmlSliceTot->at(k)); + // std::cout<at(k)<GetEntries(); i++) { + datalep->GetEntry(i); + bglep.push_back(bgTlep); + ExSgmTotlep.push_back(tmpExSgmTotlep); + // std::cout<< " bgTlep "<GetEntries(); i++) { + datahad->GetEntry(i); + bghad.push_back(bgThad); + ExSgmTothad.push_back(tmpExSgmTothad); + } +} + +TF1* AlgData::read_graph(TString dataAlg, TString cvName, TString fitName) { + DataAlg = dataAlg; + file = TFile::Open(dataAlg.Data(), "read"); + TCanvas* cv = (TCanvas*)file->Get(cvName.Data()); + TGraph* gr = (TGraph*)cv->GetListOfPrimitives()->FindObject("Graph"); + TF1* ft = (TF1*)gr->GetListOfFunctions()->FindObject(fitName.Data()); + TFormula* ftFormula = (TFormula*)ft->GetFormula(); + FtFormula = ftFormula; + return ft; +} + +void AlgData::Load_file(TString dataAlg) { + read_file(dataAlg.Data()); + FitMPVExtra = read_graph(dataAlg.Data(), "cMPVExtrabgTot", "fit_MPV"); + FitSgmExtra = read_graph(dataAlg.Data(), "cSgmExtrabgTot", "fit_sgmEx"); + FitMeanExtra1 = read_graph(dataAlg.Data(), "cMeanExtra1bgTot", "expeff"); + FitSgmExtra1 = read_graph(dataAlg.Data(), "cSgmExtra1bgTot", "expeffNeg"); + FitSlopeExtra1 = read_graph(dataAlg.Data(), "cSlopeExtra1bgTot", "fit_slp"); + FitFracExtra1 = read_graph(dataAlg.Data(), "cfracbgTot", "fit_frEx"); + Calc_ExSgmhad(); + Calc_ExSgmlep(); + Calc_maxExSlp(); +} + +TF1* AlgData::get_fit() { return Ft; } + +TFormula* AlgData::get_formula() { return FtFormula; } + +double AlgData::get_MPVExtra(double betagamma) const { return FitMPVExtra->Eval(betagamma); } + +double AlgData::get_SgmExtra(double betagamma) const { return FitSgmExtra->Eval(betagamma); } + +double AlgData::get_MeanExtra1(double betagamma) const { return FitMeanExtra1->Eval(betagamma); } + +double AlgData::get_SgmExtra1(double betagamma) const { return FitSgmExtra1->Eval(betagamma); } + +double AlgData::get_FracExtra1(double betagamma) const { return FitFracExtra1->Eval(betagamma); } + +double AlgData::get_SlopeExtra1(double betagamma) const { return FitSlopeExtra1->Eval(betagamma); } + +ROOT::Math::Interpolator* AlgData::Interpvalue(std::vector bg, std::vector Ydata, int Intertype) { + //tipo di interpolazione 0=linear;1=pol;2=cspline;4=akima + return itp1 = new ROOT::Math::Interpolator(bg, Ydata, (ROOT::Math::Interpolation::Type)Intertype); +} + +void AlgData::Load_interp() { + itpm["maxEx0Tot"] = Interpvalue(bgv, maxEx0Tot, 4); + itpm["tFfracTot"] = Interpvalue(bgv, tFfracTot, 4); + itpm["tFmpv1Tot"] = Interpvalue(bgv, tFmpv1Tot, 4); + itpm["tFsgm1Tot"] = Interpvalue(bgv, tFsgm1Tot, 4); + itpm["tFmpv2Tot"] = Interpvalue(bgv, tFmpv2Tot, 4); + itpm["tFsgm2Tot"] = Interpvalue(bgv, tFsgm2Tot, 4); + + itpm["CorrIntTot"] = Interpvalue(bgv, CorrIntTot, 4); + itpm["CorrSlpTot"] = Interpvalue(bgv, CorrSlpTot, 4); + for (unsigned int i = 0; i < tmpCorrmeanSliceTot->size(); ++i) { + TString nameCm = "CorrmeanSliceTot"; + TString nameCs = "CorrsgmSliceTot"; + Corrmean.push_back(nameCm + i); + Corrsgm.push_back(nameCs + i); + itpm[Corrmean.at(i)] = Interpvalue(bgv, CorrmeanSliceTot[i], 4); + itpm[Corrsgm.at(i)] = Interpvalue(bgv, CorrsgmSliceTot[i], 4); + } + + for (unsigned int j = 0; j < tmpCorrdglfracSliceTot->size(); ++j) { + TString nameCdf = "Corrdglfrac"; + TString nameCdm = "Corrdglmeang"; + TString nameCds = "Corrdglsgmg"; + TString nameCdmpv = "Corrdglmpv"; + TString nameCdsgl = "Corrdglsgmgl"; + Corrdglfrac.push_back(nameCdf + j); + Corrdglmeang.push_back(nameCdm + j); + Corrdglsgmg.push_back(nameCds + j); + Corrdglmpv.push_back(nameCdmpv + j); + Corrdglsgml.push_back(nameCdsgl + j); + itpm[Corrdglfrac.at(j)] = Interpvalue(bgv, CorrdglfracSliceTot[j], 4); + itpm[Corrdglmeang.at(j)] = Interpvalue(bgv, CorrdglmeangSliceTot[j], 4); + itpm[Corrdglsgmg.at(j)] = Interpvalue(bgv, CorrdglsgmgSliceTot[j], 4); + itpm[Corrdglmpv.at(j)] = Interpvalue(bgv, CorrdglmpvSliceTot[j], 4); + itpm[Corrdglsgml.at(j)] = Interpvalue(bgv, CorrdglsgmlSliceTot[j], 4); + } + + for (unsigned int n = 0; n < tmpCorrdgmeanSliceTot->size(); ++n) { + TString nameCdmg = "Corrdgmean"; + TString nameCdsgg = "Corrdgsgm"; + Corrdgmean.push_back(nameCdmg + n); + Corrdgsgm.push_back(nameCdsgg + n); + itpm[Corrdgmean.at(n)] = Interpvalue(bgv, CorrdgmeanSliceTot[n], 4); + itpm[Corrdgsgm.at(n)] = Interpvalue(bgv, CorrdgsgmSliceTot[n], 4); + } +} + +double AlgData::get_maxEx0(double betagamma) { return itpm["maxEx0Tot"]->Eval(betagamma); } +double AlgData::get_Ffrac(double betagamma) { return itpm["tFfracTot"]->Eval(betagamma); } +double AlgData::get_Fmpv1(double betagamma) { return itpm["tFmpv1Tot"]->Eval(betagamma); } +double AlgData::get_Fsgm1(double betagamma) { return itpm["tFsgm1Tot"]->Eval(betagamma); } +double AlgData::get_Fmpv2(double betagamma) { return itpm["tFmpv2Tot"]->Eval(betagamma); } +double AlgData::get_Fsgm2(double betagamma) { return itpm["tFsgm2Tot"]->Eval(betagamma); } +double AlgData::get_ClSzCorrInt(double betagamma) { return itpm["CorrIntTot"]->Eval(betagamma); } +double AlgData::get_ClSzCorrSlp(double betagamma) { return itpm["CorrSlpTot"]->Eval(betagamma); } + +std::vector AlgData::get_ClSzCorrpmean(double betagamma) { + ClSzCorrpmean.clear(); + // std::cout<<"size "<size()<size(); ++i) { + ClSzCorrpmean.push_back(itpm[Corrmean.at(i)]->Eval(betagamma)); + } + + return ClSzCorrpmean; +} + +std::vector AlgData::get_ClSzCorrpsgm(double betagamma) { + ClSzCorrpsgm.clear(); + for (unsigned int i = 0; i < tmpCorrmeanSliceTot->size(); ++i) { + ClSzCorrpsgm.push_back(itpm[Corrsgm.at(i)]->Eval(betagamma)); + } + return ClSzCorrpsgm; +} + +std::vector AlgData::get_ClSzCorrdglfrac(double betagamma) { + ClSzCorrdglfrac.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglfrac.push_back(itpm[Corrdglfrac.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglfrac; +} + +std::vector AlgData::get_ClSzCorrdglmeang(double betagamma) { + ClSzCorrdglmeang.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglmeang.push_back(itpm[Corrdglmeang.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglmeang; +} + +std::vector AlgData::get_ClSzCorrdglsgmg(double betagamma) { + ClSzCorrdglsgmg.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglsgmg.push_back(itpm[Corrdglsgmg.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglsgmg; +} + +std::vector AlgData::get_ClSzCorrdglmpvl(double betagamma) { + ClSzCorrdglmpvl.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglmpvl.push_back(itpm[Corrdglmpv.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglmpvl; +} + +std::vector AlgData::get_ClSzCorrdglsgml(double betagamma) { + ClSzCorrdglsgml.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglsgml.push_back(itpm[Corrdglsgml.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglsgml; +} + +std::vector AlgData::get_ClSzCorrdgmean(double betagamma) { + ClSzCorrdgmean.clear(); + for (unsigned int i = 0; i < tmpCorrdgmeanSliceTot->size(); ++i) { + ClSzCorrdgmean.push_back(itpm[Corrdgmean.at(i)]->Eval(betagamma)); + } + return ClSzCorrdgmean; +} + +std::vector AlgData::get_ClSzCorrdgsgm(double betagamma) { + ClSzCorrdgsgm.clear(); + for (unsigned int i = 0; i < tmpCorrdgmeanSliceTot->size(); ++i) { + ClSzCorrdgsgm.push_back(itpm[Corrdgsgm.at(i)]->Eval(betagamma)); + } + return ClSzCorrdgsgm; +} + +double AlgData::get_maxExSlp() { return maxExSlp; } + +void AlgData::Calc_maxExSlp() { + double sum = 0.0; + for (unsigned int i = 0; i < maxExSlpTot.size(); ++i) { + sum += maxExSlpTot[i]; + } + maxExSlp = sum / (double)maxExSlpTot.size(); +} + +double AlgData::get_ExSgmlep() { return ExSgmlep; } + +void AlgData::Calc_ExSgmlep() { + double sum = 0.0; + for (unsigned int i = 0; i < ExSgmTotlep.size(); ++i) { + sum += ExSgmTotlep[i]; + } + ExSgmlep = sum / (double)ExSgmTotlep.size(); +} + +double AlgData::get_ExSgmhad() { return ExSgmhad; } + +void AlgData::Calc_ExSgmhad() { + double sum = 0.0; + for (unsigned int i = 0; i < ExSgmTothad.size(); ++i) { + sum += ExSgmTothad[i]; + } + ExSgmhad = sum / (double)ExSgmTothad.size(); +} + +#endif // ALGDATA_H_INCLUDED diff --git a/DCHdigi/include/DCHdigi_v01.h b/DCHdigi/include/DCHdigi_v01.h new file mode 100644 index 0000000..b52bfb9 --- /dev/null +++ b/DCHdigi/include/DCHdigi_v01.h @@ -0,0 +1,215 @@ +/** ======= DCHdigi_v01 ========== + * Gaudi Algorithm for DCH digitization + * + * + * @author Alvaro Tolosa-Delgado, Brieuc Francois + * @date 2024-08 + * + *

Input collections and prerequisites

+ * Processor requires a collection of SimTrackerHits
+ * This code uses DD4hep length natural unit (cm), but EDM4hep data is (usually) in mm. Please be careful with units.
+ *

Output

+ * Processor produces collection of Digitized hits of Drift Chamber v2
+ * @param DCH_simhits The name of input collection, type edm4hep::SimTrackerHitCollection
+ * (default name empty)
+ * @param DCH_DigiCollection The name of out collection, type extension::DriftChamberDigiV2Collection
+ * (default name DCH_DigiCollection)
+ * @param DCH_name DCH subdetector name
+ * (default value DCH_v2)
+ * @param calculate_dndx Optional flag to calcualte dNdx information
+ * (default value false)
+ * @param fileDataAlg File needed for calculating cluster count and size
+ * (default value /eos/.../DataAlgFORGEANT.root)
+ * @param zResolution_mm Resolution (sigma for gaussian smearing) along the sense wire, in mm
+ * (default value 1 mm)
+ * @param xyResolution_mm Resolution (sigma for gaussian smearing) perpendicular the sense wire, in mm
+ * (default value 0.1 mm)
+ * @param create_debug_histograms Optional flag to create debug histograms
+ * (default value false)
+ * @param GeoSvcName Geometry service name
+ * (default value GeoSvc)
+ * @param uidSvcName The name of the UniqueIDGenSvc instance, used to create seed for each event/run, ensuring reproducibility.
+ * (default value uidSvc)
+ *
+ */ + +#ifndef DCHDIGI_V01_H +#define DCHDIGI_V01_H + +// Gaudi Transformer baseclass headers +#include "Gaudi/Property.h" +#include "k4FWCore/Transformer.h" + +// Gaudi services +#include "k4Interface/IGeoSvc.h" +#include "k4Interface/IUniqueIDGenSvc.h" + +// EDM4HEP +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/ParticleIDData.h" +#include "edm4hep/SimTrackerHitCollection.h" + +// EDM4HEP extension +#include "extension/DriftChamberDigiV2Collection.h" +#include "extension/MCRecoDriftChamberDigiV2AssociationCollection.h" + +// DD4hep +#include "DD4hep/Detector.h" // for dd4hep::VolumeManager +#include "DDSegmentation/BitFieldCoder.h" + +// STL +#include +#include + +// data extension for detector DCH_v2 +#include "DDRec/DCH_info.h" + +// ROOT headers +#include "TFile.h" +#include "TH1D.h" +#include "TRandom3.h" +#include "TVector3.h" + +// Class developed by Walaa for the CLS +#include "AlgData.h" + +/// constant to convert from mm (EDM4hep) to DD4hep (cm) + +struct DCHdigi_v01 final + : k4FWCore::MultiTransformer< + std::tuple( + const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&)> { + DCHdigi_v01(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize() override; + StatusCode finalize() override; + + std::tuple + operator()(const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&) const override; + +private: + /// conversion factor mm to cm, static to the class to avoid clash with DD4hep + static constexpr double MM_TO_CM = 0.1; + + //------------------------------------------------------------------ + // machinery for geometry + + /// Geometry service name + Gaudi::Property m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"}; + Gaudi::Property m_uidSvcName{this, "uidSvcName", "uidSvc", "The name of the UniqueIDGenSvc instance"}; + + /// Detector name + Gaudi::Property m_DCH_name{this, "DCH_name", "DCH_v2", "Name of the Drift Chamber detector"}; + + /// Pointer to the geometry service + SmartIF m_geoSvc; + + /// Decoder for the cellID + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + + /// Pointer to drift chamber data extension + dd4hep::rec::DCH_info* dch_data = {nullptr}; + + //------------------------------------------------------------------ + // machinery for smearing the position + + /// along the sense wire position resolution in mm + Gaudi::Property m_z_resolution{ + this, "zResolution_mm", 1.0, + "Spatial resolution in the z direction (from reading out the wires at both sides) in mm. Default 1 mm."}; + /// xy resolution in mm + Gaudi::Property m_xy_resolution{this, "xyResolution_mm", 0.1, + "Spatial resolution in the xy direction in mm. Default 0.1 mm."}; + + /// create seed using the uid + SmartIF m_uidSvc; + /// use thread local engine from C++ standard + inline static thread_local std::mt19937_64 m_engine; + void PrepareRandomEngine(const edm4hep::EventHeaderCollection& headers) const; + + // Operator std::normal_distribution::operator()(Generator& g) is a non-const member function and thus cannot be called for a constant object. So we defined the distribution as mutable. + // Gaussian random number generator used for the smearing of the z position, in cm! + mutable std::normal_distribution m_gauss_z_cm; + // Gaussian random number generator used for the smearing of the xy position, in cm! + mutable std::normal_distribution m_gauss_xy_cm; + + /// members with internal state (such as random engines) must be defined thread local + inline static thread_local TRandom3 myRandom; + //------------------------------------------------------------------ + // ancillary functions + + bool IsFileGood(std::string& ifilename) const { return std::ifstream(ifilename).good(); } + + /// Print algorithm configuration + void PrintConfiguration(std::ostream& io); + + /// Send error message to logger and then throw exception + void ThrowException(std::string s) const; + + int CalculateLayerFromCellID(dd4hep::DDSegmentation::CellID id) const { + //return m_decoder->get(id, "layer") + dch_data->nlayersPerSuperlayer * m_decoder->get(id, "superlayer") + 1; + return dch_data->CalculateILayerFromCellIDFields( m_decoder->get(id, "layer"), m_decoder->get(id, "superlayer") ); + } + + int CalculateNphiFromCellID(dd4hep::DDSegmentation::CellID id) const { return m_decoder->get(id, "nphi"); } + + TVector3 Convert_EDM4hepVector_to_TVector3(const edm4hep::Vector3d& v, double scale) const { + return {v[0] * scale, v[1] * scale, v[2] * scale}; + }; + edm4hep::Vector3d Convert_TVector3_to_EDM4hepVector(const TVector3& v, double scale) const { + return {v.x() * scale, v.y() * scale, v.z() * scale}; + }; + + //------------------------------------------------------------------ + // cluster calculation, developed by Walaa + + /// Flag to create to calculate cluster counting information + Gaudi::Property m_calculate_dndx{this, "calculate_dndx", false, + "Calculate number of clusters and electron per cluster"}; + + /// file with distributions to be sampled + Gaudi::Property m_fileDataAlg{ + this, "fileDataAlg", "/eos/project/f/fccsw-web/www/filesForSimDigiReco/IDEA/DataAlgFORGEANT.root", + "ROOT file with cluster size distributions"}; + + /// pointer to wrapper class, which contains the cluster size and number distributions + AlgData* flData; + + /// code developed by Walaa for calculating number of clusters and cluster size of each one + std::pair> CalculateClusters(const edm4hep::SimTrackerHit& input_sim_hit) const; + + bool IsParticleCreatedInsideDriftChamber(const edm4hep::MCParticle &) const ; + + + //------------------------------------------------------------------ + // debug information + + /// Flag to create output file with debug histgrams + Gaudi::Property m_create_debug_histos{this, "create_debug_histograms", false, + "Create output file with histograms for debugging"}; + + /// name for the file that will contain the histograms for debugging + Gaudi::Property m_out_debug_filename{this, "out_debug_filename", "dch_digi_alg_debug.root", + "name for the file that will contain the histograms for debugging"}; + /// histogram to store distance from sim hit position to the sense wire + TH1D* hDpw; + + /// histogram to store distance from digi-hit to the wire. Should be zero because digi-hit central position lies on the wire. + /// This histogram is a consistency check, because the function used to calculate the distance to the wire is different from + /// the function used to calculate the digi-hit central position from a sim-hit position + TH1D* hDww; + + /// histogram to store smearing along the wire + TH1D* hSz; + + /// histogram to store smearing perpendicular the wire + TH1D* hSxy; + + /// Create ROOT file for debug histograms + /// Does not change ROOT directory + void Create_outputROOTfile_for_debugHistograms(); +}; + +DECLARE_COMPONENT(DCHdigi_v01); + +#endif diff --git a/DCHdigi/src/DCHdigi_v01.cpp b/DCHdigi/src/DCHdigi_v01.cpp new file mode 100644 index 0000000..a981046 --- /dev/null +++ b/DCHdigi/src/DCHdigi_v01.cpp @@ -0,0 +1,525 @@ +#include "DCHdigi_v01.h" + +// STL +#include +#include + +#include "extension/MutableDriftChamberDigiV2.h" + +/////////////////////////////////////////////////////////////////////////////////////// +////////////////////// DCHdigi_v01 constructor //////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// +// -- KeyValues("name of the variable that holds the name of the collection exposed in the python steering file", {"default name for the collection"}), +DCHdigi_v01::DCHdigi_v01(const std::string& name, ISvcLocator* svcLoc) + : MultiTransformer(name, svcLoc, + { + KeyValues("DCH_simhits", {""}), + KeyValues("HeaderName", {"EventHeader"}), + }, + {KeyValues("DCH_DigiCollection", {"DCH_DigiCollection"}), + KeyValues("DCH_DigiSimAssociationCollection", {"DCH_DigiSimAssociationCollection"})}) { + m_geoSvc = serviceLocator()->service(m_geoSvcName); + m_uidSvc = serviceLocator()->service(m_uidSvcName); +} + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// initialize //////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// +StatusCode DCHdigi_v01::initialize() { + if (!m_uidSvc) + ThrowException("Unable to get UniqueIDGenSvc"); + + if( 0 > m_z_resolution.value() ) + ThrowException("Z resolution input value can not be negative!"); + m_gauss_z_cm = std::normal_distribution(0., m_z_resolution.value() * MM_TO_CM); + + if( 0 > m_xy_resolution.value() ) + ThrowException("Radial (XY) resolution input value can not be negative!"); + m_gauss_xy_cm = std::normal_distribution(0., m_xy_resolution.value() * MM_TO_CM); + + //----------------- + // Retrieve the subdetector + std::string DCH_name(m_DCH_name.value()); + if (0 == m_geoSvc->getDetector()->detectors().count(DCH_name)) { + ThrowException("Detector <<" + DCH_name + ">> does not exist."); + } + + dd4hep::DetElement DCH_DE = m_geoSvc->getDetector()->detectors().at(DCH_name); + + /////////////////////////////////////////////////////////////////////////////////// + /////////////////////////// retrieve data extension ////////////////////////// + /////////////////////////////////////////////////////////////////////////////////// + this->dch_data = DCH_DE.extension(); + if (not dch_data->IsValid()) + ThrowException("No valid data extension was found for detector <<" + DCH_name + ">>."); + + /////////////////////////////////////////////////////////////////////////////////// + + //----------------- + // Retrieve the readout associated with the detector element (subdetector) + dd4hep::SensitiveDetector dch_sd = m_geoSvc->getDetector()->sensitiveDetector(DCH_name); + if (not dch_sd.isValid()) + ThrowException("No valid Sensitive Detector was found for detector <<" + DCH_name + ">>."); + + dd4hep::Readout dch_readout = dch_sd.readout(); + // set the cellID decoder + m_decoder = dch_readout.idSpec().decoder(); + + /////////////////////////////////////////////////////////////////////////////////// + ////////////////// initialize Walaa's code for Cluster counting ///////////////// + /////////////////////////////////////////////////////////////////////////////////// + debug() << Form("Opening %s...", m_fileDataAlg.value().c_str()) << endmsg; + if (not IsFileGood(m_fileDataAlg.value())) + ThrowException("File <<" + m_fileDataAlg.value() + ">> not found."); + flData = new AlgData(); + flData->Load_file(m_fileDataAlg.value().c_str()); + flData->Load_interp(); + debug() << Form("Opening %s... Done", m_fileDataAlg.value().c_str()) << endmsg; + /////////////////////////////////////////////////////////////////////////////////// + + std::stringstream ss; + PrintConfiguration(ss); + info() << ss.str().c_str() << endmsg; + if (m_create_debug_histos.value()) { + hDpw = new TH1D("hDpw", "Distance from sim-hit to the wire, in cm", 100, 0, 1); + hDpw->SetDirectory(0); + hDww = new TH1D("hDww", "Distance from digi-hit to the wire, in cm. Should be zero because digi-hit central position lies on the wire", 100, 0, 1); + hDww->SetDirectory(0); + hSz = new TH1D("hSz", "Smearing along the wire, in cm", 100, 0, 5 * m_z_resolution.value()); + hSz->SetDirectory(0); + hSxy = new TH1D("hSxy", "Smearing perpendicular the wire, in cm", 100, 0, 5 * m_xy_resolution.value()); + hSxy->SetDirectory(0); + } + return StatusCode::SUCCESS; +} + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// operator() //////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// +std::tuple +DCHdigi_v01::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits, + const edm4hep::EventHeaderCollection& headers) const { + // initialize seed for random engine + this->PrepareRandomEngine(headers); + + debug() << "Input Sim Hit collection size: " << input_sim_hits.size() << endmsg; + + // Create the collections we are going to return + extension::DriftChamberDigiV2Collection output_digi_hits; + extension::MCRecoDriftChamberDigiV2AssociationCollection output_digi_sim_association; + + //loop over hit collection + for (const auto& input_sim_hit : input_sim_hits) { + dd4hep::DDSegmentation::CellID cellid = input_sim_hit.getCellID(); + int ilayer = this->CalculateLayerFromCellID(cellid); + int nphi = this->CalculateNphiFromCellID(cellid); + auto hit_position = Convert_EDM4hepVector_to_TVector3(input_sim_hit.getPosition(), MM_TO_CM); + + // ------------------------------------------------------------------------- + // calculate hit position projection into the wire + TVector3 hit_to_wire_vector = this->dch_data->Calculate_hitpos_to_wire_vector(ilayer, nphi, hit_position); + TVector3 hit_projection_on_the_wire = hit_position + hit_to_wire_vector; + if (m_create_debug_histos.value()) { + double distance_hit_wire = hit_to_wire_vector.Mag(); + hDpw->Fill(distance_hit_wire); + } + TVector3 wire_direction_ez = this->dch_data->Calculate_wire_vector_ez(ilayer, nphi); + + // ------------------------------------------------------------------------- + // smear the position + + // smear position along the wire + double smearing_z = m_gauss_z_cm(m_engine); + if (m_create_debug_histos.value()) + hSz->Fill(smearing_z); + + hit_projection_on_the_wire += smearing_z * (wire_direction_ez.Unit()); + if (m_create_debug_histos.value()) { + // the distance from the hit projection and the wire should be zero + TVector3 dummy_vector = this->dch_data->Calculate_hitpos_to_wire_vector(ilayer, nphi, hit_projection_on_the_wire); + hDww->Fill(dummy_vector.Mag()); + } + + // smear position perpendicular to the wire + double smearing_xy = m_gauss_xy_cm(m_engine); + if (m_create_debug_histos.value()) + hSxy->Fill(smearing_xy); + float distanceToWire_real = hit_to_wire_vector.Mag(); + + // protect against negative values + float distanceToWire_smeared = std::max(0.0, distanceToWire_real + smearing_xy); + + std::int32_t type = 0; + std::int32_t quality = 0; + float eDepError = 0; + // length units back to mm + auto positionSW = Convert_TVector3_to_EDM4hepVector(hit_projection_on_the_wire, 1. / MM_TO_CM); + auto directionSW = Convert_TVector3_to_EDM4hepVector(wire_direction_ez, 1. / MM_TO_CM); + float distanceToWire = distanceToWire_smeared / MM_TO_CM; + + extension::MutableDriftChamberDigiV2 oDCHdigihit; + oDCHdigihit.setCellID(input_sim_hit.getCellID()); + oDCHdigihit.setType(type); + oDCHdigihit.setQuality(quality); + oDCHdigihit.setTime(input_sim_hit.getTime()); + oDCHdigihit.setEDep(input_sim_hit.getEDep()); + oDCHdigihit.setEDepError(eDepError); + oDCHdigihit.setPosition(positionSW); + oDCHdigihit.setDirectionSW(directionSW); + oDCHdigihit.setDistanceToWire(distanceToWire); + // For the sake of speed, let the dNdx calculation be optional + if( m_calculate_dndx.value() ) + { + auto [nCluster, nElectrons_v] = CalculateClusters(input_sim_hit); + oDCHdigihit.setNCluster(nCluster); + // to return the total number of electrons within the step, do the following: + // int nElectronsTotal = std::accumulate( nElectrons_v.begin(), nElectrons_v.end(), 0); + // oDCHdigihit.setNElectronsTotal(nElectronsTotal); + // to copy the vector of each cluster size to the EDM4hep data extension, do the following: + for( auto ne : nElectrons_v ) + oDCHdigihit.addToNElectrons(ne); + } + + output_digi_hits.push_back(oDCHdigihit); + + extension::MutableMCRecoDriftChamberDigiV2Association oDCHsimdigi_association; + oDCHsimdigi_association.setDigi(oDCHdigihit); + oDCHsimdigi_association.setSim(input_sim_hit); + output_digi_sim_association.push_back(oDCHsimdigi_association); + + } // end loop over hit collection + + ///////////////////////////////////////////////////////////////// + return std::make_tuple( + std::move(output_digi_hits), std::move(output_digi_sim_association)); +} + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// finalize ////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// + +void DCHdigi_v01::Create_outputROOTfile_for_debugHistograms() +{ + // save current ROOT directory + TDirectory* currentDir = gDirectory; + + // save the debug histograms in a file + // file is saved and closed when going out of scope + { + auto filename = m_out_debug_filename.value().c_str(); + std::unique_ptr ofile{TFile::Open( filename, "recreate")}; + if (!ofile || ofile->IsZombie()) + { + error() << "Error: Could not open file " << filename << std::endl; + return; + } + ofile->cd(); + hDpw->Write(); + hDww->Write(); + hSxy->Write(); + hSz->Write(); + } + + // Restore previous ROOT directory + if(currentDir && ( not currentDir->IsDestructed() ) ) + currentDir->cd(); + return; +} + +StatusCode DCHdigi_v01::finalize() { + if (m_create_debug_histos.value()) + { + this->Create_outputROOTfile_for_debugHistograms(); + } + + return StatusCode::SUCCESS; +} + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// ThrowException //////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// +void DCHdigi_v01::ThrowException(std::string s) const { + error() << s.c_str() << endmsg; + throw std::runtime_error(s); +} + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// PrintConfiguration //////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// +void DCHdigi_v01::PrintConfiguration(std::ostream& io) { + io << "DCHdigi will use the following components:\n"; + io << "\tGeometry Service: " << m_geoSvcName.value().c_str() << "\n"; + io << "\tUID Service: " << m_uidSvcName.value().c_str() << "\n"; + io << "\tDetector name: " << m_DCH_name.value().c_str() << "\n"; + io << "\t\t|--Volume bitfield: " << m_decoder->fieldDescription().c_str() << "\n"; + io << "\t\t|--Number of layers: " << dch_data->database.size() << "\n"; + io << "\tCluster distributions taken from: " << m_fileDataAlg.value().c_str() << "\n"; + io << "\tResolution along the wire (mm): " << m_z_resolution.value() << "\n"; + io << "\tResolution perp. to the wire (mm): " << m_xy_resolution.value() << "\n"; + io << "\tCreate debug histograms: " << ( m_create_debug_histos.value() ? "true" : "false" ) << "\n"; + if( true == m_create_debug_histos.value() ) + io << "\t\t|--Name of output file with debug histograms: " << m_out_debug_filename.value() << "\n"; + + + return; +} + +void DCHdigi_v01::PrepareRandomEngine(const edm4hep::EventHeaderCollection& headers) const { + uint32_t evt_n = headers[0].getEventNumber(); + uint32_t run_n = headers[0].getRunNumber(); + size_t seed = m_uidSvc->getUniqueID(evt_n, run_n, this->name()); + m_engine.seed(seed); + // advance internal state to minimize possibility of creating correlations + m_engine.discard(10); + for (int i = 0; i < 10; ++i) + myRandom.Rndm(); +} + + + +/////////////////////////////////////////////////////////////////////////////////////// +/////////////////////// CalculateNClusters //////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////////// + +bool DCHdigi_v01::IsParticleCreatedInsideDriftChamber(const edm4hep::MCParticle& thisParticle) const { + auto vertex = thisParticle.getVertex(); // in mm + auto vertexRsquared = vertex[0] * vertex[0] + vertex[1] * vertex[1]; + auto vertexZabs = std::fabs(vertex[2]); + float DCH_halflengh = dch_data->Lhalf / dd4hep::mm; //2000; + float DCH_rin_squared = std::pow(dch_data->rin / dd4hep::mm, 2); //350 * 350; + float DCH_rout_squared = std::pow(dch_data->rout / dd4hep::mm, 2); //2000 * 2000; + return (vertexZabs < DCH_halflengh) && (vertexRsquared > DCH_rin_squared) && (vertexRsquared < DCH_rout_squared); +} + +std::pair > DCHdigi_v01::CalculateClusters(const edm4hep::SimTrackerHit& input_sim_hit) const { + + const edm4hep::MCParticle& thisParticle = input_sim_hit.getParticle(); + // if gamma, optical photon, or other particle with null mass, or hit with zero energy deposited, return zero clusters + if( 22 == abs(thisParticle.getPDG()) || 0 == thisParticle.getMass() || 0 == input_sim_hit.getEDep() ) + return {0., std::vector{} }; + + /// vector to accumulate the size (number of electrons) of each cluster + std::vector ClSz_vector; + //_________________SET NECESSARY PARAMETERS FOR THE CLS ALGORITHM-----WALAA_________________// + + // Parameters needed for the CL Algo //Added by Walaa/////// + ///from Createclusters.cpp//// + float Eloss = 0.0; + float EIzp = 15.8; + float ExECl1 = 0; + float cut = 1000; //controlla + float EIzs = 25.6; + float ExECl1totRec = 0; + float rndCorr(0); + const int nhEp = 10; + float hEpcut[10] = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000}; + int minE = 1000; + int maxE = 10000; + int binE = 1000; + int nhE = (maxE - minE) / binE; + float maxEx0len(0), ExSgmlen(0); + std::vector CorrpMean; + std::vector CorrpSgm; + std::vector Corrdgmean; + std::vector Corrdgsgm; + std::vector Corrdglfrac; + std::vector Corrdglmpvl; + std::vector Corrdglsgml; + std::vector Corrdglmeang; + std::vector Corrdglsgmg; + float maxEx0(0), maxExSlp(0), ExSgmlep(0), ExSgmhad(0); + float MPVEx(0), SgmEx(0), MeanEx1(0), SgmEx1(0), frac(0), Slp(0), CorrSlp(0), CorrInt(0); + + /*________________________________________________________________________________*/ + bool IsSecondaryWithinDCH = IsParticleCreatedInsideDriftChamber(thisParticle); + + // Momentum from EDM4hep, in GeV + double Momentum = sqrt((input_sim_hit.getMomentum().x * input_sim_hit.getMomentum().x) + + (input_sim_hit.getMomentum().y * input_sim_hit.getMomentum().y) + + (input_sim_hit.getMomentum().z * input_sim_hit.getMomentum().z)); + + int thisparticle_pdgid = thisParticle.getPDG(); + int electron_pdgid = 11; + constexpr double me = 0.511; // electron mass in MeV + constexpr double me_GeV = me * 1e-3; // electron mass in GeV + + /// number of clusters, with size > 1 electron + int NClp(0); + + /// number of clusters, with size = 1 electron + int NCl1(0); + + //___________________________________________________________________ + double thisparticle_mass = (thisParticle.getMass() / 1000.); // mass in GeV, required in MeV + double bg = Momentum / thisparticle_mass; + + CorrpMean = flData->get_ClSzCorrpmean(bg); + CorrpSgm = flData->get_ClSzCorrpsgm(bg); + Corrdgmean = flData->get_ClSzCorrdgmean(bg); + Corrdgsgm = flData->get_ClSzCorrdgsgm(bg); + Corrdglfrac = flData->get_ClSzCorrdglfrac(bg); + Corrdglmpvl = flData->get_ClSzCorrdglmpvl(bg); + Corrdglsgml = flData->get_ClSzCorrdglsgml(bg); + Corrdglmeang = flData->get_ClSzCorrdglmeang(bg); + Corrdglsgmg = flData->get_ClSzCorrdglsgmg(bg); + maxEx0 = flData->get_maxEx0(bg); //379.4 electron 100 gev + maxExSlp = flData->get_maxExSlp(); + + MPVEx = flData->get_MPVExtra(bg); //103.2 + SgmEx = flData->get_SgmExtra(bg); //28.5;// + MeanEx1 = flData->get_MeanExtra1(bg); //13.8;// + SgmEx1 = flData->get_SgmExtra1(bg); //9.72;// + frac = flData->get_FracExtra1(bg); //0.13;// + Slp = flData->get_SlopeExtra1(bg); //7.95;// + CorrSlp = flData->get_ClSzCorrSlp(bg); + CorrInt = flData->get_ClSzCorrInt(bg); + + double Tmax = (2.0 * me * pow(bg, 2) / + (1 + (2.0 * (1 + pow(bg, 2)) * me / thisparticle_mass) + pow(me / thisparticle_mass, 2))) * + 1e+6; + float maxEcut = cut; + if (Tmax < maxEcut) { + maxEcut = Tmax; + } + //___________________________________________________________________ + if (not IsSecondaryWithinDCH) { + // lepton PDG id goes from 11 to 16 (antiparticles have negative id) + bool IsLepton = (11 <= abs(thisparticle_pdgid)) && (16 >= abs(thisparticle_pdgid)); + if (IsLepton) { + ExSgmlep = flData->get_ExSgmlep(); + } else { + // TODO Alvaro: gamma rays treated as hadrons, is that ok? + ExSgmhad = flData->get_ExSgmhad(); + } + + + /*________________________________________________________________________________*/ + + TF1* land = new TF1("land", "landaun"); + land->SetParameter(0, 1); + land->SetParameter(1, MPVEx); + land->SetParameter(2, SgmEx); + land->SetRange(0, maxEcut); + + TF1* exGauss = new TF1("exGauss", "[0]*([1]*TMath::Gaus(x,[2],[3],true)+(1.0-[1])*TMath::Exp(-x/[4])/[4])"); + exGauss->SetParameter(0, 1); + exGauss->SetParameter(1, frac); + exGauss->SetParameter(2, MeanEx1); + exGauss->SetParameter(3, SgmEx1); + exGauss->SetParameter(4, Slp); + exGauss->SetRange(0, 90); + + float totExECl = 0.0; + float ExECl = 0.0; + double LengthTrack = input_sim_hit.getPathLength(); + double Etot_per_track = input_sim_hit.getEDep(); + LengthTrack *= 0.1; //from mm to cm + Eloss = Etot_per_track * 1.e09; //from GeV to eV + maxEx0len = maxEx0 * LengthTrack; //maxEx0 is a parameter + if (IsLepton) { + ExSgmlen = ExSgmlep * TMath::Sqrt(LengthTrack); + } else { + ExSgmlen = ExSgmhad * TMath::Sqrt(LengthTrack); + } + float maxExECl = (Eloss - maxEx0len + this->myRandom.Gaus(0, ExSgmlen)) / maxExSlp; + + if (maxExECl < EIzs) { + maxExECl = 0.0; + } //EIzs const = 25.6 + + debug() << "Eloss= " << Eloss << "EIzs= " << EIzs << "EIzp= " << EIzp << "maxExECl= " << maxExECl << "totExECl" + << totExECl << endmsg; + + // The following loop calculate number of clusters of size > 1, NClp + for (int while1counter = 0; Eloss > (EIzp + EIzs) && maxExECl > totExECl && while1counter < 1e6; while1counter++) { + ExECl = land->GetRandom(0, maxEcut); + + if (ExECl > EIzs) { + Eloss -= EIzp; + if (ExECl > (maxExECl - totExECl)) { + ExECl = maxExECl - totExECl; + } + if (ExECl > Eloss) + ExECl = Eloss; + totExECl += ExECl; + Eloss -= ExECl; + NClp++; + //CLSZ + float tmpCorr = 0.0; + for (int i = 0; i < nhEp; ++i) { + if (ExECl >= (i == 0 ? 0 : hEpcut[i - 1]) && ExECl < hEpcut[i]) { + tmpCorr = this->myRandom.Gaus(CorrpMean[i], CorrpSgm[i]); + } + } + int ClSz = TMath::Nint(ExECl * CorrSlp + CorrInt - tmpCorr); + if (ClSz < 2) { + ClSz = 2; + } + ClSz_vector.push_back(ClSz); + } + } + + debug() << "Eloss= " << Eloss << "EIzp= " << EIzp << endmsg; + + // The following loop calculate number of clusters of size 1, NCl1 + for (int while2counter = 0; Eloss >= EIzp && while2counter < 1e6; while2counter++) { + Eloss -= EIzp; + ExECl1 = exGauss->GetRandom(); + if (ExECl1 > Eloss) { + ExECl1 = Eloss; + } + ExECl1totRec += ExECl1; + NCl1++; + Eloss -= ExECl1; + // cluster size is 1 + ClSz_vector.push_back(1); + } + + } //-- end if particle is not secondary + + // if particle is a delta electron created inside the drift chamber + int NCld(0); + if (IsSecondaryWithinDCH && electron_pdgid == thisparticle_pdgid) { + // 1 delta ray cause 1 cluster NCld (d=delta) + NCld = 1; + // Ekdelta in keV + float Ekdelta = (TMath::Sqrt(Momentum * Momentum + me_GeV * me_GeV) - me_GeV) * 1e6; + + // TODO Alvaro: what is this for? + { + float tmpCl; + int tmphE = (Ekdelta - minE) / binE; + if (tmphE >= nhE) + tmphE = nhE - 1; + if (tmphE == nhE - 1) { + rndCorr = this->myRandom.Uniform(0, 1); + if (rndCorr < Corrdglfrac[tmphE]) { + tmpCl = this->myRandom.Gaus(Corrdglmeang[tmphE], Corrdglsgmg[tmphE]); + } else { + tmpCl = this->myRandom.Landau(Corrdglmpvl[tmphE], Corrdglsgml[tmphE]); + } + } else { + tmpCl = this->myRandom.Gaus(Corrdgmean[tmphE], Corrdgsgm[tmphE]); + } + + int ClSz = TMath::Nint(Ekdelta * CorrSlp + CorrInt - tmpCl); + // TODO Alvaro: should it be 1 instead? + if (ClSz < 2) { + ClSz = 2; + } + ClSz_vector.push_back(ClSz); + } + } + // conclusion: if hit caused by delta electron, NCld = 1, number of electrons ClSz=2 + + // value to be returned, sum of number of clusters size=1, size>1, and clusters caused by delta rays + int total_number_of_clusters = NCl1 + NClp + NCld; + debug() << "Ncl= " << total_number_of_clusters << " NCl1= " << NCl1 << "NClp= " << NClp << "NCld= " << NCld << endmsg; + + if( ClSz_vector.size() != std::size_t(total_number_of_clusters) ) + debug() << "Array of cluster sizes does not match total number of clusters\n"; + + // return {total_number_of_clusters, total_number_of_electrons_over_all_clusters}; + return {total_number_of_clusters, ClSz_vector}; +} diff --git a/DCHdigi/test/test_DCHdigi/check_DCHdigi_output.py b/DCHdigi/test/test_DCHdigi/check_DCHdigi_output.py new file mode 100644 index 0000000..b84d219 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/check_DCHdigi_output.py @@ -0,0 +1,41 @@ +# file: check_DCHdigi_output.py +# author: Alvaro Tolosa-Delgado, CERN 2024 +# to run: python3 check_DCHdigi_output.py +# goal: check distance hit-wire and print out a number: +# 0 : calculation of distance hit-wire was done properly +# 1 : problem with calculation of distance hit-wire +# 2 : problem with calculation of hit-projection position onto the wire + +import ROOT +import sys + +def main(): + exit_code = 0 + + # Open debug output file generated by DCHdigi alg + f = ROOT.TFile("dch_digi_alg_debug.root") + + # Retrieve the hit-wire distance distribution + hDpw = f.Get("hDpw") + hDpw.Rebin(10) + + # Retrieve the distance at which the distance distribution has its maximum + distance_hit_wire_more_frequent = hDpw.GetXaxis().GetBinCenter(hDpw.GetMaximumBin()) + + # Check if the distance hit wire has the maximum around d=0.66cm + if 0.05 < abs(distance_hit_wire_more_frequent - 0.65): + exit_code += 1 + + # Retrieve the hit-projection onto the wire to the wire distance distribution + hDww = f.Get("hDww") + + # Check the integral, excluding the first bin + hDww_integral = hDww.Integral(2, -1) + if hDww_integral != 0: + exit_code += 2 + + return exit_code + +if __name__ == "__main__": + code = main() + sys.exit(code) diff --git a/DCHdigi/test/test_DCHdigi/compact/DCH_standalone_o1_v02.xml b/DCHdigi/test/test_DCHdigi/compact/DCH_standalone_o1_v02.xml new file mode 100644 index 0000000..0ec10a3 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/compact/DCH_standalone_o1_v02.xml @@ -0,0 +1,197 @@ + + + + The compact format of the DCH subdetector, built from XLSX spreadsheet + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,superlayer:5,layer:4,nphi:11,stereosign:-1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/DCHdigi/test/test_DCHdigi/compact/elements_o1_v01.xml b/DCHdigi/test/test_DCHdigi/compact/elements_o1_v01.xml new file mode 100644 index 0000000..f35eb34 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/compact/elements_o1_v01.xml @@ -0,0 +1,884 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/DCHdigi/test/test_DCHdigi/compact/materials_o1_v02.xml b/DCHdigi/test/test_DCHdigi/compact/materials_o1_v02.xml new file mode 100644 index 0000000..189bc20 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/compact/materials_o1_v02.xml @@ -0,0 +1,601 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + CMS ECAL https://gitlab.cern.ch/geant4/geant4/-/blob/master/examples/advanced/composite_calorimeter/dataglobal/material.cms#L183 + + + + + + + + + + + + + + + + 0.98810714*2.7 + 0.011892860*10.49 + + field wire center: 50 um Al (core), 0.3 um Ag (coating) + + + + + 0.98516708*2.7 + 0.014832922*10.49 + + field wires top/bottom: 40 um Al (core), 0.3 um Ag (coating) + + + + + 0.97066175*19.28+0.029338250*19.3 + + sense wire: 20 um W (core), 0.3 um Au (coating) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/DCHdigi/test/test_DCHdigi/runDCHdigi.py b/DCHdigi/test/test_DCHdigi/runDCHdigi.py new file mode 100644 index 0000000..60ea103 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/runDCHdigi.py @@ -0,0 +1,38 @@ +# +# gaudi steering file that runs DCHdigi +# +# to execute: +# k4run runDCHdigi.py + +from Gaudi.Configuration import INFO,DEBUG +from Configurables import EventDataSvc, UniqueIDGenSvc +from k4FWCore import ApplicationMgr, IOSvc + +svc = IOSvc("IOSvc") +svc.input = [ "dch_proton_10GeV.root"] +svc.output = "dch_proton_10GeV_digi.root" + +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +geoservice.detectors = ['./compact/DCH_standalone_o1_v02.xml'] + +from Configurables import DCHdigi_v01 +DCHdigi = DCHdigi_v01("DCHdigi") +DCHdigi.DCH_simhits=["DCHCollection"] +DCHdigi.DCH_name="DCH_v2" +DCHdigi.fileDataAlg="DataAlgFORGEANT.root" +DCHdigi.calculate_dndx=True +DCHdigi.create_debug_histograms=True +DCHdigi.zResolution_mm=1 +DCHdigi.xyResolution_mm=0.1 + + +DCHdigi.OutputLevel=INFO + +mgr = ApplicationMgr( + TopAlg=[DCHdigi], + EvtSel="NONE", + EvtMax=-1, + ExtSvc=[geoservice,EventDataSvc("EventDataSvc"),UniqueIDGenSvc("uidSvc")], + OutputLevel=INFO, +) diff --git a/DCHdigi/test/test_DCHdigi/sim_steering.py b/DCHdigi/test/test_DCHdigi/sim_steering.py new file mode 100644 index 0000000..1796a21 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/sim_steering.py @@ -0,0 +1,597 @@ +from DDSim.DD4hepSimulation import DD4hepSimulation +from g4units import mm, GeV, MeV +SIM = DD4hepSimulation() + +## The compact XML file, or multiple compact files, if the last one is the closer. +SIM.compactFile = ['compact/DCH_standalone_o1_v02.xml'] +## Lorentz boost for the crossing angle, in radian! +SIM.crossingAngleBoost = 0.0 +SIM.enableDetailedShowerMode = False +SIM.enableG4GPS = False +SIM.enableG4Gun = False +SIM.enableGun = True +## InputFiles for simulation .stdhep, .slcio, .HEPEvt, .hepevt, .pairs, .hepmc, .hepmc.gz, .hepmc.xz, .hepmc.bz2, .hepmc3, .hepmc3.gz, .hepmc3.xz, .hepmc3.bz2, .hepmc3.tree.root files are supported +SIM.inputFiles = [] +## Macro file to execute for runType 'run' or 'vis' +SIM.macroFile = "scripts/vis.mac" +## number of events to simulate, used in batch mode +SIM.numberOfEvents = 10 +## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported +SIM.outputFile = "dch.root" +## Physics list to use in simulation +SIM.physicsList = None +## Verbosity use integers from 1(most) to 7(least) verbose +## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS +SIM.printLevel = 3 +## The type of action to do in this invocation +## batch: just simulate some events, needs numberOfEvents, and input file or gun +## vis: enable visualisation, run the macroFile if it is set +## qt: enable visualisation in Qt shell, run the macroFile if it is set +## run: run the macroFile and exit +## shell: enable interactive session +SIM.runType = "qt" +## Skip first N events when reading a file +SIM.skipNEvents = 0 +## Steering file to change default behaviour +SIM.steeringFile = None +## FourVector of translation for the Smearing of the Vertex position: x y z t +SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0] +## FourVector of the Sigma for the Smearing of the Vertex position: x y z t +SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0] + + +################################################################################ +## Helper holding sensitive detector and other actions. +## +## The default tracker and calorimeter sensitive actions can be set with +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.tracker=('Geant4TrackerWeightedAction', {'HitPositionCombination': 2, 'CollectSingleDeposits': False}) +## >>> SIM.action.calo = "Geant4CalorimeterAction" +## +## The default sensitive actions for calorimeters and trackers are applied based on the sensitive type. +## The list of sensitive types can be changed with +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.trackerSDTypes = ['tracker', 'myTrackerSensType'] +## >>> SIM.calor.calorimeterSDTypes = ['calorimeter', 'myCaloSensType'] +## +## For specific subdetectors specific sensitive detectors can be set based on patterns in the name of the subdetector. +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['tpc'] = "TPCSDAction" +## +## and additional parameters for the sensitive detectors can be set when the map is given a tuple +## +## >>> SIM = DD4hepSimulation() +## >>> SIM.action.mapActions['ecal'] =( "CaloPreShowerSDAction", {"FirstLayerNumber": 1} ) +## +## Additional actions can be set as well with the following syntax variations: +## +## >>> SIM = DD4hepSimulation() +## # single action by name only: +## >>> SIM.action.run = "Geant4TestRunAction" +## # multiple actions with comma-separated names: +## >>> SIM.action.event = "Geant4TestEventAction/Action0,Geant4TestEventAction/Action1" +## # single action by tuple of name and parameter dict: +## >>> SIM.action.track = ( "Geant4TestTrackAction", {"Property_int": 10} ) +## # single action by dict of name and parameter dict: +## >>> SIM.action.step = { "name": "Geant4TestStepAction", "parameter": {"Property_int": 10} } +## # multiple actions by list of dict of name and parameter dict: +## >>> SIM.action.stack = [ { "name": "Geant4TestStackAction", "parameter": {"Property_int": 10} } ] +## +## On the command line or in python, these actions can be specified as JSON strings: +## $ ddsim --action.stack '{ "name": "Geant4TestStackAction", "parameter": { "Property_int": 10 } }' +## or +## >>> SIM.action.stack = ''' +## { +## "name": "Geant4TestStackAction", +## "parameter": { +## "Property_int": 10, +## "Property_double": "1.0*mm" +## } +## } +## ''' +## +## +################################################################################ + +## set the default calorimeter action +SIM.action.calo = "Geant4ScintillatorCalorimeterAction" + +## List of patterns matching sensitive detectors of type Calorimeter. +SIM.action.calorimeterSDTypes = ['calorimeter'] + +## set the default event action +SIM.action.event = [] + +## Create a map of patterns and actions to be applied to sensitive detectors. +## +## Example: if the name of the detector matches 'tpc' the TPCSDAction is used. +## +## SIM.action.mapActions['tpc'] = "TPCSDAction" +## +SIM.action.mapActions = {} +SIM.action.mapActions['DCH_v2'] = "Geant4SimpleTrackerAction" + +## set the default run action +SIM.action.run = [] + +## set the default stack action +SIM.action.stack = [] + +## set the default step action +SIM.action.step = [] + +## set the default track action +SIM.action.track = [] + +## set the default tracker action +SIM.action.tracker = ('Geant4TrackerWeightedAction', {'HitPositionCombination': 2, 'CollectSingleDeposits': False}) + +## List of patterns matching sensitive detectors of type Tracker. +SIM.action.trackerSDTypes = ['tracker'] + + +################################################################################ +## Configuration for the magnetic field (stepper) +################################################################################ +SIM.field.delta_chord = 0.25 +SIM.field.delta_intersection = 0.001 +SIM.field.delta_one_step = 0.01 +SIM.field.eps_max = 0.001 +SIM.field.eps_min = 5e-05 +SIM.field.equation = "Mag_UsualEqRhs" +SIM.field.largest_step = 10000.0 +SIM.field.min_chord_step = 0.01 +SIM.field.stepper = "ClassicalRK4" + + +################################################################################ +## Configuration for sensitive detector filters +## +## Set the default filter for 'tracker' +## >>> SIM.filter.tracker = "edep1kev" +## Use no filter for 'calorimeter' by default +## >>> SIM.filter.calo = "" +## +## Assign a filter to a sensitive detector via pattern matching +## >>> SIM.filter.mapDetFilter['FTD'] = "edep1kev" +## +## Or more than one filter: +## >>> SIM.filter.mapDetFilter['FTD'] = ["edep1kev", "geantino"] +## +## Don't use the default filter or anything else: +## >>> SIM.filter.mapDetFilter['TPC'] = None ## or "" or [] +## +## Create a custom filter. The dictionary is used to instantiate the filter later on +## >>> SIM.filter.filters['edep3kev'] = dict(name="EnergyDepositMinimumCut/3keV", parameter={"Cut": 3.0*keV} ) +## +## +################################################################################ + +## +## default filter for calorimeter sensitive detectors; +## this is applied if no other filter is used for a calorimeter +## +SIM.filter.calo = "edep0" + +## list of filter objects: map between name and parameter dictionary +SIM.filter.filters = {} + +## a map between patterns and filter objects, using patterns to attach filters to sensitive detector +SIM.filter.mapDetFilter = {} + +## default filter for tracking sensitive detectors; this is applied if no other filter is used for a tracker +SIM.filter.tracker = None + + +################################################################################ +## Configuration for the Detector Construction. +################################################################################ +SIM.geometry.dumpGDML = "" +SIM.geometry.dumpHierarchy = 0 + +## Print Debug information about Elements +SIM.geometry.enableDebugElements = False + +## Print Debug information about Materials +SIM.geometry.enableDebugMaterials = False + +## Print Debug information about Placements +SIM.geometry.enableDebugPlacements = False + +## Print Debug information about Reflections +SIM.geometry.enableDebugReflections = False + +## Print Debug information about Regions +SIM.geometry.enableDebugRegions = False + +## Print Debug information about Shapes +SIM.geometry.enableDebugShapes = False + +## Print Debug information about Surfaces +SIM.geometry.enableDebugSurfaces = False + +## Print Debug information about Volumes +SIM.geometry.enableDebugVolumes = False + +## Print information about placements +SIM.geometry.enablePrintPlacements = False + +## Print information about Sensitives +SIM.geometry.enablePrintSensitives = False + + +################################################################################ +## Configuration for the GuineaPig InputFiles +################################################################################ + +## Set the number of pair particles to simulate per event. +## Only used if inputFile ends with ".pairs" +## If "-1" all particles will be simulated in a single event +## +SIM.guineapig.particlesPerEvent = "-1" + + +################################################################################ +## Configuration for the DDG4 ParticleGun +################################################################################ + +## direction of the particle gun, 3 vector +SIM.gun.direction = (1,1,1) + +## choose the distribution of the random direction for theta +## +## Options for random distributions: +## +## 'uniform' is the default distribution, flat in theta +## 'cos(theta)' is flat in cos(theta) +## 'eta', or 'pseudorapidity' is flat in pseudorapity +## 'ffbar' is distributed according to 1+cos^2(theta) +## +## Setting a distribution will set isotrop = True +## +#SIM.gun.distribution = 'cos(theta)' + +## Total energy (including mass) for the particle gun. +## +## If not None, it will overwrite the setting of momentumMin and momentumMax +SIM.gun.energy = '10*GeV' + +## Maximal pseudorapidity for random distibution (overrides thetaMin) +SIM.gun.etaMax = None + +## Minimal pseudorapidity for random distibution (overrides thetaMax) +SIM.gun.etaMin = None + +## isotropic distribution for the particle gun +## +## use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions +## if one of these options is not None the random distribution will be set to True and cannot be turned off! +## +SIM.gun.isotrop = False + +## Maximal momentum when using distribution (default = 0.0) +SIM.gun.momentumMax = 10000.0 + +## Minimal momentum when using distribution (default = 0.0) +SIM.gun.momentumMin = 0.0 +SIM.gun.multiplicity = 1 +#SIM.gun.particle = "mu+" +#SIM.gun.particle = "pi+" +#SIM.gun.particle = "kaon+" +SIM.gun.particle = "proton" + +## Maximal azimuthal angle for random distribution +SIM.gun.phiMax = None + +## Minimal azimuthal angle for random distribution +SIM.gun.phiMin = None + +## position of the particle gun, 3 vector +SIM.gun.position = (0.0, 0.0, 0.0) + +## Maximal polar angle for random distribution +SIM.gun.thetaMax = None + +## Minimal polar angle for random distribution +SIM.gun.thetaMin = None + + +################################################################################ +## Configuration for the hepmc3 InputFiles +################################################################################ + +## Set the name of the attribute contraining color flow information index 0. +SIM.hepmc3.Flow1 = "flow1" + +## Set the name of the attribute contraining color flow information index 1. +SIM.hepmc3.Flow2 = "flow2" + +## Set to false if the input should be opened with the hepmc2 ascii reader. +## +## If ``True`` a '.hepmc' file will be opened with the HEPMC3 Reader Factory. +## +## Defaults to true if DD4hep was build with HEPMC3 support. +## +SIM.hepmc3.useHepMC3 = True + + +################################################################################ +## Configuration for Input Files. +################################################################################ + +## Set one or more functions to configure input steps. +## +## The functions must take a ``DD4hepSimulation`` object as their only argument and return the created generatorAction +## ``gen`` (for example). +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: GeneratorAction +## ''' +## from DDG4 import GeneratorAction, Kernel +## # Geant4InputAction is the type of plugin, Cry1 just an identifier +## gen = GeneratorAction(Kernel(), 'Geant4InputAction/Cry1' , True) +## # CRYEventReader is the actual plugin, steeringFile its constructor parameter +## gen.Input = 'CRYEventReader|' + 'steeringFile' +## # we can give a dictionary of Parameters that has to be interpreted by the setParameters function of the plugin +## gen.Parameters = {'DataFilePath': '/path/to/files/data'} +## gen.enableUI() +## return gen +## +## SIM.inputConfig.userInputPlugin = exampleUserPlugin +## +## Repeat function definition and assignment to add multiple input steps +## +## +SIM.inputConfig.userInputPlugin = [] + + +################################################################################ +## Configuration for the generator-level InputFiles +################################################################################ + +## Set the name of the collection containing the MCParticle input. +## Default is "MCParticle". +## +SIM.lcio.mcParticleCollectionName = "MCParticle" + + +################################################################################ +## Configuration for the LCIO output file settings +################################################################################ + +## The event number offset to write in slcio output file. E.g setting it to 42 will start counting events from 42 instead of 0 +SIM.meta.eventNumberOffset = 0 + +## Event parameters to write in every event. Use C/F/I ids to specify parameter type. E.g parameterName/F=0.42 to set a float parameter +SIM.meta.eventParameters = [] + +## The run number offset to write in slcio output file. E.g setting it to 42 will start counting runs from 42 instead of 0 +SIM.meta.runNumberOffset = 0 + + +################################################################################ +## Configuration for the output levels of DDG4 components +################################################################################ + +## Output level for geometry. +SIM.output.geometry = 2 + +## Output level for input sources +SIM.output.inputStage = 3 + +## Output level for Geant4 kernel +SIM.output.kernel = 3 + +## Output level for ParticleHandler +SIM.output.part = 3 + +## Output level for Random Number Generator setup +SIM.output.random = 6 + + +################################################################################ +## Configuration for Output Files. +################################################################################ + +## Use the DD4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceDD4HEP = False + +## Use the EDM4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceEDM4HEP = False + +## Use the LCIO output plugin regardless of outputfilename. +SIM.outputConfig.forceLCIO = False + +## Set a function to configure the outputFile. +## +## The function must take a ``DD4hepSimulation`` object as its only argument and return ``None``. +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: None +## ''' +## from DDG4 import EventAction, Kernel +## dd = dd4hepSimulation # just shorter variable name +## evt_root = EventAction(Kernel(), 'Geant4Output2ROOT/' + dd.outputFile, True) +## evt_root.HandleMCTruth = True or False +## evt_root.Control = True +## output = dd.outputFile +## if not dd.outputFile.endswith(dd.outputConfig.myExtension): +## output = dd.outputFile + dd.outputConfig.myExtension +## evt_root.Output = output +## evt_root.enableUI() +## Kernel().eventAction().add(evt_root) +## return None +## +## SIM.outputConfig.userOutputPlugin = exampleUserPlugin +## # arbitrary options can be created and set via the steering file or command line +## SIM.outputConfig.myExtension = '.csv' +## +SIM.outputConfig.userOutputPlugin = None + + +################################################################################ +## Configuration for the Particle Handler/ MCTruth treatment +################################################################################ + +## Enable lots of printout on simulated hits and MC-truth information +SIM.part.enableDetailedHitsAndParticleInfo = False + +## Keep all created particles +SIM.part.keepAllParticles = False + +## Minimal distance between particle vertex and endpoint of parent after +## which the vertexIsNotEndpointOfParent flag is set +## +SIM.part.minDistToParentVertex = 2.2e-14 + +## MinimalKineticEnergy to store particles created in the tracking region +SIM.part.minimalKineticEnergy = 1.0 + +## Printout at End of Tracking +SIM.part.printEndTracking = False + +## Printout at Start of Tracking +SIM.part.printStartTracking = False + +## List of processes to save, on command line give as whitespace separated string in quotation marks +SIM.part.saveProcesses = ['Decay'] + +## Optionally enable an extended Particle Handler +SIM.part.userParticleHandler = "Geant4TCUserParticleHandler" + + +################################################################################ +## Configuration for the PhysicsList and Monte Carlo particle selection. +## +## To load arbitrary plugins, add a function to be executed. +## +## The function must take the DDG4.Kernel() object as the only argument. +## +## For example, add a function definition and the call to a steering file:: +## +## def setupCerenkov(kernel): +## from DDG4 import PhysicsList +## seq = kernel.physicsList() +## cerenkov = PhysicsList(kernel, 'Geant4CerenkovPhysics/CerenkovPhys') +## cerenkov.MaxNumPhotonsPerStep = 10 +## cerenkov.MaxBetaChangePerStep = 10.0 +## cerenkov.TrackSecondariesFirst = True +## cerenkov.VerboseLevel = 2 +## cerenkov.enableUI() +## seq.adopt(cerenkov) +## ph = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys') +## ph.addParticleConstructor('G4OpticalPhoton') +## ph.VerboseLevel = 2 +## ph.enableUI() +## seq.adopt(ph) +## return None +## +## SIM.physics.setupUserPhysics(setupCerenkov) +## +## # End of example +## +################################################################################ + +## Set of Generator Statuses that are used to mark unstable particles that should decay inside of Geant4. +## +SIM.physics.alternativeDecayStatuses = set() + +## If true, add decay processes for all particles. +## +## Only enable when creating a physics list not based on an existing Geant4 list! +## +SIM.physics.decays = False + +## The name of the Geant4 Physics list. +SIM.physics.list = "FTFP_BERT" + +## location of particle.tbl file containing extra particles and their lifetime information +## +## For example in $DD4HEP/examples/DDG4/examples/particle.tbl +## +SIM.physics.pdgfile = None + +## The global geant4 rangecut for secondary production +## +## Default is 0.7 mm as is the case in geant4 10 +## +## To disable this plugin and be absolutely sure to use the Geant4 default range cut use "None" +## +## Set printlevel to DEBUG to see a printout of all range cuts, +## but this only works if range cut is not "None" +## +SIM.physics.rangecut = 0.7 + +## Set of PDG IDs that will not be passed from the input record to Geant4. +## +## Quarks, gluons and W's Z's etc should not be treated by Geant4 +## +SIM.physics.rejectPDGs = {1, 2, 3, 4, 5, 6, 3201, 3203, 4101, 4103, 21, 23, 24, 5401, 25, 2203, 5403, 3101, 3103, 4403, 2101, 5301, 2103, 5303, 4301, 1103, 4303, 5201, 5203, 3303, 4201, 4203, 5101, 5103, 5503} + +## Set of PDG IDs for particles that should not be passed to Geant4 if their properTime is 0. +## +## The properTime of 0 indicates a documentation to add FSR to a lepton for example. +## +SIM.physics.zeroTimePDGs = {17, 11, 13, 15} + + +################################################################################ +## Properties for the random number generator +################################################################################ + +## If True, calculate random seed for each event basedon eventID and runID +## Allows reproducibility even whenSkippingEvents +SIM.random.enableEventSeed = True +SIM.random.file = None +SIM.random.luxury = 1 +SIM.random.replace_gRandom = True +SIM.random.seed = 123 +SIM.random.type = None + + +################################################################################ +## Configuration for setting commands to run during different phases. +## +## In this section, one can configure commands that should be run during the different phases of the Geant4 execution. +## +## 1. Configuration +## 2. Initialization +## 3. Pre Run +## 4. Post Run +## 5. Terminate / Finalization +## +## For example, one can add +## +## >>> SIM.ui.commandsConfigure = ['/physics_lists/em/SyncRadiation true'] +## +## Further details should be taken from the Geant4 documentation. +## +################################################################################ + +## List of UI commands to run during the 'Configure' phase. +SIM.ui.commandsConfigure = [] + +## List of UI commands to run during the 'Initialize' phase. +SIM.ui.commandsInitialize = [] + +## List of UI commands to run during the 'PostRun' phase. +SIM.ui.commandsPostRun = [] + +## List of UI commands to run during the 'PreRun' phase. +SIM.ui.commandsPreRun = [] + +## List of UI commands to run during the 'Terminate' phase. +SIM.ui.commandsTerminate = [] diff --git a/DCHdigi/test/test_DCHdigi/test_DCHdigi.sh b/DCHdigi/test/test_DCHdigi/test_DCHdigi.sh new file mode 100644 index 0000000..4980638 --- /dev/null +++ b/DCHdigi/test/test_DCHdigi/test_DCHdigi.sh @@ -0,0 +1,25 @@ +#!/bin/bash +# file: test_DCHdigi.sh +# author: Alvaro Tolosa-Delgado, CERN 2024 +# to run: sh + test_DCHdigi.sh +# goal: run sim-digitizer of the DCH v2, and return code printed by check_DCHdigi_output.py + +# run simulation with the drift chamber alone +ddsim --steeringFile sim_steering.py --outputFile 'dch_proton_10GeV.root' -N 10 --runType batch --random.seed 42 + +# download file for cluster counting technique +ifilename="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/IDEA/DataAlgFORGEANT.root" +# -nc, --no-clobber: skip downloads that would download to existing files +wget --no-clobber $ifilename + +# Check if the input file exists +if [[ ! -f "$(basename $ifilename)" ]]; then + echo "Error: File '$(basename $ifilename)' not found." + exit 1 +fi + +# run digitizer for position smearing and cluster counting calculation +k4run runDCHdigi.py + +# check distribution of distance from hit position to the wire +python3 check_DCHdigi_output.py