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