source: trunk/examples/extended/geometry/olap/src/OlapGenerator.cc@ 1199

Last change on this file since 1199 was 807, checked in by garnier, 17 years ago

update

File size: 6.3 KB
Line 
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//
27// $Id: OlapGenerator.cc,v 1.4 2006/06/29 17:22:56 gunter Exp $
28// GEANT4 tag $Name: $
29//
30//
31// --------------------------------------------------------------
32// OlapGenerator
33//
34// Author: Martin Liendl - Martin.Liendl@cern.ch
35//
36// --------------------------------------------------------------
37//
38#include "OlapGenerator.hh"
39#include "G4Event.hh"
40#include "G4VisExtent.hh"
41#include "G4PrimaryParticle.hh"
42#include "G4Geantino.hh"
43#include "G4ChargedGeantino.hh"
44#include "G4RunManager.hh"
45#include "G4Run.hh"
46#include "G4Event.hh"
47
48//#define OLAPDEBUG1
49
50OlapGrid::OlapGrid()
51 : count(3,0), gsize(3,3), axis(0), eventsPerRun(27)
52{
53 Reset();
54 G4cout << "OlapGrid(): " << G4endl
55 << " gsize: " << gsize[0] << " " << gsize[1]
56 << " " << gsize[2] << G4endl;
57}
58
59void OlapGrid::Next()
60{
61 G4int out = (axis+1)%3; // outer loop
62 G4int in = (axis+2)%3; // inner loop
63 if ( !(count[in] = (count[in]+1)%gsize[in]) )
64 {
65 if( !(count[out] = (count[out]+1)%gsize[out]) )
66 {
67 axis = (axis+1)%3;
68 }
69 }
70}
71
72OlapGenerator::OlapGenerator() : autoinc(false)
73{
74 if (!ext.size())
75 {
76 ext.push_back(1);
77 ext.push_back(2);
78 ext.push_back(3);
79 }
80 /*
81 ext[0]=1.;
82 ext[1]=2.;
83 ext[2]=3.;
84 */
85}
86
87
88OlapGenerator::~OlapGenerator()
89{
90}
91
92
93void OlapGenerator::GeneratePrimaries(G4Event * anEvent)
94{
95 G4int out = (grid.axis+1)%3; // outer loop
96 G4int in = (grid.axis+2)%3; // inner loop
97
98 posAB[grid.axis]=-ext[grid.axis]/2.;
99 posAB[in]=-ext[in]/2.+ext[in]/G4double(grid.gsize[in]-1)
100 *G4double(grid.count[in]);
101 posAB[out]=-ext[out]/2.+ext[out]/G4double(grid.gsize[out]-1)
102 *G4double(grid.count[out]);
103
104 posBA = posAB;
105 posBA[grid.axis]=-posBA[grid.axis];
106
107 G4ThreeVector dirAB, dirBA;
108 dirAB[grid.axis] = 1;
109 dirBA[grid.axis] = -1;
110
111 #ifdef OLAPDEBUG1
112 G4int runID = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();
113 G4int evtID = anEvent->GetEventID();
114 G4cout << "generator: "
115 << "run=" << runID << " evt=" << evtID
116 << " axis=" << grid.axis << " out=" << grid.count[out]
117 << " in=" << grid.count[in]
118 << " posAB=" << posAB
119 << " posBA=" << posBA << G4endl;
120 #endif
121
122 // now generator 2 geantinos flying in opposite direction from A->B and B->A
123// ==========================
124 G4ParticleDefinition* pd = G4Geantino::GeantinoDefinition();
125
126 // create 2 opposite positioned vertices
127 G4PrimaryVertex* vertexAB =
128 new G4PrimaryVertex(posAB,0);
129 G4PrimaryVertex* vertexBA =
130 new G4PrimaryVertex(posBA,0);
131
132 // create one geantino per vertex as primary
133 G4double mass = pd->GetPDGMass();
134 G4double energy = 1. + mass;
135 G4double pmom = std::sqrt(energy*energy-mass*mass);
136 G4double pxAB = pmom*dirAB.x();
137 G4double pyAB = pmom*dirAB.y();
138 G4double pzAB = pmom*dirAB.z();
139 G4double pxBA = pmom*dirBA.x();
140 G4double pyBA = pmom*dirBA.y();
141 G4double pzBA = pmom*dirBA.z();
142
143 G4PrimaryParticle* particleAB = new G4PrimaryParticle(pd,pxAB,pyAB,pzAB);
144 G4PrimaryParticle* particleBA = new G4PrimaryParticle(pd,pxBA,pyBA,pzBA);
145 particleAB->SetMass( mass );
146 particleBA->SetMass( mass );
147 particleAB->SetCharge(0.);
148 particleBA->SetCharge(0.);
149 vertexAB->SetPrimary( particleAB );
150 vertexBA->SetPrimary( particleBA );
151 anEvent->AddPrimaryVertex( vertexAB );
152 anEvent->AddPrimaryVertex( vertexBA );
153
154 // ok, set all counters to be ready for the next position
155 if (autoinc)
156 grid.Next();
157 /*
158 if ( !(count[in] = (count[in]+1)%grid[in]) ) {
159 if( !(count[out] = (count[out]+1)%grid[out]) ) {
160 axis = (axis+1)%3;
161 }
162 }
163 */
164
165}
166
167
168void OlapGenerator::SetExtent(const G4VisExtent & anExt)
169{
170 ext[0]=(anExt.GetXmax()-anExt.GetXmin());
171 ext[1]=(anExt.GetYmax()-anExt.GetYmin());
172 ext[2]=(anExt.GetZmax()-anExt.GetZmin());
173 for (G4int i=0; i<3; i++)
174 ext[i] += ext[i]/2000.; // enlarge by 1%
175
176 Reset();
177
178}
179
180void OlapGenerator::SetExtent(G4double e)
181{
182 ext[0]=e;
183 ext[1]=e;
184 ext[2]=e;
185}
186
187
188void OlapGenerator::Reset()
189{
190 grid.Reset();
191 /*
192 for (G4int i=0; i<3; i++) {
193 count[i]=0;
194 }
195
196 axis=0;
197 */
198}
199
200
201void OlapGenerator::SetGrid(G4int x, G4int y, G4int z)
202{
203 if (x>2)
204 grid.gsize[0]=x;
205 else
206 G4cerr << "Warning in OlapGenerator::SetGrid(): wrong x-grid parameter: " << x << G4endl;
207
208
209 if (y>2)
210 grid.gsize[1]=y;
211 else
212 G4cerr << "Warning in OlapGenerator::SetGrid(): wrong y-grid parameter: " << y << G4endl;
213
214 if (z>2)
215 grid.gsize[2]=z;
216 else
217 G4cerr << "Warning in OlapGenerator::SetGrid(): wrong z-grid parameter: " << z << G4endl;
218
219 grid.eventsPerRun = x*y + y*z + x*z;
220 Reset();
221}
Note: See TracBrowser for help on using the repository browser.