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

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

nvx fichiers dans CVS
File size: 8.8 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: G4CollisionTest.cc
30//
31// Author: Maria Grazia Pia (pia@genova.infn.it),
32//
33// Creation date: 15 April 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 "G4Proton.hh"
60#include "G4AntiNeutron.hh"
61#include "G4Neutron.hh"
62#include "G4PionPlus.hh"
63#include "G4PionMinus.hh"
64#include "G4PionZero.hh"
65#include "G4KaonMinus.hh"
66#include "G4KaonPlus.hh"
67#include "G4Gamma.hh"
68#include "G4MuonMinus.hh"
69#include "G4MuonPlus.hh"
70#include "G4NeutrinoMu.hh"
71#include "G4AntiNeutrinoMu.hh"
72#include "G4Lambda.hh"
73
74#include "G4VShortLivedParticle.hh"
75#include "G4ShortLivedConstructor.hh"
76#include "G4ParticleTable.hh"
77#include "G4ShortLivedTable.hh"
78
79#include "G4KineticTrack.hh"
80#include "G4KineticTrackVector.hh"
81#include "G4KineticTrackVector.hh"
82
83#include "G4ParticleDefinition.hh"
84
85#include "G4VCollision.hh"
86#include "G4CollisionNN.hh"
87#include "G4CollisionNNElastic.hh"
88
89int main()
90{
91 // MGP ---- HBOOK initialization
92 HepTupleManager* hbookManager;
93 hbookManager = new HBookFile("collision.hbook", 58);
94 assert (hbookManager != 0);
95
96 // MGP ---- Book histograms
97
98 // MGP ---- Book a ntuple
99 HepTuple* ntuple;
100 ntuple = hbookManager->ntuple("Collision ntuple");
101 assert (ntuple != 0);
102
103
104 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
105
106 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
107 G4ParticleDefinition* antiProton = G4AntiProton::AntiProtonDefinition();
108 G4ParticleDefinition* neutron = G4Neutron::NeutronDefinition();
109 G4ParticleDefinition* antiNeutron = G4AntiNeutron::AntiNeutronDefinition();
110
111 G4ParticleDefinition* pionPlus = G4PionPlus::PionPlusDefinition();
112 G4ParticleDefinition* pionMinus = G4PionMinus::PionMinusDefinition();
113 G4ParticleDefinition* pionZero = G4PionZero::PionZeroDefinition();
114
115 G4ParticleDefinition* kaonPlus = G4KaonPlus::KaonPlusDefinition();
116 G4ParticleDefinition* kaonMinus = G4KaonMinus::KaonMinusDefinition();
117
118 G4ParticleDefinition* lambda = G4Lambda::LambdaDefinition();
119
120 G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
121 G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
122
123 G4ParticleDefinition* theNeutrinoMu = G4NeutrinoMu::NeutrinoMuDefinition();
124 G4ParticleDefinition* theAntiNeutrinoMu = G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
125
126 G4ShortLivedConstructor ShortLived;
127 ShortLived.ConstructParticle();
128
129 // Select track 1
130
131 G4cout << "---- Set KineticTrack 1 ----" << G4endl;
132 G4int id1 = 0;
133 G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- " << G4endl;
134 G4cin >> id1;
135
136 G4ParticleDefinition* definition1;
137 switch (id1)
138 {
139 case 1:
140 definition1 = proton;
141 break;
142 case 2:
143 definition1 = antiProton;
144 break;
145 case 3:
146 definition1 = neutron;
147 break;
148 case 4:
149 definition1 = antiNeutron;
150 break;
151 case 5:
152 definition1 = pionPlus;
153 break;
154 case 6:
155 definition1 = pionMinus;
156 break;
157 case 7:
158 definition1 = kaonPlus;
159 break;
160 case 8:
161 definition1 = kaonMinus;
162 break;
163 default:
164 definition1 = proton;
165 break;
166 }
167 G4double mass1 = definition1->GetPDGMass();
168
169 // Formation time
170 G4double time1 = 0.;
171 // G4cout << "Formation time (ns)" << G4endl;
172 // G4cin >> time1;
173 // if (time1 < 0.) time1 = 0.;
174 // time1 = time1 * ns;
175
176 // Initial position
177 G4double x1 = 0.;
178 G4double y1 = 0.;
179 G4double z1 = 0.;
180 // G4cout <<"Initial position (cm) " << G4endl;
181 // G4cin >> x1 >> y1 >> z1;
182 G4ThreeVector position1(x1*cm, y1*cm, z1*cm);
183
184 // Momentum
185 G4double pzMin1;
186 G4double pzMax1;
187 G4cout << "Pz min and max (GeV)" << G4endl;
188 G4cin >> pzMin1 >> pzMax1;
189 pzMin1 = pzMin1 * GeV;
190 pzMax1 = pzMax1 * GeV;
191
192 // Select track 2
193
194 G4cout << "---- Set KineticTrack 2 ----" << G4endl;
195 G4int id2 = 0;
196 G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- " << G4endl;
197 G4cin >> id2;
198
199 G4ParticleDefinition* definition2;
200 G4KineticTrack trk2;
201 switch (id2)
202 {
203 case 1:
204 definition2 = proton;
205 break;
206 case 2:
207 definition2 = antiProton;
208 break;
209 case 3:
210 definition2 = neutron;
211 break;
212 case 4:
213 definition2 = antiNeutron;
214 break;
215 case 5:
216 definition2 = pionPlus;
217 break;
218 case 6:
219 definition2 = pionMinus;
220 break;
221 case 7:
222 definition2 = kaonPlus;
223 break;
224 case 8:
225 definition2 = kaonMinus;
226 break;
227 default:
228 definition2 = proton;
229 break;
230 }
231 G4double mass2 = definition2->GetPDGMass();
232
233 // Formation time
234 G4double time2 = 0.;
235 // G4cout << "Formation time (ns)" << G4endl;
236 // G4cin >> time2;
237 // if (time2 < 0.) time2 = 0.;
238 // time2 = time2 * ns;
239
240 // Initial position
241 G4double x2 = 0.;
242 G4double y2 = 0.;
243 G4double z2 = 0.;
244 // G4cout <<"Initial position (cm) " << G4endl;
245 // G4cin >> x2 >> y2 >> z2;
246 G4ThreeVector position2(x2*cm, y2*cm, z2*cm);
247
248 // Momentum
249 G4double pzMin2;
250 G4double pzMax2;
251 G4cout << "Pz min and max (GeV)" << G4endl;
252 G4cin >> pzMin2 >> pzMax2;
253 pzMin2 = pzMin2 * GeV;
254 pzMax2 = pzMax2 * GeV;
255
256 G4VCollision* collision = new G4CollisionNNElastic;
257 // collision->Print();
258
259 // Number of iterations
260 G4int n = 0;
261 G4cout << "---- Number of iterations ----" << G4endl;
262 G4cin >> n;
263 if (n < 1) n = 1;
264
265 G4double step1 = (pzMax1 - pzMin1) / n;
266 G4double step2 = (pzMax2 - pzMin2) / n;
267
268 G4cout << "pz min max step "
269 << pzMin1 / GeV << ", "
270 << pzMax1 / GeV << ", "
271 << step1 / GeV << ", " << G4endl
272 << "pz min max step "
273 << pzMin2 / GeV << ", "
274 << pzMax2 / GeV << ", "
275 << step2 / GeV << G4endl;
276
277 G4int i;
278 for (i=0; i<n; i++)
279 {
280 G4double pStep1 = pzMin1 + step1 * i;
281 G4ThreeVector p1(0.,0.,pStep1);
282 G4double e1 = sqrt(mass1*mass1 + p1.mag()*p1.mag());
283 G4LorentzVector p41(p1,e1);
284 G4KineticTrack trk1(definition1,time1,position1,p41);
285
286 G4double pStep2 = -(pzMin2 + step2 * i);
287 G4ThreeVector p2(0.,0.,pStep2);
288 G4double e2 = sqrt(mass2*mass2 + p2.mag()*p2.mag());
289 G4LorentzVector p42(p2,e2);
290 G4KineticTrack trk2(definition2,time2,position2,p42);
291
292 G4double sqrts = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
293
294 G4KineticTrackVector* finalState = collision->FinalState(trk1,trk2);
295 G4int nFinal = finalState->size();
296 for (G4int it=0; it<nFinal; it++)
297 {
298 G4KineticTrack* trk = (*finalState)[it];
299 G4double mass = trk->GetDefinition()->GetPDGMass();
300 G4ThreeVector p3 = trk->Get4Momentum().vect();
301 G4double p = p3.mag() / GeV;
302 G4double px = p3.x() / GeV;
303 G4double py = p3.y() / GeV;
304 G4double pz = p3.z() / GeV;
305 G4cout << mass
306 << " - p = "
307 << p
308 << ", (px,py,pz) = "
309 << px
310 << ","
311 << py
312 << ","
313 << pz
314 << G4endl;
315 }
316 }
317
318
319 hbookManager->write();
320
321 return EXIT_SUCCESS;
322}
323
324
325
326
327
Note: See TracBrowser for help on using the repository browser.