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