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.1)


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