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.3.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 "HadrontherapyAnalysisManager.hh"
 38 #include "G4SystemOfUnits.hh"                      39 #include "G4SystemOfUnits.hh"
 39 #include "HadrontherapyMatrix.hh"              << 
 40                                                << 
 41                                                << 
 42 #include "G4SteppingManager.hh"                << 
 43 #include "G4TrackVector.hh"                    << 
 44 #include "HadrontherapySteppingAction.hh"      << 
 45 #include "G4ios.hh"                            << 
 46 #include "G4SteppingManager.hh"                << 
 47 #include "G4Track.hh"                          << 
 48 #include "G4Step.hh"                           << 
 49 #include "G4StepPoint.hh"                      << 
 50 #include "G4TrackStatus.hh"                    << 
 51 #include "G4TrackVector.hh"                    << 
 52 #include "G4ParticleDefinition.hh"             << 
 53 #include "G4ParticleTypes.hh"                  << 
 54 #include "G4UserEventAction.hh"                << 
 55 #include "G4TransportationManager.hh"          << 
 56 #include "G4VSensitiveDetector.hh"             << 
 57 #include "HadrontherapyRunAction.hh"           << 
 58 #include "G4SystemOfUnits.hh"                  << 
 59 #include "HadrontherapyRBE.hh"                 << 
 60 #include <G4AccumulableManager.hh>             << 
 61                                                << 
 62                                                    40 
 63 //////////////////////////////////////////////     41 /////////////////////////////////////////////////////////////////////////////
 64 HadrontherapyDetectorSD::HadrontherapyDetector     42 HadrontherapyDetectorSD::HadrontherapyDetectorSD(G4String name):
 65 G4VSensitiveDetector(name)                     <<  43     G4VSensitiveDetector(name)
 66 {                                              <<  44 { 
 67     G4String HCname;                               45     G4String HCname;
 68     collectionName.insert(HCname="Hadrontherap     46     collectionName.insert(HCname="HadrontherapyDetectorHitsCollection");
 69     HitsCollection = NULL;                     <<  47     HitsCollection = NULL; 
 70     sensitiveDetectorName = name;                  48     sensitiveDetectorName = name;
 71                                                <<  49 
 72 }                                                  50 }
 73                                                    51 
 74 //////////////////////////////////////////////     52 /////////////////////////////////////////////////////////////////////////////
 75 HadrontherapyDetectorSD::~HadrontherapyDetecto     53 HadrontherapyDetectorSD::~HadrontherapyDetectorSD()
 76 {}                                             <<  54 { 
                                                   >>  55 }
 77                                                    56 
 78 //////////////////////////////////////////////     57 /////////////////////////////////////////////////////////////////////////////
 79 void HadrontherapyDetectorSD::Initialize(G4HCo     58 void HadrontherapyDetectorSD::Initialize(G4HCofThisEvent*)
 80 {                                                  59 {
 81                                                <<  60   
 82     HitsCollection = new HadrontherapyDetector     61     HitsCollection = new HadrontherapyDetectorHitsCollection(sensitiveDetectorName,
 83                                                <<  62       collectionName[0]);
 84 }                                                  63 }
 85                                                    64 
 86 //////////////////////////////////////////////     65 /////////////////////////////////////////////////////////////////////////////
 87 G4bool HadrontherapyDetectorSD::ProcessHits(G4     66 G4bool HadrontherapyDetectorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* )
 88 {                                                  67 {
 89                                                <<  68 
 90                                                <<  69  
 91     if (aStep -> GetPreStepPoint() -> GetPhysi <<  70   if (aStep -> GetPreStepPoint() -> GetPhysicalVolume() -> GetName() != "RODetectorZDivisionPhys") return false;
 92                                                <<  71   
 93                                                <<  72 
 94     // Get kinetic energy                          73     // Get kinetic energy
 95     G4Track * theTrack = aStep  ->  GetTrack()     74     G4Track * theTrack = aStep  ->  GetTrack();
 96     G4double kineticEnergy = theTrack->GetKine <<  75     G4double kineticEnergy =  theTrack -> GetKineticEnergy();  
                                                   >>  76 
 97     G4ParticleDefinition *particleDef = theTra     77     G4ParticleDefinition *particleDef = theTrack -> GetDefinition();
 98     //Get particle name                        <<  78     //Get particle name  
 99     G4String particleName =  particleDef -> Ge <<  79     G4String particleName =  particleDef -> GetParticleName();  
100                                                <<  80 
101     // Get particle PDG code                   << 
102     G4int pdg = particleDef ->GetPDGEncoding()     81     G4int pdg = particleDef ->GetPDGEncoding();
103                                                <<  82 
104     // Get unique track_id (in an event)           83     // Get unique track_id (in an event)
105     G4int trackID = theTrack -> GetTrackID();      84     G4int trackID = theTrack -> GetTrackID();
106                                                <<  85 
107     G4double energyDeposit = aStep -> GetTotal     86     G4double energyDeposit = aStep -> GetTotalEnergyDeposit();
108                                                <<  87 
109     G4double DX = aStep -> GetStepLength();        88     G4double DX = aStep -> GetStepLength();
110     G4int Z = particleDef-> GetAtomicNumber();     89     G4int Z = particleDef-> GetAtomicNumber();
111     G4int A = particleDef-> GetAtomicMass();       90     G4int A = particleDef-> GetAtomicMass();
112     G4StepPoint* PreStep = aStep->GetPreStepPo <<  91 
113                                                <<  92   
114     // Read voxel indexes: i is the x index, k <<  93  // Read voxel indexes: i is the x index, k is the z index
115     const G4VTouchable* touchable = aStep->Get <<  94  const G4VTouchable* touchable = aStep->GetPreStepPoint()->GetTouchable();
116     G4int k  = touchable->GetReplicaNumber(0); <<  95  G4int k  = touchable->GetReplicaNumber(0);
117     G4int i  = touchable->GetReplicaNumber(2); <<  96  G4int i  = touchable->GetReplicaNumber(2);
118     G4int j  = touchable->GetReplicaNumber(1); <<  97  G4int j  = touchable->GetReplicaNumber(1);
119                                                << 
120     G4TouchableHandle touchPreStep = PreStep-> << 
121     G4VPhysicalVolume* volumePre = touchPreSte << 
122     G4String namePre = volumePre->GetName();   << 
123                                                << 
124                                                    98 
125                                                    99     
                                                   >> 100 
                                                   >> 101 #ifdef G4ANALYSIS_USE_ROOT
                                                   >> 102     HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
                                                   >> 103 #endif
                                                   >> 104 
126     HadrontherapyMatrix* matrix = Hadrontherap    105     HadrontherapyMatrix* matrix = HadrontherapyMatrix::GetInstance();
127     HadrontherapyLet* let = HadrontherapyLet::    106     HadrontherapyLet* let = HadrontherapyLet::GetInstance();
128                                                   107     
129     G4int* hitTrack = matrix -> GetHitTrack(i, << 
130                                                << 
131                                                   108     
132     //  ******************** let ************* << 109   //  ******************** let ***************************
133     if (let)                                   << 110      if (let)
134     {                                             111     {
135         if ( !(Z==0 && A==1) ) // All but not  << 112      if ( !(Z==0 && A==1) ) // All but not neutrons 
136         {                                      << 113     {
137             if( energyDeposit>0. && DX >0. )// << 114     if( energyDeposit>0. && DX >0. )
138             {                                  << 115       {
139                 if (pdg !=22 && pdg !=11) // n << 116         if (pdg !=22) // not gamma
140                 {                              << 117     {
141                                                << 118       let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i, j, k);
142                     // Get the pre-step kineti << 119     }
143                     G4double eKinPre = aStep - << 120         else if (kineticEnergy > 50.*keV) // gamma cut
144                     // Get the post-step kinet << 121     {
145                     G4double eKinPost = aStep  << 122         let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i , j, k);
146                     // Get the step average ki << 123     }
147                     G4double eKinMean = (eKinP << 124 
148                                                << 125       }
149                     // get the material        << 126   }
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                 }                              << 
181             }                                  << 
182         }                                      << 
183     }                                             127     }
184                                                   128     
185                                                   129     
                                                   >> 130     
                                                   >> 131     //  ******************** let ***************************
                                                   >> 132   
                                                   >> 133     
186     if (matrix)                                   134     if (matrix)
187     {                                             135     {
188                                                << 
189         // Increment Fluences & accumulate ene << 
190         // Hit voxels are marked with track_id << 
191         //G4int* hitTrack = matrix -> GetHitTr << 
192         if ( *hitTrack != trackID )            << 
193         {                                      << 
194             *hitTrack = trackID;               << 
195             /*                                 << 
196              * Fill FLUENCE data for every sin << 
197              * Exclude e-, neutrons, gamma, .. << 
198              */                                << 
199             if ( Z >= 1) matrix -> Fill(trackI << 
200                                                << 
201         }                                      << 
202                                                << 
203         if(energyDeposit != 0)                 << 
204         {                                      << 
205             /*                                 << 
206              *  This method will fill a dose m << 
207              *  A method of the HadrontherapyM << 
208              *  is called automatically at the << 
209              *  It permits to store all dose/f << 
210              */                                << 
211             // if (A==1 && Z==1) // primary an << 
212             if ( Z>=1 )    //  exclude e-, neu << 
213                 matrix -> Fill(trackID, partic << 
214             /*                                 << 
215              * Create a hit with the informati << 
216              */                                << 
217             HadrontherapyDetectorHit* detector << 
218             detectorHit -> SetEdepAndPosition( << 
219             HitsCollection -> insert(detectorH << 
220         }                                      << 
221     }                                          << 
222                                                   136 
223     auto rbe = HadrontherapyRBE::GetInstance() << 137   // Increment Fluences & accumulate energy spectra
224     if (rbe->IsCalculationEnabled())           << 138   // Hit voxels are marked with track_id throught hitTrack matrix
225     {                                          << 139   G4int* hitTrack = matrix -> GetHitTrack(i,j,k); // hitTrack MUST BE cleared at every eventAction!
226         if (!fRBEAccumulable)                  << 140   if ( *hitTrack != trackID )
227         {                                      << 141   {
228             fRBEAccumulable = dynamic_cast<Had << 142                  *hitTrack = trackID;
229             if (!fRBEAccumulable)              << 143     /*
230             {                                  << 144      * Fill FLUENCE data for every single nuclide 
231                 G4Exception("HadrontherapyDete << 145      * Exclude e-, neutrons, gamma, ...
232             }                                  << 146      */
233         }                                      << 147     if ( Z >= 1) matrix -> Fill(trackID, particleDef, i, j, k, 0, true);
234   if (A>0) //protect against gammas, e- , etc  << 148         
235     {                                          << 149 
236       fRBEAccumulable->Accumulate(kineticEnerg << 150 #ifdef G4ANALYSIS_USE_ROOT
237     }                                          << 151 /*
                                                   >> 152       // Fragments kinetic energy (ntuple)
                                                   >> 153       if (trackID !=1 && Z>=1) 
                                                   >> 154       {
                                                   >> 155     // First step kinetic energy for every fragment 
                                                   >> 156      analysis -> FillKineticFragmentTuple(i, j, k, A, Z, kineticEnergy/MeV);
                                                   >> 157       }  
                                                   >> 158       // Kinetic energy spectra for primary particles 
                                                   >> 159 
                                                   >> 160       if ( trackID == 1 && i == 0) 
                                                   >> 161       {
                                                   >> 162     // First step kinetic energy for primaries only
                                                   >> 163     analysis -> FillKineticEnergyPrimaryNTuple(i, j, k, kineticEnergy/MeV);
                                                   >> 164       }
                                                   >> 165 */
                                                   >> 166 #endif
                                                   >> 167   }  
                                                   >> 168 
                                                   >> 169   if(energyDeposit != 0)
                                                   >> 170   {
                                                   >> 171 /*
                                                   >> 172  *  This method will fill a dose matrix for every single nuclide. 
                                                   >> 173  *  A method of the HadrontherapyMatrix class (StoreDoseFluenceAscii())
                                                   >> 174  *  is called automatically at the end of main (or via the macro command /analysis/writeDoseFile.
                                                   >> 175  *  It permits to store all dose/fluence data into a single plane ASCII file. 
                                                   >> 176 */      
                                                   >> 177       // if (A==1 && Z==1) // primary and sec. protons 
                                                   >> 178       if ( Z>=1 )    //  exclude e-, neutrons, gamma, ...
                                                   >> 179         matrix -> Fill(trackID, particleDef, i, j, k, energyDeposit);
                                                   >> 180       /*
                                                   >> 181        * Create a hit with the information of position is in the detector     
                                                   >> 182        */
                                                   >> 183       HadrontherapyDetectorHit* detectorHit = new HadrontherapyDetectorHit();       
                                                   >> 184       detectorHit -> SetEdepAndPosition(i, j, k, energyDeposit); 
                                                   >> 185       HitsCollection -> insert(detectorHit);
                                                   >> 186   }
238     }                                             187     }
239                                                   188 
                                                   >> 189 #ifdef G4ANALYSIS_USE_ROOT
                                                   >> 190     if(energyDeposit != 0)
                                                   >> 191     {  
                                                   >> 192   if(trackID != 1)
                                                   >> 193   {
                                                   >> 194       if (particleName == "proton")
                                                   >> 195     analysis -> SecondaryProtonEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 196 
                                                   >> 197       else if (particleName == "neutron")
                                                   >> 198     analysis -> SecondaryNeutronEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 199 
                                                   >> 200       else if (particleName == "alpha")
                                                   >> 201     analysis -> SecondaryAlphaEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 202 
                                                   >> 203       else if (particleName == "gamma")
                                                   >> 204     analysis -> SecondaryGammaEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 205 
                                                   >> 206       else if (particleName == "e-")
                                                   >> 207     analysis -> SecondaryElectronEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 208 
                                                   >> 209       else if (particleName == "triton")
                                                   >> 210     analysis -> SecondaryTritonEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 211 
                                                   >> 212       else if (particleName == "deuteron")
                                                   >> 213     analysis -> SecondaryDeuteronEnergyDeposit(i, energyDeposit/MeV);
                                                   >> 214 
                                                   >> 215       else if (particleName == "pi+" || particleName == "pi-" ||  particleName == "pi0")
                                                   >> 216     analysis -> SecondaryPionEnergyDeposit(i, energyDeposit/MeV);     
                                                   >> 217   }
                                                   >> 218     }
                                                   >> 219 #endif
240                                                   220 
241     return true;                                  221     return true;
242 }                                                 222 }
243                                                   223 
244 //////////////////////////////////////////////    224 /////////////////////////////////////////////////////////////////////////////
245 void HadrontherapyDetectorSD::EndOfEvent(G4HCo    225 void HadrontherapyDetectorSD::EndOfEvent(G4HCofThisEvent* HCE)
246 {                                                 226 {
247                                                << 227  
248     static G4int HCID = -1;                       228     static G4int HCID = -1;
249     if(HCID < 0)                                  229     if(HCID < 0)
250     {                                          << 230     { 
251         HCID = GetCollectionID(0);             << 231   HCID = GetCollectionID(0); 
252     }                                             232     }
253                                                << 233 
254     HCE -> AddHitsCollection(HCID,HitsCollecti    234     HCE -> AddHitsCollection(HCID,HitsCollection);
255 }                                                 235 }
256                                                   236 
257                                                   237