| 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