diff --git a/example/SteeringFile_IDEA_o1_v03.py b/example/SteeringFile_IDEA_o1_v03.py index 8d9713d29..c2ddf15e1 100644 --- a/example/SteeringFile_IDEA_o1_v03.py +++ b/example/SteeringFile_IDEA_o1_v03.py @@ -20,7 +20,7 @@ ## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported SIM.outputFile = "testIDEA_o1_v03.root" ## Physics list to use in simulation -SIM.physicsList = None +SIM.physicsList = "FTFP_BERT" ## Verbosity use integers from 1(most) to 7(least) verbose ## or strings: VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL, ALWAYS SIM.printLevel = 3 @@ -177,7 +177,8 @@ ## default filter for calorimeter sensitive detectors; ## this is applied if no other filter is used for a calorimeter ## -SIM.filter.calo = "edep0" +# note: do not turn on the calo filter, otherwise all optical photons will be killed! +SIM.filter.calo = "" ## list of filter objects: map between name and parameter dictionary SIM.filter.filters = { @@ -251,7 +252,7 @@ ################################################################################ ## direction of the particle gun, 3 vector -SIM.gun.direction = (1.0, 1.0, 1.0) +SIM.gun.direction = (1.0,0.1,0.1) ## choose the distribution of the random direction for theta ## @@ -285,10 +286,10 @@ SIM.gun.isotrop = False ## Maximal momentum when using distribution (default = 0.0) -SIM.gun.momentumMax = 10000.0 +#SIM.gun.momentumMax = 10000.0 ## Minimal momentum when using distribution (default = 0.0) -SIM.gun.momentumMin = 0.0 +#SIM.gun.momentumMin = 0.0 SIM.gun.multiplicity = 1 SIM.gun.particle = "e-" @@ -557,7 +558,7 @@ def Geant4Output2EDM4hep_DRC_plugin(dd4hepSimulation): ## 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 +SIM.physics.rangecut = None ## Set of PDG IDs that will not be passed from the input record to Geant4. ## @@ -621,7 +622,8 @@ def setupOpticalPhysics(kernel): opt = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys') opt.addParticleConstructor('G4OpticalPhoton') opt.VerboseLevel = 1 - opt.BoundaryInvokeSD = True + # set BoundaryInvokeSD to true when using DRC wafer as the SD + # opt.BoundaryInvokeSD = True opt.enableUI() seq.adopt(opt) @@ -653,7 +655,8 @@ def setupDRCFastSim(kernel): phys.adopt(ph) phys.dump() -SIM.physics.setupUserPhysics(setupDRCFastSim) +# turn-on fastsim if the skipScint option of the SD is set to false +#SIM.physics.setupUserPhysics(setupDRCFastSim) ################################################################################ ## Properties for the random number generator diff --git a/plugins/DRCaloFastSimModel.h b/plugins/DRCaloFastSimModel.h index 4d33d733c..3b9273ce2 100644 --- a/plugins/DRCaloFastSimModel.h +++ b/plugins/DRCaloFastSimModel.h @@ -82,42 +82,43 @@ class DRCFiberModel { DRCFiberModel()=default; ~DRCFiberModel()=default; - G4OpBoundaryProcess* pOpBoundaryProc{nullptr}; - G4OpAbsorption* pOpAbsorption{nullptr}; - G4OpWLS* pOpWLS{nullptr}; - G4bool fProcAssigned{false}; - - FastFiberData mDataPrevious{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))}; - FastFiberData mDataCurrent{FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0))}; - - G4int fSafety{2}; - G4double mNtransport{0.}; - G4double mTransportUnit{0.}; - G4ThreeVector mFiberPos{G4ThreeVector(0)}; - G4ThreeVector mFiberAxis{G4ThreeVector(0)}; - G4bool fKill{false}; - G4bool fTransported{false}; - G4bool fSwitch{true}; - G4int fVerbose{0}; + G4OpBoundaryProcess* pOpBoundaryProc = nullptr; + G4OpAbsorption* pOpAbsorption = nullptr; + G4OpWLS* pOpWLS = nullptr; + G4bool fProcAssigned = false; + + FastFiberData mDataPrevious = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0)); + FastFiberData mDataCurrent = FastFiberData(0, 0., 0., 0., G4ThreeVector(0), G4ThreeVector(0), G4ThreeVector(0)); + + G4int fSafety = 1; + G4double mNtransport = 0.; + G4double mTransportUnit = 0.; + G4ThreeVector mFiberPos = G4ThreeVector(0); + G4ThreeVector mFiberAxis = G4ThreeVector(0); + G4bool fKill = false; + G4bool fTransported = false; + G4bool fSwitch = true; + G4int fVerbose = 0; G4bool checkTotalInternalReflection(const G4Track *track) { if (!fProcAssigned) setPostStepProc(track); // locate OpBoundaryProcess only once + G4int theStatus = pOpBoundaryProc->GetStatus(); + + if (fVerbose > 1) { + std::cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID(); + std::cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus; + std::cout << " | Track status = " << track->GetTrackStatus(); + std::cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << std::endl; + } + if (track->GetTrackStatus() == fStopButAlive || track->GetTrackStatus() == fStopAndKill) return false; // accumulate step length mDataCurrent.AddStepLengthInterval(track->GetStepLength()); - G4int theStatus = pOpBoundaryProc->GetStatus(); - - if (fVerbose > 1) { - G4cout << "DRCFiberModel::checkTotalInternalReflection | TrackID = " << std::setw(4) << track->GetTrackID(); - G4cout << " | G4OpBoundaryProcessStatus = " << std::setw(2) << theStatus; - G4cout << " | StepLength = " << std::setw(9) << track->GetStepLength() << G4endl; - } - // skip exceptional iteration with FresnelReflection if (theStatus == G4OpBoundaryProcessStatus::FresnelReflection) mDataCurrent.SetOpBoundaryStatus(theStatus); diff --git a/plugins/FiberDRCaloSDAction.cpp b/plugins/FiberDRCaloSDAction.cpp index 3341c4382..37dc55d9c 100644 --- a/plugins/FiberDRCaloSDAction.cpp +++ b/plugins/FiberDRCaloSDAction.cpp @@ -47,6 +47,7 @@ namespace dd4hep */ struct DRCData { + public: bool skipScint = true; DRCFiberModel fastfiber; @@ -58,66 +59,134 @@ namespace dd4hep G4float fTimeEnd; G4float fWavlenStep; G4float fTimeStep; - static const int fArraySize = 24; + private: // from 900 nm to 460 nm - double fGraph_X[fArraySize] = { 1.37760 * CLHEP::eV, - 1.45864 * CLHEP::eV, - 1.54980 * CLHEP::eV, - 1.65312 * CLHEP::eV, - 1.71013 * CLHEP::eV, - 1.77120 * CLHEP::eV, - 1.83680 * CLHEP::eV, - 1.90745 * CLHEP::eV, - 1.98375 * CLHEP::eV, - 2.06640 * CLHEP::eV, - - 2.10143 * CLHEP::eV, - 2.13766 * CLHEP::eV, - 2.17516 * CLHEP::eV, - 2.21400 * CLHEP::eV, - 2.25426 * CLHEP::eV, - 2.29600 * CLHEP::eV, - 2.33932 * CLHEP::eV, - 2.38431 * CLHEP::eV, - 2.43106 * CLHEP::eV, - 2.47968 * CLHEP::eV, - - 2.53029 * CLHEP::eV, - 2.58300 * CLHEP::eV, - 2.63796 * CLHEP::eV, - 2.69531 * CLHEP::eV }; + const std::vector fGraph_X = { + 1.37760 * CLHEP::eV, + 1.45864 * CLHEP::eV, + 1.54980 * CLHEP::eV, + 1.65312 * CLHEP::eV, + 1.71013 * CLHEP::eV, + 1.77120 * CLHEP::eV, + 1.83680 * CLHEP::eV, + 1.90745 * CLHEP::eV, + 1.98375 * CLHEP::eV, + 2.06640 * CLHEP::eV, + + 2.10143 * CLHEP::eV, + 2.13766 * CLHEP::eV, + 2.17516 * CLHEP::eV, + 2.21400 * CLHEP::eV, + 2.25426 * CLHEP::eV, + 2.29600 * CLHEP::eV, + 2.33932 * CLHEP::eV, + 2.38431 * CLHEP::eV, + 2.43106 * CLHEP::eV, + 2.47968 * CLHEP::eV, + + 2.53029 * CLHEP::eV, + 2.58300 * CLHEP::eV, + 2.63796 * CLHEP::eV, + 2.69531 * CLHEP::eV, + 2.75520 * CLHEP::eV, + 2.81782 * CLHEP::eV, + 2.88335 * CLHEP::eV, + 2.95200 * CLHEP::eV, + 3.09960 * CLHEP::eV, + 3.54241 * CLHEP::eV, + + 4.13281 * CLHEP::eV + }; // filter efficiency of the Kodak Wratten No.9 filter - double fGraph_Y[fArraySize] = { 0.903, - 0.903, - 0.903, - 0.903, - 0.903, - 0.903, - 0.902, - 0.901, - 0.898, - 0.895, - - 0.893, - 0.891, - 0.888, - 0.883, - 0.87, - 0.838, - 0.76, - 0.62, - 0.488, - 0.345, - - 0.207, - 0.083, - 0.018, - 0. }; - - G4double wavToE(G4double wav) { return CLHEP::h_Planck * CLHEP::c_light / wav; } - - int findWavBin(G4double en) { + const std::vector fKodakEff = { + 0.903, + 0.903, + 0.903, + 0.903, + 0.903, + 0.903, + 0.902, + 0.901, + 0.898, + 0.895, + + 0.893, + 0.891, + 0.888, + 0.883, + 0.87, + 0.838, + 0.76, + 0.62, + 0.488, + 0.345, + + 0.207, + 0.083, + 0.018, + 0., + 0., + 0., + 0., + 0., + 0., + 0., + + 0. + }; + + // SiPM efficiency Hamamatsu S14160-1310PS + // TODO migrate this part to the digitization step! + // Note: Ideally, this should be part of the digitization step. + // (not the simulation step) + // But, to do this, we need to store the distribution + // of the optical photon wavelength. + // While we can develop another feature to enable this, + // let's emulate the SiPM efficiency in the simulation step for now. + // We just need a working code and without this, + // the number of Cherenkov photon will be order of magnitude higher + const std::vector fSipmEff = { + 0.02, + 0.025, + 0.045, + 0.06, + 0.0675, + 0.075, + 0.0925, + 0.11, + 0.125, + 0.14, + + 0.146, + 0.152, + 0.158, + 0.164, + 0.17, + 0.173, + 0.176, + 0.178, + 0.179, + 0.18, + + 0.181, + 0.182, + 0.183, + 0.184, + 0.18, + 0.173, + 0.166, + 0.158, + 0.15, + 0.12, + + 0.05 + }; + + public: + G4double wavToE(G4double wav) const { return CLHEP::h_Planck * CLHEP::c_light / wav; } + + int findWavBin(G4double en) const { int i = 0; for (; i < fWavBin + 1; i++) { if (en < wavToE((fWavlenStart - static_cast(i) * fWavlenStep) * CLHEP::nm)) @@ -127,7 +196,7 @@ namespace dd4hep return fWavBin + 1 - i; } - int findTimeBin(G4double stepTime) { + int findTimeBin(G4double stepTime) const { int i = 0; for (; i < fTimeBin + 1; i++) { if (stepTime < ((fTimeStart + static_cast(i) * fTimeStep) * CLHEP::ns)) @@ -139,17 +208,17 @@ namespace dd4hep // Linear interpolation function for calculating the efficiency of yellow filter // used in rejectedByYellowFilter - double getFilterEff(G4double G4energy) { + double getFilterEff(const std::vector& yarray, const G4double G4energy) const { // If the photon energy <= 1.37760 eV, than return maximum filter efficiency - if (G4energy <= fGraph_X[0]) - return fGraph_Y[0]; + if (G4energy <= fGraph_X.at(0)) + return yarray.at(0); - for(int idx = 1; idx < fArraySize; ++idx) { - if (G4energy <= fGraph_X[idx]) { - double x1 = fGraph_X[idx - 1]; - double x2 = fGraph_X[idx]; - double y1 = fGraph_Y[idx - 1]; - double y2 = fGraph_Y[idx]; + for(unsigned idx = 1; idx < yarray.size(); ++idx) { + if (G4energy <= fGraph_X.at(idx)) { + double x1 = fGraph_X.at(idx-1); + double x2 = fGraph_X.at(idx); + double y1 = yarray.at(idx-1); + double y2 = yarray.at(idx); // return linear interpolated filter efficiency return (y1 + ((y2 - y1) / (x2 - x1))*(G4energy - x1)); @@ -160,8 +229,8 @@ namespace dd4hep } // If true, then the photon is rejected by yellow filter - bool rejectedByYellowFilter(G4double G4energy, double rndVal) { - double FilterEff = getFilterEff(G4energy); // Get efficiency of filter using photon's energy + bool rejectedByYellowFilter(G4double G4energy, double rndVal) const { + const double FilterEff = getFilterEff(fKodakEff,G4energy); // Get efficiency of filter using photon's energy // filter efficiency == probability of photon accepted by filter // == Probability of random value (from uniform distribution with range 0 ~ 1) smaller than filter efficiency @@ -169,6 +238,12 @@ namespace dd4hep return (rndVal > FilterEff); } + // check sipm efficiency + // TODO migrate this to the digitization step + bool rejectedBySiPM(double G4energy, double rndVal) const { + return (rndVal > getFilterEff(fSipmEff,G4energy)); + } + DRCData() : fWavBin(120), fTimeBin(650), fWavlenStart(900.), fWavlenEnd(300.), fTimeStart(5.), fTimeEnd(70.) @@ -191,8 +266,8 @@ namespace dd4hep initialize(); InstanceCount::increment(this); - m_userData.fastfiber.fSafety = 0.; - m_userData.fastfiber.fVerbose = 0.; + m_userData.fastfiber.fSafety = 1; + m_userData.fastfiber.fVerbose = 0; } /// Define collections created by this sensitive action object @@ -249,10 +324,22 @@ namespace dd4hep // now the touchable is the tower // get local position and volumeID - dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager(); + // dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager(); + // dd4hep::VolumeID volID = volMgr.volumeID(touchable); + // note the above method started to throw warnings since the PR DD4hep#1390 + // so we switch to the copy number + + // trick: we retrieve the tower's volID by using the fact that + // in the DRconstructor.cpp, the first 32bits of the volID + // contain the information relevant to the tower + // and this 32bit integer becomes the tower's copy number + auto copyNo = touchable->GetCopyNumber(); + + // ideally the below operation should be GridDRcalo_k4geo::convertFirst32to64 + // but this is trivial, so let's avoid doing dynamic cast + dd4hep::VolumeID volID = static_cast(copyNo); G4ThreeVector global = step->GetPostStepPoint()->GetPosition(); G4ThreeVector local = touchable->GetHistory()->GetTopTransform().TransformPoint(global); - dd4hep::VolumeID volID = volMgr.volumeID(touchable); // convert G4 position to dd4hep position dd4hep::Position loc(local.x() * dd4hep::millimeter / CLHEP::millimeter, local.y() * dd4hep::millimeter / CLHEP::millimeter, local.z() * dd4hep::millimeter / CLHEP::millimeter); @@ -299,6 +386,19 @@ namespace dd4hep double timeUnit = m_userData.fastfiber.mDataCurrent.globalTime - m_userData.fastfiber.mDataPrevious.globalTime; double timeShift = timeUnit * m_userData.fastfiber.mNtransport; + // check wheter the optical photon will be rejected by the SiPM + // TODO migrate this to the digitization step + Geant4Event& evt = context()->event(); + Geant4Random& rnd = evt.random(); // Get random generator + double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1] + G4double energy = step->GetTrack()->GetTotalEnergy(); + + if (m_userData.rejectedBySiPM(energy, rndVal)) { + step->GetTrack()->SetTrackStatus(fStopAndKill); + + return false; + } + // default hit (optical photon count) Geant4DRCalorimeter::Hit* drHit = coll->find(CellIDCompare(cID)); @@ -314,7 +414,6 @@ namespace dd4hep } // everything should be in the G4 unit - G4double energy = step->GetTrack()->GetTotalEnergy(); // (approximate) timing at the end of the fiber G4double hitTime = step->GetPostStepPoint()->GetGlobalTime() + timeShift; @@ -392,15 +491,20 @@ namespace dd4hep dd4hep::DDSegmentation::VolumeID ceren = static_cast(m_segmentation->decoder()->get(cID, "c")); bool IsCeren = static_cast(ceren); - if (!IsCeren) { - Geant4Event& evt = context()->event(); - Geant4Random& rnd = evt.random(); // Get random generator - double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1] + Geant4Event& evt = context()->event(); + Geant4Random& rnd = evt.random(); // Get random generator + double rndVal = rnd.uniform(0, 1); // Get random number from uniform distribution [0, 1] + if (!IsCeren) { if (m_userData.rejectedByYellowFilter(energy, rndVal)) - return true; + return false; } + // check wheter the optical photon will be rejected by the SiPM + // TODO migrate this to the digitization step + if (m_userData.rejectedBySiPM(energy, rndVal)) + return false; + if (!hit) { hit = new Geant4DRCalorimeter::Hit(m_userData.fWavlenStep, m_userData.fTimeStep); hit->cellID = cID;