Geant4 Cross Reference |
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 "G4SteppingManager.hh" 29 #include "G4SteppingManager.hh" 30 #include "G4TrackVector.hh" 30 #include "G4TrackVector.hh" 31 #include "HadrontherapySteppingAction.hh" 31 #include "HadrontherapySteppingAction.hh" 32 #include "G4ios.hh" 32 #include "G4ios.hh" 33 #include "G4SteppingManager.hh" 33 #include "G4SteppingManager.hh" 34 #include "G4Track.hh" 34 #include "G4Track.hh" 35 #include "G4Step.hh" 35 #include "G4Step.hh" 36 #include "G4StepPoint.hh" 36 #include "G4StepPoint.hh" 37 #include "G4TrackStatus.hh" 37 #include "G4TrackStatus.hh" 38 #include "G4TrackVector.hh" 38 #include "G4TrackVector.hh" 39 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh" 40 #include "G4ParticleTypes.hh" 40 #include "G4ParticleTypes.hh" 41 #include "G4UserEventAction.hh" 41 #include "G4UserEventAction.hh" 42 #include "G4TransportationManager.hh" 42 #include "G4TransportationManager.hh" 43 #include "G4VSensitiveDetector.hh" 43 #include "G4VSensitiveDetector.hh" 44 #include "HadrontherapyRunAction.hh" 44 #include "HadrontherapyRunAction.hh" 45 #include "HadrontherapyMatrix.hh" << 45 #include "HadrontherapyAnalysisManager.hh" 46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh" 47 47 48 ////////////////////////////////////////////// 48 ///////////////////////////////////////////////////////////////////////////// 49 HadrontherapySteppingAction::HadrontherapyStep 49 HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run) 50 { 50 { 51 runAction = run; 51 runAction = run; 52 } 52 } 53 53 54 ////////////////////////////////////////////// 54 ///////////////////////////////////////////////////////////////////////////// 55 HadrontherapySteppingAction::~HadrontherapySte 55 HadrontherapySteppingAction::~HadrontherapySteppingAction() 56 { 56 { 57 } 57 } 58 58 59 ////////////////////////////////////////////// 59 ///////////////////////////////////////////////////////////////////////////// 60 void HadrontherapySteppingAction::UserStepping 60 void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep) 61 { 61 { 62 << 62 G4StepPoint* PreStep = aStep->GetPreStepPoint(); 63 // The followings are calls to usefuls inf << 63 G4StepPoint* PostStep = aStep->GetPostStepPoint(); 64 // Please, comment out them if want to use << 64 65 << 65 G4double PreStepX =PreStep->GetPosition().x(); 66 // G4Track* theTrack = aStep->GetTrack(); << 66 G4double PreStepY =PreStep->GetPosition().y(); 67 << 67 G4double PreStepZ =PreStep->GetPosition().z(); 68 G4StepPoint* PreStep = aStep->GetPreStepPo << 68 G4double parentID =aStep->GetTrack()->GetParentID(); 69 G4StepPoint* PostStep = aStep->GetPostStep << 69 G4double trackID =aStep->GetTrack()->GetTrackID(); 70 << 70 71 G4TouchableHandle touchPreStep = PreStep-> << 71 G4double PostStepX =PostStep->GetPosition().x(); 72 G4TouchableHandle touchPostStep = PostStep << 72 G4double PostStepY =PostStep->GetPosition().y(); 73 << 73 G4double PostStepZ =PostStep->GetPosition().z(); 74 //G4double PreStepX =PreStep->GetPosition( << 74 75 //G4double PreStepY =PreStep->GetPosition( << 75 // positions in the global coordinate system: 76 //G4double PreStepZ =PreStep->GetPosition( << 76 //G4ThreeVector posPreStep = PreStep->GetPosition(); 77 << 77 // G4ThreeVector posPostStep = PostStep->GetPosition(); 78 //G4double PostStepX =PostStep->GetPositio << 78 79 //G4double PostStepY =PostStep->GetPositio << 79 G4TouchableHandle touchPreStep = PreStep->GetTouchableHandle(); 80 //G4double PostStepZ =PostStep->GetPositi << 80 G4TouchableHandle touchPostStep = PostStep->GetTouchableHandle(); 81 << 81 //To get the current volume: 82 //To get the current volume: << 82 G4VPhysicalVolume* volumePre = touchPreStep->GetVolume(); 83 G4VPhysicalVolume* volumePre = touchPreSte << 83 G4VPhysicalVolume* volumePost =touchPostStep->GetVolume(); 84 //G4VPhysicalVolume* volumePost =touchPost << 84 85 << 85 //To get its name: 86 //To get its name: << 86 G4String namePre = volumePre->GetName(); 87 G4String namePre = volumePre->GetName(); << 87 G4String namePost; >> 88 >> 89 if(volumePost){ >> 90 namePost = volumePost->GetName(); >> 91 } >> 92 >> 93 >> 94 G4int eventNum = G4RunManager::GetRunManager() -> GetCurrentEvent() -> GetEventID(); >> 95 G4double eKin = aStep -> GetPreStepPoint() -> GetKineticEnergy(); >> 96 G4double PosX = aStep->GetTrack()->GetPosition().x(); >> 97 G4double PosY = aStep->GetTrack()->GetPosition().y(); >> 98 G4double PosZ = aStep->GetTrack()->GetPosition().z(); >> 99 G4String material= aStep -> GetTrack() -> GetMaterial() -> GetName(); >> 100 G4String volume= aStep->GetTrack()->GetVolume()->GetName(); >> 101 G4Track* theTrack = aStep->GetTrack(); >> 102 >> 103 if((namePre== "collimator")||(namePre== "PhysicExternalMagnet_1Down")||(namePre== "PhysicExternalMagnet_1")||(namePre== "PhysicMagnet_1Right")||(namePre== "PhysicMagnet_1Left")||(namePre== "PhysicExternalMagnet_2")||(namePre== "PhysicExternalMagnet_2Down")||(namePre== "PhysicMagnet_2Right")||(namePre== "PhysicMagnet_2Left")||(namePre== "PhysicExternalMagnet_3")||(namePre== "PhysicExternalMagnet_3Down")||(namePre== "PhysicMagnet_3Right")||(namePre== "PhysicMagnet_3Left")||(namePre== "PhysicExternalMagnet_4")||(namePre== "PhysicExternalMagnet_4Down")||(namePre=="physQuadChamberWall")||(namePre== "PhysicMagnet_4Right")||(namePre== "PhysicMagnet_4Left")||(namePre=="ExternalChamber")||(namePre=="collimatorFinal")||(namePre=="ExternalSlit")||(namePre=="PhysFourthQuad")||(namePre=="PhysThirdQuad")||(namePre=="PhysSecondQuad")||(namePre=="PhysFirstQuad")) >> 104 { >> 105 theTrack -> SetTrackStatus(fKillTrackAndSecondaries); >> 106 } >> 107 >> 108 // G4TransportationManager* tManager = G4TransportationManager::GetTransportationManager(); >> 109 //G4VPhysicalVolume* pW = tManager->GetParallelWorld ("DetectorROGeometry"); >> 110 >> 111 //G4Navigator* gNav = tManager->GetNavigator(pW); >> 112 >> 113 // G4VPhysicalVolume* currentVol = gNav->LocateGlobalPointAndSetup(aStep->GetTrack()->GetPosition()); >> 114 >> 115 // G4cout << "Step: " << currentVol->GetName() << " " << aStep->GetTrack()->GetPosition() << G4endl; >> 116 //G4cout << "G4LogicalVolume: = " << currentVol->GetLogicalVolume()->GetName() << G4endl; >> 117 //if (currentVol->GetLogicalVolume()->GetSensitiveDetector()) >> 118 // G4cout << "Sensitive Detector: " << currentVol->GetLogicalVolume()->GetSensitiveDetector()->GetName() << G4endl; >> 119 >> 120 >> 121 >> 122 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// >> 123 //// A METHOD TO RETRIEVE INFORMATIONS ABOUT SECONDARY ELECTRONS IN DIFFERENT POINTS OF FARADAY CUP //// >> 124 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// >> 125 88 126 89 127 90 // positions in the global coordinate syst << 128 91 //G4ThreeVector posPreStep = PreStep->Get << 129 //if( (aStep->GetTLocateGlobalPointAndSeturack()->GetVolume()->GetName() == "PVirtualMag") 92 //G4ThreeVector posPostStep = PostStep->Ge << 130 if((aStep->GetTrack()->GetVolume()->GetName()=="PVirtualMag") && aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") 93 << 131 { 94 //G4int eventNum = G4RunManager::GetRunMan << 132 std::ofstream WriteDataIn("new200.out", std::ios::app); 95 << 133 WriteDataIn << eKin << '\t' << " " 96 //G4double parentID =aStep->GetTrack()->Ge << 134 << eventNum << '\t' << " " 97 //G4double trackID =aStep->GetTrack()->Get << 135 << PosX << '\t' << " " 98 << 136 << PosY << '\t' << " " 99 G4double eKin = aStep -> GetPreStepPoint() << 137 << PosZ << '\t' << " " 100 << 138 //<< material << '\t' << " " 101 G4double PosX = aStep->GetTrack()->GetPosi << 139 // << volume << '\t' << " " 102 G4double PosY = aStep->GetTrack()->GetPosi << 140 << G4endl; 103 G4double PosZ = aStep->GetTrack()->GetPosi << 141 >> 142 >> 143 } >> 144 >> 145 // USEFULL METHODS TO RETRIEVE INFORMATION DURING THE STEPS >> 146 >> 147 /* G4cout << "ENERGIA: " << aStep->GetTrack()->GetKineticEnergy() >> 148 << " VOLUME " << aStep->GetTrack()->GetVolume()->GetName() >> 149 << " MATERIALE " << aStep -> GetTrack() -> GetMaterial() -> GetName() >> 150 << " EVENTO " << G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID() >> 151 << " POS " << aStep->GetTrack()->GetPosition().x() >> 152 << G4endl;*/ >> 153 >> 154 >> 155 if ((namePre=="PhysicEntranceWindow") && >> 156 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)) >> 157 { >> 158 std::ofstream WriteDataIn("finestra.out", std::ios::app); >> 159 WriteDataIn << eKin << '\t' << " " >> 160 << eventNum << '\t' << " " >> 161 << PreStepX << '\t' << " " >> 162 << PreStepY << '\t' << " " >> 163 << PreStepZ << '\t' << " " >> 164 << parentID << '\t' << " " >> 165 << trackID << '\t' << " " >> 166 << G4endl; >> 167 >> 168 //theTrack -> SetTrackStatus(fKillTrackAndSecondaries); >> 169 >> 170 } >> 171 >> 172 if ((namePre=="PhysicCup") && (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)) >> 173 { >> 174 >> 175 std::ofstream WriteDataBack("fondo.out", std::ios::app); >> 176 WriteDataBack << eKin << '\t' << " " >> 177 << eventNum << '\t' << " " >> 178 << PreStepX << '\t' << " " >> 179 << PreStepY << '\t' << " " >> 180 << PreStepZ << '\t' << " " >> 181 << parentID << '\t' << " " >> 182 << trackID << '\t' << " " >> 183 << G4endl; >> 184 104 185 105 G4String volume= aStep->GetTrack()->GetVo << 186 //theTrack -> SetTrackStatus(fKillTrackAndSecondaries); 106 G4Track* theTrack = aStep->GetTrack(); << 187 } >> 188 >> 189 //////////////////////////////////// VIRTUAL WINDOW ////////////////////////////////////////////////// >> 190 >> 191 if (((namePost=="PhysicVirtualWindow") && (namePre!="PhysicVirtualWindow")) && >> 192 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)) //To check that the particle has just entered in the current volume (i.e. it is at the first step in the volume; the preStepPoint is at the boundary): >> 193 { //(namePost=="VirtualWindow") && >> 194 >> 195 if ((PostStepX - PreStepX)>0) >> 196 { >> 197 >> 198 //To check that the particle is leaving the current volume (i.e. it is at the last step in the volume; the postStepPoint is at the boundary): >> 199 // if (PostStep->GetStepStatus() == fGeomBoundary) >> 200 >> 201 >> 202 >> 203 >> 204 std::ofstream WriteDataIn("DatiFCWindowIn.out", std::ios::app); >> 205 WriteDataIn << eKin << '\t' << " " >> 206 << eventNum << '\t' << " " >> 207 << PostStepX << '\t' << " " >> 208 << PostStepY << '\t' << " " >> 209 << PostStepZ << '\t' << " " >> 210 //<< direction.x() << '\t' << " " >> 211 //<< direction.y() << '\t' << " " >> 212 //<< direction.z() << '\t' << " " >> 213 << parentID << '\t' << " " >> 214 << trackID << '\t' << " " >> 215 << G4endl; >> 216 >> 217 } >> 218 >> 219 else //if ((PostStepX - PreStepX)<0) >> 220 { >> 221 >> 222 //G4cout<<"ecco il nome del volume nella condizione Out "<< namePre<<" "<<namePost<<G4endl; >> 223 >> 224 std::ofstream WriteDataBack("DatiFCWindowBack.out", std::ios::app); >> 225 WriteDataBack << eKin << '\t' << " " >> 226 << eventNum << '\t' << " " >> 227 << PostStepX << '\t' << " " >> 228 << PostStepY << '\t' << " " >> 229 << PostStepZ << '\t' << " " >> 230 << parentID << '\t' << " " >> 231 << trackID << '\t' << " " >> 232 << G4endl; >> 233 >> 234 } >> 235 } >> 236 >> 237 >> 238 >> 239 ///////////////////////////////////////// proton in window ////////////////////////////////////////////////// >> 240 >> 241 /* >> 242 if ((namePost=="SpectrumPVolume") && (namePre!="SpectrumPVolume") && >> 243 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton"))//&& (PreStep->GetStepStatus() == fGeomBoundary)To check that the particle has just entered in the current volume (i.e. it is at the first step in the volume; the preStepPoint is at the boundary): >> 244 { //(namePost=="VirtualVolumeWindow") && >> 245 107 246 108 //G4String material= aStep -> GetTrack() - << 247 109 //G4cout << "material " << material << G << 248 //To check that the particle is leaving the current volume (i.e. it is at the last step in the volume; the postStepPoint is at the boundary): 110 //G4String volume= aStep->GetTrack()->Get << 249 // if (PostStep->GetStepStatus() == fGeomBoundary) 111 //G4String pvname= pv-> GetName(); << 250 >> 251 >> 252 */ >> 253 /* >> 254 >> 255 std::ofstream WriteDataP("DatiFCWindowProton.out", std::ios::app); >> 256 WriteDataP << eKin << '\t' << " " >> 257 << eventNum << '\t' << " " >> 258 << PostStepX << '\t' << " " >> 259 << PostStepY << '\t' << " " >> 260 << PostStepZ << '\t' << " " >> 261 << parentID << '\t' << " " >> 262 << trackID << '\t' << " " >> 263 << G4endl; >> 264 >> 265 } >> 266 */ >> 267 >> 268 if (((namePost=="PhysicVirtualWindow") && (namePre!="PhysicVirtualWindow")) && >> 269 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton"))//&& (PreStep->GetStepStatus() == fGeomBoundary)) //To check that the particle has just entered in the current volume (i.e. it is at the first step in the volume; the preStepPoint is at the boundary): >> 270 { //(namePost=="VirtualVolumeWindow") && >> 271 112 272 113 G4String particleName = aStep->GetTrack()- << 273 >> 274 //To check that the particle is leaving the current volume (i.e. it is at the last step in the volume; the postStepPoint is at the boundary): >> 275 // if (PostStep->GetStepStatus() == fGeomBoundary) >> 276 >> 277 >> 278 >> 279 >> 280 // >> 281 >> 282 std::ofstream WriteData("DatiFCAfterWindowProton.out", std::ios::app); >> 283 WriteData << eKin << '\t' << " " >> 284 << eventNum << '\t' << " " >> 285 << PostStepX << '\t' << " " >> 286 << PostStepY << '\t' << " " >> 287 << PostStepZ << '\t' << " " >> 288 << parentID << '\t' << " " >> 289 << trackID << '\t' << " " >> 290 << G4endl; >> 291 >> 292 } >> 293 >> 294 >> 295 114 296 115 G4double momentumX = aStep->GetTrack()->G << 297 ////////////////////////////////////////////////// GUARD RING /////////////////////////////////////////////////////////////// 116 G4double momentumY = aStep->GetTrack()->G << 298 117 G4double momentumZ = aStep->GetTrack()->G << 118 299 >> 300 if ((namePre=="PhysicGuardRing") && >> 301 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary)) >> 302 { >> 303 >> 304 if ((PostStepX - PreStepX)>0) >> 305 { >> 306 >> 307 >> 308 std::ofstream WriteDataIn("DatiFCGuardRingIn.out", std::ios::app); >> 309 WriteDataIn << eKin << '\t' << " " >> 310 << eventNum << '\t' << " " >> 311 << PreStepX << '\t' << " " >> 312 << PreStepY << '\t' << " " >> 313 << PreStepZ << '\t' << " " >> 314 << parentID << '\t' << " " >> 315 << G4endl; >> 316 >> 317 } >> 318 >> 319 else >> 320 { >> 321 >> 322 >> 323 std::ofstream WriteDataBack("DatiFCGuardRingBack.out", std::ios::app); >> 324 WriteDataBack << eKin << '\t' << " " >> 325 << eventNum << '\t' << " " >> 326 << PreStepX << '\t' << " " >> 327 << PreStepY << '\t' << " " >> 328 << PreStepZ << '\t' << " " >> 329 << parentID << '\t' << " " >> 330 << G4endl; >> 331 >> 332 } >> 333 } >> 334 >> 335 //////////////////////////////////////// VIRTUAL MIDDLE /////////////////////////////////////////////////////////////// 119 336 120 G4ParticleDefinition *particleDef = theTra << 337 if (((namePost=="PhysicVirtualMiddle") && (namePre!="PhysicVirtualMiddle")) && 121 G4int pdg = particleDef ->GetPDGEncoding() << 338 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)) >> 339 { >> 340 >> 341 if ((PostStepX - PreStepX)>0) >> 342 { >> 343 >> 344 >> 345 // >> 346 >> 347 std::ofstream WriteMDataIn("DatiFCMiddleIn.out", std::ios::app); >> 348 WriteMDataIn << eKin << '\t' << " " >> 349 << eventNum << '\t' << " " >> 350 << PostStepX << '\t' << " " >> 351 << PostStepY << '\t' << " " >> 352 << PostStepZ << '\t' << " " >> 353 << parentID << '\t' << " " >> 354 << trackID << '\t' << " " >> 355 << G4endl; >> 356 >> 357 } >> 358 >> 359 else >> 360 { >> 361 >> 362 std::ofstream WriteMDataBack("DatiFCMiddleBack.out", std::ios::app); >> 363 WriteMDataBack << eKin << '\t' << " " >> 364 << eventNum << '\t' << " " >> 365 << PostStepX << '\t' << " " >> 366 << PostStepY << '\t' << " " >> 367 << PostStepZ << '\t' << " " >> 368 << parentID << '\t' << " " >> 369 << trackID << '\t' << " " >> 370 << G4endl; >> 371 >> 372 } >> 373 } >> 374 >> 375 /////////////////////////////////////// VIRTUAL BOTTOM /////////////////////////////////////////////////////////////// >> 376 >> 377 if (((namePost=="PhysicVirtualBottom") && (namePre!="PhysicVirtualBottom")) && >> 378 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary)) >> 379 { >> 380 >> 381 if ((PostStepX - PreStepX)>0) >> 382 { >> 383 >> 384 std::ofstream WriteDataIn("DatiFCBottomIn.out", std::ios::app); >> 385 WriteDataIn << eKin << '\t' << " " >> 386 << eventNum << '\t' << " " >> 387 << PostStepX << '\t' << " " >> 388 << PostStepY << '\t' << " " >> 389 << PostStepZ << '\t' << " " >> 390 << parentID << '\t' << " " >> 391 << trackID << '\t' << " " >> 392 << G4endl; >> 393 >> 394 } >> 395 >> 396 else >> 397 { >> 398 >> 399 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl; >> 400 >> 401 std::ofstream WriteDataBack("DatiFCBottomBack.out", std::ios::app); >> 402 WriteDataBack << eKin << '\t' << " " >> 403 << eventNum << '\t' << " " >> 404 << PostStepX << '\t' << " " >> 405 << PostStepY << '\t' << " " >> 406 << PostStepZ << '\t' << " " >> 407 //<< direction.x() << '\t' << " " >> 408 //<< direction.y() << '\t' << " " >> 409 //<< direction.z() << '\t' << " " >> 410 << parentID << '\t' << " " >> 411 << trackID << '\t' << " " >> 412 << G4endl; >> 413 >> 414 } >> 415 } >> 416 >> 417 >> 418 >> 419 ////////////////////////////////////////////////// Dose Calculation /////////////////////////////////////////////////////////////// >> 420 >> 421 //namePre!=("Cup") >> 422 if (((namePost=="PhysicCup") && (namePre!="PhysicCup")) && >> 423 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton") && (PreStep->GetStepStatus() == fGeomBoundary)) //|| (namePre=="Cup") && ) >> 424 { >> 425 >> 426 std::ofstream Carica("CaricaRaccolta.out", std::ios::app); >> 427 Carica << eKin << '\t' << " " >> 428 << eventNum << '\t' << " " >> 429 << PostStepX << '\t' << " " >> 430 << PostStepY << '\t' << " " >> 431 << PostStepZ << '\t' << " " >> 432 << parentID << '\t' << " " >> 433 << trackID << '\t' << " " >> 434 << G4endl; >> 435 >> 436 >> 437 } >> 438 >> 439 //namePre!=("PhysicFaradayCupBottom") >> 440 if (((namePost=="PhysicFaradayCupBottom") && (namePre!="PhysicFaradayCupBottom"))&& >> 441 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton") && (PreStep->GetStepStatus() == fGeomBoundary)) // || (namePre=="cup") && ) >> 442 { >> 443 >> 444 std::ofstream CaricaLatFC("CaricaRaccoltaLatFC.out", std::ios::app); >> 445 CaricaLatFC << eKin << '\t' << " " >> 446 << eventNum << '\t' << " " >> 447 << PostStepX << '\t' << " " >> 448 << PostStepY << '\t' << " " >> 449 << PostStepZ << '\t' << " " >> 450 << parentID << '\t' << " " >> 451 << trackID << '\t' << " " >> 452 << G4endl; >> 453 >> 454 >> 455 } >> 456 >> 457 >> 458 >> 459 if (((namePre=="PhysicFaradayCupBottom") || (namePre=="PhysicCup")) && ((aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary))) >> 460 { >> 461 >> 462 if ((PostStepX - PreStepX)>0) >> 463 { >> 464 >> 465 std::ofstream WriteDataIn("DatiFCConeBottomCupIn.out", std::ios::app); >> 466 WriteDataIn << eKin << '\t' << " " >> 467 << eventNum << '\t' << " " >> 468 << PreStepX << '\t' << " " >> 469 << PreStepY << '\t' << " " >> 470 << PreStepZ << '\t' << " " >> 471 << parentID << '\t' << " " >> 472 << G4endl; >> 473 >> 474 } >> 475 >> 476 else >> 477 { >> 478 >> 479 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl; >> 480 >> 481 std::ofstream WriteDataBack("DatiFCConeBottomCupBack.out", std::ios::app); >> 482 WriteDataBack << eKin << '\t' << " " >> 483 << eventNum << '\t' << " " >> 484 << PreStepX << '\t' << " " >> 485 << PreStepY << '\t' << " " >> 486 << PreStepZ << '\t' << " " >> 487 << parentID << '\t' << " " >> 488 << G4endl; >> 489 >> 490 } >> 491 } >> 492 >> 493 ///////////////////////////////////////////////// VIRTUAL OVER BOTTOM /////////////////////////////////////////////////////////////// 122 494 123 if(namePre == "VirtualLayer") << 495 >> 496 if ((namePre=="PhysicVirtualOverBottom") && >> 497 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary)) >> 498 { >> 499 >> 500 if ((PostStepX - PreStepX)>0) >> 501 { >> 502 >> 503 >> 504 // >> 505 >> 506 std::ofstream WriteDataIn("DatiFCOverBottomIn.out", std::ios::app); >> 507 WriteDataIn << eKin << '\t' << " " >> 508 << eventNum << '\t' << " " >> 509 << PreStepX << '\t' << " " >> 510 << PreStepY << '\t' << " " >> 511 << PreStepZ << '\t' << " " >> 512 << parentID << '\t' << " " >> 513 << trackID << '\t' << " " >> 514 << G4endl; >> 515 >> 516 } >> 517 >> 518 else >> 519 { >> 520 >> 521 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl; >> 522 >> 523 std::ofstream WriteDataBack("DatiFCOverBottomBack.out", std::ios::app); >> 524 WriteDataBack << eKin << '\t' << " " >> 525 << eventNum << '\t' << " " >> 526 << PreStepX << '\t' << " " >> 527 << PreStepY << '\t' << " " >> 528 << PreStepZ << '\t' << " " >> 529 << parentID << '\t' << " " >> 530 << trackID << '\t' << " " >> 531 << G4endl; >> 532 >> 533 } >> 534 >> 535 >> 536 } >> 537 >> 538 >> 539 ///////////////////////////////////////////// VIRTUAL LATERAL /////////////////////////////////////////////////////////////// >> 540 >> 541 >> 542 if (((namePost=="PhysicVirtualLateral") && (namePre!="PhysicVirtualLateral"))&& >> 543 (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)) >> 544 { >> 545 >> 546 >> 547 std::ofstream WriteDataIn("DatiFCLateral.out", std::ios::app); >> 548 WriteDataIn << eKin << '\t' << " " >> 549 << eventNum << '\t' << " " >> 550 << PostStepX << '\t' << " " >> 551 << PostStepY << '\t' << " " >> 552 << PostStepZ << '\t' << " " >> 553 << parentID << '\t' << " " >> 554 << trackID << '\t' << " " >> 555 << G4endl; >> 556 >> 557 >> 558 } >> 559 >> 560 >> 561 >> 562 if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){ >> 563 #ifdef G4ANALYSIS_USE_ROOT >> 564 G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition(); >> 565 G4double secondaryParticleKineticEnergy = aStep->GetTrack()->GetKineticEnergy(); >> 566 G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei >> 567 G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus" >> 568 if(particleType == "nucleus") { >> 569 G4int A = def->GetBaryonNumber(); >> 570 G4double Z = def->GetPDGCharge(); >> 571 G4double posX = aStep->GetTrack()->GetPosition().x() / cm; >> 572 G4double posY = aStep->GetTrack()->GetPosition().y() / cm; >> 573 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm; >> 574 G4double energy = secondaryParticleKineticEnergy / A / MeV; >> 575 >> 576 HadrontherapyAnalysisManager* analysisMgr = HadrontherapyAnalysisManager::GetInstance(); >> 577 analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ); >> 578 } else if(particleName == "proton") { // proton (hydrogen-1) is a special case >> 579 G4double posX = aStep->GetTrack()->GetPosition().x() / cm ; >> 580 G4double posY = aStep->GetTrack()->GetPosition().y() / cm ; >> 581 G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ; >> 582 G4double energy = secondaryParticleKineticEnergy * MeV; // Hydrogen-1: A = 1, Z = 1 >> 583 HadrontherapyAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ); >> 584 } >> 585 >> 586 G4String secondaryParticleName = def -> GetParticleName(); >> 587 //G4cout <<"Particle: " << secondaryParticleName << G4endl; >> 588 //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl; >> 589 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); >> 590 //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this. >> 591 if(secondaryParticleName == "proton") { >> 592 analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV); >> 593 } >> 594 if(secondaryParticleName == "deuteron") { >> 595 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV); >> 596 } >> 597 if(secondaryParticleName == "triton") { >> 598 analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV); >> 599 } >> 600 if(secondaryParticleName == "alpha") { >> 601 analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV); >> 602 } >> 603 if(secondaryParticleName == "He3"){ >> 604 analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV); >> 605 } >> 606 #endif >> 607 >> 608 aStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries); >> 609 } >> 610 >> 611 // Electromagnetic and hadronic processes of primary particles in the phantom >> 612 //setting phantomPhys correctly will break something here fixme >> 613 if ((aStep -> GetTrack() -> GetTrackID() == 1) && >> 614 (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") && >> 615 (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL)) 124 { 616 { 125 std::ofstream WriteDataIn("Virtual_Lay << 617 G4String process = aStep -> GetPostStepPoint() -> 126 WriteDataIn << 618 GetProcessDefinedStep() -> GetProcessName(); 127 << 619 128 << eKin <<" " // 1 << 620 if ((process == "Transportation") || (process == "StepLimiter")) {;} 129 << PosX <<" " // 2 << 621 else 130 << PosY <<" " // 3 << 622 { 131 << PosZ <<" " // 4 << 623 if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni")) 132 << momentumX <<" " // 5 << 624 { 133 << momentumY <<" " // 6 << 625 runAction -> AddEMProcess(); 134 << momentumZ <<" " // 7 << 626 } 135 << pdg << 627 else 136 //<< theTrack << '\t' << " << 628 { 137 << 629 runAction -> AddHadronicProcess(); 138 << G4endl; << 630 139 << 631 if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") ) 140 theTrack -> SetTrackStatus(fKillTrackA << 632 G4cout << "Warning! Unknown proton process: "<< process << G4endl; 141 << 633 } 142 << 634 } >> 635 } >> 636 >> 637 // Retrieve information about the secondary particles originated in the phantom >> 638 >> 639 G4SteppingManager* steppingManager = fpSteppingManager; >> 640 >> 641 // check if it is alive >> 642 //if(theTrack-> GetTrackStatus() == fAlive) { return; } >> 643 >> 644 // Retrieve the secondary particles >> 645 G4TrackVector* fSecondary = steppingManager -> GetfSecondary(); >> 646 >> 647 for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++) >> 648 { >> 649 G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); >> 650 >> 651 if (volumeName == "phantomPhys") >> 652 { >> 653 #ifdef G4ANALYSIS_USE_ROOT >> 654 G4String secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); >> 655 G4double secondaryParticleKineticEnergy = (*fSecondary)[lp1] -> GetKineticEnergy(); >> 656 >> 657 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance(); >> 658 >> 659 if (secondaryParticleName == "e-") >> 660 analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV); >> 661 >> 662 if (secondaryParticleName == "gamma") >> 663 analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV); >> 664 >> 665 if (secondaryParticleName == "deuteron") >> 666 analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV); >> 667 >> 668 if (secondaryParticleName == "triton") >> 669 analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV); >> 670 >> 671 if (secondaryParticleName == "alpha") >> 672 analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV); >> 673 >> 674 G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge(); >> 675 if (z > 0.) >> 676 { >> 677 G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber(); >> 678 G4int electronOccupancy = (*fSecondary)[lp1] -> GetDynamicParticle() -> GetTotalOccupancy(); >> 679 >> 680 // If a generic ion is originated in the detector, its baryonic number, PDG charge, >> 681 // total number of electrons in the orbitals are stored in a ntuple >> 682 analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV); >> 683 } >> 684 #endif >> 685 } 143 } 686 } 144 << 145 << 146 << 147 << 148 } 687 } >> 688 >> 689 149 690 150 691 151 692 152 693