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

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

nvx fichiers dans CVS
File size: 9.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// $Id: G4DecayTest.cc,v 1.2 2003/06/19 14:39:22 gunter Exp $
24// -------------------------------------------------------------------
25//
26// -------------------------------------------------------------------
27// GEANT 4 class file --- Copyright CERN 1998
28// CERN Geneva Switzerland
29//
30//
31// File name: G4CrossSectionTest.cc
32//
33// Author: Maria Grazia Pia (pia@genova.infn.it),
34//
35// Creation date: 15 April 1999
36//
37// Modifications:
38//
39// -------------------------------------------------------------------
40
41#include "globals.hh"
42
43#include "G4ios.hh"
44#include <fstream>
45#include <iomanip>
46#include <iostream>
47#include <assert.h>
48
49#include "CLHEP/Hist/TupleManager.h"
50#include "CLHEP/Hist/HBookFile.h"
51#include "CLHEP/Hist/Histogram.h"
52#include "CLHEP/Hist/Tuple.h"
53
54#include "Randomize.hh"
55
56#include "G4ThreeVector.hh"
57#include "G4LorentzVector.hh"
58#include "G4LorentzRotation.hh"
59
60#include "G4AntiProton.hh"
61#include "G4Proton.hh"
62#include "G4Neutron.hh"
63#include "G4AntiNeutron.hh"
64#include "G4KaonPlus.hh"
65#include "G4KaonMinus.hh"
66#include "G4PionPlus.hh"
67#include "G4PionMinus.hh"
68#include "G4PionZero.hh"
69#include "G4Gamma.hh"
70#include "G4MuonMinus.hh"
71#include "G4MuonPlus.hh"
72#include "G4NeutrinoMu.hh"
73#include "G4AntiNeutrinoMu.hh"
74#include "G4Lambda.hh"
75#include "G4LambdacPlus.hh"
76
77
78#include "G4VDecayChannel.hh"
79#include "G4DecayTable.hh"
80
81#include "G4KineticTrack.hh"
82#include "G4KineticTrackVector.hh"
83
84#include "G4VShortLivedParticle.hh"
85#include "G4ShortLivedConstructor.hh"
86#include "G4ParticleTable.hh"
87#include "G4ShortLivedTable.hh"
88
89#include "G4ParticleTable.hh"
90#include "G4VCrossSectionSource.hh"
91#include "G4XAqmTotal.hh"
92#include "G4XAqmElastic.hh"
93#include "G4XNNTotalLowE.hh"
94#include "G4XPDGTotal.hh"
95#include "G4XNNTotal.hh"
96#include "G4XPDGElastic.hh"
97
98
99int main()
100{
101 // MGP ---- HBOOK initialization
102 HepTupleManager* hbookManager;
103 hbookManager = new HBookFile("decay.hbook", 58);
104 assert (hbookManager != 0);
105
106 // MGP ---- Book histograms
107
108 HepHistogram* hSigma;
109
110 hSigma = hbookManager->histogram("CrossSection", 100,0.,10.);
111 assert (hSigma != 0);
112
113 // MGP ---- Book a ntuple
114 HepTuple* ntuple;
115 ntuple = hbookManager->ntuple("Decay ntuple");
116 assert (ntuple != 0);
117
118 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
119
120 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
121 G4ParticleDefinition* antiProton = G4AntiProton::AntiProtonDefinition();
122 G4ParticleDefinition* neutron = G4Neutron::NeutronDefinition();
123 G4ParticleDefinition* antiNeutron = G4AntiNeutron::AntiNeutronDefinition();
124
125 G4ParticleDefinition* pionPlus = G4PionPlus::PionPlusDefinition();
126 G4ParticleDefinition* pionMinus = G4PionMinus::PionMinusDefinition();
127 G4ParticleDefinition* pionZero = G4PionZero::PionZeroDefinition();
128
129 G4ParticleDefinition* kaonPlus = G4KaonPlus::KaonPlusDefinition();
130 G4ParticleDefinition* kaonMinus = G4KaonMinus::KaonMinusDefinition();
131
132 G4ParticleDefinition* lambda = G4Lambda::LambdaDefinition();
133
134 G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
135 G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
136
137 G4ParticleDefinition* theNeutrinoMu = G4NeutrinoMu::NeutrinoMuDefinition();
138 G4ParticleDefinition* theAntiNeutrinoMu = G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
139
140 G4ParticleDefinition* lambdacPlus = G4LambdacPlus::LambdacPlusDefinition();
141
142 // G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
143 // G4ParticleDefinition* deltaPlus = particleTable->FindParticle("delta+");
144 // G4ParticleDefinition* deltaPlus = particleTable->FindParticle("N(1440)0");
145
146 G4ShortLivedConstructor ShortLived;
147 G4ShortLivedTable ShortLivedTab;
148 ShortLived.ConstructParticle();
149
150 G4int DEBUG;
151 cout << "insert Debug value 0/1" << G4endl;
152 G4cin >> DEBUG;
153
154 // Select track 1
155
156 G4cout << "---- Set KineticTrack 1 ----" << G4endl;
157 G4int id1 = 0;
158 G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- " << G4endl;
159 G4cin >> id1;
160
161 G4ParticleDefinition* definition1;
162 switch (id1)
163 {
164 case 1:
165 definition1 = proton;
166 break;
167 case 2:
168 definition1 = antiProton;
169 break;
170 case 3:
171 definition1 = neutron;
172 break;
173 case 4:
174 definition1 = antiNeutron;
175 break;
176 case 5:
177 definition1 = pionPlus;
178 break;
179 case 6:
180 definition1 = pionMinus;
181 break;
182 case 7:
183 definition1 = kaonPlus;
184 break;
185 case 8:
186 definition1 = kaonMinus;
187 break;
188 // case 9:
189 // definition1 = deltaPlus;
190 // break;
191 default:
192 definition1 = proton;
193 break;
194 }
195 if(DEBUG==1)cout << definition1 << G4endl;
196 G4double mass1 = definition1->GetPDGMass();
197 G4String name1;
198 if(DEBUG==1)name1 = definition1->GetParticleName();
199 if(DEBUG==1)cout << "Particle Type: " << name1 << G4endl;
200 if(DEBUG==1)cout << "Particle Mass: " << mass1 << G4endl;
201
202
203 // Formation time
204 G4double time1 = 0.;
205
206 // Initial position
207 G4double x1 = 0.;
208 G4double y1 = 0.;
209 G4double z1 = 0.;
210 G4ThreeVector pos1(x1*cm, y1*cm, z1*cm);
211
212 // Momentum
213 G4double px1 = 0. * GeV;
214 G4double py1 = 0. * GeV;
215 G4double pz1 = 0. * GeV;
216 cout << " enter Pz (Px=Py=0.)" << G4endl;
217 G4cin >> pz1;
218 G4ThreeVector p1(px1,py1,pz1);
219 G4double e1 = sqrt(mass1*mass1 + p1.mag()*p1.mag());
220 G4LorentzVector p41(p1,e1);
221
222 G4KineticTrack mytrack;
223
224 G4int iter=0, iter_max;
225 cout << " Enter number of iteration" << G4endl;
226 G4cin >> iter_max;
227 while(iter < iter_max){
228 cout << " =================> Processing event:......" << (iter+1) << G4endl;
229 mytrack=G4KineticTrack(definition1,time1,pos1,p41);
230 if(DEBUG==1)cout << mytrack.GetDefinition()->GetParticleName() << G4endl;
231 // mytrack.SetFormationTime(time1);
232 // mytrack.SetPosition(pos1);
233 // mytrack.Set4Momentum(p41);
234
235 // G4DecayTable Dectab=G4DecayTable();
236 // definition1->SetDecayTable(&Dectab);
237 if(DEBUG==1){
238 cout << "\n " << G4endl;
239 definition1->DumpTable();
240 cout << "\n " << G4endl;}
241 G4DecayTable* TheDec=definition1->GetDecayTable();
242 // cout << " n decay channels: " << TheDec << G4endl;
243 // if(TheDec == NULL)cout << " **************** " << G4endl;
244
245 // cout << " Table not yet dumped " << G4endl;
246 // TheDec->DumpInfo();
247 // cout << " Table dumped " << G4endl;
248
249 G4KineticTrackVector* Dec = mytrack.Decay();
250 if(Dec != NULL){
251// if(DEBUG==1)TheDec->DumpInfo();
252 if(DEBUG==1)cout << " Pointer to Decay Products: "<< Dec << G4endl;
253 G4int nchan = TheDec->entries();
254 /* G4double TheBR[10];
255 cout << nchan << G4endl;
256 for(int index=nchan-1; index >=0; index--){
257 TheBR[index]=TheDec->GetDecayChannel(index)->GetBR();
258 cout << " BR: " << TheBR[index] << G4endl;} */
259 }
260
261 // get kinetictracks from kinetictracksvector
262
263 G4int nComponents = 0;
264 if (Dec != 0)
265 nComponents = Dec->entries();
266 else
267 continue;
268 G4int i;
269 G4LorentzVector Sumv=G4LorentzVector();
270 G4double pxx1=-999.;
271 G4double pyy1=-999.;
272 G4double pzz1=-999.;
273 G4double pxx2=-999.;
274 G4double pyy2=-999.;
275 G4double pzz2=-999.;
276 G4double pxx3=-999.;
277 G4double pyy3=-999.;
278 G4double pzz3=-999.;
279 for (i=0; i<nComponents; i++)
280 {
281 if(DEBUG==1)G4cout << "* Component " << i << ": " <<
282 Dec->at(i)->GetDefinition()->GetParticleName() << " " <<
283 Dec->at(i)->GetPosition() << " " <<
284 Dec->at(i)->Get4Momentum() << " " <<
285 G4endl;
286
287 switch (i) {
288 case 0:
289 pxx1=Dec->at(i)->Get4Momentum().px();
290 pyy1=Dec->at(i)->Get4Momentum().py();
291 pzz1=Dec->at(i)->Get4Momentum().pz();
292 break;
293 case 1:
294 pxx2=Dec->at(i)->Get4Momentum().px();
295 pyy2=Dec->at(i)->Get4Momentum().py();
296 pzz2=Dec->at(i)->Get4Momentum().pz();
297 break;
298 case 2:
299 pxx3=Dec->at(i)->Get4Momentum().px();
300 pyy3=Dec->at(i)->Get4Momentum().py();
301 pzz3=Dec->at(i)->Get4Momentum().pz();
302 break;}
303 Sumv=Sumv+Dec->at(i)->Get4Momentum();
304 }
305 G4double mass=Sumv.m()/1000.;
306 if(DEBUG==1)cout << "Mass:" << mass << G4endl;
307
308 ntuple->column("mass",mass);
309 ntuple->column("px1",pxx1);
310 ntuple->column("px2",pxx2);
311 ntuple->column("px3",pxx3);
312 ntuple->column("py1",pyy1);
313 ntuple->column("py2",pyy2);
314 ntuple->column("py3",pyy3);
315 ntuple->column("pz1",pzz1);
316 ntuple->column("pz2",pzz2);
317 ntuple->column("pz3",pzz3);
318 // ntuple->column("id",i);
319 ntuple->dumpData();
320
321 iter++;
322 int k;
323 }
324 hbookManager->write();
325
326 // delete source;
327
328 return EXIT_SUCCESS;
329}
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
Note: See TracBrowser for help on using the repository browser.