Geant4  10.01
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 
158 
159 
160 
161 
162 
163 if ((namePre=="PhysicEntranceWindow") &&
164  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary))
165 {
166  std::ofstream WriteDataIn("finestra.out", std::ios::app);
167  WriteDataIn << eKin << '\t' << " "
168  << eventNum << '\t' << " "
169  << PreStepX << '\t' << " "
170  << PreStepY << '\t' << " "
171  << PreStepZ << '\t' << " "
172  << parentID << '\t' << " "
173  << trackID << '\t' << " "
174  << G4endl;
175 
176  //theTrack -> SetTrackStatus(fKillTrackAndSecondaries);
177 
178  }
179 
180  if ((namePre=="PhysicCup") && (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary))
181 {
182 
183  std::ofstream WriteDataBack("fondo.out", std::ios::app);
184  WriteDataBack << eKin << '\t' << " "
185  << eventNum << '\t' << " "
186  << PreStepX << '\t' << " "
187  << PreStepY << '\t' << " "
188  << PreStepZ << '\t' << " "
189  << parentID << '\t' << " "
190  << trackID << '\t' << " "
191  << G4endl;
192 
193 
194  //theTrack -> SetTrackStatus(fKillTrackAndSecondaries);
195  }
196 
198 
199  if (((namePost=="PhysicVirtualWindow") && (namePre!="PhysicVirtualWindow")) &&
200  (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):
201 { //(namePost=="VirtualWindow") &&
202 
203  if ((PostStepX - PreStepX)>0)
204  {
205 
206 //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):
207 // if (PostStep->GetStepStatus() == fGeomBoundary)
208 
209 
210 
211 
212  std::ofstream WriteDataIn("DatiFCWindowIn.out", std::ios::app);
213  WriteDataIn << eKin << '\t' << " "
214  << eventNum << '\t' << " "
215  << PostStepX << '\t' << " "
216  << PostStepY << '\t' << " "
217  << PostStepZ << '\t' << " "
218  //<< direction.x() << '\t' << " "
219  //<< direction.y() << '\t' << " "
220  //<< direction.z() << '\t' << " "
221  << parentID << '\t' << " "
222  << trackID << '\t' << " "
223  << G4endl;
224 
225  }
226 
227  else //if ((PostStepX - PreStepX)<0)
228  {
229 
230 //G4cout<<"ecco il nome del volume nella condizione Out "<< namePre<<" "<<namePost<<G4endl;
231 
232  std::ofstream WriteDataBack("DatiFCWindowBack.out", std::ios::app);
233  WriteDataBack << eKin << '\t' << " "
234  << eventNum << '\t' << " "
235  << PostStepX << '\t' << " "
236  << PostStepY << '\t' << " "
237  << PostStepZ << '\t' << " "
238  << parentID << '\t' << " "
239  << trackID << '\t' << " "
240  << G4endl;
241 
242  }
243 }
244 
245 
246 
248 
249 /*
250 if ((namePost=="SpectrumPVolume") && (namePre!="SpectrumPVolume") &&
251  (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):
252 { //(namePost=="VirtualVolumeWindow") &&
253 
254 
255 
256 //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):
257 // if (PostStep->GetStepStatus() == fGeomBoundary)
258 
259 
260 */
261 /*
262 
263  std::ofstream WriteDataP("DatiFCWindowProton.out", std::ios::app);
264  WriteDataP << eKin << '\t' << " "
265  << eventNum << '\t' << " "
266  << PostStepX << '\t' << " "
267  << PostStepY << '\t' << " "
268  << PostStepZ << '\t' << " "
269  << parentID << '\t' << " "
270  << trackID << '\t' << " "
271  << G4endl;
272 
273  }
274 */
275 
276 if (((namePost=="PhysicVirtualWindow") && (namePre!="PhysicVirtualWindow")) &&
277  (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):
278 { //(namePost=="VirtualVolumeWindow") &&
279 
280 
281 
282 //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):
283 // if (PostStep->GetStepStatus() == fGeomBoundary)
284 
285 
286 
287 
288 //
289 
290  std::ofstream WriteData("DatiFCAfterWindowProton.out", std::ios::app);
291  WriteData << eKin << '\t' << " "
292  << eventNum << '\t' << " "
293  << PostStepX << '\t' << " "
294  << PostStepY << '\t' << " "
295  << PostStepZ << '\t' << " "
296  << parentID << '\t' << " "
297  << trackID << '\t' << " "
298  << G4endl;
299 
300  }
301 
302 
303 
304 
306 
307 
308 if ((namePre=="PhysicGuardRing") &&
309  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary))
310 {
311 
312  if ((PostStepX - PreStepX)>0)
313  {
314 
315 
316  std::ofstream WriteDataIn("DatiFCGuardRingIn.out", std::ios::app);
317  WriteDataIn << eKin << '\t' << " "
318  << eventNum << '\t' << " "
319  << PreStepX << '\t' << " "
320  << PreStepY << '\t' << " "
321  << PreStepZ << '\t' << " "
322  << parentID << '\t' << " "
323  << G4endl;
324 
325  }
326 
327  else
328  {
329 
330 
331  std::ofstream WriteDataBack("DatiFCGuardRingBack.out", std::ios::app);
332  WriteDataBack << eKin << '\t' << " "
333  << eventNum << '\t' << " "
334  << PreStepX << '\t' << " "
335  << PreStepY << '\t' << " "
336  << PreStepZ << '\t' << " "
337  << parentID << '\t' << " "
338  << G4endl;
339 
340  }
341 }
342 
344 
345  if (((namePost=="PhysicVirtualMiddle") && (namePre!="PhysicVirtualMiddle")) &&
346  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary))
347 {
348 
349  if ((PostStepX - PreStepX)>0)
350  {
351 
352 
353 //
354 
355  std::ofstream WriteMDataIn("DatiFCMiddleIn.out", std::ios::app);
356  WriteMDataIn << eKin << '\t' << " "
357  << eventNum << '\t' << " "
358  << PostStepX << '\t' << " "
359  << PostStepY << '\t' << " "
360  << PostStepZ << '\t' << " "
361  << parentID << '\t' << " "
362  << trackID << '\t' << " "
363  << G4endl;
364 
365  }
366 
367  else
368  {
369 
370  std::ofstream WriteMDataBack("DatiFCMiddleBack.out", std::ios::app);
371  WriteMDataBack << eKin << '\t' << " "
372  << eventNum << '\t' << " "
373  << PostStepX << '\t' << " "
374  << PostStepY << '\t' << " "
375  << PostStepZ << '\t' << " "
376  << parentID << '\t' << " "
377  << trackID << '\t' << " "
378  << G4endl;
379 
380  }
381 }
382 
384 
385  if (((namePost=="PhysicVirtualBottom") && (namePre!="PhysicVirtualBottom")) &&
386  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary))
387 {
388 
389  if ((PostStepX - PreStepX)>0)
390  {
391 
392  std::ofstream WriteDataIn("DatiFCBottomIn.out", std::ios::app);
393  WriteDataIn << eKin << '\t' << " "
394  << eventNum << '\t' << " "
395  << PostStepX << '\t' << " "
396  << PostStepY << '\t' << " "
397  << PostStepZ << '\t' << " "
398  << parentID << '\t' << " "
399  << trackID << '\t' << " "
400  << G4endl;
401 
402  }
403 
404  else
405  {
406 
407 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl;
408 
409  std::ofstream WriteDataBack("DatiFCBottomBack.out", std::ios::app);
410  WriteDataBack << eKin << '\t' << " "
411  << eventNum << '\t' << " "
412  << PostStepX << '\t' << " "
413  << PostStepY << '\t' << " "
414  << PostStepZ << '\t' << " "
415  //<< direction.x() << '\t' << " "
416  //<< direction.y() << '\t' << " "
417  //<< direction.z() << '\t' << " "
418  << parentID << '\t' << " "
419  << trackID << '\t' << " "
420  << G4endl;
421 
422  }
423 }
424 
425 
426 
428 
429 //namePre!=("Cup")
430  if (((namePost=="PhysicCup") && (namePre!="PhysicCup")) &&
431  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton") && (PreStep->GetStepStatus() == fGeomBoundary)) //|| (namePre=="Cup") && )
432 {
433 
434  std::ofstream Carica("CaricaRaccolta.out", std::ios::app);
435  Carica << eKin << '\t' << " "
436  << eventNum << '\t' << " "
437  << PostStepX << '\t' << " "
438  << PostStepY << '\t' << " "
439  << PostStepZ << '\t' << " "
440  << parentID << '\t' << " "
441  << trackID << '\t' << " "
442  << G4endl;
443 
444 
445 }
446 
447 //namePre!=("PhysicFaradayCupBottom")
448  if (((namePost=="PhysicFaradayCupBottom") && (namePre!="PhysicFaradayCupBottom"))&&
449  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton") && (PreStep->GetStepStatus() == fGeomBoundary)) // || (namePre=="cup") && )
450 {
451 
452  std::ofstream CaricaLatFC("CaricaRaccoltaLatFC.out", std::ios::app);
453  CaricaLatFC << eKin << '\t' << " "
454  << eventNum << '\t' << " "
455  << PostStepX << '\t' << " "
456  << PostStepY << '\t' << " "
457  << PostStepZ << '\t' << " "
458  << parentID << '\t' << " "
459  << trackID << '\t' << " "
460  << G4endl;
461 
462 
463 }
464 
465 
466 
467  if (((namePre=="PhysicFaradayCupBottom") || (namePre=="PhysicCup")) && ((aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary)))
468 {
469 
470  if ((PostStepX - PreStepX)>0)
471  {
472 
473  std::ofstream WriteDataIn("DatiFCConeBottomCupIn.out", std::ios::app);
474  WriteDataIn << eKin << '\t' << " "
475  << eventNum << '\t' << " "
476  << PreStepX << '\t' << " "
477  << PreStepY << '\t' << " "
478  << PreStepZ << '\t' << " "
479  << parentID << '\t' << " "
480  << G4endl;
481 
482  }
483 
484  else
485  {
486 
487 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl;
488 
489  std::ofstream WriteDataBack("DatiFCConeBottomCupBack.out", std::ios::app);
490  WriteDataBack << eKin << '\t' << " "
491  << eventNum << '\t' << " "
492  << PreStepX << '\t' << " "
493  << PreStepY << '\t' << " "
494  << PreStepZ << '\t' << " "
495  << parentID << '\t' << " "
496  << G4endl;
497 
498  }
499 }
500 
502 
503 
504 if ((namePre=="PhysicVirtualOverBottom") &&
505  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-")&& (PreStep->GetStepStatus() == fGeomBoundary))
506 {
507 
508  if ((PostStepX - PreStepX)>0)
509  {
510 
511 
512 //
513 
514  std::ofstream WriteDataIn("DatiFCOverBottomIn.out", std::ios::app);
515  WriteDataIn << eKin << '\t' << " "
516  << eventNum << '\t' << " "
517  << PreStepX << '\t' << " "
518  << PreStepY << '\t' << " "
519  << PreStepZ << '\t' << " "
520  << parentID << '\t' << " "
521  << trackID << '\t' << " "
522  << G4endl;
523 
524  }
525 
526  else
527  {
528 
529 //G4cout<<"ecco il nome del volume nella condizione Back "<< namePre<<" "<<namePost<<G4endl;
530 
531  std::ofstream WriteDataBack("DatiFCOverBottomBack.out", std::ios::app);
532  WriteDataBack << eKin << '\t' << " "
533  << eventNum << '\t' << " "
534  << PreStepX << '\t' << " "
535  << PreStepY << '\t' << " "
536  << PreStepZ << '\t' << " "
537  << parentID << '\t' << " "
538  << trackID << '\t' << " "
539  << G4endl;
540 
541  }
542 
543 
544 }
545 
546 
548 
549 
550 if (((namePost=="PhysicVirtualLateral") && (namePre!="PhysicVirtualLateral"))&&
551  (aStep->GetTrack()->GetDefinition()->GetParticleName() == "e-") && (PreStep->GetStepStatus() == fGeomBoundary))
552 {
553 
554 
555  std::ofstream WriteDataIn("DatiFCLateral.out", std::ios::app);
556  WriteDataIn << eKin << '\t' << " "
557  << eventNum << '\t' << " "
558  << PostStepX << '\t' << " "
559  << PostStepY << '\t' << " "
560  << PostStepZ << '\t' << " "
561  << parentID << '\t' << " "
562  << trackID << '\t' << " "
563  << G4endl;
564 
565 
566  }
567 
568 
569 
570  if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){
571 #ifdef G4ANALYSIS_USE_ROOT
572  G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
573  G4double secondaryParticleKineticEnergy = aStep->GetTrack()->GetKineticEnergy();
574  G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
575  G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
576  if(particleType == "nucleus") {
577  G4int A = def->GetBaryonNumber();
578  G4double Z = def->GetPDGCharge();
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 / A / MeV;
583 
584  HadrontherapyAnalysisManager* analysisMgr = HadrontherapyAnalysisManager::GetInstance();
585  analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ);
586  } else if(particleName == "proton") { // proton (hydrogen-1) is a special case
587  G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;
588  G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;
589  G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
590  G4double energy = secondaryParticleKineticEnergy * MeV; // Hydrogen-1: A = 1, Z = 1
591  HadrontherapyAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
592  }
593 
594  G4String secondaryParticleName = def -> GetParticleName();
595  //G4cout <<"Particle: " << secondaryParticleName << G4endl;
596  //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
597  HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
598  //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this.
599  if(secondaryParticleName == "proton") {
600  analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);
601  }
602  if(secondaryParticleName == "deuteron") {
603  analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
604  }
605  if(secondaryParticleName == "triton") {
606  analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
607  }
608  if(secondaryParticleName == "alpha") {
609  analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
610  }
611  if(secondaryParticleName == "He3"){
612  analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);
613  }
614 #endif
615 
616  aStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
617  }
618 
619  // Electromagnetic and hadronic processes of primary particles in the phantom
620  //setting phantomPhys correctly will break something here fixme
621  if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
622  (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
623  (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
624  {
625  G4String process = aStep -> GetPostStepPoint() ->
626  GetProcessDefinedStep() -> GetProcessName();
627 
628  if ((process == "Transportation") || (process == "StepLimiter")) {;}
629  else
630  {
631  if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
632  {
633  runAction -> AddEMProcess();
634  }
635  else
636  {
637  runAction -> AddHadronicProcess();
638 
639  if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
640  G4cout << "Warning! Unknown proton process: "<< process << G4endl;
641  }
642  }
643  }
644 
645  // Retrieve information about the secondary particles originated in the phantom
646 
647  G4SteppingManager* steppingManager = fpSteppingManager;
648 
649  // check if it is alive
650  //if(theTrack-> GetTrackStatus() == fAlive) { return; }
651 
652  // Retrieve the secondary particles
653  G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
654 
655  for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
656  {
657  G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
658 
659  if (volumeName == "phantomPhys")
660  {
661 #ifdef G4ANALYSIS_USE_ROOT
662  G4String secondaryParticleName = (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();
663  G4double secondaryParticleKineticEnergy = (*fSecondary)[lp1] -> GetKineticEnergy();
664 
665  HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
666 
667  if (secondaryParticleName == "e-")
668  analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
669 
670  if (secondaryParticleName == "gamma")
671  analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
672 
673  if (secondaryParticleName == "deuteron")
674  analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
675 
676  if (secondaryParticleName == "triton")
677  analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
678 
679  if (secondaryParticleName == "alpha")
680  analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
681 
682  G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
683  if (z > 0.)
684  {
685  G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
686  G4int electronOccupancy = (*fSecondary)[lp1] -> GetDynamicParticle() -> GetTotalOccupancy();
687 
688  // If a generic ion is originated in the detector, its baryonic number, PDG charge,
689  // total number of electrons in the orbitals are stored in a ntuple
690  analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
691  }
692 #endif
693  }
694  }
695 }
696 
697 
698 
699 
700 
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