- Timestamp:
- Jun 14, 2010, 3:54:58 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/examples/advanced/hadrontherapy/src/HadrontherapySteppingAction.cc
r1230 r1313 39 39 #include "G4ParticleDefinition.hh" 40 40 #include "G4ParticleTypes.hh" 41 #include "G4UserEventAction.hh" 41 42 42 43 #include "HadrontherapyAnalysisManager.hh" … … 47 48 HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run) 48 49 { 49 runAction = run;50 runAction = run; 50 51 } 51 52 … … 58 59 void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep) 59 60 { 61 /* 62 // USEFULL METHODS TO RETRIEVE INFORMATION DURING THE STEPS 63 if( (aStep->GetTrack()->GetVolume()->GetName() == "DetectorPhys") 64 && aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton") 65 //G4int evtNb = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID(); 66 { 67 G4cout << "ENERGIA: " << aStep->GetTrack()->GetKineticEnergy() 68 << " VOLUME " << aStep->GetTrack()->GetVolume()->GetName() 69 << " MATERIALE " << aStep -> GetTrack() -> GetMaterial() -> GetName() 70 << " EVENTO " << G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID() 71 << " POS " << aStep->GetTrack()->GetPosition().x() 72 << G4endl; 73 } 74 */ 75 60 76 if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){ 61 #ifdef ANALYSIS_USE62 63 64 65 66 67 G4int A = def->GetBaryonNumber();68 G4double Z = def->GetPDGCharge();69 G4double posX = aStep->GetTrack()->GetPosition().x() / cm;70 G4double posY = aStep->GetTrack()->GetPosition().y() / cm;71 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm;72 G4double energy = secondaryParticleKineticEnergy / A / MeV;73 74 HadrontherapyAnalysisManager* analysisMgr = HadrontherapyAnalysisManager::getInstance();75 analysisMgr->fillFragmentTuple(A, Z, energy, posX, posY, posZ);76 77 G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;78 G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;79 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;80 G4double energy = secondaryParticleKineticEnergy * MeV; // Hydrogen-1: A = 1, Z = 181 HadrontherapyAnalysisManager::getInstance()->fillFragmentTuple(1, 1.0, energy, posX, posY, posZ);82 83 84 85 86 87 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager:: getInstance();77 #ifdef G4ANALYSIS_USE_ROOT 78 G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition(); 79 G4double secondaryParticleKineticEnergy = aStep->GetTrack()->GetKineticEnergy(); 80 G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei 81 G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus" 82 if(particleType == "nucleus") { 83 G4int A = def->GetBaryonNumber(); 84 G4double Z = def->GetPDGCharge(); 85 G4double posX = aStep->GetTrack()->GetPosition().x() / cm; 86 G4double posY = aStep->GetTrack()->GetPosition().y() / cm; 87 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm; 88 G4double energy = secondaryParticleKineticEnergy / A / MeV; 89 90 HadrontherapyAnalysisManager* analysisMgr = HadrontherapyAnalysisManager::GetInstance(); 91 analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ); 92 } else if(particleName == "proton") { // proton (hydrogen-1) is a special case 93 G4double posX = aStep->GetTrack()->GetPosition().x() / cm ; 94 G4double posY = aStep->GetTrack()->GetPosition().y() / cm ; 95 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ; 96 G4double energy = secondaryParticleKineticEnergy * MeV; // Hydrogen-1: A = 1, Z = 1 97 HadrontherapyAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ); 98 } 99 100 G4String secondaryParticleName = def -> GetParticleName(); 101 //G4cout <<"Particle: " << secondaryParticleName << G4endl; 102 //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl; 103 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); 88 104 //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this. 89 105 if(secondaryParticleName == "proton") { 90 analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);106 analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV); 91 107 } 92 108 if(secondaryParticleName == "deuteron") { 93 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);109 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV); 94 110 } 95 111 if(secondaryParticleName == "triton") { 96 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);112 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV); 97 113 } 98 114 if(secondaryParticleName == "alpha") { 99 analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);115 analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV); 100 116 } 101 117 if(secondaryParticleName == "He3"){ 102 analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);118 analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV); 103 119 } 104 120 #endif … … 107 123 } 108 124 109 // Electromagnetic and hadronic processes of primary particles in the phantom 110 //setting phantomPhys correctly will break something here fixme 111 if ((aStep -> GetTrack() -> GetTrackID() == 1) && 112 (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") && 113 (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL)) 125 // Electromagnetic and hadronic processes of primary particles in the phantom 126 //setting phantomPhys correctly will break something here fixme 127 if ((aStep -> GetTrack() -> GetTrackID() == 1) && 128 (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") && 129 (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL)) 130 { 131 G4String process = aStep -> GetPostStepPoint() -> 132 GetProcessDefinedStep() -> GetProcessName(); 133 134 if ((process == "Transportation") || (process == "StepLimiter")) {;} 135 else 136 { 137 if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni")) 138 { 139 runAction -> AddEMProcess(); 140 } 141 else 114 142 { 115 G4String process = aStep -> GetPostStepPoint() -> 116 GetProcessDefinedStep() -> GetProcessName(); 117 118 if ((process == "Transportation") || (process == "StepLimiter")) {;} 119 else { 120 if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni")) 121 { 122 runAction -> AddEMProcess(); 123 } 124 else 125 { 126 runAction -> AddHadronicProcess(); 127 128 if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") ) 129 G4cout << "Warning! Unknown proton process: "<< process << G4endl; 130 } 131 } 143 runAction -> AddHadronicProcess(); 144 145 if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") ) 146 G4cout << "Warning! Unknown proton process: "<< process << G4endl; 132 147 } 133 134 // Retrieve information about the secondaries originated in the phantom 135 136 G4SteppingManager* steppingManager = fpSteppingManager; 137 138 // check if it is alive 139 //if(theTrack-> GetTrackStatus() == fAlive) { return; } 140 141 // Retrieve the secondary particles 142 G4TrackVector* fSecondary = steppingManager -> GetfSecondary(); 143 144 for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++) 148 } 149 } 150 151 // Retrieve information about the secondaries originated in the phantom 152 153 G4SteppingManager* steppingManager = fpSteppingManager; 154 155 // check if it is alive 156 //if(theTrack-> GetTrackStatus() == fAlive) { return; } 157 158 // Retrieve the secondary particles 159 G4TrackVector* fSecondary = steppingManager -> GetfSecondary(); 160 161 for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++) 145 162 { 146 147 148 163 G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); 164 165 if (volumeName == "phantomPhys") 149 166 { 150 #ifdef ANALYSIS_USE151 G4String secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();152 G4double secondaryParticleKineticEnergy = (*fSecondary)[lp1] -> GetKineticEnergy();153 154 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::getInstance();155 156 157 158 159 if (secondaryParticleName == "gamma")160 161 162 if (secondaryParticleName == "deuteron")163 164 165 if (secondaryParticleName == "triton")166 167 168 if (secondaryParticleName == "alpha")169 170 171 G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();172 if (z > 0.)173 174 175 176 177 178 179 167 #ifdef G4ANALYSIS_USE_ROOT 168 G4String secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); 169 G4double secondaryParticleKineticEnergy = (*fSecondary)[lp1] -> GetKineticEnergy(); 170 171 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); 172 173 if (secondaryParticleName == "e-") 174 analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV); 175 176 if (secondaryParticleName == "gamma") 177 analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV); 178 179 if (secondaryParticleName == "deuteron") 180 analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV); 181 182 if (secondaryParticleName == "triton") 183 analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV); 184 185 if (secondaryParticleName == "alpha") 186 analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV); 187 188 G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge(); 189 if (z > 0.) 190 { 191 G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber(); 192 G4int electronOccupancy = (*fSecondary)[lp1] -> GetDynamicParticle() -> GetTotalOccupancy(); 193 194 // If a generic ion is originated in the detector, its baryonic number, PDG charge, 195 // total number of electrons in the orbitals are stored in a ntuple 196 analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV); 180 197 } 181 198 #endif
Note: See TracChangeset
for help on using the changeset viewer.