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

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

nvx fichiers dans CVS
File size: 9.7 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: G4BrTest.cc,v 1.2 2003/06/19 14:39:18 gunter Exp $
24// -------------------------------------------------------------------
25//
26// -------------------------------------------------------------------
27// GEANT 4 class file --- Copyright CERN 1998
28// CERN Geneva Switzerland
29//
30//
31// File name: G4DecayTest.cc
32//
33// Author: Lorenzo Bellagamba (bellagamba@bo.infn.it),
34//
35// Creation date: 30 October 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 "G4Eta.hh"
76#include "G4LambdacPlus.hh"
77
78
79#include "G4VDecayChannel.hh"
80#include "G4DecayTable.hh"
81
82#include "G4KineticTrack.hh"
83#include "G4KineticTrackVector.hh"
84
85#include "G4VShortLivedParticle.hh"
86#include "G4ShortLivedConstructor.hh"
87#include "G4ParticleTable.hh"
88#include "G4ShortLivedTable.hh"
89
90#include "G4ParticleTable.hh"
91#include "G4VCrossSectionSource.hh"
92#include "G4XAqmTotal.hh"
93#include "G4XAqmElastic.hh"
94#include "G4XNNTotalLowE.hh"
95#include "G4XPDGTotal.hh"
96#include "G4XNNTotal.hh"
97#include "G4XPDGElastic.hh"
98
99
100int main()
101{
102 // MGP ---- HBOOK initialization
103 HepTupleManager* hbookManager;
104 hbookManager = new HBookFile("br.hbook", 58);
105 assert (hbookManager != 0);
106
107 // MGP ---- Book histograms
108
109 // MGP ---- Book a ntuple
110 HepTuple* ntuple;
111 ntuple = hbookManager->ntuple("Decay ntuple");
112 assert (ntuple != 0);
113
114 G4ShortLivedConstructor ShortLived;
115 G4ShortLivedTable ShortLivedTab;
116 ShortLived.ConstructParticle();
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 G4ParticleDefinition* eta = G4Eta::EtaDefinition();
134
135 G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
136 G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
137
138 G4ParticleDefinition* theNeutrinoMu = G4NeutrinoMu::NeutrinoMuDefinition();
139 G4ParticleDefinition* theAntiNeutrinoMu = G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
140
141 // G4ParticleDefinition* lambdacPlus = G4LambdacPlus::LambdacPlusDefinition();
142
143 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
144 G4ParticleDefinition* N1440_Plus = particleTable->FindParticle("N(1440)+");
145 G4ParticleDefinition* N1440_0 = particleTable->FindParticle("N(1440)0");
146 G4ParticleDefinition* N1535_0 = particleTable->FindParticle("N(1535)0");
147 // G4ParticleDefinition* deltaPluslus = particleTable->FindParticle("delta++");
148 G4ParticleDefinition* deltaPlus = particleTable->FindParticle("delta+");
149 G4ParticleDefinition* deltaMinus = particleTable->FindParticle("delta-");
150
151
152 G4int DEBUG;
153 cout << "insert Debug value 0/1" << G4endl;
154 G4cin >> DEBUG;
155
156 // Select track 1
157
158 G4cout << "---- Set KineticTrack 1 ----" << G4endl;
159 G4int id1 = 0;
160 G4cout << "1 p 2 ap 3 n 4 an 5 pi+ 6 pi- 7 K+ 8 K- 9 N(1440)0 10 N(1440)+ 11 delta+ 12 N(1535)0 13 delta++" << G4endl;
161 G4cin >> id1;
162
163 G4ParticleDefinition* definition1;
164 switch (id1)
165 {
166 case 1:
167 definition1 = proton;
168 break;
169 case 2:
170 definition1 = antiProton;
171 break;
172 case 3:
173 definition1 = neutron;
174 break;
175 case 4:
176 definition1 = antiNeutron;
177 break;
178 case 5:
179 definition1 = pionPlus;
180 break;
181 case 6:
182 definition1 = pionMinus;
183 break;
184 case 7:
185 definition1 = kaonPlus;
186 break;
187 case 8:
188 definition1 = kaonMinus;
189 break;
190 case 9:
191 definition1 = N1440_0;
192 break;
193 case 10:
194 definition1 = N1440_Plus;
195 break;
196 case 11:
197 definition1 = deltaPlus;
198 break;
199 case 12:
200 definition1 = N1535_0;
201 break;
202 // case 13:
203 // definition1 = deltaPlusPlus;
204 // break;
205 default:
206 definition1 = proton;
207 break;
208 }
209 if(DEBUG==1)cout << definition1 << G4endl;
210 G4double mass0 = definition1->GetPDGMass();
211 G4double width0 = definition1->GetPDGWidth();
212 G4String name0;
213 if(DEBUG==1)name0 = definition1->GetParticleName();
214 if(DEBUG==1)cout << "Particle Type: " << name0 << G4endl;
215 if(DEBUG==1)cout << "Particle Mass: " << mass0 << G4endl;
216 if(DEBUG==1)cout << "Particle Width: " << width0 << G4endl;
217
218 G4KineticTrack wig;
219 cout << "Enter number of width to be consdered for BW mass distribtion" << G4endl;
220 cout << " If Nwidth = 0 Mass = PDGMass " << G4endl;
221 double nwid;
222 G4cin >> nwid;
223
224 // Formation time
225 G4double time1 = 0.;
226
227 // Initial position
228 G4double x1 = 0.;
229 G4double y1 = 0.;
230 G4double z1 = 0.;
231 G4ThreeVector pos1(x1*cm, y1*cm, z1*cm);
232
233 // Momentum
234 G4double px1 = 0. * GeV;
235 G4double py1 = 0. * GeV;
236 G4double pz1 = 0. * GeV;
237 cout << " enter Pz (Px=Py=0.)" << G4endl;
238 G4cin >> pz1;
239
240 G4int iter=0, iter_max;
241 G4int nchan;
242 cout << " Enter number of iteration" << G4endl;
243 G4cin >> iter_max;
244 G4double Dch[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
245 G4double max=wig.BrWig(width0, mass0, mass0);
246 G4double Dmas=(2*nwid*width0)/iter_max;
247 while(iter < iter_max){
248 cout << " =================> Processing event:......" << (iter+1) << G4endl;
249 G4double mass1=(mass0-nwid*width0)+iter*Dmas;
250 cout << mass1 << G4endl;
251 G4double width1=width0;
252
253 ntuple->column("mass",mass1);
254
255 G4ThreeVector p1(px1,py1,pz1);
256 G4double e1 = sqrt(mass1*mass1 + p1.mag()*p1.mag());
257 G4LorentzVector p41(p1,e1);
258
259 G4KineticTrack mytrack(definition1,time1,pos1,p41);
260
261 ntuple->column("bwmas",mytrack.BrWig(width0,mass0,mass1));
262
263 if(DEBUG==1)cout << mytrack.GetDefinition()->GetParticleName() << G4endl;
264 G4DecayTable* TheDec=definition1->GetDecayTable();
265 nchan=TheDec->entries();
266 G4double* br=mytrack.GetActualWidth();
267 G4int nn;
268 G4double TotWd=0.;
269 // Total Width
270 for (nn=0; nn < nchan; nn++)TotWd = TotWd + br[nn];
271 ntuple->column("totwd",TotWd);
272 for (nn=0; nn < nchan; nn++){
273 // fill decay channel histos
274 switch (nn){
275 case 0:
276 ntuple->column("wd0",br[nn]);
277 ntuple->column("br0",br[nn]/TotWd);
278 break;
279 case 1:
280 ntuple->column("wd1",br[nn]);
281 ntuple->column("br1",br[nn]/TotWd);
282 break;
283 case 2:
284 ntuple->column("wd2",br[nn]);
285 ntuple->column("br2",br[nn]/TotWd);
286 break;
287 case 3:
288 ntuple->column("wd3",br[nn]);
289 ntuple->column("br3",br[nn]/TotWd);
290 break;
291 case 4:
292 ntuple->column("wd4",br[nn]);
293 ntuple->column("br4",br[nn]/TotWd);
294 break;
295 case 5:
296 ntuple->column("wd5",br[nn]);
297 ntuple->column("br5",br[nn]/TotWd);
298 break;
299 case 6:
300 ntuple->column("wd6",br[nn]);
301 ntuple->column("br6",br[nn]/TotWd);
302 break;
303 case 7:
304 ntuple->column("wd7",br[nn]);
305 ntuple->column("br7",br[nn]/TotWd);
306 break;
307 case 8:
308 ntuple->column("wd8",br[nn]);
309 ntuple->column("br8",br[nn]/TotWd);
310 break;
311 case 9:
312 ntuple->column("wd9",br[nn]);
313 ntuple->column("br9",br[nn]/TotWd);
314 break;
315 case 10:
316 ntuple->column("wd10",br[nn]);
317 ntuple->column("br10",br[nn]/TotWd);
318 break;
319 }
320 }
321 // ntuple->column("id",i);
322 ntuple->dumpData();
323
324 iter++;
325 }
326 hbookManager->write();
327
328 // delete source;
329
330 return EXIT_SUCCESS;
331}
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
Note: See TracBrowser for help on using the repository browser.