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