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.4.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 {
                                                   >>  75 }
 77                                                    76 
 78 //////////////////////////////////////////////     77 /////////////////////////////////////////////////////////////////////////////
 79 void HadrontherapyDetectorSD::Initialize(G4HCo     78 void HadrontherapyDetectorSD::Initialize(G4HCofThisEvent*)
 80 {                                                  79 {
 81                                                    80     
 82     HitsCollection = new HadrontherapyDetector     81     HitsCollection = new HadrontherapyDetectorHitsCollection(sensitiveDetectorName,
 83                                                    82                                                              collectionName[0]);
 84 }                                                  83 }
 85                                                    84 
 86 //////////////////////////////////////////////     85 /////////////////////////////////////////////////////////////////////////////
 87 G4bool HadrontherapyDetectorSD::ProcessHits(G4     86 G4bool HadrontherapyDetectorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* )
 88 {                                                  87 {
 89                                                    88     
 90                                                    89     
 91     if (aStep -> GetPreStepPoint() -> GetPhysi     90     if (aStep -> GetPreStepPoint() -> GetPhysicalVolume() -> GetName() != "RODetectorZDivisionPhys") return false;
 92                                                    91     
 93                                                    92     
 94     // Get kinetic energy                          93     // Get kinetic energy
 95     G4Track * theTrack = aStep  ->  GetTrack()     94     G4Track * theTrack = aStep  ->  GetTrack();
 96     G4double kineticEnergy = theTrack->GetKine <<  95     G4double kineticEnergy =  theTrack -> GetKineticEnergy();
                                                   >>  96     
 97     G4ParticleDefinition *particleDef = theTra     97     G4ParticleDefinition *particleDef = theTrack -> GetDefinition();
 98     //Get particle name                            98     //Get particle name
 99     G4String particleName =  particleDef -> Ge     99     G4String particleName =  particleDef -> GetParticleName();
100                                                   100     
101     // Get particle PDG code                   << 
102     G4int pdg = particleDef ->GetPDGEncoding()    101     G4int pdg = particleDef ->GetPDGEncoding();
103                                                   102     
104     // Get unique track_id (in an event)          103     // Get unique track_id (in an event)
105     G4int trackID = theTrack -> GetTrackID();     104     G4int trackID = theTrack -> GetTrackID();
106                                                   105     
107     G4double energyDeposit = aStep -> GetTotal    106     G4double energyDeposit = aStep -> GetTotalEnergyDeposit();
108                                                   107     
109     G4double DX = aStep -> GetStepLength();       108     G4double DX = aStep -> GetStepLength();
110     G4int Z = particleDef-> GetAtomicNumber();    109     G4int Z = particleDef-> GetAtomicNumber();
111     G4int A = particleDef-> GetAtomicMass();      110     G4int A = particleDef-> GetAtomicMass();
112     G4StepPoint* PreStep = aStep->GetPreStepPo    111     G4StepPoint* PreStep = aStep->GetPreStepPoint();
113                                                   112     
114     // Read voxel indexes: i is the x index, k    113     // Read voxel indexes: i is the x index, k is the z index
115     const G4VTouchable* touchable = aStep->Get    114     const G4VTouchable* touchable = aStep->GetPreStepPoint()->GetTouchable();
116     G4int k  = touchable->GetReplicaNumber(0);    115     G4int k  = touchable->GetReplicaNumber(0);
117     G4int i  = touchable->GetReplicaNumber(2);    116     G4int i  = touchable->GetReplicaNumber(2);
118     G4int j  = touchable->GetReplicaNumber(1);    117     G4int j  = touchable->GetReplicaNumber(1);
119                                                   118     
120     G4TouchableHandle touchPreStep = PreStep->    119     G4TouchableHandle touchPreStep = PreStep->GetTouchableHandle();
121     G4VPhysicalVolume* volumePre = touchPreSte    120     G4VPhysicalVolume* volumePre = touchPreStep->GetVolume();
122     G4String namePre = volumePre->GetName();      121     G4String namePre = volumePre->GetName();
                                                   >> 122     //G4double eKin = aStep -> GetPreStepPoint() -> GetKineticEnergy();
123                                                   123     
124                                                << 
125                                                   124     
126     HadrontherapyMatrix* matrix = Hadrontherap    125     HadrontherapyMatrix* matrix = HadrontherapyMatrix::GetInstance();
127     HadrontherapyLet* let = HadrontherapyLet::    126     HadrontherapyLet* let = HadrontherapyLet::GetInstance();
128                                                   127     
129     G4int* hitTrack = matrix -> GetHitTrack(i, << 128     G4int* hitTrack = matrix -> GetHitTrack(i,j,k);
130                                                   129     
131                                                   130     
132     //  ******************** let *************    131     //  ******************** let ***************************
133     if (let)                                      132     if (let)
134     {                                             133     {
135         if ( !(Z==0 && A==1) ) // All but not     134         if ( !(Z==0 && A==1) ) // All but not neutrons
136         {                                         135         {
137             if( energyDeposit>0. && DX >0. )// << 136             if( energyDeposit>0. && DX >0. )
138             {                                     137             {
139                 if (pdg !=22 && pdg !=11) // n << 138                 if (pdg !=22) // not gamma
140                 {                                 139                 {
141                                                << 140                     let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i, j, k);
142                     // Get the pre-step kineti << 
143                     G4double eKinPre = aStep - << 
144                     // Get the post-step kinet << 
145                     G4double eKinPost = aStep  << 
146                     // Get the step average ki << 
147                     G4double eKinMean = (eKinP << 
148                                                << 
149                     // get the material        << 
150                     const G4Material * materia << 
151                                                << 
152                     // get the secondary patic << 
153                     G4Step fstep = *theTrack - << 
154                     // store all the secondary << 
155                     const std::vector<const G4 << 
156                                                << 
157                     size_t SecondarySize = (*s << 
158                     G4double EnergySecondary = << 
159                                                << 
160                     // get secondary electrons << 
161                     if (SecondarySize) // calc << 
162                     {                          << 
163                         for (size_t numsec = 0 << 
164                         {                      << 
165                             //Get the PDG code << 
166                             G4int PDGSecondary << 
167                                                << 
168                             if(PDGSecondary == << 
169                             {                  << 
170                                 // calculate t << 
171                                 EnergySecondar << 
172                             }                  << 
173                         }                      << 
174                                                << 
175                     }                          << 
176                                                << 
177                     // call the let filldatas  << 
178                     let -> FillEnergySpectrum( << 
179                                                << 
180                 }                                 141                 }
                                                   >> 142                 else if (kineticEnergy > 50.*keV) // gamma cut
                                                   >> 143                 {
                                                   >> 144                     let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i , j, k);
                                                   >> 145                 }
                                                   >> 146                 
181             }                                     147             }
182         }                                         148         }
183     }                                             149     }
184                                                   150     
185                                                   151     
186     if (matrix)                                   152     if (matrix)
187     {                                             153     {
188                                                   154         
189         // Increment Fluences & accumulate ene    155         // Increment Fluences & accumulate energy spectra
190         // Hit voxels are marked with track_id    156         // Hit voxels are marked with track_id throught hitTrack matrix
191         //G4int* hitTrack = matrix -> GetHitTr    157         //G4int* hitTrack = matrix -> GetHitTrack(i,j,k); // hitTrack MUST BE cleared at every eventAction!
192         if ( *hitTrack != trackID )               158         if ( *hitTrack != trackID )
193         {                                         159         {
194             *hitTrack = trackID;                  160             *hitTrack = trackID;
195             /*                                    161             /*
196              * Fill FLUENCE data for every sin    162              * Fill FLUENCE data for every single nuclide
197              * Exclude e-, neutrons, gamma, ..    163              * Exclude e-, neutrons, gamma, ...
198              */                                   164              */
199             if ( Z >= 1) matrix -> Fill(trackI    165             if ( Z >= 1) matrix -> Fill(trackID, particleDef, i, j, k, 0, true);
200                                                   166             
201         }                                         167         }
202                                                   168         
203         if(energyDeposit != 0)                    169         if(energyDeposit != 0)
204         {                                         170         {
205             /*                                    171             /*
206              *  This method will fill a dose m    172              *  This method will fill a dose matrix for every single nuclide.
207              *  A method of the HadrontherapyM    173              *  A method of the HadrontherapyMatrix class (StoreDoseFluenceAscii())
208              *  is called automatically at the    174              *  is called automatically at the end of main (or via the macro command /analysis/writeDoseFile.
209              *  It permits to store all dose/f    175              *  It permits to store all dose/fluence data into a single plane ASCII file.
210              */                                   176              */
211             // if (A==1 && Z==1) // primary an    177             // if (A==1 && Z==1) // primary and sec. protons
212             if ( Z>=1 )    //  exclude e-, neu    178             if ( Z>=1 )    //  exclude e-, neutrons, gamma, ...
213                 matrix -> Fill(trackID, partic    179                 matrix -> Fill(trackID, particleDef, i, j, k, energyDeposit);
214             /*                                    180             /*
215              * Create a hit with the informati    181              * Create a hit with the information of position is in the detector
216              */                                   182              */
217             HadrontherapyDetectorHit* detector    183             HadrontherapyDetectorHit* detectorHit = new HadrontherapyDetectorHit();
218             detectorHit -> SetEdepAndPosition( << 184             detectorHit -> SetEdepAndPosition(i, j, k, energyDeposit); 
219             HitsCollection -> insert(detectorH    185             HitsCollection -> insert(detectorHit);
220         }                                         186         }
221     }                                             187     }
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;                                  188     return true;
242 }                                                 189 }
243                                                   190 
244 //////////////////////////////////////////////    191 /////////////////////////////////////////////////////////////////////////////
245 void HadrontherapyDetectorSD::EndOfEvent(G4HCo    192 void HadrontherapyDetectorSD::EndOfEvent(G4HCofThisEvent* HCE)
246 {                                                 193 {
247                                                   194     
248     static G4int HCID = -1;                       195     static G4int HCID = -1;
249     if(HCID < 0)                                  196     if(HCID < 0)
250     {                                          << 197     { 
251         HCID = GetCollectionID(0);             << 198         HCID = GetCollectionID(0); 
252     }                                             199     }
253                                                   200     
254     HCE -> AddHitsCollection(HCID,HitsCollecti    201     HCE -> AddHitsCollection(HCID,HitsCollection);
255 }                                                 202 }
256                                                   203 
257                                                   204