source: trunk/source/processes/hadronic/models/binary_cascade/test/G4BinaryCascadeTest.cc@ 1200

Last change on this file since 1200 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 11.9 KB
Line 
1//
2// ********************************************************************
3// * DISCLAIMER *
4// * *
5// * The following disclaimer summarizes all the specific disclaimers *
6// * of contributors to this software. The specific disclaimers,which *
7// * govern, are listed with their locations in: *
8// * http://cern.ch/geant4/license *
9// * *
10// * Neither the authors of this software system, nor their employing *
11// * institutes,nor the agencies providing financial support for this *
12// * work make any representation or warranty, express or implied, *
13// * regarding this software system or assume any liability for its *
14// * use. *
15// * *
16// * This code implementation is the intellectual property of the *
17// * GEANT4 collaboration. *
18// * By copying, distributing or modifying the Program (or any work *
19// * based on the Program) you indicate your acceptance of this *
20// * statement, and all its terms. *
21// ********************************************************************
22//
23#include "G4BinaryCascade.hh"
24#include "G4Nucleus.hh"
25#include "G4Proton.hh"
26#include "G4PionPlus.hh"
27#include "G4ThreeVector.hh"
28#include "G4LorentzVector.hh"
29#include "G4DynamicParticle.hh"
30#include "G4StateManager.hh"
31
32#include "G4LeptonConstructor.hh"
33#include "G4Gamma.hh"
34
35#include "G4ExcitationHandler.hh"
36#include "G4PreCompoundModel.hh"
37
38#define APPLYYOURSELF 1
39
40/*
41#include "ResourceMgr.hh"
42// in G4BinaryCascade.cc
43ResourceMgr absorbRes;
44ResourceMgr buildTargetListRes;
45ResourceMgr findCollisionRes;
46ResourceMgr findFirstCollisionRes;
47ResourceMgr totalRes;
48ResourceMgr thePropagatorRes;
49ResourceMgr updateRes;
50ResourceMgr scattererRes;
51void PrintResourceUsage(ResourceMgr & resMgr);
52*/
53G4V3DNucleus * global3DNucleus;
54
55
56void Init(G4int code, std::vector<G4ParticleDefinition *> & particles);
57void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax,
58 G4int &nEvents, G4int & particleCode);
59G4ThreeVector GetSpherePoint(G4double r);
60G4ThreeVector GetSurfacePoint(G4double r);
61
62int main(int argc, char * argv[])
63{
64 G4StateManager* g4State=G4StateManager::GetStateManager();
65 if (! g4State->SetNewState(G4State_Init))
66 G4cout << "error changing G4state"<< G4endl;;
67
68 G4String a=G4Proton::ProtonDefinition()->GetParticleName();
69 G4cout << " At Start " << a << G4endl;
70 G4cout.setf( std::ios::scientific, std::ios::floatfield );
71// totalRes.Start();
72 G4LeptonConstructor leptons;
73 leptons.ConstructParticle();
74 G4IonConstructor ions;
75 ions.ConstructParticle();
76 const G4ParticleDefinition* proton = G4Proton::Proton();
77 const G4ParticleDefinition* neutron = G4Neutron::Neutron();
78 const G4ParticleDefinition* pin = G4PionMinus::PionMinus();
79 const G4ParticleDefinition* pip = G4PionPlus::PionPlus();
80 const G4ParticleDefinition* pi0 = G4PionZero::PionZero();
81 const G4ParticleDefinition* deu = G4Deuteron::DeuteronDefinition();
82 const G4ParticleDefinition* tri = G4Triton::TritonDefinition();
83 const G4ParticleDefinition* alp = G4Alpha::AlphaDefinition();
84 const G4ParticleDefinition* ion = G4GenericIon::GenericIon();
85 G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
86 partTable->SetReadiness();
87
88 if(!G4StateManager::GetStateManager()->SetNewState(G4State_Idle))
89 G4cout << "G4StateManager PROBLEM! " << G4endl;
90
91 G4int A, Z, nEvents, particleCode;
92 G4double pMin, pMax;
93 std::vector<G4ParticleDefinition *> particles;
94 if(argc == 7)
95 {
96 A = atoi(argv[1]);
97 Z = atoi(argv[2]);
98 pMin = G4double(atoi(argv[3]))*MeV;
99 pMax = G4double(atoi(argv[4]))*MeV;
100 nEvents = atoi(argv[5]);
101 particleCode = atoi(argv[6]);
102 }
103 else
104 GetParam(A, Z, pMin, pMax, nEvents, particleCode);
105
106 if(pMin > pMax)
107 {
108 G4double pTmp = pMin;
109 pMin = pMax;
110 pMax = pTmp;
111 }
112
113 Init(particleCode, particles);
114
115 G4BinaryCascade bic;
116 G4ExcitationHandler * excite = new G4ExcitationHandler();
117 G4VPreCompoundModel * precompound = new G4PreCompoundModel(excite);
118
119// bic.SetDeExcitation(0);
120 bic.SetDeExcitation(precompound);
121
122
123 for(G4int event = 0; event < nEvents; ++event)
124 {
125// Get particle definition
126 G4cout << " ---- event " << event << "----------------------"<< G4endl;
127 G4double type = G4UniformRand();
128 unsigned int index = G4int(type*particles.size());
129 if(index == particles.size()) // if random gave 1!
130 --index;
131 G4ParticleDefinition * particle = particles[index];
132// position
133 G4ThreeVector pos(0, 0, -5.*fermi);
134// momentum
135 G4double p = G4UniformRand()*(pMax-pMin)+pMin;
136 G4double mass = particle->GetPDGMass();
137 G4double e = sqrt(mass*mass+p*p);
138 G4LorentzVector mom(0, p, 0, e);
139
140 G4cout << "initial p /E . Ekin, nucleus " << p << " / " << e
141 << " " << e-mass << " "
142 <<G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z, A)<< G4endl;
143
144// build the nucleus
145#ifdef APPLYYOURSELF
146 // G4Nucleus nucleus(G4double(A), G4double(Z)); // for ApplyYourself()
147 G4Nucleus nucleus(A, Z); // for ApplyYourself()
148// global3DNucleus = new G4Fancy3DNucleus;
149// global3DNucleus->Init(A, Z);
150#else
151 G4V3DNucleus * fancyNucleus = new G4Fancy3DNucleus; // for PropagateSTL()
152 fancyNucleus->Init(A, Z);
153#endif
154
155// build the projectile
156#ifdef APPLYYOURSELF
157 G4DynamicParticle * dynamic = new G4DynamicParticle(particle, mom);
158 G4double time = 0.;
159 G4Track projectile(dynamic, time, pos);
160 // needed by G4ParticleChange::Initialize(G4Track &)
161 // G4Step * step = new G4Step();
162 //step->SetStepLength(1.0*mm);
163 // projectile.SetStep(step);
164 G4HadFinalState * aFinalState = bic.ApplyYourself(projectile, nucleus);
165
166 G4double QValue = 0, Etot=0;
167 G4LorentzVector pTot = 0;
168// G4Track * second;
169 G4HadSecondary * second;
170 const G4DynamicParticle * aSec;
171
172 for(G4int isec=0;isec<aFinalState->GetNumberOfSecondaries();isec++)
173 {
174 second = aFinalState->GetSecondary(isec);
175 aSec = second->GetParticle();
176 G4cout << "SECONDARIES info ";
177 G4cout << aSec->GetDefinition()->GetBaryonNumber() << " ";
178 G4cout << aSec->GetTotalEnergy();
179 G4cout << aSec->GetMomentum()<< " ";
180 G4cout << (1-isec)*aFinalState->GetNumberOfSecondaries();
181 G4cout << G4endl;
182 QValue += aSec->GetKineticEnergy();
183 pTot += aSec->Get4Momentum();
184 if ( aSec->GetDefinition()->GetBaryonNumber() == 0 )
185 {
186 Etot += aSec->GetTotalEnergy();
187 } else {
188 Etot += aSec->GetKineticEnergy();
189 if (aSec->Get4Momentum().mag() > 0.95*GeV )
190 {
191 Etot += aSec->Get4Momentum().mag() - 939*MeV;
192 }
193
194 }
195
196
197// delete second;
198 }
199 G4cout << "QVALUE = "<<QValue<<G4endl;
200 G4cout << "Ptotal, pseudoE, exciteE = "
201 << pTot << " " << Etot << " " << e-mass-Etot <<G4endl;
202
203 aFinalState->Clear();
204// delete dynamic;
205// delete step;
206
207// change->DumpInfo();
208
209#else
210 G4double radius = fancyNucleus->GetOuterRadius();
211 G4KineticTrackVectorSTL * ktv = new G4KineticTrackVectorSTL;
212// bic.GetSpherePoint(radius, pos); // start from the surface
213 pos = GetSpherePoint(radius); // start already inside, pos.mag()< r
214 pos = GetSurfacePoint(radius*1.01); // start outside, pos.mag()= r*1.01
215// pos = G4ThreeVector(1.*fermi,0.*fermi,-15.*fermi);
216 G4KineticTrack * kt = new G4KineticTrack(particle, 0., pos, mom);
217 ktv->push_back(kt);
218 G4ReactionProductVector *result=bic.PropagateSTL(ktv, fancyNucleus);
219
220// bic.GetExcitationEnergy();
221// std::cout << bic.GetExcitationEnergy()/MeV << endl;
222
223 G4ReactionProductVector::iterator result_iter;
224 G4int isec=0;
225 for(result_iter = result->begin();result_iter != result->end();
226 ++result_iter)
227 {
228 G4ThreeVector pFinal= (*result_iter)->GetMomentum();
229
230 G4cout << " Final particle # "<< ++isec << ": "
231 << pFinal.x() << " "
232 << pFinal.y() << " "
233 << pFinal.z() << " "
234 << (*result_iter)->GetTotalEnergy() << " "
235 << G4endl;
236 }
237#endif
238 }
239/*
240 totalRes.Stop();
241 std::cout << "Absorb() resources:" << endl;
242 PrintResourceUsage(absorbRes);
243 std::cout << "FindFirstCollision() resources:" << endl;
244 PrintResourceUsage(findFirstCollisionRes);
245 std::cout << "FindCollision() resources:" << endl;
246 PrintResourceUsage(findCollisionRes);
247 std::cout << "BuildTargetList() resources:" << endl;
248 PrintResourceUsage(buildTargetListRes);
249 std::cout << "thePropagator() resources:" << endl;
250 PrintResourceUsage(thePropagatorRes);
251 std::cout << "UpdateTrackAndCollisions() resources:" << endl;
252 PrintResourceUsage(updateRes);
253 std::cout << "scatterer resources:" << endl;
254 PrintResourceUsage(scattererRes);
255 std::cout << "main() resources:" << endl;
256 PrintResourceUsage(totalRes);
257*/
258// std::cout << "Exiting" << endl;
259 std::cout.flush();
260
261 return 0;
262}
263
264/*
265
266void PrintResourceUsage(ResourceMgr & resMgr)
267{
268 double userTime = resMgr.GetUserTime();
269 double sysTime = resMgr.GetSysTime();
270 double cpuTime = userTime+sysTime;
271 std::cout << "cpu time: " << cpuTime << " (" << userTime << "+"
272 << sysTime << ")" << endl;
273}
274
275
276*/
277
278void Init(G4int code, std::vector<G4ParticleDefinition *> & particles)
279{
280 if(code == 0)
281 {
282 particles.push_back(G4Gamma::Gamma());
283 return;
284 }
285
286 G4int count = 0;
287 while(code != 0)
288 {
289 if(code & 1)
290 {
291 switch (count)
292 {
293 case 0:
294 particles.push_back(G4Proton::Proton());
295 break;
296 case 1:
297 particles.push_back(G4Neutron::Neutron());
298 break;
299 case 2:
300 particles.push_back(G4AntiProton::AntiProton());
301 break;
302 case 3:
303 particles.push_back(G4PionPlus::PionPlus());
304 break;
305 case 4:
306 particles.push_back(G4PionZero::PionZero());
307 break;
308 case 5:
309 particles.push_back(G4PionMinus::PionMinus());
310 break;
311 case 6:
312 particles.push_back(G4KaonPlus::KaonPlus());
313 break;
314 case 7:
315 particles.push_back(G4KaonZero::KaonZero());
316 break;
317 case 8:
318 particles.push_back(G4KaonMinus::KaonMinus());
319 break;
320 case 9:
321 particles.push_back(G4SigmaPlus::SigmaPlus());
322 break;
323 case 10:
324 particles.push_back(G4SigmaZero::SigmaZero());
325 break;
326 case 11:
327 particles.push_back(G4SigmaMinus::SigmaMinus());
328 break;
329 }
330 }
331 ++count;
332 code = code>>1;
333 }
334}
335
336
337void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax,
338 G4int &nEvents, G4int & particleCode)
339{
340// A and Z
341 std::cout << "Input Nucleus A and Z: ";
342 std::cin >> A >> Z;
343// initial impulse
344 std::cout << "Input min and max of |momentum|: ";
345 std::cin >> pMin >> pMax;
346// number of events
347 std::cout << "Input number of events: ";
348 std::cin >> nEvents;
349// particle code
350 std::cout << "Binary code for projectile type:" << std::endl;
351 std::cout << "1 = proton, 2 = neutron, 4 = antiproton" << std::endl;
352 std::cout << "8 = pi+, 16 = pi0, 32 = pi-" << std::endl;
353 std::cout << "64 = k+, 128 = k0, 256 = k-" << std::endl;
354 std::cout << "512 = s+, 1024 = s0, 2048 = s-" << std::endl;
355 std::cout << " 0 for photon" << std::endl;
356 std::cout << "Input projectile type (binary code): ";
357 std::cin >> particleCode;
358}
359
360
361
362
363G4ThreeVector GetSpherePoint(G4double r)
364{
365// random point INSIDE a sphere
366
367 G4double x = r*(2*G4UniformRand()-1);
368 G4double tmp = sqrt(r*r-x*x);
369 G4double y = tmp*(2*G4UniformRand()-1);
370 tmp = sqrt(r*r-x*x-y*y);
371 G4double z = tmp*(2*G4UniformRand()-1);
372 G4ThreeVector point;
373 point.setX(x);
374 point.setY(y);
375 point.setZ(z);
376 return point;
377}
378G4ThreeVector GetSurfacePoint(G4double r)
379{
380// random point oin surface of a sphere
381
382 G4double x = r*(2*G4UniformRand()-1);
383 G4double tmp = sqrt(r*r-x*x);
384 G4double y = tmp*(2*G4UniformRand()-1);
385 G4double z = -sqrt(r*r-x*x-y*y);
386 G4ThreeVector point;
387 point.setX(x);
388 point.setY(y);
389 point.setZ(z);
390 return point;
391}
Note: See TracBrowser for help on using the repository browser.