Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/iort_therapy/src/IORTSteppingAction.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/iort_therapy/src/IORTSteppingAction.cc (Version 11.3.0) and /examples/advanced/iort_therapy/src/IORTSteppingAction.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 // This is the *BASIC* version of IORT, a Gean     26 // This is the *BASIC* version of IORT, a Geant4-based application
 27 //                                                 27 //
 28 // Main Authors: G.Russo(a,b), C.Casarino*(c),     28 // Main Authors: G.Russo(a,b), C.Casarino*(c), G.C. Candiano(c), G.A.P. Cirrone(d), F.Romano(d)
 29 // Contributor Authors: S.Guatelli(e)              29 // Contributor Authors: S.Guatelli(e)
 30 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d     30 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d)
 31 //                                                 31 //    
 32 //   (a) Fondazione Istituto San Raffaele G.Gi     32 //   (a) Fondazione Istituto San Raffaele G.Giglio, Cefalù, Italy
 33 //   (b) IBFM-CNR , Segrate (Milano), Italy        33 //   (b) IBFM-CNR , Segrate (Milano), Italy
 34 //   (c) LATO (Laboratorio di Tecnologie Oncol     34 //   (c) LATO (Laboratorio di Tecnologie Oncologiche), Cefalù, Italy
 35 //   (d) Laboratori Nazionali del Sud of the I     35 //   (d) Laboratori Nazionali del Sud of the INFN, Catania, Italy
 36 //   (e) University of Wollongong, Australia   <<  36 //   (e) University of Wallongong, Australia
 37 //                                                 37 //
 38 //   *Corresponding author, email to carlo.cas     38 //   *Corresponding author, email to carlo.casarino@polooncologicocefalu.it
 39 //////////////////////////////////////////////     39 //////////////////////////////////////////////////////////////////////////////////////////////
 40                                                    40 
 41 #include "G4SystemOfUnits.hh"                  <<  41 #include <CLHEP/Units/SystemOfUnits.h>
                                                   >>  42 
 42 #include "G4SteppingManager.hh"                    43 #include "G4SteppingManager.hh"
 43 #include "G4TrackVector.hh"                        44 #include "G4TrackVector.hh"
 44 #include "IORTSteppingAction.hh"                   45 #include "IORTSteppingAction.hh"
 45 #include "G4ios.hh"                                46 #include "G4ios.hh"
 46 #include "G4SteppingManager.hh"                    47 #include "G4SteppingManager.hh"
 47 #include "G4Track.hh"                              48 #include "G4Track.hh"
 48 #include "G4Step.hh"                               49 #include "G4Step.hh"
 49 #include "G4StepPoint.hh"                          50 #include "G4StepPoint.hh"
 50 #include "G4TrackStatus.hh"                        51 #include "G4TrackStatus.hh"
 51 #include "G4TrackVector.hh"                        52 #include "G4TrackVector.hh"
 52 #include "G4ParticleDefinition.hh"                 53 #include "G4ParticleDefinition.hh"
 53 #include "G4ParticleTypes.hh"                      54 #include "G4ParticleTypes.hh"
 54 #include "G4UserEventAction.hh"                <<  55 #include "G4UserEventAction.hh"  // NOT INCLUDED IN geant4.9.3p01_version
                                                   >>  56 
                                                   >>  57 #include "IORTAnalysisManager.hh"
                                                   >>  58 
 55 #include "IORTRunAction.hh"                        59 #include "IORTRunAction.hh"
 56 #include "G4SystemOfUnits.hh"                      60 #include "G4SystemOfUnits.hh"
 57                                                    61 
 58 IORTSteppingAction::IORTSteppingAction( IORTRu <<  62 /////////////////////////////////////////////////////////////////////////////
 59 {}                                             <<  63 IORTSteppingAction::IORTSteppingAction( IORTRunAction *run)
                                                   >>  64 {
                                                   >>  65     runAction = run;
                                                   >>  66 }
 60                                                    67 
                                                   >>  68 /////////////////////////////////////////////////////////////////////////////
 61 IORTSteppingAction::~IORTSteppingAction()          69 IORTSteppingAction::~IORTSteppingAction()
 62 {}                                             <<  70 {
                                                   >>  71 }
 63                                                    72 
 64 void IORTSteppingAction::UserSteppingAction(co <<  73 /////////////////////////////////////////////////////////////////////////////
                                                   >>  74 void IORTSteppingAction::UserSteppingAction(const G4Step* aStep)
 65 {                                                  75 { 
 66                                                <<  76     /*
                                                   >>  77     // USEFULL METHODS TO RETRIEVE INFORMATION DURING THE STEPS
                                                   >>  78     if( (aStep->GetTrack()->GetVolume()->GetName() == "DetectorPhys") 
                                                   >>  79     && aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton")
                                                   >>  80     //G4int evtNb = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID();
                                                   >>  81     {
                                                   >>  82     G4cout << "ENERGIA:   " << aStep->GetTrack()->GetKineticEnergy() 
                                                   >>  83     << "   VOLUME  " << aStep->GetTrack()->GetVolume()->GetName()
                                                   >>  84     << "   MATERIALE    " <<  aStep -> GetTrack() -> GetMaterial() -> GetName()
                                                   >>  85     << "   EVENTO     " << G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID()
                                                   >>  86     << "   POS     " << aStep->GetTrack()->GetPosition().x()
                                                   >>  87     << G4endl;
                                                   >>  88     }
                                                   >>  89     */
                                                   >>  90 
                                                   >>  91     if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){
                                                   >>  92   G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
                                                   >>  93   G4double secondaryParticleKineticEnergy =  aStep->GetTrack()->
                                                   >>  94     GetKineticEnergy();     
                                                   >>  95   G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
                                                   >>  96   G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
                                                   >>  97   if(particleType == "nucleus") {
                                                   >>  98       G4int A = def->GetBaryonNumber();
                                                   >>  99       G4double Z = def->GetPDGCharge();
                                                   >> 100       G4double posX = aStep->GetTrack()->GetPosition().x() /cm;
                                                   >> 101       G4double posY = aStep->GetTrack()->GetPosition().y() /cm;
                                                   >> 102       G4double posZ = aStep->GetTrack()->GetPosition().z() /cm;
                                                   >> 103       G4double energy = secondaryParticleKineticEnergy / A /MeV;
                                                   >> 104 
                                                   >> 105       IORTAnalysisManager* analysisMgr =  IORTAnalysisManager::GetInstance();   
                                                   >> 106       analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ);
                                                   >> 107   } else if(particleName == "proton") {   // proton (hydrogen-1) is a special case
                                                   >> 108       G4double posX = aStep->GetTrack()->GetPosition().x() /cm ;
                                                   >> 109       G4double posY = aStep->GetTrack()->GetPosition().y() /cm ;
                                                   >> 110       G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
                                                   >> 111       G4double energy = secondaryParticleKineticEnergy / MeV;    // Hydrogen-1: A = 1, Z = 1
                                                   >> 112       IORTAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
                                                   >> 113   }
                                                   >> 114 
                                                   >> 115   G4String secondaryParticleName =  def -> GetParticleName();  
                                                   >> 116   //G4cout <<"Particle: " << secondaryParticleName << G4endl;
                                                   >> 117   //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
                                                   >> 118   IORTAnalysisManager* analysis =  IORTAnalysisManager::GetInstance();   
                                                   >> 119   //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this.
                                                   >> 120   if(secondaryParticleName == "proton") {
                                                   >> 121       analysis->hydrogenEnergy(secondaryParticleKineticEnergy /MeV);
                                                   >> 122   }
                                                   >> 123   if(secondaryParticleName == "deuteron") {
                                                   >> 124       analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
                                                   >> 125   }
                                                   >> 126   if(secondaryParticleName == "triton") {
                                                   >> 127       analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
                                                   >> 128   }
                                                   >> 129   if(secondaryParticleName == "alpha") {
                                                   >> 130       analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
                                                   >> 131   }
                                                   >> 132   if(secondaryParticleName == "He3"){
                                                   >> 133       analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);   
                                                   >> 134   }
                                                   >> 135 
                                                   >> 136   aStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
                                                   >> 137     }
                                                   >> 138 
                                                   >> 139     // Electromagnetic and hadronic processes of primary particles in the phantom
                                                   >> 140     //setting phantomPhys correctly will break something here fixme
                                                   >> 141     if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
                                                   >> 142       (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
                                                   >> 143       (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
                                                   >> 144     {
                                                   >> 145   G4String process = aStep -> GetPostStepPoint() -> 
                                                   >> 146       GetProcessDefinedStep() -> GetProcessName();
                                                   >> 147 
                                                   >> 148   if ((process == "Transportation") || (process == "StepLimiter")) {;}
                                                   >> 149   else 
                                                   >> 150   {
                                                   >> 151       if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni")) 
                                                   >> 152       { 
                                                   >> 153     runAction -> AddEMProcess();
                                                   >> 154       } 
                                                   >> 155       else  
                                                   >> 156       {
                                                   >> 157     runAction -> AddHadronicProcess();
                                                   >> 158 
                                                   >> 159     if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
                                                   >> 160         G4cout << "Warning! Unknown proton process: "<< process << G4endl;
                                                   >> 161       }
                                                   >> 162   }         
                                                   >> 163     }
                                                   >> 164 
                                                   >> 165     // Retrieve information about the secondary particles originated in the phantom
                                                   >> 166 
                                                   >> 167     G4SteppingManager*  steppingManager = fpSteppingManager;
                                                   >> 168 
                                                   >> 169     // check if it is alive
                                                   >> 170     //if(theTrack-> GetTrackStatus() == fAlive) { return; }
                                                   >> 171 
                                                   >> 172     // Retrieve the secondary particles
                                                   >> 173     G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
                                                   >> 174 
                                                   >> 175     for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
                                                   >> 176     { 
                                                   >> 177   G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); 
                                                   >> 178 
                                                   >> 179   if (volumeName == "phantomPhys")
                                                   >> 180   {
                                                   >> 181     G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();  
                                                   >> 182       G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();     
                                                   >> 183 
                                                   >> 184       IORTAnalysisManager* analysis =  IORTAnalysisManager::GetInstance();   
                                                   >> 185 
                                                   >> 186       if (secondaryParticleName == "e-")
                                                   >> 187     analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
                                                   >> 188 
                                                   >> 189       if (secondaryParticleName == "gamma")
                                                   >> 190     analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
                                                   >> 191 
                                                   >> 192       if (secondaryParticleName == "deuteron")
                                                   >> 193     analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
                                                   >> 194 
                                                   >> 195       if (secondaryParticleName == "triton")
                                                   >> 196     analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
                                                   >> 197 
                                                   >> 198       if (secondaryParticleName == "alpha")
                                                   >> 199     analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
                                                   >> 200 
                                                   >> 201       G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
                                                   >> 202       if (z > 0.)
                                                   >> 203       {      
                                                   >> 204     G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
                                                   >> 205     G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy(); 
                                                   >> 206 
                                                   >> 207     // If a generic ion is originated in the detector, its baryonic number, PDG charge, 
                                                   >> 208     // total number of electrons in the orbitals are stored in a ntuple 
                                                   >> 209     analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
                                                   >> 210       }
                                                   >> 211   }
                                                   >> 212     }
 67 }                                                 213 }
 68                                                   214