source: trunk/source/physics_lists/lists/src/G4HadronInelasticQLHEP.cc @ 836

Last change on this file since 836 was 825, checked in by garnier, 16 years ago

import all except CVS
File size: 15.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// $Id: G4HadronInelasticQLHEP.cc,v 1.1 2006/10/31 11:35:10 gunter Exp $
27// GEANT4 tag $Name:  $
28//
29//---------------------------------------------------------------------------
30//
31// ClassName:   G4HadronInelasticQLHEP
32//
33// Author: 11 April 2006 V. Ivanchenko
34//
35// Modified:
36// 05.07.2006 V.Ivanchenko fix problem of initialisation of HP
37//
38//----------------------------------------------------------------------------
39//
40
41#include "G4HadronInelasticQLHEP.hh"
42
43#include "G4HadronInelasticProcess.hh"
44
45#include "G4ParticleDefinition.hh"
46#include "G4ProcessManager.hh"
47
48#include "G4MesonConstructor.hh"
49#include "G4BaryonConstructor.hh"
50
51#include "G4PiNuclearCrossSection.hh"
52
53#include "G4TheoFSGenerator.hh"
54#include "G4GeneratorPrecompoundInterface.hh"
55#include "G4PreCompoundModel.hh"
56#include "G4ExcitationHandler.hh"
57#include "G4QGSMFragmentation.hh"
58#include "G4ExcitedStringDecay.hh"
59
60#include "G4CascadeInterface.hh"
61#include "G4BinaryCascade.hh"
62#include "G4HadronFissionProcess.hh"
63#include "G4HadronCaptureProcess.hh"
64#include "G4LFission.hh"
65#include "G4LCapture.hh"
66
67#include "G4NeutronHPInelastic.hh"
68#include "G4NeutronHPFission.hh"
69#include "G4NeutronHPCapture.hh"
70
71#include "G4LEAntiLambdaInelastic.hh"
72#include "G4LEAntiNeutronInelastic.hh"
73#include "G4LEAntiOmegaMinusInelastic.hh"
74#include "G4LEAntiProtonInelastic.hh"
75#include "G4LEAntiSigmaMinusInelastic.hh"
76#include "G4LEAntiSigmaPlusInelastic.hh"
77#include "G4LEAntiXiMinusInelastic.hh"
78#include "G4LEAntiXiZeroInelastic.hh"
79#include "G4LEKaonMinusInelastic.hh"
80#include "G4LEKaonPlusInelastic.hh"
81#include "G4LEKaonZeroSInelastic.hh"
82#include "G4LEKaonZeroLInelastic.hh"
83#include "G4LENeutronInelastic.hh"
84#include "G4LELambdaInelastic.hh"
85#include "G4LEProtonInelastic.hh"
86#include "G4LEPionPlusInelastic.hh"
87#include "G4LEPionMinusInelastic.hh"
88#include "G4LEOmegaMinusInelastic.hh"
89#include "G4LESigmaMinusInelastic.hh"
90#include "G4LESigmaPlusInelastic.hh"
91#include "G4LEXiMinusInelastic.hh"
92#include "G4LEXiZeroInelastic.hh"
93
94#include "G4HEAntiLambdaInelastic.hh"
95#include "G4HEAntiNeutronInelastic.hh"
96#include "G4HEAntiOmegaMinusInelastic.hh"
97#include "G4HEAntiProtonInelastic.hh"
98#include "G4HEAntiSigmaMinusInelastic.hh"
99#include "G4HEAntiSigmaPlusInelastic.hh"
100#include "G4HEAntiXiMinusInelastic.hh"
101#include "G4HEAntiXiZeroInelastic.hh"
102#include "G4HEKaonMinusInelastic.hh"
103#include "G4HEKaonPlusInelastic.hh"
104#include "G4HEKaonZeroInelastic.hh"
105#include "G4HENeutronInelastic.hh"
106#include "G4HELambdaInelastic.hh"
107#include "G4HEProtonInelastic.hh"
108#include "G4HEPionPlusInelastic.hh"
109#include "G4HEPionMinusInelastic.hh"
110#include "G4HEOmegaMinusInelastic.hh"
111#include "G4HESigmaMinusInelastic.hh"
112#include "G4HESigmaPlusInelastic.hh"
113#include "G4HEXiMinusInelastic.hh"
114#include "G4HEXiZeroInelastic.hh"
115
116#include "G4HadronProcessStore.hh"
117
118G4HadronInelasticQLHEP::G4HadronInelasticQLHEP(const G4String& name, 
119    G4int ver, G4bool qgs, G4bool bert, G4bool bic, G4bool hp)
120  : G4VPhysicsConstructor(name), verbose(ver), qgsFlag(qgs), 
121    bertFlag(bert), bicFlag(bic), hpFlag(hp), wasActivated(false)
122{
123  if(verbose > 1) G4cout << "### HadronInelasticQLHEP" << G4endl;
124  store = G4HadronProcessStore::Instance();
125  theCascade = 0;
126  theQGStringDecay = 0;
127  theQGStringModel = 0;
128  thePreEquilib = 0;
129  theHPXSecI = 0;
130  theHPXSecC = 0;
131  theHPXSecF = 0;
132}
133
134G4HadronInelasticQLHEP::~G4HadronInelasticQLHEP()
135{
136  if(wasActivated) {
137    delete theCascade;
138    delete theQGStringDecay;
139    delete theQGStringModel;
140    delete thePreEquilib;
141    delete theHPXSecI;
142    delete theHPXSecC;
143    delete theHPXSecF;
144  }
145}
146
147void G4HadronInelasticQLHEP::ConstructParticle()
148{
149  G4MesonConstructor pMesonConstructor;
150  pMesonConstructor.ConstructParticle();
151
152  G4BaryonConstructor pBaryonConstructor;
153  pBaryonConstructor.ConstructParticle();
154}
155
156void G4HadronInelasticQLHEP::ConstructProcess()
157{
158  if(wasActivated) return;
159  wasActivated = true;
160
161  if(verbose > 1) G4cout << "### HadronInelasticQLHEP Construct Process" << G4endl;
162
163  G4double minEneutron   = 0.0*GeV;
164  G4double minELEP       = 0.0*GeV;
165  G4double maxELEP       = 55.*GeV;
166  G4double minEstring    = 12.*GeV;
167  if(hpFlag) minEneutron = 19.5*MeV;
168
169  //Bertini
170  G4HadronicInteraction* theBERT = 0;
171  if(bertFlag) {
172    minELEP = 9.5*GeV;
173    theBERT = new G4CascadeInterface();
174    theBERT->SetMinEnergy(0.0);
175    theBERT->SetMaxEnergy(9.9*GeV);
176  }
177
178  //Binari
179  G4HadronicInteraction* theBIC = 0;
180  if(bicFlag) {
181    minELEP = 9.5*GeV;
182    theBIC = new G4BinaryCascade();
183    theBIC->SetMinEnergy(0.0);
184    theBIC->SetMaxEnergy(9.9*GeV);
185  }
186
187  //QGSP
188  G4TheoFSGenerator* theQGSModel = 0;
189  if(qgsFlag) {
190    maxELEP     = 25.*GeV;
191    theQGSModel = new G4TheoFSGenerator();
192    theQGStringModel  = new G4QGSModel< G4QGSParticipants >;
193    theQGStringDecay  = new G4ExcitedStringDecay(new G4QGSMFragmentation());
194    theQGStringModel->SetFragmentationModel(theQGStringDecay);
195    theCascade = new G4GeneratorPrecompoundInterface;
196    thePreEquilib = new G4PreCompoundModel(new G4ExcitationHandler);
197    theCascade->SetDeExcitation(thePreEquilib);
198    theQGSModel->SetTransport(theCascade);
199    theQGSModel->SetHighEnergyGenerator(theQGStringModel);
200    theQGSModel->SetMinEnergy(minEstring);
201    theQGSModel->SetMaxEnergy(100*TeV);
202  }
203
204  theParticleIterator->reset();
205  while( (*theParticleIterator)() ) {
206    G4ParticleDefinition* particle = theParticleIterator->value();
207    G4String pname = particle->GetParticleName();
208    if(verbose > 1) G4cout << "### HadronInelasticQLHEP:  " << pname << G4endl;
209    if(pname == "anti_lambda"  ||
210       pname == "anti_neutron" ||
211       pname == "anti_omega-"  || 
212       pname == "anti_proton"  || 
213       pname == "anti_sigma-"  || 
214       pname == "anti_sigma+"  || 
215       pname == "anti_xi-"  || 
216       pname == "anti_xi0"  || 
217       pname == "kaon-"     || 
218       pname == "kaon+"     || 
219       pname == "kaon0S"    || 
220       pname == "kaon0L"    || 
221       pname == "lambda"    || 
222       pname == "neutron"   || 
223       pname == "omega-"    || 
224       pname == "pi-"       || 
225       pname == "pi+"       || 
226       pname == "proton"    || 
227       pname == "sigma-"    || 
228       pname == "sigma+"    || 
229       pname == "xi-"       || 
230       pname == "xi0") {
231
232      G4ProcessManager* pmanager = particle->GetProcessManager();
233      G4HadronInelasticProcess* hp = 
234        new G4HadronInelasticProcess("hInelastic", particle);
235      pmanager->AddDiscreteProcess(hp);
236
237      if(pname == "proton") {
238        if(qgsFlag) {
239          Register(particle,hp,theQGSModel,"QGS");
240          hp->AddDataSet(&theXSecP);
241        } else {
242          AddHEP(particle, hp, 25.*GeV, 100.*TeV);
243        }
244        AddLEP(particle, hp, minELEP, maxELEP);
245       
246        if(bicFlag)       Register(particle,hp,theBIC,"Binary");
247        else if(bertFlag) Register(particle,hp,theBERT,"Bertini");
248
249      } else if(pname == "neutron") {
250        if(qgsFlag) {
251          Register(particle,hp,theQGSModel,"QGS");
252          hp->AddDataSet(&theXSecN);
253        } else {
254          AddHEP(particle, hp, 25.*GeV, 100.*TeV);
255        }
256        AddLEP(particle, hp, minELEP, maxELEP);
257
258        G4HadronCaptureProcess* theNeutronCapture = 
259          new G4HadronCaptureProcess("nCapture");
260        G4HadronFissionProcess* theNeutronFission = 
261          new G4HadronFissionProcess("nFission");
262        pmanager->AddDiscreteProcess(theNeutronCapture);
263        pmanager->AddDiscreteProcess(theNeutronFission);
264
265        if(hpFlag) {
266          theHPXSecI = new G4NeutronHPInelasticData;
267          theHPXSecC = new G4NeutronHPCaptureData;
268          theHPXSecF = new G4NeutronHPFissionData;
269          hp->AddDataSet(theHPXSecI);
270          theNeutronCapture->AddDataSet(theHPXSecC);
271          theNeutronFission->AddDataSet(theHPXSecF);
272          G4NeutronHPInelastic* hpi = new G4NeutronHPInelastic();
273          G4NeutronHPCapture* hpc = new G4NeutronHPCapture();
274          G4NeutronHPFission* hpf = new G4NeutronHPFission();
275          Register(particle,hp,hpi,"HP");
276          Register(particle,theNeutronCapture,hpc,"HP");
277          Register(particle,theNeutronFission,hpf,"HP");
278        }
279
280        G4HadronicInteraction* theC = new G4LCapture();
281        theC->SetMinEnergy(minEneutron);
282        Register(particle,theNeutronCapture,theC,"LCapture");
283
284        G4HadronicInteraction* theF = new G4LFission();
285        theF->SetMinEnergy(minEneutron);
286        Register(particle,theNeutronFission,theF,"LFission");
287
288        if(bicFlag) {
289          G4BinaryCascade* theB = new G4BinaryCascade();
290          theB->SetMinEnergy(minEneutron);
291          theB->SetMaxEnergy(9.9*GeV);
292          Register(particle,hp,theB,"Binary");
293        } else if(bertFlag) {
294          G4CascadeInterface* theB = new G4CascadeInterface();
295          theB->SetMinEnergy(minEneutron);
296          theB->SetMaxEnergy(9.9*GeV);
297          Register(particle,hp,theB,"Bertini");
298        }
299
300      } else if(pname == "pi-" || pname == "pi+") {
301        if(qgsFlag) {
302          Register(particle,hp,theQGSModel,"QGS");
303          hp->AddDataSet(&thePiCross);
304        } else {
305          AddHEP(particle, hp, 25.*GeV, 100.*TeV);
306        }
307        AddLEP(particle, hp, minELEP, maxELEP);
308        if(bertFlag) Register(particle,hp,theBERT,"Bertini");
309
310      } else if(pname == "kaon-"     || 
311                pname == "kaon+"     || 
312                pname == "kaon0S"    || 
313                pname == "kaon0L") {
314
315        if(qgsFlag) Register(particle,hp,theQGSModel,"QGS");
316        else        AddHEP(particle, hp, 25.*GeV, 100.*TeV);
317       
318        AddLEP(particle, hp, minELEP, maxELEP);
319        if(bertFlag) Register(particle,hp,theBERT,"Bertini");
320
321      } else {
322
323        AddHEP(particle, hp, 25.*GeV, 100.*TeV);
324        AddLEP(particle, hp, 0.0, 55.*GeV);
325      }
326
327      if(verbose > 1)
328        G4cout << "### HadronInelasticQLHEP: " << hp->GetProcessName()
329               << " added for " << pname << G4endl;
330    }
331  }
332  store->Dump(verbose);
333}
334
335void G4HadronInelasticQLHEP::AddLEP(G4ParticleDefinition* particle,
336                                    G4HadronicProcess* hp,
337                                    G4double emin,
338                                    G4double emax)
339{
340  G4HadronicInteraction* hi = 0;
341  G4String pname = particle->GetParticleName();
342
343  if(pname == "anti_lambda"      ) hi = new G4LEAntiLambdaInelastic();
344  else if(pname == "anti_neutron") hi = new G4LEAntiNeutronInelastic();
345  else if(pname == "anti_omega-" ) hi = new G4LEAntiOmegaMinusInelastic();
346  else if(pname == "anti_proton" ) hi = new G4LEAntiProtonInelastic();
347  else if(pname == "anti_sigma-" ) hi = new G4LEAntiSigmaMinusInelastic();
348  else if(pname == "anti_sigma+" ) hi = new G4LEAntiSigmaPlusInelastic();
349  else if(pname == "anti_xi-"    ) hi = new G4LEAntiXiMinusInelastic();
350  else if(pname == "anti_xi0"    ) hi = new G4LEAntiXiZeroInelastic();
351  else if(pname == "kaon-"       ) hi = new G4LEKaonMinusInelastic();
352  else if(pname == "kaon+"       ) hi = new G4LEKaonPlusInelastic();
353  else if(pname == "kaon0S"      ) hi = new G4LEKaonZeroSInelastic();
354  else if(pname == "kaon0L"      ) hi = new G4LEKaonZeroLInelastic();
355  else if(pname == "neutron"     ) hi = new G4LENeutronInelastic();
356  else if(pname == "lambda"      ) hi = new G4LELambdaInelastic();
357  else if(pname == "omega-"      ) hi = new G4LEOmegaMinusInelastic();
358  else if(pname == "proton"      ) hi = new G4LEProtonInelastic();
359  else if(pname == "pi+"         ) hi = new G4LEPionPlusInelastic();
360  else if(pname == "pi-"         ) hi = new G4LEPionMinusInelastic();
361  else if(pname == "sigma-"      ) hi = new G4LESigmaMinusInelastic();
362  else if(pname == "sigma+"      ) hi = new G4LESigmaPlusInelastic();
363  else if(pname == "xi-"         ) hi = new G4LEXiMinusInelastic();
364  else if(pname == "xi0"         ) hi = new G4LEXiZeroInelastic();
365
366  if(hi) {
367    hi->SetMinEnergy(emin);
368    hi->SetMaxEnergy(emax);
369    Register(particle,hp,hi,"LHEP");
370  } else {
371    G4cout << "### G4HadronInelasticTHEO: ERROR - no LHEP model for "
372           << pname << G4endl;
373  }
374}
375
376
377void G4HadronInelasticQLHEP::AddHEP(G4ParticleDefinition* particle,
378                                    G4HadronicProcess* hp,
379                                    G4double emin,
380                                    G4double emax)
381{
382  G4HadronicInteraction* hi = 0;
383  G4String pname = particle->GetParticleName();
384
385  if(pname == "anti_lambda"      ) hi = new G4HEAntiLambdaInelastic();
386  else if(pname == "anti_neutron") hi = new G4HEAntiNeutronInelastic();
387  else if(pname == "anti_omega-" ) hi = new G4HEAntiOmegaMinusInelastic();
388  else if(pname == "anti_proton" ) hi = new G4HEAntiProtonInelastic();
389  else if(pname == "anti_sigma-" ) hi = new G4HEAntiSigmaMinusInelastic();
390  else if(pname == "anti_sigma+" ) hi = new G4HEAntiSigmaPlusInelastic();
391  else if(pname == "anti_xi-"    ) hi = new G4HEAntiXiMinusInelastic();
392  else if(pname == "anti_xi0"    ) hi = new G4HEAntiXiZeroInelastic();
393  else if(pname == "kaon-"       ) hi = new G4HEKaonMinusInelastic();
394  else if(pname == "kaon+"       ) hi = new G4HEKaonPlusInelastic();
395  else if(pname == "kaon0S"      ) hi = new G4HEKaonZeroInelastic();
396  else if(pname == "kaon0L"      ) hi = new G4HEKaonZeroInelastic();
397  else if(pname == "neutron"     ) hi = new G4HENeutronInelastic();
398  else if(pname == "lambda"      ) hi = new G4HELambdaInelastic();
399  else if(pname == "omega-"      ) hi = new G4HEOmegaMinusInelastic();
400  else if(pname == "proton"      ) hi = new G4HEProtonInelastic();
401  else if(pname == "pi+"         ) hi = new G4HEPionPlusInelastic();
402  else if(pname == "pi-"         ) hi = new G4HEPionMinusInelastic();
403  else if(pname == "sigma-"      ) hi = new G4HESigmaMinusInelastic();
404  else if(pname == "sigma+"      ) hi = new G4HESigmaPlusInelastic();
405  else if(pname == "xi-"         ) hi = new G4HEXiMinusInelastic();
406  else if(pname == "xi0"         ) hi = new G4HEXiZeroInelastic();
407
408  if(hi) {
409    hi->SetMinEnergy(emin);
410    hi->SetMaxEnergy(emax);
411    Register(particle,hp,hi,"HEP");
412  } else {
413    G4cout << "### G4HadronInelasticQLHEP: ERROR - no HEP model for "
414           << pname << G4endl;
415  }
416}
417
418void G4HadronInelasticQLHEP::Register(G4ParticleDefinition* p, 
419                                      G4HadronicProcess* hp, 
420                                      G4HadronicInteraction* hi, 
421                                      const G4String& m)
422{
423  hp->RegisterMe(hi);
424  store->Register(hp,p,hi,m);
425  if(verbose > 1)
426    G4cout << "### QLHEP: Register new model " << m
427           << " for " << p->GetParticleName() << " and " 
428           << hp->GetProcessName()
429           << " E(GeV) " << hi->GetMinEnergy()/GeV
430           << " - " << hi->GetMaxEnergy()/GeV << G4endl;
431}
Note: See TracBrowser for help on using the repository browser.