| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include "PrimaryGeneratorAction.hh"
27 #include "PrimaryGeneratorMessenger.hh"
28
29 #include "G4Box.hh"
30 #include "G4LogicalVolume.hh"
31 #include "G4LogicalVolumeStore.hh"
32 #include "G4ParticleGun.hh"
33 #include "G4ParticleTable.hh"
34 #include "G4RunManager.hh"
35 #include "Randomize.hh"
36
37 #ifdef WITHROOT
38 #include "G4MTRunManager.hh"
39 #include "G4Run.hh"
40 #include "TFile.h"
41 #include "TTreeReader.h"
42 #include "TTreeReaderValue.h"
43 #endif
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 PrimaryGeneratorAction::PrimaryGeneratorAction
48 : G4VUserPrimaryGeneratorAction(), fPartic
49 fEnvelopeBox(nullptr) {
50 G4int n_particle = 1;
51 fParticleGun = new G4ParticleGun(n_particle)
52
53 // default particle kinematic
54 G4ParticleTable *particleTable = G4ParticleT
55 G4String particleName;
56 G4ParticleDefinition *particle =
57 particleTable->FindParticle(particleName
58 fParticleGun->SetParticleDefinition(particle
59 fParticleGun->SetParticleMomentumDirection(G
60 fParticleGun->SetParticleEnergy(30. * CLHEP:
61 fParticleGun->SetParticlePosition(G4ThreeVec
62
63 fMessenger = new PrimaryGeneratorMessenger(t
64 #ifdef WITHROOT
65 if (fReadInputFile)
66 OpenInput();
67 #endif
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 #ifdef WITHROOT
73 void PrimaryGeneratorAction::OpenInput() {
74 if (fInputFile != nullptr) {
75 delete fInputFile;
76 fInputFile = nullptr;
77 }
78 fInputFile = new TFile(fPathInputFile, "READ
79 if (fInputFile == nullptr || fInputFile->IsZ
80 G4ExceptionDescription msg;
81 msg << "Input file '" << fPathInputFile <<
82 G4Exception("PrimaryGeneratorAction::Gener
83 FatalErrorInArgument, msg);
84 }
85 auto fileDirectory =
86 dynamic_cast<TDirectory *>(fInputFile->G
87 if (fileDirectory == nullptr) {
88 G4ExceptionDescription msg;
89 msg << "Input file '" << fPathInputFile
90 << "' does not contain 'VirtualDetecto
91 G4Exception("PrimaryGeneratorAction::Gener
92 "WrongInputDirectory", FatalEr
93 }
94 fHgcalReader = new TTreeReader("HGCAL", file
95 fHgcalEventId = new TTreeReaderValue<Float_t
96 fHgcalPdgId = new TTreeReaderValue<Float_t>(
97 fHgcalPosX = new TTreeReaderValue<Float_t>(*
98 fHgcalPosY = new TTreeReaderValue<Float_t>(*
99 fHgcalMomX = new TTreeReaderValue<Float_t>(*
100 fHgcalMomY = new TTreeReaderValue<Float_t>(*
101 fHgcalMomZ = new TTreeReaderValue<Float_t>(*
102
103 return;
104 }
105 #endif
106
107 //....oooOO0OOooo........oooOO0OOooo........oo
108
109 PrimaryGeneratorAction::~PrimaryGeneratorActio
110 delete fParticleGun;
111 #ifdef WITHROOT
112 delete fInputFile;
113 #endif
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117
118 void PrimaryGeneratorAction::GeneratePrimaries
119 #ifdef WITHROOT
120 if (fReadInputFile) {
121 // get primary particles from input file
122 if (fInputFile == nullptr) {
123 OpenInput();
124 }
125 Float_t currentEventId = -1;
126 while (fHgcalReader->Next()) {
127 // First particle - set event ID
128 if (currentEventId == -1) {
129 currentEventId = **fHgcalEventId;
130 fEventCounter++;
131 // check if event should be ignored (i
132 if (fEventCounter < fStartFromEvent)
133 continue;
134 } else {
135 // check if event is to be skipped
136 if (fEventCounter < fStartFromEvent) {
137 if (currentEventId != **fHgcalEventI
138 currentEventId = -1;
139 continue;
140 } else {
141 // check if current particle belongs
142 if (currentEventId != **fHgcalEventI
143 break;
144 }
145 }
146 auto vertex = new G4PrimaryVertex(
147 G4ThreeVector(**fHgcalPosX, **fHgcalPo
148 auto particleDefinition =
149 G4ParticleTable::GetParticleTable()-
150 auto particle = new G4PrimaryParticle(pa
151 G4ThreeVector momentum(**fHgcalMomX, **f
152 particle->SetMomentumDirection(momentum.
153 particle->SetKineticEnergy(momentum.mag(
154 vertex->SetPrimary(particle);
155 anEvent->AddPrimaryVertex(vertex);
156 }
157 // if no particles are left - terminate ru
158 if (fEventCounter >= fStartFromEvent &&
159 anEvent->GetNumberOfPrimaryVertex() ==
160 G4ExceptionDescription msg;
161 msg << "Input file does not contain any
162 << anEvent->GetEventID() << ").\n";
163 msg << "Run will be aborted with " << an
164 << " events processed.\n";
165 G4Exception("PrimaryGeneratorAction::Gen
166 "EndOfInputFile", JustWarnin
167 G4RunManager::GetRunManager()->AbortRun(
168 }
169 // warn if this is the last event and some
170 // unprocessed
171 if (G4RunManager::GetRunManager()
172 ->GetCurrentRun()
173 ->GetNumberOfEventToBeProcessed()
174 G4ExceptionDescription msg;
175 msg << "Input file contains more events
176 G4Exception("PrimaryGeneratorAction::Gen
177 "UnusedEventsInInputFile", J
178 }
179 } else
180 #endif
181 {
182 // Beam position: z
183 // Use value set by UI command, or set it
184 if (fBeamZ0 == -999 * CLHEP::m) {
185 G4double worldDZ = 0;
186 if (!fEnvelopeBox) {
187 G4LogicalVolume *envLV =
188 G4LogicalVolumeStore::GetInstance(
189 if (envLV)
190 fEnvelopeBox = dynamic_cast<G4Box *>
191 }
192 if (fEnvelopeBox) {
193 worldDZ = fEnvelopeBox->GetZHalfLength
194 fBeamZ0 = -worldDZ;
195 } else {
196 G4ExceptionDescription msg;
197 msg << "World volume of box shape not
198 msg << "Perhaps you have changed geome
199 msg << "The gun will be place at the c
200 G4Exception("PrimaryGeneratorAction::G
201 "WorldVolume", JustWarning
202 }
203 }
204 // Beam position: x,y
205 switch (fBeamType) {
206 case eNone:
207 default:
208 fParticleGun->SetParticlePosition(G4Thre
209 break;
210 case eGaussian:
211 fParticleGun->SetParticlePosition(
212 G4ThreeVector(G4RandGauss::shoot(0.,
213 G4RandGauss::shoot(0.,
214 break;
215 case eFlat:
216 fParticleGun->SetParticlePosition(
217 G4ThreeVector(G4RandFlat::shoot(-fSi
218 G4RandFlat::shoot(-fSi
219 break;
220 }
221 // Particle momentum
222 if (fMomentumGaussianSpread > 0) {
223 double energy = fParticleGun->GetParticl
224 energy += G4RandGauss::shoot(0., fMoment
225 fParticleGun->SetParticleEnergy(energy);
226 }
227 fParticleGun->GeneratePrimaryVertex(anEven
228 }
229 }
230 //....oooOO0OOooo........oooOO0OOooo........oo
231