Geant4  10.03
HadrontherapySteppingAction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
29 #include "G4SteppingManager.hh"
30 #include "G4TrackVector.hh"
31 #include "HadrontherapySteppingAction.hh"
32 #include "G4ios.hh"
33 #include "G4SteppingManager.hh"
34 #include "G4Track.hh"
35 #include "G4Step.hh"
36 #include "G4StepPoint.hh"
37 #include "G4TrackStatus.hh"
38 #include "G4TrackVector.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4ParticleTypes.hh"
41 #include "G4UserEventAction.hh"
42 #include "G4TransportationManager.hh"
43 #include "G4VSensitiveDetector.hh"
44 #include "HadrontherapyRunAction.hh"
45 #include "HadrontherapyAnalysisManager.hh"
46 #include "G4SystemOfUnits.hh"
47 
49 HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run)
50 {
51  runAction = run;
52 }
53 
55 HadrontherapySteppingAction::~HadrontherapySteppingAction()
56 {
57 }
58 
60 void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep)
61 {
62 G4StepPoint* PreStep = aStep->GetPreStepPoint();
63 G4StepPoint* PostStep = aStep->GetPostStepPoint();
64 
65 G4double PreStepX =PreStep->GetPosition().x();
66 G4double PreStepY =PreStep->GetPosition().y();
67 G4double PreStepZ =PreStep->GetPosition().z();
68 G4double parentID =aStep->GetTrack()->GetParentID();
69 G4double trackID =aStep->GetTrack()->GetTrackID();
70 
71 G4double PostStepX =PostStep->GetPosition().x();
72 G4double PostStepY =PostStep->GetPosition().y();
73 G4double PostStepZ =PostStep->GetPosition().z();
74 
75 // positions in the global coordinate system:
76 //G4ThreeVector posPreStep = PreStep->GetPosition();
77 // G4ThreeVector posPostStep = PostStep->GetPosition();
78 
79 G4TouchableHandle touchPreStep = PreStep->GetTouchableHandle();
80 G4TouchableHandle touchPostStep = PostStep->GetTouchableHandle();
81 //To get the current volume:
82 G4VPhysicalVolume* volumePre = touchPreStep->GetVolume();
83 G4VPhysicalVolume* volumePost =touchPostStep->GetVolume();
84 
85 //To get its name:
86 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 
125 
126 
127 
128 
129  //if( (aStep->GetTLocateGlobalPointAndSeturack()->GetVolume()->GetName() == "PVirtualMag")
130  if((aStep->GetTrack()->GetVolume()->GetName()=="PVirtualMag") && aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")
131  {
132  std::ofstream WriteDataIn("new200.out", std::ios::app);
133  WriteDataIn << eKin << '\t' << " "
134  << eventNum << '\t' << " "
135  << PosX << '\t' << " "
136  << PosY << '\t' << " "
137  << PosZ << '\t' << " "
138  //<< material << '\t' << " "
139  // << volume << '\t' << " "
140  << G4endl;
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 
185 
186  //theTrack -> SetTrackStatus(fKillTrackAndSecondaries);
187  }
188 
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 
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 
246 
247 
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):
249 // if (PostStep->GetStepStatus() == fGeomBoundary)
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 
272 
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 
296 
298 
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 
336 
337  if (((namePost=="PhysicVirtualMiddle") && (namePre!="PhysicVirtualMiddle")) &&
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 
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 
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 
494 
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 
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))
616  {
617  G4String process = aStep -> GetPostStepPoint() ->
618  GetProcessDefinedStep() -> GetProcessName();
619 
620  if ((process == "Transportation") || (process == "StepLimiter")) {;}
621  else
622  {
623  if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
624  {
625  runAction -> AddEMProcess();
626  }
627  else
628  {
629  runAction -> AddHadronicProcess();
630 
631  if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
632  G4cout << "Warning! Unknown proton process: "<< process << G4endl;
633  }
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  }
686  }
687 }
688 
689 
690 
691 
692 
G4Track::SetTrackStatus
void SetTrackStatus(const G4TrackStatus aTrackStatus)
G4Track::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4Track::GetParentID
G4int GetParentID() const
HadrontherapyAnalysisManager::GetInstance
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
Definition: HadrontherapyAnalysisManager.cc:66
HadrontherapySteppingAction::HadrontherapySteppingAction
HadrontherapySteppingAction(HadrontherapyRunAction *)
Definition: HadrontherapySteppingAction.cc:49
G4StepPoint.hh
HadrontherapySteppingAction.hh
G4Track.hh
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4StepPoint::GetStepStatus
G4StepStatus GetStepStatus() const
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Step.hh
HadrontherapyAnalysisManager.hh
G4int
int G4int
Definition: G4Types.hh:78
G4TrackVector.hh
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4TrackStatus.hh
G4Step::GetPreStepPoint
G4StepPoint * GetPreStepPoint() const
G4ParticleTypes.hh
G4Track::GetKineticEnergy
G4double GetKineticEnergy() const
G4VSensitiveDetector.hh
HadrontherapyRunAction.hh
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4VPhysicalVolume::GetName
const G4String & GetName() const
G4StepPoint::GetPosition
const G4ThreeVector & GetPosition() const
G4UserSteppingAction::fpSteppingManager
G4SteppingManager * fpSteppingManager
Definition: G4UserSteppingAction.hh:72
cm
static constexpr double cm
Definition: G4SIunits.hh:119
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:141
G4ParticleDefinition.hh
G4Step
Definition: G4Step.hh:76
G4Track::GetTrackID
G4int GetTrackID() const
G4TransportationManager.hh
G4UserEventAction.hh
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
fGeomBoundary
Definition: G4StepStatus.hh:54
G4TrackVector
std::vector< G4Track * > G4TrackVector
Definition: G4TrackVector.hh:49
G4VTouchable::GetVolume
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4Step::GetPostStepPoint
G4StepPoint * GetPostStepPoint() const
G4ios.hh
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
HadrontherapyAnalysisManager
A class for connecting the simulation to an analysis package.
Definition: HadrontherapyAnalysisManager.hh:49
G4SteppingManager.hh
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4Step::GetTrack
G4Track * GetTrack() const
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4StepPoint::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
fKillTrackAndSecondaries
Definition: G4TrackStatus.hh:58