Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ML2SDWithParticle.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 // The code was written by :
27 // ^Claudio Andenna claudio.andenna@ispesl.it, claudio.andenna@iss.infn.it
28 // *Barbara Caccia barbara.caccia@iss.it
29 // with the support of Pablo Cirrone (LNS, INFN Catania Italy)
30 // with the contribute of Alessandro Occhigrossi*
31 //
32 // ^INAIL DIPIA - ex ISPESL and INFN Roma, gruppo collegato Sanità, Italy
33 // *Istituto Superiore di Sanità and INFN Roma, gruppo collegato Sanità, Italy
34 // Viale Regina Elena 299, 00161 Roma (Italy)
35 // tel (39) 06 49902246
36 // fax (39) 06 49387075
37 //
38 // more information:
39 // http://g4advancedexamples.lngs.infn.it/Examples/medical-linac
40 //
41 //*******************************************************//
42 
43 #include "ML2SDWithParticle.hh"
44 #include "ML2ExpVoxels.hh"
46 #include "G4SystemOfUnits.hh"
47 
49 : G4VSensitiveDetector("killer_plane"),particles(0)
50 {
51  idType=idSD_KillerPlane;
52  bStopAtPhaseSpace=true;
53  nTotalParticles=0;
54  bActive=true;
55 }
56 
58  G4int seed, G4int nMaxPart, G4String name,
59  G4String PhaseSpaceOutFile, G4bool bSave, G4bool bStop,
60  SPrimaryParticle *pData, G4double ZPos)
61 : G4VSensitiveDetector(name),particles(0)
62 {
63  accTargetZPosition=ZPos;
64  max_N_particles_in_PhSp_File=maxPartFile;
65  nMaxParticlesInRamPhaseSpace = nMaxPart;
66  idType=id;
67  primaryParticleData=pData;
68  bActive=true;
69  nParticle=0;
70  nTotalParticles=0;
71  bSavePhaseSpace=bSave;
72  bStopAtPhaseSpace=bStop;
73  if (bSavePhaseSpace)
74  {particles=new Sparticle[nMaxPart];}
75  G4String seedName;
76  char a[10];
77  sprintf(a,"%d", seed);
78  seedName=(G4String)a;
79 
80  fullOutFileData=PhaseSpaceOutFile+"_"+seedName+".txt";
81  fullOutFileData=PhaseSpaceOutFile+"_"+seedName+".txt";
82 }
84 {
85  if (particles!=0)
86  {
87  delete [] particles;
88  }
89 }
90 void CML2SDWithParticle::saveHeaderParticles()
91 {
92  std::ofstream out;
93  out.open(fullOutFileData, std::ios::out);
94  out << "Sensitive Detector-Particles"<<G4endl;
95  out << "n Total Events,\t x [mm],\t y [mm],\t z [mm],\t dirX,\t dirY,\t dirZ,\t KinEnergy [MeV],\t part Type,\t primary part type,\t nPrimaryPart" << G4endl;
96  out.close();
97 }
98 void CML2SDWithParticle::saveDataParticles(G4int nPart)
99 {
100  std::ofstream out;
101  out.open(fullOutFileData, std::ios::app);
102  static G4int nTotParticles=0;
103  for (int i=0; i< nPart; i++)
104  {
105  out << nTotParticles++ << '\t';
106 // out << particles[i].volumeName << '\t';
107  out << particles[i].pos.getX()/mm << '\t';
108  out << particles[i].pos.getY()/mm << '\t';
109  out << (accTargetZPosition + particles[i].pos.getZ())/mm << '\t'; // it translates the current z value in global coordinates to the accelerator local coordinates (only z)
110  out << particles[i].dir.getX() << '\t';
111  out << particles[i].dir.getY() << '\t';
112  out << particles[i].dir.getZ() << '\t';
113  out << particles[i].kinEnergy/MeV << '\t';
114  out << particles[i].partPDGE<< '\t';
115  out << particles[i].primaryParticlePDGE<< '\t';
116  out << particles[i].nPrimaryPart<< G4endl;
117  }
118  out.close();
119 }
120 
122 {
123  if (bActive && (CML2AcceleratorConstruction::GetInstance()->getPhysicalVolume()->GetRotation()->isIdentity()))
124  {
125  G4double energyKin= aStep->GetTrack()->GetKineticEnergy();
126  static bool bFirstTime=true;
127  if (idType==idSD_KillerPlane)
128  {
129  nTotalParticles++;
131  }
132  else
133  {
134  if (energyKin>0.)
135  {
136  particles[nParticle].volumeName="";
137  particles[nParticle].pos=aStep->GetPreStepPoint()->GetPosition();
138  particles[nParticle].dir=aStep->GetPreStepPoint()->GetMomentumDirection();
139  particles[nParticle].kinEnergy=energyKin;
140  particles[nParticle].partPDGE=aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
141  particles[nParticle].primaryParticlePDGE=primaryParticleData->partPDGE;
142  particles[nParticle].nPrimaryPart=primaryParticleData->nPrimaryParticle;
143  nParticle++;
144  nTotalParticles++;
145  if (nTotalParticles==max_N_particles_in_PhSp_File)
146  {
147  if (bFirstTime)
148  {
149  bFirstTime=false;
150  saveHeaderParticles();
151  }
152  saveDataParticles(nParticle);
153  nParticle=0;
154  bActive =false;// to stop the phase space creation
155  }
156  if (nParticle==nMaxParticlesInRamPhaseSpace)
157  {
158  if (bFirstTime)
159  {
160  bFirstTime=false;
161  saveHeaderParticles();
162  }
163  saveDataParticles(nParticle);
164  nParticle=0;
165  }
166 
167  Sparticle *particle=new Sparticle;
168  particle->dir=aStep->GetPreStepPoint()->GetMomentumDirection();
169  particle->pos=aStep->GetPreStepPoint()->GetPosition();
170  particle->kinEnergy=energyKin;
171  particle->nPrimaryPart=nTotalParticles; // to pass the id of this phase space particle
172  particle->partPDGE=aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
173  particle->primaryParticlePDGE=primaryParticleData->partPDGE;
174  particle->volumeId=-1;
175  particle->volumeName="-1";
176  }
177  if (bStopAtPhaseSpace)
178  {aStep->GetTrack()->SetTrackStatus(fStopAndKill);}
179  }
180  }
181  else
182  {
183  if (bStopAtPhaseSpace)
184  {aStep->GetTrack()->SetTrackStatus(fStopAndKill);}
185  }
186 
187  return true;
188 }
189 
191 {
192  if ((bActive) && (nParticle>0) && (CML2AcceleratorConstruction::GetInstance()->getPhysicalVolume()->GetRotation()->isIdentity()))
193  {saveDataParticles(nParticle);nParticle=0;}
194 }
void SetTrackStatus(const G4TrackStatus aTrackStatus)
G4ParticleDefinition * GetDefinition() const
const XML_Char * name
Definition: expat.h:151
static constexpr double mm
Definition: G4SIunits.hh:115
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static CML2AcceleratorConstruction * GetInstance(void)
G4int partPDGE
double getY() const
G4ThreeVector dir
int G4int
Definition: G4Types.hh:78
G4int nPrimaryPart
G4StepPoint * GetPreStepPoint() const
double getX() const
G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *ROHist)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
tuple app
Definition: demo.py:189
G4int primaryParticlePDGE
long seed
Definition: chem4.cc:68
const G4ThreeVector & GetPosition() const
bool G4bool
Definition: G4Types.hh:79
Definition: G4Step.hh:76
G4ThreeVector pos
G4double kinEnergy
double getZ() const
G4String volumeName
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4int volumeId
double G4double
Definition: G4Types.hh:76
G4Track * GetTrack() const