Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: trunk/source/processes/hadronic/models/binary_cascade/test/G4XResonanceTest.cc@ 1228

Last change on this file since 1228 was 1199, checked in by garnier, 16 years ago
nvx fichiers dans CVS
File size: 10.3 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	//
24	// -------------------------------------------------------------------
25	// GEANT 4 class file --- Copyright CERN 1998
26	// CERN Geneva Switzerland
27	//
28	//
29	// File name: G4XResonanceTest.cc
30	//
31	// Author: Maria Grazia Pia (pia@genova.infn.it),
32	//
33	// Creation date: 15 June 1999
34	//
35	// Modifications:
36	//
37	// -------------------------------------------------------------------
38
39	#include "globals.hh"
40
41	#include "G4ios.hh"
42	#include <fstream>
43	#include <iomanip>
44	#include <iostream>
45	#include <assert.h>
46
47	#include "CLHEP/Hist/TupleManager.h"
48	#include "CLHEP/Hist/HBookFile.h"
49	#include "CLHEP/Hist/Histogram.h"
50	#include "CLHEP/Hist/Tuple.h"
51
52	#include "Randomize.hh"
53
54	#include "G4ThreeVector.hh"
55	#include "G4LorentzVector.hh"
56	#include "G4LorentzRotation.hh"
57
58	#include "G4AntiProton.hh"
59	#include "G4AntiNeutron.hh"
60	#include "G4Proton.hh"
61	#include "G4Neutron.hh"
62	#include "G4PionPlus.hh"
63	#include "G4PionMinus.hh"
64	#include "G4PionZero.hh"
65	#include "G4Gamma.hh"
66	#include "G4MuonMinus.hh"
67	#include "G4MuonPlus.hh"
68	#include "G4KaonMinus.hh"
69	#include "G4KaonPlus.hh"
70	#include "G4NeutrinoMu.hh"
71	#include "G4AntiNeutrinoMu.hh"
72	#include "G4Lambda.hh"
73
74	#include "G4ParticleDefinition.hh"
75
76	#include "G4VDecayChannel.hh"
77	#include "G4DecayTable.hh"
78
79	#include "G4KineticTrack.hh"
80	#include "G4KineticTrackVector.hh"
81
82	#include "G4ParticleTable.hh"
83	#include "G4VShortLivedParticle.hh"
84	#include "G4ShortLivedConstructor.hh"
85	#include "G4ShortLivedTable.hh"
86
87	#include "G4VCrossSectionSource.hh"
88	#include "G4XAqmTotal.hh"
89	#include "G4XAqmElastic.hh"
90	#include "G4XNNTotalLowE.hh"
91	#include "G4XPDGTotal.hh"
92	#include "G4XNNTotal.hh"
93	#include "G4XPDGElastic.hh"
94
95	#include "G4VXResonanceTable.hh"
96
97	#include "G4XResonance.hh"
98	#include "G4XNNstarTable.hh"
99	#include "G4XDeltaNstarTable.hh"
100	#include "G4XDeltaDeltaTable.hh"
101	#include "G4XNDeltaTable.hh"
102	#include "G4XNDeltastarTable.hh"
103	#include "G4XDeltaDeltastarTable.hh"
104
105
106	int main()
107	{
108	// MGP ---- HBOOK initialization
109	HepTupleManager* hbookManager;
110	hbookManager = new HBookFile("xresonance.hbook", 58);
111	assert (hbookManager != 0);
112
113	// MGP ---- Book histograms
114
115	HepHistogram* hSigma;
116	hSigma = hbookManager->histogram("CrossSection", 100,0.,10.);
117	assert (hSigma != 0);
118
119	// MGP ---- Book a ntuple
120	HepTuple* ntuple;
121	ntuple = hbookManager->ntuple("CrossSection ntuple");
122	assert (ntuple != 0);
123
124	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
125
126	G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
127	G4ParticleDefinition* antiProton = G4AntiProton::AntiProtonDefinition();
128	G4ParticleDefinition* neutron = G4Neutron::NeutronDefinition();
129	G4ParticleDefinition* antiNeutron = G4AntiNeutron::AntiNeutronDefinition();
130
131	G4ParticleDefinition* pionPlus = G4PionPlus::PionPlusDefinition();
132	G4ParticleDefinition* pionMinus = G4PionMinus::PionMinusDefinition();
133	G4ParticleDefinition* pionZero = G4PionZero::PionZeroDefinition();
134
135	G4ParticleDefinition* kaonPlus = G4KaonPlus::KaonPlusDefinition();
136	G4ParticleDefinition* kaonMinus = G4KaonMinus::KaonMinusDefinition();
137
138	G4ParticleDefinition* lambda = G4Lambda::LambdaDefinition();
139
140	G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
141	G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
142
143	G4ParticleDefinition* theNeutrinoMu = G4NeutrinoMu::NeutrinoMuDefinition();
144	G4ParticleDefinition* theAntiNeutrinoMu = G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
145
146	G4ShortLivedConstructor ShortLived;
147	ShortLived.ConstructParticle();
148
149	// Select track 1
150
151	G4cout << "---- Set KineticTrack 1 ----" << G4endl;
152	G4int id1 = 0;
153	G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- " << G4endl;
154	G4cin >> id1;
155
156	G4ParticleDefinition* definition1;
157	switch (id1)
158	{
159	case 1:
160	definition1 = proton;
161	break;
162	case 2:
163	definition1 = antiProton;
164	break;
165	case 3:
166	definition1 = neutron;
167	break;
168	case 4:
169	definition1 = antiNeutron;
170	break;
171	case 5:
172	definition1 = pionPlus;
173	break;
174	case 6:
175	definition1 = pionMinus;
176	break;
177	case 7:
178	definition1 = kaonPlus;
179	break;
180	case 8:
181	definition1 = kaonMinus;
182	break;
183	default:
184	definition1 = proton;
185	break;
186	}
187	G4double mass1 = definition1->GetPDGMass();
188
189	// Formation time
190	G4double time1 = 0.;
191	// G4cout << "Formation time (ns)" << G4endl;
192	// G4cin >> time1;
193	// if (time1 < 0.) time1 = 0.;
194	// time1 = time1 * ns;
195
196	// Initial position
197	G4double x1 = 0.;
198	G4double y1 = 0.;
199	G4double z1 = 0.;
200	// G4cout <<"Initial position (cm) " << G4endl;
201	// G4cin >> x1 >> y1 >> z1;
202	G4ThreeVector position1(x1cm, y1cm, z1*cm);
203
204	// Momentum
205	G4double pzMin1;
206	G4double pzMax1;
207	G4cout << "Pz min and max (GeV)" << G4endl;
208	G4cin >> pzMin1 >> pzMax1;
209	pzMin1 = pzMin1 * GeV;
210	pzMax1 = pzMax1 * GeV;
211
212	// Select track 2
213
214	G4cout << "---- Set KineticTrack 2 ----" << G4endl;
215	G4int id2 = 0;
216	G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- " << G4endl;
217	G4cin >> id2;
218
219	G4ParticleDefinition* definition2;
220	G4KineticTrack trk2;
221	switch (id2)
222	{
223	case 1:
224	definition2 = proton;
225	break;
226	case 2:
227	definition2 = antiProton;
228	break;
229	case 3:
230	definition2 = neutron;
231	break;
232	case 4:
233	definition2 = antiNeutron;
234	break;
235	case 5:
236	definition2 = pionPlus;
237	break;
238	case 6:
239	definition2 = pionMinus;
240	break;
241	case 7:
242	definition2 = kaonPlus;
243	break;
244	case 8:
245	definition2 = kaonMinus;
246	break;
247	default:
248	definition2 = proton;
249	break;
250	}
251	G4double mass2 = definition2->GetPDGMass();
252
253	// Formation time
254	G4double time2 = 0.;
255	// G4cout << "Formation time (ns)" << G4endl;
256	// G4cin >> time2;
257	// if (time2 < 0.) time2 = 0.;
258	// time2 = time2 * ns;
259
260	// Initial position
261	G4double x2 = 0.;
262	G4double y2 = 0.;
263	G4double z2 = 0.;
264	// G4cout <<"Initial position (cm) " << G4endl;
265	// G4cin >> x2 >> y2 >> z2;
266	G4ThreeVector position2(x2cm, y2cm, z2*cm);
267
268	// Momentum
269	G4double pzMin2;
270	G4double pzMax2;
271	G4cout << "Pz min and max (GeV)" << G4endl;
272	G4cin >> pzMin2 >> pzMax2;
273	pzMin2 = pzMin2 * GeV;
274	pzMax2 = pzMax2 * GeV;
275
276	// N1440
277	G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
278	G4ParticleDefinition* N1440Zero = particleTable->FindParticle("N(1440)0");
279	G4cout << N1440Zero->GetParticleName() << " created, type is "
280	<< N1440Zero->GetParticleType() << G4endl;
281	G4ParticleDefinition* N1440Plus = particleTable->FindParticle("N(1440)+");
282	G4cout << N1440Plus->GetParticleName() << " created, type is "
283	<< N1440Plus->GetParticleType() << G4endl;
284	G4ParticleDefinition* deltaPlus = particleTable->FindParticle("delta+");
285	G4cout << deltaPlus->GetParticleName() << " created, type is "
286	<< deltaPlus->GetParticleType() << G4endl;
287
288	G4String name = N1440Plus->GetParticleName();
289
290	G4VXResonanceTable* table = new G4XDeltaNstarTable;
291
292	// Create the ResonanceCrossSection
293
294	G4VCrossSectionSource* source = new G4XResonance(definition1,definition2,
295	deltaPlus,N1440Plus,*table);
296	G4cout << "*** " << source->Name() << "*** created " << G4endl;
297	source->Print();
298
299	const G4PhysicsVector* sigmaVector = table->CrossSectionTable(N1440Plus->GetParticleName());
300
301	// Number of iterations
302	G4int n = 0;
303	G4cout << "---- Number of iterations ----" << G4endl;
304	G4cin >> n;
305	if (n < 1) n = 1;
306
307	G4double step1 = (pzMax1 - pzMin1) / n;
308	G4double step2 = (pzMax2 - pzMin2) / n;
309
310	G4cout << "pz min max step "
311	<< pzMin1 / GeV << ", "
312	<< pzMax1 / GeV << ", "
313	<< step1 / GeV << ", " << G4endl
314	<< "pz min max step "
315	<< pzMin2 / GeV << ", "
316	<< pzMax2 / GeV << ", "
317	<< step2 / GeV << G4endl;
318
319	G4int i;
320	for (i=0; i<n; i++)
321	{
322	G4double pStep1 = pzMin1 + step1 * i;
323	G4ThreeVector p1(0.,0.,pStep1);
324	G4double e1 = sqrt(mass1mass1 + p1.mag()p1.mag());
325	G4LorentzVector p41(p1,e1);
326	G4KineticTrack trk1(definition1,time1,position1,p41);
327
328	G4double pStep2 = -(pzMin2 + step2 * i);
329	G4ThreeVector p2(0.,0.,pStep2);
330	G4double e2 = sqrt(mass2mass2 + p2.mag()p2.mag());
331	G4LorentzVector p42(p2,e2);
332	G4KineticTrack trk2(definition2,time2,position2,p42);
333
334	G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
335
336
337	// source->PrintAll(trk1,trk2);
338	G4double sigma = source->CrossSection(trk1,trk2);
339
340	G4bool dummy = false;
341	G4double ppSigma = ((G4PhysicsVector*)sigmaVector)->GetValue(sqrtS,dummy);
342
343	G4cout << definition1->GetParticleName() << " " << definition2->GetParticleName()
344	<< " cross section = " << sigma / millibarn << G4endl;
345	G4cout << "pp cross section = " << ppSigma / millibarn << G4endl;
346
347	// ntuple->column("ecm",sqrtS/GeV);
348	// ntuple->column("sigma",sigma/millibarn);
349	// ntuple->dumpData();
350	}
351
352
353	hbookManager->write();
354
355	delete source;
356
357	delete table;
358
359	return EXIT_SUCCESS;
360	}
361
362

Note: See TracBrowser for help on using the repository browser.

Download in other formats: