// // // RPC-SD.cc // // // class RPC_SD - implementation // // // ************************************* // * Implements RPC Sensitive Detector * // ************************************* // // // ---------------------------------------------------------------------- // Author: Miguel Couceiro // Version: 3.0 // Date: April 2011 // ---------------------------------------------------------------------- // ---------------------------------------------------------------------- // Revisions: // // Authors: Miguel Couceiro and Ana Luísa Lopes // Version: 2.0 // Date: October 2018 // ---------------------------------------------------------------------- /* ******************** */ /* **** Directives **** */ /* ******************** */ // #ifdef DEBUG // #undef DEBUG // #endif // DEBUG #ifndef DEBUG #define DEBUG #endif // DEBUG /* ********************** */ /* **** Header Files **** */ /* ********************** */ #include #include #include #include #include #include #include #include #include #include #include #include "../include/RPC-SD.hh" /* ***************************************** */ /* **** Static variables initialization **** */ /* ***************************************** */ RPC_SD* RPC_SD::m_pSD = 0; /* ********************************* */ /* **** Constructors/Destructor **** */ /* ********************************* */ // Default Constructor RPC_SD::RPC_SD () : G4VSensitiveDetector ("") { // **** Initialize class variables // General counters m_nNumberOfAcceptedEvents = 0; m_nNumberOfRejectedEvents = 0; // Auxiliary variables m_nNumAllTracks = 0; m_nAllTracksID.clear(); m_nAllParentsID.clear(); // General source data m_nSourceID = 0; m_vSourcePosition = G4ThreeVector(0,0,0); m_vAnnihilationPosition = G4ThreeVector(0,0,0); // General photon data m_nNumPhotons = 0; m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); // Phantom related data m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); // Scanner related data m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); // **** Register Hits Collection // Primary Photons Hits m_pPrimPhotHitsCollection = 0; collectionName.insert("PrimPhotHitsCollection"); // All Photons Hits m_pAllPhotHitsCollection = 0; collectionName.insert("AllPhotHitsCollection"); // Electron Hits m_pElectHitsCollection = 0; collectionName.insert("ElectHitsCollection"); // RPC Hits m_pRPCHitsCollection = 0; collectionName.insert("RPCHitsCollection"); // **** Register the Sensitive Detector in the sensitive detector manager G4SDManager* pSDManager = G4SDManager::GetSDMpointer(); pSDManager->AddNewDetector(this); } // Copy Constructor RPC_SD::RPC_SD (RPC_SD& refRPC_SD) : G4VSensitiveDetector (dynamic_cast (refRPC_SD)) { // **** Initialize class variables // General counters m_nNumberOfAcceptedEvents = refRPC_SD.GetNumberOfAceptedEvents(); m_nNumberOfRejectedEvents = refRPC_SD.GetNumberOfRejectedEvents(); // Auxiliary variables m_nNumAllTracks = refRPC_SD.GetNumberOfAllTracks(); m_nAllTracksID.clear(); m_nAllParentsID.clear(); for (G4int i = 0 ; i < (G4int) m_nNumAllTracks ; i++) { m_nAllTracksID.push_back(refRPC_SD.GetAllTrackID(i)); m_nAllParentsID.push_back(refRPC_SD.GetAllParentTrackID(i)); } // General source data m_nSourceID = refRPC_SD.GetSourceID(); m_vSourcePosition = refRPC_SD.GetSourcePosition(); m_vAnnihilationPosition = refRPC_SD.GetAnnihilationPosition(); // Photon data m_nNumPhotons = refRPC_SD.GetNumberOfPhotonTracks(); m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); for (G4int i = 0 ; i < m_nNumPhotons ; i++) { m_nPhotParentID.push_back(refRPC_SD.GetPhotonParentID(i)); m_nPhotTrackID.push_back(refRPC_SD.GetPhotonTrackID(i)); m_nPhotNumHits.push_back(refRPC_SD.GetPhotonNumberOfHits(i)); m_nPhantomPhotNumRaylInter.push_back(refRPC_SD.GetPhantomNumberOfRayleighInteractions(i)); m_nPhantomPhotNumCompInter.push_back(refRPC_SD.GetPhantomNumberOfComptonInteractions(i)); m_nPhantomPhotNumTotalInter.push_back(refRPC_SD.GetPhantomNumberOfTotalInteractions(i)); m_fScannerPhotEntryEnergy.push_back(refRPC_SD.GetScannerPhotonEntryEnergy(i)); m_vScannerPhotEntryPos.push_back(refRPC_SD.GetScannerPhotonEntryPosition(i)); m_vScannerPhotEntryMomDir.push_back(refRPC_SD.GetScannerPhotonEntryMomentumDirection(i)); m_nScannerPhotNumRaylInter.push_back(refRPC_SD.GetScannerNumberOfRayleighInteractions(i)); m_nScannerPhotNumCompInter.push_back(refRPC_SD.GetScannerNumberOfComptonInteractions(i)); m_nScannerPhotNumPhotoInter.push_back(refRPC_SD.GetScannerNumberOfPhotoelectricInteractions(i)); m_nScannerPhotNumTotalInter.push_back(refRPC_SD.GetScannerNumberOfTotalInteractions(i)); } // **** Register Hits Collection // Primary Photons Hits m_pPrimPhotHitsCollection = 0; collectionName.insert("PrimPhotHitsCollection"); // All Photons Hits m_pAllPhotHitsCollection = 0; collectionName.insert("AllPhotHitsCollection"); // Electron Hits m_pElectHitsCollection = 0; collectionName.insert("ElectHitsCollection"); // RPC Hits m_pRPCHitsCollection = 0; collectionName.insert("RPCHitsCollection"); // **** Register the Sensitive Detector in the sensitive detector manager G4SDManager* pSDManager = G4SDManager::GetSDMpointer(); pSDManager->AddNewDetector(this); } // Constructor RPC_SD::RPC_SD (G4String name) : G4VSensitiveDetector (name) { // **** Initialize class variables // General counters m_nNumberOfAcceptedEvents = 0; m_nNumberOfRejectedEvents = 0; // Auxiliary variables m_nNumAllTracks = 0; m_nAllTracksID.clear(); m_nAllParentsID.clear(); // General source data m_nSourceID = 0; m_vSourcePosition = G4ThreeVector(0,0,0); m_vAnnihilationPosition = G4ThreeVector(0,0,0); // General photon data m_nNumPhotons = 0; m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); // Phantom related data m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); // Scanner related data m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); // **** Register Hits Collection // Primary Photons Hits m_pPrimPhotHitsCollection = 0; collectionName.insert("PrimPhotHitsCollection"); // All Photons Hits m_pAllPhotHitsCollection = 0; collectionName.insert("AllPhotHitsCollection"); // Electron Hits m_pElectHitsCollection = 0; collectionName.insert("ElectHitsCollection"); // RPC Hits m_pRPCHitsCollection = 0; collectionName.insert("RPCHitsCollection"); // **** Register the Sensitive Detector in the sensitive detector manager G4SDManager* pSDManager = G4SDManager::GetSDMpointer(); pSDManager->AddNewDetector(this); } // Destructor RPC_SD::~RPC_SD () { // **** Reset class data // Auxiliary variables m_nAllTracksID.clear(); m_nAllParentsID.clear(); // General photon data m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); // Phantom related data m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); // Scanner related data m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); } /* ***************************** */ /* **** Protected Operators **** */ /* ***************************** */ // Explicit declare the assignment operator to avoid copies of the singleton class RPC_SD& RPC_SD::operator= (const RPC_SD&) { return *this; } /* ******************************************** */ /* **** Virtual public interface functions **** */ /* ******************************************** */ // Initialize - Invoked at the beginning of each event void RPC_SD::Initialize (G4HCofThisEvent* pHCE) { // **** Get PrimPhotHits Collection ID, create the collection, and add it to the hits collection of this event static G4int PrimPhot_HCID = -1; m_pPrimPhotHitsCollection = new RPC_PrimaryPhotonHitsCollection(SensitiveDetectorName, collectionName[0]); if (PrimPhot_HCID < 0) PrimPhot_HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); pHCE->AddHitsCollection(PrimPhot_HCID, m_pPrimPhotHitsCollection); // **** Get AllPhotHits Collection ID, create the collection, and add it to the hits collection of this event static G4int AllPhot_HCID = -1; m_pAllPhotHitsCollection = new RPC_AllPhotonHitsCollection(SensitiveDetectorName, collectionName[1]); if (AllPhot_HCID < 0) AllPhot_HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[1]); pHCE->AddHitsCollection(AllPhot_HCID, m_pAllPhotHitsCollection); // **** Get ElectHits Collection ID, create the collection, and add it to the hits collection of this event static G4int Elect_HCID = -1; m_pElectHitsCollection = new RPC_ElectronHitsCollection(SensitiveDetectorName, collectionName[2]); if (Elect_HCID < 0) Elect_HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[2]); pHCE->AddHitsCollection(Elect_HCID, m_pElectHitsCollection); // **** Get RPCHits Collection ID, create the collection, and add it to the hits collection of this event static G4int RPC_HCID = -1; m_pRPCHitsCollection = new RPC_HitsCollection(SensitiveDetectorName, collectionName[3]); if (RPC_HCID < 0) RPC_HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[3]); pHCE->AddHitsCollection(RPC_HCID, m_pRPCHitsCollection); // **** Initialize class variables // Auxiliary variables m_nAllTracksID.clear(); m_nAllParentsID.clear(); // General source data const PrimaryGeneratorAction* pPrimaryGeneratorAction = dynamic_cast ((G4RunManager::GetRunManager())->GetUserPrimaryGeneratorAction()); m_nSourceID = pPrimaryGeneratorAction->GetCurrentSourceID(); m_vSourcePosition = pPrimaryGeneratorAction->GetCurrentSourcePosition(); m_vAnnihilationPosition = m_vSourcePosition; // General photon data m_nNumPhotons = 0; m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); // Phantom related data m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); // Scanner related data m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); } // EndOfEvent - Invoked in the end of each event void RPC_SD::EndOfEvent (G4HCofThisEvent*) { // **** Save all hits // Get RunAction const RPC_RunAction* pRunAction = dynamic_cast ((G4RunManager::GetRunManager())->GetUserRunAction()); // Save PrimPhot hits to ASCII file if (pRunAction->GetPrimPhotHitsSaveASCII()) { FILE* pFlPnt = pRunAction->GetPrimPhotHitsASCIIFilePointer(); G4int numHits = m_pPrimPhotHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pPrimPhotHitsCollection)[i]->WriteHitToASCIIFile(pFlPnt); } // Save PrimPhot hits to Binary file if (pRunAction->GetPrimPhotHitsSaveBinary()) { FILE* pFlPnt = pRunAction->GetPrimPhotHitsBinaryFilePointer(); G4int numHits = m_pPrimPhotHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pPrimPhotHitsCollection)[i]->WriteHitToBinaryFile(pFlPnt); } // Save AllPhot hits to ASCII file if (pRunAction->GetAllPhotHitsSaveASCII()) { FILE* pFlPnt = pRunAction->GetAllPhotHitsASCIIFilePointer(); G4int numHits = m_pAllPhotHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pAllPhotHitsCollection)[i]->WriteHitToASCIIFile(pFlPnt); } // Save AllPhot hits to Binary file if (pRunAction->GetAllPhotHitsSaveBinary()) { FILE* pFlPnt = pRunAction->GetAllPhotHitsBinaryFilePointer(); G4int numHits = m_pAllPhotHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pAllPhotHitsCollection)[i]->WriteHitToBinaryFile(pFlPnt); } // Save Elect hits to ASCII file if (pRunAction->GetElectHitsSaveASCII()) { FILE* pFlPnt = pRunAction->GetElectHitsASCIIFilePointer(); G4int numHits = m_pElectHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pElectHitsCollection)[i]->WriteHitToASCIIFile(pFlPnt); } // Save Elect hits to Binary file if (pRunAction->GetElectHitsSaveBinary()) { FILE* pFlPnt = pRunAction->GetElectHitsBinaryFilePointer(); G4int numHits = m_pElectHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pElectHitsCollection)[i]->WriteHitToBinaryFile(pFlPnt); } // Save RPC hits to ASCII file if (pRunAction->GetRPCHitsSaveASCII()) { FILE* pFlPnt = pRunAction->GetRPCHitsASCIIFilePointer(); G4int numHits = m_pRPCHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pRPCHitsCollection)[i]->WriteHitToASCIIFile(pFlPnt); fflush(pFlPnt); } // Save RPC hits to Binary file if (pRunAction->GetRPCHitsSaveBinary()) { FILE* pFlPnt = pRunAction->GetRPCHitsBinaryFilePointer(); G4int numHits = m_pRPCHitsCollection->entries(); for (G4int i = 0 ; i < numHits ; i++) (*m_pRPCHitsCollection)[i]->WriteHitToBinaryFile(pFlPnt); } // **** Clear class variables // Auxiliary variables m_nAllTracksID.clear(); m_nAllParentsID.clear(); // General photon data m_nPhotParentID.clear(); m_nPhotTrackID.clear(); m_nPhotNumHits.clear(); // Phantom related data m_nPhantomPhotNumRaylInter.clear(); m_nPhantomPhotNumCompInter.clear(); m_nPhantomPhotNumTotalInter.clear(); // Scanner related data m_fScannerPhotEntryEnergy.clear(); m_vScannerPhotEntryPos.clear(); m_vScannerPhotEntryMomDir.clear(); m_nScannerPhotNumRaylInter.clear(); m_nScannerPhotNumCompInter.clear(); m_nScannerPhotNumPhotoInter.clear(); m_nScannerPhotNumTotalInter.clear(); } // PrintAll function void RPC_SD::PrintAll (void) { // **** Print hits information G4int nEventVerboseLevel = G4EventManager::GetEventManager()->GetVerboseLevel(); if (nEventVerboseLevel >= 3) { G4cout << "\n\n\n\n"; G4cout << " Event - ID: " << G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID() << " \n"; G4cout << "\n\n"; G4cout << " Run - Number of Accepted Events: " << m_nNumberOfAcceptedEvents << " \n"; G4cout << " Run - Number of Rejected Events: " << m_nNumberOfRejectedEvents << " \n"; G4cout << "\n\n"; G4cout << " Auxiliary - Total Number of Tracks: " << m_nNumAllTracks << " \n"; for (G4int i = 0 ; i < m_nNumAllTracks ; i++) { G4cout << " Auxiliary - Track Number: " << i << " \n"; G4cout << " Auxiliary - Track ID: " << (((G4int) m_nAllTracksID.size() > i) ? m_nAllTracksID[i] : -1) << " \n"; G4cout << " Auxiliary - Track Parent ID: " << (((G4int) m_nAllParentsID.size() > i) ? m_nAllParentsID[i] : -1) << " \n"; } G4cout << "\n\n"; G4cout << " Source - Source ID: " << m_nSourceID << " \n"; G4cout << " Source - Position: " << m_vSourcePosition/mm << " mm \n"; G4cout << " Source - Annihilation Position: " << m_vAnnihilationPosition/mm << " mm \n"; G4cout << "\n\n"; G4cout << " Photon - Number of Photons: " << m_nNumPhotons << " \n"; for (G4int i = 0 ; i < m_nNumPhotons ; i++) { G4cout << "\n"; G4cout << " Photon - Photon Number: " << i << " \n"; G4cout << "\n"; G4cout << " Photon - Parent ID: " << (((G4int) m_nPhotParentID.size() > i) ? m_nPhotParentID[i] : -1) << " \n"; G4cout << " Photon - Track ID: " << (((G4int) m_nPhotTrackID.size() > i) ? m_nPhotTrackID[i] : -1) << " \n"; G4cout << " Photon - Number of Hits: " << (((G4int) m_nPhotNumHits.size() > i) ? m_nPhotNumHits[i] : -1) << " \n"; G4cout << "\n"; G4cout << " Photon - Phantom Number of Rayleigh Interactions: " << (((G4int) m_nPhantomPhotNumRaylInter.size() > i) ? m_nPhantomPhotNumRaylInter[i] : -1) << " \n"; G4cout << " Photon - Phantom Number of Compton Interactions: " << (((G4int) m_nPhantomPhotNumCompInter.size() > i) ? m_nPhantomPhotNumCompInter[i] : -1) << " \n"; G4cout << " Photon - Phantom Total Number of Interactions: " << (((G4int) m_nPhantomPhotNumTotalInter.size() > i) ? m_nPhantomPhotNumTotalInter[i] : -1) << " \n"; G4cout << "\n"; G4cout << " Photon - Scanner Number of Rayleigh Interactions: " << (((G4int) m_nScannerPhotNumRaylInter.size() > i) ? m_nScannerPhotNumRaylInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Number of Compton Interactions: " << (((G4int) m_nScannerPhotNumCompInter.size() > i) ? m_nScannerPhotNumCompInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Number of Photoelectric Interactions: " << (((G4int) m_nScannerPhotNumPhotoInter.size() > i) ? m_nScannerPhotNumPhotoInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Total Number of Interactions: " << (((G4int) m_nScannerPhotNumTotalInter.size() > i) ? m_nScannerPhotNumTotalInter[i] : -1) << " \n"; G4cout << "\n"; G4cout << " Photon - Scanner Entry Energy: " << (((G4int) m_fScannerPhotEntryEnergy.size() > i) ? m_fScannerPhotEntryEnergy[i]/keV : -1) << " keV\n"; G4cout << " Photon - Scanner Entry Position: " << (((G4int) m_vScannerPhotEntryPos.size() > i) ? m_vScannerPhotEntryPos[i]/keV : -1) << " keV\n"; G4cout << " Photon - Scanner Entry Momentum Direction: " << (((G4int) m_vScannerPhotEntryMomDir.size() > i) ? m_vScannerPhotEntryMomDir[i] : -1) << " \n"; G4cout << "\n"; G4cout << " Photon - Scanner Number of Rayleigh Interactions: " << (((G4int) m_nScannerPhotNumRaylInter.size() > i) ? m_nScannerPhotNumRaylInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Number of Compton Interactions: " << (((G4int) m_nScannerPhotNumCompInter.size() > i) ? m_nScannerPhotNumCompInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Number of Photoelectric Interactions: " << (((G4int) m_nScannerPhotNumPhotoInter.size() > i) ? m_nScannerPhotNumPhotoInter[i] : -1) << " \n"; G4cout << " Photon - Scanner Number of Total Interactions: " << (((G4int) m_nScannerPhotNumTotalInter.size() > i) ? m_nScannerPhotNumTotalInter[i] : -1) << " \n"; } G4cout << "\n\n\n\n"; } } /* ******************************************** */ /* **** Virtual protected member functions **** */ /* ******************************************** */ // ProcessHits function G4bool RPC_SD::ProcessHits (G4Step* thisStep, G4TouchableHistory*) { const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); G4Track* thisTrack = thisStep->GetTrack(); G4int trackID = thisTrack->GetTrackID(); G4int parentID = thisTrack->GetParentID(); G4String particleName = thisStep->GetTrack()->GetDefinition()->GetParticleName(); G4bool hitProcessed = true; // Verify if this is the first time of the trackID in the sensitive detector, and it is, store the trackID and the parentID if (GetAllTrackIndex(trackID) < 0) { m_nAllTracksID.push_back(trackID); m_nAllParentsID.push_back(parentID); } if ((particleName == "gamma") || (particleName == "geantino")) hitProcessed = ProcessGammaHits(thisStep); else if (particleName == "e-" || (particleName == "chargedgeantino")) hitProcessed = ProcessElectronHits(thisStep); return hitProcessed; } /* ******************************************** */ /* **** General protected member functions **** */ /* ******************************************** */ // ProcessGammaHits function G4bool RPC_SD::ProcessGammaHits (G4Step* thisStep) { const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); G4Track* thisTrack = thisStep->GetTrack(); G4int trackID = thisTrack->GetTrackID(); G4int parentID = thisTrack->GetParentID(); G4int photonTrackIndex = GetPhotonTrackIndex(trackID); const G4StepPoint* preStepPoint = thisStep->GetPreStepPoint(); const G4TouchableHandle& preTouchable = preStepPoint->GetTouchableHandle(); const G4VPhysicalVolume* prePhysicalVolume = preTouchable->GetVolume(); G4StepStatus preStepStatus = (prePhysicalVolume == 0) ? fWorldBoundary : preStepPoint->GetStepStatus(); G4String preDepthName = (preStepStatus != fWorldBoundary) ? prePhysicalVolume->GetName() : G4String("unknown"); const G4VProcess* preProcess = preStepPoint->GetProcessDefinedStep(); G4String preProcessName = (preProcess != 0) ? preProcess->GetProcessName() : G4String("unknown"); const G4int preLayerID = preTouchable->GetReplicaNumber(0); const G4StepPoint* postStepPoint = thisStep->GetPostStepPoint(); const G4TouchableHandle& postTouchable = postStepPoint->GetTouchableHandle(); const G4VPhysicalVolume* postPhysicalVolume = postTouchable->GetVolume(); G4StepStatus postStepStatus = (postPhysicalVolume == 0) ? fWorldBoundary : postStepPoint->GetStepStatus(); G4String postDepthName = (postStepStatus != fWorldBoundary) ? postPhysicalVolume->GetName() : G4String("unknown"); const G4VProcess* postProcess = postStepPoint->GetProcessDefinedStep(); G4String postProcessName = (postProcess != 0) ? postProcess->GetProcessName() : G4String("unknown"); const G4int postLayerID = postTouchable->GetReplicaNumber(0); // Just for testing #ifdef DEBUG G4cout << "\n"; G4cout << "\n"; G4cout << " eventID: " << eventID << "\n"; G4cout << " trackID: " << trackID << "\n"; G4cout << " parentID: " << parentID << "\n"; G4cout << "photonTrackIndex: " << photonTrackIndex << "\n"; G4cout << "\n"; G4cout << " preLayerID: " << preLayerID << "\n"; G4cout << " preDepthName: " << preDepthName << "\n"; G4cout << "\n"; G4cout << " postLayerID: " << postLayerID << "\n"; G4cout << " postDepthName: " << postDepthName << "\n"; G4cout << "\n"; G4cout << "\n"; #endif // DEBUG // Initialize the variables for the current trackID, if this is the first time the photon enters the scanner sensitive detector, if (photonTrackIndex < 0) { // Set the photon track index to m_nNumPhotons and then increment the m_nNumPhotons that triggered the scanner sensitive detector photonTrackIndex = m_nNumPhotons++; // Get the phantom sensitive detector G4int numPhantoms = PhantomFactory::GetPhantomFactory()->GetNumberOfPhantoms(); const Phantom_SD* thePhantomSD = (numPhantoms > 0) ? Phantom_SD::GetSensitiveDetector() : 0; // Set the annihilation position m_vAnnihilationPosition = (thePhantomSD != NULL) ? thePhantomSD->GetAnnihilationPosition() : m_vSourcePosition; // Store photon track and parent IDs, and set the photon track index to the newly created entry m_nPhotNumHits.push_back(0); m_nPhotParentID.push_back(parentID); m_nPhotTrackID.push_back(trackID); photonTrackIndex = m_nPhotTrackID.size()-1; // Store the number of interactions of each kind that the photon has suffered in the phantom G4int phantPhotTrackIndex = (thePhantomSD != NULL) ? thePhantomSD->GetPhotonTrackIndex(trackID) : -1; m_nPhantomPhotNumRaylInter.push_back((thePhantomSD != NULL) ? thePhantomSD->GetNumberOfRayleighInteractions(phantPhotTrackIndex) : -1); m_nPhantomPhotNumCompInter.push_back((thePhantomSD != NULL) ? thePhantomSD->GetNumberOfComptonInteractions(phantPhotTrackIndex) : -1); m_nPhantomPhotNumTotalInter.push_back((thePhantomSD != NULL) ? thePhantomSD->GetNumberOfTotlaInteractions(phantPhotTrackIndex) : -1); // Store the energy, position and momentum direction of the photon when it enters in the scanner m_fScannerPhotEntryEnergy.push_back(preStepPoint->GetKineticEnergy()); m_vScannerPhotEntryPos.push_back(preStepPoint->GetPosition()); m_vScannerPhotEntryMomDir.push_back(preStepPoint->GetMomentumDirection()); // Set to zero the number of each kind of interaction of the photon in the scanner m_nScannerPhotNumRaylInter.push_back(0); m_nScannerPhotNumCompInter.push_back(0); m_nScannerPhotNumPhotoInter.push_back(0); m_nScannerPhotNumTotalInter.push_back(0); } // If an interaction occurs, increment the appropriate counter if (postProcessName.contains("Rayleigh")) { m_nScannerPhotNumRaylInter[photonTrackIndex]++; m_nScannerPhotNumTotalInter[photonTrackIndex]++; } else if (postProcessName.contains("ompt")) { m_nScannerPhotNumCompInter[photonTrackIndex]++; m_nScannerPhotNumTotalInter[photonTrackIndex]++; } else if (postProcessName.contains("hot")) { m_nScannerPhotNumPhotoInter[photonTrackIndex]++; m_nScannerPhotNumTotalInter[photonTrackIndex]++; } else if (!postProcessName.contains("Transportation")) m_nScannerPhotNumTotalInter[photonTrackIndex] += ((fabs(thisStep->GetDeltaEnergy()) > 0.0) || (fabs(thisStep->GetTotalEnergyDeposit()) > 0.0) || (preStepPoint->GetMomentumDirection() != postStepPoint->GetMomentumDirection())) ? 1 : 0; // Create photon hits if (!postProcessName.contains("Transportation")) { // Get the creator process name const G4VProcess* creatorProcess = thisTrack->GetCreatorProcess(); const G4String creatorProcessName = (creatorProcess != NULL) ? creatorProcess->GetProcessName() : "primaryParticle"; // Create the primary photon hit if (creatorProcessName.contains("nnih") || creatorProcessName.contains("decay") || (creatorProcessName == "primaryParticle")) CreatePrimaryPhotonHit(thisStep); // Create the photon hit CreateAllPhotonHit(thisStep); } // Suspend the track in order to first process secondary particles if (thisStep->GetSecondary()->size() > 0) thisTrack->SetTrackStatus(fSuspend); return true; } // ProcessElectronHits function G4bool RPC_SD::ProcessElectronHits (G4Step* thisStep) { const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); const G4int eventID = thisEvent->GetEventID(); const G4Track* thisTrack = thisStep->GetTrack(); const G4int trackID = thisTrack->GetTrackID(); const G4StepPoint* preStepPoint = thisStep->GetPreStepPoint(); const G4TouchableHandle& preTouchable = preStepPoint->GetTouchableHandle(); const G4VPhysicalVolume* prePhysicalVolume = preTouchable->GetVolume(0); const G4StepStatus preStepStatus = (prePhysicalVolume == 0) ? fWorldBoundary : preStepPoint->GetStepStatus(); const G4String preDepthName = (preStepStatus != fWorldBoundary) ? prePhysicalVolume->GetName() : G4String("unknown"); const G4int preLayerID = preTouchable->GetReplicaNumber(0); const G4StepPoint* postStepPoint = thisStep->GetPostStepPoint(); const G4TouchableHandle& postTouchable = postStepPoint->GetTouchableHandle(); const G4int postHistoryDepth = postTouchable->GetHistoryDepth(); const G4VPhysicalVolume* postPhysicalVolume = postTouchable->GetVolume(0); const G4StepStatus postStepStatus = (postPhysicalVolume == 0) ? fWorldBoundary : postStepPoint->GetStepStatus(); const G4String postDepthName = (postStepStatus != fWorldBoundary) ? postPhysicalVolume->GetName() : G4String("unknown"); const G4int postLayerID = postTouchable->GetReplicaNumber(0); // Further process electrons that are at the boundary of, or were created in, an RPC layer if (preDepthName.contains("Layer") && (preStepStatus == fUndefined || preStepStatus == fGeomBoundary)) { G4int layerID = preTouchable->GetCopyNumber(0); const RPC_LayerParameters* preLayerParameters = RPC_Parameters::GetParameters()->GetModuleParameters()->GetLayer(layerID); if (preLayerParameters->IsSensitive()) { G4int photonTrackIndex = GetAllTracksPhotonIndex(trackID); if ((photonTrackIndex >= 0) && (photonTrackIndex < GetNumberOfPhotonTracks())) m_nPhotNumHits[photonTrackIndex]++; if (CreateElectronHit(thisStep)) { m_nNumberOfAcceptedEvents++; CreateRPCHit(thisStep); } else m_nNumberOfRejectedEvents++; } } if (thisStep->GetSecondary()->size() > 0) thisStep->GetTrack()->SetTrackStatus(fSuspend); return true; } // CreatePrimaryPhotonHit function G4bool RPC_SD::CreatePrimaryPhotonHit (G4Step* thisStep) { const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); const G4Track* thisTrack = thisStep->GetTrack(); const G4int trackID = thisTrack->GetTrackID(); const G4int photonTrackIndex = GetPhotonTrackIndex(trackID); const G4StepPoint* preStepPoint = thisStep->GetPreStepPoint(); const G4TouchableHandle& preTouchable = preStepPoint->GetTouchableHandle(); const G4String preDepthName = preTouchable->GetVolume(0)->GetName(); const G4StepPoint* postStepPoint = thisStep->GetPostStepPoint(); const G4TouchableHandle& postTouchable = postStepPoint->GetTouchableHandle(); const G4String postDepthName = postTouchable->GetVolume(0)->GetName(); RPC_PrimaryPhotonHit* pPrimPhotHit = new RPC_PrimaryPhotonHit; G4int layerID = -1; G4int moduleID = -1; G4int detectorID = -1; G4int headID = -1; // Get the IDs of the geometry if (preDepthName.contains("Head_Case")) headID = preTouchable->GetReplicaNumber(0); else if (preDepthName.contains("Head_Fill")) headID = preTouchable->GetReplicaNumber(1); else if (preDepthName.contains("Detector")) { detectorID = preTouchable->GetReplicaNumber(0); headID = preTouchable->GetReplicaNumber(2); } else if (preDepthName.contains("Module") && !preDepthName.contains("Layer")) { moduleID = preTouchable->GetReplicaNumber(0); detectorID = preTouchable->GetReplicaNumber(1); headID = preTouchable->GetReplicaNumber(3); } else if (preDepthName.contains("Layer")) { layerID = preTouchable->GetReplicaNumber(0); moduleID = preTouchable->GetReplicaNumber(1); detectorID = preTouchable->GetReplicaNumber(2); headID = preTouchable->GetReplicaNumber(4); } // Create the primary photon hit pPrimPhotHit->SetEventID (G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID()+1); pPrimPhotHit->SetSourceID (m_nSourceID); pPrimPhotHit->SetSourcePosition (m_vSourcePosition); pPrimPhotHit->SetAnnihilationPosition (m_vAnnihilationPosition); pPrimPhotHit->SetParentID (GetPhotonParentID(photonTrackIndex)); pPrimPhotHit->SetTrackID (GetPhotonTrackID(photonTrackIndex)); pPrimPhotHit->SetLayerID (layerID); pPrimPhotHit->SetModuleID (moduleID); pPrimPhotHit->SetDetectorID (detectorID); pPrimPhotHit->SetHeadID (headID); pPrimPhotHit->SetPhantomRayleighInteractions (GetPhantomNumberOfRayleighInteractions(photonTrackIndex)); pPrimPhotHit->SetPhantomComptonInteractions (GetPhantomNumberOfComptonInteractions(photonTrackIndex)); pPrimPhotHit->SetPhantomTotalInteractions (GetPhantomNumberOfTotalInteractions(photonTrackIndex)); pPrimPhotHit->SetScannerRayleighInteractions (GetScannerNumberOfRayleighInteractions(photonTrackIndex)); pPrimPhotHit->SetScannerComptonInteractions (GetScannerNumberOfComptonInteractions(photonTrackIndex)); pPrimPhotHit->SetScannerPhotoelectricInteractions (GetScannerNumberOfPhotoelectricInteractions(photonTrackIndex)); pPrimPhotHit->SetScannerTotalInteractions (GetScannerNumberOfTotalInteractions(photonTrackIndex)); pPrimPhotHit->SetGlobalTime (thisTrack->GetGlobalTime()); pPrimPhotHit->SetEnergy (preStepPoint->GetKineticEnergy()); pPrimPhotHit->SetEnergyTransfer (thisStep->GetDeltaEnergy()); pPrimPhotHit->SetEnergyDeposit (thisStep->GetTotalEnergyDeposit()); pPrimPhotHit->SetGlobalPosition (thisTrack->GetPosition()); pPrimPhotHit->SetIncomingMomentumDirection (preStepPoint->GetMomentumDirection()); pPrimPhotHit->SetOutgoingMomentumDirection (postStepPoint->GetMomentumDirection()); m_pPrimPhotHitsCollection->insert(pPrimPhotHit); return true; } // CreateAllPhotonHit function G4bool RPC_SD::CreateAllPhotonHit (G4Step* thisStep) { const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); const G4Track* thisTrack = thisStep->GetTrack(); const G4StepPoint* preStepPoint = thisStep->GetPreStepPoint(); const G4TouchableHandle& preTouchable = preStepPoint->GetTouchableHandle(); const G4String preDepthName = preTouchable->GetVolume(0)->GetName(); const G4StepPoint* postStepPoint = thisStep->GetPostStepPoint(); const G4TouchableHandle& postTouchable = postStepPoint->GetTouchableHandle(); const G4String postDepthName = postTouchable->GetVolume(0)->GetName(); RPC_AllPhotonHit* pAllPhotHit = new RPC_AllPhotonHit; G4int layerID = -1; G4int moduleID = -1; G4int detectorID = -1; G4int headID = -1; // Get the IDs of the geometry if (preDepthName.contains("Head_Case")) headID = preTouchable->GetReplicaNumber(0); else if (preDepthName.contains("Head_Fill")) headID = preTouchable->GetReplicaNumber(1); else if (preDepthName.contains("Detector")) { detectorID = preTouchable->GetReplicaNumber(0); headID = preTouchable->GetReplicaNumber(2); } else if (preDepthName.contains("Module") && !preDepthName.contains("Layer")) { moduleID = preTouchable->GetReplicaNumber(0); detectorID = preTouchable->GetReplicaNumber(1); headID = preTouchable->GetReplicaNumber(3); } else if (preDepthName.contains("Layer")) { layerID = preTouchable->GetReplicaNumber(0); moduleID = preTouchable->GetReplicaNumber(1); detectorID = preTouchable->GetReplicaNumber(2); headID = preTouchable->GetReplicaNumber(4); } // Create the photon hit pAllPhotHit->SetEventID (G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID()+1); pAllPhotHit->SetSourceID (m_nSourceID); pAllPhotHit->SetSourcePosition (m_vSourcePosition); pAllPhotHit->SetAnnihilationPosition (m_vAnnihilationPosition); pAllPhotHit->SetParentID (thisTrack->GetParentID()); pAllPhotHit->SetTrackID (thisTrack->GetTrackID()); pAllPhotHit->SetLayerID (layerID); pAllPhotHit->SetModuleID (moduleID); pAllPhotHit->SetDetectorID (detectorID); pAllPhotHit->SetHeadID (headID); pAllPhotHit->SetGlobalTime (thisTrack->GetGlobalTime()); pAllPhotHit->SetEnergy (preStepPoint->GetKineticEnergy()); pAllPhotHit->SetEnergyTransfer (thisStep->GetDeltaEnergy()); pAllPhotHit->SetEnergyDeposit (thisStep->GetTotalEnergyDeposit()); pAllPhotHit->SetGlobalPosition (thisTrack->GetPosition()); pAllPhotHit->SetIncomingMomentumDirection (preStepPoint->GetMomentumDirection()); pAllPhotHit->SetOutgoingMomentumDirection (postStepPoint->GetMomentumDirection()); m_pAllPhotHitsCollection->insert(pAllPhotHit); return(true); } // CreateElectronHit function G4bool RPC_SD::CreateElectronHit (G4Step* thisStep) { // Get data from RunAction const RPC_RunAction* pRunAction = dynamic_cast ((G4RunManager::GetRunManager())->GetUserRunAction()); G4bool bAcceptAllEvents = pRunAction->GetAcceptAllGeneratedEvents(); G4String sDetectionPointComputationMethod = pRunAction->GetDetectionPointComputationMethod(); G4double fActiveLayerFractionFromCathode = pRunAction->GetActiveGapFractionFromCathode(); G4String sAssignmentPointElement = pRunAction->GetAssignmentPointElement(); // Get event and track data const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); const G4Track* thisTrack = thisStep->GetTrack(); const G4int trackID = thisTrack->GetTrackID(); const G4int parentID = thisTrack->GetParentID(); // Get layer parameters from the pre step point const G4StepPoint* preStepPoint = thisStep->GetPreStepPoint(); const G4TouchableHandle& preTouchable = preStepPoint->GetTouchableHandle(); const RPC_LayerParameters* preLayerParameters = RPC_Parameters::GetParameters()->GetModuleParameters()->GetLayer(preTouchable->GetCopyNumber(0)); const G4double fLayerThickness = preLayerParameters->GetThickness()/mm; // NOTICE: fLayerThickness is in mm. // Get the IDs of the layer, module, detector and head of the sensitive layer where the electron is detected G4int layerID = preTouchable->GetCopyNumber(0); G4int moduleID = preTouchable->GetCopyNumber(1); G4int detectorID = preTouchable->GetCopyNumber(2); G4int headID = preTouchable->GetCopyNumber(4); // Get the extraction position and momentum direction in the global coordinate system G4ThreeVector vGlobalExtPos = preStepPoint->GetPosition(); G4ThreeVector vGlobalExtMomDir = preStepPoint->GetMomentumDirection(); G4ThreeVector vGlobalDetPos; // Get the layer coordinates transformation matrix and the extraction position and momentum direction in the layer's coordinate system, // and create the local point to hold the detection point const G4AffineTransform& headTransform = preTouchable->GetHistory()->GetTransform(preTouchable->GetHistoryDepth()-3); G4ThreeVector vHeadExtPos = headTransform.TransformPoint(vGlobalExtPos); G4ThreeVector vHeadExtMomDir = headTransform.TransformAxis(vGlobalExtMomDir); const G4AffineTransform& detectorTransform = preTouchable->GetHistory()->GetTransform(preTouchable->GetHistoryDepth()-2); G4ThreeVector vDetectorExtPos = detectorTransform.TransformPoint(vGlobalExtPos); G4ThreeVector vDetectorExtMomDir = detectorTransform.TransformAxis(vGlobalExtMomDir); const G4AffineTransform& moduleTransform = preTouchable->GetHistory()->GetTransform(preTouchable->GetHistoryDepth()-1); G4ThreeVector vModuleExtPos = moduleTransform.TransformPoint(vGlobalExtPos); G4ThreeVector vModuleExtMomDir = moduleTransform.TransformAxis(vGlobalExtMomDir); const G4AffineTransform& layerTransform = preTouchable->GetHistory()->GetTopTransform(); G4ThreeVector vLayerExtPos = layerTransform.TransformPoint(vGlobalExtPos); G4ThreeVector vLayerExtMomDir = layerTransform.TransformAxis(vGlobalExtMomDir); G4ThreeVector vLayerDetPos; // Create a variable to hold the electron flight time G4double fElectronFlightTime = 0.0; // Reject the event if all the events are not to be accepted and the electron was extracted outside the active region with a // linear momentum along the radial direction from the cathode to the anode if ((vLayerExtMomDir.getX() >= 0.0) && (vLayerExtPos.getX() > (fActiveLayerFractionFromCathode-1)*0.5*fLayerThickness) && !bAcceptAllEvents) return false; // Get the voltage drop and electric field of the sensitive layer where the electron was extracted to G4double fLayerVoltageDrop = preLayerParameters->GetVoltageDrop()/kilovolt; // NOTICE: fLayerVoltageDrop is in kV. G4double fLayerElectricField = (fabs(fLayerThickness) > 0) ? fLayerVoltageDrop/fLayerThickness : 0; // NOTICE: fLayerElectricField is in kV/mm, meaning that e*fGapElectricField is in keV/mm. // If the detection point computation method is set to "Automatic", set it to one of the "AnodeSideWithElectricField" or "AnodeSideWithoutElectricField", depending on the existence of an electric field in the layer if (sDetectionPointComputationMethod == "Automatic") sDetectionPointComputationMethod = (fLayerElectricField > 0) ? "AnodeSideWithElectricField" : "AnodeSideWithoutElectricField"; // Compute the detection point without electric field if (sDetectionPointComputationMethod.contains("WithoutElectricField") || (fLayerElectricField == 0)) { if (vLayerExtMomDir.getX() == 0) { vLayerDetPos = vLayerExtPos; vLayerDetPos.setX(0.5*fLayerThickness); } else { vLayerDetPos = vLayerExtPos + vLayerExtMomDir*(((vLayerExtMomDir.getX() > 0) ? 1 : -1)*0.5*fLayerThickness-vLayerExtPos.getX())/vLayerExtMomDir.getX(); fElectronFlightTime = (vLayerDetPos.getX()-vLayerExtPos.getX())/(preStepPoint->GetVelocity()*vLayerExtMomDir.getX()); } } // Compute the detection point with electric field else if (sDetectionPointComputationMethod.contains("WithElectricField")) { G4double fInitialEnergy = preStepPoint->GetTotalEnergy()/keV; // NOTICE: The (total) initial energy is in keV. G4ThreeVector vInitialMomentum = layerTransform.TransformAxis(preStepPoint->GetMomentum()/keV); // NOTICE: Geant4 gives the linear momentum vector in units of keV/c. So, when using the momentum in an equation that involves G4double fInversionTime = -vInitialMomentum.getX()/(eplus*fLayerElectricField*c_light); // the product of the linear momentum by the light speed (p*c), there is no need to multiply the Geant4 linear momentum by c. G4double fFinalTime = -1; // However, if the momentum is used in an equation in which p appears isolated, one must divide it by c. G4double fStartPosX = vLayerExtPos.getX(); G4double fEndPosX; G4double fMinX; // Set the minimum and final position along the x direction if (fInversionTime <= 0) { fMinX = vLayerExtPos.getX(); fEndPosX = 0.5*fLayerThickness; } else { fMinX = vLayerExtPos.getX() + (sqrt(pow(fInitialEnergy,2)-pow(vInitialMomentum.getX(),2)) - fInitialEnergy)/(eplus*fLayerElectricField); if (fMinX <= -0.5*fLayerThickness) fEndPosX = -0.5*fLayerThickness; else if (sDetectionPointComputationMethod.contains("OppositeSide")) { fEndPosX = fMinX; fFinalTime = fInversionTime; } else fEndPosX = 0.5*fLayerThickness; } if (fMinX > (fActiveLayerFractionFromCathode-0.5)*fLayerThickness && !bAcceptAllEvents) return false; // Compute the time needed for the electron to reach the final position if (fFinalTime < 0) { G4double fDelta = pow(vInitialMomentum.getX(),2) + eplus*fLayerElectricField*(fEndPosX-fStartPosX)*(eplus*fLayerElectricField*(fEndPosX-fStartPosX)+2*fInitialEnergy); if (fDelta < 0) return false; if ((fFinalTime = ((vInitialMomentum.getX() >= 0 ? +1 : ((fMinX > -0.5*fLayerThickness) ? +1 : -1))*sqrt(fDelta)-vInitialMomentum.getX())/(eplus*fLayerElectricField*c_light)) < 0) return false; } // Compute the final position of the electron G4double yzLayerDetPosCommonDen = sqrt(pow(fInitialEnergy,2)-pow(vInitialMomentum.getX(),2)); G4double yzLayerDetPosCommon = (asinh((eplus*fLayerElectricField*c_light*fFinalTime+vInitialMomentum.getX())/yzLayerDetPosCommonDen) - asinh(vInitialMomentum.getX()/yzLayerDetPosCommonDen))/(eplus*fLayerElectricField); G4double xLayerDetPos = vLayerExtPos.getX() + (sqrt(pow(fInitialEnergy,2) + eplus*fLayerElectricField*c_light*fFinalTime*(eplus*fLayerElectricField*c_light*fFinalTime+2*vInitialMomentum.getX())) - fInitialEnergy)/(eplus*fLayerElectricField); G4double yLayerDetPos = vLayerExtPos.getY() + vInitialMomentum.getY()*yzLayerDetPosCommon; G4double zLayerDetPos = vLayerExtPos.getZ() + vInitialMomentum.getZ()*yzLayerDetPosCommon; vLayerDetPos.setX(xLayerDetPos); vLayerDetPos.setY(yLayerDetPos); vLayerDetPos.setZ(zLayerDetPos); // Set the electron flight time fElectronFlightTime = fFinalTime; } if (sDetectionPointComputationMethod.contains("AnodeSide")) vLayerDetPos.setX(0.5*fLayerThickness); if ((preTouchable->GetSolid(0)->Inside(vLayerDetPos) == kOutside) && !bAcceptAllEvents) return false; vGlobalDetPos = layerTransform.Inverse().TransformPoint(vLayerDetPos); // Compute the assignment position G4int elementID = preTouchable->GetHistoryDepth(); if (sAssignmentPointElement == "Module") elementID -= 1; else if (sAssignmentPointElement == "Detector") elementID -= 2; const G4AffineTransform& elementTransform = preTouchable->GetHistory()->GetTransform(elementID); G4ThreeVector vElementDetPos = elementTransform.TransformPoint(vGlobalDetPos); G4ThreeVector vElementAssPos = G4ThreeVector(0.0*mm, vElementDetPos.getY(), vElementDetPos.getZ()); G4ThreeVector vGlobalAssPos = elementTransform.Inverse().TransformPoint(vElementAssPos); // **** Create the electron hit // Get the index of the primary photon hits collection that contains the last photon hit that has given rise to the current trackID G4int numPrimPhotHits = m_pPrimPhotHitsCollection->entries(); RPC_PrimaryPhotonHit* pPrimPhotHit = (*m_pPrimPhotHitsCollection)[numPrimPhotHits-1]; G4int numAllPhotHits = m_pAllPhotHitsCollection->entries(); RPC_AllPhotonHit* pAllPhotHit = (*m_pAllPhotHitsCollection)[numAllPhotHits-1]; G4int nAllPhotIndex = GetPhotonTrackIndex(pAllPhotHit->GetTrackID()); // Create the electron hit RPC_ElectronHit* pElectHit = new RPC_ElectronHit; // General data pElectHit->SetEventID (thisEvent->GetEventID()+1); // Source data pElectHit->SetSourceID (m_nSourceID); pElectHit->SetSourcePosition (m_vSourcePosition); pElectHit->SetAnnihilationPosition (m_vAnnihilationPosition); // Primary Photon Data if (pPrimPhotHit != 0) { pElectHit->SetPrimaryPhotonParentID (pPrimPhotHit->GetParentID()); pElectHit->SetPrimaryPhotonTrackID (pPrimPhotHit->GetTrackID()); pElectHit->SetPrimaryPhotonLayerID (pPrimPhotHit->GetLayerID()); pElectHit->SetPrimaryPhotonModuleID (pPrimPhotHit->GetModuleID()); pElectHit->SetPrimaryPhotonDetectorID (pPrimPhotHit->GetDetectorID()); pElectHit->SetPrimaryPhotonHeadID (pPrimPhotHit->GetHeadID()); pElectHit->SetPrimaryPhotonPhantomRayleighInteractions (pPrimPhotHit->GetPhantomRayleighInteractions()); pElectHit->SetPrimaryPhotonPhantomComptonInteractions (pPrimPhotHit->GetPhantomComptonInteractions()); pElectHit->SetPrimaryPhotonPhantomTotalInteractions (pPrimPhotHit->GetPhantomTotalInteractions()); pElectHit->SetPrimaryPhotonScannerRayleighInteractions (pPrimPhotHit->GetScannerRayleighInteractions()); pElectHit->SetPrimaryPhotonScannerComptonInteractions (pPrimPhotHit->GetScannerComptonInteractions()); pElectHit->SetPrimaryPhotonScannerPhotoelectricInteractions (pPrimPhotHit->GetScannerPhotoelectricInteractions()); pElectHit->SetPrimaryPhotonScannerTotalInteractions (pPrimPhotHit->GetScannerTotalInteractions()); pElectHit->SetPrimaryPhotonGlobalTime (pPrimPhotHit->GetGlobalTime()); pElectHit->SetPrimaryPhotonEnergy (pPrimPhotHit->GetEnergy()); pElectHit->SetPrimaryPhotonEnergyTransfer (pPrimPhotHit->GetEnergyTransfer()); pElectHit->SetPrimaryPhotonEnergyDeposit (pPrimPhotHit->GetEnergyDeposit()); pElectHit->SetPrimaryPhotonGlobalPosition (pPrimPhotHit->GetGlobalPosition()); pElectHit->SetPrimaryPhotonIncomingMomentumDirection (pPrimPhotHit->GetIncomingMomentumDirection()); pElectHit->SetPrimaryPhotonOutgoingMomentumDirection (pPrimPhotHit->GetOutgoingMomentumDirection()); } // All Photon Data if (pAllPhotHit != 0) { pElectHit->SetPhotonParentID (pAllPhotHit->GetParentID()); pElectHit->SetPhotonTrackID (pAllPhotHit->GetTrackID()); pElectHit->SetPhotonLayerID (pAllPhotHit->GetLayerID()); pElectHit->SetPhotonModuleID (pAllPhotHit->GetModuleID()); pElectHit->SetPhotonDetectorID (pAllPhotHit->GetDetectorID()); pElectHit->SetPhotonHeadID (pAllPhotHit->GetHeadID()); pElectHit->SetPhotonGlobalTime (pAllPhotHit->GetGlobalTime()); pElectHit->SetPhotonEnergy (pAllPhotHit->GetEnergy()); pElectHit->SetPhotonEnergyTransfer (pAllPhotHit->GetEnergyTransfer()); pElectHit->SetPhotonEnergyDeposit (pAllPhotHit->GetEnergyDeposit()); pElectHit->SetPhotonGlobalPosition (pAllPhotHit->GetGlobalPosition()); pElectHit->SetPhotonIncomingMomentumDirection (pAllPhotHit->GetIncomingMomentumDirection()); pElectHit->SetPhotonOutgoingMomentumDirection (pAllPhotHit->GetOutgoingMomentumDirection()); } // Electron Data pElectHit->SetElectronParentID (parentID); pElectHit->SetElectronTrackID (trackID); pElectHit->SetElectronHitID (GetPhotonNumberOfHits(nAllPhotIndex)); pElectHit->SetElectronLayerID (layerID); pElectHit->SetElectronModuleID (moduleID); pElectHit->SetElectronDetectorID (detectorID); pElectHit->SetElectronHeadID (headID); pElectHit->SetElectronExtractionTime (preStepPoint->GetGlobalTime()+fElectronFlightTime); // pElectHit->SetElectronExtractionTime (theStepPoint->GetGlobalTime()); pElectHit->SetElectronExtractionEnergy (preStepPoint->GetKineticEnergy()); pElectHit->SetElectronExtractionPosition (vGlobalExtPos); pElectHit->SetElectronDetectionPosition (vGlobalDetPos); pElectHit->SetElectronAssignmentPosition (vGlobalAssPos); pElectHit->SetElectronMomentumDirection (vGlobalExtMomDir); m_pElectHitsCollection->insert(pElectHit); return true; } // CreateRPCHit function G4bool RPC_SD::CreateRPCHit (G4Step* thisStep) { // Get event and track data const G4Event* thisEvent = G4RunManager::GetRunManager()->GetCurrentEvent(); G4int eventID = thisEvent->GetEventID(); const G4Track* thisTrack = thisStep->GetTrack(); const G4int trackID = thisTrack->GetTrackID(); const G4int parentID = thisTrack->GetParentID(); G4int numPrimPhotHits = m_pPrimPhotHitsCollection->entries(); RPC_PrimaryPhotonHit* pPrimPhotFirstHit = (numPrimPhotHits > 0) ? (*m_pPrimPhotHitsCollection)[0] : 0; RPC_PrimaryPhotonHit* pPrimPhotLastHit = (numPrimPhotHits > 0) ? (*m_pPrimPhotHitsCollection)[numPrimPhotHits-1] : nullptr; G4int numAllPhotHits = m_pAllPhotHitsCollection->entries(); RPC_AllPhotonHit* pAllPhotHit = (numAllPhotHits > 0) ? (*m_pAllPhotHitsCollection)[numAllPhotHits-1] : nullptr; G4int numElectHits = m_pElectHitsCollection->entries(); RPC_ElectronHit* pElectHit = (numElectHits > 0) ? (*m_pElectHitsCollection)[numElectHits-1] : nullptr; // If any of the pointers to the hits are null, do not process further the current hit if ((pPrimPhotFirstHit == nullptr) || (pPrimPhotLastHit == nullptr) || (pAllPhotHit == nullptr) || (pElectHit == nullptr)) return false; // Create the RPC hit RPC_Hit* pRPCHit = new RPC_Hit; // General data pRPCHit->SetEventID (eventID+1); // Phantom data pRPCHit->SetSourceID (m_nSourceID); pRPCHit->SetSourcePosition (m_vSourcePosition); pRPCHit->SetAnnihilationPosition (m_vAnnihilationPosition); pRPCHit->SetPhotonParentID (pAllPhotHit->GetParentID()); pRPCHit->SetPhotonID (pAllPhotHit->GetTrackID()); pRPCHit->SetPrimaryPhotonPhantomNumberOfRayleighInteractions (pPrimPhotLastHit->GetPhantomRayleighInteractions()); pRPCHit->SetPrimaryPhotonPhantomNumberOfComptonInteractions (pPrimPhotLastHit->GetPhantomComptonInteractions()); pRPCHit->SetPrimaryPhotonPhantomNumberOfTotalInteractions (pPrimPhotLastHit->GetPhantomTotalInteractions()); // Scanner data pRPCHit->SetPrimaryPhotonScannerEntryEnergy (pPrimPhotFirstHit->GetEnergy()); pRPCHit->SetPrimaryPhotonScannerEntryPosition (pPrimPhotFirstHit->GetGlobalPosition()); pRPCHit->SetPrimaryPhotonScannerEntryMomentumDirection (pPrimPhotFirstHit->GetIncomingMomentumDirection()); pRPCHit->SetPrimaryPhotonScannerNumberOfRayleighInteractions (pPrimPhotLastHit->GetScannerRayleighInteractions()); pRPCHit->SetPrimaryPhotonScannerNumberOfComptonInteractions (pPrimPhotLastHit->GetScannerComptonInteractions()); pRPCHit->SetPrimaryPhotonScannerNumberOfPhotoelectricInteractions (pPrimPhotLastHit->GetScannerPhotoelectricInteractions()); pRPCHit->SetPrimaryPhotonScannerNumberOfTotalInteractions (pPrimPhotLastHit->GetScannerTotalInteractions()); pRPCHit->SetPrimaryPhotonScannerDetectionTime (pPrimPhotLastHit->GetGlobalTime()); pRPCHit->SetPrimaryPhotonScannerDetectionEnergy (pPrimPhotLastHit->GetEnergy()); pRPCHit->SetPrimaryPhotonScannerDetectionEnergyTransfer (pPrimPhotLastHit->GetEnergyTransfer()); pRPCHit->SetPrimaryPhotonScannerDetectionEnergyDeposit (pPrimPhotLastHit->GetEnergyDeposit()); pRPCHit->SetPrimaryPhotonScannerDetectionPosition (pPrimPhotLastHit->GetGlobalPosition()); // Electron Data pRPCHit->SetHitID (GetPhotonNumberOfHits(GetPhotonTrackIndex(trackID))); pRPCHit->SetElectronGapID (pElectHit->GetElectronLayerID()); pRPCHit->SetElectronModuleID (pElectHit->GetElectronModuleID()); pRPCHit->SetElectronDetectorID (pElectHit->GetElectronDetectorID()); pRPCHit->SetElectronHeadID (pElectHit->GetElectronHeadID()); pRPCHit->SetElectronExtractionTime (pElectHit->GetElectronExtractionTime()); pRPCHit->SetElectronExtractionEnergy (pElectHit->GetElectronExtractionEnergy()); pRPCHit->SetElectronExtractionPosition (pElectHit->GetElectronExtractionPosition()); pRPCHit->SetElectronDetectionPosition (pElectHit->GetElectronDetectionPosition()); pRPCHit->SetElectronAssignmentPosition (pElectHit->GetElectronAssignmentPosition()); pRPCHit->SetElectronExtractionMomentumDirection (pElectHit->GetElectronMomentumDirection()); m_pRPCHitsCollection->insert(pRPCHit); return true; } /* ******************************************** */ /* **** General public interface functions **** */ /* ******************************************** */ // GetSensitiveDetector function RPC_SD* RPC_SD::GetSensitiveDetector (void) { if (m_pSD == 0) m_pSD = new RPC_SD("/Detector/ProbeVolume"); return m_pSD; } // NumberOfAceptedEvents functions G4long RPC_SD::GetNumberOfAceptedEvents (void) const { return m_nNumberOfAcceptedEvents; } void RPC_SD::ResetNumberOfAcceptedEvents (void) { m_nNumberOfAcceptedEvents = 0; } // NumberOfRejectedEvents functions G4long RPC_SD::GetNumberOfRejectedEvents (void) const { return m_nNumberOfRejectedEvents; } void RPC_SD::ResetNumberOfRejectedEvents (void) { m_nNumberOfRejectedEvents = 0; } // AllTracks functions G4int RPC_SD::GetNumberOfAllTracks (void) const { return m_nNumAllTracks; } G4int RPC_SD::GetAllTracksPhotonIndex (G4int nTrackID, G4int nPrimPhotTrackID) const { G4int nTrackIndex; G4int nOut; G4int nNumCycles = m_nAllTracksID.size(); do { nNumCycles--; nTrackIndex = GetAllTrackIndex(nTrackID); nOut = GetPhotonTrackIndex(m_nAllParentsID[nTrackIndex]); if ((nPrimPhotTrackID > 0) && (nOut >= 0) && (m_nAllParentsID[nTrackIndex] == nPrimPhotTrackID)) break; else if ((nPrimPhotTrackID == 0) && (nOut >= 0)) break; else nTrackID = m_nAllParentsID[nTrackIndex]; } while (((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nAllParentsID.size())) || (nNumCycles > 0)); return nOut; } G4int RPC_SD::GetAllTrackIndex (G4int nTrackID) const { G4int nOut = -1; for (G4int i = 0 ; i < (G4int) m_nAllTracksID.size() ; i++) { if (nTrackID == m_nAllTracksID[i]) { nOut = i; break; } } return nOut; } G4int RPC_SD::GetAllTrackID (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nAllTracksID.size())) ? m_nAllTracksID[nTrackIndex] : -1; } G4int RPC_SD::GetAllParentTrackID (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nAllParentsID.size())) ? m_nAllParentsID[nTrackIndex] : -1; } // Source functions G4int RPC_SD::GetSourceID (void) const { return m_nSourceID; } G4ThreeVector RPC_SD::GetSourcePosition (void) const { return m_vSourcePosition; } G4ThreeVector RPC_SD::GetAnnihilationPosition (void) const { return m_vAnnihilationPosition; } // PhotonTrack functions G4int RPC_SD::GetNumberOfPhotonTracks (void) const { return m_nNumPhotons; } G4int RPC_SD::GetPhotonTrackIndex (G4int nTrackID) const { G4int nOut = -1; for (G4int i = 0 ; i < (G4int) m_nPhotTrackID.size() ; i++) { if (nTrackID == m_nPhotTrackID[i]) { nOut = i; break; } } return nOut; } G4int RPC_SD::GetPhotonParentID (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhotParentID.size())) ? m_nPhotParentID[nTrackIndex] : -1; } G4int RPC_SD::GetPhotonTrackID (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhotTrackID.size())) ? m_nPhotTrackID[nTrackIndex] : -1; } G4int RPC_SD::GetPhotonNumberOfHits (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhotNumHits.size())) ? m_nPhotNumHits[nTrackIndex] : -1; } // OhantomNumberInteractions functions G4int RPC_SD::GetPhantomNumberOfRayleighInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhantomPhotNumRaylInter.size())) ? m_nPhantomPhotNumRaylInter[nTrackIndex] : -1; } G4int RPC_SD::GetPhantomNumberOfComptonInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhantomPhotNumCompInter.size())) ? m_nPhantomPhotNumCompInter[nTrackIndex] : -1; } G4int RPC_SD::GetPhantomNumberOfTotalInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nPhantomPhotNumTotalInter.size())) ? m_nPhantomPhotNumTotalInter[nTrackIndex] : -1; } // ScannerPhotonEntry functions G4double RPC_SD::GetScannerPhotonEntryEnergy (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_fScannerPhotEntryEnergy.size())) ? m_fScannerPhotEntryEnergy[nTrackIndex] : -1; } G4ThreeVector RPC_SD::GetScannerPhotonEntryPosition (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_vScannerPhotEntryPos.size())) ? m_vScannerPhotEntryPos[nTrackIndex] : -1; } G4ThreeVector RPC_SD::GetScannerPhotonEntryMomentumDirection (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_vScannerPhotEntryMomDir.size())) ? m_vScannerPhotEntryMomDir[nTrackIndex] : -1; } // ScannerNumberoOfInteractions functions G4int RPC_SD::GetScannerNumberOfRayleighInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nScannerPhotNumRaylInter.size())) ? m_nScannerPhotNumRaylInter[nTrackIndex] : -1; } G4int RPC_SD::GetScannerNumberOfComptonInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nScannerPhotNumCompInter.size())) ? m_nScannerPhotNumCompInter[nTrackIndex] : -1; } G4int RPC_SD::GetScannerNumberOfPhotoelectricInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nScannerPhotNumPhotoInter.size())) ? m_nScannerPhotNumPhotoInter[nTrackIndex] : -1; } G4int RPC_SD::GetScannerNumberOfTotalInteractions (G4int nTrackIndex) const { return ((nTrackIndex >= 0) && (nTrackIndex < (G4int) m_nScannerPhotNumTotalInter.size())) ? m_nScannerPhotNumTotalInter[nTrackIndex] : -1; }