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

Last change on this file since 1199 was 1199, checked in by garnier, 15 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
106int 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(x1*cm, y1*cm, 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(x2*cm, y2*cm, 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(mass1*mass1 + 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(mass2*mass2 + 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.