Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyDetectorSD.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/hadrontherapy/src/HadrontherapyDetectorSD.cc (Version 11.3.0) and /examples/advanced/hadrontherapy/src/HadrontherapyDetectorSD.cc (Version 10.6.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Hadrontherapy advanced example for Geant4       26 // Hadrontherapy advanced example for Geant4
 27 // See more at: https://twiki.cern.ch/twiki/bi     27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
 28                                                    28 
 29 #include "HadrontherapyDetectorSD.hh"              29 #include "HadrontherapyDetectorSD.hh"
 30 #include "HadrontherapyDetectorHit.hh"             30 #include "HadrontherapyDetectorHit.hh"
 31 #include "G4Step.hh"                               31 #include "G4Step.hh"
 32 #include "G4VTouchable.hh"                         32 #include "G4VTouchable.hh"
 33 #include "G4TouchableHistory.hh"                   33 #include "G4TouchableHistory.hh"
 34 #include "G4SDManager.hh"                          34 #include "G4SDManager.hh"
 35 #include "HadrontherapyMatrix.hh"                  35 #include "HadrontherapyMatrix.hh"
 36 #include "HadrontherapyLet.hh"                     36 #include "HadrontherapyLet.hh"
 37 #include "G4Track.hh"                              37 #include "G4Track.hh"
 38 #include "G4SystemOfUnits.hh"                      38 #include "G4SystemOfUnits.hh"
 39 #include "HadrontherapyMatrix.hh"                  39 #include "HadrontherapyMatrix.hh"
 40                                                    40 
 41                                                    41 
 42 #include "G4SteppingManager.hh"                    42 #include "G4SteppingManager.hh"
 43 #include "G4TrackVector.hh"                        43 #include "G4TrackVector.hh"
 44 #include "HadrontherapySteppingAction.hh"          44 #include "HadrontherapySteppingAction.hh"
 45 #include "G4ios.hh"                                45 #include "G4ios.hh"
 46 #include "G4SteppingManager.hh"                    46 #include "G4SteppingManager.hh"
 47 #include "G4Track.hh"                              47 #include "G4Track.hh"
 48 #include "G4Step.hh"                               48 #include "G4Step.hh"
 49 #include "G4StepPoint.hh"                          49 #include "G4StepPoint.hh"
 50 #include "G4TrackStatus.hh"                        50 #include "G4TrackStatus.hh"
 51 #include "G4TrackVector.hh"                        51 #include "G4TrackVector.hh"
 52 #include "G4ParticleDefinition.hh"                 52 #include "G4ParticleDefinition.hh"
 53 #include "G4ParticleTypes.hh"                      53 #include "G4ParticleTypes.hh"
 54 #include "G4UserEventAction.hh"                    54 #include "G4UserEventAction.hh"
 55 #include "G4TransportationManager.hh"              55 #include "G4TransportationManager.hh"
 56 #include "G4VSensitiveDetector.hh"                 56 #include "G4VSensitiveDetector.hh"
 57 #include "HadrontherapyRunAction.hh"               57 #include "HadrontherapyRunAction.hh"
 58 #include "G4SystemOfUnits.hh"                      58 #include "G4SystemOfUnits.hh"
 59 #include "HadrontherapyRBE.hh"                 << 
 60 #include <G4AccumulableManager.hh>             << 
 61                                                    59 
 62                                                    60 
 63 //////////////////////////////////////////////     61 /////////////////////////////////////////////////////////////////////////////
 64 HadrontherapyDetectorSD::HadrontherapyDetector     62 HadrontherapyDetectorSD::HadrontherapyDetectorSD(G4String name):
 65 G4VSensitiveDetector(name)                         63 G4VSensitiveDetector(name)
 66 {                                                  64 {
 67     G4String HCname;                               65     G4String HCname;
 68     collectionName.insert(HCname="Hadrontherap     66     collectionName.insert(HCname="HadrontherapyDetectorHitsCollection");
 69     HitsCollection = NULL;                         67     HitsCollection = NULL;
 70     sensitiveDetectorName = name;                  68     sensitiveDetectorName = name;
 71                                                    69     
 72 }                                                  70 }
 73                                                    71 
 74 //////////////////////////////////////////////     72 /////////////////////////////////////////////////////////////////////////////
 75 HadrontherapyDetectorSD::~HadrontherapyDetecto     73 HadrontherapyDetectorSD::~HadrontherapyDetectorSD()
 76 {}                                                 74 {}
 77                                                    75 
 78 //////////////////////////////////////////////     76 /////////////////////////////////////////////////////////////////////////////
 79 void HadrontherapyDetectorSD::Initialize(G4HCo     77 void HadrontherapyDetectorSD::Initialize(G4HCofThisEvent*)
 80 {                                                  78 {
 81                                                    79     
 82     HitsCollection = new HadrontherapyDetector     80     HitsCollection = new HadrontherapyDetectorHitsCollection(sensitiveDetectorName,
 83                                                    81                                                              collectionName[0]);
 84 }                                                  82 }
 85                                                    83 
 86 //////////////////////////////////////////////     84 /////////////////////////////////////////////////////////////////////////////
 87 G4bool HadrontherapyDetectorSD::ProcessHits(G4     85 G4bool HadrontherapyDetectorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* )
 88 {                                                  86 {
 89                                                    87     
 90                                                    88     
 91     if (aStep -> GetPreStepPoint() -> GetPhysi     89     if (aStep -> GetPreStepPoint() -> GetPhysicalVolume() -> GetName() != "RODetectorZDivisionPhys") return false;
 92                                                    90     
 93                                                    91     
 94     // Get kinetic energy                          92     // Get kinetic energy
 95     G4Track * theTrack = aStep  ->  GetTrack()     93     G4Track * theTrack = aStep  ->  GetTrack();
 96     G4double kineticEnergy = theTrack->GetKine <<  94     
 97     G4ParticleDefinition *particleDef = theTra     95     G4ParticleDefinition *particleDef = theTrack -> GetDefinition();
 98     //Get particle name                            96     //Get particle name
 99     G4String particleName =  particleDef -> Ge     97     G4String particleName =  particleDef -> GetParticleName();
100                                                    98     
101     // Get particle PDG code                       99     // Get particle PDG code
102     G4int pdg = particleDef ->GetPDGEncoding()    100     G4int pdg = particleDef ->GetPDGEncoding();
103                                                   101     
104     // Get unique track_id (in an event)          102     // Get unique track_id (in an event)
105     G4int trackID = theTrack -> GetTrackID();     103     G4int trackID = theTrack -> GetTrackID();
106                                                   104     
107     G4double energyDeposit = aStep -> GetTotal    105     G4double energyDeposit = aStep -> GetTotalEnergyDeposit();
108                                                   106     
109     G4double DX = aStep -> GetStepLength();       107     G4double DX = aStep -> GetStepLength();
110     G4int Z = particleDef-> GetAtomicNumber();    108     G4int Z = particleDef-> GetAtomicNumber();
111     G4int A = particleDef-> GetAtomicMass();      109     G4int A = particleDef-> GetAtomicMass();
112     G4StepPoint* PreStep = aStep->GetPreStepPo    110     G4StepPoint* PreStep = aStep->GetPreStepPoint();
113                                                   111     
114     // Read voxel indexes: i is the x index, k    112     // Read voxel indexes: i is the x index, k is the z index
115     const G4VTouchable* touchable = aStep->Get    113     const G4VTouchable* touchable = aStep->GetPreStepPoint()->GetTouchable();
116     G4int k  = touchable->GetReplicaNumber(0);    114     G4int k  = touchable->GetReplicaNumber(0);
117     G4int i  = touchable->GetReplicaNumber(2);    115     G4int i  = touchable->GetReplicaNumber(2);
118     G4int j  = touchable->GetReplicaNumber(1);    116     G4int j  = touchable->GetReplicaNumber(1);
119                                                   117     
120     G4TouchableHandle touchPreStep = PreStep->    118     G4TouchableHandle touchPreStep = PreStep->GetTouchableHandle();
121     G4VPhysicalVolume* volumePre = touchPreSte    119     G4VPhysicalVolume* volumePre = touchPreStep->GetVolume();
122     G4String namePre = volumePre->GetName();      120     G4String namePre = volumePre->GetName();
123                                                   121     
124                                                << 122     
                                                   >> 123     
                                                   >> 124     
125                                                   125     
126     HadrontherapyMatrix* matrix = Hadrontherap    126     HadrontherapyMatrix* matrix = HadrontherapyMatrix::GetInstance();
127     HadrontherapyLet* let = HadrontherapyLet::    127     HadrontherapyLet* let = HadrontherapyLet::GetInstance();
128                                                   128     
129     G4int* hitTrack = matrix -> GetHitTrack(i,    129     G4int* hitTrack = matrix -> GetHitTrack(i,j,k);
130                                                   130     
131                                                   131     
132     //  ******************** let *************    132     //  ******************** let ***************************
133     if (let)                                      133     if (let)
134     {                                             134     {
135         if ( !(Z==0 && A==1) ) // All but not     135         if ( !(Z==0 && A==1) ) // All but not neutrons
136         {                                         136         {
137             if( energyDeposit>0. && DX >0. )//    137             if( energyDeposit>0. && DX >0. )// calculate only energy deposit
138             {                                     138             {
139                 if (pdg !=22 && pdg !=11) // n    139                 if (pdg !=22 && pdg !=11) // not gamma and electrons
140                 {                                 140                 {
141                                                   141                     
142                     // Get the pre-step kineti    142                     // Get the pre-step kinetic energy
143                     G4double eKinPre = aStep -    143                     G4double eKinPre = aStep -> GetPreStepPoint() -> GetKineticEnergy();
144                     // Get the post-step kinet    144                     // Get the post-step kinetic energy
145                     G4double eKinPost = aStep     145                     G4double eKinPost = aStep -> GetPostStepPoint() -> GetKineticEnergy();
146                     // Get the step average ki    146                     // Get the step average kinetic energy
147                     G4double eKinMean = (eKinP    147                     G4double eKinMean = (eKinPre + eKinPost) * 0.5;
148                                                   148                     
149                     // get the material           149                     // get the material
150                     const G4Material * materia << 150                     G4Material * materialStep = aStep -> GetPreStepPoint() -> GetMaterial();
151                                                   151                     
152                     // get the secondary patic    152                     // get the secondary paticles
153                     G4Step fstep = *theTrack -    153                     G4Step fstep = *theTrack -> GetStep();
154                     // store all the secondary    154                     // store all the secondary partilce in current step
155                     const std::vector<const G4    155                     const std::vector<const G4Track*> * secondary = fstep.GetSecondaryInCurrentStep();
156                                                   156                     
157                     size_t SecondarySize = (*s    157                     size_t SecondarySize = (*secondary).size();
158                     G4double EnergySecondary =    158                     G4double EnergySecondary = 0.;
159                                                   159                     
160                     // get secondary electrons    160                     // get secondary electrons energy deposited
161                     if (SecondarySize) // calc    161                     if (SecondarySize) // calculate only secondary particles
162                     {                             162                     {
163                         for (size_t numsec = 0    163                         for (size_t numsec = 0; numsec< SecondarySize ; numsec ++)
164                         {                         164                         {
165                             //Get the PDG code    165                             //Get the PDG code of every secondaty particles in current step
166                             G4int PDGSecondary    166                             G4int PDGSecondary=(*secondary)[numsec]->GetDefinition()->GetPDGEncoding();
167                                                   167                             
168                             if(PDGSecondary ==    168                             if(PDGSecondary == 11) // calculate only secondary electrons
169                             {                     169                             {
170                                 // calculate t    170                                 // calculate the energy deposit of secondary electrons in current step
171                                 EnergySecondar    171                                 EnergySecondary += (*secondary)[numsec]->GetKineticEnergy();
172                             }                     172                             }
173                         }                         173                         }
174                                                   174                         
175                     }                             175                     }
176                                                   176                     
177                     // call the let filldatas     177                     // call the let filldatas function to calculate let
178                     let -> FillEnergySpectrum(    178                     let -> FillEnergySpectrum(trackID, particleDef, eKinMean, materialStep,
179                                                   179                                               energyDeposit,EnergySecondary,DX, i, j, k);
180                 }                                 180                 }
181             }                                     181             }
182         }                                         182         }
183     }                                             183     }
184                                                   184     
185                                                   185     
186     if (matrix)                                   186     if (matrix)
187     {                                             187     {
188                                                   188         
189         // Increment Fluences & accumulate ene    189         // Increment Fluences & accumulate energy spectra
190         // Hit voxels are marked with track_id    190         // Hit voxels are marked with track_id throught hitTrack matrix
191         //G4int* hitTrack = matrix -> GetHitTr    191         //G4int* hitTrack = matrix -> GetHitTrack(i,j,k); // hitTrack MUST BE cleared at every eventAction!
192         if ( *hitTrack != trackID )               192         if ( *hitTrack != trackID )
193         {                                         193         {
194             *hitTrack = trackID;                  194             *hitTrack = trackID;
195             /*                                    195             /*
196              * Fill FLUENCE data for every sin    196              * Fill FLUENCE data for every single nuclide
197              * Exclude e-, neutrons, gamma, ..    197              * Exclude e-, neutrons, gamma, ...
198              */                                   198              */
199             if ( Z >= 1) matrix -> Fill(trackI    199             if ( Z >= 1) matrix -> Fill(trackID, particleDef, i, j, k, 0, true);
200                                                   200             
201         }                                         201         }
202                                                   202         
203         if(energyDeposit != 0)                    203         if(energyDeposit != 0)
204         {                                         204         {
205             /*                                    205             /*
206              *  This method will fill a dose m    206              *  This method will fill a dose matrix for every single nuclide.
207              *  A method of the HadrontherapyM    207              *  A method of the HadrontherapyMatrix class (StoreDoseFluenceAscii())
208              *  is called automatically at the    208              *  is called automatically at the end of main (or via the macro command /analysis/writeDoseFile.
209              *  It permits to store all dose/f    209              *  It permits to store all dose/fluence data into a single plane ASCII file.
210              */                                   210              */
211             // if (A==1 && Z==1) // primary an    211             // if (A==1 && Z==1) // primary and sec. protons
212             if ( Z>=1 )    //  exclude e-, neu    212             if ( Z>=1 )    //  exclude e-, neutrons, gamma, ...
213                 matrix -> Fill(trackID, partic    213                 matrix -> Fill(trackID, particleDef, i, j, k, energyDeposit);
214             /*                                    214             /*
215              * Create a hit with the informati    215              * Create a hit with the information of position is in the detector
216              */                                   216              */
217             HadrontherapyDetectorHit* detector    217             HadrontherapyDetectorHit* detectorHit = new HadrontherapyDetectorHit();
218             detectorHit -> SetEdepAndPosition(    218             detectorHit -> SetEdepAndPosition(i, j, k, energyDeposit);
219             HitsCollection -> insert(detectorH    219             HitsCollection -> insert(detectorHit);
220         }                                         220         }
221     }                                             221     }
222                                                << 
223     auto rbe = HadrontherapyRBE::GetInstance() << 
224     if (rbe->IsCalculationEnabled())           << 
225     {                                          << 
226         if (!fRBEAccumulable)                  << 
227         {                                      << 
228             fRBEAccumulable = dynamic_cast<Had << 
229             if (!fRBEAccumulable)              << 
230             {                                  << 
231                 G4Exception("HadrontherapyDete << 
232             }                                  << 
233         }                                      << 
234   if (A>0) //protect against gammas, e- , etc  << 
235     {                                          << 
236       fRBEAccumulable->Accumulate(kineticEnerg << 
237     }                                          << 
238     }                                          << 
239                                                << 
240                                                << 
241     return true;                                  222     return true;
242 }                                                 223 }
243                                                   224 
244 //////////////////////////////////////////////    225 /////////////////////////////////////////////////////////////////////////////
245 void HadrontherapyDetectorSD::EndOfEvent(G4HCo    226 void HadrontherapyDetectorSD::EndOfEvent(G4HCofThisEvent* HCE)
246 {                                                 227 {
247                                                   228     
248     static G4int HCID = -1;                       229     static G4int HCID = -1;
249     if(HCID < 0)                                  230     if(HCID < 0)
250     {                                             231     {
251         HCID = GetCollectionID(0);                232         HCID = GetCollectionID(0);
252     }                                             233     }
253                                                   234     
254     HCE -> AddHitsCollection(HCID,HitsCollecti    235     HCE -> AddHitsCollection(HCID,HitsCollection);
255 }                                                 236 }
256                                                   237 
257                                                   238