| 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 // 26 //
27 // 27 //
28 28
29 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
31 31
32 #include <fstream> 32 #include <fstream>
33 #include <cstdlib> 33 #include <cstdlib>
34 34
35 #include "FCALPrimaryGeneratorAction.hh" 35 #include "FCALPrimaryGeneratorAction.hh"
36 36
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 #include "G4Event.hh" 38 #include "G4Event.hh"
39 #include "G4ParticleGun.hh" 39 #include "G4ParticleGun.hh"
40 #include "G4ParticleTable.hh" 40 #include "G4ParticleTable.hh"
41 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh"
42 #include "Randomize.hh" 42 #include "Randomize.hh"
43 #include "G4DataVector.hh" 43 #include "G4DataVector.hh"
44 #include "G4AutoLock.hh" 44 #include "G4AutoLock.hh"
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 47
48 48
49 // Migration to MT: there is a single input fi 49 // Migration to MT: there is a single input file that is read by all threads.
50 // The idea is that the events are read by a s 50 // The idea is that the events are read by a single thread and processed
51 // by all threads. Threads ask for the next ID 51 // by all threads. Threads ask for the next ID to be processed. When
52 // events are all processed we start over from 52 // events are all processed we start over from the beginning of the file
53 namespace { 53 namespace {
54 G4bool isFileRead = false; 54 G4bool isFileRead = false;
55 G4Mutex mFileRead = G4MUTEX_INITIALIZER; 55 G4Mutex mFileRead = G4MUTEX_INITIALIZER;
56 //Primary kinematics 56 //Primary kinematics
57 G4DataVector fX; 57 G4DataVector fX;
58 G4DataVector fY; 58 G4DataVector fY;
59 G4DataVector fZ; 59 G4DataVector fZ;
60 G4DataVector fCosX; 60 G4DataVector fCosX;
61 G4DataVector fCosY; 61 G4DataVector fCosY;
62 G4DataVector fCosZ; 62 G4DataVector fCosZ;
63 size_t nextEventId = 0; 63 size_t nextEventId = 0;
64 G4Mutex mNextEventId = G4MUTEX_INITIALIZER 64 G4Mutex mNextEventId = G4MUTEX_INITIALIZER;
65 65
66 size_t GetNextId() { 66 size_t GetNextId() {
67 G4AutoLock l(&mNextEventId); 67 G4AutoLock l(&mNextEventId);
68 if ( nextEventId >= fX.size() ) //file 68 if ( nextEventId >= fX.size() ) //file data are over, restart file
69 { 69 {
70 G4Exception("FCALPrimaryGenera 70 G4Exception("FCALPrimaryGeneratorAction::GeneratePrimaries","lAr002",
71 JustWarning,"Data 71 JustWarning,"Data file with kinematics is over, restart it");
72 nextEventId=0; 72 nextEventId=0;
73 } 73 }
74 return nextEventId++; 74 return nextEventId++;
75 } 75 }
76 76
77 void ReadKinematicFromFile(G4double energy 77 void ReadKinematicFromFile(G4double energy) {
78 //Only one thread shoud read input fil 78 //Only one thread shoud read input file
79 G4AutoLock l(&mFileRead); 79 G4AutoLock l(&mFileRead);
80 if ( isFileRead ) return; 80 if ( isFileRead ) return;
81 // Read Kinematics from file 81 // Read Kinematics from file
82 G4String file_name = "data-tracks/trac 82 G4String file_name = "data-tracks/tracks-80GeV.dat";
83 if (energy < 30*GeV) 83 if (energy < 30*GeV)
84 file_name = "data-tracks/tracks-20 84 file_name = "data-tracks/tracks-20GeV.dat";
85 else if (energy < 50*GeV) 85 else if (energy < 50*GeV)
86 file_name = "data-tracks/tracks-40 86 file_name = "data-tracks/tracks-40GeV.dat";
87 else if (energy < 70*GeV) 87 else if (energy < 70*GeV)
88 file_name = "data-tracks/tracks-60 88 file_name = "data-tracks/tracks-60GeV.dat";
89 else if (energy < 90*GeV) 89 else if (energy < 90*GeV)
90 file_name = "data-tracks/tracks-80 90 file_name = "data-tracks/tracks-80GeV.dat";
91 else if (energy < 150*GeV) 91 else if (energy < 150*GeV)
92 file_name = "data-tracks/tracks-12 92 file_name = "data-tracks/tracks-120GeV.dat";
93 else 93 else
94 file_name = "data-tracks/tracks-20 94 file_name = "data-tracks/tracks-200GeV.dat";
95 std::ifstream Traks_file(file_name); 95 std::ifstream Traks_file(file_name);
96 if(!Traks_file) 96 if(!Traks_file)
97 { 97 {
98 G4ExceptionDescription ed; 98 G4ExceptionDescription ed;
99 ed << "Failed to open file " << fi 99 ed << "Failed to open file " << file_name << G4endl;
100 G4Exception("FCALPrimaryGeneratorA 100 G4Exception("FCALPrimaryGeneratorAction::FCALPrimaryGeneratorAction()",
101 "lAr001",FatalExceptio 101 "lAr001",FatalException,ed);
102 } 102 }
103 G4double xx=0,yy=0,zz=0,c1=0,c2=0,c3=0 103 G4double xx=0,yy=0,zz=0,c1=0,c2=0,c3=0;
104 G4int iev = 0; 104 G4int iev = 0;
105 while(!(Traks_file.eof())) { 105 while(!(Traks_file.eof())) {
106 Traks_file >> iev >> xx >> yy >> z 106 Traks_file >> iev >> xx >> yy >> zz >> c1 >> c2 >> c3;
107 fX.push_back(xx*cm); 107 fX.push_back(xx*cm);
108 fY.push_back(yy*cm); 108 fY.push_back(yy*cm);
109 fZ.push_back(zz*cm); 109 fZ.push_back(zz*cm);
110 fCosX.push_back(c1); 110 fCosX.push_back(c1);
111 fCosY.push_back(c2); 111 fCosY.push_back(c2);
112 fCosZ.push_back(c3); 112 fCosZ.push_back(c3);
113 } 113 }
114 G4cout << "Read " << fX.size() << " ev 114 G4cout << "Read " << fX.size() << " events from file " << file_name << G4endl;
115 isFileRead= true; 115 isFileRead= true;
116 Traks_file.close(); 116 Traks_file.close();
117 return; 117 return;
118 } 118 }
119 } 119 }
120 120
121 //....oooOO0OOooo........oooOO0OOooo........oo 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122 122
123 FCALPrimaryGeneratorAction::FCALPrimaryGenerat 123 FCALPrimaryGeneratorAction::FCALPrimaryGeneratorAction() :
124 fVerbosity(0) 124 fVerbosity(0)
125 { 125 {
126 particleGun = new G4ParticleGun(); 126 particleGun = new G4ParticleGun();
127 127
128 // default Particle 128 // default Particle
129 G4String particleName; 129 G4String particleName;
130 G4ParticleTable* particleTable = G4ParticleT 130 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
131 G4ParticleDefinition* particle = particleTab 131 G4ParticleDefinition* particle = particleTable->FindParticle(particleName="e-");
132 particleGun->SetParticleDefinition(particle) 132 particleGun->SetParticleDefinition(particle);
133 133
134 // default Energy 134 // default Energy
135 particleGun->SetParticleEnergy(20*GeV); 135 particleGun->SetParticleEnergy(20*GeV);
136 } 136 }
137 137
138 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 139
140 FCALPrimaryGeneratorAction::~FCALPrimaryGenera 140 FCALPrimaryGeneratorAction::~FCALPrimaryGeneratorAction()
141 { 141 {
142 delete particleGun; 142 delete particleGun;
143 } 143 }
144 144
145 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 146
147 void FCALPrimaryGeneratorAction::GeneratePrima 147 void FCALPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
148 { 148 {
149 //this function is called at the begining of 149 //this function is called at the begining of event
150 ReadKinematicFromFile(particleGun->GetPartic 150 ReadKinematicFromFile(particleGun->GetParticleEnergy());
151 151
152 //Get next event to be processed 152 //Get next event to be processed
153 size_t nEvent = GetNextId(); 153 size_t nEvent = GetNextId();
154 particleGun->SetParticlePosition(G4ThreeVe 154 particleGun->SetParticlePosition(G4ThreeVector(fX[nEvent],fY[nEvent],fZ[nEvent]));
155 particleGun->SetParticleMomentumDirection(G4 155 particleGun->SetParticleMomentumDirection(G4ThreeVector(-1.0*fCosX[nEvent],
156 fCosY[nEvent], 156 fCosY[nEvent],
157 -1.0*fCosZ[nEvent])); 157 -1.0*fCosZ[nEvent]));
158 158
159 particleGun->GeneratePrimaryVertex(anEvent); 159 particleGun->GeneratePrimaryVertex(anEvent);
160 160
161 if (fVerbosity) 161 if (fVerbosity)
162 { 162 {
163 G4cout<< " Event "<<anEvent->GetEvent 163 G4cout<< " Event "<<anEvent->GetEventID()<< " Generated Vertex : "
164 <<anEvent->GetEventID() <<" (x,y,z 164 <<anEvent->GetEventID() <<" (x,y,z)=(" << fX[nEvent] << ","
165 <<fY[nEvent] << "," << fZ[nEvent]< 165 <<fY[nEvent] << "," << fZ[nEvent]<< ") (cosX,cosY,cosZ)=("
166 << -1.*fCosX[nEvent] << "," << fCo 166 << -1.*fCosX[nEvent] << "," << fCosY[nEvent]
167 <<"," << -1.*fCosZ[nEvent] << ")"< 167 <<"," << -1.*fCosZ[nEvent] << ")"<<G4endl;
168 } 168 }
169 169
170 } 170 }
171 171
172 //....oooOO0OOooo........oooOO0OOooo........oo 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 173
174 174
175 175