| Geant4 Cross Reference |
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 #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........oooOO0OOooo........oooOO0OOooo......
46
47 PrimaryGeneratorAction::PrimaryGeneratorAction()
48 : G4VUserPrimaryGeneratorAction(), fParticleGun(nullptr),
49 fEnvelopeBox(nullptr) {
50 G4int n_particle = 1;
51 fParticleGun = new G4ParticleGun(n_particle);
52
53 // default particle kinematic
54 G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
55 G4String particleName;
56 G4ParticleDefinition *particle =
57 particleTable->FindParticle(particleName = "e+");
58 fParticleGun->SetParticleDefinition(particle);
59 fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
60 fParticleGun->SetParticleEnergy(30. * CLHEP::GeV);
61 fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
62
63 fMessenger = new PrimaryGeneratorMessenger(this);
64 #ifdef WITHROOT
65 if (fReadInputFile)
66 OpenInput();
67 #endif
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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->IsZombie()) {
80 G4ExceptionDescription msg;
81 msg << "Input file '" << fPathInputFile << "' cannot be opened.\n";
82 G4Exception("PrimaryGeneratorAction::GeneratePrimaries()", "WrongInputName",
83 FatalErrorInArgument, msg);
84 }
85 auto fileDirectory =
86 dynamic_cast<TDirectory *>(fInputFile->Get("VirtualDetector"));
87 if (fileDirectory == nullptr) {
88 G4ExceptionDescription msg;
89 msg << "Input file '" << fPathInputFile
90 << "' does not contain 'VirtualDetector' directory.\n";
91 G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
92 "WrongInputDirectory", FatalErrorInArgument, msg);
93 }
94 fHgcalReader = new TTreeReader("HGCAL", fileDirectory);
95 fHgcalEventId = new TTreeReaderValue<Float_t>(*fHgcalReader, "EventID");
96 fHgcalPdgId = new TTreeReaderValue<Float_t>(*fHgcalReader, "PDGid");
97 fHgcalPosX = new TTreeReaderValue<Float_t>(*fHgcalReader, "x");
98 fHgcalPosY = new TTreeReaderValue<Float_t>(*fHgcalReader, "y");
99 fHgcalMomX = new TTreeReaderValue<Float_t>(*fHgcalReader, "Px");
100 fHgcalMomY = new TTreeReaderValue<Float_t>(*fHgcalReader, "Py");
101 fHgcalMomZ = new TTreeReaderValue<Float_t>(*fHgcalReader, "Pz");
102
103 return;
104 }
105 #endif
106
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108
109 PrimaryGeneratorAction::~PrimaryGeneratorAction() {
110 delete fParticleGun;
111 #ifdef WITHROOT
112 delete fInputFile;
113 #endif
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117
118 void PrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent) {
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 (is before the start position)
132 if (fEventCounter < fStartFromEvent)
133 continue;
134 } else {
135 // check if event is to be skipped
136 if (fEventCounter < fStartFromEvent) {
137 if (currentEventId != **fHgcalEventId)
138 currentEventId = -1;
139 continue;
140 } else {
141 // check if current particle belongs to the next event
142 if (currentEventId != **fHgcalEventId)
143 break;
144 }
145 }
146 auto vertex = new G4PrimaryVertex(
147 G4ThreeVector(**fHgcalPosX, **fHgcalPosY, 32.5 * CLHEP::mm), 0 * CLHEP::s);
148 auto particleDefinition =
149 G4ParticleTable::GetParticleTable()->FindParticle(**fHgcalPdgId);
150 auto particle = new G4PrimaryParticle(particleDefinition);
151 G4ThreeVector momentum(**fHgcalMomX, **fHgcalMomY, **fHgcalMomZ);
152 particle->SetMomentumDirection(momentum.unit());
153 particle->SetKineticEnergy(momentum.mag());
154 vertex->SetPrimary(particle);
155 anEvent->AddPrimaryVertex(vertex);
156 }
157 // if no particles are left - terminate run
158 if (fEventCounter >= fStartFromEvent &&
159 anEvent->GetNumberOfPrimaryVertex() == 0) {
160 G4ExceptionDescription msg;
161 msg << "Input file does not contain any more events (current event ID = "
162 << anEvent->GetEventID() << ").\n";
163 msg << "Run will be aborted with " << anEvent->GetEventID()
164 << " events processed.\n";
165 G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
166 "EndOfInputFile", JustWarning, msg);
167 G4RunManager::GetRunManager()->AbortRun();
168 }
169 // warn if this is the last event and some events are left in the file
170 // unprocessed
171 if (G4RunManager::GetRunManager()
172 ->GetCurrentRun()
173 ->GetNumberOfEventToBeProcessed() == anEvent->GetEventID() + 1) {
174 G4ExceptionDescription msg;
175 msg << "Input file contains more events than requested in the run.\n";
176 G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
177 "UnusedEventsInInputFile", JustWarning, msg);
178 }
179 } else
180 #endif
181 {
182 // Beam position: z
183 // Use value set by UI command, or set it to edge of the world volume
184 if (fBeamZ0 == -999 * CLHEP::m) {
185 G4double worldDZ = 0;
186 if (!fEnvelopeBox) {
187 G4LogicalVolume *envLV =
188 G4LogicalVolumeStore::GetInstance()->GetVolume("World");
189 if (envLV)
190 fEnvelopeBox = dynamic_cast<G4Box *>(envLV->GetSolid());
191 }
192 if (fEnvelopeBox) {
193 worldDZ = fEnvelopeBox->GetZHalfLength();
194 fBeamZ0 = -worldDZ;
195 } else {
196 G4ExceptionDescription msg;
197 msg << "World volume of box shape not found.\n";
198 msg << "Perhaps you have changed geometry.\n";
199 msg << "The gun will be place at the center.";
200 G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
201 "WorldVolume", JustWarning, msg);
202 }
203 }
204 // Beam position: x,y
205 switch (fBeamType) {
206 case eNone:
207 default:
208 fParticleGun->SetParticlePosition(G4ThreeVector(0, 0, fBeamZ0));
209 break;
210 case eGaussian:
211 fParticleGun->SetParticlePosition(
212 G4ThreeVector(G4RandGauss::shoot(0., fSigmaBeamX),
213 G4RandGauss::shoot(0., fSigmaBeamY), fBeamZ0));
214 break;
215 case eFlat:
216 fParticleGun->SetParticlePosition(
217 G4ThreeVector(G4RandFlat::shoot(-fSigmaBeamX, fSigmaBeamX),
218 G4RandFlat::shoot(-fSigmaBeamY, fSigmaBeamY), fBeamZ0));
219 break;
220 }
221 // Particle momentum
222 if (fMomentumGaussianSpread > 0) {
223 double energy = fParticleGun->GetParticleEnergy();
224 energy += G4RandGauss::shoot(0., fMomentumGaussianSpread) * energy;
225 fParticleGun->SetParticleEnergy(energy);
226 }
227 fParticleGun->GeneratePrimaryVertex(anEvent);
228 }
229 }
230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
231